EP1295223A2 - Method and apparatus for modeling three-dimensional electromagnetic scattering from arbitrarily shaped three-dimensional objects - Google Patents

Abstract

A computer program product, an apparatus, and a method for modeling equivalent surface sources on a closed, arbitrarily shaped object that result from the imposition of an arbitrary time-harmonic incident field (e.g., a radar wave) are provided. An object divided into by a plurality of patches 302 and parameters for the incident field and the properties of the object 304 are provided to the system. Next, the system produces a discretized representation of a well-conditioned boundary integral equation in the form of a well-conditioned matrix equation, modeling the interaction between the incident field and the object 306. The well-conditioned linear system is solved numerically to determine equivalent surface sources on the object 308. The equivalent sources may then be used to determine electromagnetic fields resulting therefrom 310. The invention provides improved computational efficiency as well as increased modeling accuracy by effectively reducing the numerical precision necessary.

Description

METHOD AND APPARATUS FOR MODELING THREE-DIMENSIONAL ELECTROMAGNETIC SCATTERING FROM ARBITRARILY SHAPED THREE- DIMENSIONAL OBJECTS

STATEMENT OF GOVERNMENT RIGHTS

This invention was made with U.S. Government support from the Defense Advanced Research Projects Agency, under Contract No. MDA972-95-C-0021. The U.S. Government has certain rights in this invention.

PRIORITY CLAIM

This application claims the benefit of priority to provisional application number 60/213,102, filed in the United States on June 21, 2000, and titled "Well-conditioned Boundary Integral Equations and Discretization for Three-Dimensional Electromagnetic Scattering".

FIELD OF THE INVENTION

The present invention relates generally to computer aided design and object modeling. More specifically, it relates to a computer program product, apparatus, and method for determining equivalent surface sources (or actual surface currents in the case of a metal object) on the surface of a closed, arbitrarily shaped object (e.g., a component of a vehicle such as an airplane part) that result from the imposition of an arbitrary time-harmonic incident field (e.g., a radar wave).

BACKGROUND Software utilizing boundary integral equation techniques to solve electromagnetic scattering problems (e.g., using the method of moments) has existed at least since the 1960s. Such techniques have the potential for making very accurate radar cross-section predictions.

However, practical application of the techniques to realistic problems is severely limited by poor scalability of computational resource requirements (i.e., CPU time and memory) relative to the problem size and complexity. With the relatively recent introduction of computational techniques involving high-order and "fast" mathematical methods (such as the fast multipole method (FMM)), the situation has changed drastically. Now, instead of the resource requirements for a problem described by N unknowns scaling in O(N³) time, they scale in O(NlogN) time, much like the Fast Fourier Transform (FFT) algorithm. Consequently, it has become possible to handle a much wider array of problems of practical interest. Although the methods used have improved drastically, fast methods must still be used in conjunction with iterative solvers. If the set of linear equations that arises from discretizing a boundary integral equation is ill conditioned, the number of iterations required to obtain an accurate solution is not well controlled, and may be very large. The overall efficiency of the solution process suffers accordingly. In addition, a high degree of numerical precision is generally required in the solution process, causing the need for rounding, and resulting in a loss of solution accuracy.

Unfortunately, the most widely used integral equation formulation for closed, perfect electrical conductors (PEC), namely, the combined field integral (CFIE) equation, is not, in general, well conditioned. The ill conditioning stems from the fact that the CFIE necessarily includes an integral operator component that is "hypersingular". This component is, in fact, not simply an integral operator at all, but rather the gradient of an integral operator. Whereas, ordinary integral operators tend to filter out the numerical "noise" associated with a discretization, differential operators tend to amplify it.

The eigenvalue spectrum of a hypersingular integral operator is unbounded, which leads directly to an ill-conditioned linear system. The ill effects become evident in the form of an increase in the required number of iterations as one refines the discretization of a target. Fine discretizations are often required to spatially resolve source variations or geometric detail on the subwavelength scale. This difficulty with the CFIE (as well as with some other boundary integral formulations in computational electromagnetics) is known as the "low frequency problem."

Various methods have been employed to overcome this problem. The most commonly practiced method is to multiply the CFIE matrix by a "preconditioner" matrix at every step of the iterative solver procedure. An effective preconditioner improves the condition number of a linear system, thereby reducing the number of iterations required for solution. Examples of preconditioners that have been utilized are the block-diagonal preconditioner, the incomplete LU preconditioner, and the sparse approximate inverse preconditioner. Preconditioners such as these are general purpose - in other words, they can be applied with varying degrees of success to any matrix. Unfortunately, however, none of them takes specific advantage of the analytical properties of the CFIE, so none of them is ideally suited for use with the CFIE. Instead of preconditioning the equation, it is also possible to "post-condition" it. This involves defining the unknown part of the original linear equation to be the driving term for some new linear system. If the combined linear system is properly conditioned, then the problem is effectively solved because, once the new auxiliary unknowns have been solved for, it is straightforward to determine the solution to the original linear system. The moment method approaches that use "loop-star" basis functions fall into this category. While some success has been achieved with this technique, it is cumbersome with limited extensibility, and has been applied only to low-order, moment method discretizations.

Therefore, it is desirable to provide a preconditioner method that is specifically tailored for use with the CFIE. More specifically, the preconditioner should be designed to minimize the number of iterations required for a solution. This is especially desirable when the use of a fine discretization to resolve source variations or geometric detail on a subwavelength scale results in an ill-conditioned linear equation (the "low-frequency" problem). Most specifically, it is desirable to provide a computer program product, an apparatus, and a method for modeling the electromagnetic response of a closed arbitrarily shaped three-dimensional object to an arbitrary time-harmonic incident field by means of a well-conditioned boundary integral equation (BIE) by discretizing the well-conditioned BIE to provide and solve a well-conditioned finite dimensional linear system to determine the distribution of equivalent surface sources on the object.

References:

[1] Leslie Greengard and Stephen Wandzura. Fast multipole methods. IEEE Computational Science & Engineering, 5(3):16-18, July 1998.

[2] John J. Ottusch, Mark A. Stalzer, John L. Visher, and Stephen M. Wandzura. Scalable electromagnetic scattering calculations for the SGI Origin 2000. In Proceedings SC99, Portland, OR, November 1999. IEEE.

[3] Jiming M. Song and Weng Cho Chew. The fast Illinois solver code: Requirements and scaling properties. IEEE Computational Science & Engineering, 5(3):19-23, July 1998. [4] Andrew F. Peterson. The "interior resonance" problem associated with surface integral equations of electromagnetics: Numerical consequences and a survey of remedies. Electromagnetics. 10:293-312, 1990.

[5] George C. Hsiao and Ralph E. Kleinman. Mathematical foundations for error estimation in numerical solutions of integral equations in electromagnetics. IEEE Transactions on Antennas and Propagation, 45(3):316-328, March 1997.

[6] Nagayoshi Morita, Nobuaki Kumagai, and Joseph R. Mautz. Integral Equation Methods for Electromagnetics. Artech House, Boston, 1990.

[7] Robert J. Adams. A Class of Robust and Efficient Iterative Methods for Wave Scattering Problems. Ph.D. Thesis, Virginia Polytechnic Institute and State University, Blacksburg, VA, December 1998.

[8] Robert J. Adams and Gary S. Brown. Stabilization procedure for electric field integral equations. Electronic Letters, 35(23):2015-2016, November 1999.

[9] Peter Kolm and Vladimir Rokhlin. Quadruple and octuple layer potentials in two dimensions I: Analytical Apparatus. Technical Report YALEU/DCS/RR-1176, Yale University, Department of Computer Science, March 1999.

[10] Barrett O'Neil. Elementary Differential Geometry. Academic Press, New York, 1997.

[11] John David Jackson. Classical Electrodynamics. John Wiley & Sons, New York, second edition, 1975.

[12] Milton Abramowitz and Irene A. Stegun. Handbook of Mathematical Functions. Applied Mathematics Series. National Bureau of Standards, Cambridge, 1972.

[13] Lawrence S. Canino, John J Ottusch, Mark A. Stalzer, John L. Visher, and Stephen M. Wandzura. Numerical solution of the Helmholtz equation in 2d and 3d using a high-order Nystrδm discretization. Journal of Computational Physics, 146:627-663, 1998. SUMMARY OF THE INVENTION The present invention provides an improved method for modeling electromagnetic scattering from an arbitrarily shaped three-dimensional object to an arbitrary time-harmonic incident field. It allows for the determination of the scattering with an improved degree of numerical accuracy, while also providing for increased computational efficiency.

In a first embodiment, the present invention provides a computer program product, used in a computer system having a processor, for modeling an electromagnetic response of an arbitrarily shaped three-dimensional object to an arbitrary time-harmonic incident field by means of a well- conditioned boundary integral equation (BEE). The computer program product comprises recording medium; mans, recorded on the recording medium for producing a discretized representation of a well-conditioned BIE to provide a well-conditioned, finite-dimensional linear system whose solution approximates a distribution of equivalent surface sources on a model caused by the arbitrary time-harmonic incident field; and means, recorded on the recording medium, for solving the linear system to determine values of the equivalent surface sources; whereby solution of the well-conditioned finite-dimensional linear system provides improved computational efficiency and solution accuracy. Additionally, this embodiment of the present invention may comprise means, recorded on the recording medium for modeling electromagnetic scattering from the three-dimensional object due to the time-harmonic incident field as determined from the equivalent surface sources. Preferably, the computer program product is tailored for determining the electromagnetic scattering from a time-harmonic incident field in the form of a radar wavefront and where the three-dimensional object represents a component of a vehicle.

In another, more specific embodiment, the BIE produced by the computer program product is in the form of:

(E excitation term)+ γijΑ excitation term) = {(E integral operator term)+ ^(H integral operation term)}

where, γ is a complex constant and J represents the unknown surface source distribution. In yet another, still more specific embodiment, the BIE produced by the computer program product is in the form of: kT^s (ik)n x E'"^c + iZT ^a (ik)n • H "^c = {i [- T^a (ik) ° T^τ (k) + T^p (ik) ° T^L (k) - k²T^s (ik) - T^s (k)\- K* (k)}j

In a further embodiment, the present invention provides an apparatus for modeling the electromagnetic response of an arbitrarily shaped three-dimensional object to an arbitrary time- harmonic incident field by means of a well-conditioned boundary integral equation (BIE). The apparatus comprises a computer system having a processor, a memory connected with the processor, and an input/output device connected with the processor for receiving parameters representing an arbitrarily shaped three-dimensional object and an arbitrary time-harmonic incident field; means, operating in the processor and the memory of the computer for utilizing the parameters representing the arbitrarily shaped three-dimensional object and the arbitrary time-harmonic incident field for producing a discretized representation of a well-conditioned BIE to provide a well-conditioned finite dimensional linear system whose solution approximates a distribution of equivalent surface sources on a model caused by the time-harmonic incident field; and means, operating in the processor and the memory of the computer for solving the linear system to determine values of the equivalent surface sources; whereby solution of the well-conditioned finite dimensional linear system provides improved computational efficiency and solution accuracy. Similar to the computer program product embodiments, the apparatus further comprises means, operating in the processor and the memory of the computer for modeling electromagnetic scattering from the three-dimensional object due to the time-hannonic incident field as determined from the equivalent surface sources. Again, the apparatus is preferably tailored for determining the electromagnetic scattering from a time-harmonic incident field represented in the form of a radar wavefront and where the three-dimensional object represents a component of a vehicle.

In another, more specific embodiment, the BIE produced by the apparatus is in the form of:

(E excitation term)+ ^(H excitation term) = {(E integral operator term)+ ^(H integral operation term)}

where, γ is a complex constant and J represents the unknown surface source distribution.

In yet another, still more specific embodiment, the BEE produced by the apparatus product is in the form of: kT^s (ik)n x E'"^c + iZT^a (ifc)n • H"° = { [- T (ik) o T^τ (k) + r^p (ik) ° r^L (Jfc) - k²T^s (ik) o T^s (Jfc)]- JT (£)}j

The present invention further provides a method for modeling an electromagnetic response of an arbitrarily shaped three-dimensional object to an arbitrary time-harmonic incident field by means of a well-conditioned boundary integral equation (BIE). The method operates on a computer system having a processor, a memory connected with the processor, and an input/output device connected with the processor for receiving parameters representing an arbitrarily shaped three- dimensional object and an arbitrary time-harmonic incident field. The method comprises the steps of: utilizing the parameters representing the arbitrarily shaped three-dimensional object and the arbitrary time-harmonic incident field for producing a discretized representation of a well-conditioned BIE to provide a well-conditioned finite dimensional linear system whose solution approximates a distribution of equivalent surface sources on a model caused by the arbitrary time-harmonic incident field; and solving the linear system to determine values of the equivalent surface sources; whereby solution of the well-conditioned finite-dimensional linear system provides improved computational efficiency and solution accuracy.

Additionally, the method can further comprise the step of modeling electromagnetic scattering from the three-dimensional object due to the time-harmonic incident field as determined from the equivalent surface sources. The method is preferably operated where the time-harmonic incident field is a radar wavefront and wherein the three-dimensional object represents a component of a vehicle, in order to determine the scattering of radar waves from a vehicle component due to in incident radar wave.

In another, more specific embodiment, the BIE produced by the method is in the form of: (E excitation term)+ χ(ϊL excitation term) = {(E integral operator term)+ γ(jΑ integral operation term)}

where, γ is a complex constant and J represents the unknown surface source distribution.

In yet another, still more specific embodiment, the BIE produced by the method is in the form of: kT^s (ik)n x E'"^c + iZ (ik)n ■ W"^c = {i [- T^a (ik) ° T^τ (k) + -T^p (ik) o A (k) - k²T^s (ik)

DESCRIPTION OF THE FIGURES FIG. 1 is a block diagram representing a general purpose computer system used in conjunction with the present invention;

FIG. 2(a) is an illustrative diagram showing an object having a diameter of 2.5λ (where λ represents a wavelength), and corresponding to element 300 of FIG. 3; FIG. 2(b) is an illustrative diagram showing the object of FIG. 2(a) modeled as a plurality of triangular patches, and corresponding to element 302 of FIG. 3;

FIG. 2(c) is an illustrative diagram showing the object of FIG. 2(b) with an incident field, and corresponding to element 304 of FIG. 3;

FIG. 2(d) is an illustrative diagram showing the object of FIG. 2(c) with an incident field and with electromagnetic fields resulting from equivalent surface sources due to the incident field, and corresponding to elements 306, 308, and 310 of FIG. 3;

FIG. 3 is a flow diagram demonstrating the method of the present invention in conjunction with the steps of inputting object and incident field parameters, dividing the model of the object into a plurality of patches, and the optional step of determining electromagnetic fields; and FIG. 4 is an illustrative diagram of a computer program product according to the present invention, stored on CD or DVD.

DETAILED DESCRIPTION The present invention provides a computer program product, a method, and an apparatus for generating and solving well-conditioned boundary integral equations for modeling three dimensional electromagnetic scattering of an incident wavefront from an object. It should be understood that this description is not intended to limit the invention. On the contrary, the invention is intended to cover alternatives, modifications and equivalents, which are included in the spirit and scope of the invention as defined by the appended claims. Furthermore, in the following detailed description of the present invention numerous specific details are set forth in order to provide a thorough understanding of the present invention. However, it will be obvious to one of ordinary skill in the art that the present invention may be practiced without the specific details.

Introduction

In the presence of an incident wavefront, the surface properties of an object affect scattering of the wave from the object's surface. By properly modeling the object and the incoming wavefront, it is possible to determine how waves will be scattered from the object's surface. Determination of this scattering is of particular importance with regard to designing objects having reduced radar cross-sections as well as in designing antennae or other objects that are configured to maximize electromagnetic interception.

The present invention, in a general sense, provides a computer program product, a method, and an apparatus for receiving a computer model of an object comprised of a mesh of a plurality of patches; for specifying an incident wavefront having an excitation wavelength, wavefront geometry (e.g., plane wave, etc.), and incident angle; and, based on the properties of the object and the wavefront, determining the equivalent surface sources resulting from the interaction of the wavefront with the object, and from the equivalent surface sources, modeling the scattering of energy from the surface of the object.

Invention Overview

The following overview is presented with reference to the figures in order to provide a high-level conceptual understanding of the present invention prior to discussing the mathematical details and specific advantages inherent therein. First, an overview of the data processing system of the present invention is provided, followed by a brief illustration and flow diagram explaining the operation of the present invention. Subsequent to the overview, further details of the present invention will be provided, followed by a set of numerical results demonstrating the advantages of the modified CF E of the present invention over the conventional CFIE. Finally, an appendix and glossary are provided to give the reader a convenient source for understanding some of the terms and conventions utilized herein.

Data Processing System The present invention is designed to be operated on a data processing system, which may comprise a general or special purpose computer. Having reference now to the drawings, in FIG. 1, there is shown a computer or data processing system generally designated by reference character 100. The computer 100 generally includes a central processing unit (CPU) 102 which is coupled with a memory 104. The memory 104 generally comprises both random access memory (RAM) and read-only memory (ROM) portions. The computer 100 further comprises at least one input/output device 106 connected with the CPU 102, such as a keyboard, a pointing device, a scanning device, or other available and useful devices. A display 108 is also connected with the CPU 102 to provide output to a user. Further, a storage mechanism 110 is provided to enable a user to store and retrieve data. Non-limiting examples of various storage mechanisms 110 that may be of use with the present invention include hard drives, optical drives (such as a compact disc or digital versatile disc), and magnetic drives (such as floppy disks, magnetic tapes, and Zip™ disks). Finally, a communications adapter 112 (such as a modem or network interface card) may also coupled with the computer 100 in order to allow communication with other computers 100. Note further that the data processing system may also comprise a network of computers either local to or remote from each other. In addition, the various components of the data processing system may all be local or may be distributed remotely from each other.

Although described at a major-component level, it is assumed that the computer 100 includes all of the necessary adapters to drive all of the peripherals necessary for use with the present invention. Additionally, various commercially available processors could be used in conjunction with the present invention, and could be programmed to operate the logic of the preferred embodiments.

Illustration and Flow Diagram of the Operation of the Present Invention Referring to FIG. 2, and FIG. 3, the operation of the present invention is described by means of an illustrative description and a flowchart. After starting operation of the present invention, represented by the begin 300 element in FIG. 3, an object 200 is divided into a plurality of patches, represented by element 202 in FIG. 2 and element 302 in FIG. 3. Note that methods for dividing the object into a plurality of patches are well-known in the literature, and that the patches may be of an arbitrary geometric shape (e.g., triangles, squares, etc.). Next, parameters are supplied for the incident field and the properties of the object as depicted by element 304 of FIG. 3. The incident field is depicted as element 204 in FIG. 2. The dashed lines surrounding elements 302 and 304 of FIG. 3 are intended to indicate that they are environmental to the present invention, as the patches 202 representing the object 200 and the parameters for the incident field 204 and the object 200 may simply be passed as inputs to a system incorporating the present invention. Optionally, however, the present invention may include a system for dividing the object 200 into a plurality of patches 202 and for more directly determining the associated parameters. In addition, although the object 200 is represented as a sphere and the incident field 204 is represented as a series of planar waves, these specifics are provided for illustrative convenience only. The object 200 may be any closed arbitrary shape and the incident field 204 may be any arbitrary time-harmonic field. Next, as shown by element 306 in FIG. 3, the system develops a discretized representation of a well-conditioned boundary integral equation (which takes the form of a well-conditioned matrix equation), modeling the interaction between the incident field 204 and the object 200. The well-conditioned linear system is then solved numerically to determine equivalent surface sources on the object, as depicted by element 308 in FIG. 3. After solving for the equivalent sources, the present invention may utilize the equivalent sources in order to determine electromagnetic fields resulting therefrom, as depicted by element 310 of FIG. 3. The dash-dot-dot line in element 310 is intended to indicate that it is optional. The electromagnetic scattering from the object is represented by element 206 in FIG. » 2.

By providing a well-conditioned linear system in the determination of the equivalent sources, the present invention helps to ensure that the solution to the linear system can be performed with fewer iterations and with increased numerical precision relative to conventional methods. In fact, without the use of the well-conditioned linear system, it is sometimes impractical to determine the equivalent surface sources accurately due to the limited numerical precision of one's computer.

Recording Medium In one embodiment, the present invention takes the form of a computer program product recorded on a recording medium. A disc, shown in FIG. 4 represents a recording medium 400 such as a compact disc (CD) or a digital versatile disc (DVD). Although depicted as a CD or DVD, any recording medium could be utilized, examples of which include magnetic discs and tapes, and other optical storage medium.

Further Details of the Present Invention

Well-Conditioned CFIE Formulation In accordance with features of the present invention, to provide improved computational efficiency and accuracy, modifications are made to the combined field integral equation (CFIE) for the solution of electromagnetic scattering problems in three dimensions. Unlike the conventional CFIE, the modified CFIE of the present invention is well-conditioned. It combines, a standard magnetic field integral operator with a pre-conditioned electric field integral operator, with the resulting modified CFIE being an integral equation of the second-kind that has no spurious resonances. As a result of the pre-conditioning, the number of iterations required to solve the CFIE stabilizes even as discretization of the surface is refined, whereas refinement of the discretization in conjunction with the conventional CFEE results in an ill-conditioned linear system and a growing number of iterations. The conventional CFEE, not being well-conditioned, suffers from two major drawbacks.

As mentioned previously, there are three principal reasons why well-conditioned boundary integral formulations are beneficial. The first is the fact that sparse decompositions of the integral operators of scattering theory, but not their inverses, exist, so iterative solvers must be employed. Well-conditioned systems of equations can be solved with few iterations. Second, using a fine discretization to resolve source variations or geometric detail on a sub-wavelength scale typically results in an ill-conditioned linear equation as well (also known as the "low- frequency" problem in computational electromagnetics). Third, because the process of generating a solution to the set of discretized equations can require a very wide range of values (in terms of magnitude), the resulting solutions may be quite inaccurate because a high degree of numerical precision may be required, which may exceed the accuracy permitted by the data processing system.

Only second-kind integral equations, or objects with similar spectral behavior (such as appropriately preconditioned differential equations) can be solved with fully controlled approximation error. The correct operators are the sum of a constant (or at least well- conditioned and easily invertible) operator and a compact operator.

Boundary integral operators of scattering typically violate this requirement in one of three ways: the spectrum may accumulate at zero; the operator may have an unbounded spectrum; or the operator may have small eigenvalues associated with resonances, often unphysical (e.g., "spurious resonances").

For electromagnetic scattering from perfectly conducting (PEC) surfaces, the standard boundary integral equations are the electric field integral equation (EFIE)

- n x E'"° = TJ (where inc stands for incident), (1)

and the magnetic field integral equation (MFEE)

where the integral operators Jand K are defined by TJ ≡ T(k)3

≡ ikn(x) V[VG(&; x, x') ^• J(x')]| and (3)

KJ s K(k)i ≡ - n(x) x jds' VG(k; x, x' ) x j(x' ) , (4) s where V denotes differentiation with respect to x, and n(x) is the unit normal to the surface at x. Note that the other terms follow the usual conventions: j ≡ Zn x H is the unknown surface current, E^mc and H^wc are the incident electric and magnetic fields, respectively, Z = Xμjε is the wave impedance, and G(k; x, x' ) = βxp(ikr) 14πr is the 3d Helmholtz kernel with r = |x - x'| being the distance separating field and source points. Harmonic time dependence e^~'°" is assumed.

The MFEE is a second-kind integral equation. Unfortunately, this equation is suitable for an unacceptably small class of electromagnetic problems. It is inapplicable to open surfaces, becomes ill-conditioned in the presence of geometric singularities, and suffers from spurious resonances. The EFIE has both a compact piece and a hypersingular piece (coming from the double gradient term). The spurious resonances of the MFIE may be eliminated by adding the EFIE to form a combined field integral equation (CFIE). The cost of doing so is the introduction of EFIE' s hypersingular piece, which spoils the conditioning for fine discretizations (or low frequencies).

It is known that a hypersingular integral operator and a first-kind integral operator are ideal preconditioners for each other, in the sense that the composition of the two has the spectral characteristics of a second-kind integral operator. A similar approach may be employed to analytically precondition the EEFE. In fact, as disclosed herein, the electric field integral operator J preconditions itself. Two issues are important for the successful application of this idea to closed bodies. First, only the local (or short distance) behavior of the preconditioner is important for asymptotic conditioning. Thus, the EFEE may be preconditioned by multiplying it by an electric field integral operator corresponding to an arbitrary wavenumber, real or complex; if the wavenumber has a positive imaginary part, the introduction of any additional resonances may be avoided. If the EFIE preconditioner reproduced the MFIE resonances, then the CFIE would also have them. Second, one must take care that the discretization of the product of the preconditioner and preconditioned operators preserves the correct spectral properties.

The present invention provides well-conditioned formulations for closed surfaces. Example numerical results are provided further below in order to demonstrate the advantages of the modified CFEE formulation over the conventional CFIE.

In the literature, integral operators have been constructed from the kernel for the Laplace and Helmholtz equations in 2d. Through the construction of these operators, it has been observed that the product of a first-kind operator, constructed from an undifferentiated kernel, and a hypersingular operator, constructed from a twice differentiated kernel, has the desirable spectral characteristics of a second-kind operator. Since the EFIE integral operator Jhas both of these, it may be seen that the composition of two such operators T² ≡ T o T would include a constant operator and a compact operator. Also of concern is that the product of hypersingular components would produce another hypersingular operator. However, it is easy to see that the rotation operation n x in the definition of T, as shown in equation (3) above, annihilates the component of the surface vector field normal to the surface and ensures that the product of the two hypersingular operators is identically equal to zero. Indeed, applying the hypersingular component of the second T operator to an arbitrary tangential surface vector function f(x') produces a surface gradient function, n x V(z)(x) x, x') ■ f(x') , (5)

(note the dot-product operator, ^•) which the hypersingular component of the first J operator, in turn, annihilates (for closed surfaces) by virtue of the identity

V_s - [nxV^(x)] = 0- (6) with V_s denoting the surface gradient operator on S; the identity shown in equation (6), the surface analog of the 3d identity V • [V x can be found and is valid for any sufficiently smooth function φon S. It follows immediately from equations (3), (5), and (6) that J² behaves as a second-kind integral operator.

In this specification, the spectral properties of the EFIE and MFEE integral operators and combinations thereof for the PEC sphere, a simple 3 d target for which the spectral properties of these operators are known analytically, are discussed. However, this discussion may be expanded to incorporate other object shapes. The use of the PEC sphere herein is intended to provide a non-limiting example of an object to which the present invention may be applied. A complete set of basis functions on the surface of a sphere of radius a is given by the vector spherical harmonics,

*AΘ,φ) ≡ . , * x VY_lm(θ,φ) and (7) ιφ(l + l)

V_!m(θ,φ) ≡ n x X_lm(θ,φ) , (8) defined here in terms of the scalar spherical harmonics Y,_m(θ,φ) ■

The result of applying T and K⁺ ≤ (K + — ) to each basis function is:

TOI luj = l^{~ j;} (k ^{ka^a)⁾M ^M',^/ ( ^{k^ka^a)⁾X ^uA and

where j, and H, are Riccati-Bessel and (first-kind) Riccati-Hankel functions of order /, and k is the wavenumber associated with the kernel of each integral operator. The Riccati-Bessel and Riccati-Hankel functions are defined in terms of spherical Bessel and Hankel functions j_t(x)

S,(x) ≡ xj,(x) and (11) M,(x) ≡ xh Xx) (12) Although chosen basis functions X_!m and U_/m are not eigenfunctions of the operator T(k), they are eigenfunctions of T² (k) ≡ T(k) ° T(k) :

- J, (ka) H, (ka) J, (ka) H| (ka) (13)

The operator T²^) has a bounded spectrum, since, in the limit of large /, its eigenvalues accumulate at -1/4 (a result which follows from the asymptotic properties of / and h ¹).

However, as is evident, from equations (10) and (13), the operator T²(k) also shares resonances (at the zeros of $', (ka) for the X_/m modes, and at the zeros of S, (ka) and for the U_/m modes) with the MFIE operator K⁺(k) . This is also evident from the identity T²(k) = - - K²(k) = -K^~(k) o K⁺(k), (14)

where K^~ ≡ — K . Therefore, although T² (k) is a second-kind integral operator, it is not a suitable component of a resonance-free combined field integral equation for closed bodies.

As previously mentioned, boundedness of the spectrum of the product of two EFIE operators (of the form shown by equation (3)) is assured if they have the same short-distance behavior, a condition that does not require the two operators to share the same wavenumber (propagation constant). If EFIE operators are chosen with different wavenumbers, T(k ) and T(k₂) , it is possible to simultaneously obtain a bounded product and avoid MFIE resonances.

The following analysis indicates that ik is a particularly good choice for the wavenumber in the preconditioning operator (assuming that the wavenumber k is real). The eigensystem for T (ik) T(k) on a sphere is

T C _W\ o ^ frzJM J - - J^f>^ ^{ika) _π ^m^^{'! {i} _a ^k ₎ ^a) _j ^J _; ' _{ ⁽ _k ^fa _a) ⁾ _H ^H,' _{ ⁽ _k ^fa _a) ⁾ _v . n15 <) .

It is straightforward to show (given the properties of/_/ and b,⁽¹⁾ ) that the eigenvalues of

T (ik) T(k) accumulate at — and — for the X/_m and U_/m eigenmodes, respectively, and that

T (ik) o T(k) does not share any resonances with the MFIE operator K⁺ (k) . Since J (ik) ° T(k) is a second-kind integral operator (in the sense described below) and does not share any resonances with K⁺(k) , it is possible now to develop the well-conditioned CFIE operator. The simplest form of such an operator is

T (ik)oT(k) + cK⁺(k), (16)

where is a constant to be chosen. In developing this modified CFIE operator, the EFIE portion has been modified before being added to the MFIE portion (which is already a second-kind integral operator). The same applies to the excitation side of the equation. The resulting CFIE may be written as - T(ik) (n x E'^nc) + aZn x H^mc = [τ(ik) ° T(k) + aK⁺( . (17)

The eigensystem for the CFEE operator in (16) is

If one chooses a = +1 then, as a function of the argument ka, these eigenvalues have no zeros. For = +1 , they circle the origin of the complex plane.

Other well-conditioned (modified) CFIE operators can be devised, for example, by preconditioning the MFIE part before combining it with the preconditioned EFIE part. Two example forms of alternative well-conditioned (modified) CFIE operators are as follows: and T(ik) ° T(k) + nx K⁺(ik) ° n x K k) . (20)

It has been observed that the numerical behavior of the three modified CFIE operators (in (16), (19), and (20)) is similar.

By virtue of the analysis on spheres, the modified CFEE operators in (16), (19), and (20) have been shown to be second-kind and resonance-free. However, given that the asymptotic behavior of the eigenvalues on a smooth surface stems from the short distance behavior of the kernel, it may be seen that the asymptotic behavior of the various operators on spheres should also obtain for any closed surface that can be obtained by smooth deformation of a sphere. The numerical results presented further below support this argument. Results are also presented for a cube, which, like many targets of practical interest, has geometric singularities. These suggest that the modified CFEE formulations should be well-conditioned for a wide class of closed surfaces.

Discretization and an Alternate Form of the Modified CFEE Formulation

Discretization of the modified (well-conditioned) CFIE formulation utilized in the present invention may take place by any appropriate method. Suitable methods are common in the literature, non-limiting examples of which include the Nystrδm discretization method and the method of moments. A particularly useful reference in regard to the Nystrδm discretization method, which may be used in conjunction with the present invention may be found in reference [13], in the Background, titled "Numerical solution of the Helmholtz equation in 2d and 3d using a high-order Nystrδm discretization" by Canino et al. Note also, incorporating the approach of the present invention into existing "fast" solvers is completely straightforward.

There are several ways to discretize the product operator T (k_x) o T(k₂) . The simplest and most straightforward approach, multiplying the discretized representations of the individual operators, can lead to numerical difficulties. The reason is that it is relatively difficult to make the discretized representations of the hypersingular part of each operator sufficiently accurate (especially for high-spatial-frequency eigenmodes) to numerically effect the cancellation that obtains analytically.

Effective discretizations of T (k_{) ° T(k₂) can be obtained either by discretizing the product operator directly or by reformulating the product operator to eliminate the product of hypersingular operators. The second approach has been implemented using a reformulated product operator that eliminates all instances of hypersingular operators. A short derivation of the reformulated equation is provided below.

The first step toward obtaining a more useful form of the product operator J ( c.) o T(k₂) is to separate each integral operator into its singular and hypersingular components. Introducing the abbreviations,

T₂ = T(k₂) , (22) and substituting, the following can be written:

where,

J_m ^s J --- n(x) x \ds'G(k_m ; x, x' ) J(x' ) , and (25)

T J _≡ (n(x) x, x') ■ J(x') . (26)

The product operator T ° T₂ can be expanded into four terms. Each of the two cross teπns, 7j^s o T" and T ° T₂ , can be transformed (by Stokes 's theorem) into a product of new, single- gradient integral operators on S plus a line integral around the boundary of S. The term formed by the product of the hypersingular integral operators, T" ° T₂ , reduces to a line integral. The result is further simplified by noticing that two of the three line integrals, when applied to J, can be combined into a single term whose argument is identical to the incident electric field E'"^e by virtue of equation (1).

The next step is to reformulate the excitation side of the equation, taking advantage of the fact that the incident wave obeys Maxwell's equations. By applying Stokes's theorem, the term T [n(x' ) x E'"^c (x' )] may be rewritten as the sum of a single-gradient integral operator on

V x E'"^c(x') and a line integral that exactly cancels the line integral involving Ε"^κ on the other side of the equation. A further simplification follows from Faraday's Law, V x E = iωμR .

The final result for the analytically preconditioned EFEE with reformulated integral operator product is

•zA_ιr_ι xE^,) o T₂ ^L)i

(27) where the various integral operators are defined by: τ ψ ≡ n(x) x x, x')φ(x'), (28)

T^Pφ ≡ n(x) x J n(x') x VG(k_m; x, x^<)φ(x') , (29) s

T^Ef ≡ jd_S' G(k_m, x, x') - f(x'), (30) s J_m ^rf ≡ n(x) ■ x, x') x f (x') , (31)

J_m ^sf _≡ n(x) x jds' G(k_m; x, x'X(x') , and (32) s

T*φ ≡ n(x) x x, x')φ(x') , and (33) where m=l, 2. Note that T , Tf , and T^E map scalar functions to surface vector functions, whereas J_m nd ^ do the reverse. The operator on the right hand side of equation (27) maps surface vector functions into surface vector functions.

Now, closed surfaces will be discussed, and it may be observed that J_t o T₂ behaves like a second-kind integral operator.

If S is a closed surface, the term T T₂ J vanishes and equation (27) simplifies to

- / (n x E') - Z^-T_y ^a(n ■ H') = S(k_y, k₂)J, (34) k₂ where

S₁₂ ≡ S(k„ k₂) ^T^a T₂ ^T + _Tl ^/' o T₂ ^L - k_lk₂T₁ ^s AT . (35) kγ k₂

Now, several features of S_/2 will be discussed.

First, all of the individual operators that comprise S₇₂ involve kernels with one or no gradients on the Helmholtz Green's function G. All such integral operators are bounded.

Second, the eigenvalues of the integral operator S ₂ do not accumulate at the origin. This may be demonstrated by examining its three components, T" o T₂ , Tf T₂ , and T o T₂ . The operator T" T₂ ^T is a second-kind integral operator for the transverse (divergence-free) component of J, and is identically zero for the longitudinal (irrotational) component of J. Likewise, the operator

TJ³ oT₂ is a second-kind operator for the longitudinal component of J, and is identically zero for the transverse component of J. Since any surface current distribution can be decomposed into k k longitudinal and transverse components, the sum —T" °T₂ +~T °T₂ is a second-kind k k₂ integral operator; subtracting °T₂ , a compact operator, does not change this result. As observed previously, resonance sharing may be avoided by setting k = ik and k₂ = k. In this

case, the eigenvalues of S_/2 accumulate at two points, + — , rather than only at — .

Third, the spectrum of Sy₂, after discretization, is bounded and includes accumulation points at the expected locations. However, an accurate discretization will have zero (or very small) eigenvalues wherever the EFIE operator T(k₂) has a resonance. Thus, it has to be combined with an appropriate discretization of the MFIE operator, to obtain an effective discretization of the CFEE.

Finally, it should be noted that equation (34) is manifestly insusceptible to the "low-frequency" problem that plagues the EFIE. Since the well-conditioned behavior of S_/2 comes from the k k composite operators —Tf °T₂ and — J °T₂ , both of whose prefactors have modulus unity Λ_j k₂

(assuming ° T₂ tends to zero as k — > 0 , the full operator S_/2 remains well conditioned in the limit of low frequency.

In summary, although the operators 7j o T₂ and S_/2 have identical spectral properties for closed bodies, it is easier to construct an accurate discretization for S;₂ because it is composed of less singular integral operators. Matrix representations of S_/2 have bounded spectra, but also suffer from spurious resonances inherited from the EFIE operator T(k₂) . These resonances can be eliminated by combining S/₂ with K⁺(k₂) (or the modified MFIE operators in equations (19) and (20)). The result is a well-conditioned system of linear algebraic equations. An Alternate General Formulation of the Modified CFIE The boundary integral equation used for the present invention may be formulated more generally, as follows:

(E excitation term) + ^(H excitation teπn) = {(E integral operator term) + γ (H integral operation term)}

where, γ is a complex constant,

J represents the unknown surface source distribution, and the (E excitation term) and (E integral operator term) are of a formulation selected from the group consisting of A and B, with

A representing the combination of

(E excitation term)--- T(k₁)(- n(x) x E'^∞(x)) and

(E integral operator term) ≡ T(k ) T(k₂ ) , and with B representing the combination of

(E excitation term) ≡

^J^s(^)(-n(x) x E'^∞(x)) - ^T^a(kJ)(zn( ) ■ H'^∞(x)) and

(E integral operator term) ≡

^T k_l) oT^T(k₂) + ^T k₂) o T^L(k₂) -k_lk₂T^s(k_l) o T^s(k₂) ,- and where Λ₂ the (H excitation term) and (H integral operation term) are of a formulation selected from the group consisting of A', B', and C, with A' representing the combination of

(H excitation term) ≡ Zn(x) x H'"^c(x) and

(H integral operation term) ≡ K⁺(k₂) ; B' representing the combination of

(H excitation term)--; K k_λ)(Zn(x) H'^∞(x)) and

(H integral operation term) ≡ K⁺ (k o K ⁺ (k₂) ; and C representing the combination of (H excitation term)--- n x iT (£,)(Zn(x) x H'^∞(x)) and

(H integral operation term) ≡ n x iC⁺(&₁)]o n x K⁺(k₂)\, wherein;

T(k) ≡ ikn(x) V[VG(&; x, x') ^■ J(x')] ,

K⁺(k) ≡- + K(k) , and

K(k) ≡ -n(x) x jds'VG(k; x, x') x J(x') , and s wherein J ≡ Zn x H is the unknown surface current; ζ is a complex number;

E^mc and H^mc are the incident electric and magnetic fields, respectively;

Z = ■sjμjε is the wave impedance; G(k; x, x') = exp(ikr) 14πr is the 3d Helmholtz kernel, where r = |x - x is the distance separating field and source points with harmonic time dependence e^~'^ωt assumed, and where

T^a(k)φ ≡ n(x) x jds'VG(k; x, x')φ(x') , s

T^p(k)φ ≡ n(x) x Jώ'n(x') x VG(k; x, x')φ(x') , s T^L(k)f ≡ jds G(k; x, x')f(x') , s

T^τ(k)f ≡ n(x) ■ jds'VG(k; x, x') x f (x') , and s

T^s(k)f ≡ n(x) x jds'G(k; x, x')f(x') , and where s

V and V denote differention with respect to the field point x and source point x' , respectively, and n(x) and n(x') are outward-pointing unit vectors normal to the target surface S at x and x', respectively, and i ≡ V-l ; k₂ is the propagation constant in the medium external to S; k_x is a complex number chosen to avoid resonances. This general form may be modified to provide the preferred form of the modified CFIE, as previously mentioned, (ik) is substituted for (k) in the various parts of the equation as defined above. Specifically, for example, the term

T^a(k)φ ≡ n(x) x ϊds'VG(k; x, x')φ(x') is simply changed to s

T^a(ik)φ ≡ n(x) x jds'VG(ik; x, x')φ(x') . s

The same type of change is applied to all of the other terms.

Numerical Results Demonstrating the Advantages of the Modified CFIE Over the Conventional CFIE

Below, the computational performance of the conventional CFEE (referred to below as CFIEl), - n x (n x E') + Zn x H' = [n x T(k) + K⁺(k)}j (36) is compared with the preconditioned CFEE (refeπed to below as CFIE2) kT^s (ik)n xE' + iZT^a (ik)n ■ H' - Zn x H'

= {/ [- T^a (ik) o T^τ (k) + T^β (ik) o T^L (k) - k²T^s (ik) o T^s (k)] - K⁺ (k)} (37) produced by combining equation (17) (with a = -1) and equation (34) (with k_λ = ik and k₂ = k ). Individual operators in these equations were discretized using a high-order Nystrδm scheme. In all cases, the wave impedance Z was set to unity.

Three examples are provided below to demonstrate the operation of the present invention. The first example shows how the condition number of each operator, defined as the ratio of the largest to smallest singular values, depends on the fineness of discretization. Table 1, below, lists the condition number (CN) of the matrix representing each CFIE operator as the size of the sphere decreases. In all cases, the same discretization was used, created by placing a 6-point quadrature rule on each of the 80 nearly identical patches that cover the sphere, for a total of 960 unknowns. As the sphere radius decreases, the condition number for the CFIE2 integral operator stabilizes at about 2, whereas the condition number for the CFIEl continues to grow in inverse proportion to the radius.

Table 1 : Condition number of CFEE matrices for shrinking PEC spheres

The second test compares iterative solver performance for the modified CFEE with that for the conventional CFEE. The target geometry consists of two PEC spheres, one with a radius of 2/2 , the other set at a resonant radius, namely, the first zero of S[ (2πr/λ) or r « 0.436674571. The spheres are separated by a λl\ 00 gap. The patches near the gap were subdivided by a factor of 10 to adequately resolve the currents, which vary rapidly there. Table 2 compares the iteration counts and radar cross section (RCS) errors for several discretizations. The iterations columns list the maximum number of iterations a conjugate gradient squared (CGS) route required to reach a residual error of 10^~3. A solution computed from a substantially more refined discretization provided an accuracy reference. The stated error is the root mean squared (RMS) value of the difference between the monostatic φφ RCS of the comparison solution and the reference solution at 181 angles. For identical discretizations, the two methods had about the same error. The data show a dramatic difference, however, in the iteration count behavior of the two methods in response to discretization refinements.

Table 2: Iteration count and solution error vs. discretization for two PEC spheres.

The third test also compared iterative solver performance for the two CFIE formulations. In this case, the target was a cube of size \λ (where λ is the wavelength of the radiation). Numerical results are presented for five different discretizations, the first of which was derived from a mesh (i.e., a set of patches) obtained by dividing each face into four squares. The second mesh was constructed from the first by subdividing each square into four smaller squares. The third mesh was constructed from the second by subdividing edge-touching patches in half along a line parallel to the edge; patches adjacent to two edges (i.e., corner patches) were divided into quarters. Meshes for the fourth and fifth discretizations were constructed by recursively applying the procedure by which the third mesh was constructed from the second. This process, known as patch tapering, is useful for resolving the source singularities that arise in the vicinity of geometric singularities. It also puts stress on the conventional CFEE because points near the edges get close together. Table 3 lists the maximum and average number of iterations the CGS routine required to obtain solutions for 92 independent excitations to a residual error of 10^"3. The total number of unknowns is the result of using a 9-point quadrature rule on each square or rectangular patch. The iteration count for CFIE2 grows very slowly with increasing taper depth, whereas for CFIEl it increases steadily as expected.

Table 3: Iteration count vs. taper depth for \λ PEC cube.

Appendix - Second-Kind Integral Operators The standard definition of a second-kind integral operator is an operator of the form I + K, (38) where λ is a constant, / is the identity, and K is a compact operator. In scattering theory, one encounters operators of the form

W +^ +K, (39) where λ and λ are constants and Pj and P₂ are orthogonal projection operators such that

P_λ +P₂ =I . (40) Operators of the form shown in equation (39) possess most of the desirable properties of second- kind integral operators. Although somewhat non-standard, such expressions are referred to as second-kind operators throughout this description.

Glossary of Terms and Conventions

Composition Operator - The composition operator used herein, and represented by the symbol, o , describes the following operation. Given a first integral operator, A ≡ x') , and a second integral operator, B ≡ \ds"K^B(x', x") , then the composed s operator A ° B ≡ x") .

Equivalent surface sources - According to a well-known mathematical theorem (the "equivalence principle"), the fields outside a volume V can be completely specified by a distribution of equivalent sources on the boundary of the volume S. If S corresponds to the surface of a metal target (which approximates a perfect electrical conductor), the equivalent sources correspond to true, physical surface currents.

Means - Generally, for purposes of this invention, the term refers to a computer program or subprogram, or a portion of a computer program operative for performing the task that the means is for. The means is generally considered to be language and computer platform independent, and in the form of computer instructions (or code) that may be either in the form of software or hardware.

Object - This term object refers to primarily any arbitrarily shaped closed surface for which equivalent surface sources and the resultant electromagnetic scattering are to be determined.

Objects of interest could include, as non-limiting examples, aircraft or parts thereof, automobiles, tanks, boats, missiles, and fixed structures. The equations set forth herein are developed for perfectly electrically conducting (PEC) materials or approximately perfectly conducting materials such as metals. However, extension to surface scattering with more general boundary conditions is considered routine and within the scope of the present invention. The extension to an interface between two or more dielectrics, for example, is straightforward, with the resulting operators having condition numbers slightly lower than the case described herein. Although primarily developed for closed surfaces, the approach can potentially be extended for the case of open, perfectly conducting surfaces.

Claims

CLAIMS What is claimed is:

1. A computer program product, used in a computer system having a processor, for modeling an electromagnetic response of an arbitrarily shaped three-dimensional object to an arbitrary time-harmonic incident field by means of a well-conditioned boundary integral equation (BIE), the computer program product comprising: a. a recording medium; b. means, recorded on the recording medium for producing a discretized representation of a well-conditioned BIE to provide a well-conditioned, finite- dimensional linear system whose solution approximates a distribution of equivalent surface sources on a model caused by the arbitrary time-harmonic incident field; and c. means, recorded on the recording medium, for solving the linear system to determine values of the equivalent surface sources; whereby solution of the well- conditioned finite-dimensional linear system provides improved computational efficiency and solution accuracy.

2. A computer program product as set forth in claim 1, wherein the computer program product further comprises means, recorded on the recording medium for modeling electromagnetic scattering from the three-dimensional object due to the time-harmonic incident field as determined from the equivalent surface sources.

3. A computer program product as set forth in claim 2, wherein the computer program product is tailored for determining the electromagnetic scattering from a time-harmonic incident field in the form of a radar wavefront and where the three-dimensional object represents a component of a vehicle.

4. A computer program product, as set forth in claim 1, wherein the BEE produced is in a form of: (E excitation term) + γ (H excitation term) = {(E integral operator term) + χ(H. integral operation term)}.

where, γ is a complex constant, J represents the unknown surface source distribution, and the (E excitation term) and (E integral operator term) are of a formulation selected from the group consisting of A and B, with A representing a combination of

(E excitation term) ≡ T(kJ)(- n(x) x E'"^c(x)) and (E integral operator term) ≡ T(k_x ) T(k₂ ) , and with B representing a combination of (E excitation term) ≡

ikJT^s(k_x)(→(x) x E'"^c(x)) - rh, T^a( n(x) • H'^rø(x)) and

(E integral operator term) s

; and where the (H excitation term) and (H integral operation term) are of a formulation selected from the group consisting of A', B', and C, with A' representing a combination of (H excitation term) ≡ Zn(x) x H'"^c (x) and

(H integral operation term) ≡ K⁺(k₂) ; B' representing a combination of

(H excitation term)--- ⁺(/_I)(Zn(x) x H'"^c(x)) and

(H integral operation term) ≡ K⁺(k_l) K⁺(k₂) ; and C representing a combination of

(H excitation term)≡ n x K⁺ (k_{)(Zn(x) x H'^∞(x)) and

(H integral operation term) ≡ [nx ⁺(A:₁)Jo [n x K⁺(k₂)\, wherein;

T(k) ≡ ikn(x) x jds'i, G(k; x, x')J(x) + x, x') • J(x')]}

K⁺(k) ≡ -+K(k) , and

L K(k) ≡ -n(x) x jds'VG(k; x, x')x J(x') , and s wherein J s Zn x H is the unknown surface current; ζ is a complex number;

E^mc and H"^!C are incident electric and magnetic fields, respectively;

Z = sjμ/ε is the wave impedance;

G(k; x, x') = exp(/£r)/4πris a 3d Helmholtz kernel, where r = points with harmonic time dependence e^~'°" assumed, and where

T^a(k)φ ≡ n(x) x x, x')φ(x') ,

T^β(k)φ ≡ n(x) x J -s'n(x') x VG(k; x, x')φ(x') , s

T^L(k)f ≡ jds'VG(k; x, x')f(x') , s

T^τ(k)f ≡ n(x) • fø'VG(k; x, x') x f (x') , and s

T^s(k)f ≡ n(x) x x, x')f (x') , and where

V and V denote differentiation with respect to the field point x and source point x' , respectively, and n(x) and n(x') are outward-pointing unit vectors normal to the target surface S at x and x' , respectively, and i ≡ V-l ; k₂ is the propagation constant in a medium external to S; k_x is a complex number chosen to avoid resonances.

5. A computer program product as set forth in claim 4, wherein the computer program product further comprises means, recorded on the recording medium for modeling electromagnetic scattering from the three-dimensional object due to the time-harmonic incident field as determined from the equivalent surface sources.

6. A computer program product as set forth in claim 4, wherein the BEE produced is in a form of: kT^s (ik)n x E'"^c + iZT* (ik)ιι ^■ H'"^c = { [- T^a (ik) ° T^τ (k) + T^p (ik) ° T^L (k) - k²T^s (ik) and wherein the terms of (ik) are developed by substituting (ik) for (k) in the terms of the formulation set forth in claim 4.

7. A computer program product as set forth in claim 6, wherein the computer program product further comprises means, recorded on the recording medium for modeling electromagnetic scattering from the three-dimensional object due to the time-harmonic incident field as determined from the equivalent surface sources.

8. An apparatus for modeling the electromagnetic response of an arbitrarily shaped three- dimensional object to an arbitrary time-harmonic incident field by means of a well- conditioned boundary integral equation (BIE), the apparatus comprising: a. a computer system having a processor, a memory connected with the processor, and an input/output device connected with the processor for receiving parameters representing an arbitrarily shaped three-dimensional object and an arbitrary time- harmonic incident field; b. means, operating in the processor and the memory of the computer for utilizing the parameters representing the arbitrarily shaped three-dimensional object and the arbitrary time-harmonic incident field for producing a discretized representation of a well-conditioned BIE to provide a well-conditioned finite dimensional linear system whose solution approximates a distribution of equivalent surface sources on a model caused by the time-harmonic incident field; and c. means, operating in the processor and the memory of the computer for solving the linear system to determine values of the equivalent surface sources; whereby solution of the well-conditioned finite dimensional linear system provides improved computational efficiency and solution accuracy.

9. An apparatus as set forth in claim 8, wherein the apparatus further comprises means, operating in the processor and the memory of the computer for modeling electromagnetic scattering from the three-dimensional object due to the time-harmonic incident field as determined from the equivalent surface sources.

10. An apparatus as set forth in claim 9, wherein the apparatus is tailored for determining the electromagnetic scattering from a time-harmonic incident field represented in the form of a radar wavefront and where the three-dimensional object represents a component of a vehicle.

11. An apparatus as set forth in claim 8, wherein the BIE produced is in a form of:

(E excitation term)+ ^(H excitation term) = {(E integral operator term)+ ^(H integral operation term)}.

where, γ is a complex constant,

J represents an unknown surface source distribution, and the (E excitation term) and (E integral operator term) are of a formulation selected from the group consisting of A and B, with A representing a combination of (E excitation term) ≡ T(k_l)(- n(x) x E^i∞(x)) and

(E integral operator term) ≡ T(k_λ ) T(k₂ ) , and with B representing a combination of (E excitation term) ≡

ϊkJJ" k_x)(→(x) x E'"^c(x)) - ζ^-J^α(^)(Zn(x) • H^mc(x)) and

(E integral operator term) =

o T^L(k₂) - ° T^s(k₂) ; and where the (H excitation term) and (H integral operation term) are of a formulation selected from the group consisting of A', B', and C, with A' representing a combination of (H excitation term) ≡ Zn(x) x H'"^c(x) and

(H integral operation term) ≡ K⁺ (k₂ ) ; B' representing a combination of

(H excitation term)-≡ ⁺(^)(Zn(x) x H'^∞(x)) and

(H integral operation term) ≡ K⁺ (kj) ° K ⁺ (k₂) ; and C representing a combination of

(H excitation term) ≡ n x ⁺(^)(Zn(x) x H'^∞(x)) and

(H integral operation term) ≡ n x if⁺(A:₁)Jo [n x K⁺(k₂)\, wherein;

T(k) ≡ ikn(x) + v[VG(&; x, x') • J(x')]}

K⁺(k) ≡ - + K(k) , and

K(k) ≡ -n(x) x | -sNG(it; x, x') x J(x') , and s wherein J ≡ Zn x H is the unknown surface current; ζ is a complex number;

E'"^c and H'"^c are incident electric and magnetic fields, respectively; Z = τ]μ/ε is a wave impedance;

G(k; x, x') = exp(ikr) 1 %r is a 3d Helmholtz kernel, where r - |x - x is a distance separating field and source points with harmonic time dependence e^~ assumed, and where

T^α(k)φ ≡ n(x)x jd 'VG(k; x, x')φ(x') , s T^p(k)φ ≡ n(x) x jds'n(x') x VG(k; x, x')φ(x') , s

T^τ(k)i ≡ n(x) • jds G(k; x, x*) x f (x') , and s

T^s(k)f ≡ n(x) x x, x')f (x') , and where

V and V denote differentiation with respect to the field point x and source point x' , respectively, and n(x) and n(x') are outward-pointing unit vectors normal to the target surface S at x and x' , respectively, and / ≡ V-ϊ ; k₂ is the propagation constant in a medium external to S; k_x is a complex number chosen to avoid resonances.

12. An apparatus as set forth in claim 11, wherein the apparatus further comprises means, operating in the processor and the memory of the computer for modeling electromagnetic scattering from the three-dimensional object due to the time-harmonic incident field as determined from the equivalent surface sources.

13. A computer program product as set forth in claim 11, wherein the BEE produced is in a form of: kT^s (ik)n x E^ta + iZT^a(ik)n • H'"^c = {i [- T^a (ik) ° T^τ (k) + and wherein the terms of (ik) are developed by substituting (ik) for (k) in the terms of a formulation set forth in claim 11.

14. An apparatus as set forth in claim 13, wherein the apparatus further comprises means, operating in the processor and the memory of the computer for modeling electromagnetic scattering from the three-dimensional object due to the time-harmonic incident field as determined from the equivalent surface sources.

15. A method for modeling an electromagnetic response of an arbitrarily shaped three- dimensional object to an arbitrary time-harmonic incident field by means of a well- conditioned boundary integral equation (BIE), the method operating on a computer system having a processor, a memory connected with the processor, and an input/output device connected with the processor for receiving parameters representing an arbitrarily shaped three-dimensional object and an arbitrary time-harmonic incident field, the method comprising the steps of: a. utilizing the parameters representing the arbitrarily shaped three-dimensional object and the arbitrary time-harmonic incident field for producing a discretized representation of a well-conditioned BEE to provide a well-conditioned finite dimensional linear system whose solution approximates a distribution of equivalent surface sources on a model caused by the arbitrary time-harmonic incident field; and b. solving the linear system to determine values of the equivalent surface sources; whereby solution of the well-conditioned finite-dimensional linear system provides improved computational efficiency and solution accuracy.

16. A method as set forth in claim 15, wherein the method further comprises the step of modeling electromagnetic scattering from the three-dimensional object due to the time- harmonic incident field as determined from the equivalent surface sources.

17. A method as set forth in claim 16, wherein the time-harmonic incident field is a radar wavefront and wherein the three-dimensional object represents a component of a vehicle, in order to determine the scattering of radar waves from a vehicle component due to in incident radar wave.

18. A method as set forth in claim 15, wherein the BIE produced is in a form of:

(E excitation term)+ γ (H excitation term) = {(E integral operator term)+ ^(H integral operation term)}.

where, γ is a complex, J represents an unknown surface source distribution, and the (E excitation term) and (E integral operator term) are of a formulation selected from the group consisting of A and B, with A representing a combination of

(E excitation term) ≡ T(k_l)(- n(x) x E^mc(x)) and (E integral operator term) ≡ T(k_λ ) T(k₂ ) , and with

B representing a combination of (E excitation term) ≡

ikJT^s(k )(-ιv(x) x E'"^c(x)) - ■ H"^κ(x)) and

(E integral operator term) ≡

ζ η^T^a(k_x) o T^τ(k₂) + ^T k₂) o T^L(k₂) - k_xk₂T^s(iX) ° T^s(k₂) ; and where

the (H excitation term) and (H integral operation term) are of a formulation selected from the group consisting of A', B', and C, with A' representing a combination of

(H excitation term) ≡ Zn(x) x H'"^c(x) and (H integral operation teπn) ≡ K⁺ (k₂ ) ; B' representing a combination of

(H excitation term)--- ⁺(£.)(Zn(x) x H"^,c(x)) and

(H integral operation term) ≡ K⁺(k₁) ° K⁺(k₂) ; and C representing a combination of

(H excitation term)≡ n x i ⁺(Jt₁)(Zn(x) x H"^,c(x)) and

(H integral operation term) ≡ n x ⁺(A;₁)Jo n x K⁺(k₂)\, wherein;

T(k) ≡ ika(x) x x, x') ■ J(x')]}

K⁺(k) ≡ -+K(k) , and

K(k) ≡ -n(x) x x, x') x J(x') , and wherein J == Zn x H is an unknown surface current; ζ is a complex number;

E'^∞ and H^mc are incident electric and magnetic fields, respectively;

Z = [μ/ε is a wave impedance;

G(k; x, x') = exp(ikr)/4πris a 3d Helmholtz kernel, where r = |x - x is a distance separating field and source points with harmonic time dependence e^~'°" assumed, and where

T^α(k)φ ≡ n(x) x jds'VG(k; x, x')φ(x') , s

T^β(k)φ s n(x) x jds'n(x') x VG(k; x, x')φ(x') , s

T^L(k)i = jds'VG(k; x, x')f (x') , s T^τ (k)f ≡ n(x) • jds'VG(k; x, x') x f (x') , and s

T^s(k)f ≡ n(x) x x, x')f(x') , and where

V and V denote differentiation with respect to the field point x and source point x' , respectively, and n(x) and n(x') are outward-pointing unit vectors normal to the target surface S at x and x' , respectively, and / ≡ V-ϊ ; k₂ is the propagation constant in a medium external to S; /c- is a complex number chosen to avoid resonances.

19. A method as set forth in claim 18, wherein the method further comprises the step of modeling electromagnetic scattering from the three-dimensional object due to the time- harmonic incident field as determined from the equivalent surface sources.

20. A method as set forth in claim 18, wherein the BIE produced is in the form of: kT^s (ik)n x E'"^c + iZT^a(ik)n • ET = and wherein the terms of (ik) are developed by substituting (ik) for (k) in the terms of the formulation set forth in claim 18.

21. A method as set forth in claim 20, wherein the method further comprises the step of modeling electromagnetic scattering from the three-dimensional object due to the time- harmonic incident field as determined from the equivalent surface sources.