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

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.


I. INTRODUCTION
B EAMFORMING 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/ 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 R with 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. Transverse cut of the parallel-plate channel in a normalized rotationally symmetric geodesic Luneburg lens of radius R with 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 A and 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 B on 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 A and B satisfies the variational equation where ds is 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 B along the lens aperture enables us to discretize the electric field distribution, from which we calculate the radiated field. Note that the solution 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 A k of the electric field at each k point 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: where A k is the amplitude at a given reference position, and dL k and dL k are the elementary arc lengths of the wavefront in the kth tube at the reference and final position, respectively, 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 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, ϕ 3 dB . The variable φ n can 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 k at a radial reference distance r is given by the following equation: where the angle dϕ k = (ϕ k+1 −ϕ k−1 )/2 is the linear angle subtended by the kth tube bounded by the k+1 and k−1 geodesics at the source position [see Fig. 4 Here, it should be noted that the specific value of r will not have any effect on the normalized radiation pattern, and hence it will be taken as r = 1 without loss of generality. As for dL k , this elementary arc length can be expressed as a function of the elementary distance dL k between the extreme points of the ray tube at the lens aperture where dL k is the magnitude of the elementary aperture vector being the coordinates of the target point in the aperture of the kth ray,n k is the unit vector associated with the normal direction to dL k , corresponding to the local normal to the aperture, andŝ k is the unit vector along the direction of arrival of the kth geodesic, corresponding to the local Poynting vector. The factor (ŝ k ·n k ) 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: 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 (R obs λ) radiated by the lens aperture (see Fig. 5): where σ k is the path length of the kth geodesic joining the source and the kth target point, r k = r krk is the vector from this target point to the observation point defined by (R obs , φ) (r k = |r k | andr k its associated unitary vector). E k (φ) 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: with p = 0.1, C = 1 if (ŝ k ·r k ) > 0 and C = 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 [n k ·ŝ k +n k ·r k ] 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: where a R and b R are the semiaxes of the ellipse, the parameters a and b being scale factors in the x and y directions, respectively. Here, the radius R is assumed to be 75 mm, or 7.5λ 0 at the design frequency f 0 = 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 x y-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 x y-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 x y-plane is sufficient for comparison.

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 w g = 8.64 mm, corresponding to the dimensions of the standard WR34 suitable for designs at K a -band. Full-wave numerical results are used to calibrate the model of the feed at 30 GHz, with ϕ 3 dB = 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.7 and b = 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 z as a function of the radius ρ is used: where h 0 , p, and q are 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 h 0 = 0.6321, p = 2 and q = 1.8. However, this theoretical profile has a 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 h 0 = 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 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 x y-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 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 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 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λ, h 0 = 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 K a -band [38], [39], [50]. A total of nine feeds is considered, spaced every 15 • in the x y-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.5 and 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 h 0 = 0.59, p = 2, and q = 1.6. In effect, only h 0 needs 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 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.7 and b = 1. In this case, the SLL requirement was relaxed to −15 dB, as convergence is more challenging due to the larger number  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 t 1 = 0.03, t 2 = 0.11, t 3 = 0.23, t 4 = 0.29, t 5 = 0.41, t 6 = 0.47, P 0 = 0.205, P 1 = −0.187, P 2 = 0.016, P 3 = −0.087, P 4 = −0.048, and P 5 = −0.107. The starting point for the PSO algorithm was a flat lens, corresponding to P i = 0 for i = 0, . . . , 5, with regular intervals t i+1 − t i = 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) 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 better results. Nevertheless, the main objective of the present work is achieved here, since the RT model results are in good agreement even in this case, which fully confirms that it can be used in a design procedure. In addition, scan losses are well approximated with the assumptions discussed above.
To complete the analysis of this first attempt at a compressed water drop lens design, the ray tracing and corresponding electric field distributions are reported in Fig. 14. Compared to other full-wave field distribution results reported in this article, it is clear that the strongly modulated profile affects more the propagation of the signal in the PPW section as the electric field strength shows more variation across the wavefront, similar to interference patterns, likely due to reflections at some of the bends. These effects are indeed not accounted for in the RT model and explain to some extent the small discrepancies observed. Despite the approximations, the results remain excellent and open new opportunities for the design of compact geodesic lenses.

V. CONCLUSION
In this manuscript, we have presented a fast and efficient RT approach for the analysis and design of non-rotationally symmetric geodesic lenses with arbitrary height profiles. This asymptotic method is based on a combination of geometrical optics and the scalar theory of diffraction. Geometrical optics is employed to compute the magnitude and phase of the electric field at the lens aperture. Phase information is extracted from the calculation of the ray trajectories as the shortest path between the feed and the discretized aperture via the appropriate meshing of the mean geodesic surface. The amplitude of the field is obtained by applying the ray-tube power conservation theory. The far-field pattern radiated by the lens antenna is finally calculated by means of the Kirchhoff-Fresnel diffraction formula. The radiated fields computed with this approach have been satisfactorily validated by comparison with the commercial simulator HFSS for both the trivial flat lens and the elliptically compressed Luneburg lens (this nonrotationally symmetric case is out of reach for theoretical methods previously reported in the literature). The computational CPU and memory cost of the proposed RT approach are found to be very advantageous, with reduction factors above 40 and 2 respectively, when compared to simulations using commercial software. This feature is particularly relevant when the RT method is exploited for the design of compact high-performance lenses, in particular, a modified Luneburg profile which largely improves the off-axis performance and a compressed water drop geodesic lens design with promising results while maintaining the mechanical simplicity of the original geodesic lenses. Future work will explore the use of this RT approach applied to more general geodesic surfaces to improve, for example, the scanning performance of highly compressed designs.

APPENDIX
For the optimizations carried out in this article, we have used the Python module "PySwarms," which implements a PSO. The target function that we chose to optimize with this PSO is defined as follows: where E(φ; φ s ) is the radiated electric field [see (6)] produced by a given source located at φ s and Mask(φ; φ s ) is defined in terms of the rectangular function [rect(x) has the change for |x| = 1/2] as follows: where SLL is the required SLL (in dB),φ s the angle at which the radiated field is maximum, and φ SLL the angle of approximately the first null. With this definition of the target function, we intend to carry out a "global" optimization that takes into account the radiation patterns produced by sources at different angle positions. The optimization goal is to minimize the target function and the process is stopped after 90 iterations. This was found sufficient to ensure acceptable convergence of the PSO algorithm for the cases considered in this manuscript.