Solution to the integral in Eq. (31). Different cases are presented to overcome the singularities in the original expression.
In most of the electromagnetic problems, the number of unknowns to evaluate the scattered fields grows whenever the size of the antenna, device or scenario increases or the working frequency becomes higher. In this context, the rigorous full-wave methods –e.g. Method of Moments (MoM), Fast Multipole Method (FMM) (Engheta et al., 1992), Finite-Difference Time-Domain (FDTD) (Taflove & Umashankar, 1987) or Finite-Difference Frequency-Domain (FDFD) (Rappaport & McCartin, 1991), Finite Element Method (FEM) (Kempel et al., 1998) – can not tackle the analysis of such problems beyond an upper limit determined by the computational requirements in terms of time and memory. High frequency techniques consist in the asymptotic evaluation of the Maxwell’s equations. As a consequence, they provide good accuracy when dealing with electrically large geometries meanwhile the computational needs diminish with respect to the aforementioned methods.
Within the high frequency techniques, the Geometrical Optics (GO) and the Physical Optics (PO) approximation are the most extended methods due to the successful results obtained in various fields such as Radar Cross Section (RCS), design of reflector antennas or radioelectric coverage calculation. Since the Physical Optics approximation is detailed in the following section, the Geometrical Optics is briefly summarised.
The main interest in the GO lies in the fact that incident, reflected and transmitted electromagnetic waves are studied based on the conservation of the energy flux along a ray tube between a source and an observation point. Therefore, the Geometrical Optics is usually referred to as Ray Optics. The GO comprises two different methodologies (Rossi & Gabillet, 2002): Ray Tracing (Glassner, 1989) – the starting point is the receiver or observation point and a path to the source is sought analysing the reflections on walls, buildings, mountains – and Ray Launching – multiple rays are launched from the source, so they are independently followed until an observation point or the receiver is reached. One of the common applications of the GO is the evaluation of radio electric coverage or the channel characterization in urban scenarios.
Both the GO and the PO techniques require of an additional method to compute the contribution due to the diffraction phenomenon. The GO can be complemented by means of the Geometrical Theory of Diffraction (GTD) (Keller, 1962) or the Uniform Theory of Diffraction (UTD) (Pathak & Kouyoumjian, 1974). On the other hand, the Physical Theory of Diffraction (PTD) (Ufimtsev, 1962) is applied in joint with the PO formulation.
This chapter focuses on the PO approximation and especially on its extension to dielectric and lossy materials, namely the Modified Equivalent Current Approximation (MECA) method. In section 2, the PO approximation is presented and the formulation is obtained from the equivalence principle. Then the MECA method is introduced, so the equivalent current densities are obtained at the end of section 3.1.3. In order to complete the expressions, the reflection coefficients are then calculated and the determination of the electromagnetic field levels at the observation points based on an analytical solution of the radiation integral is accomplished in section 3.2. A couple of validation examples are studied before continuing to describe the algorithms for solving the visibility problem and, specifically, the Pyramid method in section 4. Afterwards, some of the applications of MECA are listed in “
2. The Physical Optics (PO) approximation
The Physical Optics (PO) approximation (Harrington, 2001; Balanis, 1989) is a well-known high frequency technique based on the determination of the equivalent current densities induced on the surface of an illuminated perfect electric conductor (PEC) plane. The final expressions can be achieved through the equivalence principle as it will be shown in section 2.1, so once these equivalent or PO current densities are obtained, both electric and magnetic field levels can be calculated from the corresponding radiation integrals.
In order to be in compliance with the constraints introduced by the PO approximation, some aspects have to be evaluated prior to the selection of this technique to tackle any electromagnetic problem:
The geometry must be made of electrically large obstacles with a smooth variation of their surfaces.
Given a radiating source, a distinction between lit and non-lit regions must be possible to perform.
The latter implies the need of a complementary algorithm to identify line of sight (LOS) directions from a specific point of view. Some methods will be commented in section 4.
Henceforth, for the sake of simplicity, a time-harmonic variation,, of the electromagnetic fields and current densities is assumed. Likewise, the spatial dependence of those variables is not explicitly written.
2.1. PO formulation
The expressions for the electric and magnetic equivalent current densities, and respectively, can be derived from the physical equivalent (Balanis, 1989) in Fig. 1(left), where a scheme consisting of two different media is presented: a PEC on the left and a medium with permittivity ε1, permeability μ1 and conductivity σ1 on the right. Consequently, the boundary
The incident electric and magnetic fields due to external sources and in absence of any obstacle are and. On the other hand, the total fields inside the PEC are null (), while in the second medium and are calculated by adding those incident fields to the reflected ones denoted by and. Therefore the electric and magnetic induced current densities, and, at the boundary S can be obtained from the tangential components of the total fields as:
Fig. 1 (centre) shows the equivalent problem for the non-PEC medium. In order to keep the same boundary conditions, this medium has been extended replacing the PEC part and the total fields are now and at the left side of S and and at its right side.
Considering the characterization of the boundary
At this point, an additional consideration to complete the PO formulation is included: when dealing with finite geometries, the PO current density is null in the regions not illuminated by the source:
This is the reason why the distinction between shadowed and illuminated parts of the scenario is one of the aforementioned constraints to correctly apply the PO approximation.
3. The Modified Equivalent Current Approximation (MECA) method: a PO extension for penetrable and non-metallic objects
Even though several electromagnetic problems fulfil the restrictions by the PO approximation in terms of electric size or radius of curvature, the characterization of the obstacles as perfect electric conductors reduces the scope of application of this high frequency technique.
Previous works (Rengarajan & Gillespie, 1988; Hodges & Rahmat-Samii, 1993, Sáez de Adana et al., 2004) agree that the extension of the PO approximation for penetrable and non-metallic objects have to account for the reflection coefficients
where and are the unit vectors in the direction of the TE and TM components of the incident electric field respectively. Likewise, the incident magnetic field can be written in an analogous way.
The Modified Equivalent Current Approximation (MECA) method (Meana et al., 2010) is a new high frequency technique based on the evaluation of a set of equivalent currents to calculate electromagnetic field levels at any observation point. The most important features of this method are:
MECA deals with electrically large scenarios consisting of dielectric and lossy surfaces due to both electric and magnetic equivalent currents, and, are taken into account.
MECA reduces to the PO expressions when considering PEC obstacles. Therefore, it can be seen as an extension and an improvement of the classic Physical Optics approximation.
Reflection coefficients are calculated in a rigorous manner.
A constant amplitude and linear phase current density distribution on the surface of the facets is adopted to analytically solve the radiation integral.
3.1. MECA formulation
Let us suppose a skew plane wave impinging on the interface with an angle of incidence denoted by θinc and a unit incident vector. The two media are now characterised by their respective permittivities (,), permeabilities (,) and conductivities (,). The analysis of the TE and TM components will be presented separately.
3.1.1. TE component
For the incident wave, the coordinate system where
has been chosen (see Fig. 2). Substituting by the unit reflection vector, the coordinate system is defined. Then, the incident electric and magnetic fields are rewritten in terms of these unit vectors:
where the intrinsic impedance of the first medium is,.
From the Snell law, the angle of reflection θref is equal to θinc. Therefore, can be expressed as:
Inserting the TE reflection coefficient
and introducing Eq. (8) in Eq. (9), the total fields for the TE component are:
3.1.2. TM component
In Fig. 3 the same directions of incidence, reflection and transmission as well as the boundary
3.1.3. MECA equivalent currents
Adding the contributions of the TE and TM components in Eqs. (10) and (11), the total electric and magnetic fields can be written as a function of the incident electric field, the reflection coefficients, the intrinsic impedance of the first medium, the propagation vectors and the unit vector:
Going back to the boundary conditions, the MECA equivalent current densities are reckoned from the total fields at
As it can be seen and remarked in the previous expressions, the amplitude and phase of the current densities depend on the polarization and also on the angle of incidence θinc of the incident wave.
3.1.4. TE and TM reflection coefficients
In order to complete the evaluation of the MECA current densities,
where denotes any propagation vector, the subscript
The continuity condition at
where is the wavelength in the first medium and
Inserting Eqs. (16) and (17) in Eq. (15), the real and imaginary parts of the third component of the transmission propagation vector are obtained as:
Once is determined,
and the corresponding
In case of working with non-metallic obstacles, the reflection coefficients become the well-known and simpler Snell’s expressions where the calculation of
When the second medium is a PEC (), the absolute value of the third component of tends to infinite and consequently as well as. Substituting these values in Eq. (12), the total magnetic field equals at the boundary
3.2. Electromagnetic field levels
Once and can be calculated with the inclusion of the reflection coefficients in the previous section, the evaluation of the scattered fields due to this set of equivalent current densities on the surface of an obstacle requires solving a radiation integral. For arbitrary distributions and geometries, an analytical solution can not be reached. As a consequence, two assumptions are imposed (Arias et al. 2000; Lorenzo et al., 2005):
The incident wave is a plane wave.
The obstacle consists of flat triangular facets.
With the objective of distinguishing between the reflected fields at the boundary and the scattered fields at the observation points, the latter are denoted by and. Note that these do not take into account the contribution of the diffraction phenomenon.
3.2.1. Incident plane wave approximation
The position vector of a point on the surface of a flat triangular facet can be written as:
where is the position vector of the barycentre of the i-th facet and is the vector from the barycentre to the source point (see Fig. 4).
The current densities are particularised for the i-th facet, and, having constant amplitudes,and and a linear phase distribution which varies with the unit direction of the Poynting vector as:
where and are the incident electric and magnetic fields at respectively, and * denotes the complex conjugate value of the complex vector in brackets. Therefore, in a medium where the working wavelength is 1,
The reason to introduce the linear phase variation is that this technique allows employing larger facets than when assuming a constant phase distribution.
3.2.2. Application to flat triangular facets
In order to analytically evaluate the radiation integral, the observation points are supposed to be in the far field of the radiating facet – but not necessarily in the far field of the whole scenario. Most of times this constraint is not an additional restriction because applications, such us radar cross section, directly impose huge distances between the object under test and the coordinates where the fields are determined.
Mathematically, the magnitude of the vector between the source point and the observation point is approximated (Balanis, 1989) by the magnitude of the latter. The phase is taken as, where.
The formulation to compute the scattered fields can be derived from the Maxwell’s equations, e.g. using the magnetic vector potential and the electric potential vector separately and combining both solutions. The resulting expressions are written in terms of the Green’s function in unbounded media and its gradient. The far field approximation allows expressing these functions as:
Consequently, the scattered fields are:
and they can be simplified inserting the auxiliary vectors and
Considering that the expressions of and are quite similar, the following reasoning is only performed for the magnetic fields, so the contribution due to the electric current density on the i-th facet at the observation point is given by:
In case of dealing with flat triangular facets, can be expressed as a function of the vectors in Fig. 4:
with, being the position vector of the vertices of the facet and, the position vector of the barycentre. and are real coefficients. Taking into account this transformation, the integral in Eq. (29) is converted into:
where is the area of the facet and the values of and are defined as:
The formulation of the Modified Equivalent Current Approximation method is completed once all the intermediate results that have been described in this section are combined to evaluate the scattered fields at a specific observation point because of one radiating facet. Until this point, the steps to follow are summarised in:
Decompose the incident electric field into its TE and TM components.
Obtain the reflection coefficients in Eqs. (19) and (20).
Calculate the MECA electric and magnetic current densities in Eq. (13).
Particularise the expressions for the barycentre of every radiating facet.
Compute the integral with the help of Table 1.
Evaluate the auxiliary fields and introducing the expressions of the MECA current densities and the result of the radiation integral.
Calculate the scattered fields and
Afterwards, the superposition theorem is employed to add the contribution due to all the radiating facets in the geometry for the observation point
If an evaluation of the line of sight between the radiating sources and the obstacle were necessary because of its shape, this task would be faced at the beginning of the whole.
3.3. Validation schemes
In this section, two canonical examples for different constitutive parameters and geometries are presented: one assures the good behaviour in high frequency and the other consists in the analysis of electrically large surfaces. In the first scenario, the validation is accomplished in terms of the Radar Cross Section (RCS) which is calculated as:
Since the field levels due to the diffraction contribution are not significant in the selected examples, this effect can be neglected in the computation of the RCS. Consequently, the term in the numerator only takes into account the electric field levels owing to reflections on the surface.
The methods to contrast the results provided by MECA are an analytical solution taken from the references and the full-wave technique the Method of Moments (MoM) (Medgyesi-Mitschang et al., 1994).
3.3.1. High frequency behaviour
As a high frequency technique, MECA has to show an accurate behaviour when the dimensions of the obstacle are large in comparison with the working wavelength. A good example to test this is a sphere, whose monostatic RCS can be obtained theoretically for both PEC (Balanis, 1989) and non-metallic (Van-Bladel, 2007) characterisations. A frequency sweep is performed for a sphere of radius with relative permittivity and permeability of, respectively and conductivity of, being, for the lossy case. The region in Fig. 5 for corresponds to the high frequency zone where a great coincidence in the results can be observed.
3.3.2. Electrically large surfaces
In this setup, a resonant horizontal dipole is placed at a fixed position whose height is of
For the working frequency of 1800 MHz and the values of ranging from 1 to 35 metres, the magnitude of the total electric field, including the direct illumination from the source, is evaluated. In Fig. 6 MECA results are compared with those by MoM (Medgyesi-Mitschang, 1994) for the following constitutive parameters:, and. On the left, is plotted as a function of υ, while on the right the variable is the logarithm of the distance in wavelengths, , so the small deviation –less than 2dB in the worst case- can be appreciated. Consequently, a high degree of overlapping in the curves of MoM and MECA clearly demonstrates the accuracy of the high frequency technique.
4. Fast visibility algorithm for solving the visibility problem
Because MECA (also PO) predicts null equivalent current densities in shadowed regions and non-null equivalent current densities in those regions with line of sight from the source point, a distinction between non-illuminated and illuminated surfaces is mandatory prior to the evaluation of electromagnetic levels at the observation points. This distinction is known as the visibility problem and it is widely found in computer graphics (Dewey, 1988; Foley, 1992; Bittner & Wonka, 2003), e.g. virtual reality or games, but also in the context of electromagnetic problems such as radio electric coverage evaluation or tracking applications in radar.
Within the classic algorithms, the Z-buffer and the Painter’s algorithm can be cited. Assuming the screen in the XY plane, both are based on the identification of the z-coordinate of the elements in the geometry to finally project the nearest parts of the objects. Another method is the Binary Space Partitioning (BSP) (Fuch et al. 1980; Gordon & Chen, 1991) which recursively divides the space into front and back semi-spaces, creating a binary tree at the same time. Afterwards, and in accordance with the target direction, the tree is walked and a priority list of the facets can be built.
Most of the techniques to solve the visibility problem are thought to project an image onto a display unit, so the algorithm itself has not better resolution than the pixel size. Therefore, an error is introduced, although almost imperceptible to the human eye. In order to overcome it, solutions where no approximations are made can be implemented. For example, in the Trimming method (Meana et al., 2009) when a piece of surface is partially occluded by other one, only the region in shadow is trimmed and removed from the original geometry. At the end, the remaining surface is the exact part that can be seen from a specific point of view. This type of algorithms is very time consuming but its precision counteracts this disadvantage.
To speed up the computation of the existing techniques, some modifications in the phase of preprocessing can be inserted. One option to carry out this is to split the scenario into parallelepipedic or conic macrodomains. Then the facets are classified in one or other macrodomain depending on their position. As a consequence, a preliminary discrimination is accomplished by blocks, so the facets are discarded faster. In this category, the Angular Z-buffer (AZB) (Cátedra et al., 1998) or the Space Volumetric Partitioning (SVP) (Cátedra & Arriaga, 1999) is included.
Additional acceleration can be achieved when the same algorithms are developed in Graphic Processing Units (GPUs) instead of Central Process Units (CPUs) (Ricks & Kuhlen, 2010). Fortunately, this is not restricted to visibility algorithms but it is suitable for all kind of algorithms consisting of loops, matrix operations, etc. As a consequence, generic GPUs are being extensively utilised to run the codes with smaller costs in comparison with CPU parallelisation. Even more, thanks to the use of libraries, e.g. DirectX or OpenGL (Shreiner, 2004) or the interaction between Matlab® and the graphic cards by means of the free plug-in by Nvidia® or Jacket by AccelerEyes®, the implementation of the routines turns into less complex programming.
In the following section the Pyramid method (Meana et al., 2009) is described due to its suitability to deal with flat triangular facets. Briefly, a right pyramid is built so that its vertex coincides with the source point. The i-th wall is the plane containing the source point and the i-th edge of the facet under test. Thus, a point behind this facet and inside the walls of the pyramid is occluded. This is a fast operation that is accomplished by substituting the coordinates of the observation point in the equation of the plane for all the walls.
4.1. Pyramid method
The scenario consists of facets which have been sorted by their distance from the source point to their barycentres, so the closest facet, denoted by, is always seen. Additionally, the origin of the coordinate system has been moved to the coordinates of.
In order to know whether a facet occludes a generic point in the geometry, the plane in which is contained is defined:
where is the unit outward normal vector and is the barycentre of. Similarly, the three planes containing one of the edges of the and are written as:
with being the position vector of the i-th vertex () of. These planes constitute the walls of the pyramid whose base is the facet under test. There are two conditions to conclude that can not be seen from the source point:
where and the vertices are supposed to be in clockwise or counterclockwise order.
The algorithm determines the occlusion of points by a triangular facet as depicted in Fig. 8, but the Pyramid method can be applied to determine the existence of line of sight among facets. The process to perform this is based on simplifying the representation of the triangles choosing some important points like the barycentre, vertices, inner points and executing the algorithm for them.
Even though the expressions have been particularised for the case of triangular facets, the implementation can be accomplished for any polygon by including additional walls in the pyramid and substituting the upper limit of the summatory in Eq. (38).
A quite similar algorithm to the Pyramid method is the Cone method (Meana et al., 2009), but a right cone is defined instead of a right pyramid. Its radius is a mean value of the distances from the barycentre of the facet under test to each of its vertices. Analogically, any point behind that facet and inside the cone is occluded. This is a faster algorithm because only one operation comparing the cone angle and the angle between the observation point and the axis of the cone is carried out. On the other hand, the exactness diminishes.
5. Application examples
The scope of application of the Modified Equivalent Current Approximation comprises all the problems which have been analysed by means of the Physical Optics traditionally, but also the extension to dielectric and lossy materials. In the following paragraphs some of these fields where MECA could be employed are summarised.
The RCS computation is one of the most cited topics in the literature: from canonical geometries (e.g. spheres, flat plates and dihedrals) (Griesser & Balanis, 1987; Ross, 1966) for contrasting and validating results to electrically large random surfaces and complex targets (ships or airplanes) (Adana et al., 2000; Uluisik et al., 2008) with a remarkable decrease in the computational time in comparison to full wave techniques. The consideration of absorbing materials allows studying the reduction in the radar signature of aircrafts or missiles.
In order to deal with more realistic scenarios in the RCS computation, the analysis of open-ended cavities including reflections and resonances (Burkholder & Lundin, 2005) is also a line of investigation MECA can face. Likewise, MECA has proven a good accuracy when dealing with electrically large rough surfaces (Meana et al., 2010) that can model the sea surface or the orography of a rural environment. Therefore, the evaluation of the electromagnetic levels for RCS or other applications can consider not only the target itself but also its surroundings.
Another field of application is the design of reflector antennas (Boag & Letrou, 2003; Lorenzo et al., 2005) with or without dielectric radomes because of their electric size. They are usually fed by an array of antennas whose radiation pattern is determined to compute the equivalent current densities and to calculate the fields at distances much larger than the wavelength afterwards.
On the other hand, the Physical Optics approximation and its extensions can also provide some reference values to validate new algorithms or the results by some imaging or shape reconstruction techniques (Saeedfar & Barkeshli, 2006). Due to the fact that the illuminating sources are known a priori, the simulation of a set of constitutive parameters and the shape of the objects can be tackled in a fast manner.
Although the high frequency techniques have lost relevance in the recent years due to the Fast Multipole Method (FMM) –developed by Rohklin in the field of acoustic dispersion (Rohklin, 1985; Rohklin, 1990) and then extended to electromagnetic dispersion (Engheta et al., 1992) – in joint with multilevel schemes and other acceleration methodologies, the advances in the terahertz band will become notable again in the foreseeable future. This means similar applications but working at higher frequencies where the wavelength is smaller than 1 mm and, as a consequence, most of the obstacles and objects is electrically large. Some papers have already proven the validity of the asymptotic approximations at these frequencies.
In addition to the previously enumerated fields of application of the MECA method and in spite of the typical employment of Ray Tracing and Launching in joint with empirical models to deal with the radioelectric coverage evaluation (Papkelis et al., 2007), this problem must be added to that list. The reason is that it fulfils all the requisites detailed through this chapter. In the next section an application example will be shown.
5.1. Radioelectric coverage
Network planning and optimization in rural and urban environments can be studied for different radiocommunication systems (
The radio electric stations are located at their real emplacements (see Fig. 10) simulating their electromagnetic behaviour based on the available data, such as configuration parameters (e.g. power, gain), provided by the telecommunications company. As a consequence, the radiation patterns were synthesised assuming a dipole array with a different number of elements for every different antenna. Once the illumination from the sources is explicitly determined, the MECA equivalent current densities are computed and the electromagnetic field levels are calculated.
The power density at the observation points situated at 1.5 metres above the terrain is represented in Fig. 10. A threshold has been obtained by taking into account the sensibility of the receivers, so that the minimum in the colorbar has been fixed in a lower value. This means that any region in dark blue can not initiate or maintain a communication with the base station. In order to improve the coverage in specific zones, new base stations could be added and a fast evaluation can be accomplished by computing only the direct illumination due to those stations.
In this chapter “
One of the most relevant features of MECA is that the reflection coefficients are calculated from the general dispersion relation for the transversal electric and transversal magnetic components independently. Therefore any conductivity and relative permittivity and permeability are allowed. Then, they are inserted in the corresponding TE/TM analysis which was shown in Fig. 2 and Fig. 3, to be finally combined in order to obtain the MECA equivalent current densities. With the objective of determining the electromagnetic fields at the observation points analytically, the incident wave is supposed to be a plane wave which impinges on the surface and generates a current density distribution with constant amplitude and linear phase variation. Assuming a flat triangular facet, the radiation integral can be solved by parts.
The good behaviour was proven in the validation examples, where the results from the frequency sweep in the high frequency region agreed the theoretical values. Likewise, an excellent overlapping was obtained for different angles of incidence when dealing with a non-PEC electrically large surface.
Because one of the constraints to employ PO and MECA is the determination of line of sight between the source and the observation points, some algorithms to solve the visibility problem were described. The classic methods are complemented by acceleration techniques and then, they are translated into the GPU programming languages. The Pyramid method was explained as an example of fast algorithm which was specifically developed for evaluating the occlusion by flat facets. Undoubtedly, this can be employed in joint with the MECA formulation, but the Pyramid method can also be helpful in other disciplines of engineering.
Throughout the section “