Abstract
Typically, as a means to obtain a less complicated ray tracing method on a gradient-index (GRIN) medium, a normalization is done. This normalization is based upon the fact that the values of the refractive index on the surface of the GRIN medium and the value of the refractive index medium where it is immersed are the same. In this paper, a Fermat’s-ray-invariants-based ray tracing method in a non-normalized GRIN medium is presented. This method is so simple that it is no longer necessary to perform normalization to reduce ray tracing complications in GRIN media. In order to observe its operation, the method is implemented in different GRIN media with rectangular, cylindrical, and spherical symmetry. Also, the method is implemented in two different axicon GRIN lenses. MATLAB codes for the implementations are provided as supplementary files. Finally, it is possible to observe that the Fermat’s ray invariant is still preserved even outside the lens. This implies that it is not necessary to apply Snell’s law when the ray leaves the GRIN medium, because the conservation of the Fermat’s ray invariant performs the refraction as if Snell’s law had been applied.
© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
It is well known that, in a homogeneous medium, rays propagate along straight lines [1,2]. However, for a gradient index (GRIN) medium, an inhomogeneous medium where the refractive index varies from point to point [3,4], the rays are propagated along curved paths. The profile of these curves are given by the analytical [5–12] or numerical solution of the ray equation [13–15] or by creating a layered medium where Snell’s law is calculated for each layer [16]. In turn, the solutions of the ray equation depend on the geometry of the isoindicial surfaces (surfaces of constant index) of the GRIN medium. This is because the ray equation is solved in a coordinate system which is more adequate for the isoindical surfaces geometry. For example, the ray equation is solved in cartesian coordinates if the isoindicial surfaces are planes perpendicular to the optical axis; in cylindrical coordinates if its surfaces are concentric cylinders about the optical axis; and in spherical coordinates if the isoindicial surfaces are concentric spheres. These three types of media have something in common: they have at least one axis of symmetry in its coordinate system. For this reason, these media can be called symmetric GRIN media.
In the literature, it can be noticed that the ray tracing in GRIN lenses with spherical isoindical surfaces has been of great interest because of their stigmatic properties [17–20]. However, in 2021, Gómez-Correa et al. proposed an exact general method to perform ray tracing in any symmetric GRIN media [21], and it is called Physical Ray Tracing (PRT). This method is based on the conservation of Fermat’s ray invariants, i.e., this method comes from applying the Lagrangian formalism to solve the problem of Fermat’s variational principle. It should be mentioned that the Lagrangian formalism together with the Hamiltonian formalism have been developed and implemented in different disciplines of physics [22]. In recent years, the analogy between classical/quantum mechanics and optics has taken on special interest, to the point of being developed [23–26] allowing that ray tracing algorithms [18,19,27] and a method to characterize optical systems and compute aberrations [28] to be implemented. Since this method is based on Fermat’s ray invariants, we provide an explanation on Fermat’s ray invariants.
Fermat’s principle states that out of the many paths that can connect two given points P and Q, the light ray would follow that path for which the optical pathlength between the two points is an extremum [29]. Mathematically, Fermat’s principle can be expressed as
To find a solution to the Fermat’s principle problem, it is necessary that the Lagrangian satisfies the Euler-Lagrange equations,
For optical systems with spherical, rectangular and cylindrical symmetry, ray tracing is always performed assuming that the system is immersed in a medium with a constant refractive index ($n_{m}$) that is equal to the surface refractive index ($n_{Sur}$) of the GRIN system, i.e. $n_{m}=n_{Sur}$ [29]. This implies that the value of the GRIN at the boundary of the system is normalized with respect to the surrounding refractive index [30,31]. This normalization is possible because it is considered that the ray tracing inside the system does not lose generality. However, this normalization limits rays outside the GRIN medium to propagate only in a medium with refractive index equal to $n_{Sur}$.
In this paper, we propose a ray tracing method in a symmetric GRIN medium immersed in a non-normalized media in order to avoid this limitation presented in practically all the ray tracing methods [17–21]. Normalizing the refractive index of the optical system surface with respect to the medium to obtain a less complicated ray trace is not necessary because of the simplicity of this method. The proposed method is based on the Fermat’s ray invariants, and its applicability and versatility is demonstrated with examples of ray tracing in GRIN media with rectangular, cylindrical, and spherical symmetry immersed in a non-normalized media. Furthermore, the method is implemented in an interesting example given by the axicon GRIN lenses. The conservation of the invariant is so powerful that it is preserved after leaving the GRIN medium, i.e., by using this method, it is not necessary to apply Snell’s law when the ray leaves the medium; the conservation of the Fermat’s ray invariant performs the refraction as if Snell’s law had been calculated. We call this method Non-normalized Physical Ray Tracing (N-NPRT).
2. Non-normalized physical ray tracing in spherical lenses
Because of its geometry, one of the most studied GRIN lenses is the spherical GRIN lens. The GRIN of this lens is described by a function $n(r)$ that only depends on its radial coordinate $r=\sqrt {x^{2}+y^{2}+z^{2}}$, where $x$, $y$, and $z$ are Cartesian coordinates. The function $n(r)$ represents a refractive index distribution. Its maximum value is located in the center of the lens and decreases from the center to a minimum on the lens surface [3,4].
It is very well known that the Optical Lagrangian in spherical coordinates $\left (r,\theta,\phi \right )$, in the equatorial plane $\left (\theta =\pi /2\right )$, is given by
where $\phi _{r}=\textrm {d}\phi /\textrm {d}r$, and its corresponding equation of Lagrange in the variable $\phi$ isThe value of $K$ can be found at the lens surface. If $R$ is the maximum radius of the lens, $n\left (R\right )$ is the value of the refractive index of the lens surface, and $\varphi _{i}$ is the calculated value of $\varphi$ at the surface of the lens, we have
Figure 1 shows that $n\left (R\right )=n_{Sur}$, and $\varphi _{i}=\pi -\theta _{Ref}$, where $\theta _{Ref}=\sin ^{-1}\left [\frac {n_{m}}{n_{Sur}}\sin \theta _{Inc} \right ]$. Substituting these values into Eq. (8), and considering that $\sin \left (\pi \pm A\right )=\mp \sin A$, we obtain that where $\theta _{Inc}=\theta _{r}+\theta _{i}$ (see Fig. 1 and Fig. 2).Notice that Eq. (9) contains the value of the Fermat’s ray invariant in a non-normalized spherical medium. This equation is crucial, because, by using this equation, PRT in a lens embedded in a medium with $n_{m}\neq n_{Sup}$ can be performed, that is, an N-NPRT.
2.1 Non-normalized physical ray tracing stages
The N-NPRT can be implemented in any programming language by complying the following steps:
- (1) Define the values of $R$, $n\left (r\right )$, $n_{m}$, $P_{0}$, $P_{i}$, and $\Delta _{s}$. From Fig. 2, it is easy to see that $\Delta _{s}$ is the length of a straight line that has an angle of inclination $\gamma$ with respect to the horizontal axis, $P_{0}=(x_{0},y_{0})$ is the point where the incident ray originates, and $P_{i}=(x_{i},y_{i})$ indicates the point where the incident ray enters at the spherical lens.
- (2) Calculate the value of $K$ by means of Eq. (9).
- (6) The point $P_{i}$ is redefined as $P_{i}=P_{i+\Delta }$.
- (7) Repeat steps 3 to 6 until the desired point is reached.
Note that step 7 does not impose any restriction on the ray propagation, i.e., it is not confined to propagate only inside the lens. There is no restriction because the Fermat’s ray invariant is still preserved even outside the lens, i.e. it is not necessary to apply Snell’s law when the ray leaves the lens. The conservation of Fermat’s ray invariant performs the refraction as if Snell’s law had been calculated. In the most general case, as shown in Fig. 1, when the ray comes out after propagating inside the lens, the value of $\gamma$ is $\theta _{p}$, where $\theta _{p}=\gamma _{Ref}-\theta _{f}$, and $\gamma _{Ref}=\sin ^{-1}\left [\frac {n_{Sup}}{n_{m}}\sin \gamma _{Inc} \right ]$.
2.2 Implementation of the non-normalized physical ray tracing
The most famous spherical GRIN lens is the Luneburg lens. This lens has the wonderful capacity to focus the light that comes from infinity on its surface, i.e. it is a stigmatic lens [32]. This lens was introduced by Rudolf K. Luneburg in 1944. And its GRIN distribution is given by
where $0\leq r\leq R$.In 1955, S. Gutman proposed a new lens that is a generalization on the Luneburg lens [33]. This lens has the capacity to focus light within it. Its GRIN distribution is given by
where $f$ is the focal point. It is very important to say that if $0\leq f\leq R$, the Gutman lens is a stigmatic lens [33]. However, if $f>R$ the Gutman lens has a spherical aberration [34,35], i.e., the lens is not stigmatic, as shown in Fig. 3(a). Notice that, if $f=R$ the GRIN distribution (Eq. (13)) of the Gutman lens becomes a GRIN distribution of the Luneburg lens (Eq. (12)). It is also important to note that the Gutman lens is immersed in a normalized medium, that is, $n_{m}=n_{Sur}=1$.In 2021, Liuxian Zhao et al. design a Luneburg lens that reduces the spherical aberration produced by a Gutman lens with $f>R$ [36]. To minimize the aberration, a correction factor $\alpha$ that multiplies the radial coordinate is introduced in Eq. (13), as shown in Eq. (14),
The value of $\alpha$ can be determined based on an optimization function that reduce the ray aberration. This function isRay tracing in Fig. 3(b) has been done using the N-NPRT, because the lens is immersed in a non-normalized medium. The parameter values of the N-NPRT to perform ray tracing on the Modified Luneburg Lens are $R=1$, $n_{m}=1$, $n_{Sur}=1.0562$, $P_{0}=\left (-1000,0\right )$, $\Delta _{s}=0.0011$, and $n\left (r\right )$ is given by Eq. (14), as we can see in Code 1, Ref. [37]. By using all these values, we obtain an average value of the interception of the ray with the $x$-axis of $1.8766$. This means that the Focal distance is $F\approx 1.87$, as reported by Liuxian Zhao et al. in their paper [36].
3. Non-normalized physical ray tracing in axial refractive index lenses
The N-NPRT can also be implemented in other symmetric GRIN media. For example, let us consider an axial GRIN media represented by a function $n=n(y)$, i.e., the isoindical surfaces are planes that are parallel to the propagation axis. Its Optical Lagrangian is given by
where $y_{x}=\textrm {d}y/\textrm {d}x$. And its corresponding equation of Lagrange, in its alternative form [38], for the variable $y$ isNext, let us calculate the value of $K_{L}$ without normalization conditions. From Fig. 4, it is possible to show that $\varphi _{L}$ on the surface is $\varphi _{Li}=\theta _{Ref}$. Considering that $\theta _{Ref}=\sin ^{-1}\left [\frac {n_{m}}{n_{Sup}}\sin \theta _{Inc} \right ]$, then
To calculate the value of $\varphi _{L}$ we substituted Eq. (20) in Eq. (19) and we obtained that
By using the N-NPRT method described in section 2, but calculating the value of Fermat’s ray invariant this time, step 2, by means of Eq. (20), and the value of $\varphi _{L}$, step 3, by means of Eq. (21), it was possible to generate a ray tracing in an axial GRIN distribution immersed in a non-normalized medium ($n_{m}\neq n_{Sur}$) represented by
4. Modified physical ray tracing in a cylindrical fiber with parabolic GRIN
Let us consider an optical fiber with parabolic refractive index given by
Notice that the parabolic refractive index is represented by a function $n=n\left (r\right )$, and its isoindicial surfaces are concentric cylinders, i.e., this GRIN possesses a radial cylindrical profile. It is known that, for a profile of this type, its optical Lagrangian is given by
where $\phi _{r}=\textrm {d}\phi /\textrm {d}r$, and $z_{r}=\textrm {d}z/\textrm {d}r$. Its corresponding equation of Lagrange in the variable $z$ isIt is possible to find the values of the invariants at the moment the ray enters the GRIN medium at any point. Let us consider a ray propagating on the incident plane with a propagation angle $\theta _{ri}$ and enters the GRIN medium at the point $P_{i}=(x_{i},y_{i},0)$, (see Fig. 6), namely let us consider the most general possible propagation case. By using the geometry of Fig. 6,
Without loss of generality for the non-normalized system, we can make the value of the invariant $\tilde {l}$, as calculated by Eq. (31), do not depend on the value of $\vartheta$. This is possible if the incident plane is parallel to the $y$-$z$ plane, located at a height $x=r'$, that is, the ray is launched on the x-axis (at $x=r'$) in the $y$-$z$ plane (making angle $\theta _{ri}$ with the $z$-axis) [29]. This implies that $\sin {\left (\phi +\phi '\right )}=1$, then
Notice that the N-NPRT method (subsection 2.1) can be used to perform ray tracing on the cylindrical gradient described by Eq. (23). However, this method must undergo some modifications because the propagated rays in this GRIN are skew rays, and they are not contained in a plane. The modifications are:
By implementing these modifications in the N-NPRT method it is possible to obtain a ray trace for the cylindrical GRIN, as shown in Fig. 7. The N-NPRT was obtained by means of the parameters $n_{m}=1.1$, $n(r=0)=1.38$, $\Delta =0.2$, $a=5$, $\Delta _{s}=0.1$, $\theta _{ri}=45^{\circ }$, $\theta _{Inc}=45^{\circ }$, as we can see in Code 3, Ref. [42].
On the other hand, Fig. 6 shows that if the incident ray is launched on the $x$-axis ($x=r_{0}$) in the $x$-$z$ plane, (the incident plane is the $x$-$z$ plane), then $\sin {\left (\phi +\phi '\right )}=0$ and $\tilde {l}=0$. This implies that the skew rays become rays that are contained in a plane, i.e., the ray tracing problem in a cylindrical GRIN medium becomes a ray tracing problem in an axial medium, like that presented in section 3. It is important to see that the difference between the ray tracing in the meridional plane of the cylindrical GRIN and the axial GRIN is that, in the first, the gradient is given in an interval of $\left [-a,a\right ]$ whereas in the second, the gradient is in an interval of $\left [0,a\right ]$. Consequently, it is possible to realize that
and that it is also possible to perform a ray tracing in the meridional plane of the cylindrical GRIN by using the N-NPRT method described in section 3, as shown in Fig. 8.The N-NPRT shown in Fig. 8 was obtained by setting the parameters as follows: $n_{m}=1.1$, $n(r=0)=1.3$, $\Delta =0.2$, $a=5$, $\Delta _{s}=0.0027$, $\theta _{r}=19.7989^{\circ }$, $\theta _{Inc}=70.2011^{\circ }$, as we can see in Code 2, Ref. [39].
5. Axicon gradient lenses
An interesting study case is the axicon gradient lenses [43]. In 1954, H. McLeod defined an axicon as a figure of revolution that has the property that a point source when is refracted or reflected through its axis of revolution generates a continuous line of focus [44]. We could think that an axicon gradient lens is a GRIN lens that generates this type of line of focus, like the spherical GRIN lens introduced by J. R. Flores in 1999 [45]. However, the term "axicon gradient index" is not used in this same sense, this term is the name that the radial gradient medium takes when its index profile equation has a linear term in the radial distance from the axis [43,45,46], i.e., when its GRIN function can be written as
where $n_{0}$ is the maximum value of refractive index, and $\nu$ is a constant. From Fig. 9, it is possible to see that that this kind of medium, given by Eq. (37), presents a local minimum at the axis of the medium ($r=0$), and it presents a local maximum at $r=a_{0}$.Notice that, throughout this paper, the ray tracing has been performed on GRIN media that have their local maximum on the propagation axis. Even though the GRIN medium of the axicon gradient lens does not have the same distribution, however, this is not an impediment for the proposed method to work. In these types of lenses, cylindrical coordinates are used and the ray tracing analysis is performed in the meridional plane [43,46]. Using these same conditions, it is possible to generate a ray tracing of an axicon gradient lens using the N-NPRT method described in section 3. For instance, using Eq. (37) with $n_{0}=1.5$, $\nu =0.3758$, and $a_{0}=1$, it is possible to generate an N-NPRT shown in Fig. 10(a). The parameters used in this ray tracing are $a=2$, $\Delta _{s}=0.00001$, $P_{0}=\left (-600000,0\right )$, and $n_{m}=n(y_{i})$ where $y_{i}$ is the incident height of each ray, i.e., $n_{m}$ has the same value as the refractive index of the lens at the incident point, as we can see in Code 4, Ref. [47].
From Fig. 10(a), we can observed that the GRIN function given by Eq. (37) leads to a non-periodic ray tracing when the rays coming from infinity. This is because the GRIN profile has a linear term and is not paraxial and generates spherical aberration [43,46]. However, in 2020, J. Guénette et al. propose a GRIN function that generates a perfectly periodic ray tracing when the rays coming from infinity [48]. This exact GRIN funtion is given by
where the $\alpha$ parameter is expressed as [49]The GRIN distribution of Eq. (38) is very similar to the distribution given Eq. (37) as shown in Fig. 9. However, the small differences generate a perfectly periodic ray tracing, as shown in Fig. 10(b). The N-NPRT was implemented in this distribution using the same parameters for the ray tracing of the GRIN distribution of Eq. (37). The value of $\alpha$ was obtained with $n(r=0)=1.4$, i.e., the value is $\alpha =\nu =0.3758$.
6. Discussion
In specialized literature, the expression of $K$ given by Eq. (7) is known as the generalized Snell’s law for inhomogeneous media with radial symmetry [29]. However, it is difficult to see why it is considered as a generalized Snell’s law because it is in terms of the angle $\varphi$, the position vector magnitude ($r=\|\overrightarrow {r}\|$), and the value of the refractive index at $\overrightarrow {r}$. In contrast, it is easy to see in Eq. (9), which is given by terms of Snell’s law, that is, in terms of the refractive index where the ray is propagating $n_{m}$ (or $n(r)$ inside the lens), and the angle of incidence $\theta _{Inc}$ which is measured with respect to the normal of the isoindical surfaces of the GRIN medium. Furthermore, for the case of the axial GRIN medium, it can be seen Eq. (20) that $K_L$ (Eq. (19)) is a generalized Snell’s law for inhomogeneous media with axial symmetry. And, for the case of the cylindrical GRIN medium, it is evident in Eq. (30) that $\tilde {\beta }$ is a generalized Snell’s law for inhomogeneous media with cylindrical symmetry and $\tilde {l}$ is in terms of $\tilde {\beta }$ (Eq. (31)). Because of this, we affirm that the Fermat’s ray invariants are generalized Snell’s laws. Notice that, by preserving the invariants along the ray path, what is really preserved is Snell’s law inside and outside the GRIN medium. This is the reason for calculating the entire path of the rays and also why it is not necessary to calculate Snell’s law when the ray leaves the GRIN medium to a medium with constant refractive index, that is, we are always implicitly calculating Snell’s law by conserving the invariants.
7. Conclusions
An easy-to-implement method, for a ray tracing through symmetric GRIN media immersed in a non-normalized media by means of Fermat’s ray invariants was presented. This method takes advantage over other methods because it eliminates the limitation presented in practically all the ray tracing methods: the rays outside the GRIN medium only propagate in a medium with refractive index equal to that of the surface of the GRIN medium. Its applicability and versatility was observed by performing ray tracing on axicon GRIN lenses and GRIN media with rectangular, cylindrical and spherical symmetry immersed in a non-normalized medium. In addition, in this work, we observed that the Fermat’s ray invariants are generalized Snell’s laws that allow to know the path of the ray inside and outside the GRIN medium, that is, by conserving the invariants, we conserve Snell’s law at each point of the ray path, which allowed us to do the refractions when we went from a GRIN medium to a medium with a constant refractive index.
This method can be of particular interest in the area of visual optics because, by implementing some changes, this method can be used to make ray tracing in a simpler way in the human eye. It is important to remember that the ray tracing inside of the human eye is not easy because the crystalline lens, which is the internal lens of the human eye and it is a gradient-index lens, is immersed in a medium with a refractive index than different its surface refractive index [50–52].
Funding
Centro de Investigación Científica y de Educación Superior de Ensenada, Baja California.
Acknowledgments
The author would like to acknowledge A. L. Padilla-Ortiz for fruitful discussions and for her help in improving this paper. I dedicate this work to my mother and in memory of my father.
Disclosures
The author declares no conflicts of interest.
Data availability
MATLAB codes for the implementations in different coordinate systems are available in Code 1, Ref. [37], Code 2, Ref. [39], Code 3, Ref. [42], and Code 4, Ref. [47].
References
1. M. Born and E. Wolf, Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light (Cambridge University Press, 1999), 7th ed.
2. R. Jagannathan and S. A. Khan, Quantum Mechanics of Charged Particle Beam Optics: Understanding Devices from Electron Microscopes to Particle Accelerators (CRC Press, 2019).
3. C. Gómez-Reino, M. V. Pérez, and C. Bao, Gradient-Index Optics: Fundamentals and Applications (Springer, 2002).
4. E. W. Marchand, Gradient Index Optics (Academic Press, 1978).
5. E. W. Marchand, “Ray tracing in gradient-index media,” J. Opt. Soc. Am. 60(1), 1–7 (1970). [CrossRef]
6. E. W. Marchand, “Rapid ray tracing in radial gradients,” Appl. Opt. 27(3), 465–467 (1988). [CrossRef]
7. E. W. Marchand, “Ray tracing in cylindrical gradient-index media,” Appl. Opt. 11(5), 1104–1106 (1972). [CrossRef]
8. D. T. Moore, “Ray tracing in gradient-index media,” J. Opt. Soc. Am. 65(4), 451–455 (1975). [CrossRef]
9. 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]
10. 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]
11. H. A. Buchdahl, “Rays in gradient-index media: separable systems,” J. Opt. Soc. Am. 63(1), 46–49 (1973). [CrossRef]
12. J. A. Grzesik, “Focusing properties of a three-parameter class of oblate, Luneburg-like inhomogeneous lenses,” J. Electromagn. Waves Appl. 19(8), 1005–1019 (2005). [CrossRef]
13. 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]
14. A. Sharma, “Computing optical path length in gradient-index media: a fast and accurate method,” Appl. Opt. 24(24), 4367–4370 (1985). [CrossRef]
15. T. Sakamoto, “Ray trace algorithms for GRIN media,” Appl. Opt. 26(15), 2943–2946 (1987). [CrossRef]
16. D. A. Atchison and G. Smith, Optics of the Human Eye (Butterworth-Heinemann, 2002).
17. W. Liu, “Ray tracing in concentric gradient-index media: optical binet equation,” J. Opt. Soc. Am. A 39(6), 1025–1033 (2022). [CrossRef]
18. H. Ohno and T. Usui, “Points-connecting neural network ray tracing,” Opt. Lett. 46(17), 4116–4119 (2021). [CrossRef]
19. H. Ohno, “Symplectic ray tracing based on Hamiltonian optics in gradient-index media,” J. Opt. Soc. Am. A 37(3), 411–416 (2020). [CrossRef]
20. 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]
21. 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 fermat’s ray invariants,” Opt. Express 29(21), 33009–33026 (2021). [CrossRef]
22. K. Gambár and F. Márkus, “Hamilton-lagrange formalism of nonequilibrium thermodynamics,” Phys. Rev. E 50(2), 1227–1231 (1994). [CrossRef]
23. J. Evans and M. Rosenquist, ““F=ma” optics,” Am. J. Phys. 54(10), 876–883 (1986). [CrossRef]
24. S. A. Khan, “Quantum methods in light-beam optics,” Opt. Photonics News 27(12), 47 (2016). [CrossRef]
25. 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]
26. S. A. Khan, “Hamilton’s optical–mechanical analogy in the wavelength-dependent regime,” Optik 130, 714–722 (2017). [CrossRef]
27. 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]
28. A. J. Dragt, E. Forest, and K. B. Wolf, “Foundations of a lie algebraic theory of geometrical optics,” in Lie Methods in Optics, J. Sánchez Mondragón and K. B. Wolf, eds. (Springer Berlin Heidelberg, 1986), pp. 105–157.
29. V. Lakshminarayanan, A. K. Ghatak, and K. Thyagarajan, Lagrangian Optics (Springer, 2002).
30. J. R. Flores, “Estudio de elementos ópticos de gradiente de índice de simetría esférica,” Ph.D. thesis (Universidade de Santiago, 1992).
31. J. E. Gómez-Correa, “Complete gradient refractive index lens schematic human eye model,” Ph.D. thesis (Instituto Nacional de Astrofísica, Óptica y Electrónica, 2015).
32. R. K. Luneburg, Mathematical Theory of Optics (University of California Press, 1964).
33. A. S. Gutman, “Modified Luneberg lens,” J. Appl. Phys. 25(7), 855–859 (1954). [CrossRef]
34. S. P. Morgan, “General solution of the Luneberg lens problem,” J. Appl. Phys. 29(9), 1358–1368 (1958). [CrossRef]
35. J. A. Lock, “Scattering of an electromagnetic plane wave by a Luneburg lens. i. ray theory,” J. Opt. Soc. Am. A 25(12), 2971–2979 (2008). [CrossRef]
36. 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]
37. J. E. Gómez-Correa, “Numerical implementation of non-normalized physical ray tracing in a spherical GRIN medium,” figshare (2022). https://doi.org/10.6084/m9.figshare.20047073
38. G. B. Arfken and H. J. Weber, Mathematical Methods for Physics; 6th ed. (Academic Press, 2005).
39. J. E. Gómez-Correa, “Numerical implementation of non-normalized physical ray tracing in an axial GRIN medium,” figshare (2022). https://doi.org/10.6084/m9.figshare.20047091
40. 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]
41. A. Ankiewicz, “Geometric optics theory of graded index optical fibres,” Ph.D. thesis (The Australian National University, 1978).
42. J. E. Gómez-Correa, “Numerical implementation of non-normalized physical ray tracing in a cylindrical GRIN medium,” figshare (2022). https://doi.org/10.6084/m9.figshare.20047112
43. E. W. Marchand, “Axicon gradient lenses,” Appl. Opt. 29(28), 4001–4002 (1990). [CrossRef]
44. J. H. McLeod, “The axicon: A new type of optical element,” J. Opt. Soc. Am. 44(8), 592–597 (1954). [CrossRef]
45. J. R. Flores, “Gradient-index axicons with spherical symmetry,” J. Mod. Opt. 46(10), 1513–1525 (1999). [CrossRef]
46. R. M. González, J. L. nares, and C. Gómez-Reino, “Gradient-index axicon lenses: a quasi-geometrical study,” Appl. Opt. 33(16), 3420–3426 (1994). [CrossRef]
47. J. E. Gómez-Correa, “Numerical implementation of non-normalized physical ray tracing in an axicon gradient lens,” figshare (2022). https://doi.org/10.6084/m9.figshare.20317311
48. J. Guénette, M. Quémener, P. Bélanger, and S. Thibault, “Exact graded-index profile to produce an annular beam,” Opt. Commun. 472, 126037 (2020). [CrossRef]
49. M. Quémener, N. Grégoire, S. Morency, D. Côté, and S. Thibault, “Mcvd-based grin-axicon for the generation of scalable bessel–gauss beams,” Opt. Lett. 46(6), 1333–1336 (2021). [CrossRef]
50. 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]
51. 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]
52. 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]