Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group
Journal Home About Tutorials Issues in Progress Current Issue All Issues Early Posting Feature Issues

Symplectic ray tracing based on Hamiltonian optics in gradient-index media

Hiroshi Ohno

Open Access Open Access

Abstract

A method based on Hamiltonian optics for ray tracing through gradient-index (GRIN) media is proposed. The ray equation that describes light-ray paths can be written in the form of the Hamiltonian equations. Although the Hamiltonian equations can be numerically calculated using a finite-difference explicit method, deviations from the exact equations are generally inevitable at subsequent time steps. An optical Hamiltonian can be constructed of two independent terms, i.e., one term dependent on position and the other term dependent on momentum. The symplectic integrator is applicable to such a separable optical Hamiltonian system and makes the optical Hamiltonian equations form invariant at each time step of numerical calculations. Accuracies of light-ray paths calculated using the first-order symplectic ray tracing in GRIN lenses approximate those calculated on the basis of the fourth-order Runge–Kutta algorithm, which shows the promising potential of the symplectic-ray-tracing method.

© 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

1. INTRODUCTION

Light-ray paths can be controlled by refractive index fields. Since novel design techniques of gradient-index (GRIN) media based on the coordinate transformation have rapidly been developed, new types of GRIN media have been created [1–6]. In designing such GRIN media, there is a growing need for a highly efficient and practical ray-tracing method. Physical quantities such as temperature, pressure, strain, and density in media can generate GRIN fields. Several methods to measure the distributions of such physical quantities in transparent media have been rapidly developed by means of measuring deflection angles of light rays passing through the generated GRIN fields [7,8]. These measuring methods also require a highly efficient and practical ray-tracing method.

The ray equation that is derived based on Fermat’s principle can describe light-ray paths in a GRIN medium. However, because the light-ray paths in the GRIN medium are arbitrarily curved in general, the ray equation should be solved numerically except for a few cases in which analytical solutions are possible. Several ray-tracing techniques based on finite-difference explicit methods have been developed for GRIN media [9–13]. In particular, numerical integration based on the Runge–Kutta algorithm uses a standard numerical technique, is practically treatable, and has proven to be a faster and more accurate numerical method than the other methods [14].

The ray equation can be transformed into the Hamiltonian equations. Several forms of the Hamiltonian equations have been extensively studied [15]. Although the Hamiltonian equations can be numerically calculated using finite-difference explicit methods, deviations from their exact equation are inevitable at subsequent time steps. This causes the accumulation of error. The symplectic integrator that is able to preserve the equation forms has therefore been developed for celestial mechanics, nonlinear Hamiltonian systems, and quantum physics [16–19]. The symplectic integrator has also revealed its worth in optical Hamiltonian systems to increase accuracies of numerical calculations of the finite-difference time-domain (FDTD) method and wavelet transformation [20–22]. A method based on the symplectic integrator algorithm for ray tracing in GRIN media is thus proposed here. The optical Hamiltonian can be constructed of two independent terms, i.e., one term dependent on position and the other term dependent on momentum. The symplectic integrator is applicable to such a separable optical Hamiltonian system and ensures that the optical Hamiltonian equations are almost form invariant at each time step of numerical calculations. This should increase accuracy of the ray tracing. The light-ray paths in the cylindrical and spherical GRIN lenses are calculated using the symplectic ray tracing method and are compared with those calculated on the basis of the fourth-order Runge–Kutta algorithm. It is shown that the deviations of light-ray paths from ideal paths are supremely reduced by the symplectic ray-tracing method.

2. SYMPLECTIC-RAY-TRACING METHOD BASED ON HAMILTONIAN OPTICS

The path of light rays in a GRIN medium can be described by the ray equation that can be derived from the Fermat’s principle as

$$\frac{d}{{ds}}\left[ {n\left( {\boldsymbol q} \right)\frac{{d{\boldsymbol q}}}{{ds}}} \right] = \nabla n\left( {\boldsymbol q} \right),$$
where ${\boldsymbol q}$ is a position vector on the light-ray path, $n({\boldsymbol q})$ is a refractive index that depends on the position ${\boldsymbol q}$, and $ds$ is an infinitesimal element of the length along the path. Here, we define a light-ray momentum ${\boldsymbol p}$ as
$${\boldsymbol p} = n\left( {\boldsymbol q} \right)\frac{{d{\boldsymbol q}}}{{ds}}.$$
The ray equation of Eq. (1) can be transformed into the Hamiltonian equations with some Hamiltonian $H$ as
$$\frac{{d{\boldsymbol p}}}{{dt}} = - \frac{{\partial H}}{{\partial {\boldsymbol q}}},$$
$$\frac{{d{\boldsymbol q}}}{{dt}} = \frac{{\partial H}}{{\partial {\boldsymbol p}}},$$
where $t$ is a parameter and is hereafter called time even when it is not the physical time.

A point ${\boldsymbol a}$ in the phase space can be written with ${\boldsymbol q}$ and ${\boldsymbol p}$ as

$${\boldsymbol a}\left( t \right) = \left[ {\begin{array}{c}{{\boldsymbol q}\left( t \right)}\\{{\boldsymbol p}\left( t \right)}\end{array}} \right].$$
The Hamiltonian equations corresponding to Eqs. (3) and (4) can then be written with the phase space point ${\boldsymbol a}$ as
$$\frac{{d{\boldsymbol a}}}{{dt}} = {\boldsymbol J}\frac{{\partial H}}{{\partial {\boldsymbol a}}},$$
where ${\boldsymbol J}$ is a matrix defined as
$${\boldsymbol J} = \left[ {\begin{array}{cc}\textbf{0}&\textbf{1}\\{ - \textbf{1}}&\textbf{0}\end{array}} \right].$$

Another phase space point at a next time step, ${\boldsymbol A}$, can be written with an infinitesimal time $\Delta t$ as

$${\boldsymbol A} = \left[ {\begin{array}{c}{\boldsymbol Q}\\{\boldsymbol P}\end{array}} \right] = {\boldsymbol a}\left( {t + \Delta t} \right) = \left[ {\begin{array}{c}{{\boldsymbol q}\left( {t + \Delta t} \right)}\\{{\boldsymbol p}\left( {t + \Delta t} \right)}\end{array}} \right].$$
The phase space point ${\boldsymbol A}$ should also satisfy Eq. (6), which can be written as
$$\frac{{d{\boldsymbol A}}}{{dt}} = {\boldsymbol J}\frac{{\partial H}}{{\partial {\boldsymbol A}}}.$$
Using Eqs. (6) and (9), the following equation can be derived as
$$\frac{{d{A^i}}}{{dt}} = \frac{{\partial {A^i}}}{{\partial {a^l}}}{J^\textit{lk}}\frac{{\partial {A^j}}}{{\partial {a^k}}}\frac{{\partial H}}{{\partial {A^j}}} = {J^\textit{ij}}\frac{{\partial H}}{{\partial {A^j}}},$$
where the indices of $i$, $j$, $k$, and $l$ indicate components of the vectors or the matrix ${\boldsymbol J}$, and Einstein’s summation convention is used. From Eq. (10), the following equation can be derived as
$$\frac{{\partial {A^i}}}{{\partial {a^l}}}{J^\textit{lk}}\frac{{\partial {A^j}}}{{\partial {a^k}}} = {J^\textit{ij}}.$$

Equation (11) is called the symplectic conditions that ensure the phase space volume conservation, and they can be written using Eqs. (5) and (8) as

$${\left\{ {{P^i},{P^j}} \right\}_{{\boldsymbol q},{\boldsymbol p}}} = 0,$$
$${\left\{ {{Q^i},{Q^j}} \right\}_{{\boldsymbol q},{\boldsymbol p}}} = 0,$$
$${\left\{ {{Q^i},{P^j}} \right\}_{{\boldsymbol q},{\boldsymbol p}}} = {\delta ^\textit{ij}},$$
where {} denotes the Poisson bracket, $i$ and $j$ denote components of each vector, and $\delta^{ij}$ indicates the Kronecker’s delta.

The ${\boldsymbol P}$ and ${\boldsymbol Q}$ of the phase point ${\boldsymbol A}$ can be written with the ${\boldsymbol p}$ and ${\boldsymbol q}$ of the phase point ${\boldsymbol a}$ in a finite-difference explicit form using Eqs. (3) and (4) as

$${P^i} = {p^i} - \frac{{\partial H}}{{\partial {q^i}}}\left( {{\boldsymbol q},{\boldsymbol p}} \right)\Delta t,$$
$${Q^i} = {q^i} + \frac{{\partial H}}{{\partial {p^i}}}\left( {{\boldsymbol q},{\boldsymbol p}} \right)\Delta t.$$
It can be easily shown that Eqs. (15) and (16) satisfy two of the symplectic conditions, i.e., Eqs. (12) and (13). However, the rest of the symplectic conditions, i.e., Eq. (14) become
$$\begin{split} & {\left\{ {{Q^i},{P^j}} \right\}_{{\boldsymbol q},{\boldsymbol p}}}=\\& {\delta ^\textit{ij}} +\sum\limits_l {\left[ {\frac{{{\partial ^2}H}}{{\partial {p^i}\partial {p^l}}}\frac{{{\partial ^2}H}}{{\partial {q^l}\partial {q^j}}} - \frac{{{\partial ^2}H}}{{\partial{p^i}\partial {q^l}}}\frac{{{\partial ^2}H}}{{\partial {p^l}\partial {q^j}}}} \right]} \Delta {t^2}.\end{split}$$
On the right-hand side of Eq. (17), the first order of $\Delta t$ vanishes. The second term on the right-hand side of Eq. (17) is the second order of $\Delta t$ and will not vanish in general when $\Delta t$ is finite, which means that the Hamiltonian equations of Eq. (9) at the phase space point ${\boldsymbol A}$ slightly deviate from its exact form. This causes cumulative error of a light-ray path calculated numerically at subsequent time steps. The cumulative error should be improved by reducing the deviation of the Hamiltonian equation from the exact form.

On the other hand, instead of the finite-difference explicit form of Eqs. (15) and (16), the following equations in the form of the first-order symplectic integrator (i.e., semi-implicit Euler form) are therefore adopted as

$${P^i} = {p^i} - \frac{{\partial H\left( {{\boldsymbol q},{\boldsymbol p}} \right)}}{{\partial {q^i}}}\Delta t,$$
$${Q^i} = {q^i} + \frac{{\partial H\left( {{\boldsymbol q},{\boldsymbol P}} \right)}}{{\partial {P^i}}}\Delta t.$$
Note that the second term on the right-hand side of Eq. (19) is defined at the phase space of (${\boldsymbol q}$, ${\boldsymbol P}$). It can be easily shown that the two symplectic conditions represented by Eqs. (12) and (13) can be satisfied by Eqs. (18) and (19). The rest of the equations of the symplectic condition [i.e., Eq. (14)] can satisfy the following equation as
$${\left\{ {{Q^i},{P^j}} \right\}_{{\boldsymbol q},{\boldsymbol p}}} = {\delta ^\textit{ij}} + O\left( {\Delta {t^3}} \right),$$
when the following condition holds:
$$\frac{{\partial H\left( {{\boldsymbol q},{\boldsymbol p}} \right)}}{{\partial {\boldsymbol q}\partial {\boldsymbol p}}} = 0.$$
On the right-hand side of Eq. (20), the first and second order of $\Delta t$ vanishes, which means that Eq. (20) is more accurate than Eq. (17). Thus, the Hamiltonian equation is ensured to be almost form invariant under the transformation represented by Eqs. (18) and (19). This should reduce the accumulation of errors of numerical calculations along time steps. Analogously, the fourth-order symplectic integrator can also be constructed [18,23].

As one of the Hamiltonians that satisfy the condition of Eq. (21) for an optical system, the following optical Hamiltonian, ${H_{\rm opt}}$, that consists of two independent terms, i.e., one term dependent on the position, ${\boldsymbol q}$, and the other term dependent on the momentum, ${\boldsymbol p}$, is adopted as

$${H_{{\rm opt}}} = \frac{{{{\boldsymbol p}^2}}}{2} - \frac{{{n^2}\left( {\boldsymbol q} \right)}}{2}.$$

Using the optical Hamiltonian, the ray equation of Eq. (1) can be transformed with the light-ray momentum represented by Eq. (2) as

$$\frac{{d{\boldsymbol p}}}{{{{ds} / n}}} = - \frac{{\partial {H_{{\rm opt}}}}}{{\partial {\boldsymbol q}}}.$$
This equation can be further transformed as
$$\frac{{d{\boldsymbol p}}}{{dt}} = - \frac{{\partial {H_{{\rm opt}}}}}{{\partial {\boldsymbol q}}},$$
where $dt$ is defined as
$$dt = \frac{{ds}}{n}.$$

In this paper, $t$ is called time, which is not actually physical time.

The derivative of the position ${\boldsymbol q}$ with respect to $t$ can be derived using the optical Hamiltonian ${H_{\rm opt}}$ as

$$\frac{{d{\boldsymbol q}}}{{dt}} = \frac{{\partial {H_{{\rm opt}}}}}{{\partial {\boldsymbol p}}}.$$

Equations (24) and (26) are in the form of the Hamilton’s equations with canonical coordinates of (${\boldsymbol q}$, ${\boldsymbol p}$). The initial condition of the optical Hamiltonian can be derived from Eq. (2) as

$${H_{{\rm opt}}} = \frac{{{{\boldsymbol p}^2}}}{2} - \frac{{{n^2}}}{2} = 0.$$

Using the optical Hamiltonian of Eq. (22), the first-order symplectic integrator of Eqs. (18) and (19) can be written as

$${\boldsymbol P} = {\boldsymbol p} + {\left. {\nabla \left( {\frac{{{n^2}}}{2}} \right)} \right|_{\boldsymbol q}}\Delta t,$$
$${\boldsymbol Q} = {\boldsymbol q} + {\boldsymbol P}\Delta t.$$

Note that the second term on the right-hand side of Eq. (29) is set to ${\boldsymbol P}({\boldsymbol q},\;{\boldsymbol p})$ instead of ${\boldsymbol p}$. Equations (28) and (29) represent a symplectic ray tracing on the basis of Hamiltonian optics, which ensures that the Hamiltonian equations corresponding to Eqs. (24) and (26) are almost form invariant. The canonical transformation that satisfies the symplectic conditions ensures that the phase space volume is conserved. In geometrical optics, the four-dimensional phase space volume conservation is called etendue conservation, which is widely used for lighting lens designs [24,25].

3. NUMERICAL CALCULATIONS OF LIGHT-RAY PATHS USING SYMPLECTIC RAY TRACING IN GRIN LENSES

A. GRIN Lens with Cylindrical Symmetry

GRIN lenses having cylindrical symmetry and flat surfaces are widely used for several applications, such as imaging optics, photocopiers, and scanners. A typical refractive index of a cylindrical GRIN lens [26] can be written in a cylindrical coordinate system ($\rho $, $\theta $, $z$) as

$$n = {n_0}\sqrt {1 - {{\left( {A\rho } \right)}^2}} ,$$
where $\rho $ is a radius, ${n_0}$ is a refractive index of center core of the GRIN lens, and $A$ is a positive constant. In the Cartesian coordinate system ($x$, $y$, $z$), the radius $\rho $ can be written as
$$\rho = \sqrt {{x^2} + {y^2}} .$$
The analytical solution of the ray equation [i.e., Eq. (1)] for the refractive index of Eq. (30) leads to the following trajectory of a light ray initially parallel to the $z$ axis at the position of $x = {x_0}$ and $z = {0}$ in a cross-sectional $xz$-plane at $y = {0}$:
$$\frac{x}{{{x_0}}} = \sin \left( {\frac{{2\pi }}{\Lambda }z + \frac{\pi }{2}} \right),$$
where $\Lambda$ denotes a period and can be written as
$$\Lambda = \frac{{2\pi }}{A}\sqrt {1 - {{\left( {A{x_0}} \right)}^2}} .$$

In this work, the values of the parameters of $A$, ${n_0}$, and ${x_0}$ are set to ${0.5}\;{{\rm mm}^{ - 1}}$, 1.564, and 0.5 mm, respectively. These values of the parameters are taken from Ref. [26]. The analytical value of the period $\Lambda $ is 12.16733603 mm.

Figure 1 shows the light-ray paths numerically calculated using the several methods in the cross-sectional $xz$-plane. The horizontal axis indicates the $z$ axis normalized by the period $\Lambda $, and the vertical axis indicates the $x$ axis normalized by the initial $x$ coordinate of ${x_0}$. The light rays are incident on the cylindrical GRIN lens in a direction parallel to the $z$ axis. An ideal light-ray path arrives at $(x,z) = ({x_0}, \Lambda) $. Numerical calculations using the first-order explicit Euler method represented by Eqs. (15) and (16), the first-order symplectic method represented by Eqs. (28) and (29), and the fourth-order symplectic method [18,23] are, respectively, plotted by a dotted line with filled circle, a solid line with open circle, and a solid line with square. The analytical solution is also plotted by a dashed line. The time step $\Delta t$ for each method is set to 0.5 mm (i.e., $\Delta t/\Lambda = {0.04109363}$). Note that the time is not actual time and has the dimension of length in this work. It can be found from the figure that the first-order symplectic method drastically reduces cumulative error in comparison with the first-order explicit Euler method and approximates the analytical solution. The fourth-order symplectic method is found to be more accurate than the first-order symplectic method, as expected.

 figure: Fig. 1.

Fig. 1. Light-ray paths numerically calculated using several methods in a cross-sectional $xz$-plane at $y = {0}$ of a cylindrical GRIN lens. The horizontal axis indicates the $z$ axis normalized by a period $\Lambda $, and the vertical axis indicates the $x$ axis normalized by an initial $x$ coordinate of ${x_0}$. The light rays are incident on a cylindrical GRIN lens in a direction parallel to the $z$ axis. An ideal light-ray path arrives at $(x, z) = (x_0, \Lambda)$. Numerical calculations using the first-order explicit Euler method, the first-order symplectic method, and the fourth-order symplectic method are plotted by a dotted line with filled circle, a solid line with open circle, and a solid line with square, respectively. The analytical solution is also plotted by dashed line.

Download Full Size | PDF

B. GRIN Lens with Spherical Symmetry

One of the widely known GRIN lenses having spherical symmetry is the Luneburg lens [27]. The refractive index of the Luneburg lens having a radius of ${R_0}$ can be written as

$$n = \sqrt {2 - {{\left( {\frac{r}{{{R_0}}}} \right)}^2}} ,$$
where $r$ denotes a radius from the center of the Luneburg lens. A light-ray path can be solved analytically in a cross-sectional $xz$-plane at $y = {0}$ with an initial $x$ position of ${x_0}$ and an initial direction parallel to the $z$ axis as
$${\left( {\frac{x}{{{x_0}}}} \right)^2} + {\left( {\frac{z}{{{R_0}}} + \sqrt {1 - {{\left( {\frac{{{x_0}}}{{{R_0}}}} \right)}^2}} \frac{x}{{{x_0}}}} \right)^2} = 1,$$
where the center of the Luneburg lens is placed at the coordinate origin. An ideal focus point can be derived from Eq. (35) as $(x,z) = ({0},{R_0})$.

Figure 2 shows light-ray paths numerically calculated using several methods in the cross-sectional $xz$-plane at $y = {0}$. The horizontal axis indicates the $z$ axis normalized by the lens radius ${R_0}$, and the vertical axis indicates the $x$ axis normalized by the radius ${R_0}$. The light rays are incident on the Luneburg lens in a direction parallel to the $z$ axis. An ideal focus point is $(x,z) = ({0},{R_0})$. Numerical calculations using the first-order explicit Euler method, the first-order symplectic method, and the fourth-order symplectic method are plotted by a dotted line with filled circle, a solid line with open circle, and a solid line with square, respectively.

 figure: Fig. 2.

Fig. 2. Light-ray paths numerically calculated using several methods in a cross-sectional $xz$-plane at $y = {0}$ of the Luneburg lens. The horizontal axis indicates the $z$ axis normalized by the lens radius ${R_0}$, and the vertical axis indicates the $x$ axis normalized by the lens radius ${R_0}$. The light rays are incident on the Luneburg lens in the direction parallel to the $z$ axis with initial $x$ positions, ${x_0}$, of 0.15 and 0.75. An ideal focus point is $(x, z) = (0, R_0)$. Numerical calculations using the first-order explicit Euler method, the first-order symplectic method, and the fourth-order symplectic method are plotted by a dotted line with filled circle, a solid line with open circle, and a solid line with square, respectively. The analytical calculation is also plotted by dashed line.

Download Full Size | PDF

The analytical solution is also plotted by dashed line. The time step $\Delta t$ for each method is set to 0.01 ${R_0}$. The $x$ positions of ${x_0}$ of incident light rays (i.e., $x$ positions of light rays at the periphery of the Luneburg lens) are set to 0.15 ${R_0}$ and 0.75 ${R_0}$, respectively. It can be found from the figure that the first-order symplectic method drastically reduces cumulative error in comparison with the first-order explicit Euler method and approximates the analytical solution. The first-order symplectic method also approximates the fourth-order symplectic method.

4. RESULTS

Errors of numerically calculated light-ray paths using several methods for the cylindrical and spherical GRIN lenses are calculated and are compared with each other.

First, errors of light-ray paths for the cylindrical GRIN lens as mentioned in Section 3.A are calculated. The error is defined as a ratio of an absolute value of deviation of $x$ from analytical solution ${x_{\rm ana}}$ to ${x_{\rm ana}}$ at $z = \Lambda $. Figure 3 shows errors with respect to the time step $\Delta t$ normalized by $\Lambda $. The horizontal axis indicates the time step, and the vertical axis indicates the error. Numerical calculations using the first-order explicit Euler method, the first-order symplectic method, and the fourth-order symplectic method are plotted by a dotted line, circle mark, and square mark, respectively. Numerical calculation using the fourth-order Runge–Kutta method [14] is also plotted by dashed line. The time steps, $\Delta t$, for both the symplectic method and Runge–Kutta method are set to the same value where the definitions of the time in both methods are the same. It can be found from the figure that the first-order symplectic method drastically reduces cumulative error in comparison with the first-order Euler method. The first-order symplectic method also approximates the fourth-order Runge–Kutta method. The fourth-order symplectic method agrees well the fourth-order Runge–Kutta method, whereas there is a slight difference in the larger $\Delta t$ of more than 0.09 $\Lambda $.

 figure: Fig. 3.

Fig. 3. Errors with respect to time step $\Delta t$ normalized by $\Lambda $ for the GRIN lens having cylindrical symmetry as mentioned in Section 3.A. The error is defined as a ratio of an absolute value of deviation of $x$ from the analytical solution ${x_{\rm ana}}$ to ${x_{\rm ana}}$ at $z = \Lambda $. The horizontal axis indicates the time step, and the vertical axis indicates the error. Numerical calculations using the first-order explicit Euler method, the first-order symplectic method, and the fourth-order symplectic method are plotted by a dotted line, circle mark, and square mark, respectively. Numerical calculation using the fourth-order Runge–Kutta method is also plotted by dashed line.

Download Full Size | PDF

Second, errors of light-ray paths for the spherical GRIN lens (i.e., the Luneburg lens) as mentioned in Section 3.B are calculated. The error is defined as a ratio of an absolute value of deviation of $x$ from the analytical solution to the lens radius ${R_0}$ at the focal point of $z = {R_0}$. Figure 4 shows errors with respect to the time step $\Delta t$ normalized by ${R_0}$. The horizontal axis indicates the time step, and the vertical axis indicates the error. Numerical calculations using the first-order explicit Euler method, the first-order symplectic method, and the fourth-order symplectic method are plotted by a dotted line, circle mark, and square mark, respectively. Numerical calculation using the fourth-order Runge–Kutta method is also plotted by dashed line. It can be found from the figure that the first-order symplectic method drastically reduces cumulative error in comparison with the first-order Euler method. The first-order symplectic method also approximates the fourth-order symplectic method. The first-order symplectic method, the fourth-order symplectic method, and the fourth-order Runge–Kutta method are the same order of magnitude. The error of the first-order symplectic method can be less than that of the fourth-order Runge–Kutta method.

 figure: Fig. 4.

Fig. 4. Errors with respect to the time step $\Delta t$ normalized by ${R_0}$ for the GRIN lens having spherical symmetry (i.e., the Luneburg lens) as mentioned in Section 3.B. The error is defined as a ratio of an absolute value of deviation of $x$ from the analytical solution to the lens radius ${R_0}$ at the focal point of $z = {R_0}$. The horizontal axis indicates the time step, and the vertical axis indicates the error. Numerical calculations using the first-order explicit Euler method, the first-order symplectic method, and the fourth-order symplectic method are plotted by a dotted line, circle mark, and square mark, respectively. Numerical calculation using the fourth-order Runge–Kutta method is also plotted by dashed line.

Download Full Size | PDF

The first-order symplectic method represented by Eqs. (28) and (29) is simple enough to reduce the computational cost and approximates higher-order methods. Thus, it should be a promising candidate for ray tracing in GRIN media.

5. CONCLUSIONS

A method for ray tracing through GRIN media with a simple formulation and high accuracy is proposed on the basis of Hamiltonian optics. The optical Hamiltonian can be constructed of two independent terms, i.e., one term dependent on position and the other term dependent on momentum. The symplectic integrator is applicable to the separable Hamiltonian system. It can be shown that the Hamiltonian equations are ensured to be almost form invariant at each time step of numerical calculations using the symplectic integrator, which can increase accuracy of the ray tracing. The light-ray paths in GRIN lenses are calculated using the symplectic-ray-tracing method of the first and fourth order, which approximate those calculated on the basis of the fourth-order Runge–Kutta algorithm. It is shown that the deviations of light-ray paths from the ideal paths are supremely reduced by the symplectic-ray-tracing method in comparison with the explicit Euler method. The first-order symplectic ray tracing represented by Eqs. (28) and (29) is simple enough to reduce the computational cost and approximates higher-order methods. Thus, symplectic ray tracing should be a promising candidate for ray tracing in GRIN media and should enlarge the range of the method choice.

Disclosures

The author declares no conflicts of interest.

REFERENCES

1. J. B. Pendry, D. Schurig, and D. R. Smith, “Controlling electromagnetic field,” Science 312, 1780–1782 (2006). [CrossRef]  

2. U. Leonhardt, “Optical conformal mapping,” Science 312, 1777–1780 (2006). [CrossRef]  

3. Y. A. Urzhumov, N. B. Kundtz, D. R. Smith, and J. B. Pendry, “Cross-section comparisons of cloaks designed by transformation optical and optical conformal mapping approaches,” J. Opt. 13, 024002 (2011). [CrossRef]  

4. U. Leonhardt, “Notes on conformal invisibility devices,” New J. Phys. 8, 118 (2006). [CrossRef]  

5. C. Lu and Z. L. Mei, “Multi-functional lens based on conformal mapping,” Opt. Express 23, 19901–19910 (2015). [CrossRef]  

6. H. Ohno and T. Usui, “Gradient-index dark hole based on conformal mapping with etendue conservation,” Opt. Express 27, 18493–18507 (2019). [CrossRef]  

7. H. Ohno and K. Toya, “Reconstruction method of axisymmetric refractive index fields with background-oriented schlieren,” Appl. Opt. 57, 9062–9069 (2018). [CrossRef]  

8. H. Ohno and K. Toya, “Scalar potential reconstruction method of axisymmetric 3D refractive index fields with background-oriented schlieren,” Opt. Express 27, 5990–6002 (2019). [CrossRef]  

9. E. W. Marchand, “Fifth-order analysis of GRIN lenses,” Appl. Opt. 24, 4371–4374 (1985). [CrossRef]  

10. D. T. Moore, “Ray tracing in gradient-index media,” J. Opt. Soc. Am. 65, 451–455 (1975). [CrossRef]  

11. M. Sluijter, D. K. G. de Boer, and H. P. Urbach, “Ray-optics analysis of inhomogeneous biaxially anisotropic media,” J. Opt. Soc. Am. A 26, 317–329 (2009). [CrossRef]  

12. Y. Nishidate, “Closed-form analytical solutions for ray tracing in optically anisotropic inhomogeneous media,” J. Opt. Soc. Am. A 30, 1373–1379 (2013). [CrossRef]  

13. J. E. Gomez-Correa, V. Coello, A. Garaza-Rivera, N. P. Puente, and S. Chavez-Cerda, “Three-dimensional ray tracing in spherical and elliptical generalized Luneburg lenses for application in the human eye lens,” Appl. Opt. 55, 2002–2010 (2016). [CrossRef]  

14. A. Sharma, D. V. Kumar, and A. K. Ghatak, “Tracing rays through graded-index media: A new method,” Appl. Opt. 21, 984–987 (1982). [CrossRef]  

15. G. W. Forbes, “On variational problems in parametric form,” Am. J. Phys. 59, 1130–1140 (1991). [CrossRef]  

16. E. Hairer, C. Lubich, and G. Wanner, “Geometric numerical integration illustrated by the Störmer/Verlet method,” Acta Numer. 12, 399–450 (2003). [CrossRef]  

17. R. D. Ruth, “A canonical integration technique,” IEEE Trans. Nucl. Sci. 30, 2669–2671 (1983). [CrossRef]  

18. E. Forest and R. D. Ruth, “Fourth-order symplectic integration,” Physica D 43, 105–117 (1990). [CrossRef]  

19. H. Goto, K. Tatsumura, and A. R. Dixon, “Combinatorial optimization by simulating adiabatic bifurcations in nonlinear Hamiltonian systems,” Sci. Adv. 5, eaav2372 (2019). [CrossRef]  

20. T. Hirono, W. W. Lui, K. Yokoyama, and S. Seki, “Stability and numerical dispersion of symplectic fourth-order time-domain schemes for optical field simulation,” J. Lightwave Technol. 16, 1915–1920 (1998). [CrossRef]  

21. T. Hirono and Y. Yoshikuni, “Accurate modeling of dielectric interfaces by the effective permittivities for the fourth-order symplectic finite-difference time-domain method,” Appl. Opt. 46, 1514–1524 (2007). [CrossRef]  

22. H. Fan and H. Lu, “Symplectic wavelet transformation,” Opt. Lett. 31, 3432–3434 (2006). [CrossRef]  

23. H. Yoshida, “Construction of higher order symplectic integrators,” Phys. Lett. A 150, 262–268 (1990). [CrossRef]  

24. R. Winston, J. Minano, and P. Benitez, Nonimaging Optics (Academic, 2005).

25. H. Ohno, “Design of a coaxial light guide producing a wide-angle light distribution,” Appl. Opt. 56, 3977–3983 (2017). [CrossRef]  

26. B. Richerzhagen, “Finite element ray tracing: A new method for ray tracing in gradient-index media,” Appl. Opt. 35, 6186–6189 (1996). [CrossRef]  

27. R. K. Luneburg and Brown University, Mathematical Theory of Optics (Providence, 1944).

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (4)
Fig. 1.
Fig. 1. Light-ray paths numerically calculated using several methods in a cross-sectional $xz$ -plane at $y = {0}$ of a cylindrical GRIN lens. The horizontal axis indicates the $z$ axis normalized by a period $\Lambda $ , and the vertical axis indicates the $x$ axis normalized by an initial $x$ coordinate of ${x_0}$ . The light rays are incident on a cylindrical GRIN lens in a direction parallel to the $z$ axis. An ideal light-ray path arrives at $(x, z) = (x_0, \Lambda)$ . Numerical calculations using the first-order explicit Euler method, the first-order symplectic method, and the fourth-order symplectic method are plotted by a dotted line with filled circle, a solid line with open circle, and a solid line with square, respectively. The analytical solution is also plotted by dashed line.

View in Article | Download Full Size | PDF

Fig. 2.
Fig. 2. Light-ray paths numerically calculated using several methods in a cross-sectional $xz$ -plane at $y = {0}$ of the Luneburg lens. The horizontal axis indicates the $z$ axis normalized by the lens radius ${R_0}$ , and the vertical axis indicates the $x$ axis normalized by the lens radius ${R_0}$ . The light rays are incident on the Luneburg lens in the direction parallel to the $z$ axis with initial $x$ positions, ${x_0}$ , of 0.15 and 0.75. An ideal focus point is $(x, z) = (0, R_0)$ . Numerical calculations using the first-order explicit Euler method, the first-order symplectic method, and the fourth-order symplectic method are plotted by a dotted line with filled circle, a solid line with open circle, and a solid line with square, respectively. The analytical calculation is also plotted by dashed line.

View in Article | Download Full Size | PDF

Fig. 3.
Fig. 3. Errors with respect to time step $\Delta t$ normalized by $\Lambda $ for the GRIN lens having cylindrical symmetry as mentioned in Section 3.A. The error is defined as a ratio of an absolute value of deviation of $x$ from the analytical solution ${x_{\rm ana}}$ to ${x_{\rm ana}}$ at $z = \Lambda $ . The horizontal axis indicates the time step, and the vertical axis indicates the error. Numerical calculations using the first-order explicit Euler method, the first-order symplectic method, and the fourth-order symplectic method are plotted by a dotted line, circle mark, and square mark, respectively. Numerical calculation using the fourth-order Runge–Kutta method is also plotted by dashed line.

View in Article | Download Full Size | PDF

Fig. 4.
Fig. 4. Errors with respect to the time step $\Delta t$ normalized by ${R_0}$ for the GRIN lens having spherical symmetry (i.e., the Luneburg lens) as mentioned in Section 3.B. The error is defined as a ratio of an absolute value of deviation of $x$ from the analytical solution to the lens radius ${R_0}$ at the focal point of $z = {R_0}$ . The horizontal axis indicates the time step, and the vertical axis indicates the error. Numerical calculations using the first-order explicit Euler method, the first-order symplectic method, and the fourth-order symplectic method are plotted by a dotted line, circle mark, and square mark, respectively. Numerical calculation using the fourth-order Runge–Kutta method is also plotted by dashed line.

View in Article | Download Full Size | PDF

Equations (35)

Equations on this page are rendered with MathJax. Learn more.

d d s [ n ( q ) d q d s ] = ∇ n ( q ) ,
p = n ( q ) d q d s .
d p d t = − ∂ H ∂ q ,
d q d t = ∂ H ∂ p ,
a ( t ) = [ q ( t ) p ( t ) ] .
d a d t = J ∂ H ∂ a ,
J = [ 0 1 − 1 0 ] .
A = [ Q P ] = a ( t + Δ t ) = [ q ( t + Δ t ) p ( t + Δ t ) ] .
d A d t = J ∂ H ∂ A .
d A i d t = ∂ A i ∂ a l J lk ∂ A j ∂ a k ∂ H ∂ A j = J ij ∂ H ∂ A j ,
∂ A i ∂ a l J lk ∂ A j ∂ a k = J ij .
{ P i , P j } q , p = 0 ,
{ Q i , Q j } q , p = 0 ,
{ Q i , P j } q , p = δ ij ,
P i = p i − ∂ H ∂ q i ( q , p ) Δ t ,
Q i = q i + ∂ H ∂ p i ( q , p ) Δ t .
{ Q i , P j } q , p = δ ij + ∑ l [ ∂ 2 H ∂ p i ∂ p l ∂ 2 H ∂ q l ∂ q j − ∂ 2 H ∂ p i ∂ q l ∂ 2 H ∂ p l ∂ q j ] Δ t 2 .
P i = p i − ∂ H ( q , p ) ∂ q i Δ t ,
Q i = q i + ∂ H ( q , P ) ∂ P i Δ t .
{ Q i , P j } q , p = δ ij + O ( Δ t 3 ) ,
∂ H ( q , p ) ∂ q ∂ p = 0.
H o p t = p 2 2 − n 2 ( q ) 2 .
d p d s / n = − ∂ H o p t ∂ q .
d p d t = − ∂ H o p t ∂ q ,
d t = d s n .
d q d t = ∂ H o p t ∂ p .
H o p t = p 2 2 − n 2 2 = 0.
P = p + ∇ ( n 2 2 ) | q Δ t ,
Q = q + P Δ t .
n = n 0 1 − ( A ρ ) 2 ,
ρ = x 2 + y 2 .
x x 0 = sin ⁡ ( 2 π Λ z + π 2 ) ,
Λ = 2 π A 1 − ( A x 0 ) 2 .
n = 2 − ( r R 0 ) 2 ,
( x x 0 ) 2 + ( z R 0 + 1 − ( x 0 R 0 ) 2 x x 0 ) 2 = 1 ,
Select as filters
Select Topics Cancel
Top
       
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.
Privacy | Terms of Use