      PROGRAM DENSTY
************************************************************************
* Program to produce a density plot of a function Z = F(X, Y).  F is   *
* defined by a function subprogram and is passed to subroutine DNPLOT. *
* Identifiers used are:                                                *
*     F          : function to be plotted                              *
*     XMIN, XMAX : minimum and maximum X values                        *
*     YMIN, YMAX : minimum and maximum Y values                        *
*     ZMIN, ZMAX : minimum and maximum Z values                        *
*                                                                      *
* Input:  XMIN, XMAX, YMIN, YMAX                                       *
* Output: User prompts and plot of Z = F(X, Y)                         *
************************************************************************

      REAL F, XMIN, XMAX, YMIN, YMAX, ZMIN, ZMAX
      EXTERNAL F

      PRINT *, 'ENTER MINIMUM AND MAXIMUM X VALUES, THEN Y VALUES'
      READ *, XMIN, XMAX, YMIN, YMAX
      PRINT *, 'ENTER MINIMUM AND MAXIMUM VALUES OF FUNCTION'
      READ *, ZMIN, ZMAX
      CALL DNPLOT(F, XMIN, XMAX, YMIN, YMAX, ZMIN, ZMAX)

      END

    
**F*********************************************************************
*                Function whose graph is to be plotted                 *
************************************************************************

      FUNCTION F(X,Y)

      REAL F, X, Y
      
      F = EXP(-(X**2 + Y**2))
      END

**DNPLOT****************************************************************
* Subroutine to generate a density plot of a function Z = F(X, Y) for  *
* X ranging from XMIN to XMAX and Y ranging from YMIN to YMAX; Z is    *
* allowed to range from ZMIN to ZMAX.  Local identifiers used are:     *
*     FNAME  : name of file containing the output                      *
*     HLIMIT,                                                          *
*     VLIMIT : parameters:  limits on the size of the graphics WINDOW  *
*     WINDOW : two-dimensional character array -- the graphics window  *
*     LGRAY  : parameter: largest index in array GRAY                  *
*     GRAY   : array of symbols representing shades of gray            *
*     DELX   : X increment                                             *
*     DELY   : Y increment                                             *
*     DELZ   : Z increment                                             *
*     X, Y, Z: a point on the graph                                    *
*     XLOC,                                                            *
*     YLOC   : location of a point in the window                       *
*     ZSHADE : shade of gray used to represent height Z                *
*                                                                      *
* Accepts:          F, XMIN, XMAX, YMIN, YMAX, ZMIN, ZMAX              *
* Input (keyboard): FNAME                                              *
* Output (screen):  User prompt                                        *
* Output (file):    The graphics WINDOW                                *
************************************************************************

      SUBROUTINE DNPLOT (F, XMIN, XMAX, YMIN, YMAX,ZMIN, ZMAX)

      REAL F, XMIN, XMAX, YMIN, YMAX, ZMIN, ZMAX, DELX, DELY, DELZ
      INTEGER HLIMIT, VLIMIT, XLOC, YLOC, LGRAY, ZSHADE
      PARAMETER (HLIMIT = 75, VLIMIT = 45, LGRAY = 9)
      CHARACTER*1 WINDOW(0: HLIMIT, 0:VLIMIT), GRAY(0:LGRAY), FNAME*20

*    Initialize array of characters to indicate "densities of gray"
*    that in turn represent heights of function F

      DATA GRAY /'0','1','2','3','4','5','6','7','8','9'/


      PRINT *, 'ENTER NAME OF FILE TO CONTAIN THE DENSITY PLOT'
      READ '(A)', FNAME
      OPEN (UNIT = 20, FILE = FNAME, STATUS = 'NEW')
     
      DELX = (XMAX - XMIN) / REAL(HLIMIT)
      DELY = (YMAX - YMIN) / REAL(VLIMIT)
      DELZ = (ZMAX - ZMIN) / REAL(LGRAY)

* "Shade" each element of WINDOW with appropriate gray

      Y = YMIN
      DO 20 YLOC = 0, VLIMIT
         X = XMIN
         DO 10 XLOC = 0, HLIMIT
            Z = F(X, Y)
              
*           Find gray shade corresponding to Z value

            IF (Z .GE. ZMAX) THEN
               ZSHADE = LGRAYS
            ELSE
               ZSHADE = NINT((Z - ZMIN) / DELZ)
            END IF

            WINDOW(XLOC, YLOC) = GRAY(ZSHADE)
            X = X + DELX
 10      CONTINUE

         Y = Y + DELY
 20   CONTINUE

* Draw the WINDOW in the file

      DO 30 YLOC = VLIMIT, 0, -1
         WRITE(20, *)(WINDOW(XLOC,YLOC), XLOC = 0, HLIMIT)
 30   CONTINUE

      END
