Ray-Tracing Model for Generalized Geodesic-Lens Multiple-Beam Antennas

2640 IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 71, NO. 3, MARCH 2023 Ray-Tracing Model for Generalized Geodesic-Lens Multiple-Beam Antennas Qingbi Liao, Student Member, IEEE,NelsonJ.G.Fonseca ,Senior Member, IEEE, Miguel Camacho ,Member, IEEE, Ángel Palomares-Caballero ,FranciscoMesa ,Fellow, IEEE, and Oscar Quevedo-Teruel ,Fellow, IEEE Abstract—Geodesic lenses are a compelling alternative to traditional planar dielectric lens antennas, as they are low loss and can be manufactured with a simple mechanical design. However, a general approach for the design and analysis of more advanced geodesic-lens antennas has been elusive, limiting the available tools to rotationally symmetric surfaces. In this article, we present a fast and efficient implementation built on geometrical optics and scalar diffraction theory. A numerical calculation of the shortest ray path (geodesic) using an open-source library helps quantify the phase of the electric field in the lens aperture, while the amplitude is evaluated by applying ray-tube power conservation theory. The Kirchhoff-Fresnel diffraction formula is then employed to compute the far field of the lens antenna. This approach is validated by comparing the radiation patterns of a transversely compressed geodesic Luneburg lens (elliptical base instead of circular) with the ones computed using commercial full-wave simulators, demonstrating a substantial reduction in computational resources. The proposed method is then used in combination with an optimization procedure to study possible compact alternatives of the geodesic Luneburg lens with size reduction in both the transverse and vertical directions. Index Terms—Geodesic lenses, lens antennas, non-euclidean transformation optics, parallel plate waveguides (PPWs), ray tracing. Manuscript received 11 February 2022; revised 19 November 2022; accepted 20 December 2022. Date of publication 9 January 2023; date of current version 6 March 2023. The work of Qingbi Liao and Oscar Quevedo-Teruel was supported in part by the Office of Naval Research (ONR) under Grant N62909-20-1-2040. The work of Miguel Camacho and Francisco Mesa was supported in part by the MCIN/AEI /10.13039/501100011033 under Grant PID2020-116739GB-I00. (Corresponding author: Oscar Quevedo-Teruel.) Qingbi Liao is with the Filter and Antenna System, Ericsson AB, 164 83 Kista, Stockholm, Sweden (e-mail: [email protected]). Nelson J. G. Fonseca is with the Antennas and Sub-Millimeter Waves Section, European Space Agency, 2200 AG Noordwijk, The Netherlands (e-mail: [email protected]). Miguel Camacho is with the Department of Electronics and Electromagnetism, School of Physics, Universidad de Sevilla, 41012 Seville, Spain (e-mail: [email protected]). Ángel Palomares-Caballero is with the Department of Signal Theory, Telematics and Communications, University of Granada (UGR), 18071 Granada, Spain (e-mail: [email protected]). Francisco Mesa is with the Department of Applied Physics 1, ETS Ingeniería Informática, Universidad de Sevilla, 41012 Seville, Spain (e-mail: [email protected]). Oscar Quevedo-Teruel is with the Division of Electromagnetic Engineering, School of Electrical Engineering and Computer Science, KTH Royal Institute of Technology, SE 100 44 Stockholm, Sweden (e-mail: os[email protected]). Color versions of one or more figures in this article are available at https://doi.org/10.1109/TAP.2022.3233643. Digital Object Identifier 10.1109/TAP.2022.3233643 I. INTRODUCTION BEAMFORMING devices are essential to produce high-gain antenna solutions of interest for a wide range of applications, including wireless communications, automotive and surveillance radars, microwave imaging and sensing, etc. Specifically, beamforming networks, a discrete form of beamforming devices based on interconnected elementary components, are attractive when high integration is required [1]. Solutions for multiple-beam antenna design include the Blass matrix [2], [3], the Butler matrix [4], [5], [6] and the Nolen matrix [7], [8], [9], [10]. However, their complexity tends to increase exponentially with the size of the antenna aperture, and losses may be prohibitive at very high frequencies. Therefore, quasi-optical solutions are generally preferred to achieve higher gain values [11]. They are also particularly adapted to very high frequencies, including terahertz, thanks to their generally simpler designs, more tolerant to manufacturing errors. They are referred to as quasi-optical solutions because the designs are inspired by optical principles and concepts. Well-known solutions include the pillbox antenna [12], [13], [14], the Rotman lens [15], [16], [17], [18], homogeneous dielectric lenses [19], [20], [21], and shaped parallel-plate lenses [22]. One quasi-optical solution of interest is the Luneburg lens [23]. Its intrinsic symmetry provides in theory perfect collimation of the beams in any direction. This is achieved thanks to a graded-index distribution with spherical or rotational symmetry for a volumetric or planar design, respectively. A range of focusing properties can be engineered by adapting the graded-index profile [24], [25], [26], [27]. Several designs have been reported in the recent literature with applications in the microwave, terahertz, and optical domains [28], [29], [30], [31], [32], [33]. The main challenge for its implementation is the graded-index profile, which can be very lossy at millimeter-waves and frequencies above. In this sense, geodesic lenses are an interesting alternative for applications that require a planar beamformer design [34], [35], [36]. Geodesic lenses make use of the out-of-plane dimension in a parallel plate waveguide (PPW) configuration to emulate the path length differences introduced by the refractive index in planar Luneburg lenses. Their simple design is particularly suitable for millimeter-wave applications. The equivalence This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ LIAO et al.: RAY-TRACING MODEL FOR GENERALIZED GEODESIC-LENS MULTIPLE-BEAM ANTENNAS 2641 between graded-index and geodesic lenses may be seen as a first implementation of transformation optics [34], [37] and was recently used to develop a ray-tracing (RT) model of a modulated geodesic lens based on its equivalent planar lens [38], [39]. This RT approach was later adapted directly to geodesic-lens antennas in [40]. However, this approach is restricted to rotationally symmetric lenses, since the rays are calculated analytically based on the conservation of the angular momentum, as reported in [41] and [37]. Non-rotationally symmetric graded-index lenses show strong potential for reduced footprint and could lead to more degrees of freedom to adapt the radiation pattern for different feed positions [42], [43], [44], [45]. Some works have also been reported on non-rotationally symmetric geodesic lenses in the optical domain [46], [47]. In this article, we present a generalized and computationally efficient approach based on open-source RT tools to compute the radiation patterns of multiple beam geodesic-lens antennas. We compare the results obtained for a set of non-rotationally symmetric surfaces with a well-established commercial full-wave simulator and demonstrate the capabilities for the analysis of compressed geodesic Luneburg lenses, both in the transversal and vertical directions, applied to the design of multiple-beam antennas. The article is organized as follows. Section II discusses the main theoretical aspects of the RT approach used in this work. In Section III, a detailed study of the numerical implementation of the method is carried out. The accuracy and computational cost of the method are also evaluated by comparison with full-wave commercial software. Section IV illustrates the potential of the proposed technique with the design of two lens antennas: an elliptically compressed Luneburg-type geodesic lens and a modified water drop lens. Finally, some conclusions are drawn in Section V. II. RAY-TRACING MODEL From early on [34], [48], [49], it was proposed that wave propagation in geodesic lenses took the form of a transverse electromagnetic (TEM) mode that travels through the free-space channel formed between two parallel metallic surfaces, as illustrated in Fig. 1 for the case of a Rinehart–Luneburg lens of radius Rwith a transition at its end to reduce reflections. This assumes that the height of the PPW, t, is small enough compared to the wavelength to support only the fundamental TEM mode. As mentioned in Section I, the simulation of the wave propagation through the resulting parallel-plate channel and its further radiation is often performed using commercial full-wave solvers [50]. Since these simulations require long computation times and high memory resources, RT-based techniques have been considered, demonstrating to be both time-efficient and reasonably accurate to enable optimization processes [51], [52]. However, the techniques reported are limited to rotationally symmetric solutions [38], [40] or to specific lens designs [53]. Here, we propose an asymptotic numerical method using RT to compute the electric field in the aperture of the generalized geodesic lens, combined with scalar diffraction theory, to calculate the resulting electric field in the far-field region. Fig. 1. Transverse cut of the parallel-plate channel in a normalized rotationally symmetric geodesic Luneburg lens of radius Rwith a transition for reducing reflections (ρis the radial distance in the xy-plane). The black solid lines illustrated the PPW. Fig. 2. Schematic representation of the ray trajectory in a geodesic lens. z(ρ, φ) defines the geodesic surface, and the red curve represents the geodesic (shortest path) between points Aand B(the lens is expanded along the z -axis for readability). In geometrical optics, rays are defined as the trajectories orthogonal to the light’s wavefronts [41], [54]. To compute these rays in a geodesic lens, we define a curved surface as the mean surface within the PPW (see the so-called geodesic mean curve in Fig. 1). The phase delay that the wave will accumulate when traveling from a point A, on the edge of the surface and corresponding to the feed position as shown in Fig. 2, to any point Bon the aperture of the lens can be obtained applying Fermat’s principle [41], [54]. Under the assumption that only a TEM mode propagates in the PPW, the shortest path between points Aand Bsatisfies the variational equation δB A ds=0(1) where dsis the elementary arc length of the ray on a homogeneous geodesic surface (the refractive index is constant). The solution to (1) provides the path followed by the ray traveling between these two points. Sweeping Balong the lens aperture enables us to discretize the electric field distribution, from which we calculate the radiated field. Note that the solution 2642 IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 71, NO. 3, MARCH 2023 Fig. 3. Mesh model of a non-rotationally symmetric geodesic lens. The geodesic surface is obtained from a rotationally symmetric design of radius R, scaled by factors aand balong the xand ydirections, respectively. The surface includes a toroidal transition region to ensure the continuity of the slope. to (1) is generally unique. In some cases, however, there may be multiple solutions, which translates into rays intersecting in the region under study. In the case of a focal point inside the lens or at its edge, as in the Maxwell fish eye lens, the rays are all converging to a single point, corresponding to an infinite number of solutions to (1). These particular configurations may affect the accuracy of the method described in this article. The intersection of rays generally does not occur in lenses designed to collimate rays, which is of interest here for multiple-beam antenna design. For the solution of the variational equation (1), we can resort to well-established numerical procedures for the calculation of geodesic paths on an arbitrary polyhedral surface given in terms of a mesh [55], [56], [57], [58]. For example, in Fig. 3, we illustrate the discretized mesh of a non-rotationally symmetric geodesic lens which is used for the computation of geodesic paths. The length of the geodesic trajectories obtained is used to evaluate the corresponding phase of the field in the lens aperture. In this work, the shortest distance for each ray traveling through the lens is evaluated using a specialized library in Python programming language called “gdist” [59], based on [60]. The authors have modified this library to also provide the full path of the different rays along the meshed surface. Specifically, these changes in the C++ to Python interfaces enable us to extract the so-called “angle of departure” and “angle of arrival” of the rays as they depart from the source point and arrive at the target points on the aperture lens, respectively. This information is needed to obtain the amplitude Akof the electric field at each kpoint of the aperture using ray-tube power conservation theory [41], [53], [54]. Following the notation in Fig. 4, this amplitude is proposed to be computed as follows: Ak=A kdL k/dLk(2) where A kis the amplitude at a given reference position, and dL kand dLkare the elementary arc lengths of the wavefront in the kth tube at the reference and final position, respectively, Fig. 4. (a) Projections in the xy-plane of the geodesics on a rotationally symmetric Rinehart–Luneburg lens with a point source located at (R,φ s). (b) kth ray tube and geometrical variables involved in the evaluation of the field amplitude. assuming that the PPW has a constant height, so this parameter does not have an influence on the field strength. Taking the reference position very close to the source, it is assumed that the power distribution is independent of the angle of departure of the different geodesics for an isotropic point source excitation, that is, A k=1. In the case of a directive source, the power distribution does depend on the angle at which each ray is emitted from the feed. A solution often used to feed geodesic lenses is the open-ended waveguide, with an amplitude distribution that can be approximated by a Gaussian function defined as A k=10−(3ξ2/20),where ξ=(ϕk−φn)/ϕ3dB , given in terms of the ratio between the difference between the angles of departure, ϕk, and the angle of the normal to the waveguide aperture, φn, with the half-power beamwidth angle, ϕ3dB .Thevariableφncan be used as a design parameter to control the amplitude distribution in the lens aperture. The phase center also needs to be precisely defined to reduce errors when comparing the RT model with full-wave results, and may be adjusted in the model by adequate positioning of the source point. The elementary arc length dL kat a radial reference distance ris given by the following equation: dL k=rdϕk(3) where the angle dϕk=(ϕk+1−ϕk−1)/2 is the linear angle subtended by the kth tube bounded by the k+1andk−1 geodesics at the source position [see Fig. 4(b)]. Here, it should be noted that the specific value of rwill not have any effect on the normali z ed radiation pattern, and hence it will be taken as LIAO et al.: RAY-TRACING MODEL FOR GENERALIZED GEODESIC-LENS MULTIPLE-BEAM ANTENNAS 2643 Fig. 5. Geometric parameters associated with the numerical evaluation of the far-field pattern (Robs λ). r=1 without loss of generality. As for dLk, this elementary arc length can be expressed as a function of the elementary distance dLkbetween the extreme points of the ray tube at the lens aperture dLk=dLk(ˆ sk·ˆ nk)(4) where dLkis the magnitude of the elementary aperture vector dLk=[(xk−1−xk+1)ˆ x+(yk−1−yk+1)ˆ y]/2, (xk,yk)being the coordinates of the target point in the aperture of the kth ray, ˆ nk is the unit vector associated with the normal direction to dLk, corresponding to the local normal to the aperture, and ˆ skis the unit vector along the direction of arrival of the kth geodesic, corresponding to the local Poynting vector. The factor (ˆ sk·ˆ nk) is introduced to take into account that the aperture of the lens does not necessarily coincide with the wavefront. The amplitude in (2) is finally obtained as follows: Ak=A kdϕk dLk(ˆ sk·ˆ nk).(5) Knowing the electric field distribution in the aperture, we can apply Kirchhoff’s scalar diffraction theory [41] to obtain the following approximation of the far-field (Robs λ) radiated by the lens aperture (see Fig. 5): E(φ)∝ k Ek(φ)Ak e−jk0(rk+σk) rk [ˆ nk·ˆ sk+ˆ nk·ˆ rk]dLk(6) where σkis the path length of the kth geodesic joining the source and the kth target point, rk=rkˆ rkis the vector from this target point to the observation point defined by (Robs,φ) (rk=|rk|and ˆ rkits associated unitary vector). Ek(φ) is an angle-dependent factor, which is introduced to avoid the field discontinuities that a binary formulation would cause at the transition between the forward and backward directions of propagation, approximated as follows: Ek(φ)=C|ˆ sk·ˆ rk|p(7) with p=0.1, C=1if(ˆ sk·ˆ rk)>0andC=0, otherwise. A similar factor was also introduced in [53] for point source patterns in a circular array approximation of the aperture. It should be noted that the factor [ˆ nk·ˆ sk+ˆ nk·ˆ rk]in (6) reduces to (1+cos φ) for the case of a circular aperture; corresponding to the obliquity factor in the Huygens-Fresnel principle formulation used in [38], [39], and [40]. The proposed numerical method may be seen as a simplified implementation of physical optics (PO) adapted to the unique properties of geodesic lenses, and more generally PPW beamformers, enabling to treat the electric field as a scalar wave. III. MODEL VALIDATION To validate the implementation of the RT approach presented above, we consider two cases: 1) a flat lens, or a flat PPW corresponding to a trivial invisible lens [27], with an elliptical outer rim and 2) an elliptical geodesic lens obtained from the compression of a rotationally symmetric one of radius R. In both cases, the outer rim is given by the following equation: x aR2 +y bR2 =1(8) where aR and bR are the semiaxes of the ellipse, the parameters aand bbeing scale factors in the xand ydirections, respectively. Here, the radius Ris assumed to be 75 mm, or 7.5λ0at the design frequency f0=30 GHz, corresponding to the value also used to validate the equivalent planar lens RT model in [38]. The radiating aperture corresponds to the half contour on the positive side of the x-axis and the rays hitting the remaining part of the outer rim are neglected, corresponding to the typical aperture size implemented in multiple feed lens designs [38], [50]. The numerical results obtained at 30 GHz with the RT model are validated by comparison with full-wave simulation results using the commercial software HFSS, which has already demonstrated multiple times excellent agreement with the experimental results for various PPW lens antenna designs [38], [39], [50], [53] and is therefore used here as a reference for the newly developed tool. In this section, a thin PPW section, t=1 mm, is considered to minimize the numerical errors resulting from a volumetric design as the corresponding RT model is a surface one. Absorbing boundary conditions are also applied in the HFSS model to the part of the outer rim with x<0, delimiting the radiating aperture to the positive half contour. Radiation patterns are only evaluated in the xy-plane and thus the directivity cannot be estimated with the RT model. For this reason, the electric far-field patterns are normalized. Interestingly, in the case of focused lens designs, scan losses may be evaluated under the assumption that the variations of the far-field patterns are predominantly defined by the change in patterns in the xy-plane. This hypothesis has already demonstrated good results in [39] and also proves to be a good approximation in the cases studied here. The RT model may be extended in the future to account for the antenna flare and define the complete pattern, thus enabling to evaluate also directivity and gain. Here, we focus on the beamforming capabilities of the lenses and the analysis in the xy-plane is sufficient for comparison. 2644 IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 71, NO. 3, MARCH 2023 A. Elliptical Flat Lens First, a canonical flat lens case is considered. It provides a separate validation of the far-field computation, as the ray tracing is trivial in this particular case. It is also used to demonstrate the validity of the ray-tube power normalization by comparing two configurations with sufficiently different source excitations: an isotropic source and an open-ended waveguide source, specifically a rectangular waveguide with width wg=8.64 mm, corresponding to the dimensions of the standard WR34 suitable for designs at Ka-band. Full-wave numerical results are used to calibrate the model of the feed at 30 GHz, with ϕ3dB =32.5◦and a phase center displaced 0.5 mm inside the waveguide with reference to the aperture. This particular waveguide excitation will also be used in other examples in this work. The flat lens is compressed along the x-axis with scale factors a=0.7andb=1. The source is located on the x-axis, that is, φs=0◦. Note that in these particular cases, the boundary condition on the negative side of the x-axis is important for good validation. The RT model assumes only direct “line of sight” contributions from the feed to the aperture, and secondary reflections on the side walls are neglected. This is generally acceptable in geodesic lenses that focus the electric field toward the aperture but was found in these flat lens cases to degrade the comparison between our proposed model and the full-wave analysis results. A similar issue occurs in Rotman lenses [15], which also comprise a flat PPW section, and is typically solved by adapting the shape of the side walls to suppress direct reflections [61]. The HFSS results reported in this section were obtained by implementing a similar corrective action. Radiation patterns are plotted in Fig. 6(a) and (b). The results are in good agreement; in particular, the overall shape of the beam is well predicted, including the beamwidth and beam slope. Some small discrepancies are found on the beam slopes, with the HFSS results showing some residual ripples not present in the RT results, which are the consequence of the boundary effect described above. However, the agreement remains largely acceptable. These results validate the normalization of the amplitude across the lens aperture and the far-field computation approach. B. Elliptically-Compressed Rinehart–Luneburg Lens The second case of validation considers geodesic lenses. First, a reference geodesic lens corresponding to a non-compressed rotationally symmetric design is analyzed. To simplify the implementation of the model, the following analytical expression defining the height zas a function of the radius ρis used: z R=h01−ρ Rp1/q(9) where h0,p,andqare superellipse design parameters set to adjust the profile shape. As demonstrated in [38], this simple expression is a good approximation to the theoretical Rinehart–Luneburg profile with parameters h0=0.6321, p=2andq=1.8. However, this theoretical profile has a Fig. 6. Radiation patterns of a flat elliptical lens at 30 GHz with design parameters a=0.7,b=1, and R=7.5λ, fed with: (a) isotropic source and (b) rectangular open-ended waveguide. slope discontinuity at the rim of the lens, which is not suitable for practical implementation. Thus, instead, we use as a reference the profile with a toroidal bend as optimized in [38]. The corresponding parameters in this case are h0=0.56, p=2, and q=1.6 for a toroidal bend with a radius equal to 5 mm, or half-a-wavelength at the design frequency. Here, an open-ended waveguide source with two different locations is also used, φs=0◦and φs=45◦, to validate the performance of the RT model for scanned beams. The waveguide feeds are oriented so that φn=φs, corresponding to feeds that point toward the center of the rotationally symmetric lens. The boundary condition on the negative x-axis side of the HFSS model is applied here directly on the rim of the lens, as secondary reflections are found to have a marginal impact on the far field. Being the first study case with a shaped surface, intermediate results are provided here to validate the ray tracing implemented with the proposed approach. The normalized electric field in the aperture is provided in Fig. 7, both in amplitude and phase. This quantity is generally more sensitive LIAO et al.: RAY-TRACING MODEL FOR GENERALIZED GEODESIC-LENS MULTIPLE-BEAM ANTENNAS 2645 Fig. 7. (a) Amplitude and (b) phase of the aperture near-field in a rotationally symmetric Rinehart–Luneburg lens at 30 GHz with design parameters R= 7.5λ, h0=0.56,p=2, and q=1.6, including a toroidal transition of radius 5 mm at the lens rim. The lens is fed by a waveguide source at φs=0◦(the ideal phase corresponds to a reference plane wave from the circular aperture). to numerical errors resulting from the mesh and secondary effects, such as boundary conditions and discontinuities, both in HFSS and in our RT model, in line with previously reported results for other types of lenses [53]. Nevertheless, the amplitude shows good agreement in the overall variation, which is essential for a good prediction of the sidelobe level (SLL), and the phase response is well aligned with the ideal phase required to produce a plane wave, as expected for this focused design. The evaluation of the far field through the integration of this complex quantity typically reduces the effect of numerical errors. This is evidenced by the results reported in Fig. 8, including ray tracing and electric field distribution for both feed positions and compared far-field patterns in the xy-plane. The good collimation of the rays can be appreciated in Fig. 8(a) as the rays are coming out parallel on the aperture side. It is also clear from this representation that part of the geodesic surface is not covered by rays for the scanned feed position. This illustrates the energy lost by spillover. With the electric field distribution reported in Fig. 8(b), the plane wave Fig. 8. (a) Ray tracing, (b) electric field distribution, and (c) radiation patterns of a rotationally symmetric Rinehart–Luneburg lens at 30 GHz with design parameters R=7.5λ, h0=0.56,p=2, and q=1.6, including a toroidal transition of radius 5 mm at the lens rim. The lens is fed by a waveguide source at φs=0◦and φs=45◦. that excites the lens is clearly visible in both feed position cases. The patterns presented in Fig. 8(c) show excellent agreement in the shape of the main beam and the first few side lobes. These results also indicate that the scan losses due to the finite aperture size are still marginal at 45◦, in line with previously reported results for this reference lens [38]. The RT tool predicts scan losses of 0.36 dB, while HFSS provides a value of 0.67 dB at 45◦, indicating the differences are marginal and acceptable. The final step in the validation of the proposed RT model is a compressed geodesic lens with the same profile as the reference geodesic lens discussed above and scale factors a=0.5, b=1. Interestingly, this configuration provides a geodesic lens with dimensions equivalent to the half-Luneburg lens discussed in [39]. Here, the waveguide feeds are aligned with 2646 IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 71, NO. 3, MARCH 2023 Fig. 9. (a) Ray tracing, (b) electric field distribution, and (c) radiation patterns of a compressed Rinehart–Luneburg lens at 30 GHz with design parameters a=0.5,b=1,R=7.5λ, h0=0.56,p=2, and q=1.6, including a toroidal transition of radius 5 mm at the lens rim. The lens is fed by a waveguide source at φs=0◦and φs=45◦. the normal to the lens contour, as this configuration is better suited for the multiple-feed designs discussed in Section IV. The boundary condition in HFSS is also applied to the rim of the lens, as this configuration is more representative of the multiple-beam antenna designs discussed next. Similar numerical results are provided in Fig. 9, reporting the ray tracing and electric field distribution for the same two feed positions and comparing far-field patterns. Despite the strong compression factor, the elliptical geodesic lens still demonstrates reasonably good collimation of the electric field, as evidenced by the parallel rays in the aperture of the lens in Fig. 9(a), and the parallel wavefronts in Fig. 9(b). This is further validated with the numerical results compared in Fig. 9(c), where the RT model is again in excellent agreement with the full-wave simulation results. The main beam is particularly well predicted, while some small discrepancies appear in the side lobes attributed to the boundary conditions in HFSS. The levels are still in good agreement. Focusing properties are mostly preserved when the lens is compressed. However, the shoulders at about −13 dB below the maximum directivity of the on-axis beam (φs=0◦) indicate that it introduces some phase aberrations, well accounted for by the RT model. The phase aberrations are quantified in HFSS by evaluating the phase of the electric field along a line orthogonal to the direction of propagation for the waveguide source at φs=0◦and tangent to the aperture of the lens. The results are reported in Fig. 10 for both the reference and compressed geodesic lenses. In the central area, corresponding to the strongest amplitude values, the phase of the reference geodesic lens was mostly flat, in line with the expected plane wave response. In the case of the compressed lens, quadratic aberrations are clearly visible, which explains the shoulders observed in the far-field patterns. Due to the high compression factor, the pointing direction of the scanned beam deviates from the angular feed position value, φs. It remains, nevertheless, well predicted by the proposed method. Note that this beam deviation is a result of the phase distribution and not the amplitude distribution; thus the pointing angle of the waveguide feed, φn, would have only a secondary impact on this parameter. It is also noted that the scan losses are well-predicted in this configuration. As expected, they are slightly higher than in the rotationally symmetric case, as a direct consequence of the reduced projected aperture. The value predicted by the RT model is 0.96 dB, compared to a value of 1.04 dB in HFSS, in both cases for a maximum at −34◦. In terms of computational effort, the RT implementation is found to achieve on average a CPU time reduction factor of over 40 and a RAM usage reduction factor above 2 when compared to the HFSS simulations with a similar convergence level in far-field patterns. This comparative test was carried out on an Intel Core i7 laptop with 8 GB of RAM without making use of any parallelization scheme. RT methods previously reported in [38], [39], and [40] are generally faster, but restricted to rotationally symmetric lenses. The approach described here provides a generalized formulation applicable to any geodesic surface while preserving the benefits in terms of computational effort compared to full-wave methods. These results fully validate the proposed modeling approach for the design of non-rotationally symmetry geodesic lenses and confirm its potential for optimization purposes. Section IV illustrates the use of this tool for the design of multiple-beam antennas. IV. MULTIPLE-BEAM ANTENNA DESIGN EXAMPLES In this section, practical multiple-beam antenna designs are demonstrated using the developed RT model together with an optimization routine. Two solutions are discussed: one based on a conventional profile and one using a modulated profile, the so-called water drop lens profile [38]. If optimization was carried out using commercial simulators, the required computational time would become prohibitive given the large number of simulations that are usually required by the optimization algorithms, as well as the long time required for each simulation. The RT model provides significant benefits with its demonstrated accuracy and its much reduced computational requirements, both in terms of memory and CPU. In this work, we propose to search for the convenient height profile of the lens by means of a procedure that combines particle swarm optimization (PSO) with the RT approach. Other optimization LIAO et al.: RAY-TRACING MODEL FOR GENERALIZED GEODESIC-LENS MULTIPLE-BEAM ANTENNAS 2647 Fig. 10. Near-field wavefront phase for a reference (a=1) and a compressed (a=0.5) rotationally symmetric Rinehart–Luneburg lens at 30 GHz with design parameters R=7.5λ, h0=0.56,p=2, and q=1.6, including a torodial transition of radius 5 mm at the lens rim. The lens is fed by a waveguide source at φs=0◦. techniques could be used. The choice of the best optimization algorithm for this specific case study is beyond the scope of this work. The information about the employed PSO, as well as the fitness function, can be found in the Appendix. For a more realistic design, multiple feeds are implemented in the HFSS model and radiating boundary conditions are only applied to the aperture, meaning that the walls in between feeds remain fully metallic. This is expected to introduce some discrepancies between the RT and HFSS results larger than discussed in Section III addressing the validation of the model. In addition, the PPW thickness is increased to t=2 mm, in line with previously reported multiple-beam antenna prototypes at Ka-band [38], [39], [50]. A total of nine feeds is considered, spaced every 15◦in the xy-plane and symmetrically arranged around the negative x-axis. Due to obvious design symmetries, results are only reported for φs=0◦,15◦,30◦,45◦, and 60◦. A. Multiple-Beam Elliptical Rinehart–Luneburg Lens The first design is based on a compressed Rinehart–Luneburg lens profile similar to the one described in Section III. Compressing a rotationally symmetric design without further adjustments results in phase aberrations that affect the shape of the main beam. The focusing properties are nevertheless preserved, indicating that the profile is not far from an optimum when using scale factors a=0.5and b=1. The advantage of this configuration is that it has a very limited number of parameters, and the optimization is therefore straightforward. In this particular case, it was decided to limit the optimization to the boresight beam with a target SLL of −20 dB. Based on the aperture size, the sidelobe mask is further defined with φSLL =10◦. The resulting parameters are h0=0.59, p=2, and q=1.6. In effect, only h0needs to be adjusted. All remaining parameters, including the lens radius and the toroidal bend radius, were kept unchanged. The resulting profile is illustrated in Fig. 11, where it is compared to the Fig. 11. Height profiles of an optimized Luneburg-type elliptical lens with a toroidal transition of radius 5 mm at the lens rim and of a modulated water drop lens (the ideal Rinehart–Lunerburg lens profile is provided as reference). theoretical Rinehart–Luneburg profile. Note that this profile representation is before compression along the x-axis. Radiation patterns obtained with the RT model (solid lines) and the full-wave model (dashed lines) are compared in Fig. 12. To facilitate visualization, the same color is used for the corresponding beams. Thus, excellent agreement is confirmed when the two curves cannot be distinguished. This is particularly the case for the beams around boresight, where a very good agreement is found in the main beam, with differences appearing only in the low SLLs, in line with the agreement reported in Section III. This design achieved low SLL for the on-axis beam, with the main beam improvement also confirmed with the full-wave model including all beam ports. These results indicate that secondary effects, such as scattering from the structure and diffraction on the edge of the aperture, are indeed negligible and the RT model is relevant for more realistic multiple-beam antenna designs. Discrepancies between the HFSS and the RT model results are more visible on the most scanned beams, but the agreement remains fairly acceptable since the RT results appear more conservative. B. Multiple-Beam Elliptical Water Drop Lens While the previous design demonstrates the potential of non-rotationally symmetric geodesic lenses for in-plane size reduction, the height of these lenses remains relatively unchanged and may be too high for some applications requiring further integration. Here, a modulated lens profile is implemented to significantly reduce the height of the lens, while also reducing the in-plane dimensions. The geodesic surface is obtained by compressing a rotationally symmetric water drop lens using a spline profile as described in [38]. The optimized profile before compression is illustrated in Fig. 11, where it is compared with the reference profile. A height reduction by a factor of 3 is achieved. The in-plane dimensions are compressed using the scale factors a=0.7andb=1. In this case, the SLL requirement was relaxed to −15 dB, as convergence is more challenging due to the larger number 2648 IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 71, NO. 3, MARCH 2023 Fig. 12. Radiation patterns of an elliptical geodesic lens at 30 GHz with design parameters a=0.5,b=1,R=7.5λ, h0=0.59,p=2, and q=1.6, including a toroidal transition of radius 5 mm at the lens rim and nine feeds spaced every 15◦within the range φs=±60◦(only five beams are reported due to symmetries). Fig. 13. Radiation patterns of an elliptical modulated water drop lens at 30 GHz with design parameters a=0.7,b=1, R=7.5λ, and nine feeds spaced every 15◦within the range φs=±60◦(only five beams are reported due to symmetries). of variables and the more complex profile. The sidelobe mask is defined with φSLL =10◦. Using the same definition of the spline profile as in [38], corresponding to a lens diameter with enforced symmetry to have a horizontal tangent at the center of the profile, the optimized values are t1=0.03, t2=0.11, t3=0.23, t4=0.29, t5=0.41, t6=0.47, P0=0.205, P1=−0.187, P2=0.016, P3=−0.087, P4=−0.048, and P5=−0.107. The starting point for the PSO algorithm was a flat lens, corresponding to Pi=0fori=0,...,5, with regular intervals ti+1−ti=1/13 for i=0,...,5. This case corresponds to the radiation pattern reported in Fig. 6(b). Profile parameters are optimized so that the lens height remains within ±0.1R. The radiation patterns obtained with the RT model (solid lines) and with the full-wave model (dashed lines) Fig. 14. Ray tracing (left) and electric field distribution (right) of an elliptical modulated water drop lens at 30 GHz with design parameters a=0.7,b=1, R=7.5λ, and with φsvarying every 15◦from 0◦to 60◦(top to bottom). are compared in Fig. 13. Corresponding beams are also represented here with the same color. Despite the more complex profile, good agreement is still demonstrated between the RT model and the full-wave model, particularly in the main beam for beams around boresight. SLLs are a bit higher in this case but remain acceptable considering the significant size reduction of the geodesic lens. Larger discrepancies are visible in the most scanned beam (φs=60◦), indicating that this lens has a lower scanning range as a consequence of the drastic compression. Further optimization work would be required to investigate the potential of these lenses. Possibly, alternative profiles and compression techniques may lead to