C Program 'sp.for' (Seismogram Plotting) to plot seimograms previously C stored in the GSE data exchange format. C C Version: 5.00 C Date: 1996, September 30 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 This Fortran77 file consists of the following external procedures: C SP... Main program to read and plot the seismograms. C FRAME...Subroutine called by the main to plot the rectangular C frame around the seismograms and supplement it with simple C descriptions. C C Other external procedures required: C RGSE1,RGSE2... Subroutines of the Fortran 77 file 'gse.for' C (package MODEL), designed to read seimograms 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 Description of the data files: C C The data are read in by the list directed input (free format). C In the description of data files, each numbered paragraph indicates C the beginning of a new input operation (new READ statement). C If the symbolic name of the input variable is enclosed in apostrophes, C the corresponding value in input data is of the type CHARACTER, i.e. C it should be a character string enclosed in apostrophes. If the first C letter of the symbolic name is I-N, the corresponding value is of the C type INTEGER. Otherwise, the input parameter is of the type REAL and C may or may not contain a decimal point. C C Input data read from the * external unit: C The interactive * external unit may also be redirected to the file C containing the relevant data. C (1) 'SPDAT','RF',KOMP,'SRC',/ C 'SPDAT'... String with the name of the input data file containing C the data describing how to plot the seismograms (e.g., C scale, minimum and maximum values of the coordinates, C etc.). C 'RF'... String with the name of the input data file in the GSE C data exchange format, containing the seismograms to be C plotted stored . In the present version, all seismograms C must correspond to a single component. Other components C must be stored in separate files. C KOMP... Component of the seismograms to be plotted. C 'SRC'.. String with the name of the input data file containing the C coordinates of the hypocentre. If left blank (default), C 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 Default: 'SPDAT'='sp.dat', 'RF'='rf.out', KOMP=3, 'SRC'=' '. C C Input file 'SPDAT' with the data to control plotting: C (1) I1,J1,K1,I2,J2,K2,I3,J3,K3,I4,J4,K4,I5,J5,K5,I6,J6,K6,/ C Selection of seismograms to be plotted. C First, the I1th to J1th seismograms of input file 'RF', with C increment K1, are plotted. Then I2th to J2th seismograms with C increment K2, and so on, until zero Ki is encountered. C The most frequent form of this input line will be C 1,N,1,/ C where N is the number of seismograms located in input file 'RF'. C (2) TMIN,TMAX,TAXIS,TM,TN,CODE,VRED,/ C Time axis (vertical). C TMIN,TMAX... Minimum and maximum time (or reduced time), i.e. C the times corresponding to the bottom and top of the C seismogram plot. C TAXIS...Length of the vertical time axis in cm. C TM... ABS(TM) is the number of intervals along the time axis, C starting at the bottom. Must be TM.NE.0. C TM.GT.0.: the marks at the endpoints of intervals will C be supplemented with corresponding times (or reduced C times). C TM.LT.0.: the time axis will have no description. C TN... Number of subintervals to be marked in each interval. C Must be TN.GT.0. C CODE... Specifies the distribution and description of seismograms C in the plot. See the description of XMIN,XMAX,YMIN,YMAX C below. C VRED... Reduction velocity. If non-zero, the time at each C receiver is reduced by the value of RR/VRED, where RR is C the hypocentral distance. C No time reduction is applied if VRED=0. C Default: TMIN=0., TMAX=1., TAXIS=10., TM=1., TN=1., CODE=0., C VRED=0. C (3) XMIN,XMAX,XAXIS,XM,XN,YMIN,YMAX,HEIGHT,/ C Distance axis (horizontal). C XMIN,XMAX,YMIN,YMAX: C For CODE=0.: Horizonal axis is divided into MR+1 C intervals, where MR 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, XMIN,XMAX,YMIN,YMAX have no meaning. C For CODE=1.: Horizonal axis represents the profile with C endpoints (X1,X2)=(XMIN,YMIN) and (X1,X2)=(XMAX,YMAX), C situated in a horizontal plane. The seismograms are C located at the orthogonal projections of the receivers C onto the profile. If XM.GT.0., the horizontal axis is C supplemented with the values of the X1 coordinates. C For CODE=2.: Horizonal axis represents the profile with C endpoints (X1,X2)=(XMIN,YMIN) and (X1,X2)=(XMAX,YMAX), C situated in a horizontal plane. The seismograms are C located at the orthogonal projections of the receivers C onto the profile. If XM.GT.0., the horizontal axis is C supplemented with the values of both X1 and X2 C coordinates. C For CODE=3.: Horizonal axis represents the vertical C profile with endpoints X3=XMIN and X3=XMAX. C The seismograms are located at the horizontal C projections of the receivers onto the profile. C If XM.GT.0., the horizontal axis is supplemented with C the values of X3 coordinates. C YMIN and YMAX have no meaning. C For CODE=4.: Horizonal axis represents the hypocentral C distance with endpoints RR=XMIN and RR=XMAX. C The seismograms are located according to the hypocentral C distances the receivers. C If XM.GT.0., the horizontal axis is supplemented with C the values of the hypocentral distance. C YMIN and YMAX have no meaning. C XAXIS...Length of the horizontal axis in cm. C XM... ABS(XM) is the number of intervals along the horizontal C axis, starting at the left. Must be XM.NE.0. C XM.GT.0.: the marks at the endpoints of intervals will C be supplemented with the corresponding values. C XM.LT.0.: the horizontal axis will have no description. C XN... Number of subintervals to be marked in each interval. C Must be XN.GT.0. C HEIGHT..Height of the characters in cm. C Default: TMIN=0., TMAX=1., TAXIS=FLOAT(MR+1), TM=1., TN=1., C CODE=0., VRED=0., HEIGHT=0.4, where MR is the number of C seismograms to be plotted, see (1). C (4) AMP,B1,EPIC,EPS,/ C AMP... AMP=0: Maximum amplitude at each trace normalized to the C fixed value. C Otherwise: All seismograms have the same scaling or power C scaling. C B1,EPIC,EPS... Amplitude scaling parameters. EPIC and EPS are C not used if AMP=0. C Amplitude scale AA is for AMP=0: C AA=B1*XD/AMAX C where XD is the average distance between seismograms, C and AMAX is the maximum amplitude in each seismogram. C Otherwise: C AA=B1*(RR/EPIC)**EPS C where EPIC is the hypocentral distance of each receiver. C Default: AMP=0., B1 =1., EPIC=1., EPS=0. C C Input file 'RF' with the seismograms to plot: C File in the GSE data exchange format, see the description in file C 'gse.for'. In the present version, all seismograms must C correspond to a single component. Other components must be stored C in separate files. 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 '. 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 C======================================================================= C C Filenames: CHARACTER*80 FSPDAT,FILERF,FILSRC C PARAMETER (MPTS=8196) CHARACTER*1 STAT,CHAN CHARACTER*80 TEXT INTEGER KR(18),NR(50) REAL SS(MPTS) 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) CHARACTER*4 TCOM(NCOM) REAL VCOM(NCOM) DATA TCOM/'XS1 ','XS2 ','XS3 '/ C DATA KX /0,1,1,1,19/ DATA KXTEXT/' ','X1','X1','X3','HYPOCENTRAL DISTANCE'/ C C....................................................................... C C Opening data files: FSPDAT='sp.dat' FILERF='rf.out' KOMP=3 FILSRC=' ' WRITE(*,'(A)') ' Enter ''sp.dat'', ''rf.out'', KOMP, ''SRC'': ' READ(*,*) FSPDAT,FILERF,KOMP,FILSRC VCOM(1)=0. VCOM(2)=0. VCOM(3)=0. IF(FILSRC.NE.' ') THEN OPEN(5,FILE=FILSRC,STATUS='OLD') READ(5,*) (TEXT,I=1,20) READ(5,*) STAT,VCOM(1),VCOM(2),VCOM(3) CLOSE(5) END IF OPEN(5,FILE=FSPDAT,STATUS='OLD') OPEN(4,FILE=FILERF,STATUS='OLD') C C Selection of receivers: DO 1 I=1,18 KR(I)=0 1 CONTINUE READ (5,*) (KR(I),I=1,18) MR= 0 DO 3 K=1,16,3 L = KR(K) M = KR(K+1) N = KR(K+2) IF(N.EQ.0) THEN GO TO 10 END IF DO 2 J=L,M,N MR= MR+1 NR(MR)= J 2 CONTINUE 3 CONTINUE 10 CONTINUE C C C Initial values for plotting frame: TMIN=0. TMAX=1. TAXIS=10. TM=1. TN=1. CODE=0. VRED=0. READ (5,*) TMIN,TMAX,TAXIS,TM,TN,CODE,VRED XMIN=0. XMAX=1. XAXIS=FLOAT(MR+1) XM=1. XN=1. YMIN=0. YMAX=0. HEIGHT=0.4 READ (5,*) XMIN,XMAX,XAXIS,XM,XN,YMIN,YMAX,HEIGHT KOD= INT(CODE+0.5) SCY= TAXIS/(TMAX-TMIN) XD = XAXIS/FLOAT(MR+1) XA = XMAX-XMIN YA = YMAX-YMIN C C Plotting frame: CALL RGSE1(4,TEXT) IX = KOD IT = -1 IF(KOD.GE.3) IX=1 IF(XM.LT.0.) IX=0 IF(TM.LT.0.) IT=0 KX1= KX(KOD+1)+1 CALL FRAME(XAXIS,TAXIS,ABS(XM),XN,ABS(TM),TN,0,IX,IT,HEIGHT, * XMIN,XMAX,KX1,KXTEXT(KOD+1),YMIN,YMAX,2,'X2', * TMIN,TMAX,4,'TIME', * 1,' ',0,0.,0,1,' ',0,0.,0,80,TEXT,1,' ') C C Amplitude scaling: AMP=0. B1 =1. EPIC=1. EPS=0. READ (5,*) AMP,B1,EPIC,EPS IF(B1.EQ.0.) B1=1. IF(EPIC.LE.0.) EPIC=1. C C C Plotting seismograms: JR= 0 DO 39 IR=1,MR 31 CONTINUE JR= JR+1 32 CONTINUE CALL RGSE2(4,STAT,CHAN,I,XR,YR,ZR,T,TD,NPTS,MPTS,SS, * NCOM,TCOM,VCOM) IF(I.NE.KOMP) GO TO 32 IF(JR.LT.NR(IR)) GO TO 31 C RR= SQRT((XR-VCOM(1))**2+(YR-VCOM(2))**2+(ZR-VCOM(3))**2) IF(VRED.NE.0.) THEN T=T-RR/VRED END IF AA= B1*XD/999. IF(AMP.NE.0..AND.EPS.EQ.0.) AA=B1 IF(AMP.NE.0..AND.EPS.NE.0.) AA=B1*((RR/EPIC)**EPS) X = XD*FLOAT(IR) IF(KOD.EQ.1.OR.KOD.EQ.2) * X = XAXIS*((XR-XMIN)*XA+(YR-YMIN)*YA)/(XA*XA+YA*YA) IF(KOD.EQ.3) X= XAXIS*(ZR-XMIN)/XA IF(KOD.EQ.4) X= XAXIS*(RR-XMIN)/XA Y = 0. CALL PLOT(X,Y,3) C DO 34 I=1,NPTS IF (TMIN.LE.T.AND.T.LE.TMAX) THEN A = X-AA*SS(I) Y = SCY*(T-TMIN) CALL PLOT(A,Y,2) END IF T = T+TD 34 CONTINUE C Y = TAXIS CALL PLOT(X,Y,2) 39 CONTINUE C CALL PLOT(0.,0.,999) WRITE(*,'(A)') '+SP: done. ' STOP END 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 LX1= MOD(MOD(KX1(I0),256)+256,256)+1 C LX2= MOD(MOD(KX2(I0),256)+256,256)+1 C LY = MOD(MOD(KY (I0),256)+256,256)+1 C L1 = MOD(MOD(K1 (I0),256)+256,256)+1 C L2 = MOD(MOD(K2 (I0),256)+256,256)+1 C L3 = MOD(MOD(K3 (I0),256)+256,256)+1 C L4 = MOD(MOD(K4 (I0),256)+256,256)+1 C C Initialization of CALCOMP: CALL PLOTS(0,0,0) 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 Portait 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 B = YY-YY/YM/2.-HEIGH2 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 'gse.for' INCLUDE 'calcops.for' C C======================================================================= C