ENHN 389741-3-T UMH3631 THE UNIVERSITY OF MICHIGAN COLLEGE OF ENGINEERING DEPARTMENT OF ELECTRICAL ENGINEERING & COMPUTER SCIENCE $ Radiation Laboratory CD SCATTERING BY A NARROW GAP S: T.B.A. Senior, K. Sarabandi, and J.R. Natzke Radiation Laboratory Department of Electrical Engineering and Computer Science The University of Michigan Ann Arbor, MI 48109 April, 1989 McDonnell Aircraft Company St. Louis, MO 63166 Ann Arbor, Michigan

Abstract For a plane wave incident on a cavity-backed gap in a perfectly conducting 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-analytic solution previously derived, and for a narrow gap the agreement is excellent for all cavity geometries and for all material fillings that have been tested.

TABLE OF CONTENTS Page# I. INTRODUCTION 1 II. FORMULATION 2 2.1 H-Polarization 4 2.2 E-Polarization 8 III. QUASI-ANALYTICAL SOLUTION 11 IV. NUMERICAL RESULTS 19 V. CONCLUSIONS 36 REFERENCES 37 APPENDIX A: MOMENT METHOD SOLUTION OF THE COUPLED INTEGRAL EQUATIONS 38 A.1 H-Polarization 38 A.2 E-Polarization 43 A.3 Program Listings 45 APPENDIX B: PROGRAM LISTING FOR QUASI-ANALYTICAL 65 SOLUTION

1. Introduction A topic of some concern in radar cross section studies is the scattering from the gap or crack that may exist where two component parts of a target come together. Even if the crack is wholly or partially filled with a material, it can still provide a significant contribution to the overall scattering pattern of the target, and it is then necessary to develop methods for predicting the scattering. One method for doing this was described recently [1]. For a plane wave of either principal polarization incident on a narrow (kw << 1) resistive strip insert in an otherwise perfectly conducting plane, the low frequency approximations to the integral equations for the currents induced in the strip were solved in a quasi-analytic manner, leading to expressions for the far zone scattered field that are accurate for almost any resistivity R of the insert. If, instead, the insert is characterized by a surface impedance 1, the results differ only in having R replaced by ir/2 and the scattered field doubled, and this suggests that for a narrow gap backed by a cavity, the scattered field can be obtained by identifying Ti with the impedance looking into the cavity. An alternative approach is to use the equivalence principle [2] to develop coupled integral equations for the electric and magnetic currents which exist on the walls of the cavity and in the aperture, and this is the method employed here. For an incident plane wave either H- or E-polarized, 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 a variety of simple cavities are presented. For gap widths which are electrically small the results are compared with those 1

obtained using the previous method. The agreement is excellent, and confirms the utility of the original method [1] as an accurate and simple design tool. 2. Formulation The problem considered is the two-dimensional one shown in Figure 1. The plane y = 0 is perfectly conducting apart from the aperture A: - w/2 < x < w/2, which forms the entrance to a cavity whose walls S are also perfectly conducting. The cavity is filled with a homogeneous dielectric material of permittivity e1 = e~e and permeability 91 = PrWL where the quantities without subscripts refer to free space. A plane wave of either principal polarization is incident on the surface y = 0 from above, and we choose - i A -ik (x coso + y sino) H =ze (1) for H-polarization, and Ai -ik (x cos + y sin) E =ze (2) for E-polarization, where k is the propagation constant in the free space region above the surface. A time factor e-i(t is assumed and suppressed. In the far zone of the gap the scattered field can be written as - s A / 2 i (ikp- pV'4) H =z e PH ((4 OO) Xkp 2

y A / Fig. 1. Gap geometry. 3

for H-polarization, with a similar result for E-polarization, and the task is to determine the far field amplitudes PH,E (<, 4o) with particular emphasis on the case of a narrow gap (kw < 1). 2.1 H-Polarization We consider first the free space region y > 0. Using Green's theorem in conjunction with the half space Green's function G 4 {H k( 2 )2 +Ho (k J(x-x') + (y+y')2 + ) _-* A - the scattered field can be attributed to a magnetic current J = - y x E in the aperture, and the total magnetic field is then w/2 i r k Jz (x') H(1) r2 2 Hz(xy)= Hz (xy) y)+ (x) + (xx y dx' -w/2 where r -ik (x coso- y sino) Hz=e is the reflected plane wave and Y (= 1/Z) is the intrinsic admittance of free space. Hence, w/2 PH ( AI) - Jz (x')e dx'. (3) -w/2 4

and in the aperture w/2 -ikx coso kY f () (klx-x'l)dx'. Hz (x,o) = 2e - J (x') kx-x) dx (4) 2 z 0 (xx) d -w/2 We now turn to the region y < 0 occupied by the cavity. In accordance with the equivalence principle [2] it is assumed that the gap is closed with a perfect conductor, and that a magnetic current - J is placed just below, thereby ensuring the continuity of the tangential electric field in the open gap. The magnetic Hertz vector is therefore w/2 _*. Y r * (1) ( 2 1 (x,y)= 4kY ()xx') +yo dx' (5) - * A * and since J = z Jz, the magnetic field produced is w/2 (xy)=VxVx z J (x')H k +y-2dx where k1 = k J/err is the propagation constant. The electric current J = n x H on the cavity walls S and in the (closed) aperture A also implies an electric Hertz vector 11(xly) =- J(s,)H (k1l(X-X, +(y.y) ds', (6) 4 =kr S+A and the corresponding magnetic field is 5

- (2) H (x,y) =-ikY erV x n 4 H(k x-x') + (y)2 xJ (s-y ' ds' S+A A A A A where the tangential unit vector s is such that n, s, z form a right-handed system A A A A with n directed into the cavity. Clearly, J (s) = s Js (s), and in the aperture, s = x - (1) -(2) The total magnetic field is H H, and by allowing the observation point to approach the boundary of the closed cavity, we can construct an integral equation for the currents. We find J (s) = (n x z) kY J Jz (x') Ho (k /(x-x')2 + y2 dx' w/2 J (s) (n x z) J r Hz 1 (x') H o( (2 j dx' 4 1 y -W-2 lim ^ i (1).2 2 A + lm |Ss) n x I J (s') V H(x-x + (y-y')x sds (x,y)- S+-A S+A giving wr2 J (S)kY J (x') H k1 (x-x + y dx' -w/2 iki 2* ()? + -J1 Js (s') sin H k (xx + (y) ds (7) S+A where A A A (x-x') x + (Y-y') y x As siny'=z* xs', (8) (XX',) + (y-y)2 valid at all points of S and A. 6

The only remaining task is to enforce the continuity of Hz through the aperture. When the observation point is in the aperture w/2 kY f ~.~1, (~ Ix-x'l) dx' Hz (x,o) =9-r Jz (x') Ho (k Ix-x') dx' -w/2 lim A ( f 2 + Z 4 J (S')VHo )k1 (x-x') +y )xs'ds' y-40 O S+A and therefore w/2 kY h (1) 1 Hz (x,o) =4, r Jz (x') H0 (k1 Ix-x') dx' + (s) -w/2 +h Js ((s')sinY H (k1 (x-x') +y ds' (9) S+A When this is equated to the expression (4) for Hz (x,o) on the outside of the gap, we obtain w/2 w/2 -ikxcoso kY 1 2 e -= 2 H Jz (x') Ho (k Ix-x'l) dx' + -4 r Jz (x') Ho (k1 Ix-x') dx' 2 J -w/2 -w/2 1 2(1) + iJs (x) 4 Js (s') sin y' H 1 k (x-x)+y ds' (10) S+A valid for - w/2 < x < w/2. Since (7) is also valid in A, it can be used to simplify (10) by eliminating two of the integrals. The result is 7

w/2 -ikx cos% kY J (x) = 2e - Jz (x') H (k Ix-x') dx (1 -w/2 valid for x in A, and (7) and (11) constitute a pair of coupled integral equations for Jz (x) and Js (s). These are the equations that will be used, and we note the similarity of (11) and (4). When the maximum dimension of the cavity is electrically small, the Hankel function Ho can be replaced by its logarithmic approximation, and though this does not significantly simplify the numerical solution of (7) and (11), the fact that -ikx cos * e can also be replaced by unity shows that Jz (x) and Js (s) are aspect independent. The same approximation to (3) then leads to a far field amplitude which is independent of 4 and 0,0 and this is a feature of the low frequency situation. 2.2 E-Polarization The procedure is similar to that given above. For the region y > 0 Green's theorem in conjunction with the Green's function { (x-x)2 + (y-y,)2) - (k(x-x') + (y+y)2)} gives w/2 E=EEz +EZ + - x2 J J (x) Ho (k (x-x') + )dx' _., ~ A * where J - x Jx is the assumed magnetic current in the gap and 8

r -ik (x coso- y sinto) Ez=-e is the reflected plane wave. Hence w/2 k ( -ikx'coso PE (, o) = -- sin< Jx (x') e dx', (12) -w/2 iY aEz and since H = --, the tangential component of the magnetic field in the k ay aperture is w/2 Hx (xo)=- 2Y sin(o ekY 1 + 2 2 J x (x') Ho (k Ix-x'l) dx'. (13)' k a * In the region y < O0 occupied by the cavity, the field can be attributed to the - * _A A magnetic Hertz vector (5) with J = x Jx and the electric Hertz vector (6) with J = z Jz The magnetic field is therefore w/2 H (x,y)= VxV x- J (x') Ho ) k1 /(x-x') + y dx' 4k.r ir. + J H (x-xt) + (y-y) x J (s') ds', S+A and by allowing the observation point to approach the boundary of the closed cavity, we obtain the integral equation 9

w/2 Y A a r(1) 2\ I J, (s)= (ne V) - Jx (x') Ho k (x-x') + y dx' 2klLr ay -w/2 +"2- Jz (s') sinH (k1 (x-x) + (y-y')2 ds' (14) S+A where A A ^ (x-x') x + (y-y')y A siny=z. xs (15) J (x-x) + (y-y,)2 valid at all points of S and A. When the observation point is in the aperture, 1 o2 wM2 Hx (x,o)= kY e( JJx (x') H (k1 Ix-x'l) dx'+ -Jz (x) 1 o-w2 ik1lf (1), j2.. + 4 Jz (s') sin y H( k (x-x' ds' 4k (X-XZ)H + y S+A and on equating this to the expression (13) for the magnetic field on the outside of the gap, we obtain ~ 2 x2 \ w 42 -ikxcos% kY 1+ k1 a \2 2Y sino e = 2 a- Jx (x') H ) (k Ix-x'I) dx' (2 1-w/2 ~ kY(1 + 4 i I Jx (x') f (k, Ix-I) dx' ^ + JA (s) sin yH k ds (1 S+A 10

valid for x in A. This can be simplified using (14) and the result is (2 w/ -ikxcos% kY 1 a k2 Jz (x) =- 2Y sin(o e ~ - 1 + 2 J2 x (x') H() (k Ix-x'l) dx' (17) ' 2 oax -w/2 for x in A in accordance with (13), and (14) and (17) constitute a pair of coupled integral equations for Jx (x) and Jz (s). There is a third integral equation that can be developed and this has some advantages for numerical purposes. In the region y < 0 the electric field produced by the electric and magnetic Hertz vectors is kZr f (1) ----- E (x,y)= 4 J z (s') H k1 (xx) (yy,)2 ds z S+A w/2 4y J(x') () (k1 k (x-x')2 + y2 dx' z and when the boundary condition on the perfectly conducting surface is applied, we find -w/2 ay Jx (x) = i | Jx (x') ~a Hol k1 /(X-X'+y2 dxk kZLr 2 (1) ( ( ---2 - 2 (18) + 2 Jz (s') H k (x-x') +(y-y') ds (18) S+A valid on S + A. Of course, Jx (x) is non-zero only in A, and (17) and (18) are the pair of integral equations used to compute J, (x) and Jz (s). 11

3. Quasi-Analytical Solution An alternative approach was proposed by Senior and Volakis [1]. In effect, the problem which they considered is a uniform impedance insert in an otherwise perfectly conducting plane. If ir is the surface impedance, the integral equations for H- and E-polarizations are identical to (4) and (13) respectively, with H (x,o) =J, ( Hx,o) = - J (x) (19) at the insert. At low frequencies for which kw << 1 the integral equations can be simplified, and for H-polarization it is found that JJ2 () K'-cI d'= 1 +a J2 () (20) for-1 < < 1 with 2i Z a= 2. (21) kw J2 (() is a modified current in terms of which PH (. { o) = A + KH(a) (22) with 1 KH (a)= J2 (,d (23) I 1 and 12

kw e4 2 where y= 0.5772157.... is Euler's constant. We observe that PH (4, 4o) is independent of < and %O, and since KH (a) is real if a is, IPH ( 2o)| 2 (24) for real a. Similarly, for E-polarization the low frequency approximation to the integral equation is 1 a2 J3 (') n I' - CI dC'= 1 - b J3 () (25) a'2 'F -1 for-1 < C < 1 with ikw Z b =- k (26) 2 i: where the modified current J3 (0) is such that J3 (+1) = 0. In terms of J3 (r) PE ( 0o)= =- 4 (kw) sinp sin4o KE (b) (27) with 13

1 KE (b) = J3 (O) dC (28) and the angle dependence is explicit in the expression for PE (Q,,o). Computer programs were written to solve (20) and (25) by the moment method and, hence, compute KH (a) and KE (b) for all complex a and b. From an examination of the results it was found that KH (a) can be approximated as KH (a) = (- -(a + 0.15) (a + 0.29) (29) (- +Ln 2) (a+0.15) (a+0.29) +0.10a (a+0.20) for all a apart from those in the immediate vicinity of the portion -1.1 < a < 0.3 of the real axis is the complex a plane. In this region an empirical expression for KH (a) is KH (a)= — (30) -a + 2 2+0.1 and since, for other a, (30) differs from (29) by no more than 3 percent, it is sufficient to use (30) for all a. Similarly, for E-polarization the approximation is (b) 0.62 (b + 4.08) (b + 7.26) (b + 10.37) (b + 13.43) (b + 16.46) 31 E b + 1.15 (b + 4.27) (b + 7.37) (b + 10.45) (b + 13.49) (b + 16.50) valid for all b not in the immediate vicinity of the negative real axis. For positive real b, KE (b) < 1/2 and hence PE () (Po) < - (kw)2 sin( sin. (32) In their regions of validity, the estimated accuracy of (29) - (31) is about three percent. 14

To use these results to predict the scattering from a narrow gap, it was proposed that 1 be identified with the impedance looking into the gap, with n calculated using a simple transmission line (or other) model that takes into account the geometry and material filling of the gap. To show how this is done, consider a crack such as those illustrated in Figure 2. For H-polarization the cavity supports a variety of TE modes, but since the width w is small, the only mode which is not evanescent is the TEM mode, and this is the main contributor to the field in the gap. Under the assumption that this is the only mode that must be considered, the effective surface impedance r1 can be deduced from the input impedance Z1n of a parallel plate transmission line. The voltage across the gap is w/2 V = J Ex (x) dx = wEx -w/2 and since the current I is proportional to the tangential magnetic field, Ex 1 V Zln ^^IX^ (33) Hz w I w (33) For a parallel plate transmission line whose plate separation is w, the inductance and capacitance per unit length and width are L = o0grw and C = eoe1/w respectively, and the characteristic impedance is Zc = Z1 w with Z1 = Z J./ The L-shaped gap in Figure 2(b) can be viewed as two cascaded lines. The first line has length d1 and characteristic impedance Zc, whereas the second (of length 15

f Zen FZinC Zin Zin w w — w — w —I -— w 1d1 r i1d1 1\ ____ 2 I2 -] W24 W3 - W24 Fig. 2. Gap and cavity configurations. The cavity is filled with a homogeneous dielectric having relative permittivity e, and relative permeability l,.

w2) has characteristic impedance Zc = Z1d2 and is shorted. As a load its impedance is ZL = -i ZC tan klw2 (34) The junction of these lines can be modelled as a lumped parameter pi-network whose reactance and susceptance elements are [3] X = kZw d2 B ---d 1 -- 2 1- Z1 (d2+w z k w 2 B2= (d2+w (t - n 2) The input impedance of the first line cascaded with the pi-network and the second line is then ZL- iZ tan kd,(35) in- = zc, (35) Zc- iZ tan kid1 where ZL-iX (1 -iB2ZL) ZL = (1- BX) (1 - iBZL)- iBZL (36) Similarly, the T-shaped gap in Figure 2(c) can be treated as a transmission line loaded with two shorted lines in series. For the shorted lines of lengths w2 and w3, the load impedance is 17

ZL = -i Z tan kw2 + tan kw3} (37) The junction is modelled with a shunt susceptance and a series reactance in series with ZL: X3 = 2kZi wd2 k d2 3 d2 + 0.7822 where the constant was determined empirically, and the input impedance Z1i is then given by (35) with ZL- iX3 ZL = - i B3 (Z -iX3) 8) The rectangular gap is the special case d2 = 0 of either of the above structures, and for this Zin = -i Zc tan kd1. (39) Finally, for the V-shaped gap in Figure 2 (d), the inductance and capacitance per unit length of the line are functions of position, but when the coupled differential equations for the voltage and current are solved, we obtain Zn = c J, (kd,) (40) 18d

where Jo and J1 are Bessel functions. For E-polarization all of the modes are evanescent, but if we again assume that the first mode dominates in the gap, simple formulas for the surface impedance can be found. In a parallel plate waveguide of width w 1 aEz H-1 ikZ,r ay and for the lowest order mode the propagation constant is ikp where 1/2 P{( )2 (41 ) Since Ez/Hx is independent of position, a transmission line analogy can be made. The characteristic impedance of the line is -iZg/yp, which is also the impedance looking into the gap, and the results previously obtained for H-polarization are now applicable if k, is replaced by ikp and Z1 by -iZlr/p. Thus, for a rectangular gap ZR,w Zln - i- tanhkpd,, (42) and for a triangular gap ZJrw 1, (kpd 1) Zln=-i (43) P lo(kpd1) 19

where Z1n and n1 are related via (33) and 10o and 11 are modified Bessel functions [4]. Formulas for L- and T-shaped cracks can be deduced in a similar manner, but since the modes are evanescent, the shape of the lower cavity has little or no effect on the impedance. 4. Numerical Results The integral equation pairs (7), (11) and (17), (18) for H- and E-polarizations respectively were programmed for solution by the moment method, using pulse basis and point matching functions as described in Appendix A. In the case of (17), the derivative was applied to the kernel, and because of the order of the resulting singularity, the contributions from two cells on either side of the self cell were evaluated analytically, in addition to the contribution of the self cell itself. Comparison with the results of a finite element method [5] for H-polarization showed excellent agreement, and for purposes of comparison with the quasi-analytical solution, the moment method data will be regarded as exact. The computer program used to implement the expressions for the quasi-analytical solution is listed in Appendix B. Considering first the results for H-polarization, Figure 3 shows the backscattering from a rectangular air-filled gap as a function of aspect for three gap widths. The aspect variation decreases with w. It is less than 4 percent for w/A = 0.15, and since aspect independence is a feature of the quasi-analytic solution, we will henceforth confine attention to this case. It is then sufficient to take ( = (o = n/2 corresponding to normal incidence backscatter. 20

2tf -3 0 21T C.. a 2- 0~* n0 0 3 $ 1.6 IPHI 1 -0.5 -0 0 1 30 45 o0 76 90 00 in degrees Fig. 3; Modulus of th far field amplitude PH with respect to aspect 0o for a rectangular gap with * x/2 and d/. 0.2: w/X. 0.15, Q wX = 0.2, * w. 0.25. 21

2.6 IP"I l IPJ / ^ J \I 0- / i 0 0.2 0.4 0.6 0.8 1 Gap depth d/A Fig. 4. Modulus of the far field amplitude PH for a rectangular gap of varying di di = d with *. b * r2 and w/X. 0.15: exact, analytical. 22

180 120 80 arg PH o0 (deg.) *60 -120 -180 -0 0.2 0.4 0. 0. 1 Gap depth d/X Fig. 5. Argument of the far field amplitude PH for a rectangular gap of varying depth di d with 0. r/2 and w/X 0.15: | exact, -- analytical. 23

In Figures 4 and 5 the amplitude and phase of the far field amplitude PH (xi/2, i/2) are shown as a function of depth for a rectangular air-filled gap of width wA. = 0.15. We observe the cyclical behavior with zeros at d/k = 0, 0.5, 1.0,... resulting from the periodicity of the impedance looking into the gap. From (39) and (21) the corresponding a are real and vary from -oo to oo over each cycle. Over the entire range of d/k the agreement between the quasi-analytic and moment method results is excellent, but in spite of this the computed aperture impedances do not agree. This is evident from Figure 6 where I|E/Hz| is plotted as a function of x for w/X = 0.15 and d/X = 0.20. The U-shaped behavior is in accordance with the edge condition 1/2 at x = - w/2, and the data fit the curve C 1 - (with C = 860 ohms. The average value is therefore ncC/2 = 1350 ohms, compared with which (36) gives 11 I = 1 160 ohms. A similar discrepancy was found with all gap geometries. Nevertheless, the quasi-analytic solution provides an excellent approximation to the far field, and this is illustrated in Figures 7 through 10 showing IPHI for a material-filled rectangular gap and for air-filled L-, T- and V-shaped gaps. Turning now to E-polarization, Figures 11 and 12 show the amplitude and phase of PE (x/2, c/2) as functions of w/k for a rectangular air-filled gap having d/k = 0.1. The quasi-analytic and exact data diverge with increasing w/k, but the difference is less than 4 percent in amplitude and 5 degrees in phase for w/ <~ 0.20. 24

3000 -* 2500 2000 HIaI 100o (ohms) 1000 U 500 * -0.075 -O0.00 -0.025 0 0.025 0.060 O.075 x/X Fig. 6. Aperture impedance IE /HJ evaluated at - w/2 < x < w/2 and y=0 for a rectangular gap with -. o0. x/2, w/X. 0.15, and d/X 0.2: * exact, --- analyticl. 25

2.5 2 1.5,. -., i IPHI A v -.... I.... I..... I.... 0,0 0.1 0.2 0.3 0.4 0.6 Gap depth d/A Fig. 7. Modulus of the far field amplitude PH or a materalflled rectangular gap of varying depth dl d d with *0 * x/2, wX. 0.15, and,. 1: e, 2 + i1 exact, --- analytical.3 + i0.5 exact, - - - - - -analytical. 26

2 iPHI 01t 0.5 0 0.0 0.1 0.2 0.3 0.4 0.6 Gap depth d/X Fig. 8(a). Moduke of the far field amplitude PH for an air-filled L-shaped gap of varying depth di + d2 - d with *0 * x/2, w/X * 0.1 5, w2/X h 0.15, and d1/d2 3: exa - anaytical. 27

2- S m 2.5 2 1.5 IPHI: J ^ \ I 0'~~ ' ~.O~ O~.1 0.-2 0.3 O~.4 o.0.5 Gap depth d/A Fig. 8(b). Modulus of the far field amplitude, PH, for an air-filled L-shaped gap of varying depth d1 + d2 = d with = 0 = r/2, w = 0.15, and d1/d2 = 1: *&=0.05 * exact, - analytical 0.0.1 exact, analytical..0 0.0.2 0ec3 0. 0.6 w2/ =0.15 ~ exact, -—.analytical. 28

2 1.5 IPHI and di/d2 3: * exat, analytical. 29 29

2.5 2 1.5 0' ' 0', 0.0 0.1 0.2 0.3 0.4 0.5 Gap depth d/A Fig. 9(b). Modulus of the far field amplitude, PH, for an air-filled T-shaped gap of varyingdepthdi +d2=dwith= 0 =g/2,w/X=0.15,andd1/d2= 1: 30 0.0 0.1 0.2 0.3 0.4 0.6 Gap depth d/X Fig. 9(b). Modulus of the far field amplitude, PH, for an air-filled T-shaped gap of varying depth d1 + d2 = d with <) = ()o = t/2, w/. = 0.15, and dl/d2 = 1: wa/X = w3/X = 0.025 I exact, --— analytical w2/X = w3 = 0.075 ~ exact, —..analytical. 30

2.5 2 mum IPgHI 1 l 0.5 0P "i |i B i 1 i 0. 1 0.2 0.3 0.4 0.5 Gap depth d/X Fig. 10. Modulus of the far field amplitude PH for an air-filled V-shaped gap of varying depth dl. d with *. 0 x/2 and w/X 0.15: I exact, --- analyticl. 31

1.50 1.25 -1-.I* / 0.75 * / 0.50 o.so */ 0.25 0 _, ".. '.... ' 0.0 0.1 0.2 0. 0.4 Gap width w/X Fig. 1 1. Modulus of the far field amplitude PE for an air-filled rectangular gap c varying width with e = 0o = /2 and d/X 0.1: * exact, analytic 32

0 -20 -40 (do.J ) so / 0.00 0*00 0.3 WVidth W/r\ 0.4 Fig. 12. ArgUment ofthe^,, v'saynt wingh id.far e W ~nPlu e wwwjn *h A -o -. - A, for an air.llnlg.:et analytfcal 33

For a rectangular gap with w/k = 0.1 5, the quasi-analytic and exact results for PE (l IC) as a function of d/k are presented in Figure 13. The agreement is excellent, and as a consequence of the mode attenuation, the scattenng is independent of the depth for d/ 2 0.15. A similar comparison for a triangular gap is given in Figure 14. 34

0.25- 0.20 0.15-i Mall 0.10 0.06 -000-... I...... 00 0. 0. 0.1 0.2 0.3 4 Gap depth d/X Fig. 13. Modulus of the ar field amplitude PE for an airfilled rectangular gap of varying depth d, - d with * - u x/2and w/ 0.15: * exact, - analytcl. 35

0.25 0.20 0.15 IPEI O, lO 0.08 0.0 00 Gap depth d/A Fig. 14. Moduu of the far feld amtud PE for an air-fild V-shap gap o varying depth di d with o /2an w O.15: exact, -- anaytica. 36

5. Conclusions The quasi-analytic method described in [1] is based on the low frequency solution of the integral equations for a constant impedance insert in a perfectly conducting plane, and when used in conjunction with an estimate of the impedance looking into a gap, it provides a simple approximation to the far field scattering from the gap. To determine its accuracy, we have analyzed the problem of a plane wave incident on a gap backed by a cavity of arbitrary shape. The equivalence principle was used to develop coupled integral equations for the induced electric and magnetic currents, and the equations were then solved by the moment method. When the impedance looking into the cavity was determined using a transmission line model, it was found that for gap widths w/X 0.15 the quasi-analytic and moment method results for the scattered field were in excellent agreement for both polarizations and for all gap configurations that were tested. It therefore appears that the quasi-analytic method is an efficient and effective tool for predicting the scattering from the junction where two component parts of a target come together. 37

References [1] T.B.A. Senior and J.L. Volakis, "Scattering by Gaps and Cracks," IEEE Trans. Antennas Propagat., vol. AP-37, 1989. [2] R.F. Harrington, Time Harmonic Electromagnetic Fields, McGraw Hill Book Co., New York, 1961; pp. 106 et seq. [3] S. Ramo, J.R. Whinnery and T. Van Duzer, "Fields and Waves in Communication Electronics," John Wiley & Sons, New York, 1984, p. 569. [4] G.N. Watson, "A Treatise on the Theory of Bessel Functions," Cambridge Univ. Press, Cambridge, 1948, p. 77. [5] S.K. Jeng, "Aperture Admittance Matrix by Finite Element Method for Scattering by a Cavity-Backed Aperture," IEEE AP-S International Symposium, Syracuse, NY, June 1988. 38

Appendix A Moment Method Solution of the Coupled Integral Equations The integral equation pairs given by (7), (11) and (17), (18) are solved by the moment method. Using pulse basis functions in the moment method, the aperture A and the cavity walls S of Figure 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 (xi,y1), i = 1,...,N, which describe the location of each of the segments. The Hankel functions can then be expressed in terms of rotated coordinates (s,n) for the observation position and (si,ni) for each segment or source position Figure 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 (11) become J(sn) J(sn) e Jz(s) H(k (s-s) + n ds 1 -1 jAS i- As i + - $ XJS(snj'JsinYj 1 +(s-s +(-nj)2J ds9 (A.1) i-1 As. 39

J(s) =2e - kY JZ J H (k Is-sl) dsi (A.2) i1 AAsi where M are the number of segments across the aperture and (n-ni) sin 7 = - n (A.3) (S-Sj) + (n-ni) Applying point matching over the N segments of the aperture and cavity walls, N ",ik (1) [ I + - sin 1 (1)( kl R. ) ds. i. \ \.i= i ik j = i As, " r e Ii Ho k1 Rj i) ds = 0 j= 1,...,N (A.4) i-N+1 As, C ~iN+M ~ l kY Ii O(1) 1 -iks cos(o J 1 + |Ho (k Rj,) dsi =2e j = N+1,...,N+M (A.5) i = N+1 AS where Rj, = (j-s)2 +(nj-n)2 (A.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. Ii 40

in (A.4) and (A.5) are the electric currents, for i = 1,...,N, and the magnetic currents, i = N+1,...,N+M, to be determined. In matrix form, (A.4) and (A.5) become [ zj] [ ][ V] (A.7) The impedance matrix is given as zL z z [ i =.............. (A.8) where the sets of elements are as follows: ik1 (1) isin Y H1 (k1 R,) ds ji Zel a si (A.9). -1 j=i for i =1,...,N and j =1,...,N; ' k r e Ho )(k1 Rj ) ds, j i-N =z AS } (A. 10) 12 ~Er{ 9[2(sj-si) (Rj)- (2-A'(k)) (s-sj)] j = i-N for i = N+1,...,N+M and j = 1,...,N; 41

{ 0 j i+N Z j=i+N (A.11) Ze2 1 j = i+N for i = 1,...,N and j = N+1,...,N+M; kY |H( )(k Rji) ds j i Z 2 si (A. 12) m2 2 I{ 2(sj-s) RL (R,,)- (2-A'(k)) (s,-sj)] } j= i for i = N+1,...,N+M and j = N+1,...,N+M. In (A.10), the expression for A' is kl: A'(k)=2(em 2 Y+-i ) In (A.12), A' is a function of k. For the self-cells in (A.10) and (A.12), s5 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. The analytical expressions for the impedance elements of the adjacent cells are k2 (n.-n.) s.-s. Ze1 -2 2(n.-n.i)s. + - atan- I i i=j+1 (A.13) el 2 2 -I I) S ' I-nj-ni). 42

Z1 K2 r 2(Si-Sj) t (R,i) - (2-A'(k,)) s + 2|1n-nl atan |i- | ml 2 e'rn' i-N =j+l (A. 14) where si is evaluated at the endpoints of the ith segment. The adjacent cell expression for Zm2 is given by (A.14) with A'(k1) replaced with A'(k) for i = j+1. The source matrix is given by 0 j=1,...,N V. = -iks cos<o (A. 15) 1 2e 0 j=N+1,...,N+M. The currents Ii are determined by solving (A.7), given that [Zj i] is nonsingular. The aperture impedance is defined in terms of the total fields as Ex(xi,0) = Hz(xj,O) where the total electric field is equal to the magnetic current in the aperture, Ij for j = N+1,...,N+M, and from (4), the total magnetic field is now expressed as Hz(xj.,) = 2e ' - TY i I(1)( H 1o k I dx. (A. 16) for j = N+1,...,N+M. From (3), the far field amplitude at the angle is now 43

kV r -ikx co PH-2 Ij. Je dx,. (A. 17) i =N+1 Ax A.2 E-Polarization The integral equation pair given by (17) and (18) was solved in the same manner as described for the H-polarization case. The elements of the impedance matrix defined in (A.8) for the E-pol case are as follows: kZ r |Ho (k1 Rji) dsi j i Ze1={ ASi (A.18) t 1 j=i for i = 1,...,N and j= 1,...,N; I |S R H1 (1) R.j) dsi j i-N Zm1 = A Si Jl (A. 19) t -1 j=i-N for i = N+1,...,N+M and j = 1,...,N; f 0 j i+N Z9 l - 1 j i+N (A.20) for i = 1,...,N and j = N+1,...,N+M; 44

H (k Rj..) ds. J W As, i' Z m2 (A.21) -- - 2[2(S-Sj) 'n(Rj i) - (2-A'(k)) (s-sj) - H1 k Is-s l) j i ~~2 L i i i J k 1 ~-= for i = N+1,...,N+M and j = N+1,...,N+M. As for the H-polarization case, in the self-cell expressions, si is evaluated at the endpoint of the ith segment. Because of the sensitivity of the impedance element Zm2 to the 1/Rj,i term for segments near the self-cell, the adjacent cells needed to be evaluated analytically, as follows: k. kY i 1 8 m2 2 K 3k2 s. 3k As. S.-S. + 2(si-s ) (Rj,) — A'(k))s. + 21n -nl ata t inji (A.22) for i = j~1, where si is evaluated over the ith segment. The source matrix is given by 0 j=1,...,N Vj= -iks o(A.23) 2Ysin~oe = j=N+1,...,N+M Given that [Zj,] is nonsingular, the currents Ii can be calculated. 45

The aperture impedance for the E-polarization case is defined as Ez(xj.O) iH(x1iO) where the total electric field is equal to the magnetic current in the aperture, Ij for j = N+1,...,N+M, and from (13), the total magnetic field in the aperture is now expressed as -ikxj cos|o kY 1 (1) Hx(xj,O) = - 2Y sino, e +. Ij (klx -x l) dx, (A.24) ijN+1 x J ' for j = N+1,...,N+M. From (12), the far field amplitude at the specified angle < is now k Mf -ikxi cost PE( =) ksin I |f e dx. (A.25) i-N+1 Axi A.3 Program Listing The expressions for the impedance and source matrices, the aperture impedance, and the far field amplitude were programmed for solution, as shown in the program listing of GAPSCAT.FTN below. The subroutines used in the program are contained in the file GAPSUB.FTN listed below also. In running the program, the user is prompted for the polarization of the incident field, angle of incidence, angle of far field observation, and the relative permittivity e, of the gap cavity. For the relative permeability, it is assumed that r = 1, although this need not be the case. A menu is provided for the choice 46

of shapes as shown in Figure 2. The dimensions are requested according to those defined in Figure 2. An arbitrary shaped gap may also be evaluated by specifying the coordinates of its corner points. The user is also prompted for the maximum segment size Asi to be used for the pulse basis functions. A segment size of Asi/k = 0.01 was used for the results of Figures 3 to 14. The impedance matrix [ Zji ] is solved for the H- or E-polarization case using the expressions (A.9) to (A.14) or (A.18) to (A.22), respectively. The numerical integration is done for the appropriate segments using Simpson's three-point composite integration over each segment. With the source matrix [ Vj ] calculated from (A. 15) or (A.23), the electric and magnetic currents contained in [ Ij ] can then be determined. As listed, the program calculates the far field amplitude as a function of the gap depth using (A. 17) for H-polarization or (A.23) for E-polarization, where the number of iterations is specified. For one iteration, the program also outputs the aperture impedance calculated from the total fields in the aperture. 47

C C GAPSCAT.FTN C C C This FORTRAN program computes the far field scattering due C to a narrow gap of arbitrary shape in an infinite ground 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, and number of C iterations with respect to gap depth. C C OUTPUT FILES C GAPDAT Contains input data. C IMPDAT For one iteration, field or impedance C in the aperture of 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 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,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,ctemp,er,ur,Ao,A,kl complex Z(pn,pn),LHO,LH1,LH2,aLHO,Ii(pn),Vj (pn) complex eta (pn),Hi(pn),Hs(pn),E(pn),Lsca,Psca integer ipvt (pn), iretrn real rc,krho complex wk(pn),HO, H1,HOo,Hlo,ckrho logical Epol,Lossy,side,neg(50) common /prompts/ EorH,phio,phi,er,igap,wStp,dStp,w,d,noS, & q,maxC,noIter 1 format(il) 2 format(iS) 3 format(gl6.8) 4 format(a4) 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(O.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 ur-cmplx (1., 0.0) gam-0.5772157 Ao-2* (log (k/2)+gam-ci*pi/2) iprg-1 c...Setting default values EorH-1 phio-90.0 phi-90.0 er-cmplx(1.,0.0) w-0.15 d-0.2 noS-3 maxC-0.01 noIter-30 adj-0.00001 Epol-. false. Lossy-. false. side-.true. 48

c...Prompting user for input data call gaprom(iprg) if(EorH.eq. 1) Epol-.true. phio-phio*pi/180.0 phi-phi*pi/180.0 drat-dStp (1) /d kl=k*csqrt (er) A-2*(clog(kl/2)+gam-ci*pi/2) if(aImag(er).ne. 0.0) Lossy-.true. dmin-0.025 if(Epol) dmin-0.025 dmax-d if(noIter.ne. 1)then dstep- (dmax-dmin) / (noIter-1) d-dmin endif DO 700 iter-l,noIter c...Determining coordinates of corner points given gap type 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)then c L-SHAPED noS-7 dStp (1) -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) — w/2-wStp(3) q(7,2)-q(4,2) q(5, 1)-q(4,1) q(5,2)-q(4, 2)-dStp(2) q(6,1)-q(7, 1) q(6, 2)-q(5, 2) else if(igap.eq. 3)then c V-SHAPED noS-2 q(3, 1) -0.0 q(3,2) —d ad -0.75*maxC else if(igap.eq. 4)then c T-SHAPED noS-5 dStp (1) -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)-dStp (2) q(6,1) —w/2 q(6, 2)-q(5, 2) endif c...Corner points of gap at y-0 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 c******************************************************************* c************* Current segment locations (xi,yi) ************** c**+*************************************************************** N-0 do 175 1-1,noS+l c...Size (length) of lth side of gap szSd(1) -sqrt ((q(1+1, 1)-q(l, 1)) **2 & + (q(l+l, 2) -q(1, 2) ) **2) c..Angle of rotation for each side with respect to x axis if(q(l+1,1).lt. q(l,1))then psi (1) -asin( (q(1,2)-q (1+1,2) )/szSd (1)) neg(l) -.true. else psi(1)-asin ( (q(1+1,2)-q(1,2))/szSd(1)) neg(1)-.false. endif szN (1) -int (szSd (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) 49

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-1,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 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 LH2 Integral of HOo for H-pol, of dHlo/dy C 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 -1,...,N, C in the clockwise direction, starting at (x,y)-(-w/2,0). C istart-1 istop-szN(l) c...Source point is i of the Ith side, observation point is j do 230 1-1,noS+1 do 220 i-istart,istop do 210 j-l,N c Coordinate rotation for observation point sJ-m(, l)*cos(psi(l))+m(j,2)*sin(psi(l)) nj-m(j,l)*sin(psi(l))-m(j,2)*cos (psi(1)) if (neg(1))then sj —sj nj —nj endif c... Integration over ith segment LHO-czero LHl-czero LH2-czero aLHO-czero c Magnitude between midpoints Rm-lr-r'l Rm-sqrt((m(j,l)-m(i,l))**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(l))+p (ip, 2)*sin(psi (l)) ni-p(ip, 1)*sin(psi(1))-p(ip,2)*cos(psi(1)) if (neg (1)) 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-1.0 else tanf-atan((si-sj)/abs(nj-ni)) absf-abs (nj-ni) endif LH1 —kl**2/2* (nj-ni) *si & +ci*2./pi*(nj-ni)/absf*tanf-LH1 LHO0-ci/pi* (2* (si-s) *log (R) & - (2-A) *si+2.0*abs (nj-ni) *tanf) -LHO LH2-ci/pi* (2* (si-sj) *log (R) & - (2-Ao) *si+2.0*abs (nj-ni) *tanf) -LH2 if(j.eq. i) GOTO 202 50

200 continue 202 if(j.eq. i)then LHO-2*(LHO+ci/pi*(2-A)*sj) LH2-2*(LH2+ci/pi*(2-Ao)*sj) if(Epol)then krho-k*abs (sj-si) call Hankzl(krho,l,HO,H1) LH2-LH2-2./k*Hl 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, 1) *sin(psi(1))-p (ip, 2)*cos(psi (1)) if(neg ())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 krho-k*R call Hankzl(krho,2,HOo,Hlo) LHO-HO+LHO if(Epol)then LHl-kl*m(j, 2) /R*H1+LH1 if(j.eq. (i-1).or. j.eq. (i+l))then if(abs(nj-ni).eq. 0.0)then tanf-pi/2 absf-1.0 else tanf-atan((si-sj)/abs(nj-ni)) absf-abs (nj-ni) endif aLHO-ci/pi* (2* (si-sj) *log(R) c& ~ -(2.0-Ao)*si+2.0*abs(nj-ni)*tanf)-aLHO else LH2-Hlo/k/R+LH2 endif else LHl-kl* (nj-ni) /R*H+LH1 LH2-HOo+LH2 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 (1)) 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,Hl) else krho-Real (kl) *R call Hankzl(krho,2,HO,H1) endif krho-k*R call Hankzl (krho,2,HOo,Hlo) c GREEN'S FUNCTION INTEGRALS LHO-stepS/3* (4*HO+LHO) if (Epol) then LHl-stepS/3*(4*kl*m(J,2)/R*H1+LH1) if(j.eq. (i-1).or. j eq. (i+l))then LH2 —ci*8./3/pi/k**2/ (DelS)-aLHO else LH2-stepS/3* (4*Hlo/k/R+LH2) endif else LHl-stepS/3* (4*kl* (nJ-ni)/R*H1+LH1) LH2-stepS/3* (4*HOo+LH2) 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 51

if(i.le. gN.and. j.ne. i)then Z(,N+i) —ci/2*LH1 else if(i.le. gN.and. j.eq. i)then Z(j,N+i) —1. endif if(j.le. gN)then if(j.ne. i)then Z(N+, i)-czero else Z(N+J,i) —1. endif if(i.le. gN)then Z(N+j,N+i) —k*Yo/2*LH2 endif endif else C H-POL IMPEDANCE MATRIX if(j.ne. i)then Z (j, i) -ci/2*LH1 else Z(j,i) —l. endif if(i.le. gN)then Z (j,N+i) -k*Yo*er/2*LH0 endif if(j.le. gN)then if(j.ne. i)then Z(N+j,i)-czero else Z(N+j,i)-l. endif if(i.le. gN)then Z(N+, N+i)-k*Yo/2*LH2 endif endif endif 210 continue c...Incident Field (Source) matrix elements xj-m(i,l) yj-m(i,2) Vj(i)-czero if(i.le. gN)then if (Epol) then c E-POL INCIDENT FIELD Hx Vj(N+i)-2*Yo*sin(phio)*cexp(-ci*k*xj*cos(phio)) else c H-POL INCIDENT FIELD Hz Vj (N+i) -2*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 if(noIter.eq. l)then C****************************************************************** c******************** Aperture Impedance ******************** c**************************************************************** 1-1 do 500 j-1,gN Hs(j)-czero do 480 i-l,gN s-m(j, 1) *cos (psi (1))+m(,2)*sin(psi(1)) nj-m(j, 1)*sin(psi(l))-m(j,2)*cos(psi(1)) if (neg (1)) then n -n endif c...Integration over ith segment LH2-czero aLHO-czero if(j.eq. i)then c SMALL ARGUMENT APPROXIMATION integral for self-cell c and adjacent cells do 470 ip-i+l,i,-1 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 (1)) -p (ip, 2) *cos (psi (1)) if (neg (1)) then si —si 52

ni —ni endif R-sqrt ((sj-si) **2+(nj-ni) **2) if(j.eq. i.or. abs(nj-ni).eq. 0.0)then tanf-0.0 absf-1.0 else tanf-atan((si-sj)/abs(nj-ni)) absf-abs (nj-ni) endif LH2-ci/pi*(2*(si-s) *log(R) & - (2-Ao)*si+2.0*abs(nj-ni) *tanf)-LH2 if(j.eq. i) GOTO 472 470 continue 472 if(j.eq. i)then LH2-2*(LH2+ci/pi*(2-Ao)*sj) if(Epol)then krho-k*abs(sj-si) call Hankzl(krho,l,HOo,Hlo) LH2-LH2-2./k*Hlo endif endif else c SIMPSON'S THREE POINT COMPOSITE INTEGRATION do 474 ip-i+l,i,-1 c Coordinate rotation for source segment endpoints si-p (ip, 1)*cos(psi(1))+p (ip, 2)*sin (psi(1)) ni-p (ip, 1) *sin (psi (1))-p (ip, 2) *cos (psi (1)) if (neg (1)) then si —si ni —ni endif stepS-si R-sqrt ((sj-si) **2+ (nj-ni) **2) krho-k*R call Hankzl(krho,2,HOo,Hlo) if(Epol) then if(j.eq. (i-1).or. j.eq. (i+l))then if(abs(nj-ni).eq. 0.0)then tanf-pi/2 absf-1.0 else tanf-atan((si-sj)/abs(nj-ni)) absf-abs (n-ni) endif aLHO-ci/pi* (2* (si-sj) *log (R) & -(2.0-Ao) *si+2.0*abs(nj-ni)*tanf) -aLHO else LH2-Hlo/k/R+LH2 endif else LH2-HOo+LH2 endif 474 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 (1))then si —si ni —ni endif s (sS-abs(stepS-si) R-sqrt ((sj-si)**2+ (n-ni) **2) krho-k*R call Hankzl(krho,2,HOo,Hlo) if(Epol)then if(j.eq. (i-1).or. j.eq. (i+l))then LH2 —ci*8./3/pi/k**2/ (2*stepS)-aLHO else LH2-stepS/3* (4*Hlo/k/R+LH2) endif else LH2-stepS/3*(4*HOo+LH2) endif endif c SCATTERED MAGNETIC FIELD IN THE APERTURE Hs( j)k*Yo/2*Ii(N+i)*LH2+Hs () 480 continue c ELECTRIC FIELD IN THE APERTURE E (j)-I (N+J) xj-m(j,l) y -m( j,2) c INCIDENT MAGNETIC FIELDS if(Epol)then Hi (j)-2*Yo*sin(phio) *cexp(-ci*k*x*cos (phio)) else Hi (j)-2*cos (k*yj*sin (phio)) *cexp (-ci*k*xj*cos (phio)) endif c APERTURE IMPEDANCE eta () -E(j) / (Hi (j) +Hs (j)) write(2,*) m(j,l),cabs(eta(j)) 500 continue endif 53

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+l,1)*cos(phi))) Psca-Ii (N+i) *Lsca+Psca 600 continue if (Epol)then Psca —k*sin(phi)/2*Psca else Psca —k*Yo/2*Psca endif write (3,*) d, cabs(Psca) write(4,*) d,180/pi*(atan2(aImag(Psca),Real(Psca))) print *,' Exact: IPscaI - ',cabs(Psca),' arg Psca & 180/pi*(atan2(aImag(Psca),Real(Psca))) d-d+dstep 700 continue 800 call exit END 54

C C GAPSUB.FTN C C This file contains the subroutines and functions used by C GAPSCAT.FTN and ANAGAP.FTN. C C******************************************************************* SUBROUTINE GAPROM(IPRG) C********************************************************************* C C Called to prompt user for the input parameters. C integer EorH,N,noS,gN real phi,phio,w,d,maxC,q(50,2) real dStp (50),wStp(50) complex er,ctemp common /prompts/ EorH,phio,phi,er,igap,wStp,dStp,w,d,noS, & q,maxC,noIter 1 format(il) 2 format(i5) 3 format(g16.8) 6 format(2gl6.8) open (1, file-'gapdat') print *,' ' 20 write(*,25) EorH 25 format(/'Polarization of incident field:'/ & ' 1) E- or 2) H-pol [',il,']? ') read(*, 1, err-20) itemp if(itemp.eq. 2) EorH-2 30 write(*,35) phio 35 format(/,'Angle of incidence - ',gl2.6,'degrees ') read(*,3,err-20) rtemp if(rtemp.ne. 0.0) phio-rtemp 36 write(*,37) phi 37 format(/,'Angle of observation - ',gl2.6,'degrees ') read (*, 3, err-30) rtemp if(rtemp.ne. 0.0) phi-rtemp 38 write(*,39) er 39 format(/,'Relative permittivity - (',g12.6,',',g12.6,')') read(*,6,err-36) ctemp if(ctemp.ne. 0.0) er-ctemp print *,er 40 write(*,45) 45 format (/'Type of gap:'/ & ' 1) rectangular, 2) T-shape, 3) triangular,'/ & ' 4) L-shape, or 5) arbitrary? ') read(*, 1, err-20) igap if(igap.eq. 2.or. igap.eq. 4)then if(igap.eq. 2)then now-3 else now-2 endif do 47 i-l,now write(*,*) 'Enter width w',i read(*,*) wStp(i) if(i.ne. 3)then write(*,*) 'Enter depth d',i read(*,*) dStp(i) endif 47 continue w-wStp (l) d-dStp (1) +dStp (2) else if(igap.eq. 1.or. igap.eq. 3.or. & (igap.eq. 5.and. iprg.eq. 1))then 90 write(*,95) w 95 format (/,'Gap width - ',g12.6,'wavelengths ') read(*, 3,err-30) rtemp if(rtemp.gt. 0.0) w-rtemp 100 write(*,105) d 105 format (/'Gap depth - ',g12.6,'wavelengths ') read (*, 3, err-90) rtemp if(rtemp.gt. 0.0) d-rtemp if(igap.eq. 5)then 110 write(*,115) noS 115 format(/,'Number of sides (excluding aperture) - ',12) read (*,*, err30) rtemp if(rtemp.gt. 0.0) noS-int(rtemp) print *,'Enter the following coordinates, beginning ' print *,' with (-w/2,0) and going cw: ' do 127 i-l,noS+1 120 write(*,125) i 125 format(/'Corner ',il,' (x,y) ') read(*,3,err-110) q(i,l) read(*,3,err-110) q(i,2) 127 continue noIter-1 endif else GOTO 40 55

endif if(iprg.eq. 1)then 130 write(*,135) maxC 135 format(/'Max segment size - ',g12.6,' wavelengths ') read(*,3,err-110) rtemp if(rtemp.gt. 0.0) maxC-rtemp endif if(igap.ne. 5)then 140 write(*,145) noIter 145 format(/'Number of iterations - ',i3,' ') read(*, *,err-130) rtemp if(rtemp.gt. 0.0) noIter-int(rtemp) endif c...Writing input data to file GAPDAT write(1,2) EorH write(1,3) phio write(1,3) phi write(1,6) er write(1,2) igap if(igap.eq. 2.or. igap.eq. 4)then do 150 i-1,2 write(1,3) wStp(i) write(1,3) dStp(i) 150 continue else write(1,3) w write(1,3) d if(igap.eq. 5)then write(1,2) noS do 160 i-l,noS+1 write(1,6) q(i, 1),q(i,2) 160 continue endif endif if(iprg.eq. 1)then write(1,3) maxC endif write(1,2) noIter close () return end C************************www*www*w*w*w *~~wwwrwww**~*~w***w***w**** COMPLEX FUNCTION CTAN(CARG) C****************************************************************** C Calculates the tangent given a complex argument C complex ci,carg ci-cmplx (0., 1) ctan —ci*(cexp(ci*carg)-cexp(-ci*carg)) & / (cexp (ci*carg) +cexp (-ci*carg)) return end C***************************************************************** C C SUBROUTINE HANKZ1(R,N,HZERO,HONE) C C C*****************************************************************C C C Called to compute Hankel functions of the first kind C for orders one and zero. The argument is variable R C and must be positive. C C.....HANKEL FUNCTIONS ARE OF FIRST KIND —J+IY C..... N-0 RETURNS HZERO (H-zero) C..... N-1 RETURNS HONE (H-one) C..... N-2 RETURNS HZERO AND HONE C.....SUBROUTINE REQUIRES R>0 C.....SUBROUTINE ADAM MUST BE SUPPLIED BY USER C DIMENSION A(7),B(7),C(7),D(7),E(7),F(7),G(7),H(7) COMPLEX HZERO,HONE DATA A,B,C,D,E,F,G,H/1.0,-2.2499997,1.2656208,-0.3163866, &0.0444479,-0.0039444,0.00021,0.36746691,0.60559366,-0.74350384, &0.25300117,-0.04261214,0.00427916,-0.00024846,0.5,-0.56249985, &0.21093573,-0.03954289,0.00443319,-0.00031761,0.00001109, &-0.6366198,0.2212091,2.1682709,-1.3164827,0.3123951,-0.0400976, &0.0027873,0.79788456,-0.00000077,-0.0055274,-0.00009512, &0.00137237,-0.00072805,0.00014476,-0.78539816,-0.04166397, &-0.00003954,0.00262573,-0.00054125,-0.00029333,0.00013558, &0.79788456,0.00000156,0.01659667,0.00017105,-0.00249511, &0.00113653,-0.00020033,-2.35619449,0.12499612,0.0000565, &-0.00637879,0.00074348,0.00079824,-0.00029166/ IF (R.LE.0.0) GO TO 50 IF (N.LT.0.OR.N.GT.2) GO TO 50 IF (R.GT.3.0) GO TO 20 X-R*R/9.0 IF (N.EQ.1) GO TO 10 CALL ADAM(A, X,BJ) CALL ADAM(B,X,Y) BY-0.6366198*ALOG(0.5*R) *BJ+Y HZERO-CMPLX (BJ, BY) IF (N.EQ.0) RETURN 10 CALL ADAM(C,X,Y) 56

BJ-R*Y CALL ADAM(D,X,Y) BY-0.6366198*ALOG(0.5*R) *BJ+Y/R HONE-CMPLX (BJ, BY) RETURN 20 X-3.0/R IF (N.EQ.1) GO TO 30 CALL ADAM(E,X,Y) FOOL-Y/SQRT (R) CALL ADAM(F,X,Y) T R+Y BJ-FOOL*COS (T) BY-FOOL*SIN(T) HZERO-CMPLX (BJ, BY) IF (N.EQ.0) RETURN 30 CALL ADAM(G,X,Y) FOOL-Y/SQRT (R) CALL ADAM(H, X,Y) T-R+Y BJ-FOOL*COS (T) BY-FOOL*SIN(T) HONE-CMPLX (BJ, BY) RETURN 50 WRITE(6,90) N,R 90 FORMAT(32HOSICK DATA IN HANKZ1 *QUIT* N-,I2,2X,2HR-,E11.3) CALL SYSTEM END C*************************** * *********************** *********** C C SUBROUTINE ADAM (C, X, Y) C C C*******************************************************************C C C Called by subroutine HANKZ1 to compute the value of a 7th C order polynomial whose argument is X and coefficients are C contained in vector C. C DIMENSION C(7) C Y-X*C (7) C DO 10 I-1,5 Y-X* (C (7-I) +Y) 10 CONTINUE C Y-Y+C(1) RETURN END C******* ***************************** ******************************* C C SUBROUTINE CHANK(Z,N,HO,H1) C C C********************************************************************* C CALCULATES HANKEL FUNCTIONS ZEROTH AND FIRST ORDER C OF THE FIRST KIND OF COMPLEX ARGUMENT HO-JO+I*YO, H1-J1+I*Y1. C THE ACCURACY IS VERY GOOD UP TO IZI<10. ABOVE THAT THE C ACCURACY IS TO THE ORDER OF LARGE ARGUMENT EXPANSION. C C N-0 RETURNS HO C N-1 RETURNS HI C N-2 RETURNS HO AND HI C N-3 RETURNS JO AND JI C COMPLEX JO,YO,J1,Y1,Z,TERM COMPLEX I,HO,H1 I-(0.0,1.0) PI-4.*ATAN(1.) IF(N.EQ.1)GOTO10 IF(CABS(Z).GT. 12.)THEN JO-CSQRT (2/(PI*Z)) *CCOS (Z-PI/4.) YO-CSQRT (2/(PI*Z) ) *CSIN (Z-PI/4.) ELSE AONE-1. TERM- (1.0,0.0) JO-(1.0,0.0) YO-(0.0,O.O) SUM-0.0 M-0 100 M-M+1 SUM-SUM+1. /FLOAT (M) AONE-AONE (-1. ) TERM-TERM*( (Z/2)/FLOAT(M) )**2 JO-JO+AONE*TERM YO-YO+AONE*TERM*SUM IF(M.LE.10.OR. CABS(FLOAT(M)/(Z/2)).LT. 5)GOTO 100 YO-2.*(JO* (CLOG(Z/2)+0.57721)-YO)/PI ENDIF IF(N.EQ.3)THEN HO-JO ELSE HO-JO+I*YO ENDIF IF (N.EQ.O) RETURN C COMPUTATION OF J1 AND Y1 10 IF (CABS(Z).GT. 15.)THEN J1-CSQRT (2/ (PI*Z) ) *CCOS (Z-3.*PI/4.) 57

Y1-CSQRT(2/(PI*Z) ) *CSIN(Z-3.*PI/4.) ELSE AONE-1. TERM-(1.0,0.0) J- (1.0,0.0) Y1- (1.0,0.0) SUM-0.0 M-0 200 M-M+1 SUM-SUM+1. /FLOAT (M) AONE-AONE* (-1.) TERM-TERM* (Z/2)**2/(FLOAT(M) *FLOAT (M+1)) J1-J1+AONE*TERM Y1-Yl+AONE*TERM* (1./FLOAT (M+1) +2. *SUM) IF(M.LE.10.OR. CABS(FLOAT(M)/(Z/2)).LT. 5)GOTO 200 J1-J1*Z/2. Y1-(2.*J1*(CLOG(Z/2)+0.57721)-2./Z-Y1*Z/2.)/PI ENDIF IF (N.EQ.3)THEN H1IJ1 ELSE Hl=JI+I*Y1 ENDIF RETURN END C******************************************************************* C SUBROUTINE MODBESS(X, I0,I1) C C**** *************************************************************** C CALCULATES THE MODIFIED BESSEL FUNCTION OF THE C ZEROTH AND FIRST ORDER. ARGUMENT X IS POSITIVE AND REAL. C SEE PAGE 378 ABRAMOVITZ C REAL I0,I1 T-X/3.75 IF(T.LT. -1)THEN PRINT *,'ERROR' STOP ELSE ENDIF IF(T.GE. -1.OR. T.LE.1)THEN I0-1+3.5156229*T**2+3.0899424*T**4+1.2067492*T**6+ 6 0.2659732*T**8+0.0360768*T**10+0.0045813*T**12 I11-0.5+0.8789059*T**2+0.51498869*T**4+0.15084934*T**6+ & 0.02658733*T**8+0.00301532*T**10+0.00032411*T**12 I1-X*I1 ELSE I0-0.39894228+0.01328592/T+0.00225319/T**2-0.00157565/T**3 & +0.00916281/T**4-0.02057706/T**5+0.02635537/T**6 & -0.01647633/T**7+0.00392377/T**8 I0-IO*EXP (X) *SQRT (1/X) I1-0.39894228-0.03988024/T-0.00362018/T**2+0.00163801/T**3 & -0.01031555/T**4+0.02282967/T**5-0.02895312/T**6 & +0.01787654/T**7-0.00420059/T**8 ENDIF RETURN END C******************************************************************c C C C The following subroutines are standard LINPACK routines C C to perform L-U decomposition and back substitution on a C C single precision complex matrix. See CC-Memo 407 sec 2.1 C C for documentation on these routines. C C C C** **************************************************************** C C SUBROUTINE CGECO (A, LDA, N, IPVT, RCOND, Z) C C C********* ******************************************************** *C C C C NAASA 2.1.042 CGECO FTN-A 05-02-78 THE UNIV OF MICH COMP CTR C INTEGER LDA,N,IPVT(1) COMPLEX A(LDA,1), Z (1) REAL RCOND C C CGECO FACTORS A COMPLEX MATRIX BY GAUSSIAN ELIMINATION C AND ESTIMATES THE CONDITION OF THE MATRIX. C C IF RCOND IS NOT NEEDED, CGEFA IS SLIGHTLY FASTER. C TO SOLVE A*X - B, FOLLOW CGECO BY CGESL. C TO COMPUTE INVERSE(A)*C, FOLLOW CGECO BY CGESL. C TO COMPUTE DETERMINANT(A), FOLLOW CGECO BY CGEDI. C TO COMPUTE INVERSE(A), FOLLOW CGECO BY CGEDI. C C ON ENTRY C C A COMPLEX (LDA, N) C THE MATRIX TO BE FACTORED. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY A. C C N INTEGER C THE ORDER OF THE MATRIX A. 58

C C ON RETURN C C A AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS C WHICH WERE USED TO OBTAIN IT. C THE FACTORIZATION CAN BE WRITTEN A - L*U WHERE C L IS A PRODUCT OF PERMUTATION AND UNIT LOWER C TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR. C C IPVT INTEGER(N) C AN INTEGER VECTOR OF PIVOT INDICES. C C RCOND REAL C AN ESTIMATE OF THE RECIPROCAL CONDITION OF A C FOR THE SYSTEM A*X - B, RELATIVE PERTURBATIONS C IN A AND B OF SIZE EPSILON MAY CAUSE C RELATIVE PERTURBATIONS IN X OF SIZE EPSILON/RCOND C IF RCOND IS SO SMALL THAT THE LOGICAL EXPRESSION C 1.0 + RCOND.EQ. 1.0 C IS TRUE, THEN A MAY BE SINGULAR TO WORKING C PRECISION. IN PARTICULAR, RCOND IS ZERO IF C EXACT SINGULARITY IS DETECTED OR THE ESTIMATE C UNDERFLOWS. C C Z COMPLEX(N) C A WORK VECTOR WHOSE CONTENTS ARE USUALLY UNIMPORTANT. C IF A IS CLOSE TO A SINGULAR MATRIX, THEN Z IS C AN APPROXIMATE NULL VECTOR IN THE SENSE THAT C NORM(A*Z) - RCOND*NORM(A)*NORM(Z) C C LINPACK. THIS VERSION DATED 07/14/77 C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LABS. C C SUBROUTINES AND FUNCTIONS C C LINPACK CGEFA C BLAS CAXPY,CDOTC,CSSCAL,SCASUM C FORTRAN ABS,AIMAG,AMAX1,CMPLX,CONJG,REAL C C INTERNAL VARIABLES C COMPLEX CDOTC,EK, T,WK,WKM REAL ANORM,S,SCASUM,SM,YNORM INTEGER INFO, J,K,KB,KP1, L C COMPLEX ZDUM,ZDUM1,ZDUM2,CSIGN1 REAL CABS1 CABS1(ZDUM) - ABS(REAL(ZDUM)) + ABS(AIMAG(ZDUM)) CSIGN1(ZDUM1,ZDUM2) - CABS1(ZDUM1)* (ZDUM2/CABS1(ZDUM2)) C CCC Compute 1-NORM of A C ANORM - O.OEO DO 10 J - 1, N ANORM - AMAX1(ANORM,SCASUM(N,A(1,J),1)) 10 CONTINUE C CCC Factor C CALL CGEFA(A,LDA,N,IPVT,INFO) C C RCOND - 1/(NORM(A)*(ESTIMATE OF NORM(INVERSE(A)))) C ESTIMATE - NORM(Z)/NORM(Y) WHERE A*Z - Y AND CTRANS(A)*Y - E C CTRANS(A) IS THE CONJUGATE TRANSPOSE OF A. C THE COMPONENTS OF E ARE CHOSEN TO CAUSE MAXIMUM LOCAL C GROWTH IN THE ELEMENTS OF W WHERE CTRANS(U)*W - E C THE VECTORS ARE FREQUENTLY RESCALED TO AVOID OVERFLOW. C C SOLVE CTRANS(U)*W - E C EK - CMPLX(1.OE0,0.OEO) DO 20 J - 1, N 2(J) - CMPLX(0.OE0,O.OEO) 20 CONTINUE DO 100 K - 1, N IF (CABS1(Z(K)).NE. 0.OEO) EK - CSIGN1(EK,-Z(K)) IF (CABS1(EK-Z(K)).LE. CABS1(A(K,K))) GO TO 30 S - CABS1(A(K,K))/CABS1(EK-Z(K)) CALL CSSCAL(N,S,Z,1) EK - CMPLX(S,O.OEO)*EK 30 CONTINUE WK - EK - Z(K) WKM - -EK - Z(K) S - CABS1(WK) SM - CABS1(WKM) IF (CABS1(A(K,K)).EQ. O.OEO) GO TO 40 WK - WK/CONJG(A(K,K)) WKM - WKM/CONJG (A(K,K)) GO TO 50 40 CONTINUE WK - CMPLX(1.OE0,0.OEO) WKM - CMPLX(1.OEO,O.OEO) 50 CONTINUE KP1 - K + 1 IF (KP1.GT. N) GO TO 90 DO 60 J - KP1, N 59

SM - SM + CABS1(Z(J) +WKM*CONJG (A (K, J))) Z(J) - Z(J) + WK*CONJG(A(K,J)) S - S + CABS1(Z(J)) 60 CONTINUE IF (S.GE. SM) GO TO 80 T - WKM - WK WK - WKM DO 70 J - KP1, N Z(J) - Z(J) + T*CONJG(A(K,J)) 70 CONTINUE 80 CONTINUE 90 CONTINUE Z(K) - WK 100 CONTINUE S = 1.0EO/SCASUM(N,Z,1) CALL CSSCAL(N,S,Z,1) C CCC Solve CTRANS(L)*Y - V C DO 120 KB - 1, N K - N + 1 - KB IF (K.LT. N) Z(K) - Z(K) + CDOTC(N-K,A(K+1,K),1,Z(K+1),1) IF (CABS1(Z(K)).LE. 1.OEO) GO TO 110 S - 1.0EO/CABS1(Z(K)) CALL CSSCAL(N,S, Z,1) 110 CONTINUE L - IPVT(K) T - Z(L) Z(L) - Z(K) Z(K) - T 120 CONTINUE S - 1.0EO/SCASUM(N,Z,1) CALL CSSCAL(N,S,Z,1) C YNORM - 1.OEO C CCC Solve L*V - Y C DO 140 K - 1, N L - IPVT(K) T - Z(L) Z(L) - Z(K) Z(K) - T IF (K.LT. N) CALL CAXPY(N-K,T,A(K+1,K),1,Z(K+1),1) IF (CABS1(Z(K)).LE. 1.OE0) GO TO 130 S - 1.0EO/CABS1(Z(K)) CALL CSSCAL(N,S,Z,1) YNORM - S*YNORM 130 CONTINUE 140 CONTINUE S - 1.OEO/SCASUM(N,Z,1) CALL CSSCAL(N,S,Z,1) YNORM - S*YNORM C CCC Solve U*Z - V C DO 160 KB - 1, N K - N + 1 - KB IF (CABS1(Z(K)).LE. CABS1(A(K,K))) GO TO 150 S - CABS1(A(K,K))/CABS1(Z(K)) CALL CSSCAL(N,S,Z,1) YNORM - S*YNORM 150 CONTINUE IF (CABS1(A(K,K)).NE. O.OEO) Z(K) - Z(K)/A(K,K) IF (CABS1(A(K,K)).EQ. 0.0E0) Z(K) - CMPLX(1.OEO,O.OEO) T - -Z(K) CALL CAXPY(K-1,T,A(1,K),l,Z(1),1) 160 CONTINUE C MAKE ZNORM - 1.0 S - 1.0EO/SCASUM(N,Z,1) CALL CSSCAL(N,S,Z,1) YNORM - S*YNORM C IF (ANOR[M.NE. O.OEO) RCOND - YNORM/ANORM IF (ANORM.EQ. 0.OEO) RCOND - 0.OEO RETURN END C******************************************************************c C C SUBROUTINE CGEFA(A,LDA,N,IPVT,INFO) C C C******************************************************************C C C C NAASA 2.1.043 CGEFA FTN-A 05-02-78 THE UNIV OF MICH COMP CTR C INTEGER LDA,N,IPVT (1),INFO COMPLEX A(LDA,1) C C CGEFA FACTORS A COMPLEX MATRIX BY GAUSSIAN ELIMINATION. C C CGEFA IS USUALLY CALLED BY CGECO, BUT IT CAN BE CALLED C DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED. C (TIME FOR CGECO) - (1 + 9/N)*(TIME FOR CGEFA) C C ON ENTRY C 60

C A COMPLEX (LDA, N) C THE MATRIX TO BE FACTORED. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY A C C N INTEGER C THE ORDER OF THE MATRIX A C C ON RETURN C C A AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS C WHICH WERE USED TO OBTAIN IT. C THE FACTORIZATION CAN BE WRITTEN A - L*U WHERE C L IS A PRODUCT OF PERMUTATION AND UNIT LOWER C TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR. C C IPVT INTEGER(N) C AN INTEGER VECTOR OF PIVOT INDICES. C C INFO INTEGER C - 0 NORMAL VALUE. C - K IF U(K,K).EQ. 0.0. THIS IS NOT AN ERROR C CONDITION FOR THIS SUBROUTINE, BUT IT DOES C INDICATE THAT CGESL OR CGEDI WILL DIVIDE BY ZERO C IF CALLED. USE RCOND IN CGECO FOR A RELIABLE C INDICATION OF SINGULARITY. C C LINPACK. THIS VERSION DATED 07/14/77 C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LABS. C C SUBROUTINES AND FUNCTIONS C C BLAS CAXPY, CSCAL, ICAMAX C FORTRAN ABS,AIMAG,CMPLX, REAL C C INTERNAL VARIABLES C COMPLEX T INTEGER ICAMAX, J, K, KP1,L,NM1 C COMPLEX ZDUM REAL CABS1 CABS1(ZDUM) - ABS(REAL(ZDUM)) + ABS(AIMAG(ZDUM)) C CCC Gaussian elimination with partial pivoting C INFO - 0 NM1 N - 1 IF (NM1.LT. 1) GO TO 70 DO 60 K - 1, NM1 KP1 - K + 1 C C FIND L - PIVOT INDEX C L - ICAMAX(N-K+1,A(K,K),1) + K - 1 IPVT(K) - L C CCC Zero pivot implies this column already triangularized C IF (CABS1(A(L,K)).EQ. O.OEO) GO TO 40 C CCC Interchange if necessary C IF (L.EQ. K) GO TO 10 T - A(L,K) A(L,K) - A(K,K) A(K,K) - T 10 CONTINUE C CCC Compute multipliers C T - -CMPLX(1.OEO,O.OEO)/A(K,K) CALL CSCAL(N-K,T,A(K+1,K),1) C CCC Row elimination with column indexing C DO 30 J - KP1, N T - A(L,J) IF (L.EQ. K) GO TO 20 A(L,J) - A(K,J) A(K,J) - T 20 CONTINUE CALL CAXPY(N-K,T,A(K+1,K),1,A (K+1,J),1) 30 CONTINUE GO TO 50 40 CONTINUE INFO - K 50 CONTINUE 60 CONTINUE 70 CONTINUE IPVT(N) - N IF (CABS1 (A (N, N)).EQ. O.OEO) INFO - N RETURN END C C 61

C****************************************** **************************c C C SUBROUTINE CGESL(A,LDA,N,IPVT,B, JOB) C C C*************************************************************************** C C C NAASA 2.1.044 CGESL FTN-A 05-02-78 THE UNIV OF MICH COMP CTR C INTEGER LDA,N, IPVT (1), JOB COMPLEX A(LDA, 1), B (1) C C CGESL SOLVES THE COMPLEX SYSTEM C A * X B OR CTRANS(A) * X - B C USING THE FACTORS COMPUTED BY CGECO OR CGEFA. C C ON ENTRY C C A COMPLEX(LDA, N) C THE OUTPUT FROM CGECO OR CGEFA. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY A C C N INTEGER C THE ORDER OF THE MATRIX A C C IPVT INTEGER(N) C THE PIVOT VECTOR FROM CGECO OR CGEFA. C C B COMPLEX (N) C THE RIGHT HAND SIDE VECTOR. C C JOB INTEGER C - 0 TO SOLVE A*X - B, C - NONZERO TO SOLVE CTRANS(A)*X - B WHERE C CTRANS(A) IS THE CONJUGATE TRANSPOSE. C C ON RETURN C C B THE SOLUTION VECTOR X C C ERROR CONDITION C C A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS A C ZERO ON THE DIAGONAL. TECHNICALLY THIS INDICATES SINGULARITY C BUT IT IS OFTEN CAUSED BY IMPROPER ARGUMENTS OR IMPROPER C SETTING OF LDA. IT WILL NOT OCCUR IF THE SUBROUTINES ARE C CALLED CORRECTLY AND IF CGECO HAS SET RCOND.GT. 0.0 C OR CGEFA HAS SET INFO.EQ. 0 C C TO COMPUTE INVERSE(A) * C WHERE C IS A MATRIX C WITH P COLUMNS C CALL CGECO(A,LDA,N,IPVT,RCOND,Z) C IF (RCOND IS TOO SMALL) GO TO... C DO 10 J - 1, P C CALL CGESL(A,LDA,N,IPVT,C(1,J),0) C 10 CONTINUE C C LINPACK. THIS VERSION DATED 07/14/77. C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LABS. C C SUBROUTINES AND FUNCTIONS C C BLAS CAXPY,CDOTC C FORTRAN CONJG C C INTERNAL VARIABLES C COMPLEX CDOTC,T INTEGER K,KB, L,NM C NM1 - N - 1 IF (JOB.NE. 0) GO TO 50 C C JOB - 0, SOLVE A * X - B C FIRST SOLVE L*Y - B C IF (NM1.LT. 1) GO TO 30 DO 20 K - 1, NM1 L - IPVT(K) T - B(L) IF (L.EQ. K) GO TO 10 B(L) - B(K) B(K) - T 10 CONTINUE CALL CAXPY(N-K,T,A(K+1,K), 1,B(K+1),1) 20 CONTINUE 30 CONTINUE C C NOW SOLVE U*X - Y C DO 40 KB - 1, N K - N + 1 - KB B(K) - B(K)/A(K,K) T - -B(K) CALL CAXPY(K-1,T,A(1,K),1,B(1),1) 62

40 CONTINUE GO TO 100 50 CONTINUE C C JOB - NONZERO, SOLVE CTRANS(A) * X - B C FIRST SOLVE CTRANS(U)*Y - B C DO 60 K - 1, N T - CDOTC(K-1,A(1,K),1,B(1),1) B(K) - (B(K) - T)/CONJG(A(K,K)) 60 CONTINUE C C NOW SOLVE CTRANS(L)*X - Y C IF (NM1.LT. 1) GO TO 90 DO 80 KB - 1, NM1 K - N - KB B(K) - B(K) + CDOTC(N-K,A(K+1,K),1,B(K+1),1) L - IPVT(K) IF (L.EQ. K) GO TO 70 T - B(L) B(L) - B(K) B(K) - T 70 CONTINUE 80 CONTINUE 90 CONTINUE 100 CONTINUE RETURN END C C C ***** ************************************* **************** C C SUBROUTINE CAXPY(N,CA,CX,INCX,CY,INCY) C C C****************************************************************** C C C NAASA 1.1.014 CAXPY FTN-A 05-02-78 THE UNIV OF MICH COMP CTR C C CONSTANT TIMES A VECTOR PLUS A VECTOR. C JACK DONGARRA, LINPACK, 6/17/77. C COMPLEX CX(1),CY(1),CA INTEGER I,INCX,INCY,IX,IY,N C IF(N.LE.O)RETURN IF (ABS(REAL(CA)) + ABS(AIMAG(CA)).EQ. 0.0 ) RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GOTO 20 C CCC Code for unequal increments or equal increments CCC Not equal to 1 C IX - 1 IY - 1 IF(INCX.LT.O)IX - (-N+1)*INCX + 1 IF(INCY.LT.O)IY - (-N+1)*INCY + 1 DO 10 I - 1,N CY(IY) - CY(IY) + CA*CX(IX) IX - IX + INCX IY - IY + INCY 10 CONTINUE RETURN C CCC Code for both increments equal to 1 C 20 DO 30 I - 1,N CY(I) - CY(I) + CA*CX(I) 30 CONTINUE RETURN END C ' C C******************************************************************c COMPLEX FUNCTION CDOTC(N,CX,INCX,CY,INCY) C*******************************************************************c C C C NAASA 1.1.012 CDOTC FTN-A 05-02-78 THE UNIV OF MICH COMP CTR C C FORMS THE DOT PRODUCT OF TWO VECTORS, CONJUGATING THE FIRST C VECTOR. C JACK DONGARRA, LINPACK, 6/17/77. C COMPLEX CX(1),CY(1),CTEMP INTEGER I, INCX, INCY, IX, IY, N C CTEMP - (0.0,0.0) CDOTC - (0.0,0.0) IF(N.LE.O)RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GOTO 20 C CCC Code for unequal increments or equal increments CCC Not equal to 1 C IX - 1 IY - 1 IF(INCX.LT.O)IX - (-N+1)*INCX + 1 IF(INCY.LT.O)IY - (-N+1)*INCY + 1 DO 10 I - 1,N 63

CTEMP - CTEMP + CONJG(CX(IX))*CY(IY) IX - IX + INCX IY - IY + INCY 10 CONTINUE CDOTC - CTEMP RETURN C CCC Code for both increments equal to 1 C 20 DO 30 I - 1,N CTEMP - CTEMP + CONJG(CX(I))*CY(I) 30 CONTINUE CDOTC - CTEMP RETURN END C C C***************************************************************** C C SUBROUTINE CSCAL(N,CA,CX,INCX) C C C******************************************************************c C C C NAASA 1.1.019 CSCAL FTN-A 05-02-78 THE UNIV OF MICH COMP CTR C C SCALES A VECTOR BY A CONSTANT. C JACK DONGARRA, LINPACK, 6/17/77. C COMPLEX CA,CX(1) INTEGER I,INCX,N,NINCX C IF(N.LE.) RETURN IF(INCX.EQ.1)GOTO 20 C CCC Code for increment not equal to 1 C NINCX - N*INCX DO 10 I - 1,NINCX,INCX CX(I) - CA*CX(I) 10 CONTINUE RETURN C CCC Code for increment equal to 1 C 20 DO 30 I - 1,N CX(I) - CA*CX(I) 30 CONTINUE RETURN END C C C****************************************************************** C C SUBROUTINE CSSCAL(N,SA,CX,INCX) C C C****+*************************************************************c C C C NAASA 1.1.018 CSSCAL FTN-A 05-02-78 THE UNIV OF MICH COMP CTR C C SCALES A COMPLEX VECTOR BY A REAL CONSTANT. C JACK DONGARRA, LINPACK, 6/17/77. C COMPLEX CX(1) REAL SA INTEGER I,INCX,N,NINCX C IF(N.LE.O)RETURN IF(INCX.EQ.1)GOTO 20 C CCC Code for increment not equal to 1 C NINCX - N*INCX DO 10 I - 1,NINCX,INCX CX(I) - CMPLX(SA*REAL(CX(I)),SA*AIMAG(CX(I))) 10 CONTINUE RETURN C CCC Code for increment equal to 1 C 20 DO 30 I - 1,N CX(I) - CMPLX(SA*REAL(CX(I)),SA*AIMAG(CX(I))) 30 CONTINUE RETURN END C C C***************************************************************** INTEGER FUNCTION ICAMAX(N,CX,INCX) C*****************************************************************c C C C NAASA 1.1.021 ICAMAX FTN-A 05-02-78 THE UNIV OF MICH COMP CTR C C FINDS THE INDEX OF ELEMENT HAVING MAX. ABSOLUTE VALUE. C JACK DONGARRA, LINPACK, 6/17/77. C COMPLEX CX(1) REAL SMAX INTEGER I, INCX, IX,N COMPLEX ZDUM 64

REAL CABS1 CABS1(ZDUM) - ABS(REAL(ZDUM)) + ABS(AIMAG(ZDUM)) C ICAMAX - 1 IF (N.LE.1) RETURN IF(INCX.EQ.1)GOTO 20 C CCC Code for increment not equal to 1 C IX - 1 SMAX - CABS1(CX(1)) IX - IX + INCX DO 10 I - 2,N IF(CABS1(CX(IX)).LE.SMAX) GO TO 5 ICAMAX - I SMAX - CABS1(CX(IX)) 5 IX - IX + INCX 10 CONTINUE RETURN C CCC Code for increment equal to 1 C 20 SMAX - CABS1(CX(1)) DO 30 I - 2,N IF(CABS1(CX(I)).LE.SMAX) GO TO 30 ICAMAX - I SMAX - CABS1(CX(I)) 30 CONTINUE RETURN END C C C******************************************************************* REAL FUNCTION SCASUM(N,CX,INCX) C****************************************************************c** ******** *** C C C NAASA 1.1.010 SCASUM FTN-A 05-02-78 THE UNIV OF MICH COMP CTR C C TAKES THE SUM OF THE ABSOLUTE VALUES OF A COMPLEX VECTOR AND C RETURNS A SINGLE PRECISION RESULT. C JACK DONGARRA, LINPACK, 6/17/77. C COMPLEX CX(1) REAL STEMP INTEGER I,INCX, N,NINCX C SCASUM - O.OEO STEMP - O.OEO IF (N.LE.0) RETURN IF(INCX.EQ.1)GOTO 20 C CCC Code for increment not equal to 1 C NINCX - N*INCX DO 10 I - 1,NINCX,INCX STEMP - STEMP + ABS(REAL(CX(I))) + ABS(AIMAG(CX(I))) 10 CONTINUE SCASUM - STEMP RETURN C CCC Code for increment equal to 1 C 20 DO 30 I - 1,N STEMP - STEMP + ABS(REAL(CX(I))) + ABS(AIMAG(CX(I))) 30 CONTINUE SCASUM - STEMP RETURN END 65

Appendix B Program Listing for the Quasi-Analytical Solution The quasi-analytical solution proposed by Senior and Volakis [1] as described in Section 3 was programmed for solution, as listed in the program ANAGAP.FTN below. The subroutines used by this program are listed in GAPSUB.FTN in Appendix A. As with the program for the exact solution GAPSCAT, the user is prompted for the following: the polarization of the incident field, angle of incidence, angle of far field observation, the relative permittivity er of the gap cavity, the shape of the gap, and the number of iterations for calculating the far field amplitude versus gap depth. The choice of shapes and dimensions requested for the gap are according to Figure 2. The input impedance of the gap as a parallel plate waveguide is calculated according to the specified shape. For the L- and T-shaped gaps, (35) is used, given the other necessary expressions as contained in Section 3. The input impedance of the rectangular gap is given by (39), and for the Vshaped gap, (40) is used. These expressions are for the H-polarization case. As mentioned previously, for the E-polarization case, the propagation constant k1 is replaced by ikp and the characteristic impedance Z1 by -iZp.r /p, where p is given by (41). The desired effective surface impedance 11 is then calculated according to (33). For H-polarization, the far field amplitude PH is calculated from (22). PH is a function of KH(a) given by (30). The argument a is a function of the effective surface impedance as given in (21). For the E-polarization, the far field amplitude PE is calculated from (27), where KE(b) is given by (31). The argument b is a function of the effective surface impedance as given in (26). 66

c C ANAGAP.FTN C C C This FORTRAN program computes the far field scattering due C to a narrow gap of specified shape in an infinite ground 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, and number of C iterations with respect to gap depth. 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 MODBESS 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 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 logical Epol,Lossy common /prompts/ EorH,phio,phi,er,igap,wStp,dStp,w,d,noS, & q,maxC,noIter 1 format(il) 2 format(i5) 3 format(gl6.8) 4 format(al) 5 format(i3) 6 format(2g16.8) 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 ur-cmplx(1.,0.0) gam-0.5772157 iprg-2 c...Setting default values 10 EorH-1 phio-90.0 phi-90.0 er-cmplx(1., 0.0) w-0.15 d-0.5 noS-3 maxC-0.01 noIter-30 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 (l)/d if(aImag(er).ne. 0.0) Lossy-.true. dmin-0.05 67

if(Epol) dmin-0.025 dmax-d if(noIter.ne. 1)then dstep- (dmax-dmin)/ (nolter-1) d-dmin endif DO 700 iter-l,noIter c********** **********************GpIpdne ** ********************* ********************** Gap Impedance ************************ C*** *************** *** *********** ***************************** c...Complex propagation constant kl and characteristic impedance c Z1 of the T-line model if(Epol)then kl-ci*k*csqrt ((l./2/w)**2-er*ur) Zl —ci*Zo*ur/csqrt ((1./2/w) **2-er*ur) else Zl-Zo*csqrt (ur/er) kl-k*csqrt(ur*er) endif if(igap.eq. 1)then c RECTANGULAR ETA —ci*Zl*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 ( (./2/d2) **2-er*ur) Zc —ci*Zo*d2*ur/csqrt ( (1./2/d2) **2-er*ur) else Zc-Zo*d2*csqrt (ur/er) kc-k*csqrt (ur*er) endif if(igap.eq. 4)then c L-SHAPED Zi —ci*Zc*ctan(kc*w2) X-k*Zc*wl Bl-k/Zc*d2* (d2/(d2+wl)) * (.-2./pi*log(2.)) B2-k/Zc*d2* (wl/(d2+wl) ) * (.-2./pi*log (2.) ) ZL-(Zi-ci*X* (l-ci*B2*Zi)) & / ((1-B1*X) * (-ci*B2*Zi)-ci*Bl*Zi) else c T-SHAPED Zi —ci*Zc* (ctan (kc*w2) +ctan (kc*w3)) X-2*k*Zc*wl Bl-k/Zc*d2* (d2/(d2+wl)) *0.7822 ZL-(Zi-ci*X) / (-ci*Bl*(Zi-ci*X)) endif ETA-Z1* (ZL-ci*Zl*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 cHank (carg, 3, HO, H1) carg-1. else if(Epol)then rarg-Real (kl/ci) *d call ModBess (rarg, I0, I1) HO-IO H1-I1 carg-ci else rarg-Real (kl) *d call Hankzl(rarg,2,H0,H1) carg-1. endif endif ETA —ci*Zl*Real (H) /Real (HO) *carg endif c write(2,*) d,cabs(ETA) c****************************************************************** c****************** Far Field Amplitude *********************** if (Epol) then b —ci*k*w/2*Zo/ETA Ke-0.62/(b+1.15) * (b+4.08) * (b+7.26) * (b+10.37) & * (b+13.43) * (b+16.46) & /((b+4.27)*(b+7.37)*(b+10.45)*(b+13.49) & *(b+16.5)) Psca —ci*pi/4*(k*w)**2*sin(phio)*sin(phi)*Ke else a-ci*2./k/w*Zo/eta Kh —l/ (pi/2*a+0.1+log (2.) ) Ao-log (k*w/4) +gam-ci*pi/2 Psca-ci*pi*Kh/(l+Ao*Kh) endif 68

c...Outputting the far field magnitude and phase print,' d - ',d print *,' Analytical: IPscal - ',cabs(Psca), ' arg Psca - ',180/pi*(atan2(aImag(Psca),Real(Psca))) write(3,*) d,cabs(Psca) write(4,*) d,180/pi*(atan2(aImag(Psca),Real(Psca))) d-d+dstep 700 continue print *,' Again (1-yes)? read(*, 1)ians if(ians.eq. 1) GOTO 10 800 call exit END 69

UNIVERSITY OF MICHIGAN 3 9015 03524 4550II 3 9015 03524 4550