U-M Technical Report 389604-7-T Electromagnetic Scattering by a Perfectly Conducting Patch Array on a Dielectric Slab Jian-Ming Jin and John L. Volakis Radiation Laboratory Department of Electrical Engineering and Computer Science The University of Michigan Ann Arbor, MI 48109 March 22, 1989 Prepared for: General Dynamics - Fort Worth Div. P.O. Box 748 Ft. Worth, Texas 76101

REPORT SUMMARY This report describes the theoretical development and numerical implementation of a methodology for the electromagnetic characterization of an infinite and truncated periodic array of perfectly conducting patches on a dielectric slab. More specifically, it contains the following: * Derivation of the dyadic Green's functions for an infinite dielectric slab * Development of the equation for computing the currents on an infinite periodic patch array supported by an infinite slab * Application of a Conjugate Gradient (CG) - FFT method for the solution of the equation * Exact solution of the reflection by and transmission through an infinite periodic patch array on a dielectric slab * Approximate solution of the scattering by a truncated patch array on an infinite dielectric slab * Approximate solution of the scattering by a truncated patch array on a finite slab. Of particular interest is the application of the CG-FFT method and the approximate solution for the truncated array scattering. Unlike traditional matrix inversion techniques, the CG-FFT solution method allows a substantial computer

storage economy and is thus suitable for practical size geometries. The approximate solution for the truncated array scattering is compared with available exact data and is found to be reasonably accurate.

Table of Contents 1. INTRODUCTION.....................................................1 2. DYADIC GREEN'S FUNCTIONS......................................4 G eom etry.............................................................. 4 Maxwell's Equations in Spatial Form....................................6 Maxwell's Equations in Spectral Form.................................. 7 Solutions in Spectral Form.............................................. 7 Dyadic Green's Functions............................................... 9 3. FORMULATION FOR PATCH CURRENTS.......................... 12 Geometry of an Infinite Array......................................... 12 Floquet's Theorem.................................................... 15 Scattered Field........................................................ 16 Equation for Patch Currents.......................................... 17 4. NUMERICAL IMPLEMENTATION................................... 19 Discretization of Geometry............................................ 19 Discretization of Currents................................. 19 Application of Galerkin's Technique................................... 22 CG-FFT Technique................................................... 24 5. REFLECTION, TRANSMISSION AND DIFFRACTION OF A PERIODIC A R R A Y........................................................ 26 Reflection and Transmission Coefficients............................... 26 Bragg D iffractions..................................................... 27 N um erical Results...................................................... 28 1

6. SCATTERING BY A TRUNCATED ARRAY.........................37 Bidirectionally Truncated Array.......................................37 Unidirectionally Truncated Array.......................................42 Numerical Results.....................................................43 7. TRUNCATION EFFECT OF DIELECTRIC SLAB...........51 Physical Optics Approximation........................................51 Numerical Results.....................................................56 8. CONCLUSION.............................61 9. REFERENCES.......................................................63 10. APPENDIX..........................................................65

Index to Figures 2.1 Geometry of a dielectric slab............................................ 5 3.1 Geometry of a periodic patch array on a dielectric slab. (a) Top view. (b) Side view............................................................. 13 4.1 Discretization of a square, circular, and triangular patch................20 5.1 Three types of cell geometries. (a) A square patch cell. (b) A circular patch cell. (c) A cross-shaped patch cell...................................... 29 5.2 (a) Reflection and (b) transmission coefficients vs frequency for an array of square patches with various dielectric slab thicknesses.................. 30 5.3 (a) Reflection and (b) transmission coefficients vs frequency for an array of square patches for three incident angles for TM incidence.............. 32 5.4 (a) Reflection and (b) transmission coefficients vs frequency for an array of square patches for three incident angles for TE incidence............... 33 5.5 (a) Reflection and (b) transmission coefficients vs angle of incidence for an array of circular patches............................................... 34 5.6 (a) Reflected, (b) transmitted and (c) dissipated power vs frequency for an array of cross-shaped patches with various lossy dielectric slabs...... 35-36 6.1 Far field patterns in the k = 180~ half-plane for a unidirectionally truncated array having nine patches in the x direction......................... 44-45 6.2 Bistatic radar cross section for a unidirectionally truncated array having seven patches in the x direction....................................... 46 6.3 Bistatic radar cross section patterns for a 9 x 9 cross-shaped patch array at f = 5 G H z........................................................... 48 iii

6.4 Bistatic radar cross section patterns for a 9 x 9 cross-shaped patch array at f = 9GHz...............................49 6.5 Bistatic radar cross section patterns for a 9 x 9 cross-shaped patch array at f =20OGHz...............................50 7.1 Bistatic radar cross section patterns for a 26 cm x 26 cm dielectric slab at 9 GHz...................................57 7.2 Bistatic radar cross section patterns for a 26 cm x 26 cm dielectric slab at 20 GHz.................................58 7.3 Bistatic radar cross section patterns for a 9 x 9 cross-shaped patch array on a 26 cm x 26 cm dielectric slab at 9 GHz.................59 7.4 Bistatic radar cross section patterns for a 9 x 9 cross-shaped patch array on a 26 cm x 26 cm dielectric slab at 20 GHz................60

Chapter 1 Introduction The subject of electctromagnetic characterization of periodic structures has been of continuing interest over the recent years. The problem that has received most attention is that of free-standing arrays. In typical applications, however, dielectric slabs and modify their characteristic spectral response [1]-[3]. One approach for the treatment of such arrays is to first pursue individual studies of the free-standing array and dielectric slab in isolation. The composite structure can then be analyzed by cascading the generalized scattering matrix parameters found for each problem [4], [5]. Alternatively, one can treat the composite structure directly by including the effect of the dielectric slab through the use of the dyadic Green's function associated with the slab and this approach is adopted in this study. The formulation for scattering by periodic structures usually invokes Floquet's theorem leading to an equation for the patch currents upon enforcement of the required boundary condition on the array elements. Among various numerical methods for solving the resulting system of equations, the method of moments has 1

been widely used [1]-[3], [6], [7]. The major drawback of this method, however, is a need to explicitly generate and store a square impedance matrix whose size could exceeds the memory capacity of modern computing facilities. On the other hand, an iterative approach such as the conjugate gradient method when combined with the fast Fourier transform [8] avoids the explicit generation of the square matrix, allowing the study of practical size geometries. In the first part of this report, after a derivation of the pertinent dyadic Green's functions (Chapter 2), an equation is formulated for the patch currents on an infinite periodic array of perfectly conducting patches residing on a dielectric slab (Chapter 3). The Galerkin's procedure is then employed using subdomain roof-top basis functions to develop a system of equations for the patch currents. The system is subsequently solved via the conjugate gradient method in conjuction with the fast Fourier transform technique (Chapter 4). A computation of the radiation by the patch currents provides a solution for the reflection and transmission of the infinite periodic array and this is discussed in Chapter 5. In practice, an array is always of finite extent and the second part of the report deals with the characterization of the truncated array. An exact analysis of the finite array, however, is time consuming, costly and sometimes impractical to achieve for engineering design purposes. As a result, the scattering by the truncated array is considered here approximately (Chapter 6). In particular, the infinite array currents are integrated over the extent of the truncated array to evaluate the scattering pattern. The accuracy of this approximation is examined 2

and discussed in relation to some available exact data. In actual applications, both the patch array and the supporting dielectric slab are of finite extent. A method to characterize the scattering by the truncated array on a finite slab is developed in Chapter 7. In particular, a physical optics approximation is employed to compute the field diffracted by the slab edge which is then added to that contributed by the truncated patch array. 3

Chapter 2 Dyadic Green's Functions This chapter describes the derivation of the Green's functions associated with a dielectric slab. Though there are such Green's functions available in the literature, they are usually presented in an incomplete form and often contain typographical errors. It is, therefore, instructive to derive them here since their use is crucial in our study. 2.1 Geometry The geometry we consider here is illustrated in Figure 2.1 where a unit infinitesimal electric current source oriented in the x direction is located on the upper surface of a dielectric slab. The coordinate system is chosen so that the x-y plane is coincident with the lower surface of the slab. The slab is assumed to have a thickness d, a permittivity e,.e, and a permeability of the free-space. The current source is then located at (x',y',d). In the following, we begin with Maxwell's equations to find the resulting Es, Ey, and E, fields. 4

dielectric slab (x', y', d) '%,. ell, x T 0 Figure 2.1: Geometry of a dielectric slab. 5

2.2 M~axwell's Equations in Spatial Form The source-free Maxwell's equations are V xE =-jwlzoH (2.1) V xH =jwf-iOE (2.2) where an eiwt time convention is employed and w is the angular frequency. In the above, e, = El= C. inside the dielectric slab and ej = C2= 1 in air. Equations (2.1) and (2.2) can be easily transformed to the two vector wave equations V XV xE - kofiE= (2.3) V XV xH - kof4H= (2.4) and since V.- (4ej) =0 and V.- H = 0, they can be further reduced into six scalar wave equations. They are,V2E~x+ k02.ExE 0 (2.5) V2Z+k~ciHz 0 (2.6) for the z components and ( 2 C+9 2 (2.7 (eik +x = +z iwu Hz (2.8) 2 '9 ' Ieik0 + aZ) Hv.~6Hz - jwf4EoT.-Ez (2.10) for the other components. 6

2.3 M~axwell's Equations in Spectral Form To solve (2.5) and (2.6), we introduce the Fourier transform pair F(kx, ks,,z) =L J F(x, y, Z) e-k e-kYvdk3.,,dky, (.1 f00f00,, F (x, y, z) - __F(k~,kv,,z)e-1kzx'ejAYydxdy (2.12) where F represents any component of electric and magnetic fields. In the transform domain, equations (2.5) and (2.6) can be written as CIZ2+ kiz 0(2.13) 2z 0 2ff+ kiHo (2.14) ~z2 where k? = k2= e, k02-_2 within the slab, k? = k2= ko-/32 in air and /2 =2 k+k 2 Similarly (2.7)-(2.10) become Ex # 9 x+ # H (2.15) - jky 09 w/sok.f Ey #219Z z - 2 Hz(2.16) jk2,0-.-Hz w#2 Ek (2.17) /32+ Ez. (2.18) 2.4 Solutions in Spectral Form The solution of (2.13) and (2.14) must satisfy the radiation condition and as a result the spectral field components take the form E,= Ae kz ink2 < 0, for z > d (2.19) 7

Hz- Be3I2Z, ink2 <0, for z> d (2.20) E2.=Ccoskiz+Dsinkiz, for d>z>0 (2.21) Hz=Esinklz+Fcosklz, for d>z>0 (2.22) Ez= Gejk2z, link2 < 0, for z < 0 (2.23) -z He jk2Z link2 < 0, for z < 0 (2.24) where A, B, C, D, E, F, G and H are the unknown constants to be determined by applying the boundary conditions. The boundary conditions for this problem are E: - E,-=0, for z=O0and z= d (2.25) E+-E17,=0, forz=Oandz=d (2.26) H,+- H,,j= 0, for z=O0and z= d (2.27) H+ - HIf=0, for z=O0and (2.28.) H+j - H1- = - 6(x -x')6(y -y'), for z d (2.29) where the superscript "+" denotes the associated quantities on the plane z = 0 + u or z = d+ua and "-" denotes the corresponding quantities on the plane z = 0-u, or Z = d - a, in which a is an infinitesimal number. Upon application of (2.25)- (2.29) we find A =kj(k, kx Z e-Jk1)x'e-jSYee (2.30) ei k -k~x' eikvv'ejk2d (2.31) C = — kAinkek~ ~ (2.32) jsnk1d + ek2 cos A~ 8

D = j C (2.33) E = k k2e-Jkd k2 sin kd - jkl cos kid (2.34) F = -iE (2.35) G = rC (2.36) H=F (2.37) where r rk2 sin k1d - jkl cos kd (2.38) f k2 cos k1d + jkl sin k1d k2cos kd + jkl sin kd (239) k2 2 sin k1d - jkl cos kd (2.39) Substituting now the constants into (2.19)-(2.20), we find, from (2.15)-(2.16), that for z > d = _ zo jkk2kr- _ _2k + k e-jk.'e-jk'e1e-jk2(zd) (2.40) ko 2 (rk2r-jk- jkrJ2- k2 Ey - 2 (erk - A jkl - Jb) e e e * (2.41) The fields E, and E. for z > d can then be found via inverse Fourier transformation as described in (2.12). 2.5 Dyadic Green's Functions The component of the total radiated field transverse to z can be written as ET(X, y, z) = / O(x, y, zIx', y') * 7(X, y')dx'dy (2.42) 9

where G= xxGxs+yxGyx+yyGy `+xyGxv is the pertinent dyadic Green's function and J =- Jxx + Jyy is a two-dimensional radiating surface current residing on the surface z = d. In accordance with the development in the previous section, G zo _ jklk2k2 + ko2k2 472kJooJ /oo 2 r,.k2r- jkl jklF2-k2, e-jk2(z-d) eikx(x-x') eikY(Y-Y') dkdky (2.43) yX 44r2koo0o oJ o 2 erk2r,- jkk- jkr2 - k2. e-jk2(z-d) e jk (x-x') edk(y-Y')dk)dky (2.44) valid for z > d. Introducing now the definitions Zo jk k2k k okv G k2 =- Ek2erl-jk, jkr2 - k2 (2.45). _ Zo ( jkik2k(k. _ kok,.k ) ko0/2 erk2r - jkl jkir2 - k2 we may rewrite (2.43)-(2.44) as Gx = __ Ge e jk2(z-d)ek,(-)eiky(-')dkdkd (2.47) 41r J-oo - oo Gyv = (4d2 | G( e-(-)e* eJ(-Y')dk dk (2.48) 47r2 0 <-oo where we observe that G.s and G,, are not transforms of G. and G., but simply functionals introduced for convenience. The other two components of G, namely Gyy and G., can be found by a simple interchange in coordinates. They are G = 4r2 J| Gyye-jk2 (z)e((-)ejky(y-')dkxdk (2.49) GY= 2 L L Gae-jk2(z-dek (-') eik(-')dkdk (2.50) 10

where ' zo - jkl2k2 k2k2 yy ko/32 erk2Fi - jkl jkl2 - k2F Gry = Gya (2.51) (2.52) Equations (2.47)-(2.50) can be written in a compact form as - =00 _-00 - d' G = J G (ks, ky)e-jk2 (z-d) ejkx(x-x') ejk(y-')dkdky 471' Joo 2 O in which G = xxGx + yxGy + yyGyy + xyGxy. (2.53) (2.54) The above is the Green's function valid in the region z > d. The Green's function for z < 0 can also be derived in similar manner and is given by =1 1 roo koo.= G =whereG (k ky)eJ2Ze )e y')dkxdk (2.55) where '/ = Zo xY kop2 G - vY koP2 ' eyx jerkik2k2 jk2kik2 1 7i(erk2Fr - jkl) 72(jkil2 - k2)J jerkk2k ky + jkoklkky 1 71(erk2F1 - jkl) -2(jkiF2 - k2) (2.56) (2.57) (2.58) (2.59) with Y1 = jkl sin k1d + erk2 cos kid 72 = k2 sin k1d - jkl cos k1d. (2.60) (2.61) 11

Chapter 3 Formulation for Patch Currents In this chapter, Floquet's theorem is employed in conjuction with the Green's functions derived earlier to develop an equation for the electric currents generated on an infinite periodic patch array due to a plane wave excitation. 3.1 Geometry of an Infinite Array The geometry of an infinite periodic patch array is shown in Figure 3.1, where a periodic array of conducting patches with periodicities T. and Ty in the x and y directions, respectively, is residing on the surface of a dielectric slab. To be consistent with the notation used in the last chapter, the thickness of the slab is assumed to be d and the chosen coordinate system is shown in Figure 3.1. The excitation field to be considered is a plane wave, but other excitation can be employed if expressed as a superposition of plane waves through the use of the Fourier integral. Further, although the formulation to be presented is in general valid for arbitrary polarizations, particular attention is concentrated in the case of transverse electric (TE) and transverse magnetic (TM) incidences. For the TE 12

(a) _______________ d N. x 1* 0 (b) Figure 3.1: Geometry of a periodic patch array on a dielectric slab. (a) Top view. (b) Side view. 13

case, the incident electric field can be written as (, z) = (-x sin i' + y cos i)Eoeekiikozcos~ ' (3.1) while for the TM case, it is given by (p, z) = (x cos 0' cos q~ + y cos 0' sin /9 - z sin Oi)Eoei Pe3JkozcO' (3.2) where p xx + yy and k/ = X4 + ~yky in which k = ko sin 0 cos O and k= ko sin O sin qY. If the conducting patches are absent, a reflected and transmitted plane wave will result. Since these fields are required later in the formulation, we proceed to derive them for the TE and TM cases. For the TE case, the reflected field can be expressed as T'(Q, z) = (-x sin i + y cos qi)EoRTEei2kodcosieJkT e-kocos0 (3.3) in which RTE is the corresponding reflection coefficient given by RTE COS - j. - sin2 (34) cos 0i + j /e - sin2 0iCo where j = J /e- sin2 i sin(kod /e,. - sin2 i) + cos i cos(kode - sin2 i) jr - sin2 i cos(kod/er - sin2 0i) - cos Oi sin(kod er sin2 oi) The transmitted field is found as t(p, z) = (-: sin i' + y cos i)EoTTEe'kodc~eteiikpejkzcos(3.5) 14

in which TTE is the corresponding transmission coefficient given by TTE j i Ve.-sin2i( + RTE) j/ - sin2 0i cos(koder - sin2 Oi) - cos Oi sin(kod/r -sin2 0i) Similarly, for the TM case the reflected field is TE (, z) = -(x cos 0i cos qi + y cos 0i sin q>i + z sin 0i). EoRTMe j2ko d cos ' ejk e- jkoz cos' (37) where again RTM is the reflection coefficient given by RTM = e cos - j,. - sin2 iC1 (3.8) r Cos 0t + j Er-sin2 C with j Er- sin2 sl in(kod r - sin2 ) + e r cos ' cos(kod r- sin2 ') j V/e - sin2 i cos(kod jr - sin2 0i) - e cos 9i sin(kod e - sin2 Oi) The transmitted field is found as (Et, z) = (x cos 0i cos q' + y cos 9' sin t - z sin 0).E TMeJkd cods ejk' jkoz cos (39) in which TTM is the corresponding transmission coefficient given by TTM= j / - sin2'(1+R ) (3.10) j '/e'- sin2 'i cos(kod er - sin2 9i) - er cos Oi sin(kod er - sin2 Oi) 3.2 Floquet's Theorem In accordance with Floquet's theorem, the electric current supported by the array residing on the plane z = d may be deduced by that on a single patch. 15

Namely, the array current distribution can be written as 7(X, y) =(x, Y)e (xkY)(.1 where j (x, y) is a two-dimensional periodic function. It can, therefore, be expanded in a Fourier series 00 00A j X Y)j Jpq Vpq (X, Y) (3.12) p=-oo q=-oo where Jp '~J (x', y')?p*q(x', y')dx'dy' (3.13)?i/pq (X, Y) = e i(kzTpX+kyqY) (3.14) in which Sp, denotes the area of a single patch, kxp = 2 irp/Tx, kyq =2irq/Ty, and p and q are integers representing the Floquet modes. 3.3 Scattered Field From Section 2.5, the transverse part of the scattered field for z > d as produced by the current J5 is T4x) yz) =Z L L. G._J,(x', y') dx'dy' (3.15) and by employing (2.53) in conjuction with (3.11), it becomes y, z) 100 L, 1L J00 =eik2(z-d)e'kx(x-x')eiky(YY)kd 00 00 A Z pqei, PxkYqY) dx~dyI (3.16) p=-oo q=-oo 16

where k' k.,v+ k' and k'q kyq + k'. Interchanging the order of integration, we now obtain T~(x~yz =e- > jq2(z-) p=-oo q=-oo - LL~ ~:~~o~pqjk,2(z-d)ejkxxejk~yy -P kx)Skq - ky)dkxdky. (3.18) Further, by interchanging the orders of summation and integration we obtain E' xy, z E E (k'k') pqz-) yz)=0XP7 yq) Jp (3.19) p=-oo q=-oo upon evaluation of the integrals, where G~k k',)i bandbyrpaigk n k. with k' and k q, 7respectively, in (2.45)-(2.46) and (2.51)-(2.52). 3.4 Equation for Patch Currents The equation for the electric currents can be obtained by applying the boundary condition on the patch. This demands a vanishing tangential electric field on the patch, that is T T+T+ = on the patch. (3.20) Upon substitution of (3.19) into (3.20) we then obtain p=-oo q=-oo 17

for 0 < 9' < ir/2. In the case when ir/2 < O' < ir 7, + ET in (3.21) should be replaced by the transverse component of the transmitted field. We note that an alternative expression for (3.21) is 00 00 ~ (k pk' q) ipq OPpq(X,Y) =B (3.22) p=-oo q=-oo where B is given by = - [7(x) Y, d) + TT(x, y, d)] e-(kx 4Y (3.23) and can be shown to be invariant to (x, ). Specifically, for the TE case =(-X'sinq5' + ' cos 0t)Eo(1 ~ RTE)eikd ~(3.24) and for the TM case =X (5cosO' cosqt + ' cosO' sin4 'Eo (1 - RTM)eikd G.(3.25) Solving (3.22), one obtains the current distribution on the patches. 18

Chapter 4 Numerical Implementation In this chapter, we describe a numerical procedure for solving (3.22). The procedure employs the Conjugate Gradient-Fast Fourier Transform (CG-FFT) technique which appears to be the natural and best choice for the solution. 4.1 Discretization of Geometry The first step of the CG-FFT solution is to discretize the geometry. For this, we first divide the region of one period into small rectangular cells of dimension Ax x Ay. The surface of the patch is then generated based on the premise that any planar area may be approximated by a collection of the rectangular cells. Figure 4.1 shows three examples of such a discretization as applied to square, circular and triangular patches. 4.2 Discretization of Currents Roof-top functions are used below to expand the patch currents. These allow a better representation of the currents and can provide a faster converging series 19

Lr _ = 1 1 - 1 '" LLill'lil1 l1l II 1 1 1 p I I _^ _1 J I I I i Ii i I _ | | | |.. _ _ 1 __ _i_ I I I I I I 1 1 i...- t.tt X X X X X.tl ~ t.-..........tt t J~ t i l.t.......* ** -. _ _ I _ tr I _ _t _ _ _ _ _ _ _ _ _ _ _ I... *... * *XXt t t~ l t X.... i..1 _ _. 5~~~~~~~~~~~~16 1414=_1 ______.....!.S1*161 1 1 XII ___1*r L______ *to* 191:019191,1-91 ~ ~ lefe Waiteee ~ eeej ~!~~~~~ ~~~~~ ~ ~~~ ~9r l'9i~L _:::::::'.i:^. S..;- Sii*;::::111:::: y _.'.'.'.'.'. I I I _' I I I I I; _.;1 _~~~~~~~~~~~~~~~**0 I I I Ii I I~f I= 1 T =*-*2 _ _ _! *1S 1 * I*In IsI* l* T I I I I I _- S. I,S S. I S1114 " I S. I. * Sb * 11 1 1 1*|| _ ~e~e~~~~e!~~let _ I I _ _ * 1 **1\* *lllellllulllullblgliiii__*iiiii!i,*nii**2 *^^111*val11 Figure 4.1: Discretization of a square, circular, and triangular patch. _~~_* 1 ~ ~r.~rI. ~T ~*.... ~ ~f ~e~ ~ 4 ~ 1,1~4 ~ _ ~ ~ ~ eL *|sLIsL~SssS*-SSSS@_ e_ ~ S e ~ ~~* ~e~~*~eul l ~~~el- ~lle~~~lie~u |||I **S S*S*S 1 _ _ _ _ ~i e* eiXe~~~lieu lel~ui uui uuef1 1I___ Le!L _ _ e e e ~Lj__i 1{111111 - ~ i i eI*~~uuuuus ei ~uui ~ e | e ee ~1~eee Ire _~._2fi_2si~XEZ11 _ e _ * _ elee lee ~ e llle~lllll ~ e e ee ~ rrr _ _ _ _ HT {iT Inf l l elflllflllll~lllIII__;.S __ i1111 _ _ e i ~.rr rrl lelrllllllrrlllrWX^xI_____ _ss.s l 1111111|r _~~~~~~~ ~~~ _~ T@ IIIlIII*||*lui ______*fsfi___i* 1 _ _ l _ e ~T Le ~e~[~ ~[~ ~:e ~5||||f- ~~i@|X ~ e@|s Ii@i 1______ ___*,... _e e e-[e ~ 'l leD _ e~~~~ _~~iil _e ~l[X eleelee elee_ e _e el lleee _ eee mll~~~~~~~~~~~ie __ _ ~ _~e~~e- ~ ~ ere ej~e!1 _T~e~ ~~ ee~e~~ ~ ~~ e~ ~~ ~~~~e = ~ ~i e ~~et~eeiee_ ~e __ _e~~e ~~ ~~~ ~leee~ Fi:gure 4.1]:DiscretiJzation of a square, circular, and triangular patch. 20

than pulse basis. The current distribution can then be written as M/2-1 N/2-1 jx (X,y) = >J Z jxmnAm+1/2(X)Hln(y) (4.1) m=-M/2 n=-N/2 M/2-1 N/2-1 jy(X' Y) = E >j jymn~lm(x)An+i12(Y) (4.2) m=-M/2 n=-N/2 where I x- mAx I/Ax, Ix - mAxI < Ax Am (x) =(4.3) { 7Ix - mAxI < Ax/2 H~m(X) =(4.4) M and N denote the number of cells along the x and y directions, respectively, and Jmn = XJxmn +Y~dmn are the unknown constant coefficients of the expansion. Substituting now (4.1)-(4.2) into (3.13) leads to 1M/2-1 N/2-1 JxTpq = y z Jxmn Tx~ym=-/2 n=-N/2 1Am+1/2(X)ejkzpx dx J Hn(y)e-jk~yydy (4.5) = ysinc2 ~ sinc e-3pM M/2-1 N/2-1 Z Z jxmnei(2irpm/M+21rqn/N) (4.6) m=-M/2 n=-N/2 where we note that Ax/Tv = 1/M and Ay/T. = 1/N. Similarly,.1% 1 i. (pIr.iC r e-qI M/2-1 N12-1 e Z Z j~(21rpm/M+27rqn/N) (4.7) m=-M/2 n=-N/2 21

The left hand side of (3.22) then becomes E E G(kxp, kq) - jpq,pq(X Y) p=-oo q=-oo -MN E E sinc (PM) sinc () * { [ xx(kp kq) + kq)] sinc ) e + [xyGk kyq) + yyGyv(k' kq)] sinc (N) e-jw/N } rM/2-1 N/2-1 z j e-i-(2'rpm/M+2'qn/N) ej(kTpX+kyqy) m=-M/2 n=-N/2 (4.8) 4.3 Application of Galerkin's Technique From (3.22) the equations for determining the constants are M J p [-)sine2 (M) p ----oo q=-o ---- sine (q ) eN~)ej~/ +xyGx(kx'p, k1 )sinc (^) sinc2 (N) e —/l] M/2-1 N/2-1 E M jmne- (2rpm/M+2rqn/N) e(kxpx+kqy) Lm=-M/2 n=-N/2 (4.9) for x component and 1 o 00. MN E E y (kpkq)sinc2 p=-oo q=-oo +yG(k, kq)sinc (P)sinc2 (- ) M/2-1 N/2-1 Z / -ee-ji(2xrpm/M+21rqn/N) Lm=-M/2 n=-N/2 (pr sine ) ) ej(kkxpX+kyqy) (4.10) 22

for y component. To solve for Jm,(4.9)-(4.10) must now be converted to a system of equations by employing Galerkin's technique. Accordingly, the weighting functions are chosen the same as the current density expansion functions. In so doing, we multiply (4.9) by A,+112(x)llt(y) and (4.10) by L18(x)At+112 (y) (a -M/ 2,...,IM/2 - 1; t = -N/2,..., N/2 - 1) and by integrating over the patch, we obtain MN ~ [ X~kP kq)sinc4 sn2 M p=-oo q=o FM/2-1 N/2-1 I E I _mnei(21rpm/M+21rqn/N)] 1p/+2rtN [m=-M/2 n=-N/2 - x (4.11) and MN Z. yx~x(Ykx'p7 klq)sinc ~ 7sinc7) p=-oo q=-oo FM/2 -1 N/2 -11 I Z Zp M21rn/ i(21rp8/M+2,7rt/N) LM=-M/2 n=-N/2J -By. (4.12) Equations (4.1 1) and (4.12) now form a system of equations for the unknown coefficients ljmn., 23

4.4 CG-FFT Technique It can be shown that (4.11)-(4.12) can be written as M/2-1 N/2-1 [ M/2-1 N/2-1 1 A(k, kq) I Z -j(2rp'm/M+21r'n/N)| p'=-M/2 q'=-N/2 L=-M/2 n=-N/2 J *e (2-rp'slM+2-rq'tlN) = (4.13) where 00 00 A(kp, k) = f j G,(kp,, k q,,)sinc4 " ) sinc2 U=-00 V=-00 -— 700 V-007 u=-00oo V=-00oo *ej(p"r/M-q"ir/N) 00 00 Ay(k= (k,)p,, kyq.)sinc ) sinc4 (j ) u=-00oo v=-0oo0 ~e-j(p",r/M-q"./N) Ayy ^)I k I G( I'q/"snc sinc I with p" = p' + uM, q" = q' + vN. Here, we note that for large (u, v) the the summations in (4.14)-(4.17) are of the order (4.14) (4.15) (4.16) (4.17) terms of 1 (4.18) (UV)2(u2 + v2)1/2 resulting in a rapid series convergence. In accordance with the definition of the discrete Fourier transform and its inverse [9], (4.13) can be symbolically written as FFT-1 {A(k',, kq) [FFT {Jm}]} = B (4.19) 24

and it is then clear that (4.19) is suited for a conjugate gradient iterative solution eliminating a need for an explicit generation of a square impedance matrix. The CG-FFT algorithm for solving (4.19) is described as follows: Initialize the residual and search vectors: RI= B - FFT-1 {A(k ~II k' qI>) [FT{n} (4.20) Xi= FFT' 1[ (ki 1, k'qI)] [FFT {}]}(4.21) 1%= (i0,i1 (4.22) =i P3X1. (4.23) Iterate for k = 12 M Yk =FFT' q kp kjql)' [FFT {Pk}} (4.24) 1 ( 4.25) Uk=(VFk Vk j(k~l) - kI (4.26) Rk+1 = k - Cikyk (4.27) Xk+1 =FF'T1 { A (kxcpt,' 1k qI)].FFT {Rk+1}l (4.28) 1 k(4.29) (%k+, Xk+1) Pk+1 =Pk + /3kXk+l1. (4.30) Terminate at k =MN or when I IRk+l1112 < tolerance. (4.31) 11B112 Note that the quantity [A (k' I k'q) is the adjoint expression for A(k' /, k1,i) 25

Chapter 5 Reflection, Transmission and Diffraction of a Periodic Array The usual parameters characterizing an infinite periodic patch array are the reflection and transmission coefficients. In high frequency applications, high order Bragg diffractions could also be of interest. In this chapter, we describe the procedure for computing these reflected, transmitted and diffracted fields and present some numerical results. 5.1 Reflection and Transmission Coefficients Once the currents jmn are found from a solution of (4.19), the reflection and transmission coefficients can be evaluated in a straightforward manner. To show this, we refer to (3.19) which is repeated below 00 00 (,yz) = E E p(k'p,k,) -. (5.1) p=-oo q=-oo It is seen that only the Floquet mode corresponding to (p, q) = (0,0) contributes the scattered field in the direction of reflection. Thus, the transverse component 26

of this scattered field is T(y, z) = G(k, k). o[+k-k(z-d)] (5.2) where Joo = MN a ZJmn. (5.3) m,in The same Floquet mode (0,0) contributes the scattered field in the direction of transmission and this can again be found in a similar way by employing the dyadic Green's function applicable in the region z < 0. The transverse component of this scattered field is found as ET(x, y, z) = G (k,, ky). j ooe (k:++kZ) (5.4) where Joo is again given by (5.3). The total reflected and transmitted fields of the array are then obtained by adding the scattered fields to the reflected field, (3.3) or (3.7), and transmitted field, (3.5) or (3.9), associated with the dielectric slab in isolation. 5.2 Bragg Diffractions Under certain circumstance, it is possible to obtain higher order Bragg diffractions than the zeroth order (reflection and transmission) considered above. Such diffractions occur when kpq > O0 or k + k <ko (5.5) 27

which may be satisfied by certain nonzero (p, q). The direction (Opq, qpq) of a possible (p, q) diffraction is determined from the relations sin Opq cos qpq -= sin 6i cos qI + T (5.6) sin Opq sin qbpq = sin sin + (5.7) By in which A is the free-space wavelength, and the magnitude of the corresponding diffracted field in the half space z > d can be obtained from its transverse component given by ET(p, q) = G(k'P, kyq). jpq (5.8) From (5.6)-(5.7) it is seen that the diffractions appear in pairs and are symmetric with respect to the x-y plane. The diffracted field in the half space z < 0 can be found in a similar way. 5.3 Numerical Results In the following, we present some numerical results for the reflection by and transmission through an infinite periodic array associated with three types of element patch geometries illustrated in Figure 5.1. Figure 5.2 shows the reflection and transmission coefficients of a square patch array on a dielectric slab having d = 0, 0.1 and 0.2 cm for normal incidence. For the free-standing (d = 0) case the result agrees with that obtained by Cwik and Mittra [8]. The general effect of the dielectric layer is seen to shift the resonant frequency of the array. Also, note that at frequencies greater than 15 GHz higher order Bragg 28

T Lx — Tx ------ (a) U] ' T I II HF -Tx (b) Tx() (c) Figure 5.1: Three types of cell geometries. (a) A square patch cell. (b) A circular patch cell. (c) A cross-shaped patch cell. 29

1.0 0.8 0.6 0.4 0.2 0.0 1.0 0.8 0.6 0.4 0.2 0.0 I 4. 6. 8. 10. 12. 14. 16. 18. 20. 22. 24. FREQUENCY (GHz) (a) 4. 6. 8. 10. 12. 14. 16. 18. 20. 22. 24. FREQUENCY (GHz) (b) Figure 5.2: (a) Reflection and (b) transmission coefficients vs frequency for an array of square patches with various aieiectric slab thicknesses; T_, = Ty= 2.0 cm, Lx = LY=1.O CM),Er= 3.5,9t'= OsqY=0; ( ~) d =0.O cm; ( ---- )d =0.l cm; (- — ) d =0.2 cm. 30

diffractions occur and as a result the reflection and transmission coefficients are no longer complementary. Figures 5.3-5.4 show the reflection and transmission coefficients of the same patch array for three different angles of incidence. As the angle of incidence changes the frequency where the higher order Bragg diffraction occurs also changes. Figure 5.5 shows the reflection and transmission coefficients of a circular patch array as a function of the angle of incidence. The circular patch is approximately modelled as a collection of square cells. It is seen that higher order Bragg diffractions occur when Oi > 26~. The effect of loss in the dielectric slab is shown in Figure 5.6 where the reflected, transmitted and dissipated power are plotted versus frequency for a cross-shaped patch array. Results are shown for various lossy slabs and the dissipated power is computed in accordance with the relation Dissipated Power = Incident Power - Reflected Power -Transmitted Power (5.9) which is valid up to the frequency where the first order Bragg mode beyond the zeroth order appears. 31

1.0 0.8 I 0.6 0.4 0.2 0.0 1.0 0.8 0.6 0.4 I 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. FRE'QUTENCY (GHz) (a) I 0.2 0.0 4. S. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. FREQUENCY (GHz) (b) Figure 5.3: (a) Reflection and (b) transmission coefficients vs frequency for an array of square patches for three incident angles for TM incidence; T~ = = 2. 0 cm, LZ=I4Y= 1.0 cm,e,.=3.5, d =0.2 cm,k%0; ( ~) O' =Odegree;(-) 9' = 30 degrees; (- - ' =60 degrees.32

1.0I 0.8 I 0.6 0.4 0.2 0.0 1.0 0.8 0.6 0.4 0.2 0.0 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. FREQENCY(GHz) (a) I 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. FREQUENCY (GHz) (b) Figure 5.4: (a) Reflection and (b) transmission coefficients vs frequency for an array of square patches for three incident angles for TE incidence; T. = TV= 2.0 cmL = L2= 1.0 CM, f.= 3.5 d = 0.2cm qiO( ) O"=O0degree; ( ----) 9' = 30 degrees; (- - -) 0' = 60 degrees.33

1.0 0.8 H 0.6 0.4 0.2 0.0 1.0 0.8 0.6 0.4 0.2 0.0 I 0. 10. 20. 30. 40.:50. 60. 70. 80. 90. Oi (degrees) (a) 0. 10. 20. 30. 40. 50. 60. 70. 50. 90. O0 (degree) (b) Figure 5.5: (a) Reflection and (b) transmission coefficients vs angle of incidence for an array of circular patches; T_, = Tv= 2.0 cm, R = 0.625 cm, e,. = 3.5, d = 0.2 cm, O6' = 0, f = 10.4 GHz; ( ~) TM incidence;(- -)TE incidence. 34

1.0 0.8 0.6 0.4 0.2 0.0 1.0 0.8 0.6 0.4 0.2 0.0 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. ]FREQUENCY (GHz) (a) 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. FREQUENCY (GHz) (b) 35

1.0 0.8 0.6 0.4 15 0.2.isu Ca. 0.0 t --- 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. FREQUENCY (GHz) (c) Figure 5.6: (a) Reflected, (b) transmitted and (c) dissipated power vs frequency for an array of cross-shaped patches with various lossy dielectric slabs; Tr = =I 2.0 cm L =1.5 cmW=0.5 cm. d= 0.1 cm, O'= 0; ( ~) e,=4.O; ( — =4.0.jO.1;(- ).=4.0 - jO.5. 36

Chapter 6 Scattering by a Truncated Array In practical applications, an infinite periodic patch array is usually truncated but an exact solution in this case is computationally intensive and sometimes impractical because of high memory demand. In this chapter, we, instead, present an approximate solution achieved by integrating the currents J. obtained in Chapter 4 over the truncated region. Using this procedure, the scattering pattern for the finite array can be written as a product of the element pattern corresponding to a single patch and the pertinent array factor. 6.1 Bidirectionally Truncated Array To find the far zone radiated field by a single patch, one can proceed by evaluating the Fourier integrals in (2.53) using the stationary-phase formula [10]. Another much easier approach is to employ the reciprocity theorem [11]. Here we discuss both approaches. 37

Consider the integral I(1) = f, fX ) (x,y)e (zYdxdy, f > 0 (6.1) where q(x, y) is a real function and it has a simple, real stationary point (X., YB) defined by Oq O9q =0and T- = O at (x,8.(6.2) Provided f is regular near (x., y.), then an asymptotic approximation -of I(SI) for large fl is given by )eVjflq(xs,.1i) (2X 20\4a (6.3) 0JE Idet(92q/Ox.O~yj)1/2 where a = sgnd1 + sgnd2 in which (di, d2) are the eigenvalues of the matrix comprising the elements 82 q/dx8Oy,. The above stationary-phase formula can be applied for the evaluation of (2.53) which can be written as kL ii POkO e9(k2d-kx-k&PY') ~.ir(kz sin Gcos-o+ky sinesinq-k con e)dkdk, (6.4) where (r, 0, 0b) are the usual spherical coordinates describing the observation point. Comparing (6.4) to (6.1), the following identifications are possible: f = -G(k~, kY 4r.2 q =x ksin0 cos 0+ k,,sin 0sink k2 COS0. 38

The stationary point is found in accordance with (6.2) and 'is given by kxs=-ko sin 0cosqS ky=-ko sin 0sinq (6.5) and together with )11/1 jdet(&2'q/d~k.T8OkY8 1/ ko cosO' 0 ~=2 (6.6) we obtain 'jkocos9 eikor (kX8,I kyseko3 Gdsn ~(6.7).for large r. The transverse component of the far field radiated by a single patch with current J8(x, y) is then given by eT-jkO cos90io=k o iT2irr eor0 (kx3 Ikk.) ejkdc Jj- J8( YI)e63ko (sin6 cos kx'+sin6 0f s'Oy')'dx'dy'. (6.8) Note that f is not regular at 9 = ir/2 and therefore the above asymptotic expression is not valid there. We now consider the second approach. In accordance with reciprocity theorem, if El is the field due to J1 and E2 is the field due to J2, then J EJJ - J2dvJJJ7 - J7dV (6.9) holds. Assume now that a current element II 'is situated at point (r, 0, q$) and oriented in the 0 direction. This current element, if away from the dielectric surface, produces a spherical wave of the form = ek-ol..k(J.fI) (6.10) -J4rr 39

where f' is a point near the origin. We observe that (6.10) can be approximated as a plane wave with TM polarization and, therefore, the plane wave reflection coefficient (3.8) can be employed to find the reflected field near the origin. The total field at the dielectric surface near the origin is then given by = -i x kZle-iko(l -R 0 (6.11) Let now i~ denote the electric field at (r, 0, q!) produced by the patch 'current.1,. From the reciprocity theorem (6.9), we then have ~. -jkZozil -ikor( - TM)d~kodcose 47rr * j ko(8inOcos~x'+sinesinWY)~. 7,(x' y)dx'dy' 6.2 which yields ee = &3k07~o(l -RTM),j/kodccw J ej ko(sinecos.#:'+sine8inq~s')~..(' y)dx'dy'. (.3 Similarly, by orienting the element current in the direction, we find that =o ~ ~~eik (1 + RTE)eiIjd~ J jko (sin e csx'+sin 0sin4~' 7s(x', y') dx'dy'. (6.14) It can be shown that the field given by (6.13)-(6.14) is the same as that in (6.8). They all exhibit a feature of geometrical optics and are, thus, accurate provided the observation point is not near the dielectric surface where the surface wave 40

contribution may become appreciable. Therefore, in the region 0 < 9 < ir/2 we have =.koZo e-jkor -ee 4r r(1 - R TM )e jko dC Gcos COSq%[CS + coO sin qOI.) (6.15) k0Z0 e-jkortEjo o =~ -J- 4 r(1 + RTE)e kd~O~sin qOI. + cos q$I,) (6.16) where the components Ix~ are integrals given by 1X'y = JjJxye iko (x sin 0cos O+y sin 0sin -k)dxdy. (6.17) However, after employing the expansions (4.1)-(4.2) in (3.11), Ix,y can be written as summations 2 Zmixn-+Jr).e X (6.18) =YE(ymn1 + jyrnn) 2 mn ei((ko sinG COS 4k+k'g)XZmn+(k~o sinG sin b+k:V)YMn] 6.9 in which (Xmni Ymn) denotes the center point of the (mn)th cell. The above discussion applies to the radiated field in the region 0 < 0 < ir/2. The radiated field in the region ir/2 < 0 < r is obtained in a similar manner. In particular, we find that ee,4, are again given by (6.15)-(6.16) provided (1 - R TM) is replaced by TTM given in (3.10) and (1 + RTE) by TTE given in (3.6). The array factor can be easily found in the traditional manner. For an arbitrary array, it is given by AF(0, q I)e=Z [(iko sin~cos4O+kl,)xi+(ko sinG sin k+kl,)yi] (6.20) 41

where Np is the number of patches and (xi, y,) is the center point of the ith patch. In particular, for a rectangular array with N. patches in the x direction and N. patches in the y direction it is sin[(N~T./2)(ko sin 9 cos q0 + k',)] lAFO~~) -sin[(T.,/2)(ko sin 0cosq + k.,,)] sin[(N,,Ty/2)(kO sinO0 sin 0 + kg )](61 sin[(Tj~,/2)(ko sin80 sin 0 + kil,)) The scattering pattern of the truncated periodic array is now given by IE9I = leel. IAF(9,kO)J (6.22) IE,0I = IeI* I AF(0, Ok)I (6.23) and the corresponding radar cross section (RCS) can be calculated as oa(9, q!) = lim 4irr 2 E912 + jE'0j2 (6.24) r-+00 J1E012-dX 6.2 Unidirectionally Truncated Array For the special case where the array is truncated only in one direction, say x direction, we can again use the above two approaches to find the element pattern produced by a row of patches along the y~ axis. The resulting expressions are jk0 jo~ AMkds ee =-Zol -e kP(RMekdcecos OI/T(.5 eV=.-0~ e ''(+RTE)e3kodco I y / Ty (6.26) for 0 <90 < 412 and =e -z0 PTos9I/T (6.27) 42

ey = -zo e-jkOPTTEejkodosIV /Ty (6.28) for 7r/2 < 0 < ir, where p = V-T2+z2 and I,y is as given by (6.18)-(6.19) with qi, q+ = 0 or 7r. The array pattern is then computed by multiplying (6.25)-(6.26) or (6.27)-(6.28) with the array factor AF( ) sin[(NxTx/2)(ko sin 0 cos 5 + k')] AF(e )1 = sin[(T/2)(kosincos+ k)] (6.29) and the RCS is given by (6.24) with 47rr2 replaced by 2rxp. That is a(O, ) = lim 2irp IEel + EIV2 (6.30) 6.3 Numerical Results Before presenting results relating to the truncated patch array, it is first necessary to examine the accuracy of the approximate solution described above by comparison with the exact solution. While this exact solution is not available for a bidirectionally truncated array, some results do exist for arrays truncated in one direction only [12], [13]. The available data are for free-standing arrays, but this is not as relevant to the comparison. Figure 6.1 shows the approximate scattering pattern corresponding to an array having nine square patches in the x direction and being infinite periodic in the y direction. The incident field is a TM polarized plane wave with Oi = 45~ and Oi = 0. Referring to [12], the approximate results in Figure 6.1 are in good agreement with the exact given in Figures 5-9 of [12]. Another comparison is given 43

0.4 0.3 F6 0.2 0.1 0.0 4.0 3.0 9; 2.0 1.0 0.0 4.0 3.0 9- 2.0 1.0 0.0 0. 10. 20. 30. 40. 50. 60. 70. 80. 90. 0 (degrees) (a) I 0. 10. 20. 30. 40. s0. 60. 70. 80. 90. o (degrees) (b) 0. 10. 20. 30. 40. 50. 60. 70. 80. 90. o (degrees) (C) 44

6.0 4.5 3.0 -1.5 0.0I V 0. 10. 20. 30. 40. 50. 60. 70. 80. 90. o (dlegree) (d) 6.0 4.5 -3.0 - 1.5 0.0 0. 10. 20. 30. 40. 50. 60. 70. 80. 90. o (degrees) (e) Figure 6.1: Far field patterns in the4 = 1800 half-plane for a unidirectionally truncated array having niepatches in the x direction; T., = Ty= 1.0 cm, L.,= = 0.8 cm, 0' = 450, q'=0, TM incidence. (a) f = 3 GHz. (b) f = 15 GHz. (c) f = 24 GHz. (d) f = 36 GHz. (e) f = 45 GHz. 45

20. 10.. 60. 30.. 30. 60. 90. -20. -30. -40. 90. 60. 30. 0. 30. 60. 90. = 0 deg. 0 (degrees) ~ = 180 deg. Figure 6.2: Bistatic radar cross section for a unidirectionally truncated array having seven patches in the x direction; T. = Ty = 1.78 cm, Lc = 1.27 cm, Ly = 0.127 cm, Oi = 30~, qi = 180~, f = 10.8 GHz, TM incidence. The circles represent the exact data given in [13]. 46

in Figure 6.2 where the approximate scate scattering cross section, again for the TM incidence, is plotted and compared with the exact data [13] for an array having seven rectangular patches in the x direction. Note that the approximate solution agrees with the exact to within 1.0 dB in the specular direction and a reasonable agreement is observed overall. Clearly, the data in Figures 6.1 and 6.2 show that for a truncated array the termination effect can be reasonably modelled by the approximate expressions (6.22)-(6.23) provided, of course, the array size is large. Figures 6.3-6.5 present the bistatic scattering patterns in the x-z plane for a truncated array having nine cross-shaped patches in both x and y directions. The excitation is a plane wave of either TM or TE polarization with O6 = 30~ and = 180~. The results are shown for three different frequencies f = 5, 9 and 20 GHz. In view of the remarks in the previous paragraph, these are expected to be of reasonable accuracy although no comparative check has been made against the exact solutions. 47

30. 20.10. 0. -10. -20. 0. 60. 120. 180. 120. 60. 0. *=Odeg. O (degrees) *=8O deg. (a) 30. 20. 10. 0. -10. -20. 0. 60. 120. 180. 120. 60. 0. *=Odeg. O (degrees) *=8O deg. (b) Figure 6.3: Bistatic radar cross section patterns for a 9 x 9 cross-shaped patch array at f = 5 GHz; T.,, = T,= 2.0 cm, L = 1.5 cm, W = 0.5 cm, f, = 4.0, d = 0.1 cm, 0' = 300, qS = 1800. (a) TM incidence; (b) TE incidence. 48

40. 30. 20. 10. 0. -10. 0. 60. 120. 180. 120. 60. 0. * =O0deg. O (degrees) * =l18Odeg. (a) 40. 30. 20. 10. 0. -10. 0. 60. 12. 18. 120. 60. * =O0deg. O (degrees) * =l18Odeg. (b) 0. Figure 6.4: Bistatic radar cross section patterns for a 9 x 9 cross-shaped patch array at f = 9 GHz; T,, = T,= 2.0 cm, L = 1.5 cm, W = 0.5 cm, c,~ = 4.0, d = 0.1 cm, 9' = 300, qO' = 1800. (a) TM incidence. (b) TE incidence. 49

50. 40. 30. 20. 10. 0. 0. 60. 120. 180. 120. 60. 0. * =O0deg. O (degrees) = 18O deg. (a) 50. 40. 30. 20. I= s N Ac 0 10. 0. 0. 60. 120. 180. 120. 60. 0. *=Odeg. O (degrees) = 18O deg. (b) Figure 6.5: Bistatic. radar cross section patterns for a 9 x 9 cross-shaped patch array at f = 20 GHz; T,, = Ts,= 2.0 cm, L= 1.5 cm, W =0.5 cm, f,= 4.0, d = 0.1 cm, O'= 300, qS' = 1800. (a) TM incidence. (b) TE incidence. 50

Chapter 7 Truncation Effect of Dielectric Slab In Chapter 6, we dealt with a truncated periodic patch array residing on an infinite dielectric slab. If the dielectric slab is also of finite size, the diffraction caused by the edges of the slab will alter the overall scattering pattern. In this chapter, we present an approach to evaluate the diffraction contributed by the truncation of the slab. In particular, a volume integral physical optics approximation is developed to compute the diffracted field by the slab which must be added to the radiated field by the patch currents in determining the total scattering of the configuration. 7.1 Physical Optics Approximation Consider a finite planar dielectric slab of thickness d occupying the region 0 > z > d (a geometry similar to this is shown in Figure 2.1). The far zone 51

scattered field by this slab can be written as a volume integral E = j e-jkrrX r x J7(x, y,z)e(kx+ky+kzz)dxdydz (7.1) where V denotes the volume occupied by the slab, k, = ko sin 0 cos o, ky = k0o sin 0 sin q and k, = ko cos 0. In addition, J is the equivalent polarized current given by 7(x, y, z) = jWeo(e - 1)E(x, y, z) (7.2) in which E is the total electric field in the slab. In accordance with a physical optics approximation, the total field E is approximated by the field E induced in a corresponding infinite slab of the same thickness and property. For the infinite slab, E can be found analytically and for the TE case with the excitation (3.1) it is given by ' A = Asin(koz /,-sin2 + B cos ko sin2')] ~(-x sin Oi + y cos i)e(k+ ) (7.3) where A and B are constants given by j cos 0i A EoTTEe d a, A/ =- sin2 0' B = EoTTEej^kdc~ '. For the TM case with the excitation (3.2), T is given by ~ =:{ - sjin [Acos (Akoz/ - sin20' - Bsin (koz/e.- sin20 )] fr 52

(x coS O, + y sin <')- -- A sin (koz~r-sin2 0i) + B Cos (koz er-sin2 ')] Z} e(if+4) (7.4) where A and B are again the constants A- je cos 0 TMekodcos /er - sin2 0i B = EoTTMe' dc0'. Let us denote E' as ~= I(z)(k+k* (7.5) for both TE and TM cases. When (7.5) now is substituted into (7.1)-(7.2) we obtain = = k(- 1)ejkor, x x ] (z)e^Jzdz eJ i(k'+k)z+J(k'+k1)dxdy (7.6) 47rr o r Jx where S is the area occupied by the slab in the x-y plane. Introducing the definitions = d (7.7) F |= Jj )dxdy (7.8) (7.6) may be written as 4E= _ x (- )Fe-Jxrr x I. (7.9) 41rr It can be shown that r x x = (-cos 0 cos I - cos 0 sin qIy + sin lIz)0 +(sin qI, - cos bIy)$ (7.10) 53

and when this is substituted into (7.9) we finally have E9 = k(- -) lFe-Jkor(cos 0 cos qI, + cos O sin Iy - sin OI) 4irr E- = e 1) Fe-kor (- sin qI + cos qI,). 4irr (7.11) (7.12) While F depends on the specific geometry, I can be easily integrated out. In particular, for the TE case we find I- = -(As + Be) sin q' (7.13) Iy = (As + Bc) cos Oi (7.14) I= 0. (7.15) Similarly, for the TM case Is = (ac - Bs) cos q' er j/~T -sinn2O i V= - - - (ac - Bs) sin qY (7.16) (7.17) sin i Iz, = - (As + Bc). Er (7.18) In the above, s and c are given by s.= sin (koze, - sin2 0) eJ.kocOdz.o er(-am''+cosOd c = /ods(o z/.' sio')e ikcc*edz { 1 _e 'ko (fe"-"ain2 Oi'+cc* e d jk, (V/e-C 2Oi-co. ) d - 7 sin2 i -cos 0) -( r sm2 i _- cos ) (7.19) (7.20) 2ko (/r - sin2 0e + cos ) 2koEO 54

For a rectangular a x b slab, we have pa/2, b/2 II e3~(kx+k' )x+j(kyI+k'Ydd Ja/2J-b/2 - ab sinc [(ks, + k'-sifc [(k. + k.) (7.21) We note that while the above scattered field was formulated as that generated by the equivalent polarized currents, some choose to relate it to the concept of edge diffraction which may provide a more understandable interpretation. To do so, we note that F~xk k)xl ej~kyk~ki )k+ky ~. v ~[(ks +~ - (~ + ~)~.j (kx +14)x+j(k + k')]j(.2 and thus (7.6) can be rewritten as 4i7rr A r(A~~(c tX[r~ - Y j (kx+k'~)x+j(ky+klYt (.3 upon the use of Stoke's theorem, where C is the contour enclosing the area S and t is the unit vector tangential to C. The corresponding expression for the field radiated by a line source of current W(t) is - ko__ ZojkorA A j(t)e3(kxx+k~ydt(.4 and by comparison with (7.23) we observe that the edge diffraction may be attributed to an equivalent current given. by K~t) = ko~o(E t - (kA + k') - (ky + k')5 eX~ +J K)1t[(krk+ 1)2 + (k ~)2] (7.25) 55

We also note that K(t) can be further splitted into the equivalent electric and magnetic currents flowing along the edge [14, ch. 2]. 7.2 Numerical Results In the following, we present the approximate solution for the scattering by a 9 x 9 cross-shaped patch array on a 26 cm x 26 cm dielectric slab with all other parameters the same as those given in Figures 6.4 and 6.5. Figures 7.1 and 7.2 show the physical optics solutions for the bistatic scattering patterns of the slab at frequencies 9 GHz and 20 GHz, respectively. When these slab scattered fields are added to the scattered fields by the patch array (Figures 6.4-6.5), we obtain the total scattered fields of the truncated array on a finite slab. The results are shown in Figures 7.3 and 7.4. Comparing these with the corresponding data in Figures 6.4-6.5 and Figures 7.1-7.2, we observe observe that the fields radiated by the patch currents are dominant at 9 GHz while at 20 GHz the fields diffracted by the slab are dominant. 56

50. 40. 30. 20. 10. 0. -10. 0. 60. 120. 180. 120. 60. 0. * =O0deg. O (degrees) * =l18Odeg. (a) 50. 40. 30. 20. 10. 0. -l0. 0. 60. 120. 180. 120. 60. 0. * =O0deg. O (degrees) * =l18Odeg. (b) Figure 7.1: Bistatic radar cross section patterns for a 26 cm x 26 cm dielectric slab at 9 GHz; e,. = 4.0, d = 0.1 cm, 0' = 300, O~' = 1800. (a) TM incidence. (b) TE incidence. 57

60. 50. @ 40. 30. 20. 10. o.- - 0. 60. 120. 180. 120. 60. 0. *=Odeg 0 (degrees) 41=8Odeg. (a) 60. 50.40.~' 30. 20. 10. 0. 0. 60. 12. 180. 120. 60. 0. *=Odeg. O (degrees) *=8O deg. (b) Figure 7.2: Bistatic radar cross section patterns for a 26 cm x 26 cm dielectric slab at 20 GHz; c,. = 4.0, d = 0.1 cm, 0' = 300, qO" = 1800. (a) TM incidence. (b) TE incidence. 58

s0. 40.30.20. 10. 0.-10.IIII II 0. 60. 120. 180. 120. 60. 0. *=Odeg. O (degrees) *=8O deg. (a) 50. 40. 30. 20. 10. 0. -10. III _ _ _ 0. 60. 120. 180. 120. 60. 0. *=Odeg. O (degrees) *=8O deg. (b) Figure 7.3: Bistatic radar cross secti on patterns for a 9 x 9 cross-shaped patch array on a 26 cm x 26 cm dielectric slab at 9 GHz; T,, = Ts, = 2.0 cm, L= 1.5 cm, W = 0.5 cm,) 4.0, d = 0.1 cm, O' = 300,7~ 1800. (a) TM incidence. (b) TE incidence. 59

60. 50. 40. 30. 20. 10. 0. 60. 50. 40. 30. 20. 10 10. 0. 60. 120. * = 0 deg. 180. 120. 60. 0. o (degrees) * = 180 deg. (a) 0. 60. 120. 180. 120. 60. 0. * =O0deg. O (degrees) * =l18Odeg. (b) Figure 7.4: Bistatic radar cross section patterns for a 9 x 9 cross-shaped patch array on a 26 cm x 26 cm dielectric slab at 20 0Hz; T2., = TV= 2.0 cm, L = 1.5 cm, W = 0.5 cm, e, = 4.0, d = 0.1 cm, Os = 300, Os = 1800. (a) TM incidence. (b) TE incidence. 60

Chapter 8 Conclusion In this report, a spectral expression was derived for the dyadic Green's function associated with a dielectric slab. This was subsequently employed together with Floquet's theorem to formulate an equation for the patch currents on an infinite periodic array. A system of equations was then obtained by employing subdomain roof-top basis functions to expand the patch currents and applying Galerkin's technique. The conjugate gradient method in conjuction with the fast Fourier transform (CG-FFT) procedure for the solution of the system of equations was described in detail. This CG-FFT technique avoids the problem of large computer storage requirement and thus allows the treatment of large size arrays of practical interest. Based on the computed patch currents, an exact numerical solution was presented for the infinite periodic array scattering. An approximate solution for the truncated periodic array scattering was also given and its accuracy was examined by comparison with available exact numerical data for unidirectionally truncated arrays. It was found that the approximate solution is of reasonable accuracy in pre 61

dicting the scattering by the truncated array. Furthermore, the truncation effect of the dielectric slab was also studied by employing a physical optics approximation to compute the scattering by the slab. The FORTRAN computer program which was used to generate the numerical data is described and listed in the appendix. 62

References [1] S. W. Lee, "Scattering by dielectric-loaded screen," IEEE Trans. Antennas Propagat., vol. AP-19, pp. 656-665, Sept. 1971. [2] C. C. Chen, "Transmission of microwave through perforated flat plates," IEEE Trans. Microwave Theory Tech., vol. MTT-21, pp. 1-6, Jan. 1973. [3] J. P. Montgomery, "Scattering by an infinite periodic array of thin conductors on a dielectric sheet," IEEE Trans. Antennas Propagat., vol. AP-23, pp. 70-75, Jan. 1975. [4] T. A. Cwik and R. Mittra, "The cascade connection of planar periodic surfaces and lossy dielectric layers to form an arbitrary periodic screen," IEEE Trans. Antennas Propagat., vol. AP-35, pp. 1397-1405, Dec. 1987. [5] R. C. Hall, R. Mittra, and K. M. Mitzner, "Analysis of multilayered periodic structures using generalized scattering matrix theory," IEEE Trans. Antennas Propagat., vol. AP-36, pp. 511-517, April 1988. [6] B. J. Rubin and H. L. Bertoni, "Reflection from a periodically perforated plane using a subsectional current approximation," IEEE Trans. Antennas 63

Propagat., vol. AP-31, pp. 829-836, Nov. 1983. [7] S. Contu and R. Tascone, "Scattering from passive arrays in plane stratified regions," Electromagnetics, vol. 5, no. 4, pp. 285-306, 1985. [8] T. A. Cwik and R. Mittra, "Scattering from a periodic array of free-standing arbitrarily shaped perfectly conducting or resistive patches," IEEE Trans. Antennas Propagat., vol. AP-35, pp. 1226-1234, Nov. 1987. [9] W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, Numerical Recipes, chap.12, Cambridge Univ. Press, 1986. [10] L. B. Felsen and N. Marcuvitz, Radiation and Scattering of Waves, chap.4, Prentice-Hall, Inc., Englewood Cliffs, N.J., 1973. [11] R. F. Harrington, Time-Harmonic Electromagnetic Fields. chap.3, New York: McGraw-Hill, 1961. [12] W. L. Ko and R. Mittra, "Scattering by a truncated periodic array," IEEE Trans. Antennas Propagat., vol. AP-36, pp. 496-503, April 1988. [13] P. W. Grounds and K. J. Webb, "Numerical analysis of finite frequency selective surfaces," 1988 IEEE AP-S International Symposium Digest, vol. II, pp. 746-749, June 1988. [14] R. J. Marhefka, 0. M. Buyukdura, and W. Ebihara, "Radar cross section studies, phase III," The Ohio State University ElectroScience Laboratory Report No. 716621-1, April 1986. 64

Appendix This appendix contains a brief description and a list of Fortran program and subroutines for the solution of the problem considered in this report. 1. INPUT DATA There are three sets of input data describing the geometry of the problem, the numerical discretization and the incident field. All of them are contained and explained in the main program ARRAY. By looking through ARRAY, one should have no difficulty to understand and change these data. 2. OUTPUT DATA The current output format is shown by an example computation listed below. If other output format is desired, one can change the main program ARRAY or the subroutine PCPA. ****AIAY GEOMETRY INFORMATION UIITCM**** ARRAY PERIODICITY: TX * 2.000 TY - 2.000 RECTANGULAR PATCH: LI a 1.000 LY = 1.000 SLAB PSKEITTIVITY: EPS = 4.000 O. OOOJ SLAB THICKUESS: D ' 0.100 ****IICIDENT FIELD INFORMATION**** POLARIZATION: TRANSVERSE MAGNETIC (TM) WAVELENGTH: LAMBDA * 3.333 CM INCIDENT ANGLE: PHI * 180.0 DEGREES 65

THETA - 30.0 DEGREES FREQ(GHz)in 9.00 REFL(%)- 48.9904 TRANS(%) 51.1338 SCATTERING PATTERN FOR THE TRUNCATED ARRAY NO. NO. PHI(DEG) 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 1.80.0 180.0 180.0.180.0 180.0 180.0 OF PATCHES IN OF PATCHES IN THETA (DEG) 0.0000 30.00 60.00 120.0 150.0 180.0 0.0000 30.00 60.00 120.0 150.0 180.0 I-DIRECTION - I-DIRECTIONRCS(DB) 15.24 33.18 -9.256 -9.165 33.57 15.70 15.24 13.80 13.70 13.79 14.20 15.70 9 9 ELMT PATT -3.846 -4.995 -9.608 -9.516 -4.596 -3.380 -3.846 -5.283 -10.12 -10.02 -4.885 -3.380 SCATTERING PATTERI FOR THE TRUNCATED ARRAY ON A TRUNCATED SLAB A(CN)in 26.000 PHI(DEG) THETA(DEG) 0.0000 0.0000 0.0000 30.00 0.0000 60.00 0.0000 120.0 0.0000 150.0 0.0000 180.0 180.0 0.0000 180.0 30.00 180.0 60.00 180.0 120.0 180.0 150.0 180.0 180.0 B(CN)u 26.000 RCS (DB) EDGE DIFFRAC 15.22 2.096 24.23 32.05 -0.1905 -3.052 6.472 5.120 25.49 33.58 16.98 2.154 15.22 2.096 15.58 1.126 13.10 -0.3070 13.99 -8.452 15.06 -0.3098 16.98 2.153 3. PROGRAM DESCRIPTION ARRAY - PCPA - COEFF - This program is a driver for subroutine PCPA. Computes the currents on a periodic patch array residing on a dielectric slab. Calls subroutines COEFF, FARF and FAFSFF. Computes the reflection and transmission coefficients of the periodic array. 66

FARF - FAFSFF - PADFINC - GREENF - SPECTRA - SPECT1 - SPECT2 -FFTNPFFFT2 -TABLE - VNORM - PRODUCT - CPRODUCT REFLTRANS Computes the scattering pattern of a truncated patch array on an infinite dielectric slab. Computes the scattering pattern of a truncated patch array on a finite dielectric slab. Generates the coded geometry for a patch. Computes the incident field vector B as given by (3.24) and (3.25). Computes A., Ayy, Ay, and Ayy given in (4.14)-(4.17). Decides whether to call SPECT1 or SPECT2. Calls FFTN to perform FFT computation. Calls PFFFT2 to perform FFT computation. Ordinary FFT algorithm. Prime factor FFT algorithm. Finds prime factors for FFT. Calculates the euclidian norm of a vector. Multiplies two vectors. Multiplies the conjugate of a vector with a vector. Computes the reflected field by a slab in isolation. Computes the transmitted field through a slab in isolation. Note that though the program treats rectangular, circular and cross-shaped patches, other patch geometries can also be treated by a simple modification on subroutine PAD. Also note that if one has a pad size larger than 32x32, he or she needs to change the number 1024 = (32x32) into a proper number in PCPA and SPECT2. 67

PROGRAM ARRAY C THIS PROGRAM IS A DRIVER FOR SUBROUTINE PCPA. THE PROGRAM COMPUTES C C ELECTROMAGNETIC SCATTERING FROM A PERFECTLY CONDUCTING PATCH ARRAY C C ON A DIELECTRIC SLAB. C C FOR THE DETAILS OF THE FORMULATION, PLEASE REFER TO THE ATTACHED C C UNIVERSITY OF MICHIGAN TECHNICAL REPORT 389604-7-T. C C C C PCPA -- THIS SUBROUTINE HAS FOUR TASKS: C C 1: COMPUTES CURRENTS ON A PERIODIC CONDUCTING PATCH C C ARRAY RESIDING ON AN UNGROUNDED DIELECTRIC SLAB FOR A C C PLANE WAVE INCIDENCE. C C 2: COMPUTES THE REFLECTION AID TRANSMISSION COEFFICIENTS C C OF THE PERIODIC ARRAY. C C 3: COMPUTES THE SCATTERING PATTERN IN RCS FOR THE C C CORRESPONDING TRUNCATED PATCH ARRAY. C C 4: COMPUTES THE SCATTERING PATTERN II RCS FOR THE C C TRUNCATED PATCH ON THE TRUNCATED SLAB. C C***********************************************************************C C JIAINING JIN C C RADIATION LABORATORY C C DEPT. OF ELECTRICAL ENGINEERING & COMPUTER SCIENCE C C THE UNIVERSITY OF MICHIGAN C C ANN ARBOR, MICHIGAN 48109 C C MARCH 20, 1989 C COMPLEX EPS REAL LI,LY,L C$,15$$$*****$****$**$********* $$$s,$,$,$$**$***************$********$$ *****C C.....FIRST SET INPUT DATA DESCRIPTION CONCERNING THE GEOMETRY: C C C C IS - IITEGER TYPE OF PATCH C C IS-1 FOR RECTANGULAR PATCH C C IS=2 FOR CIRCULAR PATCH C C ISS3 FOR CROSS PATCH C C OTHER GEOMETRIES COULD BE TREATED BY C C MODIFYING SUBROUTINE PAD C C TX - REAL PERIODIC LENGTH II -RECTION C C TY - REAL PERIODIC LENGTH I Y-RECTIOi C C FOR IS-1 C C LI - REAL PATCH DIMENSION IN I-DIRECTION C C LY - REAL PATCH DIMENSION IN Y-DIRECTION C C FOR IS2 C C RP - REAL PATCH RADIUS C C FOR IS3 C C V - REAL WIDTH OF CROSS C C L - RtEAL LENGTH OF CROSS C C EPS - COPLI RELATIVE PERMITTIVITY OF SLAB C C D - RAL THICKNESS OF SLAB C C NI - INTIG NIUBER OF PATCHES IN THE I-DIRECTION C C MY - INTEsIG NUBER OF PATCHES IN THE Y-DIRECTIOI C C A - REAL TRUCATED SLAB DIMENSION IN I-DIRECTION C C B - REAL TRUNCATED SLAB DIMENSION IN Y-DIRECTION C C$$*$*$$$$$$$$$$$*$$$$$$$$$*$*$$$*$$$$$$$$$$$$$$$$$**$$$$$$*$*$***$$$**C IS-1 T1=2.0 TY=2.0 LI-i. 68

LYi1. RP-0. 625 V.0.56 LIl.5 EPS-(4.0,0.0) D-0.1 Am26. 0 B=26. 0 MXu9 MY=9 C...SECOND SET INPUT DATA DESCRIPTION CONCERNING THE DISCRETIZATION: C C NI - INTEGER NUMBER OF CELLS IN I-DIRECTION C C NY - INTEGER NUMBER OF CELLS IN Y-DIRECTION C C IFFT - INTEGER TYPE OF FFT ALGORITHM TO USE C C THERE ARE TWO SUBROUTINES FOR FFT COMPUTATION: C C ONE IS ORDINARY FFT ALGORITHM AND TEE OTHER IS C C THE PRIME FACTOR FFT ALGORITHM. C C FOR IFFT-1 USE THE ORDINARY FFT ALGORITHM C C NXIAND NY CAN BE ANY NUMBER FROM THE SET (8, 16, C C 32, 64, 128, 256, 512,...) C C FOR IFFT-2 USE THE PRIME FACTOR FFT ALGORITHM C C NI(=NY) CAN BE ANY NUMBER FROM THE SET (8, 9, 10,C C 12, 14, 16, 18, 18, 20, 21, 24, 28, 30, 35, 36, C C 40, 42, 45, 48, 56, 60, 63, 70, 72, 80, 84, 90, C C 105, 112, 120, 126, 140, 144, 168, 180, 210, 240,C C 252, 280, 315, 336, 360, 420, 560,..)C C TOL - REAL TOLERANCE (USUALLY.001 TO.01) C C ITMAX - INTEGER MAXIMU ITERATION NUMBER C IFFT-2 Nx16 lY-16 TOLO0.005 ITMA1=60l C...THIRD SET INPUT DATA DESCRIPTION CONCERNING THE INCIDENT FIELD: C C C C C C C C C C C C C WAVE - REAL IPOL -INTEGER PHINC - REAL THINC - REAL PHI - RBEALPH2 - REAL THi - REA TH2 - REA DPI - REAL DTN1- RMA FREE-SPACE WAVELENGTH POLARIZATION IPOLmI E-POLARIZATION (TE) (Ez=0) IPOL=2 H-POLARIZATION (TM) (Moo) INCIDENT PHI ANGLE IF KODEI, (DEGREES) INCIDENT THT ANGLE IF KODEil (DEGREES) INITIAL PHI ANGLE (DEGREES) FINAL PHI ANGLE (DEGREES) INITIAL THETA ANGLE (DEGREES) FINAL THETA ANGLE (DEGREES) INCRUMENT OF PHI ANGLE (DEGREES) INCioUENT OF THETA ANGLE (DEGREES) C C C C C C C C C C C C C C WAVNE3.33333333 IP0L=2 PHINCU18O. 0 THINC=30.0 PH1=O.0 PH2=270.0 THI=O.0 69

T12=180.O0 DPH=180.O0 DTH-30. C... CALL SUBROUTINE TO COMPUTE CURRENT DISTRIBUTION AID THEN THEE C... REFLECTION, TRAISNISSION, SCATTERING PATTERNS, ETC. CALL PCPA(WAVE,EPS,D,T1,TY,LI,LY,RP,W,L,IS,IFFT,NINY,PE1,PE2, &TEI,TH2,DPE,DTE,PHINC,TRINC,IPOL,TOL,ITNAI,NI,NY,A,B) STOP END 70

SUBROUTINE PCPA(VAVE,EPS,DD,T1,TY,L1,LY,RP,W,L,IS,IFFT,NI,NY,PH1, &PH2,TH1,TN2,DPH,DTH,PHINC,THINC,IPOL,TOL,ITEAI,MI,NY,A,BS) C THIS SUBROUTINE SOLVES THE CURRENT DISTRIBUTION ON A PERIODIC C CONDUCTING PATCH ARRAY RESIDING ON AN UN1GROUNDED DIELECTRIC SLAB C FOR A PLANE WAVE INCIDENCE. THE APPROACH IS BASED ON THE CG-FFT C METHOD. C ONCE THE CURRENT IS COMPUTED, THE SUBROUTINE COMPUTES THE C REFLECTION AND TRANSMISSION OF THE ARRAY, THE SCATTERING PATTERNS C FOR THE TRUNCATED ARRAY. C C C C C C C C. NH.11 u NI*NY PARAMETER (NM=1024) COMPLEX CUX(NM),CUY(NH),EXI(N),EY(NM),EPS,IT(NH),YT(NM),AXINM)., & AY(NM),AIXT(NM),AYT(NM),GIINH),GIY(NH),GYY (NM),GYI(NM), &RX(NH),RY(NM),PX(M),PY(NM) REAL L,LX,LY,iO,KII,IYI INTEGER SG(INM),SGX(NM),SGY(NM),SIGN,NF(4) PI=3.141592654 C...FIND PRIME FACTORS OF NI IF (IFFT.EQ.2) CALL TABLE(NX,MF,NF) CALL PAD(NE,SG,SGI,SGY,DI,DY,TX,TY,NI,NY,LI,LY,RP,V,L,IS) C...OUTPUT DATA DESCRIPTION: C C NE - INTEGER NUMBER OF CELLS IN THE PATCH C C SG - INTEGER VECTOR TAG ARRAY FOR THE PATCH C C SG(I)=O CELL OUTSIDE THE PATCH C C SG(I)inl CELL INSIDE THE PATCH C C SGI - INTEGER VECTOR SHIFTED TAG FOR ROOF-TOP FOR I CURRENT C C SGY - INTEGER VECTOR SHIFTED TAG FOR ROOF-TOP FOR Y CURRENT C NT=N1*NY CALL FINC(IPOL,VAVE,EPS,DD,DI,DY,PHINC,THINCNIT,EI,EYKOIII,KYI) C...OUTPUT DATA DESCRIPTION: C C El - COMPLEX VECTOR 1-COMPONENT OF INCIDENT FIELD C C EY - COMPLEX VECTOR Y-COMPONENT OF INCIDENT FIELD C CALL GREENIF(GXX,GIYGYY,GYX,NINY,III,KYI,KO,EPS,DD,TI,TY) C...OUTPUT DATA DESCRIPTION: C C GXX - COMPLEX VECTOR GREEN'S FUNCTION COMPONENT C C GIT - COMPLEX VECTOR GREEN'S FUNCTION COMPONENT C C GYY - COMPLEX VECTOR GREEN'IS FUNCTION COMPONENT C C GYE - CONVLEX VECTOR GREEN'S FUNCTION COMPONENT C SNITFLGAT(NT) DO Iu1,NT GII(I)GIIX(I)/SNT GYX(I)=GYI(I)/SNIT Gyy (I) GYY(I) / SNT GIY(I)-GIY(I)/SNT ENDDO 71

VNRNMVNORN(EI,EY,NT,SGX,SGY) RSSO.00 ITER-0 C... INITIALIZE CURRENT DISTRIBUTION DO I1l,NT cuX(I)-(.0,o10o) CUY (I) u(O.0,0. 0) ENDDO SIGN —i CALL SPECTRA(CUI,XT,NX,NY,SIGNMF,NF,IFFT) CALL SPECTRA(CUY,YT,NX,NY,SIGN,MF,NF,IFFT) C...OUTPUT DATA DESCRIPTION: C C IT - COMPLEX VECTOR FOURIER TRANSFORMED I-COMPONENT CURRENT C C YT - COMPLEX VECTOR FOURIER TRANSFORMED Y-COMPONENT CURRENT C CALL PRODUCT(AIT,GII,GXY,XT,YT,NT) CALL PRODUCT(AYT,GYY,GYX,YT,XT,NT) C...OUTPUT DATA DESCRIPTION: C C AXT - COMPLEX VECTOR PRODUCT TO TE TRANSFORMED C C AYT - COMPLEX VECTOR PRODUCT TO TE TRANSFORMED C SIGN-i CALL SPECTRA(AX,AXT,NX,NY,SION,KF,NF,IFFT) CALL SPECTRA(AY,AYT,NX,NY,SIGN,MF,NF,IFFT) C...OUTPUT DATA DESCRIPTION: C C AX - COMPLEX VECTOR PRODUCT AFTER TRANSFORMATION C C AY - COMPLEX VECTOR PRODUCT AFTER TRANSFORMATION C C. COMPUTE THE RESIDUAL FOR THE INTIAL CURRENTS DO 10 Iz1,NT IF(SGXI).EQ.1) THEN RXIM - EX(I)-AX(I) ELSE.RIMI u (0.0,0.0) END IF IF(SGY(I).EQ.1) TME RY(I - EY(I)-AY(I) ELSE RY(I) * (0.0,0.0) END IF 10 CONTINUE SIGNM-I CALL SPECTRA(RI,ITNI,NY,SIGNMF,N,IFFT) CALL SPECTR(RY,YT,NIX,NY,SIGN,MF,NF,IFFT) CALL CPRODCT(AXT,GIZ,GYX,IT,YT,NT) CALL CPRODUCT(AYT,GYY,GY,YT.XT,NT) CALL SPECTRA(AX,AXT,NX,NY,SIGN,MF,NF,IFFT) CALL SPECTRA(AY,AYT,NX,NY,SIQGNMF,N[F,IFFT) 72

B-I. /VNONN(AI,AY,NT,SGX,SGY) C... CONPUTE THE SEARCH VECTORS DO 20 I l,1T IF(SGI(I).EQ.1) THEN PX(I) - B*AI(I) ELSE P1(I) - (0.0,0.0) END IF IF(SGY(I).EQ.1) THEN PY(I) = B*AY(I) ELSE PYCI) = (0.0,0.0) END IF 20 CONTINUE C...BEGIN ITERATION TO FIND SOLUTION TO THE CURRENT 30 CONTINUE SIGN —I CALL SPECTRA(PX,XT,NXINY,SIGN,MF,NF,IFFT) CALL SPECTRA(PY,YT,NI,NY,SIGN,NF,NF,IFFT) CALL PRODUCT(AIT,GII,GXY,IT,YT,NT) CALL PRODUCT(AYT,GYY,GYI,YT,IT,NT) SIGN-I CALL SPECTRA(AX,AXT,NI,NY,SIGNHF,NF,IFFT) CALL SPECTBRA(AY,AYT,NI,NY,SIGN,NF,NF,IFFT) T - I./VNORN(AI,AY,NT,SGI,SGY) DO 70 I - 1,NT IF(SGX(I).EQ.I) THEN CUX(I) - CUI(I)+T*PX(I) ELSE CUI(I) a (0.0,0.0) RX(I) - (0.0,0.0) END IF IF(SGY(I).EQ.1) THE CUY(I) = CUY(I)+T*PY(I) RY(I) - RY(l)-TsAY(I) ELSE CUY(I) - (0.0,0.0) RY(I - (0.0,0.0) END IF 70 CONTINUE SIGN —I CALL SPECTRA(RX,XT,NI,NY,SIGN,HF,NF,IFFT) CALL SPECTRA(RY,YTr,NI,NY,SIGN,NF,NF,IFFT) CALL CPP.DUC(AXT, GXX,GYX,IXT,YT, NT) CALL CPIOWCT (AYT, GYY,GIXY, YT,XT,NIT) SIGoll CALL SPECTRA(AI,AXT,NX,NY,SIGN,NF,NF,IFFT) CALL SPECTRA(AY,AYT,NX,NY,SIGNHF,NF,IFFT) 73

B-I. /VNORN(AI,AY,NT,SGX,SGY) DO 50 I * l,NT IF(SGI(I).EQ.1) TREEN PX(I) - PX(I)+B*AX(I) ELSE PM() - (0.0,0.0) END IF IF(SGY(I).EQ.1) THEN PY(I) - PY(I)+B*AY(I) ELSE PY(I) = (0.0,0.0) END IF 50 CONTINUE ITER - ITERt+1 VNRKR-VNlORN(RIX,RY,NT, SGIXSGY) RSS - SQRT(VNRNR/VNRN) C...IF NOT TO MONITOR TEE ITERATIVE PROCEDURE, COMMENT TEE NEXT LINE C PRINT 140,ITER,RSS C...IF TEE CRITERIA IS SARISFIED, THEN IF(RSS.LE.TOL) TEEN PRINT 120 C...PRINT TEE INFORMATION ABOUT TEE ARRAY PRINT 200, TX, TY IF(IS.EQ.1) PRINT 220, LI, LY IF(IS.EQ.2) PRINT 240, RtP IF(IS.EQ.3) PRINT 260, L, W PRINT 280, EPS, DD C...PRINT TEE INCIDENT FIELD INFORMATION PRINT 300 IF(IPOL.EQ.1) PRINT 320 IF(IPOL.EQ.2) PRINT 340 PRINT 380, WAVE, PHINC, THINC C...IF THE CURRENT DATA ARE DESIRED, UNCOMNEN TEE FOLLOWING It LINES C... THE FIRST TASK IS COMPLETED HERB C PRINT 150 C DO Iml,NT C IF(CUI(I).EQ.(0.,0.)) PEI=O. C IF(CWI(I).NE.(0.,0.)) PHXuATAN2(AINAG(CUI(I)),REAL(CUI(I)))* C & I80./pi C IF(CUY(I).EQ.(0.,0.)) PHYO0. C IF(CWY(I).NE.(0.,0.)) PHYAITAN2(AINAG(CUY(I)),REAL(CUY(I)))* C & I80./pl C PRINT 90,CABS(CUI(I)),PHI,CABS(CUY(I)),PHY C 90 FORE&T(4413.4) C EIDDO C...CALL COUNT TO COMPUTE THE REFLECTION AND TRANSMISSION COEFFICIENTS C... THE SECOND TASK IS COMPLETED HERE CALL COEFF(CUI,cUY,SGI,SGY,NXINY,KII,IYI,K0,EPS,DD,TI,TY,THINCI A PHINC,IPOL) C....CALL FARF TO COMPUTE THE SCATTERING PATTERN FOR THE CORRESPONDING C... TRUNCATED PATCH ARRAY. TEE THIRD TASK IS COMPLETED HERE PRINT 160,MI,MY 74

CALL FAIF(CUX,CUY,SG,NI,NY,KXI,KYI,K0,EPS,DD,TX,TY,PHI,PH2,DPH, THI,TH2,DTH,NK,NY) & C...CALL FAFSFF TO COMPUTE THE SCATTERING PATTERN FOR THE TRUNCATED ON C... A TRJNCATED SLAB. THE FOURTH OR FINAL TASK IS COMPLETED HERE PRINT 170,A,BS CALL FAFSFF(CUX,CUY,SG,N1INY,PHINC,THINC,K0,EPS,DD,TX,TY,PH1, & PH2,DPH,THI,TH2,DTH,HI,NY,A,BS,IPOL) C...IF THE CRITERIA IS NOT SATISFIED, THEN ELSE IF(ITER.EQ.ITNAI) THEN PRINT 130 ELSE GO TO 30 END IF END IF 120 FORKAT( CONVERGENCE ACHIEVED, WOW!!I'//) 130 FORKAT( ITKAX EXCEEDED; NO CONVERGENCE.') 140 FORKAToITER'),I4,' RSS-',G14.4) 150 FORIIAT(' THE CURRENT DISTRIBUTION ON A PATCH'/ & ' MAG(J..I) PHASE(J..Z) HAG(J..Y) PHASE(J..)' 160 a 170 a a a FORMAT(// FORMAT(// 200 FORMAT(' 220 FORNAT(' 240 FORMAT(' 260 FORNAT(' 280 FORNAT(' 300 FORNAT( 320 FORMIAT(' 340 FORKAT( 360 FORMAT(' RETURN END SCATTERING PATTERN FOR THE TRUNCATED ARRAY')/ NO. OF PATCHES IN I-DIRECTION -',I4/ NO. OF PATCHES IN I-DIRECTION -',I4/ PHI(DEG) THETA(DEG) RCS(DB) ELNT PATTI) SCATTERING PATTERN FOR THE TRUNCATED ARRAY '/ ON A TRUNCATED SLAB ' A(CN)u',F8.3,' B(CN)u',F8.3/ PHI(DEG) THETA(DEG) RCS(DB) EDGE DIFFRAC') ****ARRAY GEOMETRY INFORMATION UNIT-CM****'// ARRAY PERIODICITY: TI in',F8.3,' TY in',F8.3) RECTANGULAR PATCH: LI in',F8.3,' LY m',F8.3) CIRCULAR PATCH: RADIUS 'l,F8.3) CROSS SHAPED PATCH: L u',F8.3,' V '),F8.3) SLAB PERMITTIVITY: EPS -',2F8.3,'3'/ SLAB THICKNESS: D 's,F8.3//) ****INCIDENT FIELD INFORMATION*s**'/) POLARIZATION: TRANSVERSE ELECTRIC (TEP)' POLARIZATION: TRANSVERSE MAGNETIC (TN)') WAVELENGTH: LAMBDA u',F8.3,' CM'/ INCIDENT ANGLE: PHI =',F8.1,' DEGREES'/ THNETA u',F8.1,' DEGREES'/) 75

SUBROUTIIE PAD(NE,SG,,SGSGY,,D,DY,TI,TY,II,NY,LX,LY,R,W,L,IS) C THIS SUBROUTINE GEIERATES A CODED GEOMETRY FOR A RECTANGULAR OR A C C CIRCULAR OR A CROSS PATCH. C C IS1 --- RECTANGULAR; IS-2 --- CIRCULAR; IS"3 -- CROSS. C IITEGER SG(*),SGX(*),SGY(*) REAL LX,LY,L IT"IX*IY DO KI1,IT SGI(x)0o SGY()=0 EIDDO HLIO. 5*LX HLY0. 5*LY DX-TX/IX DY"TY/IY HV=0.5*W HL=O.5*L IO0-O. 5*TX+O. 5*DI YO —O. 5*TY+O. 5*DY 1K0 IE-O C.....FOR RECTANGULAR PATCH IF(IS.EQ.1) THEN DO Jl1,IY Y=YO+FLOAT(J-1)*DY DO Il1,IX I"IO+FLOAT(I-1)*DI 1=1+1 IF(ABS ().LT.HLX) THEI IF(ABS (Y).LT. HLY) THEM SG(W)=1 NE-NE+1 ELSE SG()0O END IF ELSE SG(X)O0 END IF ENDDO ENIDDO GOTO 1 ELSE END IF C.....FOR CIRCULAR PATCH IF(IS.EQ.2) THEI DO Jl1,NY YmYO+FLOAT(J-1)*DY DO Iul,11 IXIO+FLOAT(I-1)*DI M1-+1 RSQIX*X+TY* IF(RSQ.LT.R*R) THEI SG(l)1l NIEIE+1 ELSE SOG()-0 END IF ENDDO 76

EIDDO GOTO 1 ELSE END IF C.....FOR CROSS PATCH IF(IS.EQ.3) THEI DO J-1,IY Y-YO+FLOAT(J-1)*DY DO I-l,11 I-XO+FLOAT(I-1 ) *DI Kw-+1 IF(ABS(I).LT.HV.AND.ABS(Y)E.LT.HL) THE ELSE IF(ABS ().LT.HL. AD.ABS(Y).LT. H) THEN SG(K)-1 ELSE SG(K)-0 END IF END IF NE"NE+1 EIDDO EIDDO ELSE END IF C.....SHIFT LABELS IN I AID Y DIRECTION I CONTINUE 1-0 DO J-I1,Y DO I-1,II IKK+1 IF(SG(!).EQ.1.AND.SG(+1).EQ. O) THEN SGI(K)-0 ELSE SGI(l)=SG(x) END IF END DO END DO DO In1,I DO J-1,NY-1 Kll(J-1)*IN+I 12"J*NI+I IF(SG(I1).EQ.1.AND.SG(12).EQ.O) THEN SGY(Il)0O ELSE SGY(r1)-SG(K1) END IF END DO END DO RETURN END 77

SUBROUTINE FINC(IFLAG,VAVE,EPS,D,DX,DY,PH,TH,NT,EI,EY,K0,KI,KY) C THis SUBrOUINE COMPUTES THE INCIDENT AND REFLECTED TANGENTIAL C C ELECTRIC FIELD ON THE PATCH. C C IFLAGI1 — TE INCIDENCE (Ez-0) C C IFLAG*2 — TM INCIDENCE (HzO0) C COMPLEX EI(*),EY(*),EPS REAL 10,II,KY COMPLEX CJ,CO,C1,C2,C3,C4,R,EITE,EYTE,EITMEYTM DATA TPI,RAD/.6283i85E+01,.1745329E-01/ CJQIMPLI(0. 11.) KO-TPI/WAVE SIP=SIX(RAD*PH) COP-COS(RAD*PH) SIT5SIN(RAD*TH) COT-COS(RAD*TH) SID-SIE(K0*D*COT) COD-COS(I0*D*COT) KXuIO*SIT*COP IYuIO*SIT*SIP C1-CSQRT(EPS-SIT*SIT) C2=CSIN(I0*D*C1) C3=CCOS(I0*D*Ci) CO0(CJ*CI*C2+COT*C3)/(CJ*CI*C3-COT*C2) IF(IFLAG.EQ.2) GOTO 10 R-(CJ*COT+CO*Cl)/(CJ*COT-CO*C1) C4-(1.+R)*CMPLI(COD,SID) EXTE=SIP*C4 EYTE=-COP*C4 DO Iinl,NT EI(I)EITE EY(I)-EYTE ENDDO RETURN 10 COu(CJ*C1*C2+EPS*COT*C~3) /(CJ*CI*c3-EPS*COT*C2) R-(EPS*COT-CJ*C1*CO)/(EPS*COT+CJ*C1*CO) C4u(1 -R)*CMPLI(COD,SID) EITM'-COT*COP*C4 EYTN'.-COT*SIP*C4 DO I-1,NT EX(I)EZXTM EY(I)=EYTN EIDDO RETURN END 78

SUBROUTINE GREENF(GXX,GXY,GYY,GYX,NX,NFY,KXI,XYIKO,EPS,D,TX,TY) C THIS SUBROUTINE COMPUTES GREEN'S FUNCTIONS G..XX, G..YX, GYY All)D C G..XY (=G-YX) FOR A GEOMETRY OF A DIELECTRIC SLAB HAVING THICK- * C -NESS D AND PERNITTIVITY EPS. THE OUTPUT IS THOSE FUNCTIONS MUL-* C TIPLIED BY EXPANSION AND TESTING FACTORS.* INTEGER P,Q REAL IXI,XYI,XXP,KYQ,X0112S COMPLEX GAM1,GAM2,SI,C0,CJ,K1,X2,Cl,C2,C3,C4,CS,C6,FAC4,EPS, & GXX(*),GYY(*),GXY(*),GYI(*) PI-3.i141592654 TPI-6.283185308 ZOu377. CJQImPLX (0. II.) HDXBTX/FLOAT(2*NX) HDY-TY/FLOAT(2*NY) IJ-3 C...BEGIN TO GENERATE ELEMENTS 1.0 DO 100 3=1,NY DO 100 I-1.NX GYY(1)u(0..0.) GXY(1)u(0..0.) GYX(K)-(0. 0.) IF(J.LE.NrY/2) THEN QuJ-1 ELSE Qu-Niy-+J END IF IF(I.LE.IX/2) THEN P.I1 -ELSE PM-N[X-1+I END IF DO 90 JV —IJ,XJ DO 90 IU —IJ,IJ KXPKXXI+TPI*FLOAT(P+IU*NX) /TX IYQ.IYI+TPI*FLOAT(Q+JV*NY)/TY BETA2IKXP*IXP+IYQ*KYQ KI=cSQRT(EPS*K[0*K0-BETA2) I1QIKPLX(REAL(11). -AIKAG(11)) 125=10*I0-BETA2 IF(125.GE.0.0) THEN I2=SQRT(125) ELSE 12m-I. 0sCJ7*SQRT(ABS (125)) END IF IF(ABS(K1WD.LE.80.) THEEN SIuC2IN(lg1*D) Gom0COS(110D) GANI=(EPS*12*SI-CJ*11*CO)/(EPS*12*CO+CJ*11*SI) GAN2U(K2*CO4+CJ*K1*Sl)/(12*SI-CJ*11*CO) ELSE GAM1n-CJ GAN20CJ KIXPTPI*FLOAT(P+IU*NX)/TX IYQ-TPI*FLOAT(Q+JV*N[Y)/TY '79

END IF IF(KIP.EQ.O.) THEN SIICI-1. ELSE SINICISI (XIP*HDI)/(KZP*HDI) END IF IF(IYQ.EQ. O.) THE SIICYsl. ELSE SIICY-SII(KYQ*HDY)/(IYQ*HDY) END IF FAC-SIICX*SIICI*SIJCY*SIICY FAC1FAC*SIICI*SINIC FAC2-FAC*SINCY*SIICY FAC3=FAC*SIICI*SIICY FAC4~CEIP(CJ*(KXP*HDX-KYQ*HDY)) IF(BETA2.EQ.O.O) THEN GIIXX(K)GXI()-ZO/(1.O+CJ*CSQRT(EPS)*GAN1)*FAC1 GYY(K) GYY(K)-ZO/(1. O+CJ*CSQRT(EPS)*GAI1)*FAC2 ELSE Zl-ZO/(IO*BETA2) Cl-CJ*Kl*12*KXP*KXP/(EPS*12*GA1I-CJ*1 ) C2-KYQ*IYQ*KO*KO/(CJ*i1*GAM2-X2) C3~CJ*I1*12*IYQ*KYQ/(EPS*K2*GANI-CJ*K1 ) C4-KXP*KXP*KO*KO/(CJ*Ki*GAN2-12) CSCJ*KI*12*KXP*KYQ/(EPS*K2*GAl1-CJ*Kl ) C6S-KP*KYQ*KO*KO/(CJ*KI*GAH2-12) GI (K) UGX(K) +Z1*(CI+C2) *FAC1 GYY (X) GYY(K)+Z1*(C3+C4) *FAC2 GY() GXY (K)+ZI*(C5-C6) *FAC3*FAC4 GYXI(X)GYX(K)+Z1*(CS-C6)*FAC3*COIJG(FAC4) END IF 90 CONTINUE 100 CONTINUE RETURN END 80

C THIS SUBROUTINE CALCULATES THE EUCLIDIAN NORM OF A VECTOR C REAL FUNCTION VNORN(A,B,N,SG1,5G2) COMPLEX A(N),B(N) INTEGER SG1l(N), SG2 (N) SUKH0. DO 10 I=1,N IF(SGI(I).EQ.1) THEN SUKNSUK+CABS (A(I) )**2 ELSE END IF IF(SG2(I).EQ.1) THEN SUN-SUK+CABS(B(I) )**2 ELSE END IF 10 CONTINUE VNOtM-SUN RETURN END SUBROUTINE PRODUCT(A,GX,GY,X,Y,N) COMPLEX A(*),GX(*),GY(*),I(*),Y(*) DO I-1,N A(I)-GI(I) *I(I)+GY(I) *Y(I) ENDDO RETURN END SUBROUTINE CPRODUCT(A,GX,GY,X,Y,N) COMPLEX A(*),GX(*),GY(*),X(*),Y(*) DO Iini,N A(I)-CONJG(GX(I) )*I(I)+CONJG(GY(I))-*Y(I) ENDDO RETURN END 81

C FOURIER TRANSFORM VIA FFT C C PARAMETER DESCRIPTION C Z...............DATA IN SPATIAL DOMAIN; SIZE N C C ZT..............DATA IN SPECTRAL DOMAIN; SIZE N C CM=== --- — C SUBROUTINE SPECTRA(Z,ZT,NZ,NZ,ISIGN,NF,NF,IFFT) INTEGER N(2),NF(4) COMPLEX Z(*),ZT(*) IF(IFFT.EQ.1) CALL SPECT1(Z,ZT,MZ,NZ,ISIGN) IF(IFFT.EQ.2) CALL SPECT2(Z,ZT,MZ,NZ,ISIGN,MF,NF) RETURN END SUBROUTINE SPECTi1(Z,ZT,NZ,NZ,ISIGN) C.....CALL THE ORDINARY FFT ALGORITHM INTEGER N(2) COMPLEX Z(*),ZT(*) N1(1) = MZ, N(2) = NZ NT - MZ*NZ C -- - - - - - - - - - - - - CI C IFORWARD FOURIER TRANSFORMI C I --- —---------- IF(ISIGN.EQ.-l) THEN DO 1=1,NT ZT(X)=Z(K) ENDDO CALL FFTN(ZT,N,2,-1) C -- - - - - - - - - - - - - CI C I IVERSE FOURIER TRANSFORMI C I-I --- —-------- ELSE DO J-1,NT Z(J)=ZT(J) ENDDO CALL FFTN(Z,N,2,1) END IF END SUBROUTINE SPECT2(Z,ZT,NZNZ,ISIGN,MF,NF) C.....CALL THE PRIME FACTOR FFT ALGORITHM INTEGER NF(4) COMPLEX Z(*),ZT(*) REAL X(1024),Y(1024) NT = MZ*NZ LPUMZ-2 82

c —..- - -—. - - ------------ —.. c I I C I FORWARD FOURIER TRAINSFOR I C I ------------------------ -I IF(ISIGI.EQ.-1) THEI DO K=1,NT X(K)=REAL(Z(K)) Y(K)=AINAG(Z(K)) EIDDO CALL PFFFT2(MZ,LP,MF,IF,X,Y,ISIGN) DO K=1,NT ZT(K)=CMPLX(X(K),Y(K)) ENDDO C --------------------------- C C I I C I INVERSE FOURIER TRANSFORM [ C I -—......... ---. ----. ---I ELSE DO K=, NT I(K)=REAL(ZT(K)) Y(K)=AIMAG(ZT(K)) EIDDO CALL PFFFT2(MZ,LP,MF,NF,X,Y,ISIGN) DO K=1,NT Z(K)-CPLI(X(K),Y(K)) ENDDO EID IF RETURN END 83

SUBROUTINE COEFF(CUI,CUY,SGI,SGY,N1INY,KII,KYI,KO,EPS,D,TX,TY,TH, & PH,IFLAG) C THIS SUBROUTINE COMPUTES REFLECTION AND TRANSMISSION COEFFICIENTS* C OF THE INFINITE ARRAY USING THE COMPUTED CURRENT INTEGER SGI(*),SGY(*) REAL 1II,KYI,K0,K2S COMPLEX GAM1,GAM2,SI,CO,CJ,K1,K2,Z1,C1,C2,C3,C4,CS,C6,FAC3,FkC4, &EPS,GII,GYY,GXY,GYX,CUI(*),CUY(*),EX,EY,SCUX,SCUY,EXT,EYT,EIR,EYR DATA TPI,RAD/.6283185E+O1,.1745329E-01/ COTZCOS(RAD*TH) COT2=COT*COT PI=3. 141592654 TP-6.283185308 ZO=377. CJ=CMPLX(0.,1.) HDX=TX/FLOAT(2*NX) HDY=TY/FLOAT(2*NY) C....COMPUTE G.XI, G..YY, G-XY AND G-YX FOR REFLECTED FIELD BETA2=IXI*KXI+IYI*XYI K1-CSQRT (EPSsKO*KO-BETA2) K11Q!PLI(REAL(K1),-AIMAG (11)) K2S510*KO-BETA2 IF(K2S.GE.0.0) THEN I2-SQRT(125) ELSE 12=-i.0*CJ*SQRT (ABS (12S)) END IF IF(ABS(K1*D).LE.80.) THEN SI-CSIN(Kl*D) CO=CCOS(11*D) GAM1U(EPS*K2*SI-CJ*11*CO)/ (EPS*12*CO+CJ*K1*SI) GAM2-(12*CO+CJ*K1*SI) / (2*SI-CJ*KI*CO) ELSE GAMi —CJ GAM2-CJ END IF Ii(BETA2.EQ.0.0) THEN GIXXn-ZO/(i. +CJ*CSQRT(EPS) *GAM1) GYY —ZO/ (1. +CJ*CSQRT(EPS) *GAM1) GXY(0. 10.) GYI=(0.,0.) ELSE Z1=ZO/ (K0*BETA2) C1-CJ*11*12*1XI*1XI/(EPS*K2*GAM1-CJ*K1) C2IKYI*KYI*10*KO/(CJ*11*GAN2-K2) C3.CJ*K1*12*IYI*IYI/(EPS*12*GAM1-CJ*Kl) C4111I*KXI*KO*KO/(CJ*K1*GAN2-12) CSUCJ*11*12*KXI*KYI/(EPS*12*GAM1-CJ*Kl) C~6uKII*1YI*10*K0/(CJ*K1*GAM2-X2) GIX=Zl*(CI+02) GYYuZ1*(C3+C4) GXYuZ1*(CS-C6) GYIZZI*(C5-C6) END IF C...COMPUTE THE SUMMATION OF THE CURRENT NTNI*N[Y SCU-(0.,0.) SCUY-(0.,0.) 84

DO 70 I - 1,NT IF(SGX(I).EQ.1) THEN SCUX = SCUI+CUX(I) ELSE END IF IF(SGY(I).EQ.1) THEN SCUY - SCU`Y+CUY(I) ELSE END IF 70 CONTINUE C....COMPUTE REFLECTED FIELD E1-(GXX*SCUI+GXY*SCUY)/NT EY=(GYX*SCUX+GYY*SCUY) /NT CALL REFL(IFLAG,K0,EPS,D,PH,TH,EIR,EYR) IF(IFLAG.EQ.i) THEN SEX-CABS (EX+EXR) **2 SEYUCABS (EY+EYR) **2 ELSE SEI-CABS (EX/COT+EIR)**2 SEY-CABS (EY/COT+EYR)**2 END IF R=SEI+SEY C....COMPUTE G-XX, G-YY, G.IY AND G-Y1 FOR TRANSMITTED FIELD IF(BETA2.EQ.0.0) THEN GXXICJ*ZO*CSQRT(EPS)/(CSQRT(EPS)*GAMl-CJ)/(CJ*SI+CSQRT(EPS)*CO) GYY-GII GIIUGXX GIY(0.,0.) GY=(0.,0.) ELSE Z1=ZO/(K0*BETA2) CI-CJ*EPS*Kl*12*12*KXI/(EPS*K2*GAMI-CJ*Kl)/(CJ*Ki*SI+EPS*12*CO) C2-CJ*K1*KO*KO*KYI/ (CJ*K1*GAM2-K2) / (2*SI-CJ*K1*CO) C3-CJ*EPS*KI*12*12*IYI/(EPS*12*GAMl-CJ*K1)/(CJ*K1*SI+EPS*12*CO) C4-CJ*I1*KO*R0*KII/ (CJ*K1*GAM2-K2) /(K2*SI-CJ*K1*CO) GXXIZ1*(C1*KXI-.C2*KYI) GYY=Z1*(C3*KYI-C4*KXI) GXY=Zl*(Cl*KYI+C2*KXI) GYX=Zl*(C1*KYI+C2*KXI) END IF C...COMPUTE TRANSMITTED FIELD EXI(GII*SCUI+GXY*SCUY) /NT EY=(GYX*SCUI+GYY*SCUY)/NT CALL TRANS(IFLAG,I0,EPS,D,PH,TH,EIT,EYT) IF(IFLAG.EQ.1) THEN SEXUCABS (EX+EXT) **2 SEY-CABS (EY+EYT) **2 ELSE SEI-CABS (EI/COT+EXT)s*2 SEY=CABS(EY/COT+EYT)**2 END IF T=SEI+SEY DIS-1.-R-T VRITE(6,90) K0/TPI*30.,R*100.,T*100. 90 FORMAT(//'FREQ(G~z)-',F6.2,' REFL(%)=',F8.4,' TRANSQ.)=',F8.4) RETURN END 85

SUBROUTINE REFL(IFLAG,K0,EPS,D,PE,TH,EI,EY) C THIS SUBROUTINE COMPUTES THE REFLECTED FIELD BY SUAB IN ISOLATION C C IFLAG1l TE INCIDENCE (Ez=O) C C IFLAG-2 TX INCIDENCE (HzO0) C COMPLEX EX,EY,EPS,CJ,CO,C1,C2,C3,C4,R,T REAL IO,KI,KY DATA TPI,RAD/.8283185E+O1,.1745329E-O1/ CJGCMPLI(O.,1.) SIP=SIN(RAD*PH) COP=COS(RAD*PH) SIT-SIN(RAD*TH) COT-COS(RAD*TH) SID=SIN(KO*D*COT) COD-COS(KO*D*COT) Cl=CSQRT(EPS-SIT*SIT) C2=CSIN(KO*D*Cl) C3=CCOS(KO*D*Cl) CO=(CJ*Cl*C2+COT*C3)/ (CJ*CI*C3-COT*C2) IF(IFLAG.EQ.2) GOTO 10 R-(CJ*COT+CO*Ci) /(CJ*COT-CO*Ci)*auPLI(COD, SID) E1=-SIP*R EY=COP*R RETURN 10 CO0(CJ*C1*C2+EPS*COT*C3) /(CJ*C1*C3-EPS*COT*C2) R-(EPS*COT-CJ*Cl*CO)/(EPS*COT+CJ*CI*CO)*CMPLI(COD,SID) EXu-COP*R EY —SIP*R RETURN END SUBROUTINE TRANS(IFLAG,KO,EPS,D,PH,TH,EI,EY) C THIS SUBROUTINE COMPUTES THE TRANSMITTED FIELD THROUGH THE SLAB IN C C ISOLATION C C IFLAG-1 - TE INCIDENCE (Ez0O) C C IFLAG-2 — TM INCIDENCE (Hz-0) C COMPLEX EX,EY,EPS,CJ,CO,Ci,C2,C3,C4,R,T REAL KO,KX,KY DATA TPI,RAD/.6283185E+01,.1745329E-01/ CjCcMPLI(0. I I.) SIP=SIN(RAD*PH) COP-COS(RAD*PH) SIT=SIN(RAD*TH) COT=COS(RAD*TH) SID=SINI(IO*D*COT) COD-COS(KO*D*COT) C1UCSQRT(EPS-SIT*SIT) C2-CSIN(I0*D*C1) C3inCCOS(I0*D*C1) COu(CJ*C1*C2+COT*C3)/ (CJ*C1*C3-COT*C2) IF(IFLAG.EQ.2) GOTO 10 R-(CJ*COT+CO*C1) /(CJ*COT-CO*Cl) T=(1.+R) /(CJ*COT*C2/Cl+C3) *CMPLI(COD, SID) E1=-SIP*T 86

EY-COP*T RETURN 10 CO.(CJ*CI*C2+EPS*COT*C3) /(CJ*Ci*C3-EPS*COT*C2) R-z(EPS*COT-CJ*C1*CO)/ (EPS*COT+CJ*Cl*CO) T-(1.+R)/(CJ*EPS*COT*C2/Ci+C3)*CNPLI(COD,SID) EI=COP*T EY=SIP*T RETURN END 87

SUBROUTINE FARF(CUX,CUY,SG,NI,NY,KII,XYI,KO,EPS,D,TX,TY,PH1I,PH2, k DPH,TH1,TH2,DTH,MI,MY) C THIS SUBROUTINE COMPUTES SCATTERED FIELD OF THE FINITE ARRAY * C ON THE INFINIT SLAB USING THE COMPUTED CURRENT* INTEGER SG(*) REAL KII,KYI,KO,KI,KY COMPLEX CO,CJ,C1,C2,C3,EPS,CUX(*),CUY(*),PHASE,RTE,RTM,SCUX,SCUY, &ETH,EPH,TTE,TTM DATA TPI,RAD/.6283185E+01,.1745329E-Ol/ Z0=377. WAVE=TPI/KO N[PHINT( (PH2-PH1)/DPH) NTH=INT( (TH2-TH1)/DTH) DO 1000 I10,NPH PH-PH1+FLOAT(I)*DPH DO 1000 J30,NTH TH-TH1+FLOAT(J)*DTH IF(TH.EQ.90.) GOTO 1000 C....COMPUTE THE REFLECTION COEFFICIENT SIP-SIN(RAD*PH) COP-COS(RAD*PH) SIT-SIN(RAD*TH) COT-COS(RAD*TH) SID-SIN(K0*D*COT) COD=COS(K0*D*COT) C1=CSQRT(EPS-SIT*SIT) C2-CSIN(K0*D*C1) C3-CCOS(KO*D*C1) CO (CJ*Cl*C2+COT*C3) /(CJ*C1*C3-COT*C2) RTE=(CJ*COT+CO*C1) /(CJ*COT-CO*CI) CO0=(CJ*Cl*C2+EPS*COT*C3) /(CJ*C1*C3-EPS*COT*C2) RTM=(EPS*COT-CJ*C1*CO) /(EPS*COT+CJ*Cl*CO) IXinK0*SIT*COP+KXI KY=IO*SIT*SIP+KYI IF(TH.LT.90.) GOTO 100 C...COMPUTE THE TRANSMISSION COEFFICIENT TTE-(1.+RTE)/ (CJ*COT*C2/C1+C3) TTM-(1.+RTM)/ (CJ*EPS*COT*C2/C1+C3) 100 CONTINUE C...INTEGRATE THE PATCH CURRENT SCUI-(0.,0.) SCUYU(0. 10.) DIUTI/NZ DY.TY/NY 10-0. 5*TI+0. 5*DX YO —0. 5*TY+0. 5*DY K=0 DO JJ-1,NY Y=YO+FLOAT(JJ-1) *DY DO IIin1,NI X*IO+FLOAT(II-1) *DI K=K+1 IF(SG(K).EQ.1) THEN 88

PHASECEXP(CJ* (X*Kx+Y*KY)) SCUI=SCUI+0.5* (CUX(K) +CUIX(K-i) )*PHASE SCUY=SCUY+0.5* (CUY(K)+CUY (K-NIX)) *PHASE ELSE END IF ENDDO ENDDO SCUZ=SCUI*DX*DY SCUY=SCUY*DX*DY C...ELEMENTARY PATTERN IF(TH.GT.90.) GOTO 200 ETH-CJ*KO*ZO*(i.-RTN)/(2.*TPI)*(COT*COP*SCUX+COT*SIP*SCUY) EPH-CJ*KO*ZO*(l. +RTE)/(2.*TPI)*(-SIP*SCUX+COP*SCUY) GOTO 300 200 ETH-CJ*KO*ZO*TTN/(2.*TPI)* (COT*COP*SCUI+COT*SIP*SCUY) EPH=CJ*KO*ZO*TTE/(2.*TPI)*(-SIP*SCUX+COP*SCUY) 300 CONTINUE C...ARRAY FACTOR IF(KX.EQ.0.) THEN FACX=FLOAT(MI) ELSE FACI=SIN (0. *FLOAT(NI)*TX*KXA)/SIN(0. 5*TX*KX) END IF IF(KY.EQ.0.) THEN FACY-FLOAT (NY) ELSE FACY-SIN(0. 5*FLOAT(NY)*TY*KY)/SIN(0. 5*TY*KY) END IF FAC=FACI*FACY C...SCATTERING CROSS SECTION SIGNA=2.*TPI* (CABS (ETH)**2+CABS(EPH)**2) *FAC*FAC SIGN&-iO.*LOG10(SIGNA/WAVE/WAVE) SIGNA1*2.*TPI* (CABS (ETH) **2+CABS (EPH) **2) SIGNAI1O. *LOGiO (SIGNAl/WAVE/WAVE) WRITE(6,90) PH,TH,SIGNA,SIGNA1 90 FORKAT4G13.4) 1000 CONTINUE RETURN END 89

SUBROUTINE FAFSFF(CUI,CUY,SG,NI,NY,PHINC,THINC,KO,EPS,D,TX,TY,PHi, & PH2,DPH,TH1,TH2,DTH,MI,NY,A,B,IPOL) C THIS SUBROUTINE COMPUTES SCATTERED FIELD OF THE FINITE ARRAY C ON A FINITE SLAB THROUGH A SUPERPOSITION OF THE FIELD RADIATED C BY THE PATCH CURRENTS AND THE FIELD DIFFRACTED BY THE SLAB EDGES C OBTAINED USING THE VOLUME INTEGRAL PHYSICAL OPTICS METHOD. * * * * INTEGER SG(* REAL KII,KYI,10,KX,KY,IX,IY,IT COMPLEX CO,CJ,C1,C2,C3,C4,EPS,CUIX(*),CUYC*),PHASE,RTE,RTM,SCUX, &SCUJY,ETH,EPH,TTE,TTM,J1,JY,JZ,FTHR,FPHRt,FTHD,FPHD,AA,BB, &CC,DD,SIC1,COC1,SIC2,COC2,SIC,COC,JX1,JY1,JZ1I,FAC,Cll DATA TPI,RAD/.6283185E+01,.1745329E-O1/ CJ=CMPLX(O. A.) ZO=377. WAVE=TPI/K0 C...CALCULATE THE REFLECTLECTION COEFFICIENTS WITH INCIDENT ANGLE SIP-SIN(RAD*PHINC) COP-COS(RAD*PHINC) SIT=SIN(RAD*THINC) COT-COS(RAD*THINC) SID-SIN(I0*D*COT) CODUCOS(X0*D*COT) 11I1K0*SIT*COP 1YI110*SIT*SIP C1-CSQRT(EPS-SIT*SIT) C1I-C1 C2-CSIN(K0*D*Cl) C3uCCOS(10*D*C1) COu(CJ*C1*C2+COT*C3)/ (CJ*CI*C3-COT*C2) RTE-(CJ*COT+CO*C1) /(CJ*COT-CO*C1) CO. (CJ*C1*C2+EPS*COT*C3) /(CJ*C1*C3-EPS*COT*C2) RTM-(EPS*COT-CJ*Ci*CO) /(EPS*COT+CJ*Cl*CO) C.CA..LCULATE THE EQUIVALENT POLARIZED CURRENTS C4-10*K0/(2.*TPI)*(EPS-1.) IF(IPOL.EQ.2) GOTO 500 AA(1.+RTE)/(C2-CJ*C1/COT*C3) *CMPLIX(COD, SID) BB —CJ*CI/COT*AA JXII-C4*SIP JYl1C4*C0P JZI-CMPLI(0..0.) GOTO 600 500 kAA(1.+RTN)/(C2-CJ*C1/COT/EPS*C3)*CMPLX(COD,SID) BB=-CJ*CI/COT./EPS*AAi JX1-C4*COP JYI1C4*SIP JZlI-C4 600 CONTINUE C...START TO COMPUTE THE SCATTERED FAR FIELD AT DIFFERENT ANGLES NPH-INT( (PH2-PH1)/DPH) NTH-INT((TH2-TH1)/DTH) DO 1000 I10,NPH PH-PHI+FLOAT(I)*DPH DO 1000 J-0,NTH TH.TH1+FLOAT(J)*DTH IF(TH.EQ.90.) GOTO 1000 90

C...COMPUTE THE REFLECTION COEFFICIENTS WITH OBSERVATION ANGLE SIP=SIN(RAD*PH) COP-.COS(RAD*PH) SIT-SIN(RAD*TH) COT-COS(RAD*TH) SID-SIN (X0*D*COT) COD-COS(X0*D*COT) Cl1CSQRT(EPS-SIT*SIT) C2=CSIN (K0*DsC1) C3=CCOS(K0*D*C1) CO0i(CJ*C1*C2+COT*C3)/(CJ*C1*C3-COT*C2) RTE-(CJ*COT+CO*Cl) /(CJ*COT-CO*CI) CO0(CJ*C1*C2+EPS*COT*C3) /(CJ*C1*C3-EPS*COT*C2) RTMU(EPS*COT-CJ*C1*CO) /(EPS*COT+CJ*Cl*CO) KX=K0*SIT*COP+KXI KY=K0*SIT*SIP+KYI C....COM PUT E THE CORRESPONDING TRANSMISSION COEFFICIENTS IF(TH.LT.90.) GOTO 100 TTE-(1.+RTE)/ (CJ*COT*C2/C1-+C3) TTM*(1.+RTM)/(CJ*EPS*COT*C2/C1+C3) 100 CONTINUE C...INTEGRATE THE PATCH CURRENTS SCUZ-(0. 10.) SCUY*(0..0.) DI-TZ/NX DY-TY/N`Y 10u-0. 5*T1+0. *DI YOE-. 5*TY+0. 5*DY 1.0 DO JJ-1,NY Y-YO+FLOAT(JJ-1) *DY DO II-1,NI 1=10+FLOAT(II-1) *DX I=1+1 IF(SG(K).EQ.1) THEN PHASE=CEXP(CJ* (X*KX+Y*XY)) SCUI-SCUI+0. 5*(CUI(K)+CUX(1-1) )*PHASE SCUY=SCUY+0.5*(CUY(K)+CUY(K-NI) )*PHASE ELSE END IF ENDDO ENDDO SCU'I-SCU'I*DX*DY SCUY-SCUY*DX*DY C...COMPU`TE THE ELEMENTARY PATTERN IF(TN.GT.90.) GOTO 200 ETE —CJ*1O*ZO*(1.-RTN)/(2. *TPI)*(COT*COP*SCUI+COT*SIP*SCUY) EPH=-CJ*KO*ZO*(1l.+RTE)/(2. *TPI)*(-SIP*SCUX+COP*SCUY) GOTO 300 200 ETH —CJ*K0*ZO*TTM/(2.*TPI) *(COT*COP*'SCUI+COT*SIP*SCUY) EPH —CJ*K0*ZO*TTE/(2. *TPI)*(-SIP*SCUI+COP*SC!JY) 300 CONTINUE ETH-ETH*CMPLX(C3,C2) EPE-EPH*CIPLX(C3,C2) C....CONMPUE THE ARRAY FACTOR 91

IF(IZ.EQ.0.) THEN FACI-FLOAT (Ml) ELSE FACXISIN (0.5*FLOAT(MI)*TI*X) /SIN(O. 5*TX*KX) END IF IF(KY.EQ.O.) THEN FACY=FLOAT (NY) ELSE FACY=SIN(0. 5*FLOAT (MY) *TY*KY) /SIN(0. 5*TY*KY) END IF FAC=FACI*FACY*CEIP(CJ*0. 5*FLOAT(MXI-) *TI*KI) k*cEIP(CJ*0. 5*FLOAT (MY-I) *TY*KY) C...COMPUTE THE FIELD RADIATED BY THE PATCH CURRENTS FTHR-FAC*ETH FPHR-FAC*EPH C...START TO COMPUTE THE DIFFRACTED FIELD BY THE SLAB IF(KI.EQ.O.) THEN I11A ELSE IX-nA*SIN(. 5*A*KX) /(0. 5*A*KI) END IF IF(KY.EQ.0.) THEN IY-B ELSE IY=B*SIN(0. 5*B*KY) / (0. *B*KY) END IF IT-IX*IY C1-Ci1 SICI-CSIN(K0*D*(CI+COT)) COCI1CCOS(KO*D*(C1+COT)) SIC2-CSIN(KO*D*(C1-COT)) COC2-CCOS(KO*D*(C1-COT)) SIC-(1.-QIPLX(COC1,SICI))/(2.*KO*(CI+COT))+ & (1.-CMPLX(COC2,-SIC2))/(2.*KO*(CI-COT)) COC-((1.-CNPLI(COCI,SICI))/(2.*KO*(CI+COT))& (i.-CMPLX(COC2,-SIC2).)/(2.*KO*(C1-COT)))*CJ IF(IPOL.EQ.2) GOTO 700 ~CC-AA*SIC+BB*COC JI-JXI*CC JYUJYI*CC iz=0. GOTO 800 700 ~CuC3*C1/EPS*(AA*COC-BB*SIC) DD —SIT/EPS*(AA*SIC+BB*COC) J1-J11*CC JYuJY1*CC JZ-JZI*DD 800 CONTINUE FTHD.(COT*COP*JI+COT*SIP*JY.-SIT*JZ) *IT FPED=(-SIP*JX+COP*JY)*IT C....FIND THE TOTAL SCATTERED FIELD ETH-FTUR+FTHD EPHWFPHR+FPHD C...SCATTERING CROSS SECTION SIGN=2. *TPI* (CABS (ETH)**2+CABS(EPH)**2) SIGNAIO0. *LOGIO(SIGM/ WAVE/WAVE) 92

SIGNA1=2.*TPI*(CABS(FTED)**2+CABS(FPHD)**2) SIGNRiIlO. *LOGGO(SIGNA1/WAVE/UAVE) VRITE(6,90) PH,TH,SIGKA,SIGNA1 90 FORNAT(4G13.4) 1000 CONTINUE RETURN END 93

C ********************************** C * FFT ROUTINE FOR COMPUTATION OF N DIMENSIONAL FOURIER TRANSFORM* C * TASKS: 1)BIT REVERSAL 2)TRIGONOMETRIC RECURRENCE 3)TRANSFORM * C *ALL DONEIINTHIS PROGRAM* C PARAMETER DESCRIPTION C C DATA..........A REAL ARRAY IN WHICH DATA ARE STORED AS IN C A MULTIDIMENSIONAL COMPLEX FORTRAN ARRAY C NDIM..........DIMENSION OF DATA AND THE FFT C NN...........INTEGER ARRAY OF LENGTH NDIM C ISIGN.........DIRECTION OF THE TRANSFORM: C I -FOREWARD FFT C -1 -INVERSE FFT TIMES THE PRODUCT OF C LENGTHS OF ALL DIMENSIONS C u SUBROUTINE FFTN(DATA,NN,NDIM,ISIGN) REAL*8 WR,WI,WPR,WPI,WTEMP,THETA DIMESION NN(NFDIM),DATA(*) NTOTI1 DO IDIM1l,NDIM NTOTUNTOT*NN(IDIM) ENDDO NPREVI1 DO IDIMU1,NDIM NNN1(IDIM) NREMN=TOT/(N*NPREV) IPl12*NPREV IP2=IPI*N IP3=IP2*NREM I2REV-1 C BIT REVERSAL SECTION DO I2=1,IP2,IP1 IF(I2.LT.I2REV)THEN DO Ii1I2,I2+IPi-2,2 DO I311I,IP3,IP2 I3REV-I2REV+I3-I2 TEMPR-DATA(I3) TEMPI=DATA(I3+1) DATA(I3)=DATA(I3REV) DATA(I3+1)=DATA(I3REV+1) DATA (I3REV)-TEMPR DATA(I3REV+1)-TEMPI ENDDO ENDDO END IF IBIT=IP2/2' 1 IF((IBIT.GE.IPI).AND.(I2REV.GT.IBIT))THEN I2REVuI2REV-IBIT IBITwIBIT/2 GO TO 1 END IF I2REV-I2REV+IBIT ENDDO C DANIELSON-LANCZOS FORMULA IFPIlIPI 2 IF(IFP1I.LT.IP2)THEN IFP2=2*IFP1 THETA=ISIGN*6. 28318530717959D0/ (IFP2/IPI) WPR=-2.DO*DSIN(O. 5D0*THETA)**2 94

WPIsDSII(THETA) WRil.DO VI=O.DO DO I3=1,IFP1,IP1 DO I1=I3,I3+IP1-2,2 DO I2=I1,IP3,IFP2 KI=I2 K2=K1+IFP1 TEMPR-SIGL(iR)*DATA(K2)-SIGL(VI)*DATAA(X2+1) TEMPImSIGL(WR)*DATA(K2+1)+SIGL(WI)*DATA(K2) DATA(K2)~DATA(K1)-TEMPR DATA(K2+1)-DATA(K1+1)-TEMPI DATA (K1)DATDATA(K1)+TEMPR DATA(Kl+1)~DATA(K1+1)+TEMPI ENDDO EIDDO WTEMPiWRi C TRIGONOMETRIC RECURRENCE WR=WR*WPR-WI*WPI+WR WI=WI*WPR+WTENP*WPI+WI EIDDO IFP1=IFP2 GO TO 2 END IF NPREV=I*IPREV EIDDO RETURN END 95

SUBROUTINE PFFFT2(L,LP,M,N,X,Y,SIGN) C THIS SUBROUTINE COMPUTES A TWO - DIMENSIONAL FFT USING C C A PRIME FACTOR ALGORITHM C C C C ALLOWED FFT DIMENSIONS FROM THE SET (2,3,4,5,7,8,9,16) C C C C L=2 L=14 L-36 L=80 L=180 L=560 C C L-3 L=15 L=40 Lz84 L=210 L=630 C C Lin4 L:16 L-42 L-90 L-240 L-720 C C L-5 L-18 L-45 L106 L-262 L-840 C C L-6 L-20 L-48 L-112 L-n280 L=1008 C C L=7 L-21 L-56 L-120 L=315 L=1260 C C L=8 L=24 L=60 L=126 L-336 L-1680 C C L=9 L=28 L=63 L-140 L=360 L=2520 C C L-10 L-30 L-70 L-144 L=420 L56040 C C L=12 L=35 L-72 L1i68 L=504 C C C C INPUT VARIABLE DESCRIPTIONS: C C C C L - I*4 ROW AND COLUMN DIMENSION OF FFT C C LP - I*4 ROW WIDTH AND COLUMN 'WIDTH OF NON-ZERO DATA C C M - I*4 NUMBER OF PRIME FACTORS C C N(1:4) - I*4 PRIME FACTORS (MAY BE IN ANY ORDER) C C XI(:L*L) - R*4 REAL PART OF INPUT C C Y(1:L*L) - R*4 IMAGINARY PART OF INPUT C C SIGN - I*4 SIGN-i FOR FORWARD AND -1 FOR INVERSE TRANSFORM C C C C OUTPUT VARIABLE DESCRIPTIONS: C C C C XI(:L*L) - R*4 REAL PART OF TRANSFORM 0 C Y(1:L*L) - R*4 IMAGINARY PART OF TRANSFORM C REAL X(*),Y(*) INTEGER N(s),MA(16),MB(16),KI(16),SIGN DATA C31,C32,CS1,C52 /.86602540,.50000000,.95105652,1.5388418/ DATA C53,C54,C55,C71 /.36327126,.56901699,-1.250000,-1.1666667/ DATA C72,C73,C74,C75 /.79015647,.05585427,.7343022,.44095855/ DATA C76,C77,C78,C81 /.34087293,.53396936,.87484229,.70710678/ DATA C92,C93,C94,C96 /.93969262, -.17364818,.76604444, -.34202014/ DATA C97,C98,C162,C163 /-.98480775, -.64278761,.38268343,1.306563/ DATA C164,C165 /.54119610,.92387953/ IF(SIGN.EQ. 1) THEN DO 1 I1i,L*L Y(I)u-Y(I) 1 CONTINUE ELSE EN IF C STEP THROUI ROW AND COLUMN TRANSFORMS DOING ROWS FIRST AND C THEN COLUMNS C C ICI12-L IC2-2*L-1 DO 2 IRC-1,2 ICI1ICI+L-1 IC2-IC2-L+l 96

C STEP THROUGH THE N STEPS C DO 3 I-l,M C COMPUTE THE PERMUTATION OUTPUT MAP FOR THE ROTATED DFTS C DO 4 II=1,1(I) JT-(II-1) *1 MB(II)=MOD(JT,N(I) )+l 4 CONTINUE C INITIALIZE THE FIRST GROUP OF POINTS AND CALL THE APPROPRIATE DFT C DO 5 IJ-0,N(I)-i ZI(MB(IJ+1))1IC2* (IJ*X+1)+1-IC2 5 CONTINUE IF ( ((IGN.EQ. -1) AND. (IRC EQ. 1)) &.OR. ((SIGN.EQ. 1) AND. (IRC.EQ. 2))) THEN ILIM-LP+2-1 ELSE ILIM-L-1 END IF IF (N (I).EQ. 2) GO TO 100 IF (NI).EQ. 3) GO TO 101 IF (N(I).EQ. 4) GO TO 102 IF (NI).EQ. 5) GO TO 103 IF (N(I).EQ. 7) GO TO ~04 IF (N (I).EQ. 8) GO TO 105 IF (NI).EQ. 9) GO TO 106 IF WNI).EQ. 16) GO TO 107 100 CONTINUE C 2 POINT DFT C DO 6 J31,K HMA(1)I(1) MA (2) =11(2) DO 7 IG0O,ILIN IO=IG*IC1 I1-MA(1)+IO I2MNA(2)+IO C*****s*2-POINIT DFT 1RA**~*s***s****s******s** TiIl(I1) X(MA(MB(l))+IO)=T1+lI(2) I(MA(MB(2))+IO)=T1-lI(2) TI-Y(I1) Y(MA(NB(1) )+IO)-T1+Y(I2) Y(MA(XB(2))+IO)inTl-Y(I2) 7 CONTINUE 1I(1)=NA(2)+IC2 1I(2)MNA(1)+IC2 6 CONTINUE GO TO 998 101 CONTINUE C 3 POINT DFT C 97

DO 8 J"1,1 Ai(1)uI(l) NA(2)=KI(2) NA(3)-KI(3) DO 9 IG50,ILIM IO.IG*IC1 II-NA(1)+IO I2-A (2)+IO I3MA (3)+IO C*******3-POINT DFT iERJAL*********************************************C Tl(X(I2)-X(I3))*C31 U1-(Y(I2)-Y(I3))*C31 R1=X(I2)+X(I3) S1Y(I2)+Y (I3) T2I(I1) -R1*C32 U2=Y(I1)-S1*C32 I(MA(MB(1) )+IO)-X(I1)+Ri Y(MA(MB(1))+IO)=Y(I1)+S1 X(NA(MB(2) )+IO)-T2+U1 I(MA(MB(3))+IO)=T2-U1 Y(NA(MB(2))+IO)-U2-T1 Y(NA(NB(3) )+IO)=U2+T1 9 COITIIUE KI(1)-MA(3)+IC2 KI(2)-NA(1)+IC2 KI(3)=NA(2)+IC2 8 CONTINUE GO TO 998 102 COITINUE C**********************************************************************C C 4 POINT DFT C DO 10 J=1,K MA(1)-KI(1) NA(2)=KI(2) MA(3)=KI(3) MA(4)KI(4) DO 11 IGmO,ILIN IO-IG*IC1 I1=NA(1)+IO I21MA(2)+IO I3=N1(3)+IO I4NMA(4)+IO C*******4-POINT DFT *ERJIL********************************************C Rl=X(Il)+X(I3) R2=1(I1)-X(I3) Sl1Y(I1)+Y(I3) S2=Y(I1)-Y(I3) R3=X(I2)+X(I4) R4=I(I2)-X(I4) S3=Y(I2)+Y(I4) S4=Y(12)-Y(I4) Z(xA(MB(1))+IO)-Rl+R3 x(MA(NB(3))+IO)-R1-R3 Y(MA(IB(l))+IO)=S1+S3 Y(MA(MB(3))+IO)=S1-S3 X(MA(MB(2))+IO)-R2+S4 I(MA(MB(4))+IO)=R2-S4 Y(MA (MB (2)) +IO)=S2-R4 98

Y(A (MB (4)) +IO) =S2+R4 C*$****$ *$**~*$************* **************************$****************C 11 CONTINUE KI (1)=A(4)+IC2 KI(2) NA(1)+IC2 KI (3)-NA(2)+IC2 KI (4)-MA(3)+IC2 10 CONTINUE GO TO 998 103 CONTINUE C 5 POINT DFT C C***********************************************************************C DO 12 Jl1,K MA(1)=KI(1) MA(2)=KI(2) NA(3)-KI(3) 1U(4)-KI(4) MA(5)-KI(5) DO 13 IG-0,ILIM IO-IG*ICI IlI-A(1)+IO I2-MA(2)+IO I3NMA(3)+IO I4nNA(4)+IO I5SA(5)+IO C*******5-POINT DFT *ERJAL*********************************************C R1lX(I2)+X(I5) R2=I(I2)-x(I5) SlY(I2)+Y(I5) S2=Y(I2)-Y(I5) R3-X(I3)+X(I4) R4X. (I3)-X(I4) S3Y (I3)+Y(I4) S4=Y(I3)-Y(I4) T1=(R2+R4)*C51 U1=(S2+S4)*C51 R2~T1-R2*C52 S2=U1-S2*C52 R4=T1-R4*C53 S4-U1-S4*C53 T1- (R1-R3)*C54 Ul"(S1-S3)*C54 T2R1+R3 U2-S1+S3 IX(A(NB(1))+IO)=X(I1)+T2 Y(MlA(B(1))+IO)-Y(I1)+U2 T2X(MA(MB(1) )+IO)+T2*C55 U2Y(AU (B( 1 ) ) +IO)+U2*C55 R1-T2+T1 R3sT2-T1 SI"U2+Ul S3-U2-U1 (MA (B (2)) +IO) =R1+S4 I(NA(NB(5))+IO)'Rl-S4 Y(MA(MB(2) )+IO)=S1-R4 Y(NA (MB(S))+IO)-S1+R4 X(NA(MB(3))+IO)-R3-S2 X(MA(MB(4))+IO)=R3+S2 Y (HA (MB (3)) +IO) "S3+R2 99

Y(nA(NB(4))+IO)sS3-R2 13 CONTINUE KI(1)-=A(5)+IC2 KI(2)-NA(1)+IC2 KI(3)-NA(2)+IC2 KI(4)-NMA(3)+IC2 KI(5)=NA(4)+IC2 12 CONTINUE GO TO 998 104 CONTINUE C 7 POINT DFT C DO 14 J-1,1 NA(1)-KI(1) MA(2)-KI(2) NA(3)-KI(3) MA(4)-KI(4) MA(5)-KI(5) MA(6)-KI(6) NA(7)-KI(7) DO 15 IGaO,ILIM IO-IG*IC1 I1-MA(1)+IO I2-PA(2)+IO I3MNA(3)+IO I4-MA(4)+IO I5-MA(5)+IO I6MNA(6)+IO I7-MA(7)+IO C*******7-POINT DFT *ERIAL*********************************************C Rl-X(I2)+X(I7) R2=X(I2)-X(I7) S1-Y(I2)+Y(I7) S2=Y(I2)-Y(I7) R3IX(I3)+X(I6) R4=X(I3)-X(I6) S3uY(I3)+Y(I6) S4Y (I3)-Y(I6) R5=X(I4)+X(I5) R6=X(I4)-X(I5) S5=Y(I4)+Y(I5) S6=Y(I4)-Y(I5) TI-R1+R3+RS UlS1+S3+S5 X(MA(MB(i)>+IO)-x(Il)+T1 Y(MA(MB(1))+IO)-Y(I1)+U1 T1-XA(M (B() )+IO) +T1*C71 U1-Y(<A(MB(1) ) +IO)+U1*C71 T2-(Rl1-R5)*C72 U2U(S1-SS)*C72 T3. (R65-3)*C73 U3-(S5-S3)*C73 T4. (R3-R1)*C74 U4=(S3-S1)*C74 R1-Tl+T2+T3 R3-T1-T2-T4 RST1-T3+T4 S1IU1+U2+U3 100

S3-U1-U2-U4 S5SU1-U3+U4 U1(S2+S4-S6)*C75 T1w(R2+R4-R6)*C75 T2-(R2+R6)*C76 U2-(S2+S6)*C76 T3 (R4+R6)*C77 U3-(S4+S6)*C77 T4-(R4-R2)*C78 U4-(S4-S2)*C78 R2-T1+T2+T3 R4-T1-T2-T4 R6-T1-T3+T4 S2-U1+U2+U3 S4-U1-U2-U4 S6.U1-U3+U4 X(MA(MB(2))+IO)=R1+S2 XI(A(MB(7))+IO)=R1-S2 Y(pA(MB(2))+IO)=Sl-R2 Y(MA(!B(7))+IO)=S1+R2 XI(A(NB(3))+IO)-=R3+S4 XI(A(B (6)) +IO)-R3-S4 Y(MA(MB(3))+IO)=S3-R4 Y(UA(MB(6))+IO)-S3+R4 X( A(MB(4) )+IO)-R5-S6 XI(A(pB(5))+IO)=R5+S6 Y(IA(MB(4))+IO)-S5+R6 Y(MA(lB(5))+IO)=S5-R6 15 CONTINUE KI (1)-A(7)+IC2 KI(2)=MA(1)+IC2 KI (3)~uMA(2)+IC2 KI(4)-MA(3)+IC2 KI(5)*!A(4)+IC2 KI(6)=NA(5)+IC2 KI(7)-MA(6)+IC2 14 CONTINUE GO TO 998 105 CONTINUE C 8 POINT DFT C DO 16 J-1,K MA(l)KI(l) MA(2)-KI(2) MA(3)-KI(3) M(4)=KI(4) KA(5)=KI(5) NA(6)-KI(6) A (7)TI1(7) MA(8)-1I(8) DO 17 IGO,ILIK IO-IG*IC1 II-nA(1)+IO I2-NA(2)+IO I3-MA(3)+IO I4=MA(4)+IO I5-A(5)+IO I6=MA(6)+IO 101

I7-MA(7)+IO I8SMA(8)+IO C*******8-POIIT DFT KERIAL*** Rl1X(I1)+I(I5) R2-I(I1)-X(I5) Sl-Y(I1)+Y(I5) S2-Y(I1)-Y(I5) R3-X(I2)+X(I8) R4-X(I2)-X(I8) S3-Y(I2)+Y(I8) S4=Y(I2)-Y(I8) RSX(I3)+X(I7) R6X1(I3)-X(I7) SS-Y(I3)+Y(I7) S6=Y(I3)-Y(I7) R7iX(I4)+X(I6) R8-X(I4)-X(I6) S7Y (I4)+Y(I6) S8-Y(I4)-Y(I6) TI=Rl+R5 T2-R1-R5 U1'S1+S5 U2-S1-S5 T3-R3+R7 R3-(R3-R7)*C81 U3=S3+S7 S3-(S3-S7)*C81 T4-R4-R8 R4-(R4+R8)*C81 U4-S4-S8 S4=(S4+S8)*C81 TSR2+R3 T6-R2-R3 U5SS2+S3 U6"S2-S3 T7=R4+R6 T8SR4-R6 U7=S4+S6 U8~S4-S6 I(MNA(B(1))+IO)=Ti+T3 IX(A(NB(5))+IO)-T1-T3 Y(NA(NB(1))+IO)-Ul+U3 Y(MA(MB(5))+IO)=U1-U3 X(HA(MB(2))+IO)-T5+U7 X (MA (B(8))+IO)=T5-U7 Y(MA(B (2))+IO)=U5-T7 Y(MA(MB(8))+IO)=U5+T7 I(NA(MB(3))+IO)=T2+U4 I(A(iB(7)) +IO)=T2-U4 Y(MA(NB(3))+IO)=U2-T4 Y(MA(MB(7)) +IO)=U2+T4 I(MA(HB(4))+IO)=T6+U8 IX(A(MB(6))+IO) T6-U8 Y(MA(MB(4))+IO) U6-T8 Y(A A(MB(6))+IO) U6+T8 17 CONTINUE KI(1)-MA(8)+IC2 KI(2)=NA(i)+IC2 KI(3)=NA(2)+IC2 102

KI (4)NA (3)+IC2 I (5) -A (4)+IC2 KI (6)=MA(5)+IC2 KI (7)MNA(6)+IC2 KI (8)-NM(7)+IC2 16 COITINUE GO TO 998 106 CONTINUE C 9 POINT DFT C DO 18 Js1,K MA(1)=KI(1) MA(2)=KI(2) NA(3)-KI(3) MA(4)-KI(4) MA(5)-KI(5) MA(6)KI(6) MA(7)KI(7) MA(8)=KI(8) MA(9)=KI(9) DO 19 IG=O,ILIM IO'IG*IC1 IlNA(1)+IO I2'MA(2)+IO I3-MA(3)+IO I4AN1(4)+IO I5'A(5)+IO I6-NA(6)+IO I7=NA(7)+IO I18NA(8)+IO I9i-A(9)+IO C*******9-POINT DFT *ERNAL*********************************************C RI —(I2)+X(I9) R2=I(I2)-X(I9) Sl=Y(I2)+Y(I9) S2=Y(I2)-Y(I9) R3=I(I3)+I(I8) R4X (I3)-X(I8) S3=Y(I3)+Y(I8) S4-Y(I3)-Y(I8) R5-X(I4)+X(I7) T-(X(I7)-1(I4))*C31 S5-Y(I4)+Y(I7) U-(Y(I7)-Y(I4))*C31 R7=X(I5)+X(I6) R8X1(I5)-X(I6) S7-Y(I5)+Y(I6) S8=Y(I5)-Y(I6) R9=X(I1)+R6 S9=Y(I1)+S5 T1-X(I1)-R5*C32 U1-Y(I1)-S5*C32 T2=(R3-R7)*C92 U2-(S3-S7)*C92 T3= (R1-R7)*C93 U3= (S1-S7)*C93 T4=(R1-R3)*C94 U4=(S1-S3)*C94 RlO=R1+R3+R7 103

S10-S1+S3+S7 RIT1+T2+T4 R3-T1-T2-T3 R7=T1+T3-T4 SlUI+U2+U4 S3=U1-U2-U3 S7=U1+U3-U4 I(MA(MB(l))+IO)=R9+R10 Y(MA(MB(1))+IO)=S9+S10 RSR9-RlO*C32 S5=S9-S1O*C32 R6=(R4-R2-R8)*C31 S6=(S4-S2-S8)*C31 T2=(R4+R8)*C96 U2=(S4+S8)*C96 T3-(R2-R8)*C97 U3=(S2-S8)*C97 T4=(R2+R4)*C98 U4-(S2+S4)*C98 R2-T+T2+T4 R4sT-T2-T3 R8-T+T3-T4 S2=U+U2+U4 S4=U-U2-U3 SSU+U3-U4 I(MA(MB(2))+IO)=R1-S2 X(MA(MB(9))+IO)=Rl+S2 Y(NA(NB(2))+IO)=Sl+R2 Y(NA(MB(9))+IO)=S1-R2 I(MA(MB(3))+IO)=R3+S4 X(MA(NB(8))+IO)=R3-S4 Y(NA(NB(3))+IO)=S3-R4 Y(IA(NB(8))+IO)=S3+R4 XI(A(MB(4))+IO)aR5-S6 X(MA(MB(7))+IO)=R5+S6 Y(KA(MB(4))+IO)-S5+R6 Y(MA(MB(7))+IO)=S5-R6 X(IA(MB(5))+IO)=R7-S8 XI(A(MB(6))+IO)=R7+S8 Y(NA(MB(5))+IO)=S7+R8 Y(NA(OB(6))+IO)=S7-R8 19 CONTINUE KI(1)=MA(9)+IC2 KI(2)=MN(1)+IC2 KI (3)=MA(2)+IC2 KI(4)-NA(3)+IC2 KI(5)NMA(4)+IC2 KI(6)=KA(s)+IC2 KI(7)T=!A(6)+IC2 KI(8)=!A(7)+IC2 KI(9)=KA(8)+IC2 18 CONTINUE GO TO 998 107 CONTINUE C 16 POINT DFT C DO 20 J1,K MA(1)=KI(1) 104

MA(2)=KI(2) MA(3)-KI(3) MA(4)-KI(4) MA(5)KXI(5) MA(6)-KI(6) A(7)-xI(7) NA(8)=XI(8) MA(9)-KI(9) NA(10)~KI(10) MA (11) -XI (11) NA(12) KI(12) NA(13)KI (13) MA(14)=KI(14) MA(15)=KI(15) MA(16)=KI(16) DO 21 IG=0,ILIM IO-IG*IC1 INlMA(1)+IO I2=NA(2)+IO I3=NA(3)+IO I4-NA(4)+IO I5=MA(5)+IO I6NMA(6)+IO I7-MA(7)+IO I8NAU(8)+IO I9,MA(9)+IO I10lOA(10)+IO Ill-nA(11)+IO I12WNA(12)+IO I13=MA(13)+IO I14iNA(14)+IO I5IaM(15)+IO I16UMA(16)+IO C*******16-POIIT DFT KERJIL********************************************C RI=X(I1)+X(I9) R2"X(I1)-X(I9) S1-Y(I1)+Y(I9) S2-Y(I1)-Y(I9) R3X,(I2)+x(I10) R4=X(I2)-X(I10) S3-Y(I2)+Y(I10) S4-Y(I2)-Y(I10) R5S-X(I3)+X(I11) R6-X(I3)-X(I11) S5SY(I3)+Y(I11) S6=Y(I3)-Y(I11) R7=X(I4)+X(I12) R8,X(I4)-X(I12) S7-Y(I4)+Y(I12) S8.Y(I4)-Y(I12) R9,X(I5)+x(I13) RlO-X(I5)-X(I13) S9,Y(I5)+Y(I13) S10,Y(IS)-Y(I13) RlllX(I6)+X(I14) R121i(I6)-X(I14) SlilY(16)+Y(I14) S12-Y(I6)-Y(I14) R13XI(I7)+X(I15) R14XI(I7)-X(I15) 105

S13-Y(I7)+Y(I15) S14=Y(I7)-Y(I15) R15=X(I8)+X(I16) R16=X(I8)-X(I16) S15SY(I8)+Y(I16) S16=Y(I8)-Y(I16) Tl-Rl+R9 T2-R1-R9 UliS1+S9 U2"S1-S9 T3-R3+R1l T4=R3-Rll U3=S3+S11 U4-S3-Sll T5-R5+R13 T6-R5-R13 U5S5+S13 U6-S5-S13 T7=R7+R15 T8=R7-R15 U7=S7+S15 U8ZS7-S15 T9= (T4+T8)*C81 T10 (T4-T8)*C81 U9 (U4+U8)*C81 U10=(U4-U8) *C81 RlTI+T5 R3-T1-T5 S1IUi+U5 S3-U1-US R5-T3+T7 R7-T3-T7 S5-U3+U7 S7=U3-U7 R9=T2+T10 R11T2-T10 S9-U2+U10 S11-U2-U10 R13"T6+T9 R15"T6-T9 S13-U6+U9 S15=U6-U9 TlR4+R16 T2fR4-R16 UlS4+S16 U2-S4-S16 T3-(R6+R14) *C81 T4 (R6-R14)*C81 U3 (S6+S14)*C81 U4i(S6-S14) *C81 T5S-R8+R12 T61R8-R12 U5SS8+S12 U6,S8-S12 T7T(T2-T6)*C162 T8=T2*C163-T7 T9sT6*C164-T7 TlO1R2+T4 Tll=R2-T4 R2sTlO+T8 106

R4-T10-T8 R6uTll+T9 R8=T11-T9 U7m(U2-U6)*C162 U8SU2*C163-U7 U9-U6*C164-U7 U10S2+U4 U11S2-U4 S2'UlO+U8 S4=U10-U8 S6=Ull+U9 S8SUll-U9 T7=(TI+T5)*C165 T8ST7-T1*C164 T9=T7-T5*C163 TlO-RO1+T3 Tl=R10O-T3 R10-TlO+T8 R12-T1O-T8 R14-T11+T9 R16-Tll-T9 U7=(U1+U5)*C165 U8-U7-C164*U1 U9=U7-C163*U5 UlOSlO+U3 U11-S10-U3 S10-UlO+U8 S12-U1O-U8 S14-U11+U9 S16-U11-U9 I(MA(MB(l))+IO)=Rl+RS I(NA(MB(9))+IO)=R1-R5 Y(NA(MB(1))+IO)=S1+S5 Y(OA(NB(9))+IO)=S1-S5 I(MA(MB(2))+IO)=R2+S10O IX(A(MB(16))+IO)=R2-S10O Y(MA(MB(2))+IO)=S2-R10O Y(IA(MB(16))+IO)=S2+R10 X(RA(MB(3))+IO)=R9+S13 XI(A(MB(15))+IO)=R9-S13 Y(HA(MB(3))+IO)=S9-R13 Y(HA(MB(15))+IO)=S9+R13 I(MA(MB(4))+IO)=R8-S16 I(MA(B (14))+IO)=R8+S16 Y(MA(MB(4))+IO)=S8+R16 Y(MA(MB(14))+IO)=S8-R16 X(HA(MB(5))+IO)=R3+S7 I(MA(MB(13) )+IO)=R3-S7 Y(NA(IB(5))+IO)=S3-R7 Y(llA(MB(13) )+IO)-S3+R7 X(MA(BB(6))+IO)=R6+S14 I(RA(NB(12) )+IO)-R6-S14 Y(,A(IB(6))+IO)=S6-R14 Y(HNA(MB(12) )+IO)-S6+R14 X(NA(MB(7) )+IO)=Rll-S15 X(MA(NB(11) )+IO)=Rll+S15S Y(MA(MB(7))+IO)=Sll+R15 Y(MA(MB(ll))+IO)=S11-R15 I(MA(MB(8))+IO)=R4-S12 I(MA(MB(10) )+IO)-R4+S12 107

Y(IAOIB(8) )+IO)uS4+R12 Y(lkA(NB (10) )+IO)-S4-RI2 21 CONTINUE KI(1)N1U(16)+IC2 KIC2)-NAL(1)+IC2 KII(3)-NA(2)+IC2 KII(4)=NA(3)+IC2 KI (5)-NA(4)+IC2 KII(6)=NA(5)+IC2 KII(7)=NA(6)+IC2 KI (8)=NA(7)+IC2 KII(9)=NA(8)+IC2 KI (10)-NA (9) +IC2 KI (11)-NA(10)+IC2 KI(12)=N&(11)+IC2 KI (13)=NA(12)+IC2 KI(14)NMA(13)+IC2 KI(15)=NA(14)+IC2 KI(16)NMA(16)+IC2 20 CONTINUE 998 CONTINUE 3 CONTINUE 2 CONTINUE IF(SIGN.EQ. 1) THEN DO 22 I-1,L*L Y(I) —Y(I) 22 CONTINUE ELSE END IF RETURN END 108

SUBROUTINE TABLE(L,N,N) C THIS SUBROUTINE RETURNS THE PRIME FACTORS OF L C INTEGER N(* INTEGER NF (69,4) DATA (IF( 1,J),J31,4) /2,0,0,0/ DATA (NF( 2,J),J=1,4) /3,0,0,0/ DATA (NF( 3,J),J=1,4) /4,0,0,0/ DATA (NF( 4,J),J-1,4) /5,0,0,0/ DATA (IF( 5,J),J=1,4) /2,3,0,0/ DATA (IF( 6,J),J-i,4) /7,0,0,0/ DATA (IF( 7,J),J-1,4) /8,0,0,0/ DATA (NF( 8,J),J-1,4) /9,0,0,0/ DATA (IF( 9,J),J-1,4) /2,5,0,0/ DATA (NF(1O,J),J=1,4) /3,4,0,0/ DATA (NFF(11,J),J-1,4) /2,7,0,0/ DATA (NF(12,J),J=i,4) /3,5,0,0/ DATA (NF(i3,J),J=1,4) /16,0,0,0/ DATA (lF(14,J),J=1,4) /2,9,0,0/ DATA (NF(15,J),J=1,4) /4,5,0,0/ DATA (NF(16,J),J=i,4) /3,7,0,0/ DATA (JF(17,J),J-1,4) /3,8,0,0/ DATA (NF(18,J),J-1,4) /4,7,0,0/ DATA (NF(19,J),J-1,4) /2,3,5,0/ DATA (NFF(20,J),J-1,4) /5,7,0,0/ DATA (NFF(21,J),J-1,4) /4,9,0,0/ DATA (NF(22,J),J=1,4) /5,8,0,0/ DATA (NF(23,J),J=1,4) /2,3,7,0/ DATA (NF(24,J),J-1,4) /5,9,0,0/ DATA (NF(25,J),J-1,4) /3,16,0,0/ DATA (NF(26,J),J=1,4) /7,8,0,0/ DATA (NF(27,J),J=1,4) /3,4,5,0/ DATA (NF(28,J),J=1,4) /7,9,0,0/ DATA (NF(29,J),J=1,4) /2,5,7,0/ DATA (NF(30,J),J=1,4) /8,9,0,0/ DATA (NF(31,J),J=1,4) /5,16,0,0/ DATA (NF(32,J),J=1,4) /3,4,7,0/ DATA (NF(33,J),J=1,4) /2,5,9,0/ DATA (NF(34,J),J-1,4) /3,5,7,0/ DATA (NF(35,J),J-1,4) /7,16,0,0/ DATA (NF(36,J),J-1,4) /3,5,8,0/ DATA (NF(37,J),Jinl,4) /2,7,9,0/ DATA (NF(38,J),J-1,4) /4,5,7,0/ DATA (NF(39,J),J-1,4) /9,16,0,0/ DATA (NF(40,J),J-1,4) /3,7,8,0/ DATA (NF(41,J),J-1,4) /4,5,9,0/ DATA (NF(42,J),J=1,4) /2,3,5,7/ DATA (NF(43,J),J-1,4) /3,5,16,0/ DATA (NF(44,J),J-1,4) /4,7,9,0/ DATA (NF(45,J),Jn1,4) /5,7,8,0/ DATA (NF(46,J),J-1,4) /5,7,9,0/ DATA (NF(47,J),Jni,4) /3,7,16,0/ DATA (NF(48,J),J-1,4) /5,8,9,0/ DATA (NF(49,-J),J-1,4) /3,4,5,7/ DATA (NF(50,J),J=1,4) /8,7,9,0/ DATA (NF(51,J),J-1,4) /5,7,16,0/ DATA (NFCS2,J),J-1,4) /2,5,7,9/ DATA (NF(53,J),J=1,4) /5,9,16,0/ DATA (NF(54,J),J-1,4) /3,5,7,8/ 109

DATA (N[F(55,J),J=1,4) /7,9,16,0/ DATA (IF(56,J),J=1,4) /4,5,7,9/ DATA (NF(57,J),J-1,4) /3,5,7,16/ DATA (NF(58,J),J-1,4) /5,7,8,9/ DATA (NF(59,J),J-1,4) /5,7,9,16/ IF(L.EQ. 2) THEN M1= X1= ELSE IF(L EQ. 3) THEN Hal 1=2 ELSE IF(L.EQ. 4) THEN Hal 1=3 ELSE IF(.EQ. 5) THEN Mini 1-4 ELSE IF(L.EQ. 6) THEE H=2 laS ELSE IF(L.EQ. 7) THEN N=2 1=6 ELSE IF(L.EQ. 8) THEN Hal 1=7 ELSE IF(L.EQ. 9) THEN Hal 1=8 ELSE IF(L.EQ. 10) THEN H=2 1=9 ELSE IF(L.EQ. 12) THEN H=2 1=10 ELSE IF(.EQ. 14) THEN H-2 1=11 ELSE IF(L.EQ. 15) THEN N=2 1=12 ELSE IF(L.EQ. 16) THEN Hal 1=13 ELSE IF(L.EQ. 18) THEN H=2 1=14 ELSE IF(L.EQ. 20) THEN H-2 1=15 ELSE IF(L.EQ. 21) THEN KH2 1=16 ELSE IF(L.EQ. 24) THEN H-2 1=17 ELSE IF(L.EQ. 28) THEN H-2 1=18 ELSE IF(L.EQ. 30) THEN 110

M'3 K-19 ELSE IF(L.EQ. 1*2 X=20 ELSE IF(L.EQ. M*2 K121 ELSE IF(L.EQ. M-2 1K22 ELSE IF(L.EQ. 1M3 K"23 ELSE IF(L.EQ. M*2 1K24 ELSE IF(L.EQ. M12 K125 ELSE IF(L.EQ. M-2 K126 ELSE IF(L.EQ. M-3 K=27 ELSE IF(L.EQ. M-2 1K28 ELSE IF(L.EQ. M-3 1K29 ELSE IF(L.EQ. M!2 K130 ELSE IF(L.EQ. M-2 K131 ELSE IF(L.EQ. 1-3 K-32 ELSE IF(L.EQ. Ml3 K=33 ELSE IF(L.EQ. M-3 K=34 ELSE IF(L.EQ. MI2 K=35 ELSE IF(L.EQ. M"3 K-36 ELSE IF(L.EQ. M-3 Ks37 ELSE IF(L.EQ. M-3 K138 ELSE IF(L.EQ. 35) THEN 36) THEN 40) THEN 42) THEN 45) THEN 48) THEN 56) THEN 60) THEN 63) THEN 70) THEi 72) THEN 80) THEN 84) THEI 90) THEN 105) THEI 112) THEN 120) THEN 126) THEN 140) THEN 144) THEI 111

Ms2 Ks39 ELSE IF(L M-3 K=40 ELSE IF(L M=3 K-41 ELSE IF(L M-4 K=42 ELSE IF(L M=3 K=43 ELSE IF(L M-3 Ks44 ELSE IF(L M-3 K-45 ELSE IF(L M=3 K=46 ELSE IF(L M=3 K=47 ELSE IF(L M-3 K-48 ELSE IF(L M-4 K149 ELSE IF(L M=3 K-50 ELSE IF(L M=3 K=51 ELSE IF(L M=4 K=52 ELSE IF(L N=3 K1=53 ELSE IF(L M=4 K=54 ELSE IF(L M=3 K=55 ELSE IF(L M=4 1=56 ELSE IF(L M=4 K=57 ELSE IF(L M=4 K158 ELSE IF(L.EQ..EQ..EQ..EQ..EQ..EQ..EQ..EQ..EQ..EQ..EQ..EQ..EQ..EQ..EQ..EQ..EQ. 168) THEN 180) THEN 210) THEN 240) THEN 252) THEN 280) THEN 315) THEN 336) THEN 360) THEN 420) THEN 504) THEN 560) THEN 630) THEN 720) THEN 840) THEN 1008) THEN 1260) THEN EQ. 1680) THEI EQ. 2520) THEN EQ. 5040) THEN 112

M14 K659 ELSE END IF DO 1 JI1,4 I(J)-IF(K,J) 1 CONTINUE RETURN END 113

IVERSITY OF MICHIGAN 3 9015 03024 4167