FAST MEMORY-SAVING HYBRID ALGORITHMS FOR ELECTROMAGNETIC SCATTERING AND RADIATION by Sunil S. Bindiganavale A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy (Electrical Engineering) in The University of Michigan 1997 Doctoral Committee: Professor John L. Volakis, Chair Assistant Professor Emad S. Ebbini Associate Professor Kamal Sarabandi Professor Thomas B.A. Senior Research Scientist Valdis V. Liepa RL-952 = RL-952

Sunil S. Bindiganavale 1997 All Rights Reserved

This thesis is dedicated to my parents, Indumati and Shamanna and my wife, Lynn ii

ACKNOWLEDGEMENTS My first sincere gratitude goes to my advisor, Prof. John Volakis, for all his help in my transformation from a 23 year old with a bachelor's and little research experience to a 28 year old with a wealth of research experience. His enthusiasm and boundless energy for research and teaching is very infectious. I am extremely greatful to the rest of my committee for their valuable feedback on all aspects of my thesis. The brilliant and lucid lectures given by Prof. Senior in EECS 631 (Fall 1992) will forever remain etched in my memory. I had the good fortune of working with a number of extremely talented and helpful people. At the beginning of my graduate program, the technical discussions with Drs. Jeff Collins, Leo Kempel and Dan Ross were invaluable.. The discussions and friendship of an ideal scholar and gentleman, Dr. JongGwvvan Yook, are deeply appreciated. The friendship of Dr. Angelos Alexanian, Dr. Steve "Mr.Cool" Robertson, Dr. Paul Siqueira and Jim Ahne is deeply appreciated. Discussions with Dr. Jian Gong and his assistance with his FEM program are appreciated. I thank all the other past and present members of the JLV group (in alphabetical order) - Lars Andersen (any discussion with him was enriching), Youssry Botros, Arik Brown (for being a great neighbor), Mike Carr, Mark Casicato, Yunus Erdeinli (especially for his generous help with meshing), Dejan Filipovic, Zhifang Li, Stephane Legault (for his knowledgeable opinions on a host of subjects), Mike Nurnberger, Tayfun Ozdemir and Mikhail Smelyanskiy. The friendship of Herte Gebretsadik is some iii

thing which is encountered only once in a lifetime. I thank her for her extremely cheerful and positive demeanor which enriched my life considerably. Jasmeet Judge was another special friend and discussions with her on anything under the sun was always entertaining. I thank the Radiation Lab and EECS administrative staff for their willingness to go out of their way while doing their job. Finally, I wish to thank my family, the people who have most influenced my life. My parents, with their endless sacrifices, ensured that I had all the time and resources for my schooling. They went through endless hardships themselves but made life a bed of roses for me. I could not ask for anything more from them. My wife, Lynn, was the epitome of love, laughter, support, help and understanding whom I could always turn to during the stressful last few months. The two and a half years since I met her were some of the happiest of my life. Her strong work ethic is most admirable. My parents-in-law, Bev and Jim Loomis, are two singularly outstanding people and their love added a new dimension to my life. iv

TABLE OF CONTENTS DEDICATION...................................... ii ACKNOWLEDGEMENTS.............................. iii LIST OF FIGURES.................................. vii LIST OF TABLES....................................xiii LIST OF APPENDICES............................... xiv CHAPTER I. INTRODUCTION.............................. 1 1.1 Objectives and Background........................ 1 1.2 Dissertation Overview.......................... 5 II. Two and three dimensional boundary integral and hybrid formulations...................................... 10 2.1 Solution of the TM integral equation for modeling two dimensional m etallic bodies.............................. 10 2.2 Finite Element - Boundary Integral solution for a two dimensional groove in a ground plane........................ 13 2.2.1 Formulation and discretization for TE incidence....... 15 2.3 Solution of the Electric Field Integral Equation for three dimensional scatterers with resistive boundary conditions.............. 18 2.4 FE-BI formulation for scattering and radiation by three dimensional cavity-backed antennas.......................... 22 2.5 An approach to generate a sparse matrix directly from the fully populated system.................................... 25 2.5.1 M ethod............................ 26 2.5.2 Considerations for 2D and 3D implementations......... 28 2.5.3 Results......................................... 29 2.5.4 Summary.................................. 32 III. The Fast Multipole Method and guidelines for using the method to calculate the RCS of large objects...................... 40 v

3.1 The Fast Multipole Method........................40 3.1.1 Exact FMM........................... 42 3.1.2 Windowed FMM...................... 45 3.1.3 Fast Far Field Algorithm................... 48 3.1.4 Logic Flow........................... 51 3.2 Guidelines for using the method to calculate the RCS of large objects 55 3.2.1 RMS Error as a function of near-group separation distance and sampling rate................................. 56 3.2.2 CPU requirements........................... 58 3.3 Summary........................................ 60 IV. Applications of the Fast Multipole Method............. 74 4.1 Computation of nose radome scattering by employing the fast multipole method to calculate the RCS of large objects.......... 74 4.1.1 Formulation.......................... 75 4.1.2 Results............................. 78 4.1.3 Summary.................................. 80 4.2 Application of fast algorithms to hybrid FE-BI solutions........88 4.2.1 Problem Definition........................... 89 4.2.2 Results.................................. 89 4.2.3 Summary..................................... 95 V. Applications of the Adaptive Integral Method (AIM)........... 97 5.1 Scattering from Plates Containing Small Features Using the Adaptive Integral Method (AIM)...................... 97 5.1.1 AIM for Planar Scatterers...................... 98 5.1.2 Results.................................... 101 5.1.3 Summary........................... 105 5.2 Computation of radiation and scattering from planar cavity-backed antennas with the Adaptive Integral Method (AIM)......... ]..06 5.2.1 Implementation........................ 107 5.2.2 Results.............................. 108 5.2.3 Summary.................................109 VI. Conclusions................................................. 125 6.1 Comparison between FMM and AIM methodologies.......... 125 6.2 Contributions of this dissertation................... 128 6.3 Recommendations for applicability...................... 129 6.4 Future work............................................. 131 6.4.1 Applications.............................. 131 6.4.2 Improved techniques..................... 131 APPENDICES...................................... 134 BIBLIOGRAPHY............................ 164 vi

LIST OF FIGURES Figure 2.1 Discretization of the cylinder surface in the moment method analysis... 11 2.2 Geometry of the groove in a ground plane.................. 13 2.3 Equivalent problem for (a) Region I, and (b) Region II.......... 14 2.4 Local coordinates for the nth edge....................... 19 2.5 Geometry of a cavity-backed annular slot antenna in a ground plane.. 23 2.6 (a) Illustration of a 15A by 2A rectangular cylinder with circular end-caps (b) Impedance matrix (row=150) (c) and (d) Regions of interest..... 33 2.7 E-pol currents on a 18A metallic strip (360 unknowns)........... 34 2.8 H-pol currents and bistatic RCS of a 12A metallic strip........... 35 2.9 H-pol bistatic RCS at normal incidence of different rectangular cylinders of width 2A with circular end-caps (a) L=15A (b) L=30A (c) L=45A... 36 2.10 E-pol currents on a 8A diameter circular cylinder.............. 37 2.11 ago iteratively refined backscatter pattern of a metallic cylinder (radius=0.5A; height=1.5A); the impedance matrix is 36% full................ 37 2.12 ao, backscatter RCS of a 3A long almond (a) Azimuth plane (b) Elevation plane....................................... 38 2.13 oag backscatter RCS pattern of a 7.5A long missile of radius 0.68A (azimuth cut)........................................ 38 2.14 Comparison of our technique and the FAFFA................... 39 3.1 Computation of the boundary integral matrix vector product using exact FM M....................................... 42 vii

3.2 Computation of the boundary integral matrix vector product using exact FMM....................................................... 43 3.3 The Translation operator for different groups on the boundary of a 50A long segment; 750 BI unknowns; 27 groups..................... 47 3.4 Computation of the boundary integral matrix vector product using windowed FMM.......................................... 48 3.5 Computation of the boundary integral matrix vector product using the FA FFA...................................... 49 3.6 Sequence of operations to be performed in the Exact FMM........ 54 3.7 Sequence of operations to be performed in the translation process of the windowed FMM......................................... 62 3.8 Sequence of operations to be performed in the FAFFA.............. 63 3.9 Code indicating the computation of the matrix-vector product in the Exact F M M....................................... 64 3.10 Code indicating the computation of the matrix-vector product in the translation phase of the windowed FMM...................... 65 3.11 Code indicating the computation of the matrix-vector product in the FAFFA 66 3.12 The process of grouping unknowns - two groups are in the near field of each other if the distance between their centers Pli, is less than dmin... 67 3.13 (a) An isosceles triangular metallic cylinder - 17.95A high with a 2.5A base (b) A rectangular metallic cylinder - 25A x 4A................... 67 3.14 Error curves for the FAFFA for a isosceles triangular metallic cylinder - 17.95A high with a 2.5A base......................... 68 3.15 Error curves for the FAFFA for a rectangular pec cylinder - 25A x 4A... 68 3.16 Comparison of bistatic patterns of an isosceles triangular pec cylinder - 17.95A high with a 2.5A base computed using theFAFFA at different sampling rates........................................... 69 3.17 Time curves for the FAFFA for a isosceles triangular pec cylinder - 17.95A high with a 2.5A base............................. 69 3.18 Time curves for the FAFFA for a rectangular pec cylinder - 25A x 4A... 70 viii

3.19 Time - error curves for the FAFFA for a isosceles triangular pec cylinder - 17.95A high with a 2.5A base............................... 70 3.20 Time - error curves for the FAFFA for a rectangular pec cylinder - 25A x 4A 71 3.21 Error and time curves for a smooth body - a circular cylinder of diameter 50A................................. 71 3.22 (a) Triangular cylinder geometry with discretization segment numbering (b) Currents on the triangular cylinder for nose-on incidence............ 72 3.23 (a) Bistatic RCS computed at nose-on incidence for the geometry of Figure 3.2.2 using the moment method and the exact FMM (b) Error curves for the exact FMM.................................. 73 4.1 The process of grouping unknowns - two groups are in the near field of each other if the distance between their centers rll is less than dmin.. 76 4.2 (a) Resistive plate geometry (b) Elevation plane backscatter RCS (~bb polarization) (c) Elevation plane backscatter RCS (08 polarization)...... 81 4.3 (a) Von Karman radome geometry (2D section) - L = 1.1 m, D = 0.6 m (b) Meshed Von Karman radome (c) Bistatic RCS (i = 0~) - 00 polarization (d) Bistatic RCS (- = 0~) - b! polarization................. 82 4.4 Backscatter RCS for a 10A long nose radome of diameter 1A........ 83 4.5 (a) and (b) Nose antenna - radome geometry - 6 is the angle of antenna inclination w.r.t x-z plane (c) aoo backscatter RCS (y-z plane)........ 84 4.6 (a) and (b) Nose antenna - radome geometry - 6 is the angle of antenna inclination w.r.t x-z plane (c) a++ backscatter RCS (y-z plane)....... 85 4.7 (a) and (b) Nose antenna - radome geometry - 6 is the angle of antenna inclination w.r.t x-z plane (c) aso backscatter RCS (x-y plane)........ 86 4.8 (a) and (b) Nose antenna - radome geometry - 6 is the angle of antenna inclination w.r.t x-z plane (c) aoo backscatter RCS (x-y plane)........ 87 4.9 Geometry of the groove recessed in a ground plane............ 90 4.10 Typical form of the FE-BI system matrix arising from the scattering/radiation problem of a groove in a ground plane........................91 4.11 Convergence curves for the hybrid algorithms for the groove of width 25A 93 ix

4.12 Scalability of the hybrid techniques to smaller problems (a) Problem geometry (b) Bistatic patterns (c) Error table................. 95 5.1 The process of transformation from the original MM grid onto the AIM grid 100 5.2 (a) Matrix build operations and (b) Matrix vector product computation in AIM........................................102 5.3 Monostatic RCS for a circular plate of diameter 2A; Comparison of the standard MM & AIM............................... 110 5.4 Monostatic RCS for a circular plate of diameter 2A with three slots computed with standard MM & AIM (a) Geometry (b) Effect of the slots on the RCS (c) 90 polarization backscatter RCS for the plate with slots (d) A4 polarization backscatter RCS for the plate with slots........... 11 5.5 Monostatic RCS for a square plate of side 4A with three holes computed with standard MM & AIM...........................112 5.6 Monostatic RCS for a circular plate of diameter 1A sampled at 0.03 A (smaller than the AIM grid spacing) due to the narrow center ridge.... 113 5.7 (a) Geometry and mesh of the grated plate (b) Geometry and mesh of the "Groove" plate without gratings (c) 00 polarization backscatter RCS computed by AIM and MM (d) +X polarization backscatter RCS computed by AIM and MM............................... 14 5.8 Electric currents (Solution coefficients) on the narrow ridge for the geometry of Figure 5.6................................ 115 5.9 Original discretization and equivalent AIM grids for the geometry of Figure 5.7......................................... 116 5.10 Error introduced by retaining only the self-cell interactions of the moment method..................................... 116 5.11 (a) Flat and (b) Curved plate with equal side lengths and discretization rates, resulting in equal number of unknowns. While the moment method yields equal solution time for both geometries, AIM would accelerate the solution for the geometry in (a) considerably more than that for the geometry in (b)............................. 117 x

5.12 (a) Geometry of a 1A square plate gridded at A/40 and (b) 4A square plate gridded at A/10. While the moment method results in equal solution times for both geometries since they have equal number of unknowns, AIM would accelerate the solution for the geometry in (a) considerably more than that for the geometry in (b) owing to the smaller FFT pad for the geometry in (a)........................................ 117 5.13 Radiation pattern from an annular slot in the b = 0~ elevation plane... 118 5.14 Radiation pattern from an annular slot in the b = 900 elevation plane..119 5.15 Bistatic scattering pattern from an annular slot; Normal incidence in the = 00 plane and observation is in the q = 900 elevation plane..........120 5.16 Input impedance of a very narrow annular slot computed with FE-BI and FE-AIM..................................... 121 5.17 Bistatic RCS at normal incidence (b = 900 plane) from a cavity-backed slot array computed with FE-BI and FE-AIM.................22 5.18 Radiation from a cavity-backed slot array computed with FE-B][ and FEAIM in the q = 900 plane........................... IL23 5.19 (a) Geometry and surface discretization of a cavity-backed patch antenna (b) Monostatic RCS at normal incidence versus frequency - cavity filling has a Er = 2.2- j0.002 and gr = 1....................... 124 A.1 Traditional Dielectric/Ferrite body formulation............... 142 A.2 The rectangular brick element.......................... 142 A.3 Geometry of a dielectric plate illuminated by a plane wave.......... 143 A.4 Backscatter RCS for a 0.2Ao x 0.2Ao x 0.025A0 dielectric/magnetic plate having cr = 7.4 - jl.ll and Lr = 1.4 - j0.672 (Sampling cell dimensions: 0.05Ao x 0.05Ao x 0.025Ao).......................... 44 A.5 Backscatter RCS for a 0.5Ao x 0.5Ao x 0.025A0 dielectric/magnetic plate having er = 7.4 - jl.11 and Lr = 1.4 - j0.672 (Sampling cell dimensions: 0.05Ao x 0.05Ao x 0.025Ao).......................... 145 A.6 Backscatter RCS for a lAo x AXo x 0.025Ao dielectric/magnetic plate having e6 = 7.4 - jl.ll and l/ = 1.4 - j0.672 (Sampling cell dimensions: 0.05Ao x 0.05Ao x 0.025Ao).......................... 146 B.1 (a) Geometry of the groove in an impedance plane (b) Impedance model (c) Equivalent scattering problem....................... 149 xi

B.2 Geometry of the gap at the junction of two semi-infinite dielectric coatings 156 B.3 E-polarization backscatter echowidth for a narrow gap in a grounded dielectric layer as in Figure 2. (a) w=0.1A, d=0.2A and 6, = 2 (b)w=0.2A, d=0.1A and ~r = 10................................ 158 B.4 H-polarization backscatter echowidth for a narrow gap in a grounded dilectric layer as in Figure 2. (a) w=0.1A, d=0.2A and ~r = 2 and (b) w=0.2A, d=0.4A and Er = 2............................... 159 B.5 E-polarization backscatter echowidth at oblique incidence (45~) for an impedance insert of normalized impedance r = 0.25 in various impedance planes as a function of the insert's width (symbols: small width approx, lines: high freq solution)............................ 160 B.6 E-polarization backscatter echowidth at oblique incidence (45~) for various three-part impedance planes as a function of the insert's width (symbols: small width approx, lines: high freq solution)..................1L61 B.7 H-polarization backscatter echowidth at oblique incidence (45~) for various impedance inserts in a metallic plane as a function of the insert's width (symbols: small width approx, lines: high freq solution).............162 B.8 H-polarization backscatter echowidth at normal incidence for various impedance inserts in an impedance plane of normalized impedance r/1 = j0.75 as a function of the insert's width (symbols: small width approx, lines: high freq solution).................................. 163 xii

LIST OF TABLES Table 4.1 CPU Times, RMS error and Storage of the hybrid algorithms...........92 4.2 CPU Times of FE-AA and FE-AIM algorithms...................93 4.3 Comparison of mesh truncation schemes........................93 4.4 Exact BI operation count of the hybrid algorithms. Nb is the number of BI unknowns. NNG is the number of near groups (groups which are treated with the exact moment method). QWIN is the width of the window in the windowed FM M................................. 94 5.1 Solution CPU time and memory requirement of the moment method... 106 5.2 Solution CPU time, memory requirement and RMS error of AIM (all entries in this table were computed with an AIM grid spacing of 0.05A)... 106 xiii

LIST OF APPENDICES Appendix A. Solution of a new reduced-unknown integral equation for modeling inhomo geneous dielectrics................... 135 A.1 Introduction................... A.2 Discretization.................. A.2.1 The Integral equation........ A.3 Results and Validation............. A.4 Conclusion.................... B. Scattering by a Narrow Groove in an Impedance Plane B.1 Formulation................ B.1.1 E-Polarization.......... B.1.2 H-Polarization.......... B.2 Integral Equation Solution........ B.2.1 E-Polarization.......... B.2.2 H-Polarization.......... B.3 Far Zone Scattered Fields and Echowidth B.4 Validation of Results........... B.5 Summary.................. 135 136 136 140 141 147 148 149:51 152 153 154 155 '155 157 xiv

CHAPTER I INTRODUCTION 1.1 Objectives and Background The main theme in this dissertation is to address and recommend techniques for alleviating the memory requirement and execution time associated with implementations of popular frequency domain numerical techniques for electromagnetic scattering and radiation. The excessive computational demands associated with several of these techniques have precluded the simulation of large structures and limited their use as design tools. Of particular interest in computational electromagnetic is the reduction of computational complexity associated with integral equation solution techniques such as the moment method [1]. Incorporation of the radiation condition in the form of the Green's function ensures that integral equation formulations are rigorous. However, the price to be paid is the associated full interaction matrix with a storage requirement of 0(N2) and a solution time which may be as high as 0(N3) for direct matrix factorization, N being the number of discretization unknowns. Integral equations have also proved to be quite accurate for terminating the finite element computational domain [2] but tend to dominate the overall computational complexity, since the finite element method has a low (0(N)) computational requirement. Consequently, any reduction in computational complexity would lead to CPU and 1

2 memory improvements to the hybrid Finite Element - Boundary Integral (FE-BI) methods with their partly full-partly sparse systems. The techniques which can be used to alleviate limitations in integral equation methods can be classified into three categories. The first encompasses techniques yielding approximate sparse matrices for direct solution, by using specially designed basis and testing functions. A sparse matrix can be generated if the field due to a single basis function is restricted to within a certain region over the surface. Such basis functions could take the form of multi-resolution basis such as wavelets [3], [4]. In [3] the electromagnetic coupling through a double-slot aperture in a planar conducting screen is considered by using a wavelets basis in a moment method integral equation solution. The same technique is employed in [4] to model planar dielectric waveguide structures. The above applications yielded highly sparse matrices but are difficult to apply to arbitrary three dimensional problems due to their complexity and uncertainty in the degree of sparsity. Also, certain wavelets like those by Daubechlies [5] destroy the original problem symmetry. Special basis functions referred to as "standing wave" basis functions were successfully used in [6] to achieve a low operation count in integral equation solutions. These basis functions are bi-directional, defined over the surface of the scatterer, and radiate collimated beams in the outward and inward directions relative to the surface. To further reduce interactions within the body, uni-directional basis functions i.e. ones that radiate primarily outward have also been proposed [7]. However, these techniques require a smooth and mostly canonical surface and so far they have not been applied to arbitrary three dimensional bodies. Another technique which relies on clumping testing functions to reduce the number of interactions [8] suffers from the same limitation. The second category, closely related to the first, is based on isolating the dominant

3 elements of the interaction matrix. This is done by either filtering out off-diagonal matrix elements [9], [10] or by defining a set of linearly independent functions over the problem geometry to account for non-self cell terms [11]. In [9] and [10] an iterative refinement procedure is also recommended to improve the initial solution which guarantees convergence under certain conditions. These techniques can be easily implemented for two-dimensional and three-dimensional scatterers but are most accurate for relatively smooth bodies. A common feature of the first two categories is that storage and operation count are still CN2 albeit C is much smaller than that for traditional implementations. The third category of techniques concentrates on speeding up the iterative solution (matrix-vector product) of the moment method linear system. The k-space (CGFFT) technique [12],[13],[14],[15] utilizes Fast Fourier Transforms (FFT) to compute the convolutional kernel arising from a regular grid. The Toeplitz system matrix results in storage of 0(N) and use of the FFT results in 0(N log N) CPU requirements, making it the method of choice for regular grids. Techniques which circumvent this limitation are the Fast Multipole Method (FMM) and the Adaptive Integral Method (AIM) which are amenable to iterative solutions only. A common feature of both AIM and FMM is that they break the interaction matrix into "near" and "far" zone sections based on the distance between testing and source locations. The near zone matrix terms are treated with the exact moment method, thus ensuring full accuracy of the self and near field interactions. Since the iterative solver requires only matrix-vector product computations, this precludes an explicit calculation and storage of the far zone interaction matrix entries resulting in substantial memory savings. Both AIM and FMM are robust, implying that they can be applied to any problem which can be solved with the moment method without any concern for the

4 geometrical shape of the scatterer. The FMM is an efficient computational procedure for calculating all interactions in some "N-body" problem (i.e. a system where each of the N particles interacts with all others). So far the FMM has been applied to a diverse set of applications ranging from the calculation of gravitational interactions between stellar bodies to interactions between biological molecules in systems such as blood hemoglobin. The N-body gravitational problem was solved using the FMM in [16] and such simulations provide a better understanding of the process by which galaxies are formed. The FMM has enabled the simulation of systems with more than 17 million particles [17]. In molecular biology, an adaptation of the fast multipole method referred to as the cell multipole method [18] was used to model power-law forces in molecular systems. The reader is referred to [19] for an excellent overview of the FMM for static fields. An investigation of these applications suggests that the FMM could be used to calculate radiating fields by replacing the particles with electromagnetic current sources. However, due to the oscillatory nature of dynamic fields, the FRMM implementation used for gravitational and molecular computations does not expedite matrix-vector multiplies for dynamic fields. This is because the number of multipoles needed to accurately represent a field depends on the size of the source distribution in wavelengths. The larger the source distribution, the more the number of multipoles needed to approximate the field with sufficient accuracy, independent of the distance between the source and observation points. The first applications of the FMM to dynamic fields are described in [20],[21]. While the algorithm described in [21] is of O(N1 5) complexity, further reduction in operation count can be achieved either by nesting [22] or by recognizing the directional characteristics of interactions between groups of radiators [23]. Guidelines for

5 the practical implementation of the FMM to electromagnetic scattering problems are given in [24]. Application of the FMM to reduce the operation count of the boundary integral in Finite Element - Boundary Integral (FE-BI) solutions is described in 125] and [26]. A comparison of the various FMM algorithms which could be hybridized with the FEM are given in [27] and [28]. AIM [29] is a variation of k-space techniques and transforms the moment method grid to a regular cartesian grid on which matrix-vector products are performed using the FFT and O(N 15) or less computational complexity. Such transformations result in certain interesting characteristics for this algorithm which distinguish it from the FMM. In this dissertation we have focussed our efforts on the second and third category of techniques. Our goals include comparative evaluation of various fast algorithms for electromagnetic simulations, their implementation in hybrid FE-BI systems, error analysis and examination of CPU and memory. An assessment of the effectiveness of these algorithms for electromagnetic simulations is provided at the end of this work. 1.2 Dissertation Overview Chapter 2 of the dissertation presents a review of the integral equation formulations required in the subsequent chapters. We present an algorithm for generating sparse matrices from the original fully populated interaction matrix. This algorithm is based on filtering the dominant elements of the interaction matrix. A feature of the algorithm is its iterative refinement scheme which guarantees convergence to the exact solution provided certain criteria are satisfied. Although yielding high degrees of sparsity, this technique is still 0(N2) and yields more accurate results for relatively smooth scatterers. This indicated the need for a more robust technique capable of

6 handling arbitrary geometries. In Chapter 3 we introduce the FMM and its variants. For the first time an error analysis of the method when applied to electromagnetic scattering problems is performed. Based on given error criteria, we set guidelines for choosing the various parameters affecting the speed and solution accuracy. Such an error analysis aids in the selection of the various parameters of the method in order to achieve a desired compromise between solution speed, memory and accuracy. Chapter 4 describes some applications of the FMM and its variants. We present the application of a version of the FMM to compute the scattering from a three dimensional aircraft nose radome-shaped geometry. The purpose of this application is to characterize the interactions between a nose radome and the antenna located at the base of the radome. Traditional moment method implementations [30] do not permit analysis of realistic size radomes owing to the overwhelming computational requirements. With the application of the FMM we are able to simulate realistic dielectric radome geometries. This chapter also describes the application of the FMM and its variants to a hybrid FE-BI solution for a groove in a ground plane. The latter method is well known for its geometrical adaptability and material generality without compromising accuracy. By using FMM to treat the BI submatrix the method can be used for large scale simulations as well. Each version of the FMM introduces a different approximation to the implementation of the boundary integral. On the basis of comparisons among execution times for given error criteria, we are able to recommend an algorithm which is the best compromise between accuracy and execution time. Since all three algorithm versions are executed on the same computing platform, a realistic interpretation of the relative computational requirements are given for the first time. This is of paramount importance since most of the algorithms

7 belonging to the FMM family have a stated computational requirement of Q(N1 5) or O(N1 33), where N is the number of unknowns. This computational requirement is asymptotic, but for smaller N, the constant multiplying the operation count assumes significant importance. Actual execution of the codes on the same platform gives information on the relative magnitude of the constants. We also examine the effect of the FMM on the hybrid system conditioning and convergence. Chapter 5 discusses the suitability of AIM to model planar scatterers with intricate geometrical details. It will be shown that AIM is extremely accurate while saving a significant amount of memory, especially for bodies requiring high tessellation rates as is the case with antennas which include small geometrical details. AIM achieves its speed up by translating the original problem grid to a new overlaid regular grid (AIM grid). We show that this regular grid can be much coarser than the original discretization (up to a factor of 5 for far field calculations) and this provides for significant speed-ups. Since the grid is uniform, the fast fourier transform can be used to compute the matrix vector products resulting in considerable CPU reduction. AIM is particularly attractive for modeling planar three-dimensional scatterers even on workstation platforms, reducing the need for elaborate parallelizat ion and domain decomposition procedures which are necessary for highly curved scatterers. The capability of AIM to model small and intricate perturbations on otherwise smooth bodies and to predict near fields accurately, even for bodies which are not electrically large, makes it attractive for many scattering and radiation computations. Chapter 5 also describes the application of AIM to a hybrid FE-BI solution for cavity backed slot antennas. Owing to their intricate construction, antennas are not easily modeled by existing algorithms unless they are simple in shape. In contrast to scattering computations, where uniform sampling can be employed everywhere, an

8 antenna geometry may necessitate use of a high discretization rate. We demonstrate this technique for the analysis of cavity-backed slot and patch antennas. Chapter 6 recommends a suitable fast algorithm based on the desired application. It also details the contributions of this dissertation and suggests future areas of investigation. Appendix 1 describes the solution by the moment method of a new, reducedunknown integral equation for scattering by an inhomogeneous dielectric/magnetic structure. The new integral equation involves only the electric field (3 scalar components) as the unknown quantity whereas traditional formulations employed both electric and magnetic fields for a total of 6 scalar components per cell [31]. The discretization of this integral equation is accomplished using rectangular brick elements. Edge-based linear shape functions are used for the expansion of the field inside the dielectric and a modified Galerkin's technique is employed for testing the integral equation. Results demonstrating the validity of the integral equation are presented. Although the integral equation results in reduced unknowns, the computational burden of O(N2) storage and O(N3) execution time makes it an ideal candidate for a fast integral algorithm implementation. While the focus of the previous chapters is on electrically large problems, Appendix 2 describes a simple, accelerated, quasi-analytical solution for an integral equation when the problem is electrically small. Based on the lines of [32]-[33], it precludes construction of the moment method system interaction matrix. We illustrate this technique with the analysis of two-dimensional scattering from a narrow groove in an impedance plane. The groove is represented by an impedance surface and hence the problem reduces to that of scattering from an impedance strip in an otherwise uniform impedance plane. On the basis of this model, appropriate inte

9 gral equations are constructed using a form of the impedance plane Green's functions involving rapidly convergent integrals. The integral equations are solved by introd(ucing a single-basis representation of the equivalent current on the narrow impedance insert. Both transverse electric and transverse magnetic polarizations are treated. The resulting solution is validated by comparison with results from the standard boundary integral method and a high-frequency solution. It is found that the presented solution for narrow impedance inserts can be used in conjunction with the high-frequency solution for the characterization of impedance inserts of any given width.

CHAPTER II Two and three dimensional boundary integral and hybrid formulations In this chapter we will review a few typical integral equation and hybrid formulations for numerical solutions. These solutions will be candidates for speed-up by fast algorithms in subsequent chapters. At the end of the chapter, we suggest a technique to directly generate a sparse matrix from a fully populated moment method matrix. We examine its application to two and three dimensional bodies and conclude by advocating the necessity for fast algorithms based on iterative solvers. 2.1 Solution of the TM integral equation for modeling two dimensional metallic bodies Consider a two-dimensional metallic cylinder with an arbitrary contour C as depicted in Figure 2.1. For plane wave incidence (E' = zJko(xcos5o+~ysinqo)), the integral equation for E-polarization (TM incidence) is constructed by enforcing the boundary condition Ez + Ez = 0, where Ez denotes the z component of the scattered field. We obtain [34] ko Zo / J()H~(2 (k3 IP — ) dl' = jko(xcoso+ysinqo), C C (2.1) 4 where -p = xx + yy and p' = x'x + y'y describe the testing and source points on the contour C. Jz indicates the surface electric current and the kernel represents the two 10

11 A It (Xn-2,Yn-2) n-2 ' ' (n+l)th n segment (Xn+, Yn+1 Figure 2.1: Discretization of the cylinder surface in the moment method analysis dimensional Green's function. We discretize the contour into N segments and denote the points joining the discretized segments by (xn, yn) with n = 0, 1, 2,..., N - 1. The parametric representation for the points on the (n + l)th segment is x = Xn + COS On y = yn + I sin On (2.2) with = tan-) Yn+ - Y (2.3) Xn+1 - Xn and 1 is the distance measured along the nth segment starting from the point (Xn, yn). The current is expanded using standard subdomain pulse basis functions as n=0 2 J,(l)= JP-1 (1 - (2.4) where Aln indicates the segment length. Using (2.4) in (2.1) and by point matching at (xm + -lm cos Om, ym Y+ -Im sin Om) we obtain the linear equations m - oZ0 N-1 in. H(2) Vm k00 Jn j H0(k.oR')dl', m = 0, 1,..., N - 1 (2.5) n=O Jo

12 with [ (m + slmn- (2 R! — XM Osm O m - Xn — I' COS On For m = n the logarithmic singularity of the kernel requires analytical integration. Thus, employing the small argument expansion for the Hankel function yields 2 Al to 2; H(Jdl' = J.H() (tol)d.A + im -- ( (n ( o Ynl) - (2.8) (2.5) is written in matrix form as and [ZR{Ji = {V} (2.9) where koZo | AlH()(LoR, n) m n Zmn = 1 \ 4 A/.[1- (ln( 4)- )] m=n (2.lo) and omn = I Xm+-an + O OM COS OA sin O n V0 d 2 2 + (m-Yn+ -sinAm - 2 sinn) (2.11) The solution of (2.9) by direct methods involves the inversion/factorization fully populated matrix N-rank matrix which could result in an execution time of 0(N3). 2i 2m- -

13 Solution by iterative methods would require a matrix-vector product involving N2 multiplications (operations) for each iteration. For large N, the computational burden can be overwhelming and this issue with techniques to circumvent it are addressed in this dissertation. 2.2 Finite Element - Boundary Integral solution for a two dimensional groove in a ground plane In the previous section, we discussed the solution of an integral equation for modeling two dimension metallic scatterers. In this section, we review the integral equation termination of a finite element computational domain referred to as the FEBI method. This hybrid method combines the rigor of the integral equation mesh termination with the adaptability of the FEM. We shall consider its application to scattering from a material-filled groove of width w and depth din a ground plane (geometry is depicted in Figure 2.2. The relative permittivity and permeability of the Y Region I (, go) Ground plane r d ^(n pec x w Region II Figure 2.2: Geometry of the groove in a ground plane material filling the groove will be denoted by cr and jT, respectively. The upper halfspace (y > 0) is denoted as region I, and the cross section of the groove (0 > y > — d) is denoted as region II. To take advantage of the FEM's geometrical adaptability and low O(N) storage requirement, it is necessary to decouple regions 1 and 2. This can

14 be accomplished by extending the ground plane over the aperture. The aperture's scattering is then modeled by introducing equivalent magnetic currents M and -M above and below the original location of the aperture as illustrated in Fig 2.3. In this j M. Ml Region I ~o Hy=O ground plane (a) r r) Region II 1 (b) Figure 2.3: Equivalent problem for (a) Region I, and (b) Region II manner the continuity of the tangential electric field is satisfied a priori. It is still necessary to ensure tangential magnetic field continuity and when this is done we generate an integral equation for the unknown equivalent magnetic currents. That is, we must enforce the equality HI(M, J M)Io = Hf (-M)Iy=O (2.12) where HI and HII are the tangential magnetic fields in regions 1 and 2, respectively and J', M' are impressed sources. In our case, HI can be expressed as an integral of M using the ground plane Green's function whereas HI" is formulated via the finite element method. To begin with, the variational equation is employed in region II SF= 0 (2.13) with F =Jj (V x H") (V x H) - k2rHI, H I. dxdy Q Cr Ot r

15 -2jkoYo Mi- Hldx (2.14) where Q denotes the cross-sectional area of region II and r is the line segment specifying the aperture. The introduction of (2.13) eliminates a need to find the Green's function associated with the internal structure of the groove and hence permits the modeling of grooves of arbitrary cross sectional geometry and material filling. The boundary condition (2.12) in conjunction with (2.13) and (2.1)4) imply a system of equations for the solution of the magnetic current M. The next section describes the formulation and discretization of the boundary integral for TE incidence and its eventual combination with the FEM system. 2.2.1 Formulation and discretization for TE incidence The groove is illuminated by the H-polarized plane wave Hi = ejko( sin o +y cos 0o) (2. 1) where ko = 27r/Ao is the free space wavenumber and oQ is the angle of incidence measured from the normal direction. Since the impressed magnetic field Hi has only a z-component, the scattered magnetic field will also be z-directed, and consequently the equivalent magnetic current M may be written as M = ZMz(x) (2.16) From Fig. 2.3(a), the magnetic field in region I due to Mz is given by HI(r) = Hic(x, y) +Hzefl(x, y) _- 2oY Mz(x')HO2)(klor- x'l)dx' (2.17) The functional of (2.14) reduces to F //{- [(F H) (H )2] -k2tr (HzI)2} dxd L. _'I d _f

16 - 2jkoYo MH'dx (2. L8) (2.17) and (2.18) need to be discretized to facilitate the coupling of fields and enforcement of boundary conditions. For discretization, the cross sectional area Q is subdivided into M small triangular or rectangular elements. Also, the line segment r is broken into L segments. With a linear filed distribution, within each element, the field within the eth element having n nodes is expanded as n Hz = E N[(x y)) = {Ne}T{* } = {*e}x{Ne} (2.19) i=l where Ni(x,y) are the shape functions and qe represent the nodal fields. Using (2.19) in (2.18) and employing a pulse basis expansion for the magnetic current kM M L F = E{be}T[KIe]{(e} - j (' s + qs2) 5 As (2.20) e=l s=l where {(fe} represents the nodal magnetic fields within the eth element, (fs1 and ()s2 denote the magnetic fields at the ends of the sth segment and Os = k(YoMz denotes the magnetic current on the sth segment. The matrix [Ke] is given by f I1 N6 ON T (ONe 1 [K'] = ] { [I x x + { { ja - kI,{ Ne}{Ne }Tdxdy (2.21) Minimization of (2.20) with respect to nodal fields yields [K]{~} + [B]{} = 0 (2.22) where {4} is a column vector representing the magnetic field at the N nodes within region II and on F while {b} = koYo{Mz} is a column vector with L elements corresponding to the number of segments on F. The sparse matrices [K] and [B] are assembled by M L [K] = ][Ke], [B] = E[Bs] (2.23) e=l s=l

17 where the elements of [Bs] are given by B8 = Jo (2.24) 2 Discretization of (2.17) by employing Galerkin's procedure yields a second set of linear equations [C]{p} = {if} - [p]{4} (2.25) where {fr} is a column vector representing the nodal magnetic fields on F. The sparse matrix [C] is assembled from [C] = S[C8] (2.26) s=1 where [C8] is a row vector whose elements are given by C = -J (2.27) 2 The excitation column vector {qCinc} is given by OsfC = —2Hznc(xs0)6 m = 1,2,...,L1 (2.28) where xs denotes the mid-point of the sth segment on F. Also [P] is a dense matrix of order L with elements given by Pst = [1- - log(0.1638ko8)J 5s s= t (2.29) Pst = -j o )(koxs - Xt\)6S6t s t (2.30) 2 0 The final system is obtained by enforcing (2.12) on (2.22) and (2.25) to obtain [K] [B] (2.31) [B T] [p { ic }

18 It should be noted that the term 6, (segment length) is retained in (2.27),(2.28),(2.29) and (refeq:ch2Pst) in order to conveniently couple (2.22) and (2.25). The subsystems [K] and [B] are sparse but [P] is fully populated and the product [P]{b} proves to be the computational bottleneck for large hybrid systems especially when the size of the boundary integral is large. This can be remedied by the application of fast algorithms for integral equations which will be discussed in the subsequent chapters. 2.3 Solution of the Electric Field Integral Equation for three dimensional scatterers with resistive boundary conditions For this application, the thin dielectric layer is modeled using the resistive boundary condition [35] n x (Ei + ES) = 7oRJ (2.32) where R is the resistivity, J is the surface current, and E' is the incident field which is a plane wave of unit amplitude given by E = (Ocosa + isina) cjk(xsinOis+y sinisini) (2.33) where ko is the free space wavenumber, a is the polarization angle and (0i, Xi) indicate the direction of incidence. The scattered field Es can be determined from J according to E5 =-jwA- Vq (2.34) where the magnetic vector potential A is given by A e -jkoR A(r) SJ(r) R dS' (2.35) with S being the surface of the body. The scalar potential q is given by (r) - = 4- c (r') dS'. (2.36) 4ir sR dS '

19 where R is the distance between observation and source points, viz. R = -r = /(x-x')2 + (y- y')2 + (z - z')2 (2.:37) The continuity equation is used to relate the surface charge density and the current Vs. J = -jr (2.:38) Enforcing (2.32) on S yields the electric field integral equation for J Etan = (jwA + V)tan + oR J r E S (2.39) To model the current, the scatterer is discretized into triangular patches. The current is then expanded in terms of vector basis functions [36] which are especially suited for triangular domains. Each basis function is associated with an interior (nonboundary) edge, and is nonzero only on the two triangles sharing that edge. Figure 2.4 shows the n" interior edge shared by triangles T: and Tny of area A+ and A- respectively. A point in the triangle pair can be designated by either the global position vector r, th n edge A'-% r+ u Figure 2.4: Local coordinates for the nth edge or local position vectors p = r - r'. The basis function f (r) for the nth edge is

20 defined as n+ + r in T+ fn(r) = 2A- PJ, r in T~n (2.40) 0, otherwise The current J on S is approximated by N J Inf (r) (2. 41) n=1 where N is the number of interior edges and the unknown coefficients In represents the current density flowing across the nth edge of the mesh shared by the T+ and Tntriangles. To solve for the basis coefficients, Galerkin's technique is applied to (2.39) giving jEi. fmdS= jJA fmdS -Jj Vs fm dS + 1o RJ -fmdSm= 1,...,N. (2.42) Using (2.40) in (2.42) yields the N x N system of linear equations, V = ZI where In is the Nth basis coefficient, Zmn is the impedance matrix whose elements are computed from T~oln \J T i T[f 1 - _ ejkR mn = AA Pm(r) ' (r) dS'dS A ~/ ff A~ m _ 2w2JJT AtA Rd' 2 7 JJTm i J + JJT ApA Pm(r) Pf(r)dS (243) where cm and Cn are the positive current reference signs for edges m and n, defined as +1 r in Tm+ cm = (2.414) -1 r in Tm

21 and +1 r' in T+ (2.45) -1 r' in T,The elements of the interaction matrix can be computed directly from (2.43). However, a more convenient way of evaluating these elements is to consider a pair of faces and compute all nine interactions between edges contained by this pair. This enables the loops for assembly of the matrix elements to be over faces, instead of edges, thus speeding up the assembly process. For an observation face p paired with a source face q, the quantity Zpq is computed for all mn edge pairs as 4 AA 2 IP p f r) (r) e dS'dS2r2 = J'eTJ e m dS'dS - R' RJ ' d + Rq JJ p +(r)p+(r)dS} (2.46) The positive current reference signs, Ce and Ce, are now assigned according to f +1, if TP is T+ C = (2.47) -1, if TP is Tm and +1, if Tg is T+ n - (2.48) -1, if Tq is Tn The integrals in (2.46) are evaluated for near and self cells by the techniques detailed in [37]. It should be noted that in (2.46) Tp = Tp+ + T- and Tq = T~ + T;~, thus computation of ZPq involves summation over four triangles. The elements of the excitation vector are given by Vm i JT / r) O (O cos a + q sin a) 2 m c s jko(xsin9i cos,+ysin i sin',i) dS. ( O dS.(2.49)

22 The N x N linear system can be solved either by direct methods such as matrix factorization (which would mean an execution time of O(N3)) or iterative methods involving an operation count of O(N2)/iteration. 2.4 FE-BI formulation for scattering and radiation by three dimensional cavity-backed antennas In this section, we review a finite element - boundary integral formulation for analyzing three dimensional cavity-backed antennas. The finite element discretization is in the form of triangular prisms. Such prisms are the element of choice for modeling planar antennas with fine detail (as small as 50th or 100th of a wavelength) as they require only surface discretization information. In contrast to tetrahedral elements [38], this eliminates the need to generate volume meshes which could be tedious and also removes the possibility of ill-conditioned systems due to degraded mesh quality. In general, for modeling planar configurations the prism element also requires lesser number of unknowns than tetrahedral elements. However, very small details and consequently dense meshes can still lead to boundary integrals with extremely large computational requirement (as described in Section 2.2.1). A technique which reduces the computational requirement considerably is examined in Chapter V. In this section, we present the key elements of a three dimensional finite element - boundary integral formulation with emphasis on the boundary integral computation. For details on the prism element the reader is referred to [39]. Consider a cavity-backed antenna recessed in a ground plane as depicted in Figure 2.5. To solve for the E-field inside and on the aperture of the cavity, it is necessary to extremize a functional, which for radiation and scattering problems may be gen

23 Annular Slot Ground plane Feed S Figure 2.5: Geometry of a cavity-backed annular slot antenna in a ground plane eralized as F(E) = J {(v X (V x E) - k2E.E}dV + JJJ E (jkoZoJi + V x r-1 Mi) dV + jkoZo fL E (H x n,)dS (2.50) o+Sf where Or and /l denote the relative tensor permittivity and permeability of the cavity filling, So represents the non-metallic portions of the aperture and Sf denotes the junction opening to the feeding structures. The volume Vs refers to the volume occupied by the impressed sources J, and Mi. Also, H denotes the magnetic field on So or Sf and n is the outer normal to these surfaces. For a unique solution of E we require knowledge of H over So and Sf. In the case of Sf, H is determined by the feed excitation while that over the non-metallic portions of the aperture is determined by the boundary integral equation H = Hg~ + 2jkoYo G(r, r') ( x E(r')) dS' (2.51) where G is the electric dyadic Green's function of the first kind such that n x G == 0 is satisfied on the metallic platform. For the cavity recessed in a ground plane, G becomes the half space dyadic Green's function G= + kVV (2.52) k 2 i 7r

24 with R = Ir - r'I and I is the unit dyad. For this problem, Hg~ is equal to the sum of the incident and reflected fields for scattering computations and zero for antenna analysis. To discretize (2.50) the volume region is subdivided using prismatic elements. The field in each prism is approximated using a linear edge-based expansion as 9 E = EV= [V] {E} (2.53) j=1 where [V]e = [{Vj}, {Vy}, {V,}] and {Ee} = {E E2,..., E}T. On the aperture, since the top and bottom faces of the prism are triangles, we have a corresponding representation for the aperture fields as 3 Es(r) = E'S (r) = [S]{ES} (2.54) i=l where [S], = [S, Sy]. To generate a linear system for the solution of E, (2.53) and (2.54) are substituted into (2.50). Subsequent minimization of the functional yields (vFvA Nv Ns Nv Ns OEe - = [A ]{E=} + L[BS]{ES} + H{Ke} + {LS} = 0 (2.55) e e=l s=1 e=l s=l where Nv and Ns indicate the number of volume and surface elements, respectively. The matrix elements are given by Aei =J/J {(V x Vi) ' r-l (V x Vj) -koVi ~ 4' Vj}dV (2.56) K = f Vi, [jkoZoJi + V x Ir-1 V x Mi] dV (2.57) B -Jj J 2kSs(r) x Ss(r')Go(r, r')dSdS' + S2ffJ [V x S-(r)][V' x Sj(r')],Go(r, r')dSdS' (2.58)

25 L = 2jkoZo 1 S -. (Hn x )dS (2.59) The boundary integral equation in (2.58) is discretized using basis functions (lefined on the top face of the prism as Si = 2A- x (r- ri) (2.60) 2Ae similar to the function defined in (2.40). Substitution of (2.60) into (2.58) gives the discretized boundary integral which is treated using the procedure outlined in Section 2.3. 2.5 An approach to generate a sparse matrix directly from the fully populated system We suggest a technique for generating sparse matrices directly from a fully populated moment method system [40]. The approach is based on the observation that the interaction between elements lying electrically far away from each other is small and if they are sufficiently far away from the self element, they can be zeroed out as a first approximation. Thus, the full matrix reduces to one that is sparse. Of course, the number of elements to be zeroed out depends on the accuracy of the desired results. This technique should work even better for large scale scattering problems since the percentage of sparsity increases as a function of the physical size of the body. We have examined the idea both for two dimensional (2D) as well as three dimensional (3D) geometries. However, the magnitude of the 2D Green's function tapers off as 1/X/r, whereas the 3D Green's function tapers off as 1/r. Consequently the far field contributions are relatively more important for two dimensional problems and in this case the inclusion of the physical optics currents as a first approximation was

26 found useful. An important component of the technique is the capability to improve upon the approximation in an iterative manner without using additional storage resources. One or two iterations beyond the first yields a good approximation of the correct result provided certain sparsity criteria are satisfied. 2.5.1 Method Consider the linear system [Z{J} = {V} (2.61) where {J} represents the unknown current vector of length N, [Z] is the N x N impedance matrix and {V} is the excitation vector. We are interested in generating a sparse matrix approximation of the original matrix [Z] for the purpose of speeding up and reducing the memory requirements for solving (2.61). We would like to write the impedance matrix [Z] as [Z] = [Z]{F}+ [Z]{1 - F} (2.62) [zS] [ZR] where F is a filter function with a finite flat-top and either a gradual (typically gaussian) or abrupt taper. Employing an abrupt taper is equivalent to applying a tolerancing criterion. Specifically, an element in the ith row of the impedance matrix is discarded if IZ(, j) I < {min(lZ(i, j), Vj) + T.F Z(i,)|} (2.63) where Z(i,j) is an element belonging to the h row and th column of the e impedance matrix and T.F. is a threshold factor. As a first iteration, an approximate solution of {J} is obtained by solving the sparse system [ZS]{J()} = {V} (2.64)

27 where [Zs] is the sparse matrix obtained from [Z] after applying the thresholding criterion (2.63). The iterative refinement technique, performs subsequent iterations in accordance with [ZS]{J(m)} = {V} - [ZRI{J(-)} (2.65) where [ZR] - [Z] - [Zs] and {J(m)} denotes the current estimate after the mrth iteration. The convergence of this iteration process can be investigated by noting that for m = 2, (2.65) can be rewritten as {j(2)} = [ZS]- ([I] [ZR][ZS]-1) {V} (2.66) Then, after a subsequent iteration {J()} = [Z]-1 ([I] - [ZR][ZS]-1 + ([ZR][ZS]-1)) {V} (2.67) and at the mth iteration jm)} = [ZS]{[I] - [ZR][ZS] 1 — + ([ZR][ZS]-)2 ([ZR]- )3 +.([ZR][ZS] -l)m }{V} (2.68) This geometric series will converge if every eigenvalue of [ZR][ZS]-1 has a magnitude less than one [41]. The result of this convergence after an infinite number of iterations would be {j(oo)} - [ZS]- ([I] + [ZR][ZS]-1) {V} (2.69) which can be written using the definition of the residual matrix as {j(oo)} = [zS]-1 ([I] + ([Z] - [ZS]) [[ZS]-1)-1 {V} (2.70) which further reduces to {J( oo)} = [ZS]-1 ([Z][ZS]-1)1 {V} (2.71) [Z]-1 { V} (2.72)

28 hence it follows that {J(~)} is the exact result from the original moment method problem. 2.5.2 Considerations for 2D and 3D implementations The governing integral equations for computing two dimensional scattering from metallic bodies are the well-known magnetic field integral equation (MFIE) and electric field integral equation (EFIE). Expansion of the current using pulse basis functions and point matching, was employed in the moment method solution of these integral equations as discussed in section 2.1. Application of this technique to strips is considered in accordance with (2.64) and (2.65). The second class of scatterers considered were two-dimensional. Owing to the operation of the derivatives on the Green's function, the discretization of the integral equation for H-polarization (TE incidence) results in a more diagonally dominant matrix than that for E-polarization (TM incidence). A single row of the impedance matrix for the cylinder whose cross-section is depicted in Figure 2.6(a) is shown in Figure 2.6(b),(c) and (d). This row, being representative of the matrix character encountered in pulse-basis point-matching solutions of 2D scattering problems, hints that any attempt to make full matrices sparse should address both polarizations separately. Owing to the strongly diagonal nature of the H-pol impedance matrix, the H-pol currents can be found by directly applying a tolerancing criterion. For the E-polarization, the physical optics currents were employed as an initial estimate and the iterative refinement technique was then used to improve on this estimate. The governing integral equation for solving 3D scattering problems is given in [36] and was discussed in section 2.3. The approach combines the advantages of triangular patch modeling and the EFIE formulation resulting in a simple and efficient

29 algorithm. The electric current is discretized in the form of special basis functions associated with each nonboundary edge of the triangulation. Galerkin's method is then used as the testing procedure to obtain the final system of linear equations. A sparse matrix is again generated by applying a magnitude tolerancing criterion to the matrix elements. It is important to note that this criterion is not suitable for generating the sparse matrix because it does not require any special rearrangement of the matrix elements before its application. That is, there is no requirement for the matrix to be diagonally dominant for a successful application of the method. 2.5.3 Results Figure 2.7 shows the E-pol currents on a 18A metallic strip after application of the following filters on the impedance matrix - (a) rectangular filter with flat halfwidth = 40 (b) Gaussian tapered filter with flat half-width = 10, a = 0.005. At normal incidence, employing a gaussian-tapered filter with a decay parameter (a) of 0.005 yields very accurate currents at the first iteration, and considerably more accurate than that obtained by the application of a rectangular filter with the same non-zero width. Even at oblique incidence, solution of the sparse system obtained by employing the Gaussian filter gives very accurate results. H-polarization currents on a 12A strip are shown in Figure 2.8. It is seen that the iteration process converges to yield almost the exact MoM currents by the fifth iteration while the bistatic RCS for normal incidence is accurate after the third iteration. The normal incidence bistatic RCS patterns for 2D rectangular cylinders with circular end-caps computed with different degrees of matrix sparsity are shown in Figure 2.9. Three different cylinder lengths (15A, 30A and 45A) were considered. The results for the 15A long cylinder were obtained by setting T.F = 0.0035 and this gave a 23% full matrix;

30 setting the T.F to 0.005 generates a 19% full matrix. This computation indicates that the maximum error in RCS prediction was less than 3dB with the 19% full matrix and the error was at angles close to grazing. Maximum error for the 23% full matrix was less (2dB or so) at angles close to grazing. We note that the degree of sparsity can be increased for physically large bodies and this is illustrated in Figures 2.9(b) and 2.9(c) where the cylinder length was increased to 30A and 45A, respectively. With the threshold factor set to 0.005, the sparse matrix is 11% full for the 30A body and 7% full for the 45A body. Even with these high sparsities, the predicted RCS is again computed very accurately everywhere except near grazing where a maximum error of 3dB is observed. We next consider the RCS analysis of circular cylinders. The H-polarization analysis of this geometry resulted in the same experiences as in the previous case. However, while computing E-polarization scattering, use of the PO currents as an initial estimate was found useful (the same estimate was used for H-polarization but it did not lead to appreciable improvement in the convergence of the iteration process). That is, for E-polarization we set J(1) equal to the physical optics current and then proceed with the computation of J(2) as described in (2.65). It should be noted, however, that the contribution of [ZR]{J(1)} can be obtained much faster by direct evaluation of the radiation integral using larger elements and without a need to store the [ZR] matrix. Figure 2.10 shows the subsequent PO and iteration currents for an 8A diameter circular cylinder, plotted as a function of the segment number on the cylinder. It is seen that the second iteration predicts currents which are near identical to the exact. In the case of E-polarization, it is again seen that the second iteration is almost coincident with the exact solution. This appreciable improvement is attributed to that the residual matrix [ZR] contains entries which are larger than

31 those in the corresponding H-polarization matrix. A typical 3D test body is a PEC cylinder and the iteratively refined backscatter pattern of a metallic cylinder of radius 0.5A and height 1.5A is shown in Figure 2.11. The matrix is 36% full and we observe that the iteration procedure converges to the reference solution. Calculations for other scatterers are shown in Figures 2.12 and 2.13. Azimuth and elevation backscatter cuts for a 3A long metallic almond are shown in Figure 2.12 and depict very good agreement with 34% full matrices. T'he solution starts to degrade at near grazing incidence, an expected phenomenon, since the traveling wave information is lost when the small magnitude matrix elements are zeroed out for achieving sparsity. A 7.5A long missile of radius 0.68A, with two fins, was also analyzed and a 23% full matrix yields a reasonably accurate solution within the first two iterations as shown in Figure 2.13. These results are quite promising and demonstrate that lower matrix sparsities can be achieved when considering larger 3D structures. In order to estimate the computational requirement of this technique, execution time and error comparison was made with a variation of the fast multipole method (FMM) known as the fast far field approximation (FAFFA) [42]. The test body was a 50A metallic strip excited at an incident angle of 60 degrees from grazing, as shown in Figure 2.14. For the FAFFA the 750 unknowns (N) were grouped into 27 groups (vN) and a near group radius of 2A was employed. For our technique we employed a gaussian-tapered filter with a decay parameter (a) of 0.005 and a flat half-width of 10. Solution times on a HP 9000/750 using a conjugate- gradient solver were about 140 sec for the unreduced matrix, 25.56 sec for the FAFFA and 21.42 sec for our technique. The average error for the FAFFA was 0.53 dB while our technique yielded an error of 0.64 dB when compared to the MoM solution using the unreduced

32 matrix. 2.5.4 Summary A simple technique for reducing the memory and computational requirements of the moment method was introduced. The solution can be improved by successive iterations which account for the far-zone interactions of the moment matrix elements. This iterative refinement algorithm is an attractive feature and was shown, using simple matrix concepts, that the iterative procedure converges to the exact value in the limit, under certain conditions. However, this algorithm works best for relatively smooth bodies and at angles of incidence close to normal on large surfaces. TFhe interactions discarded at the first iteration tend to be crucial for complex bodies, resulting in divergence of (2.68). For such bodies, the levels of sparsity achieved is very small. This lead to the investigation of more robust techniques capable of treating complex arbitrary geometries, such as techniques based on iterative solvers.

33 Segment # 50 Segment # 1, 1.5 E 1 r x E a 1~2 4 15 (a) Segment # 275 ~2kX,! Segment # 325 I I I I - H-pol -- E-pol I I — ---- - - - -- - - - - - - - 100 200 300 Segment number n (b) 400 500 600 t1 E 4 c0.0.5 a 1V 1J5 152 15 154 1 400 410 420 430 440 450 460 mt number n Segmen number n (c) (d) Figure 2.6: (a) Illustration of a 15A by 2A rectangular cylinder with circular end-caps (b) Impedance matrix (row=150) (c) and (d) Regions of interest

34.08.06.04 N.a 02.00 -.02 -.04 5.5.................................................................. x X | || X w X I F I i i II I I I; I IE ZI. I..... I..... I........I I Original (MoM) matrix - --- Filter Function - a=0.005 A Filtered Matrix I I I I I I i i 'l'l '...'' I.'....l '.......I I....I.. I... I.......... 0 30 60 90 120 150 180 210 240 270 300 330 360 n (matrix column) I I 4.5 Q 'a 3.5 2.5 1.5.5 -11 9 MoM (exact) / J - rectangularfilter.^.... v b 1 - gaussian filter- a = 0.005 J - gaussian filter- a = 0.005 g=090~ MM — (exact) ) = 30~ ~ J1 - gaussian filter - a= 0.005....I-.._..- - ( D -5 0 Position x (X) 5 10 Figure 2.7: E-pol currents on a 18A metallic strip (360 unknowns)

35 5 ---- M(exat) 4_ J- gausiaffilter-a= 0.005 X - J3-g~ssianfflter-a=O.OQ5..- f -g assianfilter- a=0.005 o...I.........I...I.... i. = -7 -5 -3 -1 1 3 5 7 Position x (X) 30... I.... 0 MoM (exact) 2 0 - _._ Iter#l. -—.. Iter#2 6 10 5 - Iter#3 0 I' / -20 %, it -30 ' I I 0 15 30 45 60 75 90 Angle 1 (deg.) Figure 2.8: H-pol currents and bistatic RCS of a 12A metallic strip

36 y A iL 00, I42X L 40 35 -30 m I 25 ' 20 15 10 C 5 0 -5 0 15 30 45 60 Observation angle (deg.) (a) 75 90 40 35 i- 30,< 25 - c 20 | 15 0 -5 I I - 40 I I ~ ~ ' ' I I I I0 x 35 — Exact 30 7% Full 30 25 20 I5 15 1| 10 75 90 0 15 30 45 60 Observation angle (deg.) (c) 0 15 30 45 60 Observation angle (deg.) (b) 75 90 Figure 2.9: H-pol bistatic RCS at normal incidence of different rectangular cylinders of width 2A with circular end-caps (a) L=15A (b) L=30A (c) L=45A

37 t Segment # 150 S\0, I Segment # 300... b, I Segment # 1 Segment # 450 2.50 2.00 -o ' 1.50 ct rI 1.00 3.50.00 0 100 200 300 400 Segment number 500 600 Figure 2.10: E-pol currents on a 8A diameter circular cylinder 15 10 m0 5 c-o -10 -15 -20 -25 0 30 60 90 Observation Angle 00, deg. Figure 2.11: o)e iteratively refined backscatter pattern of a metallic cylinder (radius=0.5A; height=1.5A); the impedance matrix is 36% full.

38 5 -0 -m -5 -o-10 c b-15 --20 --25 Observation Angle p,, deg Observation Angle 9., deg. (a) (b) Figure 2.12: ac, backscatter RCS of a 3A long almond (a) Azimuth plane (b) Elevation plane CI) u 10 -20 -20 0 30 60 90 120 150 180 Observation angle 0 (deg.) Figure 2.13: cso backscatter RCS pattern of a 7.5A long missile of radius 0.68A (azimuth cut).

39 0:..........,I....~..... - A I..... 50 40 30 U 20. 10 o -10 0 30 60 90 120 150 180 Angle 4 (deg.) Figure 2.14: Comparison of our technique and the FAFFA

CHAPTER III The Fast Multipole Method and guidelines for using the method to calculate the RCS of large objects The Fast Multipole Method (FMM) can be described succinctly as a recipe to reduce the number of multiplications and consequently accelerate the computation of matrix-vector products of an iterative solution. It also saves a significant amount of memory since only the near zone matrix entries are stored explicitly. In this chapter, we describe the methodology of the FMM and its variants. Each version of the FMM is associated with inherent approximations and a goal of this chapter is the comprehension of these approaches. At the end of this chapter, we discuss the choice of solution parameters for these algorithms based on an error study. Chapter IV discusses applications of these techniques to electromagnetic scattering simulations. 3.1 The Fast Multipole Method As stated above, the FMM is a fast method for calculating the matrix-vector products resulting from a discretization of an integral equation. For the system (2.1) the pertinent matrix vector product is obtained from the discretization of the integral equation and is given in (2.9). We examine three versions of the FMM to accelerate the matrix-vector product computation. Specifically, the exact FMM [20][21], a windowed FMM [23] and an approximate FMM [42] are examined. 40

41 The scattered field in (2.1) is given by E = - ZoJJ(H(2(kop )dl' pp E C (3.1) In the next few subsections we examine the evaluation of the integral (3.1) using the various forms of the FMM. By using consistent notation and appropriate explanation of the procedures, this exposition provides the first comparative look at the characteristics of the FMM techniques. Consequently, one is able to make conclusions on their CPU performance and accuracy compromises. In accordance with the FMM (see Figure 3.2), the N unknowns introduced for the discretization of (3.1) are subdivided into groups with each group assigned M unknowns. Thus, a total of L - groups are constructed. The key step in all FIMM procedures is to rewrite the integral (3.1) as a product of terms each being a function of p (observation point) or p' (integration point) but not both. In this manner, the evaluation of the integral is carried out by considering the group to group interactions separately from the inter-group interactions. Beyond the math, this breakdown of interactions/operations can be viewed in the context of the manager-worker model. Basically, we can view each group as managed by the center element with the workers comprising the elements of the group. Communication/interaction among the groups takes place through the managers which in turn interact with the group elements. However, this type of model is based on certain simplifications/decompositions of the original boundary integral. The decompositions reduce the direct interdependence of each group member with the other elements belonging to different groups and this is at the heart of the CPU speed-up afforded by FMM. However there are inherent approximations as part of the group decomposition process which must be understood in order to assess the accuracy of each FMM algorithm. Below we discuss the decomposition procedure employed in three of the fast multipole methods.

42 Basis elements / I Test element / I \ I P w f f --- - - - - -- Group center \ \./. Group Ei) \f/ (Test group) /. Group I Source element _ (Source group) Y / Kor/ I Ye- - -- Global origin ' ~ Group center Figure 3.1: Computation of the boundary integral matrix vector product using exact FRMM 3.1.1 Exact FMM To achieve the decomposition of (3.1) into a product of functions in 7p and '. wwe first invoke the addition theorem to rewrite the Hankel function as Q/2 H2(kolpwol l-pI) E J( ol a-P)H2)i(koPI) ejn("' - ) py, > I PI-Pl\ n=-Q/2 (3.2) where plt denotes the distance between the centers of the 1 and 1' groups as illustrated in Figure 3.1 and 3.2. Also, qI/w and Opip are the angles between the vectors p1il and p, - pi and the x-axis respectively. The source and observation points p', and p have their origin at the center of the 1' and I groups, respectively, while p' and p are measured from the origin. Typically, the semi-empirical formula Q/2 = koD + 5 ln(koD + 7r) (3.3)

43 I - /I I Group P/i,/ N I Mb At /s if 'I/I a> dr Far I I Group GrouGroup d K + K K- 1 min - - nin ' Jr- ' /- I -N' - group | P,1< d min => Near group y x Figure 3.2: Computation of the boundary integral matrix vector product using exact FMM is used to truncate the sum (3.2), where D is the diameter of the circle enclosing the groups In general, Q/2 = M, assures convergence (It will be shown that Q is the number of directions in which the radiation of the group is sampled. With M being the number of basis elements in the group, Q = 2M satisfies the Nyquist criterion for faithful replication of the source group radiation). Next we introduce the Fourier integral of the Bessel function Jn(ko pl - ) = p1 l eJ,-pn(O-0,p+/2)do Pi 27r r and in conjunction with (3.2) we can now rewrite (3.1) as ES) = - k j V()()e-kd where k = ko0(x cos 0 + y sin q$) is measured from the x-axis. In this, V,( = f Jz (p')ejk''dl is identified as the far-field pattern of the source group and Q/2 T'll ()= E Hn2)(kopl, )e-jn(n-q1+r/2) n=-Q/2 (3.4) (3.5) (3.6) (3.7)

44 is referred to as the translation operator providing the group-to-group (1 to 1') interactions. From (3.5)-(3.7), we observe that we have accomplished the decomposition of the integral (3.1) into terms which separate out the dependence on p and p'. The final evaluation of Es proceeds by discretizing the integral over q to yield the expression k0Z7 A( (3.8) 87r q=l which is the radiated fields from some location in the source group 1' to a point within the receiving group 1. Note that AX = 27'/Q indicates the angular spacing between the propagation vectors of plane waves emanating from a group. Thus Oq = qA4, q 1... Q, whereas kq = ko(x cos qOq + y sin qq). As mentioned earlier, the number of plane wave directions is set equal to twice the number of elements in the group (Q = 2M), thus, satisfying the Nyquist sampling theorem with respect to the integration over (. Given the above steps, the exact FMM procedure for carrying out the matrix vector product can be summarized as follows: 1. Compute pattern of the source group (aggregation). Mathematically, this corresponds to evaluating Vc(hq) given in (3.6). The evaluation of VI(qq) for a single source group and at a single direction requires M operations, corresponding to the number of elements in the group (the integration over the line segment is performed as a summation). Consequently for L groups and Q directions for each group, the operation count is QML. 2. The next step is to employ the translation operator to evaluate the pattern of a source group at the center of the test group. Mathematically, this operation amounts to computing the coefficient Ai(qq) = VI(qq)Tiit(qq). The evaluation of Ai(qq) involves an operation count of QL2, where again L denotes the number

45 of groups and Q is the number of directions. 3. Finally, at the receiving group, the fields are redistributed (disaggregation). Mathematically, this amounts to computing the expression Es(p) as given in (3.8). Evaluating the sum at a single point requires Q operations. Thus for L groups, each containing M unknowns, the operation count is QLM. From the above we conclude that the total operation count of the above three steps is C1QML + C2QL2. Also, the near field (by this we mean that groups in the near vicinity of each other are treated using the standard method of moments procedure) operation count is of O(NM). On choosing, Q - M, sufficient to achieve convergence, the operation count of the three steps reduces to C1MN + C2-. On setting M N,/A: implying L = /N (an optimal choice) the final operation count is N1'5 and this should be compared to the usual N2 operation count of the direct CG or LU solvers. The reduction of the operation count from O(N2) down to O(N1'5) is indeed dramatic. An appreciation of the CPU reduction can be acquired by setting, for example, N = 2000 which is a relatively small number of elements. However, further improvements can still be achieved by nesting groups leading to the multilevel FMM [22]. 3.1.2 Windowed FMM In the exact FMM, the translation operation between groups assumed isotropic radiation. However, it is suggestive that the groups would interact strongly along the line joining them and less so in other directions. Indeed, it was shown in [23] that the translation operator could be contemplated as composed of a geometrical optics (GO) term (along the line joining the source and test group) and two diffraction terms associated with the shadow boundaries of the GO term. To illustrate the

46 validity of this concept, we plot in Figure 3.3 the translation operator for different group separation distances along a segment of length 50A. For this example, the number of unknowns on the segment was 750, resulting in 27 groups. As seen, the "lit" region of the translation operator narrows as the group separation distance is increased, eventually displaying the predictable sinc function behavior for large group separation distances. The tapering off of the translation operator from a value oscillating around 2 down to zero for larger q - q115 values is characteristic of the geometrical optics plus diffraction terms in the context of traditional high frequency methods. We may also comment that this high frequency model enables the identification of a lit region even for groups which are not widely separated (for example, see Figure 3.3 for the translation operator between groups 1 and 3. The key characteristic of the windowed FMM is the exploitation of the diminished value of Tw,(q) for large q-q$. Basically, in the windowed FMM, the computation of Till$() for these angles is avoided altogether. This can be accomplished by multiplying Tili(X) with the filter (windowing) function W ( = { (q5 q - q1ii'I <i3s) Wyil'W = \ -Ce(jq-l-()2 (q - q-011'1 > ds) (3.9) where = sin 2kop (3.10) and a is a taper factor to be specified. Note also that /3s was selected to provide a larger bandpass window when pal' is smaller as dictated from high frequency analysis.

47 The discretized plane wave expansion can now be written as ES = - A koZ q W,(q)T() ()e- (3.11) 87I' q=l By taking into account only the non-zero sector of Wui/(q), the operation count of the translation process is now reduced to C3L2 N2 with the corresponding total operation count given by C1,MN +C4. Grouping the unknowns into N1/3 elements per group, results in a total operation count of O(N4/3). This should be compared with the O(N3/2) operation count of the exact FMM. The computation of the boundary integral matrix vector product by employing the windowed FMM is depicted pictorially in Figure 3.4 illustrating that the filter function has the effect of eliminating plane wave interactions at directions away from the line joining the interacting groups. 3.00,, ' ' Groups 1 and 3 -. 200 \Q=54;P, 45 X 3 l Groups and 8 Groups 1an 27! Ir 1.00 - -180 -120 -60 0 60 120 180 - z/'(de~gees) Figure 3.3: The Translation operator for different groups on the boundary of a 50A long segment; 750 BI unknowns; 27 groups

48: ~1=0 i f fft = Group ~, Group Group,Group Nb j+b d m K+l K / KWindow - -' - -- z x 'P,- > dmin - ' Far group | P /< d min Near group Figure 3.4: Computation of the boundary integral matrix vector product using windowed FMM 3.1.3 Fast Far Field Algorithm This is an approximate version of the FMM since the algorithm is based on introducing the large argument approximation of the Hankel function. That is, the approximation H - ) ( koP-l /l e-Jk PP 6-jkoPm e - o P 1 T (3.12) V kopp-j is used. As shown in Figure 3.5, pl'i is the distance between the center of the test group 1 and the center of the source group l'; pnl' is the distance between the nth source element and its group center and pim is the distance between the mth test element and its group center. The introduction of the large large argument expansion necessitates that the FM procedure be used only for groups which are very well separated. However, (3.12) allows for the immediate decoupling of the test-source element interactions, thus, enabling the computation of the matrix-vector product for far-field groups with a

49 I iI G -u,' Group Group,'Group p d K+ K / K/ ' ', min K + /I,I, /,, " m n d im m IPI min Far group | P m < d min Near group Figure 3.5: Computation of the boundary integral matrix vector product using the FAFFA reduced operation count. This is illustrated below by going through the three steps of the FMM 1. The aggregation of source elements in a single source group now involves M operations, corresponding to the number of elements in the source group. Specifically, M Vl= jje-jko0Pi'-Pnl (3.13) and since the above aggregation needs to be done for all source groups, the operation count becomes O(-M) 0(N), where - represents the total number of groups. Also this operation, being dependent only on the test group rather than the test element, needs to be repeated for all Q - test groups leading to a total operation count of 0(N 2) for aggregation. It should be noted that use of the large argument expansion, rather than the addition theorem for the Hankel function, results in the aggregation sum being a function of the test group also (V/1) unlike the exact FMM where the aggregation sum is a function of source group only (V1'(q)). Thus, the technique by which the Exact FMMl and the FAFFA reduce the operation count differ in the fact that while in the

50 exact FMM the aggregation sum is characterized by a source group (I') and a direction (q) which is not interwined with the test group direction, the aggregation sum in the FAFFA is characterized by the source group (1') and the test group (1). 2. The main advantage of FAFFA is due to the faster computation of the translation operator. We have Ail = TlVIt (3.14) where in the FAFFA the translation operator simplifies to e-jkop1, l Tltl= 0 —kp (3.15) V \/kopl,l This should be compared to the sum (3.7) for the exact FMM. Clearly, (3.15) needs to be done only at the group level and involves O( 2) operations for all possible test and source group combinations, making it the least computation — ally intensive step. 3. The disaggregation or redistribution process is again the operation kz:Z~ N/M E (p) = - ~ Aie-ijkOP (3.16) 1'=1 Since this operation involves only the source group instead of the source element, it needs to be done for each source group, implying 0(M) operations to generate a single row of the matrix-vector product. To generate M rows, corresponding to a test group, the operation count would be O(N). With - test groups, the operation count is O(N ) Consolidating the above three steps for the FAFFA algorithm we have N2 Op.count - CiNM -+ C2 (3.17) MI

51 where the first term refers to the operations associated with the near field terms. As before, M = VN and the total operation count is O(N1'5). While, the operation count for this algorithm could be further reduced down to O(N1'33) by performing the process of "interpolation" and "anterpolation" as described in [42] for very large objects, we found that the accuracy deteriorated for the considered applications. Hence only the O(N1 5) version was used. 3.1.4 Logic Flow The operation counts described in the previous section for the various algorithms are illustrated with the help of flow diagrams and sections of code from the computation of the matrix vector products for the far groups. Figures 3.6 and 3.9 depict the flow diagram and code for computing the matrix vector product in the exact FMM. It is seen in Figure 3.9 that each of the aggregation, translation and disaggregation operations consists of a single multiplication which is described below * The aggregation operation consists of the product of an entry of the trial vector represented as (Dum(J) in Figure 3.9) with an aggregation factor, represented in Figure 3.9 for the Jth element and Kith direction as SrcGc(J,K). This is given by SrcGc(J,K) = AJeik~(xco60 n) PJc where Aj is the length of the Jth discretization element, qK is the Kt th radiation direction and P'JGr is the direction vector of the Jth element, measured from the center of the group(JGr) it belongs to. The result of the aggregation operation yields a term characterized by only the source group and radiation direction (represented in Figure 3.9 as V(JGr,K)). * The translation operation involves the multiplication of the aggregation sum, V(JGr,K), with a translation factor, represented in Figure 3.9 for the IGrth test

52 group, J(rth source group and Ith radiation direction by Trans(IGr,JGr,K). K/2 This is given by Trans(IGr,JGr,K) = e-jn(K-'Gr,JcGr+)H(2) (0korIGJGr) n=-K/2 where qk4 is the KIth radiation direction, rlGr,JGr and q5IGr,JGr are the distance and angle between the IGrth and JGrth groups. The result of the translation operation yields a term dependent only on the test group and radiation direction (represented in Figure 3.9 as GrGr(IGr,K)). * The disaggregation operation involves the multiplication of the translation sum, GrGr(IGr,K), with a disaggregation factor, represented in Figure 3.9 for the Ith test element and the Kth radiation direction by TestGc(I,K). This is given by TestGc(I,K) = e-(jk0cosK ^+1 sin I)'PIGr where pIGr is the direction vector of the Ith element measured from the center of the group(IGr) it belongs to. The result of the disaggregation operation yields a term dependent only on the test element alone and is the contribution to the Ith entry of the product vector. The windowed FMM differs from the exact FMM in the translation phase and this is illustrated in Figures 3.7 and 3.10. These figures illustrate that the windowed FMM achieves its reduced operation count by eliminating some of the directions in which plane wave interaction takes place. The innermost loop in the translation phase has an operation count which is a constant (15-25 in our simulations, 20 in [23]) and is a significant reduction from the corresponding operation count in the exact FMM. The technique by which the FAFFA achieves its speed-up is depicted in Figures 3.8 and 3.11. It is seen that the FAFFA "recycles" the plane wave spectra of the source group. For a given test group, the aggregation and translation operation are performed only once for each source group, necessitating that only the disaggregation

53 operation needs to be performed for each individual element of the test group. Similar to the exact FMVI, the aggregation, translation and disaggregation process consist of a single multiplication. However the factors used in the three processes and the method by which the reduced operation count is achieved are different. * The aggregation operation again consists of the product of an entry of the trial vector with an aggregation factor, represented in Figure 3.11 for the Jth element and IGrth test group as SrcGc(IGr,J). This is given by SrcGc(IGr,J) = AJC-ijkPJGrIG'PJ,J7Gr where Aj is the length of the Jth discretization element, PJGr,IGr is the unit vector along the line joining the source and test groups while PJJGr is the vector along the line joining the source element with its group center. Thus, an aggregation sum is formed for each combination of source and test groups. * The translation operation involves the multiplication of the aggregation sum with a translation factor, represented in Figure 3.11 for the IGrth test group and JGrth source group by Trans(IGr,JGr) and given by Trans(IGr,JGr) = e '-jkJG, Again, the translation operation needs to be done for each pair of /k0OPJGr,IGr test and source groups. * The disaggregation operation involves the multiplication of the translation sum, GrGr(JGr), with a disaggregation factor, represented in Figure 3.11 for the Ith test element and the JGrth source group by TestGc(I,JGr). This is given by TestGc(I,JGr) = e-jk~JGrIGrpIGr I, where IGrI is the vector along the line joining the test element with its group center. It should be noted that to compute the interactions between a pair of groups, the aggregation and translation need to be done only once and thus the crux of the FAFFA is

54 indicated with highlighted sections of code in Figure 3.11. Aggregation Set source group counter = 1 (Start with the first source group) Set direction counter =1 (Start with the first radiation direction) Set element counter = 1 (Start with the first basis element in the source group) Perform aggregation operation for a single element of the source group for a single M M irection Operation count = 1 ) Increment element counter by one Is element counter > No number of elements in the group \(M)? Yes Increment direction counter by one Is direction counter> number \ No of radiation directions (M)9 Yes Increment source group counter by one Is source goup No counter > total number f groups (N\, Yes End aggregation.i Go to translation. Aggregation op. count (M) (M) (M) - NM V) = NM Translation Set test group counter = 1 (Start with the first test group) Set source group counter = 1 (Start with the first source group) i; Set direction counter = 1 (Start with the first radiation direction) Perform translation operation for a single N pair of source and test N - groups and for a single M direction Ops Operation count = 1 ) OP Increment direction counter by one Is direction counter > No number of radiation directions \(M)? yYes Increment source group counter by one [ s source group counter > total number No of groups Yes Increment test group counter by one Is test group No counter >total number f groupsN?/ Yes End translation. Go to disaggregation. / Translation op. count (M) (N) () = ( ) Disaggregation Set test group counter = 1 (Start with the first test group) Set direction counter = 1 (Start with the first radiation direction) *' --- Set element counter = 1 (Start with the first basis element in the test group) Perform disaggregation Ioperation for a single N element of the test - group for a single M M direction s Operation count = 1 ) Increment element counter by one Is element counter > No number of elements in the group (M)? Yes Increment direction counter by one Is direction /counter> number No of radiation directions (M)? Yes Increment test group ounter by one s \M/ st group No Yes End disaggregation. Matrix-vector product done/ Disaggregation op. count (M) (M) ()- NM N Figure 3.6: Sequence of operations to be performed in the Exact FMM

55 3.2 Guidelines for using the method to calculate the RCS of large objects Although, we demonstrated in the previous section that the fast multipole method has O(N15) or even lower computational complexity, there are several parameters which play a role in the CPU requirements and accuracy of the solution. More specifically, the grouping scheme, the sampling rate and the group size, all have an effect on the performance of the FMM. In this section we examine the role of these parameters on the accuracy and computational efficiency of the FMM on the basis of some error-criterion. Specifically, we considered the scattering by two dimensional metallic structures using a solution of the EFIE and the MFIE with pulse basis and point matching as illustrated in section 2.1 [431. In our study, the Fast Far Field Algorithm (FAFFA) also referred to as the approximate FMM was used for carrying out the matrix-vector products. Referring to Figure 3.12, the following parameters were examined with respect to the accuracy and efficiency of the FAFFA: * the near-group radius,dmin * the sampling rate * group size, and * memory requirements Our benchmark for accuracy was the RMS error given by ERRORRMS = E [RCSREF(i) - RCSFF(i)] (3.18) i=1 where RCSREF denotes the reference radar cross sections as computed by the stan — dard moment method approach without grouping, RCSFF is the RCS calculated using FAFFA and M being the number of points at which the RCS is computed.

56 We note that this error formula is among several that have been considered, but has been found to represent a reasonable measure of the accuracy. It has been employed by Schuh and Woo [44] for a study on the accuracy of RCS computations by various codes. In the following we will show error curves as a function of the aforementioned parameters for two cylinders. One is a PEC isosceles triangle having a base 2.5A in length and a height of 17.95A. The other is a rectangular PEC cylinder 25A x 4A in size. These geometries (see Figure 3.13) are analyzed at grazing (0 degrees) incidence since this excitation is most stressing for RCS calculations, particularly so for the triangular cylinder. In fact, the two geometries were selected to have different characteristics in the backscatter region. 3.2.1 RMS Error as a function of near-group separation distance and sampling rate First we look at the RMS error of the bistatic RCS pattern as a function of the near-group radius, dmin. Such error curves are shown in Figures 3.14 and 3.15 for the isosceles and rectangular cylinders, respectively. These figures depicts several curves, each corresponding to a different sampling rate (A/10, A/15, A/20, A/30 or A/40 as designated in the figure). Surprisingly, the error curves for each geometry give rise to the same conclusion. Basically, Figures 3.14 and 3.15 demonstrate that the sampling rate has a profound effect on the accuracy of the solution and the value of the error is strongly dependent on the near-group radius. For example, if dmin is set to 1.7A, a tessellation rate of 10 segments per wavelength yields an RMS error of 4.79 dB whereas a tessellation rate of 30 segments per wavelength leads to an RMS error down to 1.19 dB. As expected, a higher sampling rate leads to smaller errors. However, these errors are different and much larger than those resulting from a standard moment method implementation. Specifically, a tessellation rate of 10

57 segments per wavelength yields an RMS error of 0.1615 dB when employing the standard unreduced moment method solution but increases to 4.79 dB when FAFFA is employed with dmin = 1.7A. As depicted in Figures 3.14 and 3.15, the values of dmin plays an equally important role in achieving a given RMS error. For example, when the segment length is A/10, the error is over 4 dB when dmin = 2A but the error decreases to less than 2 dB when dmin reaches 3A for the triangular cylinder and 5A for the rectangular cylinder. This trend is, of course, logical since for larger values of dmin the FAFFA implementation resembles more and more the unreduced moment method. Having found such a strong dependence of the RMS error on the tessellation rate and the near-group distance, to proceed further with our study, it is essential to look at what may be an acceptable RMS error. This error value can then be used to determine acceptable relationships between sampling rates and dmin on the basis of a certain RMS error. We therefore refer to Figure 3.16 which shows the bistatic RCS pattern of the isosceles triangle as computed using FAFFA for different sampling rates. The shown (exact) pattern has a dynamic range of 30 dB and we observe that an RMS error over 4 dB is associated with large and unacceptable RCS deviations (over 10 dB) from the reference RCS values. However, when the RMS error is about 1 dB or so, it is clear that the FAFFA results have only small deviations from the exact curve which is actually tracked very well by the FAFFA results even in low RCS regions. Consequently, we may select the RMS error value of 1 dB as the threshold level in determining whether a given choice of FAFFA parameters lead to acceptable accuracy or not.

58 3.2.2 CPU requirements Returning now to Figures 3.14 and 3.15, we observe that an RMS error of 1 dB or less is obtained for the triangular cylinder if dmin > 4A when the sampling rate is 10 segments per wavelength, but dmin must be at least 7A for the rectangular cylinder with the same sampling rate. Unfortunately, these values for dmin are much larger than IA which is the value often suggested by other users of the fast multipole method and its variations. It is of course important to keep the near-group radius as small as possible for CPU and memory reduction. However, the sampling rate must be increased to 30 segments per wavelength (for the isosceles triangle) if we are to reduce dmin to 2A, but clearly the increased sampling is counter to the benefits from a smaller near-group distance. Thus in Figures 3.17 and 3.18, we show curves of the CPU time as a function of the near-group distance for different sampling rates. Figure 3.17 refers to the triangular cylinder whereas Figure 3.18 corresponds to the rectangular cylinder, and both show a similar dependence on dmin and the sampling rate. For both cylinders the CPU time increases almost linearly with dmin and is, of course, higher for denser tessellations. These curves can be used in conjunction with those in Figures 3.14 and 3.15 to obtain the CPU time to determine dmin for a given RMS error and sampling rate. Alternatively, Figures 3.14 and 3.15 can be consolidated with the data in Figures 3.17 and 3.18 to obtain a new set of curves which explicitly give the CPU time as a function of a desired/given RMS error for the selected sampling rate. These CPU curves as a function of a specified RMS error are shown in Figures 3.19 and 3.20, corresponding to the triangular and rectangular cylinders, respectively. The trends in each of these Figures are identical, demonstrating that the dependence between CPU time, sampling, near group-distance and RMS error are likely the same regardless of geometry considerations.

59 Basically, Figures 3.19 and 3.20 show that the CPU time increases quadratically with sampling rate. Clearly, a lower sampling rate is more attractive in terms of CPU time since it leads to smaller systems and therefore faster convergence rate. That is, it is better to use a larger C1min with lower sampling rates rather than a higher sampling rate with a lower near-group window radius. For example, to achieve an RMS error of 1 dB for the triangular cylinder, the CPU time is about 10 sec when using a tessellation rate of 10 segments per wavelength and from Figure 3.14 the corresponding dmin is 3.75A. However, to achieve the same error using a tessellation rate of 30 segments per wavelength, the corresponding CPU time is more than five times higher but the required near-group window radius is less than 2.5 wavelengths. It may seem that the higher number of degrees of freedom associated with higher segmentation rates leads to unequivocally higher memory requirements as well. However, this is not always the case with the FAFFA implementation because the value of dmin plays a major role on the bandwidth of the dense section of the matrix. For example, in the case of the rectangular cylinder, sampling at A/10 would imply 580 unknowns while sampling at A/30 would imply 1740 unknowns. To achieve a RMS error of 3.75 dB, it would require a near group window of 3.5A for sampling at A/10 and a window of 1S for sampling at A/30. Execution time considerations are of the order of 20 seconds for the first case and 75 seconds for the latter case. An examination of the grouping of the unknowns would reveal that the memory requirement for the latter case is smaller. Sampling at A/30 and employing a near group radius of 1A would mean a storage requirement of 42 elements per row. Sampling at A/10 and employing a near group radius of 3.5A would mean a storage requirement of 72 elements per row, which is higher than the previous case. Consequently, although Figures 3.19 and 3.20 suggest that a lower

60 sampling rate has the least CPU requirements in satisfying a given error criterion, the memory requirements will be higher when compared to those associated with a higher sampling error and the same error criterion. Memory considerations will of course be an issue for very large simulations and the choice of using a high sampling or not will depend on the available computing resources. The geometries and incidence angles selected in Figure 3.13 are among the most stressing for RCS calculations and result in relatively large error for small near-group window radius. For a smooth body such as the 50A diameter circular cylinder shown in Figure 3.21, the RMS error is small even with a small near-group window radius. However, the relations between near-group window radius, sampling rate, memory and execution time expounded in the previous section continue to hold for such smooth bodies also. It should be noted that while the FAFFA employs the far-zone Green's function which results in larger error, corresponding results for the exact FMM shown in Figure 3.22 and 3.23 show much greater accuracy. It is seen that the exact FMM is capable of replicating even the surface currents on the triangular cylinder for nose-on incidence with remarkable accuracy. (Figure 3.2.2). 3.3 Summary In this section we examined the CPU requirements of the FAFFA as a function of the different parameters affecting its performance, including the near-group window radius, sampling rate and error. We presented curves which show the required sampling rates, near-group window radius and CPU time as a function of a given error criterion. Based on these curves, one concludes that * The error is strongly dependent on dmin and the sampling rate. This error dependence on the sampling is particularly inherent to the FAFFA solution

61 method and should not be confused with similar errors associated with the unreduced moment method implementations. With higher sampling, the spatial extent of each group is reduced and thus the process of aggregation and disaggregation of group elements into group centers induces less error. * In lieu of CPU time efficiency, it is best to employ lower sampling and larger near-group window radius. * Stringent memory requirements behooves the adoption of high sampling and small near group radius but at the expense of high execution time.

62 Figure 3.7: Sequence of operations to be performed in the translation process of the windowed FMM

63 Is source element counter > total number of elements in the group (M)? Yes Set ITestGrNew = 0 (Indicates that the test group is handled for aggregation and translation and therefore only disaggregation needs to be computed for each test element)! _ ~ ~ I I --- — I Increment test element counter by one I Perform translation for a single source group & test group pair (1 operation) Total op. count (M + ) (N2) + (N) A g iM T al t M D Aggregation2 Translation2 Disaggregation2 Aggregation Translation Disaggregation I!.A. Is test element counter > total \ No number of elements in the group(M)? Yes Figure 3.8: Sequence of operations to be performed in the FAFFA

64 c NGr - Number of groups; IQAngles - Number of radiation directions c DPhi - Angular spacing between radiation directions c NElGr(NGr) - Array containing number of elements in each group c GetGlobal - function which gets the global element # given a group c number and local element number c DMin - Minimum distance beyond which two group are treated as far c groups c Distance - Function which returns the distance between groups c Dum - Vector multiplying the matrix - in the CG algorithm c will be the search and residual vectors c AX - Array which is the matrix-vector product c SrcGc,Trans & TestGc - aggregation, translation & disaggregation factors c respectively described in the text IQAngles = 2*NGr DPhi = 2*Pi/IQAngles c Aggregation - Operation count O(NM) do JGr =l,NGr do K = 1,IQAngles do JEl =l,NElGr(JGr) J = GetGlobal(JGr,JEl) if (S.eq.'*') then _N_ V V(JGr,K) = V(JGr,K) + Dum(J)*conjg(SrcGc(J,K)) Mv else V(JGr,K) = V(JGr,K) + Dum(J)*SrcGc(J,K) endif enddo enddo enddo Ic Translation - Operation count O(N^2/M) do IGr = 1,NGr do JGr = l,NGr if (Distance(IGr,JGr).gt.DMin) then do K = 1,IQAngles if (S.eq.'*') then N J GrGr(IGr,K) = GrGr(IGr,K) + conjg(Trans(IGr,JGr,K))*V(JGr,K) N N M else M M GrGr(IGr,K) = GrGr(IGr,K) + Trans(IGr,JGr,K) * V(JGr,K) endif enddo endif enddo enddo Ic Disaggregation - Operation count O(NM) ] do IGr = l,NGr do K = 1,IQAngles do IEl =,NElGr(IGr) I = GetGlobal(IGr,IEl) if (S.eq.'*') then N _ AX(I) = AX(I) + GrGr(IGr,K)*conjg(TestGc(I,K)) M else m AX(I) = AX(I) + GrGr(IGr,K)*TestGc(I,K) endif enddo enddo enddo Figure 3.9: Code indicating the computation of the matrix-vector product in the Exact FMM

65 c The symbols in this section of code represent the c same quantities as in the code for the exact FMM c Translation - Operation count O(NA2/M^2) do IGr = 1,NGr do JGr = 1,NGr if (Distance(IGr,JGr).gt.DMin) then do K = 1,IQAngles if (S.eq.'*') then [ if (abs(Trans(IGr,JGr,K)).eq.0) then continue else GrGr(IGr,K) = GrGr(IGr,K) + conjg(Trans(IGr,JGr,K)) * V(JGr,K) < | endif else \A | if (abs(Trans(IGr,JGr,K)).eq.0)then i continue else GrGr(IGr,K) = GrGr(IGr,K) + Trans(IGr,JGr,K) * V(JGr,K) [endif endif enddo endif ienddo enddo 1.10: Code indicating the computation of the matrix-vector product in the translation phase of the windowed FMM N M

66 c C C C C C: c c c The symbols in this section of code represent the same quantities as in the code for the exact FMM ITestGrNew is a counter which is set to 1 if a test group is "new". This means that since the aggregation and translation operations involve only the test group rather than the test element, they need to be done only once (for the first element in the test group). For the rest of the elements in the test group only the disaggregation operation needs to be done (corresponds to ITestGrNew = 0). do IGr = 1,NGr ITestGrNew = 1 do IE = 1,NEIGr(IGr) I = GetGlobal(IGr,IEl) do JGr= 1,NGr if (Distance(IGr,JGr).gt.Dmin) then if (ITestGrNew.eq. 1) then I V = (0.,0.) do JEl = 1,NElGr(JGr) J = GetGlobal(JGr,JEl) M if (S.eq.'*') then V = Dum(J)*conjg(SrcGc(IGr,J)) + V else V = Dum(J)*SrcGc(IGr,J) + V endif enddo if (S.eq.'*') then GrGr(JGr) = conjg(Trans(IGr,JGr))*V else GrGr(JGr) = Trans(IGr,JGr)*V endif lendif. if (S.eq.'*') then AX(I) = AX(I) + GrGr(JGr)*conjg(TestGc(I,JGr)) else AX(I) = AX(I) + GrGr(JGr)*TestGc(I,JGr) endif endif enddo ITestGrNew = 0 enddo enddo Figure 3.11: Code indicating the computation of the matrix-vector product in the FAFFA

67 < dmin = Near group Group K Figure 3.12: The process of grouping unknowns - two groups are in the near field of each other if the distance between their centers Pil is less than dmin y 18k 2.5hX (a) y 4X4[ ----~ --- -1 25; (b) Figure 3.13: (a) An isosceles triangular metallic cylinder - 17.95A high with a 2.5A base (b) A rectangular metallic cylinder - 25A x 4A

68 5.0 4.0 3.0 2.0 I I I I I I..... — 0 — Sampling- /10 - -& - Sampling - XJ15.............. Sampling - V20 — v — Sampling - /30 ------ Sampling - /40 "0 1d Cz t-. w 1.0.0 Figure 3.14: Error cur high with K I I I I I I I I I I... 0 1 2 3 4 5 6 7 8 9 Radius for near-group treatment, dmin (X) ves for the FAFFA for a isosceles triangular metallic cylinder - 17.95A a 2.5A base "0 - U.t:P 6.0 5.0 4.0 3.0 2.0 1.0 ' I I I I ' I. ' I I I ' I ' I ' — 0 --- Sampling - V10 - -n - Sampling- h/20 4\: ' - ' \a ':\ \ A,, \q, ', a,. n 0 1 2 3 4 5 6 7 8 9 10 11 12 13 Radius for near-group treatment, dmin (X) Figure 3.15: Error curves for the FAFFA for a rectangular pec cylinder - 25A x 4A

69 20 F I.. I I.. I. I I..... I.... I... I I...... 10 - I ~3. P4.,4 CT(rA 5. 0 It -lo -,,",4 1: I I1: I. I1 AA 44 ^ I.I I, I;,,, -,, FAFAAA=.7-amln As0RSEr~ 1, '",,,...,l~, ' ~ i: ' It MoM (E~act)...... FAFFA-dm-in:-l.7X —sampling VI0-RMS Err —4.79 A FAFFA-d min-l.7.-.sampling M15-RMS Err=-2.61 FAFFA-dmin=l.7)3.-sampling X/30-RMS Err=-l.19 -20 -30 -40............... l l | 0 30 60 90 120.......... 150 180 Observation angle ~ (deg.) Figure 3.16: Comparison of bistatic patterns of an isosceles triangular pec cylinder - 17.95A high with a 2.5A base computed using theFAFFA at different sampling rates GO 0 <u 0 CA3 O >-^ O 0 lN 0 0 c~' 0 cS D. - <Uj u ~I0 C) 225 200 175 150 125 100 75 50 25 n ~-o-_-__ Sambling-_A0 - I' 0'. --- - ' - Sampling - X/15 /-..-A..Sampling - V20 ----- -Sampling - V30 --— ~> ---- Sampling - J40 A/............I&....... /~........................ E-l - - - -... -..-/~ E3. --- — ~_ _ i_ ~ ~ ~ ~..... -- - A - - "" V 0 1 2 3 4 5 6 7 8 9 Radius for near-group treatment, dmin (X) Figure 3.17: Time curves for the FAFFA for a isosceles triangular pec cylinder - 17.95A high with a 2.5A base

70 a) -\ O O O C3 0 O' o cr * —4 0, E c) 250 225 200 175 150 125 100 75 50 25 n...................... I I I 'I I I' I '' I...4' I ' I.. -''.oZ~!..' ---- Sampling - J10 - -- - Sampling - 20 '.." - -. --- —- Sampling - /30 -. - / - - -' I. I I I. I I v 0 1 2 3 4 5 6 7 8 9 10 11 12 13 Radius for near-group treatment, dmin () curves for the FAFFA for a rectangular pec cylinder - 25A x 4A Figure 3.18: Time cn Q c - Ct O o. 0u.a 0 E ~,. r 225 200 175 150 125 100 75 50 25 n -.... -.......... -— o --- Sampling - 10 -- -- Sampling- X/15 - -....... --- Sampling- X/20 -. — V — Sampling- V30 ' ------- Sampling- X/40 -- - - - - - - --- - —, ( - - - 4- -,,,,, - 1, ^ _ "" - "- - - -A,-..-.-....................................0 1.0 2.0 3.0 4.0 5.0 RMS Error in dB Figure 3.19: Time - error curves for the FAFFA for a isosceles triangular pec cylinder - 17.95A high with a 2.5A base

71 0 0n 0.0 0) 1 -011 S4:250.225:20-0 175 150 125 100 75 50 25 0 I -T- I I I El -a - - Sampli..................A....... Sampli ling - VIO 7 ling -X./20 ling -X.30 - 0, - - -- -- -- - - - - -0 cb -.0 ------— 0 I I I I I 0 1 2 3 4 5 6 RMS Error in dB Figure 3.20: Time - error curves for the FAFFA for a rectangular pec cylinder - 25A x 4A - I I 1.5 I I I I i 1.0 [ K >Sampling - -10 - -- Sampling - I20 — l-Sampling -X/20 K >' — ---- - - - - - - - - A =-I I- - 0 0% Ck 0 A0 S4 15 —..........'..... ---- Sampling - A/10 — l — Sampling - A120 10 - -a ---Sampling - V30 -El — - 5~~-E.5 -n 0 5 10 15 20 25 Radius for near-gmup treatment, d.j (X) 3 0 5 10 15 20 25 Radius for near-gmup trvatment d,,. (X) 30 Figure 3.21: Error and time curves for a smooth body - a circular cylinder of diameter 50A

72 y Segment # 50 x 18~ Segme #410 2.5 X Segmen 1 Segment # 770 (a) I 2 -- Moment Method ) 1 [ Exact FMM 3 1 I 0 100 200 300 400 500 600 700 Segment number (b) Figure 3.22: (a) Triangular cylinder geometry with discretization segment numbering (b) Currents on the triangular cylinder for nose-on incidence

73 20 C —/ PQ u *4-.c, sr 10 0 -10 0 30 60 90 120 150 180 Observation angle t (deg.) (a).20.15 "0 -/ 01.10 ' I. I ' I ' I ' I ' I ' I. I ' I ' I ' -E- Sampling - X/10 ------- Sampling - V20 --- - - - - - - - - - - - - - - - - - -.05.00 0 1 2 3 4 5 6 7 8 9 - - - 10 11 Radius for near-group treatment, dmin (X) (b) Figure 3.23: (a) Bistatic RCS computed at nose-on incidence for the geometry of Figure 3.2.2 using the moment method and the exact FMM (b) Error curves for the exact FMM

CHAPTER IV Applications of the Fast Multipole Method In this chapter, we present two applications of the FMM. First, to model node radome geometries of realistic electrical sizes, wre employ a version of the FMM to accelerate and reduce memory requirements of the EFIE formulation presented in Section 2.3. Next, we discuss the savings in solution time and memory when the FMM is used to compute the boundary integral matrix-vector products in a hybrid FE-BI solution. 4.1 Computation of nose radome scattering by employing the fast multipole method to calculate the RCS of large objects Nose radomes serve as enclosures for antennas and are generally pointed to reduce aerodynamic drag. The performance of a radome is generally described by parameters [45] which include the far-field radiation pattern, power transmittance, boresight error and boresight error slope. Approximate methods for treating the propagation of the plane wave through a radome include high frequency techniques which typically consider the radome to be locally plane and omit guided waves as well as interactions between the radome sections. Also, treatments of the higher order interactions between the radome and the antenna have so far been of little attention. Clearly, a more accurate approach is to employ an exact analysis method such as 74

75 the moment method technique [30] which is capable of including all first and higher order phenomena. However, due to the traditional 0(N2) storage and 0(N3) CPU requirements the method of moments approach can only be used for small radome structures. Recently though techniques have been introduced which can reduce the CPU requirements down to 0(N1'5) or less for large scale simulations. The fast multipole method (FMM) introduced in section III is one of these techniques and accomplishes the CPU reduction by lumping the far zone moment method elements into groups. The groups are subsequently interacted (rather than the elements) to achieve the purported CPU reduction. In the FMM implementation, the near-zone and self-cell elements are evaluated unaltered and thus the accuracy of the original method of moments formulation is retained. 4.1.1 Formulation To make use of the FMM to speed-up the solution it is necessary to employ iterative methods. In this case it has been shown that the FMM can reduce the operation count down to O(N1 5)/iteration, a substantial reduction. A detailed study of the parameters affecting the accuracy and solution time of this method is given in section 3.2 [24]. To employ the approximate version of the FMM, the unknowns are divided into groups with M unknowns in each group and thus the number of groups will then be M. For large source to observation distances, the kernel in (2.46) is approximated by using the large argument expansion as 3-jkoR. -jkorl/1 6-ejkoilrji -jko- ll r-e (4.1) R ruil where rv/j is the distance between the center of the test group I and the center of the source group 1'; rjlt is the distance between the jth source element and its group center and r1i is the distance between the ith test element and its group center (see figure

76 4.1). It is important to note that for large source to observation distances integration Group - ri < dmin Near group Group 1 >. N, i/ ' \ Group 1min/ ri > dmi. n Far group J Group K Figure 4.1: The process of grouping unknowns - two groups are in the near field of each other if the distance between their centers rl,l is less than dmin over the triangular domains is accomplished by mid-point integration; thus the source and observation co-ordinates in (2.46) are replaced by the centroid of the triangular domain of integration. Also, while the magnitude terms represented by the basis functions are more easily approximated because of their non-oscillatory nature, the phase term in the kernel needs to be accurately computed. The decoupling of testsource element interactions in the kernel as in (4.1) enables the computation of the matrix-vector product for far-field groups with a reduced operation count as detailed in the following sequence, where we have considered only the oscillatory term in the kernel. 1. For each test group, the aggregation of source elements in a single source group involves Al operations, corresponding to the number of elements in the source group. The aggregation operation corresponds to M bll = E I-ko-rl^ -r'j (4.2) j=l

77 2. Since the above aggregation operation needs to be done for all source groups the operation count becomes O(-M) 0 O(N), where N represents the total number of groups. Also this operation, being dependent only on the test group rather than the test element, needs to be repeated for - test groups leading to a total operation count of 0(N-) for aggregation. 3. The next step would be a translation operation corresponding to e-jkorll C' = -bl'l (4.3) rill Since this needs to be done only at the group level, it involves O( 2) operations for all possible test and source group combinations and is the least computationally intensive step. 4. The final step in the sequence would be the process of disaggregation corresponding to the operation N/M Iil = E Clle-jkorl'r (4.4) 1/=1 Conceptually, this process is the converse of aggregation. Since this operation involves only the source group instead of the source element it needs to be done for each source group thus implying an O(-) operation to generate a single row of the matrix-vector product. To generate M rows corresponding to a test group the operation count would be O(N). With - test groups, the operation count would be of O(N-). 5. The near field operation count being of O(NM) and the far field being O( ) gives a total operation count of N2 Op.count ~ C1NM + C2 M(4.5) M( 4.5

78 Typically, we can set the elements in each group, M = VW and as a result the total operation count is O(N1'5). 4.1.2 Results The formulation was first validated for resistive plates. The validating code employed was based on the resistive boundary condition and representation of the fields employing the Stratton-Chu integral equation [46]. Figure 4.2(b) shows the backscatter RCS of a 1.16A x 0.85A plate in the ~ = 0~ plane for the )( polarization. Results for a metal plate and a resistive plate with normalized resistivity of 2.12 - j0.2 are shown and comparison between the two sets of data is very good. Figure 4.2(c) shows the corresponding plots for the 00 polarization. To validate the typical nose radome geometry, available scattering data for the Von Karman radome were used for reference. Although, this is a BOR structure, it should be remarked that our formulation was not specialized to this class of geometries. The generating curve for this radome (see figure 4.3) is given by [47] r /= D {cos- (1- - ) - 2 sin [2cos- - (- (4.6) with r = 2 + 2 (4.7) and D is the diameter of the radome base whereas L is the length of the radome. The thickness of the dielectric shell was h = 0.05m and the dielectric constant was Cr = 4.0. The resistive sheet condition (R = k o(e-)h) gives an equivalent resistivity of -jl.061. Discretization of this geometry at twelve points/wavelength, results in 187 nodes, 362 triangular facets and 549 edges (unknowns) as shown in figure 4.3(b). The bistatic RCS with nose-on incidence is shown in figure 4.3(c) and (d). Also shown in these figures is the comparison with a surface integral formulation [30]. As

79 seen there is very good agreement between the two codes even though the RCS levels are very low at angles away from forward scattering. The incorporation of the fast algorithm alleviates the limitation on the size of the bodies analyzed as can be inferred from the results shown in figure 4.4. In this figure a metallic nose radome 10A long with a circular base of diameter of 1A is analyzed. Again, a body of revolution was the geometry of choice because of the availability of validation data from CICERO [48]. Discretization of this radome results in 3120 nodes and 6204 triangular facets. Comparison of the asatte RCS fom the from the two codes (see figure 4.4) reveals very good agreement for the c0 polarization while the 4 q polarization shows some deviation particularly at the low RCS levels. The near group radius employed in the FMM solution was 2A and the number of unknowns was 9324. Solution times were 12 seconds/iteration for the FMM radome code while the conjugate gradient solver employing the full, unreduced system would need an estimated 60 seconds/iteration. Estimates from scaling smaller problems indicate that LU decomposition would need a solution time of about five and a half hours while requiring an unrealistic 0.67 GB of memory for storage of the full matrix. An important advantage of integral equation analysis for the nose radome problem is that it avoids modeling the free space between the radome and nose antenna (see figures 4.5 - 4.8) thus reducing the computational cost. Results which include the interactions between the nose-antenna and the radome are given in figures 4.5 - 4.8. Unlike the Von Karman radome which was generated by an algebraic equation, the radome in figures 4.5 - 4.8 was generated from user specified elliptical or circular cross sections, which were then interpolated with cubic splines to generate the cross section at any intermediate point along the axis. Figure 4.7 depicts the (o-Q backscatter crosssection of the dielectric nose radome alone and in the presence of a circular plate at

80 the base of the radome (nose antenna). The antenna is inclined with respect to the plane perpendicular to the radome axis. The effect of the antenna on the RCS is pronounced at incidence angles close to specular directions of the antenna plate and this was to be expected. This is particularly seen in figure 4.7 for the backscatter RCS patterns perpendicular to the plane with respect to which the antenna is inclined. Specifically the large peak observed in figure 4.7 at q = 1500 when the antenna is inclined at 600 with respect to the x-z plane is due to the large specular return at that angle of incidence. 4.1.3 Summary A version of the fast multipole method was employed to compute the electromagnetic scattering from electrically large nose radome structures. The application of this technique while preserving to a great extent the accuracy of the moment method, significantly alleviates the limitations on the size of the bodies analyzed. In this work we analyzed nose radome shaped structures which were composed of metal or thin dielectric shells and in the case of the latter the resistive boundary condition was employed. The hybridization of a technique such as the finite element method will enable the treatment of a wider class of materials and could thus handle radomes with frequency selective surfaces. An application of such a hybrid technique is discussed in the next section.

81 y lx AV/ 7 7 77 /,V/ /V / / ///zzz 7/zzz z 1 777777ZZ ZZZZzz// 1.16 k (a) 0.85 X 20. 10 0 o U r> w -10 U o -20 ct c3 -30 -40 20 ^ 10 o O ) -20 30 -30 -- 9,"..,... - I... I ~ I ' '- R-0 (Thiswork) R=0O (Naor & Senior) - - - R=2.12-j0.2 (This work). v R=2.12-j0.2 (Naor & Senior) - / '\ - 'I\ i I I II 0 15 30 45 60 Observation angle 0 (deg.) (b). * -40 75 90 0 15 30 45 60 Observation angle e (deg.) (c) 75 L....~..-. 90 Figure 4.2: (a) Resistive plate geometry (b) Elevation plane backscatter RCS (0b polarization) (c) Elevation plane backscatter RCS (00 polarization)

82 Y (a) (b) J * *, *. I...... S E2 c/:).iM C4 E -ec cI.M 1 - coi FMM Radome code * Arvas & Ponnapalli (AP-T May 1989) ~ I * {}.........1 n L n - 0 30 60 90 120 150 Observation angle e (deg.) (c) 180 0 30 60 90 120 Observation angle H (deg.) (d) 150 180 Figure 4.3: (a) Von Karman radome geometry (2D section) - L = 1.1 m, D = 0.6 m (b) Meshed Von Karman radome (c) Bistatic RCS (i = 0~) - 99 polarization (d) Bistatic RCS (0i = 0~) - /s4 polarization

83 25 ' 25 CICERO - 0po ' a FMM Radome code - O0 pol ------ CICERO - pol 15 v FMM Radome codea. J pol 5 1 I | \\ W // x V 5\ v ---,,,I [ 8 I I I,I, v 0 15 30 45 60 75 90 Observation angle (deg.) Figure 4.4: Backscatter RCS for a 10A long nose radome of diameter 1A

84 I// \\\\ / / \ \ / / I \ I I \ \ I \ / i \ I l \ I I \ I - - - -' I 1 I i I I Y 0 ~Lx 6 (a) (b) 1-1 1-1 e~s C8 cn s ~s 6 it co cd IA a ci co 0m 15 10 5 0 -5 -10 -15 -20 -25 -30 -35 - Resistive radome only - E=4 - - Antenna - 0 deg inclination - - - Antenna - 30 deg inclination T Antenna - 60 deg inclination 0 Antenna - 90 deg inclination 0 30 60 90 120 150 180 Observation angle 0 (deg.) (c) Figure 4.5: (a) and (b) Nose antenna - radome geometry - 6 is the angle of antenna inclination w.r.t x-z plane (c) aoQ backscatter RCS (y-z plane)

85 X 'I \\ 'I \ 'I \\ / II I I I I I I I+ I Y x 5 (a) (b) I-6 C( Q, 9 15 10 5 0 -5 -10 -15 -20 -25 - Resistive radome only - er=4 - - Antenna - 0 deg inclination - - Antenna - 30 deg inclination Antenna - 60 deg inclination Antenna - 90 deg inclination 180 6 is the angle of antenna incli(y-z plane) -30 -35 0 30 60 90 120 150 Observation angle 0 (deg.) (c) Figure 4.6: (a) and (b) Nose antenna - radome geometry - nation w.r.t x-z plane (c) oaw backscatter RCS

86 / / I L4,x, / i\,\ / I / I I I I i I Y x ~ (a) (b) 15 10:, 5 1-/ ri 1 -10 I -15 4 -20 M -25 -30 -35 Resistive radome only - e=4:\ / v/ --- Antenna - 0 deg inclination --- Antenna- 30 deg inclination Antenna - 60 deg inclination Antenna - 90 deg inclination 150 180 0 30 60 90 120 Observation angle ) (deg.) (c) Figure 4.7: (a) and (b) Nose antenna - radome geometry - 6 is the angle of antenna inclination w.r.t x-z plane (c) oas backscatter RCS (x-y plane)

87 / I I I I! I I ht' ---.. - -' (a) (b) 15 | | | | 10 '..\_\. -- Resistive radomeonly.- e=4 0 Antenna - deg inclination 5 I \ — Antenna-30 deg inclination / / Antenna - 60 deg inclination -10 -o - \ / o Antenna - 90 deg inclination -15. -20 m -25 -30 -35......... 0 30 60 90 120 150 180 Observation angle ~ (deg.) (c) Figure 4.8: (a) and (b) Nose antenna - radome geometry - is the angle of antenna inclination w.r.t x-z plane (c),09 backscatter RCS (x-y plane)

88 4.2 Application of fast algorithms to hybrid FE-BI solutions The Finite Element-Boundary Integral (FE —BI) method reviewed in Chapter II has been quite popular and extensively applied to many scattering and radiation problems. The method [49] combines the geometrical adaptability and material generality of the FEM with the rigorous boundary integral (BI) for truncating the mesh. Nevertheless, although "exact", the FE-BI leads to a partly full and partly sparse system which is computationally intensive for large boundary integrals. When the boundary is rectangular or circular, the FFT can be used to reduce the memory and CPU requirements down to N log N [49] [50]. However, in general, the boundary integral is not convolutional and in that case the CPU requirements will be O(Nb), where Nb denotes the unknowns on the boundary. The application of the fast multipole method enables the computation of the boundary matrix-vector product using O(NbJ5) operations per iteration [51]. Two-level schemes can also be employed to reduce the operation count down to O(Nb'33) [552]. In this section, we apply the three different versions of the Fast Multipole Method (FMM), introduced in section 3.1, to reduce the storage and computational requirements of the boundary integral when the latter is used to terminate the finite element mesh. By virtue of its low operation count, the application of the FMM results in substantial speed-up of the boundary integral portion of the code independent of the boundary shape. Each version of the FMM is associated with inherent approximations and a goal of this chapter is the evaluation of these approaches in terms of CPU requirements and accuracy.

89 4.2.1 Problem Definition As an application of the proposed study we consider the scattering by a cavitybacked aperture shown in Figure 2.2. The FE-BI formulation for this problem was already outlined in section 2.2 and results in the system [Ant] [A] {f } 0 (48) [Ax] [B] {_\} {,inc} whose typical form is given in Figure 4.10. For H,-incidence, the vector {(} represents the magnetic field at the nodes within the groove and on the boundary whereas {4} is proportional to the magnetic current (or tangential electric field) on the boundary cells. By virtue of the finite element method, the matrices [AInt] and [Ax] are sparse and thus the corresponding matrix-vector products are implemented using O(N) operations. However [B] is a full sub-matrix and thus O(Nb2) operations are needed to perform the product [B]{b} with Nb denoting the number of nodes on the cavity aperture. Consequently, in an iterative solution, this matrix-vector product becomes the computational bottleneck. To reduce the operation count we can employ the FMM procedure for implementing the products [B]{4}. We employ the three different versions of the FMM discussed in Chapter III to accelerate the boundary integral matrix-vector product computation. 4.2.2 Results A computer code based on the three FMM formulations and utilizing the conjugate gradient solver was implemented and executed on a HP 9000/750 workstation with a peak flop rate of 23.7 MFLOPS. The geometry analyzed was a rectangular groove of depth 0.35A with a material filling of 6r = 4 and [r,=1. It was illuminated by a plane wave at normal incidence and BI subsystems of different sizes were

90 J i M. A i Y // ( go4,o) (Boundary/aperture) d j i.|(;r ii._pec x w FEM domain Figure 4.9: Geometry of the groove recessed in a ground plane achieved by considering three different groove widths: 25A, 35A and 50A. Although the boundary integral unknowns were about a seventh of the total unknown count, the number of boundary integral multiplications was more than 90% of the total number of multiplications for performing the entire matrix-vector product. This is because the bandwidth of the FEM subsystem was only about five, while that of the fully populated BI subsystem was equal to the number of BI unknowns. Thus, the overall computational requirement is dominated by the BI, hence necessitating use of fast integral algorithms like FMM to alleviate it. Table 4.1 compares the execution time and RMS error [24] of the standard FE-BI to the FE-FMMExact, FE-FAFFA and the FE-Windowed FMM (FE-WFMM). The data reveal that the FE-FMMExact offers almost a 50% savings in execution time with almost no compromise in accuracy. While the FE-FAFFA is the fastest of the three algorithms, the RMS error was substantially higher (>1 dB). If the maximum tolerable RMS error 1 is set at 1 dB 'To our experience, a 1 dB RMS error criterion for patterns with large dynamic range as encountered in these simulations gives a calculated pattern that is in excellent agreement with the exact. Typically, a 3 dB error criterion can lead to large deviations between calculated and exact data.

91 0 a X o. %. L System 80 5... gi B:L: 200 60 80 100 120. column Figure 4.10: Typical form of the FE-BI ystem matrix arising from the scattering/radiation 80 t \ l 0 20 40 60 80 100 120 column problem of a groove in a ground plane. [24], the FE-WFMM is the most attractive option since it meets the error criterion and is only slightly slower than the FE-FAFFA. The BI memory requirement for the FE-WFMM is higher than that for the FE-FMMExact but lower than that for the FE-FAFFA. This stems from the fact that the FE-FMMExact provides the best representation for the matrix elements which lie electrically far apart, thus enabling the use of a lower threshold distance and consequently a less populated near zone interaction matrix. The solution CPU time for the FE-WFMM competes favorably with a hybrid algorithm obtained by hybridizing the Finite Element Method with the Adaptive Integral Method [29] (FE-AIM) as seen from the data in Table 4.2 for reference purposes. Table 4.2 also shows execution times for the Finite Element - Artificial Absorber (FE-AA) [53] technique. The latter uses a fictitious artificial

92 absorber to terminate the Finite Element Mesh and leads to a completely sparse system. Consequently the CPU time per iteration is smaller than either the FE-BI or its variants. However, due to the larger system condition, a greater number of iterations was necessary to achieve convergence (For the linear system Ax = b, convergence is defined as LL < 10-4, where r is the residual vector). Thus, the total time required lIbl " ' for convergence is greater than the standard FE-BI itself. The advantages of a fast integral mesh truncation over the standard FE-BI and an absorber termination are summarized in Table 4.3. Of interest is the comparison of the residual error as a function of the number of iterations in the CG solver. Such a comparison is shown in Figure 4.11 and is seen that the curves for the FE-BI and the FE-FMMExact overlap to graphical accuracy whereas the FE-WFMM shows a very small deviation from the exact result. Thus, the hybridization of the FMM does not have any adverse effect on the condition of the FE-BI system. Solution CPU Time (Minutes, seconds) Groove Total BI Width Unknowns Unknowns FE-BI (CG) FE-FAFFA FE-FMMExact FE-WFMM 25A 2631 375 (9,33) (4,21) (6,20) (5,08) 35A 3681 525 (17,33) (6,56) (11,32) (8,23) 50A 5256 750 (47,21) (16,33) (28,20) (18,12) I /. - % RMS error (dB) Groove FE- FE- FEWidth FAFFA FMMExact WFMM 25A 1.12 0.0752 0.6218 35A 1.2 0.1058 0.721 50A 1.36 0.1123 0.843 Storage of BI (KB)______ Groove FE- FE- FE- FEWidth BI (CG) FAFFA FMMExact WFMM 25A 1072 277 166 185 35A 2102 458 275 381 50A 4291 784 470 613 Table 4.1: CPU Times, RMS error and Storage of the hybrid algorithms

93 10 10-1 o 0 a) ' 10-2 "o a) 10-3 -4 0 100 200 300 400 500 Iteration number 600 Figure 4.11: Convergence curves for the hybrid algorithms for the groove of width 25A Solution CPU Time (Minutes,seconds) Groove Width FE-AA FE-AIM 25A (12,1) (4,48) 35A (18,56) (9,27) 50A (49,52) (25,20) Table 4.2: CPU Times of FE-AA and FE-AIM algorithms Truncation scheme BI Fast BI Absorber Mesh termination Can be on Can be on Typically 0.3A surface the body the body away Storage O(N2) O(N1'5) or O(N1 33) O(N) Time/iteration Highest Lower than BI Smallest Total solution Typically less Smallest Highest time than Absorber Table 4.3: Comparison of mesh truncation schemes

94 Algorithm Operation count for BI Operations (multiplications) computation (as implemented) required for BI computation for the 50A groove FE-BINb 562500 FE-FAFFA (NNG + 2)Nb 5 + (1 - 2NNG)Nb - NNGN 136890 FE-FMMExact (NNG + 6)NA-5 - 2NNGNb 180356 FE-WFMM (4 + NNG + QwIN)Nb' - NNGQWINNbU667 153780 Table4.4: Exact BI operation count of the hybrid algorithms. Nb is the number of BI unknowns. NNG is the number of near groups (groups which are treated with the exact moment method). QWIN is the width of the window in the windowed FMM. Table 4.4 gives the exact BI operation count rather than merely stating its order. The knowledge of the constants associated with each exponent of Nb enables us to compare the requirements of two algorithms which might have the same order of operation count. In Table 4.4, NNG is the number of near groups (groups which are treated with the exact moment method procedure owing to their electrical proximity) which depends on the algorithm and the problem geometry. NNG is smallest for the FE-FMMExact and largest for the FE-FAFFA, due to the use of the far-zone Green's function in the latter. Table 4.4 also gives the number of multiplications in a single BI matrix-vector product for the 50A groove. For the FE-FMMEXact the number of multiplications is reduced by a factor of three over the FE-BI. However, the actual CPU time is reduced by a smaller factor due to computational overhead for the various function calls. The performance of the hybrid algorithms at a more stressing angle of incidence is depicted in Figure 4.12. We observe that the RMS error follows the same trend as for normal incidence illumination. For this example calculation the width of the groove was 10A. However, even for this smaller size aperture the scalability of the speed-up is maintained. The employed near group radius was 1A and thus the matrix vector products for groups separated by a distance less than a wavelength was computed

95 using the exact method of moments procedure. Smaller near-group distances can be employed to reduce the CPU time even further and near group distances down to 0.3A have been found to yield sufficiently accurate results. 30. -— i —. - FE -BI DMin = 1 x 20 0 FE - Exact FMM 0 FE- Windowed FMM O ~ FE-FAFFA CA 0: -lo. i_ pec -20 w=10 X3 _________________________ (a) 0 30 60 90 120 150 180 Observation angle (deg.) (b), I o I,. I --- Groove Width FE-FAFFA 10A 0.7627 RMS error (dB) FE-Exact FMM FE-WFMM 0.1621 0.3291 I -I (c) Figure 4.12: Scalability of the hybrid techniques to smaller problems (a) Problem geometry (b) Bistatic patterns (c) Error table 4.2.3 Summary We examined the computational performance of the finite element-boundary integral method using the fast multipole method for evaluating the BI matrix-vector products. Three different versions of the fast multipole method were used and it was shown that a considerable reduction in CPU time could be achieved with even further reductions for larger problem sizes. The FE-WFMM provided the best compromise

96 in terms of speed and accuracy.

CHAPTER V Applications of the Adaptive Integral Method (AIM) 5.1 Scattering from Plates Containing Small Features Using the Adaptive Integral Method (AIM) The Adaptive Integral Method is another algorithm which reduces the computationally complexity of moment method solutions. In the case of AIM, the CPU reduction is achieved by mapping the original MM discretization onto a rectangular grid and exploiting the Toeplitz property of the Green's function on this grid. That is, the Fast Fourier Transform (FFT) is invoked to compute the matrix-vector products in the iterative solver. For an arbitrary three dimensional body, a three dimensional FFT is required and as can be understood, this calculation is very time consuming. For planar scatterers the dimensionality of the FFT is reduced by one, thereby significantly accelerating the solution. In this chapter, we examine the benefits of AIM when the body is not electrically large, but is highly tessellated owing to its intricate construction, thus leading to a large unknown count. We show that significant savings in CPU and memory can be achieved by AIM and examine its accuracy for near field and far field computations. 97

98 5.1.1 AIM for Planar Scatterers In this section, we describe the application of AIM to planar scatterers. Following the discretization procedure outlined in Section 2.3, we begin with the linear system [Z]{I} = {V} (5.1) with [Z] being the elements interaction matrix, whereas {I} is the vector of the unknown coefficients and {V} is the excitation vector. The matrix [Z] is fully populated, demanding 0(N2) storage, and each [Z]{I} matrix-vector product requires O(N2) multiplications. Fast algorithms such as FMM and AIM are used to reduce the operation count from N2 down to NV, where a < 1.5. Both algorithms work on approximating the far zone interactions. In the case of AIM, the CPU reduction is achieved by first splitting the matrix as [Z] = [near] + [zfar] (5.2) based on a threshold distance referred to as the near-zone radius. The matrix [Znear] contains the interactions between elements separated less than the threshold distance, whereas [Zfar] contains the remaining interactions. The elements of [Znear] are evaluated with the exact MM while those of [Zfar] and the product [Zfar]{I} are evaluated in an approximate manner as prescribed by the AIM procedure [29]. Application of AIM requires that the whole geometry be enclosed in a regular rectangular grid. Basically, the fields of each interior edge is re-expressed using a new expansion based on delta sources located at the nodes of the uniform AIM grid as depicted in Figure 5.1. For the mth edge, this new expansion has the form M2 m = Z ^ -(x - Xmq)(y - ymq)[Aq x + AyqY] (5.3) q=l

99 where rq are the position vectors of M2 points on the square surrounding the center of the edge and 6(x) is the usual Dirac delta function. The coefficients AxY are suitably chosen so that the new expansion is equivalent to the original representation using triangular elements. A similar expansion is used for the divergence of the basis functions M2 E=m (x- Xmq)6(Y- Ymq)Ahq (5.4) q=1 To find a relation between the AxY and In coefficients, we equate moments of the two expansions up to order M. Specifically, we set M, = Fm (5.5) ql,q2 qljq2 where OO r00 Mm = / m(x - x)q(y - ya)q2dxdy for 0 < ql,q2 < M * -00 -00 M2 = (xmq - a)q1(Ymq - ya)2 [Ax + Amqy] with q = ql + q2 (5.6) q=1 r00 r00 Fmq2 = J fm(x - xy)(y - y)2dxdy (5.7) ql q2 -o -oo Similarly, by equating moments of V,. J3 with the new expansion (5.4), we establish a relation between Ad and In. That is, we set Dl2= Hl,q2 (5.8) where roo o M2 Dlq2 = d (x-a )(y -a y)2dxdy = (xmq- xa ) (mq- ya)2Adq (5.9) qj q2 Pm (l q2mq o o 'q=l Hq2 = J Vs * fm(X -a)q(y - )dxdy (5.10) (5.5) and (5.8) give three M2 x M2 systems yielding the equivalence coefficients as the solution. This process is depicted pictorially in Figure 5.1.

100 Original MM discretization AIM grid (section) Am2xAm1 Y. Y *< *Iz Z —X * * _A 5m 9 A m7 t, t AIM Representation for edge * * "in" (Note: Each of the * * / g * * * * original basis functions is *;/ * * defined on triangle pairs) Original MM discretization 2 M =3 > AIM coefficients = M =9 Figure 5.1: The process of transformation from the original MM grid onto the AIM grid Were we to use the equivalent expansions to represent the currents everywhere, the resulting impedance matrix will be of the form 3 [ IM = [A]i[][A] (5.11) i=1 In this, [A]i are the sparse matrices containing the coefficients of the expansion (5.3 ) and (5.4) whereas [G] is the Toeplitz matrix whose elements are the free space Green's function evaluated at the grid points. It has been shown in [29] that [ZJAIM] is not of sufficient accuracy for modeling the interactions between the nearby current elements. To take advantage of the Toeplitz structure of [G] and sparsity of [A] we can still use [Ztoa] to represent the far element interactions. However, we will retain the exact interaction matrix elements for the near element interactions. That is, we rewrite [Z"AIM] as [zitotal vnear ar 512 [ JAIM [Z]AIM + [Z]A (5.12) Comparing this to (5.2) and setting [Z]far [Z]M we can rewrite the original [Z] matrix as [z] - ([Z]r - [Z]a) + [Z]total [Z '~ ([Z] — ZAiM) — [ZAiM (5.13)

101 or 3 [Z] -_ [S] + i[A] i[G][A]T (5.14) i=l where [S] = [Z]nea- [Z]ne4 is a sparse matrix corresponding to the difference between the near field interactions computed by MM and AIM. The Toeplitz property of the Green's function, defined on the regular grid, enables use of the FFT to accelerate the computation of the matrix-vector product. The sequence of operations involved in the construction of the coefficient and Green's function matrices are indicated in Figure 5.2(a); those for the matrix-vector product execution are outlined in Figure 5.2(b). In the computation of the matrix-vector product, the initial step of transforming the currents from the original MM grid onto the uniform AIM grid is comparable to the grouping operation of the FMM. While the FMM relies on grouping to reduce the number of scattering centers, the sequence of operations in AIM can be interpreted as a realignment of scattering centers onto a regular grid. Although, this process does not reduce the number of scattering centers, the regularity of their location enables use of the FFT for fast computation of matrixvector products. 5.1.2 Results When examining the merits of a fast integral algorithm such as AIM, of importance is the memory and CPU requirements, both contrasted to the delivered accuracy. Although approximate analytical expressions have been derived in [29] for some of these parameters, these refer to implementations involving cubical grids and the three-dimensional FFT. Our goal in this chapter is to assess the accuracy of AIM in treating small details within an aperture/surface and to provide the reader with quantitative measures on the performance of AIM when implemented with the two dimensional FFT. The near-zone radius or threshold distance has a dramatic impact

102 Construct cartesian grid which extends over the object Calculate moments of the original basis functions (RWGs)................................................................................................... Solve for equivalence coefficients Calculate Toeplitz Green's function over the cartesian grid......................................................................................... Compute the near-field matrix (interactions between elements closer than threshold) (a) Compute I = A I (Performs transformation of current from MM grid to AIM grid) Compute 1= I (Compute the Fourier transform of the discrete current distribution)................................................................................................... Compute V = G I (Far-field Fourier Transform - G is the Green's function) Compute V =,l V (Inverse transform to AIM grid) far A Compute V = AV (Far field back to MM representation) Compute v 7 r 1 (Near field with exact MM).................................................................................................. Add to form total field far nvar V(b)= V +V (b) Figure 5.2: (a) Matrix build operations and (b) Matrix vector product computation in AIM on the CPU requirements since it controls the non-zero element population of the system matrix. In the case of AIM, because of the inherent mapping to a uniform grid, we are highly interested in examining its suitability to model small and fine details embedded in much larger scale structures.. The calculations for the plate con

103 figurations given next are intended to address this issue by examining the method's performance for a number of representative and practical situations. All of the included results were generated using single precision arithmetic on an HP9000/C-110 workstation with a rated peak speed of 47 Mflops (the level 4 optimization option was also used). In all cases, a third order (M=3) multipole expansion was used with a grid spacing of 0.05A. Figures 5.3-5.7 depict the 00 and q$ polarization radar cross section patterns (< = 0~ cut) as calculated by AIM for the different threshold distances indicated on the figures. The first circular plate has no holes and was used to validate the method. From the pattern comparisons, it is clear that AIM recovers the exact result very well. As given in Table 5.1 and 5.2, AIM achieves this with at least a factor of five less memory than the traditional MM. even though the geometries are still rather small to demonstrate the full impact of AIM. Also, Table 5.2 shows that a near zone radius of 0.3A is sufficient to maintain good accuracy (below one dB in RMS error [24]). The advantage of AIM is more pronounced when gaps are inserted into the plate's surfaces and this is the primary reason that one may prefer AIM over other fast integral methods for planar structures. As depicted in Figures 5.4 and 5.5, AIM maintains its accuracy for the same threshold criterion even though the gaps/slots have a dominant effect on the RCS pattern as shown in Figure 5.4. In the case of narrow slots (or thin ridges in the plates) of width 0.03A, the memory requirements of the traditional MM increase quickly due to the higher element density. For the geometry in Figure 5.6, AIM yields memory saving of 79% and the CPU time is reduced by a factor of 12 while retaining the monostatic pattern accuracy to within a tenth of a dB. This is achieved by using a uniform AIM grid density of 20 points

104 per linear wavelength even though the cell density of the original plate mesh is much greater due to the narrow slot. One may assume that this change in grid density will affect the near zone field. However, our observations indicate that the surface current is equally accurate. For the configuration in Figure 5.6 the average current density error is 7.3% for a threshold distance of 0.2A and 6% for a threshold distance of 0.4A. The currents for the geometry in Figure 5.6 along the center narrow strip are plotted and compared in Figure 5.8. These results demonstrate the important feature that the near zone threshold criterion is not affected by the specific geometrical details, leading to tremendous memory savings. Moreover, the accuracy of the results provide a convincing argument that AIM can efficiently handle highly irregular and resonant (i.e. antenna) geometries as well as smooth scatterers. At the same time, the convergence rate of the AIM system is unaffected indicating that the system condition is unchanged. This is of critical importance for fast iterative solutions, since an increase in the iteration count would annul the faster computation of the matrix-vector product. Figure 5.7 shows the monostatic RCS pattern for a grating structure which acts as a "polarization filter". The thin ridges in the grating cause a strong specular A A return for the 4) polarization (almost 10 dB above the return in the absence of the gratings) as is evident from the results in Figure 5.7(d). Of importance is that the MM triangular mesh in Figure 5.7 required a cell size of 0.02A per linear dimension because of the narrow grating. However, the overlaid rectangular AIM grid could be selected to have a much coarser discretization. More specifically, we chose grid spacings of 0.05A and 0.1A for the AIM grid and, thus, computational requirements of AIM were much lower. For the 0.1A grid spacing the solution time was reduced from 2.75 minutes down to only 12 sees at the expense of some accuracy (fraction of a dB).

105 To further increase in accuracy, we employed a 0.05A grid spacing and as shown in Figure 5.7(b) the AIM curve is now indistinguishable from the reference MM result (within 0.1 dB,). From Tables 5.1 and 5.2, the AIM computational and memory requirements are 8 times and 9 times less, respectively, without loss of accuracy. This is a significant observation and we have found that both the convergence rate and condition of the AIM system remains essentially unchanged from the original moment method system. The original discretization for the geometry in Figure 5.7 and the equivalent AIM grids are pictorially depicted in Figure 5.9. It should be noted that even though the size of the discretization is very small, retaining the selfcell term alone in the moment method system introduces huge error (Figure 5.10), thus emphasizing the importance of non-self terms. 5.1.3 Summary The performance of AIM is much improved when applied to scattering from flat complex scatterers and scatterers with high discretization rates. Thus, the reduction of solution time is considerably more for the geometries depicted in Figure 5.11(a) and 5.12(a) than for the geometries in Figure 5.11(b) and 5.12(b). A memory reduction of 5 to 10 times over traditional MM was observed without compromise in accuracy when using a threshold radius of 0.2A. This CPU reduction is achieved without resorting to parallelization or optimization techniques (as is known AIM is particularly amenable to such improvements). More importantly, the AIM algorithm is capable of modeling very small details in large bodies with a high degree of accuracy, while simultaneously saving considerable memory. This is of importance when modeling broadband antennas (spirals or log-periodics) and gratings which are both large in overall size but can contain features as small as A/100 in size. Application of AIM for analysis of cavity-backed antennas is described in the next section.

106 Discretization Geometry Facets Edges Unknowns MM memory (MB) MM solution time 99 pol (0 = 0~ inc.) Figure 5.3 586 908 850 5.51 32 secs Figure 5.4 554 890 772 4.54 29 secs Figure 5.5 1130 1806 1584 19.14 4 mins 50 sees Figure 5.6 1036 1667 1441 15.84 4 mins Figure 5.7 1038 1957 1157 10.21 2 mins 45 secs Table 5.1: Solution CPU time and memory requirement of the moment method AIM Data Geometry Threshold Non-Zeros Memory Solution time RMS Error(dB) (A) in Near Z (MB) 00 pol (0 = 0~ inc.) 00 pol |x pol 0.3 59928 0.68 23 secs 0.1718 0.0755 Figure 5.3 0.4 100182 1.14 25 secs 0.1490 0.0693 0.7 257390 2.94 28 secs 0.0728 0.0490 Figure 5.4 0.4 79030 0.9 21 secs 0.0728 0.0583 0.6 157994 1.8 27 secs 0.0721 0.0520 Figure 5.5 0.7 283774 3.24 3 mins 32 secs 0.8017 0.5185 Figure 5.6 0.2 296250 3.39 20 secs 0.1063 0.0949 0.4 649556 7.43 31 secs 0.0548 0.0632 Figure 5.7 0.2 120220 1.37 18 secs 0.0469 0.0469 Table 5.2: Solution CPU time, memory requirement and RMS error of AIM (all entries in this table were computed with an AIM grid spacing of 0.05A) 5.2 Computation of radiation and scattering from planar cavity-backed antennas with the Adaptive Integral Method (AIM) Several cavity-backed antennas contain small features and details which may necessitate high discretization. This could take the form of very narrow slots which may be a fiftieth or hundredth of a wavelength in width. Discretization of such geometries could lead to very large numerical systems even if the size of the antenna is not electrically very large. To efficiently treat such systems, the properties of an algorithm based on an iterative solver, should include the following

107 * It is of paramount importance that the "threshold" distance (distance beyond which interactions are treated as of the far zone variety) is as small as possible. * It should be capable of characterizing small perturbations in an otherwise smooth surface. * It should be capable of modeling near fields accurately. * If the algorithm incorporates a process by which very small discretization details can be "mapped" onto a different domain which is less dense than the original, computation of the matrix vector product in this domain would simulate the effect of a reduced number of unknowns. Figures 5.7 and 5.8 depict two planar configurations analyed by the AIM from which it can be gleaned that all the above criteria are met. Unlike FMM, which carries out the matrix vector product on the original moment method discretization, the ability of AIM to map the small details onto a sparse grid and still retain accuracy makes it the method of choice to analyse such antennas. For efficient modeling of the cavity we employ FEM with its low O(N) storage and execution time. Triangular prisms are used for discretization of the cavity volume for the reasons described in Section 2.4. 5.2.1 Implementation The FE-BI formulation for three dimensional cavity-backed antennas using prismatic elements is described in Section 2.4. Substitution of (2.60) in (2.58) gives a discretized boundary integral of the form in (2.46(). The near and far zone terms are treated as outlined in 5.1.1. The FEM matrix and the near zone interactions of AIM are stored in a sparse storage format, thus affecting significant savings in memory.

108 5.2.2 Results Figure 5.13 shows the radiation pattern for an annular slot computed in the elevation plane, A = 5~. The reference FE-BI solution [39] is contrasted with computations of BI using AIM (indicated as FE-AIM). It is seen that for this example, the threshold distance in AIM can be reduced to 0.25A without significant loss of accuracy. This enables the reduction of matrix entries stored in the near field portion by a factor of three resulting in a corresponding savings in memory as indicated in the tabulation of the near-zone non-zero entries. Figure 5.14 shows the radiation pattern for the same antenna in the X = 90~ elevation plane. The normal direction in this plane, reveals the characteristic separation between co-polarization and crosspolarization levels for the annular slot at observation angles close to normal in the elevation plane. From this figure, it is gleaned that the threshold distance in AIM can be reduced down to even 0.15A if an average error of a dB could be tolerated. From the computation of near-zone matrix entries, such a threshold would result in a factor of five saving in memory. Figure 5.15 shows a scattering cross-section for the same slot but at a frequency of 0.73 GHz (at which the antenna is electrically even smaller) instead of the previous 1 GHz. It should be noted that for a threshold of 0.4A (larger than the diameter of the BI contour) the near-zone and far-zone entries for AIM cancel each other in accordance with (5.13), thus yielding a very small error (0.00086 dB) in comparison to the FE-BI solution. A quantity of vital importance in antenna computations is input impedance. Figure 5.16 depicts the input impedance of a very narrow probe-fed annular slot, computed using FE-BI and FE-AIM. The probe is placed at y = 0. It is seen that evaluation of the boundary integral with AIM enables the reduction of the near-zone non-zeros by more than half. Computation of input impedance demands very high accuracy and the threshold distance was

109 held constant at 10.5 cm (corresponding to 0.35A at 1 GHz and 0.49A at 1.4 GHz - the corresponding diameter of the entire BI contour varying from 0.513A to 0.718A). While, Figures 5.13-5.16 demonstrate the ability of AIM to translate very fine details such as a narrow slot onto a coarser equivalent grid, Figure 5.17 and 5.18 indicate the importance of a low threshold distance in modeling cavity-backed antenna arrays. Figure 5.17 and 5.18 indicate that for an average error of less than a dB in scattering and radiation patterns it is possible to reduce the number of non-zeros in the near-zone part of the impedance matrix by a factor of six, resulting in substantial saving in memory. This is a consequence of employing a threshold distance of 10 cm, which is about a fifth of the cavity diameter. It is necessary to note that employing such a threshold distance results in a majority of the interactions between different slots being treated with the AIM procedure. This is of paramount importance in modeling antenna arrays and spiral antennas. While Figures 5.13-5.17 compare spatial domain FE-BI and FE-AIM solutions, Figure 5.19 compares the spatial domain FE-AIM solution with a spectral-domain FE-BI solution presented in [54] for the scattering by a cavity-backed patch antenna. 5.2.3 Summary AIM, with its low threshold distance, and ability to translate to an equivalent grid is capable of saving a significant amount of memory and solution time for bodies which are finely discretized even though they may not be electrically large. Its accuracy is preserved even while performing radiation computations thus making it the method of choice for analyzing antennas with intricate details.

110 Plate orientation reference YAz - O Z X 2 z -- --- A0 Polarization reference / tO Polarization X 25,15 a c) -25 1-5.1 -25 00 Polarization..... I... -I....I..... I.... I..... -. MM -E AIM - Thbreshold = 0.3/ 'A AIM - Threshold = 0.4. ABi - Threshold = 0.74,..... I..... I..... I..... I.... 25 20 15 U 0rs at Q -5 -5 -— MM - \ AIM - Threshold = 0.3X A AlM - Theshold = 0.4k * AIM - Threshold = 0.7X..... I I.. I..... I..... I..... -10 0 15 30 45 60 Observation Angle 0 (deg.) 75 90 0 15 30 45 60 Observation Angle 0 (deg.) 75 90 Figure 5.3: Monostatic RCS for a circular plate of diameter 2A; Comparison of the standard MM & AIM

111 Plate orientation reference y z -x 25:: 15 a-15 -15 -25 1.6 X (a) 00 Polarization 0 15 30 45 60 Observation Angle 0 (deg.) (b) 75 90 25 a15 la 5 -5 a -5.-15 -25 () Polarization - 1 MM I 0 15 30 45 60 75 90 0 15 30 45 60 75 90 Observation Angle u (deg.) z Observation Angle o (deg.) (c).. (d) Polarization referenc I Y e / - x Figure 5.4: Monostatic RCS for a circular plate of diameter 2A with three slots computed with standard MM & AIM (a) Geometry (b) Effect of the slots on the RCS (c) 00 polarization backscatter RCS for the plate with slots (d) H0 polarization backscatter RCS for the plate with slots

112 Plate orientation reference y 3. 2 "A 3 X z x 1.94X Polarization reference / y 00 Polarization 4 Polarization x 30..... ' --- —-....... 30 25 o ALM - Threshold 0.7?k 25 mm:E20 El AINI% - ~11reshold =0.7k. 20 015 -5 0 15 30 45 60 75 90 0 15 30 45 60 75 90,Observation Angle e (deg.) Observation Angle 0 (deg.) Figure 5.5: Monostatic RCS for a square plate of sidle 4A with three holes computed with standard MM & AIM

113 I1 Plate orientation reference Z X Z X o.o in Z 0.03? ( Smaller than the AIM grid size) - Polarization reference y of 00 Polarization 0 fi -10; -20 9 cl an a-30 -40 0 15 30 45 60 75 90 0 15 30 45 60 Observation Angle 0 (deg.) Observation Angle e (deg.) 75 90 Figure 5.6: Monostatic RCS for a circular plate of diameter 1A sampled at 0.03 A (smaller than the AIM grid spacing) due to the narrow center ridge

114 1 0.05A li/ lX 0.04 ) (a) Polarization reference 00 Polarization x )) Polarization AZ.................................... 0 I -5 Ie -10 cV -15 u, -20 -25 -30 -, — IMM (For Fig (a)) ~- - - lAIM - Grid spacing = 0.053, (For Fig'a)t - - - - -A - - AIM - Grid spacing = 0.1, (For Fig (a)J \ MM ForFig(b)) 15 10 m eo la L, C4 O cn 'o u a U! co m eu CO 5 0 -5.,... I..... I..... I........ I..... MM (For Fig (a)) -- - AIM - Grid spacing = (.05, (For Fig (a) - --- AIM - (Gid spacing = 0.1, (For Fig (a)) /\ \ / _ " \\ / —..... I..... I..... I..... I.. ~.. I,,.,. -10 -15 0 15 30 45 60 Observation Angle 0 (deg.) (c) 75 90 0 15 30 45 60 Observation Angle e (deg.) (d) 75 90 Figure 5.7: (a) Geometry and mesh of the grated plate (b) Geometry and mesh of the "Groove" plate without gratings (c) 00 polarization backscatter RCS computed by AIM and MM (d) bb polarization backscatter RCS computed by AIM and MM

115 Plate orientation reference YA Z X -- - - ik 0.6 '. Currents at the centroid of the shaded triangular patches on the narrow ridge 0.03 k ( Smaller than the AIM grid size) -- MM -* AIM - Threshold = 0.4 k a Cu.a co eo 0 f3 ^3 4.5 4.3 4.1 3.9 3.7 3.5 3.3 3.1 2.9 2.7 2.5 -0 x-4 K1O I I I I I-~ ---~ — - -- -- -- -~- - -~ ---- --- -~ I I I' I I ---— ~- ---- - -hh h1- --— r- — r - I I I' I -— r- -r I I I' I -— r -— r - ' —r --- —-— r --- —-— r -~~ ~ I I --- ------ I — - - - I- - I- - - - - - - - - - - - - -I I- - - - - - I I I t.3 - -0.2 -0.1 0 0.1 y in wavelengths 0.2 0.3 Figure 5.8: Electric currents (Solution coefficients) on the narrow ridge for the geometry of Figure 5.6

116 0.05 A I77 4 -4-. 0.1 x 0.1k 0 0 0 0 0 0 o,~1O o o O o o o o o o O o o O o o o o o O o o o o o o 0 'o 0 0 0 0 0 'o 0 0 0 00 C, 0 00 C, 0 0 0a 0 a 0 00 0 0 0.05 iA 0 00 O O O O O O O O O O O O o O O O O 0 0 0 0 0 0000000000000000000 O0 0 0 0 0 0 00a000000000000000000 0 0 0 0 0 000000000000000000000 0 0 0 0 0 000000000000000000000 0 0 0 0 0 000000000000000000000 o o o o o verlaidoooAIMoogridso o o o o o ooooooooooooooooooooo oo oo oo oo oo oo oo oo oo oo o o o o o o oo oo oo oo oo oo oo ooooo oo oo oo oo oo oo oo oo oo oo o o o o o o oo oo oo oo oo oo oo ooooo oo oo oo oo oo oo oo oo oo oo o o o o o o oo oo oo oo oo oo oo ooooo oo oo oo oo oo oo oo oo oo oo o o o o o o oo oo oo oo oo oo oo ooooo oo oo oo oo oo oo oo oo oo oo o o o o o o oo oo oo oo oo oo oo ooooo oo oo oo oo oo oo oo oo oo oo o o o o o o oo oo oo oo oo oo oo ooooo oo oo oo oo oo oo oo oo oo oo o w217o~whw~w~tao -w - 4 - ix 0.04A Figure 5.9: Original discretization and equivalent AIM grids for the geometry of Figure 5.7 Standard moment method (0 pol). --- —--- Standard moment method (O pol) A Self-cell only (0 pol) [Avg. Error = 36.81 dB] v Self-cell only (O pol) [Avg. Error = 31.96 dB] 10 "0i C)4 U C/) cn Q Q 0 4 --- 4 — 0 -10 -20 -30 -40 -50 0 15 30 45 60 75 90 Observation angle 0 (deg.) Figure 5.10: Error introduced by retaining only the self-cell interactions of the moment method

117 B -0.5 --1. -1.5 - ~4' -2 2 2 0 0 -2 2 -2 4 -4 y x (a) (b) Figure 5.11: (a) Flat and (b) Curved plate with equal side lengths and discretization rates, resulting in equal number of unknowns. While the moment method yields equal solution time for both geometries, AIM would accelerate the solution for the geometry in (a) considerably more than that for the geometry in (b) Both geometries contain 1681 nodes, 3200 elements, 4800 edges Gridded at 2/40 Gridded at 7110 x (a) x (b) Figure 5.12: (a) Geometry of a 1A square plate gridded at A/40 and (b) 4A square plate gridded at A/10. While the moment method results in equal solution times for both geometries since they have equal number of unknowns, AIM would accelerate the solution for the geometry in (a) considerably more than that for the geometry in (b) owing to the smaller FFT pad for the geometry in (a)

118 y X a = 12.35 cm b = 0.75 cm p = 7.7 cm Frequency = 1 GHz 0 FE-BI FE-AIM (0.53 A) FE-AIM (0.35 X) FE-AIM (0.25.) Non-zeros in near Z matrix 4656 4656 2245 1537 -5 X -10 IQ -15 en 9 -20 0 o i -25 Cu ( -30 N -35 -40 -45 I Elevation plane pattern at )0=5... I V -,..I.I....I.I....,.I....,.I.. I L L- L Average Error (dB) Co-pol X-pol 0.4402 0.2998 -5 0...................... -90 -75 -60 -45 -30 -15 0 1. -- FE-BI (Co-pol) (deg.) - - FE-BI (X-pol) FE-AIM (Threshold = 0.53A) (Co-pol) FE-AIM (Threshold = 0.53X) (X-pol) -- FE-AIM (Threshold = 0.35X) (Co-pol) - - FE-AIM (Threshold = 0.35X) (X-pol) FE-AIM (Threshold = 0.25X) (Co-pol) FE-AIM (Threshold = 0.25k) (X-pol) 5 30 45 60 75 90............. FE-AIM (0.35 2) FE-AIM (0.25 A) 0.6343 0.5334 Figure 5.13: Radiation pattern from an annular slot in the - = 0~ elevation plane

119 Non-zeros in near Z matrix FE-BI 4656 y Average Error (dB) Co-pol X-pol 1.0524 1.0373 0.6142 0.9105 a = 12.35 cm: b = 0.75 cm p =7.7 cm Frequency = 1 GHz Frequency = 1 GHz FE-AIM (0.15 X) FE-AIM (0.25 X) FE-AIM (0.35 X) 915 1537 2245 FE-AIM (0.15 ) FE-AIM (0.25 3) FE-AIM (0.35 3) 0.876 0.2506 (a) 7:$ 1 -m ~' 0 o.-.Cl ao -4 Ct - o 7 Z -I 0 i * I I * I.. I. I I I I I 0 -10 -20 - -30 -40 -50 -5 - - - -90 -75 -60 -45 -30 -15 0 15 30 45 60 75 90 Elevation angle 0 (deg.) - FE-BI (0 pol) (Co-pol) --- FE-BI (0 pol) (X-pol) --- FE-AIM (0 pol) (Co-pol) (Threshold = 0.15 X) - - - FE-AIM (0 pol) (X-pol) (Threshold = 0.15 X) * FE-AIM (0 pol) (Co-pol) (Threshold = 0.25 X) ~ FE-AIM () pol) (X-pol) (Threshold = 0.25 X) FE-AIM (0 pol) (Co-pol) (Threshold = 0.35 X) FE-AIM () pol) (X-pol) (Threshold = 0.35 A) (b) Figure 5.14: Radiation pattern from an annular slot in the b = 900 elevation plane

120 Non-zeros in near Z matrix FE-BI 4656./ 4...-. i/ ^ iii'':,,ii..... '.. i..'' x b a = 12.35 cm... b = 0.75 cm p = 7.7 cm Frequency = 0.73 GHz FE-AIM (0.25 X) FE-AIM (0.33 X) FE-AIM (0.4 X) 2184 3272 4656 FE-AIM (0.25 X) FE-AIM (0.33 X) FE-AIM (0.4 X) Average Error (dB) Co-pol 1.3124 1.0666 0.00086 (a) E 4 -T3 ~ a3 0 -10 -20 -30 -40 -50 -60 -70 -80 FE-AIM (O pol) (Threshold = 0.53 k) 0 15 30 45 0 15 30 45 60 75 90 Elevation angle 0 (deg.) (b) Figure 5.15: Bistatic scattering pattern from an annular slot; Normal incidence in the - = 00 plane and observation is in the f; = 900 elevation plane

121 Probe / y x a = 12.35 cm b = 0.75 cm p = 7.7 cm Near Z non-zeros BI: 4656 AIM: 2245 Near-zone threshold I held constant at 10.5 cm 80(4 60( 40( 120( - 204 -20v!! \ I II I II ~~~ ---+ - --- -, - /- -/ —7- ~ + -r- -+^ - -- - I I I I '1 I I I I I I I I 'l- - -" - I Ie I m i \ i i I,, I --- - -- ---- ------- -- -- I -- -----. - I I I I ~I I I I + — I I I --...... -— I ---. —. —.'.' —. — --------- J Ii I " I I I I I I I I I I I I 1 Ir 1 \1 1 -l - - - - - - - - - - - - - - - - -I I I I I I - I I I I I I -- -— I — ------— I- - ---- -----------—, I I I I I I I I I ~I I I I I \ I I I, I I I I I \ I 1 - ' I I I I Lines - FE-BI Dots - FE-AIM I -41 li I 19 0.95 1 1.05 Frequency (GHz) 1.0667 1.1333 1.2 1.25 1.3 1.35 1.4 1.1 1.15 1.2 1.25 Frequency (GHz) FE-BI 71.74 +j259.04 248.51+j422.88 708.59-jl70.31 238.47-j259.27 109.77-jl46.96 67.63 - j71.61 49.7 - j16.48 1.3 1.35 1.4 FE-AIM 71.46 +j268.51 258.03 +j453.73 710.39 - j269.71 206.96 - j266.05 95.48 -jl142.82 59.84- j 67.16 44.42-jl2.95 Figure 5.16: Input impedance of a very narrow annular slot computed FE-AIM with FE-BI and

122 Total Edges: 5529 Non-PEC: 983 Surface Edges: 2110 (PEC) 496 (Slot) 248 (Non-Pec) Inner radius of each slot = 7.325 cm Cavity diameter = 49.4 cm Cavity depth = 1.5 cm Cavity filling Er= 2 Slot width = 1.5 cm Frequency = 1 GHz (wavelength = 30 cm) 10 0 I I I. I. i.10 -20 o -30 I AIM (Non-zeros in near, = 5008) -BI (Non-zeros in near Z = 30876) Average error = 0.8 dB I. -40 -50so -60 I 1 -70 -80 I- — ~ I _e -L-J ~ II _C ~ _~ __ _L _' _ 0 15 30 45 60 Observation angle 0 (deg.) FE-BI --- FE-AIVI (10 cm theshold) 75 90 Figure 5.17: Bistatic RCS at normal incidence (- = 900 array computed with FE-BI and FE-AIM plane) from a cavity-backed slot

123 Total Edges: 5529 Non-PEC: 983 Surface Edges: 2110 (PEC) 496 (Slot) 248 (Non-Pec) Inner radius of each slot = 7.325 cm Cavity diameter = 49.4 cm Cavity depth = 1.5 cm Cavity filling %= 2 Slot width = 1.5 cm Frequency = 1 GHz (Lambda = 30cm) Co-pol (FE-BI) 0 -10 1-1 -20 u * -30 It -40 -50.. -. I.. I.. I I. I. I.. I,, I.. I.. I. -90 -75 -60 -45 -30 -15 0 15 30 45 60 75 Observation angle (deg.) Figure 5.18: Radiation from a cavity-backed slot array computed with FE-BI and FE-AIM in the < = 900 plane

124 z = Patch,;- /E,.- -\,,-,> '.s' @'......../.. " - x - '.''-....~.?......:...... -..,.:~.:?.S,. /::........-.,-...... -........... '-. *.,Lx'. %,,.:, 'X.,.\,;' f L 0.5 cn (a) E -3 (XQ (Xa ~) 0 -10 -20 -30 -40 -50 -60 1 2 3 4 5 Frequency (GHz) (b) Figure 5.19: (a) Geometry and surface discretization of a cavity-backed patch antenna (b) Monostatic RCS at normal incidence versus frequency - cavity filling has a ~r = 2.2- jO.002 and Ur = 1

CHAPTER VI Conclusions In this dissertation we developed and implemented fast, memory-saving algorithms for integral equations. For each algorithm, comparisons of memory requirements and execution time were made with "exact" techniques. A large portion of the dissertation was devoted to a discussion on fast integral algorithms employing iterative solvers - FMM and AIM and their performance in conjunction with FEM. It is necessary to elucidate the differences between them in order to make a decision on their suitability for scattering or radiation problems. Cognizance of similarities and differences between the two algorithms also suggests future areas of research. We also list accomplishments made during the course of the dissertation work. 6.1 Comparison between FMM and AIM methodologies * FMM and AIM are employed in conjunction with iterative solvers only. This follows from the fact that iterative solvers require only the product of the system matrix and a, trial vector rather than the explicit system matrix itself. * The fundamental principle behind the two algorithms is the same. Both algorithms work on the basis that interactions between far away source and test) elements can be treated using approximations. To preserve accuracy, the inter 125

126 actions between elements electrically close to each other continue to be treated exactly, requiring careful evaluation of numerical singularities when the source and testing points collocate. It is necessary to note that discarding interactions between elements which are electrically far apart will lead to large and significant errors. * Only the near zone interactions are generated explicitly by both the techniques. These interactions are stored in a sparse storage format, thus effecting significant savings in memory. Storage required for calculation of terms used in the matrix-vector product for elements which are electrically far apart is negligible in comparison to the 0(N2) storage required by the moment method. This is an aspect of crucial importance since most commercial solvers for large-scale applications [55], [56] do not save memory. * FMM uses a series expansion for the Green's function to uncouple source and test points and enable grouping of unknowns. It does not affect the basis functions in any way. Grouping of unknowns effectively reduces scattering centers on the body. AIM, however, makes use of a new delta basis and forces equality of moments between the original and new basis. The new delta basis simulate the effect of the original basis and source and test points are not decoupled. A succinct way of restating the above is that FMM reduces the operation count (multiplications in the matrix-vector product of the iterative solution) by effectively reducing the scattering centers on the body. AIM does not reduce the scattering centers but regularizes the scattering centers by translating them to the regular grid enabling the use of the FFT, thus reducing the operation count. The regular AIM grid can be less dlense than the original moment method grid,

127 thus reducing the CPU time and memory even further. An analogy that ilustrates the difference between FMM and AIM is the process of counting beads. If the solution process is compared to the counting process, two possible techniques can accelerate the solution process - either grouping the beads with equal number of beads in each group or lining them up in a regular matrix format with rows and columns. The process of grouping and regular arrangement is akin to the calculation of moments, multipole expansion, discrete plane wave expansion etc. in FMM and AIM. However, once this is done the counting process (analogous to solution) is accelerated as only the number of groups needs to be counted (as in FMM) or the number of rows and columns need to be counted (as in AIM) thus precluding the necessity to count each bead (as in standard moment method). * Accuracy in FMM is controlled by the number of terms in the series expansion for the Green's function. In AIM, accuracy is controlled by the order of moments. Theoretically, both of these quantities should be infinite but are truncated for numerical purposes. * The AIM procedure of employing equivalent basis functions results in smaller values for threshold distance (electrical distance beyond which interactions between elements are treated approximately) in comparison to the FMM process of decoupling source and test interactions. However, since the basis functions are unaffected by the FMM procedure it is possible to take advantage of this as detailed in the later subsection on future work. * Interactions between elements which are electrically close and far from each other cannot be separated a priori in AIM. The Green's function matrix is

128 Toeplitz only if the entire structure is considered. Hence, it becomes necessary, while considering elements within the threshold distance, to find the difference between interactions computed with moment method and AIM as indicated in (5.13). FMM computes these interactions with just the standard moment method. 6.2 Contributions of this dissertation The goals of this dissertation were a better understanding and comparative evaluation of fast algorithms for integral equation solutions, their implementation in hybrid systems, analysis of error, CPU and memory requirement with an objective of assessing their effectiveness. Below, we list our contributions in detail and indicate the relevant refereed journal publications. Chapter 2 described an algorithm for generating sparse matrices from the original moment method linear system. This technique yields high degree of sparsity for smooth bodies but the sparsity is highly unpredictable for arbitrary scatterers suggesting the need for algorithms based on iterative solvers. Chapter 3 gives a complete description and understanding of the FMM as applied to electromagnetics problems [28]. Based on an error study we were the first to suggest parameters for the FMM such as threshold distance and sampling rate [24]. Chapter 4 describes two applications of the FMM. The overwhelming computational requirement of the moment method had precluded the simulation of realistic size nose radome geometries, and we were able to alleviate this to a considerable degree using a variant of the FMM. For the first time, FMM was used to reduce the computational requirement of the boundary integral in FE-BI solutions [26]. This makes the BI termination (dense system) computationally competitive with a termi

129 nation like an artificial absorber (sparse system) while preserving the accuracy and rigor associated with the BI termination [27]. Chapter 5 describes special features of AIM [57]. We demonstrated that this technique is particularly amenable to model planar bodies, especially those which might require high discretization owing to their intricate construction [58]. For the first time, we showed that the regular AIM grid can be less dense than the original moment method grid thus providing the effect of a reduced number of unknowns in the computation of the matrix-vector product. Chapter 5 also describes the performance of AIM for radiation simulations by considering a cavity-backed slot antenna. To our knowledge, this is the most efficient technique for treating narrow cavity-backed slots. Appendix 1 [59] describes a special technique for the discretization and solution of an integral equation for modeling inhomogeneous dielectrics. This integral equation reduces the number of unknowns at each volume location from six to three. This reduction is at the cost of a higher degree of singularity in the kernel which is treated with a special discretization procedure referred to as the modified Galerkin's technique. Appendix 2 [60] describes a simple, accelerated, quasi-analytical solution for an integral equation when the problem is electrically small. Using a single-basis representation of the equivalent current, we compute the scattering from a narrow groove in an impedance plane, precluding the need for setting up a numerical system of equations. 6.3 Recommendations for applicability Having developed and examined the parameters that control the accuracy, speed and memory requirement of various solution algorithms it is of paramount importance

130 to recommend a category of applications for which they are most suitable. An integral equation-based FMM formulation is most suitable for simulating electrically large conducting and thin dielectric scatterers with normal discretization (normal discretization rates for moment method analysis of scattering problems is 8 to 10 points per linear wavelength). A hybrid FE-BI formulation incorporating FMM combines the twin virtues of the FEM with its ability to treat inhomogeneous material sections, low 0(N) storage and computational requirement and O(N1 33) or lower computational requirement of FMM while evaluating the boundary integral. The adaptability of the FEM to inhomogeneous material sections removes the material restrictions imposed by a stand-alone integral equation formulation. This technique is thus suitable for treating electrically large composite material bodies efficiently. The AIM is particularly effective in modeling planar or nearly planar scatterers, since the dimensionality of the FFT is reduced to two from three, which would be required for an arbitrarily shaped three dimensional scatterer. The ability to translate to an equivalent grid, which can be less dense than the original discretization makes this algorithm particularly attractive for modeling highly discretized scatterers. Such discretization is necessary for intricate details on the scatterer. Threshold distance becomes a very important criterion when analyzing cavitybacked antennas such as slots and spirals since they may not be electrically large but generally require very high discretization. AIM, with its small threshold distance, is particularly attractive for evaluating the BI matrix-vector products, in FE-BI formulations.

131 6.4 Future work Fast algorithms have rejuvenated integral equation analysis in electromagnetics and have rendered several previously unattractive solution techniques more tractable. Here we list possible future areas of research, classified from the mundane to the esoteric. A couple of these suggestions could lead to future doctoral dissertations. Suggestions are given both for improving the techniques discussed, as well as applications of the currently existing techniques. 6.4.1 Applications * Applications in inverse scattering: [61] and [62] are excellent references which deal with ultrasonic focusing through inhomnogeneous media by inverse scattering. During the application of this technique, the forward scattering problem involving the 2D Green's function needs to be solved. The techniques outlined in Chapter 3 can considerably enhance the speed of these simulations. * Scattering and radiation from cavity-backed antennas on curved platforms: [63] and [64] describe a FE-BI solution for analyzing cavity-backed antennas on curved platforms. It employs the Green's function of a cylindrical body. For large-scale simulations or to enhance the speed of existing simulations the techniques outlined in Chapter 3 and 5 can be used. 6.4.2 Improved techniques * Faster, more robust iterative solvers: In this dissertation, we used the conjugate gradient (CG) and biconjugate gradient solvers (BiCG) while solving the linear systems iteratively. Recently, other solvers like the GMRES and Flexible GMRES [6)5] have been under investigation which might perform better than

132 CG and BiCG. GMRES, like CG and BiCG, requires a matrix-vector product, the computation of which can be accelerated by using the techniques described in this dissertation. Thus, these solvers could be very easily ported into the computer programs which implement the solution techniques described in this dissertation. * Improvements in the AIM algorithm: The order of moments in the AIM algorithm is inversely proportional to the grid spacing. A complete investigation of AIM moments could lead to a more efficient application of AIM. In our simulations, we used a radix 2 FFT which requires rounding off the size of the FFT pad to the nearest power of two. This rounding off can be wasteful especially with increasing FFT pad size. Use of FFTs with a different radix could lead to a more efficient AIM algorithm. The AIM grid extends over the entire dimensions of the body being analyzed. This may prove wasteful if the body is slender in two perpendicular directions like an aeroplane. Domain decomposition could be used to reduce the portion of grid not enclosing a part of the body. * Use of wavelets as basis functions in FMM solutions: The fast multipole method achieves its CPU reduction by employing a multipole expansion for the Green's function. It does not affect the basis functions in any way. Employing wavelet bases could lead to a two-fold CPU reduction. Use of wavelets produces the sparse moment-method matrix and concurrent use of the FMM results in further speed-up of the matrix vector products. Thus, use of multiresolution basis along with the FIMM would produce an operation count lower than that associated with a stand-alone version of either of these techniques. The FMM has been shown to be superior to wavelets [66] with regards to memory and CPU time,

133 but in this hybrid scheme, the wavelet basis would simply complement the FMM.

APPENDICES 134

135 APPENDIX A Solution of a new reduced-unknown integral equation for modeling inhomogeneous dielectrics A.1 Introduction The modeling of inhomogeneous dielectrics via integral equation methods has traditionally been accomplished by introducing equivalent volume electric and magnetic currents as shown in Figure A.1. For a dielectric with non-trivial permittivity and permeability this type of modeling implies six scalar unknowns at each volume location. As a result, the implementation of the resulting integral equation is computationally intensive and has excessive storage requirements. Recently [31], it was demonstrated that any inhomogeneous dielectric material, regardless of its permittivity and permeability profile can be modeled by a single electric or magnetic current density. Alternatively, either the electric or magnetic fields within the dielectric can be used as the unknown quantities. However, the simplest form of this reduced-unknown integral equation involves derivatives of the unknown quantities and, thus, a higher basis function is required for discretizing the resulting integral equation. Hence linear (or higher order) shape functions must be employed for the field expansion.

136 A.2 Discretization A.2.1 The Integral equation The volume electric-field integral equation to be discretized is of the form E'(r) = E(r)-/J [VGo(r, r') x I] * { V' x E(r') +V X [(er- 1)E(r')] } d (A.1) where I = xx + yy + zI, Go(r, r') is the free-space Green's function, Ei is the incident field and Vd is the dielectric volume. From the properties of the identity dyad, equation (A.1) becomes E'(r)= E(r) + J VGo(r, r') x { V x E(r') +V' x [(er- 1)E(r')] dv' (A.2) The second term in (A.2) can be written as V> x [(r - 1)E(r')] = (Er - 1)V x E(r') + V'er x E(r') (A.3) giving Ei(r) E(r) + e {re 7r —} JJJ VGo(r, r') x [V x Ee(r')] dv' e==l Pre N + E ///v VGo(r, r') x [V'Ere x Ee(r')] dv' (A.4) e=-1 In this iPre and Ere are the relative permeability and permittivity within the eth cell, Ve is the volume of the eth cell, N denotes the number of cells used in the discretization of Vd and the subscript e indicates the field in the eth cell. Note that in deriving (A.4) from (A.2) we replaced the constitutive parameters over the cthe cell by their average value (Erel're) [67]. We further remark that the second volume integral

137 reduces to a surface integral at the interface of p)iecewise homogeneous regions. For the application considered here the second integral of (A.4) vanishes and will be dropped hereon. For simplicity, we choose to discretize the volume Vd as a collection of rectangular brick elements. Within each brick element the fields are represented using the edgebased expansion functions borrowed from traditional finite element analysis methods [47], [68]. Edge-based expansion functions maintain field continuity across dielectric interfaces and avoid explicit specification of the fields at corners. Referring to the brick element shown in Figure A.2, we introduce the expansion (using local element coordinates) 4 4 Ee(r') = xEW (y', z' ) + y W (z' x')q j=l j=l 4 +Z S We (x', y')Oe$ j=1 = {}{W } (A.5) The explicit expressions for the shape functions are [69] (b-y')(t- z') _ = y'(t- z') (b -y W') ' y'z' bt bt (t - z')(a - x') z'(a - x') Wyl1 - Wy2 -ta ta ta (A.6) =(t - z')x' ' zx 3 ta 4 ta e (a-x')(b-y') - x'(b-y') wz3 = a-w = ab z ab We shall hereon refer to the shape functions by we, where 4 4 4 We(r) = xEWxj (r) +y Wy (r) +z Wz (r) j=l j=l j=l 12 =E we (A.7) n=l

138 To construct a linear system of equation from (A.4), we follow the Galerkin's procedure, wherein the weighting functions are the same as the shape functions in form, but the domain of integration for the testing functions may be arbitrary. Specifically, by testing (A.4) with the mth shape function of the pth element, we have JJ WP(r). E(r)dv = [P (r) Ee(r')dv 6(p-e) + c — Jr J lre JJfvpWP (r) JJJf VGo(r,r') x [V'x Ee(r')] dv'dv (A.8) where 6(p - e) is the Kronecker delta function, and WP are defined in (A.7). After assembly, equation (A.8) gives rise to a system of equations which can be written as {Zmn}{n} = {Vm} The integrals which must be evaluated in constructing this system are of the form VP = p(r) W E(r) dv ZP = WP (r) W. (r') dv] (p- q)+ eq - u, [///Up ] { 1rq fv WP (r)' ff VGo(r r') x [V x W (r')] dv' (A.9) where the superscripts refer to the brick element number and the subscripts are associated with the (local) element edges. For the self-cells (p = q) the kernel of the integrand in equation (A.9) has a - singularity which must be carefully evaluated. To do so, we rewrite the corresponding impedance matrix elements (by adding and subtracting an appropriate singular term) as ZPq = [ WP(r) WW(r') dv &(p- ) + E - lrq [JJJ~] { '}

139 JJ/q (' x Wq(r')) JIJ [W (r)-WP (r')] x VGo(r, r')dv dv' - rq irq {'r - Iq1 }fffs (v' xW (r')) {W (r) x V ' vf Go(r,r')dv} dv' (A.10) t (A.10) Clearly, the integrand of the second integral in (A.10) is now well behaved when p = q. The singularity has been transferred to the last integral which involves only the scalar Green's function as the integrand. We observe that on choosing Vp to be a small spherical volume, the last integral in (A.10) over Vp can be evaluated in closed form. Provided we consistently choose Vp to be spheres, this testing procedure does not introduce any approximations even though the discretization cells are maintained as rectangular bricks. The described testing procedure was employed by Tsai et al.[70] who also selected Vp to be a sphere having a volume equal to the pth rectangular brick cell. We perform an analytical integration of the free space Green's function over the spherical volume to get fff, r) dv = k- e d (s o - dcos(kod)) + sin(kod)e-jk ( + ja) - sin(kod)e-k~d ( +d) } (A. 1) where a is the radius of a sphere equal in volume to the rectangular brick and d is the distance between the center of the sphere and the source point at r'. Also, V'JJ Go(r, r)d v = [kod cos(kod) - sin(kod)] ja [kodcos(kod) - sin(kod)] Go(r,)dv kd2 eko a k+ d2Ceakoa {:(x-xC) + y(y-ycc)f+ (z - z) } (A.12) where (xc, Yc, zc) are the coordinates of the sphere's origin. Note that as d -+ 0 all the terms in equations (A.11) and (A.12) converge to a finite quantity.

140 A.3 Results and Validation The implementation of the presented numerical solution is lengthy as is usually the case with most three-dimensional numerical solutions. Herewith, we present some results for the reduced-unknown formulation. 'The results are radar cross section plots corresponding to ferrite plates. We choose plates because of the availability of reference data and codes for these simpler structures. In our implementations, the excitation was the plane wave Ei(r) = [(a. i) i0 + (a. 'i) i] e-jk.r (A.13) where a = i1 cos a + Xi sin a is the polarization vector and the propagation vector is given by ki =- ko (sin Oi cos Xix + sin 0' sin qiy + cos oi) (A.14) whereas Oi and qi are the usual unit vectors in the spherical coordinate system (see Figure A.3). Given the field components in each cell, the far zone scattered field is found from jk0-3kor N Es(r) = {4 E }re/ - -} ek' [r x V' x Ee(r)] dv' (A.15) e=-1 Pre where r is the far zone distance and r is the radial unit vector along the direction of observation. The radar cross section (RCS) 2IE (r)I2 a = lim 47rr E2(r) 2 (A.16) r —oo lEir) 2 was calculated for several ferrite/dielectric plates of different sizes. Herewith we present RCS results for three plate configurations. These are given in Figures A.4-A.6 for Ef and E1 incidence and correspond to square plates of size 0.2Ao x 0.2Ao x 0.025Ao, 0.5Ao x 0.5Ao x 0.025Ao and 1Ao x 1Ao x 0.025Ao, respectively. As seen, the agreement

141 between results based on the reduced-unknown integral equation and those based on a more traditional formulation [71], [72] is excellent. A.4 Conclusion In this work we discretized and implemented a new, reduced-unknown integral equation for scattering by a dielectric layer. Results demonstrating the validity of the integral equation were presented. The singular elements of the impedance matrix were evaluated analytically by resorting to modified Galerkin's testing, employing spheres for the testing elements and rectangular bricks for discretizing the volume. As a result the resulting matrix was asymmetric. When the system is solved using iterative algorithms it may be desirable to employ a testing procedure which leads to a symmetric matrix. This may lengthen the evaluation of the matrix elements but could lead to substantial solution speed-ups. Application of fast algorithms such as the Fast Multipole Method (FMM) and Adaptive Integral Method (AIM) could result in substantial solution speed-ups.

142 Boundary integral formulation Volume integral formulation Figure A.1: Traditional Dielectric/Ferrite body formulation 5 1 a Figure A.2: The rectangular brick element

143 Hi +t E i I / I/77//, i /////,///, / i X H — K-g*~ A A A\ E 0 Figure A.3: Geometry of a dielectric plate illuminated by a plane wave

144 -25.0 -30.0 \ A \ A \ A \ A A Ia r- " C9 6 — C 4 C --- IS. N A \" A N, -35.0 -40.0 -45.0 o CGFFT (E-0 pol) [7] - - - New Intgral eqn.(E-0 pol) A CGFFT (E-~ pol) [7] - - - New Intgral eqn.(E-4 pol) \ A \\ a \ A \ A \ a ' -&A&A I - 0.00 15.00 30.00 45.00 60.00 75.00 90.00 Incidence angle, 0i(degrees) Figure A.4: Backscatter RCS for a 0.2Ao x 0.2Ao x 0.025Ao dielectric/magnetic plate having Er = 7.4 - jl.11 and jt = 1.4 - jO.672 (Sampling cell dimensions: 0.05Ao x 0.05Ao x 0.025Ao)

145 -5.00............ -15.0 - '. - i, -25.0 -.... t:> -35.0 o CGFFT (E-0 pol) [7] -45.0. --- —- New Integral eqn. (E-0 pol) ' * CGFFT (E-0 pol) [7], - New Integral eqn. (E-~ pol) -55.0. I.00 15.00 30.00 45.00 60.00 75.00 90.00 Incidence angle,Oi(degrees) Figure A.5: Backscatter RCS for a 0.5A0 x 0.5A0 x 0.025Ao dielectric/magnetic plate having cr = 7.4 - jl.ll and /tr = 1.4 - jO.672 (Sampling cell dimensions: 0.05Ao x 0.05Ao x 0.025Ao)

146 5.00.. | | | |.. -5.00.500* -25.0, -35.0 - s CGFFI (E-0 pol) [7] -. ---— 4- New Integral eqn. (E-0 pol) * -35.0 '* * CGFFT (E-0 pol) [7] ', ' New Integral eqn. (E-0 pol) -55.0...00 15.00 30.00 45.00 60.00 75.00 90.00 Incidence angle,0i(degrees) Figure A.6: Backscatter RCS for a 1Ao x 1Ao x 0.025Ao dielectric/magnetic plate having Er = 7.4 - jl.ll and ur = 1.4 - jO.672 (Sampling cell dimensions: 0.05Ao x 0.05Ao x 0.025Ao)

147 APPENDIX B Scattering by a Narrow Groove in an Impedance Plane A topic of some concern in radar cross section studies is the scattering from gaps or cracks that may exist where two component parts of a target come together. Even if the crack is entirely or partially filled with some material, it can still provide a significant contribution to the overall scattering pattern of the target, and it is therefore necessary to develop methods for predicting its scattering. The scattering by a groove or an impedance insert in a ground plane has already been considered by a variety of techniques. Integral equation [73]-[74] and finite element-boundary integral solutions [75]-[76] have been effectively used in this respect. Also, in the case of narrow grooves or impedance inserts in a ground plane, closed form solutions have been obtained which were found quite accurate for widths 0.15A or less [32]-[33]. We present a similar solution to the scattering from narrow grooves in an otherwise uniform impedance plane. The solution is useful for computing the scattering by grooves in coated ground planes and is intended to supplement high frequency solutions which are suitable for large width impedance inserts [77],[78]. As in the solutions given in [79], the groove is simulated by an impedance insert, thus, forming a three-part impedance plane (see Figure B.1). The solution of the scattering by the illustrated three-part impedance surface is obtained by introducing appropriate equivalent currents over the extent of the middle impedance strip

148 and then constructing an integral equation for the solution of these currents. This integral equation is obtained by making use of the impedance plane Green's functions which are conveniently expressed in terms of rapidly converging semi-infinite integrals. An important aspect of our solution is the introduction of a one-term basis expansion for the equivalent current, appropriate for small width impedance inserts. The coefficient of the single basis function is found in closed form involving integrals which are functions of the insert's width but independent of the impedance characterizing the insert. Unfortunately, these integrals are also functions of the background impedance, thus, precluding their tabulation for different insert widths. In the following sections we first proceed with the construction of the integral equations for the three-part impedance plane for both TE and TM incidences. The solution of the integral equations are then considered on the assumption of a small width impedance insert. The validity of the results are examined by comparison with a closely related BIM solution [80] and a third order high frequency solution. We conclude by establishing the bounds of the proposed small-width approximations. B. 1 Formulation Consider a narrow filled groove of width w and depth d situated in an otherwise uniform impedance plane as shown in Figure B.1. The groove is illuminated by the plane wave Et(orHI) = ^Jko(xcoso+Ysino) (B.1) for E- (or H) polarization, where ko = 2wx/Ao is the free space wavenumber and qo is the angle of incidence. We shall consider the E- and H-polarizations separately.

149 ( Co, go) Impedance surface (Normalized impedance) \ d:(e':):. -pec x w (a) 1 Th (b) Satisfy the IBC YXYXE =-1] ZoyxH J M 11 11 111 (c) Figure B.1: (a) Geometry of the groove in an impedance plane (b) Impedance model (c) Equivalent scattering problem B.1.1 E-Polarization For TM-incidence, the rectangular groove can be approximated by a strip of impedance [32] j Zo' 7r - tanh(pkod) (B.2) P in which Zo is the free space intrinsic impedance, \(A)- 2 (B.3) d denotes the groove's depth, w is the groove's width, whereas e' and ti are the relative permeability and permittivity of the material filling the groove. This impedance was derived by assuming a single propagating mode in the groove and is consequently accurate when this condition is satisfied irrespective of the groove's depth.

150 The scattering by the subject impedance insert can be represented by introducing the equivalent electric and magnetic currents J= x H M= E x (B.4) over the extent of the groove. These currents are assumed to radiate in the presence of a uniform impedance plane satisfying the boundary condition y x y x E= -11ZoY x H (B.5) where r1 is the associated normalized impedance (see Figure B.1) and (E, H) denote total fields. One can think of (B.5) as being imposed at y = 0-, but at y = 0+ E and H must satisfy the boundary condition y x x E = -Zoy x H -w/2 < x < w/2, y = 0+ (B.6) over the extent of the groove as dictated by the original problem. Using (B.4) in (B.6) we obtain that (assuming E = zEz) M = -rZoJz, Ez = - ZoJz (B.7) i.e, J and M are linearly related over the impedance insert, and thus we need only consider the solution of one of these, say Jz. To solve for Jz, we construct an integral equation by enforcing (B.6). First we decompose Ez(Ez = Z: E) as Ez,= E + E; (B.8) where EZ denotes the scattered field and is equal to that radiated by J and M in the presence of the uniform impedance plane (see Figure B.l(c)). On using (B.7), Es can be expressed as Ez = -jkoZo w/2 Jz(x')GIE(X, Y/X', Y')dx' Jw/2 j Jw/2 -Zo- w/2 Jz(x')GiE(, yx, ylx' ydx (B.9) ay J-w/2

151 where GIE is the electric Green's function of the uniform impedance plane (of normalized impedance r1/). It is given by [81] GIE(, y/X', y') — {H( (o - x')2 + (y - y)2 + QE (B.10) where QE H (o x-x + + y')2 — 2 -2 a'e- vH (ko(x - x')2 + (y + y - )2) dv (B.11) and a' = ko. On substituting (B.9) into (B.8) and then into the second of (B.7) we obtain the integral equation jkoxcoso /2 J-w/2 +rZo- /2 JZ (x')GIE(x, Y -O 0/x', 0)dx (B.12) Oy J-w/2 to be enforced at y 0+ for the solution of Jz(x). B.1.2 H-Polarization For TE-incidence, the rectangular gap can be approximated by a strip of normalized impedance [32] jZop, tan(pkod) (B.13) where p= e (B.14) As before, equivalent currents are introduced and the application of the impedance boundary condition to be satisfied over the groove yields the relations (assuming H= HZ) Mz = i ZoJJ, Ex =: yrZoJx (B.15)

152 where E = -E + EX (B.16) The corresponding tangential scattered E field is given by jo (k+ 2 + w/2 ) 2 J( )GIH(x,y y)dx 0 fw/2 -rlZo- J J(x')GIH(x, y/x, y')dx (B.17) dy J-w/2 where GIH is the Green's function of the uniform impedance surface (of normalized surface impedance /71). It is given by [81] GIH(X, Y/x',y') H (ko (x - x')2 + (y - y)2 - QH (B.18) where QH = H (k -) + ( + 2) -2 j f'e-H VH (k ( -)2 + (y + y'-V)2) dv (B.19) and /3 = k0o71. Substituting (B.17) into the second of (B.15) yields the integral equation 3 k c2 2 \ w/2 sin co ko +z2 J (x')GIH(X, Y - 0+x', O)dx' x J-w/2 o fw/2 d+yy J-w/2 x)GH(X, - x', +rJ (x) (B.20) for the solution of Jx(x). B.2 Integral Equation Solution Typically, the solution of (B.12) and (B.20) can be accomplished numerically using standard techniques. Such a numerical solution is usually the only alternative

153 for moderately sized grooves, but if the groove is small (kow <K 1) then certain simplifications are possible. Once again, we shall consider the E- and H-polarizations separately. B.2.1 E-Polarization It can be shown [82] that HI, and consequently Jz, is of 0(1) near each of the impedance junctions (see Figure B.l(b)) provided neither il nor t are zero. Thus, for kow < 1, we can assume that HI is nearly constant at y = 0+ over the extent of the groove. Based on this, we set Jz(x) = Xe (B.21) where Xe is a constant to be found. Substituting (B.21) into (B.12) and point matching at x = 0 yields Xe I= h l + (Z + 7 I) (B.22) where the integrals I1, 12 and 13 are given by IkoZo W/2 Ho(2)(koll)d' 2 J-w/2 2 -- 7 [n 4)-1} (B.23) -ko =o wf/2 GE(O' O/', O)dx' (B.24) Zo w/2 13 = kfL ' GIE(Y/X) dx' (B.25) 2 0y J-w/2 y-/,O where 7 = 1.781 and GIE(x, y/x', y') = j aceSC H(2) (ko (x - x' + (y + y - jv)2) dv (B.26)

154 In carrying out the integrals in (B.24) and (B.25), GIE is computed numerically. Because of the exponential factor e-~', the infinite integral defining GIE converges very rapidly. For our implementation GCE, was evaluated by subdividing the range of integration in 6 intervals and employing 4 point Gaussian quadrature integration in each interval. When the argument of the Hankel -function is small, to avoid numerical difficulties, the integration from 0 to 0.1 is carried out analytically by introducing the small argument expansion of the Hankel function [83]. B.2.2 H-Polarization From diffraction theory, it can again be shown that Hz is of 0(1) and thus we can approximate J, by W)- = Xh (B.27) where Xh is a constant to be found. Substituting (B.27) into (B.20) and point matching at x == 0 yields sin fo0 1 = (2 + 1) (B.28) The integrals 11 and I2 are given by 1 o2 Gw/2 I1= 2 ( + 2) j G /X,0) dx (B.29) 2k0 X J-w/2 X-+O 12 = L GH(0,y/XO ) d' (B.30) oy J-w/2 ' y-o where GI y Y/XY ) = | /3 e-bL Jo H(2) (ko(x- 2 + (y + y' - j)2) dv (B.31) These can be evaluated numerically without difficulty since GI converges rapidly due to the presence of the exponential factor e-:". Again, it is necessary to introduce

155 the small argument expansion of the Hankel function and carry out the integration from 0 to 0.1 analytically as was done in [83]. B.3 Far Zone Scattered Fields and Echowidth The far zone scattered fields are computed from (B.12) and (B.20) on using the large argument approximation of the Hankel function. By evaluating the integrals via the stationary phase method, we obtain the simple formulae Es = -wZo 0 - (kop — -) 711 sin (B. EZ = -wZoX e -j(kp —4 si (B.32) 2wp 1 +?r1 sin ( E=WZOXh -j(kop —2)^ (B.33) V 27 7/1 + sin b The corresponding echowidth is given by ( xz0 \)2 a ko (7/i sin w)2 l sin E-pol (r-sn H-pol (B.34) B.4 Validation of Results The derived echowidth expressions areare based on a low frequency solution of the exact integral equation. They are, thus, expected to be valid for small groove widths and it is, therefore, of interest to examine their accuracy limitations as the width of the groove increases. For this validation, we used the solution based on the BIM and a high frequency solution. The BIM reference solution was that presented by Moore and Ling [80] and applies to an isolated conductor-backed dielectric gap, such as that shown in Figure B.2. To compare with this solution the dielectric coated conductor was approximated by the normalized impedance Il == i tan (kov/\rl-d) (B.35) E r

156 \ / (!o ) Pec backed Dielectric (represented by equivalent impedance l) / Dielectric pec Gap ( represented by equivalent impedance s ) Figure B.2: Geometry of the gap at the junction of two semi-infinite dielectric coatings where Or and jir are the permittivity and permeability of the dielectric layer and d is the thickness of the dielectric. Also the gap impedances are computed from (B.2) and (B.13). In Figures B.3 and B.4 we compare our small width approximation with the BIM data given in [80] for E and H polarizations, respectively. Each figure displays two curves, one corresponding to a gap of width O.1Ao and the other to a gap of width w = 0.2.Ao. The other gap parameters are given in the figures. As can be expected, the agreement between the two solutions is good for those gaps in the smaller thickness coatings. However, for the 0.2Xo gap in a d = 0.4Ao thick coating, our small width and impedance approximations are no longer expected to be valid. Not surprisingly, the H-polarization echowidth for this geometry as computed via the small width approximation is 2-3 dB off from the BIM solution. The greater disagreement near grazing is due to the inaccurate simulation of the dielectric coating by an impedance surface. Having validated our solution we now proceed with an assessment of its accuracy range and limitation. Figures B.5-B.8 show the backscatter echowidth as a function

157 of the insert's width for different values of 71/ and i. The small width and high frequency solutions are compared in these figures. Of course, by its derivation, the high frequency solution becomes more accurate as w increases, whereas the presented small width approximation does the same as wu decreases. Consequently, the two solutions will have a certain intermediate range of widths where agreement is likely to be expected. Figures B.5 and B.6 show E-polarization curves at some oblique (45~) incidence.. Similarly, Figures B.7 and B.8 display some H-polarization curves at oblique and normal incidences. From these figures, we observe that the presented approximate solutions are indeed in good agreement for 0.15Ao < w < 0.3A0, and in many cases even outside this range. The conclusion that can be reached from Figures B.5-B.8 is that the combination of the simple closed form high frequency and small width approximations provide an accurate evaluation of the scattering by the three part impedance plane for all insert widths. Specifically, the high frequency solution can be used for w > 0.25Ao, whereas the small width approximation is suitable for smaller w. B.5 Summary Integral equations for the analysis of the scattering by a groove in an impedance plane were constructed by representing the groove as an impedance surface. The integral equations for the three-part impedance plane were solved on the assumption of a small width impedance insert. The coefficient of the single basis function for the equivalent current across the insert was explicitly given in terms of integrals which were functions of the insert's width and the background impedance. The accuracy and limitations of the derived small width approximations were evaluated by 1The reference high frequency solution is a modification of that given in [77] using the impedance junction diffraction coefficients [82]

158 comparison with the boundary integral method and a high frequency solution. From these comparisons, we inferred that the presented solution was accurate for insert widths up to 0.25Ao and could be extended up to 0.4A0 for near normal incidences. Thus, the presented approximate small-width approximation and the high frequency solution complement each other for the evaluation of the scattering by a three part impedance plane for all insert widths. 0.000 I I I I I I r"0 /-"0 0-, U -10.0 -20.0 -30.0 -maI - idt I - I -d-2I - ' I M [16](*- -. 1= Small width approx. d 2.. — =.o \Small width approx.(w=.,d=.A\Er=) \ \ \ l E BIM [16](w=O.2X,d=0.2X,=2) \ \.... Small width approx.(w=0.2=kd0.1,r=10) \ A BIM [16](w=0.2k,d=0.1),er=10). AR flf -4U.U -- 1 0.00 15.00 30.00 45.00 60.00 75.00 90.00 Observation Angle (90 - )) [deg] Figure B.3: E-polarization backscatter echowidth for a narrow gap in a grounded dielectric layer as in Figure 2. (a) w=0.1A, d=0.2.A and ~r = 2 (b)w=0.2A, d=0.1A and er = 10

159 0.000 -10.0 9 I I I I I r ~ 1 A AAA A A A AA AAA - A A r-1 "Ic 6-. "C$ -20.0 -30.0 'Xd=0.2',-r=2):2) 'Xhd=0.4XEr=2):2) 'A A I X, A I 'A A I El I -40.0 -50.0L 0.04 Small width approx.(w=O.1 El BIM [9](w=O.1X,d=O.2X,Ir= Small width approx.(w=O.2 A BIM [9](w=0.2X,d=0.4X,E r 0 r r I. i 15.00 30.00 45.00 60.00 75.00 90.00 Observation Angle (4) [deg] Figure B.4: H-polarization backscatter echowidth for a narrow gap in a grounded dilectric layer as in Figure 2. (a) w=O.1A, d=O.2A and r, = 2 and (b) w=O.2A, d=O.4A and (fr = 2

160 0 rn=0.25 qr=-j0.25.. --- ——. =0.25 rl=-j0.25 A n=0.25 n_=-i.0 r-" 73 6 — 4-) tQ CIl 5.00 0.000 -5.00 -10.0 -15.0 -20.0 -25.0 -30.0 *.'1 J ""., r- =0.25 n1i=-jo.5 0<> =0.25 r1=-j0.75 ---- =0.25 1=-jo.75 A A Aoo o o -... A A - - -. -.. -^ A- -0 <> A o' - A - A, I, I I I. I -35.0 -40.0 ).( 0DO I., 1 0.05 0.10 0.15 0.20 0.25 0.30 Width[X] 0.35 Figure B.5: E-polarization backscatter echowidth at oblique incidence (45~) for an impedance insert of normalized impedance r = 0.25 in various impedance planes as a function of the insert's width (symbols: small width approx, lines: high freq solution)

161 o Ti=0.15 71=-j0.25.-.. n=0.15 n1=-j0.25 5.00 A rl=j. 15 rl=-jO.5.. --- - r=jO.15rlq=-j0.5 0.000 O rq=jO.25 r1=-jO.75 = -10.0 A a 0-15.00. -20.0 > A...-.. -35.0 - o -40.0. -. 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 Width[X] Figure B.6: E-polarization backscatter echowidth at oblique incidence (450) for various three-part impedance planes as a function of the insert's width (symbols: small width approx, lines: high freq solution)

162 r-" "10 L-.-j 10 4-A tb cu " —4 (1) C13 4 —j 7. P-4 0.J= C.) 4 5.00 0.000 -5.00 -10.0 -15.0 -20.0 -25.0 -30.0.0 rj=O.25rj1=O A ri=0. 15 i1 =0.. ---- - n=0.15rj1=0 <c- ri=O.5 i1,=O. --- — i=O.5ri1,=O -i KE~ ---------------- - --------------------------- 0 (DAA — A A A - A A -35.0 [-40. 0 -- 0.00 (D A 4, I I - I I I I 0.05 0.10 0.15 0.20 0.25 0.30 0.35 WidthR]~ Figure B.7: H-polarization backscatter echowidth at oblique incidence (450) for various impedance inserts in a metallic plane as a function of the insert's width (symbols: small width approx, lines: high freq solution)

163 5.00 ' i i0.000 - 0'O 0 A -' - 2 0.0' ^ | ' A-O.0-, A -1 - ) -10.0 - -15.0 --20.0 - A.0 r=.25 ill=j0.75 - - - r =.25 1,1=j0.75 -25.0 - A rj=j0.25 r,=j0.75 A A -30.0 ------ i=j0.25 ll=jO.75 -35.0 0 9l=-jO.25 i,=jO.75 ri=-jO.25 ii=jO.75 -40.0 - - 0.00 0.10 0.20 0.30 0.40 0.50 Width[X] Figure B.8: H-polarization backscatter echowidth at normal incidence for various impedance inserts in an impedance plane of normalized impedance 7/1 = j0.75 as a function of the insert's width (symbols: small width approx, lines: high freq solution)

BIBLIOGRAPHY 164

165 BIBLIOGRAPHY [1] R.F. Harrington, Field Computation by moment methods. IEEE press: New York, 1993. [2] J.M. Jin and J.L. Volakis, "TE scattering by an inhomogeneously filled aperture in a thick conducting plane," IEEE Transactions on Antennas and Propagation, vol. 38, pp. 1280-1286, August 1990. [3] B.Z. Steinberg and Y. Leviatan, "On the use of wavelet expansions in the method of moments," IEEE Transactions on Antennas and Propagation, vol. 41, pp. 610-619, May 1993. [4] K. Sabetfakhri and L. Katehi, "Analysis of integrated millimeter-wave and submillimeter-wave waveguides using orthonormal wavelet expansions," IEEE Transactions on Microwave Theory and Techniques, vol. 42, no. 12, pp. 659-667, 1994. [5] I. Daubechies, "Orthonormal bases of compactly supported wavelets," Commun. Pure Appl. Math., vol. 41, pp. 909-996, 1988. [6] F. Canning. "Improved impedane matrix localization method," IEEE Transactions on Antennas and Propagation, vol. 41, pp. 659-667, 1993. [7] A. Boag and R. Mittra, "Complex multipole beam approach to electromagnetic scattering problems," IEEE Transactions on Antennas and Propagation, vol. 42, no. 3, pp. 366-372, 1994. [8] R. Kastner and G. Nocham, "Matrix thinning with reduced field testing (RFT)," Microwave and Optical Technology Letters, vol. 9, no. 6, pp. 329-333, 1995. [9] J.L. Volakis and S.V. Krestianinov, "ISOMOM: A new method for reducing storage and execution time of moment-method solutions in electromagnetics," Microwave and Optical Technology Letters, vol. 8, pp. 184-186, March 1995. [10] S. S. Bindiganavale and J. L. Volakis, "ISOMOM: A sparse moment method technique for a wide class of 2D and 3D scattering problems," in 1995 IEEE Int. Symp. on Antennas and Propagation Digest, pp. 1536-1539, June 1995. [11] G.K. Gothard and S.M. Rao, "A new technique to generate sparse matrix using the method of moments - Application to two-dimensional problems," in 1995 IEEE Int. Symp. on Antennas and Propagation Digest, June 1995. [12] N.N.Bojarski, "k-space formulation of the electromagnetic scattering problem," Tech. Rep. AFAL-TR-71-5, USAF, March 1971.

166 [13] 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 Transactions on Antennas and Propagation, vol. 34, no. 5, 1986. [14] C.Y.Shen, K.J.Glover, M.I.Sancer, and A.D.Varvatsis, "The discrete Fourier transform method of solving differential integral equations in scattering theory," IEEE Transactions on Antennas and Propagation, vol. 37, pp. 1032-1049, August 1989. [15] K. Barkeshli and J.L. Volakis, "On the implementation of the conjugate gradient fourier transform method for scattering by planar plates," IEEE Transactions on Antennas and Propagation, vol. 32, no. 2, pp. 20-26, 1990. [16] J. Barnes and P. Hut, "A hierarchical O(N log N) force-calculation algorithm," Nature, vol. 324, no. 4, pp. 446-449, 1986. [17] M.S. Warren and J.K. Salmon, "Astrophysical N-body simulations using hierarchical tree data structures," in Proceedings Supercomputing 92, pp. 570-6, 1992. [18] H-.Q. Ding, N. Karasawa, and W. A. Goddard III, "Atomic level simulations on a million particles: The cell multipole method for Coulomb and London nonbond interactions," J. Chem. Phys., vol. 97, no. 6, pp. 4309-4315, 1992. [19] L. Greengard, The rapid evaluation of potential fields in particle systems. The MIT Press, 1988.. [20] V. Rokhlin, "Rapid solution of integral equations of scattering theory in two dimensions," Journal of Computational Physics, vol. 86, no. 2, pp. 414-439, 1990. [21] R. Coifman, V. Rokhlin, and S. Wandzura, "The fast multipole method for the wave equation: A pedestrian prescription," IEEE Antennas and Propagation Magazine, vol. 35, no. 3, pp. 7-12, 1993. [22] J. Song and W. Chew, "Multilevel fast multipole algorithm for soving combined field integral equation of electromagnetic scattering," Microwave and Optical Technology Letters, vol. 10, pp. 14-19, 1995. [23] R.J. Burkholder and D.H. Kwon, "High-frequency asymptotic acceleration of the fast multipole method," Radio Science, vol. 31, no. 5, pp. 1199-1206, 1996. [24] S. S. Bindiganavale and J. L. Volakis, "Guidelines for using the fast multipole method to calculate the RCS of large objects," Microwave and Optical Technology Letters, vol. 11, no. 4, pp. 190-194, 1996. [25] N. Lu and J. M. Jin, "Application of the fast multipole method to finite elementboundary integral solution of scattering problems," IEEE Transactions on Antennas and Propagation, vol. 44, no. 6, pp. 781-786, 1996. [26] S. S. Bindiganavale and J. L. Volakis, "A hybrid FEM-FMM technique for electromagnetic scattering," IEEE Transactions on Antennas and Propagation, vol. 45, no. 1, pp. 180-181, 1997.

167 [27] S. S. Bindiganavale and J. L. Volakis, "Comparison of fast integral mesh truncation schemes for hybrid FE-BI systems," lEE Electronics Letters, vol. 33, no. 11, pp. 924 -925, 1997. [28] S.S. Bindiganavale and J.L. Volakis, "Comparison of three FMM techniques for solving hybrid FE-BI systems," IEEE Antennas and Propagation Magazine, Accepted for publication. [29] E. Bleszynski, M. Bleszynski, and T. Jaroszewicz, "AIM: Adaptive integral method for solving large-scale electromagnetic scattering and radiation problems," Radio Science, vol. 31, no. 5, pp. 1225-1251, 1996. [30] E. Arvas and S. Ponnapalli, "Scattering cross section of a small radome of arbitrary shape," IEEE Transactions on Antennas and Propagation, vol. 37, pp. 655-658, May 1989. [31] J.L. Volakis, "Alternative field representations and integral equations for modeling inhomogeneous dielectrics," IEEE Transactions on Microwave Theory and Techniques, vol. 40, pp. 604-608, March 1992. [32] T.B.A. Senior and J.L. Volakis, "Scattering by gaps and cracks," IEEE Transactions on Antennas and Propagation, vol. 37, pp. 744 —750, June 1989. [33] K. Barkeshli and J.L. Volakis, "Electromagnetic scattering from an aperture formed by a rectangular cavity recessed in a ground plane," Journal of Electromagnetic Waves and Applications, 1991. [34] J.L. Volakis, EECS 633 Coursepack. University of Michigan, 1994. [35] T.B.A. Senior and J.L. Volakis, Approximate boundary conditions in electromagnetics. IEE press: London, 1995. [36] S. Rao, D. Wilton, and A. Glisson, "Electromagnetic scattering by surfaces of arbitrary shape," IEEE Transactions on Antennas and Propagation, vol. 30, no. 3, pp. 409-418, 1982. [37] D.R. Wilton, S.M. Rao, A.W. Glisson, D.H. Schaubert, O.M. Al-Bundak, and C.M. Butler, "Potential integrals for uniform and linear source distributions on polygonal and polyhedral domains," IEEE Transactions on Antennas and Propagation, vol. 32, pp. 276-281, March 1984. [38] J. Gong, J.L. Volakis, and A.C. Woo, "A hybrid finite element-boundary integral method for the analysis of cavity-backed antennas of arbitrary shape," IEEE Transactions on Antennas and Propagation, vol. 42, no. 9, pp. 1233-1242, 1994. [39] J. Gong, J.L. Volakis, and H.T.G. Wang, "Efficient finite element simulation of slot antennas using prismatic elements," Radio Science, vol. 31, no. 6, pp. 1837-1844, 1996. [40] S. S. Bindiganavale and J. L. Volakis, "Computation of nose radome scattering by employing the fast multipole method," in 1996 IEEE Int. Symp. on Antennas and Propagation Digest, July 1996.

168 [41] Major. Paul Skinner, "Personal communication," Air Force Institute of Technology, 1994. [42] C. Lu and W. Chew, "Fast far field approximation for calculating the RCS of large objects," Microwave and Optical Technology Letters, vol. 8, no. 5, pp. 238-241, 1995. [43] J.J.H. Wang, Generalized moment methods in electromagnetics. Wiley-Interscience, 1991. [44] M.J. Schuh and A.C. Woo, "Code scaling," IEEE Transactions on Antennas and Propagation, 1995. [45] Y.T. Lo and S.W. Lee ed., Antenna handbook. Van Nostrand: New York, 1988. [46] M. Naor and T.B.A. Senior, "Scattering by resistive plates," Tech. Rep. 018803-1-T, Radiation Laboratory, University of Michigan, September 1981. [47] J.D. Walton ed., Radome engineering handbook. Marcel-Dekker: New York, 1970. [48] J.M Putnam and L.N. Medgyesi-Mitschang, "Combined field integral equation formulation for axially inhomogeneous bodies of revolution," Tech. Rep. QA003, McDonnell Douglas research laboratories, St.Louis, MO 63166, December 1987. [49] J.L. Volakis, J. Gong, and A. Alexanian, "A finite element boundary integral method for antenna RCS analysis," Electromagnetics, vol. 14, no. 1, pp. 83-85, 1994. [50] J.-M. Jin and J.L. Volakis, "TM scattering by a inhomogeneously filled aperture in a thick conducting plane," IEE Proceedings, Part H, vol. 137, pp. 153-159, June 1990. [51] S. S. Bindiganavale and J. L. Volakis, "A hybrid FEM-FMM technique for electromagnetic scattering," in 12th Annual Review of Progress in Applied Computational Electromagnetics (ACES) Digest, Naval Postgraduate school, Monterey CA, pp. 563 — 570, March 1996. [52] N. Lu and J. M. Jin, "Application of the fast multipole method to finite elementboundary integral solution of scattering problems," in 12th Annual Review of Progress in Applied Computational Electromagnetics (ACES) Digest, Naval Postgraduate school, Monterey CA, pp. 1182-1189, March 1996. [53] J.-M. Jin, J.L. Volakis, and V.V. Liepa, "Fictitious absorber for truncating finite element meshes in scattering," IEE Proceedings, Part H, vol. 139, no. 5, pp. 472-476, 1992. [54] A.C. Polycarpou, M.R. Lyons, J. Aberle, and C.A. Balanis, "Analysis of arbitrary shaped cavity-backed patch antennas using a hybritization of the finite element and spectral domain methods," in 1996 IEEE Int. Symp. on Antennas and Propagation Digest, July 1996. [55] 0.0. Storaasli, "Performance of NASA equation solvers on computational mechanics application,"' in AIAA Paper No.96-1505, April 1996. [56] R.C. Young, "Fast matrix solver," Multipath Corporation, 1996.

169 [57] S.S. Bindiganavale, J.L. Volakis, and H. Anastassiu, "Scattering from plates containing small features using the adaptive integral method (AIM)," IEEE Transactions on Antennas and Propagation, submitted for publication. [58] H.T. Anastassiu, M. Smelyanskiy, S. Bindiganavale, and J.L. Volakis, "Scattering from relatively flat surfaces using the adaptive integral method (AIM)," Radio Science, accepted for publication. [59] S.S. Bindiganavale and J.L. Volakis, "Solution of a new reduced-unknown integral equation for modeling inhomogeneous dielectrics," Journal of Electromagnetic Waves and Applications, vol. 8, no. 11, pp. 1531-1541, 1994. [60] S.S. Bindiganavale and J.L. Volakis, "Scattering by a narrow groove in an impedance plane," Radio Science, vol. 31, no. 2, pp. 401-408, 1996. [61] O.S. Haddadin and E.S. Ebbini, "Ultrasonic focusing through inhomogeneous media by application of the inverse scattering problem," Journal of Acoustic Society of America, accepted for publication. [62] O.S. Haddadin and E.S. Ebbini, "Multiple frequency distorted born iterative method for tomogrpahic imaging," IEEE transactions on ultrasonics, ferroelectrics, and frequency control, submitted for publication. [63] L.C. Kempel and J.L. Volakis, "Scattering by cavity-backed antennas on a circular cylinder," IEEE Transactions on Antennas and Propagation, vol. 42, no. 9, pp. 1268 -1279, 1994. [64] L.C. Kempel, J.L. Volakis, and R.J. Sliva, "Radiation by cavity-backed antennas on a circular cylinder," IEE Proceedings, Part H, vol. 142, pp. 233-239, 1995. [65] Y. Saad, Iterative methods for sparse linear systems. PWS Publishing Company: Boston, 1996. [66] R.L. Wagner and W.C. Chew, "A study of wavelets for the solution of electromagnetic integral equations," IEEE Transactions on Antennas and Propagation, vol. 43, no. 8, pp. 802-810, 1995. [67] D.H.Schaubert, D.R.Wilton, and A.W.Glisson, "A tetrahedral modeling method for electromagnetic scattering by arbitarily shaped inhomogeneous dielectric bodies," IEEE Transactions on Antennas and Propagation, vol. 32, pp. 77-85, January 1984. [68] A.Bossavit and I.Mayergoyz, "Edge elements for scattering problems," IEEE Transactions on Magnetics, vol. 25, pp. 2816-2821, July 1989. [69] J.M.Jin and J.L.Volakis, "Electromagnetic scattering by and transmission through a three-dimensional slot in a thick conducting plane," IEEE Transactions on Antennas and Propagation, vol. 39, pp. 543-550, April 1991. [70] C.T.Tsai, H.Massoudi, C.H.Durney, and M.F.Iskander, "A procedure for calculating fields inside arbitarily shaped inhomogeneous dielectric bodies using linear basis functions with the moment method," IEEE Transactions on Microwave Theory and Techniques, vol. 34, pp. 1131-1139, November 1986.

170 [71] T.J.Peters and J.L.Volakis, "Application of a conjugate gradient fft method to scattering from thin planar material plates," IEEE Transactions on Antennas and Propagation, vo]L. 36, pp. 518-526, April 1988. [72] C.Y.Shen, "Application of the discrete fourier transform method to thin dielectric structures,'" IEEE Transactions on Antennas and Propagation, vol. 37, pp. 1277-1283, October 1989. [73] R.F.Harrington and J.R.Mautz, "Generalized network formulation for aperture problems," IEEE Transactions on Antennas and Propagation, vol. 24, pp. 870-873, November 1976. [74] K. Barkeshli and J.L. Volakis, "TE scattering by a two-dimensional groove in a ground plane using higher order impedance boundary conditions," IEEE Transactions on Antennas and Propagation, vol. 38, pp. 1280-1286, August 1990. [75] S.K.Jeng, "Scattering from a cavity-backed slit in a ground plane - TE case," IEEE Transactions on Antennas and Propagation, vol. 38, pp. 1523-1529, October 1990. [76] K.W.Whites, E.Michielssen, and R.Mittra, "Approximating the scattering by a material-filled 2-D trough in an infinite plane using the impedance boundary condition," IEEE Transactions on Antennas and Propagation, vol. 41, pp. 146-153, February 1993. [77] M.I.Herman and J.L.Volakis, "High frequency scattering by a double impedance wedge," IEEE Transactions on Antennas and Propagation, vol. 36, pp. 664-678,1988. [78] R.G.Rojas, H.C.Ly, P.H.Pathak, and R.Tiberio, "Electromagnetic plane wave diffraction by a three-part thin, planar dielectric/magnetic slab," Radio Science, vol. 26, no. 5, pp. 1267-1280, 1988. [79] K. Barkeshli and J.L. Volakis, "Scattering from narrow rectangular filled grooves," IEEE Transactions on Antennas and Propagation, vol. 39, pp. 804-810, 1991. [80] J.Moore and H.Ling, "Scattering by gaps in coated structures," Journal of Electromagnetic Waves and Applications, vol. 7, no. 3, p. 1993, 1993. [81] K.Sarabandi, "Scattering from dielectric structures above impedance surfaces and resistive sheets," IEEE Transactions on Antennas and Propagation, vol. 40, no. 1, pp. 68-78, 1992. [82] T.B.A. Senior, "Diffraction by half plane junctions," Tech. Rep. RL 892, Radiation Laboratory, University of Michigan, March 1993. [83] M.A. Ricoy and J.L. Volakis, "Integral equations with reduced unknowns for the simulation of two-dimensional composite structures," IEEE Transactions on Antennas and Propagation, vol. 37, pp. 362-372, March 1989.