C
C Program PTSPS to display point values in Post Script
C
C Version: 6.20
C Date: 2008, June 5
C
C
C.......................................................................
C                                                    
C Description of data files:
C
C Input data read from the standard input device (*):
C     The data are read by the list directed input (free format) and
C     consist of a single string 'SEP':
C     'SEP'...String in apostrophes containing the name of the input
C             SEP parameter or history file with the input data.
C     No default, 'SEP' must be specified and cannot be blank.
C
C                                                     
C Input data file 'SEP':
C     File 'SEP' has the form of a SEP
C     parameter file.  The parameters, which do not differ from their
C     defaults, need not be specified in file 'SEP'.
C Names of the input and output files:
C     GRD='string' ... Name of the input ASCII file with the grid values
C             to be plotted.
C             For general description of files with gridded data refer
C             to file forms.htm.
C             Default: GRD='grd.out'
C     PS='string' ... Name of the output PostScript file.
C             If non-blank PSHEX is given, output file PS contains just
C             the prolog part of the output PostScript file, see PSHEX.
C             If the value of parameter N4 is greater than 1, PS is
C             the output filename for the first spatial grid.
C             The filenames for the subsequent spatial grids are
C             generated in such a way that all digits contained within
C             the filename are considered a number which is increased
C             by 1 for each new time slice.
C             Default: PS='grd.ps'
C     PSHEX='string' ... Optional name of the output file containing
C             hexadecimal grid data and the PostScript trailer.
C             If PSHEX is blank (default), file PS contains the whole
C             PostScript file.
C             Otherwise, file PS contains just the prolog, setup, and
C             object definition parts, whereas the grid data are written
C             to file PSHEX.
C             The separation into two files is designed to facilitate
C             batch or manual editing of the PostScript definitions.
C             Resulting PostScript file is then obtained by merging
C             PS, possible new definitions and PSHEX.
C             If the value of parameter N4 is greater than 1, PSHEX is
C             the output filename for the first spatial grid.
C             The filenames for the subsequent spatial grids are
C             generated in such a way that all digits contained within
C             the filename are considered a number which is increased
C             by 1 for each new time slice.
C             Default: PSHEX=' '
C     GRDPS2='string', GRDPS3='string', ... , GRDPS9='string' ... Names
C             of the optional input ASCII files with the grid values to
C             be plotted together with data of GRD.
C             For example, data GRD may control hue (colour) and data
C             GRDPS2 may control brightness or saturation.
C             For general description of files with gridded data refer
C             to file forms.htm.
C             Defaults: GRDSP2=' ', ... , GRDPS9=' '
C Data specifying grid dimensions:
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     N4=positive integer... Number of time slices (snapshots).
C             Individual time slices from the input file 'FILEIN' are
C             separated into different files.  The name(s) of the first
C             file(s) is(are) given by input parameters 'PS' and
C             'PSHEX' described above.  Next filenames are generated in
C             such a way that all digits contained within the filenames
C             are considered as numbers which are increased by 1 for
C             each new time slice.
C             Default: N4=1
C     D4=real... Time interval between two consecutive time slices.
C             Useful only if explicitly specifying time-dependent colour
C             scaling of snapshots (see input parameters VCIRC4, VSAT4
C             and VBRI4 below).
C             Default: D4=1.
C Additional data specific to this program:
C (a) Dimensions and layout:
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     UNIT='string'... All lengths controlling the size and position of
C             the plot are assumed to be expressed in the units given
C             by the string.  The units also influence the default
C             paper size, plot size and margins.  Allowed values:
C             UNIT='cm': centimetres (default),
C             UNIT='in': inches (1in=2.54cm),
C             UNIT='pt': big points (1in=72pt).
C             Points (pt) are useful to generate plots for conversion to
C             bitmap forms, e.g., using GhostScript.
C             Default: UNIT='cm'
C     XSIGN=real... Determines the sign of the default value of HSIZE.
C             Default: XSIGN=1.
C     HSIZE=real... Size (in UNITs) of the image, corresponding to the
C             X1 axis (horizontal before a possible rotation).
C             If negative, the values will be displayed from the right
C             to the left.
C             Default: HSIZE=SIGN( 16.0,XSIGN) for UNIT='cm',
C                      HSIZE=SIGN(  6.5,XSIGN) for UNIT='in',
C                      HSIZE=SIGN(NH*N1,XSIGN) for UNIT='pt'.
C     YSIGN=real... Determines the sign of the default value of VSIZE.
C             Default: YSIGN=1.
C     VSIZE=real... Size (in UNITs) of the image, corresponding to the
C             X2 axis (vertical before a possible rotation).
C             If negative, the values will be displayed from the top to
C             the bottom.
C             Default (proportional display):
C               VSIZE=SIGN(HSIZE*((N3-1)/NH+1)*N2/(NH*N1),YSIGN), i.e.,
C               VSIZE=SIGN(HSIZE*N2/N1,YSIGN) for N3=1.
C     HOFFSET=real... Distance (in UNITs) 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 for UNIT='cm',
C                      HOFFSET=1.0 for UNIT='in',
C                      HOFFSET=0.0 for UNIT='pt'.
C     VOFFSET=real... Distance (in UNITs) 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     HEIGHT=real... Height of the paper in a portrait position.
C             Default: HEIGHT=29.7 for UNIT='cm',
C                      HEIGHT=11.0 for UNIT='in',
C                      HEIGHT=((N3-1)/NH+1)*N2 for UNIT='pt'.
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 (b) Range of displayed values:
C     VMIN=real, VMAX=real... Values less than or equal to VMIN,
C             or greater than 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=999999000000.
C     R=real, G=real, B=real... Colour of the undefined
C             values.
C             Defaults: R=0.80, G=0.80, B=0.80 (light grey)
C (c) Colour interpretation of values contained within file GRD:
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=(GMAX1-GMIN1+VPLUS)*VSIGN,
C               where GMINi and GMAXi are the minimum and maximum grid
C               values defined in file 'GRDi'.
C               If the above value of VCIRC is zero (e.g., if
C               GMAX1=GMIN1), the default value of VCIRC is changed to
C               VCIRC=1.
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=GMIN1, CREF=4/6 (blue).
C               where GMIN1 is the minimum defined grid value.
C               Default VMIN, VMAX, VPLUS, VSIGN, VCIRC, VREF and CREF
C               render the central value (GMIN1+GMAX1)/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               yellow.
C     VCIRC4=real... If VCIRC is specified, VCIRC4 enables to modify it
C             in dependence on progressing time levels I4=1,2,...,N4:
C               VCIRC(I4)=VCIRC*EXP(-VCIRC4*(I4-1)*D4)
C             Default:  VCIRC4=0.
C (d) Saturation interpretation of values contained within file GRDPS2:
C     VSAT=real... The extent of values corresponding to the whole
C             saturation range from 0 (white) to 1 (fully saturated
C             colours).
C             Default:  VSAT=GMAX2-GMIN2 (or VSAT=1 for GMIN2=GMAX2).
C     VWHITE=real... Value corresponding to white (saturation=0).
C             Default:  VWHITE=GMIN2 (strictly speaking, =GMAX2-VSAT).
C     VSAT4=real... If VSAT is specified, VSAT4 enables to modify it in
C             dependence on progressing time levels I4=1,2,...,N4:
C               VSAT(I4)=VSAT*EXP(-VSAT4*(I4-1)*D4)
C             Default:  VSAT4=0.
C (e) Brightness interpretation of values contained within file GRDPS3:
C     VBRI=real... The extent of values corresponding to the whole
C             brightness range from 0 (black) to 1 (bright colours).
C             Default:  VBRI=GMAX3-GMIN3 (or VBRI=1 for GMIN3=GMAX3).
C     VBLACK=real... Value corresponding to black (brightness=0).
C             Default:  VBLACK=GMIN3 (strictly speaking, =GMAX3-VBRI).
C     VBRI4=real... If VBRI is specified, VBRI4 enables to modify it in
C             dependence on progressing time levels I4=1,2,...,N4:
C               VBRI(I4)=VBRI*EXP(-VBRI4*(I4-1)*D4)
C             Default:  VBRI4=0.
C (f) Form of the output PostScript file:
C     SHOWPAGE=integer... PostScript command 'showpage' at the end of
C             file directs printer to print and delete the picture.
C             This is usually what we want.  However, sometimes we may
C             wish to overlay the picture with another one.  In this
C             case, we wish to remove the 'showpage'.
C             SHOWPAGE=0: Two last lines of the file containing strings
C               'showpage' and '%%EOF' are not written.
C             SHOWPAGE=1: The lines are written.
C             Default:  SHOWPAGE=1
C
C=======================================================================
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
      INTEGER IRAM(MRAM)
      EQUIVALENCE (IRAM,RAM)
C
C
C.......................................................................
C
      CHARACTER*2 UNIT
      REAL UNITPT,HEIGHT,OFFSET,WIDTH
C
C     UNIT... One of: 'cm', 'in', 'pt'.
C     UNITPT...Size of the length unit, in which input data controlling
C             the size and position of the plot are expressed, in big
C             points (pt).  E.g., UNITPT=72./2.54 corresponds to
C             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
      INTEGER MFILE,LU
      PARAMETER (MFILE=9,LU=10)
      CHARACTER*80 FSEP,FPS,FPTS,FPSHEX,FILE(MFILE)
      CHARACTER*255 TEXT
      CHARACTER*1 HEX1(0:15)
      CHARACTER*2 HEX2(0:255)
      CHARACTER*6 FGRDPS
      INTEGER LUIN(MFILE)
      REAL UNDEF
      INTEGER NFILE,IFILE,N1,N2,N3,N31,N32,N4,J,J0,I,I0,I1,I2,I31,I32,I4
      INTEGER ISHOW
      REAL XSIGN,YSIGN
      REAL VMIN,VMAX,VPLUS,VSIGN,VCIRC,VREF,CREF,ROTATE,RED,GREEN,BLUE
      REAL GMIN(MFILE),GMAX(MFILE),G,G0,GSTEP
      REAL BBOX1,BBOX2,BBOX3,BBOX4,BB1,BB2,BB3,BB4
      REAL BB1P,BB2P,BB3P,BB4P,BB2DEF,BB4DEF,AUX,VAUX,C,S
      DATA LUIN /1,2,3,4,5,6,7,8,9/
      DATA HEX1
     * /'0','1','2','3','4','5','6','7','8','9','A','B','C','D','E','F'/
C
      UNDEF=UARRAY()
C
C-----------------------------------------------------------------------
C
C     Reading name of SEP file with input data:
      WRITE(*,'(A)') '+PTSPS: Enter input filename: '
      FSEP=' '
      READ(*,*) FSEP
      WRITE(*,'(A)') '+PTSPS: Working ...           '
C
C     Reading all data from the SEP file into the memory:
      IF (FSEP.NE.' ') THEN
        CALL RSEP1(LU,FSEP)
      ELSE
C       PTSPS-01
        CALL ERROR('PTSPS-01: SEP file not given')
C       Input file in the form of the SEP (Stanford Exploration Project)
C       parameter or history file must be specified.
C       There is no default filename.
      ENDIF
C
C     Reading input parameters from the SEP file:
      CALL RSEP3T('PTS',FPTS ,'pts.out')
      CALL RSEP3T('PS',FPS,'grd.ps')
      CALL RSEP3I('SHOWPAGE',ISHOW,1)
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)
      CALL RSEP3I('N4',N4,1)
C
C     Reading points:
      OPEN(LU,FILE=FPTS,STATUS='OLD')
      READ(LU,*) (TEXT,I=1,20)
      NRAMP=0
      KQ=5
      NQ=KQ
   20 CONTINUE
        IF(NRAMP+KQ.GT.MRAM) THEN
C         PTSPS-03
          CALL ERROR('PTSPS-03: Too small array RAM')
        END IF
        TEXT='$'
        DO 21 I=NRAMP+1,NRAMP+KQ
          RAM(I)=0.
   21   CONTINUE
        READ(LU,*,END=29) TEXT,(RAM(I),I=NRAMP+1,NRAMP+KQ)
C       WRITE(*,*) TEXT,(RAM(I),I=NRAMP+1,NRAMP+KQ)
        IF(TEXT.EQ.'$') THEN
          GO TO 29
        END IF
C       Number of storage locations in RAM used for the points
        NRAMP=NRAMP+NQ
      GO TO 20
   29 CONTINUE
      CLOSE(LU)
      MRAMP=NRAMP
      WRITE(*,*) NRAMP
C
C     Loop over N3
      NN3=1
      NSTART=1
      NEND=0
      I=NSTART
  150 CONTINUE
        IF(NN3.GT.N3) GO TO 200
C       write(*,*)RAM(I+2)
        IF(RAM(I+2).EQ.NN3) THEN
          NEND=NEND+KQ
          I=I+KQ
          GO TO 150
        END IF
C
C       Calculating minimum and maximum values GMIN(*) and GMAX(*):
        CALL RSEP3R('VMIN',VMIN,-999999.)
        CALL RSEP3R('VMAX',VMAX, 999999000000.)
        NFILE=1
        KOLUMN=4
        DO 33 IFILE=1,NFILE
          GMINI=VMAX
          GMAXI=VMIN
          DO 31 I=NSTART-1+KOLUMN,NEND,NQ
            G=RAM(I)
C           WRITE(*,*) G
            IF(G.NE.UNDEF) THEN
              IF(G.LE.VMIN.OR.VMAX.LE.G) THEN
                RAM(I)=UNDEF
              ELSE
                GMINI=AMIN1(GMINI,G)
                GMAXI=AMAX1(GMAXI,G)
              END IF
            END IF
   31     CONTINUE
C
          GMIN(IFILE)=GMINI
          GMAX(IFILE)=GMAXI
   33   CONTINUE
C
C       Recalling the plotting unit and setting default dimensions:
        CALL RSEP3T('UNIT',UNIT,'cm')
        CALL LOWER(UNIT)
        IF(UNIT.EQ.'cm') THEN
          UNITPT=72./2.54
          HEIGHT=29.7
          OFFSET=2.5
          WIDTH=16.0
        ELSE IF(UNIT.EQ.'in') THEN
          UNITPT=72.
          HEIGHT=11.0
          OFFSET=1.0
          WIDTH=6.5
        ELSE IF(UNIT.EQ.'pt') THEN
          UNITPT=1.
          HEIGHT=FLOAT(N32*N2)
          OFFSET=0.0
          WIDTH=FLOAT(N31*N1)
        ELSE
C         PTSPS-02
          CALL ERROR('PTSPS-02: Unrecognized plotting units')
C         Allocated plotting units are UNIT='cm', UNIT='in', UNIT='pt'.
        END IF
C
C.......................................................................
C
C       Recalling the data for the plotting area:
        CALL RSEP3R('XSIGN'  ,XSIGN,1.)
        CALL RSEP3R('YSIGN'  ,YSIGN,1.)
        AUX=HEIGHT
        CALL RSEP3R('HEIGHT' ,HEIGHT,AUX)
        CALL RSEP3R('HSIZE'  ,BB3,SIGN(WIDTH,XSIGN))
        CALL RSEP3R('HOFFSET',BB1,OFFSET)
C       Default height of the figure
        BB4DEF=ABS(BB3)*FLOAT(N2)/FLOAT(N1)
        CALL RSEP3R('VSIZE'  ,BB4,SIGN(BB4DEF,YSIGN))
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 plotting units (e.g. centimetres) to points:
        BB1P=BB1*UNITPT
        BB2P=BB2*UNITPT
        BB3P=BB3*UNITPT
        BB4P=BB4*UNITPT
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*UNITPT/2.
          BBOX2=BBOX2-HEIGHT*UNITPT/2.
          BBOX3=BBOX3-HEIGHT*UNITPT/2.
          BBOX4=BBOX4-HEIGHT*UNITPT/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*UNITPT/2.
          BBOX2=BBOX2+HEIGHT*UNITPT/2.
          BBOX3=BBOX3+HEIGHT*UNITPT/2.
          BBOX4=BBOX4+HEIGHT*UNITPT/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(*):
C
C       Parameters controlling colours:
        CALL RSEP3R('D4'    ,D4    ,1.)
        TIME=FLOAT(I4-1)*D4
        CALL RSEP3R('VPLUS' ,VPLUS ,0.)
        CALL RSEP3R('VSIGN' ,VSIGN ,1.)
        CALL RSEP3R('VCIRC4',VAUX  ,0.)
        VAUX=EXP(VAUX*TIME)
        AUX=(GMAX(1)-GMIN(1)+VPLUS)*VSIGN*VAUX
        IF(AUX.EQ.0.) AUX=1.
        CALL RSEP3R('VCIRC' ,VCIRC ,AUX)
        VCIRC=VCIRC/VAUX
        CALL RSEP3R('VREF'  ,VREF  ,GMIN(1))
        CALL RSEP3R('CREF'  ,CREF  ,2./3.)
        IF(NFILE.GE.2) THEN
          CALL RSEP3R('VSAT4' ,VAUX  ,0.)
          VAUX=EXP(VAUX*TIME)
          AUX=(GMAX(2)-GMIN(2))*VAUX
          IF(AUX.EQ.0.) AUX=1.
          CALL RSEP3R('VSAT'  ,VSAT  ,AUX)
          VSAT=VSAT/VAUX
          CALL RSEP3R('VWHITE',VWHITE,GMAX(2)-AUX)
        END IF
        IF(NFILE.GE.3) THEN
          CALL RSEP3R('VBRI4' ,VAUX  ,0.)
          VAUX=EXP(VAUX*TIME)
          AUX=(GMAX(3)-GMIN(3))*VAUX
          IF(AUX.EQ.0.) AUX=1.
          CALL RSEP3R('VBRI'  ,VBRI  ,AUX)
          VBRI=VBRI/VAUX
          CALL RSEP3R('VBLACK',VBLACK,GMAX(3)-AUX)
        END IF
C
C       Writing PostScript prolog:
C       FPS(1:9)='pts000.ps'
        FPS(6:6)=CHAR(ICHAR('0')+NN3-(NN3/10)*10)
        FPS(5:5)=CHAR(ICHAR('0')+NN3/10)
        OPEN(LU,FILE=FPS)
        WRITE(*,'(''+'',79('' ''))')
        WRITE(*,'(2A)') '+Writing: ',FPS(1:MIN0(LEN(FPS),70))
        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: (grdps)',
     *  '%%Creator: grdps',
     *  '%-----------------------------------------------------------',
     *  '% Default UNDEFINED procedure:'
C520    WRITE(LU,'(A)')
C520 *  '/UNDEFINED {0 0 0.80 sethsbcolor} bind def'
        CALL RSEP3R('R',RED  ,0.80)
        CALL RSEP3R('G',GREEN,0.80)
        CALL RSEP3R('B',BLUE ,0.80)
        WRITE(LU,'(A,3(F4.2,1X),A)')
     *  '/UNDEFINED {',RED,GREEN,BLUE,'setrgbcolor} bind def'
        WRITE(LU,'(A)')
     *  '% Default VALUEtoCOLOR procedure:'
        WRITE(LU,'(3(A,G13.6)/A,G13.6,A/A,G13.6,A/A)')
     *  '/VCIRC  ',VCIRC ,' def % Range of grid values:',
     *                                            GMIN(1),' to',GMAX(1),
     *  '/VREF   ',VREF  ,' def',
     *  '/CREF   ',CREF  ,' def',
     *  '/VRED VREF CREF VCIRC mul sub def'
        IF(NFILE.GE.2) THEN
          WRITE(LU,'(3(A,G13.6)/A,G13.6,A)')
     *    '/VSAT   ',VSAT  ,' def % Range of grid values:',
     *                                            GMIN(2),' to',GMAX(2),
     *    '/VWHITE ',VWHITE,' def'
        END IF
        IF(NFILE.GE.3) THEN
          WRITE(LU,'(3(A,G13.6)/A,G13.6,A)')
     *    '/VBRI   ',VBRI  ,' def % Range of grid values:',
     *                                            GMIN(3),' to',GMAX(3),
     *    '/VBLACK ',VBLACK,' def'
        END IF
        WRITE(LU,'(A)')
     *    '/VALUEtoCOLOR{'
        IF(NFILE.LE.1) THEN
          WRITE(LU,'(A)')
     *      ' VRED sub VCIRC div dup truncate sub dup 0 lt {1 add} if'
        ELSE
          WRITE(LU,'(A,I1,9A))')
     *      ' ',NFILE,' -1 roll',
     *      ' VRED sub VCIRC div dup truncate sub dup 0 lt {1 add} if'
        END IF
        IF(NFILE.LE.1) THEN
          WRITE(LU,'(A)')
     *    ' 1'
        ELSE
          WRITE(LU,'(A,I1,9A))')
     *      ' ',NFILE,' -1 roll',
     *      ' VWHITE sub VSAT div'
        END IF
        IF(NFILE.LE.2) THEN
          WRITE(LU,'(A)')
     *      ' 1'
        ELSE
          WRITE(LU,'(A,I1,9A))')
     *      ' ',NFILE,' -1 roll',
     *      ' VBLACK sub VBRI div'
        END IF
        DO 51 IFILE=4,NFILE
          WRITE(LU,'(A)')
     *    ' pop'
   51   CONTINUE
        WRITE(LU,'(A)')
     *  ' sethsbcolor} bind def',
     *  '%-----------------------------------------------------------',
     *  '% User-defined VALUEtoCOLOR procedure may be inserted here:'
C       IF(FPSHEX.NE.' ') THEN
C         CLOSE(LU)
C         WRITE(*,'(''+'',79('' ''))')
C         WRITE(*,'(2A)') '+Writing: ',FPSHEX(1:MIN0(LEN(FPSHEX),70))
C         OPEN(LU,FILE=FPSHEX)
C       END IF
C
C       Writing the plotting procedure:
        WRITE(LU,'(A)')
     *  '%-----------------------------------------------------------',
     *  '%%EndProcSet',
     *  '%%EndProlog',
     *  '%-----------------------------------------------------------',
     *  '%%BeginSetup',
     *  '% Numerical values describing the image size and position:'
        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*UNITPT,' def'
        WRITE(LU,'(A,F8.1,A)')        '/ROTATE',ROTATE,' def'
        WRITE(LU,'(A)')
     *  '% Scaling of grid values:'
        WRITE(LU,'(A)') '/DX BB3 N1 div 1.5 div def'
        WRITE(LU,'(A)') '/DY BB4 N2 div 1.5 div def'
        WRITE(LU,'(A)') '/HX DX -1.0 div def'
        WRITE(LU,'(A)') '/HY DY -1.0 div def'
        WRITE(LU,'(A)') '/P {'
        WRITE(LU,'(A)') '  VALUEtoCOLOR'
        WRITE(LU,'(A)') '  exch SCALE1 div BB1 add'
        WRITE(LU,'(A)') '  exch SCALE2 div BB2 add moveto'
        WRITE(LU,'(A)') '  HX HY rmoveto DX 0.0 rlineto 0.0 DY rlineto
     *DX neg 0.0 rlineto'
        WRITE(LU,'(A)') 'closepath fill'
        WRITE(LU,'(A)') '} bind def'
        WRITE(LU,'(A)') '%%EndSetup',
     *  '%-----------------------------------------------------------',
     *  '%%BeginObject: (grdps)',
     *  'PAPERSIZE  2 div dup translate ROTATE rotate',
     *  'PAPERSIZE -2 div dup translate',
     *  '/SCALE1 N1 BB3 div def',
     *  '/SCALE2 N2 BB4 div def'
C       Plot frame
        WRITE(LU,'(A)') 'BB1 BB2 moveto BB3 0.0 rlineto 0.0 BB4 rlineto
     *BB3 neg 0.0 rlineto'
        WRITE(LU,'(A)') 'closepath stroke'
        DO 30 I=NSTART,NEND,KQ
          WRITE(LU,'(3(G13.6),A)') RAM(I),RAM(I+1),RAM(I+3),' P'
   30   CONTINUE
C
C       Writing PostScript trailer:
        WRITE(LU,'(A)')
     *  'PAPERSIZE  2 div dup translate ROTATE neg rotate',
     *  'PAPERSIZE -2 div dup translate',
     *  '%%EndObject'
        IF(ISHOW.NE.0) THEN
          WRITE(LU,'(A)')
     *    'showpage',
     *    '%%EOF'
        END IF
        CLOSE(LU)
C       write(*,*)NEND,NN3
        I=NEND+1
        NSTART=NEND+1
        NN3=NN3+1
      GO TO 150
  200 CONTINUE
C
      WRITE(*,'(A)') '+PTSPS: Done.                 '
      STOP
      END
C
C=======================================================================
C
      INCLUDE 'error.for'
C     error.for
      INCLUDE 'sep.for'
C     sep.for
      INCLUDE 'forms.for'
C     forms.for
      INCLUDE 'length.for'
C     length.for
C
C=======================================================================
C                                                                 
C C To be done: C Loop over NN3=1,N3 C I=1,N3 C NSTART, NEND C NSTART=1 C NN3=1 C YYY Loop over I=NSTART,NRAMP,KQ C READ() RAM(I) C IF(RAM(I+2).NE.NN3) GO TO XXX C NEND=NEND+KQ C XXX Generate PostScript file C GMIN,GMAX for RAM(I) I=NSTART,NEND,KQ C NSTART=NEND C GO TO YYY C Done: C NVRTX -> NRAMP C