C
C PROGRAM FRESAN C ************** C CHARACTER*80 MTEXT,PSTEXT,ITEXT,FILEIN,FILEOU,FILE1,FILE2,FILE3 COMPLEX IMAG,CTIME,CAUX,CAX,CAY,CAZ,CAX1,CAY1,CAZ1,CSOUR,AAUX DIMENSION ISEL(100),IEP(100),III(100),AMFR(1000),AMFI(1000), 1AAUX(3),P(3),POL(3),POL1(3) COMMON/PLOT1/RESP(30000),IFREQ(100),AMPL(100),DST(100),MCOMP COMMON/SOUR/ROS COMMON/SH/MSH,NSH C C C************************************************** C LIN=5 LOU=6 LU2=1 LU7=2 LU6=3 FILEIN='fres.dat' FILEOU='fres.out' FILE1='lu2.dat' FILE2='lu7.dat' FILE3='lu6.dat' WRITE(*,'(2A)') ' (FRESAN) SPECIFY NAMES OF INPUT AND OUTPUT', 1' FILES LIN, LOU, LU2, LU7, LU6: ' READ(*,*) FILEIN,FILEOU,FILE1,FILE2,FILE3 IF(FILE1.EQ.' ') GO TO 99 IF(FILE2.EQ.' ') LU7=0 IF(FILE3.EQ.' ') LU6=0 OPEN(LIN,FILE=FILEIN,FORM='FORMATTED',STATUS='OLD') OPEN(LOU,FILE=FILEOU,FORM='FORMATTED') OPEN(LU2,FILE=FILE1,FORM='FORMATTED',STATUS='OLD') IF(LU7.NE.0)OPEN(LU7,FILE=FILE2,FORM='FORMATTED') IF(LU6.NE.0)OPEN(LU6,FILE=FILE3,FORM='FORMATTED',STATUS='OLD') C C************************************************** C ITEXT='FRESAN' PSTEXT=' ' IPRINT=0 JPRINT=0 XSHIFT=3. YSHIFT=6. READ(LIN,*)ITEXT READ(LIN,*)IPRINT,JPRINT,PSTEXT WRITE(LOU,105)ITEXT WRITE(LOU,111)IPRINT,JPRINT,PSTEXT REWIND LU2 IF(LU7.NE.0)REWIND LU7 C NEPIC=0 NSPR=30000 IMAG=(0.,1.) PI=4.*ATAN(1.) C C READ INPUT PARAMETERS C MCOMP=0 MABS=0 MSOUR=0 MSELEC=0 MEPIC=0 MFR=0 MPLOT=0 READ(LIN,*)MCOMP,MABS,MSOUR,MSELEC,MEPIC,MFR,MPLOT WRITE(LOU,101)MCOMP,MABS,MSOUR,MSELEC,MEPIC,MFR,MPLOT C IF(MPLOT.NE.0)THEN IF(IPRINT.LT.0)THEN CALL PLOTN(PSTEXT,0) IPRINT=-IPRINT END IF CALL PLOTS(LDUM1,LDUM2,7) END IF C AROT=0. READ(LIN,*)RFR,GM,QRED,FREQM,AROT WRITE(LOU,100)RFR,GM,QRED,FREQM,AROT CSRT=COS(AROT) SNRT=SIN(AROT) IF(ABS(RFR).LT.0.00001)RFR=1. IF(ABS(QRED).LT.0.00001)QRED=1. C C FREQUENCY FILTERING DATA C NFREQ=0 READ(LIN,*)NFREQ,FL,FR,FD WRITE(LOU,107)NFREQ,FL,FR,FD IF(NFREQ.GT.0)FD=1./FD/FLOAT(NFREQ) NF1=INT(FL/FD-0.5) NF=INT(FR/FD+0.5)-NF1 FL=FD*FLOAT(NF1)+FD FR=FL+FD*FLOAT(NF)-FD FDA=6.2831853*FD C C SELECTION OF RECEIVERS AND ELEMENTARY WAVES C IF(MEPIC.EQ.0)GO TO 1 READ(LIN,*)NEPIC,(IEP(I),I=1,NEPIC) WRITE(LOU,101)NEPIC,(IEP(I),I=1,NEPIC) 1 IF(MSELEC.EQ.0)GO TO 2 READ(LIN,*)NSELEC,(ISEL(I),I=1,NSELEC) WRITE(LOU,101)NSELEC,(ISEL(I),I=1,NSELEC) 2 CONTINUE C C READING FROM LU2,GENERAL QUANTITIES C READ(LU2,104)MTEXT READ(LU2,106)NDST,KSH,ILOC READ(LU2,100)XSOUR,YSOUR,ZSOUR,TSOUR,RSTEP,ROS IF(NDST.NE.1)READ(LU2,100)(DST(I),I=1,NDST) IF(NDST.EQ.1)READ(LU2,100)XREC,YREC,XPRF,XATAN WRITE(LOU,104)MTEXT WRITE(LOU,101)NDST,KSH,ILOC WRITE(LOU,100)XSOUR,YSOUR,ZSOUR,TSOUR,RSTEP,ROS IF(NDST.NE.1)WRITE(LOU,100)(DST(I),I=1,NDST) IF(NDST.EQ.1)WRITE(LOU,100)XREC,YREC,XPRF,XATAN IF(NDST.EQ.1)DST(1)=XPRF NRDST=NDST-NEPIC C C INITIAL PARAMETERS FOR RADIATION PATTERNS C AND FOR FREQUENCY DEPENDENT FACTORS C IF(MSOUR.NE.0)CALL SOURCE(LIN,LOU,0,0,MSOUR,P,POL,AMSOUR,PHSOUR) IF(MFR.NE.0)CALL FDEP(0,G,FL,FD,NF,AMFR,AMFI) IF(LU6.NE.0)CALL FDQI(LU6,LOU,0,MCOMP,FL,FD,NF,AROT,AMFR,AMFI) IF(LU6.NE.0)THEN NFREQ=0 NF1=INT(FL/FD-0.5) FR=FL+FD*FLOAT(NF)-FD FDA=6.2831853*FD END IF C C C WRITING ON LU7 C WRITE(LU7,104)MTEXT WRITE(LU7,104)ITEXT WRITE(LU7,152)XSOUR,YSOUR,ZSOUR,TSOUR,RSTEP,FL,FD WRITE(LU7,106)NRDST,NF,NF1,MCOMP,ILOC C C CHECK OF THE DIMENSIONS OF THE RESP FIELD C INITIAL SPECIFICATION OF THE RESP FIELD C NSUM=2*NDST*NF IF(NSUM.GT.NSPR)WRITE(LOU,153) IF(NSUM.GT.NSPR)GO TO 99 DO 4 I=1,NSUM 4 RESP(I)=0. C C S T A R T O F S U M M A T I O N C O F F R E Q U E N C Y R E S P O N S E C **************************************** C C LOOP OVER TERMINATION POINTS OF ALL ELEMENTARY WAVES C AT ALL RECEIVERS C 6 CONTINUE READ(LU2,103,END=84)NC,IDIST,TIME,CAX,CAY,CAZ,TAST IF(NC.LT.0)READ(LU2,109)CAX1,CAY1,CAZ1 READ(LU2,110)(P(I),I=1,3) READ(LU2,110)(POL(I),I=1,3) IF(NC.LT.0)READ(LU2,110)(POL1(I),I=1,3) CAUX=CAX*CSRT+CAY*SNRT CAY=-CAX*SNRT+CAY*CSRT CAX=CAUX CAUX=CAX1*CSRT+CAY1*SNRT CAY1=-CAX1*SNRT+CAY1*CSRT CAX1=CAUX NC1=IABS(NC) IF(MSELEC.EQ.0)GO TO 8 C C SELECTION OF WAVES C DO 7 J=1,NSELEC IF(ISEL(J).EQ.NC1)GO TO 6 7 CONTINUE 8 IF(MEPIC.EQ.0)GO TO 10 C C SELECTION OF RECEIVERS C DO 9 J=1,NEPIC IF(IEP(J).EQ.IDIST)GO TO 6 9 CONTINUE 10 CONTINUE C C IF(MSOUR.NE.0)THEN CALL SOURCE(LIN,LOU,1,1,MSOUR,P,POL,AMDIR,PHDIR) CSOUR=AMDIR*CEXP(IMAG*PHDIR) END IF IF(MSOUR.EQ.0)CSOUR=(1.,0.) AAUX(1)=CAX*CSOUR AAUX(2)=CAY*CSOUR AAUX(3)=CAZ*CSOUR IF(NC.LT.0)THEN IF(MSOUR.NE.0)THEN CALL SOURCE(LIN,LOU,1,1,MSOUR,P,POL1,AMDIR,PHDIR) CSOUR=AMDIR*CEXP(IMAG*PHDIR) END IF IF(MSOUR.EQ.0)CSOUR=(1.,0.) AAUX(1)=AAUX(1)+CAX1*CSOUR AAUX(2)=AAUX(2)+CAY1*CSOUR AAUX(3)=AAUX(3)+CAZ1*CSOUR END IF IF(MCOMP.EQ.0)K=3 IF(MCOMP.EQ.1)K=1 IF(MCOMP.EQ.2)K=2 RW=REAL(AAUX(K)) CW=AIMAG(AAUX(K)) IF(MFR.NE.0)THEN CSOUR=CSQRT(CAX*CAX+CAY*CAY+CAZ*CAZ) RW=REAL(CSOUR) CW=AIMAG(CSOUR) END IF AR=SQRT(RW*RW+CW*CW) IF(ABS(RW).LT..000001.AND.ABS(CW).LT..000001)THEN PH=0. ELSE PH=ATAN2(CW,RW) END IF IF(IPRINT.EQ.1)WRITE(LOU,102)NC,IDIST,TIME,AR,PH,AMDIR,TAST C C ABSORPTION EFFECTS C IF(MABS.EQ.0)GO TO 67 AUX=QRED*TAST ABS1=.5*AUX IF(MABS.LE.2)GO TO 67 ABS2=1.+ALOG(FREQM/RFR) PH=PH+2.*AUX*FREQM 67 CONTINUE C C FREQUENCY DEPENDENT EFFECTS C KFR=0 IF(MABS.LE.1.AND.MFR.EQ.0)GO TO 74 KFR=1 DO 73 JF=1,NF AMFR(JF)=1. 73 AMFI(JF)=0. C C FREQUENCY DEPENDENT EFFECTS: SUBROUTINE FDEP C IF(MFR.NE.0)CALL FDEP(1,G,FL,FD,NF,AMFR,AMFI) C C PROCESSING OF FREQUENCY-DEPENDENT QI AMPLITUDES C IF(LU6.NE.0)CALL FDQI(LU6,LOU,1,MCOMP,FL,FD,NF,AROT,AMFR,AMFI) C C FREQUENCY DEPENDENT EFFECTS:CAUSAL ABSORPTION C IF(MABS.LE.1.OR.MABS.EQ.3)GO TO 74 DO 75 JF=1,NF FRJ=FD*FLOAT(JF+NF1) ABS4=2.*AUX*FRJ*ALOG(FRJ/RFR) CAUX=(0.,0.) IF(MABS.EQ.2)CAUX=CEXP(-IMAG*ABS4) AUXR=REAL(CAUX) AUXI=AIMAG(CAUX) AUX1=AMFR(JF)*AUXR-AMFI(JF)*AUXI AUX2=AMFR(JF)*AUXI+AMFI(JF)*AUXR AMFR(JF)=AUX1 AMFI(JF)=AUX2 75 CONTINUE C C FFR COMPUTATION C 74 CONTINUE SPH=SIN(PH) CPH=COS(PH) ARE=AR*CPH AIM=AR*SPH C C PRINT OF REAL AND IMAGINARY PARTS OF THE COMPLEX AMPLITUDE C IF(IPRINT.EQ.2)WRITE(LOU,154)ARE,AIM CTIME=TIME IF(MABS.EQ.0)GO TO 83 CTIME=CTIME+ABS1*IMAG IF(MABS.EQ.3)CTIME=CTIME-ABS2*AUX/PI 83 TR=REAL(CTIME) TI=AIMAG(CTIME) TID=TI*FDA TRD=TR*FDA TIE=TID*FLOAT(NF1) TRE=TRD*FLOAT(NF1) AUX=EXP(-TIE) ARC=AUX*COS(TRE) ARS=AUX*SIN(TRE) AA1=ARC*ARE-ARS*AIM AA2=ARC*AIM+ARS*ARE AUX=EXP(-TID) AR=AUX*COS(TRD) AI=AUX*SIN(TRD) C C AR, AI - CONSTANT FACTORS FOR RECURENT RESPONSE COMPUTATION C AA1, AA2 - REAL AND IMAGINARY PARTS OF THE AMPLITUDE C CORRESPONDING TO THE LOWEST CONSIDERED FREQUENCY C IF(IPRINT.EQ.3)WRITE(LOU,100)TR,TI,AR,AI,AA1,AA2 C C LOOP FOR FREQUENCIES C N=2*(IDIST-1)*NF-1 DO 82 JF=1,NF IF(JF.NE.1)THEN A11=AR*AA1-AI*AA2 AA2=AR*AA2+AI*AA1 AA1=A11 END IF AUX=AA1 AUXX=AA2 IF(KFR.EQ.0)GO TO 81 AUX=AA1*AMFR(JF)-AA2*AMFI(JF) Auxx=AA1*AMFI(JF)+AA2*AMFR(JF) 81 CONTINUE N=N+2 RESP(N)=RESP(N)+AUX RESP(N+1)=RESP(N+1)+AUXX 82 CONTINUE C C END OF LOOP FOR FREQUENCIES C GO TO 6 C C E N D O F C O M P U T A T I O N O F C F R E Q U E N C Y R E S P O N S E C ***************************************** C C PRINTING AND STORING THE FREQUENCY RESPONSE C 84 CONTINUE DO 95 IR=1,NDST K=2*IR*NF J=K-2*NF+1 A=0.1E-20 DO 91 I=J,K 91 A=AMAX1(A,ABS(RESP(I))) WRITE(LU7,160)DST(IR),A IF(IPRINT.EQ.1)WRITE(LOU,108)DST(IR),A A=A/99999.1 92 DO 93 I=1,12 III(I)=RESP(J)/A IF(J.GE.K)GO TO 94 93 J=J+1 WRITE(LU7,161)(III(I),I=1,12) IF(IPRINT.EQ.2)WRITE(LOU,161)(III(I),I=1,12) GO TO 92 94 WRITE(LU7,161)(III(K),K=1,I) IF(IPRINT.EQ.2)WRITE(LOU,161)(III(K),K=1,I) 95 CONTINUE REWIND LU7 REWIND LU2 C C PRINTING OF FREQUENCY DEPENDENT AMPLITUDE CURVES C IF(JPRINT.EQ.0)GO TO 350 READ(LIN,*)MTAB,(IFREQ(I),I=1,MTAB) WRITE(LOU,101)MTAB,(IFREQ(I),I=1,MTAB) DO 310 I=1,MTAB IFF=IFREQ(I) AMAX=0. IF(IFF.GT.NF)GO TO 310 FREQ=FD*FLOAT(NF1+IFF) DO 305 IR=1,NDST A=RESP(2*(IR-1)*NF+2*IFF-1) B=RESP(2*(IR-1)*NF+2*IFF) AMPL(IR)=SQRT(A*A+B*B) IF(AMPL(IR).GT.AMAX)AMAX=AMPL(IR) 305 CONTINUE B=AMAX/999.1 DO 306 IJ=1,NDST IF(ABS(B).GT..0000001)III(IJ)=AMPL(IJ)/B IF(ABS(B).LT..0000001)III(IJ)=0 306 CONTINUE WRITE(LOU,170)FREQ,AMAX WRITE(LOU,171)(III(K),K=1,NDST) 310 CONTINUE 350 CONTINUE C C PLOTTING OF FREQUENCY DEPENDENT AMPLITUDE CURVES C IF(MPLOT.EQ.0)GO TO 99 CALL ZPLOT(NF1,NF,FD,NDST,MPLOT,LIN,LOU) CALL PLOT(0.,0.,999) C C 100 FORMAT(8F10.5) 101 FORMAT(16I5) 102 FORMAT(2I3,F10.5,1X,2(1X,E12.6),3F10.6) 103 FORMAT(2I3,F12.7,6E12.6,F10.6) 104 FORMAT(A) 105 FORMAT(1X,A) 106 FORMAT(26I3) 107 FORMAT(I5,6F10.5) 108 FORMAT(F10.3,E12.5) 109 FORMAT(6E12.6) 110 FORMAT(3F10.5) 111 FORMAT(2I5,1X,A) 152 FORMAT(5F10.5,2E15.7) 153 FORMAT(1X,'DIMENSION OF RESP EXCEEDED') 154 FORMAT(3E16.8) 160 FORMAT(F10.3,E12.5) 161 FORMAT(12I6) 170 FORMAT(/,1X,'FREQUENCY=',F10.5,2X,'(HZ), MAXIMUM AMPLITUDE=', 11E15.6) 171 FORMAT(20I4) C C 99 CONTINUE STOP END C C SUBROUTINE FDQI(LU6,LOU,NREAD,MCOMP,FL,FD,NF,AROT,AMFR,AMFI) C COMPLEX W(2) DIMENSION AMFR(1000),AMFI(1000),E1(3),E2(3) C IF(NREAD.EQ.0)READ(LU6,101)FL,FD,NF IF(NREAD.EQ.0)RETURN C IF(MCOMP.EQ.0)M=3 IF(MCOMP.EQ.1)M=1 IF(MCOMP.EQ.2)M=2 READ(LU6,100)(E1(K),K=1,3),(E2(K),K=1,3) IF(M.NE.3.AND.AROT.NE.0.)THEN CSRT=COS(AROT) SNRT=SIN(AROT) EAUX=E1(1)*CSRT+E1(2)*SNRT E1(2)=-E1(1)*SNRT+E1(2)*CSRT E1(1)=EAUX EAUX=E2(1)*CSRT+E2(2)*SNRT E2(2)=-E2(1)*SNRT+E2(2)*CSRT E2(1)=EAUX END IF E1M=E1(M) E2M=E2(M) DO 1 I=1,NF READ(LU6,100)FF,W AMFR(I)=REAL(W(1)*E1M+W(2)*E2M) AMFI(I)=AIMAG(W(1)*E1M+W(2)*E2M) 1 CONTINUE C 100 FORMAT(16E15.5) 101 FORMAT(2F10.5,I5) C RETURN END C C SUBROUTINE FDEP(N,G,FL,DF,NF,AMFR,AMFI) DIMENSION AMFR(1000),AMFI(1000) IF(N.EQ.0)RETURN DO 1 I=1,NF AMFR(I)=1. AMFI(I)=0. 1 CONTINUE RETURN END C C SUBROUTINE ZPLOT(NF1,NF,FD,NDST,MPLOT,lin,lou) CHARACTER*80 MTEXT COMMON/PLOT1/RESP(30000),IFREQ(100),AMPL(100),DST(100),MCOMP C IRUN=0 SHIFT=5. C 2 READ(LIN,*)MTEXT WRITE(LOU,105)MTEXT SC=1. READ(LIN,*)XMIN,XMAX,XLEN,DTICX,SC WRITE(LOU,100)XMIN,XMAX,XLEN,DTICX,SC IF(ABS(XLEN).LT..00001)GO TO 99 READ(LIN,*)AMIN,AMAX,ALEN,DTICA WRITE(LOU,100)AMIN,AMAX,ALEN,DTICA NLOG=0 NTICX=1 NTICA=1 NDX=0 NDA=0 READ(LIN,*)NLOG,NTICX,NTICA,NDX,NDA WRITE(LOU,101)NLOG,NTICX,NTICA,NDX,NDA XMER=XLEN/(XMAX-XMIN) AMER=ALEN/(AMAX-AMIN) CALL PLOT(5.,5.,-3) C C LOOP FOR PLOTS C 1 READ(LIN,*) MTAB,(IFREQ(I),I=1,MTAB) WRITE(LOU,101)MTAB,(IFREQ(I),I=1,MTAB) IF(MTAB.EQ.0)GO TO 2 IF(IRUN.EQ.1)CALL PLOT(XLEN+SHIFT,0.,-3) IRUN=1 CALL BORDER(XLEN,DTICX,ALEN,DTICA,SC,MTEXT,0,XMIN,XMAX, 1AMIN,AMAX,NTICX,NTICA,NDX,NDA) C C LOOP FOR FREQUENCIES C DO 10 I=1,MTAB IFF=IFREQ(I) AMAXIM=0. IF(IFF.GT.NF)GO TO 10 FREQ=FD*FLOAT(NF1+IFF) DO 5 IR=1,NDST A=RESP(2*(IR-1)*NF+2*IFF-1) B=RESP(2*(IR-1)*NF+2*IFF) AMPL(IR)=SQRT(A*A+B*B) IF(AMPL(IR).GT.AMAXIM)AMAXIM=AMPL(IR) IF(NLOG.EQ.1)AMPL(IR)=ALOG10(AMPL(IR)) 5 CONTINUE IF(NLOG.EQ.1)AMAXIM=ALOG10(AMAXIM) IEXP=3 DO 8 IR=1,NDST XNEW=(DST(IR)-XMIN)*XMER IF(XNEW.LT.0..OR.XNEW.GT.XLEN)GO TO 8 A=AMPL(IR) IF(MPLOT.EQ.2)A=A/AMAXIM YNEW=(A-AMIN)*AMER IF(IR.EQ.1)GO TO 11 IF(YNEW.LT.0..OR.YNEW.GT.ALEN)GO TO 6 IF(YOLD.GE.0..AND.YOLD.LE.ALEN)GO TO 7 IF(YOLD.LT.0.)YB=0. IF(YOLD.GT.ALEN)YB=ALEN IEXP=2 4 AUX1=(DST(IR)-DST(IR-1))*XMER AUX2=YNEW-YOLD XOLD=XNEW-(YNEW-YB)*AUX1/AUX2 CALL PLOT(XOLD,YB,IEXP) IF(IEXP.EQ.3)GO TO 12 IEXP=3 GO TO 7 6 IF(YOLD.GT.ALEN)GO TO 12 IF(YNEW.LT.0.)YB=0. IF(YNEW.GT.ALEN)YB=ALEN IEXP=3 GO TO 4 11 IF(YNEW.LT.0..OR.YNEW.GT.ALEN)GO TO 12 7 CALL PLOT(XNEW,YNEW,IEXP) 12 YOLD=YNEW IEXP=2 8 CONTINUE U=(4.5+1.2*FLOAT(I))*SC UU=(4.5+1.8*FLOAT(I))*SC CALL NUMBER(U,ALEN+.2*SC,.3*SC,FREQ,0.,1) IF(MPLOT.EQ.2)CALL NUMBER(UU,ALEN+.7*SC,.3*SC,AMAXIM,0.,5) 10 CONTINUE CALL SYMBOL(.2*SC,ALEN+.2*SC,.3*SC,'FREQUENCIES(HZ)=',0.,16) ELM=.45*SC T=.5*(XLEN-6.3*SC) CALL SYMBOL(T,-1.6*SC,ELM,'DISTANCE IN KM',0.,14) T=.5*(ALEN-7.65*SC) U=-(1.6+.4*NDA)*SC CALL SYMBOL(U,T,ELM,'A M P L I T U D E',90.,17) IF(MCOMP.EQ.0) 1CALL SYMBOL(.2*SC,ALEN+1.5*SC,ELM,'VERTICAL',0.,8) IF(MCOMP.EQ.1) 1CALL SYMBOL(.2*SC,ALEN+1.5*SC,ELM,'X-COMPONENT',0.,11) IF(MCOMP.EQ.2) 1CALL SYMBOL(.2*SC,ALEN+1.5*SC,ELM,'Y-COMPONENT',0.,11) IF(MPLOT.EQ.2) 1CALL SYMBOL(.2*SC,ALEN+.7*SC,.3*SC,'REDUCTION FACTOR=',0.,17) GO TO 1 C C THE END OF THE LOOP FOR FREQUENCIES C 100 FORMAT(8F10.5) 101 FORMAT(16I5) 104 FORMAT(A) 105 FORMAT(1X,A) 99 CONTINUE RETURN END C C======================================================================= C INCLUDE 'source.for' C source.for INCLUDE 'border.for' C border.for INCLUDE 'error.for' C error.for INCLUDE 'length.for' C length.for INCLUDE 'sep.for' C sep.for INCLUDE 'calcops.for' C calcops.for C C======================================================================= C