!___________________________________________________________________________________
! This module provides some useful high-level 2D and 3D Fortran 90 plotting routines
! based on the DISLIN package, specialized for Lahey Fortran 95.
! written by Aleksandar Donev (donev@pa.msu.edu), last modifed 2/26/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 Graphics
  USE DISLIN
  IMPLICIT NONE

  ! An array of pointers 
  TYPE :: array_ptr
    REAL, DIMENSION(:), POINTER :: p
  END TYPE

  PUBLIC :: InitGraphics,AddLegend,Plot2D,Plot3D,Plot2DMatrix,Plot3DMatrix,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 :: file,file_type
    CHARACTER(*), OPTIONAL :: x_label,y_label,z_label
    CHARACTER(*), DIMENSION(:), OPTIONAL :: 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)
    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)
  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,pos)
    CHARACTER(*), DIMENSION(:) :: legends
    INTEGER, OPTIONAL :: pos
    CHARACTER(size(legends)*LEN(legends(1))) :: Leg   
    INTEGER :: i,n,length
    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(.NOT.PRESENT(pos)) THEN 
      CALL LEGEND(Leg,3)
    ELSE
      CALL LEGEND(Leg,pos)
    END IF  
  END SUBROUTINE AddLegend            
  
!___________________________________________________________________________________
  ! This subroutine plots up to 10 2D curves on a graph
  ! with specified line-types, symbol-types and colors.
  ! Again, my_func can perform any additional level 2 routines (after GRAF) needed.
  ! If specified, axis determines the scaling.
  !___________________________________________________________________________________   
  SUBROUTINE Plot2D(x1,y1,spec1,x2,y2,spec2,x3,y3,spec3,x4,y4,spec4,x5,y5,spec5, &
    x6,y6,spec6,x7,y7,spec7,x8,y8,spec8,x9,y9,spec9,x10,y10,spec10,my_func,axis)

    REAL, DIMENSION(:), TARGET, OPTIONAL :: x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,y1,y2,y3,y4,y5,y6,y7,y8,y9,y10
    CHARACTER(3), OPTIONAL :: spec1,spec2,spec3,spec4,spec5,spec6,spec7,spec8,spec9,spec10
    INTEGER, EXTERNAL, OPTIONAL :: my_func
    REAL, DIMENSION(4), OPTIONAL :: axis

    INTEGER :: my_result
    INTEGER, PARAMETER :: N=10
    TYPE(array_ptr), DIMENSION(N) :: x_arrays,y_arrays
    REAL, DIMENSION(N) :: mins_x,mins_y,maxs_x,maxs_y
    REAL :: min_x,max_x,min_y,max_y,len_x,len_y
    INTEGER :: curve_count,i,sym_type
    LOGICAL :: PresentXY,PresentSpec
    CHARACTER(3) :: this_spec
    curve_count=0

    ! Add as many curves to the plot as there are arguments (up to 10)
    ! using array pointers to point to the arrays being plotted
    ! NOTE: This is not very efficient, but is purely Fortan 90 like
    ! and it is easy to understand and modify.
    DO i=1,N
      CALL TestPresenceXY(PresentXY,i)
      IF (PresentXY) THEN
        CALL AddCurve(i)
      END IF
    END DO

    ! Find the scaling from the maximum range of values in the data
    ! or simply use axis if present:
    IF(.NOT.PRESENT(axis)) THEN
      DO i=1,curve_count
        mins_x(i)=MINVAL(x_arrays(i)%p)
        maxs_x(i)=MAXVAL(x_arrays(i)%p)
        mins_y(i)=MINVAL(y_arrays(i)%p)
        maxs_y(i)=MAXVAL(y_arrays(i)%p)
      END DO
      min_x=MINVAL(mins_x(1:curve_count))
      max_x=MAXVAL(maxs_x(1:curve_count))
      min_y=MINVAL(mins_y(1:curve_count))
      max_y=MAXVAL(maxs_y(1:curve_count))
    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
    IF(PRESENT(my_func)) my_result=my_func(2)

    ! This big loop sets the appropriate symbol or line type and color
    ! and then plots the corresponding curve:
    DO i=1,curve_count
      CALL TestPresenceSpec(PresentSpec,i,this_spec)
      IF(PresentSpec) THEN
        SELECT CASE(this_spec(1:1))
          CASE('S')
            CALL INCMRK(-1)
            SELECT CASE(this_spec(2:2))
              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(this_spec(2:2))
              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(this_spec(3:3))
          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
      CALL CURVE(x_arrays(i)%p,y_arrays(i)%p,SIZE(x_arrays(i)%p))
      NULLIFY(x_arrays(i)%p,y_arrays(i)%p)
    END DO

    CONTAINS

    !__________________________________
    ! This internal subroutine tests
    ! to see how many curves were given:
    ! (Fortran 90 does not have varargin)
    !__________________________________ 
    SUBROUTINE TestPresenceXY(TestPresentXY,k) 	  
      LOGICAL :: TestPresentXY
      INTEGER :: k
      SELECT CASE(k)
        CASE(1)
          TestPresentXY=PRESENT(x1).AND.PRESENT(y1)
        CASE(2)
          TestPresentXY=PRESENT(x2).AND.PRESENT(y2)
        CASE(3)
          TestPresentXY=PRESENT(x3).AND.PRESENT(y3)
        CASE(4)
          TestPresentXY=PRESENT(x4).AND.PRESENT(y4)
        CASE(5)
          TestPresentXY=PRESENT(x5).AND.PRESENT(y5)
        CASE(6)
          TestPresentXY=PRESENT(x6).AND.PRESENT(y6)
        CASE(7)
          TestPresentXY=PRESENT(x7).AND.PRESENT(y7)
        CASE(8)
          TestPresentXY=PRESENT(x8).AND.PRESENT(y8)
        CASE(9)
          TestPresentXY=PRESENT(x9).AND.PRESENT(y9)
        CASE(10)
          TestPresentXY=PRESENT(x10).AND.PRESENT(y10)
        CASE DEFAULT
          TestPresentXY=.FALSE.
      END SELECT
    END SUBROUTINE TestPresenceXY

    !__________________________________
    ! Same as the previous routine
    ! but for the format specification
    !__________________________________
    SUBROUTINE TestPresenceSpec(TestPresentSpec,k,Spec)
      LOGICAL :: TestPresentSpec
      INTEGER :: k
      CHARACTER(3) :: Spec
      SELECT CASE(k)
        CASE(1)
          TestPresentSpec=PRESENT(spec1)
          Spec=spec1
        CASE(2)
          TestPresentSpec=PRESENT(spec2)
          Spec=spec2
        CASE(3)
          TestPresentSpec=PRESENT(spec3)
          Spec=spec3
        CASE(4)
          TestPresentSpec=PRESENT(spec4)
          Spec=spec4
        CASE(5)
          TestPresentSpec=PRESENT(spec5)
          Spec=spec5
        CASE(6)
          TestPresentSpec=PRESENT(spec6)
          Spec=spec6
        CASE(7)
          TestPresentSpec=PRESENT(spec7)
          Spec=spec7
        CASE(8)
          TestPresentSpec=PRESENT(spec8)
          Spec=spec8
        CASE(9)
          TestPresentSpec=PRESENT(spec9)
          Spec=spec9
        CASE(10)
          TestPresentSpec=PRESENT(spec10)
          Spec=spec10
        CASE DEFAULT
          TestPresentSpec=.FALSE.
          Spec='SSR'
      END SELECT
    END SUBROUTINE TestPresenceSpec

    !__________________________________
    ! This routine adds a curve to 
    ! the plot "buffer" by assigning an
    ! array pointer to it:
    !__________________________________
    SUBROUTINE AddCurve(k)
      INTEGER :: k
      curve_count=curve_count+1
      SELECT CASE(k)
        CASE(1)
          x_arrays(curve_count)%p=>x1
          y_arrays(curve_count)%p=>y1
        CASE(2)
          x_arrays(curve_count)%p=>x2
          y_arrays(curve_count)%p=>y2
        CASE(3)
          x_arrays(curve_count)%p=>x3
          y_arrays(curve_count)%p=>y3
        CASE(4)
          x_arrays(curve_count)%p=>x4
          y_arrays(curve_count)%p=>y4
        CASE(5)
          x_arrays(curve_count)%p=>x5
          y_arrays(curve_count)%p=>y5
        CASE(6)
          x_arrays(curve_count)%p=>x6
          y_arrays(curve_count)%p=>y6
        CASE(7)
          x_arrays(curve_count)%p=>x7
          y_arrays(curve_count)%p=>y7
        CASE(8)
          x_arrays(curve_count)%p=>x8
          y_arrays(curve_count)%p=>y8
        CASE(9)
          x_arrays(curve_count)%p=>x9
          y_arrays(curve_count)%p=>y9
        CASE(10)
          x_arrays(curve_count)%p=>x10
          y_arrays(curve_count)%p=>y10
      END SELECT
    END SUBROUTINE AddCurve

  END SUBROUTINE Plot2D

  !___________________________________________________________________________________
  ! This subroutine plots 2D curves on a graph
  ! with specified line-types, symbol-types and colors.
  ! This time the data are all stored in the same matrix as COLUMNS!
  ! Again, my_func can perform any additional level 2 routines (after GRAF) needed.
  ! If specified, axis determines the scaling.
  !___________________________________________________________________________________   
  SUBROUTINE Plot2DMatrix(x,y,type_spec,line_spec,color_spec,my_func,axis)

    REAL, DIMENSION(:,:) :: x
    REAL, DIMENSION(SIZE(x,1),SIZE(x,2)) :: y
    CHARACTER(SIZE(x,2)), OPTIONAL :: type_spec,line_spec,color_spec
    INTEGER, EXTERNAL, OPTIONAL :: my_func
    REAL, DIMENSION(4), OPTIONAL :: axis

    INTEGER :: curve_count, n_points
    INTEGER :: i,my_result,sym_type
    REAL :: min_x,max_x,min_y,max_y,len_x,len_y
    CHARACTER(SIZE(x,2)) :: types,lines,colors
    curve_count=SIZE(x,2)
    n_points=SIZE(x,1)
    types="" ; lines="" ; colors=""
    
    IF(PRESENT(type_spec)) types=type_spec
    IF(PRESENT(line_spec)) lines=line_spec   
    IF(PRESENT(color_spec)) colors=color_spec  
      
    ! Find the scaling from the maximum range of values in the data
    ! or simply use axis if present:
    IF(.NOT.PRESENT(axis)) THEN
      min_x=MINVAL(x)
      max_x=MAXVAL(x)
      min_y=MINVAL(y)
      max_y=MAXVAL(y)
    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
    IF(PRESENT(my_func)) my_result=my_func(2)

    ! This big loop sets the appropriate symbol or line type and color
    ! and then plots the corresponding curve:
    DO i=1,curve_count
        SELECT CASE(types(i:i))
          CASE('S')
            CALL INCMRK(-1)
            SELECT CASE(lines(i:i))
              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(lines(i:i))
              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(colors(i:i))
          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
      CALL CURVE(x(:,i),y(:,i),n_points)
    END DO

  END SUBROUTINE Plot2DMatrix
  
  !___________________________________________________________________________________
  ! This subroutine can be used to plot 3D data with given specifications
  ! It is almost equivalent to Plot2D
  !___________________________________________________________________________________  
  SUBROUTINE Plot3D(x1,y1,z1,spec1,x2,y2,z2,spec2,x3,y3,z3,spec3,x4,y4,z4,spec4,x5,y5,z5,spec5,x6,y6,z6,spec6,my_func,axis,view)

    REAL, DIMENSION(:), TARGET, OPTIONAL :: x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6,z1,z2,z3,z4,z5,z6
    CHARACTER(3), OPTIONAL :: spec1,spec2,spec3,spec4,spec5,spec6
    INTEGER, EXTERNAL, OPTIONAL :: my_func
    REAL, DIMENSION(6), OPTIONAL :: axis
    REAL, DIMENSION(3), OPTIONAL :: view 
    
    INTEGER :: my_result
    INTEGER, PARAMETER :: N=6
    TYPE(array_ptr), DIMENSION(N) :: x_arrays,y_arrays,z_arrays
    REAL, DIMENSION(N) :: mins_x,mins_y,maxs_x,maxs_y,mins_z,maxs_z
    REAL :: min_x,max_x,min_y,max_y,min_z,max_z,len_x,len_y,len_z
    INTEGER :: curve_count,i,sym_type
    LOGICAL :: PresentXYZ,PresentSpec
    CHARACTER(3) :: this_spec
    curve_count=0

    DO i=1,N
      CALL TestPresenceXYZ(PresentXYZ,i)
      IF (PresentXYZ) THEN
        CALL AddCurve(i)
      END IF
    END DO

    IF(.NOT.PRESENT(axis)) THEN
      DO i=1,curve_count
        mins_x(i)=MINVAL(x_arrays(i)%p)
        maxs_x(i)=MAXVAL(x_arrays(i)%p)
        mins_y(i)=MINVAL(y_arrays(i)%p)
        maxs_y(i)=MAXVAL(y_arrays(i)%p)
        mins_z(i)=MINVAL(z_arrays(i)%p)
        maxs_z(i)=MAXVAL(z_arrays(i)%p)
      END DO
      min_x=MINVAL(mins_x(1:curve_count))
      max_x=MAXVAL(maxs_x(1:curve_count))
      min_y=MINVAL(mins_y(1:curve_count))
      max_y=MAXVAL(maxs_y(1:curve_count))
      min_z=MINVAL(mins_z(1:curve_count))
      max_z=MAXVAL(maxs_z(1:curve_count))
    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
    IF(PRESENT(my_func)) my_result=my_func(3)

    DO i=1,curve_count
      CALL TestPresenceSpec(PresentSpec,i,this_spec)
      IF(PresentSpec) THEN
        SELECT CASE(this_spec(1:1))
          CASE('S')
            CALL INCMRK(-1)
            SELECT CASE(this_spec(2:2))
              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(this_spec(2:2))
              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(this_spec(3:3))
          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
      CALL CURV3D(x_arrays(i)%p,y_arrays(i)%p,z_arrays(i)%p,SIZE(x_arrays(i)%p))
      NULLIFY(x_arrays(i)%p,y_arrays(i)%p,z_arrays(i)%p)
    END DO

    CONTAINS

    SUBROUTINE TestPresenceXYZ(TestPresentXYZ,k)
      LOGICAL :: TestPresentXYZ
      INTEGER :: k
      SELECT CASE(k)
        CASE(1)
          TestPresentXYZ=PRESENT(x1).AND.PRESENT(y1).AND.PRESENT(z1)
        CASE(2)
          TestPresentXYZ=PRESENT(x2).AND.PRESENT(y2).AND.PRESENT(z2)
        CASE(3)
          TestPresentXYZ=PRESENT(x3).AND.PRESENT(y3).AND.PRESENT(z3)
        CASE(4)
          TestPresentXYZ=PRESENT(x4).AND.PRESENT(y4).AND.PRESENT(z4)
        CASE(5)
          TestPresentXYZ=PRESENT(x5).AND.PRESENT(y5).AND.PRESENT(z5)
        CASE(6)
          TestPresentXYZ=PRESENT(x6).AND.PRESENT(y6).AND.PRESENT(z6)
        CASE DEFAULT
          TestPresentXYZ=.FALSE.
      END SELECT
    END SUBROUTINE TestPresenceXYZ

    SUBROUTINE TestPresenceSpec(TestPresentSpec,k,Spec)
      LOGICAL :: TestPresentSpec
      INTEGER :: k
      CHARACTER(3) :: Spec
      SELECT CASE(k)
        CASE(1)
          TestPresentSpec=PRESENT(spec1)
          Spec=spec1
        CASE(2)
          TestPresentSpec=PRESENT(spec2)
          Spec=spec2
        CASE(3)
          TestPresentSpec=PRESENT(spec3)
          Spec=spec3
        CASE(4)
          TestPresentSpec=PRESENT(spec4)
          Spec=spec4
        CASE(5)
          TestPresentSpec=PRESENT(spec5)
          Spec=spec5
        CASE(6)
          TestPresentSpec=PRESENT(spec6)
          Spec=spec6
        CASE DEFAULT
          TestPresentSpec=.FALSE.
          Spec='SSR'
      END SELECT
    END SUBROUTINE TestPresenceSpec

    SUBROUTINE AddCurve(k)
      INTEGER :: k
      curve_count=curve_count+1
      SELECT CASE(k)
        CASE(1)
          x_arrays(curve_count)%p=>x1
          y_arrays(curve_count)%p=>y1
          z_arrays(curve_count)%p=>z1
        CASE(2)
          x_arrays(curve_count)%p=>x2
          y_arrays(curve_count)%p=>y2
          z_arrays(curve_count)%p=>z2
        CASE(3)
          x_arrays(curve_count)%p=>x3
          y_arrays(curve_count)%p=>y3
          z_arrays(curve_count)%p=>z3
        CASE(4)
          x_arrays(curve_count)%p=>x4
          y_arrays(curve_count)%p=>y4
          z_arrays(curve_count)%p=>z4
        CASE(5)
          x_arrays(curve_count)%p=>x5
          y_arrays(curve_count)%p=>y5
          z_arrays(curve_count)%p=>z5
        CASE(6)
          x_arrays(curve_count)%p=>x6
          y_arrays(curve_count)%p=>y6
          z_arrays(curve_count)%p=>z6
      END SELECT
    END SUBROUTINE AddCurve
  END SUBROUTINE Plot3D
  
  !___________________________________________________________________________________
  ! This subroutine plots 3D curves on a graph
  ! with specified line-types, symbol-types and colors.
  ! This time the data are all stored in the same matrix as COLUMNS!
  ! Again, my_func can perform any additional level 3 routines (after GRAF) needed.
  ! If specified, axis determines the scaling and view the ANGLE of view.
  !___________________________________________________________________________________   
  SUBROUTINE Plot3DMatrix(x,y,z,type_spec,line_spec,color_spec,my_func,axis,view)

    REAL, DIMENSION(:,:) :: x
    REAL, DIMENSION(SIZE(x,1),SIZE(x,2)) :: y,z
    CHARACTER(SIZE(x,2)), OPTIONAL :: type_spec,line_spec,color_spec
    INTEGER, EXTERNAL, OPTIONAL :: my_func
    REAL, DIMENSION(6), OPTIONAL :: axis
    REAL, DIMENSION(3), OPTIONAL :: view    

    INTEGER :: curve_count, n_points
    INTEGER :: my_result,i,sym_type
    REAL :: min_x,max_x,min_y,max_y,min_z,max_z,len_x,len_y,len_z
    CHARACTER(SIZE(x,2)) :: types,lines,colors
    curve_count=SIZE(x,2)
    n_points=SIZE(x,1)
    types="" ; lines="" ; colors=""
    
    IF(PRESENT(type_spec)) types=type_spec
    IF(PRESENT(line_spec)) lines=line_spec   
    IF(PRESENT(color_spec)) colors=color_spec  
      
    ! Find the scaling from the maximum range of values in the data
    ! or simply use axis if present:
    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_y=axis(5)
      max_y=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
    IF(PRESENT(my_func)) my_result=my_func(2)

    ! This big loop sets the appropriate symbol or line type and color
    ! and then plots the corresponding curve:
    DO i=1,curve_count
        SELECT CASE(types(i:i))
          CASE('S')
            CALL INCMRK(-1)
            SELECT CASE(lines(i:i))
              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(lines(i:i))
              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(colors(i:i))
          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
      CALL CURV3D(x(:,i),y(:,i),z(:,i),n_points)
    END DO

  END SUBROUTINE Plot3DMatrix
  
  !___________________________________________________________________________________  
  ! This subroutine plots a 3D or 2D surface (color) plot from given data
  ! NOTE: Only one data set can be given (since this is 3D, this makes sense)
  !___________________________________________________________________________________    
  SUBROUTINE SurfPlot(x,y,z,spec,my_func,color_range,axis,view)
    REAL, DIMENSION(:), OPTIONAL :: x,y
    REAL, DIMENSION(:,:), OPTIONAL :: z
    CHARACTER(3), OPTIONAL :: spec
    INTEGER, EXTERNAL, OPTIONAL :: my_func
    REAL, DIMENSION(2), OPTIONAL :: color_range
    REAL, DIMENSION(6), OPTIONAL :: axis
    REAL, DIMENSION(3), OPTIONAL :: view

    INTEGER :: my_result
    REAL :: min_x,max_x,min_y,max_y,min_z,max_z,len_x,len_y,len_z
    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)) CALL ZSCALE(color_range(1),color_range(2))
    IF(PRESENT(spec)) this_spec=spec

    ! 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:
      CASE('2')
        ! 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-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)
        CALL CRVMAT(z,size(x),size(y),i,i)
      ! 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
  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

END MODULE Graphics