C PDSUBS: A collection of source-field subroutines for pulse C basis functions and Dirac delta testing functions, for C 2D problems involving waves incident at an oblique C angle (possibly) with respect to the cylinder axis C C See Appendix B of the text "Computational Methods for C Electromagnetics, IEEE Press, 1998 for a complete C description of the appropriate equations. C C Compiled by A. F. Peterson, December 1997 C C C---------------------------------------------------------- C LIST OF SUBROUTINES AND FUNCTIONS: C C HZPKT - HZ-FIELD PRODUCED BY KT-CURRENT C HZPJT - HZ-FIELD PRODUCED BY JT-CURRENT C HZPKZ - HZ-FIELD PRODUCED BY KZ-CURRENT C EZPJZ - EZ-FIELD PRODUCED BY JZ-CURRENT C EZPKT - EZ-FIELD PRODUCED BY KT-CURRENT C FIM1 - INTEGRAND USED BY 'EZPKT' C FIM2 - INTEGRAND USED BY 'EZPKT' C EZPJT - EZ-FIELD PRODUCED BY JT-CURRENT C HTPJT - HT-FIELD PRODUCED BY JT-CURRENT C HTPJZ - HT-FIELD PRODUCED BY JZ-CURRENT C FIM3 - INTEGRAND USED BY 'HTPJZ' C FIM4 - INTEGRAND USED BY 'HTPJZ' C HTPKT - HT-FIELD PRODUCED BY KT-CURRENT C HTPKZ - HT-FIELD PRODUCED BY KZ-CURRENT C ETPJZ - ET-FIELD PRODUCED BY JZ-CURRENT C FIM5 - INTEGRAND USED BY 'ETPJZ' C ETPJT - ET-FIELD PRODUCED BY JT-CURRENT C ETPKT - ET-FIELD PRODUCED BY KT-CURRENT C ETPKZ - ET-FIELD PRODUCED BY KZ-CURRENT C H02 - ZERO-ORDER HANKEL FUNCTION C BJ1 - FIRST-ORDER BESSEL FUNCTION C BY1 - FIRST-ORDER NEUMANN FUNCTION C K0 - ZERO-ORDER MODIFIED BESSEL FUNCTION 2ND KIND C K1 - FIRST-ORDER MODIFIED BESSEL FUNCTION 2ND KIND C ROM1D - ROMBERG NUMERICAL INTEGRATION C C (routines for HZPJZ and EZPKZ are not included since those C values are always zero) C---------------------------------------------------------- C C Format for calling statements: C C COMPLEX FUNCTION HZPKT(X,Y,W,PHI,KZ) C COMPLEX FUNCTION HTPKT(X,Y,W,PHI,PSI,KZ) C COMPLEX FUNCTION HZPJT(X,Y,W,PHI,KZ) C COMPLEX FUNCTION EZPKT(X,Y,W,PHI,KZ) C COMPLEX FUNCTION HTPJZ(X,Y,W,PHI,PSI,KZ) C COMPLEX FUNCTION HTPJT(X,Y,W,PHI,PSI,KZ) C COMPLEX FUNCTION EZPJZ(X,Y,W,PHI,KZ) C COMPLEX FUNCTION ETPJT(X,Y,W,PHI,PSI,KZ) C COMPLEX FUNCTION ETPJZ(X,Y,W,PHI,PSI,KZ) C COMPLEX FUNCTION EZPJT(X,Y,W,PHI,KZ) C COMPLEX FUNCTION ETPKT(X,Y,W,PHI,PSI,KZ) C COMPLEX FUNCTION HZPKZ(X,Y,W,PHI,KZ) C COMPLEX FUNCTION ETPKZ(X,Y,W,PHI,PSI,KZ) C COMPLEX FUNCTION HTPKZ(X,Y,W,PHI,PSI,KZ) C C C Input data: (all real-valued) C C X,Y distance between center of source cell and C observation point in x and y directions C in units of free-space wavelength (the total C free-space wavelength, NOT the wavelength in C the x-y plane) C C eg: X = x(observer) - x(source) C C W width of source strip in free-space wavelengths C C PHI polar angle in radians giving the outward C normal vector of the source strip-cell, C where the outward normal vector is defined C C U(X)*COS(PHI)+U(Y)*SIN(PHI), C C and where U(X) is a unit vector in the x- C direction and U(Y) is a unit vector in the C y-direction. C C PSI polar angle in radians giving the outward C normal vector of the observer strip-cell, C where required, defined in a similar manner C as PHI C C KZ phase constant along the cylinder axis defined C so that K**2 = KT**2 + KZ**2, where C K = 2 * pi (radians/wavelength) C and the intrinsic impedance ETA = 376.73 C C Note: For normally-incident waves, use KZ = 0.0 C C Output: Complex-valued field component in units of volts or C amperes per wavelength C C C---------------------------------------------------------- C COMPLEX FUNCTION HZPKT(X,Y,W,PHI,KZ) C C HZPKT: RETURN THE HZ-FIELD PRODUCED AT (X,Y) BY AN INFINITE C STRIP OF CONSTANT UNITY-AMPLITUDE KT-CURRENT DENSITY. C THE STRIP IS OF LENGTH 'W' AND IS CENTERED AT THE C ORIGIN, SO THAT ITS OUTWARD NORMAL VECTOR MAKES A C POLAR ANGLE 'PHI' RELATIVE TO THE X-AXIS. 'KZ' IS C THE LONGITUDINAL WAVENUMBER. C C A. F. PETERSON FEB 26 1986 C C USES EXTERNAL ROUTINES 'H02' AND 'K0' TO PROVIDE THE ZERO- C ORDER HANKEL AND MODIFIED BESSEL FUNCTIONS OF THE SECOND KIND. C COMPLEX H02 REAL K,KZ,KT,K0 DATA K/6.283185308/,TWOOPI/.636619772/,ETA/376.73/ C WO2=0.5*W X0=WO2*SIN(PHI) X1=-X0 Y0=-WO2*COS(PHI) Y1=-Y0 KT=SQRT(ABS(K*K-KZ*KZ)) C C R0 IS DISTANCE FROM THE TAIL END OF THE SOURCE TO OBS PT C R1 IS DISTANCE FROM THE HEAD END OF THE SOURCE TO OBS PT C R0=SQRT((X-X0)**2+(Y-Y0)**2)*KT R1=SQRT((X-X1)**2+(Y-Y1)**2)*KT C IF(K.GT.KZ) THEN HZPKT=CMPLX(0.,-.25*KZ/(K*ETA))*(H02(R0)-H02(R1)) ELSE HZPKT=CMPLX(.25*KZ*TWOOPI/(K*ETA),0.)*(K0(R0)-K0(R1)) ENDIF RETURN END C C ---------------------------------------- C COMPLEX FUNCTION HTPKT(X,Y,W,PHI,PSI,KZ) C C HTPKT: RETURN THE HT-FIELD AT (X,Y) PRODUCED BY AN INFINITE STRIP C OF KT-CURRENT-DENSITY. THE STRIP IS OF CROSS-SECTIONAL C LENGTH 'W' AND IS CENTERED AT THE ORIGIN, SO THAT ITS C OUTWARD NORMAL VECTOR MAKES A POLAR ANGLE 'PHI' WITH THE C X-AXIS. THE TRANSVERSE H-FIELD IS THAT PARALLEL TO A C SIMILAR STRIP WHOSE NORMAL VECTOR MAKES A POLAR ANGLE C 'PSI' RELATIVE TO THE X-AXIS. 'KZ' IS THE LONGITUDINAL C WAVENUMBER. C COMPLEX H02 REAL K,KT,KZ,K0,K1 C DATA K/6.283185308/,TWOOPI/.636619772/,ETA/376.73/ WO2=0.5*W KT=SQRT(ABS(K*K-KZ*KZ)) R=SQRT(X*X+Y*Y) WOT=W*0.1 IF(R.LT.WOT) THEN C C SOURCE AND OBSERVATION CELLS COINCIDE C R1=TWOOPI*ALOG(KT*W*0.16380498) R2=KT*WO2 IF(K.GT.KZ) THEN HTPKT=(K*K*W*CMPLX(1.,-R1)-2.*KT*CMPLX(BJ1(R2), + -BY1(R2)))/(-4.*K*ETA) ELSE HTPKT=(K*K*W*CMPLX(0.,-R1)-2.*KT*CMPLX(0.,TWOOPI*K1(R2)) + )/(-4.*K*ETA) ENDIF ELSE C C SOURCE AND OBSERVATION CELLS DIFFER C SPH=SIN(PHI) CPH=COS(PHI) SPS=SIN(PSI) CPS=COS(PSI) DX1=X-WO2*SPH DY1=Y+WO2*CPH DX2=X+WO2*SPH DY2=Y-WO2*CPH R1=SQRT(DX1*DX1+DY1*DY1) R2=SQRT(DX2*DX2+DY2*DY2) IF(K.GT.KZ) THEN HTPKT=(COS(PSI-PHI)*K*K*W*H02(KT*R)+ + KT*((-SPS*DX1+CPS*DY1)/R1*CMPLX(-BJ1(KT*R1), + BY1(KT*R1))-(-SPS*DX2+CPS*DY2)/R2* + CMPLX(-BJ1(KT*R2),BY1(KT*R2))))/(-4.*K*ETA) ELSE HTPKT=(COS(PSI-PHI)*K*K*W*CMPLX(0.,TWOOPI*K0(KT*R)) + +KT*((-SPS*DX1+CPS*DY1)/R1*CMPLX(0.,-TWOOPI* + K1(KT*R1))-(-SPS*DX2+CPS*DY2)/R2* + CMPLX(0.,-TWOOPI*K1(KT*R2))))/(-4.*K*ETA) ENDIF ENDIF RETURN END C C ---------------------------------------- C COMPLEX FUNCTION HZPJT(X,Y,W,PHI,KZ) C C HZPJT: RETURN THE HZ-FIELD PRODUCED AT (X,Y) BY AN INFINITE C STRIP OF CONSTANT UNITY-AMPLITUDE JT-CURRENT DENSITY. C THE STRIP IS OF LENGTH 'W' AND IS CENTERED AT THE C ORIGIN, SO THAT ITS OUTWARD NORMAL VECTOR MAKES A C POLAR ANGLE 'PHI' RELATIVE TO THE X-AXIS. 'KZ' IS C THE LONGITUDINAL WAVENUMBER. C C USE DUALITY AND ROUTINE 'EZPKT' C C FOR THE CASE X=Y=0, RETURN THE VALUE OF HZ AN INFINITESIMAL C DISTANCE OUTSIDE THE SOURCE CELL C COMPLEX EZPKT REAL KZ HZPJT=-EZPKT(X,Y,W,PHI,KZ) RETURN END C C ------------------------------------------------- C C ZERO-ORDER HANKEL FUNCTION (2ND KIND) FROM ABR & STEGUN PP 369 C COMPLEX FUNCTION H02(A) IF(A.GT.3.000) GOTO 5 BJ=(A/3.00)**2 BJ=1.0+BJ*(-2.2499997+BJ*(1.2656208+BJ*(-.3163866+BJ* & (.0444479+BJ*(-.0039444+BJ*.00021))))) BY=(A/3.00)**2 BY=2.0/3.1415926*ALOG(A/2.)*BJ+.36746691+BY*(.60559366+BY & *(-.74350384+BY*(.25300117+BY*(-.04261214+BY*(.00427916- & BY*.00024846))))) GOTO 10 5 BJ=3.00/A F0=.79788456+BJ*(-.00000077+BJ*(-.00552740+BJ*(-.00009512 & +BJ*(.00137237+BJ*(-.00072805+BJ*.00014476))))) T0=A-.78539816+BJ*(-.04166397+BJ*(-.00003954+BJ*(.00262573 & +BJ*(-.00054125+BJ*(-.00029333+BJ*.00013558))))) BY=SQRT(A) BJ=F0*COS(T0)/BY BY=F0*SIN(T0)/BY 10 H02=CMPLX(BJ,-BY) RETURN END C C ------------------------------------------------- C FUNCTION BJ1(X) C C BESSEL FUNCTION OF ORDER 1 FOR REAL ARGUMENT C C POLYNOMIAL APPROX FROM P. 370 ABRAMOWITZ AND STEGUN C C A. F. PETERSON UIUC ECE DEPT. FEB 26, 1986 C DATA A1/.00001109/, A2/-.00031761/, A3/.00443319/ DATA A4/-.03954289/, A5/.21093573/, A6/-.56249985/ DATA A7/.5/,C1/-.00020033/ DATA C2/.00113653/, C3/-.00249511/, C4/.00017105/ DATA C5/.01659667/, C6/.00000156/, C7/.79788456/ DATA D1/-.00029166/, D2/.00079824/, D3/.00074348/ DATA D4/-.00637879/, D5/.0000565/, D6/.12499612/ DATA D7/-2.35619449/, PI/3.141592654/ IF(X) 4,1,1 1 IF(X-3) 2,2,3 2 Y=(X/3.)**2 BJ1=X*(A7+Y*(Y*(Y*(Y*(Y*(A1*Y+A2)+A3)+A4)+A5)+A6)) RETURN 3 Y=3./X F1=C7+Y*(Y*(Y*(Y*(Y*(C1*Y+C2)+C3)+C4)+C5)+C6) T1=X+D7+Y*(Y*(Y*(Y*(Y*(D1*Y+D2)+D3)+D4)+D5)+D6) V=1./SQRT(X) BJ1=V*F1*COS(T1) RETURN 4 WRITE(*,*) 'ERROR IN ARGUMENT OF BJ1 ',X BJ1=0. RETURN END C C --------------------------------------------------- C FUNCTION BY1(X) C C NUEMANN FUNCTION OF ORDER 1 FOR REAL ARGUMENT C C POLYNOMIAL APPROXIMATION FROM ABRAMOWITZ AND STEGUN, PP. 370 C C A. F. PETERSON UIUC ECE DEPT. FEB 26 1986 C DATA A1/.00001109/, A2/-.00031761/, A3/.00443319/ DATA A4/-.03954289/, A5/.21093573/, A6/-.56249985/ DATA A7/.5/ DATA B1/.0027873/, B2/-.0400976/ DATA B3/.3123951/, B4/-1.3164827/, B5/2.1682709/ DATA B6/.2212091/, B7/-.6366198/, C1/-.00020033/ DATA C2/.00113653/, C3/-.00249511/, C4/.00017105/ DATA C5/.01659667/, C6/.00000156/, C7/.79788456/ DATA D1/-.00029166/, D2/.00079824/, D3/.00074348/ DATA D4/-.00637879/, D5/.0000565/, D6/.12499612/ DATA D7/-2.35619449/, PI/3.141592654/ IF(X) 4,4,1 1 IF(X-3) 2,2,3 2 Y=(X/3.)**2 BJ=X*(A7+Y*(Y*(Y*(Y*(Y*(A1*Y+A2)+A3)+A4)+A5)+A6)) BY1=(1./X)*((2./PI)*X*ALOG(X/2.)*BJ+B7+Y*(Y*(Y*(Y*(Y*( +B1*Y+B2)+B3)+B4)+B5)+B6)) RETURN 3 Y=3./X F1=C7+Y*(Y*(Y*(Y*(Y*(C1*Y+C2)+C3)+C4)+C5)+C6) T1=X+D7+Y*(Y*(Y*(Y*(Y*(D1*Y+D2)+D3)+D4)+D5)+D6) V=1./SQRT(X) BY1=V*F1*SIN(T1) RETURN 4 WRITE(*,*) 'ERROR IN ARGUMENT OF BY1 ',X BY1=0. RETURN END C C -------------------------------------------- C REAL FUNCTION K0(X) C C K0: RETURN ZERO-ORDER MODIFIED BESSEL FUNCTION OF SECOND C KIND FOR REAL ARGUMENT -- BASED UPON POLYNOMIAL C APPROXIMATION FROM ABRAMOWITZ AND STEGUN, PP. 379 C C A. F. PETERSON UIUC ECE DEPT FEB 25 1986 C DATA A1/.0000074/,A2/.0001075/,A3/.00262698/ DATA A4/.0348859/,A5/.23069756/,A6/.4227842/ DATA A7/-.57721566/,B1/.00053208/,B2/-.0025154/ DATA B3/.00587872/,B4/-.01062446/,B5/.02189568/ DATA B6/-.07832358/,B7/1.25331414/,C1/.0045813/ DATA C2/.0360768/,C3/.2659732/,C4/1.2067492/ DATA C5/3.0899424/,C6/3.5156229/,C7/1./ IF(X) 4,4,1 1 IF(X-2.) 2,2,3 2 Y=0.25*X*X T=Y*0.284444444 BI=T*(T*(T*(T*(T*(C1*T+C2)+C3)+C4)+C5)+C6)+C7 K0=-ALOG(0.5*X)*BI+A7+Y*(Y*(Y*(Y*(Y*(A1*Y+A2)+A3)+A4)+A5)+A6) RETURN 3 Y=2./X V=1./SQRT(X) K0=V*EXP(-X)*(Y*(Y*(Y*(Y*(Y*(B1*Y+B2)+B3)+B4)+B5)+B6)+B7) RETURN 4 WRITE(*,*) 'NEGATIVE OR ZERO ARGUMENT IN K0 ',X K0=0. RETURN END C C ----------------------------------------------------- C REAL FUNCTION K1(X) C C K1: RETURN FIRST-ORDER MODIFIED BESSEL FUNCTION OF SECOND C KIND FOR REAL ARGUMENT -- BASED UPON THE POLYNOMIAL C APPROXIMATION FROM ABRAMOWITZ AND STEGUN, PP. 379 C C A. F. PETERSON UIUC ECE DEPT FEB 25 1986 C DATA A1/-.00004686/,A2/-.00110404/,A3/-.01919402/ DATA A4/-.18156897/,A5/-.67278579/,A6/.15443144/ DATA A7/1./,B1/-.00068245/,B2/.00325614/ DATA B3/-.00780353/,B4/.01504268/,B5/-.03655620/ DATA B6/.23498619/,B7/1.25331414/,C1/.00032411/ DATA C2/.00301532/,C3/.02658733/,C4/.15084934/ DATA C5/.51498869/,C6/.87890594/,C7/.5/ IF(X) 4,4,1 1 IF(X-2.) 2,2,3 2 Y=.25*X*X T=Y*0.284444444 BI=X*(T*(T*(T*(T*(T*(T*C1+C2)+C3)+C4)+C5)+C6)+C7) K1=ALOG(0.5*X)*BI+(Y*(Y*(Y*(Y*(Y*(A1*Y+A2)+A3)+A4)+A5)+A6)+A7)/X RETURN 3 Y=2./X K1=EXP(-X)*(Y*(Y*(Y*(Y*(Y*(Y*B1+B2)+B3)+B4)+B5)+B6)+B7)/SQRT(X) RETURN 4 WRITE(*,*) 'NEGATIVE OR ZERO ARGUMENT IN K1 ',X K1=0. RETURN END C C----------------------------------------------------------------- C COMPLEX FUNCTION EZPKT(X,Y,W,PHI,KZ) C C EZPKT: RETURN THE EZ-FIELD AT (X,Y) PRODUCED BY AN INFINITE C STRIP OF CONSTANT UNITY-AMPLITUDE KT-CURRENT DENSITY. C THE STRIP IS OF LENGTH 'W' AND IS CENTERED AT THE C ORIGIN, SO THAT ITS OUTWARD NORMAL VECTOR MAKES A C POLAR ANGLE 'PHI' RELATIVE TO THE X-AXIS. 'KZ' IS C THE LONGITUDINAL WAVENUMBER. C C A. F. PETERSON FEB 25 1986 C C EXTERNAL FUNCTIONS 'BJ1', 'BY1', AND 'K1' REQUIRED FOR THE FIRST C ORDER BESSEL AND MODIFIED BESSEL FUNCTIONS OF REAL ARGUMENTS. C 'ROM1D' PERFORMS A NUMERICAL INTEGRATION OVER 'BY1' THROUGH C THE FUNCTION 'FIM1', AND OVER 'K1' THROUGH THE FUNCTION 'FIM2' C C REMARKS: COMMON BLOCK 'INTE1' IS USED TO PASS DATA TO C SUBROUTINES 'FIM1'AND 'FIM2' FOR NUMERICAL INTEGRATION C REAL K,KZ,KT,K1,RO(15) COMMON /INTE1/ U,V,CP,SP,KT EXTERNAL FIM1,FIM2 C DATA K/6.283185308/,TWOOPI/.636619772/ R=SQRT(X*X+Y*Y) wot=w*.1 WO2=W*.5 IF(R.LT.wot) THEN C C SOURCE AND OBSERVATION POINTS COINCIDE C C FOR THE SELF-TERM, THIS RETURNS THE VALUE OF EZ C AN INFINITESIMAL DISTANCE OUTSIDE THE SOURCE C CELL, WHERE 'OUTSIDE' IS DEFINED BY THE NORMAL VECTOR. C EZPKT=(0.5,0.) ELSE C C SOURCE AND OBSERVATION POINTS DIFFER C KT=SQRT(ABS(K*K-KZ*KZ)) U=X V=Y CP=COS(PHI) SP=SIN(PHI) IF(K.GT.KZ) THEN EZPKT=CMPLX(BJ1(KT*R)*W*(X*CP+Y*SP)/R, + -ROM1D(FIM1,-WO2,WO2,.01,1.E-5,15,2,RO,IER)) + *CMPLX(0.,-0.25*KT) ELSE c if(kt*r.gt.k) then c c if kt*r is large, the modified bessel function k1 becomes very c small since it goes as exp(-kt*r) --- this means many of the c matrix elements are insignificant and can be neglected; we c certainly don't need to waste time numerically integrating c something that small --- use single-point integration c ezpkt=cmplx(0.25*kt*twoopi*w* + k1(kt*r)*(x*cp+y*sp)/r,0.) else EZPKT=CMPLX(0.25*KT*TWOOPI*ROM1D(FIM2, + -WO2,WO2,.01,1.E-5,15,2,RO,IER),0.) endif ENDIF ENDIF RETURN END C C ---------------------------- C FUNCTION FIM1(X) C C SEE DOCUMENTATION FOR FUNCTION 'EZPKT' C COMMON /INTE1/ U,V,CP,SP,RKT DU=U+X*SP DV=V-X*CP R=SQRT(DU*DU+DV*DV) FIM1=BY1(RKT*R)*(DU*CP+DV*SP)/R RETURN END C C ---------------------------- C FUNCTION FIM2(X) C C SEE DOCUMENTATION FOR FUNCTION 'EZPKT' C REAL KT,K1 COMMON /INTE1/ U,V,CP,SP,KT DU=U+X*SP DV=V-X*CP R=SQRT(DU*DU+DV*DV) FIM2=K1(KT*R)*(DU*CP+DV*SP)/R RETURN END C C----------------------------------------------------------------- C COMPLEX FUNCTION HTPJZ(X,Y,W,PHI,PSI,KZ) C C HTPJZ: RETURN THE HT-FIELD AT (X,Y) PRODUCED BY AN INFINITE C STRIP OF JZ-CURRENT DENSITY. THE STRIP IS OF CROSS- C SECTIONAL WIDTH 'W' AND IS CENTERED AT THE ORIGIN, C SO THAT ITS OUTWARD NORMAL VECTOR MAKES A POLAR ANGLE C 'PHI' WITH RESPECT TO THE X-AXIS. THE TRANSVERSE C COMPONENT OF THE H-FIELD IS THAT PARALLEL TO A C SIMILAR STRIP WITH OUTWARD NORMAL VECTOR DEFINED BY C THE POLAR ANGLE 'PSI' RELATIVE TO THE X-AXIS. C 'KZ' IS THE LONGITUDINAL WAVENUMBER. C C A. F. PETERSON MARCH 6, 1986 C C EXTERNAL FUNCTIONS 'BJ1', 'BY1', AND 'K1' REQUIRED FOR FIRST- C ORDER BESSEL, NEUMANN, AND MODIFIED BESSEL FUNCTIONS OF REAL C ARGUMENTS. 'ROM1D' PERFORMS A NUMERICAL INTEGRATION OVER C 'BY1' THROUGH THE FUNCTION 'FIM3', AND OVER 'K1' THROUGH THE C FUNCTION 'FIM4' C C COMMON BLOCK 'INTE2' IS USED TO PASS DATA TO SUBROUTINES C 'FIM3' AND 'FIM4' FOR NUMERICAL INTEGRATION C REAL K,KZ,KT,RO(15) COMMON /INTE2/ U,V,CPH,SPH,CPS,SPS,KT EXTERNAL FIM3,FIM4 C DATA K/6.283185308/,TWOOPI/.636619772/ R=SQRT(X*X+Y*Y) WOT=W*0.1 IF(R.LT.WOT) THEN C C SOURCE AND OBSERVATION CELLS COINCIDE C C FOR THIS CASE (X=Y=0) RETURN THE HT-FIELD JUST OUTSIDE C THE SOURCE (WHERE 'OUTSIDE' IS DEFINED BY THE NORMAL) C HTPJZ=(0.5,0.) ELSE C C SOURCE AND OBSERVATION CELLS DIFFER C WO2=W*0.5 KT=SQRT(ABS(K*K-KZ*KZ)) U=X V=Y CPH=COS(PHI) SPH=SIN(PHI) CPS=COS(PSI) SPS=SIN(PSI) IF(K.GT.KZ) THEN HTPJZ=CMPLX(BJ1(KT*R)*W*(X*CPS+Y*SPS)/R, + -ROM1D(FIM3,-WO2,WO2,.01,1.E-5,15,2,RO,IER)) + *CMPLX(0.,-0.25*KT) ELSE HTPJZ=CMPLX(0.25*KT*TWOOPI*ROM1D(FIM4,-WO2, + WO2,.01,1.E-5,15,2,RO,IER),0.) ENDIF ENDIF RETURN END C C ------------------- C FUNCTION FIM3(S) C C SEE DOCUMENTATION FOR 'HTPJZ' C REAL KT COMMON /INTE2/ U,V,CPH,SPH,CPS,SPS,KT DU=U+S*SPH DV=V-S*CPH R=SQRT(DU*DU+DV*DV) FIM3=BY1(KT*R)*(DU*CPS+DV*SPS)/R RETURN END C C ------------------- C FUNCTION FIM4(S) C C SEE DOCUMENTATION FOR 'HTPJZ' C REAL KT,K1 COMMON /INTE2/ U,V,CPH,SPH,CPS,SPS,KT DU=U+S*SPH DV=V-S*CPH R=SQRT(DU*DU+DV*DV) FIM4=K1(KT*R)*(DU*CPS+DV*SPS)/R RETURN END C C ---------------------------------------- C COMPLEX FUNCTION HTPJT(X,Y,W,PHI,PSI,KZ) C C HTPJT: RETURN HT-FIELD AT (X,Y) PRODUCED BY AN INFINITE C STRIP OF CONSTANT UNITY-AMPLITUDE JT-CURRENT C DENSITY. THE STRIP IS OF LENGTH 'W' AND IS C CENTERED AT THE ORIGIN, SO THAT ITS OUTWARD C NORMAL VECTOR MAKES A POLAR ANGLE 'PHI' WITH C RESPECT TO THE X-AXIS. THE COMPONENT OF THE C TRANSVERSE HT-FIELD IS THAT PARALLEL TO A SIMILAR C STRIP WITH AN OUTWARD NORMAL VECTOR THAT MAKES C A POLAR ANGLE 'PSI' RELATIVE TO THE X-AXIS. C 'KZ' IS THE LONGITUDINAL WAVENUMBER. C C A. F. PETERSON UIUC ECE DEPT FEB 27 1986 C C USES EXTERNAL FUNCTIONS 'H02' AND 'K0' TO PROVIDE THE C ZERO-ORDER HANKEL AND MODIFIED BESSEL FUNCTIONS OF THE C SECOND KIND C COMPLEX H02 REAL K,KZ,KT,K0 DATA K/6.283185308/,TWOOPI/.636619772/ C R=SQRT(X*X+Y*Y) WOT=W*0.1 IF(R.LT.WOT) THEN C C SOURCE AND OBSERVATION PTS COINCIDE C HTPJT=(0.,0.) C ELSE C C SOURCE AND OBSERVATION PTS DIFFER C KT=SQRT(ABS(K*K-KZ*KZ)) R1=0.25*W*KZ*SIN(PSI-PHI) IF(K.GT.KZ) THEN HTPJT=R1*H02(KT*R) ELSE HTPJT=CMPLX(0.,TWOOPI*R1)*K0(KT*R) ENDIF ENDIF RETURN END c c ------------------------------------------------------- c function rom1d(f,a,b,rerr,aerr,maxp,minp,r,ier) c c returns estimate of the integral of 'f' over the interval c (a,b) -- Romberg integration c c 'f' is an external function of the form 'f(x)' ('f' must c be declared 'external' in the calling program) c 'rerr' is the desired relative accuracy (rerr must be between c 0.00001 and 1.0) c 'aerr' is the absolute error desired -- if the absolute value of c the estimate falls below 'aerr', the algorithm terminates c (useful termination if the integrand converges to zero c and the relative accuracy does not stabilize) c 'maxp' specifies the maximum number of function evaluations c to be (1+2**(maxp-1)) c 'minp' specifies the minimum number of function evaluations c to be (1+2**(minp-1)) -- this may be necessary if the c integrand oscillates several times within the interval. c 'r' is an array of length 'maxp' (workspace provided by c main routine). On return, contains the last row of the c Romberg extrapolation array. c 'ier' on return is set to zero if accuracy met c on return is 1 if rerr is out of range c on return is 2 if maxp is less than 1 or minp .gt. maxp c on return is 3 if accuracy is not met c c author: A. F. Peterson UIUC ECE Dept c c last date revised: Jan 7, 1987 c real r(maxp) if((maxp.lt.1) .or. (minp.gt.maxp)) then ier=2 rom1d=0. go to 30 endif if((rerr.lt.0.00001).or.(rerr.gt.1.)) then ier=1 rom1d=0. go to 30 endif ier=0 sum=0.5*(f(a)+f(b)) r(1)=sum*(b-a) n=1 c c refine estimate of integral using trapezoid rule c do 20 k=2,maxp rlast=r(1) del=(b-a)/float(n) x=a-0.5*del do 10 i=1,n 10 sum=sum+f(x+del*float(i)) n=n*2 r(k)=0.5*del*sum c c update Romberg extrapolation array c w=4. do 15 j=2,k r(k+1-j)=(w*r(k+2-j)-r(k+1-j))/(w-1) w=4.*w 15 continue c rom1d=r(1) if(k.lt.minp) go to 20 c c check integral estimate: c c if change in estimate is within desired range, or if the c estimate falls below the absolute error desired, return c if(abs(rom1d).lt.aerr) return if(abs((rom1d-rlast)/rom1d).lt.rerr) return 20 continue ier=3 c c error c 30 write(*,*) 'error in rom1d -- ier = ',ier return end c c ------------------------------------------------------- c COMPLEX FUNCTION EZPJZ(X,Y,W,PHI,KZ) C C EZPJZ: RETURNS THE EZ-FIELD AT (X,Y) PRODUCED BY AN INFINITE C STRIP OF CONSTANT UNITY-AMPLITUDE JZ-CURRENT-DENSITY. C THE STRIP IS OF CROSS-SECTIONAL LENGTH 'W' AND IS C CENTERED AT THE ORIGIN, ORIENTED SO THAT ITS OUTWARD C NORMAL VECTOR MAKES A POLAR ANGLE 'PHI' WITH THE C X-AXIS. 'KZ' IS THE LONGITUDINAL WAVENUMBER. C C A. F. PETERSON FEB 25 1986 C C EXTERNAL ROUTINES 'H02' AND 'K0' REQUIRED FOR ZERO-ORDER C HANKEL AND MODIFIED BESSEL FUNCTIONS OF REAL ARGUMENT C C REMARKS: DUE TO THE SIMPLE 'ONE-POINT' APPROXIMATION USED FOR C THE INTEGRATION, 'PHI' IS NOT USED. IT IS INCORPORATED C HERE TO FACILITATE LATER IMPROVEMENTS. C COMPLEX H02 REAL K,KZ,KT,K0,KT2 DATA K/6.283185308/,TWOOPI/0.636619772/,ETA/376.73/ C R=SQRT(X*X+Y*Y) KT2=K*K-KZ*KZ KT=SQRT(ABS(KT2)) wot=w*0.1 R1=-W*KT2*.25*ETA/K IF(R.LT.wot) THEN C C SOURCE AND OBSERVATION POINTS COINCIDE C R=TWOOPI*ALOG(KT*W*0.16380498) IF(K.LT.KZ) THEN EZPJZ=CMPLX(0.,-R)*R1 ELSE EZPJZ=CMPLX(1.,-R)*R1 ENDIF ELSE C C SOURCE AND OBSERVATION POINTS DIFFER C IF(K.GT.KZ) THEN EZPJZ=H02(R*KT)*R1 ELSE EZPJZ=CMPLX(0.,TWOOPI*R1*K0(R*KT)) ENDIF ENDIF RETURN END C C ---------------------------------------- C COMPLEX FUNCTION ETPJT(X,Y,W,PHI,PSI,KZ) C C ETPJT: RETURN THE ET-FIELD AT (X,Y) PRODUCED BY AN INFINITE STRIP C OF JT-CURRENT-DENSITY. THE STRIP IS OF CROSS-SECTIONAL C LENGTH 'W' AND IS CENTERED AT THE ORIGIN, SO THAT ITS C OUTWARD NORMAL VECTOR MAKES A POLAR ANGLE 'PHI' WITH THE C X-AXIS. THE TRANSVERSE E-FIELD IS THAT PARALLEL TO A C SIMILAR STRIP WHOSE NORMAL VECTOR MAKES A POLAR ANGLE C 'PSI' RELATIVE TO THE X-AXIS. 'KZ' IS THE LONGITUDINAL C WAVENUMBER. C C USE DUALITY AND ROUTINE 'HTPKT' C COMPLEX HTPKT REAL KZ DATA ETA2/141925.4929/ C ETPJT=HTPKT(X,Y,W,PHI,PSI,KZ)*ETA2 RETURN END C C ------------------------------------------------------- C COMPLEX FUNCTION ETPJZ(X,Y,W,PHI,PSI,KZ) C C ETPJZ: RETURN THE ET-FIELD AT (X,Y) PRODUCED BY AN INFINITE STRIP C OF JZ-CURRENT-DENSITY. THE STRIP IS OF CROSS-SECTIONAL C LENGTH 'W' AND IS CENTERED AT THE ORIGIN, SO THAT ITS C OUTWARD NORMAL VECTOR MAKES A POLAR ANGLE 'PHI' WITH THE C X-AXIS. THE TRANSVERSE E-FIELD IS THAT PARALLEL TO A C SIMILAR STRIP WHOSE NORMAL VECTOR MAKES A POLAR ANGLE C 'PSI' RELATIVE TO THE X-AXIS. 'KZ' IS THE LONGITUDINAL C WAVENUMBER. C C EXTERNAL FUNCTIONS 'BJ1' AND 'BY1' REQUIRED FOR FIRST- C ORDER BESSEL AND NEUMANN FUNCTIONS OF REAL C ARGUMENTS. 'ROM1D' PERFORMS A NUMERICAL INTEGRATION OVER C 'BY1' THROUGH THE FUNCTION 'FIM5' C C COMMON BLOCK 'INTE2' IS USED TO PASS DATA TO SUBROUTINE C 'FIM5' FOR NUMERICAL INTEGRATION C REAL K,KZ,KT,RO(15) COMMON /INTE2/ U,V,CPH,SPH,CPS,SPS,KT EXTERNAL FIM5 C DATA K/6.283185308/,ETA/376.73/ R=SQRT(X*X+Y*Y) WO2=W*0.5 wot=w*0.1 IF(R.LT.wot) THEN C C SOURCE AND OBSERVATION CELLS COINCIDE C ETPJZ=(0.,0.) C ELSE C C SOURCE AND OBSERVATION CELLS DIFFER C KT=SQRT(ABS(K*K-KZ*KZ)) U=X V=Y CPH=COS(PHI) SPH=SIN(PHI) CPS=COS(PSI) SPS=SIN(PSI) IF(K.GT.KZ) THEN ETPJZ=CMPLX(BJ1(KT*R)*W*(Y*CPS-X*SPS)/R, + -ROM1D(FIM5,-WO2,WO2,.01,1.0E-5,15,2,RO,IER)) + *CMPLX(0.,0.25*KT*ETA*KZ/K) ELSE WRITE(*,*) 'OUT OF VISIBLE RANGE -- ERROR' STOP ENDIF ENDIF RETURN END C C ------------------- C FUNCTION FIM5(S) C C SEE DOCUMENTATION FOR 'ETPJZ' C REAL KT COMMON /INTE2/ U,V,CPH,SPH,CPS,SPS,KT DU=U+S*SPH DV=V-S*CPH R=SQRT(DU*DU+DV*DV) FIM5=BY1(KT*R)*(DV*CPS-DU*SPS)/R RETURN END C C ---------------------------------------- C COMPLEX FUNCTION EZPJT(X,Y,W,PHI,KZ) C C EZPJT: RETURN THE EZ-FIELD AT (X,Y) PRODUCED BY AN INFINITE STRIP C OF JT-CURRENT-DENSITY. THE STRIP IS OF CROSS-SECTIONAL C LENGTH 'W' AND IS CENTERED AT THE ORIGIN, SO THAT ITS C OUTWARD NORMAL VECTOR MAKES A POLAR ANGLE 'PHI' WITH THE C X-AXIS. 'KZ' IS THE LONGITUDINAL WAVENUMBER. C C USE DUALITY AND ROUTINE 'HZPKT' C COMPLEX HZPKT REAL KZ DATA ETA2/141925.4929/ C EZPJT=HZPKT(X,Y,W,PHI,KZ)*ETA2 RETURN END C C ---------------------------------------- C COMPLEX FUNCTION ETPKT(X,Y,W,PHI,PSI,KZ) C C ETPKT: RETURN THE ET-FIELD PRODUCED AT (X,Y) BY AN INFINITE C STRIP OF CONSTANT UNITY-AMPLITUDE KT-CURRENT DENSITY. C THE STRIP IS OF LENGTH 'W' AND IS CENTERED AT THE C ORIGIN, SO THAT ITS OUTWARD NORMAL VECTOR MAKES A C POLAR ANGLE 'PHI' RELATIVE TO THE X-AXIS. C THE TRANSVERSE E-FIELD IS THAT PARALLEL TO A C SIMILAR STRIP WHOSE NORMAL VECTOR MAKES A POLAR ANGLE C 'PSI' RELATIVE TO THE X-AXIS. 'KZ' IS THE LONGITUDINAL C WAVENUMBER. C C USE DUALITY AND ROUTINE 'HTPJT' C COMPLEX HTPJT REAL KZ ETPKT=-HTPJT(X,Y,W,PHI,PSI,KZ) RETURN END C C ---------------------------------------- C COMPLEX FUNCTION HZPKZ(X,Y,W,PHI,KZ) C C HZPKZ: RETURN THE HZ-FIELD PRODUCED AT (X,Y) BY AN INFINITE C STRIP OF CONSTANT UNITY-AMPLITUDE KZ-CURRENT DENSITY. C THE STRIP IS OF LENGTH 'W' AND IS CENTERED AT THE C ORIGIN, SO THAT ITS OUTWARD NORMAL VECTOR MAKES A C POLAR ANGLE 'PHI' RELATIVE TO THE X-AXIS. 'KZ' IS C THE LONGITUDINAL WAVENUMBER. C C USE DUALITY AND ROUTINE 'EZPJZ' C COMPLEX EZPJZ REAL KZ DATA ETA2/141925.4929/ C HZPKZ=EZPJZ(X,Y,W,PHI,KZ)/ETA2 RETURN END C C ---------------------------------------- C COMPLEX FUNCTION ETPKZ(X,Y,W,PHI,PSI,KZ) C C ETPKZ: RETURN THE ET-FIELD AT (X,Y) PRODUCED BY AN INFINITE STRIP C OF KZ-CURRENT-DENSITY. THE STRIP IS OF CROSS-SECTIONAL C LENGTH 'W' AND IS CENTERED AT THE ORIGIN, SO THAT ITS C OUTWARD NORMAL VECTOR MAKES A POLAR ANGLE 'PHI' WITH THE C X-AXIS. THE TRANSVERSE E-FIELD IS THAT PARALLEL TO A C SIMILAR STRIP WHOSE NORMAL VECTOR MAKES A POLAR ANGLE C 'PSI' RELATIVE TO THE X-AXIS. 'KZ' IS THE LONGITUDINAL C WAVENUMBER. C C USE DUALITY AND ROUTINE 'HTPJZ' C COMPLEX HTPJZ REAL KZ C C *** NOTE: FOR THE SELF-TERM, THIS ROUTINE PRODUCES THE ET- C FIELD AN INFINITESIMAL DISTANCE OUTSIDE THE SOURCE, AS C DEFINED BY THE NORMAL VECTOR -- THIS VALUE IS -0.5 C C TO GET THE ET-FIELD JUST INSIDE THE SOURCE, C SUBSTITUTE THE VALUE +0.5 FOR THE SELF-TERM C ETPKZ=-HTPJZ(X,Y,W,PHI,PSI,KZ) RETURN END C C ------------------------------------------------------- C COMPLEX FUNCTION HTPKZ(X,Y,W,PHI,PSI,KZ) C C HTPKZ: RETURN THE HT-FIELD AT (X,Y) PRODUCED BY AN INFINITE STRIP C OF KZ-CURRENT-DENSITY. THE STRIP IS OF CROSS-SECTIONAL C LENGTH 'W' AND IS CENTERED AT THE ORIGIN, SO THAT ITS C OUTWARD NORMAL VECTOR MAKES A POLAR ANGLE 'PHI' WITH THE C X-AXIS. THE TRANSVERSE H-FIELD IS THAT PARALLEL TO A C SIMILAR STRIP WHOSE NORMAL VECTOR MAKES A POLAR ANGLE C 'PSI' RELATIVE TO THE X-AXIS. 'KZ' IS THE LONGITUDINAL C WAVENUMBER. C C USE DUALITY AND ROUTINE 'ETPJZ' C COMPLEX ETPJZ REAL KZ DATA ETA2/141925.4929/ C HTPKZ=ETPJZ(X,Y,W,PHI,PSI,KZ)/ETA2 RETURN END