C
C PROGRAM SYNT - SYNTHETIC SEISMOGRAM COMPUTATION C C ************************************************************** C EXTERNAL SIGNAL C CHARACTER*80 FILEIN,FILEOU,FILE1,FILE2 CHARACTER*80 TITLE,MPRINT,IPRINT,PSTEXT DIMENSION NDIS(100) COMPLEX S(2048),REF(2048),SS(2048) COMMON/SPECTR/S,REF,DF,NI,NT,NT2,NTM,ILS,IRS,NFS,NFILT,NWIN COMMON/SIG/F(2048),TO,DT,NPTS,NO,NSIG,NPL1,NPL2,NPL3, 1NPL4,IPR1,NDER,NINT COMMON/WIND/FLO,FLEFT,FRO,FRIGHT,FEXP,ILF,IRF,NFW COMMON/AUXIL/FF(2048),GG(2048),IS(2048),AUX(4),IAUX(4) COMMON/PLOT1/ARESP(6),ASPECT(6),ASYNT(6),JRESP(4),JSPECT(4),JSYNT( 14),NOO,RMIN,RMAX,GLR,TRMIN,TRMAX,GLTR,IR,ITR,NRED,NSH,AREDUC,TITLE 2,NPLOT,MCOMP C C C************************************************** C LIN=5 LOU=6 LU7=1 LU8=2 FILEIN='synfan.dat' FILEOU='synfan.out' FILE1='lu7.in' FILE2='lu8.out' WRITE(*,'(2A)') ' SPECIFY NAMES OF INPUT AND OUTPUT FILES', 1' LIN, LOU, LU7, LU8: ' READ(*,*) FILEIN,FILEOU,FILE1,FILE2 IF(FILE1.EQ.' ') GO TO 99 IF(FILE2.EQ.' ') LU8=0 OPEN(LIN,FILE=FILEIN,FORM='FORMATTED',STATUS='OLD') OPEN(LOU,FILE=FILEOU,FORM='FORMATTED') OPEN(LU7,FILE=FILE1,FORM='FORMATTED',STATUS='OLD') IF(LU8.NE.0)OPEN(LU8,FILE=FILE2,FORM='FORMATTED') C C************************************************** C C TITLE='SYNFAN' PSTEXT=' ' IPR1=0 IPR2=0 READ(LIN,*)TITLE READ(LIN,*)IPR1,IPR2,PSTEXT WRITE(LOU,111)TITLE WRITE(LOU,103)IPR1,IPR2,PSTEXT REWIND LU7 IF(LU8.NE.0)REWIND LU8 C C READ INPUT PARAMETERS C NSIG=0 NPTS=400 NT=1024 NWIN=0 NFILT=0 NPLOT=0 NNPLOT=0 NSTOP=0 NSH=7 NDER=0 NINT=0 READ(LIN,*)NSIG,NPTS,NT,NWIN,NFILT,NPLOT,NSTOP,NSH,NDER,NINT WRITE(LOU,101)NSIG,NPTS,NT,NWIN,NFILT,NPLOT,NSTOP,NSH, 1NDER,NINT IF(NPLOT.NE.0)THEN IF(IPR1.LT.0)THEN CALL PLOTN(PSTEXT,0) IPR1=-IPR1 END IF CALL PLOTS(LDUM1,LDUM2,7) READ(LIN,*)NNPLOT,(NDIS(I),I=1,NNPLOT) WRITE(LOU,101)NNPLOT,(NDIS(I),I=1,NNPLOT) END IF NT2=NT/2 NTM=NT2+1 READ(LU7,110)MPRINT READ(LU7,110)IPRINT READ(LU7,152)XSOUR,YSOUR,ZSOUR,TSOUR,RSTEP,FL,DF READ(LU7,102)NDST,NFS,ILS,MCOMP,ILOC DT=1./(DF*FLOAT(NT)) WRITE(LOU,111)MPRINT WRITE(LOU,111)IPRINT WRITE(LOU,152)XSOUR,YSOUR,ZSOUR,TSOUR,RSTEP,FL,DF WRITE(LOU,101)NDST,NFS,ILS,MCOMP WRITE(LU8,110)MPRINT WRITE(LU8,110)IPRINT WRITE(LU8,110)TITLE WRITE(LU8,152)XSOUR,YSOUR,ZSOUR,TSOUR,RSTEP,DT,DF WRITE(LU8,101)NDST,NT,MCOMP,ILOC IRS=ILS+NFS IRSS=IRS+1 NI=10 IF(NT.EQ.2048)NI=11 IF(NT.EQ.512)NI=9 IF(NT.EQ.256)NI=8 IF(NT.EQ.128)NI=7 READ(LIN,*)ARESP,JRESP READ(LIN,*)ASPECT,JSPECT READ(LIN,*)ASYNT,JSYNT WRITE(LOU,106)ARESP,JRESP WRITE(LOU,106)ASPECT,JSPECT WRITE(LOU,106)ASYNT,JSYNT IF(NPLOT.EQ.0)GO TO 7 CALL PLOT(3.,5.,-3) C C INPUT SIGNAL C 7 CALL SIGNAL(LIN,LOU) C C SPECTRUM OF THE INPUT SIGNAL C CALL FT(LIN,LOU) IF(NSTOP.EQ.1.AND.NPLOT.NE.0)CALL PLOT(0.,0.,999) IF(NSTOP.EQ.1)STOP C C LOOP FOR RECEIVERS C DO 40 IDIST=1,NDST NPR1=0 NPR2=0 NPR3=0 READ(LU7,160)DST,AA WRITE(LOU,160)DST,AA NFF=2*NFS DO 49 I=1,NTM FF(I)=0. SS(I)=(0.,0.) 49 REF(I)=(0.,0.) READ(LU7,161)(IS(I),I=1,NFF) DO 41 I=1,NFS A=AA*FLOAT(IS(2*(I-1)+1))/99999.1 B=AA*FLOAT(IS(2*I))/99999.1 FF(I+ILS)=SQRT(A*A+B*B) 41 REF(I+ILS)=CMPLX(A,B) IF(NNPLOT.EQ.0)GO TO 42 DO 43 I=1,NNPLOT IF(IDIST.EQ.NDIS(I))GO TO 45 43 CONTINUE GO TO 42 45 CONTINUE READ(LIN,*)NPR1,NPR2,NPR3,TMIN,TMAX,TLEN WRITE(LOU,170)NPR1,NPR2,NPR3,TMIN,TMAX,TLEN ASYNT(1)=TMIN ASYNT(2)=TMAX ASYNT(3)=TLEN CALL PLOT(3.,0.,-3) CALL SYMBOL(-0.25*TLEN,0.,0.05*TLEN,'X(KM)=',90.,6) CALL NUMBER(-0.25*TLEN,0.3*TLEN,0.05*TLEN,DST,90.,2) 42 CONTINUE C C PRINTING AND PLOTTING OF FREQUENCY RESPONSE C IF(IPR2.LE.1.AND.NPR1.EQ.0)GO TO 50 CALL REDUC(AREDUC,1,IRSS,FF,GG,IS) IF(IPR2.LE.1)GO TO 48 WRITE(LOU,114)DST,DF,AREDUC,IRSS WRITE(LOU,107)(IS(I),I=1,IRSS) 48 IF(NPR1.EQ.0)GO TO 50 NOO=1 NRED=NPR1 IF(NPR1.EQ.1)CALL ZPLOT(0.,DF,ILS,IRSS,GG) IF(NPR1.EQ.2)CALL ZPLOT(0.,DF,ILS,IRSS,FF) 50 CONTINUE C C AMPLITUDE SPECTRUM OF SYNTHETIC SEISMOGRAM C DO 18 I=1,NTM 18 FF(I)=0. DO 11 I=ILF,IRF SS(I)=REF(I)*S(I) A=REAL(SS(I)) B=AIMAG(SS(I)) 11 FF(I)=SQRT(A*A+B*B) DO 12 I=2,NT2 IJ=NT-I+2 12 SS(IJ)=CONJG(SS(I)) IF(IPR2.LE.2.AND.NPR2.EQ.0)GO TO 60 CALL REDUC(AREDUC,1,IRSS,FF,GG,IS) IF(IPR2.LE.2)GO TO 52 WRITE(LOU,115)DST,DF,AREDUC,IRSS WRITE(LOU,107)(IS(I),I=1,IRSS) 52 IF(NPR2.EQ.0)GO TO 60 NOO=5 NRED=NPR2 IF(NPR2.EQ.1)CALL ZPLOT(0.,DF,ILS,IRSS,GG) IF(NPR2.EQ.2)CALL ZPLOT(0.,DF,ILS,IRSS,FF) 60 CONTINUE C C GOING BACK TO TIME DOMAIN. SYNTHETIC SEISMOGRAMS C CALL FCOOLR(NI,SS,-1.0) DEL=1.0/FLOAT(NT) DO 20 I=1,NT 20 FF(I)=DEL*REAL(SS(I)) CALL REDUC(AREDUC,1,NT,FF,GG,IS) WRITE(LU8,162)DST,TO,AREDUC,NT WRITE(LU8,107)(IS(I),I=1,NT) C IF(IPR2.EQ.0.AND.NPR3.EQ.0)GO TO 70 IF(IPR2.EQ.0)GO TO 65 WRITE(LOU,116)DST,TO,DT,AREDUC,NT WRITE(LOU,108)(IS(I),I=1,NT) 65 IF(NPR3.EQ.0)GO TO 70 NOO=6 NRED=NPR3 IF(AREDUC.EQ.0.)GO TO 66 TSTART=TMIN-TO IPOM=IFIX(TSTART/DT)+1 67 IF(IPOM.GT.0)GO TO 61 IPOM=IPOM+NT GO TO 67 61 IF(IPOM.LE.NT)GO TO 62 IPOM=IPOM-NT GO TO 61 62 CONTINUE IIAUX=NT-IPOM DO 63 J=1,NT IF(J.GE.IPOM)JAUX=J-IPOM+1 IF(J.LT.IPOM)JAUX=J+IIAUX+1 63 FF(JAUX)=GG(J)*AREDUC TLEN=ASYNT(2)-ASYNT(1) NFIN=IFIX(TLEN/DT)+1 CALL REDUC(AREDUC,1,NFIN,FF,GG,IS) 66 CONTINUE IF(NPR3.EQ.1)CALL ZPLOT(TMIN,DT,1,NFIN,GG) IF(NPR3.EQ.2)CALL ZPLOT(TMIN,DT,1,NFIN,FF) 70 CONTINUE 40 CONTINUE C C END OF THE LOOP FOR RECEIVERS C 99 CONTINUE IF(NPLOT.NE.0)CALL PLOT(0.,0.,999) C 101 FORMAT(16I5) 102 FORMAT(26I3) 103 FORMAT(2I5,1X,A) 106 FORMAT(6F10.5,4I5) 107 FORMAT(20I4) 108 FORMAT(1X,20I4) 110 FORMAT(A) 111 FORMAT(1X,A) 114 FORMAT(/1X,'FREQUENCY RESPONSE-AMPLITUDES X(KM)=',F10.5/,2X, 1'DF=',F10.5,2X,'MAX=',E16.7,2X,'N=',I5) 115 FORMAT(/1X,'AMPLITUDE SPECTRUM OF SYNTHETIC SEISMOGRAM X(KM)=', 1F10.5/2X,'DF=',F10.5,2X,'MAX=',E16.8,2X,'N=',I5) 116 FORMAT(/1X,'SYNTHETIC SEISMOGRAM X(KM)=',F10.5/2X,'T0=', 1F10.5,2X,'DT=',F10.5,2X,'MAX=',E16.8,2X,'N=',I5) 152 FORMAT(5F10.5,2E15.7) 154 FORMAT(3E16.8) 160 FORMAT(F10.3,E12.5) 161 FORMAT(12I6) 162 FORMAT(2F10.3,1E12.5,I5) 170 FORMAT(3I5,3F10.5) REWIND LU7 REWIND LU8 STOP END C C SUBROUTINE BORD CHARACTER*80 TITLE COMMON/PLOT1/ARESP(6),ASPECT(6),ASYNT(6),JRESP(4),JSPECT(4),JSYNT( 14),NOO,RMIN,RMAX,GLR,TRMIN,TRMAX,GLTR,IR,ITR,NRED,NSH,AREDUC,TITLE 2,NPLOT,MCOMP IF(NOO.GT.1)GO TO 1 RMIN=ARESP(1) RMAX=ARESP(2) GLR=ARESP(3) TRMIN=ARESP(4) TRMAX=ARESP(5) GLTR=ARESP(6) IR=JRESP(1) ITR=JRESP(2) NDECX=JRESP(3) NDECY=JRESP(4) GO TO 3 1 IF(NOO.EQ.2.OR.NOO.GE.6)GO TO 2 RMIN=ASPECT(1) RMAX=ASPECT(2) GLR=ASPECT(3) TRMIN=ASPECT(4) TRMAX=ASPECT(5) GLTR=ASPECT(6) IR=JSPECT(1) ITR=JSPECT(2) NDECX=JSPECT(3) NDECY=JSPECT(4) GO TO 3 2 RMIN=ASYNT(1) RMAX=ASYNT(2) GLR=ASYNT(3) TRMIN=ASYNT(4) TRMAX=ASYNT(5) GLTR=ASYNT(6) IR=JSYNT(1) ITR=JSYNT(2) NDECX=JSYNT(3) NDECY=JSYNT(4) 3 CONTINUE CALL SYMBOL(-0.12*GLR,0.,.07*GLR,TITLE,90.,80) CALL PLOT(0.,0.,3) IF(TRMAX.GE.0.) 1CALL NUMBER(0.025*GLR,-0.2*GLR,0.05*GLR,TRMAX,90.,NDECY) IF(TRMAX.LT.0.) 1CALL NUMBER(0.025*GLR,-0.25*GLR,0.05*GLR,TRMAX,90.,NDECY) CALL PLOT(0.,0.,3) ZR=(RMAX-RMIN)/FLOAT(IR) ZTR=(TRMAX-TRMIN)/FLOAT(ITR) A=GLR/FLOAT(IR) B=GLTR/FLOAT(ITR) DO 100 J=1,ITR F=B*FLOAT(J) TG=TRMAX-ZTR*FLOAT(J) CALL PLOT(F,0.,2) CALL PLOT(F,0.2,2) IF(TG.GE.0.) 1CALL NUMBER(F+0.025*GLR,-0.2*GLR,0.05*GLR,TG,90.,NDECY) IF(TG.LT.0.) 1CALL NUMBER(F+0.025*GLR,-0.25*GLR,0.05*GLR,TG,90.,NDECY) CALL PLOT(F,0.,3) 100 CONTINUE CALL NUMBER(F+0.1*GLR,-.025*GLR,0.05*GLR,RMIN,90.,NDECX) IF(NOO.EQ.1) 1CALL SYMBOL(F+0.25*GLR,0.,0.05*GLR,'FREQUENCY RESPONSE',90.,18) IF(NOO.EQ.2) 1CALL SYMBOL(F+0.25*GLR,0.,0.05*GLR,'INPUT SIGNAL',90.,12) IF(NOO.EQ.3)CALL SYMBOL(F+0.25*GLR,0.,0.05*GLR, 1'AMPLITUDE SPECTRUM OF INPUT SIGNAL',90.,34) IF(NOO.EQ.4)CALL SYMBOL(F+0.25*GLR,0.,0.05*GLR, 1'AMPLITUDE SPECTRUM, WINDOW APPLIED',90.,34) IF(NOO.EQ.5) 1CALL SYMBOL(F+0.25*GLR,0.,0.05*GLR,'AMPLITUDE SPECTRUM',90.,18) IF(NOO.EQ.5)CALL SYMBOL(F+0.35*GLR,0.,0.05*GLR, 1'OF SYNTHETIC SEISMOGRAM',90.,23) IF(NOO.EQ.6)THEN CALL SYMBOL 1 (F+0.25*GLR,0.,0.05*GLR,'SYNTHETIC SEISMOGRAM',90.,20) IF(MCOMP.EQ.0) 1 CALL SYMBOL(F+0.35*GLR,0.,0.05*GLR,'VERTICAL',90.,8) IF(MCOMP.EQ.1) 1 CALL SYMBOL(F+0.35*GLR,0.,0.05*GLR,'X-COMPONENT',90.,11) IF(MCOMP.EQ.2) 1 CALL SYMBOL(F+0.35*GLR,0.,0.05*GLR,'Y-COMPONENT',90.,11) END IF IF(NOO.EQ.7)CALL SYMBOL(F+0.25*GLR,0.,0.05*GLR, 1'INPUT SIGNAL AFTER WINDOWING',90.,28) IF(NOO.EQ.7)NOO=2 CALL PLOT(F,0.,3) DO 200 J=1,IR C=A*FLOAT(J) CALL PLOT(F,C,2) CALL PLOT(F-0.2,C,2) ZG=RMIN+ZR*FLOAT(J) CALL NUMBER(F+0.1*GLR,C-.08*GLR,0.05*GLR,ZG,90.,NDECX) CALL PLOT(F,C,3) 200 CONTINUE IF(NOO.NE.2.AND.NOO.NE.6) 1CALL SYMBOL(F+0.18*GLR,C-0.25*GLR,0.04*GLR,'F(HZ)',90.,5) IF(NOO.EQ.2.OR.NOO.EQ.6) 1CALL SYMBOL(F+0.18*GLR,C-0.3*GLR,0.04*GLR,'TIME(S)',90.,7) CALL PLOT(F,C,3) DO 300 J=1,ITR D=B*FLOAT(ITR-J) CALL PLOT(D,C,2) CALL PLOT(D,C-0.2,2) 300 CALL PLOT(D,C,2) IF(NRED.EQ.1) 1CALL SYMBOL(D-0.05*GLR,0.,0.03*GLR,'REDUCTION FACTOR= ',90.,18) IF(NRED.EQ.1) 1CALL NUMBER(D-0.05*GLR,0.7*GLR,0.03*GLR,AREDUC,90.,5) CALL PLOT(D,C,3) DO 400 J=1,IR E=A*FLOAT(IR-J) CALL PLOT(D,E,2) CALL PLOT(D+0.2,E,2) 400 CALL PLOT(D,E,2) RETURN END C C SUBROUTINE ZPLOT(FMIN,DF,NMIN,NMAX,F) CHARACTER*80 TITLE DIMENSION F(1) COMMON/PLOT1/ARESP(6),ASPECT(6),ASYNT(6),JRESP(4),JSPECT(4),JSYNT( 14),NOO,RMIN,RMAX,GLR,TRMIN,TRMAX,GLTR,IR,ITR,NRED,NSH,AREDUC,TITLE 2,NPLOT,MCOMP CALL BORD GMR=GLR/(RMAX-RMIN) GMTR=GLTR/(TRMAX-TRMIN) C DO 1 II=NMIN,NMAX VO=GMR*(FMIN-RMIN+DF*FLOAT(II-1)) IF((VO+0.0001).GE.0.)GO TO 2 1 CONTINUE GO TO 99 C 2 IMAX1=0 WO=GMTR*(TRMAX-F(II)) VOS=VO WOS=WO IF(WO.GT.0.)GO TO 3 IMAX1=1 CALL PLOT(0.,VO,3) GO TO 4 3 CALL PLOT(WO,VO,3) 4 CONTINUE II=II+1 C DO 10 I=II,NMAX VO=GMR*(FMIN-RMIN+DF*FLOAT(I-1)) WO=GMTR*(TRMAX-F(I)) IF(VO.GT.(GLR+0.0001))GO TO 11 IF(WO.LT.0.)GO TO 12 IF(IMAX1.EQ.0)GO TO 13 VV=VOS-WOS*(VO-VOS)/(WO-WOS) CALL PLOT(0.,VV,2) CALL PLOT(WO,VO,2) IMAX1=0 GO TO 9 13 CALL PLOT(WO,VO,2) IMAX1=0 GO TO 9 12 IF(IMAX1.EQ.0)GO TO 15 CALL PLOT(0.,VO,2) IMAX1=1 GO TO 9 15 VV=VOS-WOS*(VO-VOS)/(WO-WOS) CALL PLOT(0.,VV,2) CALL PLOT(0.,VO,2) IMAX1=1 GO TO 9 11 IF(IMAX1.EQ.0)GO TO 16 CALL PLOT(0.,RMAX,2) GO TO 99 16 IF(WO.LT.0.)GO TO 18 17 WW=WOS+(WO-WOS)*(GLR-VOS)/(VO-VOS) CALL PLOT(WW,GLR,2) GO TO 99 18 VV=VOS-WOS*(VO-VOS)/(WO-WOS) IF(VV.GT.GLR)GO TO 17 CALL PLOT(0.,VV,2) CALL PLOT(0.,GLR,2) GO TO 99 9 VOS=VO WOS=WO 10 CONTINUE CALL PLOT(GMTR*TRMAX,VOS,2) CALL PLOT(GMTR*TRMAX,GLR,2) C 99 CONTINUE CALL PLOT(GLTR+FLOAT(NSH),0.,-3) RETURN END C C SUBROUTINE SIGNAL(LIN,LOU) CHARACTER*80 TITLE COMPLEX S(2048),REF(2048) COMMON/SIG/F(2048),TO,DT,NPTS,NO,NSIG,NPL1,NPL2,NPL3, 1NPL4,IPR,NDER,NINT COMMON/SPECTR/S,REF,DF,NI,NT,NT2,NTM,ILS,IRS,NFS,NFILT,NWIN COMMON/AUXIL/FF(2048),GG(2048),IS(2048),AUX(4),IAUX(4) COMMON/PLOT1/ARESP(6),ASPECT(6),ASYNT(6),JRESP(4),JSPECT(4),JSYNT( 14),NOO,RMIN,RMAX,GLR,TRMIN,TRMAX,GLTR,IR,ITR,NRED,NSH,AREDUC,TITLE 2,NPLOT,MCOMP C NPL1=0 NPL2=0 NPL3=0 NPL4=0 IF(NPLOT.EQ.0)GO TO 3 READ(LIN,*) NPL1,NPL2,NPL3,NPL4 WRITE(LOU,100)NPL1,NPL2,NPL3,NPL4 3 READ(LIN,*)IAUX,AUX,TO WRITE(LOU,104)IAUX,AUX,TO DO 1 I=1,NPTS F(I)=0. 1 GG(I)=0. C IF(NSIG)30,20,30 30 READ(LIN,*)(IS(I),I=1,NPTS) WRITE(LOU,120)NSIG WRITE(LOU,102)(IS(I),I=1,NPTS) A=10.**(-NSIG) DO 32 I=1,NPTS F(I)=A*FLOAT(IS(I)) IS(I)=0 32 CONTINUE GO TO 60 20 CALL SFUN 60 CONTINUE C CALL REDUC(AREDUC,1,NPTS,F,GG,IS) IF(IPR.EQ.0)GO TO 61 WRITE(LOU,121)TO,DT,AREDUC,NPTS WRITE(LOU,107)(IS(I),I=1,NPTS) 61 IF(NPL1.EQ.0)GO TO 10 NOO=2 NRED=NPL1 IF(NPL1.EQ.1)CALL ZPLOT(TO,DT,1,NPTS,GG) IF(NPL1.EQ.2)CALL ZPLOT(TO,DT,1,NPTS,F) C 10 CONTINUE RETURN 100 FORMAT(16I5) 102 FORMAT(2X,16I5) 104 FORMAT(4I5,5F10.5) 107 FORMAT(20I4) 120 FORMAT(2X,'SIGNAL. NSIG=',I3) 121 FORMAT(/,1X,'INPUT SIGNAL T0=',F10.5,2X,'DT=',F10.5,2X, 1'MAX=',E16.8,2X,'N=',I5) END C C SUBROUTINE FT(LIN,LOU) CHARACTER*80 title COMPLEX S(2048),REF(2048),SAUX COMMON/SIG/F(2048),TO,DT,NPTS,NO,NSIG,NPL1,NPL2,NPL3, 1NPL4,IPR,NDER,NINT COMMON/SPECTR/S,REF,DF,NI,NT,NT2,NTM,ILS,IRS,NFS,NFILT,NWIN COMMON/AUXIL/FF(2048),GG(2048),IS(2048),AUX(4),IAUX(4) COMMON/WIND/FLO,FLEFT,FRO,FRIGHT,FEXP,ILF,IRF,NFW COMMON/PLOT1/ARESP(6),ASPECT(6),ASYNT(6),JRESP(4),JSPECT(4),JSYNT( 14),NOO,RMIN,RMAX,GLR,TRMIN,TRMAX,GLTR,IR,ITR,NRED,NSH,AREDUC,TITLE 2,NPLOT,MCOMP C IF(NWIN.EQ.0)GO TO 1 READ(LIN,*)FLO,FLEFT,FRIGHT,FRO,FEXP WRITE(LOU,101)FLO,FLEFT,FRIGHT,FRO,FEXP 1 CONTINUE C C PREPARE FOR FCOOLR DO 10 I=1,NPTS 10 S(I)=CMPLX(F(I),0.0) IF(NT.EQ.NPTS)GO TO 11 NPTS1=NPTS+1 DO 12 I=NPTS1,NT 12 S(I)=(0.0,0.0) 11 CONTINUE CALL FCOOLR(NI,S,1.0) DO 21 I=1,NT A=REAL(S(I)) B=AIMAG(S(I)) 21 FF(I)=SQRT(A*A+B*B) C C DERIVATIVE (FOR NDER.NE.0) OR INTEGRAL (FOR NINT.NE.0) C OF THE INPUT SIGNAL C IF(NDER.NE.0.OR.NINT.NE.0)THEN DO 22 I=2,NT AUX1=FLOAT(I-1)*DF IF(NDER.NE.0)SAUX=CMPLX(0.,-AUX1) IF(NINT.NE.0)SAUX=CMPLX(0.,1./AUX1) S(I)=S(I)*SAUX 22 CONTINUE END IF C C PLOTTING AMPLITUDE SPECTRUM OF THE INPUT SIGNAL C CALL REDUC(AREDUC,1,IRS,FF,GG,IS) IF(IPR.EQ.0)GO TO 25 WRITE(LOU,120)DF,AREDUC,IRS WRITE(LOU,107)(IS(I),I=1,IRS) 25 CONTINUE IF(NPL2.EQ.0)GO TO 28 NOO=3 NRED=NPL2 IF(NPL2.EQ.1)CALL ZPLOT(0.,DF,1,IRS,GG) IF(NPL2.EQ.2)CALL ZPLOT(0.,DF,1,IRS,FF) 28 CONTINUE C C WINDOWING OF AMPLITUDE SPECTRUM C ILF=ILS+1 IRF=IRS-1 NFW=IRF-ILF+1 IF(NWIN.EQ.0.AND.NFILT.EQ.0)GO TO 30 IF(NWIN.NE.0)CALL FENSTR IF(NFILT.NE.0)CALL FILTER DO 31 I=1,IRS A=REAL(S(I)) B=AIMAG(S(I)) 31 FF(I)=SQRT(A*A+B*B) CALL REDUC(AREDUC,1,IRS,FF,GG,IS) IF(IPR.EQ.0)GO TO 32 WRITE(LOU,121)DF,AREDUC,IRS WRITE(LOU,107)(IS(I),I=1,IRS) 32 CONTINUE IF(NPL3.EQ.0)GO TO 30 NOO=4 NRED=NPL3 IF(NPL3.EQ.1)CALL ZPLOT(0.,DF,1,IRS,GG) IF(NPL3.EQ.2)CALL ZPLOT(0.,DF,1,IRS,FF) 30 CONTINUE C C INPUT SIGNAL, FREQUENCY WINDOW AND FILTER APPLIED C DO 40 I=1,NT 40 REF(I)=S(I) DEL=1./FLOAT(NT) CALL FCOOLR(NI,REF,-1.0) DO 41 I=1,NT 41 FF(I)=DEL*REAL(REF(I)) CALL REDUC(AREDUC,1,NPTS,FF,GG,IS) IF(IPR.EQ.0)GO TO 42 WRITE(LOU,122)TO,DT,AREDUC,NPTS WRITE(LOU,107)(IS(I),I=1,NPTS) 42 CONTINUE IF(NPL4.EQ.0)GO TO 50 NRED=NPL4 NOO=7 IF(NPL4.EQ.1)CALL ZPLOT(TO,DT,1,NPTS,GG) IF(NPL4.EQ.2)CALL ZPLOT(TO,DT,1,NPTS,F) 50 CONTINUE RETURN C 101 FORMAT(8F10.5) 107 FORMAT(20I4) 120 FORMAT(/1X,'AMPLITUDE SPECTRUM OF THE INPUT SIGNAL'/, 12X,'DF=',F10.5,2X,'MAX=',E16.8,2X,'N=',I5) 121 FORMAT(/1X,'AMPLITUDE SPECTRUM OF INPUT SIGNAL,WINDOW APPLIED'/ 1,2X,'DF=',F10.5,2X,'MAX=',E16.8,2X,'N=',I5) 122 FORMAT(/1X,'INPUT SIGNAL, FREQUENCY WINDOW APPLIED'/2X,'T0=', 1F10.5,2X,'DT=',F10.5,2X,'MAX=',E16.8,2X,'N=',I5) END C C C C SUBROUTINE FENSTR COMPLEX S(2048),REF(2048) COMMON/SPECTR/S,REF,DF,NI,NT,NT2,NTM,ILS,IRS,NFS,NFILT,NWIN COMMON/WIND/FLO,FLEFT,FRO,FRIGHT,FEXP,ILF,IRF,NFW C ILO=FLO/DF+1.1 IRO=FRO/DF+1.1 IF(IRO.GT.IRS)IRO=IRS IF(ILO.LT.ILS)ILO=ILS IF(IRO.GT.NTM)IRO=NTM ILEFT=FLEFT/DF+1.1 IRIGHT=FRIGHT/DF+1.1 IF(IRIGHT.GT.IRO)IRIGHT=IRO IF(ILEFT.LT.ILO)ILEFT=ILO IF(IRIGHT.GT.NTM)IRIGHT=NTM NLEFT=ILEFT-ILO NRIGHT=IRO-IRIGHT IF(NLEFT.GT.0)DLEFT=3.14159/FLOAT(NLEFT) IF(NRIGHT.GT.0)DRIGHT=3.14159/FLOAT(NRIGHT) C DO 1 I=1,NTM J=NT+2-I FIF=I-IRIGHT FAF=I-ILO IF(I.GE.ILEFT.AND.I.LE.IRIGHT)GO TO 1 IF(I.LE.ILO.OR.I.GE.IRO)S(I)=(0.0,0.0) IF(I.GT.ILO.AND.I.LT.ILEFT) 1S(I)=S(I)*(0.5*(1.-COS(DLEFT*FAF)))**FEXP IF(I.GT.IRIGHT.AND.I.LT.IRO) 1S(I)=S(I)*(0.5*(COS(DRIGHT*FIF)+1.))**FEXP IF(I.EQ.1.OR.I.EQ.NTM)GO TO 1 S(J)=CONJG(S(I)) 1 CONTINUE ILF=ILO+1 IRF=IRO-1 NFW=IRF-ILF+1 RETURN END C C SUBROUTINE SFUN COMMON/AUXIL/FF(2048),GG(2048),IS(2048),AUX(4),IAUX(4) COMMON/SIG/F(2048),TO,DT,NPTS,NO,NSIG,NPL1,NPL2,NPL3, 1NPL4,IPR,NDER,NINT C IF(IAUX(1).NE.1)GO TO 10 C C GABOR SIGNAL C DO 1 I=1,NPTS T=TO-AUX(4)+DT*FLOAT(I-1) G=6.2831853*T*AUX(1) G1=G/AUX(2) G2=G1*G1 F(I)=0. IF(G2.LT.20.)F(I)=EXP(-G2)*COS(G+AUX(3)) 1 CONTINUE RETURN 10 CONTINUE C IF(IAUX(1).NE.2)GO TO 20 C C BERLAGE SIGNAL C DO 2 I=1,NPTS T=TO-AUX(4)+DT*FLOAT(I-1) G=6.2831853*T*AUX(1) G1=AUX(2)*T F(I)=0. IF(G1.LT.20..AND.T.GT.0.)F(I)=EXP(-G1)*SIN(G)*T**AUX(3) 2 CONTINUE RETURN 20 CONTINUE IF(IAUX(1).NE.3)GO TO 30 C C MULLER SIGNAL C G1=FLOAT(IAUX(2)) G2=G1/(G1+2) G3=G1+2 G1=3.14159*G1 G3=3.14159*G3 DO 3 I=1,NPTS T=TO-AUX(2)+DT*FLOAT(I-1) F(I)=0. IF(T.GT.0.AND.T.LE.AUX(1)) 1F(I)=SIN(G1*T/AUX(1))-G2*SIN(G3*T/AUX(1)) 3 CONTINUE RETURN 30 CONTINUE IF(IAUX(1).NE.4)GO TO 40 C C RICKER SIGNAL C DO 4 I=1,NPTS T=TO-AUX(2)+DT*FLOAT(I-1) G1=AUX(1)*T G2=2.*G1 G3=G1*G1 F(I)=0. IF(G3.LT.20.)F(I)=(1.-0.5*G2*G2)*EXP(-G3) 4 CONTINUE RETURN 40 CONTINUE IF(IAUX(1).NE.5)GO TO 50 C C BOX-CAR FUNCTION C DO 5 I=1,NPTS T=TO+DT*FLOAT(I-1) F(I)=0. IF(T.GE.AUX(1).AND.T.LE.AUX(2))F(I)=AUX(3) 5 CONTINUE RETURN 50 CONTINUE IF(IAUX(1).NE.6) GO TO 60 C C RAMP FUNCTION C DO 6 I=1,NPTS F(I)=0. T=TO+DT*FLOAT(I-1) IF(T.LT.AUX(3).AND.T.GE.AUX(2))F(I)=AUX(4) IF(T.LE.AUX(2).AND.T.GE.AUX(1)) 1F(I)=AUX(4)*(T-AUX(1))/(AUX(2)-AUX(1)) 6 CONTINUE RETURN 60 CONTINUE IF(IAUX(1).NE.7)GO TO 70 C C TRIANGLE FUNCTION C DO 7 I=1,NPTS F(I)=0. T=TO+DT*FLOAT(I-1) IF(T.LE.AUX(2).AND.T.GE.AUX(1)) 1F(I)=AUX(4)*(T-AUX(1))/(AUX(2)-AUX(1)) IF(T.LE.AUX(3).AND.T.GE.AUX(2)) 1F(I)=AUX(4)*(AUX(3)-T)/(AUX(3)-AUX(2)) 7 CONTINUE RETURN 70 CONTINUE IF(IAUX(1).NE.8)GO TO 80 C C ONE SAMPLE FUNCTION C DO 8 I=1,NPTS F(I)=0. IF(I.EQ.IAUX(2))F(I)=AUX(1) 8 CONTINUE RETURN 80 CONTINUE IF(IAUX(1).NE.9)GO TO 90 C C SQUARED SIN SIGNAL C DO 9 I=1,NPTS T=TO-AUX(2)+DT*FLOAT(I-1) G1=3.1415926*T/AUX(1) F(I)=0. IF(T.GT.0.AND.T.LE.AUX(1))THEN G2=SIN(G1) F(I)=G2*G2 END IF 9 CONTINUE RETURN 90 CONTINUE C RETURN END C C C SUBROUTINE FCOOLR(K,D,SN) C FAST FOURIER TRANSFORM OF N = 2**K COMPLEX DATA POINTS C REPARTS HELD IN D(1,3,...2N-1), IMPARTS IN D(2,4,...2N). C SUBROUTINE FCOOLR(K,D,SN) DIMENSION INU(15),D(4096) 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 SUBROUTINE REDUC(AREDUC,N1,N2,F,G,IS) DIMENSION F(1),G(1),IS(1) DO 1 I=N1,N2 IS(I)=0 1 G(I)=0. AREDUC=0. DO 2 I=N1,N2 2 IF(ABS(F(I)).GT.AREDUC)AREDUC=ABS(F(I)) IF(AREDUC.EQ.0.)GO TO 5 B=1./AREDUC C=999.1*B GO TO 6 5 B=0. C=0. 6 CONTINUE DO 3 I=N1,N2 G(I)=B*F(I) 3 IS(I)=IFIX(C*F(I)) RETURN END C C SUBROUTINE FILTER COMPLEX S(2048),REF(2048) COMMON/SPECTR/S,REF,DF,NI,NT,NT2,NTM,ILS,IRS,NFS,NFILT,NWIN RETURN END C C======================================================================= C INCLUDE 'error.for' C error.for INCLUDE 'calcops.for' C calcops.for C C======================================================================= C