## Abstract

Ray tracing is an important technique for predicting optical system performance. In the field of transformation optics, the Hamiltonian equations of motion for ray tracing are well known. The numerical solutions to the Hamiltonian equations of motion are affected by the complexities of the inhomogeneous and anisotropic indices of the optical device. Based on our knowledge, no previous work has been conducted on ray tracing for transformation optics with extreme inhomogeneity and anisotropicity. In this study, we present the use of 3D reverse ray tracing in transformation optics. The reverse ray tracing is derived from Fermat’s principle based on a sweeping method instead of finding the full solution to ordinary differential equations. The sweeping method is employed to obtain the eikonal function. The wave vectors are then obtained from the gradient of that eikonal function map in the transformed space to acquire the illuminance. Because only the rays in the points of interest have to be traced, the reverse ray tracing provides an efficient approach to investigate the illuminance of a system. This approach is useful in any form of transformation optics where the material property tensor is a symmetric positive definite matrix. The performance and analysis of three transformation optics with inhomogeneous and anisotropic indices are explored. The ray trajectories and illuminances in these demonstration cases are successfully solved by the proposed reverse ray tracing method.

© 2015 Optical Society of America

## 1. Introduction

Anisotropic transformation optics can manipulate electromagnetic fields in many interesting ways. Ample studies have been conducted from the infrared spectrum to the optical spectrum such as that on an invisibility cloak [1–5] and on ray optics to enhance an imaging system [6]. Ray tracing was employed in these studies to predict the optical performance of transformation optics. In optical design, predicting optical performance is important. Ray tracing calculations are performed to describe the geometric aspects of imaging by propagating rays through an optical system [4,7]. These ray properties can be described by dispersion relations. A well-known ray tracing method is based on the Hamiltonian equations of motion. The dispersion relations are in Hamiltonian form and expressed in nonlinear ordinary differential equations (ODEs). Hamiltonian ray tracing has already been applied to traditional biaxial material [8] and transformation optics [1,2,4,6,9]. Another approach to ray tracing is by using the Newtonian equations of motion, which has been employed to analyze the performance of cloak grating [3,5]. All of these approaches are based on a set of ODEs.

Three reasons are cited for developing another ray tracing method that is not dependent on solving ODEs. First, because of the stiff or ill-conditioned [10,11] nature of the Hamiltonian equations of motion arising from the inhomogeneity and anisotropicity of the material, the uses of Hamiltonian ray tracing are limited [12]. Standard solvers such as the Euler method or the explicit Runge–Kutta methods may exhibit instability or other issues when solving stiff ODEs [10,11]. The computational effort depends on the complexity of the equations. Second, if the analytical expression for the material property tensor is unknown, the Hamiltonian equations of motion cannot be derived. Third, Hamiltonian ray tracing handling transformation optics with extreme inhomogeneity and anisotropicity have not been conducted yet. To conduct ray tracing in any possible transformed device, a method that does not require a full numerical solution to the Hamiltonian equations of motion is necessary.

In this study, we apply 3D reverse ray tracing, whose main principle is the Fermat principle [13], in transformation optics for the first time. An eikonal function of the investigated domain in transformation optics can be obtained with the application of the sweeping method [14,15]. To obtain the illuminance of a device, we further treated the wave vectors of the rays as the gradient of that eikonal function map. Illuminance is an important indicator for designing nonimaging optical systems [16]. As the rays are traced back from the target to the source, reverse ray tracing provides an efficient approach to investigate the illuminance of a region of interest compared with traditional ray tracing methods [17, 18]. Only the rays in the points of interest need to be traced. Moreover, in transformation optics, the reverse ray tracing can be performed even if the analytical expression of the material property tensor of the investigated system is unknown.

The rest of this paper is organized as follows. In Section 2, the reverse ray tracing algorithm is described. Three examples are presented in Sections 3 and 4 to demonstrate the accuracy and utility of our proposed method. A comparison of ray trajectories and illuminance calculated by reverse ray tracing and Hamiltonian ray tracing is shown in Section 3. In Section 4, two examples, which cannot be easily solved using the Hamiltonian equations of motion, are traced successfully by our proposed reverse ray tracing.

## 2. Solver and geodesic

The ray properties are described by the Hamilton–Jacobi equation. The dispersion equation in the transformed space is a Hamilton–Jacobi equation, which can be expressed as

*T*(X) is the eikonal function [8]. The surfaces of the constant eikonal function can be treated as the wavefronts. The matrix

**M**(X) is the material property tensor of simulated optical device at location X(

*x*,

*y*,

*z*).

**M**(X) is the required symmetric positive definite 3 × 3 matrix for a given material and its determinant is denoted by det(

**M**(X)). The wave vector

**K**, which is the gradient of the eikonal function, can be expressed as ∇

*T*(X).

Table 1 illustrates the algorithm of our reverse ray tracing, including solving the eikonal function, ray trajectories, and illuminance. We modified the sweeping method to obtain the eikonal function in 3D (Step 2 in Table 1). Figure 1 shows an illustration of triangulation in reverse ray tracing. Rectangular grids were adopted in this study because of a faster convergence speed. The device is divided into rectangular grids, as shown in Fig. 1(a), and the grid of the *yz*-plane is shown in Fig. 1(b). Originally, the eikonal function of each node depends on 24 neighboring triangles (8Δ∈xy plane, 8Δ∈xz plane, and 8Δ∈zy plane). In this study, eight tetrahedrons [19] (such as CNLW) are introduced to find the minimum root of the eikonal function before sweeping all 24 triangles. If the eikonal function of each node is obtained by eight tetrahedrons, sweeping 24 triangles is not necessary. When considering the eikonal value at X* _{C}* node, the real roots,

*T*(X

*), of Eq. (1) are solved in each tetrahedron and each triangle around the X*

_{C}*node. The wave vector, ∇*

_{C}*T*(X

*), of the tetrahedron CNLW can be expressed as*

_{C}*x*,

*y*,

*z*) is in a Cartesian coordinate,

**X**

*is the vector from node X*

_{NC}*to node X*

_{N}*, and $\left|{X}_{NC}\right|$ is the distance between node X*

_{C}*and node X*

_{C}*. Each value of*

_{N}*T*corresponds to one characteristic direction,

**d**. The characteristic direction, which is indicated by a black arrow in Fig. 1(a), depends on the local property tensor and can be expressed as

Our goal is to choose the minimum of the roots that satisfies the causality condition, wherein the characteristic direction, **d**, should pass through X* _{C}* within the tetrahedron for 3D case and within the triangle for 2D case. The whole algorithm for solving eikonal function is summarized as follows: first, at each X

*node, the minimum of real roots of Eq. (1) from eight neighboring tetrahedrons is determined. If no root satisfies the causality condition, the same process is done for 24 neighboring triangles. If no root satisfies the causality condition in certain triangles, the group velocity, |∇*

_{C}

_{K}*J*|, is introduced to determine the eikonal value. The group velocity is proportional to

**M**(X). The eikonal value of triangle ΔCSA is given by the following equation [14]:

**v**

_{g}

^{SC}and

**v**

_{g}

^{AC}denote the group velocities along edges SC and AC. In the cases of transformation optics, the group velocity is $\left|{\nabla}_{K}J\right|$and is proportional to

**M**(X). Therefore,$\left|{\text{X}}_{AC}\right|/{\text{v}}_{g}^{AC}$, the travel time between nodes A and C, for example, can be modified into$\sqrt{{\text{X}}_{AC}M{({\text{X}}_{C})}^{-1}{\text{X}}_{AC}^{T}}$. This form is useful for calculating wave propagation between two nodes in 3D.

Finally, minimum of *T*(X* _{C}*) depending on eight tetrahedrons and 24 triangles is determined. The convergence criterion for all fast sweeping iterations is when the maximum difference between two consecutive iterations is less than 10

^{−9}. The final eikonal value,

*T*(X

*), of each node should satisfy Eq. (1). Most of*

_{i}*T*(X

*) are determined in Step 3 in practice. After all eikonal values in the ray tracing domains are obtained, the ray trajectories based on back propagation [20] can be calculated through the following expression:*

_{i}*γ*is the discrete pathway (geodesic) for

_{k,}*k*step and

*τ*is the step size in back propagation. As shown in Fig. 1(b), the starting location of back propagation is marked with a green dot. The location of the source marked by the black dot is the initial point of the wavefront and also the reference point. Initialization assigns infinite values to eikonal function,

*T*, for all nodes except the initial point. The

*T*of the initial point is zero. Sweeping starts from the node that is nearest to the initial point; for example, the sweeping direction is from the blue to the red nodes, as shown in Fig. 1(b).

The idea of obtaining wave vectors from the gradient of the eikonal function in an investigated domain allows us to find the illuminance of a device. In the sweeping method, upwind singularity at the source exists in a point-source problem [14]. However, the illuminance in the observation plane is usually far from the source. Therefore, the calculation accuracy of the illuminance is not affected by the upwind singularity. The illuminance from a device can be obtained through the following steps: first, the wave vector is viewed as a gradient of the eikonal function, ∇*T*(X* _{i}*). The wave vector at the point X

*is*

_{i}**K**= [∇

*(X*

_{x}T*), ∇*

_{i}*(X*

_{y}T*), ∇*

_{i}*(X*

_{z}T*)]. The electric field vector*

_{i}**E**is orthogonal to the Poynting vector and is defined by the following equation:

*J*is the Hamilton–Jacobi equation expressed in Eq. (1). The corresponding magnetic field vector

**H**in the transformed space is defined as

Moreover, for calculating the transmission coefficient at the interface between transformation optics and air, the wave vectors of rays before and after refraction should be determined first and can be expressed as follows [21]:

**K**

_{1}is the incident wave vector in the transformed space;

**K**

_{2}is the wave vector after refraction in free space; and

**n**is the unit normal to the interface at the point X

*. The transmitted electromagnetic field and Fresnel equations are based on the process proposed in reference [8]. The detailed process of calculating the Fresnel coefficients is described in Appendix A1.*

_{i}In summary, the process of reverse ray tracing involves three steps: first, finding the eikonal function map of the device, *T*(X); second, determining the gradient of the eikonal function, ∇*T*(X), which represents the wave vector; and third, taking the wave vectors into the Fresnel transmission calculation. The illuminance is obtained and the geodesic of transformation optics can be traced by back propagation. The process is detailed in Table 1. The approach in this study can be applied to any transformation optics with positive definite material property tensor. All simulation programs are written in C language and the random number generator is from the GSL C library [22].

## 3. Simulation accuracy

The accuracy of trajectory and illuminance prediction calculated by reverse ray tracing are demonstrated in this section. The case of a wave collimator for converting cylinder waves into plane waves [23] was used as an example. This device serves as an optical element to reshape target illuminance and the exit surface-emitting ray angle [24].

The coordinate transformation of this device is illustrated in Fig. 2. The transformation maps a semi-circular domain with radius *a* into a rectangular domain with width of 2 × *W* and height of *H*. The transformation equation from the original space (*x, y, z*) to the transformed space (*x’*, *y’*, *z’*) can be expressed as

*x*= 0,

*y*= 0) of the coordinate along the

*z*-axis. As no transformation occurs along the

*z*-axis, we focus only on the ray tracing in the

*xy*-plane. In the simulation, the dimensions of the collimator are

*W*= 4 mm and

*H*= 2 mm along the

*x*- and

*y*-axes. Although this device is a 2D case, we calculated it in 3D to verify our 3D algorithm. The dimension along the

*z*-axis was assumed to be 0.4 mm in the calculation. To study the influence of the grid size on the ray tracing accuracy, three grid points, 51 × 51, 101 × 101, and 301 × 301, along the

*x*-axis and

*y*-axis were examined. We compared the ray trajectories and illuminances at the target calculated by Hamiltonian ray tracing and reserve ray tracing, respectively.

#### 3.1 Eikonal function and trajectories

The contour plot of the eikonal function, *T*(X), of the wave collimator is illustrated in Fig. 3. The number of grid points along the *x*- and *y*-axes is 301 × 301. The contour indicates the wavefronts that change from circular to nearly planar shapes. The black lines are ray trajectories obtained by reverse ray tracing. Figure 4 shows the ray trajectories emitted from the source with emitting angles of 10°, 30°, 45°, 60°, and 80°. The trajectories calculated by the reserve ray tracing approach with grid point numbers of 51 × 51, 101 × 101, and 301 × 301 are plotted as scatters. The analytical and Hamiltonian solutions are shown in black broken lines and red solid lines, respectively, for comparison. All trajectories obtained by analytical, Hamiltonian, and reverse ray tracing approaches are highly similar. However, a large deviation exists in the case of reverse ray tracing with a grid number of 51 × 51, especially for the 80° emitting angle. Error estimate of the sweeping method is *Ο*(*h*log(*h*)), where *h* is the grid size [25]. Accuracy decreases as grid size increases. The *O*(*h*log(*h*)) along the *x*-axis is 0.13 for the 51 × 51 grid points and 0.09 for the 101 × 101 grid points. For the 51 × 51 grid points, the absolute difference of the wavefront between two given adjacent nodes along the *x*-axis, which is |∇*T _{i,j,k}* –∇

*T*

_{i}_{± 1}

*|, is not significantly larger than the error estimate (0.13). A requirement of reverse ray tracing is that the absolute difference of the wavefront between two given adjacent nodes should be larger than the global error estimate. If the wavefront difference between two given adjacent points is larger than the global error estimate, then the accuracy of trajectories is less affected by a decrease in the number of grid points in the reverse ray tracing.*

_{,j,k}#### 3.2 Wave vector and illuminance

Furthermore, we also compare the illuminances of the tested device predicted by Hamiltonian ray tracing and reverse ray tracing. The number of grid points along the *x*- and *y*-axes is 101 × 101. The starting point of reverse ray tracing is [*W*Sin(θ_{lam}), H, 0], where θ_{lam} is the emitting angle. Angular distribution of θ_{lam} is Lambertian, which is generated by a random number generator [26]. After the map of *T*(X* _{i}*) of each node is obtained, the map of the wave vector ∇

*T*(X

*) is obtained for Fresnel transmission calculation. The Fresnel transmission is based on Appendix A1. Figure 5 shows the illuminances on a target, which is 0.001 mm above the exit surface of the wave collimator. Illuminance is the total luminous flux incidence per unit area. The illuminances obtained by Hamiltonian ray tracing and reverse ray tracing are plotted as a solid line and scatterplots, respectively. After transformation, the illuminance does not obey Lambert’s cosine law. Distribution of*

_{i}*W*Sin(θ

_{lam}) on the exit surface of the device along the

*x*-axis is uniform. The transformation converts a cylinder wave into a plane wave. The illuminances obtained by Hamiltonian and reverse ray tracing match well, and no significant difference between the illuminances from the two methods exists in different sizes of grid mesh.

## 4. Numerical examples

This section demonstrates the application of reverse ray tracing in two cases where the Hamiltonian ray tracing cannot be successfully solved by standard ODE solvers. The device in the first example can modify emission direction [27]. The issue of Hamiltonian ray tracing in Example 1 is illustrated in Appendix A2. The device in the second example can modify the emission direction [27] and polarization state [28].

#### 4.1 Example 1 – A 3D highly directive emitting device

An optical device that can perform directive emission [29] is adopted in this example. This device is used for uniform illumination and conversion of Lambertian emissions into highly directive emissions. This special case is not easily solved by Hamiltonian ray tracing. However, we successfully find the geodesics and illuminance distribution of this device by reverse ray tracing. The geometric configuration, transformed space, and material property tensors of this device are described in Appendix A2. The transformation is between the initial spherical coordinate (*r’, θ’*, *ϕ’*) and the transformed cylindrical coordinate (*ρ*, *ϕ*, *z*). A sphere with radius *R* is transformed to a cylinder slab with radius *W* and thickness *H*. The reverse ray tracing is performed in the Cartesian coordinates in the transformed space. The dimensions of the simulated device is 8 mm (2*W*) × 8 mm (2*W*) × 2 mm (*H*) along the *x*-, *y-*, and *z*-axes, respectively. The number of grid points is 101 × 101 × 101. The Lambertian point source is used in the simulation and its location is at [*x*, *y*, *z*] = [0, 0, 0]. The uniform initial locations of back propagation are distributed on the exit surface of the device. The number of simulated rays is 10^{6}. Figure 6(a) shows the trajectories and the eikonal function, *T*(X), of the emitting rays. Because the device is cylindrically symmetrical, only a cross section in the *xz*-plane with *y* equal to 0 is illustrated. The ray trajectories obtained by reverse ray tracing are plotted in blue lines and the analytical trajectories are plotted in solid black scatters. The trajectories obtained by reverse ray tracing match the analytical solution well except when the ray-emitting angle is larger than 80°. The contour plot of the eikonal function is illustrated in the background. The wavefronts, which can be represented by the contour of the eikonal function, are symmetrical and parallel to the *xy*-plane, indicating that the ray trajectories are converted into directive emission.

Based on Lambert’s cosine law [30], a practical approximation for the irradiance distribution of a Lambertian source is given by

where*θ*is the viewing angle and

*I*is the illuminance on axis at distance

_{0}‖∮|*r*from the emission source. According to the transformed equation [Appendix A2, Eq. (A6)], the viewing angles at the exit surface of the device approach 0. Illuminance at the distance

*r*from the emission source should be proportional to the thickness of the device. The distribution of ray locations of the exit surface is uniform. The Fresnel coefficient of each ray approaches 0.9 because the emission angle at the exit surface is smaller than 5°. According to these factors, the illuminance on the target should be uniform. According to our simulation, the illuminance in the

*xy*-plane on a target is indeed uniform, as shown in Fig. 6(b). The prediction of wave vectors based on our reserved ray tracing approach is correct.

#### 4.2 Example 2 – A 3D highly directive emitting device with polarization conversion

The optical device demonstrated in this section provides three major functionalities in illumination systems: the conversion of Lambertian illumination into uniform illumination; polarization control; and directive emission. A transformation optical device with these three functions was derived from references [28, 29]. Hamiltonian ray tracing is not successful in this case. The material property tensor of the device can be expressed as

_{1}is the Jacobian transformation matrix of the first transformation, which is the same as the transformation of the case in Section 4.1. It represents the transformation from the initial spherical coordinate (

*r’, θ’*,

*ϕ’*) to the transformed cylindrical coordinate (

*ρ*,

*ϕ*,

*z*). This transformation makes the luminance intensity distribution of the target proportional to the distance between the target plane and emitter.

*R*is the radius in spherical coordinates before transformation, and

*W*and

*H*are the radius and the thickness, respectively, of the cylinder slab in cylindrical coordinates after transformation. The second transformation, Λ

_{2}, rotates the polarization state of the rays by angle α under a cylindrical coordinate. In this case, the value of α is π/2.

For the simulation conditions, the device dimensions are 8 mm (2*W*) × 8 mm (2*W*) × 2 mm (*H*) along the *x*-, *y*-, and *z*-axes, respectively. The number of grid points is 101 × 101 × 101 and the number of simulated rays is 6 × 10^{5}. The initial locations of back propagation are uniformly distributed on the exit surface of the device because *W*Sin(θ_{lam}) is a uniform distribution. Angular distribution of θ_{lam} is Lambertian. Figure 7(a) shows the trajectories and the eikonal function, *T*(X), of the emitting rays. The blue lines are the ray trajectories obtained by reverse ray tracing and the solid black scatters are obtained from the analytical solution. The trajectories obtained by reverse ray tracing match the analytical solution well. The contour plot of the eikonal function is also illustrated. As the device is cylindrically symmetrical, only two cross sections in the planes with azimuth angles ϕ = 0° and ϕ = 60° are shown. The wavefronts, which can be represented by the contour of the eikonal function, approximate plane waves and are parallel to the *xy*-plane. Figure 7(b) shows the polar plot of the normalized ray numbers versus the ray emission angle at the exit surface of the device. Most of the ray-emitting angles are less than 5° except for the emitting angles of rays at the edge of the device, indicating that Lambertian light source transforms into directive emission. Figure 7(c) shows ray numbers with respect to their rotated angles. The legends in this figure represent the initial emission angle of each ray from the source. The rotation of emission rays indicates the rotation of the polarization state. The rotated angle is the difference between the azimuthal ϕ of each ray at the exit surface of the device and at the location near the source. Over 99.9 percent of rays rotate almost 90°, and only 0.03 percent of rays, which have large initial emission angles propagating from the point source, rotate less than 5°.

## 5. Conclusion

In this paper, we propose an alternative ray tracing method that is not dependent on finding solutions to ordinary differential motion equations for inhomogeneous and anisotropic materials. This method is useful in transformation optics that have indices of symmetric positive definite matrices and does not depend on an analytical expression of the material property tensor of the investigated system. In particular, it can be applied to special cases of transformation optics where Hamiltonian ray tracing is unusable. The sweeping method was employed to obtain the eikonal function of a system. To acquire the illuminance of the system, the wave vectors of the rays were further treated as the gradient of that eikonal function map. The rays are traced back from the target to the source; thus, reverse ray tracing provides an efficient approach to investigate the illuminance of a region of interest compared with traditional ray tracing methods. The device that converts a cylinder wave into a plane wave was adopted to demonstrate the accuracy of the reverse ray tracing. The illuminances obtained by Hamiltonian ray tracing and reverse ray tracing match. However, the absolute difference of the wavefront between two given adjacent nodes should be larger than the global error estimate. The second and third cases cannot be easily solved by Hamiltonian ray tracing. Instead, we solved both cases successfully by using the reverse ray tracing method. We expect that reverse ray tracing is another viable alternative to ray tracing based on ODEs.

## Appendix

#### A1. Reflection and transmission coefficients

In this section, we show the equations for calculating the transmission coefficient at the interface between transformation optics and air discussed in Section 2. It starts with the electric field vector, **E,** in the transformed space. Direction of time-averaged Poynting vector is proportional to ∇_{K}*J*, where **K** is the wave vector and *J* is a Hamiltonian–Jacobi equation. The electric field vector **E** is orthogonal to the Poynting vector and the corresponding magnetic field vector **H** field are given as Equations (6) and (7). **K**_{1} and **K**_{2} obtained from Equation (8) are the incident wave vector in the transformed space and the transmitted wave vector after refraction, respectively. Medium with material property tensors ε = μ is singly refracting. The transmitted electric field **E*** _{T}* and the corresponding magnetic field vector

**H**

*can be expressed as*

_{T}**n**is the unit normal to the interface and the indices

*s*and

*p*denote the two polarization states. To apply the boundary condition, two orthogonal vectors are defined aswhere the vectors

**t**

_{s}and

**t**

_{p}are tangential to the interface. The four Fresnel coefficients,

*a*,

_{Ts}*a*,

_{Tp}*a*, and

_{Rs}*a*, can be calculated by a linear matrix equation given as

_{Rp}**H***is the complex conjugate of

**H**. Finally, the intensity transmittance factor T [8] can be expressed as

*i*represents the incident wave.

#### A2. Coordinate transformation of Section 4.1, Example 1

In this case, a spherical wave is transformed to a plane wave. For the detailed transformation derivation, please refer to the work of Zhang [29]. In the present paper, we list the transformation process and material property tensors of the device medium. A sphere with radius *R* is transformed into a cylinder slab with radius *W* and thickness *H*. Figure 8 shows the transformation between the initial spherical coordinate (*r’, θ’*, *ϕ’*) and the transformed cylindrical coordinate (*ρ*, *ϕ*, *z*), which is given as

*R*is the radius of the sphere in the initial coordinate, and

*W*and

*H*are the radius and the thickness, respectively, of the cylinder slab in the transformed coordinate. The material property tensors in the Cartesian coordinates in the transformed space can be calculated using the following:

*Q*is the coordinate transformation matrix from the cylindrical coordinate system to the Cartesian coordinate system,

*ε*is the property tensor expressed in cylindrical coordinate, and

*ε”*is the property tensor expressed in Cartesian coordinates.

Here, the ODE solving issue in Hamiltonian ray tracing is illustrated. The material property tensor **M**(X) for the *xz* plane can be expressed as

As the emitting angle of a ray is 21°, the initial wave vector in the original space can be expressed as K_{o} = [Sin(θ), 0, Cos(θ)] = [0.36, 0,0.93]. The corresponding wave vector, K_{t}, in the transformed space is obtained using the following equations:

**n**is a unit normal to the device surface and J is the Hamilton–Jacobi equation given in Eq. (1) in Section 2; where K

_{t}= [K

_{x,}K

_{y,}K

_{z}] = [0.36, 0, 144i] = βx + iαz at [x, y, z] = [0,0,0.01]. The imaginary part of the wave vector represents the amplitude attenuation of the wave while the real part represents the directions of the propagation wave. The initial condition of the ODEs is [x, y, z, K

_{x}, K

_{y}, K

_{z}] = [0, 0, 0.01, 0.36, 0, 0]. The resulting transformed wave vector does not exactly satisfy the Hamilton–Jacobi equation. Moreover, the Jacobian matrix of the ODEs can also be critical to reliability and efficiency. In this case, the condition number of the Jacobian matrix in the solving process is large (>10

^{5}); hence, the Jacobian matrix is close to singular or badly scaled. Figure 9 shows the ray trajectories calculated using both Hamiltonian ray tracing and reverse-ray tracing. The emitting angle of the ray is 21°. The ODEs in Eq. (A9) were solved using the commands ode23s and ode45 implemented in MATLAB. The ray trajectory calculated using the Hamiltonian ray tracing cannot be solved even with a smaller step size in the integration. The real ray tracing approach [31] may also be introduced for solving this issue.

## Acknowledgment

This study was financially supported by the National Science Council of Taiwan under grants NSC 101-2221-E-006-210-MY2 and NSC 100-2120-M-009-010.

## References and links

**1. **B. Zhang, Y. Luo, X. Liu, and G. Barbastathis, “Macroscopic Invisibility Cloak for Visible Light,” Phys. Rev. Lett. **106**(3), 033901 (2011). [CrossRef] [PubMed]

**2. **X. Chen, Y. Luo, J. Zhang, K. Jiang, J. B. Pendry, and S. Zhang, “Macroscopic invisibility cloaking of visible light,” Nat. Commun. **2**, 176 (2011). [CrossRef] [PubMed]

**3. **J. C. Halimeh and M. Wegener, “Photorealistic rendering of unidirectional free-space invisibility cloaks,” Opt. Express **21**(8), 9457–9472 (2013). [CrossRef] [PubMed]

**4. **D. Schurig, J. B. Pendry, and D. R. Smith, “Calculation of material properties and ray tracing in transformation media,” Opt. Express **14**(21), 9794–9804 (2006). [CrossRef] [PubMed]

**5. **J. C. Halimeh, R. Schmied, and M. Wegener, “Newtonian photorealistic ray tracing of grating cloaks and correlation-function-based cloaking-quality assessment,” Opt. Express **19**(7), 6078–6092 (2011). [CrossRef] [PubMed]

**6. **N. Kundtz and D. R. Smith, “Extreme-angle broadband metamaterial lens,” Nat. Mater. **9**(2), 129–132 (2010). [CrossRef] [PubMed]

**7. **A. S. Glassner, ed., *An introduction to ray tracing* (Academic Press Ltd., 1989), p. 327.

**8. **M. Sluijter, D. K. G. de Boer, and J. J. M. Braat, “General polarized ray-tracing method for inhomogeneous uniaxially anisotropic media,” J. Opt. Soc. Am. A **25**(6), 1260–1273 (2008). [CrossRef] [PubMed]

**9. **A. Akbarzadeh and A. J. Danner, “Generalization of ray tracing in a linear inhomogeneous anisotropic medium: a coordinate-free approach,” J. Opt. Soc. Am. A **27**(12), 2558–2562 (2010). [CrossRef] [PubMed]

**10. **L. Shampine and C. Gear, “A User’s View of Solving Stiff Ordinary Differential Equations,” SIAM Rev. **21**(1), 1–17 (1979). [CrossRef]

**11. **P. Deuflhard, “A modified Newton method for the solution of ill-conditioned systems of nonlinear equations with application to multiple shooting,” Numerische Mathematik **22**(4), 289–315 (1974). [CrossRef]

**12. **J. A. Ogilvy, “A layered media model for ray propagation in anisotropic inhomogeneous materials,” Appl. Math. Model. **14**(5), 237–247 (1990). [CrossRef]

**13. **U. Leonhardt and T. G. Philbin, “Transformation Optics and the Geometry of Light,” in *Progress in Optics*, E. Wolf, ed. (Elsevier, 2009), Chap. 2, pp. 69–152.

**14. **J. Qian, Y.-T. Zhang, and H.-K. Zhao, “A Fast Sweeping Method for Static Convex Hamilton–Jacobi Equations,” J. Sci. Comput. **31**(1-2), 237–271 (2007). [CrossRef]

**15. **A. Capozzoli, C. Curcio, A. Liseno, and S. Savarese, “Two-dimensional fast marching for geometrical optics,” Opt. Express **22**(22), 26680–26695 (2014). [CrossRef] [PubMed]

**16. **R. J. Koshel, “Introduction and Terminology,” in *Illumination Engineering* (John Wiley & Sons, Inc., 2013), pp. 1–30.

**17. **R. Leutz and H. P. Annen, “Reverse ray-tracing model for the performance evaluation of stationary solar concentrators,” Sol. Energy **81**(6), 761–767 (2007). [CrossRef]

**18. **K. E. Moore, R. F. Rykowski, and S. Gangadhara, “Reverse radiance: a fast accurate method for determining luminance,” Proc. SPIE **8485**, 84850I (2012).

**19. **J.-M. Mirebeau, “Anisotropic Fast-Marching on Cartesian Grids Using Lattice Basis Reduction,” SIAM J. Numer. Anal. **52**(4), 1573–1599 (2014). [CrossRef]

**20. **S. Jbabdi, P. Bellec, R. Toro, J. Daunizeau, M. Pélégrini-Issac, and H. Benali, “Accurate Anisotropic Fast Marching for Diffusion-Based Geodesic Tractography,” Int. J. Biomed. Imaging **2008**, 320195 (2008). [CrossRef] [PubMed]

**21. **M. M. Crosskey, A. T. Nixon, L. M. Schick, and G. Kovačič, “Invisibility cloaking via non-smooth transformation optics and ray tracing,” Phys. Lett. A **375**(18), 1903–1911 (2011). [CrossRef]

**22. **B. Gough, *GNU Scientific Library Reference Manual, *3^{rd}. ed. (2009).

**23. **D.-H. Kwon and D. H. Werner, “Transformation optical designs for wave collimators, flat lenses and right-angle bends,” New J. Phys. **10**(11), 115023 (2008). [CrossRef]

**24. **H.-C. Chen, J.-Y. Lin, and H.-Y. Chiu, “Rectangular illumination using a secondary optics with cylindrical lens for LED street light,” Opt. Express **21**(3), 3201–3212 (2013). [CrossRef] [PubMed]

**25. **H. K. Zhao, “A fast sweeping method for Eikonal equations,” Math. Comput. **74**(250), 603–628 (2005). [CrossRef]

**26. **R. J. Koshel, “Sampling, Optimization, and Tolerancing,” in *Illumination Engineering* (John Wiley & Sons, Inc., 2013), pp. 251–297.

**27. **J. J. Zhang, Y. Luo, S. Xi, H. S. Chen, L. X. Ran, B. I. Wu, and J. A. Kong, “Directive emission obtained by coordinate transformation,” Prog. Electromagnetics Res. **81**, 437–446 (2008). [CrossRef]

**28. **D. H. Kwon and D. H. Werner, “Polarization splitter and polarization rotator designs based on transformation optics,” Opt. Express **16**(23), 18731–18738 (2008). [CrossRef] [PubMed]

**29. **J. J. Zhang, Y. Luo, S. Xi, H. S. Chen, L. X. Ran, B. I. Wu, and J. A. Kong, “Directive emission obtained by coordinate transformation,” Prog. Electromagnetics Res. **81**, 437–446 (2008). [CrossRef]

**30. **I. Moreno, M. Avendaño-Alejo, and R. I. Tzonchev, “Designing light-emitting diode arrays for uniform near-field irradiance,” Appl. Opt. **45**(10), 2265–2272 (2006). [CrossRef] [PubMed]

**31. **V. Vavrycuk, “Real ray tracing in anisotropic viscoelastic media,” Geophys. J. Int. **175**(2), 617–626 (2008). [CrossRef]