C 
C Program GRDPS to Display GRiD values in Post Script
C
C Version: 5.10
C Date: 1997, October 25
C
C Coded by Ludek Klimes
C             Department of Geophysics, Charles University Prague
C             Ke Karlovu 3,  121 16  Praha 2,  Czech Republic
C             E-mail: klimes@seis.karlov.mff.cuni.cz
C
C.......................................................................
C
C                                                    
C Main input data read from the * input device:
C     The data are read in by the list directed input (free format), by
C     a single READ statement.  The following items are read:
C     'SEP'...String in apostrophes containing the name of the input
C             SEP parameter file with the data specifying grid
C             dimensions.  The file may also contain additional data
C             specific to this program.
C     'DAT'...String in apostrophes containing the name of the input
C             SEP parameter file with data additional to the data
C             already contained in file FSEP.  These data may be, e.g.,
C             the data specific to this program.
C             'FDAT'=' ' implies no additional data.
C             Description of files SEP and DAT
C     'FILEIN'..String in apostrophes containing the name of the input
C             ASCII file with the grid values.
C     'BBOX'..String in apostrophes containing the name of the output
C             PostScript file or its prolog part, see 'FILEPS'.
C     'FILEPS'..String in apostrophes which may contain the name of the
C             output file containing hexadecimal grid data and the
C             PostScript trailer.  If 'FILEPS' is blank, file 'FBBOX'
C             contains the whole PostScript file.  Otherwise, file
C             'FBBOX' contains just the prolog, setup, and object
C             definition parts, whereas the grid data are written to
C             file 'FILEPS' in order to facilitate batch or manual
C             editing of the PostScript definitions.  Resulting
C             PostScript file is then obtained by merging 'FBBOX',
C             possible new definitions and 'FILEPS'.
C     /...    Input data line must be terminated by a slash.
C     Defaults: 'SEP'='grd.h', 'DAT'=' ', 'FILEIN'='grd.out',
C             'BBOX'='grd.ps', 'FILEPS'=' '.
C
C                                                     
C Both data files SEP and DAT have the same form of the SEP
C (Stanford Exploration Project) parameter files:
C     All the data are specified in the form of PARAMETER=VALUE, e.g.
C     N1=50, with PARAMETER directly preceding = without intervening
C     spaces and with VALUE directly following = without intervening
C     spaces.  The PARAMETER=VALUE couple must be delimited by a space
C     or comma from both sides.
C     The PARAMETER string is not case-sensitive.
C     PARAMETER= followed by a space resets the default parameter value.
C     All other text in the input files is ignored.  The file thus may
C     contain unused data or comments without leading comment character.
C     Everything between comment character # and the end of the
C     respective line is ignored, too.
C     The PARAMETER=VALUE couples may be specified in any order.
C     The last appearance takes precedence.
C Note that the division of the input data into two files FSEP and
C FDAT is, in principle, arbitrary.  The data may be arbitrarily mixed
C or concentrated to a single file, specifying the name of the other
C input file blank.  Data of file FDAT are read after FSEP and take
C thus precedence.
C
C Data specifying grid dimensions (SEP parameter file form):
C These data are common to nearly all programs dealing with SEP-like
C gridded data formats.
C     N1=positive integer... Number of gridpoints along the X1 axis.
C             Default: N1=1
C     N2=positive integer... Number of gridpoints along the X2 axis.
C             Default: N2=1
C     N3=positive integer... Number of gridpoints along the X3 axis.
C             Default: N3=1
C
C Additional data specific to this program (SEP parameter file form):
C     NH=positive integer... Layout of 2-D sections of a 3-D grid.
C             Individual 2-D sections of dimensions N1*N2 are merged
C             into rows NH sections long.  Resulting (N3-1)/NH+1 rows
C             are merged together to form the final image.
C             Default: NH=INT(SQRT(N3))
C     VMIN=real, VMAX=real... Values less than or eqal to VMIN,
C             or greater than or or equal to VMAX will be deemed
C             undefined.  For example, if VMIN=0 is set when plotting
C             velocities, zero velocities are rendered as undefined
C             values, that is usually desirable.
C             Defaults: VMIN=-999999, VMAX=999999.
C     VPLUS=real, VSIGN=real, VCIRC=real... VCIRC is the extent of
C             values corresponding to the whole colour circle RGB.
C             Negative VCIRC corresponds to the opposite direction
C             around the colour circle (BGR).
C             Undefined values are displayed in gray.
C             Defaults: VPLUS=0, VSIGN=1,
C               VCIRC=SIGN(GMAX-GMIN+VPLUS,VSIGN),
C               where GMIN and GMAX are the minimum and maximum defined
C               grid values.
C             Hint: Specify VPLUS=1 when plotting integer values.
C             Hint: To plot velocities with blue minimum velocity,
C               increasing through green, yellow and red again to blue,
C               specify VMIN=0 and VSIGN=-1 with default VREF and CREF.
C     VREF=real, CREF=real... Value VREF corresponds to colour CREF,
C             where Red=0, Yellow=1/6, Green=2/6, Cyan=3/6, Blue=4/6,
C             Magenta=5/6, Red=1, Yellow=7/6, Green=8/6, etc.
C             Default:  VREF=GMIN, CREF=4/6 (blue).
C               where GMIN is the minimum defined grid value.
C               Default VMIN, VMAX, VPLUS, VSIGN, VCIRC, VREF and CREF
C               render the central value (GMIN+GMAX)/2 in yellow.
C             Hint: To plot the values oscillating around 0, you may
C               wish to set VREF=0 and CREF=0.166667 to render 0 in
C               yelow.
C     HSIZE=real... Size (in cm) of the image, corresponding to the X1
C             axis (horizontal before a possible rotation).
C             If negative, the values will be displayed from the right
C             to the left.
C             Default: HSIZE=16.0
C     VSIZE=real... Size (in cm) of the image, corresponding to the X2
C             axis (vertical before a possible rotation).
C             If negative, the values will be displayed from the top to
C             the bottom.
C             Default: VSIZE=ABS(HSIZE)*N2/N1  (proportional display)
C     HOFFSET=real... Distance (in cm) of the image from the leftmost
C             paper edge (before a possible rotation).  Controls the
C             horizontal position of the figure.
C             Default: HOFFSET=2.5
C     VOFFSET=real... Distance (in cm) of the image from the bottom
C             paper edge (before a possible rotation).  Controls the
C             vertical position of the figure.
C             Default:
C               if VSIZE.LE.HEIGHT-2*2.5:     VOFFSET=HEIGHT-2.5-VSIZE
C               otherwise if VSIZE.LE.HEIGHT: VOFFSET=(HEIGHT-VSIZE)/2.
C               otherwise:                    VOFFSET=2.5
C             Paper size of HEIGHT=29.7 (in cm) is given by the
C               PARAMETER Fortran statement and may be changed
C               (e.g., to 29.74 for a 11 inch paper).
C     ROTATE=real... Enables to rotate the image by angle specified in
C             degrees (positive counterclockwise).  The image is rotated
C             around the centre of the square paper of size HEIGHT.
C             If applied, the user will probably wish to specify the
C             value of ROTATE=90.
C             Parameters HSIZE,VSIZE,HOFFSET,VOFFSET apply to the image
C             before rotation.
C             Attention: BoundingBox is incorrect if ROTATE is not
C               multiple of 90 degrees.
C             Default: ROTATE=0.
C
C Length units:
C     All length units controlling the size and position of the plot are
C     assumed to be centimetres in the description above.
C     Changing the plotting units
C
C=======================================================================
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
      INTEGER MGRID
      PARAMETER (MGRID=MRAM)
      REAL GRID(MGRID)
      EQUIVALENCE (GRID,RAM)
      INTEGER IGRID(MGRID)
      EQUIVALENCE (GRID,IGRID)
C
C.......................................................................
C                                                   
C     Important program parameters:
      REAL UINCH,HEIGHT,OFFSET,WIDTH
      PARAMETER (UINCH=2.54,HEIGHT=29.7,OFFSET=2.5,WIDTH=16.0)
C
C     UINCH...Number of length units, in which input data controlling
C             the size and position of the plot are expressed, per inch.
C             UINCH=2.54 corresponds to plotting in cm.
C     HEIGHT..Anticipated height of the paper sheet.
C     OFFSET..Left margin, and top or bottom margin for low or high
C             plots, respectively.
C     WIDTH...Default width of the plot.
C
C.......................................................................
C
      CHARACTER*80 FSEP,FDAT,FILEIN,FBBOX,FILEPS
      CHARACTER*1 HEX1(0:15)
      CHARACTER*2 HEX2(0:255)
      INTEGER LU
      REAL UNDEF
      PARAMETER (LU=1)
      PARAMETER (UNDEF=-999999.)
      INTEGER N1,N2,N3,N31,N32,J,J0,I,I0,I1,I2,I31,I32
      REAL VMIN,VMAX,VPLUS,VSIGN,VCIRC,VREF,CREF,ROTATE
      REAL GMIN,GMAX,G,G0,GSTEP,BBOX1,BBOX2,BBOX3,BBOX4,BB1,BB2,BB3,BB4
      REAL BB1P,BB2P,BB3P,BB4P,BB2DEF,BB4DEF,AUX,C,S
      DATA HEX1
     * /'0','1','2','3','4','5','6','7','8','9','A','B','C','D','E','F'/
C
C
C.......................................................................
C
C     Reading main input data:
      WRITE(*,'(A)') ' Enter FSEP,FDAT,FILEIN,FBBOX,FILEPS: '
      FSEP='grd.h'
      FDAT=' '
      FILEIN='grd.out'
      FBBOX='grd.ps'
      FILEPS=' '
      READ(*,*) FSEP,FDAT,FILEIN,FBBOX,FILEPS
      WRITE(*,'(A)') '+                                    '
C
C     Reading all the data from files FSEP and FDAT to the memory
C     (SEP parameter file form):
      CALL RSEP1(LU,FSEP)
      CALL RSEP1(LU,FDAT)
C
C     Recalling the data specifying grid dimensions
C     (arguments: Name of value in input data, Variable, Default):
      CALL RSEP3I('N1',N1,1)
      CALL RSEP3I('N2',N2,1)
      CALL RSEP3I('N3',N3,1)
C     Transforming 2-D arrays to N1*N2*1 and 1-D arrays to N1*1*1
      IF(N1.EQ.1) THEN
        N1=N2
        N2=N3
        N3=1
      END IF
      IF(N1.EQ.1) THEN
        N1=N2
        N2=1
      END IF
      IF(N2.EQ.1) THEN
        N2=N3
        N3=1
      END IF
C     Determining the layout of 2-D slices of 3-D arrays
      CALL RSEP3I('NH',N31,INT(SQRT(FLOAT(N3))+.001))
      N32=(N3-1)/N31+1
      IF(N1*N2*N31*N32.GT.MGRID) THEN
C       GRDPS-01
        PAUSE 'Error GRDPS-01: Too small array RAM(MRAM)'
        STOP
C       Array RAM(MRAM) allocated in include file 'ram.inc' is too small
C       to contain the grid values to be plotted.  You may wish to
C       increse the dimension MRAM in file 'ram.inc'.
C       ram.inc
      END IF
C
C.......................................................................
C
C     Reading input grid values:
      CALL RARRAY(LU,FILEIN,'FORMATTED',.TRUE.,UNDEF,N1*N2*N3,GRID)
C
C     Rearranging 3-D input N1*N2*N3 grid into (N1*N31)*(N2*N32) grid:
C     Extending the grid to N1*N2*N31*N32 points
      DO 1 I=N1*N2*N3+1,N1*N2*N31*N32
        GRID(I)=UNDEF
    1 CONTINUE
C     Transposing the N1*N2*N31*N32 grid to N1*N31*N2*N32 grid
      DO 5 I32=0,N2*N31*(N32-1),N2*N31
        DO 4 I0=0,N2*N31-1
          J0=I0
    2     CONTINUE
            I2=J0/N31
            I31=J0-I2*N31
            J0=I31*N2+I2
          IF(J0.LT.I0) GO TO 2
          I=(I32+I0)*N1
          J=(I32+J0)*N1
          DO 3 I1=1,N1
            I=I+1
            J=J+1
            AUX=GRID(I)
            GRID(I)=GRID(J)
            GRID(J)=AUX
    3     CONTINUE
    4   CONTINUE
    5 CONTINUE
C
      N1=N1*N31
      N2=N2*N32
C
C.......................................................................
C
C     Recalling the data for the plotting area:
      CALL RSEP3R('HSIZE'  ,BB3,WIDTH)
      CALL RSEP3R('HOFFSET',BB1,OFFSET)
C     Default height of the figure
      BB4DEF=ABS(BB3)*FLOAT(N2)/FLOAT(N1)
      CALL RSEP3R('VSIZE'  ,BB4,BB4DEF)
C     Default vertical position of the figure
      IF(ABS(BB4).LE.HEIGHT-2.*OFFSET) THEN
        BB2DEF=HEIGHT-OFFSET-ABS(BB4)
      ELSE IF(ABS(BB4).LE.HEIGHT) THEN
        BB2DEF=(HEIGHT-ABS(BB4))/2.
      ELSE
        BB2DEF=OFFSET
      END IF
      CALL RSEP3R('VOFFSET',BB2,BB2DEF)
      IF(BB3.LT.0.) THEN
        BB1=BB1-BB3
      END IF
      IF(BB4.LT.0.) THEN
        BB2=BB2-BB4
      END IF
      CALL RSEP3R('ROTATE',ROTATE,0.)
C
C     Transformation from centimetres to points:
      BB1P=BB1*72./UINCH
      BB2P=BB2*72./UINCH
      BB3P=BB3*72./UINCH
      BB4P=BB4*72./UINCH
C     Bounding box:
      BBOX1=AMIN1(BB1P,BB1P+BB3P)
      BBOX2=AMIN1(BB2P,BB2P+BB4P)
      BBOX3=AMAX1(BB1P,BB1P+BB3P)
      BBOX4=AMAX1(BB2P,BB2P+BB4P)
      IF(ROTATE.NE.0.) THEN
        C=COS(ROTATE*3.14159/180.)
        S=SIN(ROTATE*3.14159/180.)
        BBOX1=BBOX1-HEIGHT/2.
        BBOX2=BBOX2-HEIGHT/2.
        BBOX3=BBOX3-HEIGHT/2.
        BBOX4=BBOX4-HEIGHT/2.
        AUX  =C*BBOX1-S*BBOX2
        BBOX2=S*BBOX1+C*BBOX2
        BBOX1=AUX
        AUX  =C*BBOX3-S*BBOX4
        BBOX4=S*BBOX3+C*BBOX4
        BBOX3=AUX
        BBOX1=BBOX1+HEIGHT/2.
        BBOX2=BBOX2+HEIGHT/2.
        BBOX3=BBOX3+HEIGHT/2.
        BBOX4=BBOX4+HEIGHT/2.
        AUX  =AMIN1(BBOX1,BBOX3)
        BBOX3=AMAX1(BBOX1,BBOX3)
        BBOX1=AUX
        AUX  =AMIN1(BBOX2,BBOX4)
        BBOX4=AMAX1(BBOX2,BBOX4)
        BBOX2=AUX
      END IF
C
C.......................................................................
C
C     Calculating minimum and maximum values GMIN and GMAX:
      CALL RSEP3R('VMIN',VMIN,-999999.)
      CALL RSEP3R('VMAX',VMAX, 999999.)
      GMIN= 999999.
      GMAX=-999999.
      DO 11 I=1,N1*N2
        G=GRID(I)
        IF(G.NE.UNDEF) THEN
          IF(G.LE.VMIN.OR.VMAX.LE.G) THEN
            GRID(I)=UNDEF
          ELSE
            GMIN=AMIN1(GMIN,G)
            GMAX=AMAX1(GMAX,G)
          END IF
        END IF
   11 CONTINUE
C
C     Rescaling range GMIN...GMAX to 0...254 (undefined values are 255):
      IF(GMIN.NE.GMAX) THEN
        GSTEP=(GMAX-GMIN)/254.
      ELSE
        GSTEP=1.
      END IF
      G0=GMIN-GSTEP/2.
      DO 12 I=1,N1*N2
        IF(GRID(I).NE.UNDEF) THEN
          G=(GRID(I)-G0)/GSTEP
          IGRID(I)=INT(G)
        ELSE
          IGRID(I)=255
        END IF
   12 CONTINUE
C
C     Preparing hexadecimal coding table:
      I=0
      DO 22 I2=0,15
        DO 21 I1=0,15
          HEX2(I)(1:1)=HEX1(I2)
          HEX2(I)(2:2)=HEX1(I1)
          I=I+1
   21   CONTINUE
   22 CONTINUE
C
C     Parameters controlling colours:
      CALL RSEP3R('VPLUS',VPLUS,0.)
      CALL RSEP3R('VSIGN',VSIGN,0.)
      CALL RSEP3R('VCIRC',VCIRC,SIGN(GMAX-GMIN+VPLUS,VSIGN))
      CALL RSEP3R('VREF' ,VREF ,GMIN)
      CALL RSEP3R('CREF' ,CREF ,2./3.)
C
C     Writing PostScript prolog:
      OPEN(LU,FILE=FBBOX)
      WRITE(LU,'(A/A,4I6,/(A))')
     *'%!PS-Adobe-3.0',
     *'%%BoundingBox:',INT(BBOX1+.5),INT(BBOX2+.5),
     *                 INT(BBOX3+.5),INT(BBOX4+.5),
     *'%%EndComments',
     *'%%BeginProlog',
     *'%%BeginProcSet: (d2gridps)',
     *'%%Creator: d2gridps',
     *'%-----------------------------------------------------------',
     *'% Default UNDEFINED procedure:',
     *'/UNDEFINED {0 0 0.80 sethsbcolor} bind def',
     *'% Default VALUEtoCOLOR procedure:'
      WRITE(LU,'(3(A,F14.6)/A,F15.6,A/A,F15.6,A/A)')
     *'/VCIRC',VCIRC,' def % Range of grid values:',GMIN,' to',GMAX,
     *'/VREF',VREF,' def',
     *'/CREF',CREF,' def',
     *'/VRED VREF CREF VCIRC mul sub def'
      WRITE(LU,'(A)')
     *'/VALUEtoCOLOR{VRED sub VCIRC div',
     *'              dup truncate sub dup 0 lt {1 add} if',
     *'              1 1 sethsbcolor} bind def',
     *'%-----------------------------------------------------------',
     *'% User-defined VALUEtoCOLOR procedure:'
      IF(FILEPS.NE.' ') THEN
        CLOSE(LU)
        OPEN(LU,FILE=FILEPS)
      END IF
C
C     Writing the plotting procedure:
      WRITE(LU,'(A)')
     *'%-----------------------------------------------------------',
     *'%%EndProcSet',
     *'%%EndProlog',
     *'%-----------------------------------------------------------',
     *'%%BeginSetup',
     *'% Numerical values describing the image:'
      WRITE(LU,'(A,I6,A)')          '/N1',N1,' def'
      WRITE(LU,'(A,I6,A)')          '/N2',N2,' def'
      WRITE(LU,'(A,F8.1,A,F8.3,A)') '/BB1',BB1P,' def %',BB1,'cm'
      WRITE(LU,'(A,F8.1,A,F8.3,A)') '/BB2',BB2P,' def %',BB2,'cm'
      WRITE(LU,'(A,F8.1,A,F8.3,A)') '/BB3',BB3P,' def %',BB3,'cm'
      WRITE(LU,'(A,F8.1,A,F8.3,A)') '/BB4',BB4P,' def %',BB4,'cm'
      WRITE(LU,'(A,F8.1,A)')        '/PAPERSIZE',HEIGHT*72./UINCH,' def'
      WRITE(LU,'(A,F8.1,A)')        '/ROTATE',ROTATE,' def'
      WRITE(LU,'(A,F14.6,A)')       '/GMIN',GMIN,' def'
      WRITE(LU,'(A,F14.6,A)')       '/GMAX',GMAX,' def'
      WRITE(LU,'(A)')
     *'%%EndSetup',
     *'%-----------------------------------------------------------',
     *'%%BeginObject: (d2gridps)',
     *'PAPERSIZE  2 div dup translate ROTATE rotate',
     *'PAPERSIZE -2 div dup translate',
     *'/SCALE1 N1 BB3 div def',
     *'/SCALE2 N2 BB4 div def',
     *'/VSTEP GMAX GMIN sub 254 div def',
     *'/RGB 3 string def',
     *'N1 N2 8 [ SCALE1 0 0 SCALE2 BB1 SCALE1 mul neg',
     *'                            BB2 SCALE2 mul neg ]',
     *'{currentfile 1 string readhexstring pop 0 get',
     *' dup 255 eq {pop UNDEFINED}',
     *' {VSTEP mul GMIN add VALUEtoCOLOR} ifelse currentrgbcolor',
     *' 255 mul .5 add cvi RGB exch 2 exch put',
     *' 255 mul .5 add cvi RGB exch 1 exch put',
     *' 255 mul .5 add cvi RGB exch 0 exch put RGB}',
     *'                            bind false 3 colorimage'
C
C     Writing output hexadecimal values:
      WRITE(LU,'(40A2)') (HEX2(IGRID(I)),I=1,N1*N2)
C
C     Writing PostScript trailer:
      WRITE(LU,'(A)')
     *'PAPERSIZE  2 div dup translate ROTATE neg rotate',
     *'PAPERSIZE -2 div dup translate',
     *'%%EndObject',
     *'showpage',
     *'%%EOF'
      CLOSE(LU)
C
      STOP
      END
C
C=======================================================================
C
      INCLUDE 'sep.for'
C     sep.for
      INCLUDE 'forms.for'
C     forms.for
      INCLUDE 'length.for'
C     length.for
C
C=======================================================================
C