C Program 'ss.for' (Synthetic Seismograms) to read or generate and C filter the source time function. It may store the filtered source C time function and its Hilbert transform in the GSE data exchange C format, or read the frequency-damain response function and generate C synthetic seismograms 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 SS... Main program. C SIGNAL..Subroutine to generate Gabor's or Mueller's signal of C given parameters. 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 PLOPN...Simple subroutine to initiate plotting. C PLTIM...Simple subroutine to supplement the signal plots with the C time lables. C FCOOLR..Subroutine for the fast Fourier transform of N=2**K C complex data points. C C Other external procedures required: C WGSE1,WGSE2,WGSE3... Subroutines of the Fortran 77 file 'gse.for' C (package MODEL), designed to write 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) 'SSDAT','RF','SSGSE','SSLOG',KCOMP,/ C 'SSDAT'... String with the name of the input data file containing C the data describing the parameters of the source time C function and frequency-domain filter. C 'RF'... String with the name of the input data file with the C frequency-domain response functions at individual C receivers. Must be specified for positive KCOMP in C data 'SSDAT', but is not taken into account if KCOMP=0. C 'SSGSE'... String with the name of the output data file in the GSE C data exchange format, containing the seismograms at the C receivers (for KCOMP positive), or the filtered source C time function and its Hilbert transform (if KCOMP=0). C 'SSLOG'... String with the name of the output log file. C KCOMP...KCOMP=0: Output file 'SSGSE' will contain the filtered C source time function and its Hilbert transform, no C synthetic seismograms are generated. C KCOMP=1: Output file 'SSGSE' will contain the X1 component of C the synthetic seismograms. C KCOMP=2: Output file 'SSGSE' will contain the X2 component of C the synthetic seismograms. C KCOMP=3: Output file 'SSGSE' will contain the X3 component of C the synthetic seismograms. C Default: 'SSDAT'='ss.dat', 'RF'='rf.out', 'SSGSE'='ss.gse', C 'SSLOG'='sslog.out', KCOMP=0. C C Input file 'SSDAT' with the data describing the source time function: C (1) KSGNL,KEXP,MPTS,NPTS,DER,HILB,/ C KSGNL...Type of the source time function: C KSGNL=0: Digitized time function specified point by point. C KSGNL=1: Gabor signal. C KSGNL=2: Mueller signal. C KEXP... Exponent controling the frequency-domain filter, see (2). C MPTS... Number of points of the time functions at the output check C plot. The corresponding MPTS-1 time intervals are all C together scaled to 10.23cm. The signal is approximately C centred. MPTS does not apply to the spectra at the output C check plot, NPTS/2 points of each spectrum is plotted. C MPTS=0: No output check plot is generated. C NPTS... Number of the time samples for the fast Fourier transform. C NPTS is rounded up to the nearest power of 2. C DER,HILB... The source time function is DER-th derivative and C HILB-th Hilbert transform of the given signal. C Default: KEXP=1, MPTS=0, NPTS=MSS, DER=0, HILB=0, C where MSS is the array dimension decladred in the code. C (2) FMIN,FL,FR,FMAX,TL,TD,TRED,VRED,/ C FMIN,FL,FR,FMAX... Parameters of the frequency filter to be C applied to the source time function. The filter is zero C for frequencies F smaller than FMIN or greater than FMAX. C The filter is 1 between FL and FR. C Between FMIN and FL, cosine tappering C ( 0.5-0.5*cos(pi*(F-FMIN)/(FL-FMIN) )**KEXP C is used. C Between FR and FMAX, cosine tappering C ( 0.5-0.5*cos(pi*(F-FMAX)/(FR-FMAX) )**KEXP C is used. C TL... Reference time of the given signal. C TD... Time interval to digitize the source time function and C seismograms. C TRED,VRED... Specification of the time window for synthetic C seismograms. C VRED.EQ.0: Seismogram is centred in the time interval of C length (NPTS-1)*TD according to the travel times of the C first and last arrivals at the receiver. C VRED.NE.0: Time of the first sample of the time window is C T=TRED+R/VRED, where R is the hypocentral distance. C Default: FMIN=0, FL=0, FR=FMAX, FMAX=0.5/TD, TL=0, TRED=0, VRED=0. C (3) PAR1,PAR2,PAR3,PAR4,...,/ C Parameters of the given signal. C KSGNL=0: Digitized time function specified point by point: C PAR1,PAR2,PAR3,PAR4,PAR5,PAR6,... ...time function C digitized with step TD. PAR1 corresponds to time TL. C KSGNL=1: Gabor signal: C S(t)=PAR4*exp(-(2*pi*PAR1*(t-TL)/PAR2)**2) C *cos(2*pi*PAR1*(t-TL)+PAR3) C PAR1... Prevailing frequency. C PAR2... Relative width of the signal. C PAR3... Phase. C PAR4... Amplitude of the envelope. C KSGNL=2: Mueller signal: C For 0.LE.(t-TL).LE.PAR2/(2*PAR1): C S(t)=A*(sin(2*pi*PAR1*(t-TL)) C -sin(2*pi*PAR1*(t-TL)*(PAR2+2)/PAR2) C * PAR2 / (PAR2+2) C -(1-cos(2*pi*PAR1*(t-TL)/PAR2)) C * sin(pi*PAR2) / (PAR2+2) ) C Otherwise: C S(t)=0 C PAR1... Reference frequency. C PAR2... Relative width of the signal. If PAR2 is integer, C the number of local extrema of the signal. C PAR3... No meaning. C PAR4... Maximum amplitude of the signal. Determines A. 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) X1,X2,X3,TMIN,TMAX,AMAX,/ 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 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 C Output log file 'SSLOG': C Contains information on the input data, source time function, C synthetic seismograms. Often may be discarded. C C Output check plot: C 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 If KCOMP.GT.0, contains plots of the synthetic seismograms at the C 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 spektra. C NPTS points are plotted for the synthetic seismograms. C C======================================================================= C C Filenames: CHARACTER*80 FSSDAT,FILERF,FSSGSE,FSSLOG C INTEGER MSS PARAMETER (MSS=2048) PARAMETER (UNDEF=-999999.) CHARACTER*80 TEXT1 REAL S1(2,MSS),S2(2,MSS),PSGNL(10),SS(MSS) 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 C....................................................................... C C Opening data files: FSSDAT='ss.dat' FILERF='rf.out' FSSGSE='ss.gse' FSSLOG='sslog.out' KCOMP=0 WRITE(*,'(A)') * ' Enter 4 filenames (ss.dat, rf.out, ss.gse, sslog.out): ' READ(*,*) FSSDAT,FILERF,FSSGSE,FSSLOG,KCOMP OPEN(5,FILE=FSSDAT,STATUS='OLD') OPEN(6,FILE=FSSLOG) C Input data for the source signal WRITE(6,'(A)') ' Input data for source signal' KEXP=1 MPTS=0 NPTS=MSS DER=0. HILB=0. READ (5,*) KSGNL,KEXP,MPTS,NPTS,DER,HILB WRITE(6,'(A,5I4,2F8.3)') * ' ***',KSGNL,KEXP,KCOMP,MPTS,NPTS,DER,HILB FMIN=0. FL=0. FR=UNDEF FMAX=UNDEF TL=0. TRED=0. VRED=0. READ (5,*) FMIN,FL,FR,FMAX,TL,TD,TRED,VRED IF(FMAX.EQ.UNDEF) THEN FMAX=0.5/TD END IF IF(FR.EQ.UNDEF) THEN FR=FMAX END IF WRITE(6,'(A,10F8.3)') ' ***',FMIN,FL,FR,FMAX,TL,TD,TRED,VRED 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 PAUSE 'Error in SS: Array dimension MSS less than NPTS.' STOP END IF IF (KSGNL.LE.0) THEN DO 12 I=1,NPTS S1(1,I)=0. 12 CONTINUE READ (5,*) (S1(1,I),I=1,NPTS) 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 TL= TL-TD*FLOAT(N1) ELSE READ (5,*) (PSGNL(I),I=1,10) WRITE(6,'(A,10F8.3)') ' ***',(PSGNL(I),I=1,10) CALL SIGNAL(KSGNL,NPTS,TL,TD,S1,PSGNL) END IF DO 20 I=1,NPTS S1(2,I)= 0. 20 CONTINUE C IF(KCOMP.GT.0) THEN OPEN(4,FILE=FILERF,STATUS='OLD') END IF OPEN(7,FILE=FSSGSE) C C Plotting the input signal: WRITE(6,'(/A,I2)') ' Source signal No.',KSGNL WRITE(6,'(2A/2A)') ' * Left-hand Left-hand Right-han', * 'd Right-hand Non-zero Maximum ', * ' tip hill-side hill-sid', * 'e tip range amplitude' CALL PLSIG(MPTS,1,1,MPTS,NPTS,TL,TD,S1,A,I,J) C C Spectrum of the input signal: CALL FCOOLR(KPTS,S1,1.) FD= 1./TD/FLOAT(NPTS) C C Ampl. 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) CALL PLSIG(MPTS,3,3,NPTS/2,NPTS/2,0.,FD,S2(2,1),A,I,J) 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) 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,TL,TD,S2,A,N1,N2) CALL PLSIG(MPTS,6,6,MPTS,NPTS,TL,TD,S2(2,1),A,N,J) C Legenda: 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) IF(KCOMP.LE.0) THEN C Writing the source time function and its Hilbert transform: CALL WGSE1(7,' ') DO 51 I=1,NPTS SS(I)=S2(1,I) 51 CONTINUE CALL WGSE2(7,' ',' ',0,0.,0.,0.,TL,TD,NPTS,SS,0,TCOM,VCOM) DO 52 I=1,NPTS SS(I)=S2(2,I) 52 CONTINUE CALL WGSE2(7,' ',' ',0,0.,0.,0.,TL,TD,NPTS,SS,0,TCOM,VCOM) CALL WGSE3(7) WRITE(*,'(A)') '+SS: Source time function generated.' STOP END IF IF(MPTS.GT.0) CALL PLOT(0.,0.,999) IF (N1.GE.NPTS.OR.N.GE.NPTS) THEN PAUSE * 'Error in SS: Too small number NPTS of time samples for FFT.' STOP END IF N1= MIN0(N1,I) N2= MIN0(N2,J) IF (KCOMP.GT.3) THEN PAUSE 'Error in SS: KOMP greater than 3.' STOP END IF C C....................................................................... C C Headlines of files: WRITE(6,'(/A)') * ' Beginning of the input file with frequency response' READ (4,*) TEXT1 WRITE(6,'(2A)') ' ***',TEXT1 READ (4,*) (VCOM(I),I=1,3) WRITE(6,'(A,10F8.3)') ' ***',(VCOM(I),I=1,3) READ (4,*) FMINIM,A,NF WRITE(6,'(A,2E12.5,I4)') ' ***',FMINIM,A,NF IF (NPTS.NE.INT(1./A/TD+.5)) THEN PAUSE 'Error in SS: Inconsistent time and frequency steps.' STOP END IF MINIM= INT(FMINIM/FD+1.5) MAXIM= MINIM+NF-1 IF (NF1+1.LT.MINIM) THEN PAUSE * 'Error in SS: Missing low frequencies in response function.' STOP END IF IF (MAXIM+NF2.LT.NPTS/2) THEN PAUSE * 'Error in SS: Missing high frequencies in response function.' STOP END IF CALL WGSE1(7,TEXT1) WRITE(6,'(/A)') ' Synthetic sections at receivers' WRITE(6,'(2A/2A)') ' * Coordinates of the receiver ', * ' First Last Upper ', * ' X Y Z ', * ' arrival arrival amplitude' WRITE(6,'(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 C Loop for the receivers: NUMS= 1 DO 79 IREC=1,999999 X=UNDEF TMIN= 999999. TMAX=-999999. AMAX=0. READ (4,*,END=90) X,Y,Z,TMIN,TMAX,AMAX IF(X.EQ.UNDEF) THEN GO TO 90 END IF WRITE(6,'(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. 58 CONTINUE N = MIN0(NPTS,MAXIM+1) DO 59 I=N,NPTS S2(1,I)= 0. S2(2,I)= 0. 59 CONTINUE C IF(KCOMP.EQ.1) THEN READ (4,'(12F6.0)') * (S2(1,I),S2(2,I),AUX,AUX,AUX,AUX,I=MINIM,MAXIM) ELSE IF(KCOMP.EQ.2) THEN READ (4,'(12F6.0)') * (AUX,AUX,S2(1,I),S2(2,I),AUX,AUX,I=MINIM,MAXIM) ELSE IF(KCOMP.EQ.3) THEN READ (4,'(12F6.0)') * (AUX,AUX,AUX,AUX,S2(1,I),S2(2,I),I=MINIM,MAXIM) END IF 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 65 CONTINUE CALL FCOOLR(KPTS,S2,-1.) N = INT((TMAX+TMIN)/(TD+TD)) IF(VRED.GT.0) N=NPTS/2+ * INT((TRED+SQRT((VCOM(1)-X)**2+(VCOM(2)-Y)**2 * +(VCOM(3)-Z)**2)/VRED)/TD) TMIN= TL+TD*FLOAT(N) N = MOD(N,NPTS) IF(N.NE.0) THEN DO 66 I=1,N J = NPTS-N+I S2(2,J)= S2(1,I) 66 CONTINUE K = NPTS-N END IF DO 68 I=1,K J = N+I S2(2,I)= S2(1,J) 68 CONTINUE CALL PLSIG * (MPTS+1,NUMS,IREC,NPTS,NPTS,TMIN,TD,S2(2,1),AMAX,N1,N2) IF(N1.NE.NPTS) THEN NUMS= NUMS+1 TMIN= TMIN+TD*FLOAT(N1) N2= NPTS-N2 N = N2-N1 DO 69 I=1,N N1= N1+1 SS(I)=S2(2,N1) 69 CONTINUE CALL WGSE2 * (7,' ',' ',KCOMP,X,Y,Z,TMIN,TD,N,SS,NCOM,TCOM,VCOM) ELSE CALL WGSE2(7,' ',' ',KCOMP,X,Y,Z,0.,TD,0,SS,NCOM,TCOM,VCOM) END IF ELSE CALL WGSE2(7,' ',' ',KCOMP,X,Y,Z,0.,TD,0,SS,NCOM,TCOM,VCOM) END IF 79 CONTINUE C C End of computation: 90 IF(MPTS.GT.-1) CALL PLOT(0.,0.,999) CALL WGSE3(7) WRITE(*,'(A)') '+SS: done. ' STOP END C C======================================================================= C SUBROUTINE SIGNAL(KSGNL,NPTS,TL,TD,S,PAR) C REAL S(2,NPTS),PAR(*) C GO TO (10,20),KSGNL C C Gabor's signal 10 T = -TD*FLOAT(NPTS/2) TL= TL+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) 11 T = T+TD RETURN C C Mueller's signal: 20 N2= IFIX(PAR(2)/PAR(1)/TD/2.)+1 N1= (NPTS-N2)/2 TL= TL-TD*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 21 I=1,N1 21 S(1,I)= 0. T = 0. F = 0. N2= N1+N2 N1= N1+1 DO 22 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))) 22 T = T+TD IF(PAR(4).NE.0.) F=F/PAR(4) DO 23 I=N1,N2 23 S(1,I)= S(1,I)/F N2= N2+1 DO 24 I=N2,NPTS 24 S(1,I)= 0. RETURN END C C======================================================================= C SUBROUTINE PLSIG(KPLOT,NUMS,NUM,MPTS,NPTS,TL,TD,S,AMP,N1,N2) C REAL S(2,NPTS) C IF(MPTS.GT.NPTS) THEN WRITE(*,'(2(A,I6))') ' MPTS=',MPTS,' NPTS=',NPTS PAUSE 'Error in SS: MPTS greater than NPTS' STOP 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+TD*FLOAT(N1) T2= TL+TD*FLOAT(NPTS-N2-1) T3= TL+TD*FLOAT(N3-1) T4= TL+TD*FLOAT(N4-1) A = T2-T1 WRITE(6,'(I4,5F11.3,E11.3)') NUM,T3,T1,T2,T4,A,AMP IF(KPLOT.LE.0) RETURN C C Plotting the signal: NUM1= MOD(NUMS-1,14)+1 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 SUBROUTINE PLOPN 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 SUBROUTINE PLTIM(T3,T4,T,B) 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 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 'gse.for' INCLUDE 'calcops.for' C C======================================================================= C