1C*GRIMG3 -- gray-scale map of a 2D data array, using dither
2C+
3      SUBROUTINE GRIMG3 (A, IDIM, JDIM, I1, I2, J1, J2,
4     1                   BLACK, WHITE, PA, MODE)
5      INTEGER IDIM, JDIM, I1, I2, J1, J2, MODE
6      REAL    A(IDIM,JDIM)
7      REAL    BLACK, WHITE
8      REAL    PA(6)
9C--
10C 2-Sep-1994 - moved from GRGRAY [TJP].
11C-----------------------------------------------------------------------
12      INCLUDE 'grpckg1.inc'
13      INTEGER  I,IX,IX1,IX2,IY,IY1,IY2,J
14      REAL     DEN,VALUE,BW
15      REAL     XXAA,XXBB,YYAA,YYBB,XYAA,XYBB,YXAA,YXBB,XYAAIY,YXAAIY
16      INTEGER  M, IAA, ICC, JRAN, ILAST, JLAST, IXSTEP, IYSTEP
17      REAL     RAND, RM, FAC, FACL
18      PARAMETER (M=714025, IAA=1366, ICC=150889, RM=1.0/M)
19      PARAMETER (FAC=65000.0)
20      INTRINSIC MOD, NINT, REAL, LOG
21C-----------------------------------------------------------------------
22C
23      IF (MODE.LT.0 .OR. MODE.GT.2) RETURN
24C
25C Initialize random-number generator (based on RAN2 of Press et al.,
26C Numerical Recipes)
27C
28      JRAN = 76773
29C
30      IX1 = NINT(GRXMIN(GRCIDE))+1
31      IX2 = NINT(GRXMAX(GRCIDE))-1
32      IY1 = NINT(GRYMIN(GRCIDE))+1
33      IY2 = NINT(GRYMAX(GRCIDE))-1
34      DEN = PA(2)*PA(6)-PA(3)*PA(5)
35C
36C Calculate constants.
37C
38      BW   = ABS(BLACK-WHITE)
39      FACL = LOG(1.0+FAC)
40      XXAA = (-PA(6))*PA(1)/DEN
41      XXBB = PA(6)/DEN
42      XYAA = (-PA(3))*PA(4)/DEN
43      XYBB = PA(3)/DEN
44      YYAA = (-PA(2))*PA(4)/DEN
45      YYBB = PA(2)/DEN
46      YXAA = (-PA(5))*PA(1)/DEN
47      YXBB = PA(5)/DEN
48C
49C Choose step size: at least 1/200 inch, assuming the line-width
50C unit is 1/200 inch.
51C
52      IXSTEP = MAX(1,NINT(GRWIDT(GRCIDE)*GRPXPI(GRCIDE)/200.0))
53      IYSTEP = MAX(1,NINT(GRWIDT(GRCIDE)*GRPYPI(GRCIDE)/200.0))
54C
55C Draw dots.
56C
57      ILAST = 0
58      JLAST = 0
59      DO 120 IY=IY1,IY2,IYSTEP
60          XYAAIY = XXAA-XYAA-XYBB*IY
61          YXAAIY = YYAA+YYBB*IY-YXAA
62          DO 110 IX=IX1,IX2,IXSTEP
63              I = NINT(XYAAIY+XXBB*IX)
64              IF (I.LT.I1.OR.I.GT.I2) GOTO 110
65              J = NINT(YXAAIY-YXBB*IX)
66              IF (J.LT.J1.OR.J.GT.J2) GOTO 110
67              IF (I.NE.ILAST .OR. J.NE.JLAST) THEN
68                  ILAST = I
69                  JLAST = J
70                  VALUE = ABS(A(I,J)-WHITE)/BW
71                  IF (MODE.EQ.0) THEN
72C                     -- "linear"
73                      CONTINUE
74                  ELSE IF (MODE.EQ.1) THEN
75C                     -- "logarithmic"
76                      VALUE = LOG(1.0+FAC*VALUE)/FACL
77                  ELSE IF (MODE.EQ.2) THEN
78C                     -- "square root"
79                      VALUE = SQRT(VALUE)
80                  END IF
81              END IF
82              JRAN = MOD(JRAN*IAA+ICC, M)
83              RAND = JRAN*RM
84              IF (VALUE.GT.RAND) CALL GRDOT0(REAL(IX),REAL(IY))
85  110     CONTINUE
86  120  CONTINUE
87C-----------------------------------------------------------------------
88       END
89