!___________________________________________________________________________________ 
!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 9/19/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 FunGraphics
  USE DISLIN
  USE SimpleGraphics
  IMPLICIT NONE
  
  PRIVATE
  PUBLIC :: FunPlot2D, FunPlot3D, FunPlotContour, FunVectorPlot2D

CONTAINS

  SUBROUTINE FunPlot2D(f_x,x_range,num_points,plot_spec,axis)  
    
   INTERFACE 
     FUNCTION f_x(x) RESULT(y)
       REAL, INTENT(IN) :: x
       REAL :: y
     END FUNCTION f_x
   END INTERFACE    
   REAL, DIMENSION(2), INTENT(IN) :: x_range
   INTEGER, OPTIONAL, INTENT(IN) :: num_points
   CHARACTER(LEN=3), DIMENSION(:), OPTIONAL, INTENT(IN) :: plot_spec 
   REAL, DIMENSION(4), OPTIONAL, INTENT(IN) :: axis    

   REAL, ALLOCATABLE, DIMENSION(:) :: x, y
   REAL :: dx 
   INTEGER :: i,n_points,n_curves
   CHARACTER(LEN=3), DIMENSION(5) :: this_spec
   LOGICAL :: this_plot

   IF(PRESENT(plot_spec)) THEN
     n_curves=MIN(5,SIZE(plot_spec))
     DO i=1,n_curves
       this_spec(i)=plot_spec(i)
     END DO  
   ELSE
     n_curves=1
     this_spec(1)="STR" ! Just a default 
   END IF  
       
   IF(PRESENT(num_points)) THEN
     n_points=num_points
   ELSE
     n_points=10 ! Default
   END IF

   ALLOCATE(x(n_points),y(n_points))       
   dx=(x_range(2)-x_range(1))/(n_points-1)
   DO i=1,n_points
     x(i)=x_range(1)+(i-1)*dx
     y(i)=f_x(x(i))
   END DO

   this_plot=.TRUE.
   IF(PRESENT(axis)) THEN
     DO i=1,n_curves    
       CALL Plot2D(x=x,y=y,plot_spec=this_spec(i),axis=axis,new_plot=this_plot)
       IF(i==1) this_plot=.FALSE.
     END DO  
   ELSE
     DO i=1,n_curves    
       CALL Plot2D(x=x,y=y,plot_spec=this_spec(i),new_plot=this_plot)
       IF(i==1) this_plot=.FALSE.
     END DO  
   END IF
   DEALLOCATE(x,y)
                 
  END SUBROUTINE FunPlot2D

  SUBROUTINE FunPlot3D(f_xy,xy_range,num_points,plot_spec,axis,view) 
    
   INTERFACE 
     FUNCTION f_xy(x,y) RESULT(z)
       REAL, INTENT(IN) :: x,y
       REAL :: z
     END FUNCTION f_xy
   END INTERFACE
       
   REAL, DIMENSION(4), INTENT(IN) :: xy_range
   INTEGER, DIMENSION(2), OPTIONAL, INTENT(IN) :: num_points
   CHARACTER(LEN=3), OPTIONAL, INTENT(IN) :: plot_spec 
   REAL, DIMENSION(6), OPTIONAL, INTENT(IN) :: axis 
   REAL, DIMENSION(3), OPTIONAL, INTENT(IN) :: view

   REAL, ALLOCATABLE, DIMENSION(:) :: x, y
   REAL, ALLOCATABLE, DIMENSION(:,:) :: z
   REAL :: dx,dy
   INTEGER, DIMENSION(2) :: n_points
   INTEGER :: i,j
   CHARACTER(LEN=3) :: this_spec
   REAL, DIMENSION(6) :: this_axis    

   IF(PRESENT(plot_spec)) THEN
     this_spec=plot_spec
   ELSE
     this_spec="3SR" ! Just a default 
   END IF    
       
   IF(PRESENT(num_points)) THEN
     n_points=num_points
   ELSE
     n_points=(/10,10/) ! Default
   END IF

   ALLOCATE(x(n_points(1)),y(n_points(2)),z(n_points(1),n_points(2)))       
   dx=(xy_range(2)-xy_range(1))/(n_points(1)-1)
   dy=(xy_range(4)-xy_range(3))/(n_points(2)-1)

   x=(/(xy_range(1)+(i-1)*dx,i=1,n_points(1))/) ! Use implied DO
   y=(/(xy_range(3)+(j-1)*dy,j=1,n_points(2))/) ! Implied DO loop
   
   DO i=1,n_points(1)
     DO j=1,n_points(2)
       z(i,j)=f_xy(x(i),y(j))
     END DO
   END DO

   IF(.NOT.PRESENT(axis)) THEN
     IF(PRESENT(view)) THEN
       CALL SurfPlot(x=x,y=y,z=z,plot_spec=this_spec,view=view)
     ELSE
       CALL SurfPlot(x=x,y=y,z=z,plot_spec=this_spec)
     END IF
   ELSE
     IF(PRESENT(view)) THEN
       CALL SurfPlot(x=x,y=y,z=z,plot_spec=this_spec,axis=axis,view=view)
     ELSE
       CALL SurfPlot(x=x,y=y,z=z,plot_spec=this_spec,axis=axis)
     END IF    
  END IF   
  DEALLOCATE(x,y,z) 
                   
  END SUBROUTINE FunPlot3D

  SUBROUTINE FunPlotContour(f_xy,xyz_range,num_points,plot_spec,axis) 
    
   INTERFACE 
     FUNCTION f_xy(x,y) RESULT(z)
       REAL, INTENT(IN) :: x,y
       REAL :: z
     END FUNCTION f_xy
   END INTERFACE
       
   REAL, DIMENSION(6), INTENT(IN) :: xyz_range
   INTEGER, DIMENSION(3), OPTIONAL, INTENT(IN) :: num_points
   CHARACTER(LEN=3), OPTIONAL, INTENT(IN) :: plot_spec 
   REAL, DIMENSION(4), OPTIONAL, INTENT(IN) :: axis 

   REAL, ALLOCATABLE, DIMENSION(:) :: x, y,labels
   REAL, ALLOCATABLE, DIMENSION(:,:) :: z
   REAL :: dx,dy,dl
   INTEGER, DIMENSION(3) :: n_points
   INTEGER :: i,j
   CHARACTER(LEN=3) :: this_spec
   REAL, DIMENSION(4) :: this_axis    

   IF(PRESENT(plot_spec)) THEN
     this_spec=plot_spec
     this_spec(1:1)="C"
   ELSE
     this_spec="CNR" ! Just a default 
   END IF    
       
   IF(PRESENT(num_points)) THEN
     n_points=num_points
   ELSE
     n_points=(/10,10,10/) ! Default
   END IF

   IF(PRESENT(axis)) THEN
     this_axis=axis
   ELSE
     this_axis=(/xyz_range(1:4)/) ! Default
   END IF   

   ALLOCATE(x(n_points(1)),y(n_points(2)),labels(n_points(3)),z(n_points(1),n_points(2)))       
   dx=(xyz_range(2)-xyz_range(1))/(n_points(1)-1)
   dy=(xyz_range(4)-xyz_range(3))/(n_points(2)-1)
   dl=(xyz_range(6)-xyz_range(5))/(n_points(3)-1)

   x=(/(xyz_range(1)+(i-1)*dx,i=1,n_points(1))/) ! Use implied DO
   y=(/(xyz_range(3)+(j-1)*dy,j=1,n_points(2))/) ! Implied DO loop
   labels=(/(xyz_range(5)+(i-1)*dl,i=1,n_points(3))/)
     
   DO i=1,n_points(1)
     DO j=1,n_points(2)
       z(i,j)=f_xy(x(i),y(j))
     END DO
   END DO

   CALL SurfPlot(x=x,y=y,z=z,plot_spec=this_spec,axis=(/this_axis,xyz_range(5:6)/),contours=labels)
     
  DEALLOCATE(x,y,z) 
                   
  END SUBROUTINE FunPlotContour

  SUBROUTINE FunVectorPlot2D(f_xy,xy_range,num_points,plot_spec,axis,zoom) 
    
   INTERFACE 
     FUNCTION f_xy(x,y) RESULT(z)
       REAL, INTENT(IN) :: x,y
       REAL, DIMENSION(2) :: z
     END FUNCTION f_xy
   END INTERFACE
       
   REAL, DIMENSION(4), INTENT(IN) :: xy_range
   INTEGER, DIMENSION(2), OPTIONAL, INTENT(IN) :: num_points
   CHARACTER(LEN=3), OPTIONAL, INTENT(IN) :: plot_spec 
   REAL, DIMENSION(6), OPTIONAL, INTENT(IN) :: axis 
   REAL, OPTIONAL, INTENT(IN) :: zoom

   REAL, ALLOCATABLE, DIMENSION(:) :: x, y
   REAL :: dx,dy, abs_z, scaling
   REAL, DIMENSION(2) :: z
   INTEGER, DIMENSION(2) :: n_points
   INTEGER :: i,j, vec_spec, this_color
   CHARACTER(LEN=3) :: this_spec
   REAL, DIMENSION(6) :: this_axis  
   REAL :: min_x,max_x,min_y,max_y,min_z,max_z,len_x,len_y,len_z  

   IF(PRESENT(zoom)) THEN
     scaling=zoom
   ELSE
     scaling=1.0
   END IF
        
  IF(PRESENT(plot_spec)) THEN
     this_spec=plot_spec
  ELSE
     this_spec="11g" ! Just a default 
  END IF    

  vec_spec=0
  i=IACHAR(this_spec(1:1))-48 
  IF(.NOT.((i>=0).AND.(i<=9))) i=1
  vec_spec=vec_spec+1000*i
  i=IACHAR(this_spec(2:2))-48 
  IF(.NOT.((i>=0).AND.(i<=9))) i=1
  vec_spec=vec_spec+100*i
  vec_spec=vec_spec+21
              
   IF(PRESENT(num_points)) THEN
     n_points=num_points
   ELSE
     n_points=(/10,10/) ! Default
   END IF

   ! Again we determine the range of the axis:  
   IF(.NOT.PRESENT(axis)) THEN 
     min_x=xy_range(1) 
     max_x=xy_range(2)
     min_y=xy_range(3)
     max_y=xy_range(4) 
   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  

   SELECT CASE(this_spec(3:3)) 
    CASE('V') 
       CALL SETVLT("VGA") 
     CASE('G') 
       CALL SETVLT("GREY") 
     CASE('g') 
       CALL SETVLT("RGREY")
       CALL SETCLR(0) 
     CASE('r') 
       CALL SETVLT("RRAIN") 
       CALL SETCLR(0)
     CASE DEFAULT 
       CALL SETVLT("RAIN") 
   END SELECT 
   
   IF(PRESENT(axis)) THEN
      CALL GRAF3(min_x,max_x,min_x,len_x/10,min_y,max_y,min_y,len_y/10, & 
      axis(5),axis(6),axis(5),(axis(6)-axis(5))/10)
   ELSE
     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) 
   END IF  
   CALL TITLE()
              
   ALLOCATE(x(n_points(1)),y(n_points(2)))       
   dx=(xy_range(2)-xy_range(1))/(n_points(1)-1)
   dy=(xy_range(4)-xy_range(3))/(n_points(2)-1)

   x=(/(xy_range(1)+(i-1)*dx,i=1,n_points(1))/) ! Use implied DO
   y=(/(xy_range(3)+(j-1)*dy,j=1,n_points(2))/) ! Implied DO loop

   DO i=1,n_points(1)
     DO j=1,n_points(2)
       z=f_xy(x(i),y(j))
       abs_z=SQRT(z(1)**2+z(2)**2)
       IF(PRESENT(axis)) THEN                            
         IF( (abs_z<MIN(axis(5),axis(6))) .OR. (abs_z>MAX(axis(5),axis(6))) ) CYCLE
         this_color=MODULO(INT(256*ABS((abs_z-axis(5))/(axis(6)-axis(5)))),256)
         CALL SETCLR(this_color)           
       END IF  
       ! Unfortunately, there is no VECTOR routine that accepts
       ! user coordinates, so I use this silly  trick to call FIELD:
       z=scaling*z
       CALL FIELD(x(i:i),y(j:j),x(i:i)+z(1:1),y(j:j)+z(2:2),1,vec_spec)
     END DO
   END DO

  DEALLOCATE(x,y) 
                   
  END SUBROUTINE FunVectorPlot2D
         
ENDMODULE FunGraphics