NUFFT for the Efficient Spectral Domain MoM Analysis of a Wide Variety of Multilayered Periodic Structures

IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 67, NO. 10, OCTOBER 2019 6551 NUFFT for the Efficient Spectral Domain MoM Analysis of a Wide Variety of Multilayered Periodic Structures Miguel Camacho ,Student Member, IEEE, Rafael R. Boix ,Member, IEEE, and Francisco Medina ,Fellow, IEEE Abstract— In this paper the Method of Moments (MoM) in the spectral domain is used for the analysis of multilayered structures containing periodic arrays of either patches or apertures. The patches and apertures may have many different geometries including complex surfaces limited by two parallel lines and two arbitrary curves, circular and elliptic rings, circular and elliptic arcs, and circular and elliptic sectors. Basis functions accounting for edge singularities are used in the approximation of the electric/magnetic current density on the patches/apertures, which enables a fast convergence of MoM with respect to the number of basis functions. Since the 2-D Fourier transforms of the basis functions cannot be obtained in closed-form, these Fourier transforms are efficiently computed by means of the nonuniform fast Fourier transform (NUFFT) algorithm. Results have been obtained for frequency-selective surfaces (FSSs), and for the elements used in the design of both reflectarray and metasurface antennas. The results obtained indicate that the software based on the NUFFT is only 15% slower than the standard spectral domain MoM software used for structures in which the 2-D Fourier transform of the basis functions is analytical, and between 50 and 80 times faster than CST. Index Terms— Fourier transforms, moment of methods (MoM), multilayered media, periodic structures. I. INTRODUCTION THE efficient analysis of multilayered periodic structures is of key importance for the design of microwave devices ranging from frequency-selective surfaces (FSSs) to reflectarray/transmitarray antennas and metasurface (MTS) leakywave antennas. FSSs paved the way by enabling the control of the transmission/reflection properties of very thin screens that Manuscript received March 25, 2019; revised May 15, 2019; accepted May 26, 2019. Date of publication June 6, 2019; date of current version October 4, 2019. This work was supported in part by the Engineering and Physical Sciences Research Council (EPSRC), U.K., through the EPSRC Center for Doctoral Training in Metamaterials under Grant EP/L015331/1, in part by Junta de Andalucía under Project P12-TIC-1435, and in part by the Spanish Ministry of Science, Innovation and Universities with EU FEDER funds under Project CICYT TEC2016-75103-C2-1-R and Project CICYT TEC2017-84724-P. (Corresponding author: Rafael R. Boix.) M. Camacho is with the Electromagnetic and Acoustic Materials Group, Department of Physics and Astronomy, University of Exeter, Exeter EX4 4QL, U.K., and also with the School of Engineering and Applied Sciences, University of Pennsylvania, Philadelphia, PA 19104-6391 USA (e-mail:[email protected]). R. R. Boix and F. Medina are with the Microwaves Group, Department of Electronics and Electromagnetism, College of Physics, University of Seville, 41012 Seville, Spain (e-mail: [email protected]; [email protected]). Color versions of one or more of the figures in this article are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TAP.2019.2920361 work as bandpass or band-reject filters when made of apertures or patches respectively [1], [2]. Although these structures were strictly periodic, researchers were able to extend the methods developed for the analysis of FSSs to the analysis of nonperiodic structures with a spatially smooth variation, thanks to the use of the local periodicity assumption (LPA). One example is the design of reflectarray antennas, which are made of periodic arrangements of printed patches of different dimensions on a grounded dielectric, usually illuminated by a horn antenna [3]. By optimizing the geometry of each element in the antenna, one can compensate for the phase difference encountered between the waves scattered by the different elements. If the geometry of the elements varies slowly along the antenna, then the LPA allows for the analysis of each element as if it was surrounded by a periodic environment [3]. The LPA was validated experimentally, and was an enabling breakthrough for the fast and accurate design of reflectarray antennas [4]. More recently, the LPA has also been applied for the design of printed patch holographic or MTS antennas, based on the leaky-wave phenomena [5], [6]. These antennas are made of very sub-wavelength printed patches on a thin grounded substrate, which are treated through equivalent surface impedance concepts [7]. The very-subwavelength patterning of the surface also allows for a reliable use of the LPA. When the equivalent surface impedance is properly modulated along the surface of the antenna, one can achieve a largely controllable conversion between surface and space waves [5]–[7], whose radiation can be spatially-tailored to obtain the desired radiation diagram. One of the numerical methods most frequently used for the analysis of multilayered periodic structures is the method of moments (MoM) in the spectral domain, which was introduced by Mittra and his coworkers [2]. The spectral domain MoM is very efficient provided that the 2-D Fourier transform of the basis functions used to approximate the magnetic current density on the apertures, or the electric current density on the patches, can be obtained in closed form. Two types of basis functions have been customarily used in the application of the spectral domain MoM: subsectional basis functions (such as rooftops, piecewise sinusoidal, etc.) and entire domain basis functions. Subsectional basis functions have the capability to be adaptable to a large variety of geometries, but their main drawback is that they lead to large MoM matrices to be inverted. Chan and Mittra [8] developed a fast Fourier 0018-926X © 2019 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. Authorized licensed use limited to: Universidad de Sevilla. Downloaded on May 22,2020 at 15:06:31 UTC from IEEE Xplore. Restrictions apply. 6552 IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 67, NO. 10, OCTOBER 2019 transform (FFT) approach for the efficient computation of the spectral domain MoM matrix entries when subsectional basis functions are used, but this approach is restricted to cases where the patches or apertures fit into a uniform rectangular grid limited by the contour of the periodic unit cell. Entire domain basis functions accounting for the edge singularity condition are specially suitable for the application of MoM to multilayered periodic structures because they ensure a fast convergence of MoM with respect to the number of basis functions, and therefore, lead to small MoM matrices that can be easily inverted [9]–[11]. The nice properties of the edge singularity basis functions are due to the fact that they constitute a complete orthonormal set for the kernel of the integral equations to be solved as pointed out by Lerer and Schuchinsky [12]. Edge singularity basis functions have been proposed for simple canonical geometries such as rectangles [10], [12], parallelograms [12], circles [7], [13] and ellipses [7]. However, they have also been extended to more sophisticated shapes such as rings [14], split rings [11] or even more general surfaces limited by two parallel lines and two arbitrary curves [15]. The problem is that the 2-D Fourier transform of the edge singularity basis functions can be only obtained for the simplest geometries (rectangles, parallelograms, circles, and ellipses) [12], which prevents the use of these basis functions in the spectral domain MoM analysis of periodic structures with more complex shapes unless sophisticated mathematical tricks are used for each particular shape [11], [14]. In this paper the authors present a novel implementation of the spectral domain MoM for the analysis of multilayered periodic structures in which the Fourier transforms of the basis functions are numerically computed by means of the nonuniform fast Fourier transform (NUFFT) algorithm [16], [17]. This novel implementation has the relevant advantage that it can deal with edge singularity basis functions for which the 2-D Fourier transforms of the basis functions cannot be obtained in closed-form, and therefore, it enables the efficient MoM analysis of a large class of multilayered periodic structures. Also, the analysis of all these structures is carried out within a CPU time which is only slightly larger than that required in the cases where the 2-D Fourier transforms of the basis functions are available in closed-form. In previous papers, specially tailored mathematical tricks (interpolations in terms of Chebyshev polynomials, truncation of infinite series based on the properties of Bessel functions, etc.) have been used for the spectral domain MoM analysis of some particular periodic structures with edge-singularity basis functions for which 2-D Fourier transforms are not available in closed-form [11], [14]. The uniqueness of the spectral domain MoM with NUFFT as proposed in this paper is that this latter approach can be efficiently applied to a wide variety of different geometries with similar success. In particular, in this paper the new approach is applied to the analysis of periodic structures where the patches/apertures are surfaces limited by two parallel lines and two arbitrary curves [15], and to the cases where the patches/apertures are circular or elliptic rings [14], circular or elliptic arcs [11], and circular or elliptic sectors. The use of the NUFFT has made it possible to cover all these different cases since it does not require an equispaced grid for the computation of the samples used in the determination of the Fourier transform. This is in contrast to the standard FFT, which requires an equispaced uniform grid for the computation of the samples, and therefore, is more limited in terms of the number of different geometries that can be handled [8]. In this paper the spectral domain MoM with NUFFT has been applied to the determination of the transmission properties of FSSs made of apertures, to the determination of the phase curves of periodic structures used in the design of linearly and circularly polarized reflectarray antennas under the LPA, and finally, to the determination of the iso-frequency dispersion curves and the equivalent surface reactance of patches used in MTS antennas. The results obtained have been compared with previously published results and with results provided by the commercial software CST, and good agreement has been found in all cases. Our spectral domain MoM with NUFFT software has proven to be between 50 and 80 times faster than CST. II. NUMERICAL PROCEDURE Fig. 1(a) and (b) show a periodic array of patches printed on a two-layered substrate backed by a conducting plane. Similarly, Fig. 1(c) and (d) show a periodic array of apertures in a conducting screen, embedded in a two-layered substrate. The patches, the backing conducting plane and the conducting screen containing the apertures are assumed to be perfect electric conductors (PECs) of negligible thickness. Each of the two layers of the substrate of both periodic structures has a thickness hi(i=1,2), a complex permittivity εi=ε0εr,i(1−jtanδi)(i=1,2), and a permeability μ0.Note that although the results presented in this paper are for periodic structures on one- or two-layered substrates, the derivations carried out in the rest of this Section for the structures of Fig. 1(a)–(d) can be easily extended to deal with general substrates containing an arbitrary number of layers. This only requires one to use adequate spectral domain dyadic Green’s functions, which can be obtained for generic multilayered substrates by means of the recurrent algorithm described in [18]. Also, while the periodic structure of Fig. 1(c) and (d) acts as a bandpass FSS, the structure of Fig. 1(a) and (b) cannot be used as a reject band FSS unless the conducting plane is substituted by an interface with a semi-infinite region of free space [1], [2]. This latter possibility can also be easily contemplated in the formulation provided below by simply changing the spectral domain dyadic Green’s function of the problem. We will assume that the periodic unit cell of the periodic structures of Fig. 1(a)–(d) may have either Nepatches [Ne=2 in Fig. 1(b)] or Neapertures [Ne=2 in Fig. 1(d)]. The possible shapes for these patches and apertures are shown in Fig. 2(a)–(g), and include surfaces limited by two parallel lines and two arbitrary curves (these surfaces encompass many different types of polygons, the barrel shape, the diabolo shape, etc.) of the type shown in Fig. 2(a), circular and elliptic rings [Fig. 2(b) and (c)], circular and elliptic arcs [Fig. 2(d) and (e)], and circular and elliptic sectors [Fig. 2(f) and (g)]. A plane wave is assumed to obliquely impinge on the two periodic structures of Fig. 1(a)–(d), its incidence direction Authorized licensed use limited to: Universidad de Sevilla. Downloaded on May 22,2020 at 15:06:31 UTC from IEEE Xplore. Restrictions apply. CAMACHO et al.: NUFFT FOR THE EFFICIENT SPECTRAL DOMAIN MoM ANALYSIS 6553 Fig. 1. (a) and (b) Side and top views of a periodic array of patches printed on a two-layered substrate backed by a conducting plane, respectively. (c) and (d) Side and top views of a periodic array of apertures embedded in a two-layered substrate, respectively. In both cases, a plane wave impinges on the periodic structures. being given by the angular spherical coordinates θinc and ϕinc. In the following, a time dependence of the type ejωtwill be assumed and suppressed throughout. A. MoM in the Spectral Domain In this subsection we describe the spectral domain MoM solution of the scattering problems posed in Fig. 1(a)–(d). Also, we briefly indicate how the solution of these scattering problems is used in the design of FSSs, reflectarray antennas and MTS antennas. In order to determine the fields scattered by the periodic structure of Fig. 1(a) and (b), we need to determine the surface electric current density excited on the patches of the metallized interface z =0, J(x,y), by the impinging plane wave. This electric current density is the solution of the following electric Fig. 2. Patches or apertures appearing in the unit cells of the periodic structures of Fig. 1. (a) Surface limited by two lines parallel to the x-axis, and by two arbitrary curves x=l− p(y)and x=l+ p(y). (b) Circular and (c) elliptic rings. (d) Circular and (e) elliptic arcs. (f) Circular and (g) elliptic sectors. field integral equation (EFIE): ˆz×Eap(x,y, z =0)++∞  m=−∞ +∞  n=−∞Pmn GEJ p (x−x,y−y, z =0, z =0)·J(x,y)dxdy=0 (x,y)∈P00 (1) where Pmn (m,n=...,−1,0,1,...)is the metallized portion of the z =0 plane within the mn-th periodic unit cell, Authorized licensed use limited to: Universidad de Sevilla. Downloaded on May 22,2020 at 15:06:31 UTC from IEEE Xplore. Restrictions apply. 6554 IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 67, NO. 10, OCTOBER 2019 Eap(x,y, z )is the electric field that would be generated in all space by the plane wave impinging on the multilayered structure of Fig. 1(a) and (b) in the absence of the patches, and GEJ p(x−x,y−y, z =0, z =0)is a 2 ×2 matrix representing the dyadic Green’s function that links the transverse (to z ) electric field at the plane z =0, Et(x,y, z =0), with the electric current density J(x,y)that originates this electric field. Also, in order to determine the fields scattered by the periodic structure of Fig. 1(c) and (d), we need to determine the magnetic current density excited on the apertures of the conducting screen located at z =0, M(x,y), by the impinging plane wave. This magnetic current density is the solution of the following integral equation: Jaa(x,y)++∞  m=−∞ +∞  n=−∞Amn ×GJM a(x−x,y−y, z =0, z =0)·M(x,y)dxdy=0 (x,y)∈A00 (2) where Amn is the set of apertures existing in the z =0plane within the mn-th unit cell, Jaa(x,y)is the surface electric current density that would be induced at the conducting plane z =0 of Fig. 1(c) and (d) by the impinging wave if the apertures were not present, and GJM a(x−x,y−y, z =0, z =0) is a 2 ×2 matrix representing the dyadic Green’s function that links the surface electric current density (discontinuity in the tangential magnetic field) existing at the z =0 plane with the magnetic current density M(x,y)that originates this electric current density. Since J(x,y)and M(x,y)are Floquet-periodic functions of xand y, in order to solve the integral equations of (1) and (2), we only need to determine J(x,y)and M(x,y)within one periodic unit cell, e.g., the domain C00 ={0≤x≤a;0≤ y≤b}. Looking for a solution of the integral equations, we expand J(x,y)and M(x,y)in P00 and A00 respectively in terms of known basis functions as follows: J(x,y)= Ne  p=1 2  q=1 N1q  r=1 N2q  s=1 cpq rs Jpq rs (x,y)(3) M(x,y)= Ne  p=1 2  q=1 N1q  r=1 N2q  s=1 dpq rs Mpq rs (x,y). (4) In (3)/(4) the superindex p(p=1,...,Ne) indicates the number of the patch/aperture in P00/A00 for a total of Nepatches/apertures in Fig. 1(b) and (d). The superindex q(q=1,2) indicates the component of J(x,y)/M(x,y)which is being approximated by the basis function Jpq rs (x,y)/Mpq rs (x,y). And since J(x,y)/M(x,y) depends on two coordinates (which can be cartesian coordinates, polar coordinates [7], [11], [14] or stretched polar coordinates suitable for patches and apertures with elliptical shape [7], [12]), the subindex r(r=1,...,N1q)isswepttofit the dependence of J(x,y)/M(x,y)on one of the coordinates, while the subindex s(s=1,...,N2q) is swept to fit the dependence of J(x,y)/M(x,y)on the other coordinate. Thus, the total number of basis functions per patch (aperture) turns out to be Nb=N11 N21 +N12 N22. When (3) and (4) are substituted in (1) and (2), and Galerkin’s version of MoM is applied, the following systems of equations are obtained for the unknown coefficients cpq rs and dpq rs of (3) and (4) Ne  p=1 2  q=1 N1q  r=1 N2q  s=1 ij,pq kl,rs cpq rs =eij kl (5) (i=1,...,Ne;j=1,2;k=1,...,N1j;l=1,...,N2j) Ne  p=1 2  q=1 N1q  r=1 N2q  s=1 ij,pq kl,rs dpq rs =fij kl (6) (i=1,...,Ne;j=1,2;k=1,...,N1j;l=1,...,N2j). If we invoke Parseval’s identity for 2-D Fourier transforms, the MoM matrix entries ij,pq kl,rs and ij,pq kl,rs of (5) and (6) can be expressed in the spectral domain as double infinite summations given by [2] ij,pq kl,rs =ab +∞  m=−∞ +∞  n=−∞ Jij,d kl (kxm,kyn)∗t · GEJ,c p(kx=kxm,ky=kyn, z =0, z =0)· Jpq,d rs (kxm,kyn) (7) ij,pq kl,rs =ab +∞  m=−∞ +∞  n=−∞ Mij,d kl (kxm,kyn)∗t · GJM,c a(kx=kxm,ky=kyn, z =0, z =0)· Mpq,d rs (kxm,kyn) (8) where kxm =k0sin θinc cos ϕinc +2πm/a(k0=ω√ε0μ0= 2π/λ0)andkyn =k0sin θinc sin ϕinc +2πn/b. In (7) and (8) the 2 ×2 matrices  GEJ,c p(kx,ky, z =0, z = 0)and  GJM,c a(kx,ky, z =0, z =0)stand for the continuous 2-D Fourier transforms of GEJ p(x,y, z =0, z =0)and GJM a(x,y, z =0, z =0). It can be easily shown that  GJM,c a(kx,ky, z =0, z =0) = GEJ,c a(kx,ky, z =0, z =0)−1·0−1 10 (9) where  GEJ,c a(kx,ky, z =0, z =0)is the continuous 2-D Fourier transform of the dyadic Green’s function relating the transverse (to z ) electric field and the electric current density at the z =0 plane of Fig. 1(c). The matrix  GEJ,c a(kx,ky, z =0, z =0)is not to be confused with  GEJ,c p(kx,ky, z =0, z =0)since z =0inFig.1(c)is in between the two layers and the two-layered medium is limited by two semi-infinite half-spaces up and down, while z =0 in Fig. 1(a) is above the two layers and the two-layered medium is limited by one semi-infinite halfspace up and by a conducting plane down. Both matrices  GEJ,c a(kx,ky, z =0, z =0)and  GEJ,c p(kx,ky, z =0, z =0) can be obtained in closed-form as explained in [18]. Authorized licensed use limited to: Universidad de Sevilla. Downloaded on May 22,2020 at 15:06:31 UTC from IEEE Xplore. Restrictions apply. CAMACHO et al.: NUFFT FOR THE EFFICIENT SPECTRAL DOMAIN MoM ANALYSIS 6555 The spectral domain basis functions  Jpq,d rs (kxm,kyn) and  Mpq,d rs (kxm,kyn)(p=1,...,Ne;q=1,2;r= 1,...,N1q;s=1,...,N2q) of (7) and (8) stand for the discrete 2-D Fourier transforms of Jpq rs (x,y)and Mpq rs (x,y) respectively, and are given by  Jpq,d rs (kxm,kyn)=1 ab P00 Jpq rs (x,y)e−j(kxmx+kyn y)dxdy (10)  Mpq,d rs (kxm,kyn)=1 ab A00 Mpq rs (x,y)e−j(kxmx+kyn y)dxdy. (11) Finally, the coefficients eij kl and fij kl of the two systems of equations (5) and (6) can be obtained in the spectral domain as eij kl =−ab Jij,d kl (kx0,ky0)∗t·Eap(x,y, z =0) ×e−jk0(sin θinc cosϕincx+sin θinc sinϕinc y)(12) fij kl =−ab Mij,d kl (kx0,ky0)∗t·Jaa(x,y) ×e−jk0(sin θinc cosϕincx+sin θinc sinϕinc y)(13) where the factor e−jk0(sin θinc cosϕinc x+sin θinc sinϕinc y)has been explicitly included in (12) and (13) to absorb the dependence of Eap(x,y, z =0)and Jaa(x,y)on xand y. Note that both Eap(x,y, z =0)and Jaa(x,y)include a phase factor of the type e+jk0(sin θinc cosϕincx+sin θinc sinϕinc y)(nonvanishing for oblique incidence, i. e., θinc = 0), which has to be removed in (12) and (13) by means of the annihilating phase factor e−jk0(sin θinc cosϕincx+sin θinc sinϕinc y)since the coefficients eij kl and fij kl are not dependent on xand y. As mentioned above, the periodic structure of Fig. 1(c) and (d) is a bandpass FSS. Once the integral equation (2) is solved by means of MoM and M(x,y)is obtained, the reflection and transmission coefficients of the mn-th Floquet harmonic of the FSS, Rmn and Tmn, can be obtained in terms of M(x,y)as explained elsewhere (e. g., see [2, eqs. (75)–(78)] and [19, eqs. (24)–(27)]). As commented in Section I of the paper, the design of reflectarray antennas is customarily carried out by means of the LPA. In the frame of the LPA, the phase shift in each reflectarray element is computed as if the element were surrounded by a periodic environment. In this sense, the configurations of Fig. 1(a) and (b) represent the typical periodic structures that are used for the characterization of the elements of a reflectarray. Let us assume that one of these periodic structures contains the element used in the design of a linearly polarized antenna, and let us assume that the value of the operating frequency does not allow the excitation of grating lobes when the plane wave impinges on the periodic structure (which implies that max(a,b)<λ 0/(1+sin θinc)[1]). Under these conditions, the reflectarray element surrounded by a periodic environment can be characterized by means of a 2 ×2 linear polarization reflection matrix, RLP, defined as [20] Eref x Eref y=RLP ·Einc x Einc y=Rxx Rxy Ryx Ryy ·Einc x Einc y(14) where Einc xand Einc yare the complex amplitudes of the transverse (to z ) components of the electric field of the impinging plane wave, and Eref xand Eref yare the complex amplitudes of the transverse (to z ) components of the electric field of the reflected m=n=0 Floquet harmonic (the only reflected propagating harmonic in the absence of grating lobes). Once the integral equation (1) is solved, the quantities Eref xand Eref y can be readily obtained in terms of J(x,y). In reflectarray antenna design, the phases of the matrix coefficients Rxx and Ryy of (14) are crucial for the determination of the element dimensions that provide the adequate phase shift for a prescribed radiation pattern, and the magnitudes of Rxy and Ryx are useful to give an estimation of the cross-polarization level introduced by the element. If the reflectarray to be designed is a circularly polarized antenna, in the absence of grating lobes each reflectarray element has to be characterized by means of a 2 ×2 circular polarization reflection matrix, RCP, which can be defined as [11] Eref RHCP Eref LHCP=RCP ·Einc RHCP Einc LHCP =RRHCP,RHCP RRHCP,LHCP RLHCP,RHCP RLHCP,LHCP ·Einc RHCP Einc LHCP (15) where Einc RHCP and Einc LHCP are the complex amplitudes of the right-hand circular polarization (RHCP) and left-hand circular polarization (LHCP) components of the electric field of the impinging wave, and Eref RHCP and Eref LHCP are the RHCP and LHCP components of the electric field of the reflected m=n=0 Floquet harmonic. These latter components can be readily obtained in terms of J(x,y)once the integral equation (1) is solved. As in the case of linear polarization antennas, the phases of the matrix coefficients RRHCP,RHCP and RLHCP,LHCP of (15) are the key parameters that make it possible to adjust the dimensions of each reflectarray element to obtain the required phase shift. The periodic structures of Fig. 1(a) and (c) may support bound (nonleaky) surface wave modes that are attenuated along the z -direction (for these modes, the phase constant is larger than k0). In order to obtain the propagation constant, kρsw, of the fundamental homogeneous Floquet mode m= n=0 that propagates in a direction making an angle αwith the positive xaxis, we have to redefine kxm and kyn in (7) and (8) as kxm =kρsw cosα+2πm/a(16) kyn =kρsw sinα+2πn/b.(17) With this definition of kxm and kyn, the values of kρsw for fixed values of αand k0(i.e., for fixed values of frequency) can be found as the solutions of the transcendental equations [7], [21] detij,pq kl,rs (kρsw,α,k0)=0 (18) detij,pq kl,rs (kρsw,α,k0)=0 (19) where ij,pq kl,rs and ij,pq kl,rs are the MoM matrices (7) and (8) of the linear systems of equations (5) and (6) respectively. In practical applications in which the periodic structure of Fig. 1(a) and (b) acts as an impedance surface or an Authorized licensed use limited to: Universidad de Sevilla. Downloaded on May 22,2020 at 15:06:31 UTC from IEEE Xplore. Restrictions apply. 6556 IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 67, NO. 10, OCTOBER 2019 MTS made of sub-wavelength elements [7], [22], focus is concentrated on the isofrequency dispersion curves of the fundamental Floquet mode m=n=0. In these dispersion curves, kρsw =kρsw(α) is represented in polar coordinates for fixed values of frequency [7], [22]. The dispersion curves also provide important information about the phase and group velocities of the surface wave under study. The design of aperture MTS antennas is done through the concept of impenetrable equivalent surface reactance. By spatially modulating such reactance, one can achieve leakage from the aforementioned surface waves. The implementation of the optimal modulation requires the knowledge of the dependence of the reactance values on a pair of geometrical parameters of the element in the unit cell such as size, angle of rotation, eccentricity in the case of an elliptic element, etc. This dependence is usually represented as a colormap that provides the design curves from which one can extract the desired values of the modulated reactance [6], [7]. The equivalent impenetrable reactance of the MTS periodic structure is defined in terms of kρsw as [7] Xs=ζ0kρsw k02 −1 (20) where ζ0=√μ0/ε0is the free-space impedance. B. Basis Functions for the Electric and Magnetic Current Densities The choice of basis functions for J(x,y)/M(x,y)in the patches/apertures of Fig. 1(a)–(d) is crucial for an efficient implementation of the spectral domain MoM of Subsection II-A. As commented in Section I, the basis functions accounting for edge singularities ensure a fast convergence of MoM with respect to the number of basis functions, and therefore, lead to an efficient implementation of MoM since the size of the MoM matrices to be inverted (i.e., the coefficient matrices appearing in the linear systems of equations (5) an (6)) turns out to be small. In this subsection we define the edge singularity basis functions for all the geometries shown in Fig. 2(a)–(g). It can be verified that the components of the basis functions parallel to the edges of the patches/apertures become singular at the points of these edges. Also, the components of the basis functions that are normal to the edges of the patches/apertures become zero at the points of these edges, their derivative along the normal direction being singular at these points. These are physical constraints to be expected for the electric/magnetic current densities existing on the patches/apertures. The results presented in Section III will confirm that the edge singularity basis functions chosen for the geometries of Fig. 2(a)–(g) all lead to a fast convergence of MoM with respect to the number of these basis functions. Let Bpq rs (x,y)represent any of the two basis functions Jpq rs (x,y)and Mpq rs (x,y).Let{x,y, z }be a shifted system of coordinates centered at the point (x=xoc,y=yoc, z =0) of the domain C00 of Fig. 1(b) and (d) [see the new system of coordinates {x,y, z }in Fig. 2(a)–(g)]. The relation between the new cartesian coordinates and the original cartesian coordinates of Fig. 1(b) and (d) is given by x=x−xoc (21) y=y−yoc (22) z = z .(23) For patches and apertures which fit the geometrical shape shown in Fig. 2(a), we propose to use the following edge singularity basis functions introduced in [15] Bp1 rs (x,y) =1 lp(y) Tr−1x−ap(y) lp(y) 1−x−ap(y) lp(y)2Us−1y Lp1−y Lp2 ×x−ap(y) lp(y) dlp(y) dy+dap(y) dyˆx+ˆy (l− p(y)<x<l+ p(y);−Lp≤y≤+Lp) (r=1,...,N11;s=1,...,N21)(24) Bp2 rs (x,y) =1 Lp Ur−1x−ap(y) lp(y)1−x−ap(y) lp(y)2 × Ts−1y Lp 1−y Lp2ˆxl− p(y)≤x≤l+ p(y);−Lp<y<+Lp (r=1,...,N12;s=1,...,N22)(25) where Tr−1(·)(Ts−1(·))andUr−1(·)(Us−1(·)) are Chebyshev polynomials of first and second kind respectively. The functions ap(y)and lp(y)of (24) and (25) are given by ap(y)=l+ p(y)+l− p(y) 2(26) lp(y)=l+ p(y)−l− p(y) 2(27) where xand ycan be obtained in terms of xand yas shown in (21) and (22), and where x=l− p(y)and x=l+ p(y) (−Lp≤y≤+Lp)are the equations of the two curves that limit the left side and the right side of the patch/aperture of Fig. 2(a). Let us now introduce polar coordinates ρand ϕfor the shifted system of coordinates {x,y, z }. According to (21) and (22), these coordinates are related to the cartesian coordinates of Fig. 1(b) and (d) by means of ρcos ϕ=x−xoc (28) ρsin ϕ=y−yoc.(29) Authorized licensed use limited to: Universidad de Sevilla. Downloaded on May 22,2020 at 15:06:31 UTC from IEEE Xplore. Restrictions apply. CAMACHO et al.: NUFFT FOR THE EFFICIENT SPECTRAL DOMAIN MoM ANALYSIS 6557 For the circular ring of Fig. 2(b), we propose the use of the edge singularity basis functions [14], [23] Bp1 rs (ρ,ϕ) = Tr−12 ρ2p−ρ1pρ−ρ2p+ρ1p 2 1−2 ρ2p−ρ1pρ−ρ2p+ρ1p 22 ×ejs−N21+1 2ϕ(−sin ϕˆx+cos ϕˆy) (ρ1p<ρ<ρ2p;0≤ϕ<2π) (r=1,...,N11;s=1,...,N21;N21 odd)(30) Bp2 rs (ρ,ϕ) =Ur−12 ρ2p−ρ1pρ−ρ2p+ρ1p 2 ×1−2 ρ2p−ρ1pρ−ρ2p+ρ1p 22 ejs−N22+1 2ϕ ×(cos ϕˆx+sin ϕˆy)( ρ1p≤ρ≤ρ2p;0≤ϕ<2π) (r=1,...,N12;s=1,...,N22;N22 odd). (31) For the case of the elliptical ring of Fig. 2(c), we will assume that b2p/a2p=b1p/a1p=η<1 (i.e., we will assume that the axial ratios between the minor semi-axis and the major semi-axis are the same for the two ellipses limiting the rings). Under this assumption, we are going to define stretched polar coordinates αand βfor the shifted system of coordinates {x,y, z }, which are related to the coordinates xand yof Fig. 1(b) and (d) by means of αcos β=x−xoc (32) ηα sin β=y−yoc.(33) The edge singularity basis functions Bpq rs chosen for the elliptical ring of Fig. 2(c) can be written in terms of α and βas Bp1 rs (α,β) = Tr−12 a2p−a1pα−a2p+a1p 2 1−2 a2p−a1pα−a2p+a1p 22 ×ejs−N21+1 2β(−sin βˆx+ηcos βˆy) (a1p<α<a2p;0≤β<2π) (r=1,...,N11;s=1,...,N21;N21 odd)(34) Bp2 rs (α,β) =Ur−12 a2p−a1pα−a2p+a1p 2 ×1−2 a2p−a1pα−a2p+a1p 22 ejs−N22+1 2β ×(η cos βˆx+sin βˆy)(a1p≤α≤a2p;0≤β<2π) (r=1,...,N12;s=1,...,N22;N22 odd). (35) In the case of the elliptical arc of Fig. 2(e), we will assume that b2p/a2p=b1p/a1p=η<1 as in the case of the elliptical ring of Fig. 2(c). Under that assumption, the edge singularity basis functions we have chosen for the elliptical arc are Bp1 rs (α,β) = Tr−12 a2p−a1pα−a2p+a1p 2 1−2 a2p−a1pα−a2p+a1p 22 ×Us−12 β2p−β1pβ−β2p+β1p 2 ×1−2 β2p−β1pβ−β2p+β1p 22 ×(−sin βˆx+ηcos βˆy)(a1p<α<a2p;β1p≤β≤β2p) (r=1,...,N11;s=1,...,N21)(36) Bp2 rs (α,β) =Ur−12 a2p−a1pα−a2p+a1p 2 ×1−2 a2p−a1pα−a2p+a1p 22 × Ts−12 β2p−β1pβ−β2p+β1p 2 1−2 β2p−β1pβ−β2p+β1p 22(η cos βˆx+sin βˆy) (a1p≤α≤a2p;β1p<β <β 2p) (r=1,...,N12;s=1,...,N22)(37) where β1p=arctan tan ϕ1p η(38) β2p=arctan tan ϕ2p η,(39) and where ϕ1pand ϕ2pare the angles defined in Fig. 2(e). The coordinates αand βof (36) and (37) are the stretched polar coordinates defined in (32) and (33). The basis functions for the circular arc of Fig. 2(d) can be obtained as a particular case of (36) and (37) when η=1, a1p=ρ1p,a2p=ρ2p,β1p=ϕ1p,β2p=ϕ2p,α=ρand β=ϕ(where ρand ϕare the polar coordinates defined in (28) and (29)). In the case of the elliptical sector of Fig. 2(g), the basis functions proposed are Bp1 rs (α,β) = Tr−1α a0p 1−α a0p2α a0p πs β2p−β1p−1 ×sin πs 2+πs β2p−β1pβ−β2p+β1p 2 ×1−2 β2p−β1pβ−β2p+β1p 22 ×(−sin βˆx+ηcos βˆy) (0<α<a0p;β1p≤β≤β2p) (r=1,...,N11;s=1,...,N21)(40) Authorized licensed use limited to: Universidad de Sevilla. Downloaded on May 22,2020 at 15:06:31 UTC from IEEE Xplore. Restrictions apply. 6558 IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 67, NO. 10, OCTOBER 2019 Bp2 rs (α,β) =Ur−1α a0p1−α a0p2α a0p πs β2p−β1p−1 × cos πs 2+πs β2p−β1pβ−β2p+β1p 2 1−2 β2p−β1pβ−β2p+β1p 22 ×(η cos βˆx+sin βˆy)(0<α ≤a0p;β1p<β <β 2p) (r=1,...,N12;s=1,...,N22)(41) where the coordinates αand βhave been defined in (32) and (33), and where β1pand β2phave been defined in (38) and (39). As in the case of Fig. 2(d) in relation with (36) and (37), the basis functions for the circular sector of Fig. 2(f) can be obtained as a particular case of (40) and (41) when η=1, a0p=ρ0p,β1p=ϕ1p,β2p=ϕ2p,α=ρand β=ϕ (ρand ϕare the polar coordinates defined in (28) and (29)). In the particular case of the circular sector, the factor (α/a0p)[((πs)/(β2p−β1p))−1]and the sinusoidal functions of [(πs)/(2)+(πs)/(β2p−β1p)×(β−(β2p+β1p)/(2))]have been included in (40) and (41) to accommodate the behavior of the current density (electric and magnetic) far from the edges of the circular sector in accordance with the magnetic wall model of a microstrip circular sector proposed in [24], just as it was done in [13] for a circular microstrip patch. C. NUFFT of the Basis Functions Equations (7) and (8) show that the computation of the MoM matrix entries in the spectral domain requires the knowledge of the 2-D discrete Fourier transforms  Jpq,d rs (kxm,kyn) and  Mpq,d rs (kxm,kyn). The Fourier transforms of the edge singularity basis functions introduced in (24), (25), (30), (31), (34) to (37), (40) and (41) cannot be obtained in closedform, and this is an important drawback since the spectral domain MoM loses all its efficiency when the 2-D Fourier transforms of the basis functions are not available in closedform. In this subsection we show how the 2-D discrete Fourier transforms of all the edge singularity basis functions introduced in Subsection II-B can be numerically obtained in an efficient way by means of the NUFFT algorithm. Let Fpq rs (x,y)represent any of the two components of Jpq rs (x,y)and Mpq rs (x,y). In the case of the basis functions for the patch/aperture of Fig. 2(a), the discrete Fourier transform of Fpq rs (x,y), Fpq,d rs (kxm,kyn), would be given by (see (10) and (11))  Fpq,d rs (kxm,kyn) =e−j(kx0xoc+ky0yoc) ab +Lp −Lp ×l+ p(y) l− p(y) Fpq rs (x,y)×e−j(kx0x+ky0y) ×e−jm2π(x+xoc) a+n2π(y+yoc) bdxdy.(42) In order to numerically carry out the integral of (42), we are going to sample the integrand in the values of xand ygiven by y ip =−Lp+yp 2+ypi(i=0,...,Ny−1)(43) x ij,p=l− p(y i)+x ip 2+x ip j (i=0,...,Ny−1;j=0,...,Nx−1)(44) where yp=2Lp/Nyand x ip =2lp(y ip)/Nx(lp(y) was defined in (27)). Then, the integral of (42) can be approximately computed by means of the equation  Fpq,d rs (kxm,kyn)≈ Ny−1  i=0 Nx−1  j=0 e−j(kx0xoc+ky0yoc) ab x ipyp ×Fpq rs x ij,p,y ipe−jkx0x ij,p+ky0y ip ×e−jm2π(x ij,p+xoc) a+n2π(y ip+yoc) b. (45) In practice, the numerical computation of ij,pq kl,rs and ij,pq kl,rs requires to truncate the infinite summations of (7) and (8) within the interval −M/2≤m,n<+M/2. For these particular values of mand n, (45) can be rewritten as  Fpq,d rs (m,n)≈ NT−1  k=0 γke−j(mμk+nνk)(−M/2≤m,n<+M/2) (46) where γk(i,j)=e−j(kx0xoc+ky0yoc) ab x ipypFpq rs x ij,p,y ip ×e−jkx0x ij,p+ky0y ip(k=0,...,NT−1)(47) μk(i,j)=2πx ij,p+xoc a(k=0,...,NT−1)(48) νk(i,j)=2πy ip +yoc b(k=0,...,NT−1)(49) and where k(i,j)=i+Nyj(i=0,...,Ny−1;j=0,..., Nx−1),andNT=NxNy. The expression (46) is ready for the use of the NUFFT of type 1 as it appears in [17, Eqn. (1)], which can be efficiently implemented by means of the algorithm shown in [17, page 448]. A FORTRAN code for this algorithm is available at the internet link [25]. In the case of the patches/apertures of Fig. 2(b), (d), and (f),  Fpq,d rs (kxm,kyn)would be given by  Fpq,d rs (kxm,kyn) =e−j(kx0xoc+ky0yoc) ab ρ max,p ρ min,p ×ϕ max,p ϕ min,p Fpq rs (ρ,ϕ)×e−j(kx0ρcos ϕ+ky0ρsin ϕ) ×e−jm2π(ρcos ϕ+xoc) a+n2π(ρsin ϕ+yoc) bdϕ ρdρ (50) Authorized licensed use limited to: Universidad de Sevilla. Downloaded on May 22,2020 at 15:06:31 UTC from IEEE Xplore. Restrictions apply. CAMACHO et al.: NUFFT FOR THE EFFICIENT SPECTRAL DOMAIN MoM ANALYSIS 6559 where (ρ min,p,ρ max,p)=(ρ1p,ρ2p)in the case of Fig. 2(b) and (d), (ρ min,p,ρ max,p)=(0,ρ0p)in the case of Fig. 2(f), (ϕ min,p,ϕ max,p)=(0,2π) in the case of Fig. 2(b), and (ϕ min,p,ϕ max,p)=(ϕ1p,ϕ2p)in the case of Fig. 2(d) and (f). If the variables ρand ϕof the integrand of (50) are sampled, it is possible to obtain an approximate expression for (50) of the type shown in (46). The expressions of γk,μk and νkwould be in this case γk(i,j)=e−j(kx0xoc+ky0yoc) ab ρ ip ρpϕpFpq rs ρ ip,ϕjp  ×e−jkx0ρ ip cos ϕjp+ky0ρ ip sin ϕjp(k=0,...,NT−1) (51) μk(i,j)=2πρ ip cos ϕjp +xoc a(k=0,...,NT−1)(52) νk(i,j)=2πρ ip sin ϕjp +yoc b(k=0,...,NT−1)(53) where ρ ip =ρ min,p+ρp 2+ρpi(i=0,...,Nρ−1)(54) ϕjp =ϕ min,p+ϕp 2+ϕpj(j=0,...,Nϕ−1)(55) with ρp=(ρ max,p−ρ min,p)/Nρ,ϕp=(ϕ max,p− ϕ min,p)/Nϕ,andNT=NρNϕ. Finally, in the case of the patches/apertures of Fig. 2(c), (e) and (g),  Fpq,d rs (kxm,kyn)would be given by  Fpq,d rs (kxm,kyn) =e−j(kx0xoc+ky0yoc) ab α max,p α min,p ×β max,p β min,p Fpq rs (α,β)e−jkx0αcos β+ky0ηαsin β ×e−jm2παcos β+xoc a+n2πηαsin β+yoc bdβ ×ηα dα(56) where (α min,p,α max,p)=(a1p,a2p)in the case of Fig. 2(c) and (e), (α min,p,α max,p)=(0,a0p)in the case of Fig. 2(g), (β min,p,β max,p)=(0,2π) in the case of Fig. 2(c), and (β min,p,β max,p)=(β1p,β 2p)in the case of Fig. 2(e) and (g). If the variables αand βof the integrand of (56) are sampled, it is possible to obtain an approximate expression of (56) in terms of these samples of the type shown in (46). In this particular case, the expressions for γk,μkand νkwould be γk(i,j)=e−j(kx0xoc+ky0yoc) ab ηα  ip αpβ pFpq rs α ip,βjp ×e−jkx0α ip cos βjp+ky0ηα ip sin βjp(k=0,...,NT−1) (57) μk(i,j)=2πα ip cos βjp +xoc a(k=0,...,NT−1)(58) νk(i,j)=2πηα  ip sin βjp +yoc b(k=0,...,NT−1)(59) where α ip =αmin,p+αp 2+α pi(i=0,...,Nα−1)(60) βjp =β min,p+β p 2+β pj(j=0,...,Nβ−1)(61) with α p=(α max,p−α min,p)/Nα,β p=(β max,p− β min,p)/Nβ,andNT=NαNβ. III. NUMERICAL RESULTS AND VALIDATIONS In all the results presented in this Section we have used M=100, which means that 10 000 terms have been retained in the computation of the double infinite summations of (7) and (8). Also, the number of samples employed for an accurate determination of the NUFFT has been roughly NT≈10 000 (around 100 samples per coordinate). Finally, the number of basis functions used per patch/aperture for the convergence of the results ranges from Nb=4toNb=18. This implies that small MoM matrices have to be inverted, and therefore, clearly shows the advantage gained when using edge singularity basis functions in the application of the spectral domain MoM. Fig. 3(a) shows results for the transmission through bandpass FSSs made of arrays of apertures with bow tie or barrel shape. In this case the transmission refers to the ratio between the total power transmitted and the incident power, even in the presence of grating lobes, which means that Transmission = m n Tmnk2 xm+k2 yn<k2 0 (62) In Fig. 3(a) the results obtained with the spectral domain MoM described in Section II are compared with the results published in [15, Fig. 4], and with results obtained with the commercial software CST. Excellent agreement is found among the three sets of results when four basis functions of the type shown in (24) and (25) are used in the approximation of the magnetic current density on the apertures. We can confirm that our MoM software is around 80 times faster than CST in a laptop computer with an Intel Core i7-4790 processor at 3.6 GHz, four cores and 32 GB of RAM. Fig. 3(b) shows the convergence of the results of Fig. 3(a) as we increase the number of basis functions used in the application of MoM. The results obtained with three basis functions are nearly coincident with those obtained for four basis functions, which indicates the fast convergence provided by the edge singularity basis functions. Fig. 4(a) shows our spectral domain MoM-NUFFT results for the transmission of the fundamental Floquet harmonic m=n=0 through a bandpass FSS consisting of a periodic array of circular ring apertures. In the generation of these MoM-NUFFT results we have used the basis functions of (30) and (31) to approximate the magnetic current density on the rings. In Fig. 4(a) our MoM-NUFFT results are compared with the results plotted in [26, Fig. 2], and with results provided by CST. Excellent agreement is noticed between our results and CST results. Good agreement is found between our results and those of [26] for frequencies below the first Wood’s anomaly (occurring for (fb)/c=1), but the results of [26] are Authorized licensed use limited to: Universidad de Sevilla. Downloaded on May 22,2020 at 15:06:31 UTC from IEEE Xplore. Restrictions apply.