A Conjugate Gradient FFT Method for the Computation of the Scattering by Thin Planar Material Plates April 1988 T. J. Peters and J. L. Volakis Radiation Laboratory Dept. of Electrical Engineering and Computer Science The University of Michigan 1301 Beal Avenue Ann Arbor, MI 48109-2122 Technical Repdrt No. 389604-2-T Prepared for General Dynamics Fort Worth Division P. 0. Box 748 Fort Worth, Texas 76101

Abstract A conjugate gradient FFT formulation and implementation of the scattering by planar material plates is presented. The plates considered are of arbitrary material composition and periphery. A substantial part of this investigation concentrated on the development of efficient and higher accuracy FFT algorithms. This resulted in the generation of a scattering code 3 to 4 times faster than the traditional and, in addition, a higher convergence was achieved with the incorporation of basic functions for the representation of the plate's current density. ii

TABLE OF CONTENTS Page# CHAPTER I. INTRODUCTION 1 Motivation 1 Literature Review 1 Objectives 2 CHAPTER II. DERIVATION OF THE CONJUGATE GRADIENT METHOD 4 Introduction 4 Derivation 5 Typical Operator Equations 11 CHAPTER III. DERIVATION AND COMPUTATION OF A DISCRETE FOURIER TRANSFORM USING HIGHER ORDER INTEGRATION AND PRIME FACTORIZATION 13 Introduction 13 Derivation of a DFT 13 Higher Order Integration 18 Derivation of a Prime Factor FFT 26 Extension to Two Dimensions 32 Test Results 34 CHAPTER IV. SCATTERING FROM A SINGLE PLANAR PLATE 37 Orientation 37 Integral Equations 38 Surface Generator 42 Perfectly Conducting Plate 43 Thin Dielectric Plate 45 Thin Dielectric and Magnetic Plate 47 Calculation of Radar Cross Section 52 CHAPTER V. RESULTS 54 CHAPTER VI. FUTURE WORK 65 iii

CHAPTER I INTRODUCTION L.1 Motivation The material plate scattering problem is one of the most commonly occurring scattering phenonena which does not have a known analytic solution. Many practical electrically thin structures have edges and corners which provide a significant contribution to the scattered fields. For electrically large surfaces the scattering contributions due to edges and corners are fairly localized so that computation of the surface currents and scattered fields from plates will give much insight on the scattering behavior of more arbitrary shapes. One on the prime reasons for studying the scattering characteristics of material plates is to understand how the material distribution effects the scattered fields. Not only is it important to know the scattered fields but it is also desirable to know how to change the material composition of a scatterer to achieve a certain scattering characteristic. 1.2 Literature Review At the present time efficiency is still a prime concern for numerical solution of scattering by material plates with above resonant dimensions. This efficiency is 1

2 measured in terms of minimizing computational factors such as time, storage and cost. The purpose of this study is to develop a method which has the potential to solve a wide class of scattering problems in an efficient manner. Newman et al [1] and Naor et al [2] treat rectangular material plates, however, the basis functions are inappropriate for curved perimeters and the boundary conditions employed are not valid at edge on incidence with H-polarization. Good results have been obtained for perfectly conducting plates by Glisson et al [3] and Rao et al [4] using triangular cells and linear basis functions. Unfortunately, these methods rely on matrix solutions which may become restrictive for a body of above resonant dimensions. Studies on perfectly conducting wire and plate scatterers by Sarkar et al [5] and work by Christodoulou et al [6] on meshes combine the method of conjugate gradients [7]-[10] with the FFT [11]-[12] to solve operator equations directly without storing a large matrix. 1.3 Objectives The major goal of this study is to develop an efficient numerical modeling technique to solve for the induced current and the scattered field from a material plate illuminated by a plane wave. The technique involves combining a conjugate gradient (CG) method with a fast Fourier transform (FFT) to solve a set of coupled convolution integral equations. Achieving an efficient ailgorithm requires a good understanding of both a conjugate gradient method and the computation of the FFT. Chapter II gives a complete derivation of a conjugate gradient method. Chapter III derives the discrete Fourier transform (DFT) from the continous Fourier transform and describes how computation may be accomplished using various integration formulas combined with the FFT. Chapter IV describes the single material plate which includes the zero thickness perfectly conducting plate, the electrically thin dielectric plate, and

3 the combination dielectric and magnetic plate. Chapter V gives some preliminary numerical results and chapter VI describes some future work which will be an extension of this study.

CHAPTER II DERIVATION OF A CONJUGATE GRADIENT METHOD Introduction The conjugate gradient method seeks the solution to the operator equation Az =b (2.1) where A is an operator which is assumed to be non-singular. In order to derive the solution technique the concept of an inner product must be discussed. Let (f, g) denote the inner product of the two functions f = fr + jfi and g = g, + jgi. If the functions are continuous, then given a region s E R2 and real functions fr,i, g7,i E R2 a suitable inner product could be given as (f,g) = ff(x,y)g*(x,y)w(x,y)ds (2.2) where * denotes the complex conjugate and w(x, y) is a real weight function. It is most common to define w(x, y) as 1 for differential operators W(x, y) = N-i M-1 (2.3) 1 Z S{(x - maxr, y - nay) for integral operators n=O m=O 4

5 If the functions are discrete n x 1 vectors such that f = [fo fi f2i.. fN-l]T and g = [go 91 92..* N-1]T a suitable inner product could be given as N-i (f,g) = g fn. (2.4) n=O Closely related to the inner product is the Euclidean norm of a function defined by If112 = (,). - (2.5) Another concept that must be introduced is the adjoint operator of A denoted by Aa which satisfies the relation (Ap,q) = (p,Aaq) (2.6) for any p and q. Derivation With an understanding of these concepts, the following observation may be made. The solution of (2.1) will yield the same solution as the minimum of some quadratic functional. That functional may be written as I(z) = (b-Az,b-Az) = (Az, Az) - (b, Az) - (Az, b) + (b, b) = (Az, Az)- (Az, b)*- (Az, b) + (b, b) = 2 [(Az,Az)-Re(Az,b)+ (b,b)] = 2 [(z, AaAz) - Re(z, Aab) + I(b, b). (2.7) The method of conjugate method of conjugate gradients [7-10] is a general solution technique for minimizing I. Making the substitutions B = AaA, h = Aay and noting that the

6 solution is independent of the multiplicative factor 2 and the constant (b, b) the solution z = A-lb = B-lh also minimizes the function F(z) given by F(z) = 2(Bz, z) - Re(h, z). (2.8) If A is non-singular then B will have the following properties [13]. 1. Hermitian B = Ba 2. positive definite (Bz, z) > 0 z: 0 3. eigenvalues are positive and real and eigenvectors are orthogonal The solution is assumed to have the recursive form Zk+l = Zk + (akPk (2.9) where ak is a real positive constant and pk is called the search vector. This yields an expression for the residual rk+l defined by rk+-1 = h- Bzk+ h - B(zk + akpk) = rk - kBpk (2.10) Substituting (2.9) into (2.8) F(k+l) = F(zk) + alk(Bpk,pk) - akRe(rk,pk). (2.11) To minimize this function with respect to ak set - F =0 (2.12) aak and solve for Re(rk,pk) (Bpk,pk) (2.13)

7 Substituting (2.13) into (2.11) 1 (Re(rk,pk))2 F(zk+l) = F(zk) ( (k, pk ) (2.14) 2 (Bpkpk) The function decreases at each step if Re(rk,pk) #7 0 and B is positive definite. However, I lb - Az l2 may not decrease at each step since it is not the function to be minimized. This local minimization does not guarantee that the solution is obtained in a finite number of steps. The crux of the conjugate gradient method is that global minimization in a finite number of steps may be achieved if the search vectors Pk are chosen correctly. In order to choose pk it is informative to view the minimization from a geometric point of view. The vector pk is a direction in n-space. The residual vector rk+l is proportional to the difference between the exact solution and the k + 1 approximation. Thus r+l = h - Bzk+l B(z - Zk+l) (2.15) where z is the true solution z = [z' z2 3...,n]T. (2.16) If Pk is chosen such that (rk+l,Pk) = 0 then the minimization occurs in an ndimensional plane and since (rk+,pk) = 0 = (rk - kBpk,pk) = rk- Re(rk,pk )BP) (Bpkpk) = (rk,pk)- Re(rk, pk) (2.17) the inner product (rkpk) is real. Thus (rk+l,Pk) is recognized as the equation of an n-dimensional plane given by (B(z- Zk+1),pk) = 0 (2.18)

8 which may be rewritten as z1 + + 7z2 + 73z3 +...+ nZ =/1 (2.19) which establishes the relationship between a linear equation and an n-dimensional plane. This idea may be extended to a system of equations by noting that the solution of a system of m < n equations in n unknowns is the same as finding the intersection of m n-dimensional planes. With this in mind it can be surmised that if the search directions were generated such that (rk+l, p) = 0 for s = 1,.. k then at the kth iterative step the function would be minimized over the intersection of k planes. Another interpretation is that at each iteration a least squares solution is found to a system of m equations in n unknowns. The iterative nature of the solution stems from the fact that at each step the order of the leasts squares solution is increased and a more accurate solution is obtained. Criteria for chosing pk may be found by expanding the product of the residual and all previous search directions. (rk+l,Pk-1) = (rk,pk-1) - ak(Bpk,pk-1) (rk+lPk-2) = (rk-1,Pk-2) - ak-l(BPk- 1,Pk-2) - Ck(Bpk, Pk-2) (rk+lPk-3) = (rk-2,Pk-3) - ak2 (BpkP2,pk-3) - ak-1.(Bpk-l,pk-3) - ok (Bpk,pk-3) n-1 (rk+l,Pk-n) = -E ak. (Bpk-s,Pk-n) (2.20) s=0 Observation of (2.20) indicates that if the method is to reach an exact solution in n iterations, the condition (Bpk,ps) = 0 s = 0,..., k-1 (2.21) must be enforced. Equation (2.21) is recognized as the condition for the vectors pk to be "B-orthogonal" or "mutually conjugate". Vectors with this property may

9 be generated by the Gram-Schmidt process. Given a set of linearly independent vectors vl,..., v,n a set of "B-orthogonal" vectors may be generated as Pi = Vi P2 = V2 -bllpl P3= V3 -b2lP- b22P2 k Pk+l = Vk+l-E I kPs (2.22) s=1 where Ask = (Bvk+l,) (2.23) ljksI = /Q, ~ \' (2.23) (Bpp)'P) Note that this method requires k matrix vector products for the kth iteration so it is inherently very inefficient for large k. However, if the vectors vk are chosen to be the residuals rk then all 3ks = 0 for s = 1,..., k -1 which leaves only one matrix vector product per iteration. To show this, the relationship between residuals must be established. Taking the adjoint of (2.22) and multiplying by the residual r, yields k (rc, pk+1) = (r, rk+l) - E1 3k(rc, p). (2.24) s=1 When c = k+ 1 the products (r, pa) = 0 for s < c by virtue of (2.21) which implies (rk+l, pk+l) = (rk+, rk+l) (2.25) and for c = k + 2,...,n (rc,pk+l) = (rc, rk+) = 0. (2.26) Therefore, the residual vectors are orthogonal such that for k = s (r, rk) = 0. (2.27) With (2.27) established take the inner product of rk+l with (2.10) and observe that - -(rk+lr,+ - r,) = (Bark+l,ps) = (Brk+lps) = 0 (2.28) <-<

10 for s = 1,...,k- 1 and (rk+,rk+l) = Brk+l,pk) (2.29) (fk for s = k. Thus k, = 0 fors = 1,...,k- 1. The computations will be minimized if the search directions Pk are normalized with respect to (rk, rk). The algorithm, which will be referred to by this author as the "nested operator" algorithm, is expressed below. initialize rl = h-BXl (2.30) f0 = (2.31) d (ri) (2.31) pi = /orl (2.32) for k 1,...,n ak = (2.33) (Bpk,pk) 2.33) Xk+l = Xk + kPk (2.34) rk+l rk - akBpk (2.35) A = ( 1 ) (2.36) (rk+l, rk+ ) k+The algorithm is terminated at = n or whenk+l (2.37) The algorithm is terminated at k = n or when Ihl1112 < tolerance. (2.38) Excluding the initialization, this algorithm requires one matrix vector product per iteration. The number of multiplications and divisions per iteration is n2 + 5n + 2. However, the matrix B may or may not be explicitly known. In this case the residual expressed in (2.10) must be replaced by Aar with r = b - Az. The modified algorithm referred to by this author as the "split operator" algorithm is

11 initialize rl = y-Axi (2.39) (0o Aa, - (2.40) (Aarl, Aarl) pi = /oAarl (2.41) for k 1,...,n ak (2.42) (Apk, Apk (2.42) Xk+1 = k + akpk (2.43) rk+l = rk - kApk (2.44),k = (2.45) (Aa/rk+l, Aark+l (2.45) Pk+l = Pk + kArk+l (2.46) The algorithm is terminated at k = n or when 1rk+12 < tolerance. (2.47) IYIl 2 In this case the magnitude of the residual will decrease at each step since it represents the actual function being minimized. Excluding the initialization, this algorithm requires two matrix vector products per iteration. The number of multiplications and divisions per iteration is 2n2 + 5n + 2. Typical operator equations The most common operators will be of integral or differential type, or a matrix derived from one or the other. If the operator A is composed of Fredholm equations of the second kind lf(x, y) + J [f,(x', y')li(X, X', y, y') + fy(x', y')2(xx,, y')] ds' = h(x, y) (2.48) fy(x, y) + f [f(x', yt)3(x, x, y, y') + f(x', y')(x, x', y, y)] ds' = h(x, y) (2.49)

12 then the adjoint operator A"f would be defined by the relations;f(x, y) + IJ [fX(x', y')l(x', x, y'y) + fy(x', y')2(x', x,,y'y)] ds' = h(x, y) (2.50) f(x, y) + [f(x', (',3 x(,, y' (, y ) + f(x', y'),(x',x,y',y)] ds' = hy(xy). (2.51) If the operator A is composed of first and second order differential operators [X+ x f(x, y) + aay fy(x, y) = h(x, y) (2.52) xy2(x Y) + [ + 2 A f(xY) = y h,(x,y) (2.53) then the adjoint operator Aaf would be defined by the relations - 2 02 [-x2 ] - f (x, y+ ) + xyf(x ) =, (xy) (2.54) OxOy f(x, y) + [ - ay; fy(x,Z y) = hy(x, y) (2.55) If the operator A is a general complex n x n matrix then the adjoint is the complex conjugate of the transpose defined by A = (AT)*. (2.56)

CHAPTER III DERIVATION AND COMPUTATION OF A DISCRETE FOURIER TRANSFORM USING HIGHER ORDER INTEGRATION AND PRIME FACTORIZATION 3.1 Introduction The discrete Fourier transform (DFT) has been studied extensively by mathematicians, engineers and scientists for many years. The bibliography by Heideman et al [12] contains over 2000 reference papers concerning computation and application of the DFT. Most of the work has been concerned with increasing the speed of the computation. The purpose of this study is to develop a algorithm which is not only fast but also more accurate. 3.2 Derivation of a DFT The DFT and inverse DFT are approximate representations of the continuous transform pair z(f) = j z(x)e-J2~fxdx (3.1) boo 13

14 The transform pair is assumed to be valid. The interested reader may consult Champeney [14] for a discussion of the sufficient and necessary conditions for transformability. Assume that z(x) is a complex function of bounded support such that z(x) =,r(x) + jZi(x) for min < x < Xax (3 3) t0 otherwise. Under this restriction it may be shown that there exists fmin, fmax and S such that Iz(f)l < 6 for f > fmax and f < fmi,. In a more practical sense z(f) may be defined by -() r(f) + ji(f) for fmn, < f < frmax z(f) = (3.4) 0 otherwise. Let the spatial and spectral domains be segmented into N uniform cells of widths Ax = Xmax - min (3.5) Af = fmafmin (3.6) N Since the sample points are equally spaced, the cell width Ax implies a resolution limit on Ifmax - fminl and Af implies a resolution limit on IXmax - Xmin, These limits may be found by finding the lowest order interpolation sinusoid which passes through all the sample points. VA A A A A 7C0-r=~= B Figure 3.1. Interpolation usoids for spatial or spectral domain. Figure 3.1. Interpolation sinusoids for spatial or spectral domain.

15 Observation of figure 3.1 indicates that for curves A, B, and C there is no ambiguity in terms of differenciating each curve given the set of sample values. However, curves C and D are indistinguishible based on the sample values and furthermore, there exist an infinite number of higher order interpolation sinusoids which pass through all the sample points. Curve C denotes the highest order interpolation sinusoid which may be distinguished from any other lower order sinusoid. The resolution limits are then determined by the spatial or spectral frequency of curve C as fma - fmin - (3.7) Ax max - Xmin = A- (3.8) Af which combined with (3.5) and (3.6) may be used to define a quantity called the space-bandwidth product defined by AxAf = N (3.9) N r Note that the spatial sample interval Ax and the spectral sample interval Af are dependent in a reciprocal manner. This is undesirable since the error in the forward and inverse Fourier transforms will be dependent on Ax and Af respectively. Assume that the error in the forward and inverse DFT is considered minimial if!Ax < 61 and Af < 6b Making Ax smaller will increase N and make the forward Fourier transform more accurate but Af remains constant. In order to increase the accuracy of the inverse Fourier transform, extra cells must be appended in the spatial domain in order to increase N without decreasing Ax. This is commonly refered to as'padding' and in general if the number of non-zero samples is M then the total length N of the data set should satisfy N > 2M. It is convienient to make the integration formula independent of the location of the spatial data along the x axis so making the change of variables x = u + xo

16 and defining the quantities s(u) = z(u + xO) (3.10)'(f) = z(f)ej2rfxo (3.11) allows the transform integrals to be written as J Umax (f) = s(u)e-ij2urfdu (3.12) Umin fjmax. - fmax s(u) = f/ s(f)e2 Ufdf = [J (f)e-'df (3.13) Jfmn Jmin Assuming an even number of cells the spatial domain may be segmented as where the variables are discretized as u = pAu p =,..., N-1 (3.14) N N f vAf v=- +1.... (3.15) 2 2 UO UM-1 UNr-1 U Figure 3.2. Space limited function with N cells and similiarly, the spectral domain may be segmented as \l(/)l fo f Figure 3.3. Band limited functior with N cells

17 The integrals are approximated using an open quadrature rule of the form uemax N-1 ] ~s(u)e-j21fudu = () = s() (p) a(v)WP (3.16) ~min p=O fmagn ff *s (f)e-2rfu df = (p) = [ (+1 ()Ov(p)wV (3.17) where W = e-i=. If the indices of the sequences in (3.18) and (3.19) are interpreted as storage locations the vectors s and s would appear in memory as s(0) s(1) s(2)... s(N)... s(N-1) ~~~2 ~ (3.18) s(2N)... s(-2) s(-1) s(0) s(1) s(2)... s(N) In order to directly overwrite s onto s the negative index s elements may be shifted N elements to the right which yields s(0) s(1) s(2)... s(f) s( +1)... s(N-1) 5(0) s(1) s(2) -... (N) 5(-N+1)... S(-1) This shift has no effect on the complex exponential since Wp(QN) = Wpq. A generalized non-reciprocal DFT pair may then be defined by N-1 S (q) = sn(p)ap(q)Wpq (3.20) p=o N-1 * s(p) = [ sE n(q)/q(p)WPq (3.21) q=O If midpoint integration is used, crp(q) = Ax and 3p(q) = Af yielding the conventional reciprocal DFT pair N-1 (q) = E S(p)WPq (3.22) p=o s(p) = 1 - [eS(q)WPq (3.23) Although the forward and inverse transform may be generalized as above, in many cases it is sufficient to generalize only the forward transform. A common illustration of this is the evaluation of a convolution. The inverse transform of a convolution will in many cases be much smoother than any of its functional components. transformed

18 3.3 Higher Order Integration The midpoint integration used to derive the conventional DFT pair implies that the entire integrand is constant over each cell. However, the cell width is chosen based on the function to be transformed not the integrand in the transform. Even if the function is constant over the cell width, the product of the function and the complex exponential becomes quite oscillatory with increasing frequency. Midpoint integration applied to a single cell yields 1A j2 z(x)e-j2dfxdx t hz(xo)e-j2'fx~. (3.24) Jil The assumption on the integrand is that over the cell [Zr(X) + jzi(x)]e-j2if = ZrO + jzio= constant. (3.25) The system of equations Zr(x) cos(2irfx) - zi(x) sin(2rfx) = ZrO (3.26) Zr(x) sin(27rfx) - zi(x) cos(27rfx) -zio (3.27) may be formed and solved for the unknown real and imaginary parts of z(x) as Zr(x) = ZrO cos(2rfx) - zo sin(27rfx) (3.28) z,(x) = zo sin(2rfx) + zio cos(27rfx). (3.29) As these equations indicate, the implied basis functions are dependent on the frequency. Thus, the approximation that the integrand is constant becomes less valid as the frequency increases. To illustrate this, let h = 2 such that fmax1 = 10 and z= 1.0 +jl.0. 2h

19 1.3 1.1o - - 0.8 o0.6 0.4 -0.3 f=3 -0-f=6 \ -0.5 - * f=10 = I 0.9 I.. 0.024 -0.018 -0.012 -0.006 0.000 0.006 0.012 0.018 0.024 CELL WIDTH Figure 3.4. Real part of z(x) assumed at each f.

20 3.3.1 Weighted Open Newton-Coates Integration The simple integration used to obtain (3.24) and (3.25) does not exploit the fact that part of the integrand of (3.12) and (3.13) is known. The following formulas were derived in a manner analogous to the open Newton-Coats forms[15]. The formulas are derived for equally spaced nodes and are of the weighted open NewtonCoates (WONC) type. The weights are complex functions of f defined as a(f) = ar(f) + ja'(f). The nodes are defined by k = o + kh k=0,...,n-1 (3.30) where Xo = a + h (3.31) where ar(f) = e21fUn Ln(u)e-23fudu (3.32) and Ln(t) denotes the Lagrange polynomials defined by -1 (t - t,) L (t)= (t - (3.33) i=o (tn - t) i3n This study will consider the 1 and 3 point formulas. Making the definition V = irfh, the 1 point formula is given as 1 z(x)e-2rfdx = al(b)z(xo)e-j2 o (3.34) where a() = hsin() (3.35) The composite 1 point rule is given by b N-1 | z(x)e-j2fzdx = a1(() z(xk)e f (3.36) Ja?k=O

21 and the 3 point formula yields b z(x)e~-ifxdx = al(i)z(xo)e-32rfxO + a2(/)z(xl)e-ji2 fx+ Ca3()Z(X2)e-j32f (3.37) where,a.3 = -16h [(342 _ 2) sin(50) + (15-2 - 2) sin(O) + 44,(cos(50) + 2 cos(4,))] (3.38) aib = 6, [(342 - 2) cos(5b)- (15,2 - 2) cos(4) - 44(sin(5) - 2 sin(k))] (3.39) a2(/) = - [(2-542)sin(30)-6 cos(30)] (3.40) aC3() = Ca(f) (3.41) The composite 3 point rule may be written as ff -1 X-1 -o 3 3 z(x)e-j2 fxdx = c1(q() Z Z(X3k)e-j2 fx3k+ a12()-)e-j2 fh Z(x3k+l)e-j2rf:3k /a k=O k=O -1 +-Ct3(?P)e-j4'rfh E Z(X3k+2)e —j2f 3k (3.42) k=O 3.3.2 Weighted Open Finite Difference Integration It is possible to derive integration formulas which use not only function values but also derivative values [16]. The derivatives may then be approximated using finite differences. This integration formula will be referred to as a weighted open finite difference (WOFD) type. s(z)e-2fdx = [aos(xo)'(xo) + s'2s'(xo)] e- 2rfo (3.43) where h = b - a and xo =. The coefficients may be calculated by requiring the formula to be exact for polynomials up to and including x for n = 0,...,2. Thus,

22 this formula has the same order as the 3-pt WONC formula. An alternative and easier derivation of the coefficents is to integrate the interpolating polynomial over the interval. An interpolation scheme which uses both function and derivative values is called Birkoff interpolation. The single point integration formula is a special case of Birkoff interpolation which expands the unknown function in a Taylor polynomial around each point. s(x) = S(Xo) + (x - Xo)s'(o) + -(x- X)2S(Xo) (3.44) By inspection, the weight coefficients are expressed as the integrals ao = ej27fzo b e-ji2fzdx (3.45) a a1 = ei2-fzo~ (x- xo)e- rfxdx (3.46) 1 b 2 Ja which yield ao = hli (3.48) ai = j2h272 (3.49) a = 2-= h3 a2 = ii73 (3.50) where sin(= ) (3.51), cos(4) - sin())'2 = (3.52) - [(2 - 2) sin(+) + 2 cos()] (3. 73 3 3j ---— V (3.53) Extending the single point formula to a composite formula with N points yields [b s(x fN-1 N-1 $ s(X)e-j27rfdx =cO E s(Xn)e-j2fzxn +- C1 Z't(zxn)e-j2rfxn Jza n=O n=O N-1 +al2 E s"(x,)e-j2nf)n. (3.54) n=0

23 The derivatives may be approximated using the center difference derivative rules'(Xn) 1 [S(Xn+l) - (Xn-_)] (3.55) 2h S"(zn) I [s(Xn+) - 2s(Xn) + s(Xn-1)]. (3.56) The Fourier transform of the first and second derivatives is given by -00 2I- [00 j0 s'(x)e-j2'ffxdx 1; h [s(x + h) - s(x - h)] e-j2-f xdx [ej2fh - e -i2h] s(x)e-j2rfz 2h 00-o = j sin(2rfh) J s(x)e-3i2fx (3.57) r0 s'(x)e-2 dx t h2 [s(x + h)-2s(x) +s(x- h)] e-j2rfdx 2h [e2fh - 2 + e-32rh] f s(x)e-j2~f 2h oo [1 - cos(27rfh)] f s(x)e32f (3.58) The integral may be written as /00 [ 1 sn2-1 _ rl s(x)e-32'fxdx [ao + j- a, sin(2) - -2 2[1- cos(24p)] ] s(xn)e-j27fzx +(Io + IM_-)e-J2rxo (3.59) where Io and IM-1 are given by I0o = 1 [s(x-i)- 2- s()- (x 2)]] e2f +al [s'(xo) 2hs(,_ )+ )]- e()+~c1 [S'(Xo) - 2h[5(xl) - (XM-1)+ ]] +a2 [s"(XM_) - 2 [s(x) - 2s(xM-1) + S(X- 2)]] e +Q2 ["(XMo)- h[s(xM1)- 2S(ZM) + S(ZM1)]] elf (3.61) JM-I +1 [,(rM-i)- [S(XM) -S(M-2)] e -jr -1) [h~

24 and represent the end point corrections which replace the midpoint difference formulas by forward and backward differences, respectively. These correction terms may be simplified to Io = a2s"(xo) + act'(xo) + 22 [2a2(2-e)- 2] s(xo) 2h2 [h-1 -+-2(32],. hl -22 [ha, + 2a2] s(x1) (3.62) IM-1 [a22s"(XM-1)+ aCls(XM_-)] e-j2(M-1)L +22 [hal - 2a2(1 - ej2 )] e-j2Mk (XM_i) + [ha, - 2a2] e-j(M-)S(XM-2) (3.63) by excluding all points outside the interval. The derivatives at the initial points may be approximated with forward difference rules 3'(xo) - [s(x2) - 4s(x) + 3s(xo)] (3.64) s"(xo) >. [s(x2) - 2s(Xl) + s(xo)]. (3.65) Similarily, the end point derivatives may be approximated using the backward difference rules 8'(XM-1) 2h [3s(XM-.1) - 45(XM-2) + S(XM-3)] (3.66) 1 S"(XM-1 ); 2 [S(M-1)- - 2s(xM-2) + S(XM-3)]. (3.67) Making the assumption that s(xn) = 0 for n > M - 1 (3.59) may be rewritten as 0oo N-1 f s(x)e-2fxdx = C s(x,)e-j2"f'" -00 n=O + [C25(XO) + (C3(Xl) + 45S(X2)] e —j2x0o + [5s(x1M-l) + C6S(xM- 2) + C75(xM-3)] e-j2 x~(3.68) where, letting 0 = MBb h = h [71-72 sin(21b)-j 73[1- cos(2b)]] (3.69) h[-il 1 1,_~~~~~~~~~~~~3.9 2 12 IrrL-~~~~~

25 2 4[672 sin(2Q) + 73[3 - cos(2b)]] h -j [672[3 + cos(2b)] + 73 sin(2b)] (3.70) 24 (3 = -8 [73-j672] (3.71) 4 = 2[73 - j672] (3.72) (5 = h [6[3 sin(2(0 - b)) + sin(20)]72 + [3 cos(2(0 - b)) - cos(20)]73] +jh [6[3 cos(2(0 - b)) + cos(20)]72 - [3 sin(2(0 - b)) - sin(20)]73] (3.73) -6 = - [372 sin(2( - +y(2(-))+7cos(2-)) +j[372 cos(2(0 - )) - 73 sin(2(O - 3))]] (3.74) h 7 = [672 sin(2(0 - 4)) + 73 cos(2(0 - )) +j [672 cos(2(0 - 4)) - y3 sin(2(0 - 4))]]. (3.75) 3.3.3 Operation Count Table 3.1 shows the operation count comparison for the conventional FFT verses the Ipt and 3pt WONC formulas and the 2nd order Ipt WOFD formula. Type real mult. real add. FFT clN c2N Ipt (WONC) clN + 2N c2N 3pt (WONC) 3clN + 6N c2N + 2N lpt (WOFD) ciN + 24N c2N + 14N Table 3.1. Operation count As the next section will show for 60 < N < 504 the constants cl and c2 approximately range over 3 < c1 < 5 and 18 < c2 < 27. For functions with

26 sharp end point discontinuities the 3pt WONC formula would be more efficient. However, if the function approches zero in a smooth manner at the end points, then the end point corrections of the Ipt WOFD may be neglected and it becomes as efficient as the lpt WONC. 3.4 Derivation of a Prime Factor FFT The FFT is a fast method of computing a summation of the form N-1 z(q) = E Z(p)WP. (3.76) p=o where W = e-jN and N = N1N2N3*- Nm. The basic strategy of computation is to map each one-dimensional index p and q onto a N-dimensional map which allows the summation to be accomplished in a more efficient manner. WP may be factored as Wp = W1 qW2w.. Wmqm (3.77) and the single summation becomes a nested set of m Nk point summations Nm-1 1V3-1 N2-1i Nii(q) = E E E E (p) W1 w222 W3 W.. PWmm (3.78) Pm=O P3=0 LP2=O Lp = where Wk = W^. The indices p and q are refered as the input and output maps respectively. One possible mapping would be to let each factor Nk represent the base of a number system and express the input and output maps in mixed radix notation. This mapping attributed to Cooley and Tukey [17] expresses p and q as P = Pm + NmPm-, + Nm-1Nmpm-2 + Nm-2Nm-lNmPm-3 +... + N2... NmPl q = ql + Niq2 + NN2q3 +... +N1.. Nm-lqm (3.79) Note that the order of the radix expression for the output map is the reverse of the input map. Performing the multiplication pq yields N N N N Pq = N lq + NP2q2 + P3q3... PmQm + nN2puqv + apr,q (3.80) iv! l2 ly3 N^

27 where n is an integer, a is real and u $ v, r $ s. The first m terms of (3.96) represent m Nk point DFT's. The next term may be ignored since WnN2puqv = 1 for all n,p,q. The last term represents the cross product terms which are the'twiddle factors' in [18]. It is possible to avoid these twiddle factors by using a different map construction. An understanding of the prime factor algorithm requires a few introductory terms from number theory. A congruence relation is defined by assuming that a is congruent to b Modulo N if a and b yield the same remainder when divided by N. A congruence is defined by a b (mod N) (3.81) A set integers is called mutually prime if the greatest common divisor between them is 1. The prime factor mapping is based on the following theorem. The Chinese Remainder Theorem [19] states that given a set of mutually prime integers { N, N2, N3,..., N} the system of congruences x = rk (mod N) k=l,...,m (3.82) has a unique solution x (mod N). Two basic maps were suggested by Goode [20] and called the Sino and Ruritanian maps, respectively. Both of these expansions use the same map for input and output. The Sino map may be constructed as follows. Let N N N N p = L1 -pi + L2p + L-p + +L... + L p (mod N) (3.83) N1 N2 N3 Nm N N N N q = L1 -ql + L2 q2 + L3-P3 + Lm (mod N) (3.84) then there is a one to one correspondence between p and {pltP2, P3 *,Pm } and between q and {qlq2q3... qm} if the integers L, L2, L3,...,Lm are chosen such that LkN - 1 (mod Nk) k = 1,...,m (3.85) Nk,

28 Unfortunately, a numerical implementation of the DFT using this map will either require the solution of this set of congruences for each N or require auxiliary storage of a precomputed set. This is somewhat undersirable since speed is of prime importance. An alternative called the Ruritanian map is a special case of the Sino map where Lk = 1 for k = 1,..., m. N N N N P = -P+ P2 + 3P3 +... + Pm (mod N) (3.86) N1 N2 N3 Nm N N N N q = Nql + 2-q2 + -3 +... + - qm (mod N) (3.87) Since W is periodic in N W(U)N(V)N = Wuv and the exponent does not have to be evaluated Module N. Performing the multiplication pq yields N2 N2 N2 N2 Pq = 2Plql + 2JP2q2 + 3P3q3 +... — Pmqm + nN2 (3.88) 1 2 3 2 where n is an integer. IVP may now be factored as N f X N WV = (W1Q) N (W^ ) (3P33)N3.... (Wmqm) Nm (3.89) This factorization is not in the desired form because each W.pkqk is raised to the power -. However, it may be shown [21]-[23] that if N is mutually prime to Nk the effect of the exponent is to permute the output Nk point sequence. Replace the output sequence s(ni) i =0,...k-1 (3.90) by the permuted sequence N s(n,) r = i- (mod Nk) for i = 0,...k- 1 (3.91) Nk It is not apparent which way to sequence through the indices. As an example consider N = 2 * 3 * 5 = 30. The index sequences are given in table 3.2.

29 stage 1 stage 2 stage 3 Pi P2 P3 P P P2 P3 P P P2 P3 P O O O 1 O 0 0 1 O O O 1 1 0 0 16 0 1 0 11 0 0 1 7 0 1 1 17 0 2 0 21 0 0 2 13 I 1 1 2 1 0 1 22 0 0 3 19 0 2 2 3 1 1 1 2 0 25 1 2 2 1 1 2 1 19 1 1 0 26 0 0 3 19 0 0 2 13 1 1 1 2 J 0 3 4 0 1 2 23 1 1 2 8 0 1 4 5 0 2 2 3 1 1 3 14 I 1 1 20 3 4 1 1 4 A2 0 2 0 21 1 1 3 14 0 2 0 21 1..2 0 6 1 2 3 I24 0 2 1 27 O 0 1 7 0 0 4 25 0 2 2 3 0 1 22 0 1 4 5 0 2 3 9 0 1 2 23 0 2 4 1.5A 2 4 15 L1 L 2 1 0 16 1 0 0 16 0 2 3 9 1 1 0 26 1 0 1 22 1 2 3 24 1 2 0 6 1 0 2 28 0 0 4 25 0 0 1 7 1 0 3 4 1 0 J 1 0 0 1 1 17 1 0 4 10 0 1 0 11 0 2 1 27 0 1 0 11 1 1 0 2fi 1 0 2 28 0 1 1 17 0 2 1 27 1 1 2 8 0 1 2 23 iL 2 1 12 1 2 2 18 0 1 3 29 0 0 2 13 0 0 3 19 0 4 5 L 0 2 28 0 1 3 29 1 2 0 6 0 1 3 29 0 2 3 9 1 2 1 12 LLL 3 14 1 0 4 10 1 2 2 18 0 2 4 15 1 1 4 20 1 2 3 24 1 2 4 30 1 2 4 30 1 2 4 30 Table 3.2. Index sequence for N = 2 * 3 * 5.

30 Table 3.3 shows the possible factors for N chosen from the list of {2, 3,4, 5, 7, 8, 9, 16}. Let ak = multiplications for kth factor DFT (3.92) /k = additions for kth factor DFT (3.93) then for N = NN2N3 * N, m N total multiplications = ~- ak (3.94) k=1 Nk total additions = - /1 (3.95) k=1 Nk

31 N Factors Mult. add. N Factors Mult. Add. 2 2 0 4 80 5 16 260 1284 3 3 4 12 84 3 4.7 304 1536 4 4 0 16 90 2 5 9 380 1996 5 5 10 34 105 3 5 - 7 590 2214 6 2 3 8 36 112 7-16 396 2188 7 7 16 72 120 3-.58 460 2076 8 8 4 52 126 2 7-9 568 2780 9 9 20 88 140 4-.57 600 2952 10 2 5 20 88 144 9 16 500 2740 12 3-4 16 96 168 37 -8 692 3492 14 2.7 32 172 180 4.59 760 3704 15 3 -5 50 162 210 2.3 - 57 1180 4848 16 16 20 148 240 3 516 1100 4812 18 2*9 40 212 252 4.7-9 1136 6064 20 4.5 40 216 280 5 7 8 1340 6604 21 3.7 76 300 315 5 7 9 2050 8462 24 3 8 44 252 336 3.7 16 1636 7908 28 4.7 64 400 360 5 8 9 1700 8308 30 2-.35 100 384 420 3 4 57 2360 10536 35 5-7 150 598 504 8 7-9 2524 13388 36 4.9 80 496 560 57- 16 3100 14748 40 5 8 100 532 630 2 5 -79 4100 21964 42 2 37 152 684 720 5-.916 3940 18596 45 5.9 190 746 840 3 5 78 5140 23172 48 3.16 124 636 1008 7. 916 5804 29548 56 7 8 156 940 1260 4- 5-79 8200 38888 60 3.4.5 200 1104 1680 3. 5 7 16 11540 50964 63 7 9 284 1264 2520 5-.7.89 17660 84076 70 2 5 7 300 1588 5040 5 7.9 16 39100 182012 72 8 9 196 1172 Table 3.3. Real multiplications and additions

32 3.5 Extension to Two Dimensions The two dimensional fourier transform and inverse transform are defined as Z(fy,f) = z(x, y)e-2 = z(fxz+fYY)dx, dy (3.96) J-OO J -00 z(x, y) = (f, fy)eJ2r(fz+fYY)dfx, dfy (3.97) -oo00 J-00 Making the change of variables x = u+xo and y = v+yo and defining the quantities s(u,v) = z(u + o, v + yo) (3.98) (ff, fy) = Z(f, fy)ej2r(fX~O+fl~) (3.99) the transform pair may be written as (Af,7 f) = i L| s(u, v)e-j2(fuu+fvv) du,dv (3.100) s(u, v) = s(fu, fv)e'2r(fuu+"V)df dfv. (3.101) -90 -oo Since the exponential kernal of the integral is separable, the two-dimensional DFT may be decomposed into a series of one-dimensional DFT's. Because the data set is padded, a significant savings in computation time may be achieved by avoiding any unnecessary one-dimensional DFT's of sets of zeros. The two-dimemsional DFT is decomposed as shown in figure 3.5.

/ v /v B 33 VN-1 VN-1 --------- A - VM-1 VM_ I I I I I I I * * * * 1.. * 0 * * 0E * 1** * I* I* I ~I~U ~~N _ vo ~L. E I I s vo [I~~~~ ~U - - ~_N I1 N f f Uo UM-1 UN-l f~ f. fu- f-i fv C v D f- |- VN.fv1 ~ ~ ~ E ~ ~ ~ ~ ~ ~ VM-1 ~0~ ~0~ ~0 0 0 0 0 0 0 ~ f 2 0 0 0 f 1O1e *1 1~1~1@1 1 1 1 1 10101010 f V f vI~ ~ ~ ~ ~ ~ ~ ~ I I I fu Yu0 ~ ~ ~ ~ ~ ~ ~~~ /...1-~ _ l f2/1-'~ - Ju fU ff 2 f, f2 f v E VN_ 1 A- Two-Dimensional Data ~__ _ ~ > 1 1 1 - | B- Data after Yu VM-1 VM ~ * * * * C- Data after, * 0 0 0 0 0- ~ ~,~- - - - - 0 0 0 D- Data after''-1 " *~| * | * ~0!!- * * r E-Data after F-1 u* 0 u0 - Uo uM-1 Un -1 Figure 3.5. Two-dimensional DFT Operations.

34 3.6 Test Results The formulation was tested with a function z(x) = z(x) for -a < x < a (3102) 0 otherwise which is expanded as a finite sum of Chebyshev polynomials Tk( (x)) where 4 zr i(x) = E ri )), (3.103) k=O and ((x) =. Test Coefficients 0 0.0 0.0 0 0.0 0.0 1 0.4 0.4 2 0.8 0.2 3 0.5 0.5 4 1.2 0.3 Table 3.4. Test coefficients.

35 2.9 II I I 2.5 real part -0 Imaginary part i.7 0 1.3 0.9 UL 0.5 I-. LU 0.1 ui -0.3 -0.7 -1.1 -1.5 -3.000 -2.250 -1.500 -0.750 0.000 0.750 1.500 2.250 3.000 X Figure 3.6. Test function.

36 100.0 90.0 unmodified FFT -0- pt WONC 80.0 -' 3pt WONC'1 1pt WOFD 70.0 60.0 50.0 LU 40.0 30.0 20.0 10.0 0.0 -10.0 a I a -29 -22 -14 -7 0 8 15 23 30 FREQUENCY Figure 3.7. Percent error.

CHAPTER IV SCATTERING FROM A SINGLE PLANAR PLATE 4.1 Orientation The orientation of the plate is given as in figure 4.1. Hi / /ki / t f W/ / I I9 x Figure 4.1. Plate with plane wave incidence. 37

38 The incident electric field has unity amplitude and may be decomposed into the orthogonal components E' = [(a * 0) + (a. * )]e-j(k"') = [E,,o + Eyoy + Ezoz] h' (4.1) where ki*r = -ko [sin(,) cos(qi)x + sin(Oi) sin(oi)y]. (4.2) Exo = cos(ai) cos(i,) cos(q,) - sin(a,) sin(qi) (4.3) Eyo = cos(ai) cos(,) sin(qi) + sin(a,) cos(5i) (4.4) Ezo = -cos(ai)sin(O,) (4.5) and h' = e-j(i'f). The corresponding magnetic field components are computed from Hi = k, x Ei = [HoxA + Hyy + Hoz]h where,xo = [sin(ai) cos(,) cos(qi) + cos(at) sin(qi)] (4.6) o o = +-[sin(ai) cos(i) sin(qi) - cos(a,) cos(oi)] (4.7) Hzo = - sin(ai)sin(Oi). (4.8) 5o The E-polarization case is defined at a, = 90~ and H-polarization occurs at a, = 0~. 4.2 Integral Equations The scattered fields may be written in terms of vector potentials as EB =- V x F- ~(V x V x A- J) (4.9) H = V x A - k (V x Vx - m) (4.10) where Jem are electric and magnetic polarization currents and A and F satisfy the inhomogeneous vector Helmholtz equations V2A+ko2A = -Je (4.11) 2F + ko2F = -Jm (4.12)

39 which have the solutions in integral form given by A = fJJ Je(P)G(R, R')dv' (4.13) F = JJj Jm(R)G(R,')dv' (4.14) with e-jkolIR-R' G(R, R') l- (4.15) The incident fields induce electric and magnetic polarization currents inside the material body. The total internal fields are expressed in terms of the currents as ET = /3Je and HT = /3Jm where pe 1)ko (4.16) (~t - 1) ko -j /3m (4.17) These currents represent the source of the scattered field and are the solutions of the volume integral equations formed by relating the internal total, scattered and incident fields as T_ Es = Es (4.18) fT_ -H = Hi. (4.19) If the plate is electrically thin, the internal field components are assumed to have a constant z variation such that - = 0. It is convenient for computational purposes to express the integral equations in terms of tangential variables such that the normal components E. and H. are found from (4.9) and (4.10) directly and the tangential components are found after substituting (4.11) and (4.12) into (4.9) and (4.10). The volume integral equations may then be implicitly written as Z o 2 + )A, + A + Fz = Et (4.20) e A [aA + (ko + 2) Ay] - Fz= Et (4.21)

40 a a -z (a2 a2 rpeJ_ - Fx + F -j ) y2 Az = Ez (4.22) 3pmJ + j koZo [ + Fr+ a y- A = H} (4.23) J kozo y z + (ko + y) Fy + Az = Hy (4.24) 1lJz- + A - Ay - ( + 2) Fz = Hz. (4.25) Since the volume current has no z variation, the volume integral may be reduced to a surface integral. The volume currents may then be replaced by an equivalent normalized electric surface current Ke = rJ.eZol and a magnetic surface current Km = rJm. Assuming that the coordinates are normalized with respect to Ao (4.20)-(4.25) reduce to w K + iL [J1K: + L2KV + 6Km] dS'- E= (4.26) w1K' + JjF [I2KC + X3K - x5Kz] dS' = (4.27) w2K - Jj [L 6Km - X5sK + X4Kz] dS' = E (4.28) w3Krm + IJs ['ilK + 22Kym - 6Kz] dS' = H (4.29) w3Km + Fjs K [ + 3Kym + sKze dS' = Hy (4.30) w2Kz + Jfs [6Ke - s5Ky- 4Km] dS' = H (4.31) where q~ 1 2 j~- - 2 \ 4ar 021 02 1 4r2+ z)G T2 J= (4.32) = j1 (2 a2 G 21 (z02 a2\ 42-y2 3 4 G 2 F = 2 y2 e (4.33) 2ir ay2/ 2ir \0X2 ay2/ 0 0 5 =-G 6 = -G (4.34) Oa: Oy and -0 = -3 r efj,(1) 1 j W= -J(er - 1)27rr = 2r + (4.35) Er 11 [4' j[[(e - 1)- + (e) )2] w2 = -j (,. - 1)2irr = 27rr (e- 1)2 + (eU)2 J (4.36) [e-2g('~-l)+(e)]

41 W3 -_ _ =[ ] (4.37) (/it - 1)27rr 27rt [(14 - 1)2 + ( ))2 ) W4 = [r = ~ [14' - (1_ 1) + (1 ]] (4.38) fi(tr - 1)27rr 27rr' (14- 1)2 + ( ")2 The two dimensional fourier transform and inverse transform are defined as (fAfy) = J i| g(z y)e-3j2(f'x+fvY)ddx dy (4.39) -00 -00 g(x, y) = (f, fy)ej2,(f+fY)dff, dfy. (4.40) Joo -00 The two dimensional green's function has a fourier transform given by e-j2V/X2+ -j d for fx2+ f2 <1 e 2,r -3 ~~~~~~(4.41) 47r 2 + y2 d for +f > 2ir for f _2>1 where d = |(lf,2 + fy2 - 11) —. Let X denote the fourier transform and X* signify complex conjugation. The fourier transforms of (4.32)-(4.34) are - = Ci(1 for f2+ fy2<1 (4.42) 0 forf2+f2>1 0 for f2 + f2 < 1 i 0for1 (4.43) 1 forf +fy> 1 (2 for f2 + f2 < 1'2r = (4.44) 0 for f2+ f2 > 1 =2i O for f2 + f2 < 1 (4.45) C2 forf~+f/2>l I3r = 3 (4.46) 0 forf I+fS~>1 [ = 0 for/f2+f/2<1 ~I/^3i f (4.47).~a for f + fy2 < 1 = -— orf+f<1 (4.48) 0 forf,2+f,2>1

42 i 0 for f2+f2<1 /4i = -. (4.49) 1 4 for f + f > 1 0 for (4.50) 0 forf2+f2> 1 _ 0 for f+2 + f2 1 5ij = forf f 1 (4.51) C5 for f + f > 1 T =6r = < 6 forf2+f < 1 (4.52) 0 for fx +f >1 = 0 forf^+f2<1 | O for fx2 + fY2 < 1 6F6i - _ (4.53) 46 forf2 + fy 1 where S 1 = (1-f-2)d (4.54) 2 = -fxfyd (4.55) (3 = (1-fy)d (4.56) =4 = -(f? + f)d (4.57) 5 = fd (4.58) ~6 = fyd (4.59) The adjoint operators are computed from TIa = Tn*(-fx, -fy). 4.3 Surface Generation The surface generation is based on the premise that any planar area may be approximated by a collection of square cells of dimension As. Arbitrary perimeter geometries may be generated by passing a spuare arroy of cell centroid locations through a series of constraints which determine whether the cell is inside or out

43 side the perimeter. This information is then stored in a square tag array which represents a digital code for the geometry. section ofl a p llanaprfct.yi: ii a 11 1111 c p y a ||||I|L |||L i1 i ^ ^i I l i* l x'..1.1,.':..1.1,.' I::.:..LA.A _ _.........:.::... 1:.': _ Figure 4.2. Digital Generation of a Square, Circular, or Triangular The program CGFCON solves for the surface current and the backscatter cross III = _!....!..... 1.....__ the present time the program is set up to allow circular, square and equilateral triangular perimeters. The unknown surface currents are found by solving a pair of coupled integral equations by a combined conjugate gradient FFT method. _ _ _ _ s. 3 3- 3L]3 3L 3 i _..... j J. iele~e 1! l e I _,_. 1 v _ _ _ _..., _ _ _ _ _: _ _ i _ _ _ i,,_ _ _ _ 2V ___________________:I' I |I II of oupedintgrl euaion b a omine cnjuat grdint FTmetlo...,

44 4.4.1 Formulation and Implementation The formulation for E-pol or H-pol results in two coupled integral equations which are defined and valid on the plate. Assume the coordinates are normalized to A0 and the current is normalized to Zo1. The integral equations are fS [TlKx + 2K;] dS Ed (4.60) IffT [2K +;I3K] dS' = E (4.61) The plate is divided up into square cells of side A. The algorithm is given by solve the 2 coupled equations for Ke and Ky initialize, - F,2, + 2,3Kf] (4.62) q _= F-1 [1,2 +'"2,3 r] (4.63) 30 = (Iql2 + I2)- (4.64) P:,y = oqx, (4.65) for k-1,... qk SF-1 i2,3P] (4.66) C (qk l2 +jkl2)'1 (4.67) Ke'k+l Kk + apk (4.68) -.,y — Y K Ik kP~,Y r + rk - k (4.69) — Y xy -- tkqxy k+1 = -1 [ak +l k+l] (4.70) -- [ 1,2,x -- I 2,3 - y 47 Ok (Iqk+12 +I qk+l 12) (4.71) k+1 = k l+ kk+1 (4.72) y:,Y Py +/ k qxy terminate when (|rfc+l|12 -+ I,~l( 2I); (rk I'2 + Irk+ Iy ) < tolerance (4.73) (IE'I2 + IE1); I

45 The Algorithm is programmed by separating all parameters in to their real and imaginary parts. No complex numbers are explicitly used. 4.5 Thin Dielectric Plate The program CGFDIE solves for the equivalent surface current and the backscatter cross section of a planar dielectric plate illuminated by a plane wave. The plate is assumed to be thin such that the thickness r satisfies the relation r << Ap where Ap is the wavelength inside the plate. At the present time the program is set up to allow circular, square and equilateral triangular perimeters. The unknown surface currents for the E-pol case are found by solving a pair of coupled integral equations by a combined conjugate gradient FFT method. The, H-pol case requires the solution of an additional equation. 4.5.1 Formulation The formulation for E-pol two coupled integral equations which are defined and valid on the plate. The H-pol case requires an additional solution for the normal component of current. Assume the coordinates are normalized to A0 and the current is normalized to Zo1. The integral equations are given by WIKe + JIsj [XlK' + TI2KY] dS' = E (4.74) w Ky + J'j [2K; + i3KV] dS' = E (4.75) w2Ke- - ILs 4KdS' = Ez (4.76) 4.5.2 Implementation: E Polarization

46 The plate is divided up into square cells of side A. The algorithm for E-pol is given as initialize rxy Eu~ - wt h'~ _ ~wK"-,-1 X1,2K__ + X2,3 (4.77) r - 1K -,,2 O (4.77) 0* 1 - a 1y'4= Fwi + F [i2F + (4.78) o = (q~12+ jq~12)l (4.79) Py = /oqxy (4.80) for k= 1,... qky = wIps + _F- [i1,2P + I2,3Py] (4.81) c(k = (Iqk12+ Iqk2)-1 (4.82) Ke,k+l = Kek akk (4.83) k+1 = r, - (4.84) rxIY,y (4.84) qk+ = + +- k[=k1 12r + 3k+1 ] (4.85) k = (lqk+1[2 + Iqk+1 2)1 (4.86) px Pxy k+1 (4.87) terminate when (Jrk-+l 12 + Irk+l 12) (21 (Ir k+2 y + 1I < tolerance (4.88) (IEIJ2 + IEyIl2) The Algorithm is programmed by separating all parameters in to their real and imaginary parts. No complex numbers are explicitly used. 4.5.3 Implementation: H Polarization The H-pol case requires compution of the normal component of current in addition to the previous equations. The algorithm for the normal component of

47 current is initialize E - w2K'e, - -1 ['k48~] (4.89) qZ= w2rz +-1 [4ri] (4.90) o0 = Iq ~ -2 (4.91) P = /3oq (4.92) for k= 1,... k = -w k l 2PZ +. T-1[I4] (4.93) ak = Iq -2 (4.94) j^-ek+i- ^ *+a^(4.95) je,fk+l = KeCk + CkPz(495) k+1 k k (4.96) r-kq (4.96) k+1 = *rk+1 + -1 [ k+l (4.97):w2r z (4.97) Pk = I -2 (4.98) _+1 k I (4.99) terminate when Ir,+ll Ir'I < tolerance (4.100) Each Algorithm is programmed by separating all parameters in to their real and imaginary parts. No complex numbers are explicitly used. The more refined algorithm is given as 4.6 Thin Dielectric and Magnetic Plate The program CGFDAM solves for the equivalent surface current and the backscatter cross section of a planar material plate composed of dielectric and magnetic

48 material illuminated by a plane wave. The plate is assumed to be thin such that the thickness r satisfies the relation r << Ap where Ap is the wavelength inside the plate. At the present time the program is set up to allow circular, square and equilateral triangular perimeters. The unknown surface currents are found by solving two sets of coupled integral equations by a combined conjugate gradient FFT method. 4.6.1 Formulation Assume that the coordinates are normalized with respect to Ao and the electric currents are normalized by Zo1'. The E-pol equations are wiK + /Fj [I1K; + X/2Ky + Ii6Kz] dS' = E (4.101) twiK +, ['2K + [, 3Ky - Xl5Km] dS' = E (4.102) 4w4K + JUS ['6K - 5 - - 4Km ] dS' = H= (4.103) W2K.,m + jIs 1K + 2Ky] dS' = H (4.104) W2Km + JsF[2K, + 3K] dS' = (4.105) and the H-pol equations are w2 K + J/s [\'1K + L2Ky' - T6Kz] dS' = H (4.106) w2Kn + Jjs [X2Kxm + X3Ky +,5Ke] dS' = H (4.107) w2Kz - /s [I 6Km -5- I5Kym + I4Kz] dS' Ez (4.108) wlK+J[s, [IKxe+ 2K] dS' E' (4.109) wlKy +J [t2K + 3Ky] dS' = (4.110)

49 4.6.2 Implementation: E-Polarization The plate is divided up into square cells of side A. The algorithm for E-pol is solve the 3 coupled equations for KI, Ky and Kz initialize 1 i ut(-e,l__ ~ —-1 1l,2__xe'. T2,4Ky "-3,5z r,y= Ex,- wy1K - - [, + K + yK (4.111) 11 r' = H_-w 2Kz - -1 [13 sKK_' -- ~Kz (4.112) Qz,=y = _ + I F-1 [i1,2rl + 2,4r~,5z] (4.113) q~ = w r + F-1-' [ 2F@ r-Ir - 3Tr] (4.114) fo, = (.q. +jq+~1 -a Q - (4.115) pl~ = q~ (4.116) (4.117) fork = 1,... qkvy WPy F + 1, + 2,4PY ~ T3,5P] (4.118) ~~~~1k =' k /3,P qz = Wz - JF16PZ] (4.119) W2Pz' 3px -- 5 y z k = (lqi' + Iq + Iql2) (4.120) Kxek+ Ke + Ctkpk, (4.121),y xl;,y Kzmk+l Kmjk + akp (4.122) rIk+ = k -ak (4.123) -, y,z - sq,Z k+1 = k+1 + F-1 ak+1 a'+1 ~ l a zk+1l] qxIV w1rz,y + F [L1,2r + X2,4`l +y',5r (4.124) k+1 = w +l + -1 [k+1 a k+l a -k+1] (4.125) q- 1 2 z Trx - y - =,412 / = (lqk+1l2 + Iqk+ 2 + Iq1k+l12) (4.126) k,+Z = P + k (4.127) px41,Z 23~,YZPz Qx:y,Z

50 terminate when (I+1 12 + [+i12 + I2+r 12)1 \(jr'j + vI |~'I + Iz+ )< tolerance (4.128) (IE+I2 + IElI2 + IH I2)1 solve the 2 coupled equations for K.' and K' initialize rly = t - w3Kml-F-1 [,2K,1 + K, l] (4.129) x x,- x,y [-12 - [- I 2,4 yj O 20 * 1 +.F-1 -a a n_1 qy = 3 r+,y+F [a,2r +- I2,4r] (4.130) /o = (Iq12+ qO2)1 (4.131) Px, = foq,y (4.132) for k- 1,..., W3py + -1 [ 1,2p +2,4Py] (4.133) Ck - (Iqkl+ Iq2)-1 (4.134) tk+l - Kmk + kPxy (4.135) XY x,y ry = r, - k (4.136) — "~Y rx kq~,y k+1 = k -1 [!- a -k+l + (a4.k+137) q= Wr, + 11[Xl,2 X + XF2,4'- ] (4.137) fzk - q[k+l12 Ik+112 /k = (q+12 + Iqk+11 (4.138) Pk1 =- k (4.139) terminate when l+l 12 + I I+112) (r+12 + Ir ) < tolerance (4.140) (IH1I2 + IHyI2) Implimenting this algorithm is most efficient without use of any complex numbers.

51 4.6.3 Implementation: H-Polarization The algorithm for the H-pol case is solve the 3 coupled equations for K', Ky and Kz initialize =,- y F- [1,2K1 + I2,4Ky'1 T i3,5Kz1] (4.141) - x,y = Ei -4Kre + F1 [i3Km,1 - sK l + 1Ke' ] (4.142) 0 *1, + r'-l -]-''' I2,4 IFl3,5Z] (4.143) q~,y,~'-1 tlar a 2 1 = warl - F-1' [ 1r - rl + ]3,r] (4.144) o = (Iq + Iql +I (4.145) p = - /3oo,,, (4.146) Px,y,z (4.146) forc = 1,..., k'-1 [ +1,2. 2,4Py Ic 3,5Pz] (4.147) q W4(P T - +5Pl + 6Pz] (4.148) ak = il2l1+lIq + IJq (4.149) Km,k+l - Kmk + akpk (4.150) ec+ik Kk P+ kp (4.151) rk+l k k r,xyz (4.152) qIk +1 + _ -_l [i —I a.]k+1 ~i'3 5'z. k (4.153) =+ +l - H1 [ +1'i - r V 6+1] (4.154). qk+l -12 + Ik+12 + -_k+1 12_ W4 -k = % q (4.155) terminate when (k+ k+ l 2 k+ l 2)— I + IHI + IE"2) < tolerance (4.157) (IH|+2 + IH -Ii+|2)i 9 I z

52 solve the 2 coupled equations for K, and Ke initialize 1O w) i e,1 l - Y-1 [ 7e,1 i,] r = y t- l _-l [P1,2l +~ 2,4lyJl] (4.158) Pr = E,y (4.161) for k = 1,... qi = W,,p + - [i1l,2r + 2,4Py] (4.162) /3 = (jqk02 + Iqk12)2 (4.163) P,y+l= Keok + (4.164) +f1 = rk- (4.165) q,1 W +l + -1 [a 2 +li +2,4r y+l] (4.166) px, y + -k (lq+12 + Iqky+2) (4.167),y+1 pk + K kP+l (4.168) terminate when +rk +l 2k k,yI y 1 < tolerance (4.169) E8(R) = jkG(R)[R x N(, k) - ZTN(8, s)] (4.170) where N1'"m(0, 5) has the 8 and q- components N'm = cos(9)[cos(~)S'm (0, k) + sin(q )S'm (0, ~)]-sin(0)S^em(0, ~) = sm n-s()Sm(,) -(q+cos( + COS()Sm(0) (4.171) T+ Y

53 and S~^,m( ) = Jj K (ek')dS'. (4.172) If the current is constant over the square cell and the coordinates are normalized to Ao then ge~ (0, ) = F(0, <S) y z ker(p, q)e32-7rsin()[co()x(P)+sin()y(q)] (4.173) p=l q=1 where F(0 I\ 1 sin [rA sin(9) cos(q)] sin [rA sin(9) sin(Q)] F(Or = 7r2 sin(9) cos(o) sin(O) sin(q) ( For the electric current normalized to Zo' the backscatter cross section is given by I2 P12 a = lim47rR 2 - R-+oo J 12 = I= Tr + NgI + INJW - lNo2]. (4.175)

CHAPTER V RESULTS The conjugate gradient FFT algorithms given in chapter V were implemented using the prime factorization described in chapter III. This allowed an efficient solution of all plate shapes and sizes. In addition, the FFT employed was specialized to exploit the sparcity of the data due to the zero padding. An additional increase in the convergence rate was observed combining the higher order integration formulas with the prime factor FFT. The total computational time was reduced by a factor of 3 to 4 over the authors previous programs using a conventional Cooley-Tukey mixed radix FFT. Figures 5.1-5.4 show the magnitude of the induced surface current on a square perfectly conducting plate illuminated by a plane wave at normal and edge-on incidence. Figures 5.5-5.8 compare computed backscatter patterns with experimental mesurements for a perfectly conducting triangular plate and a RAM material square plate. Figures 5.9-5.10 compare computed results with a semi-analytical method by Chu et al [24] for a circular dielectric plate. Figure 5.11 compares the computed backscatter pattern of a square constant 1 ohm resistive plate with a quadratically tapered plate with 1 ohm at the center and 5Zo at the edges. 54

55 Z-AXIS Y-AXIS Figure 5.1. IgK\l:perfectly conducting plate,s = 2Ao, A = 2\, tol =.007, iter = 500,n = 256,0 = 0~, = 0~, scale = 1.0.

56 Z-AXIS I-AXIS IW~~~~~~X-AXIS Figure 5.2. IKli:perfectly conducting plate, = 2Ao, o = toI =.007,iter = 500,n = 256, 8 00, = 0~, scale = 6.1.

57 Z-AXISI Yi' -AXIS YlC X-AXIS Figure 5.3. |IKl:perfectly conducting plate,s = 2Ao, A = 2oA, tol =.007, iter = 500,n = 256,0 = 90~, 4 = 0~,scale = 1.0.,, Xae-A1.0.

58 Z-AXIS XXI ~~~~~~~-AXx1 Figure 5.4. IlKI:perfectly conducting plate,s = 2Ao, = #o, tol =.007, iter = 500,n = 256,0 = 90~, = 0~, scale = 0.8..007,~~~~~~~~~~~~~~~~ ite 0,n=26 0,~=0,sae=08

59 20.00 10.00 0.00 a dB -10.00- - E-Pol Measured 0 E-Pol CG-FFT -20.00 -30.00 -40.00 1II 90.00 67.50 45.00 22.50 0.00 22.50 45.00 67.50 90.00 = -90 Oi(deg) qi = 90 Figure 5.5. Perfectly conducting equilateral triangular plate: y - zplane, s = 2\o, tol =.05, A = ^Ao,n = 128.

60 20.00 10.00 0.00 g ldB -lo oo H0.00-PolMeasured 0| **** H-Pol CG-FFT -20.00 -30.00 - -40.00 1 1 1 1 90.00 67.50 45.00 22.50 0.00 22.50 45.00 67.50 90.00 qi = -90 Oi(deg) Xi = 90 Figure 5.6. Perfectly conducting equilateral triangular plate: y - zplane, s = 2Ao, tol =.05, A = ko, n = 128.

61 20.00 " -- E-Pol Measured 10.00'\ *0 * * E-Pol CG-FFT 0.00 z dB -lo.oo - dB -10.00 -20.00 -30.00 -40.00 I 1 1 0.00 11.25 22.50 33.75 45.00 56.25 67.50 78.75 90.00 Gi(deg) Figure 5.7. Square plate: s = 2Ao,r =.0254Ao, er = 7.4 -jl.ll,pr = 1.4 - jO.672,tol =.001,4A = \ A, n=128.

62 20.00 >\. - H-Pol Measured 10.00 * H-Pol CG-FFT 0.00 - AD dB -o 10.00o 0 -20.00 -30.00 -40.00 0.00 11.25 22.50 33.75 45.00 56.25 67.50 78.75 90.00 Oi(deg) Figure 5.8. Square plate: s = 2Ao,r =.0254Ao, cr = 7.4 -jl.ll,r = 1.4 - jO.672,tol =.001,A = A0Xo, n=128.

63 10.00 \0.00 E-Pol Chu-Weil-Willis [24] 0.00 \ * * E-Pol CG-FFT -10.00 4 dB -20.00 -30.00 -40.00 -50.00 I'''' 0.00 11.25 22.50 33.75 45.00 56.25 67.50 78.75 90.00 9i(deg) Figure 5.9. Circular plate: ro = A0o, =.01AO, r, = 2.0 - jlO.0,tol =.001,A =: Ao, n=128.

64 10.00 0\.00 -- H-Pol Chu-Weil-Willis [24] 0.00 - * H-Pol CG-FFT -10.00 a dB -20.00 -30.00 -40.00 -50.00 L I 0.00 11.25 22.50 33.75 45.00 56.25 67.50 78.75 90.00 Oi(deg) Figure 5.10. Circular plate: ro = Ao,r =.OlAo, Er = 2.0 -jl.O,tol =.001,A = 2 Ao, n=128.

65 23.0 17.3 11.5 5.0f dB 0 -11.4 -17.1 -22.8 -28.5 - constant 1 ohm -34.3 — 0- quadratic 1 ohm at center 5x377 ohms at edges -40.0 0.0 11.3 22.5 33.8 45.0 56.3 67.5 78.8 90.0 9i(deg) Figure 5.11. Square resistive plate:s = 20o,r =.0lAo, tol =.001,A = 25Ao, n=72.

CHAPTER VI FUTURE WORK There are basically three projects under consideration for the near future which are extensions of the ideas presented in the previous chapters. They are 1 scattering by multiple layer planar plates 2 scattering by single and multiple layer non-planar plates 3 synthesis of material tapers for different plate configurations Project 1 is a direct extension of the techniques presented in the previous chapters. Project 2 concerns the development of a method which combines the versatility of a finite element method with the efficiency of a conjugate gradient method to solve scattering by non-planar plates. Project 3 assumes the material distribution of a plate is an unknown and solves the scattering problem for not only currents but also constitutive parameters such as e and p subject to certain constraints on the radiation pattern. 66

BIBLIOGRAPHY 67

68 BIBLIOGRAPHY [1] E. H. Newman and M. R. Schrote, "An open integral formulation for electromagnetic scattering by material plates",IEEE Trans. Antennas Propagat., vol. AP-32,No. 7,July, 1984, pp. 672-678. [2] M. Naor and T. B. A. Senior, "Scattering by resistive plates",Technical report no. 018803-1-T, Radiation Lab., Dept. of EECS, The University of Michigan, Sept., 1981. [3] A. W. Glisson and D. R. Wilton, "Simple and efficient numerical methods for problems of electromagnetic radiation and scattering from surfaces",IEEE Trans. Antennas Propagat., vol. AP-28,No. 5,Sept., 1980, pp. 593-603. [4] S. M. Rao, D. M. Wilton,and A. W. Glisson, "Electromagnetic scattering by surfaces of arbitrary shape",IEEE Trans. Antennas Propagat., vol. AP-30,No. 3,May, 1982, pp. 409-418. [5] T. K. Sarkar,E. Arvas and S. M. Rao, "Application of FFT and the conjugate gradient method for the solution of electromagnetic radiation from electrically large and small conducting bodies",IEEE Trans. Antennas Propagat., vol. AP-34, May, 1986, pp. 635-640. [6] C. G. Christodoulou and J. F. Kauffman, "On the electromagnetic scattering from infinite rectangular grids with finite conductivity ",IEEE Trans. Antennas Propagat., vol. AP-34, Feb., 1986, pp. 144-154. [7] M. R. R. Hnets and E. Steifel, "Method of conjugate gradients for solving linear systems", J. Res. Nat. Bur. Standard., vol. 49, no. 6, Dec., 1952,

69 pp.409-436. [8] M. R. Hestenes, "Conjugate direction methods in optimization", New York, Springer-Verlag,1980. [9] J. W. Daniel, "The conjugate gradient method for linear and nonlinear operator equations",Siam J. Numer. Anal., vol. 4, no. 1, 1967, pp. 10-26. [10] J. W. Golub, and C. F. VanLoan, "Matrix computions", John Hopkins University Press, Baltimore, MD, 1983, pp. 363-379. [11] IEEE, "Special issue on the fast fourier transform and its application to digital filtering and spectral analysis",IEEE Trans. on Audio and Electroacoustics, Vol. AU-15, No.2, June, 1967. [12] M. T. Heideman and C. S. Burrus, "A Bibliography of Fast Transform and Convolution Algorithms II", Technical Report No. 8402, Dept. of Elec. Engin., Rice University, Houston, Texas, Feb. 24, 1984. [13] R. A. Horn, and C. A. Johnson, "Matrix analysis", Cambridge University Press, London, 1985, ch. 4,7. [14] D. C. Champeney, "A Handbook of Fourier Theorems", Camgridge University Press, Cambridge Cb2 lrp, 1987. [15] P. J. Davis and P. Rabinowitz, "Methods of Numerical Integration", 2nd Ed. Academic Press,pp. 51-198, 1984. [16] R. W. Hamming, "Numerical Methods for Scientists and Engineers", 2nd Ed. McGraw-Hill, pp. 277-316., 1973. [17] J. W. Cooley and J. W. Tukey, "An Algorithm for the Machine Calculation of Complex Fouricl Series", Math. of Comput., Vol. 19, pp. 297-301. April 1985.

70 [18] E. O. Brigham, "The Fast Fourier Transform", Prentice-Hall, Englewood Cliffs, NJ, pp. 191-195, 1974. [19] D. T. Long, "Elementary Introduction to Number Theory", 3rd Ed., PrenticeHall, Englewood Cliffs, NJ,pp. 117-123, 1987. [20] I. J. Good, "The Relationship Between Two Fast Fourier Transforms", IEEE Trans. Computers, March 1971, pp. 310-317. [21] C. Temperton, "Implementation of a Self-Sorting In-Place Prime Factor FFT Algorithm", J. Comput. Phys., 58, pp. 283-299, 1985. [22] J. H. Rothweiler, "Implementation of the In-Order Prime Factor Transform for Variable Sizes", IEEE Trans. Acoust., Speech, Signal Processing, vol. ASSP-30, no. 1, Aug. 1982, pp. 105-107. [23] C. S. Burrus and P. W. Eschenbacher, "An In-Place, In-Order Prime Factor FFT Algorithm", IEEE Trans. Acoust., Speech, Signal Processing, vol. ASSP-29, no. 4, Aug. 1981, pp. 806-817. [24] T. M. Willis and II. Weil, "Disc scattering and absorption by an improved computational method", Applied Optics, Vol. 26, No. 18, pp. 3987-3995, 1987.