Open Access
Add to BibSonomy
Abstract
Ray tracing in gradient-index (GRIN) media has been thoroughly studied and several ray tracing methods have been proposed. Methods are based on finding the ray path given a known GRIN. In recent decades, the inverse problem, which consists of finding the GRIN distribution for a given light ray path, has been gaining attention. Given that it is not an easy task, the methods proposed in the literature vary in degrees of difficulty. In this work, an alternative method is presented to derive symmetric GRIN distributions whose implementation can be considered the simplest to date. Since it is based on invariants, which result from the symmetries of the system as stated by Fermat’s principle, it is an exact numerical method, i.e., the physical system is not approximated. The robustness of the method permits the reconstruction of the GRIN distribution from a ray propagating in three-dimensions. In order to demonstrate its operation, different known symmetric GRIN media are reconstructed using rays that propagate in two and three dimensions.
© 2023 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
Gradient index (GRIN) media are characterized by a inhomogeneous index of refraction distribution that changes spatially [1,2]. These GRIN media can be classified depending on the geometry of their isoindicial surfaces, i.e., it depends on the geometry of their surfaces in which the refractive index is constant. Taking the isoindicial surfaces into consideration, and taking advantage of the fact that a new frontier in freeform optics in the form of freeform GRIN media has just been defined in 2021 [3], we propose that GRIN media can be classified into two groups:
- 2 Freeform GRIN media having refractive index distributions that may vary arbitrarily in three spatial dimensions, which means that there are not symmetries associated to a curvilinear coordinate system [3].
The first group contains the most common three GRIN media types studied in the literature: the axial, cylindrical, and spherical GRIN medium. In the first and second types their isoindicial surfaces are perpendicular or parallel planes and concentric cylinders to the optical axis, respectively. In the latter case, the isoindicial surfaces are concentric spheres around the center point [1]. The second group contains the more general kind of GRIN media, where the corresponding isoindicial surfaces are freeform. Within the first group we can also consider a modified elliptic GRIN media that has been used to create a very approximate model of the human crystalline lens [4]. An example of the freeform group has been used very successfully to model the physiological characteristics of the human lens [5,6].
In a GRIN medium, geometrical-light-propagation follows curved trajectories [7,8]. These curves are described by the solution to the differential equation of the light rays, which is given by
where $\textbf {r}$ is a position vector of a typical point on a ray, $s$ is the length of the ray measured from a fixed point on it, and $n$ is the GRIN distribution of the medium [9,10]. There are limited number of examples of GRIN distributions for which closed analytical solutions to Eq. (1) have been found [11–18], while for the most cases, numerical methods have to be employed [19–21]. Recently, it has been shown that the precise trajectory of rays in a symmetric GRIN medium can be determined without the necessity of obtaining an analytical or numerical solution for the ray equation [22,23], and it is called Physical Ray Tracing (PRT). This method is based on the conservation of quantities inherited from the symmetries of the GRIN distribution as stated by Fermat’s principle. We refer to these quantities as Fermat’s ray invariants. This leads to an exact numerical method to be implemented in a straightforward manner. The term ’exact’ refers to the fact that no approximation is applied to the physical system apart from that of numerical origin. Several ideas for such methods arise from an analogy between classical/quantum mechanics and optics which has regained significant interest and has been actively developed for a few decades [24–28]. This growing fascination stems from the realization that there are striking parallels between certain fundamental concepts and mathematical formalisms of classical and quantum mechanics and those found in the study of light and its behavior.All the methods mentioned above address the direct problem: to find the ray propagation in a given GRIN distribution, and it has been thoroughly studied. However, the inverse problem, which consists in finding the GRIN distribution for a given light ray propagation, has barely been studied, as it is not an easy task [29–37]. The inverse problem has even been dealt with using neural networks [38,39]. Recently, GRIN distributions with spherical symmetry have been gaining interest because of their stigmatic properties [38–42]. For instance, a method based on the optical Binet equation has been developed to solve the inverse problem [43].
In this paper, a method to find the general solution to the inverse problem for symmetric GRIN media is presented. This method is based on Fermat’s invariants, just as in the PRT method [22,23]. A remarkable advantage of this method is that it does not require any information about solutions to any differential equation, but the the ray trajectory and the corresponding invariants. Also, due to the physical nature of the principle on which the method is based, it is straighforward to implement in two and three dimensions. We call this method Physical GRIN Reconstruction (PhysGRIN). To demonstrate its operation, the PhysGRIN method is implemented in different relevant symmetric GRIN distributions.
2. Fermat’s ray invariants
Fermat’s principle states that “the light ray between two points $A$ and $B$ is the curve for which the optical path length (OPL) attains an extreme value” [1,7]. Mathematically, Fermat’s principle can be succinctly stated as
where $\mathcal {L}(q_{j},\dot {q_{,}},s)=n({q})$ is the optical Lagrangian, $q_{j}$ $(j=1,2,\ldots )$ are generalized coordinates associated to an arbitrary curvilinear orthogonal coordinate system and are parameterized by $s$, i.e., $q_{j}(s)$, and $\dot {q_{j}}=\text {d}q_{j}/\text {d}s$.In this variational problem, the extremum condition is expressed as
and the solution must satisfy the Euler-Lagrange equations, which are given byGiven a symmetry of the GRIN that manifests itself in the invariance of the Lagrangian with respect to a particular generalized coordinate, say $q_j$, Eq. (4) becomes
In the case of symmetric GRIN media, specifically with axial, spherical and cylindrical symmetry, it has been shown that Fermat’s ray invariant is proportional to the following product
where $n(\textbf {r})$ is the GRIN distribution, $f$ is a function of $\varphi$, which is the angle between the vector $\nabla n(\textbf {r})$ and the vector $\textbf {s}$ tangent to the ray, as shown in Fig. 1 [7,22,23,46,47]. Notice that $\nabla n(\textbf {r})$ is orthogonal to the isoindicial surface at point $\textbf {r}$.If the geometry of the isoindicial surfaces and the ray path are known, then $\varphi$ is given by
where $\hat {\nabla } n(\textbf {r}):={\nabla } n(\textbf {r})/|{\nabla } n(\textbf {r})|$, and ${}\hat {(.)}{}$ denotes a unit vector.Equations (7) and (8) are the basis of our PhysGRIN method, as will be seen in the following sections.
3. Spherical gradient-index medium reconstruction
It is well-known that a spherical symmetric GRIN medium distribution is described by a function $n=n(r)$, when written in spherical coordinates $\left (r,\theta,\phi \right )$, and that the associated Optical Lagrangian, at the $x$-$y$ plane ($\theta =\pi /2$), is given by
where $\phi _{r}=\textrm {d}\phi /\textrm {d}r$. In this case, the Lagrangian does not depend on $\phi$, which leads to the following Euler-Lagrange equationSubstituting Eq. (9) into Eq. (11), and using the geometry shown in Fig. 2, it is possible to find, after some algebra, that
where $K$ is the Fermat’s ray invariant associated to the spherical symmetry of the GRIN (Eq. (12) is also known as the generalized Snell’s law for inhomogeneous media with spherical symmetry [41]), and the angle $\varphi$ can be written, according to Eq. (8), as where $\hat {\textbf {r}}$ is the unit position vector. In this case, Eq. (8) becomes Eq. (13) when spherical symmetry is considered, since ${\nabla }n(\textbf {r})$ and $\textbf {r}$ are always parallel [7].If the direct problem is addressed, that is, if the GRIN distribution is known, $K$ can be calculated on the surface of the lens using Eq. (12), which is given by
where $\varphi _{1}$ is the value of $\varphi$ at the surface of the lens (at $P_{1}\left (x_{1},y_{1}\right )$, see Fig. 2), $R$ is the radius of the lens, and $n\left (R\right )$ is the refractive index of the lens surface.However, to address the inverse problem, $K$ must be found independently from the GRIN distribution. This can be done, as it has recently been shown, as follows [23]:
where $n_{e}$ is the external refractive index in which the GRIN distribution is immersed, and $\theta _{Inc}=\theta _{0}+\theta _{1}$ (see Fig. 2).Now, from Eqs. (12) and (15), it is possible to obtain the GRIN distribution as a function of the given path of the ray, that is,
This is one of the main results of this work. It establishes the solution to the inverse problem; once a ray path is given, Eq. (16) assigns a value to the GRIN for each point on the path through the values of $r$ and $\varphi$ by means of Eq. (13). Hence, this process leads to the reconstruction of the GRIN distribution $n(r)$ throughout the ray propagation path.
Let us first analyze the case in which the reconstruction is performed using a ray path that is known to be the result of propagation in a spherical symmetric medium that enters and leaves the GRIN. This situation can occur, for instance, when experimental data are obtained from propagation in a GRIN medium of which only its spherical symmetry is known. Such a path has a point of minimal distance to the center of the GRIN. This point will be called $P^{*}$, whose distance to the origin is denoted by $r^{*}$ (see Fig. 2), which satisfies in general $r^{*}\geq 0$. Then, the reconstruction will not necessarily span the entire interval $0<r<R$, but the subinterval $r^{*}<r<R$. It can be shown that $P^{*}$ is reached when $\varphi =\pi /2$, yielding $K=r^{*}n\left (r^{*}\right )$. In summary, the GRIN distribution is only reconstructed on an interval given by
In order to establish a systematic procedure of GRIN reconstruction, let us assume that such a ray is discretized. This can be the case, as mentioned above, when the ray path is known from experimental data. Of course, a continuous parametrization of a ray path can be discretized on demand, so this procedure can be applied as well in this case. Let us take a discretized ray inside a spherical GRIN lens as the one shown in Fig. 3, and let $P_{i}=\left (x_{i},y_{i}\right )$ be the points of the ray path, where $i$ is the index of the discretization. Then, the values of $\varphi$ and $r$ at the point $P_{i}$, which are necessary to reconstruct the GRIN (see Eq. (16)), are given by
and where $\hat {\textbf {r}}_{i}=\frac {\left (x_{i},y_{i}\right )}{\sqrt {x^2_{i}+y^2_{i}}}$ and $\hat {\textbf {s}}_{i}=\frac {\left (x_{i+1}-x_{i},y_{i+1}-y_{i}\right )}{\sqrt {\left (x_{i+1}-x_{i}\right )^{2}+\left (y_{i+1}-y_{i}\right )^{2}}}$. Notice that the value of $\hat {\textbf {s}}_{i}$ is calculated using two neighboring points $P_{i}$ and $P_{i+1}$ (see Fig. 3). Notice that there is also an alternative way of obtaining the angle $\varphi$ through geometric analysis, as it is shown in Supplement 1.By substituting Eqs. (18) and (19) into Eq. (16), the value of the GRIN distribution at the point $P_{i}$, is given by
Eq. (20) represents the discretized version of Eq. (16), where the invariant $K$ is calculated using Eq. (15).
3.1 GRIN reconstruction of the Luneburg lens
We proceed to show three instances of the reconstruction of the GRIN distribution for the spherical Luneburg lens. Each reconstruction requires one single ray path which is obtained by solving the Eikonal equation for the Luneburg GRIN distribution $n\left (r\right )=\sqrt {2-r^{2}}$, as given by [48], and each ray possesses a different set of initial conditions, they are presented in Table 1. The resulting ray is given by
Table 1. Parameters of the GRIN distribution
Through Eqs. (20)–(22), the GRIN distribution of the Luneburg lens can be reconstructed. The parameters required to trigger such a process, and the reconstruction interval for each of the test rays (see Fig. (4)), are shown in Table 1. The reconstruction of the GRIN distribution $n(r)$ is shown in Fig. 5.
Fig. 4. Three different examples of rays supported by the Luneburg lens. Despite its appearance, Ray 3 does not coincide with the propagation axis (see Table 1).
Fig. 5. Three different reconstructions corresponding to three distinct rays. The interval of reconstruction depends on $r^{*}$; the minimum radius of the points on the ray path. Ray 1 has the largest $r^{*}$, leading to a shorter interval of reconstruction, whereas Ray 3 has the smallest $r^{*}$, leading to the largest reconstruction interval.
As can be noticed from Table 1 and Fig. 5, the reconstruction interval depends on the path of the ray; it will be longer for rays that traverse through more isoindicial surfaces. For instance, Ray 3 has the longest reconstruction interval while Ray 1 has the shortest interval. Indeed, since the isoindicial surfaces are concentric spheres, the ray that passes closer to the optical axis, will traverse through more surfaces.
3.2 GRIN reconstruction error
There are several metrics to evaluate the overall accuracy of a reconstruction. In this work we will use a local measure as well as a global one. The absolute error (AE) is a pointwise metric calculated as the absolute value of the difference between target $n(r_{i})$ and reconstruction $\hat {n}(r_{i})$:
The global measure is known as the normalized root mean square error (NRMSE), and is defined as follows:
The AE for all three rays is shown in Fig. 6(a). However, let’s focus our attention on Ray 3, which generates the longest reconstruction interval of the spherical GRIN. This ray generates a maximum AE ($\text {AE}_{Max}$) of $7.07047x10^{-5}$, for $m=10,000$, which means that the error is generated in the fifth decimal position.
Fig. 6. (a) The absolute error (AE) for the reconstruction of the three rays of Fig. 5. (b) AE and (c) the normalized root mean square error (NRMSE) in decibels (dB) for the Ray 3 with different values of $m$.
Of course, the reconstruction process depends on the number of points that are generated between the points $P_{1}$ and $P^{*}$, that is, it depends on the value of $m$. For this reason, the AE for Ray 3 for different values of $m$ is presented in Fig. 6(b). Notice that for $m=70$ we get $AE_{Max}=9.99442\text {x}10^{-3}$. This means that even for small number of points in the discretization, the approximation remains within $10^{-3}$ error (precision). However, the AE presents a local measure of the error. To get a more general idea we have measured the value of the NRMSE for different values of $m$. Figure 6(c) shows that for small $m$, there are small values of NRMSE with an order of magnitude of $10^{-3}$. For larger $m$, NRMSE with order of magnitude of $10^{-5}$ can be achieved.
The analysis above shows that the GRIN reconstruction method we developed is effective and accurate. Furthermore, in this section we show that the method is easy to implement and easy to translate to any programming language.
The method works with a similar precision for any symmetric GRIN media as long as the Fermat’s ray invariant can be calculated. For this reason, in the next section 4, let us consider a medium described by a cylindrical GRIN distribution.
4. Cylindrical gradient-index media reconstruction
The cylindrical GRIN distribution is represented by a function $n=n(\rho )$ where $\rho =\sqrt {x^{2}+y^{2}}$ is the radial distance in cylindrical coordinates $(\rho,\phi,z)$. Notice that $\rho$ is contained in the $x$-$y$ plane for each $z$, i.e., the vector $\boldsymbol {\rho }=(x,y,0)$, in this distribution, is perpendicular to the $z$ axis. The optical Lagrangian of this system is given by
where $\phi _{\rho }=\textrm {d}\phi /\textrm {d}\rho$, and $z_{\rho }=\textrm {d}z/\textrm {d}\rho$. Its corresponding Euler-Lagrange equations in the variable $z$ and $\phi$ areFig. 7. The green plane shows the geometrical parameters of the incident ray, whereas the adjacent plane shows the parameters of the transmitted ray at the point $P_{1}$. The inset shows the projection on the plane $x$-$y$ of the latter planes.
These invariants are the manifestations of the translational and rotational invariance of the cylindrical GRIN medium. The parameter $\tilde {l}$ is usually referred to as the skewness parameter [7]. In the special case when $\tilde {l}=0$, the ray paths are confined into a plane. These confined rays are called meridional rays.
4.1 Meridional rays
Now, employing the PhysGRIN method, we consider the reconstruction of a cylindrical symmetric GRIN using a meridional ray. In this case, only the Fermat’s ray invariant $\tilde {\beta }$ has to be considered, and it can be found in a way that does not depend on the value of the GRIN distribution [23], i.e.,
where $\theta _{cInc}=\frac {\pi }{2}-\theta _{r}$ (see Fig. 7).From Eqs. (28) and (31), we obtain
The latter establishes the reconstruction of the cylindrical GRIN distribution, once a ray path is given.
To implement the method, let us consider the following cylindrical GRIN distribution
Fig. 8. Propagation of a ray path in a meridional plane with the following parameters: $a=5$, $\Delta =0.2$, $\theta _{r_{1}}=33.6901^{o}$, $n_{e}=n_{1}$, and $\theta _{cInc}=56.3099^{o}$.
Using Eqs. (30)-(32) together with the ray path given in Eq. (34) and its parameters shown in Fig. 8, it is possible to reconstruct the cylindrical GRIN distribution. Notice that in this case, Eq. (30) becomes $\varphi _{ci}=\cos ^{-1}\left (\hat {{\textbf {x}}}_{i}\cdot \hat {\textbf {s}}_{i}\right )$ where $\hat {{\textbf {x}}}_{i}=(x_i/|x_{i}|, 0, 0)$. This process is shown in Fig. 9(a). The reconstruction interval is given by the maximum and minimum heights reached by the meridional ray. In this case, both heights are equal. Following the procedure shown in Section 3.2, the reconstruction was performed with $m=912$ points, which generates a value of $\text {AE}_{Max}=3.9916\text {x}10^{-4}$, as shown in Fig. 9(b), and a value of the $\text {NRMSE}=2.2242\text {x}10^{-4}$. To optimize the reconstruction means that the unreconstructed area, in yellow Fig. 9(a), should be null, which is achieved for a maximal incidence angle that occurs when $\tilde \beta =n_1 \sqrt {1-2\Delta }$ and coincides with the core-cladding refractive index [7].
Fig. 9. (a)-(b) Reconstruction of the GRIN distribution at the meridional plane for the ray shown in Fig. 8. The red dotted lines mark the boundary of the region where $n(x)$ is reconstructed. (c) The AE of the reconstruction as a function of $x$. The orders of magnitude of the AE and NMRSE are $\sim 10^{-4}$, which are not noticeable at the scale of (a) and (b).
The PhysGRIN method for meridional rays may be applied to an axial GRIN distribution, since this kind of media are described by
Eqs. (36) and (37) correspond to Eqs. (31) and (30) of the meridional case, respectively. For this reason, an axial GRIN distribution is not reconstructed in this work.
The methods described in this section possess broad and versatile applicability, as they can be employed effectively in both discrete and layered distributions. This approach is particularly relevant as it can be applied to cases where distributions exhibit these specific characteristics [50–53]. As emphasized above, the key factor lies in the invariant $K_L$ which is independent of the specific index distribution, ensuring the methods’ applicability and efficiency in diverse scenarios.
4.2 Helical rays
The second case of a ray propagating in a cylindrical GRIN is that of a ray that is not confined to a meridional plane, i.e., the ray enters at the point given by $x=a'$, $y=0$, $z=0$, and making an angle $\theta _{r}$ with the $z$-axis. This condition is a special case of a skew ray and is called helical ray [7]. The ray path is given by
In previous works [22,23] it has been pointed out that in order to propagate non-meridional rays by the PRT method, two invariants have to be known. However, it is not the case with the PhysGRIN method, as both invariants are set by the refractive index according to Eqs. (28) and (29), that is, one invariant is enough to reconstruct the GRIN distribution. For this reason, we will focus our attention only on $\tilde {\beta }$, just as in the meridional case.
Taking the helical ray in Eq. (38), with $a'=4$, $a=5$, $\tilde {\beta }=0.7071$, $n_{1}=n_{e}$, and $\Delta =0.2$; and Eqs. (30)–(32), the GRIN is reconstructed only at $\rho _{i}=a'$, as shown in Fig. 11. At first glance, the method apparently does not work in this case, when the ray is not a plane curve and lies in three dimensions. However, the considered ray path is a circular helix of radius $a'$, which implies that it is a ray that propagates over the isoindicial surface with value $n(a')$. In other words, a helical ray will always travel over the isoindicial surface where it impinges the GRIN distribution. Note that if we take the value of $\rho =a'$ in Eq. (33), we obtain that $n(a')=1.1903$, which is the same value that was obtained with our method using the helical ray.
Fig. 10. The helical ray enters the cylindrical waveguide at the $x-$ axis and exactly parallel to the $yz$ plane. From the $xy$ perspective, the entering ray appears tangent to the circle, so this selection of parameters yields a helical trajectory.
Fig. 11. Reconstruction of the GRIN distribution of the ray path that enters at a tangential plane and corresponds to the helical trajectory shown in Fig. 10. Notice that the reconstruction yields a hollow cylindric surface. The reconstruction in (a) 3D, (b) 2D, and (c) the profile as function of $\rho$.
The two cases presented in this section 4, show rays propagating in the $+z$ direction. In each case, the GRIN distribution has been calculated accordingly. It is possible to calculate the radial distribution by selecting a ray propagating in a plane transverse to the $z$ axis, i.e., $x$-$y$ plane. If we use this ray, it would be very easy to find the distribution using the same method for spherical GRIN medium that was presented in section 3.
5. Remarks on symmetric and freeform GRIN reconstruction
Now we discuss some remarks on the existence of a solution to the corresponding inverse problem of reconstructing GRIN media given a ray path.
For simplicity sake, let us begin our analysis with the spherical GRIN. In this case, the ray path must satisfy the restrictions imposed by the invariant quadratures, given by Eq. (12) in this case, which requires the ray paths to be planar curves. Then we can assume that the given ray path propagates in the $x$-$y$ plane without loss of generality. Consequently, the ray path can be parametrized as $\mathbf {r}(t)=(x(t), y(t))$. This leads as well to a parametrization of the radius $r$ of each point on the path; $r(t)=\sqrt {x(t)^2+y(t)^2}$.
Let us assume that the given ray path has at least two different points, $\mathbf {r}(t_1)$ and $\mathbf {r}(t_2)$, with the same radius, namely $r(t_1)=r(t_2)$. If the same refractive index is assigned to these points, that is $n(r(t_1))=n(r(t_2))$, and the same occurs for any other points having the same radius, then the symmetry holds, since all points with the same value of $r$ must lie in a isoindicial surface. A special case of this is when the parametrization $r(t)=\sqrt {x(t)^2+y(t)^2}$ assigns to distinct values of $t$, distinct values of $r$, that is, $r(t)$ is a one to one function. On the other hand, if different values of the GRIN are assigned to points with the same $r$, namely, $n(r(t_1))\neq n(r(t_2))$, then the spherical symmetry is broken, i.e., the GRIN becomes a freeform GRIN media. In this case, our method cannot be used and this problem should be addressed following a point-wise reconstruction, since the GRIN would be a function of more than one variable.
This reasoning can easily be generalized for skew rays propagating in symmetric GRIN media that depend on one single generalized coordinate $q$ (just as the examples addressed in this work); provided that the given ray path meets the invariant quadratures given by the symmetries of the system. Then, the inverse problem can be solved yielding a symmetric GRIN as a solution, provided that the same values of the index of refraction are assigned to points of a given isoindicial surface. In any other case, point-wise reconstruction methods should be used.
6. Conclusions
We have proposed a new method that solves the inverse problem of geometrical optics, which we have called Physical GRIN (PhysGRIN) method. Its introduction offers a valuable tool to reconstruct symmetric GRIN media since it only requires the invariants of the system and one single ray path. For these reasons, it is straightforward to implement an exact numerical method whose accuracy only depends on the numerical treatment. The method has an advantage over others since it does not require neither the knowledge of an explicit solution to any equation nor the use of machine learning.
The PhysGRIN method can be applied in other areas of physics, since it relies in the establishment of the invariants given by the symmetries of the Lagrangian system, as stated by Noether’s theorem. Direct applications can be given in Classical Mechanics and Acoustics, since both have the same mathematical structure as geometrical optics [54–59]. The PhysGRIN method could even be applied to gravitational symmetric systems reconstructing, for instance, the effective potential of a gravitational lens [60]. Also, this method may be generalized towards addressing freeform GRIN reconstruction. This shall be relevant in visual optics for characterizing the GRIN distribution of the crystalline lens, since it has an asymmetric GRIN distribution, i.e. its GRIN distribution is freeform. [4–6,61]. In conclusion, our work makes a significant contribution to the field by offering tools and techniques that can be applied across a wide range of situations and physical systems, thereby facilitating the analysis and comprehension of complex phenomena in diverse domains.
Funding
Instituto Nacional de Astrofísica, Óptica y Electrónica.
Acknowledgments
The authors acknowledge the support of the Instituto Nacional de Astrofísica, Óptica y Electrónica in this publication.
Disclosures
The authors declare no conflicts of interest.
Data availability
No data were generated or analyzed in the presented research.
Supplemental document
See Supplement 1 for supporting content.
References
1. M. V. Pérez and C. Bao, Gradient-Index Optics: Fundamentals and Applications (Springer, 2002).
2. E. W. Marchand, Gradient Index Optics (Academic Press, 1978).
3. D. H. Lippman, N. S. Kochan, T. Yang, G. R. Schmidt, J. L. Bentley, and D. T. Moore, “Freeform gradient-index media: a new frontier in freeform optics,” Opt. Express 29(22), 36997–37012 (2021). [CrossRef]
4. J. E. Gómez-Correa, S. E. Balderas-Mata, B. K. Pierscionek, and S. Chávez-Cerda, “Composite modified Luneburg model of human eye lens,” Opt. Lett. 40(17), 3990–3993 (2015). [CrossRef]
5. A. Jaimes-Nájera, J. E. Gómez-Correa, V. Coello, B. K. Pierscionek, and S. Chávez-Cerda, “Single function crystalline lens capable of mimicking ciliary body accommodation,” Biomed. Opt. Express 11(7), 3699–3716 (2020). [CrossRef]
6. A. Jaimes-Nájera, J. E. Gómez-Correa, V. Coello, B. K. Pierscionek, and S. Chávez-Cerda, “A single-function model for the eye’s crystalline lens,” Opt. Photonics News 31(12), 54 (2020). [CrossRef]
7. V. Lakshminarayanan, A. K. Ghatak, and K. Thyagarajan, Lagrangian Optics (Springer, 2002).
8. A. J. Dragt, E. Forest, and K. B. Wolf, “Foundations of a lie algebraic theory of geometrical optics,” in Lie Methods in OpticsK. B. Sánchez Mondragón and E. Wolf, eds. (Springer Berlin Heidelberg, 1986), pp. 105–157.
9. M. Born and E. Wolf, Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light (Cambridge University, 1999), 7th ed.
10. R. Jagannathan and S. A. Khan, “Quantum Mechanics of Charged Particle Beam Optics: Understanding Devices from Electron Microscopes to Particle Accelerators (CRC, 2019).
11. E. W. Marchand, “Ray tracing in gradient-index media,” J. Opt. Soc. Am. 60(1), 1–7 (1970). [CrossRef]
12. E. W. Marchand, “Rapid ray tracing in radial gradients,” Appl. Opt. 27(3), 465–467 (1988). [CrossRef]
13. E. W. Marchand, “Ray tracing in cylindrical gradient-index media,” Appl. Opt. 11(5), 1104–1106 (1972). [CrossRef]
14. D. T. Moore, “Ray tracing in gradient-index media,” J. Opt. Soc. Am. 65(4), 451–455 (1975). [CrossRef]
15. W. Streifer and K. B. Paxton, “Analytic solution of ray equations in cylindrically inhomogeneous guiding media. 1: Meridional rays,” Appl. Opt. 10(4), 769–775 (1971). [CrossRef]
16. K. B. Paxton and W. Streifer, “Analytic solution of ray equations in cylindrically inhomogeneous guiding media. part 2: Skew rays,” Appl. Opt. 10(5), 1164–1171 (1971). [CrossRef]
17. H. A. Buchdahl, “Rays in gradient-index media: separable systems,” J. Opt. Soc. Am. 63(1), 46–49 (1973). [CrossRef]
18. W. Liu, Y. Ren, and F. Liu, “Analytical ray tracing based on hamilton principal function and conjugate variable pairs,” J. Opt. Soc. Am. A 40(1), 78–84 (2023). [CrossRef]
19. A. Sharma, D. V. Kumar, and A. K. Ghatak, “Tracing rays through graded-index media: a new method,” Appl. Opt. 21(6), 984–987 (1982). [CrossRef]
20. A. Sharma, “Computing optical path length in gradient-index media: a fast and accurate method,” Appl. Opt. 24(24), 4367–4370 (1985). [CrossRef]
21. T. Sakamoto, “Ray trace algorithms for GRIN media,” Appl. Opt. 26(15), 2943–2946 (1987). [CrossRef]
22. J. E. Gómez-Correa, A. L. Padilla-Ortiz, A. Jaimes-Nájera, J. P. Trevino, and S. Chávez-Cerda, “Generalization of ray tracing in symmetric gradient-index media by fermats ray invariants,” Opt. Express 29(21), 33009–33026 (2021). [CrossRef]
23. J. E. Gómez-Correa, “Geometrical-light-propagation in non-normalized symmetric gradient-index media,” Opt. Express 30(19), 33896–33910 (2022). [CrossRef]
24. J. Evans and M. Rosenquist, ““F=ma” optics,” Am. J. Phys. 54(10), 876–883 (1986). [CrossRef]
25. S. A. Khan, “Quantum methods in light-beam optics,” Opt. Photonics News 27(12), 47 (2016). [CrossRef]
26. D. Ambrosini, A. Ponticiello, G. S. Spagnolo, R. Borghi, and F. Gori, “Bouncing light beams and the Hamiltonian analogy,” Eur. J. Phys. 18(4), 284–289 (1997). [CrossRef]
27. S. A. Khan, “Hamilton’s optical mechanical analogy in the wavelength-dependent regime,” Optik 130, 714–722 (2017). [CrossRef]
28. X. Hui, Z. He, and Y. Liu, “Design of cycloidal rays in optical waveguides in analogy to the fastest descent problem,” Opt. Continuum 2(8), 1763–1771 (2023). [CrossRef]
29. W. Liu, H. Hu, F. Liu, and H. Zhao, “Manipulating light trace in a gradient-refractive-index medium: a Lagrangian optics method,” Opt. Express 27(4), 4714–4726 (2019). [CrossRef]
30. U. Leonhardt, “Optical conformal mapping,” Science 312(5781), 1777–1780 (2006). [CrossRef]
31. S. Nemoto and T. Makimoto, “Refractive-index distribution for a prescribed ray path,” J. Opt. Soc. Am. 69(3), 450–454 (1979). [CrossRef]
32. J. C. Mi nano, “Refractive-index distribution in two-dimensional geometry for a given one-parameter manifold of rays,” J. Opt. Soc. Am. A 2(11), 1821–1825 (1985). [CrossRef]
33. F. Borghero and G. Bozis, “A two-dimensional inverse problem of geometrical optics,” J. Phys. A: Math. Gen. 38(1), 175–184 (2005). [CrossRef]
34. F. Borghero and G. Bozis, “Two solvable problems of planar geometrical optics,” J. Opt. Soc. Am. A 23(12), 3133–3138 (2006). [CrossRef]
35. F. Borghero and T. Kotoulas, “Refractive-index distributions generating as light rays a given family of curves lying on a surface,” J. Opt. Soc. Am. A 28(2), 278–283 (2011). [CrossRef]
36. P. S. J. Russell and T. A. Birks, “Hamiltonian optics of nonuniform photonic crystals,” J. Lightwave Technol. 17(11), 1982–1988 (1999). [CrossRef]
37. F. Borghero and F. Demontis, “Three-dimensional inverse problem of geometrical optics: a mathematical comparison between Fermat’s principle and the eikonal equation,” J. Opt. Soc. Am. A 33(9), 1710–1722 (2016). [CrossRef]
38. H. Ohno and T. Usui, “Neural network gradient-index mapping,” OSA Continuum 4(10), 2543–2551 (2021). [CrossRef]
39. H. Ohno and T. Usui, “Design method of gradient-index optics in homogeneous medium using neural network gradient-index mapping,” Opt. Rev. 30(3), pp355–360 (2023). [CrossRef]
40. H. Ohno, “Symplectic ray tracing based on Hamiltonian optics in gradient-index media,” J. Opt. Soc. Am. A 37(3), 411–416 (2020). [CrossRef]
41. R. K. Luneburg, Mathematical Theory of Optics (University of California, 1964).
42. J. E. Gómez-Correa, V. Coello, A. Garza-Rivera, N. P. Puente, and S. Chávez-Cerda, “Three-dimensional ray tracing in spherical and elliptical generalized Luneburg lenses for application in the human eye lens,” Appl. Opt. 55(8), 2002–2010 (2016). [CrossRef]
43. W. Liu, “Ray tracing in concentric gradient-index media: optical binet equation,” J. Opt. Soc. Am. A 39(6), 1025–1033 (2022). [CrossRef]
44. J. W. Blaker and M. A. Tavel, “The Application of Noether’s Theorem to Optical Systems,” Am. J. Phys. 42(10), 857–861 (1974). [CrossRef]
45. D. Neuenschwander, Emmy Noether’s Wonderful Theorem (Johns Hopkins University, 2017).
46. A. Ankiewicz and C. Pask, “Geometric optics approach to light acceptance and propagation in graded index fibres,” Opt. Quantum Electron. 9(2), 87–109 (1977). [CrossRef]
47. A. Ankiewicz, “Geometric optics theory of graded index optical fibres,” Ph.D. thesis, The Australian National University (1978).
48. J. A. Grzesik, “Focusing properties of a three-parameter class of oblate, Luneburg-like inhomogeneous lenses,” Journal of Electromagnetic Waves and Applications 19(8), 1005–1019 (2005). [CrossRef]
49. J. M. Gordon, “Spherical gradient-index lenses as perfect imaging and maximum power transfer devices,” Appl. Opt. 39(22), 3825–3832 (2000). [CrossRef]
50. K. Xu, Y. Chen, T. A. Okhai, and L. W. Snyman, “Micro optical sensors based on avalanching silicon light-emitting devices monolithically integrated on chips,” Opt. Mater. Express 9(10), 3985–3997 (2019). [CrossRef]
51. A. G. Leal-Junior, C. Marques, A. Frizera, and M. J. Pontes, “Multi-interface level in oil tanks and applications of optical fiber sensors,” Opt. Fiber Technol. 40, 82–92 (2018). [CrossRef]
52. A. G. Leal-Junior, A. Frizera, C. Marques, and M. J. Pontes, “Optical fiber specklegram sensors for mechanical measurements: A review,” IEEE Sens. J. 20(2), 569–576 (2020). [CrossRef]
53. K. Xu, “Silicon electro-optic micro-modulator fabricated in standard cmos technology as components for all silicon monolithic integrated optoelectronic systems*,” J. Micromech. Microeng. 31(5), 054001 (2021). [CrossRef]
54. H. Goldstein, Classical Mechanics (Addison-Wesley, 2011), 2nd ed.
55. P. C. Deshmukh, Foundations of Classical Mechanics (Cambridge, 2019).
56. L. Zhao, T. Horiuchi, and M. Yu, “Broadband acoustic collimation and focusing using reduced aberration acoustic Luneburg lens,” J. Appl. Phys. 130(21), 214901 (2021). [CrossRef]
57. L. Zhao, T. Horiuchi, and M. Yu, “Acoustic waveguide based on cascaded Luneburg lens,” JASA Express Lett. 2(2), 024002 (2022). [CrossRef]
58. C. Hu, S. Xue, Y. Yin, Z. Hao, Y. Zhou, and H. Chen, “Acoustic super-resolution imaging based on solid immersion 3D Maxwell’s fish-eye lens,” Appl. Phys. Lett. 120(19), 192202 (2022). [CrossRef]
59. J. Hyun, C.-S. Park, J. Chang, W.-H. Cho, and M. Kim, “Gradient-index phononic crystals for omnidirectional acoustic wave focusing and energy harvesting,” Appl. Phys. Lett. 116(23), 234101 (2020). [CrossRef]
60. W. Xiao, S. Tao, and H. Chen, “Mimicking the gravitational effect with gradient index lenses in geometrical optics,” Photonics Res. 9(7), 1197–1203 (2021). [CrossRef]
61. C. J. Sheil and A. V. Goncharov, “Accommodating volume-constant age-dependent optical (avocado) model of the crystalline grin lens,” Biomed. Opt. Express 7(5), 1985–1999 (2016). [CrossRef]
