Analysis of Electromagnetic Scattering from Large Arrays of Cylinders via a Hybrid of the Method of Auxiliary Sources (MAS) with the Fast Multipole Method (FMM)

Abstract

The Method of Auxiliary Sources (MAS) is an established technique for the numerical solution of electromagnetic (EM) scattering and radiation problems. This paper presents a hybrid of MAS with the Fast Multipole Method (FMM), which provides a strategy for reducing the computational cost and for solving large-scale problems without notable accuracy loss (and in a reasonable time). The hybrid MAS-FMM scheme is applied to the problem of EM scattering from an arbitrarily large array of lossless/lossy dielectric cylinders. Numerical results are presented to verify the MAS and MAS-FMM schemes, as well as to illuminate the improvements stemming from the proposed hybridization (especially the ones regarding the associated complexity and computational cost). A few concluding remarks offer a summary of this work, along with a list of possible future extensions.

1. Introduction

The Method of Auxiliary Sources (MAS) is an established numerical method in Computational Electromagnetics (CEM), and it has been applied extensively to a large variety of scattering and radiation problems [1]. The standard MAS can be classified in the so-called Generalized Multipole Techniques (GMTs) [2]. Among frequency-domain numerical methods, the MAS avoids time-consuming numerical integrations; its algorithmic implementation is essentially straightforward, and it is characterized by relatively low computational cost [3]. However, despite these advantageous features, certain issues related to the implementation of the MAS undermine its applicability and effectiveness [4,5]. These include difficulties associated with the selection of the location(s) and number(s) of the auxiliary sources, occasionally severe ill-conditioning effects of the associated matrices, as well as moderately slow convergence of the computed fields and possible divergence of the currents involved [4,5,6,7,8,9,10].

Despite the fact that the MAS is often considered to have a lower computational cost compared to the Method of Moments (MoM) and the Finite-Element Method (FEM), its complexity remains high (i.e., the complexity is calculated asymptotically as the number of unknowns approaches infinity). For this reason, the relative advantages of MAS may be weakened in applications requiring increased numerical stability and accuracy, which inevitably involve large numbers of unknowns. One possible alternative to reduce the MAS complexity and computational cost would be to somehow group the interactions between auxiliary sources and matching points and subsequently to approximate them by simpler formulas. This precisely constitutes the fundamental concept behind the Fast Multipole Method (FMM) and related strategies [11]. Exploiting such an alternative, hybridizations between MAS and FMM were proposed in [12] for perfectly conducting 2-D scatterers, in [13] for internal problems corresponding to single Helmholtz resonators, and in [14] for perfectly conducting 3-D scatterers.

In this paper, a MAS-FMM scheme is presented for the problem of electromagnetic scattering from and propagation through a large array of parallel lossless/lossy dielectric (circular) cylinders. An array of this type can stand as a two-dimensional (2D) model of a forest or an orchard. Admittedly, such an approximation may look naive at a first glance; however, exact simulation of EM propagation through vegetation is a formidable task. First of all, the actual geometry itself is almost arbitrary, including trunks, branches, and foliage with random orientation and variable dielectric properties. Moreover, the large number of scatterers involved implies that standard integral-equation methods are not applicable, due to the exorbitant number of unknowns required. Therefore, modeling tree trunks as infinite dielectric cylinders, where, additionally, the number of unknowns is somehow compressed, is a tractable way to go, producing fairly accurate and sufficiently realistic results for engineering purposes. Besides, problems involving arrays of dielectric cylinders occur in other areas pertaining to applied electromagnetics, including in optics for the analysis of modes excited in photonic fibers and gratings [15,16,17].

The simpler case of arrays of perfectly electric conducting (PEC) cylinders was presented and evaluated in [18]. In this work, pertaining to lossless/lossy dielectric cylinders, a standard MAS formulation is used as a starting point. Then, according to the FMM concept, far-field interactions are approximated as perturbations of a central interaction between the centroids of properly selected clusters encompassing neighboring auxiliary sources and matching points. The accuracy and efficiency of the proposed methods are examined and assessed. The credibility of both the MAS and MAS-FMM schemes is highlighted via representative numerical examples. Regarding the computational cost required for analyzing large arrays of cylinders via MAS and MAS-FMM, substantial savings are indeed achieved when applying the latter. Furthermore, the procedure for proper selection of the basic parameters affecting the performance of the hybrid scheme is also discussed. Finally, a few concluding remarks are outlined, and possible extensions of the work conducted so far are pointed out.

2. Problem Description

As already stated, the problem at hand shown in Figure 1 involves an arbitrarily large array of parallel circular cylinders, which are excited by either a plane wave or an infinite line source parallel to their axes. The cylinders are located at the nodes of a $N_x\times N_y$ lattice with arbitrary coordinates given by ${\overrightarrow{r}}_{\mathrm{cyl}}^{\left(n_x,n_y\right)}=\left(x_{\mathrm{cyl}}^{\left(n_x,n_y\right)},y_{\mathrm{cyl}}^{\left(n_x,n_y\right)}\right)$, for $n_x=1,2,\dots ,N_x$ and $n_y=1,2,\dots ,N_y$. Obviously, the number of cylinders is $N_x{N}_y$, each of which is uniquely identified by the pair $\left(n_x,n_y\right)$. The radius of each cylinder is denoted by $a_{\mathrm{cyl}}^{\left(n_x,n_y\right)}$. Non-uniform and non-orthogonal lattices of cylinders can be treated in a similar manner and with only minor modifications, as discussed below. The incident electromagnetic field, either plane wave or that radiated by the infinite line source discussed above, is assumed to be a continuous wave of frequency $f$ and wavelength $\lambda =c/f$, where $c$ is the speed of light. A time convention $\exp \left(j2\pi ft\right)$ is assumed and suppressed throughout this work.

The aforementioned cylinders can be modeled either as PEC or as lossless/lossy dielectric, which are closer to real-world tree trunks. The former case has been studied in [18]. The cases of lossless and lossy dielectric cylinders are studied below. For brevity, only the case of transverse magnetic TM_z incidence is presented in this work. Given the polarization of the incident field, the auxiliary sources are selected to be infinite line sources of impulsive electric currents, which are aligned parallel to the axes of the cylinders under study. Therefore, the field components $E_z$, $H_x$, and $H_y$ are the only non-zero ones considered hereinafter [19]. The pertinent formulation can be readily modified for the case of transverse electric TE_z incidence, in which the non-zero field components are $E_x$, $E_y$, and $H_z$ [20].

The array of cylinders considered so far seems rather suitable for predicting electromagnetic fields inside and close to orchards. With some small compromise, propagation through or inside forests can also be predicted with ease. This is further discussed in the concluding section of this paper.

3. Formulation

3.1. Standard MAS Formulation

According to standard MAS, a finite number of auxiliary sources are placed on fictitious (auxiliary) circular curves, both inside and outside the boundaries of the cylinders comprising the array under study, for the description of the unknown EM fields in the regions outside and inside these boundaries, respectively. The locations of the auxiliary curves are selected to ensure convergent solutions (e.g., see [4,5,6,7,8,9,10,19,20,21,22,23]). The convergence in the strict sense can be investigated analytically only for certain canonical geometries; nevertheless, when analytical investigations are not possible, the convergence is most often examined numerically [i.e., by carefully checking the pertinent boundary condition(s) errors and the numerical stability of the solutions]. Given the selected locations of the auxiliary sources and the matching points, the currents of the auxiliary sources are solved for in order to satisfy the boundary conditions of the electric and/or the magnetic fields in a point-matching fashion.

For the description of the unknown scattered field in the air between cylinders and outside the lattice of cylinders, $N^{\left(n_x,n_y\right)}$ auxiliary sources are uniformly distributed at a circle of radius $a_{\mathrm{in}}^{\left(n_x,n_y\right)}<a_{\mathrm{cyl}}^{\left(n_x,n_y\right)}$ inside each cylinder. Similarly, for the description of the unknown field induced inside each cylinder, $N^{\left(n_x,n_y\right)}$ auxiliary sources are uniformly distributed at a circle of radius $a_{\mathrm{out}}^{\left(n_x,n_y\right)}>a_{\mathrm{cyl}}^{\left(n_x,n_y\right)}$ outside the cylinder. The auxiliary sources are placed at

$$\begin{gathered} \begin{array}{c}{\overrightarrow{r}}_{\mathrm{in},n}^{\left(n_x,n_y\right)}=\\ \left(x_{\mathrm{cyl}}^{\left(n_x,n_y\right)}+a_{\mathrm{in}}^{\left(n_x,n_y\right)}\cos \left(2\pi \frac{n-1}{N^{\left(n_x,n_y\right)}}\right),y_{\mathrm{cyl}}^{\left(n_x,n_y\right)}+a_{\mathrm{in}}^{\left(n_x,n_y\right)}\sin \left(2\pi \frac{n-1}{N^{\left(n_x,n_y\right)}}\right)\right)\end{array} \\ n=1,2,\dots ,N^{\left(n_x,n_y\right)} \end{gathered}$$

(1)

$$\begin{gathered} \begin{array}{c}{\overrightarrow{r}}_{\mathrm{out},n}^{\left(n_x,n_y\right)}=\\ \left(x_{\mathrm{cyl}}^{\left(n_x,n_y\right)}+a_{\mathrm{out}}^{\left(n_x,n_y\right)}\cos \left(2\pi \frac{n-1}{N^{\left(n_x,n_y\right)}}\right),y_{\mathrm{cyl}}^{\left(n_x,n_y\right)}+a_{\mathrm{out}}^{\left(n_x,n_y\right)}\sin \left(2\pi \frac{n-1}{N^{\left(n_x,n_y\right)}}\right)\right)\end{array} \\ n=1,2,\dots ,N^{\left(n_x,n_y\right)} \end{gathered}$$

(2)

The total number of auxiliary sources for all cylinders is given by

$$N_{\mathrm{T}}=2{\displaystyle \sum _{n_x=1}^{N_x}{\displaystyle \sum _{n_y=1}^{N_y}{N}^{\left(n_x,n_y\right)}}}$$

(3)

When the cylinders are identical (i.e., of equal radius $a_{\mathrm{cyl}}$), one can simply assume $N^{\left(n_x,n_y\right)}=N$ and (3) yields $N_{\mathrm{T}}=2N_x{N}_{y}N$.

Then, the electric field in the air (i.e., outside the cylinders) is expressed as superposition of the incident field $E_z^{\mathrm{inc}}$, and the fields generated by the sets of auxiliary sources inside the cylinders as follows:

$$E_z\left(\overrightarrow{r}\right)=E_z^{\mathrm{inc}}\left(\overrightarrow{r}\right)-\frac{k_0{\zeta }_0}{4}{\displaystyle \sum _{n_x=1}^{N_x}{\displaystyle \sum _{n_y=1}^{N_y}{\displaystyle \sum _{n=1}^{N^{\left(n_x,n_y\right)}}{w}_{\mathrm{in},n}^{\left(n_x,n_y\right)}{H}_0^{\left(2\right)}\left(k_0\left|\overrightarrow{r}-{\overrightarrow{r}}_{\mathrm{in},n}^{\left(n_x,n_y\right)}\right|\right)}}}$$

(4)

where $w_{\mathrm{in},n}^{\left(n_x,n_y\right)}$ are the unknown (complex) currents pertaining to the inner set of auxiliary sources associated with the cylinder denoted by $\left(n_x,n_y\right)$ and $\left|\overrightarrow{r}-{\overrightarrow{r}}_{\mathrm{in},n}^{\left(n_x,n_y\right)}\right|$ are the distances between their positions and the arbitrary observation point $\overrightarrow{r}=\left(x,y\right)$ in the air. The special function $H_{\nu }^{\left(2\right)}$ denotes the Hankel function of the second kind and order $\nu$. In (4), $k_0=2\pi f\sqrt{{\mu }_0{\epsilon }_0}$ (or $k_0=2\pi f/c$) is the wavenumber and ${\zeta }_0=\sqrt{{\mu }_0/{\epsilon }_0}$ is the intrinsic impedance of vacuum. The respective magnetic field is given by

$$\begin{array}{c}{H}_x\left(\overrightarrow{r}\right)=H_x^{\mathrm{inc}}\left(\overrightarrow{r}\right)\\ +\frac{j{k}_0}{4}{\displaystyle \sum _{n_x=1}^{N_x}{\displaystyle \sum _{n_y=1}^{N_y}{\displaystyle \sum _{n=1}^{N^{\left(n_x,n_y\right)}}{w}_{\mathrm{in},n}^{\left(n_x,n_y\right)}\frac{\left(\overrightarrow{r}-{\overrightarrow{r}}_{\mathrm{in},n}^{\left(n_x,n_y\right)}\right)\cdot {\overrightarrow{e}}_y}{\left|\overrightarrow{r}-{\overrightarrow{r}}_{\mathrm{in},n}^{\left(n_x,n_y\right)}\right|}{H}_1^{\left(2\right)}\left(k_0\left|\overrightarrow{r}-{\overrightarrow{r}}_{\mathrm{in},n}^{\left(n_x,n_y\right)}\right|\right)}}}\end{array}$$

(5)

$$\begin{array}{c}{H}_y\left(\overrightarrow{r}\right)=H_y^{\mathrm{inc}}\left(\overrightarrow{r}\right)\\ -\frac{j{k}_0}{4}{\displaystyle \sum _{n_x=1}^{N_x}{\displaystyle \sum _{n_y=1}^{N_y}{\displaystyle \sum _{n=1}^{N^{\left(n_x,n_y\right)}}{w}_{\mathrm{in},n}^{\left(n_x,n_y\right)}\frac{\left(\overrightarrow{r}-{\overrightarrow{r}}_{\mathrm{in},n}^{\left(n_x,n_y\right)}\right)\cdot {\overrightarrow{e}}_x}{\left|\overrightarrow{r}-{\overrightarrow{r}}_{\mathrm{in},n}^{\left(n_x,n_y\right)}\right|}{H}_1^{\left(2\right)}\left(k_0\left|\overrightarrow{r}-{\overrightarrow{r}}_{\mathrm{in},n}^{\left(n_x,n_y\right)}\right|\right)}}}\end{array}$$

(6)

where ${\overrightarrow{e}}_x$ and ${\overrightarrow{e}}_y$ are the unit vectors in the $x$ and $y$ directions, respectively. The dot between vectors denotes the inner product. Similarly, the fields inside each cylinder are expressed as superimpositions of the fields generated by the respective set of auxiliary sources outside the cylinder

$$E_z^{\left(n_x,n_y\right)}\left(\overrightarrow{r}\right)=-\frac{k\zeta }{4}{\displaystyle \sum _{n=1}^{N^{\left(n_x,n_y\right)}}{w}_{\mathrm{out},n}^{\left(n_x,n_y\right)}{H}_0^{\left(2\right)}\left(k\left|\overrightarrow{r}-{\overrightarrow{r}}_{\mathrm{out},n}^{\left(n_x,n_y\right)}\right|\right)}$$

(7)

$$H_x^{\left(n_x,n_y\right)}\left(\overrightarrow{r}\right)=\frac{jk}{4}{\displaystyle \sum _{n=1}^{N^{\left(n_x,n_y\right)}}{w}_{\mathrm{out},n}^{\left(n_x,n_y\right)}\frac{\left(\overrightarrow{r}-{\overrightarrow{r}}_{\mathrm{out},n}^{\left(n_x,n_y\right)}\right)\cdot {\overrightarrow{e}}_y}{\left|\overrightarrow{r}-{\overrightarrow{r}}_{\mathrm{out},n}^{\left(n_x,n_y\right)}\right|}{H}_1^{\left(2\right)}\left(k\left|\overrightarrow{r}-{\overrightarrow{r}}_{\mathrm{out},n}^{\left(n_x,n_y\right)}\right|\right)}$$

(8)

$$\begin{array}{c}{H}_y^{\left(n_x,n_y\right)}\left(\overrightarrow{r}\right)=\\ -\frac{jk}{4}{\displaystyle \sum _{n=1}^{N^{\left(n_x,n_y\right)}}{w}_{\mathrm{out},n}^{\left(n_x,n_y\right)}\frac{\left(\overrightarrow{r}-{\overrightarrow{r}}_{\mathrm{out},n}^{\left(n_x,n_y\right)}\right)\cdot {\overrightarrow{e}}_x}{\left|\overrightarrow{r}-{\overrightarrow{r}}_{\mathrm{out},n}^{\left(n_x,n_y\right)}\right|}{H}_1^{\left(2\right)}\left(k\left|\overrightarrow{r}-{\overrightarrow{r}}_{\mathrm{out},n}^{\left(n_x,n_y\right)}\right|\right)}\end{array}$$

(9)

where $w_{\mathrm{out},n}^{\left(n_x,n_y\right)}$ are the unknown (complex) currents pertaining to the outer set of auxiliary sources associated with the cylinder denoted by $\left(n_x,n_y\right)$ and $\left|\overrightarrow{r}-{\overrightarrow{r}}_{\mathrm{out},n}^{\left(n_x,n_y\right)}\right|$ are the distances between their positions and the arbitrary observation point $\overrightarrow{r}=\left(x,y\right)$ inside the said cylinder. In (7)–(9), $k=2\pi f\sqrt{{\mu }_0\epsilon }$ is the wavenumber and $\zeta =\sqrt{{\mu }_0/\epsilon }$ is the intrinsic impedance associated with the dielectric medium of the cylinders. Lossy homogeneous media with relative permittivity ${\epsilon }_r$ and conductivity $\sigma$ have complex dielectric constant $\epsilon =\left[{\epsilon }_r-j\sigma /\left(2\pi f{\mu }_0\right)\right]{\epsilon }_0$. When electric losses are absent ($\sigma =0$), the last expression simplifies to $\epsilon ={\epsilon }_r{\epsilon }_0$. In the case of PEC cylinders, the fields inside them vanish and the analysis is simplified to that presented in [18]. The formulation and results of [18] have been regenerated using the generalized formulation and the code developed in the present paper. However, these results are not included here to avoid unnecessary overlapping between the two papers.

Next, the unknown currents $w_{\mathrm{in},n}^{\left(n_x,n_y\right)}$ and $w_{\mathrm{out},n}^{\left(n_x,n_y\right)}$ are obtained by enforcing the boundary conditions of the tangential electric and the tangential magnetic field at $N^{\left(m_x,m_y\right)}$ discrete matching points on each cylinder, hereby denoted by $\left(m_x,m_y\right)$ with $m_x=1,2,\dots ,N_x$ and $m_y=1,2,\dots ,N_y$ (so as to distinguish between counters pertaining to auxiliary sources and matching points), which are placed at

$$\begin{gathered} \begin{array}{c}{\overrightarrow{r}}_{\mathrm{cyl},m}^{\left(m_x,m_y\right)}=\\ \left(x_{\mathrm{cyl}}^{\left(m_x,m_y\right)}+a_{\mathrm{cyl}}^{\left(m_x,m_y\right)}\cos \left(2\pi \frac{m-1}{N^{\left(m_x,m_y\right)}}\right),y_{\mathrm{cyl}}^{\left(m_x,m_y\right)}+a_{\mathrm{cyl}}^{\left(m_x,m_y\right)}\sin \left(2\pi \frac{m-1}{N^{\left(m_x,m_y\right)}}\right)\right)\end{array} \\ m=1,2,\dots ,N^{\left(m_x,m_y\right)} \end{gathered}$$

(10)

Standard boundary conditions imply that the tangential components of the electric and magnetic fields must be continuous across the boundaries of all cylinders. Thus, the following equations are readily obtained:

$$\begin{gathered} \begin{array}{c}-\frac{k_0{\zeta }_0}{4}{\displaystyle \sum _{n_x=1}^{N_x}{\displaystyle \sum _{n_y=1}^{N_y}{\displaystyle \sum _{n=1}^{N^{\left(n_x,n_y\right)}}{w}_{\mathrm{in},n}^{\left(n_x,n_y\right)}{H}_0^{\left(2\right)}\left(k_0\left|{\overrightarrow{r}}_{\mathrm{cyl},m}^{\left(m_x,m_y\right)}-{\overrightarrow{r}}_{\mathrm{in},n}^{\left(n_x,n_y\right)}\right|\right)}}}\\ +\frac{k\zeta }{4}{\displaystyle \sum _{n=1}^{N^{\left(n_x,n_y\right)}}{w}_{\mathrm{out},n}^{\left(n_x,n_y\right)}{H}_0^{\left(2\right)}\left(k\left|{\overrightarrow{r}}_{\mathrm{cyl},m}^{\left(m_x,m_y\right)}-{\overrightarrow{r}}_{\mathrm{out},n}^{\left(n_x,n_y\right)}\right|\right)}=-E_z^{\mathrm{inc}}\left({\overrightarrow{r}}_{\mathrm{cyl},m}^{\left(m_x,m_y\right)}\right)\end{array} \\ {m}_x=1,2,\dots ,N_x, m_y=1,2,\dots ,N_y \end{gathered}$$

(11)

$$\begin{gathered} \begin{array}{c}-\frac{j{k}_0}{4}{\displaystyle \sum _{n_x=1}^{N_x}{\displaystyle \sum _{n_y=1}^{N_y}{\displaystyle \sum _{n=1}^{N^{\left(n_x,n_y\right)}}{w}_{\mathrm{in},n}^{\left(n_x,n_y\right)}\frac{\left({\overrightarrow{r}}_{\mathrm{cyl},m}^{\left(m_x,m_y\right)}-{\overrightarrow{r}}_{\mathrm{cyl}}^{\left(m_x,m_y\right)}\right)\cdot \left({\overrightarrow{r}}_{\mathrm{cyl},m}^{\left(m_x,m_y\right)}-{\overrightarrow{r}}_{\mathrm{in},n}^{\left(n_x,n_y\right)}\right)}{a_{\mathrm{cyl}}^{\left(m_x,m_y\right)}\left|{\overrightarrow{r}}_{\mathrm{cyl},m}^{\left(m_x,m_y\right)}-{\overrightarrow{r}}_{\mathrm{in},n}^{\left(n_x,n_y\right)}\right|}}}}\\ \times H_1^{\left(2\right)}\left(k_0\left|{\overrightarrow{r}}_{\mathrm{cyl},m}^{\left(m_x,m_y\right)}-{\overrightarrow{r}}_{\mathrm{in},n}^{\left(n_x,n_y\right)}\right|\right)\\ +\frac{jk}{4}{\displaystyle \sum _{n=1}^{N^{\left(n_x,n_y\right)}}{w}_{\mathrm{out},n}^{\left(n_x,n_y\right)}\frac{\left({\overrightarrow{r}}_{\mathrm{cyl},m}^{\left(m_x,m_y\right)}-{\overrightarrow{r}}_{\mathrm{cyl}}^{\left(m_x,m_y\right)}\right)\cdot \left({\overrightarrow{r}}_{\mathrm{cyl},m}^{\left(m_x,m_y\right)}-{\overrightarrow{r}}_{\mathrm{out},n}^{\left(n_x,n_y\right)}\right)}{a_{\mathrm{cyl}}^{\left(m_x,m_y\right)}\left|{\overrightarrow{r}}_{\mathrm{cyl},m}^{\left(m_x,m_y\right)}-{\overrightarrow{r}}_{\mathrm{out},n}^{\left(n_x,n_y\right)}\right|}}\\ \times H_1^{\left(2\right)}\left(k\left|{\overrightarrow{r}}_{\mathrm{cyl},m}^{\left(m_x,m_y\right)}-{\overrightarrow{r}}_{\mathrm{out},n}^{\left(n_x,n_y\right)}\right|\right)=-H_t^{\mathrm{inc}}\left({\overrightarrow{r}}_{\mathrm{cyl},m}^{\left(m_x,m_y\right)}\right)\end{array} \\ {m}_x=1,2,\dots ,N_x, m_y=1,2,\dots ,N_y \end{gathered}$$

(12)

The Equations of (11) and (12) constitute a linear system of $N_{\mathrm{T}}$ equations with $N_{\mathrm{T}}$ unknowns (i.e., the aforesaid currents) of the form $Z\cdot w=-b$, where $Z$ is the interaction matrix, $w$ is the vector containing the unknown currents, and $b$ is the excitation vector. For convenience, these can be written in the form of sub-matrices similar to those provided in [18]. From the solution of this system, one can directly compute the electric and magnetic fields at any observation point using the expressions of (4)–(9).

As already discussed, non-uniform arrays of cylinders can be treated using the method at hand by zeroing the contributions pertaining to certain pairs $\left(n_x,n_y\right)$. This approach is rather simple and straightforward and is omitted for brevity. It is further stressed that non-orthogonal cylinder lattices are also within the capabilities of both schemes, as long as the centers of the cylinders are arbitrary, with the proviso that the cylinders do not overlap. Therefore, the MAS/MAS-FMM schemes of this work are essentially capable of coping with uniform/non-uniform and orthogonal/non-orthogonal cylinder arrays.

After solving the system for the currents $w_{\mathrm{in},n}^{\left(n_x,n_y\right)}$ and $w_{\mathrm{out},n}^{\left(n_x,n_y\right)}$, the boundary conditions utilized for solving the problem are satisfied exactly only at the collocation points. As an indicative error metric, one can adopt the (normalized) boundary-condition error, which is given by the ratio of the tangential electric or magnetic field difference across the boundary over the corresponding incident field. Obviously, this error is expected to be zero (or very close to zero due to the finite accuracy of computers) at the collocation points and non-zero between them, reaching local maxima at or very close to the respective mid-points. The occurrence of fairly small error peaks and means (small compared to what may be expected empirically), which decrease steadily as $N_{\mathrm{T}}$ is increased, is a strong indication that the numerical solutions are numerically stable and trustworthy.

In order to validate the MAS model discussed so far, extensive numerical tests were conducted and evaluated. In all cases in which the auxiliary sources were placed according to the findings of [4,5,6,7,8,9,10,18], the errors in the boundary conditions were indeed found to diminish steadily with increasing $N_{\mathrm{T}}$. Moreover, the solutions obtained were compared to those of other methods in order to independently check the validity of the results. Exemplary cases are presented below. At this point, it is noted that $N^{\left(n_x,n_y\right)}$ should be selected so that the distance between neighboring sources belonging to each auxiliary curve is comparable to the respective distance from the boundary, which is ${\delta }_{\mathrm{in}}^{\left(n_x,n_y\right)}=a_{\mathrm{cyl}}^{\left(n_x,n_y\right)}-a_{\mathrm{in}}^{\left(n_x,n_y\right)}$ for the inner set and ${\delta }_{\mathrm{out}}^{\left(n_x,n_y\right)}=a_{\mathrm{out}}^{\left(n_x,n_y\right)}-a_{\mathrm{cyl}}^{\left(n_x,n_y\right)}$ for the outer one. Therefore, as a quite safe choice for good conditioning of the system and overall behavior of the numerical solutions, one can start with $N^{\left(n_x,n_y\right)}$ close to $2\pi a_{\mathrm{in}}^{\left(n_x,n_y\right)}/{\delta }_{\mathrm{in}}^{\left(n_x,n_y\right)}$ or $2\pi a_{\mathrm{out}}^{\left(n_x,n_y\right)}/{\delta }_{\mathrm{out}}^{\left(n_x,n_y\right)}$, and progressively increase them until some error/stability criterion is met. One should also bear in mind that severe ill-conditioning problems typically occur when $N^{\left(n_x,n_y\right)}$ becomes quite larger than $2\pi a_{\mathrm{in}}^{\left(n_x,n_y\right)}/{\delta }_{\mathrm{in}}^{\left(n_x,n_y\right)}$ or $2\pi a_{\mathrm{out}}^{\left(n_x,n_y\right)}/{\delta }_{\mathrm{out}}^{\left(n_x,n_y\right)}$. These issues can be identified and isolated using higher-precision arithmetic, various system solvers, and alternative routine implementations.

3.2. Hybrid MAS-FMM Formulation

The standard MAS model discussed so far is quite simple and capable of treating arbitrarily large arrays of cylinders. However, as $N_{\mathrm{T}}$ gets larger, the computational cost of the matrix formation ensuing from the system of (11) and (12) becomes excessively high. As a first improvement, one can simply utilize the large-argument approximations $H_0^{\left(2\right)}\left(u\right)~\sqrt{2/\left(\pi u\right)}\exp \left(-ju+j\pi /4\right)$ and $H_1^{\left(2\right)}\left(u\right)~\sqrt{2/\left(\pi u\right)}\exp \left(-ju+j3\pi /4\right)$ as $u\to \infty$, which can be used for computing the interactions involving internal auxiliary sources and matching points pertaining to different cylinders (i.e., for $n_x\ne m_x$ and/or $n_y\ne m_y$). Note that the auxiliary sources located outside each cylinder do not contribute to interactions between cylinders and, therefore, their fields should be approximated using the aforesaid large-argument expressions only in cases of electrically large cylinders. Roughly speaking, these simple approximations alone can save up to about 70% of the execution time required for computing $Z$ [18]. Nevertheless, as discussed below, the computational cost for filling $Z$ can be further reduced.

Although often less costly than MoM, the complexity of MAS is still of order $O\left(N_{\mathrm{T}}^2\right)$ for filling $Z$ and $O\left(N_{\mathrm{T}}^3\right)$ for solving $Z\cdot w=-b$ via conventional methods (e.g., LU decomposition), meaning that the relative advantages of MAS, at least with respect to the execution time, may become less pronounced for very large $N_{\mathrm{T}}$. To reduce the complexity and the computational cost, the interactions between fairly distant auxiliary sources and matching points can be properly grouped together into clusters as in [18], which is exactly the fundamental idea behind the FMM. Direct implementations of FMM reduce the complexity for obtaining $Z\cdot w=-b$ from $O\left(N_{\mathrm{T}}^2\right)$ to $O\left(N_{\mathrm{T}}^{3/2}\right)$, while more sophisticated implementations (e.g., recursive clustering schemes and preconditioning) can lead to even lower complexity down to $O\left(N_{\mathrm{T}}\log N_{\mathrm{T}}\right)$ [11].

The FMM essentially consists in the decomposition $Z=Z^{\prime }+Z^{″}$, where $Z^{\prime }$ contains the near-field part of $Z$ and $Z^{″}$ contains the respective far-field part. The computation of the latter is accelerated via the decomposition of the form $Z^{″}=V^{\mathrm{T}}\cdot T\cdot \tilde{V}$ (the superscript ‘T’ denotes the transpose matrix/vector), in which $V$, $T$, and $\tilde{V}$ are all sparse matrices. The cluster grouping is accomplished via $V$, whereas $\tilde{V}$ performs the necessary disaggregation. The translation matrix $T$ essentially contains the interactions between clusters. The aforesaid decomposition is based on the addition theorem for Hankel functions (e.g., see Section 10.23 in [24])

$$H_{\nu }^{\left(2\right)}\left(k_0\left|{\overrightarrow{r}}_0+\overrightarrow{\delta }\right|\right)\exp \left(j\nu \left({\chi }_r-{\chi }_0\right)\right)={\displaystyle \sum _{l=-\infty }^{\infty }{H}_{l+\nu }^{\left(2\right)}\left(k_0\left|{\overrightarrow{r}}_0\right|\right)J_l\left(k_0\left|\overrightarrow{\delta }\right|\right)\cos \left(l\eta \right)}$$

(13)

where $\eta$ is the angle between ${\overrightarrow{r}}_0$ and $\overrightarrow{\delta }$, whereas ${\chi }_r$ and ${\chi }_0$ are the angles of the vectors ${\overrightarrow{r}}_0+\overrightarrow{\delta }$ and ${\overrightarrow{r}}_0$ measured from the $x$ axis, respectively. Here, ${\overrightarrow{r}}_0$ connects the center of a cluster of auxiliary sources with the center of a cluster of matching points, while $\overrightarrow{\delta }={\overrightarrow{\delta }}_{\mathrm{MP}}-{\overrightarrow{\delta }}_{\mathrm{AS}}$ so that the vector ${\overrightarrow{r}}_0+\overrightarrow{\delta }$ connects the location of any auxiliary source with any matching point. In (13), $H_{l+\nu }^{\left(2\right)}(\cdot )$ and $J_l(\cdot )$ denote Hankel (of the second kind) and Bessel functions of order $l$ respectively. As long as $J_l$ practically vanishes for $l>k_0\left|\overrightarrow{\delta }\right|$, the summation in (13) can be truncated so as to include only the terms with $\left|l\right|\le L$, where $L~k_0\max \left\{\left|\overrightarrow{\delta }\right|\right\}$. This parameter can be readily estimated from the radii of the cylinders yielding the rough upper bound $k_0\max \left\{\left|\overrightarrow{\delta }\right|\right\}<2k_0\max \left\{a_{\mathrm{cyl}}^{\left(n_x,n_y\right)}\right\}$. Alternatively, it could be heuristically selected via trial and error. Furthermore, one can proceed using

$$J_l\left(u\right)\cos \left(l\eta \right)=\frac{j^{-l}}{2\pi }{\displaystyle {\int }_0^{2\pi }\exp \left(ju\cos \xi \right)\cos \left(l\left(\xi -\eta \right)\right)d\xi }$$

(14)

which yields the following integral:

$$J_l\left(k_0\left|\overrightarrow{\delta }\right|\right)\cos \left(l\eta \right)=\frac{j^{-l}}{2\pi }{\displaystyle {\int }_0^{2\pi }\exp \left(j{k}_0\stackrel{⌢}{\xi }\cdot \overrightarrow{\delta }\right)\cos \left(l\gamma \right)d\xi }$$

(15)

where $\stackrel{⌢}{\xi }$ is the unit vector pointing towards the direction of integration, which is normal to the unit circle, and $\gamma$ is the angle between ${\overrightarrow{r}}_0$ and $\stackrel{⌢}{\xi }$. The integral in (15) can be computed numerically over a finite set of $P$ uniformly distributed angles on the unit circle. Then, using (13)–(15), one can write

$$\begin{array}{c}{H}_{\nu }^{\left(2\right)}\left(k_0\left|{\overrightarrow{r}}_0+\overrightarrow{\delta }\right|\right)\exp \left(j\nu \left({\chi }_r-{\chi }_0\right)\right)\approx \\ {\displaystyle {\int }_0^{2\pi }\exp \left(j{k}_0\stackrel{⌢}{\xi }\cdot \left({\overrightarrow{\delta }}_{\mathrm{MP}}-{\overrightarrow{\delta }}_{\mathrm{AS}}\right)\right){\displaystyle \sum _{l=-L}^L\frac{j^{-l}}{2\pi }{H}_{l+\nu }^{\left(2\right)}\left(k_0\left|{\overrightarrow{r}}_0\right|\right)\cos \left(l\gamma \right)}d\xi }\end{array}$$

(16)

The integral in (16) can be estimated numerically by the summation

$$\begin{array}{c}{H}_{\nu }^{\left(2\right)}\left(k_0\left|{\overrightarrow{r}}_0+\overrightarrow{\delta }\right|\right)\approx \\ \exp \left(-j\nu \left({\chi }_r-{\chi }_0\right)\right){\displaystyle \sum _{p=1}^P{V}_{\mathrm{AS}}\left({\stackrel{⌢}{\xi }}_p\right)T\left({\stackrel{⌢}{\xi }}_p\right){\tilde{V}}_{\mathrm{MP}}\left({\stackrel{⌢}{\xi }}_p\right)\Delta \xi }\end{array}$$

(17)

where $\Delta \xi =2\pi /P$,$V_{\mathrm{AS}}\left({\stackrel{⌢}{\xi }}_p\right)=\exp \left(-j{k}_0{\stackrel{⌢}{\xi }}_p\cdot {\overrightarrow{\delta }}_{\mathrm{AS}}\right)$ and ${\tilde{V}}_{\mathrm{MP}}\left({\stackrel{⌢}{\xi }}_p\right)=\exp \left(j{k}_0{\stackrel{⌢}{\xi }}_p\cdot {\overrightarrow{\delta }}_{\mathrm{MP}}\right)$ are the elements of the one-dimensional matrices $V$ and $\tilde{V}$, respectively, whereas $T\left({\stackrel{⌢}{\xi }}_p\right)$ are the elements of the matrix $T$ obtained from

$$T\left({\stackrel{⌢}{\xi }}_p\right)={\displaystyle \sum _{l=-L}^L\frac{j^{-l}}{2\pi }{H}_{l+\nu }^{\left(2\right)}\left(k_0\left|{\overrightarrow{r}}_0\right|\right)\cos \left(l{\gamma }_p\right)}$$

(18)

Note that (18) can be used for the computation of both electric-field ($\nu =0$) and magnetic-field ($\nu =1$) interactions. When the distances between cylinders are significantly larger than their radii, (18) can be simplified to

$$T\left({\stackrel{⌢}{\xi }}_p\right)\approx {\displaystyle \sum _{l=-L}^L\frac{j^{-l}}{2\pi }\sqrt{\frac{2}{\pi k_0\left|{\overrightarrow{r}}_0\right|}}\exp \left(-j{k}_0\left|{\overrightarrow{r}}_0\right|+j\frac{\left(2l+2\nu +1\right)\pi }{4}\right)\cos \left(l{\gamma }_p\right)}$$

(19)

This approximation further accelerates the procedure for filling the matrix $T$ and can achieve additional cost reductions. The number of terms in the summation of (17) should be selected in the order of $k_0\max \left\{\left|\overrightarrow{\delta }\right|\right\}$. In practice, one can simply start with $P=L$ and increase (independently) $P$ and/or $L$ so that the results of interest are found to be numerically stable and trustworthy—this point is further discussed below. When the distances between cylinders are sufficiently larger than their radii—note that this condition holds for tree trunks in orchards/forests and for many other practical purposes—the clusters coincide with the sets of auxiliary sources and matching points associated with each cylinder. As a result, $Z^{\prime }$ contains the sub-matrices corresponding to $n_x=m_x$ and $n_y=m_y$, with its other entries set to zero. With regard to the non-diagonal entries of $Z^{″}$, these are computed with the aid of (17) and (18), whereas its diagonal ones can be set to zero. This sort of matrix decomposition yields sparse matrices and facilitates the efficient storage and the fast inversion for the linear system at hand, leading to significant savings with respect to the overall computational cost.

4. Numerical Results and Discussion

Representative numerical results are provided here for the validation of the numerical schemes discussed so far and, also, for their comparison from the viewpoints of accuracy, complexity, and computational cost. Though the formulation presented in this work was written down explicitly and thoroughly tested for plane-wave and line-source excitations, the results that follow regard only the former case, for the sake of simplicity and for making direct comparisons with results from the existing literature and others obtained via other methods.

In order to validate the MAS model described above, exhaustive numerical tests were conducted. In all cases, at least when the auxiliary sources were properly placed as described above, the errors in the boundary conditions of the tangential electric/magnetic fields were found to be fairly small (see below) and, also, to diminish steadily as the numbers of auxiliary sources (and matching points) was increased. Moreover, the numerical solutions were compared to those of analytical [25] and numerical methods [26,27,28,29,30,31,32], in order to independently check their validity. A few exemplary cases are presented below.

First, the simple case of an isolated dielectric cylinder (comprising a trivial 1×1 lattice) was examined. The parameters were $\epsilon =5{\epsilon }_0$, $a_{\mathrm{cyl}}=0.2$ m, $a_{\mathrm{in}}=0.75a_{\mathrm{cyl}}$, and $a_{\mathrm{out}}=1.25a_{\mathrm{cyl}}$. Given that only one cylinder was considered in this case, the MAS-FMM hybrid was not applied, and the tests were run for MAS in its original, naive form. The incident field was a plane wave with frequency $f=1$ GHz, impinging from the left (i.e., coming from the direction ${\varphi }_0={180}^{\circ }$). Several numerical tests have shown that the numerical stability and the convergence behavior of the MAS solutions were just as anticipated from [4,5,6,7,8,9,10]. Moreover, the computed near and far fields were found to be virtually identical to those of the exact solution [25] (for sufficiently large $N$). The relative difference between the electric-field magnitude as obtained from MAS with $N=40$ and the respective analytical solution are depicted in Figure 2, as a function of the observation angle $\varphi$ at a distance $r=10a_{\mathrm{cyl}}$ from the center of the cylinder at hand. Apparently, the relative differences are symmetric and small, as they do not exceed 0.07%. The symmetric behavior of the difference between numerical and exact solution is a strong indication that the computations have not been corrupted by any significant numerical noise and roundoff errors. Furthermore, in order to examine the credibility of the developed code for coping with lossy dielectrics, the same problem was solved after changing the dielectric constant of the cylinder to the complex-valued one $\epsilon =\left(5-j3\right){\epsilon }_0$. Again, the numerical results were found to be virtually identical to those of the respective analytical solution, as also shown in Figure 2. In this latter case, the relative differences are lower than 0.01%. The respective peak (max) and average (mean) errors associated with the continuity of the electric and magnetic fields are reported in

Table 1. Peak (max) and mean boundary-condition errors for various exemplary cases of lattices of dielectric cylinders.

Case	Parameters	Max Errors (%)		Mean Errors (%)
		E-Field	H-Field	E-Field	H-Field
1 × 1 lattice $N_x=1$ $N_y=1$	$f=1$ GHz $\epsilon =5{\epsilon }_0$ $a_{\mathrm{cyl}}=0.2$ m $a_{\mathrm{in}}=0.75a_{\mathrm{cyl}}$ $a_{\mathrm{out}}=1.25a_{\mathrm{cyl}}$ $N=40$	0.023	0.154	0.007	0.059
1 × 1 lattice $N_x=1$ $N_y=1$	$f=1$ GHz $\epsilon =\left(5-j3\right){\epsilon }_0$ $a_{\mathrm{cyl}}=0.2$ m $a_{\mathrm{in}}=0.75a_{\mathrm{cyl}}$ $a_{\mathrm{out}}=1.25a_{\mathrm{cyl}}$ $N=40$	0.029	0.253	0.010	0.087
1 × 5 lattice $N_x=1$ $N_y=5$	$f=0.3$ GHz $\epsilon =5{\epsilon }_0$ $a_{\mathrm{cyl}}=0.1\lambda$ $a_{\mathrm{in}}=0.75a_{\mathrm{cyl}}$ $a_{\mathrm{out}}=1.25a_{\mathrm{cyl}}$ $d=0.75\lambda$ $N=40$	0.023	0.152	0.007	0.057
5 × 5 lattice $N_x=5$ $N_y=5$	$f=1$ GHz $\epsilon =5{\epsilon }_0$ $a_{\mathrm{cyl}}=0.2$ m $a_{\mathrm{in}}=0.75a_{\mathrm{cyl}}$ $a_{\mathrm{out}}=1.25a_{\mathrm{cyl}}$ $d=2$ m $N=40$	0.027	0.189	0.005	0.039

. These percentage errors correspond to the magnitude of the field difference (across the boundary) normalized to the magnitude of the incident plane wave. To further illustrate the very good agreement between numerical and exact solutions, Figure 3 exhibits the electric-field magnitude as obtained from MAS and the respective exact solutions for the two cases examined so far.

As another example, the 1 × 5 lattice studied in [29] was examined. The parameters pertaining to the five identical cylinders taken from [29] were $\epsilon =5{\epsilon }_0$, $a_{\mathrm{cyl}}=0.1\lambda$, $a_{\mathrm{in}}=0.75a_{\mathrm{cyl}}$, and $a_{\mathrm{out}}=1.25a_{\mathrm{cyl}}$. The distance between adjacent cylinder centers was set to $d=0.75\lambda$ and the incident field was a plane wave with frequency $f=0.3$ GHz, impinging from the direction ${\varphi }_0={180}^{\circ }$. Again, the results shown herein were obtained with $N=40$. To justify this choice, the convergence of the numerical solutions is exhibited in Figure 4, which shows the logarithm (with base 10) of the computed peak electric-field and magnetic-field boundary-condition errors for $10\le N\le 100$. Such logarithmic curves are often used for convergence/stability assessments when applying frequency-domain numerical methods (e.g., see [33]). Obviously, the errors decrease rapidly and steadily, as one could expect from the findings and discussions in [4,5,6,7,8,9,10]. It is particularly stressed that no ill-conditioning issues were encountered in the range $10\le N\le 100$ and, therefore, the results of Figure 4 (and Figure 5 below) are believed to be free of numerical noise and roundoff errors. The respective mean boundary-condition errors were found to exhibit a similar behavior and are not shown here for the sake of brevity. The computed radar cross-section or (RCS), which is analogous to the square of the scattered electric-field magnitude in the far-field region, is shown in Figure 5 as a function of the observation angle $\varphi$ at a distance $r=1000\lambda$ from the center of the lattice. For clarity, the RCS is shown in dBm (i.e., in decibels normalized to 1 m). Evidently, the results of MAS/MAS-FMM in Figure 5 are in very close agreement with those obtained from [29], which are shown as rhombic dots. In particular, the forward and backward RCS tabulated in [29] almost coincide with those of MAS/MAS-FMM, as long as they differ by (only) about 1–1.5%. Certain other RCS values estimated from ([29], Figure 5 and Figure 6) are also shown in Figure 5. The results from MAS and MAS-FMM are virtually identical at the scale of the plot. The latter were obtained with the FMM parameters set to $L=P=8$. These values were selected after extensive numerical experimentation regarding the convergence behavior of the series involved in (16)–(19) and the numerical stability of the MAS-FMM solutions. Specifically, for the 5 × 5 lattice case, the upper bound $2k_0\max \left\{a_{\mathrm{cyl}}^{\left(n_x,n_y\right)}\right\}$ of the estimator $L~k_0\max \left\{\left|\overrightarrow{\delta }\right|\right\}$ is (roughly) 8.3, which is consistent with the choice $L=P=8$.

Next, the MAS and MAS-FMM schemes are compared from the aspect of the computational cost (in terms of the execution time and the memory usage). To this end, the total execution time for forming and solving the respective linear system for the MAS currents is shown in Figure 6 as a function of the total number of unknowns. The results depicted in Figure 6 were obtained for various numbers of cylinders comprising lattices of increasing size from 2 × 2 and up to 20 × 20, which required 32,000 auxiliary sources (and matching points) for its adequate modeling (without altering $N$, so as for the error to remain unchanged). The execution times have been measured on the same personal computer and are indicative of the associated computational cost. As it is obvious, the MAS-FMM scheme is indeed notably faster than the naïve MAS scheme, an advantage that becomes even more pronounced as the array becomes larger. From the results of Figure 6, the MAS curve exhibits a polynomial behavior of order 3, whereas the MAS-FMM curve increases with a much lower rate (order 2 or slightly lower). Apparently, the savings in the execution time reach and even exceed 80% for $N_{\mathrm{T}}>\mathrm{25,000}$. Important savings in the memory consumption were also documented during the exhaustive runs conducted and performed. While small lattices (e.g., 2 × 2 and 4 × 4) seemed not to increase the memory usage during the execution of the conventional MAS code, larger lattices appeared to be quite demanding in terms of memory; namely, the memory usage reached 16% for the 8 × 8 lattice, 30% for the 12 × 12, 60% for the 16 × 16, and 95% for the 20 × 20. On the other hand, the memory usage during the execution of the MAS-FMM code did not exceed 12–13% for all these lattices. Besides, numerous runs for various large lattices have shown that these savings in the execution time and memory consumption are absolutely representative and feasible and should be anticipated for lattices with dimensions pertaining to actual orchards or forests and for frequencies in the UHF band.

Finally, the MAS-FMM scheme proposed and examined in this paper is compared to the well-established package COMSOL using FEM [34]. The aim of this comparison is twofold: to validate the results of MAS-FMM in a quite demanding case and to further illustrate the applicability of the MAS-FMM scheme to large-scale problems. At this point, it is worth mentioning that the application of standard, general-purpose FEM schemes requires mesh generation for the whole domain at hand, which may completely consume the available computer resources even for moderately large arrays of cylinders, especially when the distances between their centers are much larger than their radii. On the contrary, MAS-FMM can analyze very large arrays of cylinders, without severe limitations on the size of the array. For the aforesaid comparison, a 5 × 5 lattice of cylinders was tested. The parameters pertaining to the cylinders are contained in Table 1. The computed RCS (again in dBm) is shown in Figure 7, in which certain negative values are out of scale and are not shown. Both methods predict the rather anticipated rapidly oscillating patterns, which are very close to each other. Note, in particular, that the positive peaks agree remarkably well (taking into account the nature of the oscillatory pattern for the examined electrically large lattice). The evident discrepancies in the negative peaks (nulls) are rather unimportant, as long as they correspond to very small RCS values. The distance between cylinder centers in the lattice at hand was selected to be relatively short in order to facilitate convergence of the FEM/COMSOL solution, used as reference. However, this selection inevitably resulted in suboptimal performance of the FMM module in our algorithm, verified by a few discrepancies between MAS-FMM and FEM, as, e.g., in the central lobe (at $\varphi$ = 180°) of the RCS pattern shown in Figure 7. Such discrepancies are not present in analogous comparisons between FEM and pure MAS; therefore, they are clearly attributed to FMM, which is certainly more accurate in larger geometries, not easily computable by FEM. Other factors affecting accuracy will be investigated in a future work on MAS-FMM and its use for near- and far-field computations. Regarding the execution time, COMSOL required about 81 s, whereas MAS-FMM required about 8 s. The differences became more dramatic for the 10 × 10 lattice, for which COMSOL required about 40 min and MAS-FMM only 1 min. It is noted that COMSOL failed to analyze the case of the 20 × 20 lattice (due to lack of memory), which, however, was successfully treated by MAS-FMM (in about 15 min).

5. Conclusions and Prospect

The MAS-FMM scheme proposed and examined in this paper is a hybrid numerical technique combining the general approach of MAS together with the grouping/clustering concept of FMM. The aim of this combination is to achieve computational efficiency without notably compromising numerical accuracy, especially for large-scale problems that involve many scatterers. Many numerical experiments have shown the efficacy of the proposed hybridization and revealed huge savings in the associated computational cost. Though the scheme presented and evaluated here is for 2D problems only, 3D extensions and implementations are possible, and even more promising from the aspect of the potential computational savings.

Regarding the possible extensions of this work, these can be summarized as follows:

• The treatment of random lattices of cylinders for the deterministic or stochastic analysis of electromagnetic propagation through complex vegetation environments (like forests) is quite simple and straightforward. It is particularly stressed that macroscopic, stochastic approaches may be best suited for characterizing complex vegetation environments. For this purpose, one can start from the assumption of proper statistical distributions for the radii of the cylinders and the distances between them (depending on the specific characteristics of the orchard/forest or other environment of interest), proceed using some certain strategy for applying the MAS-FMM scheme of this work repeatedly, and finally compute distributions or moments for the statistical characterization of the propagation environment of interest.
• The generalization of the proposed MAS-FMM scheme to cope with 3D scatterers. Though not easy to manipulate and present in a comprehensive manner, this can be accomplished via lengthy but straightforward modifications [1].
• The systematic assessment of the proposed MAS-FMM scheme from the aspect of the associated complexity and computational cost and the relevant comparisons with other established methods of Computational Electromagnetics. To this end, rigorous cost metrics, such as polynomials involving the total number of unknowns, could be obtained, and subsequently utilized to provide a solid basis for cost comparisons. Furthermore, apart from the cost metrics usually used pertaining to the execution time and memory usage, one possible measure of the computational cost could be the number of machine cycles consumed by the deployed routines (i.e., matrix filling, system solving, calculation of currents and/or fields, etc.). Such an analysis is beyond the scope of the present paper and is left for a future investigation.