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