On the fast approximation of Green's functions in MPIE formulations for planar layered media

IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 50, NO. 9, SEPTEMBER 2002 2185 On the Fast Approximation of Green’s Functions in MPIE Formulations for Planar Layered Media N. V. Shuley, Member, IEEE, R. R. Boix, Member, IEEE, F. Medina, Senior Member, IEEE, and M. Horno Abstract—The numerical implementation of the complex image approach for the Green’s function of a mixed-potential integralequation formulation is examined and is found to be limited to low values of 0 (in this context 0 =2 0 , where is the distance between the source and the field points of the Green’s function and 0 is the free space wavelength). This is a clear limitation for problems of large dimension or high frequency where this limit is easily exceeded. This paper examines the various strategies and proposes a hybrid method whereby most of the above problems can be avoided. An efficient integral method that is valid for large 0 is combined with the complex image method in order to take advantage of the relative merits of both schemes. It is found that a wideoverlappingregionexistsbetweenthetwotechniquesallowing a very efficient and consistent approach for accurately calculating the Green’s functions. In this paper, the method developed for the computation of the Green’s function is used for planar structures containing both lossless and lossy media. Index Terms—Antenna, complex images, Green’s function, numerical analysis, scatterer. I. INTRODUCTION IN THE numerical modeling of printed planar structures used in monolithic integrated microwave and millimeter structures, it is generally accepted that the method of moments (MoM) [1] is one of the most efficient and rigorous methods for the analysis of small-to-medium-sized structures (up to several wavelengths). MoM formulations in either the spatial or spectral domain involve the conversion of an operator based integro-differential equation into a matrix equation that is subsequently passed to the computer for numerical processing. The individual entries of the matrix require spatial or frequency integrations involving the Green’s function in the appropriate domain and suitable basis and testing functions. Multiple numerical integration is usually required for this step of the MoM procedure that subsequently makes the filling time a long time-consuming process. This is particularly true for problems that involve electrically large geometries or many frequencies, such as is required in ground penetration radar (GPR) studies, which are interested in the short pulse scattering from buried objects (e.g., [2] and [3]) or other problems involving large Manuscript received September 28, 1998. The work of N. V. Shuley was supportedinpartbytheMinistryof Scienceand Education,Madrid,Spainunder Grant SAB95-0424. N. V. Shuley is with the School of Computer Science and Electrical Engineering, University of Queensland, Brisbane, Qld. 4072, Australia. R. R. Boix and F. Medina are with the Department of Electronics and Electromagnetism, University of Seville, 41012-Seville, Spain. M. Horno, deceased, was with the Department of Electronics and Electromagnetism, University of Seville, 41012-Seville, Spain. Publisher Item Identifier 10.1109/TMTT.2002.802333. transverse distances, such as in the mutual coupling of microstrip arrays of more than a few elements. Other typically large geometry structures such as traveling-wave antennas or even relatively small arrays of resonant elements could benefit from the fast evaluation of the Green’s functions that would considerably speed up matrix fill times. For problems involving planar stratified media, one recent trend that is becoming widely acceptable is to implement a discrete complex image approach [4]–[9]. Originally proposed and formulated by Fang et al. [4], this method is essentially an alternative means of evaluating the time-consuming Sommerfeld integrals that are integral representations of the Green’s functions. The main advantage of a complex image scheme is that the numerical evaluation of the Sommerfeld integral in conjunction with an interpolation scheme is completely avoided. The method is particularly attractive for mixed-potential integral-equation (MPIE) formulations since all the Green’s functions for this formulation are typically cast in terms of scalar functions involving Hankel transforms. The advantage of the method becomes apparent when it is realized that the Hankel transform for a spherical wave can be expressed,in closed form, as a summation of cylindrical waves via the Sommerfeld identity [11]. In addition, the MPIE formulation benefits from the fact that the singularity of both the scalar and vector potentials are of the order of and, therefore, less singular than spatial-domain electric-field integral-equation (EFIE) formulations. Furthermore, the moment integrals associated with the singular term are known in closed form for certain basis functions. Unfortunately, the numerical implementation for the complex image approach is generally limited to low values of , where is the distance between the source and field points of the Green’s function and is the free-space wavelength. As previously mentioned, this is a clear limitation for problems of large dimension or high frequency where this limit is easily exceeded. This paper examines the various strategies and proposes a hybrid method whereby most of the above problems can be avoided. An efficient integral method that is valid for large is combined with the complex image method (CIM) in order to take advantage of the relative merits of both schemes. It is found that a wide overlapping region exists between the two techniques allowing a very efficient approach for accurately calculating the Green’s functions. Section II reviews the CIMs currently in use and explains the difficulties for large for both lossless and lossy structures. Section III discusses the corresponding integral method and Section IV presents some results, which combine Sections II and III to illustrate the technique. 0018-9480/02$17.00 © 2002 IEEE Authorized licensed use limited to: Universidad de Sevilla. Downloaded on July 06,2020 at 15:09:04 UTC from IEEE Xplore. Restrictions apply. 2186 IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 50, NO. 9, SEPTEMBER 2002 Fig. 1. Multilevel planar geometry. Without loss of generality, illustrations have been confined to the scalar potential Green’s function . This component containsalltheessentialfeaturesoftheSommerfeldintegral:slowly convergent oscillatory behavior with surface wave poles in the integrand, branch point singularities, and a nonzero quasi-static limit.The worst-casescenarioalsooccurswhensource andfield points are assumed to be coplanar, which corresponds to the boundary condition used in forming the MPIE. This condition is assumed throughout this paper. Again, without loss of generality, attention has been focused on the microstrip structure for illustrative purposes. II. CIM We begin by considering a multilayered medium, where . is the source point and is the field point. Let , where be a function representing either the scalar Green’s function or any of the vector potential entries of the Green’s dyad for the multilayered medium of Fig. 1. The function can be written in terms of its zeroth-order Hankel transform as follows: (1) where the Sommerfeld integration path (SIP) is passing above the poles and branch point at of (Fig. 2) [6]. The function can be obtained in a relatively straightforward way for the multilayered geometry of Fig. 1 by means of iterative algorithms [6]. However, once is known, the Sommerfeld integral of (1) cannot be obtained in closed form and, therefore, it is not possible to obtain an exact analytical expression for . The idea underlying the complex image approach is to obtain an accurate approximation of in such a way that the Sommerfeld integrals of (1) can be determined in closed form and, thus, an accurate analytical approximation of can be obtained accordingly. As we shall see later in this paper, the problem arising (a) (b) Fig. 2. (a) Two-level complex integration paths in the complex k -plane. SIP is the Sommerfeld integration path; the same path is also applicable to lossy structures. (b) Integration contour in the k -plane corresponding to (a). from this idea is that it is impossible to build accurate approximations of both and for all values of and , respectively. In the following, we will show how the CIM has been applied inthispaperfortheparticularcaseof .Wewilldistinguish between the two cases. Case 1: It is assumed that and lie on a plane parallel to the interfaces inside the th layer. In this case, for every value of , we approximate as where and for scalar potential for vector potential (2) Here, is the quasi-static limit or the asymptotic limit of as becomes large. are the poles of and are the corresponding residues at the poles. is the function to be approximated Authorized licensed use limited to: Universidad de Sevilla. Downloaded on July 06,2020 at 15:09:04 UTC from IEEE Xplore. Restrictions apply. SHULEY et al.: FAST APPROXIMATION OF GREEN’S FUNCTIONS IN MPIE FORMULATIONS FOR PLANAR LAYERED MEDIA 2187 by an exponential series and is given by (3) Case 2: and are at the interface between the th and th layers. Again, . This is nearly alwaysthe case required for single-layer microstrip-based problems without vias. In this case, the above formulas still hold with the exception of the quasi-static limit , which should be changed to (4) The complex image formulation as outlined above is essentiallythemethodofChowet al.[5]foramicrostripstructure and later extended by Aksun and Mittra [8]. In both of these papers, the scalar and vector potential Green’s functions in the spectral domain were transformed to the spatial domain by first analytically extracting from the spectral domain, as in the above formulation, the surface-wave poles and quasi-static images. The remaining terms are then handled by performing a Prony, least square Prony, or matrix pencil of functions (MPOF) [12] exponential fit to the remaining terms that constitute the complex images. Aksun [9] then generalized the CIM to multilayer planar structures as formulated here. However, in an effort to preserve generality, where a priori information of the surface-wave poles and quasi-static information may not be known explicitly, neither the poles, nor the quasi-static terms were extracted, but a method of proceeding directly to the entire Green’s function was considered. This generalization considered a two-step sampling scheme, which addresses the problem of trying to simultaneously preserve both the low behavior (predominantly far-field information) and the slow convergent asymptotic behavior of the spectrum (predominately near-field information), although Kipp and Chan [18] were the first to report on such a split sampling scheme in order to provide for a more effective utilization of the exponential approximation. Given that regularly spaced data is required by any of the exponential fitting algorithms, this strategy represents a good compromise between the conflicting requirements as noted above. The application of the two-step exponential fitting algorithm of Aksun [9] to (3) requires two sets of samples of the known function . The first set of samples is obtained along a deformed path in the -plane passing above the surface-wave poles and branch point of [see Fig. 2(a)]. Once the surface-wave poles and branch point are safely passed, the second set of samples is obtained along a second path over the real axis of the -plane. The combined path is thus equivalent to the SIP as required by the Hankel transform. As in [9], this paper uses the MPOF algorithm [12] to determine the unknown coefficients (amplitudes and exponents) of (3). This requires that the samples of are taken in terms of a real variable , which must be chosen in such a way that the complex function of , [appearing in the exponents of (3)] is a linear function of . Following Aksun [9], the transformations linking and in the paths and of the complex -plane are given by (5) Full guidelines for choosing , —the points that delimit the respective paths—may be found in [9]. Corresponding plots of and in the complex -plane are shown in Fig. 2(b). The equations for and used in this paper (5) are the same as those of Aksun, except for the fact that the wavenumber of the open layer has been used instead of the wavenumber of the source layer. This is an important difference, as discussed by Kipp and Chan [10]. Had we used in the expressions for and instead of using , then and would have had a branch point at while would still haveasinglebranchpointat [10]. Moreover, thereisno conceptualdifficultyinusing the sameequationsfor treating the lossy substrate case. The path is the same. Whereas if was usedinstead of , thecomplexdielectricwouldforcea change in path. This difference reinforces the notion of having the function to be approximated and the exact function sharing the same branch point. Analytical continuation requires the function to be approximated to be analytic not only on the integration path in the complex -plane, but also in all regions of in the neighborhood of the path.The paths in the -and -planesare shown in Fig. 2. Sampling and application of the MPOF algorithm then converts the Green’s function to a sum of exponentials, which may be directly inverted to the spatial domain via the Sommerfeld identity [11] (6) where is an arbitrary complex constant. If we are dealing with small , the above procedures work well. However, the far-field Green’s function is dominated by cylindrical surface waves and the expression provided for the Green’sfunction bythe CIMcontains aseriesofquasi-spherical waves of the type shown in the left-hand-side term of (6), which do not have a clear physical meaning. Thus, we expect problems with the CIM for large . We now illustrate some of the above observations with the canonical set of data of [5]. Here, wehaveused thetwo-stepmethodbothwith andwithoutextraction of the surface-wave poles. We also include in the plot the effect of using the free-space (called “case A” in the figures) as compared to the wavenumber of the dielectric layer Authorized licensed use limited to: Universidad de Sevilla. Downloaded on July 06,2020 at 15:09:04 UTC from IEEE Xplore. Restrictions apply. 2188 IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 50, NO. 9, SEPTEMBER 2002 Fig. 3. G for the microstrip data of [5] using the two-step MPOF algorithm.  =12 : 6 , t =1 : 0 mm, f =10 GHz. Case A corresponds to using k 0 k in (2)–(4) and k in (5), whereas case B corresponds to using k " 0 k instead of k 0 k in (2)–(4) and k p  instead of k in (5). For case A curves: N =35 samples on C , N =20 samples on C . For case B curves: N = 100 samples on C , N = 300 samples on C . With/without pole extraction refers to if the surface wave term(s) is/are extracted prior to exponential approximation. Fig. 4. Same data as for Fig. 3. f =50 GHz. (called “case B”) as per the above discussion (Fig. 3). The difference in sampling levels required to achieve correspondence of the curves at low frequency is shown in the caption. Figs. 3 and4arefordifferentfrequencydataof [5].It isto benotedthat, for Fig. 3, there is only one TM surface wave propagating, while in Fig. 4, there are three surface waves. The “exact” result using the method of averages [17] together with straight quadrature is shown for comparison. While both data sets agree very well with [5] for low (apart from a normalizing factor), there is clearly a region around where the image method fails. The type of failure is different for cases A and B as evidenced by the figures, which seem to confirm the findings of [10]; in particular, the data of [10, Fig. 10]. However, the maximum horizontal distance before divergence of the two results seem to differ by an order of magnitude compared to that reported in [10]. In order to explain this, one should note that, whereas in Figs. 3 and 4, the source and field points are on the boundary between the substrate and free space, while in [10, Fig. 10], the source and field points are inside the substrate between dielectric layers. Thus, the errors resulting from using the wrong branch point instead of the right one will be amplified. An attempt to increase the number of samples on so as to improve the far-field result brought about inconsistent and number-of-samples-dependent results, sometimes with an even worse result for an increased number of samples. This is not really surprising when the surface-wave term is not extracted from (2) and (3). Generally, the exponential fitting procedure has a problem with terms involving a behavior, i.e., in cases where the surface wave is dominating, as well as the branch-point problems of case B previously mentioned. By realizing that the problem in the MPOF algorithm exists when the surface waves are dominant, it is instructive to extract the surface-wave contribution before commencing the exponential approximation. This process makes use of the Hankel transform pair defined by the identity [16] (7) It might be thought that prior extraction of the surface-wave contributionwouldsignificantlyimprovethecompleximageapproach in the far field, as these terms are the dominating ones here. In Figs. 3 and 4, it is shown that this is true in case A, but is not true in case B. Since case B is the case treated by Aksun [9], thisseemstojustifyAksun’sapproachinavoidingtheextraction of the surface-wave terms. Determination of the surface-wave contribution is a time-consuming task (because the poles and associated residues have to be found) and, in case B of Figs. 3 and 4, the benefit of obtaining the surface-wave contribution is lost due to the failure of the MPOF algorithm for large . Itisalsoinstructivetoexaminethedependencyonthenumber of samples in the two regions of Fig. 2(a). Ideally, all results should be independent on the sampling scheme employed to make the exponential series approximation. However, due to the limitations of the numerical algorithm, this is rarely the case. A large number of numerical tests show that the results are quite stable with different numbers of samples ( ) on the part of the spectrum. This is understandable since once the quasi-static contribution is removed the two-step process makes the spectrum decay quickly to zero with very little fine detail. On the other hand, rapid changes in the spectrum can occur at low values of , which correspond to the far-field region of the spatial domain. It, therefore, seems plausible to increase the number of samples ( ) in the region for increased accuracy in the far field. Fig. 5 shows that good results for large can be obtained for a low number of samples. It should also be noted that the number of samples required to achieve good results for when in these figures is substantially less than those mentioned in [9]. This is probably because, in the current paper, the quasi-static limit and the surface-wave Authorized licensed use limited to: Universidad de Sevilla. Downloaded on July 06,2020 at 15:09:04 UTC from IEEE Xplore. Restrictions apply. SHULEY et al.: FAST APPROXIMATION OF GREEN’S FUNCTIONS IN MPIE FORMULATIONS FOR PLANAR LAYERED MEDIA 2189 Fig. 5. Convergence study of the number of samples used in the far-field region. Microstrip structure:  =12 : 6 , t =1 : 0 mm, f =50 GHz. Case A with extraction of all surface-wave poles. Fig. 6. Comparison between the CIM and quadrature (full line) for G for two levels of loss. The CIM (dotted line) uses case A with all poles extracted. f =10 GHz,  =4 : 4 , t =10 mm. N =35 samples on C , N =20 samples on C . term are extracted before applying the CIM. Furthermore, the right branch point is used in the approximation of the spectral Green’s function, which are steps that are not considered in [9]. However, when , no matter how many samples are used in the application of the CIM proposed in this paper, the results for tend to be very poor. Once the stable number of samples has been ascertained, the exponential series approximation is stable with frequency data, as shown in Fig. 6 for both the lossless and lossy cases. All data of this plot using the CIM were calculated using the same number of samples . Typical results for the CIM involving lossy dielectrics have not yet been treated in the literature. As mentioned previously, the -plane contour does not need to be modified if the case A formulation is followed. Fig. 6 shows some typical data for increasing loss tangent using the CIM as compared to quadrature. Again, the failure of the CIM for large and increasing loss tangent is apparent. Once again, it appears that the attenuated surface wave cannot be properly reconstructed. For the lossy case, the surface waves are decaying at while the spherical wave decays at , the surface wave will thus predominate at smaller for increasing loss. Care needs to be exercised when using the CIM in high loss situations where these surface waves control the quality of the approximation. III. INTEGRAL METHOD A. Lossless Case It appears from the above section that the only difficulty with the CIM is its inability to deal with the surface-wave behavior at large . An integral technique that is attractive for large behavior was proposed in [13]. Although the technique deteriorates for decreasing , this exactly suits our purposes. This approach considers the imaginary axis as part of an alternative path for the evaluation of the Sommerfeld integral. If Cauchy’s theorem is applied to the contour of the -plane, shown in Fig. 7(a), it can be shown that the Sommerfeld integral of (1) can be rewritten as (8) where is the modified Bessel function of zeroth order and aretheresiduesatthepoles of . Although the result contains an infinite integral in thefirstterm,themodifiedBesselfunction decaysexponentially as its argument increases, thus ensuring a quick computation of the infinite integral with increasing . The second term is an integral over the finite interval and is readily computed. Althoughthetechniquealsorequiresthepolesandtheirresidues to be evaluated, the total numerical effort required is actually inversely proportional to . As noted by the originators of this technique, the various terms become dominant according to the distance from the source. Thus, the terms that correspond tothedeficiencyofthecompleximageapproach,may bereadily identified. B. Lossy Case Including losses implies that the poles now become complex and move downward off the real axis (and toward the imaginary axis for increasing loss). This means that they would no longer be captured in the lossless path noted above. However, a set of two contours and [see Fig. 7(b)] can be used for the computation of the Green’s function in this case. This new set of two contours makes it possible to account for the complex poles of the Green’s function, as explained in [14]. Evaluation of the Green’s function using the two new contours requires that the traditional Sommerfeld branch cut (which divides the complex -plane into a proper and an improper sheet) be modified to a new cut. This new cut is defined as . Ordinarily, using different branch Authorized licensed use limited to: Universidad de Sevilla. Downloaded on July 06,2020 at 15:09:04 UTC from IEEE Xplore. Restrictions apply. 2190 IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 50, NO. 9, SEPTEMBER 2002 (a) (b) Fig. 7. Integration contours in the k -plane for the integral method of Section III. (a) Lossless case. (b) Lossy case. cuts means that the leaky-wave poles, which were all on the lower Riemann sheet in the case of the Sommerfeld cuts, can now be a problem if the contour is not chosen correctly. However, it can be shown that this particular path avoids the leaky-wave poles provided that any leaky-wave pole is not close to the negative imaginary axis. Thus, with this choice of integration path, any lossy poles present are enclosed, which was not the case of the lossless formulation. When Cauchy’s theorem is applied to the two closed contours and of Fig. 7(b), it is possible to obtain an expression for the Sommerfeld integral of (1), which is a generalization of (8) This expression is given by (9) Fig. 8. Comparison of the calculation for G using the CIM (case A with full surface-wave pole-extraction full line) and the integral method (dotted line) of Section III,  =12 : 6 , t =1 : 0 mm, N =35 samples on C , N =20 samples on C . IV. COMBINATION OF COMPLEX IMAGES AND INTEGRAL METHOD It is suggested that, for all values of , a combination of the CIM and the integral method will be a very efficient means of calculating any Green’s functions required in a mixed-po- tential formulation. As has been observed, the CIM can experience problems in the region of large . On the other hand, the integral method becomes very time consuming for small . Clearly there must be a region of overlap for the two approaches that will allow a certain amount of flexibility in an automatic selectionscheme.Twobeneficialresultsofthestudyoftheoverlap procedure are, firstly, that we are able to check the relative accuracy of each of the techniques against each other in respective regions and, secondly, since the terms in the integral approach have a direct correspondence to certain sections of the field, we are able to readily identify the source of any inaccuracies in the computational process. Fig. 8. shows the combination of the two techniques for the lossless case. The integral method has been computed for an upper limit of 10 for the modified Bessel term. If required, a higher truncation level gives a more accurate result for lower , but at increased computational effort. It should be noted thatthepointofdepartureoftheCIMdefinestheupperendpoint of a very wide overlapping region between the two methods and the lower end point can be changed by selection of an appropriate upper limit of the first integral in (5). As there are two or more surface waves interacting in the 30- and 50-GHz plots, the ripple structure is more pronounced in those plots. Fig. 9 shows some new results for the CIM for a lossy structure for the same parameters of Fig. 6. Here, the “exact” quadrature curve has been repeated for convenience. Losses in the system tend to flatten out the characteristic “plateau” of the curve and the corresponding failure of the MPOF algorithm is less dramatic, but as discussedearlier. The surface wave also experiences a decaying mechanism and eventually terms that are withoutattenuationdominate,whichleadtoaninterferencephenomena at large , shown here around for the Authorized licensed use limited to: Universidad de Sevilla. Downloaded on July 06,2020 at 15:09:04 UTC from IEEE Xplore. Restrictions apply. SHULEY et al.: FAST APPROXIMATION OF GREEN’S FUNCTIONS IN MPIE FORMULATIONS FOR PLANAR LAYERED MEDIA 2191 Fig. 9. Comparison between the integral method of Section III (dotted line) and quadrature (solid line) repeated from Fig. 6. f =10 GHz,  =4 : 4 , t =10 mm. The truncation level for the first integral of (9) has been set at 10 k . curve. This interaction has been interpreted as a Zenneck–surface-wave interaction and is due to the interaction of a pole and branch point [15]. Finally, it is also instructive to check the claim that the surface-wave contributions indeed dominate in the far field. We chose a value of corresponding to a nominal distance at which CIM experiences difficulties. We evaluated the contribution of the surface-wave term(s) compared to the total field, as formulated in (8) and(9). The percentageerror between the approximate (surface-wave field) and “exact” [complete (8)] was 4.2% for the data of Fig. 3 and 4.8% for Fig. 4, both reducing for larger . For the medium-loss case of Fig. 9 (), the error increased to 26%, only falling to 3.5% by . Although somewhat frequency dependent, such error levels indicate that, for high-accuracy evaluation of the Green’s function, the surface-wave terms are, in fact, not completely dominant when the CIM begins to fail. The organizational strategy in a generalized code construction would, therefore, be to compute the Green’s function using the CIM for small and then implement the integral method for intermediate distances. The combination of the CIM and integral method of Section III for the lossless case is shown directly in Fig. 8 and the lossy case may be inferred from a combination of Figs. 6 and 9. At large , the surface waves (both cylindrical and Zenneck waves) are dominating and the Green’s function may be set without any integration at all, as explained in[15].Theoverlapbetweenthese regionsmaybechosensomewhat arbitrarily given the large overlap between the strategies. Althoughthispaperhasconcentratedonthescalarpotential , similar observations may be made for the other components of the vector potential dyadic . V. CONCLUSION For the computation of the Green’s functions using an MPIE approach in a planar multilayered structure, it is not sufficient to rely solely on the CIM for the cases where large transverse distances are involved. This paper has demonstrated a strategy of combining an integral and the CIM for very rapid and accurate evaluation of the Green’s function. The method has also been extended to the case of lossy media. However, it must be stated that the price paid for this combined approach is ultimately an increase in CPU time over the direct application of the CIM alone in spite of the very efficient numerical integration called for in the approach, as discussed in this paper. ACKNOWLEDGMENT This study was performed while author N. V. Shuley was a Visiting Academic with the Faculty of Physics, Department of ElectronicsandElectromagnetics,UniversityofSeville,Seville, Spain. REFERENCES [1] R. F. Harrington, Field Computation by Moment Methods. Melbourne, FL: Krieger, 1983. [2] S. Vitebskiy and L. Carin, “Moment-method modeling of short-pulse scattering from the resonances of a wire buried inside a lossy, dispersive half-space,” IEEE Trans. Antennas Propagat., vol. 43, pp. 1302–1312, Nov. 1995. [3] S. Vitebskiy, K. Sturgess, and L. Carin, “Short-pulse plane-wave scattering from buried perfectly conducting bodies of revolution,” IEEE Trans. Antennas Propagat., vol. 44, pp. 143–151, Feb. 1996. [4] D. G. Fang, J. J. Yang, and G. Y. Delisle, “Discrete image theory for horizontal electric dipoles in a multilayered medium,” Proc. Inst. Elect. Eng., pt. H, vol. 135, pp. 297–303, Oct. 1988. [5] 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. Microwave Theory Tech., vol. 39, pp. 588–592, Mar. 1991. [6] G. Dural and M. I. Aksun, “Closed-form Green’s functions for general sourcesandstratifiedmedia,”IEEETrans.MicrowaveTheoryTech.,vol. 43, pp. 1545–1551, July 1995. [7] K. A. Michalski and J. R. Mosig,“Discrete compleximage mixedpotential integral equation analysis of microstrip patch antennas with vertical probe feeds,” Electromagnetics, vol. 15, pp. 377–392, 1995. [8] M.I. Aksun andR.Mittra, “Derivationofclosed-form Green’s functions for a general miocrostrip geometry,” IEEE Trans. Microwave Theory Tech., vol. 40, pp. 2055–2061, Nov. 1992. [9] M. I. Aksun, “A robust approach for the derivation of closed-form Green’s function for a general microstrip geometry,” IEEE Trans. Microwave Theory Tech., vol. 44, pp. 651–658, May 1996. [10] R. A. Kipp and C. H. Chan, “Complex image method for sources in bounded regions of multilayer structures,” IEEE Trans. Microwave Theory Tech., vol. 42, pp. 860–865, May 1994. [11] W. C. Chew, Waves and Fields in Inhomogeneous Media, ser. Electromagn. Waves. Piscataway, NJ: IEEE Press, 1995. [12] T. K. Sarkar and Pereira, “Using the matrix pencil method to estimate theparameters ofasumof complexexponentials,”IEEEAntennasPropagat. Mag., vol. 37, pp. 48–55, Feb. 1995. [13] J. R. Mosig and F. E. Gardiol, “Analytical and numerical techniques in the Green’s function treatment of microstrip antennas and scatterers,” Proc. Inst. Elect. Eng., pt. H, vol. 130, pp. 175–182, Mar. 1983. [14] J. R. Mosig and A. A. Melcon, “Green’s functions in layered media: Imaginary axis integration and asymptotic behavior,” in Proc. IEEE AP-S Conf., Baltimore, MD, July 1996, pp. 416–419. [15] Y. Brand, A. A. Melcon, J. R. Mosig, and R. C. Hall, “Large distance behavior of stratified media spatial Green’s functions,” in Proc. IEEE AP-S Conf., Montreal, QC, Canada, July 1997, pp. 2334–2337. [16] I. S. Gradshteyn and I. M. Ryzhik, Tables of Integrals, Series and Products. New York: Academic Press, 1980. [17] J. R. Mosig, Numerical Techniques for Microwave and Millimeter-Wave Passive Structures, T. Itoh, Ed. New York: Wiley, 1989, ch. 3. [18] R. A. Kipp and C. H. Chan, “A complex image method for the vertical component of the magnetic potential of a horizontal dipole in layered media,” in Proc. IEEE AP-S Conf., vol. 2, Seattle, WA, June 1994, pp. 1366–1369. Authorized licensed use limited to: Universidad de Sevilla. Downloaded on July 06,2020 at 15:09:04 UTC from IEEE Xplore. Restrictions apply. 2192 IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 50, NO. 9, SEPTEMBER 2002 N. V. Shuley (S’79–M’85) was born in Melbourne, Australia, in 1951. He received the B.E.Sc. and M.E.Sc. degrees in electrical engineering from the University of New South Wales, N.S.W., Australia, in 1973 and 1975, respectively, and the Ph.D. degree in electrical engineering from the Chalmers University of Technology, Göteborg, Sweden, in 1985. From 1977 to 1978, he was an RF Engineer with Microwave Associates, Dunstable, U.K., where he was mainly concerned with the design and development of microwave solid-state power generation. From 1979 to 1985, he was a Research and Teaching Assistant and later a Post-Doctoral Scientist, until May 1988, with the Chalmers University of Technology. In May 1988, he joined the Department of Electrical Engineering, University of Queensland, Brisbane, Qld., Australia, as a Senior Lecturer. He has carried out consultancy work for Nordic industry since 1981 and is an accredited consultant for the European Space Agency. His current research interests include microwave antennas, dichroic surfaces, and analytical and computer-aided electromagnetics. R. R. Boix (M’96) was born in Melilla, Spain, in 1962. He received the Licenciado and Doctor degrees in physics from the University of Seville, Seville, Spain, in 1985 and 1990, respectively. Since 1985, he has been with the Electronics and Electromagnetics Department, University of Seville, where he became an Associate Professor in 1994. During the summers of 1991 and 1992, he was a Visiting Scholar with the Electrical Engineering Department, University of California at Los Angeles (UCLA). During the summer of 1996, he was a Visiting Scholar with the Electrical and Computer Engineering Department, Syracuse University, Syracuse, NY. His current research interest are focused on the analysis of the effects of complex substrates on the performance of planar passive microwave circuits, planar periodic transmission lines, frequency-selective surfaces, and printed circuit antennas. Dr. Boix is on the Editorial Board of the IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES. F. Medina (M’90–SM’01) was born in Puerto Real, Cádiz, Spain, in November 1960. He received the Licenciado and Doctor degrees from the University of Seville, Seville, Spain, in 1983 and 1987, respectively, both in physics. From 1986 to 1987, he spent the academic year withthe Laboratoire deMicroondesde l’ENSEEIHT, Toulouse, France. From 1985 to 1989, he was a ProfesorAyudante(AssistantProfessor)withtheDepartment of Electronics and Electromagnetism, University of Seville, and since 1990, he has been a Profesor Titular (Associate Professor) of electromagnetism. He is also currently Head of the Microwaves Group, University of Seville. His research interest includes analytical and numerical methods for guidance, resonant and radiating structures, passive planar circuits, and the influence on these circuits of anisotropic materials. Dr. Medina was a member of both the Technical Program Committee (TPC) of the 23rd European Microwave Conference, Madrid, Spain, 1993, and the TPC of ISRAMT’99, Malaga, Spain. He is on the Editorial Board of the IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES. He has been a reviewer for other IEEE and Institution of Electrical Engineers (IEE), U.K., publications. M. Horno was born in Torre del Camp, Jaén, Spain, and died in September 1998, in Seville, Spain. He received the Licenciado and Doctor degrees from the University of Seville, Seville, Spain, in 1969 and 1972, respectively, both in physics. In October 1969, he joined the Department of Electronics and Electromagnetism, University of Seville, where he became an Assistant Professor in 1970, Associate Professor in 1975, and Professor in 1986. His main fields of interest included boundary value problems in electromagnetic theory, wave propagation through anisotropic media, and microwave integrated circuits. During his final years, he was engaged in the analysis of planar transmission lines embedded in anisotropic materials, multiconductor transmission lines, and planar antennas. He was a member of the Electromagnetism Academy, Massachusetts Institute of Technology (MIT), Cambridge. Authorized licensed use limited to: Universidad de Sevilla. Downloaded on July 06,2020 at 15:09:04 UTC from IEEE Xplore. Restrictions apply.