IMPLICIT REAL *8(A-B,D-H,O-Z),COMPLEX *16(C) REAL *8 KE C********************************************************************** C* NAME PASELAG * C* * C* TYPE MAIN * C* * C* PURPOSE TO COMPUTE SECONDARY PHASELAG USING BOTH SPHERICAL * C* AND PLANE EARTH THEORIES * C* * C* AUTHOR D.GRAY 1977/78 * C* MODIFIER T.MAHMOOD OCT.1981 * C* * C* EXTERNALS TAUS,PHASE,PHASE2,PHASE0,DCBRT * C* * C* PARAMETERS ICODE =A CODE TO DETERMINE THE NEXT STEP IN THE * C* THE COMPUTATION * C* =98 IF THERE IS ANY SET OF DATA UNDER DIFFERENT * C* TITLE * C* =99 IF THERE IS NO MORE DATA * C* =0 IF THE DATA IS IN THE SAME SET * C* FREQ =FREQUENCY(KHZ) * C* V VELOCITY(METERS/SEC) * C* AVCON =CONDUCTIVITY OF EARTH OR WATER (MHOS/METER) * C* E1 =PERMITTIVITY OF ATMOSPHERE (E.S.U) * C* E2 =PERMITTIVITY OF EARTH OR WATER (E.S.U) * C* A =RADIUS OF EARTH (METRES) * C* ALF =VERTICAL LAPSE RATE OF PERMITTIVITY OF AIR * C* FF =FACTOR TO BE MULTIPLIED TO THE READ-IN FREQUENCY* C* FS =FACTOR TO BE MULTIPLIED TO THE READ-IN CONDUCT * C* FE2 =FACTOR TO BE MULTIPLIED TO THE READ-IN E2 * C* W =ANGULAR FREQUENCY =2*PI*FREQ (RADDIANS/SEC) * C* CDE =THE CONDUCTIVITY AND PERMITTIVITY PARAMETERS * C* FOR A VERTICAL DIPOLE SOURCE * C* KE =THE MODULUS OR THE AMPLITUDE OF CDE * C* CTS =TAU S ; SOLUTION OF RICCATI EQUATION * C* DISTM =DISTANCE (METERS) * C* PLSPHR=PHASELAG COMPUTED USING SPHERICAL THEORY * C* PLPLNE=PHASELAG COMPUTED USING PLANE THEORY * C* KEYSPH=NUMBER OF TERMS USED TO COMPUTE PLSPHR * C* KEYPLN=NUMBER OF TERMS USED TO COMPUTE PLPLNE * C* DPLSPH=INCREAMENT OF PLSPHR WITH DISTANCE * C* DIFPLG=DIFFERENCE BETWEEN PLANE AND SPHERICAL PHASELAG * C* NTESTS=TESTING THE IDENTICALITY OF THE TWO PHASELAGS * C* =8888888 IF THE DIFFERENCE IS LESS THAN 10 METERS* C* =11111 IF OTHERWISE * C* ERROR =FRACTIONAL ERROR;SPHERICAL PHASELAG:DISTANCE * C* * C* LANGUAGE FORTRAN * C* * C* REFRENCES JOHLER ET AL (1956) ;NBS CIRCULAR 573 * C* HOWE (1960) : NBS JOURNAL,VOL.64B,#2 * C* LIMITATION THIS PROGRAM IS NOT SUITABLE FOR LOW GROUND CONDUCT. * C* 0.001 FOR 100KHZ FREQ. IS THE LIMIT AND IS HIGHER WITH * C* HIGHER FREQUENCIES * C* * C* NOTES 1. THE FACTOR CARD IS USED SO THAT ONLY ONE CARD * C* (FACTOR) NEEDS TO BE CHANGED IF OTHER INPUT PARAME- * C* TERS ARE TO BE CHANGED.OTHERWISE1.0 IS USED FOR ALL * C* 2.IN THE OUTPUT: PLANE PHASELAG SHOULD BE USED IF IT * C* IS GREATER THAN THE SPHERICAL PHASELAG * C********************************************************************** DIMENSION TITLE(20),CTS(2500),TABLE(219),DISTM(219),PLSPHR(219), & KEYSPH(219),PLPLNE(219),KEYPLN(219),DPLSPH(219), & DIFPLG(219),NTESTS(219),ERROR(219) IR =5 IW =6 E1 =1.000676D0 A =6.367390D+6 ALF=0.75D0 V =2.9979250D+08 C C.....READ IN DATA :1) HEADING, 2) FACTOR CARD, 3) PARAMETERS C 10 CONTINUE READ(IR,1001) (TITLE(I),I=1,20) READ(IR,1002) FF,FS,FE2 20 READ(IR,1003) ICODE,FREQ,AVCON,E2 IF(ICODE.EQ.98) GO TO 10 IF(ICODE.EQ.99) GO TO 100 FREQ =FREQ*FF AVCON=AVCON*FS E2 =E2*FE2 PI =4.D0*DATAN(1.D0) W =2.D0*PI*FREQ SMC =(AVCON*PI*V*V*4.0D-7)/W KE =(DCBRT((V*ALF)/(W*A*E1*E1))*DSQRT(SMC*SMC+E2*E2))/ & DSQRT(DSQRT((E2-E1)*(E2-E1)+SMC*SMC)) PSI =DATAN(E2/SMC)-0.5D0*DATAN((E2-E1)/SMC) PP =PI*0.75D0-PSI CPPI =DCMPLX(0.D0,PP) CALL TRAPS(0,0,100,0,0) CDE =DCMPLX(KE,0.D0)*CDEXP(CPPI) C C.....PRINTING THE HEADING AND THE PARAMETERS USED IN THE COMPUTATION C WRITE(IW,1004)(TITLE(L),L=1,20),FREQ,AVCON,E2 C C.....COMPUTATION OF TAU S C CALL TAUS(CDE,CTS) WRITE(6,1005) ND =219 FYP=0.D0 C C.....COMPUTATION OF SECONDARY PHASELAG C DO 50 K=1,ND IF(K.LT.20) DIST=K*1000.D0 IF(K.GE.20) DIST=20000.D0+(K-20)*10000.D0 DISTM(K) =DIST CALL PHASE(DIST,FREQ,V,E1,ALF,A,CDE,CTS,FYP,FYCT,FYT,KEY) PLSPHR(K)=FYT KEYSPH(K)=KEY ERROR(K) =DISTM(K)/PLSPHR(K) CALL PHASE2(DIST,FREQ,V,AVCON,E1,E2,ALF,A,FYLAG2,KEY2) PLPLNE(K)=FYLAG2 KEYPLN(K)=KEY2 DIFPLG(K)=PLPLNE(K)-PLSPHR(K) NTESTS(K)=11111 IF(DABS(DIFPLG(K)).LT.10.D0) NTESTS(K)=8888888 IF(K.GT.1) GO TO 30 DPLSPH(K)=0.D0 GO TO 40 30 DPLSPH(K)=PLSPHR(K)-PLSPHR(K-1) 40 FYP=FYCT WRITE(6,1006)DISTM(K),PLSPHR(K),DPLSPH(K),PLPLNE(K),DIFPLG(K), & NTESTS(K),KEYSPH(K),KEYPLN(K),ERROR(K) 50 CONTINUE GO TO 20 1001 FORMAT(1X,A3,19A4) 1002 FORMAT(3(F10.5,10X)) 1003 FORMAT(I3,F10.1,2F10.5) 1004 FORMAT('1',1H1/10X,A3,19A4//10X,10H FREQUENCY,F13.1/10X, & 13H CONDUCTIVITY,F9.5/10X,13H PERMITTIVITY,F10.5/) 1005 FORMAT(///5X,'DISTANCE',6X,'PHASE LAG',6X,'INCREMENT',5X, & 'PHASE LAG',5X,'DIFFERENCE',4X,'TEST',4X,'NO OF TERMS', & 4X,'FRACTIONAL'/5X,'(METER)',7X,'(SPHERE)',8X,'(SPHERE)', & 6X,'(PLANE)',6X,'(PLN-SPH)',13X,'(SP)',3X,'(PL)', & 7X,'ERROR'//) 1006 FORMAT(5X,F9.1,6X,F6.1,9X,F6.1,8X,F6.1,8X,F7.1,4X,I7, & 3X,I4,2X,I4,F13.1) 100 STOP END SUBROUTINE TAUS(CDE,CTS) C********************************************************************** C* NAME TAUS * C* * C* TYPE SUBROUTINE * C* * C* PURPOSE TO COMPUTE TAU S USING POWER SERIES: * C* EITHER ASCENDING OR DESCENDING SERIES * C* * C* AUTHOR D.GRAY 1977/78 * C* MODIFIER T.MAHMOOD OCT.1981 * C* * C* EXTERNALS NONE * C* * C* CALLING CALL TAUS(CDE,CTS) * C* * C* PARAMETERS ALL EXPLAINED IN THE MAIN * C* INPUT :CDE * C* OUTPUT:CTS * C* * C* LANGUAGE FORTRAN * C* * C* REFRENCES JOHLERS ET AL (1956) : NBS CIRCULAR 573 * C* HOWE (1960) : NBS JOURNAL,VOL.64B,#2 * C********************************************************************** IMPLICIT REAL *8(A-B,D-H,O-Z),COMPLEX *16(C) DIMENSION CTS(2500), CB(500), CD(500) IW=6 PI=4.D0*DATAN(1.D0) DO 28 N=1,2500 CTS(N)=DCMPLX(0.D0,0.D0) CTSOLD=DCMPLX(0.D0,0.D0) C C.... START OF DESCENDING (RECURRENCE) SERIES C KEY2=1 S =DFLOAT(N-1) Y =0.375D0*PI*(4.D0*S+1.D0) TAU =((Y*Y/2.D0)**(1.D0/3.D0))*(1.D0-(7.D0/(48.D0*Y*Y))+ & (35.D0/(288.D0*Y*Y*Y*Y))) IF(N .EQ. 1) TAU=0.808616516D0 IF(N .EQ. 2) TAU=2.57809613D0 IF(N .EQ. 3) TAU=3.82571528D0 IF(N .EQ. 4) TAU=4.89182029D0 IF(N .EQ. 5) TAU=5.85130097D0 IF(N .EQ. 6) TAU=6.73731638D0 UR =TAU*DCOS(PI/3.D0) UIM =TAU*DSIN(PI/3.D0) DMIN =1.0D-08 DMAX =TAU*5.D0 CD(1) =DCMPLX(UR,UIM) CD(2) =-0.5D0/(CD(1)*CDE) CTS(N)=CD(1)+CD(2) DO 22 J=3,500 KEY =J CSUM=-DFLOAT(J-3)*CD(J-2)/(DFLOAT(J-1)*CDE*CDE) DO 21 K=3,J CSUM=CSUM+CD(K-1)*CD(J-K+2) 21 CONTINUE CD(J) =-0.5D0*CSUM/CD(1) CTS(N)=CTS(N)+CD(J) DM =CDABS(CD(J)) C C..... IF DIVERGENT -- GOES TO ASCENDING SERIES C IF(DM .GT. DMAX) GO TO 23 IF(DM .LT. DMIN) GO TO 28 22 CONTINUE CTSOLD=CTS(N) C C..... START OF ASCENDING (RECURRENCE) SERIES C 23 S =DFLOAT(N-1) KEY2=2 X =0.375D0*PI*(4.D0*S+3.D0) TAU0=((X*X/2.D0)**(1.D0/3.D0))*(1.D0+(5.D0/(48.D0*X*X))- & (5.D0/(36.D0*X*X*X*X))) IF(N .EQ. 1) TAU0=1.85575708D0 IF(N .EQ. 2) TAU0=3.24460762D0 IF(N .EQ. 3) TAU0=4.38167124D0 IF(N .EQ. 4) TAU0=5.38661378D0 IF(N .EQ. 5) TAU0=6.30526301D0 IF(N .EQ. 6) TAU0=7.16128272D0 BMIN =1.0D-08 BMAX =5.D0*TAU0 UR =TAU0*DCOS(PI/3.D0) UIM =TAU0*DSIN(PI/3.D0) CB(1) =DCMPLX(UR,UIM) CB(2) =-CDE CB(3) =DCMPLX(0.D0,0.D0) CTS(N)=CB(1)+CB(2) DO 25 J=4,500 KEY =J CSUM=DCMPLX(0.D0,0.D0) DO 24 K=3,J CSUM=CSUM+CB(K-2)*CB(J-K+1) 24 CONTINUE CB(J) =DFLOAT(J-3)*CDE*CDE*CSUM/DFLOAT(J-1) CTS(N)=CTS(N)+CB(J) BM =CDABS(CB(J)) C C..... IF DIVERGENT -- GIVES WARNING C IF(BM.GT.BMAX) GO TO 26 IF(BM.LT.BMIN) GO TO 28 25 CONTINUE 26 WRITE(IW,2001) N IF(CDABS(CTSOLD) .NE. 0.D0) GO TO 27 GO TO 28 27 CTS(N)=CTSOLD KEY2 =1 C C..... IF BOTH SERIES DIVERGENT -- USES DESCENDING SERIES C WRITE(IW,2002) N 28 CONTINUE 2001 FORMAT( 1H ,8H TAU S (,I4,20H ) DOES NOT CONVERGE ) 2002 FORMAT(1H , 8H TAU S (, I4,27H ) DESCENDING SERIES USED ) RETURN END SUBROUTINE PHASE (DIST,FREQ,V,E1,ALF,A,CDE,CTS,FYP,FYCT,FYT,KEY) C********************************************************************** C* NAME PHASE * C* * C* TYPE SUBROUTINE * C* * C* PURPOSE TO COMPUTE SECONDARY PHASELAG USING SPHERICAL * C* EARTH THEORY * C* * C* AUTHOR D.GRAY 1977/78 * C* MODIFIER T.MAHMOOD * C* * C* EXTERNALS DCBRT * C* * C* CALLING CALL PHASE(DIST,FREQ,V,E1,ALF,A,CDE,CTS,FYP, * C* FYCT,FYT,KEY) * C* * C* PARAMETERS ALL EXPLAINED IN THE MAIN * C* INPUT :DIST,FREQ,V,E1,ALF,A,CDE,CTS,FYP * C* OUTPUT:FYCT,FYT,KEY * C* * C* LANGUAGE FORTRAN * C* * C* REFRENCES JOHLER ET AL (1956) : NBS CIRCULAR 573 * C********************************************************************** IMPLICIT REAL *8(A-B,D-H,O-Z),COMPLEX *16(C) REAL *8 K1 DIMENSION CTS(2500) PI =4.D0*DATAN(1.D0) W =2.D0*PI*FREQ K1 =(W/V)*DSQRT(E1) DBIG =100000.D0 DCORR=1.D0/(K1*DIST) FRC =DCBRT(K1*A)*DCBRT(ALF*ALF) FWC =DSQRT(FRC*2.D0*PI*DIST/A) FYI =0.D0 MLT =0 CSUM =DCMPLX(0.D0,0.D0) DO 31 N=1,2500 KEY =N CFR2 =((FRC*CTS(N)*DIST)/A)+(DIST/A)*ALF*0.5D0+0.25D0*PI CFR2C=CFR2*DCMPLX(0.D0,1.D0) CALL TRAPS(0,0,100,0,0) CSUM =CSUM+CDEXP(CFR2C)/(CTS(N)*2.D0-(1.D0/CDE**2)) CFR =CSUM*FWC FRR =DREAL(CFR)-DCORR*DCORR FRI =DIMAG(CFR)+DCORR FYCT =DATAN2(FRI,FRR) IF(DABS(FYCT-FYI).LT.0.000001D0) MLT=MLT+1 IF(MLT.GT.5) GO TO 32 FYI=FYCT 31 CONTINUE 32 IF(FYCT.LT.0.D0) FYCT=FYCT+2.D0*PI IF(DIST.LT.DBIG) GO TO 33 IF(FYCT.LT.FYP) FYCT=FYCT+2.D0*PI IF(FYCT.LT.FYP) FYCT=FYCT+2.D0*PI IF(FYCT.LT.FYP) FYCT=FYCT+2.D0*PI 33 FYT=FYCT*V/W RETURN END SUBROUTINE PHASE2 (DIST,FREQ,V,AVCON,E1,E2,ALF,A,FYLAG2,KEY2) C********************************************************************** C* NAME PHASE2 * C* * C* TYPE SUBROUTINE * C* * C* PURPOSE TO COMPUTE SECONDARY PHASELAG USING PLANE EARTH THEORY * C* * C* AUTHOR D.GRAY 1977/78 * C* MODIFIER T.MAHMOOD * C* * C* EXTERNALS DCBRT * C* * C* CALLING CALL PHASE2(DIST,FREQ,V,AVCON,E1,E2,ALF,A,YFLAG2,KEY2 * C* PARAMETERS ALL EXPLAINED IN THE MAIN * C* INPUT :DIST,FREQ,V,AVCON,E1,E2,ALF,A * C* OUTPUT:FYLAG2,KEY2 * C* * C* LANGUAGE FORTRAN * C* * C* REFRENCES JOHLER ET AL (1956) : NBS CIRCULAR 573 * C********************************************************************** IMPLICIT REAL *8(A-B,D-H,O-Z),COMPLEX *16(C) REAL *8 K1, KE, MUO PI =4.D0*DATAN(1.D0) MUO=4.0D-07*PI W =2.D0*PI*FREQ K1 =W*DCBRT(E1)/V AK =AVCON*MUO*V*V/W CK2=(W/V)*CDSQRT(DCMPLX(E2,AK)) IF(AK.NE.0.D0) GO TO 41 PSIE=0.5D0*PI GO TO 42 41 PSIE =DATAN(E2/AK)-0.5D0*DATAN((E2-E1)/AK) 42 KE =DCBRT(V*ALF/(E1*E1*W*A))*DSQRT(E2*E2+AK*AK)/((E2-E1)* & (E2-E1)+AK*AK)**0.25D0 PP =0.75D0*PI-PSIE CPPI =DCMPLX(0.D0,PP) CALL TRAPS(0,0,100,0,0) CDELE=DCMPLX(KE,0.D0)*CDEXP(CPPI) CRHOC=-DIST*DCBRT(K1*A*ALF*ALF)/(A*2.D0*CDELE*CDELE) CRHO1=CRHOC*DCMPLX(0.D0,1.D0) CFUNC=1.D0-K1*K1/(CK2*CK2)+(K1/CK2)**4 CALL TRAPS(0,0,100,0,0) CY =(CDSQRT(PI*CRHO1)*CDEXP(-CRHO1))*DCMPLX(0.D0,1.D0)+ & DCMPLX(1.D0,0.D0) CADD =DCMPLX(1.D0,0.D0) DO 43 I=1,1000 KEY2=I CADD=CADD*(-4.D0*CRHO1)/DFLOAT(4*I-2) CY =CY+CADD IF(CDABS(CADD) .LT. 1.0D-7*CDABS(CY)) GO TO 44 43 CONTINUE 44 DCORR=1.D0/(K1*DIST) FZR =DREAL(CY*CFUNC)-DCORR*DCORR FZI =DIMAG(CY*CFUNC)+DCORR FYCT =(DATAN(FZI/FZR))/W FYLAG2 =FYCT*V 45 IF(FYLAG2.GT.0.D0) GO TO 46 FYLAG2=FYLAG2+V/(2.D0*FREQ) GO TO 45 46 RETURN END SUBROUTINE PHASE0 (DIST,FREQ,V,E1,FYT) C********************************************************************** C* NAME PHASE0 * C* * C* TYPE SUBROUTINE * C* * C* PURPOSE TO COMPUTE SECONDARY PHASE LAG FOR A SPECIAL CASE: * C* FREE SPACE ( Z=0 ) * C* * C* AUTHOR D.GRAY * C* MODIFIER T.MAHMOOD * C* * C* EXTERNALS NONE * C* * C* CALLING CALL PHASE0(DIST,FREQ,V,E1,FYT) * C* * C* PARAMETERS ALL EXPLAINED IN THE MAIN * C* INPUT :DIST,FREQ,V,E1 * C* OUTPUT:FYT * C* * C* LANGUAGE FORTRAN * C* * C* REFRENCES JOHLER ET AL (1956) : NBS CIRCULAR 573 * C* * C* NOTE THIS ROUTINE IS NOT USED AT PRESENT IN THE MAIN PROG. * C********************************************************************** IMPLICIT REAL *8(A-B,D-H,O-Z),COMPLEX *16(C) REAL *8 K1 PI =4.D0*DATAN(1.D0) W =2.D0*PI*FQ K1 =W*DSQRT(E1)/V DCORR=1.D0/(K1*DIST) FRR =-DCORR*DCORR+1.0 FRI =DCORR FYCT =(DATAN(FRI/FRR))/W FYT =FYCT*V 51 IF(FYT.GT.0.D0) GO TO 52 FYT=FYT+V/(2.D0*FREQ) GO TO 51 52 RETURN END FUNCTION DCBRT(X) C********************************************************************** C* NAME DCBRT(X) * C* * C* TYPE DOUBLE PRECISION FUCTION * C* * C* PURPOSE TO COMPUTE THE CUBE ROOT OF A REAL NUMBER * C* * C* AUTHOR D.GRAY 1977/78 * C* * C* PARAMETERS X=THE NUMBER TO BE COMPUTED * C* * C* LANGUAGE FORTRAN * C********************************************************************** IMPLICIT REAL *8(A-H,O-Z) IF(X.GE.0.D0) DCBRT=X**(1.D0/3.D0) IF(X.LT.0.D0) DCBRT=-(DABS(X))**(1.D0/3.D0) RETURN END