C
C Program SS (Synthetic Seismograms) to read or generate and filter the C source time function. It may store the filtered source time function C and its Hilbert transform in the GSE data exchange format, or read the C frequency-domain response function and generate synthetic seismograms C in the GSE data exchange format. C C Version: 5.30 C Date: 1999, May 29 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 SS... Main program. C SS C SIGNAL..Subroutine to generate Gabor or Mueller signal of C given parameters. C SIGNAL C PLSIG...Subroutine to determine the maximum amplitude of the given C function, detect zeros beyond and behind the signal, write C the parameters of the signal to the output log file, and C eventually plot the signal. C PLSIG C PLOPN...Simple subroutine to initiate plotting. C PLOPN C PLTIM...Simple subroutine to supplement the signal plots with the C time labels. C PLTIM C FCOOLR..Subroutine for the fast Fourier transform of N=2**K C complex data points. C FCOOLR C C Other external procedures required: C WGSE1,WGSE2,WGSE3... Subroutines of the Fortran 77 file 'gse.for' C (package MODEL), designed to write 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 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) 'SEP',/ C 'SEP'...String with the name of the input SEP data file C containing the data describing the parameters of the C source time function and frequency-domain filter. C Description of file SEP C Default: 'SEP'='ss.h'. C C C Input file 'SEP': C It has the form of SEP parameter file. C Name of the input file: C RF='string'... String with the name of the input data file with C the frequency-domain response functions at individual C receivers. Is not taken into account if SS=' '. C The file is generated by program 'greenss.for'. C Description of file RF C Default: RF='rf.out' (mostly sufficient) C Names of output files: C SIGGSE='string'... String with the name of the output data file in C the GSE data exchange format, containing the filtered C source time function and its Hilbert transform. C Generated just if SIGGSE.NE.' '. C For the GSE data exchange format, refer to the description C in file 'gse.for'. C Default: SIGGSE=' ' C SS='string'... String with the name of the output data file in the C GSE data exchange format, containing the seismograms at C the receivers. C Generated just if SS.NE.' '. C For the GSE data exchange format, refer to the description C in file 'gse.for'. C Default: SS='ss.gse' C SSLOG='string'... String with the name of the output log file of C program SS. C Do not specify blank name. C Description of file SSLOG C Default: SSLOG='sslog.out' (mostly sufficient) C SIGPLOT='string'... String with the name of the output PostScript C file with the sketches of C 1 Input signal, C 2 Spectrum of the input signal, C 3 Spectrum of the filter, C 4 Spectrum of the filtered signal, C 5 Filtered signal, C 6 Hilbert transform of the filtered signal. C Generated just if SIGPLOT is not blank and MPTS is C positive. C Description of check plots C Default: SIGPLOT=' ' C SSPLOT='string'... String with the name of the output PostScript C file with the sketches of the synthetic seismograms at the C receivers. C Generated just if both SSPLOT and SS are not blank. C Description of check plots C Default: SSPLOT=' ' C Time step and time interval for the Fast Fourier Transform: C DT=real... Time interval to digitize the source time function and C seismograms. C Default: DT=1. C NFFT=integer... Number of the time samples for the fast Fourier C transform. C NFFT must be an integer power of 2 (NFFT is rounded up to C the nearest power of 2 in this program but not in other C programs which may work with it). C Default: NFFT=MSS, where MSS is the array dimension C declared in the code. C TRED=real, SSVRED=real... Specification of the time window for C the calculation of synthetic seismograms. Usually need C not be specified. C SSVRED.EQ.0: Seismogram is centred in the time interval of C length (NPTS-1)*DT according to the travel times of the C first and last arrivals at the receiver (default). C SSVRED.NE.0: Time of the first sample of the time window C is T=TRED+R/SSVRED, where R is the hypocentral distance. C Default: TRED=0., SSVRED=0. (mostly sufficient) C Data describing the filtration of the source time function: C DER=real, HILB=real... The source time function is DER-th C derivative and HILB-th Hilbert transform of the given C signal. C Default: DER=0. C Default: HILB=0. C FMIN=real, FLOW=real, FHIGH=real, FMAX=real... Parameters of the C frequency filter to be applied to the source time C function. The filter is zero C for frequencies F smaller than FMIN or greater than FMAX. C The filter is 1 between FLOW and FHIGH. C Between FMIN and FLOW, cosine tapering C ( 0.5-0.5*cos(pi*(F-FMIN)/(FLOW-FMIN) )**KEXP C is used. C Between FMIN and FLOW, cosine tapering C ( 0.5-0.5*cos(pi*(F-FMAX)/(FHIGH-FMAX) )**KEXP C is used. C Default: FMIN=0. C Default: FLOW=0. C Default: FMAX=0.5/DT C Default: FHIGH=FMAX C KEXP=integer... Exponent controlling the cosine frequency-domain C filter. Usually need not be specified because the default C is the most common option. C Default: KEXP=1 (mostly sufficient) C Input data to control optional plotting: C MPTS=integer... Number of points of the time functions at the C output check plot. The corresponding MPTS-1 time C intervals are all together scaled to 10.23cm. The signal C is approximately centred. MPTS does not apply to the C spectra at the output check plot, NFFT/2 points of each C spectrum is plotted. C MPTS=0: No output check plot of the source time function C is generated. C MPTS.LT.0: No output plot is generated, including the C seismograms. C Default: MPTS=0 C Data describing the source time function: C KSIG=integer... Type of the source time function: C KSIG=0: Digitized time function specified point by point. C KSIG=1: Gabor signal. C * KSIG=2: Hermite-Gaussian (Ricker) signal. C KSIG=3: Kuepper (Mueller) signal. C * KSIG=4: Rayleigh signal. C * KSIG=5: Generalized Berlage signal. C * KSIG=6: ? C * Only KSIG=0,1,3 is implemented in this version. C Default: KSIG=0 C SIGDIG=string... Name of the file containing the digitized time C function specified point by point and terminated by a C slash. The file is read by the list-directed input. C Required just for for KSIG=0 C Default: SIGDIG=' '. C SIGT=real... Reference time of the given signal. C Used for all signals. C Default: SIGT=0. (mostly sufficient) C SIGF=positive real... Reference frequency. C Required for all analytical signals (KSIG=1,2,3,4,5,6). C Default: SIGF=1. C SIGW=positive real... Relative width of the signal. C Often the width of the signal envelope expressed in the C reference half-periods 1/(2*SIGF). C Very roughly defined for non-causal functions. C Required for all analytical signals (KSIG=1,2,3,4,5,6). C Default: SIGW=4. C SIGA=real... Amplitude of the maximum of the signal or its C envelope. C Used for all analytical signals (KSIG=1,2,3,4,5,6). C Default: SIGA=1. (mostly sufficient) C SIGPH=real... Phase in radians. C Used for KSIG=1,3,4,5,6. C Default: SIGPH=0. C PAR5=real... C Used for KSIG=5,6. C PAR6=real... C Used for KSIG=5. C Detailed description of the source time functions: C KSIG=0: Time function digitized with step DT, specified point by C point and terminated by a slash is read from the file C given by input parameter SIGDIG. The first sample C corresponds to time SIGT. C KSIG=1: Gabor signal: C S(t)=SIGA*exp(-(2*pi*SIGF*(t-SIGT)/SIGW)**2) C *cos(2*pi*SIGF*(t-SIGT)+SIGPH) C Gabor signal is not causal but has all derivatives C continuous. C SIGF... Prevailing frequency. C SIGW... Relative width of the signal. In the interval of C SIGW half-periods 1/(2*SIGF) the signal envelope C exceeds 8.48 per cent of its maximum. C SIGA... Amplitude of the envelope. C SIGPH...Phase. C KSIG=2: Hermite-Gaussian (Ricker) signal: C * for future extension, not implemented in this version * C S(t)=A*(-1)**n*Hn(SIGF*(t-SIGT))*exp(-(SIGF*(t-SIGT))**2) C where Hn(x) is the Hermite polynomial of order n=SIGW, C (-1)**n*Hn(x)*exp(-x**2)=(d/dx)**n[exp(-x**2)] C and the scaling factor is C A=SIGA*(n/2)!/n! C Then C S(t)=A*(1/SIGF*d/dt)**n[exp(-(SIGF*(t-SIGT))**2)] C Note that Ricker suggested the Hermite-Gaussian signal C with n=1, 2 and 3 for the approximation of particle C diplacement, velocity and acceleration, respectively, C at the receivers very distant from a point source. C Most widely used is the special case for SIGW=2, C originally suggested by Ricker for particle velocity, C S(t)=SIGA*(2*(SIGF*(t-SIGT))**2-1)*exp(-(SIGF*(t-SIGT))**2) C Hermite-Gaussian signal is not causal but has all C derivatives continuous. C SIGF... Reference frequency. The mean frequency is C F0=SIGF*(SIGW/2)!/(((SIGW-1)/2))!/pi and the C spectrum has maximum at FM=SIGF*SQRT(SIGW/2)/pi. C Note that for even SIGW C (SIGW/2)!/(((SIGW-1)/2))! C =2*4*...*SIGW/(1*3*...*(SIGW-1))/SQRT(pi) C and for odd SIGW C (SIGW/2)!/(((SIGW-1)/2))! C =3*5*...*SIGW/(2*4*...*(SIGW-1))*SQRT(pi)/2 C SIGW... Order n of the Hermite-Gaussian signal. The best C known option is SIGW=2 (Ricker signal for C far-field particle velocity). C SIGW is rounded to the nearest integer. C KSIG=3: Kuepper (Mueller) signal: C For 0.LE.(t-SIGT).LE.SIGW/(2*SIGF): C S(t)=SIGA*A*(sin(2*pi*SIGF*(t-SIGT)) C -sin(2*pi*SIGF*(t-SIGT)*(SIGW+2)/SIGW) C * SIGW / (SIGW+2) C -(1-cos(2*pi*SIGF*(t-SIGT)/SIGW)) C * sin(pi*SIGW) / (SIGW+2) ) C with A=1/SMAX, where SMAX is the maximum absolute value C of the signal without multiplication factors SIGA*A. C Outside the interval: C S(t)=0 C Note that the standard Kuepper signal is defined just for C integer SIGW and consists only of the two sine terms. C The last term (with cosine) has been added to make the C signal continuous also for non-integer SIGW. C Kuepper signal is causal, has continuous the first C derivative, and for integer SIGW also the second C derivative. C SIGF... Reference frequency. C SIGW... Relative width of the signal, i.e. the duration of C the signal in half-periods 1/(2*SIGF). C SIGW rounded to the next greater integer is the C number of local extrema of the signal. C SIGA... Maximum amplitude of the signal. Determines A. C SIGPH...Not applicable. C KSIG=4: Rayleigh signal: C * for future extension, not implemented in this version * C S(t)=SIGA*(COS(SIGPH)-SIN(SIGPH)*SIGF*(t-SIGT)) C /(1+(SIGF*(t-SIGT))**2) C Note that the Rayleigh signal is the Dirac distribution C D(t-SIGT) shifted by i/SIGF in the complex plane, C S(t)=Re(exp(i*SIGPH)*D(t-SIGT+i/SIGF)) C Rayleigh signal is not causal but has all derivatives C continuous. C SIGF... Reference frequency. C SIGW... Not applicable. C SIGA... Amplitude of the envelope. C SIGPH...Phase. C KSIG=5: Generalized Berlage signal: C * for future extension, not implemented in this version * C For t.LE.SIGT: C S(t)=0 C For t.GT.SIGT: C S(t)=SIGA*A(t)*(t-SIGT)**SIGW*exp(-PAR5*(t-SIGT)) C *sin(2*pi*SIGF*(t-SIGT)+SIGPH) C where A(t) normalizes the maximum of the envelope to 1 C and also enables to generalize the Berlage signal by C means of choosing PAR6 positive, C A(t)=(SIGW/(PAR5*(1+SIGW*PAR6*(t-SIGT))*TMAX**2) C *exp(PAR5*TMAX) C where time t=SIGT+TMAX corresponds to the maximum of the C envelope, C TMAX=(SQRT(1+4*SIGW*SIGW*PAR6/PAR5)-1)/(2*SIGW*PAR6) C Note that the limiting case for PAR6=0 is the standard C Berlage signal, with A(t) constant, C A(t)=PAR5/SIGW*exp(SIGW) C TMAX=SIGW/PAR5 C Note that the limiting case for SIGW=+infinity is given by C S(t)=SIGA C *exp(-1/PAR6/(t-SIGT)+2*SQRT(PAR5/PAR6)-PAR5*(t-SIGT)) C *sin(2*pi*SIGF*(t-SIGT)+SIGPH) C with maximum of the envelope in time t=SIGT+TMAX, C TMAX=1/SQRT(PAR5*PAR6) C The generalized Berlage signal is causal, has continuous C derivatives of orders less than INT(SIGW), C and for SIGPH=0 or SIGPH=pi even of orders less than C INT(SIGW)+1. C SIGF... Prevailing frequency. C SIGW... Controls the onset of the signal. Derivatives of C of orders less than INT(SIGW), and for SIGPH=0 or C SIGPH=pi even of orders less than INT(SIGW)+1, are C continuous. For PAR6=0 (standard Berlage), C distance TMAX beetween the begining of the signal C and the maximum of the envelope is SIGW/PAR5. C SIGW=999 or greater is understood as infinity. C SIGA... Amplitude of the envelope. C SIGPH...Phase. C PAR5... Controls the effective duration of the signal. C PAR6... Zero for standard Berlage signal. Positive values C enable to decrease distance TMAX beetween the C begining of the signal and the maximum of the C envelope. To decrease TMAX from SIGW/PAR5 to C K/PAR5, select PAR6=PAR5*(1/K-1/SIGW)/K. C KSIG=6: ? C * for future extension, not implemented in this version * C For 0.LE.(t-SIGT): C S(t)=SIGA*exp(-(1/TRED-2+TRED)*pi*PAR5/SIGW) C *sin(2*pi*SIGF*(t-SIGT)+SIGPH) C Where: C TRED=4*SIGF*(t-SIGT)/PAR5 C Otherwise: C S(t)=0 C SIGF... Prevailing frequency. C SIGW... Controls the relative width of the signal. C SIGA... Amplitude of the envelope. C SIGPH...Phase. C PAR5... Distance beetween the begining of the signal and C the maximum of the envelope is SIGW/4 prevailing C periods 1/SIGF. C Reasonable value for SIGPH=0 is PAR5=3. C C C Input file 'RF' with the response functions: C (1) 'TEXT1',/ C 'TEXT1'... Text string in apostrophes. The first 34 characters C will be passed to the header of the output GSE file. C (2) XS1,XS2,XS3,/ C XS1,XS2,XS3... Coordinates of the hypocentre. C (3) FMINIM,FD,NF,/ C FMINIM..The lowest frequency at the digitized response function. C FD... The frequency step. C NF... The number of discrete frequencies. C (4) For each receiver (4.1) and (4.2): C (4.1) 'REC',X1,X2,X3,TMIN,TMAX,AMAX,/ C 'REC'...Name of the receiver (6 characters). C X1,X2,X3... Coordinates of the receiver. C TMIN,TMAX...Travel times of the first and last arrivals at the C receiver. C AMAX... Maximum absolute value over the real and imaginary part of C all 3 components of the response function. C (4.2) Written only if TMIN.LE.TMAX (Attention: Formatted input!): C 6*NF integer numbers, FORMAT(12I6): C (IPR(I,IF),I=1,6,IF=1,NF) C IPR... IPR(I,IF)=IFIX(99999.1*SPR(I,IF)/AMAX), where C SPR(*,IF) is the response function at the frequency C F=FMIN+(IF-1)*FD. C SPR(1,IF) is the real part of the X1 component. C SPR(2,IF) is the imaginary part of the X1 component. C SPR(3,IF) is the real part of the X2 component. C SPR(4,IF) is the imaginary part of the X2 component. C SPR(5,IF) is the real part of the X3 component. C SPR(6,IF) is the imaginary part of the X3 component. C (5) / or end of file. C C C Output file 'SSGSE' with the seismograms or source time function: C File in the GSE data exchange format, see the description in file C 'gse.for'. C Description of format GSE C C C Output log file 'SSLOG': C Contains information on the input data, source time function, C synthetic seismograms. Often may be discarded. C C C Output check plots: C File SIGPLOT (if MPTS.GT.0) contains plots of: C 1 Input signal, C 2 Spectrum of the input signal, C 3 Spectrum of the filter, C 4 Spectrum of the filtered signal, C 5 Filtered signal, C 6 Hilbert transform of the filtered signal. C File SSPLOT (if SS.NE.' ') contains plots of the synthetic C seismograms at the receivers. C Horizontal size of each function is 10.23cm, vertical scaling is C 1cm per the maximum absolute amplitude of the function. C MPTS points are plotted for input signal, filtered signal and its C Hilbert transform. C NPTS/2 points are plotted for the spectra. C NPTS points are plotted for the synthetic seismograms. C C======================================================================= C C C EXTERNAL RSEP1,RSEP3T,RSEP3I,RSEP3R,WSEPR,WGSE1,WGSE2,WGSE2C,WGSE3 EXTERNAL ERROR,SIGNAL,PLSIG,SYMBOL,FCOOLR C C Common block /RAMC/: INCLUDE 'ram.inc' C ram.inc C C Allocation of working arrays: INTEGER MSS PARAMETER (MSS=MRAM/9) REAL S1(2,MSS),S2(2,MSS),S3(2,MSS),S4(2,MSS),SS(MSS) EQUIVALENCE (S1,RAM ) EQUIVALENCE (S2,RAM(2*MSS+1)) EQUIVALENCE (S3,RAM(4*MSS+1)) EQUIVALENCE (S4,RAM(6*MSS+1)) EQUIVALENCE (SS,RAM(8*MSS+1)) C C----------------------------------------------------------------------- C C Filenames: CHARACTER*80 FILDAT,FILLOG,FILRF,FILGSE,FILSIG,FILPS PARAMETER (LU4=1,LU5=2,LU6=3,LU7=4) C PARAMETER (UNDEF=-999999.) CHARACTER*80 TEXT1 REAL PSGNL(10) C C Receiver name: CHARACTER*6 REC C C Source coordinates transferred through the GSE file: INTEGER NCOM PARAMETER (NCOM=3) CHARACTER*80 LINE CHARACTER*3 TCOM(NCOM) REAL VCOM(NCOM) DATA TCOM/'XS1','XS2','XS3'/ C C....................................................................... C C Opening data files: FILDAT='ss.h' WRITE(*,'(A)') '+SS: Enter input filename: ' READ(*,*) FILDAT WRITE(*,'(A)') '+SS: Working... ' C CALL RSEP1(LU5,FILDAT) CALL RSEP3T('SSLOG',FILLOG,'sslog.out') OPEN(LU6,FILE=FILLOG) C C Input data for the source signal CALL RSEP3I('KSIG ',KSGNL ,0) CALL RSEP3I('KEXP ',KEXP ,1) CALL RSEP3I('MPTS ',MPTS ,0) CALL RSEP3I('NFFT ',NPTS ,MSS) CALL RSEP3R('DER ',DER ,0.) CALL RSEP3R('HILB ',HILB ,0.) CALL RSEP3R('DT ',DT ,1.) CALL RSEP3R('FMIN ',FMIN ,0.) CALL RSEP3R('FLOW ',FL ,0.) CALL RSEP3R('FHIGH',FR ,FMAX) CALL RSEP3R('FMAX ',FMAX ,0.5/DT) CALL RSEP3R('SIGT ',SIGT ,0.) CALL RSEP3R('TRED ',TRED ,0.) CALL RSEP3R('SSVRED',VRED ,0.) N = NPTS-1 NPTS= 1 DO 1 KPTS=1,15 NPTS= NPTS+NPTS IF (N-1.LE.0) GO TO 10 N= N/2 1 CONTINUE 10 CONTINUE IF (NPTS.GT.MSS) THEN C SS-01 CALL ERROR('SS-01: Array dimension MSS less than NPTS.') END IF IF (KSGNL.LE.0) THEN CALL RSEP3T('SIGDIG',FILSIG,' ') IF (FILSIG.EQ.' ') THEN C SS-02 CALL ERROR('SS-02: No filename with source time function') C For input parameter KSIG=0 (default value), input parameter C SIGDIG must contain the name of the file with the digitized C source time function. Refer to the description of the C input data. END IF OPEN(LU5,FILE=FILSIG,STATUS='OLD') DO 12 I=1,NPTS S1(1,I)=0. 12 CONTINUE READ (LU5,*) (S1(1,I),I=1,NPTS) CLOSE(LU5) DO 13 I=NPTS,1,-1 IF(S1(1,I).NE.0.) GO TO 14 S1(1,I)=0. 13 CONTINUE 14 CONTINUE N = I N1= (NPTS-N)/2 N2= NPTS-N-N1 J = NPTS DO 15 I=1,N2 S1(1,J)= 0. J = J-1 15 CONTINUE DO 16 I=1,N K = J-N1 S1(1,J)= S1(1,K) J = J-1 16 CONTINUE DO 17 I=1,N1 S1(1,I)= 0. 17 CONTINUE SIGT= SIGT-DT*FLOAT(N1) ELSE CALL RSEP3R('SIGF' ,PSGNL(1),1.) CALL RSEP3R('SIGW' ,PSGNL(2),4.) CALL RSEP3R('SIGPH',PSGNL(3),0.) CALL RSEP3R('SIGA' ,PSGNL(4),1.) CALL RSEP3R('SIG5' ,PSGNL(5),0.) CALL RSEP3R('SIG6' ,PSGNL(6),0.) CALL SIGNAL(KSGNL,NPTS,SIGT,DT,S1,PSGNL) END IF DO 20 I=1,NPTS S1(2,I)= 0. 20 CONTINUE C C Plotting the input signal: CALL RSEP3T('SIGPLOT',FILPS,' ') WRITE(LU6,'(/A,I2)') ' Source signal No.',KSGNL WRITE(LU6,'(2A/2A)') ' * Left-hand Left-hand Right-hand', * ' Right-hand Non-zero Maximum ', * ' tip hill-side hill-side', * ' tip range amplitude' CALL PLSIG(MPTS,1,1,MPTS,NPTS,SIGT,DT,S1,A,I,J,FILPS) C C Spectrum of the input signal: CALL FCOOLR(KPTS,S1,1.) FD= 1./DT/FLOAT(NPTS) C C Amplitude spectrum, frequency window: A = 2./FLOAT(NPTS) DO 38 I=1,NPTS/2 S2(1,I)= SQRT(S1(1,I)*S1(1,I)+S1(2,I)*S1(2,I)) F = FD*FLOAT(I-1) IF (F-FMIN) 31,31,32 31 S2(2,I)= 0. GO TO 38 32 IF (F-FL) 33,34,34 33 S2(2,I)= A*(.5+.5*COS(3.14159*(F-FL)/(FMIN-FL)))**KEXP GO TO 38 34 IF (F-FR) 35,35,36 35 S2(2,I)= A GO TO 38 36 IF (F-FMAX) 37,31,31 37 S2(2,I)= A*(.5+.5*COS(3.14159*(F-FR)/(FMAX-FR)))**KEXP 38 CONTINUE IF (DER) 39,41,39 39 FDA= 6.283185*FD F = 0. DO 40 I=2,NPTS/2 F = F+FDA 40 S2(2,I)= S2(2,I)*F**DER 41 CONTINUE DO 42 I=NPTS/2+1,NPTS S2(1,I)= 0. 42 S2(2,I)= 0. CALL PLSIG(MPTS,2,2,NPTS/2,NPTS/2,0.,FD,S2,A,I,J,FILPS) CALL PLSIG(MPTS,3,3,NPTS/2,NPTS/2,0.,FD,S2(2,1),A,I,J,FILPS) C C Filtration of the input signal: A = 1.570796*(DER+HILB) C = COS(A) S = -SIN(A) DO 47 I=1,NPTS A = S1(1,I)*C-S1(2,I)*S S1(2,I)= S1(1,I)*S+S1(2,I)*C 47 S1(1,I)= A 48 DO 49 I=1,NPTS S1(1,I)= S1(1,I)*S2(2,I) S1(2,I)= S1(2,I)*S2(2,I) 49 S2(1,I)= S2(1,I)*S2(2,I) CALL PLSIG(MPTS,4,4,NPTS/2,NPTS/2,0.,FD,S2,A,NF1,NF2,FILPS) DO 50 I=1,NPTS S2(1,I)= S1(1,I) 50 S2(2,I)= S1(2,I) CALL FCOOLR(KPTS,S2,-1.) CALL PLSIG(MPTS,5,5,MPTS,NPTS,SIGT,DT,S2,A,N1,N2,FILPS) CALL PLSIG(MPTS,6,6,MPTS,NPTS,SIGT,DT,S2(2,1),A,N,J,FILPS) C Legend: IF(FILPS.NE.' '.AND.MPTS.GT.0) THEN CALL SYMBOL(-1.4,-3.,0.8,'1',0.,1) CALL SYMBOL( 0.0,-3.,0.3,'INPUT SIGNAL',0.,12) CALL SYMBOL( 5.6,-3.,0.3,'RIGHT:',0.,6) CALL SYMBOL( 7.6,-3.,0.4,'MAXIMUM AMPLITUDE',0.,17) CALL SYMBOL(-1.4,-4.,0.8,'2',0.,1) CALL SYMBOL( 0.0,-4.,0.3,'SPECTRUM OF THE INPUT SIGNAL',0.,28) CALL SYMBOL(-1.4,-5.,0.8,'3',0.,1) CALL SYMBOL( 0.0,-5.,0.3,'SPECTRUM OF THE FILTER',0.,22) CALL SYMBOL(-1.4,-6.,0.8,'4',0.,1) CALL SYMBOL( 0.0,-6.,0.3, * 'SPECTRUM OF THE FILTERED SIGNAL',0.,31) CALL SYMBOL(-1.4,-7.,0.8,'5',0.,1) CALL SYMBOL( 0.0,-7.,0.3,'FILTERED SIGNAL',0.,16) CALL SYMBOL(-1.4,-8.,0.8,'6',0.,1) CALL SYMBOL( 0.0,-8.,0.3, * 'HILBERT TRANSFORM OF THE FILTERED SIGNAL',0.,40) END IF CALL RSEP3T('SIGGSE',FILSIG,' ') IF(FILSIG.NE.' ') THEN C Writing the source time function and its Hilbert transform: OPEN(LU7,FILE=FILSIG) CALL WGSE1(LU7,' ') DO 51 I=1,NPTS SS(I)=S2(1,I) 51 CONTINUE CALL WGSE2(LU7,' ',' ',0,0.,0.,0.,SIGT,DT,NPTS,SS) DO 52 I=1,NPTS SS(I)=S2(2,I) 52 CONTINUE CALL WGSE2(LU7,' ',' ',0,0.,0.,0.,SIGT,DT,NPTS,SS) CALL WGSE3(LU7) CLOSE(LU7) WRITE(*,'(A)') '+SS: Source time function generated. ' STOP END IF IF(FILPS.NE.' '.AND.MPTS.GT.0) CALL PLOT(0.,0.,999) IF (N1.GE.NPTS.OR.N.GE.NPTS) THEN C SS-04 CALL ERROR * ('SS-04: Too small number NPTS of time samples for FFT.') END IF N1= MIN0(N1,I) N2= MIN0(N2,J) C C....................................................................... C C Opening output file with seismograms: CALL RSEP3T('SS',FILGSE,'ss.gse') IF (FILGSE.EQ.' ') THEN IF(FILSIG.EQ.' ') THEN WRITE(*,'(A)') '+SS: Nothing to do. ' ELSE WRITE(*,'(A)') '+SS: Signal generated. ' END IF STOP END IF OPEN(LU7,FILE=FILGSE) CALL RSEP3T('SSPLOT',FILPS,' ') C C Opening input files with the response function: CALL RSEP3T('RF',FILRF,'rf.out') IF (FILRF.EQ.' ') THEN C SS-03 CALL ERROR('SS-03: No file with response function specified') END IF OPEN(LU4,FILE=FILRF,STATUS='OLD') C C Headlines of files: WRITE(LU6,'(/A)') * ' Beginning of the input file with frequency response' READ (LU4,*) TEXT1 WRITE(LU6,'(2A)') ' ***',TEXT1 READ (LU4,*) (VCOM(I),I=1,3) WRITE(LU6,'(A,10F8.3)') ' ***',(VCOM(I),I=1,3) READ (LU4,*) FMINIM,A,NF WRITE(LU6,'(A,2E12.5,I4)') ' ***',FMINIM,A,NF IF (NPTS.NE.INT(1./A/DT+.5)) THEN C SS-06 CALL ERROR('SS-06: Inconsistent time and frequency steps.') END IF MINIM= INT(FMINIM/FD+1.5) MAXIM= MINIM+NF-1 IF (NF1+1.LT.MINIM) THEN C SS-07 CALL ERROR * ('SS-07: Missing low frequencies in response function.') END IF IF (MAXIM+NF2.LT.NPTS/2) THEN C SS-08 CALL ERROR * ('SS-08: Missing high frequencies in response function.') END IF CALL WGSE1(LU7,TEXT1) WRITE(LU6,'(/A)') ' Synthetic sections at receivers' WRITE(LU6,'(2A/2A)') ' * Coordinates of the receiver ', * ' First Last Upper ', * ' X Y Z ', * ' arrival arrival amplitude' WRITE(LU6,'(2A/2A)') ' * Left-hand Left-hand Right-han', * 'd Right-hand Non-zero Maximum ', * ' tip hill-side hill-sid', * 'e tip range amplitude' C C Preparing source coordinates for output in the GSE file: DO 55 I=1,NCOM CALL WSEPR(LINE,TCOM(I),VCOM(I)) CALL WGSE2C(LINE) 55 CONTINUE C C Loop for the receivers: NUMS= 1 DO 79 IREC=1,999999 X=UNDEF TMIN= 999999. TMAX=-999999. AMAX=0. READ (LU4,*,END=90) REC,X,Y,Z,TMIN,TMAX,AMAX IF(X.EQ.UNDEF) THEN GO TO 90 END IF WRITE(LU6,'(I4,5F11.3,E11.3)') IREC,X,Y,Z,TMIN,TMAX,AMAX IF(TMIN.LE.TMAX) THEN C Zero range in frequency domain: N = MINIM-1 DO 58 I=1,N S2(1,I)= 0. S2(2,I)= 0. S3(1,I)= 0. S3(2,I)= 0. S4(1,I)= 0. S4(2,I)= 0. 58 CONTINUE N = MIN0(NPTS,MAXIM+1) DO 59 I=N,NPTS S2(1,I)= 0. S2(2,I)= 0. S3(1,I)= 0. S3(2,I)= 0. S4(1,I)= 0. S4(2,I)= 0. 59 CONTINUE C READ(LU4,'(12F6.0)') (S2(1,I),S2(2,I),S3(1,I),S3(2,I), * S4(1,I),S4(2,I),I=MINIM,MAXIM) A = AMAX/99999. DO 65 I=MINIM,MAXIM B = A*(S1(1,I)*S2(1,I)-S1(2,I)*S2(2,I)) S2(2,I)= A*(S1(1,I)*S2(2,I)+S1(2,I)*S2(1,I)) S2(1,I)= B B = A*(S1(1,I)*S3(1,I)-S1(2,I)*S3(2,I)) S3(2,I)= A*(S1(1,I)*S3(2,I)+S1(2,I)*S3(1,I)) S3(1,I)= B B = A*(S1(1,I)*S4(1,I)-S1(2,I)*S4(2,I)) S4(2,I)= A*(S1(1,I)*S4(2,I)+S1(2,I)*S4(1,I)) S4(1,I)= B 65 CONTINUE CALL FCOOLR(KPTS,S2,-1.) CALL FCOOLR(KPTS,S3,-1.) CALL FCOOLR(KPTS,S4,-1.) N = INT((TMAX+TMIN)/(DT+DT)) IF(VRED.GT.0) N=NPTS/2+ * INT((TRED+SQRT((VCOM(1)-X)**2+(VCOM(2)-Y)**2 * +(VCOM(3)-Z)**2)/VRED)/DT) TMIN= SIGT+DT*FLOAT(N) N = MOD(N,NPTS) IF(N.LT.0) THEN N=N+NPTS END IF DO 66 I=1,N J = NPTS-N+I S2(2,J)= S2(1,I) S3(2,J)= S3(1,I) S4(2,J)= S4(1,I) 66 CONTINUE K = NPTS-N DO 68 I=1,K J = N+I S2(2,I)= S2(1,J) S3(2,I)= S3(1,J) S4(2,I)= S4(1,J) 68 CONTINUE CALL PLSIG * (MPTS+1,NUMS,IREC,NPTS,NPTS,TMIN,DT,S2(2,1),AMAX,N1,N2,FILPS) IF(N1.LT.NPTS) THEN C Non-zero signal NUMS= NUMS+1 T0= TMIN+DT*FLOAT(N1) N2= NPTS-N2 N = N2-N1 DO 71 I=1,N N1= N1+1 SS(I)=S2(2,N1) 71 CONTINUE CALL WGSE2(LU7,REC,' ',1,X,Y,Z,T0,DT,N,SS) ELSE C Zero signal CALL WGSE2(LU7,REC,' ',1,X,Y,Z,0.,DT,0,SS) END IF CALL PLSIG * (MPTS+1,NUMS,IREC,NPTS,NPTS,TMIN,DT,S3(2,1),AMAX,N1,N2,FILPS) IF(N1.LT.NPTS) THEN C Non-zero signal NUMS= NUMS+1 T0= TMIN+DT*FLOAT(N1) N2= NPTS-N2 N = N2-N1 DO 72 I=1,N N1= N1+1 SS(I)=S3(2,N1) 72 CONTINUE CALL WGSE2(LU7,REC,' ',2,X,Y,Z,T0,DT,N,SS) ELSE C Zero signal CALL WGSE2(LU7,REC,' ',2,X,Y,Z,0.,DT,0,SS) END IF CALL PLSIG * (MPTS+1,NUMS,IREC,NPTS,NPTS,TMIN,DT,S4(2,1),AMAX,N1,N2,FILPS) IF(N1.LT.NPTS) THEN C Non-zero signal NUMS= NUMS+1 T0= TMIN+DT*FLOAT(N1) N2= NPTS-N2 N = N2-N1 DO 73 I=1,N N1= N1+1 SS(I)=S4(2,N1) 73 CONTINUE CALL WGSE2(LU7,REC,' ',3,X,Y,Z,T0,DT,N,SS) ELSE C Zero signal CALL WGSE2(LU7,REC,' ',3,X,Y,Z,0.,DT,0,SS) END IF ELSE CALL WGSE2(LU7,REC,' ',1,X,Y,Z,0.,DT,0,SS) CALL WGSE2(LU7,REC,' ',2,X,Y,Z,0.,DT,0,SS) CALL WGSE2(LU7,REC,' ',3,X,Y,Z,0.,DT,0,SS) END IF 79 CONTINUE C C End of computation: 90 IF(FILPS.NE.' '.AND.MPTS.GT.-1.AND.NUMS.GT.1) CALL PLOT(0.,0.,999) CALL WGSE3(LU7) CLOSE(LU7) CLOSE(LU6) CLOSE(LU4) WRITE(*,'(A)') '+SS: Done. ' STOP END C C======================================================================= C C C SUBROUTINE SIGNAL(KSGNL,NPTS,SIGT,DT,S,PAR) REAL S(2,NPTS),PAR(*) C EXTERNAL ERROR C GO TO (10,1,30,1,1,1) KSGNL 1 CONTINUE C SS-09 CALL ERROR('SS-09: Only signals KSIG=1,3 allowed') C C Gabor signal 10 CONTINUE T = -DT*FLOAT(NPTS/2) SIGT= SIGT+T A = 6.283185*PAR(1) B = A*A/PAR(2)/PAR(2) DO 11 I=1,NPTS S(1,I)=0. IF(B*T*T.LT.70.) S(1,I)= COS(A*T+PAR(3))*EXP(-B*T*T) IF(PAR(4).NE.0.) S(1,I)=S(1,I)*PAR(4) T = T+DT 11 CONTINUE RETURN C C Kuepper (Mueller) signal: 30 CONTINUE N2= IFIX(PAR(2)/PAR(1)/DT/2.)+1 N1= (NPTS-N2)/2 SIGT= SIGT-DT*FLOAT(N1) A = 6.283185*PAR(1) B = PAR(2)/(2.+PAR(2)) C = A/B D = SIN(3.141593*PAR(2))/(2.+PAR(2)) E = A/PAR(2) DO 31 I=1,N1 S(1,I)= 0. 31 CONTINUE T = 0. F = 0. N2= N1+N2 N1= N1+1 DO 32 I=N1,N2 S(1,I)= SIN(A*T)-B*SIN(C*T)+D*COS(E*T)-D F = AMAX1(F,ABS(S(1,I))) T = T+DT 32 CONTINUE IF(PAR(4).NE.0.) F=F/PAR(4) DO 33 I=N1,N2 S(1,I)= S(1,I)/F 33 CONTINUE N2= N2+1 DO 34 I=N2,NPTS S(1,I)= 0. 34 CONTINUE RETURN C C Generalized Berlage signal: 50 CONTINUE N2= IFIX(ALOG(1000.)/PAR(5)/DT)+1 N1= (NPTS-N2)/2 SIGT= SIGT-DT*FLOAT(N1) A = 6.283185*PAR(1) B=2.*SQRT(PAR(5)/PAR(6)) TMAX = 2.*PAR(2)*PAR(2)*PAR(6)/PAR(5) IF(TMAX.LT.0.000001) THEN TMAX=1. ELSE TMAX=(SQRT(1.+2.*B)-1.)/B ENDIF TMAX=TMAX*PAR(2)/PAR(5) DO 51 I=1,N1+1 S(1,I)= 0. 51 CONTINUE T = 0. F = 0. N1= N1+1 DO 52 I=N1+1,NPTS T=T+DT S(1,I)=0. IF(PAR(2).LE.998.) THEN TRED=1.+PAR(2)*PAR(6)*T*TMAX*TMAX TRED=PAR(2)*T/(PAR(5)*TRED) C=PAR(5)*T IF(C.LT.70.) S(1,I)= SIN(A*T+PAR(3))*EXP(-C)*TRED**PAR(2) ELSE IF(PAR(6).LE.0.) THEN C SS-10 CALL ERROR('SS-10: Signal parameter PAR6 not positive') ENDIF B=2.*SQRT(PAR(5)/PAR(6)) C=1./(PAR(6)*T)-B+PAR(5)*T IF(C.LT.70.) S(1,I)= SIN(A*T+PAR(3))*EXP(-C) ENDIF IF(PAR(4).NE.0.) S(1,I)=S(1,I)*PAR(4) 52 CONTINUE RETURN C C Signal No.6: 60 CONTINUE N2= IFIX(PAR(2)/PAR(1)/DT/2.)+1 N1= (NPTS-N2)/2 SIGT= SIGT-DT*FLOAT(N1) A = 6.283185*PAR(1) B = 4.*PAR(1)/PAR(5) DO 61 I=1,N1+1 S(1,I)= 0. 61 CONTINUE T = 0. F = 0. N1= N1+1 DO 62 I=N1+1,NPTS T=T+DT S(1,I)=0. TRED=B*T C=1./TRED-2.+TRED C=C*3.141593*PAR(5)/PAR(2) IF(C.LT.70.) S(1,I)= SIN(A*T+PAR(3))*EXP(-C) IF(PAR(4).NE.0.) S(1,I)=S(1,I)*PAR(4) 62 CONTINUE RETURN C END C C======================================================================= C C C SUBROUTINE PLSIG(KPLOT,NUMS,NUM,MPTS,NPTS,TL,DT,S,AMP,N1,N2,FILPS) REAL S(2,NPTS) CHARACTER*(*) FILPS C EXTERNAL ERROR,PLTIM,PLOTN,PLOT,NUMBER PARAMETER (LU6=3) C IF(MPTS.GT.NPTS) THEN C SS-11 WRITE(*,'(2(A,I6))') ' MPTS=',MPTS,' NPTS=',NPTS CALL ERROR('SS-11: MPTS greater than NPTS') END IF C C Maximum amplitude: AMP= 0.000001 DO 1 I=1,NPTS 1 AMP= AMAX1(AMP,ABS(S(1,I))) EPS= 0.002*AMP C C Zeros beyond and behind the signal: DO 2 N1=1,NPTS IF(ABS(S(1,N1)).GT.EPS) GO TO 3 2 CONTINUE N1= NPTS RETURN 3 N1= N1-1 DO 4 N2=1,NPTS I = NPTS-N2+1 IF(ABS(S(1,I)).GT.EPS) GO TO 5 4 CONTINUE 5 N2= N2-1 C C Writing the parameters of the signal: N3= (NPTS-MPTS)/2+1 N4= N3+MPTS-1 T1= TL+DT*FLOAT(N1) T2= TL+DT*FLOAT(NPTS-N2-1) T3= TL+DT*FLOAT(N3-1) T4= TL+DT*FLOAT(N4-1) A = T2-T1 WRITE(LU6,'(I4,5F11.3,E11.3)') NUM,T3,T1,T2,T4,A,AMP IF(FILPS.EQ.' '.OR.KPLOT.LE.0.) RETURN C C Plotting the signal: NUM1= MOD(NUMS-1,14)+1 CALL PLOTN(FILPS,(NUMS-1)/14) IF(NUM1.EQ.1.AND.NUMS.NE.1) CALL PLOT(0.,0.,999) IF(NUM1.EQ.1) CALL PLOPN CALL PLOT(0.,-2.,-3) ccc A = -0.8*FLOAT(NUM/10)-1.428 A = -0.8*AINT(ALOG10(FLOAT(NUM)+.5)+1.)-1.428 CALL NUMBER(A,-0.4,0.8,FLOAT(NUM),0.,-1) CALL PLTIM(T3,T4,T3,-.3) CALL PLTIM(T3,T4,T1,-.5) CALL PLTIM(T3,T4,T2,-.5) CALL PLTIM(T3,T4,T4,-.3) CALL NUMBER(11.016,-0.200,0.4,AMP,0.,6) CALL PLOT(10.23,0.00,3) CALL PLOT( 0.00,0.00,2) A = 10.23/FLOAT(MPTS-1) X = 0. DO 11 I=N3,N4 Y = S(1,I)/AMP CALL PLOT(X,Y,2) 11 X = X+A RETURN END C C======================================================================= C C C SUBROUTINE PLOPN C EXTERNAL PLOTS,PLOT C CALL PLOTS(0,0,0) * CALL PLOT( 0. , 0. ,3) * CALL PLOT(21.0, 0. ,2) * CALL PLOT(21.0,29.7,2) * CALL PLOT( 0. ,29.7,2) * CALL PLOT( 0. , 0. ,2) CALL PLOT(5.38,29.7,-3) RETURN END C C======================================================================= C C C SUBROUTINE PLTIM(T3,T4,T,B) C EXTERNAL PLOT,NUMBER C A = (T-T3)/(T4-T3) IF(A.LT.-0.01) RETURN IF(A.GT. 1.01) RETURN A = 10.23*A CALL PLOT(A, 0.2,3) CALL PLOT(A,-0.2,2) A = A-0.457 CALL NUMBER(A,B-0.1,0.2,T,0.,2) RETURN END C C======================================================================= C C C C Fast Fourier transform of N = 2**K complex data points C SUBROUTINE FCOOLR(K,D,SN) C DIMENSION INU(15),D(*) C LX=2**K Q1=LX IL=LX SH=SN*6.28318530718/Q1 DO 10 I=1,K IL=IL/2 10 INU(I)=IL NKK=1 DO 40 LA=1,K NCK=NKK NKK=NKK+NKK LCK=LX/NCK L2K=LCK+LCK NW=0 DO 40 ICK=1,NCK FNW=NW AA=SH*FNW W1=COS(AA) W2=SIN(AA) LS=L2K*(ICK-1) DO 20 I=2,LCK,2 J1=I+LS J=J1-1 JH=J+LCK JH1=JH+1 Q1=D(JH)*W1-D(JH1)*W2 Q2=D(JH)*W2+D(JH1)*W1 D(JH)=D(J)-Q1 D(JH1)=D(J1)-Q2 D(J)=D(J)+Q1 20 D(J1)=D(J1)+Q2 DO 29 I=2,K ID=INU(I) IL=ID+ID IF(NW-ID-IL*(NW/IL)) 40,30,30 30 NW=NW-ID 29 CONTINUE 40 NW=NW+ID NW=0 DO 6 J=1,LX IF(NW-J) 8,7,7 7 JJ=NW+NW+1 J1=JJ+1 JH1=J+J JH=JH1-1 Q1=D(JJ) D(JJ)=D(JH) D(JH)=Q1 Q1=D(J1) D(J1)=D(JH1) D(JH1)=Q1 8 DO 9 I=1,K ID=INU(I) IL=ID+ID IF(NW-ID-IL*(NW/IL)) 6,5,5 5 NW=NW-ID 9 CONTINUE 6 NW=NW+ID 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