An efficient approach for the computation of 2-D Green's functions with 1-D and 2-D periodicities in homogeneous media

IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 56, NO. 12, DECEMBER 2008 3733 An Efficient Approach for the Computation of 2-D Green’s Functions With 1-D and 2-D Periodicities in Homogeneous Media Ana L. Fructos, Rafael R. Boix, Member, IEEE, Francisco Mesa, Member, IEEE, and Francisco Medina, Senior Member, IEEE Abstract—This paper presents an algorithm for the acceleration of the series involved in the computation of 2-D homogeneous Green’s functions with 1-D and 2-D periodicities. The algorithm is based on an original implementation of the spectral Kummer-Poisson’s method, and it can be applied to the efficient computation of a wide class of infinite series. In the algorithm the number of asymptotic terms retained in Kummer’s transformation is externally controlled so that any of the series that has to be accelerated is split into one series with exponential convergence and another series with algebraic convergence of arbitrarily large order. Numerical simulations have shown that there is an “optimum” number of asymptotic terms retained in Kummer’s transformation for which the CPU time needed in the summation of the series is minimized. The CPU times required by Ewald’s method for the evaluation of 2-D Green’s functions with 1-D and 2-D periodicities have been compared with those required by the present algorithm, and the algorithm has been found to be between 1.2 and 3 times faster than Ewald’s method when working in “optimum” operation conditions. Index Terms—Convergence of numerical methods, Green’s functions, periodic structures, series. I. INTRODUCTION THE application of the method of moments [1] to the solution of periodic electromagnetic problems requires the computation of periodic Green’s functions. Thus, the 2-D homogeneous Green’s function with 1-D periodicity has been used in the determination of the scattering matrix of inductive obstacles in rectangular waveguide [2]. And that same 2-D Green’s function with 1-D periodicity has also been used in the analysis of the scattering of plane waves from 1-D periodic arrays of lossy strips in free space [3]. Recently, several authors have shown that the nonperiodic Green’s functions of multilayered media can be expressed in closed form as linear combinations of spherical and cylindrical waves in homogeneous media [4], [5]. According to this, the 3-D Green’s functions of multilayered media with 2-D periodicity can be easily obtained in terms Manuscript received September 14, 2007; revised March 06, 2008. Current version published December 30, 2008. This work was supported in part by the Spanish Ministerio de Educación y Ciencia and European Union FEDER funds (project TEC2007-65376), and in part by Junta de Andalucía (project TIC-253). A. L. Fructos, 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]). F. Mesa is with the Microwaves Group, Department of Applied Physics 1, ETS de Ingeniería Informática, University of Seville, 41012-Seville, Spain (e-mail: [email protected]). Digital Object Identifier 10.1109/TAP.2008.2007281 of both 3-D homogeneous Green’s functions with 2-D periodicity and 2-D homogeneous Green’s functions with 2-D periodicity [6], [7]. This approach for the computation of periodic Green’s functions in multilayered media has been followed in [7] as a mean to compute the reflection and transmission properties of frequency selective surfaces embedded in multilayered media. The previous examples show that the development of efficient and accurate algorithms for the computation of the 2-D homogenous Green’s functions with 1-D and 2-D periodicities is very useful as an intermediate step for studying the scattering by periodic structures both in homogeneous and multilayered media. Periodic Green’s functions in homogeneous media can be written either as spatial infinite series or as spectral infinite series when Poisson’s formula is used [8]. When losses are negligible, both the spatial and the spectral series are slowly converging, and a large number of analytical and numerical methods have been applied to accelerate the convergence of these series. Among the analytical methods, we should mention the spatial Kummer’s transformation [9] (which does not ensure convergence for complex phase shifts [8]), the spectral Kummer’s transformation [10] (which does not make it possible the analytic extraction of the Green’s functions singularity [8]), the spectral Kummer’s transformation combined with Poisson’s formula [6], [8], [11]–[13], the Veisoglu’s transformation [8], [14], [15] (which performs very well for the computation of periodic Green’s functions with 1-D periodicity but cannot be easily extended to the computation of Green’s functions with 2-D periodicity [8]), the perfectly matched layers method [16], [17] (which has proven to be powerful in the computation of homogeneous and multilayered periodic Green’s functions with 1-D periodicity but has not been extended yet to the computation of periodic Green’s functions with 2-D periodicity), Ewald’s method [8], [15], [18]–[21], and the lattice sums method [22]–[26] (which turns out to be very convenient when the periodic Green’s function has to be computed a large number of times, but has only been implemented for square lattices in the case of periodic Green’s functions with 2-D periodicity [23]). Among the numerical methods employed for the acceleration of the series involved in the computation of periodic Green’s functions, we should mention the Shanks’s transformation [27], the -algorithm [28], the Chebyschev-Toeplitz algorithm [29], the Levin’s T transform [30], and the recently developed summation-by-parts algorithm [31]. As commented in [8], these latter numerical acceleration methods can be directly applied to the spatial and spectral series as well as to the fast converging series resulting 0018-926X/$25.00 © 2008 IEEE Authorized licensed use limited to: Universidad de Sevilla. Downloaded on June 09,2020 at 17:00:47 UTC from IEEE Xplore. Restrictions apply. 3734 IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 56, NO. 12, DECEMBER 2008 from the application of some analytical acceleration methods (such as the Kummer’s transformation methods, the perfectly matched layers method, the Ewald’s method, and the lattice sums method). Several papers have carried out a comparative study among the different analytical and numerical methods for the acceleration of the series [8], [25], [32], and all of them seem to indicate that when the Green’s functions have to be evaluated for a moderate number of pairs of source and observation points, Ewald’s method is the best method in most scenarios because of its versatility and good compromise between accuracy and efficiency (only when the Green’s functions have to be evaluated for a large number of pairs of source and observation points, the lattice sums method is preferred). In the present paper an original algorithm is reported to accelerate the series involved in the computation of 2-D homogeneous Green’s functions with 1-D and 2-D periodicity. The algorithm is based on a novel implementation of the spectral Kummer-Poisson’s method. As in [9], [10], the use of an asymptotic expansion of arbitrarily large order in the application of Kummer’s transformation makes the algorithm very efficient. However, unlike [9] and [10], the algorithm basically relies on Poisson’s formula rather than on the specific mathematical expression of the terms of the series. This makes the algorithm very versatile because the same philosophy underlying its application to the computation of the single series involved in problems with 1-D periodicity can easily be extended to the computation of the double series involved in problems with 2-D periodicity, even though the axes of the periodic lattice are arbitrarily skewed. In the results section of the paper, we compare the algorithm with Ewald’s method, and we show that our algorithm can be made more efficient than Ewald’s method when asymptotic terms of sufficiently high order are retained in Kummer’s transformation. Further comparisons with other acceleration methods have not been carried out since the comparison with Ewald’s method can be taken as a good benchmark [8], [25], [32]. The paper is organized as follows. Section II shows all the mathematical derivations concerning the application of the new algorithm to the computation of the 2-D Green’s function with 1-D periodicity (Section II-A) and 2-D periodicity (Section II-B) in homogeneous media. In Section III we show the fast convergence rate of the series arising from the application of the spectral Kummer-Poisson’s method in the new algorithm, and we also compare the CPU times required by the new algorithm with those required by Ewald’s method. Conclusions are summarized in Section IV. II. DESCRIPTION OF THE NUMERICAL METHOD A. The 1-D Green’s Function With 1-D Periodicity In the following, a time dependence of the type will be assumed and suppressed throughout. Let us consider a one-dimensional array of line sources that are parallel to the axis and periodically located in a nonmagnetic homogeneous lossy medium along the direction (see Fig. 1) characterized by and . The periodic 2-D Green’s function Fig. 1. Planar array of line sources with 1-D periodicity along the x direction.  + nd ^ x =( x + nd )^ x + y ^ y is the 2-D vector pointing at the line sources ( n =... ; 0 1 ; 0 ; 1 ; ...) and  = x ^ x + y ^ y is the 2-D vector for the observation point. generated by the array of line sources is given by the well-known spatial series (1) where is the host medium complex wavenumber, is the period, is the zeroth-order Hankel function of the second kind [33], is the phase per period imposed either by a wave emitted by the array (positive sign) or by a wave incident on the array (negative sign), and are the angular spherical coordinates that indicate the propagation direction of these waves (see [3, eq. (1a)]), and is defined as with being the coordinates of the observation point, and the coordinates of the line sources . The th term of the series (1) decays like as , and therefore, the convergence of (1) is extremely slow when losses are negligible [8]. In this particular case, Poisson’s formula provides an alternative series for the computation of with a faster convergence rate. This alternative series is given by [8] (2) where and . According to [8], the series obtained in (2) is a spectral representation of the series (1). The th term of the series (2) has an exponential decay of the type when , but has a decay of the type when [8]. This means that the convergence of the spectral series is also slow when . A better alternative to speed up the convergence rate of (1) was reported in [15], [21], where the authors made use of Ewald’s method [18] to express the 2-D Green’s function with 1-D periodicity as (3) Authorized licensed use limited to: Universidad de Sevilla. Downloaded on June 09,2020 at 17:00:47 UTC from IEEE Xplore. Restrictions apply. FRUCTOS et al.: COMPUTATION OF 2-D GREEN’S FUNCTIONS 3735 In the above equation is the th-order exponential integral function defined in [33, p. 228] and where stands for the complementary error function [33, p. 297]. It can be shown that the two series of (3) over the integer index exhibit extremely fast Gaussian convergence [21]. The parameter appearing in (3) is the splitting parameter that simultaneously controls the convergence of the two series over [21]. Although Capolino et al. have suggested an expression for (see [21, eq. (34)]), in this paper we have used a simplified expression given by (4) with . By comparison with (4), the expression for provided in [21, eq. (34)] includes a third argument dependent on , which is the number of terms externally imposed to achieve convergence in the series over in (3). Our numerical simulations have shown that this does not bring any advantage in the computation of in terms of increased accuracy and/or reduced CPU time. When is computed with an accuracy of eight significant digits and is computed via (4), it has been found that the value of typically ranges between and , and it is strongly dependent on the ratio . In this paper, we propose to accelerate the computation of the series (1) by applying Kummer’s transformation to its spectral version (2), and subsequently applying Poisson’s formula to the resulting series containing the asymptotic expression of the th term of (2) as (see [11, eq. (1)]) and ([13, eq. (4)]). When applying Kummer’s transformation, instead of retaining one single term in the asymptotic expansion of the th term of (2) as [8], [12], in this work we retain several terms [9], [10], which considerably improves the efficiency of the spectral Kummer-Poisson’s method. Specifically, the series in (2) will be written as a sum of the following two fast converging series: (5) where (6) (7) In the above equations is the zeroth-order modified Bessel function [33], are real positive coefficients larger than that have to be defined, and are complex coefficients that have to be computed in terms of and the coefficients as shown below. Note that (the th term of the series ) exponentially decays as since exponentially decreases with increasing [33]. Also, we are about to show that a judicious choice of the coefficients makes it possible to achieve a fast decay of (the th term of the series )as . When we take and in (5) to (7), these equations reduce to the Kummer-Poisson’s formula reported in ([12, eq. (13)]) and ([8, eq. (20)]). In this particular case, shows the following decay for large (see [8, eq. (20)]) (8) The decay shown above can be made much faster by taking in (6). In order to show this, let us carry out an asymptotic expansion of for large (9) where it has been made use of the series expansion of and that of the exponential function [33]. The coefficients of (9) are supplied in the Appendix for . An analogous asymptotic expansion can be carried out for when to give (10) From (9) and (10), the term in square brackets in (6) can be written for large as (11) with (12) According to (11) and (12), if we wish in (6) decays as fast as possible for , we simply have to choose the coefficients in such a way that they make vanish , i.e., the coefficients must Authorized licensed use limited to: Universidad de Sevilla. Downloaded on June 09,2020 at 17:00:47 UTC from IEEE Xplore. Restrictions apply. 3736 IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 56, NO. 12, DECEMBER 2008 be chosen so that they satisfy the following system of linear equations: (13) When the coefficients fulfill the above system of equations, the terms in the series (6) and (7) for large decay like (14) (15) where we have used the asymptotic expansion of the modified Bessel function for large arguments [33]. Equation (14) and (15) show that the series (6) and (7) will have a very fast convergence provided the values of the coefficients and are sufficiently large. However, all these coefficients cannot be chosen arbitrarily large for several reasons. First of all, we have to bear in mind that, for a predetermined accuracy, the reduction in the number of terms of the series that have to be computed does not necessarily imply a CPU time reduction in the evaluation of the series [10], and CPU time reduction is the key issue in the efficient computation of . Concerning the choice of , if we start from and increase the value of , the CPU times required to compute (6) and (7) decrease. However, for a certain threshold value of , the CPU times involved in the computation of the series reach a minimum, and after that value of , they begin to grow since the evaluation of every term of the series demands an increasing number of operations and this has an influence on the overall computational cost (see [10, Fig. 1]). In next section we will show that a value of ranging between 5 and 7 is usually an optimum choice that minimizes the CPU time consumption. The choice of the coefficients is also critical. If these coefficients are too large, the series (7) converges very quickly owing to the fast exponential decay of the modified Bessel functions but the series (6) converges slowly because a relatively large value of is required to make reach its asymptotic limit. And the opposite occurs when the coefficients are too small. Numerical simulations have shown that it is possible to reach a convenient tradeoff between the number of terms that have to be computed in each series when the following optimum values of are chosen: (16) Unfortunately, when the operation frequency is high (more specifically, when ) and the coefficients are chosen as suggested in (16), the series (6) and (7) converge to very large nearly equal numbers of opposite sign and a severe accuracy loss may appear in the final result for . This situation is very similar to that observed when using the optimum splitting parameter in the application of Ewald’s method [20], [21]. As it happens with Ewald’s method, by increasing the values of beyond , one avoids that and are two nearly equal numbers of opposite sign, and, therefore, an accurate Fig. 2. Spatial array of line sources with 2-D periodicity.  + m a + n a = ( x + ma + na )^ x +( y + ma + na )^ y is the 2-D vector pointing at the line sources ( m;n =... ; 0 1 ; 0 ; 1 ; ...) and  = x ^ x + y ^ y is the 2-D vector for the observation point. result for is obtained. Thus, for high frequencies, or equivalently for large period , the choice of the coefficients should be different than optimum. In particular, we suggest choosing (17) B. The 2-D Green’s Function With 2-D Periodicity Now, let us consider a 2-D array of line sources that are parallel to the axis and are periodically located in the nonmagnetic homogeneous lossy medium of complex permittivity and permeability (see Fig. 2). Let and be the primitive vectors of the periodic 2-D lattice, which are not necessarily orthogonal (i.e., the axes of the periodic array may be skewed). The spatial series for the periodic 2-D Green’s function generated by this 2-D array of line sources is given by (18) where is the projection on the plane of the wavenumber vector of a wave emitted by the array (positive sign) or a wave incident on the array (negative sign), and are the angular spherical coordinates that indicate the propagation direction of these waves, and is given by with and . By analogy with the series (1), the series (18) converges slowly when losses are negligible (in fact, the terms of (18) decay like as both and when ). If the 2-D version of Poisson’s formula is applied to (18), one obtains a spectral series with a faster convergence rate in the lossless case. This series is given by (19) where is the area of the unit cell of the 2-D lattice and Authorized licensed use limited to: Universidad de Sevilla. Downloaded on June 09,2020 at 17:00:47 UTC from IEEE Xplore. Restrictions apply. FRUCTOS et al.: COMPUTATION OF 2-D GREEN’S FUNCTIONS 3737 The terms of the series (19) decay like as both and . Therefore, although the convergence of (19) is faster than that of (18) in the lossless case, it is still slow. Ewald’s method can also be applied to compute .In fact, if we follow the approach adopted in ([21, Sect. IV]), Ewald’s method leads to the following expression: (20) Similarly to what happens with the two series in (3), the two series of (20) over the integer indexes and exhibit Gaussian convergence. The parameter of (20) is the splitting parameter that simultaneously controls the convergence of the two double series over and . By following a procedure completely similar to those reported in ([21, Sect. V]) and ([8, Appendix I]), it is possible to find a value of that not only minimizes the total number of terms required for convergence in the calculation of the two series of (20) over and ,but also avoids cancellation errors arising from the addition of very large nearly-equal numbers of opposite sign. This value of turns out to be (21) where has been chosen in this paper as in (4). Next, the rationale presented in the previous subsection for the efficient computation of is reused for the efficient computation of . In fact, we have applied Kummer’s transformation to the spectral series (19) and then we have applied Poisson’s formula to the asymptotic remainder (see [13, eq. (4)]). We have retained several terms in the asymptotic expansion of the terms of (19) as both and in order to improve the efficiency of the spectral Kummer-Poisson’s method. As a result of these derivations, we have obtained that can be expressed as the following sum of two fast converging double series: (22) where (23) (24) In (23) and (24) are real positive coefficients larger than that have to be defined (these coefficients play the role of the coefficients in (6) and (7)), and are complex coefficients that have to be computed in terms of and the coefficients as shown below (the coefficients play the role of the coefficients in (6) and (7)). As it happens with and exponentially decays as and , and a judicious choice of the coefficients makes it possible to achieve a fast decay of as and . In order to derive the value of the coefficients that optimize the convergence of (23), we need to rewrite the term in square brackets of (23) in the following way: (25) where and are polynomials in the variable of degrees and , respectively. According to (25), if we wish in (23) decays as fast as possible as and , we have to choose the coefficients in such a way that the coefficients of the polynomial are cancelled out and the degree of the polynomial is reduced as much as possible. Equation (25) indicates that this is achieved when the coefficients are enforced to satisfy the following system of linear equations: (26) where the coefficients can be obtained in terms of and . It turns out that the solution to the system of linear (26) can be obtained in closed form, and is given by (27) When the coefficients are chosen as shown in (27), the terms of the series (23) and (24) for large and decay like (28) (29) Authorized licensed use limited to: Universidad de Sevilla. Downloaded on June 09,2020 at 17:00:47 UTC from IEEE Xplore. Restrictions apply. 3738 IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 56, NO. 12, DECEMBER 2008 where the asymptotic expansion of the modified Bessel function for large arguments [33] has been used once again. In principle, (28) and (29) seem to indicate that the larger the values of and , the better the convergence of the series (23) and (24). However, as discussed in previous subsection, the values of and have to be optimized so as to minimize the CPU time involved in the computation of within a predetermined accuracy. The results section shows that the optimum value of that minimizes the CPU time employed in the computation of ranges between and , and depends on the required accuracy. Concerning the coefficients , their choice has to ensure a good balance between the required accuracy and the minimum number of terms that have to be evaluated in the series of (23) and (24) to compute and . Also, the coefficients should be chosen so as to avoid that and are large nearly-equal numbers of opposite sign since this leads to cancellation errors in the computation of via (22). All these constraints are approximately satisfied when the coefficients are chosen in the following way [see also (16) and (17)]: (30) III. NUMERICAL RESULTS In this section, we will compare the convergence rates of the series (6) and (7) [(23) and (24)] with those of the series (1) and (2) [(18) and (19)]. All these comparisons will be carried out for the case where the host medium of the line sources of Fig. 1 (Fig. 2) is lossless, since this is the worst case scenario for the convergence of the series (1) [(18)]. Also in this section, we will compare the CPU times employed in the computation of via (5) [(22)] with those employed by Ewald’s method via (3) [(20)]. A. The 2-D Green’s Function With 1-D Periodicity In order to study the convergence rates of the series (1), (2), (6) and (7), we define the following function of the integer index : (31) where stands for the th term of the aforementioned series. In Fig. 3, we plot versus for the series (6) considering different values of . Note that all the plotted curves asymptotically tend to straight lines of negative slope , which is in agreement with the asymptotic value given for in (14) when (note that the plot is logarithmic). Fig. 3 shows that for , a small number of terms is required in the series of (6) to compute with a high accuracy (in fact, the terms suffice to compute with an accuracy of eight significant digits). However, when the single-term Kummer-Poisson’s method is used as in [8] and [12], a large number of terms is needed to obtain with high accuracy (thus, whereas only the terms are Fig. 3. Plots of f ( n ) as a function of n for the series (6). Results are presented for several values of M . d =0 : 6 ; 1 x =0 : 5 d; 1 y =0 ; = = 4 ; =0 and tan  =0 . required to obtain with an accuracy of four significant digits, it is necessary to include the terms for an accuracy of eight significant digits). The study carried out in Fig. 3 is not repeated for the series (7) because the convergence of this series is basically independent of and is mainly controlled by the smallest value of . In Fig. 4(a)–(c), is plotted for the series (1), (2), (6), and (7) considering three different values of . Whereas the convergence of (1) and (7) is basically independent of , the convergence of (2) and (6) becomes faster with increasing as a consequence of the exponential factor that appears in the th term of (2) and (6) as . Note that the convergence of (1) is extremely slow for all values of since we have assumed that the host medium is lossless. The convergence of (2) is very slow when but it considerably improves when for the reasons commented above. By comparison with the convergence of (1) and (2), the convergence of (6) and (7) is always very fast [in the cases studied in Fig. 4(a)–(c), the terms needed to evaluate the two series with an accuracy of eight significant digits range between and ]. As increases beyond the values studied in Fig. 4(a)–(c), the convergence rate of (2) becomes comparable to those of (6) and (7). In fact, numerical simulations have shown that when , (2) is more efficient than (5) to (7) for the computation of since the number of operations carried out in (2) is considerably smaller than that carried out in (6) and (7) when the three series have a similar convergence rate [10]. We have carried out numerical simulations for the case where the lossless host medium used in Fig. 4(a)–(c) is substituted by a lossy host medium , and these simulations indicate that in the lossy case the convergence of (1) improves but the convergence of (2), (6), and (7) does not substantially differ from that shown in Fig. 4(a)–(c). These latter results are not shown for the sake of saving space. Numerical simulations have shown that whereas the convergence of the series over of (3) is generally faster than the convergence of the series (6) and (7) (see [21, Figs. 3–5]), the convergence of the series over of (3) is generally slower. Therefore, in order to compare the relative efficiency of (3) and (5) in Authorized licensed use limited to: Universidad de Sevilla. Downloaded on June 09,2020 at 17:00:47 UTC from IEEE Xplore. Restrictions apply. FRUCTOS et al.: COMPUTATION OF 2-D GREEN’S FUNCTIONS 3739 Fig. 4. Plots of f ( n ) as a function of n for the series of (1), (2), (6), and (7). Results are presented for several values of 1 y . d =0 : 6 ; 1 x =0 : 5 d;  = = 4 ; =0 ; tan  =0 and M =7 . (a) 1 y =0 , (b) 1 y =0 : 01 d , (c) 1 y =0 : 05 d . the computation of , we need to compare the CPU times required by both expressions to obtain within a given accuracy. This comparison is carried out in Tables I(a) and (b) where we show the CPU times required by (2), (3), and (5) to determine with an accuracy of both four and eight significant digits. For the CPU time required by the spectral series (2) has not been included because it is too long. According to Tables I(a) and (b), whereas the CPU time required by Ewald’s method is basically independent of , the CPU times required by (2) and (5) usually decreases with increasing TABLE I CPU TIMES FOR THE COMPUTATION OF G VIA (2), (3) AND (5). G IS COMPUTED WITH ACCURACIES OF BOTH 4AND 8SIGNIFICANT FIGURES.ALL CPU TIMES ARE NORMALIZED TO THE MINIMUM VALUE OBTAINED IN EACH TABLE. 1 x =0 : 5 d; = = 4 ; =0 AND tan  =0 .(A)ACCURACY OF FOUR SIGNIFICANT FIGURES,(B)ACCURACY OF EIGHT SIGNIFICANT FIGURES . This result is to be expected since Fig. 4(a)–(c) has already shown that the convergence of the series (2) and (6) becomes faster as increases. Tables I(a) and (b) also show that both (3) and (5) require longer CPU times as increases. In order to explain this latter behavior, we have to consider that in the case we have used the optimum values of and in (3) and (5) [i.e., the values shown in ([21, eq. (31)]) and (16)], but in the case the values employed for and are other than “optimum” to avoid cancellation errors arising from the addition of two large nearly-equal numbers of opposite sign [21], [20]. Tables I(a) and (b) show that the CPU times required by the single-term Kummer-Poisson’s method of [8], [12] are usually several times longer than those required by Ewald’s method. However, the CPU times required by the Kummer-Poisson’s method described in Section II-A are comparable to those of Ewald’s method when , and are typically between 1.2 and 2 times shorter than those of Ewald’s method when (except for the case in Table I(a)). Tables I(a) and (b) do not show results for because the CPU times obtained in those cases have always been found longer than those obtained for . B. The 2-D Green’s Function With 2-D Periodicity In order to study the convergence rates of the series (18), (19), (23), and (24), the following function of the integer index is introduced: (32) Authorized licensed use limited to: Universidad de Sevilla. Downloaded on June 09,2020 at 17:00:47 UTC from IEEE Xplore. Restrictions apply. 3740 IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 56, NO. 12, DECEMBER 2008 Fig. 5. Plots of f ( n ) as a function of n for the series (23). Results are presented for several values of M and two different lattices. 1 x =1 y =0 : 4 ; =  = = 4 and tan  =0 . (a) Orthogonal lattice: a =0 : 6^ x ; a =0 : 6^ y  . (b) Nonorthogonal lattice: a =(0 : 6^ x +0 : 3^ y ) ; a =(0 : 2^ x +0 : 4^ y )  . where stands for the th term of the aforementioned series (note that has been chosen in the definition of ). Fig. 5(a)–(b) plots versus for the series (23) in the cases of an orthogonal lattice and a nonorthogonal lattice. Different values of are considered in the figures. Note that the envelopes of all the plotted curves are parallel straight lines of negative slope , which is in agreement with the asymptotic value given for in (28) (note that the plot is logarithmic). Fig. 5(a)–(b) shows that the number of terms needed to ensure convergence in the computation of is substantially larger for than for when a high accuracy is demanded (in fact, for an accuracy of eight significant digits, at least the terms and are needed in the case , and only the terms and for ). The study carried out in Fig. 5(a)–(b) is not repeated for the series (24) because its convergence is basically controlled by the smallest value of . Fig. 6(a)–(b) shows results of for (18), (19), (23), and (24) when the periodic lattice is both orthogonal and nonorthogonal. As happens with (1), the convergence of the series (18) is also extremely slow in the lossless case. The convergence of (19) is better than that of (18), but it is still very slow. However, the convergence of (23) and (24) is much Fig. 6. Plots of f ( n ) as a function of n for the series of (18), (19), (23), and (24). Results are presented for two different lattices. 1 x =1 y =0 : 4 ; =  = = 4 ; tan  =0 and M =6 . (a) Orthogonal lattice: a =0 : 6^ x ; a =0 : 6^ y  . (b) Nonorthogonal lattice: a =(0 : 6^ x +0 : 3^ y ) ; a =(0 : 2^ x +0 : 4^ y )  . faster (in fact, the terms and suffice to evaluate the series (23) with eight significant digits, and the terms and suffice to evaluate the series (24) with the same accuracy). When losses are considered in the host medium , the convergence of (18) improves but the convergence of (19), (23), and (24) does not appreciably differ from that shown in Fig. 6(a)–(b). Finally, Tables II(a) and (b) present the CPU times required by (20) and (22) to obtain with an accuracy of both four and eight significant digits. When four significant digits are demanded, the CPU time required by Kummer-Poisson’s method becomes comparable to that of Ewald’s method when , and when eight significant digits are demanded, the two methods become comparable when . The optimum performance of the Kummer-Poisson’s method described in Section II-B seems to be obtained for in Table II(a) and for in Table II(b). In these optimum conditions, the CPU times required by the Kummer-Poisson’s method are typically between 1.5 and 3 times shorter than those required by Ewald’s method. Note that the CPU time ratio between the Kummer-Poisson’s method of Section II and Ewald’s method turns out to be slightly less favorable in the computation of than in the computation of [compare Tables I(a) and (b) with Tables II(a) and (b)]. This is attributed to the fact Authorized licensed use limited to: Universidad de Sevilla. Downloaded on June 09,2020 at 17:00:47 UTC from IEEE Xplore. Restrictions apply. FRUCTOS et al.: COMPUTATION OF 2-D GREEN’S FUNCTIONS 3741 TABLE II CPU TIMES FOR THE COMPUTATION OF G VIA (20) AND (22). G IS COMPUTED WITH ACCURACIES OF BOTH 4AND 8S IGNIFICANT FIGURES.A LL CPU TIMES ARE NORMALIZED TO THE MINIMUM VALUE OBTAINED IN EACH TABLE. 1 x =1 y =0 : 4 ; =  = = 4 ,AND tan  =0 .O RTHOGONAL LATTICE: a =0 : 6^ x n ; a =0 : 6^ y n : NONORTHOGONAL LATTICE: a =(0 : 6^ x +0 : 3^ y ) n ; a =(0 : 2^ x +0 : 4^ y ) n  .( A)ACCURACY OF FOUR SIGNIFICANT DIGITS AND (B)A CCURACY OF EIGHT SIGNIFICANT DIGITS that whereas (13) has to be numerically solved, the solution to the system of (26) is given by (27) in closed form. IV. CONCLUSION The authors have presented a new algorithm to accelerate the computation of the series involved in the obtaining of the 2-D homogeneous Green’s functions with 1-D and 2-D periodicities. The algorithm is based on a novel implementation of the spectral Kummer-Poisson’s method [13]. In this novel implementation the number of terms retained in Kummer’s transformation can be varied in such a way that every slowly converging series is split into two fast converging series, one with exponential convergence and the other with algebraic convergence of arbitrarily large order. The CPU times required by the new algorithm have been compared with those required by Ewald’s method when the 2-D periodic Green’s functions are computed with an accuracy of both four and eight significant digits. We have found that when a sufficiently large number of terms is retained in the application of Kummer’s transformation, the new algorithm is typically between 1.2 and 2 times faster than Ewald’s method in the computation of the Green’s function with 1-D periodicity, and between 1.5 and 3 times faster than Ewald’s method in the computation of the Green’s function with 2-D periodicity. It should be pointed out that whereas other sophisticated algorithms based on Kummer’s transformation [10], [9] are heavily dependent on the mathematical expression of the terms of the series that have to be accelerated, the new algorithm presented in this paper is very general and works equally well for problems with 1-D periodicity and 2-D periodicity (also involving lattices with arbitrarily skewed axes in this latter case). In fact, although the algorithm has been used in this paper for the summation of cylindrical waves, it can also be used for the summation of spherical waves. The authors are currently working on the extension of the algorithm to the computation of 3-D homogenous Green’s functions with 1-, 2-, and 3-D periodicities, and this will be the topic of a future work. APPENDIX The coefficients that appear in (9) to (13) are given by the following expressions in the range : REFERENCES [1] R. F. Harrington, Field Computation by Moment Methods. New York: Macmillan, 1968. [2] Y. Leviatan, P. G. Li, A. T. Adams, and J. Perini, “Single-post inductive obstacle in rectangular waveguide,” IEEE Trans. Microw. Theory Tech., vol. 31, pp. 806–812, Oct. 1983. [3] R. E. Jorgenson and R. Mittra, “Oblique scattering from lossy strip structures with one-dimensional periodicity,” IEEE Trans. Antennas Propag., vol. 38, pp. 212–219, Feb. 1990. [4] Y. L. Chow, J. J. Yang, D. G. Fang, and G. E. Howard, “A closed-form spatial Green’s function for the thick microstrip substrate,” IEEE Trans. Microw. Theory Tech., vol. 39, pp. 588–592, Mar. 1991. [5] A. G. Polimeridis, T. V. Yioultsis, and T. D. Tsiboukis, “A robust method for the computation of Green’s functions in stratified media,” IEEE Trans. Antennas Propag., vol. 55, pp. 1963–1969, Jul. 2007. [6] R. M. Shubair and Y. L. Chow, “Efficient computation of the periodic Green’s function in layered dielectric media,” IEEE Trans. Microw. Theory Tech., vol. 41, pp. 498–502, Mar. 1993. [7] Y. Yu and C. H. Chan, “Efficient hybrid spatial and spectral techniques in analyzing planar periodic structures with nonuniform discretizations,” IEEE Trans. Microw. Theory Tech., vol. 48, pp. 1623–1627, Oct. 2000. [8] G. Valerio, P. Baccarelli, P. Burghignoli, and A. Galli, “Comparative analysis of acceleration techniques for 2-D and 3-D Green’s functions in periodic structures along one and two directions,” IEEE Trans. Antennas Propag., vol. 55, pp. 1630–1642, Jun. 2007. Authorized licensed use limited to: Universidad de Sevilla. Downloaded on June 09,2020 at 17:00:47 UTC from IEEE Xplore. Restrictions apply.