Fast Far Field Computation of Single and Dual Reflector Antennas

The physical optics (PO) method has been widely used for the analysis of the electromagnetic behavior of single and dual reflector antennas. An extensive work has been done by the authors of this paper in order to increase the speed for obtaining far field patterns from single and dual geometries and also in order to increase the accuracy of the method. This paper reviews these contributions and improves the existing published work with the physical interpretation of the radiation from a single patch and the computer implications when using acceleration techniques such as OpenMP.


Introduction
In engineering electromagnetic reflector antennas are commonly used for long-range radio communications, in applications such as satellite communications, radiolinks, and space exploration.In these applications it is necessary to accurately estimate the radiation pattern of the antenna to fulfill requirements such as gain, sidelobe levels, and crosspolar radiation.
The electromagnetic scattering by metallic surfaces has been extensively treated in the literature with special emphasis on the reflector antenna application case.Usually, the radiated integral of the current on the metallic surface is numerically solved with the aid of some mathematical and physical approximations.
The earliest method was based on series representation [1] in order to develop the radiated field.This mathematical approximation was then adapted to the case of an offset paraboloid using Jacobi polynomial series method [2].
The second type of mathematical approximation, based on numerical integration, was presented in [3], where it is demonstrated that the convergence depends on the choice of the integration grid coordinate system.In [4] popular numerical integration methods applied to the reflector antenna problem are compared.It was shown that the most accurate one to develop the surface integral consisted of using Gauss-Zernike polynomial integration for the radial coordinate and the trapezoidal rule for the angular coordinate along the circumference.
Other very popular physical approximations are based on geometrical optics with aperture integration (GO + AI) and physical optics (PO).These techniques are compared in [5] for offset parabolic antenna.An alternative consists of combining GO with the geometrical theory of diffraction (GTD) to represent the reflector rim.In [6] it is demonstrated that the PO solution agrees well with the GO + GTD for the copolar component.
In [7] the integration across the main reflector of a Cassegrain antenna is transferred to the aperture of the feed.When using several reflector surfaces, the interaction among surfaces must be taken into account.For example, in [8] a ray tracing technique is described in order to model the propagated field through a multireflector system that is numerically specified.
In [9] the authors presented a comparison between using direct numerical integration and using triangular or rectangular patches.In that book chapter, the methods that used patches demonstrated to be faster than those using numerical integration.In Section 2, a brief summary of this book chapter is included.This method can be described as an approximation of the magnetic field integral equation (MFIE).This approximation consists of neglecting the interactions between the patches on the same surface.This is a good approximation when the surfaces are flat enough as in reflector antennas.
In Section 3, a new physical interpretation of the triangular patch current integration is provided.
In Section 4, the case of dual reflector geometries is addressed.Computing the PO currents induced by the feed on the subreflector and those induced by the subreflector on the main reflector surface, in addition to the feed fields, provides an estimation of the radiated fields including the spillover effect.The blockage introduced by the subreflector can be simulated by a second iteration consisting of including the radiation of a new set of blocking currents across the subreflector.The PO blocking currents are derived from the fields produced by the main reflector into the subreflector.Using the same approximation of the MFIE as with single reflector antennas and assuming that the incident field only exists on the subreflector, the process can be reduced to invert a small matrix, which is possible to perform in modern computers.This method was introduced in [10] and improved in [11][12][13].
In Section 5 some techniques to improve the speed of the computation are presented.The first [14] is based on precomputing the complex exponentials and on using OpenMP [15].And the second [16] is based on using memory hierarchy techniques.
Finally, in Section 6 some comparative results are shown.Two computer programs were developed for educational purposes as described in [17,18].The final version is downloadable from the web for free since 2006 (http://www.com.uvigo.es).

Physical Optics for Reflector Antennas
In this section a method based on magnetic field integral equation (MFIE) [19,20] will be used to estimate the fields radiated by a reflector antenna.
MFIE model is obtained from equivalent currents on the surface [21], by a similar way as classical high frequency methods such as physical optics (PO) but is more accurate than PO because it calculates the coupling among surfaces.In this paper, MFIE model is modified in order to adapt it to the particular problem with smooth (almost flat) and electrically big surfaces [12].

MFIE Field Formulation.
The current (  ) on a PEC with smooth and closed surface can be obtained using where   is the observation point,    is the source point,   (  ) is the incident magnetic field at point   , symbol − ∫ means the principal value of the integral, and (  ,    ) is free space Green function, defined as Taking into account (2), (1) can be rewritten as where   =   −    and   = |  −    |.After obtaining equivalent currents, the radiated magnetic field at points  not belonging to the surface can be obtained by where  =  −    ,  = | −    | and R = /.

PO Field Formulation.
Physical optics (PO) is a simplification of (3) when the surface is flat enough to consider that the currents (   ) are on the same plane as the vector position   pointing to another point on the surface.In this case, their cross-product (   )×  is colinear with the normal vector n(  ) and their cross-product tends to zero.With this approximation, (3) is simplified as and the radiated magnetic field at points not belonging to the surface can be obtained with physical optics by considering When the observation point is at a distance  = || ≫ |   | for every source point, then the distance can be simplified as follows: where the term r ⋅    is important only for the phase calculation but it can be neglected in the amplitude terms.With this approximation, (6) for points far from the surface can be obtained as If the distance  = || ≫ , then (8) can be simplified by In the next subsections, ( 8) is used to obtain the radiated fields from a small patch at any point on other surfaces.

Patch Discretization.
The radiated field at an observation point   can be obtained from ( 8) by dividing the surface in subdomains (patches) small enough that we can consider that the point   is far from each patch.The magnetic field   at the point   can be represented as where   is the magnetic field radiated by the patch  at   .The field radiated by each patch will be where   represents the surface of the considered patch .PO approximation can be used when the patch surface is flat enough to consider that the whole patch has the same  normal vector n .Inserting this vector into (11) produces the following expression: 2.4.Physical Approximation.Equation ( 11) is difficult to compute due to the existence of the incident field inside the integral.But if the patch is flat and small, the incident field on this patch can be considered locally as a plane wave.With this approximation, the amplitude and the direction of the field are assumed to be constant on the patch, and the phase will vary according to the direction of the propagation of the incident plane wave p .Consider that Then, Including ( 13) and ( 14) into (12), where Then, the scattered far field will now be obtained as 2.5.Triangular Patches.Let us consider a triangular flat patch as seen in Figure 1.The triangle is defined by three points  1 ,  2 , and  3 . 0 is the point where  0 is calculated which is assumed to be at the barycenter  0 = ( 1 +  2 +  3 )/3.If we define V  =   −  , the local coordinates of the triangle can be described with two variables  and V as The normal vector and the area of the triangle are related to the cross-product of V 12 and V 13 as follows: Using this coordinate system, ( 16) is then given by whose solution is where Equation ( 21) has the following singular values: Form small values of  and , a second order series approximation can be used: Finally, for directions close to the main lobe or for very small patches, the expression (25) could be approximated as (26)

Physical Interpretation
In this section the relationship between this integral representation and other methods such as geometrical optics is explained.
3.1.Vertex Contributions.Equation ( 21) can be reformulated in order to clarify its physical meaning: Each of the three terms of ( 27) can be interpreted as the contribution of a vertex to the radiation integral   .Taking into account the definitions of V 0 it can be stated that So, expression (27) can be transformed into where each phase term  V 0 ⋅(r  −p  ) is due to the displacement from the barycenter to the vertex  and the amplitude terms depend on the vectors from the vertex  to the other two vertices, according to the values of  and  given by ( 22) and (23): 3.2.Reflection Point.The theory of geometrical optics states that the contributors for the radiated field are the reflection and the diffraction points.The diffraction comes from the border of the object, and if the surface is divided in patches, then each patch creates also diffraction fields that can be compensated with the diffraction generated by the surrounding patches.But ( 27) and (29) do not represent separately diffraction and reflection terms, as it has been done in [22].However, it is possible to study the reflection effect, that appears in the direction r = p or r = p − 2n ⋅ (p  ⋅ n).
If the triangle is big enough compared with the wavelength, then the main effect is the reflection, as it can be seen in Figure 2 where the incident field forms 30 degrees with the plane of the triangle.It can be observed that when the size of the triangle increases, then the radiated field tends to the geometrical optics effect, represented by a narrow beam at 30 ∘ .

Discretization Effect as a Function of the Size of the Object.
Methods like MoM require discretization less than /10 in order to converge to the real solution.With PO the size of the triangles depends on the required accuracy in the angular domain, as it can be seen in the next example.
Consider a circle of radius  on the plane  with a plane wave normal to the circle as incident field.The scattered field for  = 0 will be with a maximum value   ()| max =   (0) =  2 .Adapting (27) to this example, where the phase reference point will be at the center of the circle and V 12 =  ⋅ ρ and V 13 =  ⋅ ρ+Δ , this leads to   =  cos   sin  and   =  cos(  + Δ) sin .Then, the integral will be where Δ = 2/.The limit of this expression tends to the integral: Now, we are going to determine the number of angular sectors   that are necessary to obtain an accurate solution when the radius  varies.In Figures 3, 4, and 5 the convergence can be seen when the radius is 2, 10, and 50.
The convergence behaviour shown in these figures is similar because in each case the radiation pattern is calculated in a margin according to the size of the circle, where the radiation is in a margin of 40 dB less than the maximum.Greater areas imply lower beam widths, and the number of triangles required for the same accuracy is almost the same.
But if the same accuracy is required in the whole  angular view, then the number of triangles must be increased in the same rate as the radius is increased.

PO for Dual Reflector Antennas
When there is more than one surface, it is not possible to eliminate the radiation of a patch over the patches on the other surfaces.Applying this situation to a dual reflector antenna, the radiation between both surfaces must be calculated.

MFIE Approximation for Dual Reflector Antennas.
If the patches are small enough to use (26), then the amplitude and also the phase at each patch can be considered constant.So the current at patch  will be a vector   .This simplifies the integral assuming that   ≈ (  −   ) =   and |  | =   , so (3) can be obtained as where   is the area of th-triangle.
The equivalent current in (34) is due to incident (PO) currents and to the coupling among surfaces.We can define and now the equivalent current can be expressed as where  and  are for points on different surfaces.
When the currents on both surfaces are obtained, the scattered far fields can be calculated at a direction r using where  is intended for all the main reflector and subreflector patches and 4.2.Iterative Solution for MFIE Currents.The solution to (37) can be obtained with an iterative method or with a matrix method.From this point,  will be for reflector patches and  for main reflector patches.
With the iterative method, the currents can be obtained using an iterative process in which the first step is using PO currents: where This iterative process converges in few steps, usually less than six.

Matrix Formulation.
In order to obtain a matrix formulation, we need to define the local coordinates on the patches.A flat patch can be defined with two unitary and orthogonal vectors û1 and û2 for th-patch.Using triangular patches, these vectors can be obtained as The currents defined in (37) can be expressed in terms of û1 and û2 as and now (37) can be divided into two equations: can also be divided into two components.After some manipulations, (43) can be expressed as The mutual coupling between the surfaces can be expressed as This equation system can be expressed with matrix formulation: where In (47), the first two subscripts of  are for rows and the third and fourth subscripts are for the columns.With dual reflector antennas, we can assume that the incident field on main reflector is negligible.Then, in order to solve (46) the only matrix to be inverted is the one with subreflector currents: where and the components of the main reflector currents are calculated directly with (50).
Finally, all the currents must be transformed to the general coordinate system, using (42).

Computer Acceleration
The methods described here for obtaining the radiation pattern of single and dual reflector antennas are really fast, even for electrically big surfaces.The results can be observed in the next section.But there are some new computer techniques that can reduce even more the computer time.The first is related to the size of the memory.The second consists of precalculating the exponential functions.And the third is related to working with several computer cores at the same time.

Memory Size.
With single reflector antennas, almost the whole computer time is spent for obtaining the integral (from ( 21) to ( 23)), mainly with the three exponential terms.These exponentials must be computed for all the triangles and all the directions, so it can be implemented with a matrix expression or in a loop expression.With matrices, the programs can use acceleration techniques such as BLAS [23] but require more size in the computer memory.
Using loops offers two possibilities.The first is to obtain for each direction the contribution of the scattering of all the patches and then adding their contributions.The second one consists of obtaining for each patch the scattering in all the directions and then adding these vectors to obtain the radiation in each direction.
In real examples, the number of triangles is usually larger than the number of directions, so the first loop solution needs more memory than the second loop solutions.But in any case, memory requirements are much lower than the matrix implementation.

Complex Exponential.
Most of the execution time is used in the calculation of complex exponentials.
For single reflector antennas, the exponential term  r⋅  and those that appear in the calculation of the integral ( 21) must be computed for each direction and each triangle.In order to reduce the computation time, the exponentials can be precomputed previously using the property that the complex exponential is periodical with period 2.Using the function Int to extract the integer part of the real numbers, the exponential must be computed in three steps: = cos pre (arg) + sin pre (arg) .
When using a precomputed table of 3000 values, the execution time is reduced by a factor of 2.5 appearing as a slight quantization error.
For dual reflector antennas, the same problem appears for calculating the scattering fields.But obtaining   in (36) spends much more time because this exponential must be computed for each patch on the subreflector radiating at each patch on the main reflector.The proposed method to reduce the computer time of the factor (1 − /) ⋅  − / 2 has the following steps: where  min is the minimum value of  and  max the maximum value.The term  2 is calculated directly from the components  2  +  2  +  2  , so we avoid to use the square root function.When using a precomputed table of 12000 values, the execution time is reduced by a factor of 3.4 appearing as a slight quantization error.

OpenMP.
In [13] OpenMP has been used to parallelize the code and memory-hierarchy-based optimization techniques to reduce the computer time of the code.Using these techniques, the computer time can be reduced in a factor very close to the number of cores of the CPU (usually 4 or 8).
The main problem that appears using OpenMP is the order of the loops.If the directions of the observation points are used as outer loop, then each core can compute the scattered field created by all the triangles in one direction, and at the end, it should store the result in a position of the output vector.But if the index of the triangles is used as the outer loop, then each core must compute the radiation in all directions and then use the reduction method to add all the results.Unfortunately, the reduction method is not well implemented for vectors in OpenMP, and each core must wait for the others to write their results.

Results
In this section, some results are shown to present the last improvements in the physical optics technique.

Single Reflector Antenna.
As an example, an antenna of 0.5 m of radius, 0.25 m of offset height, and 0.5 m of focal length was studied.The frequency is 100 GHz and it has a cos- type feed with taper 12 dB.The scattered field is obtained for  = /2 and  from 0 to 1 ∘ or 5 ∘ .
The convergence result is obtained raising the number of sectors   until the difference will be negligible.Then, the difference in directivity is accounted for in natural units, using The integral is approximated by the area as stated in (26), because with this small angular margin, it provides less error than the other approximations, and it is the fastest (about 30% when using (25)).These two differences can be observed in Figure 6 for the area and in Figure 7 for the approximation (25) of the integral.This difference with the convergence is measured in sectors with width 0.1 ∘ .With   = 42 (area of the triangles of 8 2 ), the difference is below 1%.With   = 24 (area of the triangles of 25 2 ), the difference is below 5%.And with   = 13 (area of the triangles of 86 2 ), the difference is below 25%.In Figure 8 the first case is represented, and Figure 9 represents the third case.
The fastest version of the antenna analysis has been obtained using C language and OpenMP directives, for using all the cores of the computer, and (53).With these improvements, the case with   = 42 is obtained in 0.023 seconds in a quad-core 2.6 GHz computer.

Dual Reflector Antenna.
The antenna used for this example was a Cassegrain reflector at 10 GHz.The main reflector had a radius of 0.5 m and the subreflector had a radius of 0.1162 m.The number of rings used was 25 for the main reflector and 10 for the subreflector.In order to maximize the coupling between the surfaces, the antenna had no offset height.
The results show that the iterative method converges in 6 steps, so 12 near field calculations were required.Using the matrix formulation (48), the Matlab code spent 4.61 seconds, because of the matrix inversion.Using (37) 12 times, the result was obtained in 3.73 seconds with Matlab.If the matrix   is first computed and then used 12 times, in Matlab the spent time was reduced to 1.27 seconds.Finally, if the code is developed in C with OpenMP and the technique (57), the spent time was reduced up to 0.18 seconds.

Conclusion
A review of previous achievements regarding the analysis of single and dual reflector antennas by physical optics is presented.It includes the PO method as a simplification of the more general MFIE formulation.The discretization of the surfaces in triangles is addressed including a new physical interpretation of the meaning of the triangle integral.The method has been applied to single and dual reflector antennas.In the dual reflector case a new contribution has been reported to improve the method accuracy by using mutual coupling between reflectors from the MoM formulation.
Finally, some improvements in the calculation speed of the radiation pattern have been presented.The use of OpenMP has been shown to be of special interest to reduce the computation time.

Figure 6 :
Figure 6: Difference using (26): at (a) the difference and at (b) the pattern.

Figure 7 :
Figure 7: Difference using (25): at (a) the difference and at (b) the pattern.

Figure 8 :Figure 9 :
Figure 8: Difference when   = 42: at (a) the difference and at (b) the pattern.