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.50 C Date: 2000, November 26 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 Description of data files: C C Input data read from the standard input device (*): C The data are read by the list directed input (free format) and C consist of a single string 'SEP': C 'SEP'...String in apostrophes containing the name of the input C SEP parameter or history file with the input data. C No default, 'SEP' must be specified and cannot be blank. C C C Input data file 'SEP': C File 'SEP' has the form of the SEP C parameter file. The parameters, which do not differ from their C defaults, need not be specified in file 'SEP'. 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 (NFFT-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=real... 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 SMALL=positive real... Amplitudes in absolute value not exceeding C SMALL times the maximum absolute amplitude of the C seismogram are assumed zero. C Default: SMALL=0.002 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) X1SRC,X2SRC,X3SRC,/ C X1SRC,X2SRC,X3SRC... 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 Endpoints of the time (or frequency) interval plotted and the C leftmost and rightmost samples of amplitudes in absolute C value grater than SMALL times the maximum absolute C amplitude of the seismogram are labelled with the C corresponding time (or frequency) values. C MPTS points are plotted for input signal, filtered signal and its C Hilbert transform. C NFFT/2 points are plotted for the spectra. C NFFT points are plotted for the synthetic seismograms. 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 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*5 TCOM(NCOM) REAL VCOM(NCOM),KEXP DATA TCOM/'X1SRC','X2SRC','X3SRC'/ C C----------------------------------------------------------------------- C C Reading name of SEP file with input data: WRITE(*,'(A)') '+SS: Enter input filename: ' FILDAT=' ' READ(*,*) FILDAT C C Reading all data from the SEP file into the memory: IF (FILDAT.NE.' ') THEN CALL RSEP1(LU5,FILDAT) ELSE C SS-12 CALL ERROR('SS-12: SEP file not given') C Input file in the form of the SEP (Stanford Exploration Project) C parameter or history file must be specified. C There is no default filename. ENDIF C WRITE(*,'(A)') '+SS: Working... ' C 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('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('FMAX ' ,FMAX ,0.5/DT) CALL RSEP3R('FHIGH' ,FR ,FMAX) CALL RSEP3R('KEXP ' ,KEXP ,1.0) 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 C SS-05 CALL ERROR * ('SS-05: Too large number NFFT of time samples for FFT') C NFFT should not exceed 2**15. Check the input data. C To use greater NFFT, edit subroutine FCOOLR at the end of this C file. 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 NFFT 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 IF(FILPS.EQ.' ') THEN WRITE(*,'(A)') '+SS: Nothing to do. ' ELSE WRITE(*,'(A)') '+SS: Signal plotted. ' END IF 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' TEXT1=' ' 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))) CALL RSEP3R('SMALL',SMALL,0.002) EPS= SMALL*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