C
C Program SP (Seismogram Plotting) to plot seismograms previously stored C in the GSE data exchange format. C C Version: 5.20 C Date: 1998, November 4 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 Contents: C Description of data files C List of external procedures C C....................................................................... C C Description of data files: 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 data consist of character strings C which have to be enclosed in apostrophes. Note that the C interactive * external unit may be piped or redirected to a file. C The following items are read: C (1) 'SEP' / C 'SEP'...String in apostrophes containing the name of the main C input data file in the SEP format. C File 'SEP' contains the name of the GSE file with the C seimogram to be plotted, names of 3 output PostScript C files with the plots of 3 components and the data C describing how to plot the seismograms (e.g., scale, C minimum and maximum values of the coordinates, etc.). C /... It is recommended to terminate these data by a slash. C Default: 'SEP'='ss.h'. C C C Data file 'SEP' has the form of the SEP (Stanford Exploration C Project) parameter file: 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 Names of the input files: C SS='string'... String with the name of the input data file in the C GSE data exchange format, containing the seismograms to be C plotted. C Description of file GSE C Default: SS='ss.gse' C SS1='string', SS2='string', ..., SS9='string'... Strings with the C additional, optional names of the input data files in the C GSE data exchange format, containing the seismograms to be C plotted simultaneously with seismograms given by parameter C SS into the same figures. The order of plottting is SS, C SS1, SS2, ..., SS9, considering just nonblank filenames. C The seismograms are plotted in colours given by parameters C KOLOR, KOLOR1, KOLOR2, ..., KOLOR9, respectively. Refer C to file calcops.rgb for the C definitions of the colours. The frame and description C remain in colour number 1, which is probably black. C Defaults: SS1=' ', SS2=' ', SS3=' ', ..., SS9=' ' C SRC='string'... String with the name of the input data file C with the coordinates of the hypocentre. If left blank C (default), the hypocentral coordinates are deemed zero. C The hypocentral coordinates are used only for eventual C travel-time reduction on a given velocity, or for C amplitude power scaling. Otherwise, this filename need C not be specified. This file also need not be specified if C the seismograms are generated by programs 'greenss.for' or C 'ss.for', because those programs enter the hypocentral C coordinates directly to the comments of the data section C in GSE data exchange file. The hypocentral coordinates C from the GSE file take precedence of the ones from the C 'SRC' file. C Description of file SRC C Default: 'SRC'=' ' C REC='string'... String with the name of the input data file C with the list of the receivers. C If REC=' ' seismograms at all receivers will be plotted, C otherwise only seismograms at the receivers listed in the C file will be plotted. C This file also enables to determine the number NREC of all C receivers, in which the seimograms may be plotted, for C scaling purposes. If REC=' ', NREC=999 for scaling. C In most cases, file REC will be the same as for the C calculation of seismograms. C Description of file REC C Default: 'REC'=' ' C SPHILI='string'... String with the name of the input data file C with the list of the travel times to be highlighted. C If SPHILI=' ', nothing is highlighted. C Description of file SPHILI C Default: 'SPHILI'=' ' C Names of output files: C SP1='string'... String with the name of the output PostScript file C with the plot of the first component of the synthetic C seismograms. C Default: SP1='ss1.gse' C SP2='string'... String with the name of the output PostScript file C with the plot of the second component of the synthetic C seismograms. C Default: SP2='ss2.gse' C SP3='string'... String with the name of the output PostScript file C with the plot of the third component of the synthetic C seismograms. C Default: SP3='ss3.gse' C Data to control plotting: C Colours: C KOLOR=positive integer, KOLOR1=positive integer, C KOLOR2=positive integer, ..., KOLOR9=positive integer... Colours C to plot seismograms of files SS,SS1, SS2, ..., SS9, C respectively. Colour indices correspond to dummy C parameter INP of CalComp subroutine C NEWPEN. C The colours corresponding to the indices may be defined C or changed in file calcops.rgb. C Default: KOLOR=1, KOLOR1=2, KOLOR2=3, ..., KOLOR9=10 C Data for the time axis (vertical): C SPTMIN=real... Minimum time (or reduced time), i.e. C the time corresponding to the bottom of the C seismogram plot. C Default: SPTMIN=0. C SPTMAX=real... Maximum time (or reduced time), i.e. C the time corresponding to the top of the C seismogram plot. C Default: SPTMAX=1. C SPTLEN=real... Length of the vertical time axis in cm. C Default: SPTLEN=10. C SPTDIV=real... ABS(SPTDIV) is the number of intervals along the C time axis, starting at the bottom. Must be SPTDIV.NE.0. C SPTDIV.GT.0.: the marks at the endpoints of intervals will C be supplemented with corresponding times (or reduced C times). C SPTDIV.LT.0.: the time axis will have no description. C Default: SPTDIV=1. C SPTSUB=real... Number of subintervals to be marked in each C interval. Must be SPTSUB.GT.0. C Default: SPTSUB=1. C SPVRED=real... Reduction velocity. If non-zero, the time at each C receiver is reduced by the value of RR/SPVRED, where RR is C the hypocentral distance. C No time reduction is applied if SPVRED=0. C Default: SPVRED=0. C Data for the distance axis (horizontal): C KODESP=integer... Specifies the distribution and description of C seismograms in the plot. See the description of SPXMIN, C SPXMAX,SPYMIN,SPYMAX below. C Default: KODESP=0 C SPXMIN=real, SPXMAX=real, SPYMIN=real, SPYMAX=real: C For KODESP=0: Horizontal axis is divided into NREC+1 C intervals, where NREC is the number of seismograms to be C plotted, see (1). The I-th seismogram is plotted at the C endpoint of the I-th interval. Horizontal axis has no C description, SPXMIN,SPXMAX,SPYMIN,SPYMAX have no C meaning. C For KODESP=1: Horizontal axis represents the profile with C endpoints (X1,X2)=(SPXMIN,SPYMIN) and C (X1,X2)=(SPXMAX,SPYMAX), C situated in a horizontal plane. The seismograms are C located at the orthogonal projections of the receivers C onto the profile. If SPXDIV.GT.0., the horizontal axis C is supplemented with the values of the X1 coordinates. C For KODESP=2: Horizontal axis represents the profile with C endpoints (X1,X2)=(SPXMIN,SPYMIN) and C (X1,X2)=(SPXMAX,SPYMAX), C situated in a horizontal plane. The seismograms are C located at the orthogonal projections of the receivers C onto the profile. If SPXDIV.GT.0., the horizontal axis C is supplemented with the values of both X1 and X2 C coordinates. C For KODESP=3: Horizontal axis represents the vertical C profile with endpoints X3=SPXMIN and X3=SPXMAX. C The seismograms are located at the horizontal C projections of the receivers onto the profile. C If SPXDIV.GT.0., the horizontal axis is supplemented C with the values of X3 coordinate. C SPYMIN and SPYMAX have no meaning. C For KODESP=4: Horizontal axis represents the hypocentral C distance with endpoints RR=SPXMIN and RR=SPXMAX. C The seismograms are located according to the hypocentral C distances the receivers. C If SPXDIV.GT.0., the horizontal axis is supplemented C with the values of the hypocentral distance. C SPYMIN and SPYMAX have no meaning. C Default: SPXMIN=0., SPXMAX=1., SPYMIN=0., SPYMAX=0. C SPXLEN=real... Length of the horizontal axis in cm. C Default: SPXLEN=FLOAT(NREC+1), where NREC is the number of C seismograms to be plotted. C SPXDIV=real... ABS(SPXDIV) is the number of intervals along the C horizontal axis, starting at the left. C Must be SPXDIV.NE.0. C SPXDIV.GT.0.: the marks at the endpoints of intervals will C be supplemented with the corresponding values. C SPXDIV.LT.0.: horizontal axis will have no description. C Default: SPXDIV=1. C SPXSUB=real... Number of subintervals to be marked in each C interval. Must be SPXSUB.GT.0. C Default: SPXSUB=1. C Amplitude scaling: C NORMSP=integer... Type of amplitude scaling: C NORMSP=0: Maximum amplitude at each trace normalized to C the given value. C Simple and universal option. No other amplitude scaling C parameter usually needs to be specified for this option, C although the input parameters SPAMP and SPAMPi enable C for possible amplitude scaling. C Amplitude scale AA is C AA=SPAMPi*XD/AMAX C where AMAX is the maximum amplitude in each seismogram C and XD is the average distance between seismograms, i.e. C XD=SPXLEN/(NREC+1), where NREC is the number of plotted C seismograms. The receiver file given by input parameter C REC is thus required to determine NREC. C Otherwise: All seismograms have the same scaling or power C scaling. C Amplitude scale AA is C AA=SPAMPi*(RR/SPDIST)**SPOWER C where SPDIST is the hypocentral distance of the receiver C under consideration. C Default: NORMSP=0. C SPAMP=real... Amplitude scale for all 3 components. C Default: SPAMP=1. C SPAMP1=real, SPAMP2=real, ..., SPAMP9=real... Amplitude scales C SPAMP individually set for optional input GSE files C given by parameters SS1, SS2, ..., SS9, respectively. C Defaults: SPAMP1=SPAMP, SPAMP2=SPAMP, ..., SPAMP9=SPAMP C SPAMPX1=real, SPAMPX2=real, SPAMPX3=real... Amplitude scale C multiplicative factors for individual components. C SPAMP will be multiplied by SPAMPX1 for component 1, C by SPAMPX2 for component 2, by SPAMPX3 for component 3. C Default: SPAMPX1=1, SPAMPX2=1, SPAMPX3=1 C SPOWER=real... Exponent of the amplitude power scaling. C Need not be specified if power scaling is not required. C Has no meaning if NORMSP=0. C Default: SPOWER=0. C SPOWER1=real, SPOWER2=real, ..., SPOWER9=real... Exponents of the C amplitude power scaling for optional input GSE files C given by parameters SS1, SS2, ..., SS9, respectively. C Defaults: SPOWER1=SPOWER, SPOWER2=SPOWER, ..., C SPOWER9=SPOWER C SPDIST=real... Reference hypocentral distance for the amplitude C power scaling. C Need not be specified if power scaling is not required. C Has no meaning if NORMSP=0 or SPOWER=0. C Default: SPDIST=1. C SPEXP=real, SPEXPT=real... Exponential scaling of seismograms with C respect to time. The amplitude scale is multiplied by the C factor of EXP(SPEXP*(t-SPEXPT)). C Example 1: Exponential scaling with SPEXP=pi*F/Q, where C F is the dominant frequency and Q is the quality factor, C may roughly compensate for attenuation when plotting C late arrivals. C Example 2: Exponential scaling with SPEXP equal to half C the sum of the positive Lyapunov exponents may C compensate for large geometrical spreading when plotting C late arrivals. C Default: SPEXP=0., SPEXPT=0. C SPEXP1=real, SPEXP2=real, ..., SPEXP9=real... Exponential scaling C set individually for optional input GSE files given by C parameters SS1, SS2, ..., SS9, respectively. C Defaults: SPEXP1=SPEXP, SPEXP2=SPEXP, ..., SPEXP9=SPEXP C Other data to control plotting: C SPCHRH=real... Character height in cm. C Default: SPCHRH=0.4 C SPHIWI=real... Width of the highlighted area when plotting travel C times. C Default: SPHIWI=XD, C where XD is the average distance between seismograms, C i.e. XD=SPXLEN/(NREC+1), where NREC is the number of C plotted seismograms, see the amplitude scaling. C C C Input file GSE with the seismograms to plot: C File in the GSE data exchange format, see the description in file C 'gse.for'. C The 'sp.for' program is looking in the comment lines of the C waveform identification section for the hypocentral coordinates C identified by strings 'XS1 ', 'XS2 ' and 'XS3 '. If they are C present, they overwrite the values from the optional file SRC. C Description of format GSE C C C Input file SRC with the hypocentral coordinates: C (1) / C None to 20 character strings terminated by a slash. C (2) 'SRCNAME',XS1,XS2,XS3,/ C 'SRCNAME'... Name of the source. Not used. C XS1,XS2,XS3... Coordinates of the hypocentre. C Default: XS1=0., XS2=0., XS3=0. C C C Input file REC containing receiver coordinates: C (1) Several strings terminated by / (a slash). C The simplest way is to submit just the /. C (2) Several times (2.1): C (2.1) 'NAMER(IR)',X1REC(IR),X2REC(IR),X3REC(IR),/ C 'NAMER(IR)'... String reserved for the name of the receiver. C No meaning here, anything in apostrophes may be submitted. C X1REC(IR),X2REC(IR),X3REC(IR)... Coordinates of the receiver. C The coordinates need not be present in the file. It may C thus be comfortable to omit them if preparing the list for C the selection of particular receivers for seismogram C plotting. C Default: X1REC(IR)=0, X2REC(IR)=0, X3REC(IR)=0. C (3) / (a slash). C C C Input file SPHILI containing travel times to highlight: C (1) Several strings terminated by / (a slash). C The simplest way is to submit just the /. C (2) Several times (2.1): C (2.1) 'NAMER(IR)',X1REC(IR),X2REC(IR),X3REC(IR),TT,TTERR,K1,K2,/ C 'NAMER(IR)',X1REC(IR),X2REC(IR),X3REC(IR)... Same meaning as in C the receiver file above. C TT... Travel time. If given, it will be plotted as horizontal C line of width SPHIWI and colour K1. It is plotted after C the frame and the corresponding error bar, but before C seismograms. C TTERR...Travel time error. If given, it will be plotted as C a solid rectangle of width SPHIWI and colour K2. It is C plotted after the frame but before the corresponding C travel time and before seismograms. C K1... Colour to plot the travel time. The travel time is not C plotted if K1 is not positive. C K2... Colour to plot the error bar. The error bar is not C plotted if K2 is negative. If K2=0 (white), the contour C of the rectangle if plotted in colour 1 (usually black). C Default: X1REC(IR)=0, X2REC(IR)=0, X3REC(IR)=0.,K1=1,K2=0 C (3) / (a slash). C C....................................................................... C C This Fortran77 file consists of the following external procedures: C SP... Main program to read and plot the seismograms. C SP C FRAME...Subroutine called by the main to plot the rectangular C frame around the seismograms and supplement it with simple C descriptions. C FRAME C C Other external procedures required: C RGSE1,RGSE2... Subroutines of the Fortran 77 file 'gse.for' C (package MODEL), designed to read seismograms in the GSE C data exchange format. C PLOTS,PLOT,SYMBOL,NUMBER... CALCOMP plotting subroutines. For C example, Fortran 77 routines of file 'calcops.for' C (package MODEL) may be used to generate seismogram plots C in the PostScript files. C C======================================================================= C C C C Common block /RAMC/: INCLUDE 'ram.inc' C ram.inc C C Allocation of working array: INTEGER MPTS PARAMETER (MPTS=MRAM) REAL SS(MPTS) EQUIVALENCE (SS,RAM) C C----------------------------------------------------------------------- C C Input and output data filenames and logical unit numbers: INTEGER LU PARAMETER (LU=1) CHARACTER*80 FILSEP,FILESS(0:9),FILSRC,FILREC,FILHIL,FILEPS(3) C REAL UNDEF PARAMETER (UNDEF=-999999.) INTEGER KOLOR(0:9) REAL SPAMP(0:9),SPAMPX(3),SPOWER(0:9),SPEXP(0:9),XPTS(4),YPTS(4) C CHARACTER*1 CHAN,TEXT CHARACTER*80 COMLIN C C Data specifying labels in the plot: CHARACTER*20 KXTEXT(5) INTEGER KX(5) C C Source coordinates transferred through the GSE file: INTEGER NCOM PARAMETER (NCOM=3) REAL VCOM(NCOM),DEF C C List of receivers for plotting: INTEGER MREC,NREC PARAMETER (MREC=1000) CHARACTER*6 STAT,REC(MREC) C DATA KX /0,1,1,1,19/ DATA KXTEXT/' ','X1','X1','X3','HYPOCENTRAL DISTANCE'/ C C....................................................................... C C Main input data: FILSEP='ss.h' WRITE(*,'(A)') '+SP: Enter input filename [''ss.h'']: ' READ(*,*) FILSEP WRITE(*,'(A)') '+SP: Working... ' C C Reading all the data from input SEP parameter file: CALL RSEP1(LU,FILSEP) C C Input and output filenames: CALL RSEP3T('SS' ,FILESS(0),'ss.gse') CALL RSEP3T('SS1',FILESS(1),' ') CALL RSEP3T('SS2',FILESS(2),' ') CALL RSEP3T('SS3',FILESS(3),' ') CALL RSEP3T('SS4',FILESS(4),' ') CALL RSEP3T('SS5',FILESS(5),' ') CALL RSEP3T('SS6',FILESS(6),' ') CALL RSEP3T('SS7',FILESS(7),' ') CALL RSEP3T('SS8',FILESS(8),' ') CALL RSEP3T('SS9',FILESS(9),' ') CALL RSEP3T('SP1',FILEPS(1),'ss1.ps') CALL RSEP3T('SP2',FILEPS(2),'ss2.ps') CALL RSEP3T('SP3',FILEPS(3),'ss3.ps') C C Hypocentral coordinates: VCOM(1)=0. VCOM(2)=0. VCOM(3)=0. CALL RSEP3T('SRC',FILSRC,' ') IF(FILSRC.NE.' ') THEN OPEN(LU,FILE=FILSRC,STATUS='OLD') READ(LU,*) (TEXT,I=1,20) READ(LU,*) STAT,VCOM(1),VCOM(2),VCOM(3) CLOSE(LU) END IF C C Number and list of receivers: NREC=999 CALL RSEP3T('REC',FILREC,' ') IF(FILREC.NE.' ') THEN OPEN(LU,FILE=FILREC,STATUS='OLD') READ(LU,*) (TEXT,I=1,20) DO 1 I=1,MREC REC(I)='$$$$$$' READ(LU,*,END=2) REC(I) IF(REC(I).EQ.'$$$$$$') THEN GO TO 2 END IF 1 CONTINUE C SP-01 CALL ERROR('SP-01: Array dimension MREC small') 2 CONTINUE NREC=I-1 CLOSE(LU) END IF C C Colours: CALL RSEP3I('KOLOR ',KOLOR(0),1) CALL RSEP3I('KOLOR1',KOLOR(1),2) CALL RSEP3I('KOLOR2',KOLOR(2),3) CALL RSEP3I('KOLOR3',KOLOR(3),4) CALL RSEP3I('KOLOR4',KOLOR(4),5) CALL RSEP3I('KOLOR5',KOLOR(5),6) CALL RSEP3I('KOLOR6',KOLOR(6),7) CALL RSEP3I('KOLOR7',KOLOR(7),8) CALL RSEP3I('KOLOR8',KOLOR(8),9) CALL RSEP3I('KOLOR9',KOLOR(9),10) C C Initial values for plotting frame: C Time axis: CALL RSEP3R('SPTMIN',SPTMIN, 0.) CALL RSEP3R('SPTMAX',SPTMAX, 1.) CALL RSEP3R('SPTLEN',SPTLEN,10.) CALL RSEP3R('SPTDIV',SPTDIV, 1.) CALL RSEP3R('SPTSUB',SPTSUB, 1.) CALL RSEP3R('SPVRED',SPVRED, 0.) C Distance axis: CALL RSEP3I('KODESP',KODESP, 0 ) CALL RSEP3R('SPXMIN',SPXMIN, 0.) CALL RSEP3R('SPXMAX',SPXMAX, 1.) CALL RSEP3R('SPYMIN',SPYMIN, 0.) CALL RSEP3R('SPYMAX',SPYMAX, 0.) CALL RSEP3R('SPXLEN',SPXLEN,FLOAT(NREC+1)) CALL RSEP3R('SPXDIV',SPXDIV, 1.) CALL RSEP3R('SPXSUB',SPXSUB, 1.) C Characters: CALL RSEP3R('SPCHRH',SPCHRH, 0.4) C SCY= SPTLEN/(SPTMAX-SPTMIN) XD = SPXLEN/FLOAT(NREC+1) XA = SPXMAX-SPXMIN YA = SPYMAX-SPYMIN C C Amplitude scaling: CALL RSEP3I('NORMSP',NORMSP,0) CALL RSEP3R('SPAMP ',SPAMP(0),1.) CALL RSEP3R('SPAMP1',SPAMP(1),SPAMP(0)) CALL RSEP3R('SPAMP2',SPAMP(2),SPAMP(0)) CALL RSEP3R('SPAMP3',SPAMP(3),SPAMP(0)) CALL RSEP3R('SPAMP4',SPAMP(4),SPAMP(0)) CALL RSEP3R('SPAMP5',SPAMP(5),SPAMP(0)) CALL RSEP3R('SPAMP6',SPAMP(6),SPAMP(0)) CALL RSEP3R('SPAMP7',SPAMP(7),SPAMP(0)) CALL RSEP3R('SPAMP8',SPAMP(8),SPAMP(0)) CALL RSEP3R('SPAMP9',SPAMP(9),SPAMP(0)) CALL RSEP3R('SPAMPX1',SPAMPX(1),1.) CALL RSEP3R('SPAMPX2',SPAMPX(2),1.) CALL RSEP3R('SPAMPX3',SPAMPX(3),1.) CALL RSEP3R('SPDIST' ,SPDIST,1.) CALL RSEP3R('SPOWER' ,SPOWER(0),0.) CALL RSEP3R('SPOWER1',SPOWER(1),SPOWER(0)) CALL RSEP3R('SPOWER2',SPOWER(2),SPOWER(0)) CALL RSEP3R('SPOWER3',SPOWER(3),SPOWER(0)) CALL RSEP3R('SPOWER4',SPOWER(4),SPOWER(0)) CALL RSEP3R('SPOWER5',SPOWER(5),SPOWER(0)) CALL RSEP3R('SPOWER6',SPOWER(6),SPOWER(0)) CALL RSEP3R('SPOWER7',SPOWER(7),SPOWER(0)) CALL RSEP3R('SPOWER8',SPOWER(8),SPOWER(0)) CALL RSEP3R('SPOWER9',SPOWER(9),SPOWER(0)) CALL RSEP3R('SPEXP ',SPEXP(0),0.) CALL RSEP3R('SPEXP1',SPEXP(1),SPEXP(0)) CALL RSEP3R('SPEXP2',SPEXP(2),SPEXP(0)) CALL RSEP3R('SPEXP3',SPEXP(3),SPEXP(0)) CALL RSEP3R('SPEXP4',SPEXP(4),SPEXP(0)) CALL RSEP3R('SPEXP5',SPEXP(5),SPEXP(0)) CALL RSEP3R('SPEXP6',SPEXP(6),SPEXP(0)) CALL RSEP3R('SPEXP7',SPEXP(7),SPEXP(0)) CALL RSEP3R('SPEXP8',SPEXP(8),SPEXP(0)) CALL RSEP3R('SPEXP9',SPEXP(9),SPEXP(0)) CALL RSEP3R('SPEXPT',SPEXPT,0.) C C Higlighting given areas (e.g., travel times with error bars): CALL RSEP3T('SPHILI',FILHIL,' ') CALL RSEP3R('SPHIWI',SPHIWI,XD) C C....................................................................... C C Loop over 3 seismogram components: DO 99 KOMP=1,3 IF(FILEPS(KOMP).NE.' ') THEN C C Initialization of CALCOMP: CALL PLOTN(FILEPS(KOMP),0) CALL PLOTS(0,0,0) C C Plotting frame: CALL NEWPEN(1) IX = KODESP IT = -1 IF(KODESP.GE.3) IX=1 IF(SPXDIV.LT.0.) IX=0 IF(SPTDIV.LT.0.) IT=0 KX1= KX(KODESP+1)+1 CALL FRAME(SPXLEN,SPTLEN,ABS(SPXDIV),SPXSUB, * ABS(SPTDIV),SPTSUB,0,IX,IT,SPCHRH, * SPXMIN,SPXMAX,KX1,KXTEXT(KODESP+1), * SPYMIN,SPYMAX,2,'X2', * SPTMIN,SPTMAX,4,'TIME', * 1,' ',0,0.,0,1,' ',0,0.,0,1,' ',1,' ') C C Higlighting given areas (e.g., travel times with error bars): IF(FILHIL.NE.' ') THEN OPEN(LU,FILE=FILHIL,STATUS='OLD') READ(LU,*) (TEXT,I=1,20) C Loop for areas to highlight: 10 CONTINUE STAT='$$$$$$' T0=UNDEF TD=UNDEF K1=1 K2=0 READ(LU,*,END=19) STAT,XR,YR,ZR,T0,TD,K1,K2 IF(STAT.EQ.'$$$$$$') THEN GO TO 19 END IF IF(T0.EQ.UNDEF) THEN K1=-1 K2=-1 END IF IF(TD.EQ.UNDEF) THEN K2=-1 END IF C Determining the receiver: IF(FILREC.NE.' ') THEN C Loop for receivers DO 11 IREC=1,NREC IF(REC(IREC).EQ.STAT) GO TO 12 11 CONTINUE GO TO 10 END IF 12 CONTINUE IF(SPVRED.NE.0.) THEN RR=SQRT((XR-VCOM(1))**2+(YR-VCOM(2))**2+(ZR-VCOM(3))**2) T0=T0-RR/SPVRED END IF IF(KODESP.EQ.0) THEN IF(FILREC.EQ.' ') THEN C SP-02 CALL ERROR('SP-02: No receiver file specified') C For KODESP=0, filename REC must be specified in the C input data. END IF X=XD*FLOAT(IREC) ELSE IF(KODESP.EQ.1.OR.KODESP.EQ.2) THEN X=SPXLEN*((XR-SPXMIN)*XA+(YR-SPYMIN)*YA)/(XA*XA+YA*YA) ELSE IF(KODESP.EQ.3) THEN X=SPXLEN*(ZR-SPXMIN)/XA ELSE X=SPXLEN*(RR-SPXMIN)/XA END IF C Plotting the highlight: IF(K2.GE.0) THEN XPTS(1)=X-0.5*SPHIWI XPTS(2)=X+0.5*SPHIWI XPTS(3)=X+0.5*SPHIWI XPTS(4)=X-0.5*SPHIWI YPTS(1)=SCY*(T0-TD-SPTMIN) YPTS(2)=SCY*(T0-TD-SPTMIN) YPTS(3)=SCY*(T0+TD-SPTMIN) YPTS(4)=SCY*(T0+TD-SPTMIN) IF(K2.GT.0) THEN CALL NEWPEN(K2) CALL FILL(XPTS,YPTS,4) ELSE CALL NEWPEN(1) CALL PLOT(XPTS(1),YPTS(1),3) CALL PLOT(XPTS(2),YPTS(2),2) CALL PLOT(XPTS(3),YPTS(3),2) CALL PLOT(XPTS(4),YPTS(4),2) END IF END IF IF(K1.GT.0) THEN CALL NEWPEN(K1) CALL PLOT(X-0.5*SPHIWI,SCY*(T0-SPTMIN),3) CALL PLOT(X+0.5*SPHIWI,SCY*(T0-SPTMIN),2) END IF GO TO 10 19 CONTINUE C Closing highlighting file CLOSE(LU) END IF C C Plotting seismograms: C Loop for GSE files DO 49 ISS=0,9 IF(FILESS(ISS).NE.' ') THEN CALL NEWPEN(KOLOR(ISS)) C C Opening input GSE file with the seismograms: OPEN(LU,FILE=FILESS(ISS),STATUS='OLD') CALL RGSE1(LU,TEXT) IREC=0 C C Loop for seismograms 20 CONTINUE 21 CONTINUE C C Selecting the component: 22 CONTINUE CALL RGSE2 * (LU,STAT,CHAN,I,XR,YR,ZR,T0,TD,NPTS,MPTS,SS) IF(NPTS.LE.-1) THEN C End of the GSE file GO TO 40 END IF IF(I.NE.KOMP) GO TO 22 C C Determining the receiver: IF(FILREC.EQ.' ') THEN IREC=IREC+1 ELSE C Loop for receivers DO 31 IREC=1,NREC IF(REC(IREC).EQ.STAT) GO TO 32 31 CONTINUE GO TO 21 END IF 32 CONTINUE C C Plotting the seismogram: 33 CONTINUE CALL RGSE2C(COMLIN,*34) CALL RSEP2(COMLIN) GO TO 33 34 CONTINUE DEF=VCOM(1) CALL RSEP3R('XS1',VCOM(1),DEF) DEF=VCOM(2) CALL RSEP3R('XS2',VCOM(2),DEF) DEF=VCOM(3) CALL RSEP3R('XS3',VCOM(3),DEF) RR=SQRT((XR-VCOM(1))**2+(YR-VCOM(2))**2+(ZR-VCOM(3))**2) C IF(SPEXP(ISS).NE.0.) THEN DO 35 I=1,NPTS SS(I)=SS(I)*EXP(SPEXP(ISS)*(T0+TD*FLOAT(I-1)-SPEXPT)) 35 CONTINUE END IF IF(SPVRED.NE.0.) THEN T0=T0-RR/SPVRED END IF C IF(NORMSP.EQ.0) THEN AA=ABS(SS(1)) DO 36 I=2,NPTS AA=AMAX1(ABS(SS(I)),AA) 36 CONTINUE AA=SPAMP(ISS)*SPAMPX(KOMP)*XD/AA ELSE IF(SPOWER(ISS).EQ.0.) THEN AA=SPAMP(ISS)*SPAMPX(KOMP) ELSE AA=SPAMP(ISS)*SPAMPX(KOMP) * *((RR/SPDIST)**SPOWER(ISS)) END IF END IF IF(KODESP.EQ.0) THEN IF(FILREC.EQ.' ') THEN C SP-03 CALL ERROR('SP-03: No receiver file specified') C For KODESP=0, filename REC must be specified in the C input data. END IF X=XD*FLOAT(IREC) ELSE IF(KODESP.EQ.1.OR.KODESP.EQ.2) THEN X=SPXLEN*((XR-SPXMIN)*XA+(YR-SPYMIN)*YA)/(XA*XA+YA*YA) ELSE IF(KODESP.EQ.3) THEN X=SPXLEN*(ZR-SPXMIN)/XA ELSE X=SPXLEN*(RR-SPXMIN)/XA END IF Y = 0. CALL PLOT(X,Y,3) C DO 38 I=1,NPTS T = T0+TD*FLOAT(I-1) IF (SPTMIN.LE.T.AND.T.LE.SPTMAX) THEN A = X-AA*SS(I) Y = SCY*(T-SPTMIN) IF (T.LT.SPTMIN+TD) THEN C Straight line before the seismogram CALL PLOT(X,Y,2) END IF CALL PLOT(A,Y,2) END IF 38 CONTINUE C C Straight line after the seismogram CALL PLOT(X,Y,2) Y = SPTLEN CALL PLOT(X,Y,2) GO TO 20 40 CONTINUE C End of loop for receivers C C Closing GSE file CLOSE(LU) END IF 49 CONTINUE C End of loop for GSE files C CALL PLOT(0.,0.,999) END IF 99 CONTINUE WRITE(*,'(A)') '+SP: Done. ' STOP END C C======================================================================= C C C SUBROUTINE FRAME(XX,YY,XM,XN,YM,YN,MM,IX,IY,HEIGHT, * XL1,XR1,LX1,KX1,XL2,XR2,LX2,KX2,YL,YR,LY,KY, * L1,K1,M1,A1,N1,L2,K2,M2,A2,N2,L3,K3,L4,K4) C CHARACTER*(*) KX1,KX2,KY,K1,K2,K3,K4 REAL XX,YY,XM,XN,YM,YN,HEIGHT,XL1,XR1,XL2,XR2,YL,YR,A1,A2 INTEGER MM,IX,IY,LX1,LX2,LY,L1,M1,N1,L2,M2,N2,L3,L4 C C Input: C XX... Length of the horizontal axis in cm. C YY... Length of the vertical axis in cm. C XM... The number of intervals along the horizontal axis, C starting at the left. The intervals are marked with C long ticks and are supplemented with the numerical C values if IX is positive. XM must be positive. C XN... Number of subintervals to be marked in each interval C with short ticks. XN must be positive. C YM... The number of intervals along the vertical axis, C starting at the bottom. The intervals are marked with C long ticks and are supplemented with the numerical C values if IY is positive. YM must be positive. C YN... Number of subintervals to be marked in each interval C with short ticks. YN must be positive. C MM... If negative, the top line of the frame is not plotted. C IX... IX=0: No labeling of the horizontal axis. C IX=1: Labeling of the horizontal axis with the first C variable only. C IX=2: Labeling of the horizontal axis with both variables. C IY... IY=0: No labeling of the vertical axis. C IY=1: Labeling of the vertical axis. C HEIGHT... Height of the characters in cm. C XL1,XR1... Values of the first variable along the horizontal axis, C corresponding to left-hand and right-hand corners. C LX1... Length of string KX1. C KX1... First label of the horizontal axis. String to be written C below the horizontal axis. C XL2,XR2... Values of the second variable along the horizontal C axis, corresponding to left-hand and right-hand corners. C LX2... Length of string KX2. C KX2... Second label of the horizontal axis. String to be written C below the horizontal axis. C YL,YR...Values of the variable along the vertical axis, C corresponding to left-hand and right-hand corners. C LY... Length of string KY. C KY... Label of the vertical axis. String to be written to the C left of the vertical axis. C L1... Length of string K1. C K1,A1...String and the number to be written above the frame, C starting 0.5*HEIGHT above the left top corner. C M1... Number A1 is written only if M2 is positive. C N1... Number of decimal places of A1. C L2... Length of string K1. C K2,A2...String and the number to be written above the frame, C ending 0.5*HEIGHT above the right top corner. C M2... Width of A2 in characters. A2 is written only if M2 is C positive. C N2... Number of decimal places of A2. C L2... Length of string K2. C L3... Length of string K3. C K3... String to be written below the frame, starting 5.3*HEIGHT C below the left bottom corner. C L4... Length of string K4. C K4... String to be written below the frame, starting 7.0*HEIGHT C below the left bottom corner. C C----------------------------------------------------------------------- C HEIGH0=.215*HEIGHT HEIGH2=.500*HEIGHT C C Initial values for plotting frame: I0= 0 JX= IABS(IX)-1 MX= INT(XM*XN+0.001) NX= INT(XN+0.5) MY= INT(YM*YN+0.001) NY= INT(YN+0.5) XD= XX/XM/XN YD= YY/YM/YN XH1= (XR1-XL1)/XM XH2= (XR2-XL2)/XM YH = (YR -YL )/YM C C Plotting border 29.7*21.0 cm (landscape): C CALL PLOT( 0. , 0. ,3) C CALL PLOT(29.7, 0. ,2) C CALL PLOT(29.7,21.0,2) C CALL PLOT( 0. ,21.0,2) C CALL PLOT( 0. , 0. ,2) C C Plotting frame XX*YY cm centred on the A4 sheet: C Landscape * X = (29.7-XX)/2. * Y = (21.0-YY)/2. C Portrait X = (21.0-XX)/2. Y = (29.7-YY)/2. C Leaving 2 cm belts for description of axes IF(IX.GE.0) THEN X=X+2.5*HEIGHT END IF IF (IY.NE.0) THEN Y=Y+2.5*HEIGHT END IF C Shifting the origin and plotting the frame CALL PLOT(X,Y,-3) CALL PLOT(XX,0.,2) CALL PLOT(XX,YY,2) IF(MM.GE.0) THEN I = 2 ELSE I = 3 END IF CALL PLOT(0.,YY,I) CALL PLOT(0.,0.,2) C C Description of axes: C C Bottom horizontal axis: X = 0. X1= XL1-XH1 X2= XL2-XH2 DO 16 I=I0,MX IF (MOD(I,NX).NE.0) THEN A = 0.1 ELSE A = 0.2 IF (JX.GE.0) THEN X1= X1+XH1 X2= X2+XH2 IF(IX.GE.0.OR.MOD(I,MX).NE.0) THEN C Determination of the number of decimal places: ccc J = INT(.99-ALOG10(AMAX1(ABS(XH1),ABS(XH2),0.001))) DO 11 J=0,5 IF(AMOD(ABS(XH1)+0.000001,0.1**J).LT.0.000002) THEN GO TO 12 END IF 11 CONTINUE 12 CONTINUE C J is the preferable number of decimal places. JMAX=MAX0(INT(ALOG10(ABS(X1)+0.5*0.1**J))+1,1) IF(X1.LT.0.) JMAX=JMAX+1 C JMAX is the number of digits left to the decimal point. IF (JX.GT.0) THEN DO 13 J=J,5 IF(AMOD(ABS(XH2)+0.000001,0.1**J).LT.0.000002) THEN GO TO 14 END IF 13 CONTINUE 14 CONTINUE C J is the preferable number of decimal places. JMAX=MAX0(INT(ALOG10(ABS(X2)+0.5*0.1**J))+1,JMAX) IF(X1.GE.0..AND.X2.LT.0.) JMAX=JMAX+1 C JMAX is the number of digits left to the decimal point. END IF JMAX=INT(XX/XM/HEIGHT)-JMAX-1 C JMAX is the maximum number of decimal places. J=MIN0(J,JMAX) IF(J.LE.0) THEN J=-1 END IF C J is the number of decimal places. C B = X-( .5*FLOAT(1+JX)*AINT(ALOG10(ABS(X1)+.5))+0.285 * -FLOAT(JX)+.5*FLOAT(J+1) )*HEIGHT IF(X1.LT.0.) THEN B=B-HEIGHT END IF CALL NUMBER(B,-1.8*HEIGHT,HEIGHT,X1,0.,J) IF (JX.GT.0) THEN B = X-HEIGHT*AINT(ALOG10(ABS(X2)+.5))+.715*HEIGHT IF(X2.LT.0.) B=B-HEIGHT CALL NUMBER(B,-3.3*HEIGHT,HEIGHT,X2,0.,J) END IF END IF END IF END IF IF (MOD(I,MX).NE.0) THEN CALL PLOT(X,0.,3) CALL PLOT(X,A ,2) END IF X = X+XD 16 CONTINUE IF (JX.EQ.0) THEN B = XX-XX/XM/2.-HEIGH2*FLOAT(LX1)+HEIGH0 CALL SYMBOL(B,-3.3*HEIGHT,HEIGHT,KX1,0.,LX1) ELSE IF (JX.GT.0) THEN CALL SYMBOL(XX+2.715*HEIGHT,-1.8*HEIGHT,HEIGHT,KX1,0.,LX1) CALL SYMBOL(XX+2.715*HEIGHT,-3.3*HEIGHT,HEIGHT,KX2,0.,LX2) END IF C C Right-hand vertical axis: Y = 0. M = MY DO 23 I=1,M Y = Y+YD A = 0.1 IF (MOD(I,NY).EQ.0) THEN A = 0.2 END IF CALL PLOT(XX,Y,3) CALL PLOT(XX-A,Y,2) 23 CONTINUE C C Top horizontal axis: IF (MM.GE.0) THEN X = XX M = MX-1 DO 33 I=1,M X = X-XD IF (MOD(I,NX).NE.0) THEN A = 0.1 ELSE A = 0.2 END IF CALL PLOT(X,YY,3) CALL PLOT(X,YY-A,2) 33 CONTINUE END IF C C Left-hand vertical axis: Y = 0. Y1= YL-YH DO 45 I=I0,MY IF (MOD(I,NY).NE.0) THEN A = 0.1 ELSE A = 0.2 IF (IY.NE.0) THEN Y1= Y1+YH C C Determination of the number of decimal places: ccc J = INT(.99-ALOG10(AMAX1(ABS(YH),0.001))) DO 41 J=0,5 IF(AMOD(ABS(YH)+0.000001,0.1**J).LT.0.000002) THEN GO TO 42 END IF 41 CONTINUE 42 CONTINUE IF(J.LE.0) THEN J=-1 END IF ccc IF(J.LE.0.OR.IY.GT.0) THEN ccc J=IY ccc END IF C J is the number of decimal places. C B = -( AINT(ALOG10(ABS(Y1)+.5))+2.785+FLOAT(J) )*HEIGHT IF(Y1.LT.0.) THEN B=B-HEIGHT END IF CALL NUMBER(B,Y-HEIGH2,HEIGHT,Y1,0.,J) END IF END IF IF (I.NE.0) THEN CALL PLOT(0.,Y,3) CALL PLOT(A,Y ,2) END IF Y = Y+YD 45 CONTINUE IF (IY.NE.0) THEN A = -HEIGHT*FLOAT(LY)-1.285*HEIGHT IF (YM-FLOAT(MY/NY).GE.0.25) THEN B = YY-HEIGHT ELSE B = YY-YY/YM/2.-HEIGH2 END IF CALL SYMBOL(A,B,HEIGHT,KY,0.,LY) END IF C C Writing texts: CALL SYMBOL(HEIGH0,YY+HEIGH2,HEIGHT,K1,0.,L1) B = HEIGHT*FLOAT(L1)+1.215*HEIGHT IF(M1.GT.0) THEN CALL NUMBER(B,YY+HEIGH2,HEIGHT,A1,0.,N1) END IF B = XX-HEIGHT*FLOAT(L2+M2)-.785*HEIGHT CALL SYMBOL( B,YY+HEIGH2,HEIGHT,K2,0.,L2) B = XX-HEIGHT*FLOAT(M2) +.215*HEIGHT IF(M2.GT.0) THEN CALL NUMBER(B,YY+HEIGH2,HEIGHT,A2,0.,N2) END IF CALL SYMBOL(HEIGH0,-5.3*HEIGHT,HEIGHT,K3,0.,L3) CALL SYMBOL(HEIGH0,-7.0*HEIGHT,HEIGHT,K4,0.,L4) C RETURN END C C======================================================================= C INCLUDE 'error.for' C error.for INCLUDE 'sep.for' C sep.for INCLUDE 'gse.for' C gse.for INCLUDE 'length.for' C length.for INCLUDE 'calcops.for' C calcops.for C C======================================================================= C