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.
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 cross-polar 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  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 .
The second type of mathematical approximation, based on numerical integration, was presented in , where it is demonstrated that the convergence depends on the choice of the integration grid coordinate system. In  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  for offset parabolic antenna. An alternative consists of combining GO with the geometrical theory of diffraction (GTD) to represent the reflector rim. In  it is demonstrated that the PO solution agrees well with the GO + GTD for the copolar component.
In  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  a ray tracing technique is described in order to model the propagated field through a multireflector system that is numerically specified.
In  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  and improved in [11–13].
In Section 5 some techniques to improve the speed of the computation are presented. The first  is based on precomputing the complex exponentials and on using OpenMP . And the second  is based on using memory hierarchy techniques.
Finally, in Section 6 some comparative results are shown.
2. Physical Optics for Reflector Antennas
MFIE model is obtained from equivalent currents on the surface , 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 .
2.1. 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
After obtaining equivalent currents, the radiated magnetic field at points not belonging to the surface can be obtained by where , and .
2.2. 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 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 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.
2.3. 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 . 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 . Consider that Then,
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 , , and . is the point where is calculated which is assumed to be at the barycenter . If we define , the local coordinates of the triangle can be described with two variables and as
The normal vector and the area of the triangle are related to the cross-product of and 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
3. 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 it can be stated that
So, expression (27) can be transformed into where each phase term 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 . However, it is possible to study the reflection effect, that appears in the direction or .
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 .
3.3. Discretization Effect as a Function of the Size of the Object
Methods like MoM require discretization less than 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 will be with a maximum value .
Adapting (27) to this example, where the phase reference point will be at the center of the circle and and , this leads to and . Then, the integral will be where .
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 , , and .
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.
4. 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.
4.1. 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 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.
4.3. 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 and for th-patch. Using triangular patches, these vectors can be obtained as
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).
5. 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.
5.1. 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  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.
5.2. Complex Exponential
Most of the execution time is used in the calculation of complex exponentials.
For single reflector antennas, the exponential term 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 . Using the function Int to extract the integer part of the real numbers, the exponential must be computed in three steps:
When using a precomputed table of 3000 values, the execution time is reduced by a factor of 2.5 appearing asa 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 has the following steps: where is the minimum value of and the maximum value. The term is calculated directly from the components , 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.
In  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.
In this section, some results are shown to present the last improvements in the physical optics technique.
6.1. 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 and from 0 to or .
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 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 . With (area of the triangles of ), the difference is below . With (area of the triangles of ), the difference is below . And with (area of the triangles of ), the difference is below . 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 is obtained in 0.023 seconds in a quad-core 2.6 GHz computer.
6.2. 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.
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.
This work was supported by the Spanish Government Grants TEC2011-28683-C02-02 and CONSOLIDER-INGENIO CSD2008-00068.
V. Galindo-Israel and R. Mittra, “A new series representation for the radiation integral with application to reflector antennas,” IEEE Transactions on Antennas and Propagation, vol. 25, no. 5, pp. 631–641, 1977.View at: Google Scholar
R. Mittra, Y. Rahmat-Samii, V. Galindo-Israel, and R. Norman, “An efficient technique for the computation of vector secondary pattems of offset reflectors,” IEEE Transactions on Antennas and Propagation, vol. 27, no. 3, pp. 294–304, 1979.View at: Google Scholar
J. R. Parkinson and M. J. Mehler, “Convergence of physical optics integral by Ludwig's technique,” Electronics Letters, vol. 22, no. 22, pp. 1161–1162, 1986.View at: Google Scholar
Y. Rahmat-Samii, “A comparison between GO/aperture-field and physical optics methods for offset reflector,” IEEE Transactions on Antennas and Propagation, vol. 32, no. 3, pp. 301–306, 1984.View at: Google Scholar
N. C. Albertsen and K. Pontoppidan, “Analysis of subreflectors for dual reflector antennas,” IEE Proceedings H, vol. 131, no. 3, pp. 205–213, 1984.View at: Google Scholar
A. M. Arias-Acuña, J. O. Rubiños-López, I. Cuiñas-Gómez, and A. García-Pino, “Electromagnetic scattering of reflector antennas by fast physical optics algorithms,” in Recent Res. Devel. Magnetics, vol. 1, pp. 43–63, 2000.View at: Google Scholar
F. Obelleiro, J. M. Taboada, J. L. Rodríguez, J. O. Rubiños, and A. M. Arias, “Hybrid moment-method physical-optics formulation for modeling the electromagnetic behavior of on-board antennas,” Microwave and Optical Technology Letters, vol. 27, no. 2, pp. 88–93, 2000.View at: Google Scholar
M. Graña, M. Arias, O. Rubiños, A. García-Pino, and J. A. Martínez-Lorenzo, “Iterative physical optics solution for MFIE on dual reflector antennas,” in Electromagnetic Theory Symposium (EMTS '07), pp. 1–3, Ottawa, Canada, 2007.View at: Google Scholar
M. Graña, M. Arias, O. Rubiños, and A. García-Pino, “Asymptotic MFIE for optimizing dual reflectors pattern,” in Proceedings of the 2nd European Conference on Antennas and Propagation (EuCAP '07), pp. 1–5, Edinburgh, UK, 2007.View at: Google Scholar
M. Graña-Varela, M. Arias, O. Rubiños, and A. García-Pino, “Rapid dual reflector shaping using ant colony optimization, fast iterated PO and asymptotic MFIE,” in Proceedings of the 3rd European Conference on Antennas and Propagation (EuCAP '09), pp. 2731–2735, Berlin, Germany, March 2009.View at: Google Scholar
M. Graña-Varela, M. Arias-Acuña, O. Rubiños, and A. García-Pino, “Shaped solid contour beam reflector antennas that adjust to varying operating conditions with just a few mechanical actuators,” IEEE Antennas and Wireless Propagation Letters, vol. 8, pp. 839–842, 2009.View at: Publisher Site | Google Scholar
H. Gómez-Sousa, J. A. Martínez-Lorenzo, O. Rubiños-López et al., “Strategies for improving the use of the memory hierarchy in an implementation of the modified current approximation (MECA) method,” ACES, vol. 25, no. 10, pp. 841–852, 2010.View at: Google Scholar
I. Cuiñas, A. M. Arias, A. G. Pino, and A. Ramos, “Educational software for single- and dual-reflector antennas,” Computer Applications in Engineering Education, vol. 7, no. 1, pp. 23–29, 1999.View at: Google Scholar
N. Morita, N. Kumagai, and J. R. Mautz, Integral Equation Methods for Electromagnetics, Artech House, Norwood, Mass, USA, 1987.
R. Harrington, “Boundary integral formulations for homogeneous material bodies,” Journal of Electromagnetic Waves and Applications, vol. 3, pp. 1–15, 1989.View at: Google Scholar
C. A. Balanis, Engineering Electromagnetics, John Wiley & Sons, New York, NY, USA, 1st edition, 1989.
Y. M. Wu, L. J. Jiang, and W. C. Chew, “An efficient method for computing highly oscillatory physical optics integral,” Progress in Electromagnetic Research, vol. 27, pp. 211–257, 2012.View at: Google Scholar