MODULE Function_F
!-------------------------------------------------------------------------
! Module that contains the function whose graph is to be plotted
!-------------------------------------------------------------------------

CONTAINS

  FUNCTION F(X,Y)

    REAL :: F
    REAL, INTENT(IN) :: X, Y

    F = EXP(-(X**2 + Y**2))

  END FUNCTION F

END MODULE Function_F


MODULE GraphingRoutines
!-------------------------------------------------------------------------
! Module that contains the subroutines for producing various kinds
! of graphs:
!    ScatterPlot
!    DensityPlot
!    . . .
!-------------------------------------------------------------------------

CONTAINS

  !-DensityPlot------------------------------------------------------------
  ! Subroutine to generate a density plot of a function Z = F(X, Y) for
  ! X ranging from X_Min to X_Max and Y ranging from Y_Min to Y_Max; Z is
  ! allowed to range from Z_Min to Z_Max.  Local identifiers used are
  !   FileName    : name of file containing the output
  !   Horiz_Limit,
  !   Vert_Limit  : limits on the size of the graphics window (constants)
  !   Window      : two-dimensional character array -- the graphics window
  !   MaxGraySub  : parameter: largest subscript in array Gray
  !   Gray        : array of symbols representing shades of gray
  !   DeltaX      : X increment
  !   DeltaY      : Y increment
  !   DeltaZ      : Z increment
  !   X, Y, Z     : a point on the graph
  !   X_Loc,
  !   Y_Loc       : location of a point in the window
  !   Shade_for_Z : shade of gray used to represent height Z
  !
  ! Accepts:          F, X_Min, X_Max, Y_Min, Y_Max, Z_Min, Z_Max
  ! Input (keyboard): FileName
  ! Output (screen):  User prompt
  ! Output (file):    The graphics window
  !------------------------------------------------------------------------
  
  SUBROUTINE DensityPlot (X_Min, X_Max, Y_Min, Y_Max, Z_Min, Z_Max)
  
    USE Function_F

    IMPLICIT NONE
    REAL, INTENT(IN) ::  X_Min, X_Max, Y_Min, Y_Max, Z_Min, Z_Max
    REAL :: DeltaX, DeltaY, DeltaZ, X, Y, Z
    INTEGER, PARAMETER ::  Horiz_Limit = 75, Vert_Limit = 45, MaxGraySub = 9
    INTEGER :: X_Loc, Y_Loc, Shade_for_Z
    CHARACTER(1), DIMENSION(0: Horiz_Limit, 0:Vert_Limit) :: Window
    CHARACTER(1), DIMENSION(0:MaxGraySub) :: &
                     Gray = (/"0","1","2","3","4","5","6","7","8","9"/)
    CHARACTER(20) :: FileName
  
    WRITE (*, '(1X, A)', ADVANCE = "NO") &
          "Enter name of file to contain the density plot: "
    READ *, FileName
    OPEN (UNIT = 20, FILE = FileName, STATUS = "NEW")
  
    DeltaX = (X_Max - X_Min) / REAL(Horiz_Limit)
    DeltaY = (Y_Max - Y_Min) / REAL(Vert_Limit)
    DeltaZ = (Z_Max - Z_Min) / REAL(MaxGraySub)
  
    ! "Shade" each element of Window with appropriate gray
  
    Y = Y_Min
    DO Y_Loc = 0, Vert_Limit
       X = X_Min
       DO X_Loc = 0, Horiz_Limit
          Z = F(X, Y)
  
          ! Find gray shade corresponding to Z value
  
          IF (Z > Z_Max) THEN
             Shade_for_Z = MaxGraySub
          ELSE
             Shade_for_Z = NINT((Z - Z_Min) / DeltaZ)
          END IF
   
          Window(X_Loc, Y_Loc) = Gray(Shade_for_Z)
          X = X + DeltaX
       END DO
  
       Y = Y + DeltaY
    END DO
  
    ! Draw the Window in the file
  
    DO Y_Loc = Vert_Limit, 0, -1
       WRITE(20, *) (Window(X_Loc,Y_Loc), X_Loc = 0, Horiz_Limit)
    END DO
  
  END SUBROUTINE DensityPlot
  
END MODULE GraphingRoutines


PROGRAM Density_Plot_of_a_Function
!-----------------------------------------------------------------------
! Program to produce a density plot of a function Z = F(X, Y).  F is
! defined by an external function subprogram and is passed to 
! subroutine DensityPlot imported from module GraphingRoutines.  
! Identifiers used are:
!   F            : function to be plotted
!   X_Min, X_Max : minimum and maximum X values
!   Y_Min, Y_Max : minimum and maximum Y values
!   Z_Min, Z_Max : minimum and maximum Z values
!
! Input:  X_Min, X_Max, Y_Min, Y_Max
! Output: User prompts and density plot of Z = F(X, Y)
!-----------------------------------------------------------------------

  USE GraphingRoutines, ONLY: DensityPlot

  IMPLICIT NONE

  INTERFACE
    FUNCTION F(X, Y)
      REAL :: F
      REAL, INTENT(IN) :: X, Y
    END FUNCTION F
  END INTERFACE

  REAL :: X_Min, X_Max, Y_Min, Y_Max, Z_Min, Z_Max

  WRITE (*, '(1X, A)', ADVANCE = "NO") &
        "Enter minimum and maximum x values, then y values: "
  READ *, X_Min, X_Max, Y_Min, Y_Max
  WRITE (*, '(1X, A)', ADVANCE = "NO") &
        "Enter minimum and maximum values of function: "
  READ *, Z_Min, Z_Max
  CALL DensityPlot(X_Min, X_Max, Y_Min, Y_Max, Z_Min, Z_Max)

END PROGRAM Density_Plot_of_a_Function