!___________________________________________________________________________________ 
!This module provides some useful high-level 2D and 3D Fortran 90 plotting routines 
!based on the DISLIN package. 
!written by Aleksandar Donev (donev@pa.msu.edu), last modifed 11/2/2000 
!NOTE: This file is not heavily commented. For more information, contact me. 
!Feel free to use all or any part of this program. 
!___________________________________________________________________________________ 

MODULE SimpleGraphics
 USE DISLIN
 IMPLICIT NONE

 PUBLIC :: InitGraphics,AddLegend,Plot2D,Plot3D,SurfPlot,EndGraphics
 PRIVATE

 CONTAINS
 !___________________________________________________________________________________
 ! This subroutine must be called before plotting begins
 ! The function my_func can perform any additional level 1 routines needed 
 !___________________________________________________________________________________
 SUBROUTINE InitGraphics(file,file_type,plot_title,x_label,y_label,z_label,my_func)
   CHARACTER(*), OPTIONAL, INTENT(IN) :: file,file_type
   CHARACTER(*), OPTIONAL, INTENT(IN) :: x_label,y_label,z_label
   CHARACTER(*), DIMENSION(:), OPTIONAL, INTENT(IN) :: plot_title
   INTEGER, EXTERNAL, OPTIONAL :: my_func
   
   INTEGER :: my_result,i
 
   IF(PRESENT(file)) CALL SETFIL(file)
   IF(PRESENT(file_type)) CALL METAFL(file_type)
 
   CALL SCRMOD("REVERS")
   IF(PRESENT(my_func)) my_result=my_func(0) ! Level 0 additional routines
   CALL DISINI
   CALL PAGERA
   CALL COMPLX
   CALL NOCHEK
   CALL MIXALF    
   CALL SETMIX('_',"IND")
   CALL SETMIX('^',"EXP")
   CALL SETMIX('$',"RES")
   CALL SETMIX('/',"LEG")    
    
   IF(PRESENT(x_label)) CALL NAME(x_label,'X') 
   IF(PRESENT(y_label)) CALL NAME(y_label,'Y') 
   IF(PRESENT(z_label)) CALL NAME(z_label,'Z')
   IF(PRESENT(plot_title)) THEN
     DO i=1,SIZE(plot_title) 
       CALL TITLIN(plot_title(i),i) 
     END DO 
   END IF     
   IF(PRESENT(my_func)) my_result=my_func(1) ! Level 1 additional routines 
 END SUBROUTINE InitGraphics 
 
 !___________________________________________________________________________________  
 ! This subroutine will add a legend to the plot at the given position (see LEGEND) 
 ! The legend should be stored in an array of strings of the right size 
 ! It should be called AFTER a plotting routine 
 !___________________________________________________________________________________ 
  SUBROUTINE AddLegend(legends,position) 
   CHARACTER(*), DIMENSION(:), INTENT(IN) :: legends 
   INTEGER, OPTIONAL, INTENT(IN) :: position 
   CHARACTER(10*LEN(legends(1))) :: Leg ! For some reason SIZE does not work  
   INTEGER :: i,n,length,leg_pos 
   n=SIZE(legends)
   length=LEN(legends(1)) 
      
   CALL LEGINI(Leg,n,length) 
   CALL FRAME(0) 
   CALL MIXLEG 
   CALL LEGCLR 
   CALL LEGTIT("") 
   DO i=1,n      
     CALL LEGLIN(Leg,legends(i),i) 
   END DO   
      
   IF(PRESENT(position)) THEN  
     leg_pos=position 
   ELSE 
     leg_pos=3 
   END IF   
   CALL LEGEND(Leg,leg_pos)
   
 END SUBROUTINE AddLegend             

 !___________________________________________________________________________________ 
 ! This subroutine plots 2D curves on a graph 
 ! with specified line-types, symbol-types and colors. 
 ! ALL arguments are optional
 ! This time the x and y data are both stored in rank-one arrays-vectors!  
 ! The logical indicator new indicates whether this is a new plot or an old one.
 ! Again, my_func can perform any additional level 2 routines (after GRAF) needed. 
 ! If specified, axis determines the scaling. 
 ! Otherwise, the scaling is chosen so as to fit the data range of x and y.
 !___________________________________________________________________________________    
 SUBROUTINE Plot2D(x,y,plot_spec,new_plot,my_func,axis) 
 
   REAL, DIMENSION(:), INTENT(IN), OPTIONAL :: x 
   REAL, DIMENSION(SIZE(x)), INTENT(IN), OPTIONAL :: y 
   CHARACTER(LEN=3), OPTIONAL, INTENT(IN) :: plot_spec 
   LOGICAL, OPTIONAL, INTENT(IN) :: new_plot
   INTEGER, EXTERNAL, OPTIONAL :: my_func 
   REAL, DIMENSION(4), OPTIONAL, INTENT(IN) :: axis 
 
   INTEGER :: i,my_result,sym_type 
   REAL :: min_x,max_x,min_y,max_y,len_x,len_y 
   CHARACTER(LEN=1) :: type_spec,line_spec,color_spec  
   LOGICAL :: plot_new
   
   IF(PRESENT(new_plot)) THEN
     plot_new=new_plot
   ELSE
     plot_new=.TRUE.
   END IF  

   IF(plot_new) THEN   
     ! Find the scaling from the maximum range of values in the data 
     ! or simply use axis if present: 
     IF(.NOT.PRESENT(axis)) THEN 
       IF(PRESENT(x).AND.PRESENT(y)) THEN
		       min_x=MINVAL(x) 
		       max_x=MAXVAL(x) 
		       min_y=MINVAL(y) 
		       max_y=MAXVAL(y) 
		     ELSE
		       min_x=0.0
		       min_y=0.0
		       max_x=1.0
		       max_y=1.0
		     END IF    
     ELSE 
       min_x=axis(1) 
       max_x=axis(2) 
       min_y=axis(3) 
       max_y=axis(4) 
     END IF 
     len_x=max_x-min_x 
     len_y=max_y-min_y 
     CALL GRAF(min_x-len_x/20,max_x+len_x/20,min_x,len_x/10,min_y-len_y/20,max_y+len_y/20,min_y,len_y/10) 
     CALL TITLE
   END IF
      
   IF(PRESENT(my_func)) my_result=my_func(2) 
 
   ! This big chunk sets the appropriate symbol or line type and color 
   ! and then plots the corresponding curve.
   ! We don't change the settings if this is not a new plot and there is no specification
   ! Therefore, if this is not a new plot and there is no specification, DISLIN defaults 
   ! will be used, not mine. So be wise!
   IF(plot_new. OR. (.NOT.plot_new.AND.PRESENT(plot_spec))) THEN
				   IF(PRESENT(plot_spec)) THEN
				     type_spec=plot_spec(1:1)
				     line_spec=plot_spec(2:2)
				     color_spec=plot_spec(3:3)
				   ELSE
				     type_spec='S'
				     line_spec='C'
				     color_spec='B'
				   END IF       
       SELECT CASE(type_spec) 
         CASE('S') 
           CALL INCMRK(-1) 
           SELECT CASE(line_spec) 
             CASE('S') 
               sym_type=0 
             CASE('C') 
               sym_type=1 
             CASE('T') 
               sym_type=2 
             CASE('D') 
               sym_type=5 
             CASE('s') 
               sym_type=16 
             CASE('c') 
               sym_type=21 
             CASE('t') 
               sym_type=20 
             CASE('d') 
               sym_type=19 
             CASE DEFAULT 
               sym_type=1 
           END SELECT 
         CALL MARKER(sym_type) 
         CASE('L') 
           CALL INCMRK(0) 
           SELECT CASE(line_spec) 
             CASE('-') 
               CALL SOLID 
             CASE('.') 
               CALL DOTL 
             CASE(':') 
               CALL DASHL 
             CASE('|') 
               CALL CHNDOT 
             CASE DEFAULT 
               CALL DOTL 
           END SELECT 
         CASE DEFAULT 
           CALL INCMRK(1) 
           CALL MARKER(1) 
           CALL DOTL 
       END SELECT 
       CALL SETVLT("SMALL") 
       SELECT CASE(color_spec) 
         CASE('R') 
           CALL SETCLR(2) 
         CASE('B') 
           CALL SETCLR(4) 
         CASE('G') 
           CALL SETCLR(3) 
         CASE('Y') 
           CALL SETCLR(5) 
         CASE DEFAULT 
           CALL SETCLR(2) 
       END SELECT 
   END IF

   IF(PRESENT(x).AND.PRESENT(y)) CALL CURVE(x,y,SIZE(x))
 
 END SUBROUTINE Plot2D  

  SUBROUTINE Plot3D(x,y,z,plot_spec,new_plot,my_func,axis,view) 
 
   REAL, DIMENSION(:), INTENT(IN), OPTIONAL :: x 
   REAL, DIMENSION(SIZE(x)), INTENT(IN), OPTIONAL :: y,z 
   CHARACTER(LEN=3), OPTIONAL, INTENT(IN) :: plot_spec 
   LOGICAL, OPTIONAL, INTENT(IN) :: new_plot
   INTEGER, EXTERNAL, OPTIONAL :: my_func 
   REAL, DIMENSION(6), OPTIONAL, INTENT(IN) :: axis 
   REAL, DIMENSION(3), OPTIONAL :: view   
 
   INTEGER :: i,my_result,sym_type 
   REAL :: min_x,max_x,min_y,max_y,min_z,max_z,len_x,len_y,len_z
   CHARACTER(LEN=1) :: type_spec,line_spec,color_spec  
   LOGICAL :: plot_new
   
   IF(PRESENT(new_plot)) THEN
     plot_new=new_plot
   ELSE
     plot_new=.TRUE.
   END IF  

   IF(plot_new) THEN   
     ! Find the scaling from the maximum range of values in the data 
     ! or simply use axis if present: 
     IF(.NOT.PRESENT(axis)) THEN 
       IF(PRESENT(x).AND.PRESENT(y)) THEN
		       min_x=MINVAL(x) 
		       max_x=MAXVAL(x) 
		       min_y=MINVAL(y) 
		       max_y=MAXVAL(y)
		       min_z=MINVAL(z) 
		       max_z=MAXVAL(z)		        
		     ELSE
		       min_x=0.0
		       min_y=0.0
		       min_z=1.0
		       max_x=1.0
		       max_y=1.0
		       max_z=1.0
		     END IF    
     ELSE 
       min_x=axis(1) 
       max_x=axis(2) 
       min_y=axis(3) 
       max_y=axis(4) 
       min_z=axis(5) 
       max_z=axis(6)       
     END IF 
     len_x=max_x-min_x 
     len_y=max_y-min_y
     len_z=max_z-min_z 
    IF(PRESENT(view)) CALL VIEW3D(view(1),view(2),view(3),"ANGLE")
    CALL GRAF3D(min_x-len_x/20,max_x+len_x/20,min_x,len_x/10,min_y-len_y/20,max_y+len_y/20,min_y,len_y/10, &
       min_z-len_z/20,max_z+len_z/20,min_z,len_z/10)
     CALL TITLE
   END IF
      
   IF(PRESENT(my_func)) my_result=my_func(2) 
 
   ! This big chunk sets the appropriate symbol or line type and color 
   ! and then plots the corresponding curve.
   ! We don't change the settings if this is not a new plot and there is no specification
   ! Therefore, if this is not a new plot and there is no specification, DISLIN defaults 
   ! will be used, not mine. So be wise!
   IF(plot_new. OR. (.NOT.plot_new.AND.PRESENT(plot_spec))) THEN
				   IF(PRESENT(plot_spec)) THEN
				     type_spec=plot_spec(1:1)
				     line_spec=plot_spec(2:2)
				     color_spec=plot_spec(3:3)
				   ELSE
				     type_spec='S'
				     line_spec='C'
				     color_spec='B'
				   END IF       
       SELECT CASE(type_spec) 
         CASE('S') 
           CALL INCMRK(-1) 
           SELECT CASE(line_spec) 
             CASE('S') 
               sym_type=0 
             CASE('C') 
               sym_type=1 
             CASE('T') 
               sym_type=2 
             CASE('D') 
               sym_type=5 
             CASE('s') 
               sym_type=16 
             CASE('c') 
               sym_type=21 
             CASE('t') 
               sym_type=20 
             CASE('d') 
               sym_type=19 
             CASE DEFAULT 
               sym_type=1 
           END SELECT 
         CALL MARKER(sym_type) 
         CASE('L') 
           CALL INCMRK(0) 
           SELECT CASE(line_spec) 
             CASE('-') 
               CALL SOLID 
             CASE('.') 
               CALL DOTL 
             CASE(':') 
               CALL DASHL 
             CASE('|') 
               CALL CHNDOT 
             CASE DEFAULT 
               CALL DOTL 
           END SELECT 
         CASE DEFAULT 
           CALL INCMRK(1) 
           CALL MARKER(1) 
           CALL DOTL 
       END SELECT 
       CALL SETVLT("SMALL") 
       SELECT CASE(color_spec) 
         CASE('R') 
           CALL SETCLR(2) 
         CASE('B') 
           CALL SETCLR(4) 
         CASE('G') 
           CALL SETCLR(3) 
         CASE('Y') 
           CALL SETCLR(5) 
         CASE DEFAULT 
           CALL SETCLR(2) 
       END SELECT 
   END IF

   IF(PRESENT(x).AND.PRESENT(y).AND.PRESENT(z)) CALL CURV3D(x,y,z,SIZE(x))
 
 END SUBROUTINE Plot3D 
   
 !___________________________________________________________________________________   
 ! This subroutine plots a 3D or 2D surface (color) plot from given data 
 ! The function should not be called more than once--this crowds the plots
 ! and is rarely useful, so it is disabled.
 ! The data is given on a regular grid in x and y (array vectors)
 ! in the form of a rank-2 matrix z=f(x,y).
 ! my_func can perform any additional level 3 routines (after GRAF3) needed. 
 ! The color_range for the color pallette can be given explicitly.
 ! If specified, axis determines the scaling. 
 ! Otherwise, the scaling is chosen so as to fit the data range of x and y. 
 ! The view point (camera) can also be given in the form of (theta,phi,r).
 !___________________________________________________________________________________     
 SUBROUTINE SurfPlot(x,y,z,plot_spec,my_func,color_range,axis,view,contours) 
   REAL, DIMENSION(:), INTENT(IN) :: x,y 
   REAL, DIMENSION(SIZE(X),SIZE(Y)), INTENT(IN) :: z 
   CHARACTER(3), OPTIONAL, INTENT(IN) :: plot_spec 
   INTEGER, EXTERNAL, OPTIONAL :: my_func 
   REAL, DIMENSION(2), OPTIONAL, INTENT(IN) :: color_range 
   REAL, DIMENSION(6), OPTIONAL, INTENT(IN) :: axis 
   REAL, DIMENSION(3), OPTIONAL, INTENT(IN) :: view 
   REAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: contours
 
   INTEGER :: my_result, this_color,n_labels 
   INTEGER, DIMENSION(:), ALLOCATABLE :: these_labels
   REAL :: min_x,max_x,min_y,max_y,min_z,max_z,len_x,len_y,len_z
   REAL :: this_label, this_range(2) 
   INTEGER :: i,sym_type 
   CHARACTER(3) :: this_spec="3RS" 
    
   ! Again we determine the range of the axis:  
   IF(.NOT.PRESENT(axis)) THEN 
     min_x=MINVAL(x) 
     max_x=MAXVAL(x) 
     min_y=MINVAL(y) 
     max_y=MAXVAL(y) 
     min_z=MINVAL(z) 
     max_z=MAXVAL(z) 
   ELSE 
     min_x=axis(1) 
     max_x=axis(2) 
     min_y=axis(3) 
     max_y=axis(4) 
     min_z=axis(5) 
     max_z=axis(6) 
   END IF 
   len_x=max_x-min_x 
   len_y=max_y-min_y 
   len_z=max_z-min_z 
 
   IF(PRESENT(view)) CALL VIEW3D(view(1),view(2),view(3),"ANGLE") 
   IF(PRESENT(color_range)) THEN
     CALL ZSCALE(color_range(1),color_range(2)) 
     this_range=color_range
   ELSE  
     this_range=(/min_z,max_z/)
   END IF   
   IF(PRESENT(plot_spec)) this_spec=plot_spec 
   IF(PRESENT(contours)) THEN
     n_labels=SIZE(contours)
     ALLOCATE(these_labels(n_labels))
     these_labels=contours
   ELSE
     n_labels=1
     ALLOCATE(these_labels(n_labels))   
     these_labels=(min_z+max_z)/2
   END IF
       
   ! We set the appropriate color pallete: 
   SELECT CASE(this_spec(3:3)) 
    CASE('V') 
       CALL SETVLT("VGA") 
     CASE('G') 
       CALL SETVLT("GREY") 
     CASE('g') 
       CALL SETVLT("RGREY") 
     CASE('r') 
       CALL SETVLT("RRAIN") 
     CASE DEFAULT 
       CALL SETVLT("RAIN") 
   END SELECT 
   SELECT CASE(this_spec(1:1)) 
     ! This part plots a 2D color plot: 
     ! It is highly improvized and WRONG when 
     ! an axis is supplied and they are different from the minmax range of xy
     ! This is because CRVMAT does not accept x and y data:
     CASE('M') 
       ! In this case the second character determines the resolution increase: 
       i=IACHAR(this_spec(2:2))-48 
       IF(.NOT.((i>=1).AND.(i<=9))) i=1 
       CALL AUTRES(SIZE(x),SIZE(y))      
       CALL GRAF3(min_x,max_x,min_x,len_x/10,min_y,max_y,min_y,len_y/10, & 
         min_z-len_z/20,max_z+len_z/20,min_z,len_z/10) 
       CALL TITLE 
       IF(PRESENT(my_func)) my_result=my_func(3) 
       CALL CRVMAT(z,SIZE(x),SIZE(y),i,i)  ! THIS IS WRONG!--use for demos only
     CASE('C')
       SELECT CASE(this_spec(2:2)) 
         CASE('N')
           CALL LABELS("NONE","CONTUR")    
         CASE DEFAULT
           CALL LABELS("FLOAT","CONTUR")     
       END SELECT   
       CALL GRAF3(min_x,max_x,min_x,len_x/10,min_y,max_y,min_y,len_y/10, & 
         this_range(1),this_range(2),this_range(1),(this_range(2)-this_range(1))/10)       
       ! CALL GRAF(min_x-len_x/20,max_x+len_x/20,min_x,len_x/10,min_y-len_y/20,max_y+len_y/20,min_y,len_y/10)
       CALL TITLE 
       IF(PRESENT(my_func)) my_result=my_func(2) 
       DO i=1,n_labels
         this_label=contours(i)
         this_color=MODULO(INT(256*(this_label-this_range(1))/(this_range(2)-this_range(1))),256)
         CALL SETCLR(this_color)
         CALL CONTUR(x,SIZE(x),y,SIZE(y),z,this_label) 
       END DO                        
     ! This part plots a 3D surface plot with the given specifications:         
     CASE DEFAULT 
       CALL GRAF3D(min_x-len_x/20,max_x+len_x/20,min_x,len_x/10,min_y-len_y/20,max_y+len_y/20,min_y,len_y/10, & 
         min_z-len_z/20,max_z+len_z/20,min_z,len_z/10) 
       CALL TITLE 
       IF(PRESENT(my_func)) my_result=my_func(3) 
       SELECT CASE(this_spec(2:2)) 
         CASE('M') 
           CALL SURFCE(x,SIZE(x),y,SIZE(y),z) 
         CASE('C') 
           CALL SHDMOD("SMOOTH","SURFACE") 
           CALL SURSHD(x,SIZE(x),y,SIZE(y),z) 
         CASE DEFAULT 
           CALL SHDMOD("SMOOTH","SURFACE") 
           CALL SURSHD(x,SIZE(x),y,SIZE(y),z) 
           CALL SURFCE(x,SIZE(x),y,SIZE(y),z) 
       END SELECT 
   END SELECT 
   DEALLOCATE(these_labels)

 END SUBROUTINE SurfPlot 
 
 !___________________________________________________________________________________   
 ! This function needs to be called to end the plotting 
 ! my_func now specifies the level 0 routines at the end 
 !___________________________________________________________________________________     
 SUBROUTINE EndGraphics(my_func) 
   INTEGER, EXTERNAL, OPTIONAL :: my_func 
   INTEGER :: my_result 
   IF(PRESENT(my_func)) my_result=my_func(-1) 
   CALL DISFIN 
 END SUBROUTINE EndGraphics 
 
ENDMODULE SimpleGraphics