389741-5-T SCATTERING BY A NARROW GAP IN AN IMPEDANCE PLANE J.R. Natzke Radiation Laboratory Department of Electrical Engineering and Computer Science The University of Michigan Ann Arbor, Ml 48109 November, 1989 McDonnell Aircraft Company St. Louis, MO 63166 389741-5-T = RL-2571

ABSTRACT For a plane wave incident on a cavity-backed gap in an impedance plane, the coupled integral equations for the induced currents have been solved numerically and the far field scattering computed. The results are compared with a quasi-analytical solution previously derived and modified for the impedance plane. For narrow gaps of widths less than 0.15k, the agreement is within 12 percent for H-polarization and 14 percent for Epolarization for the cavity geometries considered, limited to small surface impedances of the plane. Excellent agreement is obtained when the material filling of the gap is lossy.

TABLE OF CONTENTS 1. Introduction 2. Theoretical Development of the Integral Equations 2.1. H-polarization 2.2. E-polarization 3. Application of the Quasi-Analytical Solution 4. Numerical Results 5. Conclusions References 1 2 2 8 10 16 28 29 Appendix A: Appendix B: Evaluation of the Half Space Green's Function of an Impedance Plane Moment Method Solution of the Coupled Integral Equations 30 B.1 H-Polarization B.2 E-Polarization B.3 Program Listings 34 34 38 39 Appendix C: Program Listing for the Quasi-Analytical Solution 52

1. INTRODUCTION A case of interest in radar cross section studies is the scattering from gaps and cracks in planar surfaces. Results were recently derived [1,2] for the narrow gap in a perfectly conducting ground plane. Of equal concern is the scattering from a gap of similar geometry but in a ground plane coated with a material having arbitrary dielectric properties. The development of the solution in this paper is for planar surfaces large in extent compared to the width of the gap. As in [2], a set of coupled integral equations for the electric and magnetic currents which exist on the walls of the gap cavity and in aperture of the gap are developed by employing the equivalence principle [3]. Since a thinly coated conducting surface can be described by the impedance boundary condition [4], the half-space Green's function for an impedance plane has been developed and is applied. For plane wave incidence the integral equations are derived for a cavity of arbitrary shape filled with a homogeneous material. The equations are solved by the moment method and data for several cavities are presented. A quasi-analytical solution was derived in [5] for the far field scattering from a uniform resistive or impedance insert in a perfectly conducting plane. In [2], the solution was applied to a narrow gap in the conducting ground plane. Accurate results were obtained by defining the surface impedance of the gap as the input impedance of the gap modeled as a shorted transmission line. Modifications are made to the quasi-analytical solution of [5] such that it can predict the scattering from a gap in an impedance plane with a small surface impedance. The modified quasi-analytical solution is applied to the gap in the impedance plane for gap widths which are electrically small, and the results are 1

compared with those obtained using the coupled integral equations. The limitations are determined for which this quasi-analytical solution provides an accurate design tool. 2. THEORETICAL DEVELOPMENT OF THE INTEGRAL EQUATIONS The gap geometry under consideration is the two-dimensional one shown in Fig. 1. The plane y = 0 is an impedance plane of infinite extent for jxj > w/2 with a surface impedance a'. The plane y = 0 for Ixl < w/2 contains the aperture A of the gap cavity whose walls S are perfectly conducting. The cavity is filled with a homogeneous dielectric material of permittivity el = cre and permeability Ii- = |rW, where the quantities without subscripts refer to free space. The incident plane wave is -i -i A -ik(x cos4o + y sin) ( H,E =ze (1) for H- and E-polarizations respectively, where k is the propagation constant in the free space region above the surface. A time factor e-i(tt is assumed and suppressed. 2.1 H-Polarization The equivalence principle is applied to the gap for the region y > 0 by shorting the gap with a perfect electric conductor in A and placing a magnetic current = - y x E(x,0) over it. The plane y = 0 is then an impenetrable surface with mixed boundary conditions: that of a perfectly conducting surface for lxi < w/2 and that of an impedance surface for Ixl > w/2. For the tangential magnetic field, the impedance boundary condition is given by [4] 2

y - w k.-.- -.- c-..I..[,_ x A Figure 1. Narrow gap of arbitrary shape in an infinite ground plane. -Zin iw —w-. Zin I (a) (b) Figure 2. Rectangular and triangular gap geometries and dimensions. 3

H=0 (2) ay H where Z is the free space impedance. By applying Green's theorem, an integral equation for the total magnetic field can be written as Hz(x,y) = H,(xy) + JJx(x') -G(x,y;x',y') + ikYJz(x') G(x,y;x',0) dx' A 2y' —O -w/2 oo + J + IfH,(x',O) -jG(x,y;x',y') - G(x,y;x',O) A (x',y') dx' - W!2 y' ---O where G is the appropriate free space Green's function. This is certainly not a desirable integral equation to solve with the path of infinite extent over Ixl > w/2. However, this path of integration over the y = 0 plane can be eliminated by applying the half space Green's function for an impedance plane, which also satisfies (2). The half space Green's function for the impedance plane is G = 4 {Ho {k(x-x) + (yy,)2J 1 r () ek[(x-x')cosa+ly+y'lsina] (3 -oG= i k(Xx+ Ho - J F (3) where H(1) is the zeroth order Hankel function of the first kind and TY- sino H() = Y + sina is the reflection coefficient for an H-polarized plane wave incident on an impedance boundary at the angle a. The path C is in the complex a plane, defined in Fig. 3. Note that as il approaches zero, (3) becomes the Green's function for a perfectly conducting plane for the H-polarization case. For the 4

numerical solution of (3), a more quickly converging integral gives G=-Ho{k/(x-x')2 + (yk-y)2j +k (x-x')2 +(y+y')2 - 2 J e Ho k (x-x')2 (y+y'+iv)2 dv (4) 0 where kir TH - 'H- Z Equation (4) was derived using a transform technique presented in [6]. For such a Green's function to be applied, the continuity of the impedance boundary condition over the entire y = 0 plane is necessary. This is accomplished by impressing a magnetic current source equivalent to fYHz over a = d + i a" C ~l - -- Ia - c 2 2 a' 7E Figure 3. Path of integration in the complex a plane. 5

A. Applying then the half space Green's function of (4), Hz(x,y) = H;(x,y) - H (x,y) + ikY JJ;(x') G(x,y;x',0) dx' A ik Hz(X',O) G(x,y;x',O) dx' (5) A where Hr T) Y - sin4() -ik(x cosoo-y sin(o) Hz x,y) e z lY + sin(o is the reflected plane wave, and Y = 1/Z. The second integral in (5) is the correction term for the additional impressed current, ensuring the original boundary condition for the total field at the shorted gap. That is, it removes the contribution of the scattered field due to the surface rj over A which was not part of the original system. Observing the field in the aperture, 2 sin(o)o -ikxcos(o)o H(xO) s in e s ikxcoso+ ikY [J(x') - iJs(x') G(x,O;x',0) dx' (6) A where Js(x') = [y x H(x',0)] * from (5). As shown in Fig. 1, s is tangential to the surface of the cavity wall, and in the aperture, s = x. For the region y < 0 occupied by the cavity, a magnetic current -J is placed just below A to ensure the continuity of the tangential electric field in the open gap. The expression for the magnetic field in the shorted cavity is the same as that constructed in [2] since the region is independent of the impedance boundary. The integral equation for the currents on the cavity walls is [2,Eq. (7)], 6

ik E= ( r JSJ(x) H(k1 (x-x) +y2 dx' A + JS(s') sinf ( H kl (x-x') + (y-y') ds' (7) 2 S+A where k1 = k er and n' (x-x') + (y-y')9 (8) sin p' =z z~ (8) 2 2 (x-) + (y-y)2 valid at all points of S and A. To ensure the continuity of Hz through the aperture, the field expression for (7) [2,Eq. (9)] is matched to (6), resulting in Js(x) = sin e-ikxcos + ikY J(x') -i Js(x') G(x,0;x',O) dx' (9) JY + sinoo A A valid for x in A, and (7) and (9) constitute a pair of coupled integral equations for JZ and JS In the far field, the large argument expression for the Green's function of (3) is found by saddle point integration, allowing C to become the path of steepest descent about o = i/2. The scattered magnetic field is then given as H= 2 (k-/4) (10) Z nkp where = -kY 1 [JZ(X')- J(')] ekx' dx' (11) A Thus the far field amplitude is given for the currents JZ and Js over A determined from (7) and (9). 7

2.2 E-Polarization The impedance boundary condition for the total tangential electric field is [4] a+ i -- Ez = 0 (12) The half-space Green's function for an impedance plane with an E-polarized incident plane wave satisfying the same boundary condition is G = {H{k(x-x + (y-y)+ ra) eik[(x-x')cosa+l+ ysina da (13 where ( Y - csco E = / Y + csca Note again that as il approaches zero, (13) becomes the Green's function for a perfectly conducting plane under E-polarized illumination. For a quickly convergent integral, the same expression is used as that given in (4), with TH replaced with kZ 1E - Using the same approach as for the H-polarization case, the total field for the region y > 0 is Ez(xy) = Ez(xy) + Ez(xy) - xJ;(x') -G(xy;x',y') dx' a y y' —0 A I- E(x',0) —G(x,y;x',y') dx' (14) JY' y' —0 A where J is the assumed equivalent magnetic current on A and 8

r lY - csc(o e -ik(x cosco - y sino) Ez(x,y) = #fly + c e zv+ csca0 is the reflected plane wave. The second integral in (14) accounts for the scattering due to the added impressed current equivalent to Ez(x',O) for the impedance plane. The dependence of this current on the surface impedance is shown more explicitly using [4] the boundary condition Ez = - qHx.,Y aEFrom Hx = iay, the tangential component of the magnetic field in the aperture is 2 Y -ikx cos%, X' ly) Y + csco iY ay [Jx(X') +TJz(X) (xy;x',')l dx. (15) A In the region y < 0 occupied by the cavity, the field source is the equivalent magnetic current JX in the aperture and the contribution from the electric current Jz on the walls of the cavity. The integral equation constructed for the currents is [2,Eq. (15)] Y A ) 2 2 J() =2 (n V) a JX(X) H k1 (x-x') + y dx' A ik1(1) I 2 2 + 2 Jz(s') sin3 H k1 k (x-x')2 (y-y)2 ds' (16) S+A where (x-x') + (y-y')y A sin = ~ xs (17) (x-x)2+ (y-y')2 valid at all points of S and A. When (15) is matched to the field expression for (16) with the observation point in the aperture, 9

J) 2 Y -ikx cos z l + csco e +- - lim k y-O ayu [;() + Jz(x')] aG(xy;x',y') 'dx' (18) A for x in A. Equations (16) and (18) are a pair of coupled integral equations for the currents Jx and Jz. A third integral equation was developed in [2] for the electric field in the cavity of the gap. When the boundary condition on the perfectly conducting surface was applied, the expression for the currents on A and S came to be Jx =(X') H1, (x -(Xx') 'H+ y2) dx' A kZLt r 2 + kZ2 JZ(s) Ho (k x-x (yY) ds' (19) S+A where Jx(x) is non-zero only in A. Equations (18) and (19) are the pair of integral equations used to compute JX and JZ. The scattering in the far field for E-polarization has a similar expression as that in (10), giving the far field amplitude k sin I[J;() + Jz(x')] ek dx'. (20) +(090) 2;q e A 3. APPLICATION OF THE QUASI-ANALYTICAL SOLUTION Consider an impedance insert of width w and characterized by q in an impedance plane f. For the H-polarization case, as in Section 2.1, the application of Green's theorem for the region y > 0 gives 10

w/2 Hz(x,y) = Hz(x,y) - Hz(x,y) + i- J (1 - T) Hz(x',O) G(x,y;x',O) dx' (21) -w/2 where the incident and reflected fields are those described in (5) and G is the Green's function of (4), satisfying (2) in the plane y = 0. The integral of (21) is the scattering due to the impedance insert, and this clearly vanishes for i1 = A. The Green's function of (4) can be written in the form, for y' = 0, G(x,y;x',0) (22) where the normalized half space impedance Green's function is LH(x,y;x',01|) g(x,y;x',O|l) = 1 -...2. for which the function LH is the integral of (4). We note that as f approaches zero, LH approaches zero, and as fr approaches infinity, LH approaches H(l). Thus Igj varies from one to zero for these limits. To remove the coordinate dependence, g is averaged over -w/2 < x < w/2, for x' and y set to zero. By curve fitting the numerically generated magnitude and phase, this average is found to be approximated analytically to within 5 percent as gH(WI) = e e (23) where c(w) = 0.245 + 1.267w/k V(wi) = 1- e'(0'098+1 760w/k) i/Z for 0.02 < w/X < 0.15 and f/Z < 2. 1 1

Now taking the observation point to be on the insert, (21) can be approximated by 2 sin^o -ikx cos% Hz(X,) iY + sinko e w/2 -- -) gH (li) J Hz(X',0) Ho (k Ix-x'l) dx' (24) -W/2 where the averaging of g has enabled us to remove the fl dependence from the integral. The assumption that g is constant over w should hold for small w and small 'f. The integral equation (24) can be written in the form [5] 1 k. coso -i J(C) Ho k I-'l d'= e2 +aJl() (25) where j()_ikw iq-~ T y + sino(XO 1 2 Z 2 sn gH(wIi) HZ(xO) and a 2 Z 1 (26) -kw T- gH (WI (26)) Referring to (11) and the integral of (21), the far field amplitude written in terms of J1 is PH( ) w= Wf) Y+ 1 sino fJ1(') e2 dC' (27) -1 A simplification is made using [5] J l() = [1 + Hw (TIY + 1) jY+sin~ PoH('I) J2(0 ) where the integral equation that the modified current J2 satisfies is [2,5] 12

1 1 -JJ 2() - d = 1 + aJ2(C) -1 for-1 < C < 1 and kw A=n + Y+-is 4 2 where y = 0.5772157.... is Euler's constant. The far field amplitude (27) becomes sino - 1 PH( -,si) = 1 A + 1 (28) qY +1 gH(wlIi) TjY + sinoo [ KH(a) with 1 KH(a) =J,2() d, -1 for which an approximate expression is [2] 1 K(a) (29) H -+ n. 2 +0.1 2 with a given in (26). Thus with the modifications to a and PH, the same quasianalytical expressions are used to solve for the scattering of the impedance insert in an impedance plane. Similarly, for E-polarization, the total field for the region y > 0 is Ez(x,y) = E'(x,y) + Ez(x,y) w/2 +J [E,(x',0h) - Ez(x',OIi)1 G(xy;x,y')I dx', (30) -w/2 and the correlation of the field with the respective impedance surface is shown explicitly. The corresponding expression for the tangential magnetic field over the insert is 13

2 Y -ikx cos4o (X')' Tly + csc)o w/2 li m ( - T)Hx(x',O)G(x,y;x',y')l dx' (31) kZ y —o ay Jaiy O' ---Y)I d -w/2 where G is (13). When (13) is expressed in the form of (4), the partial derivatives in (31) render the source and image Hankel functions zero for y' = 0, and G is given by (22) with LE (x,y;x',I0i) g(x,y;x',0l ) = Ho (kx-xI)2+y2j where LE is the integral of (4) with TH replaced with tE. For this case we note that as ln approaches zero, LE approaches H(), and as Tr approaches infinity, LE approaches zero. Thus Igi varies from one to zero for these limits. The average over w of the normalized half space impedance Green's function is approximated analytically within 10 percent as gE(WI) = e( e (32) where d(w) = 1.558- 4.226w/? i(wlT) = - (0.380 - 0.80w/) [1 - e-2586/Zj for 0.025 < w/X < 0.15 and 1/Z < 2. 14

The tangential magnetic field on the insert is then approximated as 2 Y -ikx cos0o H(x)- y + CSC e i 2 I2 (1) - (nI - ) g(WI) tk2 + a2)J J Hx(x',O) H )(k Ix-x') dx. (33) Z DX2-w/2 Given this approximation, (33) can be equated to [5] 2 1 kw + 3( ') j J Ho1Cw C-= e- bJ( - ) ac2 2-2i132 3-d' ~ with 2i l-i fly +csc3 3()w z 2Y gE(WIn) Hx(x,O) (34) and ikw Z 1 b =- (35) 2 T- gE(wT) (35) Referring to (20) and using (32), the far field amplitude for the E-polarization case is P i( ) (kw)2 sin 1 1 (b) (36) PE(') -4 flY +1 iY+csc4% gE(WIi) KE(b) 36) with KE approximated by [2,5] K (b) 0.62 (b + 4.08)(b + 7.26)(b + 10.37)(b + 13.43)(b + 16.46) (37) E( b + 1.15 (b + 4.27)(b + 7.37)(b + 10.45)(b + 13.49)(b + 16.50) As in [2], the analogy is drawn from the impedance insert to the narrow gap by equating the surface impedance ri to the input impedance of the gap modeled as a transmission line. The expressions for various gap and cavity 15

configurations are contained in [2], and the far field amplitude PH and PE can then be calculated accordingly by (28) and (36), respectively. 4. NUMERICAL RESULTS The integral equation pairs (7), (9) and (18), (19) for H- and Epolarizations respectively were programmed for solution by the method of moments, using pulse basis and point matching functions. The half space impedance Green's function of (4) and its derivatives were evaluated analytically in handling the singularities of the integration and numerically otherwise as described in Appendix A. The application of the method of moments is described in Appendix B, which also contains the program listing used for generating the results. The quasi-analytical expressions of Section 3 were programmed for solution, as listed in Appendix C. The upper limit of w/k = 0.15 was determined for the applicability of the quasi-analytical solution in [2], and this limit is maintained for the results listed here. In Figs. 4 and 5 the magnitude and phase of the far field amplitude PH(x/2, n/2) are shown as a function of depth for a rectangular air-filled gap of width w/k = 0.15, comparing the method of moments (MoM) and quasianalytical (QA) solutions. For each method it was verified that as fr approaches zero, the results approach that of the gap in a perfectly conducting plane in [2]. For f/Z = 0.1 and 0.5 the difference between the peak amplitudes of each method are within 12 percent, and the phase curves show excellent agreement. As expected, IPHI is non-zero as d approaches zero, corresponding to the scattering from a perfectly conducting strip in an impedance plane. Consistent results were verified for ff/Z = 0.5 and w/k as small as 0.025. As observed in [2], a cyclical behavior exists with increasing gap depth resulting from the periodicity of the impedance looking into the gap. 16

2.5 2 1.5 -0IP I I 0.0 0.1 0.2 0.3 0.4 Gap depth d/X Figure 4. Modulus of the far field amplitude PH for a rectangular gap of varying depth d1 d with = = 0=1/2 and w/k =0.15: I/Z = 0.1 * MoM, QA q/Z = 0.5 * MoM, - - - - QA 0.5 17

arg PH 0o (deg.) -60 -120 -0.0 0.1 0.2 0.3 0.4 Gap depth d/A Figure 5. Argument of the far field amplitude PH for a rectangular gap of varying depth di = d with 4 = 0- = i/2 and w/A. = 0.15: Tf/Z = 0.1 * MoM, QA Ti/Z = 0.5 0 MoM, - - - - QA 18

The far field amplitude backscatter response to the rectangular gap with d/k = 0.15 is contained in Fig. 6, showing excellent agreement between the MoM and QA solutions for all 4. Figure 7(a) shows the far field amplitude of the rectangular gap filled with a lossless dielectric having ~r = 2.5 for /Z = 0.5 and 1.0. It was observed that as the relative permittivity of the gap filling was increased from 1, the prediction by the QA method improved, bringing the difference at the peaks within 14 percent for the f/Z = 1.0 curve shown in the figure. The agreement is improved when some loss is introduced in the dielectric filling, as shown in Fig. 7(b) for ~r = 3 + i0.5. For the V-shaped gap of Fig. 2(b), the far field amplitude is presented in Fig. 8 for varying depth. Similar results are expected for gaps of arbitrary geometries, given the appropriate surface impedance rT of the gap necessary for the QA solution. For the E-polarization case, the magnitude and phase are plotted in Figs. 9 and 10 for the rectangular air-filled gap for Ij/Z = 0.1 and 0.3. The curves reveal the evanescent nature of the fields in the gap cavity, giving constant values for d/k > 0.1. In Fig. 9, the greatest deviation occurs as d approaches zero, for which the MoM predicts a decrease in magnitude. For f/Z = 0.3, the difference between the two methods is within 14 percent. The QA method breaks down then for rI/Z > 0.4 since the amplitude continues to increase for all d. Figure 10 shows that the phase information is lost in using the QA method. A plot of the backscatter from an air-filled rectangular gap with d1/X = 0.2 is contained in Fig. 11, and Fig. 12 shows the far field amplitude of the scattering from an air-filled V-shaped gap. 19

1.5 - IPHI 1 -/ v /./ 0.5 -0 0 15 30 45 60 75 0 in degrees Figure 6. Modulus of the far field amplitude PH for the backscatter from a rectangular gap of depth d 1/ = 0.15 with o =, w/k = 0.15: I/Z = 0.1 MoM, QA i/Z = 0.5 * MoM, - - - - QA W 20

2.5 IPHI Gap depth d/X Figure 7(a). Molulus of the far field amplitude PH for the material-filled rectangular gap of varying depth d1 = d with Er = 2.5,,r = 1, = = 0 = i7/2, and w/X = O. 15: T/Z = 0.5 /Z = 1.0 0 * MoM, QA MoM, - - - - QA 21

2.5 IPHI 2 1.5 1 -0.5 -0- - 0.0 Figure 7(b). 0.1 0.2 0.3 0.4 Gap depth d/X 0.5 Modulus of the far field amplitude P H for the material-filled rectangular gap of varying depth d1 = d with cr = 3 + i0.5, gr = 1, = 0 = 7c/2, and w/ = 0.15: j/Z = 0.5 i/Z= 1.0 0 22 MoM, QA MoM, - - - - QA

2.5 1.5 -IPHI 0.5 0.0 0.1 0.2 0.3 0.4 0.5 Gap depth d/X Figure 8. Modulus of the far field amplitude PH for an air-filled V-shaped gap of varying depth d1 = d with 4 = 0 = 7=/2 and w/X = 0.15: /Z = 0.1 I MoM, - QA i/Z = 0.5 * MoM, - - - - QA 23

0.20 0.15 IPEI 0.10 0.05 -0.00- ', ' -, * - 0.0 0.1 0.2 0.3 0.4 Gap depth d/A Figure 9. Modulus of the far field amplitude PE for an air-filled rectangular gap of varying depth d1 = d with = = ~0 = x/2 and w/X = 0.15: i/Z = 0. MoM, -- QA r-/Z = 0.3 * MoM, - - - - QA 24

0 -30 - " OZ.-, 0 0 0 0 0 o 0 0 a 0 0 a 0 0 0 0 e -60 -irg PE (deg.)9 -120 --120 - a 0 a 00 a 0 a 0m 0 0.IIl / J ^ / / J' -1fi. I I -aY I. I. I. 0.0 0.0 0.1 0O2 Gap depth d/A 0.3 0.4 Figure 10. Argument of the far field amplitude PE for an air-filled rectangular gap of varying depth d1 = d with ( = 0o = 7/2 and w/; = 0.15: T/Z = 0.1 /Z = 0.3 0 MoM, -- QA MoM, - -- QA 25

IPEI J4 lop/ 0.10-,^ 0.05 0.00 0 15 30 45 60 75 0 in degrees Figure 11. Modulus of the far field amplitude PE for the backscatter from a rectangular gap of depth d /k = 0.2 with % = X, w/k = 0.15: l/Z = 0.1 * MoM, QA j/Z = 0.3 * MoM, - - - - QA 26

0.26 0.20 -0.15 IPEI " 0.10 0.05 0.00- - 0.0 0.1 0.2 0.3 0.4 Gap depth d/A Figure 12. Modulus of the far field amplitude PE for an air-filled V-shaped gap of varying depth d1 = d with e = %0 = X/2 and w/X = 0.15: x/Z = 0.1 * MoM, ---- QA r/Z = 0.3 ~ MoM, - - - - QA 27

5. CONCLUSIONS The quasi-analytical solution of [2] was based on the low frequency approximation of the integral equations for a constant impedance insert in a perfectly conducting plane applied to a cavity-backed gap in the same such plane. This same solution has been applied to a gap in an impedance plane by modifying the expressions with coefficients dependent on the surface impedance of the plane and the width of the gap. For comparison, a solution was derived by the equivalence principle in conjunction with the half space impedance Green's function, and the exact integral equations were solved by the method of moments. For an air-filled gap, the quasi-analytical solution for the far-field amplitude is within 12 percent for an I/Z < 0.5 for H-polarization and within 14 percent for an ff/Z < 0.3 for E-polarization, assuming w/X < 0.15. The results are improved when the gap is material-filled with a relative permittivity greater than 2.5, within 14 percent for f/Z < 1 for H-polarization. The quasianalytical solution gives excellent agreement with the method of moments results when the gap is filled with a lossy dielectric material. Thus the quasianalytical solution is a simple and accurate method for determining the scattering from a narrow gap in an impedance plane for surface impedane for surfac impedances of small but significant values. 28

REFERENCES [1] T. B. A. Senior, K. Sarabandi, and J. R. Natzke, "Scattering by a Narrow Gap," to be published in IEEE Trans. Antennas and Propagation, August, 1990. [2] T. B. A. Senior, K. Sarabandi, and J. R. Natzke, "Scattering by a Narrow Gap," Radiation Laboratory Report 389741-3-T, The University of Michigan, Ann Arbor, April, 1989. [3] R. F. Harrington, Time Harmonic Electromagnetic Fields, McGraw-Hill Book Co., New York, 1961. [4] T.B.A. Senior, "Impedance Boundary Conditions for Imperfectly Conducting Surfaces," Applied Science Res., sec. B, vol. 8, p. 418, 1960. [5] T. B. A. Senior and J. L. Volakis, "Scattering by Gaps and Cracks," IEEE Trans. Antennas and Propagation, vol. AP-37, 1989. [6] K. Sarabandi, "Scattering from Variable Resistive and Impedance Sheets," Radiation Laboratory Report RL-863, The University of Michigan, Ann Arbor, March 13, 1989. 29

APPENDIX A Evaluation of the Half Space Green's Function of an Impedance Plane The half space Green's function for an impedance plane is given in (3) and, considering the integral, in a more quickly convergent form in (4). For the currents and observation in the gap, y' = y = 0, and (4) becomes., (k lxx) -V (1H ) (A.1) G(x,O;x) = H (k x-x ') - lim e H ( + (y+iv) dv (A1) 0 The Hankel function is well defined, but the integral must be analyzed carefully, since its integrand contains a singularity at v = Vo = lx-x'l in the limit as y approaches zero. Using the small argument approximation of the Hankel function, the integral about this point is vO+Av LH(X'O;x',OI) = lim J te [' + i - (x-')2 + (y+iv)2 dv (A.2) Vo-AV where A'=1 + - Y+y+t -. (A.3) Assuming Av is small such that the exponential term is nearly constant over 2Av, (A.2) is evaluated analytically in the limit as LH = THe (L + L2) (A.4) where L= = 2A' Av and 30

L2 = { (Vo+Av)[tmr(2vo Av + Av2) + ib] - (Vo-AV) Ut(2vo Av- AV2) L2 = &m(2v2Vo + Av -4Av + V, { 2vo. -Ai).j 2vo - Av Over the other portions of the path of integration, LH is evaluated numerically. For the case x = x', we have vo = 0, and the evaluation of (A.1) can be done as follows. In the application of the moment method, the path of the integral in (9) is segmented such that the currents are assumed to be constant over each segment Ax. The currents can then be removed from the integral, and with the self-cell, the integration of LH over Ax is xo+Ax/2 Av J limO i H e H A' + i2 (x 2 (y+i dv y-F20 -H + xo-Ax/2 0 + xH e Ho (x,-x')2 + (y+iv)2dv dx', Av where x0 is the midpoint of the self-cell. Evaluating the first integral within the brackets analytically, the self-cell expression of LH becomes x +Ax/2 SH(AXIi) = LH + 2 lim | eTe'H H k (x-x') + (y+iv) dv dx' (A.5) y — oH0O Xo Av where Av LH = H e 2 (L3 + L4) (A.6) with L=-andAA and 31

4 i 2 Ax-2Av) +i iu+ Av + } 4 Ax OAx - 2Av) i 2 + 4 The double integral in (A.5) is evaluated numerically. The self-cell expression of the Hankel function in (A.1) is well known and is given in Appendix B as needed. For the E-polarization case, the derivatives of the Green's function must be considered. From (18), and using the form in (4) for (13), lim,G(x,yx,y)= )im k2= i It2e' kxxd (A.7) li ' ', G(xy;x I= limO k TE Ho (kX-X)2 + (y+iv)2 dv (A.7) y' —oay ay' E2-y 2 0 The order of the differential operator is reduced by applying the integral shown in (18) over the segment Ax. Given the endpoints X1,X2 of the segment, the integral LE of (A.7) becomes X2( J J e EV H(1 k (x-x,)2 + (y+iv)2J dv dx' x1 o + k JTEe ( 2) i H1 (xx2) +(y+iv) J 0 H (x-x2) +(y+iv) (x-x1) H ({2(x-x)2+(y+iv) ) 1dv, (A.8) (X-X) () 2 2 (x-x) )2+(y+iv)2 dv, (A.8) (X_xl2+(y+iv)2 in the limit as y approaches zero, where k2 has been factored out. For vo not equal to zero, the singularity in the first integral expression of (A.8) is handled in the same way as that for the H-polarization case, giving 32

L~ E Vo (A.9) EL = 'Ee o(L1 +L2) (A.9) with L1 and L2 given above. The integration is done numerically over the remaining path, as well as with respect to x'. Let L'E(x,AX'IT) denote the second integral expression in (A.8). The contribution from the singularity is, using the small argument approximation of the first order Hankel function, 1 E -E Vo2 2o2-A (X-X2) LE = e ~k(x-x2) Av + r- 2v 2A + ic 21 E k 2 k 2Vo2 +AV TVM. 2V -Av (X-X) - e E k(X-X) Av + 1 -1 (A. 10) 1 nk 2v +Av + x-x i where Vol =IX-x11 Vo2 =Ix-21 Considering the integrations of (A.8) over the self-cell, Xo+AX/2 SE(Ax) = LE + 2 lim tE e H k (x-x')2 + (y+iv)2 dv dx xo Av AV + L lim e Ax H1 )k(Ax/2)2+(y+iv)2 dv (A.11) E k y-yoJ E 2 2 1 J Av (Ax/2)2+(y+iv)2 where the analytical expressions are Av OX tE 2 LE =T E e (L3+L4) (A.12) and Ix TE -2 2 Ax + 2(1 LE T e + 2A~ 2I (A. 13) L E k 4 k Ax - 2Av J' 33

APPENDIX B Moment Method Solution of the Coupled Integral Equations The integral equation pairs given by (7),(9) and (18),(19) are solved by the moment method [2,Appendix A]. Using pulse basis functions in the moment method, the aperture A and the cavity walls S shown in Fig. 1 are segmented into N cells of size As. The magnetic and electric currents are assumed to be constant over each of these segments. When the integrations of the coupled equations are taken over each segment, the current expressions can be removed as constants from the integrals. With the contour of integration discretized, the (x',y') coordinates become (x.,y.), i = 1,...,N, which describe the location of each of the segments. The Green's functions can then be expressed in terms of rotated coordinates (s,n) for the observation position and (sinj) for each segment or source position since the integration is with respect to A the tangential vector s as shown in Fig. 1. The expressions for the numerical solution of the coupled equations are developed in the following sections for the H- and E-polarization cases. Applying point matching, the magnetic and electric currents in the aperture and on the cavity walls are determined, and the far field amplitude is calculated. A.1 H-Polarization For the discretized contour of integration, (7) and (9) become J5(sn) = E ~r J Jz(Si) Ho k1 (s-s) + n2 ds i= 1 As. ik (1)( 22 + 2 Js(si,ni) sinai H Lk1 (s-s) + (n-nj) dsJ (B.1 ) i=1 As. 34

2 sinf.l -iks cos(o Js(s) fly + sino e 2 Jz(i)- Js(i)] I[Ho (k Is-si,)-LH(s,0;s,0) ds. (B.2) =1 Asi where M are the number of segments across the aperture and (n-ni) sin = (-. (B.3) 2 2 (s-si) + (n-n.) Applying point matching over the N segments of the aperture and cavity walls, ik1 (1) -1. + sin H (k1 Rji ) ds i =N+1 As. kY (1) i= 1 As. kY+M 2 fo rN (cos i [ j = 2 +'[HO (k R,) - LH(sk,0;sX,0)J dsi =0 Rj = (sj) + (n-n)2. (B.6) The coordinate (sj,nj) is the observation position at the midpoint of the jth segment. Hence, for i,j = 1,...,M,N+1,...,N+M, the segments are located in the aperture, and for i,j = M+1,...,N, the segments are located on the cavity walls. l\ in (B.4) and (B.5) are the electric currents, for i = 1,...,N, and the magnetic currents, i = N+1,...,N+M, to be determined. 35

In matrix form, (B.4) and (B.5) become [Z ][|, V, (B.7) The impedance matrix is given as Z, Z, el ml Z (B.8) e2 m2 where the sets of elements are as follows: ik1 (1) 2 sin1ji H1 (k1 Rj) ds, ji Z1 Asi (B.9) -1 j=i for i= 1,...,N and j= 1,...,N; E ~Er Ho (k1 Rj) ds j i-N m Asi (B.10) l S i 1 A1 k 2(s -S)[i (R ) 1+A;] j=i-N for i =,...N+ +M and j = 1,...,N; - - F Ho (k Rj) - LH(sj,O;si,0) dsi j i+N As. Ze2 (B.11) for i,V..)j = i+N 1 - 2 (s1 +2A- ) - SH(ASni(R)1 j = i+N for i=1,...,Nandj=N+1,...,N+M; 36

[HO (k R..)- LH(SO;SO)] ds j i j As Z = (B. 12) kY 2 [2(sj-sj)i/ tm(Rj) - 1 + A'J SH(AsiIT) 1 = for i = N+1,...,N+M and j = N+1,...,N+M. In (B.10), the expression for A; is i2 k A' = 1 +- 2 9 and in (B.11) and (B.12), SH is given by (A.5). For the self-cells in (B.10) to (B.12), si is taken to be the endpoint of the ith segment. The self-cell expressions were derived analytically, and a numerical integration is applied to the other segments. In the case of the V-shaped gap, the adjacent cells needed to be evaluated in the vicinity of y = -d, for Rj,i less than one cell size, and these expressions can be found in Appendix A of [2]. The source matrix is given by 0 j= 1,...,N V = 2 si no -iksj cosO (B.13) {e j= N+1,...,N+M 1y + sin~e and the currents Ii are determined by solving (B.7), given that [Zj,] is nonsingular. From (11), the far field amplitude at the angle p is now kY 1 f -ikxi cos4 pH() N- e dx. (B.14) i= 1 Ax. 37

B.2 E-Polarization The integral equation pair given by (18) and (19) was solved in the same manner as described for the H-polarization case. The elements of the impedance matrix defined in (B.8) for the E-polarization case are as follows: kZ! (1) -1r Ho (k R.ji) dsi ji Ze As (B.15) kr 2 (sH-sj)1 R) dt(Rj.) - 1 + ~ im 1 -1 j = i-N for i =N+1,...,N+M and j = 1,...,N; -2= L,E0(s,O;s,O0|j) dsi + LE(siAsii) j i+N Z.e2 J (B.17) -1 -1kY S(As ) j =i+N for i = 1,...,N and j = N+1,...,N+M; kY [ | LE(sj,O;si,O0|T) ds. + L'E(s,ASjTi) j i z E M2 LAS. (B.18) -2 SE(ASjiI) j=i fori = N+1,...,N+M andj = N+1,...,N+M. In the self-cell expression of (B.15), s; is 38

evaluated at the endpont of the ith segment. In (B.17) and (B.18), SE is given by (A.11). The source matrix is 0 j=1,...,N V. = 2Y -iksj cos (B. 19) T y+csc- e o j=N+1,...,N+M. Given that [Zjj] is nonsingular, the currents Ij can be determined, and the far field amplitude is calculated from the expression k sin (1iN +T I) |je dx.. (B.20) i= 1 Ax B.3 Program Listing The expressions for the impedance and source matrices and the far field amplitude were programmed for solution, as contained in the program listing of IMP.FTN below. The subroutines called by the program in addition to those listed are contained in the file GAPSUB.FTN of Appendix A of [2]. The numerical integration is done for the appropriate segments using Simpson's three-point composite integration over each segment. A segment size of As(/k = 0.01 was used for the results of Figs. 4 to 12. The numerical integration implemented for LH and LE is the Gauss-quadrature technique, and convergence was verified for ir/Z < 2. As listed, the program calculates the far field amplitude as a function of the gap depth using (B.14) for H-polarization and (B.20) for E-polarization. 39

c C IMP.FTN C C C This FORTRAN program computes the far field scattering due C to a narrow gap of arbitrary shape in an infinite impedance C plane. The moment method is applied to solve the currents C of two coupled equations. C C INPUT The user is prompted from the subroutine C GAPROM for the polarization and angle of the C incident field, angle of far field observation, C relative permittivity of gap filling, shape and C dimensions of gap, segment size, number of C iterations with respect to gap depth, and C normalized surface impedance. C C OUTPUT FILES C GAPDAT Contains input data. C AMPDAT Contains the magnitude of the far field. C PHADAT Contains the phase of the far field. C C SUBROUTINES C HANKZ1 Computes the Hankel functions of the first C kind of orders zero and one. C CHANK Computes the Hankel functions of the first C kind of orders zero and one given a C complex argument. C IMPGF Calculates the value of the impedance C plane Green's function. C SIMPGF Calculates the value of the self-cell C for the impedance plane Green's function. C CGECO Factors a complex matrix and estimates the C condition of the matrix. C CGESL Solves the complex set of linear C equations [A] [x] = [b]. C C integer pn parameter(pn=500) integer EorH,N,noS,gN,szN(50) real pi,k,phi,phio,w,d,maxC,etab,q(50,2) real dStp(50),wStp(50),p(pn,2),m(pn, 2),szSd(50),psi(50) real posiX,posiY,spaceX,spaceY,stepX real ni,si,nj,sj,xj,yj,DelS,R,Rm,tanf complex czero,ci,er,ur,Ao,A,kl,cselfO,cselfl,ctempO,ctempl complex Z(pn,pn),LHO,LH1, LHOo,LHlo,aLHO,Ii(pn),Vj (pn) complex eta(pn),Hi(pn),Hs(pn),E(pn),Lsca,Psca complex GFOa(pn),GFOb(pn),GFla(pn),GFlb(pn),GFO(pn),GF1(pn) integer ipvt (pn),iretrn real rc,krho complex wk(pn),HO,H1,HOo,Hlo,ckrho,iHOo,iHlo logical Epol,Lossy,side,neg(50) common /prompt/ EorH,phio,phi,er,ur,igap,wStp,dStp,w,d,noS, & q,maxC,noIter,etab 1 format(il) 2 format(i5) 3 format(g16.8) 4 format(a4) 6 format(2gl6.8) open (1,file='gapdat') open(3,file='ampdat') open(4,file='phadat') c...Declaring constant values czero=cmplx(0.0,0.0) ci=cmplx(0.0,1.) pi=4.0*atan(1.0) k=2*pi Zo=sqrt(4.e-07*pi/8.854e-12) Yo=l./Zo gam=0.5772157 Ao=2*(log(k/2)+gam-ci*pi/2) iprg=l c...Setting default values EorH=l phio=90.0 phi=90.0 er=cmplx(l.,0.0) ur=cmplx(1.,0.0) w=0.15 d=0.2 noS=3 maxC=0.01 noIter=30 etab=l. adj=0.00001 Epol=.false. Lossy=.false. 40

....Prompting user for input data call gaprom(iprg) if(EorH.eq. 1) Epol=.true. phio=phio*pi/180.0 phi=phi*pi/l80.0 drat=dStp (1) /d kl=k*csqrt (er*ur) A=2* (clog (kl/2) +gam-ci*pi/2) if(aImag(er).ne. 0.0) Lossy=.true. dmin=. 025 if(Epol) dmin=0.02 dmax=d if(nolter.ne. l)then dstep= (dmax-dmin) / (noIter-1) d=dmin endif DO 700 iter=l,nolter *... Determining coordinates if(igap.eq. 1)then c RECTANGULAR noS=3 q (3,1) =w/2 q(3, 2)=-d q(4,1)=-w/2 q(4,2)=-d else if(igap.eq. 2) c L-SHAPED noS=7 dStp(l) =d*drat dStp(2) =d* (1-drat q (3,1) =w/2 q(3, 2)=-dStp (1) q(8,1)=-w/2 q (8,2) =-dStp (1) q (4, 1) =w/2+wStp (2 q(4,2)=-dStp(1) q (7,1) =-w12-wStp 4 q(7,2)=q(4, 2) q(5, 1)=q(4, 1) q (5,2)=q(4, 2)-dSt q(6,1)=q(7, 1) q(6,2)=q(5, 2) else if(igap.eq. 3) c V-SHAPED noS=2 q(3,1)=0.0 q(3, 2)=-d adj=0.75*maxC else if(igap.eq. 4) c T-SHAPED noS=5 dStp (l)=d*drat dStp (2) =d* (1-drat q(3,1)=w/2 q (3, 2) =-dStp (1) q (4, 1) =w/2+wStp (2 q(4,2)=-dStp(1) q(5, 1)=q(4, 1) q (5,2)=q(4, 2)-dSt q(6,1)=-w/2 q(6,2)=q(5,2) endif *... Corner points of gap at q(1,1)=-w/2 q(1,2)=0.0 q(2,1)=w/2 q(2,2)=0.0 q (noS+2, 1) =-w/2 q (noS+2, 2) =0.0 of corner points given gap type c******************************** ***** *** ****~***** ***~ c************* Current segment locations (xi,yi) * c***~****~* * *~**~*~****** N=0 do 175 l=1,noS+1 c...Size (length) of ith side of gap szSd (1) =sqrt ((q (1+1,1) -q (1, 1)) **2 & +(q(l+1,2) q (l,2))**2) c..Angle of rotation for each side with respect to x axis if(q(1+1,1).lt. q(1,1))then psi (1) =asin ( (q (1,2) -q (1+1,2)) IszSd (1)) neg(l)=.true. else psi (1) =asin ( (q (1+1,2) -q (1, 2)) /szSd (1)) neg (1)=. false. endif szN (1) =int (s zSd (1) /maxC) +1 N=N+szN (1) spaceX= (q (1+1, 1) -q (1, 1) ) /szN (1) spaceY= (q (1+1, 2) -q (1, 2) ) /szN (1) posiX=q(1,1) posiY=q(1,2) 41

c...ENDPOINTS of each segment are p, MIDPOINTS are m in (x,y) c coordinates do 170 i=N-szN(l)+l,N p(i,l)=posiX p(i,2)=posiY m(i, 1)=posiX+spaceX/2.0 m (i,2)=posiY+spaceY/2.0 posiX=posiX+spaceX posiY=posiY+spaceY 170 continue 175 continue p (N+1, 1)=-w/2 p (N+1, 2)=0.0 c...Number of segments in the aperture gN=szN (1) c...Number of current coefficients to be calculated NgN=N+gN print *, ' d = ',d,' N = ',N,' gN = ',gN c...Initializing matrices to zero do 190 j=1,NgN do 180 i=l,NgN Z (j,i)=czero 180 continue Vj(j)=czero 190 continue c****************************************************************** c******************* Impedance, Source, ********************** c****************** and Current Matrices ********************** c*************************************************************** C C VARIABLES: C HOo,HO Hankel function of zero order in free C space and in material er, respectively. C Hlo,H1 Hankel function of first order in free C space and in material er, respectively. C iHOo, Evaluation of the integral expression of C iHlo the half space Green's function. C Green's Function integrals: C LHO Integral of HO. C LH1 Integral of H1 for H-pol, of dHl/dy for C E-pol. C LHOo Integral of HOo and iHOo for H-pol and C iHOo for E-pol. C LHlo Integral of iHlo for E-pol. C aLHO Analytical integral of HOo for evaluation C of adjacent cells for LH2, E-pol case. C C The integration is done one side at a time, for j=1,...,N, C in the clockwise direction, starting at (x,y)=(-w/2,0). C istart=l istop=szN (1) c...Source point is i of the lth side, observation point is j do 230 l=l,noS+l do 220 i=istart,istop do 210 j=l,N c Coordinate rotation for observation point sj=m(j,l)*cos(psi(l))+m(j,2)*sin(psi(l)) nj=m(j, 1) *sin(psi (1)) -m(j, 2) *cos (psi (1)) if(neg(l))then sj=-sj nj=-nj endif c...Integration over ith segment LHO=czero LHl=czero LHOo=czero LHlo=czero c Magnitude between midpoints Rm=lr-r' Rm=sqrt ((m(j, 1)-m (i, 1)) **2+ (m (j, 2)-m(i, 2) )**2) if(j.eq. i.or. Rm.le. adj)then c SMALL ARGUMENT APPROXIMATION integral for self-cell c and adjacent cells do 200 ip=i+l,i,-l c Coordinate rotation for source segment points si=p (ip, 1) *cos (psi (1)) +p (ip, 2) *sin (psi (1)) ni=p (ip, 1)*sin(psi(l))-p(ip, 2) *cos (psi (1)) if(neg(l))then si=-si ni=-ni endif R=sqrt ((sj-si) **2+ (nj-ni) **2) if(j.eq. i.or. abs(nj-ni).eq. 0.0)then tanf=pi/2 absf=l.0 else tanf=atan((si-sj)/abs(nj-ni)) absf=abs(nj-ni) endif LHl=-kl**2/2*(nj-ni)*si & +ci*2./pi*(nj-ni)/absf*tanf-LH1 LHO=ci/pi*(2* (si-sj)*log(R) & -(2-A)*si+2.0*abs(nj-ni)*tanf)-LHO 42

LHOo=ci/pi* (2* (si-sj) *log(R) & -(2-Ao)*si+2.0*abs(nj-ni)*tanf)-LHOo if(j.eq. i) GOTO 202 200 continue 202 if(j.eq. i)then LHO=2* (LHO+ci/pi* (2-A) *sj) if(i.le. gN.and. j.le. gN)then if(i.eq. 1.and. iter.eq. 1)then iHlo=abs (sj-si) call simpGF(0.,etab,Epol,iHOo,iHlo) cselfO=iHOo cselfl=iHlo endif if (Epol)then LHlo=cselfO+cselfl else LHOo=2*(LHOo+ci/pi*(2-Ao)*sj)-cselfO endif endif endif else c SIMPSON'S THREE POINT COMPOSITE INTEGRATION do 204 ip=i+l,i,-l c Coordinate rotation for source segment endpoints si=p (ip, 1) *cos (psi (1)) +p (ip, 2) *sin (psi(1)) ni=p(ip,l)*sin(psi(l))-p (ip,2)*cos(psi (1)) if(neg(1))then si=-si ni=-ni endif stepS=si c HANKEL FUNCTION evaluation at endpoints of segment R=sqrt ((sj-si) **2+ (nj-ni) **2) if (Lossy) then ckrho=kl*R call cHank(ckrho,2,HO, Hl) else krho=Real (kl) *R call Hankzl(krho,2,HO,H1) endif if(i.le. gN.and. j.le. gN)then krho=k*R call Hankzl(krho,2,HOo,Hlo) if(i.eq. 1.and. iter.eq. 1)then call impGF(krho,etab,Epol,ctempO,ctempl) if(i.eq. ip)then GFOa(j)=ctempO iHOo=GFOa(j) GFla (j)=ctempl iHlo=GFla(j) else GFOb (j)=ctemp0 iHOo=GFOb(j) GFlb(j)=ctempl iHlo=GFlb(j) endif else if(i.eq. ip)then iHOo=GFOa(abs (-i)+1) iHlo=GFla (abs (j-i) +1) else iHOo=GFOb (abs (j-i) +1) iHlo=GFlb (abs (j-i) +1) endif endif endif LHO=HO+LHO if (Epol) then LHl=kl*m(j, 2)/R*Hl+LH1 if(i.le. gN.and. j.le. gN)then LHlo=-LHlo-iHlo LHOo=iHOo+LHOo endif else LH1=kl* (nj-ni)/R*Hl+LH1 LHOo=(HOo-iHOo)+LHOo endif 204 continue c Coordinate rotation for source segment midpoints si=m (i,1) *cos (psi (1)) +m (i,2) *sin (psi (1)) ni=m (i, 1) *sin (psi (1)) -m (i, 2) *cos (psi (1)) if(neg(l))then si=-si ni=-ni endif stepS=abs(stepS-si) DelS=2*stepS c HANKEL FUNCTION evaluation at midpoint of segment R=sqrt ((sj-si)**2+(nj-ni)**2) if (Lossy) then ckrho=kl*R call cHank(ckrho,2,HO,H1) else krho=Real(kl)*R 43

call Hankzl(krho,2,HO,Hl) endif lf(i.le. gN.and. j.le. gN)then krho=k*R call Hankzl (krho,2,HOo,Hlo) if(i.eq. 1.and. Iter.eq. l)then call impGF(krho,etab,Epol,ctempO,ctempl) GFO (j) =ctemp0 iHOo=GFO(j) GF1 (j)=ctempl lHlo=GFl(j) else iHOo=GFO (abs (j-l) +1) iHlo=GF1 (abs (j-i)+l) endif endif c GREEN' S FUNCTION INTEGRALS LHO=stepS/3* (4*HO+LHO) if(Epol)then LH1= stepS/3* (4*kl*m(j,2) /R*Hl+LH1) if(i.le. gN.and. j.le. gN)then LHOo=stepS/3*(4*iHOo+LHOo) LHlo=LHOo+LHlo endif else LH1=stepS/3* (4*kl* (nj-ni) /R*H1+LH1) if(i.le. gN.and. j.le. gN)then LHOo=stepS/3*(4*(HOo-iHOo)+LHOo) endif endif endif if(i.ne. j.and. Rm.le. adj)then LHO=-LHO endif if (Epol) then c E-POL IMPEDANCE MATRIX Z (j, i) =k*Zo*ur/2*LHO if(i.le. gN.and. j.ne. i)then Z (j, N+i) =-ci/2*LH1 else if(i.le. gN.and. j.eq. i)then Z(j,N+i)=-l. endif if(j.le. gN)then if(j.ne. i)then if(i.le. gN)then Z(N+j,i)=-k/2*etab*LHlo else Z (N+j, i) =czero endif else Z(N+j,i)=-l.-k/2*etab*LHlo endif if(i.le. gN)then Z(N+j,N+i)=-k*Yo/2*LHlo endif endif else C H-POL IMPEDANCE MATRIX if(j.ne. i)then Z(j, i)=ci/2*LH1 else Z (j,l)=-1. endif if(i.le. gN)then Z(j,N+i)=k*Yo*er/2*LHO endif if(j.le. gN)then if(j.ne. i)then if(i.le. gN)then Z (N+j, i) =-k*etab/2*LHOo else Z (N+j, i) =czero endif else Z (N+j,i)=1.-k*etab/2*LHOo endif if(i.le. gN)then Z (N+j,N+i)=k*Yo/2*LHOo endif endif endif 210 continue c...Incident Field (Source) matrix elements xj=m ( 1,) Vj(i)=czero if(i.le. gN)then if(Epol)then c E-POL INCIDENT FIELD Hx Vj (N+i ) =2*sin (phio) / (etab*sin (phio) +1) & *Yo*cexp(-ci*k*xj*cos(phio)) else c H-POL INCIDENT FIELD Hz 44

Vj(N+i) =2*sin (phio) / (etab+sin(phio)) & *cexp(-ci*k*xj*cos(phio)) endif endif 220 continue istart=istop+l istop=istop+szN (1+1) 230 continue c...Calling subroutines to calculate the current matrix call CGECO(Z,pn,NgN,ipvt,rc,wk) call CGESL(Z,pn,NgN,ipvt,Vj,0) print *,' The condition number is ',rc do 310 i=l,NgN c CURRENT MATRIX Ii(i)=Vj(i) 310 continue c************************************************************************ c****************** Far Field Amplitude *********************** c***************************************************************** Psca=czero DelX=w/gN/2 do 600 i=l,gN c...Simpson's three point composite integration over each c segment in the aperture Lsca=DelX/3*(cexp(-ci*k*p(i,1)*cos(phi)) & +4*cexp (-ci*k*m (i, 1) *cos (phi)) & +cexp(-ci*k*p(i+1, 1)*cos(phi))) if (Epol) then Psca- (Ii (N+i) +etab*Zo*Ii (i)) *Lsca+Psca else Psca= (Yo*Ii (N+i) -etab*Ii (i)) *Lsca+Psca endif 600 continue if (Epol) then Psca=-k/2*sin(phi) / (etab+l)*Psca else Psca=-k/2/(etab+l)*Psca endif write(3,*) d,cabs(Psca) write(4,*) d,180/pi* (atan2 (aImag (Psca),Real (Psca)) ) print *,' Exact: IPscal = ',cabs(Psca),' arg Psca = & 180/pi*(atan2(aImag(Psca),Real(Psca)) ) d=d+dstep 700 continue 800 call exit END C******************************************************************* SUBROUTINE IMPGF(KRHO,ETAB,EPOL,LRE,LRM) C****************************************** real pi,krho,k,nu,nuo,numax real*8 rts(64,64),coef(64,64) complex ci,carg,croot,Ao complex re,rm,eint,mint,Lre,Lrm, HO,H1 logical Epol,self,second common /data/ rts,coef call gausq pi=3.141593 ci=cmplx(0.,1.) k=2*pi gam=0.5772157 Ao=l+ci*2/pi*(gam+alog(k/2)) self=. false. second=. false. iroot=64 ih=0 alf=k/etab bet=k*etab nuo=krho/k c...Determining size of del nu and max nu value necessary if (Epol) then dnu=alog (0.95) / (-alf) numax=alog (0. 00001) / (-alf) ih=2 else dnu=alog (0. 95) / (-bet) numax=alog(0.00001 ) / (-bet) endif dnumax=sqrt ((0. 1/k) **2+nuo**2) -nuo if(dnu.gt. dnumax) dnu=dnumax rnumax=sqrt ( (12./k)**2+nuo**2) if(nuo.eq. 0.0)then self=. true. second=.true. endif 45

if(dnu.ge. nuo.and..not.(self))then dnu=nuo second=.true. endif if(numax.le. nuo+dnu) numax=2* (nuo+dnu) if(rnumax.lt. numax) numax=rnumax c...Gauss-Quadrature Lre=O.O Lrn=O.O do 20 m=1,2 if(.not.(m.eq. 1.and. second))then if(m.eg. 1)then a=O.0 b=nuo-dnu else if(m.eq. 2)then a=nuo+dnu b=numax else a=b b=numax endif mint=0.0 eint=0.0 do 70 i=l,iroot srts=rts (iroot, i) scoef=coef (iroot, i) nu= ((b-a) *srts+b+a) /2 croot=nuo**2-nu**2 carg=k*csqrt (croot) if(m.eq. l)then rarg=Real(carg) call Hankzl (rarg, ih, HO, H1) else call cHank(carg,ih,HO,H1) endif if (Epol) then rm=alf*exp(-alf*nu)*nuo*Hl/carg re=alf*exp(-alf*nu)*HO else re=bet*exp(-bet*nu)*HO endif eint=scoef*re+eint mint=scoef*rm~mint 70 continue Lre= (b-a) /2*eint+Lre Lrm (b-a)/2*mint+Lrm endif 20 continue c...Singularity evaluation, nuo - Ixj-xil if(.not.(self))then if (Epol)then mint=alf *exp (-alf*nuo) & *(nuo*dnu+ci/pi/k**2 & * (slog((2*nuo-dnu)/ (2*nuo+dnu)) +ci*pi)) eint alf*exp (=alf*nuo) & *(2*Ao*dnu & +ci/pi*((nuo+dnu)*(alog(2*dnu*nuo+dnu**2) & +ci*pi) - (nuo-dnu) *alog (2*dnu*nuo-dnu**2) & -4*dnu+nuo* (alog( (2*nuo+dnu)/ (2*nuo-dnu)) & -ci*pi))) else eint=bet*exp (-bet*nuo) & *(2*Ao*dnu & +ci/pi*((nuo+dnu)*(alog(2*dnu*nuo+dnu**2) & +ci*pi)-(nuo-dnu)*alog(2*dnu*nuo-dnui*2) & -4*dnu+nuo* (alog( (2*nuo~dnu)/ (2*nuo-dnu)) & -ci*pi))) endif endif Lre=Lre+eint Lrr=Lrm+mint return end c**************************************~*************** SUBROUTINE SIMPGF(KRHO,ETAB,EPOL,LRE,LRM) C*** * **** *************~******~**** **************** real pi,k,nu,nuo,numax,numin,krho real*8 rts (64, 64),coef (64, 64) complex cicarg,croot,Ao,H0,H1,H(64) complex re,rm,eint,mint,Lre,Lrm,ssLre,ssLrm logical Epol,self,second common /data/ rts,coef call gausq pi=3.141593 ci=cmplx(0.,1.) k=2*pi gam=0.5772157 Ao=l+ci*2/pi*(gam+alog(k/2)) self=.false. second=. false. iroot=64 46

numi=21 ih=0 alf =k/etab bet=k*etab nuo=krho/k c... Determining size of del fu and max fu value necessary if (Epol) then dnu=alog(0. 95)/(-aif) numax~alog(0.00001)/(-alf) ih=2 else dnu=alog(0. 95) / (-bet) numax=alog(0.0000 1)/ (-bet) endif dnumax=sqrt((0.1/k)**2+nuo**2)-nuo if(dnu.gt. dnumax) dnu=dnumax rnumax=sqrt ((12./k)**2+nuo**2) if(nuo.eq. 0.0)then self=.true. second=.true. dx=real (Lrm) if(dnu.ge. dx) dnu=1.05*dx/2 numin=dnu endif if(dnu.ge. nuo.and..not.(self))then dnu=nuo second=.true. endif if(numax.le. nuo+dnu) numax=2*(nuo+dnu) if(rnumax.lt. numax) numax=rnumax c...Self-cell evaluation if(self)then if(Epol)then mint=alf*exp ( alf*dnu/2) *2* (-dx/2*dnu+ci/pi/k**2 & *alog((dx+dnu)/(dx-dnu))) eint=alf*exp (-alf*dnu/2) *2* ((Ao-ci/pi*2) *dx*dnu & +ci/pi*dx*(dnu*alog(dx**2-dnu**2)-2*dnu & +dx*alog( (dx+dnu)/ (dx-dnu))) & +ci/pi/2*((dx**2-dnu**2)* (alog( (dx-dnu)/ (dx+dnu)) & -ci*pi)+dx*dnu+ci*dx**2*pi)) else eint=bet*exp (-bet*dnu/2) *2* ((Ao-ci/pi*2) *dx*dnu & +ci/pi*dx*(dnu*alog(dx**2-dnu**2)22*dnu & +dx*slog ((dx+dnu)/ (dx-drlu))) & +ci/pi/2* ((dx**2-dnu**2)* (alog( (dx-dnu)/ (dx+dnu)) & -ci*pi)+dx*dnu+ci*dx**2*pi)) endif ssLre=eint ssLrm-mint endif rt=0.0 13 if(rt.eq. 1.3) numm=numm+2 c...Integration over self-cell do 30 mm=l,numm if(mm.gt. 1)then nuo=nuo+dx/ (numm-1) if(abs(numin-nuo).le. 0.000001)then rt=1.3 nuo=0 print *,'ACK!!!!!! Changing numim.' goto 13 endif if (Epol) then dnu=alog (0. 95) / (-alf) else dnu-alog(0.95)/(-bet) endif dnumax=sqrt((0.1/k)**2+nuo**2)-nuo if(dnu.gt. dnumax) dnu=dnumax if(nuo.lt. numin)then second=.true. self=.true. else if(dnu.ge. (nuo-numin))then dnu-nuo-numin second=.true. endif endif if(numax.le. nuo+dnu) numax=2*(nuo+dnu) if(rnumax.lt. numax) numax=rnumax endif c Gauss-Quadrature Lre=0.0 Lrm=0.0 do 20 m=1,2 if(.not.(m.eq. 1.and. second))then if(m.eq. l)then a=numin b=nuo-dnu else if(m.eq. 2)then 47

a=nuo+dnu if(nuo.lt. numin) a=numin b=numax else a=numax b=2*numax endif mint =0.0 eint=0.0 do 70 i=l,iroot srts=rts (iroot, i) scoef=coef (iroot, i) nu= ((b-a) *srts+b+a)/2 croot= ((nuo) **2-nu**2) carg=k*csqrt (croot) if(m.eq. 1)then rarg=real (carg) call Hankzl (rarg, ih, HO, H1) else call cHank (carg, ih,HO,H1) endif if (Epol) then rm=alf*exp (-alf*nu) *2*nuo*H1/carg re=alf*exp (-alf*nu) *HO else re=bet*exp (-bet*nu) *H0 endif eint =scoef* re+eint mint =s coe f* rm+mint 70 continue Lre= (b-a)/2*eint+Lre Lrm= (b-a) /2*mint+Lrm endif 20 continue c...Singularity evaluation, nuo = Ix-xij if(.not. (self))then if (Epol) then mint =alf*exp (-alf*nuo) & *2* (nuo*dnu+ci/pi/k**2 & *(alog((2*nuo-dnu)/ (2*nuo+dnu))+ci*pi)) eint=alf*exp (-alf*nuo) & *(2*Ao*dnu & +ci/pi* ((nuo+dnu) * (alog (2*dnu*nuo+dnu**2) & +ci*pi)- (nuo-dnu) *alog (2*dnu*nuo-dnu**2) & -4*dnu+nuo* (alog((2*nuo+dnu) / (2*nuo-dnu)) & -ci*pi))) else eint=bet*exp (-bet*nuo) & *(2*Ao*dnu & +ci/pi* ((nuo+dnu) * (alog (2*dnu*nuo+dnu**2) & +ci*pi)- (nuo-dnu) *alog (2*dnu*nuo-dnu**2) & -4*dnu+nuo* (alog((2*nuo+dnu) / (2*nuo-dnu)) & -ci*pi))) endif Lre=Lre+eint Lrm=Lrm+mint endif H (mm)=Lre self=. false. second=. false. 30 continue Lre=O.O do 40 i=1,numm-2,2 Lre=H (i) +4*H (i+1) +H (i+2) +Lre 40 continue Lre=ssLre+2* (dx/(numm-1)/3*Lre) Lrm=ssLrm-Lrm return end C SUBPROGRAM DATAINT CONTAINS INTEGRATION DATA C SUBROUTINE GAUSQ REAL*8 RTS (64, 64), COEF (64, 64) C COMMON /DATA/ RTS,COEF C C C --- —-------------------------------------------— + C FIXED POINTS FOR GAUSSIAN QUADRATURE C+ --- —------------------------------------— + C C ---N=64 C RTS (64,1) =.999305041735772D0 RTS(64,2)=.996340116771955D0 RTS(64,3)=.991013371476744D0 RTS(64,4)=.983336253884625D0 RTS(64,5)=.973326827789910D0 RTS (64, 6)=. 961008799652053D0 RTS(64,7)=.946411374858402D0 RTS(64,8)=.929569172131939D0 48

RTS(64, 9). 91052213707850200 RTS(64,10)=.88931544599511400 RTS(64,11)=.86599939815409200 RTS(64,12)=.84062929625258000 RTS(64,13)=.81326531512279700 RTS(64,14)=.78397235894334100 RTS(64,15)=.75281990726053100 RTS(64,16)=.719881850171610D0 RTS(64,17)=.68523631305423300 RTS(64,18)=.648965471254657DO RTS(64,19)=.61115535517239300 RTS(64,20)=.571895646202634D0 RTS(64,21)=.531279464019894D0 RTS(64,22)=.489403145707052D0 RTS(64,23)=.446366017253464D0 RTS(64,24)=.402270157963991D0 RTS(64,25)=.357220158337668D0 RTS(64,26)=.311322871990210D0 RTS(64,27)=.264687162208767D0 RTS(64,28)=.217423643740007Do RTS(64,29)=.169644420423992D0 RTS(64,30)=.12146281929612000 RTS(64,31)=.07299312178779900 RTS(64,32)=.02435029266342400 RTS(64,64)=-.999305041735772D0 RTS (64, 63)=-.996340116771955D0 RTS(64,62)=-.99101337147674400 RTS(64,61)=-.983336253884625D0 RTS(64,60)=-.973326827789910D0 RTS(64,59)=-.96100879965205300 RTS (64,58)=-. 946411374858402D0 RTS(64,57)=-.92956917213193900 RTS(64,56)=-.910522137078502D0 RTS(64,55)=-.889315445995114D0 RTS(64,54)=-.865999398154092D0 RTS(64,53)=-.840629296252580D0 RTS(64,52)=-.813265315122797D0 RTS(64,51)=-.783972358943341D0 RTS(64,50)=-.752819907260531D0 RTS(64,49)=-.71988185017161000 RTS (64, 48)=-.685236313054233D0 RTS(64,47)=-.648965471254657D0 RTS (64,46)=-.611155355172393D0 RTS(64,45)=-.571895646202634D0 RTS(64,44)=-.53127946401989400 RTS(64,43)=-.48940314570705200 RTS(64,42)=-.446366017253464D0 RTS(64,41)=-.402270157963991D0 RTS(64,40)=-.357220158337668D0 RTS(64,39)=-.311322871990210D0 RTS(64,38)=-.264687162208767D0 RTS(64,37)=-.217423643740007D0 RTS(64,36)=-.169644420423992D0 RTS(64,35)=-.121462819296120DO RTS(64,34)=-.07299312178779900 RTS(64,33)=-.02435029266342400 COEF(64,1)=.00178328072169600 COEF(64,2)=.004147033260562D0 COEF(64,3)=.006504457968978D0 COEF(64,4)=.008846759826363D0 COEF(64,5)=.01116813946013100 COEF(64,6)=.01346304789671800 COEF(64,7)=.01572603047602400 COEF(64,8)=.01795171577569700 COEF(64,9)=.02013482315353000 COEF(64,10)=.02227017380838300 COEF(64,11)=.02435270256871000 COEF(64,12)=.02637746971505400 COEF(64,13)=.02833967261425900 COEF(64,14)=.030234657072402D0 COEF(64,15)=.03205792835485100 COEF(64,16)=.033805161837141D0 COEF(64,17)=.03547221325688200 COEF(64,18)=.03705512854024000 COEF (64,19)=.038550153178615DO COEF(64,20)=.039953741132720D0 COEF(64,21)=.041262563242623D0 COEF(64,22)=.042473515123653D0 COEF(64,23)=.043583724529323D0 COEF(64,24)=.04459055816375600 COEF(64,25)=.04549162792741800 COEF(64,26)=.04628479658131400 COEF(64,27)=.04696818281621000 COEF(64,28)=.04754016571483000 COEF(64,29)=.04799938859645800 COEF(64,30)=.04834476223480200 COEF (64,31)=. 048575467441503D0 COEF(64,32)=.048690957009139D0 COEF(64,64)=.001783280721696D0 COEF(64,63)=.00414703326056200 COEF(64,62)=.00650445796897800 COEF(64,61)=.00884675982636300 COEF(64,60)=.011168139460131D0 C C C 49

COEF (64,59)= COEF(64,58)= COEF(64,57)= COEF (64, 56) = COEF(64,55)= COEF(64,54)= COEF(64,53)= COEF(64,52)= COEF(64,51)= COEF(64,50)= COEF(64,49)= COEF(64,48)= COEF(64,47)= COEF (64,46)= COEF(64,45)= COEF(64,44)= COEF(64,43)= COEF(64,42)= COEF (64,41)= COEF(64,40)= COEF (64,39)= COEF(64,38)= COEF(64,37)= COEF (64,36)= COEF(64,35)= COEF(64,34)= COEF(64,33)=.013463047896718D0.015726030476024D0.017951715775697D0.02013482315353000.02227017380838300.024352702568710D0.026377469715054D0.02833967261425900.03023465707240200.03205792835485100.033805161837141D0.03547221325688200.03705512854024000.03855015317861500.039953741132720D0 0 41262563242623D0.04247351512365300.043583724529323D0.04459055816375600.045491627927418D0.04628479658131400.04696818281621000.04754016571483000.04799938859645800.04834476223480200.04857546744150300.04869095700913900 C C ---N=64 C RTS(64,1)=.999305041735772D0 RTS(64,2)=.996340116771955D0 RTS(64,3)=.991013371476744D0 RTS(64,4)=.98333625388462500 RTS(64,5)=. 973326827789910D0 RTS (64,6)=.961008799652053D0 RTS(64,7)=.946411374858402D0 RTS(64,8)=.929569172131939D0 RTS(64,9)=.910522137078502DO RTS(64,10)=.88931544599511400 RTS(64,11)=.86599939815409200 RTS(64,12)=.84062929625258000 RTS(64,13)=.81326531512279700 RTS(64,14)=.78397235894334100 RTS(64,15)=.752819907260531D0 RTS(64,16)=.71988185017161000 RTS(64,17)=.685236313054233D0 RTS(64,18)=.648965471254657D0 RTS(64,19)=.61115535517239300 RTS(64,20)=.57189564620263400 RTS(64,21)=.53127946401989400 RTS(64,22)=.489403145707052D0 RTS(64,23)=.44636601725346400 RTS(64,24)=.40227015796399100 RTS(64,25)=.35722015833766800 RTS(64,26)=.31132287199021000 RTS(64,27)=.26468716220876700 RTS(64,28)=.21742364374000700 RTS(64,29)=.16964442042399200 RTS(64,30)=.12146281929612000 RTS(64,31)=.072993121787799D0 RTS(64,32)=.02435029266342400 RTS(64,64)=-.999305041735772D0 RTS(64,63)=-.99634011677195500 RTS(64,62)=-.991013371476744D0 RTS(64,61)=-.983336253884625D0 RTS (64, 60)=-. 97332682778991000 RTS(64,59)=-.96100879965205300 RTS(64,58)=-.946411374858402D0 RTS(64,57)=-.929569172131939D0 RTS(64,56)=-.910522137078502DO RTS(64,55)=-.88931544599511400 RTS(64,54)=-.86599939815409200 RTS(64,53)=-.84062929625258000 RTS(64,52)=-.81326531512279700 RTS(64,51)=-.78397235894334100 RTS(64,50)=-.75281990726053100 RTS(64,49)=-.719881850171610D0 RTS (64, 48) =-.685236313054233D0 RTS(64,47)=-.648965471254657D0 RTS(64,46)=-.611155355172393D0 RTS(64,45)=-.571895646202634D0 RTS(64,44)=-.53127946401989400 RTS(64,43)=-.48940314570705200 RTS(64,42)=-.446366017253464D0 RTS(64,41)=-.402270157963991D0 RTS(64,40)=-.357220158337668D0 RTS(64,39)=-.31132287199021000 RTS(64,38)=-.26468716220876700 RTS(64,37)=-.217423643740007D0 RTS(64,36)=-.1696444204 23992D0 RTS(64,35)=-.12146281929612000 RTS(64,34)=-.072993121787799D0 RTS(64,33)=-.024350292663424D0 C C 50

c COEF (64,1)= COEF(64,2)= COEF(64,3)= COEF(64,4)= COEF(64,5)= COEF(64,6)= COEF(64,7)= COEF(64,8)= COEF(64,9)= COEF(64,10)= COEF(64,11)= COEF(64,12)= COEF(64,13)= COEF (64,14)= COEF(64,15)= COEF(64,16)= COEF (64,17)= COEF(64,18)= COEF(64,19)= COEF(64,20)= COEF(64,21)= COEF(64,22)= COEF(64,23)= COEF(64,24)= COEF(64,25)= COEF(64,26)= COEF(64,27)= COEF (64,28)= COEF(64,29)= COEF(64,30)= COEF(64,31)= COEF(64,32)= COEF (64, 64)= COEF(64,63)= COEF(64,62)= COEF(64,61)= COEF(64,60) = COEF(64,59)= COEF(64,58)= COEF(64,57)= COEF(64,56)= COEF(64,55)= COEF(64,54)= COEF(64,53)= COEF(64,52)= COEF(64,51)= COEF(64,50)= COEF(64,49) = COEF(64,48)= COEF(64,47)= COEF(64,46)= COEF(64,45)= COEF(64,44)= COEF(64,43)= COEF(64,42)= COEF(64,41)= COEF(64,40)= COEF(64,39)= COEF (64,38) = COEF(64,37)= COEF(64,36)= COEF(64, 35)= COEF(64,34)= COEF(64,33)= RETURN END.001783280721696D0.004147033260562D0.006504457968978D0.008846759826363D0.011168139460131D0.013463047896718D0.015726030476024D0.017951715775697DO.020134823153530D0.0222701738088383D0.024352702568710D0.026377469715054D0.028339672614259D0.030234657072402D0.032057928354851D0.033805161837141D0.035472213256882D0.037055128540240D0.038550153178615D0.039953741132720DO.041262563242623D0.042473515123653D0.043583724529323D0.044590558163756D0.045491627927418D0.046284796581314DO.046968182816210D0.047540165714830D0.047999388596458D0.048344762234802D0.048575467441503D0.048690957009139D0.001783280721696D0.004147033260562D0.006504457968978D0.008846759826363D0.011168139460131D0.013463047896718D0.015726030476024D0.017951715775697D0.020134823153530D0.0222701738088383D0.024352702568710D0.026377469715054D0.028339672614259D0.030234657072402D0.032057928354851D0.033805161837141D0.035472213256882D0.037055128540240D0.038550153178615D0.039953741132720D0.041262563242623D0.042473515123653D0.043583724529323DO.044590558163756D0.045491627927418D0.046284796581314D0.046968182816210D0.047540165714830D0.047999388596458D0.048344762234802D0.048575467441503D0.048690957009139D0 C 51

APPENDIX C. Program Listing for the Quasi-Analytical Solution The quasi-analytical solution presented in Section 3 was programmed for solution, as listed in the program IMPQA.FTN below. The subroutines used by this program are listed in GAPSUB.FTN in Appendix A of [2]. For H-polarization, the far field amplitude PH is calculated from (28), and for E-polarization, PE is calculated from (34). The parameters a and b in (26) and (33), respectively, are dependent on the effective surface impedance. of the gap, which is calculated from the formulas in [2] according to the specified shape of the cavity and the field polarization. 52

c C IMPQA.FTN C C C This FORTRAN program computes the far field scattering due C to a narrow gap of specified shape in an infinite impedance C plane. The far field amplitude is calculated given the C input impedance of the gap, calculated from it's equivalent C transmission line model. C C INPUT The user is prompted from the subroutine C GAPROM for the polarization and angle of the C incident field, angle of far field observation, C relative permittivity of gap filling, shape and C dimensions of gap, segment size, number of C iterations with respect to gap depth, and C normalized input impedance. C C OUTPUT FILES C GAPDAT Contains input data. C IMPDAT Effective surface impedance of the gap. C AMPDAT Contains the magnitude of the far field. C PHADAT Contains the phase of the far field. C C SUBROUTINES C HANKZ1 Computes the Hankel functions of the first C kind of orders zero and one. C CHANK Computes the Hankel functions of the first C kind of orders zero and one given a C complex argument. C MODBES Computes the modified Bessel functions of C the first kind of orders zero and one. C C FUNCTION C CTAN Calculates the tangent of a complex C argument. C C integer EorH,N,noS,gN,szN(50) real pi,k,phi,phio,w,d,maxC, q (50,2) real dStp(50),wStp(50), I0,I1,krho,psi complex czero, ci, ctemp,er,ur, kl, kc, carg, ctan, eta complex Zi, Z1,Zc,ZL,X,B1,B2,HO,H1,a,b,Ke,Kh,Ao,Psca,ge,gh logical Epol,Lossy common /prompt/ EorH,phio,phi,er,ur,igap,wStp,dStp,w,d,noS, & q, maxC, noIter, etab 1 format (il) 2 format (i5) 3 format (gl6.8) 4 format (al) 5 format (i3) 6 format (2gl6.8) open (1, file='gapdat') open (2, file=' impdat' ) open (3, file='ampdat' ) open (4, file='phadat') c...Declaring constant values czero=cmplx (0.0,0.0) ci=cmplx (0.0,1.) pi=4.0*atan (1.0) k=2*pi Eo=1.0 Ho=1.0 Zo=sqrt (4.e-07*pi/8.854e-12) Yo=1./Zo gam=0. 5772157 iprg=2 c... Setting default values 10 EorH=l phio=90.0 phi=90.0 er=cmplx (1., 0.0) ur=cmplx (1., 0.0) w=0.15 d=0.5 noS=3 maxC=0.01 noIter=30 etab=1. adj=0. 00001 Epol=. false. Lossy=. false. side=. true. c...Prompting user for input data 15 call gaprom (iprg) if (igap.eq. 5) GOTO 15 if(EorH.eq. 1) Epol=.true. phio=phio*pi/180.0 phi=phi*pi/180.0 drat=dStp (1) /d 53

if(almag(er).ne. 0.0) Lossy=.true. dmin=0.025 if(Epol) dmin=0.02 dmax=d if(nolter.ne. ))then dstep= (dmax-dmin) / (nolter-1) d=dmin endif DO 700 iter=l,noIter ***** ********** * * ** ***Gap Impedance******************** * C* ***** * * t** *** * * Gp mpdace *** * * ***** * *** * * c...Complex propagation constant ki and characteristic impedance c Z1 of the T-line model if (Epol) then kl=ci*k*csqrt((l./2/w)**2-er*ur) Z1=-ci*Zo*ur/csqrt((1.121w)**2-er*ur) else Zl=Zo*csqrt (ur/er) kl=k*csgrt(ur*er) endif if(igap.eq. 1)then c RECTANGULAR ETA=-ci*Z1*ctan (kl*d) else if(igap.eq. 2.or. igap.eq. 4)then w=wStp (1) w2=wStp(2) w3=wStp(3) dl=d*drat d2=d*(l-drat) c Propagation constant and characteristic impedance of c the arms of the T- or L-shaped gaps if (Epol) then kc=ci*k*csqrt((l./2/d2)**2-er*ur) Zc=-ci*Zo*d2*ur/csqrt ((1./21/d2) **2-er*ur) else Zc=Zo*d2*csgrt (ur/er) kc=k*csqrt (ur*er) endif if(igap.eq. 4)then c L-SHAPED Zi —ci*Zc*ctan(kc*w2) X=kl *Zc*wl B1=kl/Zc*d2*(d2/(d2+wl))*(l. 2./pi*log(2.)) B2=kl/Zc*d2* (wl /(d2+wl) ) * (1.-2. /pi*log (2.)) ZL= (Zi-ci*X* (1 ci*B2*Zi)) & /((l-B1lX)*(l-ci*B2*Zi)-ci*B1*Zi) else c T-SHAPED Zi =ci*Zc* (ctan (kc*w2) +ctan (kc*w3)) X=kl*Zc*wl B1=kl/Zc*d2*(d2/(d2+w)) *0.7822 ZL=(Zi-ci*X)/(lci*B1* (Zi-ci*X)) endif ETA=Z1* (ZL-ci*Z1*wl*ctan (kl*dl)) & / (Zl*wl-ci*ZL*ctan (kl*dl)) else if(igap.eq. 3)then c TRIANGULAR if(Lossy)then carg=kl*d call c~ank(carg,3,HO,H1) ETA=-ci*Zl*H1/HO else if (Epol) then rarg=Real (kl/ci) *d call ModBes(rarg, IO,11) HO=IO H1=I1 carg=ci else rarg=Real (kl) *d call Hank zl (rarg, 2,HO,H1) carg=1. endif ETA=-ci*Z1*Real(Hi)/Real(HO)*carg endif endif c write(2,*) d,cabs(ETA) c~** ******** * * *** ****~****~**** ****** ** **** ****** * c****************** Far Field Amplitude ******************** c*** *********t~t****** Rsdp=l./ (etab+l) if (Epol) then theta= (-.380-0.8*w)* (l-exp(-2.586*sqrt (etab))) dw=1.558-4.266*w ge=exp(-dw*etab)*cexp(ci*theta) b=-ci*k*w/2*Zo/ (ETA-etab*Zo) /ge Ke=0.62/(b+1.15)*(b+4.08)*(b+7.26)*(b+10.37) & * (b+1 3. 43) * (b+1 6. 4 6) & /((b+4.27)*(b+7.37)*(b+10.45)*(b+13.49) & *(b+16.5)) 54

Ri=sin(phio)/(sin(phio)*etab+1) Psca=-ci*pi/4*(k*w)**2*sin(phi)*Rsdp/ge*Ri*Ke else psi=l-exp (-1 (0.098+1. 76*w) *etab) cw=0.245+1.267*w gh=exp (-cw*etab) *cexp (ci*psi) a=ci*2/k/w*Zo/(ETA-etab*Zo)/gh Kh=-l./(pi/2*a+0.1+log(2.)) Ao=log(k*w/4)+gam-ci*pi/2 Ri=sin (phio) / (etab+sin (phio)) Psca=ci*pi*Rsdp/gh*Ri*Kh/(l+Ao*Kh) endif c... Outputting the far field magnitude and phase print *,' d = ',d print, Analytical: IPscaj = ',cabs(Psca), & arg Psca = ',180/pi*(atan2(almag(Psca),Real(Psca))) write (3, *) d, cabs (Psca) write(4,*) d,180/pi*(atan2(almag(Psca),Real(Psca))) d=d+dstep 700 continue print *,' Again (l=yes)? read(*,l)iana if(ians.eq. 1) GOTO 10 800 call exit END 55