1. Introduction

SuperNEC (SNEC) is an object-oriented version of the FORTRAN program NEC-2. This technical manual is therefore very similar to the NEC-2 manual. In fact, many sections have been copied verbatim from the NEC-2 manual. The reason this document has been produced (as opposed to referring the reader to any readily available NEC-2 manual) is because new theory will be added to SNEC and some features that exist in NEC-2 will not be implemented in SNEC. The easiest way of keeping an up to date reference of the theory implemented in SNEC is to start with the an electronic version of the NEC-2 manual and modify the document as the SNEC code evolves. Many extensions have been made to SNEC in prototype form. These features include hybridisation with UTD, fast iterative solvers, MBPE amongst other features. When these extensions are formally incorporated into the SNEC program, then the theory behind the extension will be added to this manual.

2. The Integral Equation for Free Space

The SNEC program uses an electric-field integral equation (EFIE) to model the electromagnetic response of general structures. The EFIE can be used to model both wire structure and surfaces. Surfaces are represented by wire grids and have had reasonable success for far-field quantities but with variable accuracy for surface fields. The EFIE and its derivation are outlined in the following sections.

2.1. The electric field integral equation (EFIE)

The form of the EFIE used in NEC follows from an integral representation for the electric field of a volume current distribution .

(1 )

where

and the time convention is .

is the identity dyad . When the current distribution is limited to the surface of a perfectly conducting body, equation (1) becomes :

(2)

where

is the surface current density.

The observation point is restricted to be off the surface so that .

If approaches as a limit, equation (2) becomes :

(3)

where the now represents a principal value integral. It is a principle value integral since is now unbounded.

An integral equation for the current induced on by an incident field can be obtained from equation (3) and the boundary condition for :

(4)

where

is the unit normal vector of the surface at

is the field due to the induced current.

Substituting equation (3) for yields the integral equation :

(5)

The vector integral in equation (5) can be reduced to a scalar integral equation when the conducting surface is that of a cylindrical thin wire, thereby making the solution much easier. The assumptions applied for a thin wire, known as the thin-wire approximation, are as follows:

1. Transverse currents can be neglected relative to axial currents on the wire.

2. The circumferential variation in the axial current can be neglected.

3. The current can be represented by a filament on the wire axis.

4. The boundary condition on the electric field need be enforced in the axial direction only.

These widely used approximations are valid as long as the wire radius is much less than the wavelength and much less than the wire length. An alternate kernel for the EFIE, based on an extended thin-wire approximation in which condition 3 is relaxed, is also included in NEC for wires having too large a radius for the thin-wire approximation [ 11]. From assumptions 1, 2 and 3, the surface current on a wire of radius can be replaced by a filamentary current I where

where

= distance parameter along the wire axis at

= unit vector tangent to the wire axis at

Equation (5) then becomes

(6)

where the integration is over the length of the wire. Enforcing the boundary condition in the axial direction reduces equation (6) to the scalar equation,

(7)

Since is now the point at on the wire axis while is a point at s on the wire surface and the integrand is bounded.

3. Numerical Solution

The integral equations (7) is solved numerically in NEC by a form of the method of moments. An excellent general introduction to the method of moments can be found in R.F. Harrington’s book, “Field Computation by Moment Methods” [ 13]. A brief outline of the method follows.

The method of moments applies to a general linear-operator equation,

(8)

where f is an unknown response, e is a known excitation, and L is a linear operator (an integral operator in the present case). The unknown function f may be expanded in a sum of basis functions, , as :

(9)

A set of equations for the coefficients are then obtained by taking the inner product of equation (9) with a set of weighting functions

(10)

Due to the linearity of L, substituting equation (9) for f yields,

(11)

This equation can be written in matrix notation as :

(12)

where

(13)

The solution is then

(14)

For the solution of equation (7), the inner product is defined as

(15)

where the integration is over the structure surface. Various choices are possible for the weighting functions and basis functions . When , the procedure is known as Galerkin’s method. In NEC the basis and weight functions are different, being chosen as a set of delta functions

(16)

with a set of points on the conducting surface. The result is a point sampling of the integral equations known as the collocation method of solution. Wires are divided into short straight segments with a sample point at the centre of each segment.

The choice of basis functions is very important for an efficient and accurate solution. In NEC, the support of is restricted to a localised subsection of the surface near . This choice simplifies the evaluation of the inner-product integral and ensures that the matrix will be well conditioned. For finite N, the sum of cannot exactly equal a general current distribution so the functions should be chosen as close as possible to the actual current distribution. The basis functions used in NEC are explained in the following sections.

3.1. Current expansion on wires

Wires in SNEC are modelled by short straight segments with the current on each segment represented by three terms - a constant, a sine, and a cosine. This expansion was first used by Yeh and Mei [ 14] and has been shown to provide rapid solution convergence [ 15],[ 16]. It has the added advantage that the fields of the sinusoidal currents are easily evaluated in closed form. The amplitude of the constant, sine, and cosine terms are related such that their sum satisfies physical conditions on the local behaviour of current and charge at the segment ends. This differs from AMP where the current was extrapolated to the centres of the adjacent segments, resulting in discontinuities in current and charge at the segment ends. Matching at the segment ends improves the solution accuracy, especially at the multiple-wire junctions of unequal length segments where AMP extrapolated to an average length segment, often with inaccurate results.

The total current on segment number j in NEC has the form :

(17)

where is the value of at the centre of segment j and is the length of segment j. Of the three unknown constants , and two are eliminated by local conditions on the current leaving one constant, related to the current amplitude, to be determined by the matrix equation. The local conditions are applied to the current and to the linear charge density, q, which is related to the current by the equation of continuity :

(18)

At a junction of two segments with uniform radius, the obvious boundary conditions are that the current and charge are continuous at the junction. At a junction of two or more segments with unequal radii, the continuity of current is generalised to Kirchoff’s current law that the sum of currents into the junction is zero. The total charge in the vicinity of the junction is assumed to distribute itself on individual wires according to the wire radii. T.T. Wu and R.W.P. King [ 17] have derived a condition for the linear charge density on a wire at a junction that neglects local coupling effects. Using their derivation, , can be determined by :

(19)

where

Q is related to the total charge in the vicinity of the junction and is constant for all wires at the junction.

At a free wire end, the current may be assumed to go to zero. On a wire of finite radius, however, the current can flow into the end cap and hence be nonzero at the wire end. In one study of this effect, a condition relating the current at the wire end to the current derivative was derived [ 18]. For a wire radius , this condition is :

(20)

where J_o and J₁ are Bessel functions of order 0 and 1 respectively. The unit vector is normal to the end cap. Hence, is +1 if the reference direction, , is toward the end, and -1 if is away from the end.

Thus, for each segment two equations are obtained from the two ends. At free ends the following holds :

(21)

When the end of a segment connects to other segments, a second equation is obtained :

(22)

Two additional unknowns and are associated with the junctions but can be eliminated by Kirchoff’s current equation at each junction. The boundary condition equations provide the additional equation-per-segment to completely determine the current function of equation (17) for every segment.

To apply these condition, the current is expanded as a sum of basis functions chosen so that they satisfy the local conditions on current and charge in any linear combination. A typical set of basis functions and their sum on a four segment wire are shown in Figure 1. For a general segment I in Figure 2 the i^th basis function has a peak on segment i and extends onto every segment connected to i, going to zero with zero derivative at the outer ends of the connected segments.

Figure 1: Current basis function and sim on a four segment wire.

Figure 2: Segments covered by the i^th basis function.

The general definition of the i^th basis function is given below. For the junction and end conditions described above, the following definitions apply for the factors in the segment end conditions:

(23)

and

(24)

The condition of zero current at a free end may be obtained by setting to zero.

The portion of the i^th basis function on segment i is then :

(25)

If and , end conditions are :

(26)

(27)

If and , end conditions are :

(28)

(29)

If and , end conditions are :

(30)

(31)

Over segments connected to end 1 of segment i, the i^th basis function is :

(32)

End conditions are :

(33)

(34)

(35)

Over segments connected to end 2 of segment i, the i^th basis function is :

(36)

End conditions are :

(37)

(38)

(39)

Equations (25), (32) and(36), defining the complete basis function, involve unknown constants. Of these, unknowns are eliminated by the end conditions in terms of and which can then be determined from the two Kirchoff’s current equations:

(40)

(41)

The complete basis function is then defined in terms of one unknown constant. In this case was set to -1 since the function amplitude is arbitrary, being determined by the boundary condition equations. The final result is given below:

(42)

(43)

(44)

(45)

(46)

(47)

For and :

(48)

(49)

(50)

(51)

(52)

For and :

(53)

(54)

(55)

(56)

For and :

(57)

(58)

(59)

(60)

For all cases :

(61)

(62)

where the sum of is over segments connected to end 1 of segment i, and the sum for is over segments connected to end 2. If the complete basis function is :

(63)

When a segment end is connected to a ground plane, the end condition on both the total current and the last basis function is :

(64)

This replaces the zero current condition at a free end. This condition does not require a separate treatment, however, but is obtained by computing the last basis function as if the last segment is connected to its image segment on the other side of the surface.

It should be noted that in AMP, the basis function has unit value at the centre of segment i and zero value at the centres of connected segments although it does extend onto the connected segments. As a result, the amplitude of is the total current at the centre of segment i must be computed by summing the contributions of all basis functions extending onto segment i.

3.2. Evaluation of the fields

The current on each wire segment has the form :

(65)

where

The solution requires the evaluation of the electric field at each segment due to this current. Three approximations of the integral equation are used: a thin-wire form for most cases, an extended thin-wire form for thick wires, and a current element approximation for large interaction distances. In each case, the evaluation of the field is greatly simplified by the use of formulas for the fields of the constant and sinusoidal current components.

The accuracy of the thin-wire approximation for a wire of radius a and length depends on and . Studies have shown that the thin-wire approximation leads to errors of less than 1% for greater than 8 [ 11]. Furthermore, in the numerical solution of the EFIE, the wire is divided into segments less than about in length to obtain an adequate representation of current distribution thus restricting to less than about 0.08. The extended thin-wire approximation is applicable to shorter and thicker segments, resulting in errors less than 1% for greater than 2.

For the thin-wire kernel, the source current is approximated by a filament on the segment axis while the observation point is on the surface of the observation segment. The fields are evaluated with the source segment on the axis of a local cylindrical coordinate system as illustrated in Figure 3.

Figure 3: Current-filament geometry for the thin-wire kernel.

Using

(66)

(67)

the and components of the electric field at P due to a sinusoidal current filament of arbitrary phase

(68)

are :

(69)

(70)

For a current that is constant over the length of the segment with strength I, the fields are :

(71)

(72)

These field expressions are exact for the specified currents. The integral over of is evaluated numerically in NEC.

Substituting sine and cosine currents and evaluating the derivatives yields the following equations for the fields. For

(73)

(74)

(75)

For a constant current of strength :

(76)

(77)

Despite the seemingly crude approximation, the thin-wire kernel does accurately represent the effect of wire radius for wires that are sufficiently thin. The accuracy range was studied by Poggio and Adams [ 11] where an extended thin-wire kernel was developed for wires that are too thick for the thin-wire approximation,

The derivation of the extended thin-wire kernel starts with the current on the surface of the source segment with surface density.

(78)

where is the radius of the source segment. The geometry for evaluation of the fields is shown in Figure 4. A current filament of strength is integrated over with

(79)

(80)

Thus, the z component of the field of the current tube is

(81)

Figure 4: Current geometry for the extended thin-wire kernel.

For the component of field, the change in the direction of must be considered. The field in the direction is :

(82)

where

(83)

The integral over in equations (81) and (82) cannot be evaluated in closed form. Poggio and Adams, however, have evaluated them as a series in powers of a² [ 11]. The first term in the series gives the thin-wire kernel. For the extended thin-wire kernel, the second term involving a² is retained with terms of order neglected. As with the thin-wire kernel, the field observation point is on the segment surface. Hence, when evaluating the field on the source segment, .

The field equations with the extended thin-wire approximation are given below. For a sinusoidal current of equation (68) :

(84)

(85)

For a constant current of strength :

(86)

(87)

The term is the series approximation of :

(88)

where

(89)

Neglecting terms of order :

(90)

(91)

(92)

The term G₂ is the series approximation of :

(93)

To order :

(94)

(95)

Equation (86) makes use of the relation :

(96)

while equation (87) follows from :

(97)

When the observation point is within the wire , a series expansion in rather than is used for and. For this simply involves inter-changing and in equations (90) and (91). Then for , with :

(98)

(99)

the expressions for G₁, G₂ and their derivatives are :

(100)

(101)

(102)

(103)

(104)

Special treatment of bends in wires is required when the extended thin-wire kernel is used. The problem stems from the cancellation of terms evaluated at and in the field when segments are part of a continuous wire. The current expansion in NEC results in a current having a continuous value and derivative along a wire without junctions. This ensures that for two adjacent segments on a straight wire, the contributions to the component of electric field at of the first segment exactly cancel the contributions from , representing the same point, for the second segment. For a straight wire of several segments, the only contributions to with either the thin-wire or extended thin-wire kernel come from the two wire ends and the integral of along the wire. For the component of field or either component at a bend, while there is not complete cancellation, there may be partial cancellation of large end contributions.

The cancellation of end terms makes necessary a consistent treatment of the current on both sides of a bend for accurate evaluation of the field. This is easily accomplished with the thin-wire kernel since the current filament on the wire axis is physically continuous around a bend. However, the current tube assumed for the extended thin-wire kernel cannot be continuous around its complete circumference at a bend. This was found to reduce the solution accuracy when the extended thin-wire kernel was used for bent wires.

To avoid this problem in NEC, the thin-wire form of the end terms in equation (69) through (72) is always used at a bend or change in radius. The extended thin-wire kernel is used only at segment ends where two parallel segments join, or at free ends. The switch from extended thin-wire form to the thin-wire form is made from one end of a segment to the other rather than between segments where the cancellation of terms is critical.

When segments are separated by a large distance, the interaction may be computed with sufficient accuracy by treating the segment current as an infinitesimal current element at the segment centre. In spherical co-ordinates, with the segment at the origin along the axis, the electric field is :

(105)

(106)

The dipole moment M for a constant current I on a segment of length is :

(107)

For a current with

(108)

while for a current

(109)

Use of this approximation saves a significant amount of time in evaluating the interaction matrix elements for large structures. The minimum interaction distance at which it is used is selected by the user in NEC. A default distance of one wavelength is set, however.

For each of the three methods of computing the field at a segment due to the current on another segment, the field is evaluated on the surface of the observation segment. Rather than choosing a fixed point on the segment surface, the field is evaluated at the cylindrical coordinates with the source segment at the origin. If the centre point on the axis of the observation segment is at , then :

(110)

where is the radius of the observation segment. Also, the component of tangent to the observation segment is computed as :

(111)

Inclusion of the factor , which is the cosine of the angle between and , is necessary for accurate results at bends in thick wires.

3.3. The matrix equation for current

For a structure having wire segments, the order of the matrix in equation (12) is The matrix equation has the form :

(112)

I is then the column vector of segment basis function amplitudes. The elements of E are the left-hand side of equation (7) evaluated at the segment centres.

A matrix element in matrix represents the electric field at the centre of segment i due to the jth segment basis function, centred on segment j.

3.4. Solution of the matrix equation

The matrix equation

(113)

is solved in NEC by Gauss elimination [ 19]. The basic step is factorisation of the matrix into the product of an upper triangular matrix and a lower triangle matrix where

(114)

The matrix equation is then

(115)

from which the solution, I, is computed in two steps as

(116)

and

(117)

Equation (116) is first solved for F by forward substitution, and equation (117) is then solved for I by backward substitution.

The major computational effort is factoring G into L and U. This takes approximately multiplication steps for a matrix of order N compared to N³ for inversion of G by the Gauss-Jordan method. Solution of equations (116) and (117), making use of the triangular properties of L and U, takes approximately as many multiplications as would be required for multiplication of by the column vector E. The factored matrices L and U are saved in NEC since the solution for induced current may be repeated for a number of different excitations. This, then, requires only the repeated solution of equations (116) and (117).

Computation of the elements of the matrix G and solution of the matrix equation are the two most time-consuming steps in computing the response of a structure, often accounting for over 90% of the computation time. This may be reduced substantially by making use of symmetries of the structure, either symmetry about a plane, or symmetry under rotation.

In rotational symmetry, a structure having M sectors is unchanged when rotated by any multiple of 360/M degrees. If the equations for all segments in the first sector are numbered first and followed by successive sectors in the same order, the matrix equation can be expanded in submatrices in the form :

(118)

If there are equations in each sector, and are element column vectors of the excitations and currents in sector i. is a submatrix of order containing the interaction fields in sector 1 due to currents in sector i. Due to symmetry, this is the same as the fields in sector k due to currents in sector , resulting in the repetition pattern shown. Thus, only matrix elements in the first row of submatrices need be computed, reducing the time to fill the matrix by a factor of 1/M.

The time to solve the matrix equation can also be reduced by expanding the excitation sub-vectors in a discrete Fourier series as :

(119)

(120)

where

(121)

where

indicates the conjugate of the complex number.

Examining a component in the expansion,

(122)

it is seen that the excitation differs from sector to sector only by a uniform phase shift. This excitation of a rotationally symmetric structure results in a solution having the same form as the excitation, i.e.,

(123)

It can be shown that this relation between solution and excitation holds for any matrix having the form of that in equation (118). Substituting these components of E and I into equation (118) yields the following matrix equation of order relating to :

(124)

The solution for the total excitation is then obtained by an inverse transformation,

(125)

The solution procedure, then, is first to compute the M submatrices and Fourier-transform these to obtain

(126)

The matrices of order are then each factored into upper and lower triangular matrices by the Gauss elimination method. For each excitation vector, the transformed sub-vectors are then computed by equation (120) and the transformed current sub-vectors are obtained by solving the M equations,

(127)

The total solution is then given by equation (125).

The same procedure can be used for structures that have planes of symmetry. The Fourier transform is then replaced by even and odd excitations about each symmetry plane. All equations remain the same with the exception that the matrix with elements , given by equation (121), is replaced by the following matrices:

For one plane of symmetry :

(128)

For two orthogonal planes of symmetry :

(129)

and for three orthogonal symmetry planes :

(130)

For either rotational or plane symmetry, the procedure requires factoring of M matrices of order rather than one matrix of order . Each excitation then requires the solution of the M matrix equations. Since the time for factoring is approximately proportional to the cube of the matrix order and the time for solution is proportional to the square of the order, the symmetry results in a reduction of factor time by and in solution time by . The time to compute the transforms is generally small compared to the time for matrix operations since it is proportional to a lower power of . Symmetry also reduces the number of locations required for matrix storage by since only the first row of submatrices need be stored. The transformed matrices, can replace the matrices as they are computed.

NEC includes a provision to generate and factor an interaction matrix and save the result of a file. A later run, using the file, may add to the structure and solve the complete model without unnecessary repetition of calculations. This procedure is called the Numerical Green’s Function (NGF) option since the effect is as if the free space Green’s function in NEC were replaced by the Green’s function for the structure on the file. The NGF is particularly useful for a large structure, such as a ship, on which various antennas will be added or modified. It also permits taking advantage of partial symmetry since a NGF file may be written for the symmetric part of a structure, taking advantage of the symmetry to reduce computation time. Unsymmetric parts can be added in a later run.

For the NGF solution the matrix is partitioned as

(131)

where is the interaction matrix for the initial structure, is the matrix for the added structure, and and represent mutual interactions. The current is computed as

(132)

(133)

after the factored matrix has been read from the NGF file along with other necessary data.

Electrical connections between the new structure and the old (NGF) structure require special treatment. If a new wire connects to an old wire the current basis function for the old wire segment is changed by the modified condition at the junction. The old basis function is given zero amplitude by adding a new equation having all zeros except for a one in the column of the old basis function. A new column is added for the corrected basis function.

4. Effect of a Ground Plane

In the integral equation formulation used in NEC, a ground plane changes the solution in three ways:

(1) by modifying the current distribution through near-field interaction

(2) by changing the field illuminating the structure

(3) by changing the re-radiated field

Effects (2) and (3) are easily analysed by plane-wave reflection as a direct ray and a ray reflected from the ground. The re-radiated field is not a plane wave when it reflects from the ground, but, as can be seen from reciprocity, plane-wave reflection gives the correct far-zone field. Analysis of the near-field interaction effect is, however, much more difficult.

In Section 2.1, the kernels of the integral equations are free-space Green’s functions, representing the E or H field at a point due to an infinitesimal electric current element at . When a ground is present the free space Green’s function must be replaced by Green’s functions for the ground problem. The solution for the fields of current elements in the presence of a ground plane was developed by Arnold Sommerfeld [ 20]. While this solution has been used directly in integral-equation computer codes, excessive computation time greatly limits its use. Numerous approximations to the Sommerfeld solution have been developed that require less time for evaluation but all have limited applicability.

The NEC code has three options for grounds. The most accurate for lossy grounds uses the Sommerfeld solution for interaction distances less than one wavelength and an asymptotic expansion for larger distances. To keep the solution time reasonable, a grid of values of the Sommerfeld solution is generated and interpolation is used to find specific values. The solution for a perfectly conducting ground is much simpler since the ground may be replaced by the image of the currents above it. The third option models a lossy ground by a modified image method using the Fresnel plane-wave reflection coefficients. While specular reflection does not accurately describe the behaviour of near fields, the approximation has been found to provide useful results for structures that are not too near to the ground [ 21] and [ 22]. The attraction of this method is its simplicity and speed of computation which are the same as for the image method for perfect ground.

4.1. The Sommerfeld/Norton method

The Sommerfeld/Norton ground option in NEC originated with the code WFLLL2A [ 23] which uses numerical evaluation of the Sommerfeld integrals for ground fields when the interaction distance is small and uses Norton’s asymptotic approximations [ 24] for larger distances. Since evaluation of the Sommerfeld integrals is very time consuming, a table of Sommerfeld integrals is generated and a bivariate interpolation scheme is used to obtain the field value for integration over the current distribution. This method greatly reduces the required computation time. The method used in NEC to evaluate the field over ground is described below, and the numerical evaluation of the Sommerfeld integrals to fill the interpolation grid is discussed in Section 4.2.

The electric field above an air-ground interface due to an infinitesimal current element of strength also above the interface, with parameters shown in Figure 5, is given by the following expressions:

(134)

(135)

(136)

(137)

(138)

(139)

(140)

(141)

where the superscript indicates a vertical (V) or horizontal (H) current element and the subscript indicates the cylindrical component of the field vector. The horizontal current element is along the x-axis.

Figure 5: Coordinates for evaluating the field of a current element over ground.

G₂₂ and G₂₁ are the free space and image Green’s functions :

(142)

(143)

where

(144)

(145)

and U₂₂ and V₂₂ are Sommerfeld integrals involving the zero order Bessel function,

(146)

(147)

where

(148)

(149)

In NEC we need to compute the fields due to current filaments with arbitrary length and orientation by combining the field components in equations (134) through (138) and integrating over current distributions composed of constant, sine, and cosine components. Direct numerical integration over the segments is difficult due to singularities in the fields. has a singularity while , and each have singularities. The derivatives in the field expressions result in singularities with a triplet-like behaviour in the field components parallel to the current filament. The resulting cancellation makes accurate numerical integration near the singularity very difficult.

A free-space field has a similar singularity, but as discussed in Section 3.2, the integral over a straight filament may be evaluated in closed form for a sinusoidal current with free-space wavelength and involves only a numerical integration of for a constant current. The dominant singular component of the ground field may be integrated in the same way. The terms involving in equation (134) through (138) are, in fact, the field of the current element in free space, and their integral is obtained from the free-space routines in NEC.

The remaining terms represent the field due to ground and are singular at . The singularities in and result from the failure of the integrals in equations (130) and (131) to converge without the exponential and Bessel functions as r and go to zero. The singular behaviour of and as r and go to zero may be found by setting since the dominant contributions to the integrals for small and come from much greater than and . Here, however, we only replace by and use the integrals.

(150)

(151)

where

(152)

which have the correct singular behaviour and can be combined with the terms. The field components due to ground (equation (134) through (138) without the terms) may then be written as :

(153)

(154)

(155)

(156)

(157)

where

(158)

(159)

In equation (153) through (157) the dominant singular component has been subtracted out of and combined with . The integral for converge without the exponential or Bessel function factors and remains finite as and go to zero. The derivatives of in the field expressions have singularities, but this is much less of a problem for numerical integration than the previous singularity. The singularity could be taken out of also, but instead, a term is taken out that results in the final terms in equations (153) through (157) being the image field multiplied by . The integral over the current filament of these image terms is evaluated by the free-space equations leaving only the and terms to be integrated numerically. still has a singularity, but that is no worse than the derivatives of . With the thin-wire approximation, is never less than the wire radius so the integration is not difficult in practical cases.

The components left for numerical integration over the current distribution are then,

(160)

(161)

(162)

(163)

(164)

Since the integrals in equations (158) and (159) cannot be evaluated in closed form, the following terms must be evaluated by numerical integration over :

(165)

(166)

(167)

(168)

(169)

(170)

Where

(171)

(172)

Evaluating these integrals over for each point needed in the numerical integration over the current distribution is slow in even the fastest computers. Hence, an interpolation technique is used for the remaining field components as was done in the code SOMINT for the total field due to ground. Since the integrals depend only on and a grid of values of generated for the field components of equations (160) through (163) and bivariate interpolation is used to obtain values for integration over a current distribution.

To facilitate interpolation in the region of the singularity, the components are dived by a function having a similar singularity and interpolation is performed on the ratio. The field components of equations (160) through (163) are divided by for all values of to remove the singularity and the free-space phase before interpolation. The factors or are also omitted until after interpolation to avoid introducing the dependence. The surfaces to which interpolation is applied are then :

(173)

(174)

(175)

(176)

After interpolation on the smoothed surfaces the results are multiplied by the omitted factors to give the correct values.

With the singularity removed, interpolation may be used for arbitrarily small values of and . The values for in the interpolation grid must be found as limits for approaching zero, however, since the integrals do not converge in this case. When and approach zero the dominant contributions in equations (165) through (170) come from large . Hence, the singular behaviour can be found by setting and equal to . First, however, it is necessary to approximate D₁ and D₂ for as :