1. INTRODUCTION

Symplectic ray tracing was born out of Hamiltonian optics, first described by William Rowan Hamilton in his 1828 paper, Theory of Systems of Rays [1]. Indeed, Hamilton’s formalism was then applied to dynamics more generally [2], where it is now a fundamental component of undergraduate courses in both classical and quantum mechanics [3,4]. While several volumes have been published specifically on Hamiltonian optics [5–7], Hamilton’s analysis is seldom exploited to its full potential, typically being neglected in favor of an equivalent Lagrangian framework, which has indeed proved useful for ray tracing in gradient-index (GRIN) media [8–10].

Nonetheless, by adopting the Hamiltonian perspective, we may construct numerical techniques for ray tracing in GRIN media that are both accurate and computationally inexpensive. Since the calculation of ray trajectories in GRIN lenses typically requires the solution of a hyperbolic nonlinear partial differential equation whose solution is analytical only in certain cases [11], numerical methods prove to be indispensable for nonlinear ray tracing [12].

Unfortunately, however, symplectic ray tracing has not been given much consideration over the past two decades, despite numerous advances being made in the underlying theory [13–15]. Furthermore, most papers examining the application of symplectic numerical techniques to optical problems require that the reader be familiar with differential geometry [16]. Additional publications demand a decent knowledge of group theory to fully appreciate a proposed Hamiltonian alternative to Seidel aberrations for the creation of freeform optical systems [17,18].

This paper instead aims to provide a useful working knowledge of symplectic integration techniques, thereby enabling these numerical methods to be used for GRIN ray tracing without requiring prior familiarity with differential geometry or group theory. In particular, we also include second-order methods in our analysis, which are often omitted in other publications, presumably on the assumption that fourth-order methods give more accurate results [16,19]. The remainder of this paper has six additional sections plus an appendix. Section 2 provides a brief overview of the mathematical preliminaries necessary for Hamiltonian optics. Section 3 then describes the construction of a symplectic numerical method via the splitting technique. Section 4 compares the accuracy of symplectic and nonsymplectic methods when ray tracing a selection of lenses with spherically symmetric index profiles. The computation times required by the most accurate methods identified in Section 4 are examined when rendering a test image in Section 5. A more detailed discussion takes place in Section 6. Section 7 concludes the paper. Appendix A contains numerical recipes for each integration method used in Section 4.

2. OVERVIEW OF HAMILTONIAN OPTICS

Numerous treatments of Hamiltonian optics begin by deriving an optical Lagrangian as a consequence of Fermat’s principle before performing a Legendre transform on this Lagrangian to obtain the associated optical Hamiltonian [20–22]. While interested readers are directed to [23] for a more thorough derivation, we dispense with this procedure for the sake of brevity, choosing instead to derive the necessary optical Hamiltonian directly from the eikonal equation [24]

Now, differentiating ${\textbf p}$ with respect to $s$

Recalling $|\nabla S{|^2} = {n^2}$ from Eq. (1), once again

Finally, multiplying both sides by $n$, we get

Defining ${\rm d}t = {\rm d}s/n$, we now reformulate Eq. (1) as a system of six ordinary differential equations

However, to avoid confusion, some caution is warranted. Although the parameter $t$ is analogous to time in a dynamical system and is even referred to as time in existing literature [13,25], it does not represent physical time. Nevertheless, by making this change of variables, the optical momenta simply become the direction cosines for their respective axes. Furthermore, by choosing to parameterize our system in terms of $t$, we transform Eq. (1) from a partial differential equation into a system of Hamilton’s equations for some optical Hamiltonian $H$. Thus, we may now write

Solving for $H$ in Eq. (12) is then a straightforward exercise that provides the optical Hamiltonian

The merits of Hamilton’s formulation of geometric optics become more apparent by examining the dissimilarities in the behavior of symplectic and nonsymplectic numerical methods. Figure 1 illustrates how a ray’s position and optical momentum calculated during each iteration are handled differently by each method. A nonsymplectic method will update a ray’s optical momentum and position using data obtained exclusively from the previous iteration, while a symplectic method initially calculates the ray’s optical momentum. Then, with this updated momentum, it obtains the ray’s position during the same iteration. Hence, the error of a symplectic method is bounded for a greater number of iterations when compared to an equivalent nonsymplectic one [14,15,26]. A useful consequence of this result is that symplectic numerical methods are typically more accurate than their nonsymplectic counterparts, particularly when using a relatively large step between successive iterations [14,15].

 

Fig. 1. Differences between symplectic and nonsymplectic methods, with ${\textbf p}$ and ${\textbf q}$ representing a ray’s optical momentum and position vectors, respectively. By using the momentum calculated during the same iteration to update the ray’s position, symplectic methods typically produce more accurate results.

Download Full Size | PDF

3. CONSTRUCTING A SYMPLECTIC NUMERICAL METHOD

The application of symplectic methods to ray tracing in GRIN media is made possible by the fact that Eq. (13) is of the form $H = P({p_x},{p_y},{p_z}) + N(x,y,z)$, making it separable. This allows the terms dependent on the ray’s position and optical momentum to be treated independently. Symplectic numerical methods for separable Hamiltonians are constructed via splitting, giving two separate iterative schemes to calculate the position and optical momentum of a given ray. Constructing the relevant iterative schemes first requires us to solve Hamilton’s equations in matrix form. First, we combine the position and momentum vectors for a given ray into a single matrix ${\textbf z} = ({\textbf p},{\textbf q})$ and define the Liouville operator ${D_H}$ [27] to act on ${\textbf z}$, so

Treating ${D_H}$ as a vector operator, the solution to Eq. (15) takes the form of a matrix exponential

Recalling that Eq. (13) is separable, if we further define

However, this does not exclude the possibility of negative coefficients in the case where the sums of $\{{c_i}\}$ and $\{{d_i}\}$ would otherwise be greater than unity. Forest and Ruth [28], for instance, provide a concrete example of a fourth-order method with negative coefficients and is considered in Appendix A. Next, carrying out first-order Taylor series expansions of $\exp (t{c_i}{D_P})$ and $\exp (t{d_i}{D_N})$, we get

In the simplest case of a first-order method (i.e., taking ${c_i} = {d_i} = 1$), if we know the state of a given ray ${{\textbf z}_n}$ at $t = n$, we may calculate ${{\textbf z}_{n + 1}}\,\,{\rm at}\,\,{t_{n + 1}} = {t_n} + \Delta t$, where $\Delta t$ is some small, finite difference between subsequent iterations, by multiplying Eq. (19) by ${{\textbf z}_n}$ before substituting the result into Eq. (20). We then obtain the symplectic Euler method [15]

Higher-order symplectic methods can then be constructed by choosing appropriate values for ${c_i}$ and ${d_i}$. However, finding the optimal values for these coefficients often involves a large system of equations that may not have an exact solution [28–30].

Table 1. Lenses Used in These Numerical Experiments and Their Respective Index Profiles^a

View Table | View all tables in this article

4. RESULTS OF NUMERICAL EXPERIMENTS

To test the validity of symplectic methods to solve optical problems, simulated traces with a monochromatic light source of wavelength $\lambda = 589.3\;{\rm nm}$ were carried out within a selection of spherical, radially symmetric lenses, each of radius ${R_0} = 1\;{\rm mm}$, all of which have been extensively studied in existing optical literature. The Luneburg [31] and Maxwell fish-eye [32] lenses are classical GRIN elements that have found applications as directional antennae for radio and microwave devices [33,34] and more recently at optical wavelengths [35]. While the Luneburg lens will direct a collimated bundle of rays to a single focus on its surface, the Maxwell fish-eye is described as an absolute optical instrument, theoretically capable of ideal point-to-point imaging in a three-dimensional space [36].

The Eaton lens [37] and the optical black hole (also referred to as a concentrator lens) [25] differ from the Luneburg and Maxwell fish-eye lenses with refractive index profiles that are singular at their centers, providing a significant challenge to numerical methods that may be used to traces these lenses. Thankfully, however, we may take advantage of analogies that exist between such singular GRIN lenses and systems examined in other fields; the ray path within the Eaton lens, for instance, is effectively identical to the trajectory of a body in the Kepler problem of celestial mechanics [38]. The Eaton lens is gradually finding applications as a beam deflector [39], while the optical black hole could prove useful for laboratory astrophysics [40].

Table 1 provides further details about each of the four lenses used during numerical experiments where six numerical methods were chosen to trace each lens. However, for the Eaton lens, the symplectic Chin-Chen 4A method [41] developed for celestial mechanics was also tested because of the Eaton lens’ behavior as an optical analogue for celestial mechanics. Table 2 provides information on the order and symplecticity of each method used, with their associated numerical recipes presented in Appendix A.

Table 2. Numerical Methods Used to Ray Trace Each Lens^a

View Table | View all tables in this article

Since the Luneburg and Maxwell fish-eye lenses directed all rays to a single focal point, it was possible to trace multiple rays and examine their behavior simultaneously. However, as the Eaton lens and optical black hole direct each individual ray to a unique focus, only a single ray was traced in these lenses. The initial ray heights for the single rays in the Eaton lens and optical black hole were ${y_0} = 0.58{R_0}$ and ${y_0} = 0.94{R_0}$, respectively. These values were chosen so each numerical method would be challenged to maintain its accuracy in the vicinity of the singularity and avoid producing erroneous unphysical results.

Figure 2 depicts numerical ray traces through the Luneburg and Maxwell fish-eye lenses, with the first and third columns depicting the traces in their entirety while the second and fourth columns present a detailed view of the focal points. The step-size values range from $\Delta t = \lambda $ to $\Delta t = 100\lambda $, where $\lambda = 589.3\;{\rm nm}$. The exact solution is depicted in each case by a continuous blue line. As one might expect, the disparity between the numerical methods decreases with a reduction in $\Delta t$ and an increase in method order $m$. For the Luneburg lens, all methods with the exception of the Euler method are indistinguishable from each other and the exact solution for $\Delta t = \lambda$. Significant discrepancies between each of the methods are not readily visible until $\Delta t = 100\lambda$, where each of the methods can be easily identified, though significant overlap is still observed between RK2, RK4, and Forest and Ruth’s method. For Maxwell’s fish-eye, however, both the Euler method and the symplectic Euler method deviate noticeably at $\Delta t = \lambda$, with both methods disappearing from the inset when $\Delta t$ is increased to $100\lambda$. Each of the second- and fourth-order methods then follow suit, accumulating greater errors also deviating significantly from the expected focal point as $\Delta t$ reaches $100\lambda$. RK2, in particular, appears to significantly underestimate the focal distance, while the velocity Verlet method displays the opposite behavior, overestimating the focal distance, thought not to the same extent. RK4 and Forest and Ruth’s method both appear to underestimate the focal distance by approximately same amount, indicating unsurprisingly that any differences between symplectic and nonsymplectic methods appear to diminish with an increase in method order for the same choice of step size.

 

Fig. 2. Numerical traces through Luneburg (left) and Maxwell fish-eye (right) lenses, each of radius ${R_0}$. The region shown in each right-hand column is marked on each trace in the left-hand column by a black rectangle and depicts the focal point in greater detail.

Download Full Size | PDF

For both the Luneburg and Maxwell fish-eye lenses, Fig. 4 indicates that symplectic methods tend to exhibit lower focus errors than their nonsymplectic counterparts of the same order, with the focus error being defined as the difference between the numerical and exact focus along the lens perimeter, divided by the lens circumference, normalized in terms of the lens radius ${R_0}$. The dashed horizontal line represents the Airy disk radius, providing a benchmark for diffraction-limited imaging. Consequently, any error values below this line may be considered exact. Most noticeably, the symplectic Euler method vastly outperforms the standard Euler method by several orders of magnitude for the Luneburg lens and approximately one order of magnitude for the Maxwell fish-eye. The initial conditions presented in Fig. 4 for the Luneburg and Maxwell fish-eye lenses are associated with the rays in Fig. 2 whose optical momentum experiences the greatest rate of change. Differences between the fourth-order methods are much less drastic, with both RK4 and Forest and Ruth’s method producing errors on the same order of magnitude. However, for second-order methods, the results are much more surprising; the velocity Verlet method is the only method for which the numerical trace is diffraction limited for all step sizes examined. It not only outperforms RK2, but also RK4 and Forest and Ruth’s method, particularly when $\Delta t \gt 5\lambda$ for the Luneburg lens and when $\Delta t = 2.5\lambda $ or $25\lambda$ for the Maxwell fish-eye, suggesting that the increased computational expense incurred by using higher-order methods is not necessarily justified by greater accuracy. This is especially true of the Luneburg lens when $\Delta t = \lambda$, where the velocity Verlet method is exact to within machine precision.

A similar analysis is presented for the Eaton lens and optical black hole in Fig. 3, where the symplectic methods once again generally appear to offer superior performance to nonsymplectic methods, particularly at larger step sizes. For the Eaton lens with $\Delta t = 10\lambda$, the symplectic Euler method not only outperforms the standard Euler method (which is only contained within the inset for $\Delta t = \lambda$) but also proves more accurate than RK2. For $\Delta t = 100\lambda$, both RK4 and the Chin-Chen 4A method appear to offer the most accurate solutions, suggesting that the use of methods developed outside optics also could be of benefit for numerical ray tracing. Although it does not dominate all other methods to the same extent as it did with the Luneburg lens and Maxwell fish-eye, the velocity Verlet method still appears to perform somewhat well, particularly when compared to RK2, with the velocity Verlet method focusing just outside the inset. In contrast, Forest and Ruth’s method performs somewhat poorly here, requiring a smaller step size to match the results produced by RK4 and the velocity Verlet method. Nonetheless, Forest and Ruth’s method does mimic their behavior at least in a qualitative sense by overestimating the focusing ability of the Eaton lens. The Chin-Chen 4A method instead follows trajectories similar to all other methods, where the focusing ability of the lens is underestimated. Thus, the combination of an overestimating method with an underestimating one would perhaps form the basis for a successful predictor–corrector scheme, although this would presumably come at an elevated computational cost, especially if the methods chosen are not symplectic.

 

Fig. 3. Numerical traces through the Eaton lens (left) and optical black hole (right), each of radius ${R_0}$. The region shown in each right-hand column is marked on each trace in the left-hand column by a black rectangle and depicts the focal point in greater detail.

Download Full Size | PDF

In contrast to the retroreflecting behavior of the Eaton lens, the singular index profile of the optical black hole directs rays along an infinite inward spiral. However, the shaded gray circle in Fig. 3 provides a practical stopping condition and represents the event horizon of a Schwarzschild black hole with radius $0.67{R_0}$. Despite the necessary truncation of the ray paths, some surprising results may still be observed. For instance, RK4 exhibits uncharacteristically poor performance across all step sizes. Moreover, it is surpassed by all symplectic methods and even the standard Euler method when $\Delta t = 100\lambda$. Unlike the Eaton lens, however, RK2 performs similarly to the velocity Verlet method across all step sizes, providing a marked difference in its ray tracing capabilities for the two singular-index lenses. Finally, the symplectic Euler method outperforms the nonsymplectic Euler method by a considerable margin, particularly for $\Delta t \ge 10\lambda$, further emphasising the gains in accuracy that may be made through the use of symplectic numerical methods, particularly if a low-order method is desired. Furthermore, as was the case with the Eaton lens, the combination of methods that both overestimate and underestimate the focal distance (such as the velocity Verlet and Forest and Ruth’s methods, for instance) could once again form a useful predictor–corrector method, enabling a larger step size to be used.

On the whole, symplectic methods demonstrate more consistent performance across both singular lenses. Although the higher-order methods do not necessarily dominate nonsymplectic methods to the same extent as they did for the Luneburg and Maxwell fish-eye lenses, they tend to at least perform as well as the more established nonsymplectic algorithms at smaller step sizes. This is especially true of the optical black hole, where Forest and Ruth’s method exhibits the lowest error of all methods for $\Delta t = \lambda$, while RK4’s focus error remains almost constant for all chosen step sizes, as depicted in Fig. 4. Differences in the error between the symplectic and standard Euler methods for both the Eaton lens and optical black hole follow a similar general trend of decreasing with the step size, although the error values tend to oscillate less than those of the Luneburg lens and Maxwell fish-eye. A noteable exception, however, is Forest and Ruth’s method with $\Delta t = 2.5\lambda$ for the optical black hole, where the error increases once again for $\Delta t = \lambda$. By constrast, for the Eaton lens, RK4 still appears to be the most accurate choice for most step sizes, being the only method that offers a diffraction-limited image for a step size greater than $10\lambda$. Nonetheless, differences in the accuracy between RK4 and Forest and Ruth’s method for the Eaton lens become almost negligible when $\Delta t = \lambda$. Additionally, the Chin-Chen 4A method appears to underperform at smaller step sizes, with the velocity Verlet method making it almost redundant by offering similarly accurate solutions while demanding less computational resources.

 

Fig. 4. Error in ray tracing the Luneburg, Maxwell fish-eye, Eaton, and optical black hole lenses for numerous step sizes. The error is defined as the distance between the exact and numerical focus along the lens perimeter, divided by the lens circumference, normalized in terms of the lens radius ${R_0}$. The dashed horizontal black line represents the diffraction limit in each instance. Values below this line may be considered exact.

Download Full Size | PDF

Nevertheless, performing calculations at quadruple rather than double precision would enable us to examine any trends further and ascertain whether or not RK4 continues to hold its own versus Forest and Ruth’s method for ray tracing within the Eaton lens and also check the constancy of its focus error when tracing the optical black hole. Likewise, any differences between RK2 and the velocity Verlet method within the optical black hole could also be considered in greater detail. Indeed, the impressive performance of the velocity Verlet method within the Luneburg and Maxwell fish-eye lenses could also be further tested by higher-precision numerical experiments.

5. APPLYING SYMPLECTIC METHODS TO IMAGE RENDERING

Although it may be some time before the GRIN lenses considered in the previous section are incorporated into consumer optical devices, symplectic ray tracing finds a more immediate application in nonlinear ray tracing for rendering computer images. Considering the consistent performance of both the velocity Verlet method and Forest and Ruth’s method in the previous section, we now examine their suitability for ray tracing a black hole test image, with RK4 given as a nonsymplectic comparison. The ray trace program, written in Python, was originally made available on GitHub [42], requiring some modifications to remove dependencies on deprecated functions and the inclusion of velocity Verlet and Forest and Ruth’s methods within the existing program.

Figure 5 shows 1080 p samples of the black hole test images. In a Cartesian coordinate system whose origin is at the center of the black hole’s event horizon, the observer is located at $(x,y,z) = (0,{r_S}, - 20{r_S})$, where ${r_S}$ is the black hole’s Schwarzschild radius (i.e., the radius of its event horizon). All image rendering took place on a Dell Inspiron 5570 computer with an Intel Core i7 CPU and 8 GB of DDR4 RAM. Each ray was then traced from the observer’s location for a total of 250 iterations with a step size of $0.16{r_S}$. The rendering workload was balanced as equally as possible between four threads, with each thread running in parallel, assigned to one of the four CPU cores. At first glance, no significant differences can be seen between each of the three images in Fig. 5. However, Fig. 6 presents a normalized pixel-wise difference map between both symplectic methods and RK4. Some noticeable differences between the velocity Verlet method and RK4 may then be observed, especially in the vicinity of the black hole’s event horizon and at the edge of its accretion disk. By contrast, any discrepancies between Forest and Ruth’s method and RK4 are substantially more difficult to distinguish. The radius of the event horizon in pixels for each image was then calculated, with the results displayed in Table 3. We see that the resolutions for RK4 and Forest and Ruth’s method are identical, each resulting in an angular resolution per pixel of 64.0041 arcseconds for our observer. The angular resolution per pixel for the velocity Verlet method is marginally lower at 64.4041 arcseconds, although Fig. 5 has already demonstrated that this reduction in resolution is practically negligible.

 

Fig. 5. Test images (1080 p) of a Schwarzschild black hole, showing no significant difference in image quality between the three numerical methods shown. The observer is located at $(x,y,z) = (0,{r_S}, - 20{r_S})$, where the origin of this coordinate system is at the center of the black hole’s event horizon and ${r_S}$ represents its Schwarzschild radius.

Download Full Size | PDF

 

Fig. 6. Normalized pixel-wise difference map for Forest and Ruth’s method and the velocity Verlet method when compared to RK4. Any differences between RK4 and the velocity Verlet method are significantly more noticeable than those of Forest and Ruth’s method.

Download Full Size | PDF

Table 3. Angular Resolution Per Pixel for Each of the Three Numerical Methods Used to Render 1080p Images as Observed from $(x,y,z) = (0,{r_S}, - 20{r_S})$^a

View Table | View all tables in this article

 

Fig. 7. Plot of the necessary trace time versus image height for images with 16:9 and 4:3 aspect ratios. Both Forest and Ruth’s method and the velocity Verlet method require less time than RK4 to ray trace the same test image.

Download Full Size | PDF

Symplectic methods further demonstrate their advantage in Fig. 7, where both the velocity Verlet and Forest and Ruth’s methods provide modest reductions in the computation times required to trace the test images at various resolutions with both 16:9 and 4:3 aspect ratios. The same computer, step size, number of iterations, and observer position used to render the 1080p images in Fig. 5 were again used to render each image here. Owing to a greater number of pixels overall, each 4:3 ratio image took longer to render than a 16:9 image at the same approximate resolution. Independent of the aspect ratio, we notice that RK4 consistently requires more time to trace the same image when compared to Forest and Ruth’s method and the velocity Verlet method, with RK4 taking over 30 min to trace the largest image, yet the velocity Verlet method needs only half that time to carry out the same task. Furthermore, the difference in the time required by these two methods appears to grow with a decrease in image resolution. Forest and Ruth’s method also offers a modest reduction in the time required to trace the test image when compared to RK4. However, the fact that its image quality is virtually identical to that of the velocity Verlet method makes it a somewhat redundant choice for our task. Nonetheless, this still proves that symplectic methods prove to be a more optimal choice for nonlinear ray tracing, especially, as in our case, where the program may be executed in parallel, since parallel implementations of symplectic methods should not need extensive modification of the serial code.

Although the use of the Python language for this ray tracer made its modification straightforward, additional speed increases could be achieved by using a just-in-time Python compiler such as Cython or Numba [43,44], or even rewriting the program completely in a compiled language such as C, with multiprocessing aspects being handled by well-established parallel computing APIs such as OpenMP or MPI [45,46]. The ray tracer would presumably also stand to benefit greatly from GPU acceleration. Since only a CPU was available for the numerical experiments presented here, performance benchmarking on a GPU-enabled system would provide an interesting comparison with our present results. Symplectic numerical methods could also be tested against fast marching methods [47], another type of algorithm commonly employed to solve the eikonal equation. Nevertheless, the existing Python ray tracer could still serve as a valuable educational resource to introduce students to topics that include nonlinear ray tracing, numerical integration, and high-performance computing.

6. DISCUSSION

In general, symplectic methods offer a consistently accurate performance when ray tracing GRIN elements and often outperform popular nonsymplectic numerical techniques such as RK4. Although optics currently lags behind fields such as celestial mechanics and molecular dynamics in its adoption of symplectic methods, the velocity Verlet method appears to provide a useful template for the development of symplectic numerical methods specifically for ray tracing. Indeed, more complex optical elements such as natural and artificial eyes, for which analytical ray traces may not exist, could indeed benefit from the application of tailor-made symplectic methods. In addition to standard GRIN elements, freeform GRIN (F-GRIN) optics offer another degree of freedom in optical design, taking advantage of a lack of rotational invariance as well as a GRIN profile that may or may not be continuous [48]. A selection of bespoke symplectic numerical methods could accelerate the design and prototyping of useful F-GRIN lenses. The ability to accurately trace arbitrary F-GRIN elements would hopefully, in turn, lead to developments in fabrication and metrology, addressing in particular a current lack of useful testing methods for conventional GRIN and F-GRIN elements alike.

Moreover, each of the symplectic methods considered in this paper are explicit methods, meaning each subsequent step relies on position and optical momentum data previously obtained. Several implicit symplectic methods, often renowned for their increased numerical stability, remain to be studied. However, these methods typically require a root-solving algorithm to calculate the relevant data during each iteration, substantially increasing the computational cost required in their implementation. Nevertheless, certain index profiles could exist for which implicit methods may be used without the need for a root solver, thus eliminating any extra computational expenses incurred. Numerical methods, symplectic or otherwise, with adaptive step sizes and their application to GRIN ray tracing is another area that has yet to be explored.

7. CONCLUSION

This work emphasizes the value of symplectic numerical techniques for ray tracing GRIN media, demonstrating their increased accuracy when compared to popular nonsymplectic methods. The velocity Verlet method shows particular promise, performing better than expected when ray tracing the Luneburg lens, Maxwell fish-eye, and optical black hole. The RK4 method still performs well when tracing the Eaton lens, though the performance of Forest and Ruth’s method is equally adequate for $\Delta t = \lambda$. In all, the results presented here make a strong case for the development of specialized symplectic methods for GRIN optics where analytical ray traces may not exist. However, it is also clear that the focus error for each ray is related to that ray’s initial height. Finding an expression, empirical or otherwise, to describe the rate of change of the focus error with respect to the ray height could indeed be an area for further study. Nonetheless, in addition to future applications in optical design, symplectic methods greatly reduce the computational cost of nonlinear ray tracing, which may also be of some use to those outside the scientific community, particularly for visual effects in computer games, movies, and other forms of entertainment media.