# Fast Far Field Computation of Single and Dual Reflector Antennas

Authors: Marcos Arias-Acuña; Antonio García-Pino; Oscar Rubiños-López
Journal: Journal of Engineering
 (2013)

---

## Abstract

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.


---

## Body

## 1. 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 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 [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 Section3, a new physical interpretation of the triangular patch current integration is provided.In Section4, 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–13].In Section5 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 Section6 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).


## 2. 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].


### 2.1. MFIE Field Formulation

The currentJ-(r-S) on a PEC with smooth and closed surface can be obtained using
(1)J-(r-S)=2n^(r-S)×(H-i(r-S)+⨍S′J-(r-S′)×∇′g(r-S,r-S′)dS′),
where r-S is the observation point, r-S′ is the source point, H-i(r-S) is the incident magnetic field at point r-S, symbol ⨍ means the principal value of the integral, and g(r-S,r-S′) is free space Green function, defined as
(2)g(r-S,r-S′)=e-jk|r-S-r-S′|4π|r-S-r-S′|.Taking into account (2), (1) can be rewritten as
(3)J-(r-S)=2n^(r-S)×(H-i(r-S)+⨍S′J-(r-S′)×R-S1+jkRSRSe-jkRS4πRS2dS′),
where R-S=r-S-r-S′ and RS=|r-S-r-S′|.After obtaining equivalent currents, the radiated magnetic field at pointsr- not belonging to the surface can be obtained by
(4)H-s(r-)=14π∫S′‍J-(r-S′)×R^1+jkRR2e-jkRdS′,
where R-=r--r-S′, R=|r--r-S′| and R^=R-/R.


### 2.2. PO Field Formulation

Physical optics (PO) is a simplification of (3) when the surface is flat enough to consider that the currents J-(r-S′) are on the same plane as the vector position R-S pointing to another point on the surface. In this case, their cross-product J-(r-S′)×R-S is colinear with the normal vector n^(r-S) and their cross-product tends to zero. With this approximation, (3) is simplified as
(5)J-PO(r-S)=2n^(r-S)×H-i(r-S)
and the radiated magnetic field at points not belonging to the surface can be obtained with physical optics by considering
(6)H-s(r-)=14π∫S′‍J-PO(r-S′)×R^1+jkRR2e-jkRdS′.When the observation point is at a distancer=|r-|≫|r-S′| for every source point, then the distance can be simplified as follows:
(7)R=|r--r-S′|=r2+r-S′2-2r-·r-S′≈r-r^·r-S′,
where the term r^·r-S′ 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
(8)H-s(r-)=14π(1+jkr)·e-jkrr2[∫S′‍J-PO(r-S′)ejkr^·r-S′dS′]×r^.If the distancer=|r-|≫λ, then (8) can be simplified by
(9)H-s(r^)=j2λe-jkrr[∫S′‍J-PO(r-S′)ejkr^·r-S′dS′]×r^.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 pointr-k can be obtained from (8) by dividing the surface in subdomains (patches) small enough that we can consider that the point r-k is far from each patch. The magnetic field H-k at the point r-k can be represented as
(10)H-k=∑i‍H-ik,
where H-ik is the magnetic field radiated by the patch i at r-k. The field radiated by each patch will be
(11)H-ik=14π(1+jkrik)·e-jkrikrik2[∫Si‍J-PO(r-S′)ejkr^ik·r-S′dSi]×r^ik,
where Si represents the surface of the considered patch i.PO approximation can be used when the patch surface is flat enough to consider that the whole patch has the same normal vectorn^i. Inserting this vector into (11) produces the following expression:
(12)H-ik=12π(1+jkrik)·e-jkrikrik2[n^i×∫Si‍H-i(r-S′)ejkr^ik·r-S′dSi]×r^ik.


### 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^i. Consider that
(13)r-S′=r-i+r-".
Then,
(14)H-i(r-S′)≈H-i0·e-jkr-"·p^i.Including (13) and (14) into (12),
(15)H-ik≈12π(1+jkrik)·e-jkrikrik2[n^i×H-i0·Ii]×r^ik,
where
(16)Ii=∫Si‍ejk(r^ik-p^i)·r-"dSi.Then, the scattered far field will now be obtained as(17)H-(r^)=jλe-jkrr[∑i‍n^i×H-i0·Ii·ejkr^·r-i]×r^.


### 2.5. Triangular Patches

Let us consider a triangular flat patch as seen in Figure1. The triangle is defined by three points P-1, P-2, and P-3. P-0 is the point where H-i0 is calculated which is assumed to be at the barycenter P0=(P1+P2+P3)/3. If we define v-mn=P-n-P-m, the local coordinates of the triangle can be described with two variables u and v as
(18)r-"=v-01+u·v-12+v·v-13.Figure 1
Triangular points and vectors.The normal vector and the area of the triangle are related to the cross-product ofv-12 and v-13 as follows:
(19)v-12×v-13=2Ai·n^i.Using this coordinate system, (16) is then given by
(20)Ii=2Ae-j((α+β)/3)∫01‍∫01-u‍ej(αu+βv)dvdu,
whose solution is
(21)Ii(r^k)=2Ae-j((α+β)/3)[αejβ-βejα+β-α(α-β)αβ],
where
(22)α=kv-12·(r^ik-p^i),(23)β=kv-13·(r^ik-p^i).Equation (21) has the following singular values:
(24)α=0⇒Ii(r^k)=2Ae-j(β/3)[1+jβ-ejββ2],β=0⇒Ii(r^k)=2Ae-j(α/3)[1+jα-ejαα2],α=β⇒Ii(r^k)=2Ae-j(2α/3)[ejα(1-jα)-1α2],α=β=0⇒Ii(r^k)=A.Form small values ofα and β, a second order series approximation can be used:
(25)Ii(r^k)≈A(1-α2+β2-αβ36).Finally, for directions close to the main lobe or for very small patches, the expression (25) could be approximated as
(26)Ii(r^k)≈A.


## 2.1. MFIE Field Formulation

The currentJ-(r-S) on a PEC with smooth and closed surface can be obtained using
(1)J-(r-S)=2n^(r-S)×(H-i(r-S)+⨍S′J-(r-S′)×∇′g(r-S,r-S′)dS′),
where r-S is the observation point, r-S′ is the source point, H-i(r-S) is the incident magnetic field at point r-S, symbol ⨍ means the principal value of the integral, and g(r-S,r-S′) is free space Green function, defined as
(2)g(r-S,r-S′)=e-jk|r-S-r-S′|4π|r-S-r-S′|.Taking into account (2), (1) can be rewritten as
(3)J-(r-S)=2n^(r-S)×(H-i(r-S)+⨍S′J-(r-S′)×R-S1+jkRSRSe-jkRS4πRS2dS′),
where R-S=r-S-r-S′ and RS=|r-S-r-S′|.After obtaining equivalent currents, the radiated magnetic field at pointsr- not belonging to the surface can be obtained by
(4)H-s(r-)=14π∫S′‍J-(r-S′)×R^1+jkRR2e-jkRdS′,
where R-=r--r-S′, R=|r--r-S′| and R^=R-/R.


## 2.2. PO Field Formulation

Physical optics (PO) is a simplification of (3) when the surface is flat enough to consider that the currents J-(r-S′) are on the same plane as the vector position R-S pointing to another point on the surface. In this case, their cross-product J-(r-S′)×R-S is colinear with the normal vector n^(r-S) and their cross-product tends to zero. With this approximation, (3) is simplified as
(5)J-PO(r-S)=2n^(r-S)×H-i(r-S)
and the radiated magnetic field at points not belonging to the surface can be obtained with physical optics by considering
(6)H-s(r-)=14π∫S′‍J-PO(r-S′)×R^1+jkRR2e-jkRdS′.When the observation point is at a distancer=|r-|≫|r-S′| for every source point, then the distance can be simplified as follows:
(7)R=|r--r-S′|=r2+r-S′2-2r-·r-S′≈r-r^·r-S′,
where the term r^·r-S′ 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
(8)H-s(r-)=14π(1+jkr)·e-jkrr2[∫S′‍J-PO(r-S′)ejkr^·r-S′dS′]×r^.If the distancer=|r-|≫λ, then (8) can be simplified by
(9)H-s(r^)=j2λe-jkrr[∫S′‍J-PO(r-S′)ejkr^·r-S′dS′]×r^.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 pointr-k can be obtained from (8) by dividing the surface in subdomains (patches) small enough that we can consider that the point r-k is far from each patch. The magnetic field H-k at the point r-k can be represented as
(10)H-k=∑i‍H-ik,
where H-ik is the magnetic field radiated by the patch i at r-k. The field radiated by each patch will be
(11)H-ik=14π(1+jkrik)·e-jkrikrik2[∫Si‍J-PO(r-S′)ejkr^ik·r-S′dSi]×r^ik,
where Si represents the surface of the considered patch i.PO approximation can be used when the patch surface is flat enough to consider that the whole patch has the same normal vectorn^i. Inserting this vector into (11) produces the following expression:
(12)H-ik=12π(1+jkrik)·e-jkrikrik2[n^i×∫Si‍H-i(r-S′)ejkr^ik·r-S′dSi]×r^ik.


## 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^i. Consider that
(13)r-S′=r-i+r-".
Then,
(14)H-i(r-S′)≈H-i0·e-jkr-"·p^i.Including (13) and (14) into (12),
(15)H-ik≈12π(1+jkrik)·e-jkrikrik2[n^i×H-i0·Ii]×r^ik,
where
(16)Ii=∫Si‍ejk(r^ik-p^i)·r-"dSi.Then, the scattered far field will now be obtained as(17)H-(r^)=jλe-jkrr[∑i‍n^i×H-i0·Ii·ejkr^·r-i]×r^.


## 2.5. Triangular Patches

Let us consider a triangular flat patch as seen in Figure1. The triangle is defined by three points P-1, P-2, and P-3. P-0 is the point where H-i0 is calculated which is assumed to be at the barycenter P0=(P1+P2+P3)/3. If we define v-mn=P-n-P-m, the local coordinates of the triangle can be described with two variables u and v as
(18)r-"=v-01+u·v-12+v·v-13.Figure 1
Triangular points and vectors.The normal vector and the area of the triangle are related to the cross-product ofv-12 and v-13 as follows:
(19)v-12×v-13=2Ai·n^i.Using this coordinate system, (16) is then given by
(20)Ii=2Ae-j((α+β)/3)∫01‍∫01-u‍ej(αu+βv)dvdu,
whose solution is
(21)Ii(r^k)=2Ae-j((α+β)/3)[αejβ-βejα+β-α(α-β)αβ],
where
(22)α=kv-12·(r^ik-p^i),(23)β=kv-13·(r^ik-p^i).Equation (21) has the following singular values:
(24)α=0⇒Ii(r^k)=2Ae-j(β/3)[1+jβ-ejββ2],β=0⇒Ii(r^k)=2Ae-j(α/3)[1+jα-ejαα2],α=β⇒Ii(r^k)=2Ae-j(2α/3)[ejα(1-jα)-1α2],α=β=0⇒Ii(r^k)=A.Form small values ofα and β, a second order series approximation can be used:
(25)Ii(r^k)≈A(1-α2+β2-αβ36).Finally, for directions close to the main lobe or for very small patches, the expression (25) could be approximated as
(26)Ii(r^k)≈A.


## 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:
(27)Ii(r^k)=-2A×[e(-α-β)/3α·β+e(2α-β)/3(β-α)·(-α)+e(-α+2β)/3(α-β)·(-β)].Each of the three terms of (27) can be interpreted as the contribution of a vertex to the radiation integral Ii. Taking into account the definitions of v-0n it can be stated that
(28)v-01=P-1-P-1+P-2+P-33=-v-12-v-133,v-02=P-2-P-1+P-2+P-33=2v-12-v-133,v-03=P-3-P-1+P-2+P-33=-v-12+2v-133.So, expression (27) can be transformed into
(29)Ii(r^k)=-2A[ev-01·(r^ik-p^i)a12·a13+ev-02·(r^ik-p^i)a23·a21+ev-03·(r^ik-p^i)a32·a31],
where each phase term ev-0n·(r^ik-p^i) is due to the displacement from the barycenter to the vertex n and the amplitude terms depend on the vectors from the vertex n to the other two vertices, according to the values of α and β given by (22) and (23):
(30)anm=kv-nm·(r^ik-p^i).


### 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^i or r^=p^i-2n^·(p^i·n^).If the triangle is big enough compared with the wavelength, then the main effect is the reflection, as it can be seen in Figure2 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∘.Figure 2
Reflection versus length of one triangle.


### 3.3. 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 radiusR on the plane XY with a plane wave normal to the circle as incident field. The scattered field for ϕ=0 will be
(31)Ii(θ)=∫0R‍∫02π‍ejkρcosϕ·sinθρdϕdθ=2πRJ1(kRsinθ)ksinθ
with a maximum value Ii(θ)|max=Ii(0)=πR2.Adapting (27) to this example, where the phase reference point will be at the center of the circle and v-12=R·ρ^ϕ and v-13=R·ρ^ϕ+Δϕ, this leads to αn=Rcosϕnsinθ and βn=Rcos(ϕn+Δϕ)sinθ. Then, the integral will be
(32)Ii(θ)=sinΔϕk2sin2θ×∑n=0N-1‍[ejkRsinθcos((n+1)Δϕ)-1cos((n+1)Δϕ)·(cos(nΔϕ)-cos((n+1)Δϕ))-ejkRsinθcos(nΔϕ)-1cos(nΔϕ)·(cos(nΔϕ)-cos((n+1)Δϕ))],
where Δϕ=2π/N.The limit of this expression tends to the integral:(33)Ii(θ)=1k2sin2θ∫02π‍ejkRsinθcosϕ(1-jkRsinθ)-1cos2ϕdϕ.Now, we are going to determine the number of angular sectorsNϕ that are necessary to obtain an accurate solution when the radius R varies. In Figures 3, 4, and 5 the convergence can be seen when the radius is 2λ, 10λ, and 50λ.Figure 3
Convergence whenR=2λ.Figure 4
Convergence whenR=10λ.Figure 5
Convergence whenR=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.


## 3.1. Vertex Contributions

Equation (21) can be reformulated in order to clarify its physical meaning:
(27)Ii(r^k)=-2A×[e(-α-β)/3α·β+e(2α-β)/3(β-α)·(-α)+e(-α+2β)/3(α-β)·(-β)].Each of the three terms of (27) can be interpreted as the contribution of a vertex to the radiation integral Ii. Taking into account the definitions of v-0n it can be stated that
(28)v-01=P-1-P-1+P-2+P-33=-v-12-v-133,v-02=P-2-P-1+P-2+P-33=2v-12-v-133,v-03=P-3-P-1+P-2+P-33=-v-12+2v-133.So, expression (27) can be transformed into
(29)Ii(r^k)=-2A[ev-01·(r^ik-p^i)a12·a13+ev-02·(r^ik-p^i)a23·a21+ev-03·(r^ik-p^i)a32·a31],
where each phase term ev-0n·(r^ik-p^i) is due to the displacement from the barycenter to the vertex n and the amplitude terms depend on the vectors from the vertex n to the other two vertices, according to the values of α and β given by (22) and (23):
(30)anm=kv-nm·(r^ik-p^i).


## 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^i or r^=p^i-2n^·(p^i·n^).If the triangle is big enough compared with the wavelength, then the main effect is the reflection, as it can be seen in Figure2 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∘.Figure 2
Reflection versus length of one triangle.


## 3.3. 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 radiusR on the plane XY with a plane wave normal to the circle as incident field. The scattered field for ϕ=0 will be
(31)Ii(θ)=∫0R‍∫02π‍ejkρcosϕ·sinθρdϕdθ=2πRJ1(kRsinθ)ksinθ
with a maximum value Ii(θ)|max=Ii(0)=πR2.Adapting (27) to this example, where the phase reference point will be at the center of the circle and v-12=R·ρ^ϕ and v-13=R·ρ^ϕ+Δϕ, this leads to αn=Rcosϕnsinθ and βn=Rcos(ϕn+Δϕ)sinθ. Then, the integral will be
(32)Ii(θ)=sinΔϕk2sin2θ×∑n=0N-1‍[ejkRsinθcos((n+1)Δϕ)-1cos((n+1)Δϕ)·(cos(nΔϕ)-cos((n+1)Δϕ))-ejkRsinθcos(nΔϕ)-1cos(nΔϕ)·(cos(nΔϕ)-cos((n+1)Δϕ))],
where Δϕ=2π/N.The limit of this expression tends to the integral:(33)Ii(θ)=1k2sin2θ∫02π‍ejkRsinθcosϕ(1-jkRsinθ)-1cos2ϕdϕ.Now, we are going to determine the number of angular sectorsNϕ that are necessary to obtain an accurate solution when the radius R varies. In Figures 3, 4, and 5 the convergence can be seen when the radius is 2λ, 10λ, and 50λ.Figure 3
Convergence whenR=2λ.Figure 4
Convergence whenR=10λ.Figure 5
Convergence whenR=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.


## 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 k will be a vector J-k. This simplifies the integral assuming that R-S≈(r-k-r-n)=r-nk and |r-nk|=rnk, so (3) can be obtained as
(34)J-k=2n^k×H-ki+12πn^k×∑n‍J-n×r-nkrnk2(jk+1rnk)e-jkrnkAn,
where An is the area of nth-triangle.The equivalent current in (34) is due to incident (PO) currents and to the coupling among surfaces. We can define
(35)J-kPO=2n^k×H-ki,(36)c-nk=r-nke-jkrnk2πrnk2(jk+1rnk)
and now the equivalent current can be expressed as
(37)J-k=J-kPO+n^k×∑n‍AnJ-n×c-nk,
where k and n are for points on different surfaces.When the currents on both surfaces are obtained, the scattered far fields can be calculated at a directionr^ using
(38)H-s(r^)=j2λe-jkrr∑i‍AiJ-iejkr-i·r^,
where i is intended for all the main reflector and subreflector patches and
(39)E-s(r^)=ηH-s(r^)×r^.


### 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, k will be for reflector patches and n 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:(40)J-k0=J-kPO=2n^k×H-ki,J-n1=J-nPO+n^n×∑k‍AkJ-k0×c-kn=2n^n×H-ni+n^n×∑k‍AkJ-k0×c-kn,J-kN=J-kPO+n^k×∑n‍AnJ-kN×c-nk,J-nN+1=J-nPO+n^n×∑k‍AkJ-kN×c-kn
where c-kn=-c-nk.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 vectorsu^1k and u^2k for kth-patch. Using triangular patches, these vectors can be obtained as
(41)v-12=P-2-P-1;v-13=P-3-P-1,n^k=v-12×v-13|v-12×v-13|,u^1k=v-12|v-12|;u^2k=n^k×u^1k.The currents defined in (37) can be expressed in terms of u^1k and u^2k as
(42)J-k=J1ku^1k+J2ku^2k
and now (37) can be divided into two equations:
(43)J1k=J-k·u^1k=J-kPO·u^1k-∑n‍u^2k·(AnJ-n×c-nk),J2k=J-k·u^2k=J-kPO·u^2k+∑n‍u^1k·(AnJ-n×c-nk).J
-
n can also be divided into two components. After some manipulations, (43) can be expressed as
(44)J1kPO=J1k+∑n‍J1n(u^2k×u^1n)·c-nk+∑n‍AnJ2n(u^2k×u^2n)·c-nk,J2kPO=J2k-∑n‍J1n(u^1k×u^1n)·c-nk-∑n‍AnJ2n(u^1k×u^2n)·c-nk.The mutual coupling between the surfaces can be expressed as(45)J1kPO=J1k+∑n‍AnJ1n(u^2k×u^1n)·c-nk+∑n‍AnJ2n(u^2k×u^2n)·c-nk,J1nPO=J1n+∑k‍AkJ1k(u^2n×u^1k)·c-kn+∑k‍AkJ2k(u^2n×u^2k)·c-kn,J2kPO=J2k-∑n‍AnJ1n(u^1k×u^1n)·c-nk-∑n‍AnJ2n(u^1k×u^2n)·c-nk,J2nPO=J2n-∑k‍AkJ1k(u^1n×u^1k)·c-kn-∑k‍AkJ2k(u^1n×u^2k)·c-kn.This equation system can be expressed with matrix formulation:(46)[0k×k¯¯M2k2n¯¯1AkIk×k¯¯M2k1n¯¯M2n2k¯¯0n×n¯¯M2n1k¯¯1AnIn×n¯¯1AkIk×k¯¯M1k2n¯¯0k×k¯¯M1k1n¯¯M1n2k¯¯1AnIn×n¯¯M1n1k¯¯0n×n¯¯]·[AkJ2k¯AnJ2n¯AkJ1k¯AnJ1n¯]=[J1kPO¯J1nPO¯J2kPO¯J2nPO¯],
where
(47)M1k1n¯¯=-(u^1k×u^1n)·c-nkM1n1k¯¯=-(u^1n×u^1k)·c-kn=M1k1n¯¯TM1k2n¯¯=-(u^1k×u^2n)·c-nkM2n1k¯¯=(u^2n×u^1k)·c-kn=-M1k2n¯¯TM2k1n¯¯=(u^2k×u^1n)·c-nkM1n2k¯¯=-(u^1n×u^2k)·c-kn=-M2k1n¯¯TM2k2n¯¯=(u^2k×u^2n)·c-nkM2n2k¯¯=(u^2n×u^2k)·c-kn=M2k2n¯¯T.
In (47), the first two subscripts of M 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:
(48)[J2k¯J1k¯]=[M1¯¯M2¯¯M2¯¯M3¯¯]-1·[J1kPO¯J2kPO¯],
where
(49)M1¯¯=-M2k2nAn¯¯·M1n2kAk¯¯-M2k1nAn¯¯·M2n2kAk¯¯,M2¯¯=Ik×k¯¯-M2k2nAn¯¯·M1n1kAk¯¯-M2k1nAn¯¯·M2n1kAk¯¯,M3¯¯=-M1k2nAn¯¯·M1n1kAk¯¯-M1k1nAn¯¯·M2n1kAk¯¯
and the components of the main reflector currents are calculated directly with (50). 								(50)J1n¯=-M2n2k¯¯·AkJ2k¯-M2n1k¯¯·AkJ1k¯,J2n¯=-M1n2k¯¯·AkJ2k¯-M2n1k¯¯·AkJ1k¯.
Finally, all the currents must be transformed to the general coordinate system, using (42).


## 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 k will be a vector J-k. This simplifies the integral assuming that R-S≈(r-k-r-n)=r-nk and |r-nk|=rnk, so (3) can be obtained as
(34)J-k=2n^k×H-ki+12πn^k×∑n‍J-n×r-nkrnk2(jk+1rnk)e-jkrnkAn,
where An is the area of nth-triangle.The equivalent current in (34) is due to incident (PO) currents and to the coupling among surfaces. We can define
(35)J-kPO=2n^k×H-ki,(36)c-nk=r-nke-jkrnk2πrnk2(jk+1rnk)
and now the equivalent current can be expressed as
(37)J-k=J-kPO+n^k×∑n‍AnJ-n×c-nk,
where k and n are for points on different surfaces.When the currents on both surfaces are obtained, the scattered far fields can be calculated at a directionr^ using
(38)H-s(r^)=j2λe-jkrr∑i‍AiJ-iejkr-i·r^,
where i is intended for all the main reflector and subreflector patches and
(39)E-s(r^)=ηH-s(r^)×r^.


## 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, k will be for reflector patches and n 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:(40)J-k0=J-kPO=2n^k×H-ki,J-n1=J-nPO+n^n×∑k‍AkJ-k0×c-kn=2n^n×H-ni+n^n×∑k‍AkJ-k0×c-kn,J-kN=J-kPO+n^k×∑n‍AnJ-kN×c-nk,J-nN+1=J-nPO+n^n×∑k‍AkJ-kN×c-kn
where c-kn=-c-nk.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 vectorsu^1k and u^2k for kth-patch. Using triangular patches, these vectors can be obtained as
(41)v-12=P-2-P-1;v-13=P-3-P-1,n^k=v-12×v-13|v-12×v-13|,u^1k=v-12|v-12|;u^2k=n^k×u^1k.The currents defined in (37) can be expressed in terms of u^1k and u^2k as
(42)J-k=J1ku^1k+J2ku^2k
and now (37) can be divided into two equations:
(43)J1k=J-k·u^1k=J-kPO·u^1k-∑n‍u^2k·(AnJ-n×c-nk),J2k=J-k·u^2k=J-kPO·u^2k+∑n‍u^1k·(AnJ-n×c-nk).J
-
n can also be divided into two components. After some manipulations, (43) can be expressed as
(44)J1kPO=J1k+∑n‍J1n(u^2k×u^1n)·c-nk+∑n‍AnJ2n(u^2k×u^2n)·c-nk,J2kPO=J2k-∑n‍J1n(u^1k×u^1n)·c-nk-∑n‍AnJ2n(u^1k×u^2n)·c-nk.The mutual coupling between the surfaces can be expressed as(45)J1kPO=J1k+∑n‍AnJ1n(u^2k×u^1n)·c-nk+∑n‍AnJ2n(u^2k×u^2n)·c-nk,J1nPO=J1n+∑k‍AkJ1k(u^2n×u^1k)·c-kn+∑k‍AkJ2k(u^2n×u^2k)·c-kn,J2kPO=J2k-∑n‍AnJ1n(u^1k×u^1n)·c-nk-∑n‍AnJ2n(u^1k×u^2n)·c-nk,J2nPO=J2n-∑k‍AkJ1k(u^1n×u^1k)·c-kn-∑k‍AkJ2k(u^1n×u^2k)·c-kn.This equation system can be expressed with matrix formulation:(46)[0k×k¯¯M2k2n¯¯1AkIk×k¯¯M2k1n¯¯M2n2k¯¯0n×n¯¯M2n1k¯¯1AnIn×n¯¯1AkIk×k¯¯M1k2n¯¯0k×k¯¯M1k1n¯¯M1n2k¯¯1AnIn×n¯¯M1n1k¯¯0n×n¯¯]·[AkJ2k¯AnJ2n¯AkJ1k¯AnJ1n¯]=[J1kPO¯J1nPO¯J2kPO¯J2nPO¯],
where
(47)M1k1n¯¯=-(u^1k×u^1n)·c-nkM1n1k¯¯=-(u^1n×u^1k)·c-kn=M1k1n¯¯TM1k2n¯¯=-(u^1k×u^2n)·c-nkM2n1k¯¯=(u^2n×u^1k)·c-kn=-M1k2n¯¯TM2k1n¯¯=(u^2k×u^1n)·c-nkM1n2k¯¯=-(u^1n×u^2k)·c-kn=-M2k1n¯¯TM2k2n¯¯=(u^2k×u^2n)·c-nkM2n2k¯¯=(u^2n×u^2k)·c-kn=M2k2n¯¯T.
In (47), the first two subscripts of M 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:
(48)[J2k¯J1k¯]=[M1¯¯M2¯¯M2¯¯M3¯¯]-1·[J1kPO¯J2kPO¯],
where
(49)M1¯¯=-M2k2nAn¯¯·M1n2kAk¯¯-M2k1nAn¯¯·M2n2kAk¯¯,M2¯¯=Ik×k¯¯-M2k2nAn¯¯·M1n1kAk¯¯-M2k1nAn¯¯·M2n1kAk¯¯,M3¯¯=-M1k2nAn¯¯·M1n1kAk¯¯-M1k1nAn¯¯·M2n1kAk¯¯
and the components of the main reflector currents are calculated directly with (50). 								(50)J1n¯=-M2n2k¯¯·AkJ2k¯-M2n1k¯¯·AkJ1k¯,J2n¯=-M1n2k¯¯·AkJ2k¯-M2n1k¯¯·AkJ1k¯.
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 [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.


### 5.2. Complex Exponential

Most of the execution time is used in the calculation of complex exponentials.For single reflector antennas, the exponential termejkr^·r-i 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:
(51)cospre(n)=cos(nN·4π);sinpre(n)=sin(nN·4π),(52)arg=[αλ-Int(αλ)]·N2+N2,(53)ejkα=cospre(arg)+sinpre(arg).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 obtainingc-nk 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-j/kr)·e-jkr/r2 has the following steps: (54)a=Nrmax2-rmin2;b=a·rmin2;r2=rmin2+na;rn=k·r2,(55)realpre(n)=cosrn-sinrn/rnr2;imagpre(n)=-sinrn+cosrn/rnr2,(56)arg=Int(r2·a-b),(57)e-jkrr2(1-jkr)=realpre(arg)+j·imagpre(arg),
where rmin is the minimum value of r and rmax the maximum value. The term r2 is calculated directly from the components rx2+ry2+rz2, 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.


### 5.3. 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.


## 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 [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.


## 5.2. Complex Exponential

Most of the execution time is used in the calculation of complex exponentials.For single reflector antennas, the exponential termejkr^·r-i 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:
(51)cospre(n)=cos(nN·4π);sinpre(n)=sin(nN·4π),(52)arg=[αλ-Int(αλ)]·N2+N2,(53)ejkα=cospre(arg)+sinpre(arg).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 obtainingc-nk 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-j/kr)·e-jkr/r2 has the following steps: (54)a=Nrmax2-rmin2;b=a·rmin2;r2=rmin2+na;rn=k·r2,(55)realpre(n)=cosrn-sinrn/rnr2;imagpre(n)=-sinrn+cosrn/rnr2,(56)arg=Int(r2·a-b),(57)e-jkrr2(1-jkr)=realpre(arg)+j·imagpre(arg),
where rmin is the minimum value of r and rmax the maximum value. The term r2 is calculated directly from the components rx2+ry2+rz2, 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.


## 5.3. 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.


## 6. 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-q 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 sectorsNr until the difference will be negligible. Then, the difference in directivity is accounted for in natural units, using
(58)Difference(%)=mean(|DNr-Dref|)max(|DNr-Dref|)·100.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.Difference using (26): at (a) the difference and at (b) the pattern.
(a)
(b)Difference using (25): at (a) the difference and at (b) the pattern.
(a)
(b)This difference with the convergence is measured in sectors with width0.1∘. With Nr=42 (area of the triangles of 8λ2), the difference is below 1%. With Nr=24 (area of the triangles of 25λ2), the difference is below 5%. And with Nr=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.Difference whenNr=42: at (a) the difference and at (b) the pattern.
(a)
(b)Difference whenNr=13: at (a) the difference and at (b) the pattern.
(a)
(b)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 Nr=42 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 c-nk 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.


## 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-q 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 sectorsNr until the difference will be negligible. Then, the difference in directivity is accounted for in natural units, using
(58)Difference(%)=mean(|DNr-Dref|)max(|DNr-Dref|)·100.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.Difference using (26): at (a) the difference and at (b) the pattern.
(a)
(b)Difference using (25): at (a) the difference and at (b) the pattern.
(a)
(b)This difference with the convergence is measured in sectors with width0.1∘. With Nr=42 (area of the triangles of 8λ2), the difference is below 1%. With Nr=24 (area of the triangles of 25λ2), the difference is below 5%. And with Nr=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.Difference whenNr=42: at (a) the difference and at (b) the pattern.
(a)
(b)Difference whenNr=13: at (a) the difference and at (b) the pattern.
(a)
(b)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 Nr=42 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 c-nk 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.


## 7. 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.



---

*Source: 140254-2013-03-24.xml*