IMPLICIT REAL *8(A-B,D-H,O-Z),COMPLEX *16(C) REAL *8 K1,KE DIMENSION CTS(625),TFY(219) C*********************************************************************** C* NAME FASELAG * C* * C* TYPE MAIN * C* * C* PURPOSE TO COMPUTE THE SECONDARY PHASELAG(IN METERS) USING * C* SPHERICAL EARTH THEORY FOR VARIOUS FREQUENCIES, * C* CONDUCTIVITIES AND DIELECTRIC CONSTANTS AS ENTERED * C* * C* AUTHOR P.BRUNAVS 1977 * C* MODIFIER T. MAHMOOD OCT. 1981 * C* * C* EXTERNALS TSZERO,TSINTY,QUAD,ASYMP,TRAPS * C* NOTE:TRAPS IS A BUILT-IN ROUTINE TO OVERCOME * C* EXPONENT UNDERFLOW * C* * C* PARAMETERS * C* ICODE= A CODE TO INDICATE THE END OF DATA;GREATER THAN 0 * C* FRQ=FREQUANCY(KHZ) * C* SG =CONDUCTIVITY(MHOS/METER) * C* V =VELOCITY(METERS/SEC) * C* A =RADIUS OF THE EARTH(METER) * C* ALF=VERTICAL LAPSE RATE OF THE PERMITTIVITY OF THE ATMOSPHERE * C* D =DISTANCE(METER) * C* W =ANGULAR FREQUENCY=2*PHI*FRQ (RADIANS/SEC) * C* E1 =PERMITTIVITY OF THE AIR AT THE SURFACE OF THE EARTH (E.S.U)* C* =1.000676 (ASSUMED) * C* E2 =PERMITTIVITY (DIELECTRIC CONSTANT) OF THE EARTH (E.S.U) * C* =15 OR LESS FOR THE EARTH * C* =81 FOR THE SEA WATER * C* K1 =WAVE NUMBER OF THE ATMOSPHERE AT THE SURFACE OF THE EARTH, * C* =W/V*SQRT(E1) (RADIANS/METER) * C* CDE=THE CONDUCTIVITY AND PERMITTIVITY PARAMETER,FOR A VERTICAL * C* DIPOLE SOURCE * C* KE =THE MODULUS OR AMPLITUDE OF CDE * C* CTS=TAU S ;SOLUTION OF RICATTI EQUATION * C* IS =SUBSCRIPT ASSOCIATED WITH CTS * C* AMA=THE ABSOLUTE VALUE OF THE LAST TERM USED IN POWER SERIES * C* LFA=CODES USED TO DETERMINE WHETHER THE SERIES USED TO COMPUTE * C* CTS CONVERGE OR NOT ; * C* =1 DOES NOT CONVERGE,>0.1 * C* =2 DOES NOT CONVERVE * C* =3 DOES NOT CONVERGE AFTER FULL ITERATION, CTS*10**-7< >0.1* C* =4 CONVERGE AFTER FULL ITERATION, 32 POWER SERIES(DESCENDING OR ASCENDING) FAIL TO CONVERGE, * C* HANKEL FUNCTION IS USED * C* VK3 =THE ORDER OF HANKEL FUNCTION=1/3 OR 2/3 * C* CKNDM=COMPLEX I * C* CZH =THE QUANTITY OF THE HANKEL FUNCTION TO BE EVALUATED * C* CHNK=EVALUATED HANKEL FUNCTION * C* FYM =SECONDARY PHASELAG (METERS) * C* TFY =SECONDARY PHASELAG (METERS) * C* MR =NUMBER OF TERMS USED IN THE COMPUTATION OF THE PHASELAG * C* (RESIDUE SERIES) * C* * C* LANGUAGE FORTRAN * C* * C* REFERENCES JOHLER ET AL (1956) : NBS CIRCULAR 573 * C* HOWE (1960) : NBS JOURNAL,VOL.64B,#2 * C* WALTERS ET AL (1962) : NBS JOURNAL,VOL.66D,#1 * C*********************************************************************** DATA IR,IP,IW/5,6,6/ C C.....NMX=MAXIMUM NUMBER USED IN TS MATRIX C.....ITAUPR=1 IF PRINTOUT OF TAU'S REQUIRED C READ (IR,1001) NMX,ITAUPR IF (NMX.GT.620) NMX= 620 C C.....READ IN DATA C 101 READ (IR,1002) ICODE,FRQ,SG,E2 IF(ICODE.GE.1)GO TO 200 CC1 =DCMPLX(1.D0,0.D0) CC2 =DCMPLX(0.D0,1.D0) C1 =DCMPLX(0.5D0,-0.5D0*DSQRT(3.D0)) C2 =DCMPLX(-0.5D0,-0.5D0*DSQRT(3.D0)) PI =4.D0*DATAN(1.D0) V =2.997925D+8 E1 =1.000676D0 A =6.36739D+6 ALF =0.75D0 W =2.D0*PI*FRQ K1 =(W/V)*DSQRT(E1) SMC =(SG*PI*V*V*4.0D-7)/W KE =(((V*ALF)/(A*W*E1*E1))**4*(E2*E2+SMC*SMC)**6 & /((E2-E1)**2+SMC*SMC)**3)**(1.D0/12.D0) PSI =DATAN(E2/SMC)-0.5D0*DATAN((E2-E1)/SMC) PP =PI*0.75D0-PSI CPPI =DCMPLX(0.D0,PP) CDE =DCMPLX(KE,0.D0)*CDEXP(CPPI) FRC =(K1*A*ALF*ALF)**(1.D0/3.D0) DARG =DATAN2(DIMAG(CDE),DREAL(CDE))*180.D0/PI ER3 = 0.1D0 KP =1 KOUT =0 N =0 C C.....COMPUTATION OF TAU S(CTS) C 102 N =N+1 IS = N-1 AMT= ER3 CTS(N)= DCMPLX(0.D0,0.D0) CTSST = DCMPLX(0.D0,0.D0) CDIFT = DCMPLX(0.D0,0.D0) IF(KP.EQ.2) GO TO 105 C C.....COMPUTATION OF CTS USING ASCENDING SERIES C CALL TSZERO(IS,CDE,CTSV,AMA,LFA) CTSST=CTSV AMT =AMA GO TO (103,200,104,114,114),LFA 103 KOUT=1 104 KP =2 C C.....COMPUTATION OF TAU S(CTS) USING DESCENDING SERIES C 105 CALL TSINTY(IS,CDE,CTSV,AMA,LFA) GO TO (106,200,107,114,114),LFA 106 IF (KOUT.EQ.1) GO TO 200 107 KP=1 IF(AMA.GT.AMT) GO TO 108 KP=2 CTSST=CTSV AMT =AMA 108 CTSS =CTSST TIM =0.0001D0 TRE =0.0001D0 J =25 ITS =0 C C.....COMPUTATION OF TAU S(CTS) USING HANKEL FUNCTIONS IF BOTH ASCENDING C.....AND DECENDING SERIES DO NOT CONVERGE C 109 ITS =ITS+1 CTSI =CTSS DO 112 IK=1,3 IF(IK.EQ.2) CTSI=CTSS+DCMPLX(TRE,0.D0) IF(IK.EQ.3) CTSI=CTSS+DCMPLX(0.D0,TIM) CZH =(CDSQRT(CTSI*2)**3*DCMPLX(0.D0,-1.D0))/3.D0 CRUTS=CDSQRT(-CTSI*2) KIND =1 CKNDM=DCMPLX(0.D0,1.D0) C C..... IF(KIND.EQ.2) CKNDM=-CKNDM C C..... EVALUATION OF HANKEL FUNCTIONS C DO 110 KV=1,2 VK3=DFLOAT(KV)/3.D0 IF(N.LE.J) CALL QUAD(CKNDM,VK3,CZH,CHNK) IF(N.GE.(J+1)) CALL ASYMP(CKNDM,VK3,CZH,CHNK,IRT) IF(KV.EQ.1) CHNK1=-CHNK 110 CONTINUE CHNK2=-CHNK CZH=DCONJG(CZH) DO 111 KV=1,2 VK3=DFLOAT(KV)/3.D0 IF(N.LE.J) CALL QUAD(CKNDM,VK3,CZH,CHNK) IF(N.GE.(J+1)) CALL ASYMP(CKNDM,VK3,CZH,CHNK,IRT) IF(KV.EQ.1)CHNK1=CHNK1-DCONJG(CHNK)*C1 111 CONTINUE CZH =DCONJG(CZH) CHNK2=CHNK2+DCONJG(CHNK)*C2 CHTS =(CHNK1/CHNK2)*DCMPLX(0.5D0,0.5D0*DSQRT(3.D0))/CRUTS IF(IK.EQ.1) CDEH=CHTS IF(IK.EQ.2) CAA =CHTS-CDEH 112 CONTINUE CBB =CHTS-CDEH CC =CDE-CDEH DTRE =(DIMAG(CC/CBB)/DIMAG(CAA/CBB))*TRE DTIM =(DIMAG(CC/CAA)/DIMAG(CBB/CAA))*TIM CDTSI=DCMPLX(DTRE,DTIM) CTSS =CTSS+CDTSI IF(CDABS(CDTSI).LT.1.0D-7) GO TO 113 IF(ITS.LT.10) GO TO 109 LFA =1 113 LFA =LFA+32 CTSV =CTSS 114 CTS(N)=CTSV IF(ITAUPR.LT.1) GO TO 116 IF (N.GT.200) GO TO 116 C C.....PRINTING THE MAIN PARAMETERS IN THE COMPUTATION OF TAU S(CTS) C IF (N.EQ.1) WRITE (IW,1003) FRQ,SG,E2,E1,KE,CDE IF(N.GE.4) CDIFT=(CTS(N)-CTS(N-3)+(CTS(N-2)-CTS(N-1))*3.D0)*1.0D+6 IF(LFA.LE.5) GO TO 115 CTINC=CTS(N)-CTSST TMOD =CDABS(CTS(N)) TARG =DATAN2(DIMAG(CTS(N)),DREAL(CTS(N)))*180.D0/PI C C.....PRINTING THE VALUES OF TAU S(CTS) AND OTHERS IF TAU S(CTS) C.....ARE DETERMINED USING HANKEL FUNCTIONS C WRITE(IW,1004) IS,CTS(N),CDIFT,LFA,AMT,CTINC,TMOD,TARG,IRT GO TO 116 C C.....PRINTING THE VALUES OF TAU S(CTS) AND OTHERS IF TAU S(CTS) C.....ARE DETERMINED USING EITHER ASCENDING OR DESCENDING SERIES C 115 WRITE (IW,1004) IS,CTS(N),CDIFT,LFA 116 IF(N.LT.NMX) GO TO 102 C C.....PRINTING THE MAIN PARAMETERS IN THE COMPUTATION OF THE PHASELAG C WRITE(IW,1005) FRQ,SG,E2,E1,CDE,KE,DARG C C.....COMPUTATION OF SECONDARY PHASELAG C FYP=0.0D0 DBIG=50000.D0 DO 120 K=8,219 IF(K.LT.20) D=K*1000.D0 IF(K.GE.20) D=(K-20)*10000.D0+20000.D0 DCR =1.D0/(K1*D) FWC =DSQRT(2.D0*PI*FRC*D/A) CSUM=DCMPLX(0.D0,0.D0) FYI =0.D0 MLT =0 DO 117 M=1,NMX CT =CTS(M) CFX =CT*DCMPLX(0.D0,FRC*D/A)+DCMPLX(0.D0,0.5D0*ALF*D/ & A+0.25D0*PI) CALL TRAPS(0,0,500,0,0) CSUM=CSUM+CDEXP(CFX)/(CT*DCMPLX(2.D0,0.D0)-CC1/CDE**2) CFR =CSUM*FWC*DCMPLX(1.D0-DCR**2,DCR) FY =DATAN2(DIMAG(CFR),DREAL(CFR)) MR =M IF(DABS(FY-FYI).LT.0.000001D0) MLT=MLT+1 IF(MLT.GT.5) GO TO 118 117 FYI=FY 118 IF(FY.LT.0.D0) FY=FY+PI*2.D0 IF(D.LT.DBIG) GO TO 119 IF(FY.LT.FYP) FY=FY+PI*2.D0 IF(FY.LT.FYP) FY=FY+PI*2.D0 IF(FY.LT.FYP) FY=FY+PI*2.D0 119 FYM =FY*V/W CALL TRAPS(0,0,500,0,0) ERA =CDABS(CFR*CDEXP(DCMPLX(0.D0,K1*D)))*W*2.D0 FYB =FYM-DCR*V/W C C..... PRINTING THE DISTANCE AND SECONDARY PHASELAG C WRITE(IW,1006) D,FYB,FYM,ERA,MR FYP =FY TFY(K)=FYM 120 CONTINUE C C.....PRINTING THE SECONDARY PHASELAG(ONLY);SUITABLE FOR OTHER FORM C.....OF STORAGE(E.G. TAPE OR DISK) C WRITE(IP,1007) (TFY(I),I=8,219) GO TO 101 1001 FORMAT(5X,2I5) 1002 FORMAT(I2,F10.1,F10.5,F6.1) 1003 FORMAT(1H1,10X,F12.1,F10.5,F4.0,F10.6,3D16.8/) 1004 FORMAT(I4,2F14.8,2X,2F9.2,I4,D11.3,2F14.8,F13.6,F11.6,I3) 1005 FORMAT(1H1,10X,F12.0,F12.5,F7.1,F10.6/10X,2D17.8,2F14.8///) 1006 FORMAT(10X,F9.1,2F12.1,D17.8,I8) 1007 FORMAT(8F10.3) 200 STOP END SUBROUTINE TSINTY(IS,CDE,CTSV,AMA,LFA) C********************************************************************** C* NAME TSINTY * C* * C* TYPE SUBROUTINE * C* * C* PURPOSE TO COMPUTE TAU S(CTS) USING DESCENDING POWER SERIES * C* * C* AUTHOR P.BRUNAVS 1977 * C* MODIFIER T.MAHMOOD OCT.1981 * C* * C* EXTERNALS NONE * C* * C* CALLING CALL TSINTY(IS,CDE,CTSV,AMA,LFA) * C* * C* PARAMETERS ALL EXPLAINED IN THE MAIN * C* INPUT : IS,CDE * C* OUTPUT: CTSV,AMA,LFA * C* * C* LANGUAGE FORTRAN * C* * C* REFERENCES JOHLER 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 CA(70),TU(6) TF(Y)=((Y*Y/2.D0)**(1.D0/3.D0))*(1.D0-(7.D0/(48.D0*Y*Y)) & +(35.D0/(288.D0 *Y*Y*Y*Y))) TU(1)=0.808616516D0 TU(2)=2.57809613D0 TU(3)=3.82571528D0 TU(4)=4.89182029D0 TU(5)=5.85130097D0 TU(6)=6.73731638D0 PI =4.D0*DATAN(1.D0) N =IS+1 TAU =TF((DFLOAT(IS)*4.D0+1.D0)*PI*0.375D0) IF(N.LE.6) TAU=TU(N) CTAUC=DCMPLX(0.5D0,0.5D0*DSQRT(3.0D0))*TAU ER =TAU*1.0D-9 ER2 =TAU*1.0D-7 EMX =1.0D0 ER3 =0.1D0 CAD =-0.5D0/CTAUC CDE2 =1.D0/CDE**2 CA(1)=CTAUC CA(2)=CAD/CDE CTSV =CA(1)+CA(2) LFA =1 DO 22 I=3,50 DCI =(DFLOAT(I)-3.D0)/(DFLOAT(I)-1.D0) ME =I-1 CSUM=DCMPLX(0.D0,0.D0) DO 21 M=2,ME 21 CSUM =CSUM+CA(M)*CA(ME-M+2) CA(I)=(CSUM-CA(I-2)*CDE2*DCI)*CAD CTSV =CTSV+CA(I) AMA =CDABS(CA(I)) IF(AMA.GT.EMX)GO TO 24 IF(AMA.LT.ER )GO TO 23 22 CONTINUE IF(AMA.GE.ER3) GO TO 24 LFA=3 IF(AMA.LT.ER2) LFA=4 GO TO 24 23 LFA=5 24 RETURN END SUBROUTINE TSZERO(IS,CDE,CTSV,AMA,LFA) C********************************************************************** C* NAME TSZERO * C* * C* TYPE SUBROUTINE * C* * C* PURPOSE TO COMPUTE TAU S(CTS) USING ASCENDING POWER SERIES * C* * C* AUTHOR P.BRUNAVS 1977 * C* MODIFIER T.MAHMOOD OCT.1981 * C* * C* EXTERNALS NONE * C* * C* CALLING CALL TSZERO(IS,CDE,CTSV,AMA,LFA) * C* * C* PARAMETERS ALL EXPLAINED IN THE MAIN * C* INPUT : IS,CDE * C* OUTPUT:CTSV,AMA,LFA * C* * C* LANGUAGE FORTRAN * C* * C* REFRENCES JOHLER 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 CA(70),TU(6) TF(X)=((X*X/2.D0)**(1.D0/3.D0))*(1.D0+(5.D0/(48.D0*X*X)) & -(5.D0/(36.D0*X*X*X*X))) TU(1)=1.85575708D0 TU(2)=3.24460762D0 TU(3)=4.38167124D0 TU(4)=5.38661378D0 TU(5)=6.30526301D0 TU(6)=7.16128272D0 N =IS+1 PI =4.D0*DATAN(1.D0) TAU =TF((DFLOAT(IS)*4.D0+3.D0)*PI*0.375D0) IF(N.LE.6) TAU=TU(N) CTAUC=DCMPLX(0.5D0,0.5D0*DSQRT(3.D0))*TAU ER =TAU*1.0D-9 ER2 =TAU*1.0D-7 EMX =1.D0 ER3 =0.1D0 CDES =CDE**2 CA(1)=CTAUC CA(2)=-CDE CA(3)=DCMPLX(0.D0,0.D0) CTSV =CA(1)+CA(2)+CA(3) LFA = 1 DO 32 I=4,50 DCI=(DFLOAT(I)-3.D0)/(DFLOAT(I)-1.D0) ME =I-2 CSUM=DCMPLX(0.D0,0.D0) DO 31 M=1,ME 31 CSUM=CSUM+CA(M)*CA(ME-M+1) CA(I)=CSUM*CDES*DCI CTSV =CTSV+CA(I) AMA =CDABS(CA(I)) IF(AMA.GT.EMX)GO TO 34 IF(AMA.LT.ER )GO TO 33 32 CONTINUE IF(AMA.GE.ER3) GO TO 34 LFA=3 IF(AMA.LT.ER2) LFA=4 GO TO 34 33 LFA= 5 34 RETURN END COMPLEX FUNCTION CSX(CKNDM,VK3,CZH,T) C********************************************************************** C* NAME CSX * C* * C* TYPE COMPLEX FUNCTION (*16) * C* * C* PURPOSE TO COMPUTE THE VALUE OF THE FIRST TERM OF THE HANKEL * C* FUNCTION * C* * C* AUTHOR P.BRUNAVS 1977 * C* MODIFIER T.MAHMOOD OCT.1981 * C* * C* EXTERNALS NONE * C* * C* CALLING X=CSX(CKNDM,VK3,CZH,T) * C* * C* PARAMETERS ALL EXPLAINED IN THE MAIN AND SUBROUTINE QUAD * C* * C* LANGUAGE FORTRAN * C* * C*REFERENCES WALTERS ET AL (1962) :NBS JOURNAL,VOL.66D,#1 * C********************************************************************** IMPLICIT REAL *8(A-B,D-H,O-Z),COMPLEX *16(C) CALL TRAPS(0,0,500,0,0) CSX=CDEXP((CZH*DCMPLX(DSIN(T),0.D0)-DCMPLX(T,0.D0)* & DCMPLX(VK3,0.D0))*CKNDM) RETURN END COMPLEX FUNCTION CSU(CKNDM,VK3,CZH,VV) C********************************************************************** C* NAME CSU * C* * C* TYPE COMPLEX FUNCTION (*16) * C* * C* PURPOSE TO COMPUTE THE VALUE OF THE SECOND TERM OF THE HANKEL * C* FUNCTION * C* * C* AUTHOR P.BRUNAVS 1977 * C* MODIFIER T.MAHMOOD OCT.1981 * C* * C* EXTERNALS NONE * C* * C* CALLING CALL CSU(CKNDM,VK3,CZH,VV) * C* * C* PARAMETERS ALL EXPLAINED IN THE MAIN AND SUBROUTINE QUAD * C* * C* LANGUAGE FORTRAN * C* * C* REFRENCES WALTERS ET AL (1962) : NBS JOURNAL,VOL.66D,#1 * C********************************************************************** IMPLICIT REAL *8(A-B,D-H,O-Z),COMPLEX *16(C) PI =4.D0*DATAN(1.D0) FV =(VV-1.D0/VV)*0.5D0 PV =-PI*VK3 U =VV**VK3 CU =DCMPLX(U,0.D0) CPV=CKNDM*PV CEZ=CZH*FV CSU=DCMPLX(0.D0,0.D0) IF((DREAL(CEZ)+670.D0).LT.0.D0) GO TO 50 CALL TRAPS(0,0,500,0,0) CM =(CKNDM*CDEXP(CEZ))/VV CALL TRAPS(0,0,500,0,0) CSU=CM/CU+(CM*CDEXP(CPV))*CU 50 RETURN END SUBROUTINE QUAD(CKNDM,VK3,CZH,CHNK) C********************************************************************** C* NAME QUAD * C* * C* TYPE SUBROUTINE * C* * C* PURPOSE EVALUATION OF DEFINITE INTERGRAL OF F(X)*DX FROM * C* X=A TO X=B USING GAUSS-LEGENDRE QUADRATURE OF ORDER 48 * C* TO EVALUATE THE HANKEL FUNCTIONS OF ORDER VK3 * C* * C* AUTHOR P.BRUNAV 1977 * C* MODIFIER T.MAHMOOD OCT.1981 * C* * C* EXTERNALS CSU,CSX * C* * C* CALLING CALL QUAD(CKNDM,VK3,CZH,CHNK) * C* * C* PARAMETERS ALL EXPLAINED IN THE MAIN + * C* R=GAUSSIAN ABSCISSAS TABULATED BY KOPAL * C* W=GAUSSIAN WEIGHT TABULATED BY KOPAL * C* INPUT :CKNDM,VK3,CZH * C* OUTPUT:CHNK * C* * C* LANGUAGE FORTRAN * C* * C* REFRENCES WALTERS ET AL (1962) : NBS JOURNAL,VOL.66D,#1 * C* F.KOPAL (1955) : NUMERICAL ANALYSIS(J.WILEY & SON,N.Y) * C********************************************************************** IMPLICIT REAL *8(A-B,D-H,O-Z),COMPLEX *16(C) DIMENSION R(24),W(24) EXTERNAL CSU,CSX PI=4.D0*DATAN(1.D0) DATA N/24/ DATA R/.03238 01709 62869, .09700 46992 09463, .16122 23560 68892, & .22476 37903 94689, .28736 24873 5546 , .34875 58862 9216 , & .40868 64819 9072 , .46690 29047 5096 , .52316 09747 2223 , & .57722 47260 8397 , .62886 73967 7651 , .67787 23796 3266 , & .72403 41309 2382 , .76715 90325 1574 , .80706 62040 2944 , & .84358 82616 2439 , .87657 20202 7425 , .90587 91367 1557 , & .93138 66907 0655 , .95298 77031 6043 , .97059 15925 4625 , & .98412 45837 2283 , .99353 01722 6635 , .99877 10072 5243 / DATA W/.06473 76968 12684, .06446 61644 35950, .06392 42385 84648, & .06311 41922 86254, .06203 94231 59893, .06070 44391 65894, & .05911 48396 98396, .05727 72921 00403, .05519 95036 99984, & .05289 01894 85194, .05035 90355 53854, .04761 66584 92490, & .04467 45608 56694, .04154 50829 43465, .03824 13510 65831, & .03477 72225 64770, .03116 72278 32798, .02742 65097 08357, & .02357 07608 39324, .01961 61604 57356, .01557 93157 22944, & .01147 72345 79235, .00732 75539 01276, .00315 33460 52306/ CZ =DCMPLX(0.D0,0.D0) C C.....EVALUATION OF THE FIRST TERM OF THE HANKEL FUNCTION C XMD=0.25D0*PI SCL=0.25D0*PI XMI=0.75D0*PI CQX=CZ DO 61 K=1,N S=R(K)*SCL 61 CQX=CQX+(CSX(CKNDM,VK3,CZH,XMD+S)+CSX(CKNDM,VK3,CZH,XMD-S)+ & CSX(CKNDM,VK3,CZH,XMI+S)+CSX(CKNDM,VK3,CZH,XMI-S))* & DCMPLX(W(K),0.D0) CQX=CQX*SCL C C.....EVALUATION OF THE SECOND TERM OF THE HANKEL FUNCTION C XMD=0.5D0 SCL=0.5D0 CQU=CZ DO 62 K=1,N S=R(K)*SCL 62 CQU=CQU+(CSU(CKNDM,VK3,CZH,XMD+S)+CSU(CKNDM,VK3,CZH,XMD-S))* & DCMPLX(W(K),0.D0) CQU =CQU*SCL CHNK=(CQX-CQU )/DCMPLX(PI,0.D0) RETURN END SUBROUTINE ASYMP(CKNDM,VK3,CZH,CHNK,IRT) C********************************************************************** C* NAME ASYMP * C* * C* TYPE SUBROUTINE * C* * C* PURPOSE EVALUATION OF HANKEL FUNCTION USING ASSYMTOTIC * C* EXPANSION * C* * C* AUTHOR P.BRUNAVS 1977 * C* MODIFIER T.MAHMOOD OCT.1981 * C* * C* EXTERNAL NONE * C* * C* CALLING CALL ASYMP(CKNDM,VK3,CZH,CHNK,IRT) * C* * C* PARAMETERS ALL EXPLAINED IN THE MAIN * C* INPUT :CKNDM,VK3,CZH * C* OUTPUT:CHNK,IRT * C* * C* LANGUAGE FORTRAN * C* * C* REFRENCES UNKNOWN * C********************************************************************** IMPLICIT REAL *8(A-B,D-H,O-Z),COMPLEX *16(C) PI =4.D0*DATAN(1.D0) CEPX =(CZH-DCMPLX((0.5D0*VK3+0.25D0)*PI,0.D0))*CKNDM CM =CDSQRT(DCMPLX(2.D0/PI,0.D0)/CZH) CZ8 =DCMPLX(0.125D0,0.D0)/CZH V42 =VK3**2*4.D0 CSUMA=DCMPLX(1.D0,0.D0) CT1A =DCMPLX(1.D0,0.D0) DO 71 L=1,50 FL =DFLOAT(L) VNA =(V42-(FL*2.D0-1.D0)**2)/FL CT2A=(CZ8*CKNDM**L)*VNA IF(CDABS(CT2A).GT.0.99D0) GO TO 72 CT2A =CT1A*CT2A CSUMA=CSUMA+CT2A DCRIT=CDABS(CT2A)/CDABS(CSUMA) IF(DCRIT.LT.1.0D-11) GO TO 73 CT1A=CT2A KEA =L 71 CONTINUE 72 IRT =3 GO TO 74 73 CALL TRAPS(0,0,500,0,0) CHNK=CSUMA*CM*CDEXP(CEPX) IRT =5 74 RETURN END