GO Analysis of GRIN Lens Antennas by Combining in a Single ODE, Field and Wavefront-Curvature Transport to the Ray Tracing

A novel very efficient algorithm based on geometrical optics (GO) is presented for the analysis of graded index (GRIN) lens antennas, namely, dielectric inhomogeneous lenses with a 3-D arbitrary varying refractive index. A family of curved ray paths are traced starting from a set of points defined on the lens input interface, which is illuminated by a feed antenna, up to a corresponding set of points on the output interface, i.e., the lens radiating aperture. The ray tracing is numerically performed in combination with the field transportation along the ray by exploiting an additional wavefront-curvature transport equation, thus providing a single independent vector ordinary differential equation (ODE) for each ray. This scheme allows a complete parallelization of the algorithm by assigning any ray to a different thread. Crossing the lens input–output interfaces at the two end points of the curved ray is accomplished by augmenting the refraction Snell’s law and the Fresnel transmission coefficients with a novel compact closed-form dyadic formula for the transmitted wavefront curvature. The ODE output provides the field aperture distribution on the output lens interface, permitting the far-field calculation as a radiation integral. To this end, an ad hoc quadrature rule is exploited, which requires an aperture sampling rate corresponding to the slowly varying part of the integrand, while the rapid phase variation is accounted for analytically, thus resulting in an extremely efficient and frequency-independent scheme. The effectiveness and accuracy of the algorithm are shown by resorting to a couple of well-known analytical benchmarks and to two more general examples with real-life antennas: a spaceborne weather radar lens antenna and a feed antenna with a zero-focal lens.

capacity and reduced latency.This requirement has spurred the industry to transition to millimeter-wave (mm-wave) frequencies, which offer broader bandwidth solutions.Nonetheless, a significant obstacle arises in the form of substantial freespace path losses.To tackle this issue, there is a demand for antenna technologies capable of delivering high directivity and supporting multiple beams aimed at various users.
Graded index (GRIN) lens-based antennas can emerge as a competitive technology when compared to other conventional technologies such as homogeneous lens antennas, reflector antennas, and phased arrays.This is attributed to their rapid and cost-effective development using additive manufacturing, reduced complexity, and the greater flexibility they offer in the design process.This flexibility in design allows for the creation of slimmer and more compact implementations, even though it may involve the use of materials with higher maximum permittivity.Beam scanning can also be accomplished by altering the position or orientation of the source relative to the lens [1] or by customizing the lens distribution for multifocal capabilities.The resurgence of interest in GRIN lenses is linked to the newfound ability to manufacture them using 3-D printers [2], [3], [4], [5], [6], [7], [8].The two most commonly used GRIN lenses are the spherical Luneburg lens (LL) [9] and the Maxwell fish-eye (MFE) lens [10], for which the refractive index is available in an analytical closed form.Nevertheless, their large and bulky design can be a disadvantage in certain applications.Therefore, innovative flat cylindrical GRIN lenses have been developed using various techniques, including transformation optics (TO) [11], [12], ray-tracing optimization [13], and analytical formulas rooted in geometrical optics (GO) [14], [15].All of these methodologies are tailored to specific applications, often constrained by geometric and material limitations, making them less universally applicable.TO requires the existence of both electric and magnetic properties, in conjunction with anisotropy.When TO is rigorously applied, it provides the material parameters and the distribution of the electric and magnetic fields, as well as the Poynting vector within a GRIN lens.Nevertheless, achieving these conditions in practical applications can be a formidable task.In nearly all cases of actual TO device implementation, it becomes necessary to simplify the constituent parameter expressions by approximating them in a feasible manner, often leading to the omission of the magnetic response and, at times, anisotropy.This simplification can introduce notable errors in the original field.Consequently, the availability of efficient tools for rapid analysis of GRIN devices has become increasingly crucial.A full-wave electromagnetic simulator can offer valuable insights into wave propagation within the GRIN lens and the resulting far-field radiation patterns.However, as the electrical size and/or permittivity gradient increase, simulations become computationally demanding in terms of both time and resources.As a consequence, they cannot be effectively employed for the optimization-based design of GRIN lenses, which would necessitate numerous iterations involving varying GRIN lens geometries and refractive indices.To overcome this challenge, one can resort to GO as a very efficient analysis tool.As a special case of the method of characteristics, GO provides a solution of an asymptotic (high frequency) approximation of the wave equation by tracing rays (eikonal equation) and then by calculating the field along the ray trajectories (transport equation).Many previous works deal with the ray-tracing implementation, without addressing the transport equation [16], [17], [18], [19], which ultimately results in an incomplete GO analysis tool.However, recent research (e.g., [20], [21]) has introduced complete GO-based algorithms, including the field calculation.These algorithms not only provide analytical insight but also enable optimization for GRIN lenses.Nonetheless, there remains a shortage of efficient and accurate numerical tools for comprehensively analyzing wave propagation within inhomogeneous materials, particularly when the lens aperture is significant relative to the wavelength, rendering full-wave electromagnetic simulation tools unfeasible.
In this research, we introduce a GO analysis algorithm that combines in a single ordinary differential equation (ODE) the ray tracing and the field calculation with the capacity to analyze a wide range of GRIN lens types while accommodating arbitrary material distributions in space.This algorithm paves the way for novel and inventive designs within the framework of GRIN lenses.Our proposed approach brings several features compared to the methodology outlined in [20] and [22].These innovations encompass the following.
1) More Accurate Field Amplitude Calculation From Wavefront Curvature: In contrast to the method presented in [20], our study employs an additional transport equation for precisely computing the wavefront curvature along the ray.This enables the exact determination of the ray-tube cross-sectional area, which, in turn, permits precise field transport calculations without the need for approximations based on tracking neighboring rays.This improvement enhances the overall accuracy of our analysis.It has to be said that our approach encounters difficulties when dealing with caustics, where the ODE solver of the wavefront transport breaks down.
In contrast, the method in [20], although less accurate, is more robust and allows field calculation through caustics, facilitating inverse design by optimization.Therefore, we are working on a more efficient direct retrieval method for the refractive index in ray design that will eliminate caustic-related problems.This will be the focus of our future work and preliminary results are presented in [23].
2) Increased Ray Processing Independence and Complete GO Parallelization: As a special case of the method of characteristics, GO solves electromagnetic problems along ray trajectories.The calculation of ray trajectories is independent of each other, allowing for easy parallelization by assigning each ray to a separate thread.However, in [20] and [24], the calculation of field transport along the ray requires the estimation of the spreading factor by tracking the distance between a bundle of, at the very least, three rays.This process interconnects different threads, complicating the parallelization.Alternatively, to mitigate this issue, additional dummy rays can be traced alongside each primary ray, as done in [20] and [22], incurring an increased computational burden.The algorithm developed in this research enables the independent processing of individual rays, including field calculations.This feature makes the entire GO analysis fully parallelizable and efficient, not limited to just ray tracing.This enhancement has the potential to significantly boost computational efficiency and expedite the analysis procedure.3) More Efficient Calculation of Wavefront Curvature of the Transmitted Wave: The calculation of the transmitted wavefront curvature, essential for the initial condition in the wavefront-curvature transport equation at the input lens interface, is accomplished through a novel and concise dyadic expression introduced in this work.4) Far-Field Computation Utilizing an Ad Hoc Method: Rather than relying on the fast Fourier transform algorithm derived from the plane wave spectrum, our research adopts an alternative approach to compute the far field.This different technique may offer potential advantages in terms of either accuracy or computational efficiency, especially when the lens output interface is far from being flat.This alternative approach may offer the advantages with respect to an FFT applied directly to the aperture ray-field, especially in terms of estimate of cross-polar components.A careful analysis is provided in [25].These innovative aspects of our work contribute to advancements in the analysis of GRIN lenses, offering insights and potential enhancements in both design and optimization.The remainder of this article is divided into two sections.Section II delves into the development of the differential equations for ray tracing and transport (Section II-A) and their transformation into a unified vector ODE format, well-suited for array programming and efficient parallelization (Section II-B).The utilization of general external astigmatic rays outside the lens region to establish boundary value problems at the input interface is discussed in Section II-C.Section II-D explores the crossing of rays at interfaces between different inhomogeneous media, providing the wave parameters for the transmitted wave, which serve as the initial boundary values for the ODE.Following this, Section II-E presents a comprehensive algorithm for rapidly computing the far field, starting from the field aperture distribution at the output interface.The entire algorithm is summarized with a pseudocode in Section II-F.In Section III, we demonstrate the Fig. 1.GRIN lens antenna geometrical arrangement.A ray is originated at the source point F, hits the lens input interface S in at A, and then travels inside the inhomogeneous volume V on a curved trajectory up to B on the output interface S out .effectiveness and accuracy of the presented algorithm through four numerical examples: the MFE lens (Section III-A), the LL (Section III-B), an unconventional application of a circular scanning GRIN lens antenna presented in [6] (Section III-C), and the zero-focal lens described in [7] (Section III-D), for which a comparison against measurements is available.

II. FORMULATION
The GRIN lens antenna is assumed as a block of isotropic dielectric material with an arbitrary shape and a varying refractive index n(r) [or relative dielectric permittivity ε(r) = n 2 (r)].The refractive index is defined at every point r ∈ V , where V is the lens volume, and it is assumed at least twice differentiable, i.e., n(r) ∈ C 2 (V ).In our assumption, the refractive index is assumed as continuous and differentiable, and no attempt of modeling practical realizations with discretized stepped refractive index profiles is done, which is beyond the purpose of this present article.The boundary of V comprises an input S in and an output S out surface.In the lens intentional operation, the impinging wave enters through S inc , crosses V , and exits through S out (Fig. 1).

A. GO in Inhomogeneous Media
GO [26], [27], [28] is effective and accurate when the size of the lens, the curvature of its surface, and the refractive index changing rate are all large compared to the wavelength λ.According to GO, the time harmonic (e jωt intended and suppressed) electromagnetic field is written as with k = 2π/λ denoting the free-space wavenumber, E 0 (r) the slowly varying field polarization and strength (namely, the leading term of the Luneburg-Kline asymptotic series expansion [29], [30]), and ψ(r) the eikonal function, corresponding to the optical length.The eikonal function satisfies the eikonal equation obtained by introducing (1) in the Helmholtz wave equation for the electric field and collecting the leading asymptotic terms with respect to inverse powers of k.Instead, the first higher order asymptotic terms provide the electric field transport equation where d/ds = ŝ • ∇, with ŝ = ∇ψ/n.The argument dependence (r) is hereinafter omitted for the sake of compactness.
In this article, we consider the case of lossless materials for which n and ε are real.However, the present formulation can be easily extended to the case of small losses, for which ε = ε ′ − jε ′′ becomes complex (ε ′′ ≪ ε ′ ), by assuming n = √ ε ′ and by adding an extra term −(kε ′′ /2n)E 0 in the transport equation ( 3), accounting for weak absorption [28].
Instead of solving the eikonal equation ( 2) by calculating the eikonal in a volume grid of points, an effective modeling scheme is to resort to rays, i.e., a set of curves perpendicular to wavefronts.By taking the gradient of (2), the ray trajectory equation is derived whose solution provides the curved ray trajectory in the inhomogeneuos medium r(s), parametrized by the curvilinear abscissa s.The unit vector tangent to the ray is given by and, according to (2), satisfies ∇ψ = nŝ.The numerical solution of (4) requires the definition of initial conditions r(s = 0) = r A and ŝ(s = 0) = ŝA , i.e., a starting point r A and a launching direction ŝA .
Once the ray is traced, it is claimed in some literature that one can use (3) to calculate E 0 .However, this is not practical because the eikonal Laplacian ∇ 2 ψ = n∇ • ŝ + dn/ds is not available from the single ray trajectory.Such a quantity requires additional information on the wavefront or the rays in the neighborhood of the traced ray r(s).As a matter of fact, the eikonal Laplacian expresses the divergence of the rays which causes the variation of the cross-sectional area d of the ray tube and the spreading of the transported power.The latter is dictated by the power conservation (Poynting theorem) ∇ • (|E 0 | 2 ∇ψ) = 0 applied to a ray tube (Fig. 2).Indeed, according to differential geometry, the divergence of the unit vector normal to a surface (i.e., the wavefront) is related to the surface curvature by ∇ • ŝ = κ 1 + κ 2 , with κ 1,2 = 1/ρ 1,2 denoting the wavefront principal curvatures; i.e., the reciprocal of the principal radii of curvature ρ 1,2 (Fig. 3).
To solve this impairment in [20] and [24], a bundle of two additional paraxial rays are traced around the center ray to estimate the ray-tube cross-sectional area from the distance between the rays.However, a more accurate, efficient, and elegant solution consists in resorting to the wavefrontcurvature approach in [31] (previously developed for seismic waves [32]).The wave-field object is augmented by introducing the local wavefront-curvature dyadic with X1,2 the two orthogonal unit vectors tangent to the wavefront in the direction of the principal planes (Fig. 3).
The curvature dyadic allows a complete second-order local representation of the eikonal at any point along a ray, including a displacement transverse to the ray.Namely, the Hessian of the eikonal is given by in which the curvature dyadic provides the components transverse to the ray, while the longitudinal term is calculated from the refractive index profile; in fact, according to (4) The eikonal Hessian obeys to the transport equation [31] which can be obtained by taking the Hessian of (2).Therefore, by using (9), one can calculate the wavefront curvature along the ray and use it in (3) to calculate the Laplacian of the eikonal as the trace of the Hessian ∇ 2 ψ = tr{∇∇ψ}.

B. Ray Trajectory and
Equations ( 3)-( 5) and ( 9) can be combined in single vectorial, first-order, ODE in which we have chosen the optical length t = ψ A + s 0 nds to parameterize the ray, by introducing the change of variable dt = nds, and is a 15 × 1 column vector collecting all the ray and wavefield information.Namely, in the first three entries, we have arranged the ray point y (1:3)  = r, in the second three entries the eikonal gradient (i.e., the ray direction) y (4:6)  = ∇ψ = nŝ, in the third three entries the electric field y (7:9)  = E 0 , and in the last six entries the Hessian of the eikonal y (10:15)  = vec{∇∇ψ}.It is also worth noting that the definition of dyads according to the dyadic notation Ā = bb is straightforwardly implemented as a matrix product between a column and a row vector b • b T .Furthermore, in (12), F(y) is a 15 × 1 column vector function of the vector argument y, which is defined as follows.The first three entries are given by the second three entries by the third three entries by F (7:9) (y) = − tr vec −1 y (10:15) 2n 2 y (7:9)   − y (7:9)  • ∇n n and the last three entries by Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
Equations ( 14)-( 17) also involve the refractive index n, its gradient ∇n, and its Hessian ∇∇n all evaluated at the ray point r = y (1:3) .Those quantities can be easily rewritten in terms of the relative permittivity ε = n 2 , its gradient ∇ε = 2n∇n, and its Hessian ∇∇ε = 2∇n∇n + 2n∇∇n, if preferred.Finally, in (12), the initial conditions y A are assigned at the starting value of the parameter t = ψ A = ψ(r A ) by specifying the coordinates y (1:3) A = r A of the starting point A; then, from the knowledge of the ray direction, the electric field, and the wavefront curvature at A, the rest of the vector can be completed.Namely, y (4:6) A = n(r A )ŝ(r A ), y (7:9) A = E 0 (r A ), and y (10:15) A = vec{n(r A ) Q(r A ) + D(r A )}, with D defined as in (8).Hence, the solution of ( 12) permits the ray tracing and the wave-field propagation along the ray, by considering one ray at a time, therefore allowing a parallelization of the code.In addition, (12) calculates four quantities at a matrix, which is efficient for those matrix-based calculators such as MAT-LAB.Among the general purpose, ready-to-use, and numerical tools available for the solution of ODEs, we have adopted MATLAB's ode45 without encountering any numerical issues.As a matter of fact, the ODE is not stiff within the GO validity range ∇n/n, ∇|E 0 |/|E 0 |, and ∇ 2 ψ ≪ k, i.e., when the refractive index and the field amplitude exhibit a slow variation compared to the wavelength and far form caustics.That kind of routine accepts as input t A , y A and the pointer to a function calculating F in ( 14)-( 17) and gives as output a matrix containing a sampled version of the solution.In the present case, we obtain a 16-row matrix with the values of t in the first row Each column contains, below the parameter value in the first row, the corresponding wave-field information vector y(t).The number of columns is variable and depends on the ray length and the required accuracy, which is set via an error parameter; the higher the accuracy, the smaller the parameter step, and the larger the number of parameter samples and the computation time.Along with the initial conditions, the ODE solver routine needs a stopping criterion; in addition to an upper limit in the parameter domain t ∈ [ψ A , t max ], one can define an event that is a logical condition, which dictates the solution endpoint parameter value.By choosing as event r(t) = y (1:3) (t) ∈ V , if the starting point A ∈ S in , then the endpoint is found at the point B with coordinates r B ∈ S out , for a parameter value t = ψ B .Hence, all the wave-field information at the exit point B is available in the last column of the output matrix (18).

C. Astigmatic Ray Illumination From the External Region
In some cases, the wave illuminating the lens is originated from a source, which is external to the lens.The region around the lens r / ∈ V is assumed to be air (free space), i.e., a homogeneous medium with n = ε = 1.Here, t = s, ∇n = ∇∇n = 0, and the system of ODEs ( 12) admits the well-known GO analytical solution.A ray starting at the point F, with coordinates r F , in the direction ŝF , continues on a straight path r(s) = r F +sŝ F .The wavefront radii of curvature Fig. 4. Interface between two materials #1 and #2 with respective refractive indices n 1 and n 2 .The unit vector normal to the interface pointing toward #1 is denoted by n.A ray impinges on the interface from medium #1, with impinging direction ŝi , and hits the interface at R. Here, the ray abruptly tilts to the transmitted direction ŝt , because of refraction when n 1 (R) ̸ = n 2 (R).The electric field is also discontinuous across the interface and the transmitted electric field vector can be expressed in terms of the incident field through the dyadic transmission coefficient.The wavefront curvature also changes across the interface; the transmitted wave curvature is affected by the refractive index step and the local interface curvature but also by the local gradient of the refractive indices, for inhomogeneous materials.
increase linearly ρ i = ρ Fi + s along the ray, in the constant directions Xi (i = 1, 2), so that The electric field direction also remains constant, while its amplitude is modulated by a closed-form spreading factor, thus yielding to the well-known closed-form expression for the GO field ( 1) In conclusion, for a known source, by using standard GO expressions (19) and ( 20), the source rays can be analytically traced from the source up to the point A ∈ S in where the wave-field quantities ψ(r A ), E 0 (r A ), and Q(r A ) are calculated in closed form.

D. Discontinuity of the Refractive Index at an Interface
When a ray crosses the input (output) interface S in (S out ), a point A(B) is encountered where the refractive index n (or the relative permittivity ε) might be discontinuous or nondifferentiable, and thus, the ODE ( 12) is not applicable.To calculate the ray direction, the field, and the wave curvature beyond an interface, one has to exploit continuity conditions.
In a general arrangement (Fig. 4), a ray crosses the interface S between two inhomogeneous media at the refraction point R ∈ S, with coordinates r R .
The ray travels from the medium #1 to the medium #2 beyond the interface, characterized by n 1 (r) and n 2 (r), respectively, and crosses the interface at R. The interface surface S is characterized at R by the normal unit vector n pointing outward from #2 to #1 and by the curvature dyadic C = Û1 Û1 /R 1 + Û2 Û2 /R 2 , with R 1,2 the interface surface principal curvature radii in the principal planes along Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
the tangent unit vectors Û1,2 .The incident, reflected, and transmitted fields are written at R in the GO format (1).Then, the eikonal continuity is imposed at R and the neighbor points on the interface, for the various wavefronts.For the transmitted wave, the eikonal continuity reads Next, the eikonal gradient continuity provides the well-known Snell's refraction law Finally, the eikonal Hessian continuity provides the expression of the transmitted wave wavefront-curvature dyad in which, as well as in ( 22), all quantities are evaluated at r A and Di 1 ( Dt 2 ) is defined as in ( 8) with ŝi and n 1 (ŝ t and n 2 ).Equation ( 23) represents the generalization of the equation in [33] for the case of an interface between two homogeneous materials and reduces to it when Di 1 = Dt 2 = 0.In addition, by imposing the electric and magnetic field tangent component continuity at the interface, the well-known dyadic Fresnel transmission coefficient is obtained: Therefore, by using ( 21)-( 24), the ray and wave quantities can be transported through the interface.This is needed to derive the initial condition to initialize the ODE (12) from the impinging field at A and to calculate the field just outside the lens output interface S out at the point B, from the ray and wave-field parameters calculated by the ODE at the curved ray endpoint.It is worth noting that, in certain inhomogeneous lens designs, the refractive index at the interface is unity to avoid reflection; in such cases, ŝt = ŝi and T = 1 (no refraction as well).However, the jump discontinuity in the refractive index gradient still introduces a discontinuity in the wavefrontcurvature dyad.

E. Efficient Calculation of the Far Field
Once rays are traced from the source to the output interface S out and the field is transported across the interface, aperture field samples E t (r m ) = E t 0 (r m )e − jkψ(r m ) are available on the lens aperture at the sample points r m ∈ S out .The far field in the observation direction û is then calculated resorting to the radiation integral where a spherical factor e − jkr /r is intended and suppressed, η is the free-space impedance, and the aperture equivalent magnetic and electric currents are defined as By using ( 26) in ( 25), the radiation integral appears in the form [34] with 0 a slowly varying part of the integrand and ϕ = ψ − û • r the phase function of a rapidly oscillating complex exponential factor.If the aperture samples are sufficiently dense to accurately follow the variation of the slowly varying part (i.e., to follow the variation of the output interface geometry and its normal n, and the transmitted wave eikonal ψ, direction ŝt , and electric field E t 0 ), then one can generate on S out a triangular surface mesh from the sample points r m with a Delaunay triangulation [35].Ready-to-use routines for this operation are available; they take as input the list of N points r m (m = 1, 2, . . ., N ) arranged in an N × 3 matrix and give as output an N t × 3 matrix d, which is a lookup table where any row corresponds to one of the N t triangles and contains in the three columns the indices of the three points that are the three triangle vertices, in a counterclockwise sense.As far as the pth triangle ( p = 1, 2, . . ., N t ), the ith vertex (i = 1, 2, 3) is v p,i = r m , with m = d pi , i.e., the pth row and ith column entry of d.Analogously, F p,i = F(v m ).Triangle sides are readily obtained as l p,i = v p,i+1 − v p,i , with i + 1 intended "mod 3" (i.e., i = 4 → i = 1, see Fig. 5).Hence, by assuming a linear interpolation for G and f on each triangle, the integral ( 27) is evaluated as the sum of the integrals over each triangle, which is calculated in a closed form [25], finally leading to the quadrature rule In ( 28) Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.+ e − jδ p,i δ p,i+1 δ 2 p,i − e jδ p,i−1 δ p,i+1 δ 2 p,i−1 (29) with S p = 1/2|l p,i+1 × l p,i | denoting the pth triangle's area and δ p,i = k[ψ(v p,i+1 ) − ψ(v p,i )− û • l p,i ] the phase gradient component along the ith edge of the pth triangle (note that 3 i=1 δ p,i = 0, and i − 1 is also intended "mod 3"; i.e., i − 1 = 0 → i − 1 = 3).Indeed, unlike in [25], in the present case, under the far-field observation and linear eikonal interpolation assumption, the phase function is exactly linear over each triangle with a constant gradient.Therefore, the required spatial sampling rate is independent of the frequency and can be much larger than the wavelength, as only the slowly varying part of the integrand needs to be reconstructed.Similar to [25] and [36], an iterative scheme might be arranged where the mesh is refined at each step by tracing new rays to add aperture samples, but it is beyond the scope of this present article.Though ( 29) is not defined when observing in the far-field caustic directions on the triangle edge Keller's diffraction cones (δ p,i = 0, i = 1, 2, 3), it exhibits removable singularities there, and it is intended to assume the limit values B p,i δ p,i = 0 = S p j , and . Finally, when observing at the three Keller's cone intersection, i.e., at the forward scattering far-field caustic direction of the triangle, B p,i (δ p,i = δ p,i+1 = δ p,i−1 = 0) = S p /3. Evaluate the source ray-field ψ A , ŝi (r A ), E i 0 (r A ), and Qi (r A ) with ( 19)- (20) 5:

F. GRIN Lens Overall Analysis Algorithm
The formulation presented in Sections II-A -II-E provides a complete tool for the analysis of GRIN lens antennas, which is here schematically summarized in Algorithm 1. First, the lens geometry, the refractive index function n(r), and the source type and location are defined.Then, a loop starts on a given number of rays N r launched from the source to illuminate the input surface S in of the lens.The ray from the source to the surface point A ∈ S in with coordinates r A is traced analytically, and the impinging ray eikonal ψ A , its direction ŝi (r A ), the source illuminating field E i 0 (r A ), and wavefront curvature Qi (r A ) are calculated according to (19)  and (20).The ray-field transmitted beyond the interface is determined by using ( 21)-( 24); namely, ŝt (r A ), E t 0 (r A ), and Qt (r A ) are calculated and then collected in y A , thus allowing the ODE ( 12) initialization.Next, the ODE is numerically solved, thus tracing the curved ray inside the lens up to the endpoint B ∈ S out .The ODE output provides the endpoint location r B , and all the quantities characterize the ray-field impinging on the output interface, i.e., ψ B , ŝi (r B ), E i 0 (r B ), and Qi (r B ). Finally, the quantities are transmitted through the output interface via ( 21)-( 24), thus calculating the aperture field sample E t 0 (r B ) and the transmitted ray direction ŝt (r B ), which are stored.Once the ray-tracing loop is completed, the aperture distribution is available.Hence, the surface is meshed by a Delaunay triangulation and the far-field radiation is calculated resorting to the quadrature rule (28).It is worth mentioning that the ray-tracing loop can be easily and very efficiently parallelized because rays are independent of each other, even assigning any ray to a thread.

III. NUMERICAL EXAMPLES
In this section, we demonstrate the effectiveness of the proposed algorithm through some examples.In the first two examples, we test the numerical solution against well-known analytical benchmarks.Next, two more general and realistic cases are provided.

A. MFE Lens
As a first example, we consider an MFE lens [37], which is a sphere of dielectric with radius a and radial refractive index variation n(r ], where the reference system origin lies at the sphere center.A unit momentum (J 0 = 1 Am) dipole source J 0 = J 0 x is placed at F ≡ A = −a ẑ, i.e., on the sphere "south pole." As it is shown in [37], rays emerging from the point source travel across the sphere and converge at the "north pole" B ≡ a ẑ following circular trajectories.Namely, a ray launched in the direction ŝt A = sin θ 0 (cos φ 0 x + sin φ 0 ŷ) + cos θ 0 ẑ, i.e., with an angle θ 0 from the z-axis, traces an arc of the circumference with center at C = −a cot θ 0 (cos φ 0 x + sin φ 0 ŷ) and radius R = a csc θ 0 .In Fig. 6, a section of the MFE in a vertical plane is shown; the background is colored proportionally to Fig. 9. Ray tracing in an LL.Rays are launched by a dipole source at the "south pole" and shown in a vertical plane.The background color is proportional to the refractive index n with a surrounding medium n = 1 beyond the LL surface (black line).Numerically traced rays (white solid lines) follow the expected elliptic trajectories and focus to infinity.One ellipse corresponding to the theoretical ray trajectory for θ 0 = π/5 = 36 • is traced for comparison (black dashed line).The observation scan where the field is calculated for comparison (see Fig. 10) is also shown (red dashed line).Fig. 10.Amplitude of the electric field on the output interface (red dashed line in Fig. 9) in the LL illuminated by the unit momentum electric dipole at the "south pole."The scan is parametrized by ρ ∈ [−a, a].The proposed algorithm calculation (green solid line) and analytical formula (red dashed line), when observing in the dipole plane φ 0 = 0 or in the orthogonal plane φ 0 = π/2, are almost identical.The ray trajectories calculated with the proposed algorithm are shown (white solid lines) and precisely follow the analytical circular trajectories, one of which is reported for comparison (dashed black line) corresponding to the launching direction θ 0 = 2π /5 = 72 • .It is worth mentioning that at the source point F, the source GO field is singular, as F is a caustic where ∇ 2 ψ = ∞.Also, since on the MFE surface n = 1, then ŝt = ŝi (no refraction) and, with the source at the surface, Qt = Qi = ∞ because the discontinuity introduced by (23) becomes negligible; therefore, the fifth step of the algorithm is trivial, being the transmitted ray field equal to the incident one.
However, to numerically manage this condition, a tiny sphere with radius r 0 = 10 −6 a is considered in the MFE around the source, and the ODE starting points A are taken on its surface, slightly displaced from F, thus avoiding the singularity.The source ray-field quantities at r A = −a ẑ + r 0 (r 2 +a 2 ) 2 −(2za) 2 of the circular wavefronts in Fig. 7, along three ray trajectories: the central straight ray along the z-axis (θ 0 = 0), the abovementioned ray (θ 0 = 2π /5), and the extremal ray traveling on the sphere surface (θ 0 = π/2).It is observed a perfect correspondence between the curvature provided by the proposed algorithm (green solid line) and the analytical benchmark (red dashed line).Finally, the electric field is calculated along the x-axis (φ 0 = 0) and y-axis φ 0 = π/2 in the range ρ ∈ [−a, a] (red dashed line), and its amplitude is compared against the analytical closed-form formula in Fig. 8.The ray-tracing numerical prediction (green solid) and the analytical formula (red dashed) are superposed.In the two scan planes, the field profile is different because of the dipole radiation pattern, which is omnidirectional in the φ 0 = π/2 plane (orthogonal to the dipole), whereas it exhibits nulls in the φ 0 = 0 plane (along the dipole direction).

B. Luneburg Lens
As a second example, we consider an LL [29], which is also a radially symmetric spherical GRIN lens, with refractive index radial variation n(r ) = (2 − (r/a) 2 ) 1/2 .In this configuration, which is the most famous but only one special case of the more general kind, rays emanating by a source on the sphere surface travel across the LL along elliptic trajectories and exit from the other side all parallel, i.e., focusing to infinity.We adopt the same reference system, notation, and illumination of the previous MFE example.Here, a ray launched from the unit dipole at the "south pole" in the direction ŝt A , i.e., with an angle θ 0 from the z-axis, traces an arc of the ellipse with center at the origin, major axis along sin(θ 0 /2)(cos φ 0 x + sin φ 0 ŷ) + cos(θ 0 /2)ẑ, i.e., with an angle θ 0 /2 from the z-axis, and semiaxes a(1 ± cos θ 0 ) 1/2 .In Fig. 9, a section of the LL in a vertical plane is shown; the background is colored proportionally to the refractive index (from blue n = 1 to yellow n = √ 2), with a surrounding medium n = 1 beyond the LL surface (black line).The rays calculated with the proposed algorithm are shown (white solid lines) and precisely follow the analytical elliptic trajectories, one of which is reported for comparison (dashed black line) corresponding to the launching direction θ 0 = π/5 = 36 • .The electric field is calculated on the LL output surface (red dashed line) both in the xz (φ 0 = 0) and yz (φ 0 = π/2) planes, and  √ a 2 −ρ 2 in Fig. 10.Again, the raytracing numerical prediction (green solid) and the analytical formula (red dashed) are superposed and the field profile in the two scan planes varies because of the dipole radiation pattern.It is worth noting that along these observation scans, the wavefront is cylindrical Qi (r B ) = ρ 2 a 2 √ a 2 −ρ 2 ρ ρ in the internal side of the output surface, whereas it becomes flat Qt (r B ) = 0 (plane wave) on the external side, according to the jump prescribed by (23).In Fig. 11, the curvatures κ 1,2 , i.e., the eigenvalues of the curvature dyad Qi (r B ), along the output interface (red dashed line in Fig. 9) are plotted by comparing the numerical solution of the ODE and the analytical benchmark of this canonical shape, revealing a perfect matching.

C. Circular Beam-Scanning Lens
As a third example, we consider the offset inhomogeneous lens designed in [6] for spaceborne weather radar.The lens optimization exploits both the interface shapes and the refractive index profile to focus on the infinity of the rays arising from an offset source.To permit a conical scanning, the lens geometry and refractive index are circularly symmetric, i.e., interface surfaces are described by z = f 1,2 (ρ) and n(ρ, z).The feed, located at r F ≡ 5/6a ŷ, radiates a x linearly polarized spherical wave with an amplitude pattern cos 4.45 (θ ), pointed toward the lens center, with a 29 • tilt.The 3-D view of the lens geometry is shown in Fig. 12(a), while in Fig. 12(b), a 2-D view in the focus plane is reported where the background is proportional to the refractive index.The reader is addressed to [6] for the details of the exact geometry.The proposed algorithm is used to calculate the far field radiated by the antenna at various frequencies ( f = 20, 40, and 60 GHz) when the radius of the lens is a = 9 cm, and the results are compared against those provided by the full-wave simulation software CST Microwave Studio.The analysis was performed with 186 rays starting from the focal feed point, hitting the input surface and propagating through the inhomogeneous lens up to a point on the output surface.The ray trajectories and the resulting aperture triangular mesh are shown in Fig. 12.It is noteworthy that the aperture mesh consists of 338 triangles, with an average area of 0.28λ 2 , 1.12λ 2 , and 2.53λ 2 and an average side length of 0.72λ, 1.45λ, and 2.18λ, at 20, 40, and 60 GHz, respectively.Despite the coarse sampling, especially at the highest frequency f = 60 GHz, the quadrature formula (28), unlike FFT-based formulas, allows accurate evaluation of the radiation integral on a frequency-independent mesh over a nonflat aperture.Indeed, the normalized farfield pattern of the weather radar lens antenna is shown in Fig. 13 as predicted by the proposed GO algorithm (upper row) and compared to the corresponding full-wave CST simulation (lower row) for the different frequencies.It is apparent that GO correctly reproduces the main beam and the first sidelobes, whereas diffraction effects occurring at the lens rim, which are responsible for conical far sidelobes, are not present, while they appear in the full-wave simulation.It is rather apparent in Fig. 12(a) that the rays emitted from the feed do not reach the far end of the aperture, which is in a GO shadow region where no GO ray arrives.Therefore, the GO prediction of the aperture distribution provides no field in this area.Augmenting GO through diffracted rays, such as by using the uniform theory of diffraction, may enhance precision and describe the field in such a shaded part of the radiating aperture.Nevertheless, its contribution in the far field does not significantly impact radiation in the main lobe.In fact, such diffraction lobes are much lower than the main beam and become negligible when the frequency increases, as expected.The lens radiates a pencil beam in the θ = 33 • direction, which becomes more directive for higher frequencies, as the lens' electric size increases.It is worth noting that the ray tracing is performed once (in less than 1 s, with a nonparallelized, interpreted MATLAB code) and the various radiation patterns at the various frequencies are calculated later on (in approximately 41 s for each frequency point).On the other hand, the CST analysis requires 12M/62M/180M cells and 50/138/253 min at 20/40/60 GHz on a 64-bit operating system and 512-GB RAM.
Fig. 14 shows the far-field gain in the horizontal cut on the focal plane φ = π/2.In this picture, the accuracy of the main beam prediction of the GO (black solid line) can be better and quantitatively appreciated in terms of realized gain for the different frequencies, against the CST results (red dasheddotted line).It is evident that the gain peak value is well estimated and the pattern shape and sidelobes are accurately predicted except for the abovementioned rim diffraction effect.
The GO algorithm also permits the estimation of the power efficiency of the lens antenna system, which is useful for antenna design and assessment.Indeed, by integrating the Poynting vector flux on various surfaces, one can track the power spreading in the system.The flux integral can be easily calculated as W p,i S p 3 (30) in which W p,i = [(n|E 0 | 2 /2η)ŝ • n] r=v p,i is the normal component of the Poynting vector calculated at the sample point r m = v p,i where the mth ray intersects the integrating surface S, which corresponds to the ith vertex of the pth triangle in the Delaunay mesh of S, and S p the pth triangle area, analogously to the radiation integral calculation.Applying (30) to the input S in and output S out surfaces by considering both the impinging E i 0 and the transmitted E t 0 field, and to the farfield sphere, it is calculated that the spillover efficiency is e so ≃ 79%, reflection efficiency e r ≃ 97% both at S in and S out so that the total loss efficiency is e loss ≃ 74%.These efficiencies are frequency-independent, provided that the feed pattern is assumed independent of the frequency.Conversely, the aperture efficiency depends on the frequency because the aperture phase error in the main beam direction increases for increasing frequencies, indeed e ap = 52%, 37%, and 21% at f = 20, 40, and 60 GHz, respectively.

D. Zero-Focal Lens
To illustrate a comparison between the proposed GO algorithm prediction and measurement results, we consider a final example involving the zero-focal GRIN lens presented in [15].This lens, which was fabricated as described in [7], has a focus on the input surface, hence its name.The analytical expression for the radially varying refractive index profile has been obtained through the inverse truncated Abel transform and it is given by n(ρ) = n 0 /cosh(πρ/2d), in which n 0 = cosh(πa/2d) ≃ 2 is the maximum refractive index value at the center (ρ = 0), while d = 9 cm and a = 7.5 cm are the lens thickness and radius, respectively.A 3-D perspective of the lens geometry together with the ray trajectories (white lines) is shown in Fig. 15(a), while a section in the vertical plane is reported in Fig. 15(b) where the background color corresponds to the refractive index.
The feed is in direct contact with the input surface at the focus, thus maximizing reflection and spillover efficiency, and radiates a spherical wave with an amplitude pattern U (θ ) = cos 3 (θ ).The GRIN lens has been fabricated through additive manufacturing, specifically utilizing the fused filament fabrication (FFF) technique, commonly known as 3-D printing.The refractive index profile was divided into cylindrical layers, which were synthesized with hollow structures with different filling fractions.These correspond to intermediate refractive index values between air (n = 1) and dielectric material PREPERM (n = 2.24).For a comprehensive understanding of the fabrication process, please refer to [7] for detailed information.Pictures of the lens prototype display both the rear view (Fig. 16     Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

IV. CONCLUSION
A novel and efficient GO algorithm has been introduced for a GRIN lens illuminated by a general astigmatic source.The algorithm involves tracing curved ray paths from points on the lens input interface and combining ray tracing with field and wavefront-curvature transport equations.This results in a single independent vector ODE for each ray, enabling parallelization of the algorithm by assigning different rays to separate threads.The astigmatic ray behavior at the lens interfaces is achieved by augmenting Snell's law transmission coefficients with a compact closed-form dyadic formula for the transmitted wavefront curvature.The ODE solution provides samples of the field aperture distribution on the output lens interface, enabling far-field calculations using a novel ad hoc fast and efficient quadrature rule.Several examples have demonstrated the accuracy and efficiency of the algorithm in general cases.

Fig. 2 .Fig. 3 .
Fig. 2. Ray-tube geometry for a curved ray in an inhomogeneous medium.The ray-tube cross-sectional area d varies along the ray according to ray divergence.
Wavefront-Curvature Transport in a Single ODE Before proceeding further, let us introduce the operator vec{•}, which reshapes the triangular part of a 3×3 symmetric matrix into a 6 × 1 column vector vec vec −1 {•}, which rebuilds the 3 × 3 symmetric matrix from the 6 × 1 column vector vec −1

Fig. 5 .
Fig.5.Arrangement of the pth triangle in the Delaunay mesh.The vertexes are pointed by the vectors v p,i , i = 1, 2, 3, in the general reference system, and they are counted in a positive counterclockwise sense with respect to the outgoing normal n.Each side corresponds to the vector l p,i = v p,i+1 − v p,i .

Fig. 6 .
Fig.6.Ray tracing in an MFE.Rays are launched by a dipole source at the "south pole" and shown in a vertical plane.The background color is proportional to the refractive index n, with a surrounding medium n = 1 beyond the MFE surface (black line).Numerically traced rays (white solid lines) follow the expected circular trajectory and converge at the "north pole."One circle corresponding to the theoretical ray trajectory for θ 0 = 2π/5 = 72 • is traced for comparison (black dashed line).The observation scan where the field is calculated for comparison (see Fig.8) is also shown (red dashed line).

Fig. 7 .
Fig.7.Wavefront curvature κ along the rays in the MFE as calculated by the proposed algorithm (green solid line) and analytical value (red dashed line).Three different rays are considered with initial launching directions θ 0 = 0 (straight central ray along the z-axis), θ 0 = 2π/5 (ray compared to the analytical circular trajectory in Fig.6), and θ 0 = π/2 (extremal ray on the MFE surface).

1 :
Define Lens Geometry, n(r), source type and location //ray-tracing loop 2: for m = 1 to N do 3:Launch the mth ray from the source and find the respective hit point A ∈ S in location r A 4:

Fig. 8 .
Fig. 8. Amplitude of the electric field on the linear scan ρ ∈ [−a, a] (red dashed line in Fig. 6) in the MFE illuminated by the unit momentum electric dipole at the "south pole."The proposed algorithm calculation (green solid line) and analytical formula (red dashed line), when observing in the dipole plane φ 0 = 0 or in the orthogonal plane φ 0 = π/2, are almost identical.

Fig. 12 .
Fig.12.Geometry of the circular beam-scanning GRIN lens antenna designed in[6] for a spaceborne weather radar.(a) Three-dimensional view of the antenna system.Ray trajectories start at the feed point and hit the input surface (orange straight lines), then travel through the lens (white curved lines), and finally reach the output surface.At any endpoint B ∈ S out , the transmitted ray direction ŝt is drawn (orange arrow).The triangular Delaunay mesh for the radiation integral is also depicted (black lines).(b) 2-D section of the same antenna in the focal plane is also shown with the background color corresponding to the refractive index n(ρ, z).
and Qt (r A ) = ( 1 − ŝt A ŝt A )/r 0 .Such values are used to initialize the ODE.Along ray trajectories, the numerical solution of the ODE provides the eikonal Hessian as ∇∇ψ = vec −1 {y (10:15) }, which allows the calculation of the wavefront-curvature dyad Q = ( 1 − ŝŝ) • ∇∇ψ • ( 1 − ŝŝ)/n.The two nonzero eigenvalues of Q are the numerically calculated local wavefront curvatures κ 1,2 ; they are equal to each other and they are compared with the expected analytical value κ = − 2z √

Fig. 14 .
Fig. 14.Far-field normalized pattern of the weather radar inhomogeneous lens antenna in Fig. 12 for f =20 (right column), 40 (central column), and 60 GHz (left column), as predicted by GO (top row) and CST full-wave simulation (bottom row), in the uv plane.

Fig. 15 .
Fig. 15.Geometry of the zero-focal GRIN lens [7].Three-dimensional view of the antenna system (left).Ray trajectories start at the feed point on the input surface and travel through the lens (white curved lines) up to the output surface.At any endpoint B ∈ S out , the transmitted ray direction ŝt is drawn (orange arrow).A 2-D section of the same antenna in the focal plane is also shown (right) with the background color corresponding to the refractive index n(ρ).Prototype of the zero-focal GRIN lens made with 3-D printing (bottom).(a) Rear view.(b) Front view.

Fig. 17 .
Fig. 17.Far-field pattern of the zero-focal lens antenna at f = 18 GHz.Normalized far field as predicted by (a) presented GO algorithm and (b) experimental results in the uv plane.(c) Directivity pattern in the E-plane (φ = π/2 cut); GO prediction (black solid line) and measurements (red dashed line).
(a)), where the open waveguide feed is visible, and the frontal view (Fig. 16(b)), which reveals the layered structure and the distinct geometries of each layer.The normalized far-field pattern at f = 18 GHz in the uv plane, as predicted by the proposed GO algorithm, is shown in Fig. 17 (a) along with the corresponding experimental results in Fig.17 (b).The same comparison in terms of antenna directivity is reported in Fig.

17
(c) in the φ = π/2 (E-plane) cut.It is evident that, as well as in the previous example, the GO algorithm accurately reproduces the main beam.