Ray Tracing Method of Gradient Refractive Index Medium Based on Refractive Index Step

a large range of refractive index gradient changes, such as around the shock layer of an aircraft. Abstract: For gradient refractive index media with large refractive index gradients, traditional ray tracing methods based on reﬁned elements or spatial geometric steps have problems such as low tracing accuracy and efﬁciency. The ray tracing method based on refractive index steps proposed in this paper can effectively solve this problem. This method uses the refractive index step to replace the spatial geometric step. The starting point and the end point of each ray tracing step are on the constant refractive-index surfaces. It avoids the problem that the traditional tracing method cannot adapt to the area of sudden change in the refractive index and the area where the refractive index changes sharply. Therefore, a suitable distance can be performed in the iterative process. It can achieve high-efﬁciency and precise ray tracing in areas whether the refractive index changes slowly or sharply. According to the comparison of calculation examples, this method can achieve a tracing accuracy of 10 − 5 mm. The speed and precision of ray tracing are better than traditional methods.


Introduction
The optical phenomenon in the gradient refractive index medium, as a ubiquitous objective physical phenomenon in nature, has been noticed for a long time. The Maxwell fisheye lens, the Wood lens, the Luneburg lens, the Gradient index rod lens, etc., can all be regarded as historical signs of the research on gradient index optics [1]. Due to the change in the refractive index of the medium, the original path of the light propagating in it is changed, which results in the deflection or phase change of the light and it also results in the image offset, blur, and jitter. Therefore, ray tracing in gradient refractive index media has become a research hotspot.
The ray tracing through graded refractive index media is based on the solution of the ray formula: d ds n dr ds = ∇n (1) where s is the distance along the ray path, r is the position of the ray, and n is the index of refraction. Regarding Formula (1), Lucian proposed a universal numerical solution in 1968 [2]. Since then, the Runge-Kutta method, Euler method, and Taylor series expansion method have been applied to solve Formula (1), which have laid the foundation of gradient refractive index media ray tracing for a long time [3,4]. Most ray tracing methods are based on the three methods which have been listed on references [5][6][7]. In the ray tracing process, 2 of 9 these methods are all iteratively based on the refined elements [8][9][10] or spatial geometric steps [4,11,12]. When ray tracing is performed in a gradient refractive index medium with a gentle refractive index gradient, the accuracy of the ray tracing can be controlled by changing the geometric step. However, when encountering a gradient refractive index field with a sharp change (such as the shock wave of a high-speed aircraft [13][14][15]), it is necessary to adjust the geometric step according to the refractive index gradient aimed to be self-adaptive to the local refraction index gradient and grid geometry scale [7,16]. Even so, it is still possible to cross the region where the refractive index changes sharply in one iterative step. This situation often causes the refraction effect of the refractive index layer whose refractive index changes sharply to be ignored. In order to ensure the accuracy of the entire ray tracing path, the geometrical step-based ray tracing method (RTSGS) generally needs to set a smaller geometric tracing step. It will decrease efficiency in the tracing area where the refractive index changes slowly. Therefore, it will increase useless calculations [17]. This paper proposes a ray tracing method based on refractive index steps (RTRIS). In this method, the refractive index step ∆n is used to replace the geometric step ∆s of the traditional ray tracing method. The advantage of this method is that each iterative process takes equal refractive index step ∆n as the iterative step. In areas where the refractive index changes slowly, a relatively large distance will be performed in the iterative process to achieve high-efficiency ray tracing. In areas where the refractive index changes sharply or suddenly, the iterative step will reduce automatically to capture the sudden changes in the refractive index effectively. It can achieve refined tracing and improve ray tracing accuracy. The application of this method avoids the problem of being unable to capture regions of sudden changes in the refractive index after applying the step size regulation function. It can also effectively avoid the problem of insufficient tracing accuracy or low efficiency caused by unreasonable geometric step setting. By using the proposed method, the efficiency of ray tracing can be achieved while the accuracy of ray tracing can be improved. Thus, the method realizes the adaptive process of the gradient index medium in the ray tracing process.

The Method of Ray Tracing Based on Refractive Index
Step According to the differential geometry thought, the gradient index medium is divided into many small areas by the constant index surface. The medium between two constantindex surfaces is regarded as an equal refractive index medium. The value of the refractive index for this medium is replaced as the average of the two constant-index surfaces. Those constant-index surfaces are regarded as refractive surfaces. Snell's law is used to calculate the direction vector of the refracted ray after the ray passes through the constant index surface. So, Snell's law is the foundation of RTSGS [15].

Space Snell's Law
Snell's law in vector form is much more useful in the process of ray tracing in threedimensional space [18]. According to the reference [6], we can calculate Snell's law in vector form.
The plane Snell's law can usually be expressed as: In this formula, n and n are the refractive index of the medium on the incident side and refracted side, respectively. I and I are the incident angle and refraction angle of the ray, respectively.
Combined with Figure 1, Formula (2) can be modified as: 3 of 9 where A 0 is the unit vector along the incident ray, A 0 is the unit vector along the refracted light, and N is the unit vector along the normal.
In this formula, n and n are the refractive index of the medium on the incident side and refracted side, respectively. I and I  are the incident angle and refraction angle of the ray, respectively.
Combined with Figure 1, Formula (2) can be modified as: where 0 A is the unit vector along the incident ray, ' 0 A is the unit vector along the refracted light, and N is the unit vector along the normal. The refracted ray vector with length n and the incident ray vector with length n are written as  A and A , respectively: Formula (4) shows that vector  A -A and vector N are collinear vectors. Rewrite Formula (4) as: where P is vector deflection constant. Product both sides of Formula (4) with N [6]: when ' nn  , the vectors A' -A and N are parallel as shown in Figure 1b. When nn   , the two vectors are antiparallel as shown in Figure 1c.
When the refractive index of the two media and the incident angle of the light are known, cos nI  can be replaced with: Therefore, the refraction law of vector form is [6]:

Ray Tracing in the Gradient Refractive Index Medium
The ray tracing method in the gradient refractive index medium proposed in this paper is based on the refractive index step. The tracing process is shown in Figure 2. In this figure, The refracted ray vector with length n and the incident ray vector with length n are written as A and A, respectively: Formula (4) shows that vector A -A and vector N are collinear vectors. Rewrite Formula (4) as: where P is vector deflection constant. Product both sides of Formula (4) with N [6]: when n > n, the vectors A -A and N are parallel as shown in Figure 1b. When n < n, the two vectors are antiparallel as shown in Figure 1c.
When the refractive index of the two media and the incident angle of the light are known, n cos I can be replaced with: Therefore, the refraction law of vector form is [6]:

Ray Tracing in the Gradient Refractive Index Medium
The ray tracing method in the gradient refractive index medium proposed in this paper is based on the refractive index step. The tracing process is shown in Figure 2. In this figure, n 0 , n 1 , n 2 , n 3 are the equal refractive index surfaces, and n 3 = n 2 + t = n 1 + 2 × t = n 0 + 3 × t (t is the refractive index step). The medium between n 0 and n 1 is regarded as an equal refractive index medium, and the value of refractive index is n 0 + n 1 2 . Similarly, the refractive index of the medium between n 1 and n 2 is n 1 + n 2 2 . The refractive index of the medium between n 2 and n 3 is n 2 + n 3 2 . P 0 , P 1 , P 2 , P 3 are the intersection points of ray L and constant refractive-index surfaces n 0 , n 1 , n 2 , n 3 . When the ray L enters the medium from point P 0 along direction e p 0 p 1 , its actual propagation path is L 1 . Since the medium between n 0 and n 1 is regarded as an equal refractive index medium, the Line segment L 1 is approximately regarded as the equivalent propagation path of ray between n 0 and n 1 In the same way, the actual ray propagation path L 2 , L 3 , is approximately replaced with L 2 , L 3 , respectively. tive index is 01 2 nn  . Similarly, the refractive index of the medium between 1 n and 2 n is 12 2 nn  . The refractive index of the medium between 2 n and n 3 is 2 nn  23 . 0 P , 1 P , 2 P , 3 P are the intersection points of ray L and constant refractive-index surfaces 0 n , 1 n , 2 n , 3 n . When the ray L enters the medium from point 0 P along direction 01 pp e , its actual propagation path is ' L 1 . Since the medium between 0 n and 1 n is regarded as an equal refractive index medium, the Line segment 1 L is approximately regarded as the equivalent propagation path of ray between 0 n and 1 n . In the same way, the actual ray propaga- L is approximately replaced with 2 L , 3 L , respectively. The constant-index surface i n is regarded as refractive surface. According to Formula (8), the direction vector of the refracted ray can be expressed as: P  In the same way: P  So far, the iterative formula, which is based on the refractive index step for ray tracing in gradient index media, can be obtained. According to the iterative formula, we can calculate the OPL (optical path length): The constant-index surface n i is regarded as refractive surface. According to Formula (8), the direction vector of the refracted ray can be expressed as: where e p 1 p 2 and e p 2 p 3 are the direction vectors of the incident ray and refractive ray, and the N 12 is the normal of constant refractive-index surface. At this point, we can get the direction vector of the refracted ray. Since the refractive index step is used instead of the geometric step in this solution, the coordinate of the next refraction point P 2 cannot be determined directly according to e p 1 p 2 . Therefore, we propose to solve the coordinate of point P 2 by solving the combined formula of straight line P 1 P 2 and the constant refractive-index surface n 2 : In the same way: e p 2 p 3 = e p 1 p 2 + P 23 N 23 (11) So far, the iterative formula, which is based on the refractive index step for ray tracing in gradient index media, can be obtained. According to the iterative formula, we can calculate the OPL (optical path length): where S P i P i + 1 is the distance of P i and P i + 1 .

Algorithm Verification
For most of the non-uniform refractive index medium, the analytical solution of the ray trajectory cannot be derived from the ray formula. However, the analytical solution can be acquired for some special refractive index distributions [1]. We compare the analytical solutions of two representative refractive index distributions with the numerical solutions obtained by the proposed algorithm. All the following ray tracing processes take the + z axis as the optical axis direction, and in order to compare the tracing accuracy of RTRIS and RTSGS, the iterative steps of the two tracing methods are the same in the tracing process.

Example 1
The distribution function of the refractive index field is: where α is the medium distribution constant. The linear function distribution of the axial gradient refractive index is the simplest form of axial gradient refractive index distribution. Under the condition of paraxial, the analytical solution of the ray formula is: where p 0 , q 0 are the optical direction cosine of the ray. When x = 0, the constant a = n 0 is the refractive index of the incident surface. Since the refractive index of glass and optical fiber is generally around 1.5, so take n 0 = 1.5 and α = − 0.0005. When the incident point coordinate is (0, 0, 0) and incident ray direction vector is (0, 0.0001, 1), Formula (14) can be simplified as: When tracing length is 2000 mm, the optical path analytical solution is: Use the algorithm proposed in this paper to trace the gradient index medium with a linear function of refractive index distribution from 0 to 2000 mm, we can obtain the optical path L1 = 2.000000027 × 10 3 . It can be found that the order of tracing accuracy can reach 10 −6 , and the accuracy can completely meet the demand of engineering accuracy.
In Figure 3, the RTRIS represents the error between the optical path obtained by the ray tracing method based on the refractive index step and the actual optical path from 0 to 2000 mm. In addition, the RTSGS represents the error between the optical path obtained by the ray tracing method based on the geometric step and the actual optical path. It is showing that when the tracing distance is in the range of 0 to 400 mm, the errors of the two tracing methods are almost the same, and the error value is small. When the tracing distance is longer than 400 mm, the error growth rate of RTRIS increases rapidly, and its error curve is similar to exponential curve. Finally, when the tracing distance is 2000 mm, the error value reaches 1.47 × 10 −5 mm. For RTSGS, when the tracing distance continues to increase, the error value increases slightly, and the growth rate is small. The error curve is flat. When tracing distance reaches 2000 mm, the error value is only 2.2 × 10 −6 mm, which is only 0.15 times as much as RTRIS. The primary source of error is the accumulation of computer errors and fitting errors of the constant refractive-index surfaces. It is showing that the ray tracing method based on the gradient refractive index step is more accurate and stable than that based on the space step for the gradient refractive index media distribution in accordance with Formula (11). is flat. When tracing distance reaches 2000 mm, the error value is only , which is only 0.15 times as much as RTRIS. The primary source of error is the accumulation of computer errors and fitting errors of the constant refractive-index surfaces. It is showing that the ray tracing method based on the gradient refractive index step is more accurate and stable than that based on the space step for the gradient refractive index media distribution in accordance with Formula (11).

Example 2
The distribution function of refractive index field is: The distribution function of the refractive index field shown in Formula (17) is one of the typical distribution forms that can be obtained by ion diffusion technology. When 1.5 n  0 and   = 0 .00049 , the refractive index distribution along the z axis is shown in Figure 4. It can be seen that in the range of 0 to 1000 mm, the refractive index changes smoothly, and the gradient of the refractive index is small. As the tracing distance increases, the refractive index changes more and more rapidly, and the refractive index decreases rapidly. The gradient of the refractive index rises sharply. The refractive index distribution is extremely uneven.

Example 2
The distribution function of refractive index field is: The distribution function of the refractive index field shown in Formula (17) is one of the typical distribution forms that can be obtained by ion diffusion technology. When n 0 = 1.5 and α = − 0.00049, the refractive index distribution along the z axis is shown in Figure 4. It can be seen that in the range of 0 to 1000 mm, the refractive index changes smoothly, and the gradient of the refractive index is small. As the tracing distance increases, the refractive index changes more and more rapidly, and the refractive index decreases rapidly. The gradient of the refractive index rises sharply. The refractive index distribution is extremely uneven. is flat. When tracing distance reaches 2000 mm, the error value is only 2.2 10 mm  , which is only 0.15 times as much as RTRIS. The primary source of error is the accumulation of computer errors and fitting errors of the constant refractive-index surfaces. It is showing that the ray tracing method based on the gradient refractive index step is more accurate and stable than that based on the space step for the gradient refractive index media distribution in accordance with Formula (11).

Example 2
The distribution function of refractive index field is: The distribution function of the refractive index field shown in Formula (17) is one of the typical distribution forms that can be obtained by ion diffusion technology. When 1.5 n  0 and   = 0 .00049 , the refractive index distribution along the z axis is shown in Figure 4. It can be seen that in the range of 0 to 1000 mm, the refractive index changes smoothly, and the gradient of the refractive index is small. As the tracing distance increases, the refractive index changes more and more rapidly, and the refractive index decreases rapidly. The gradient of the refractive index rises sharply. The refractive index distribution is extremely uneven.  Under paraxial conditions, the analytical solution of the ray equation is: When the incident point coordinate is (0, 0, 0) and the incident ray direction vector is (0, 0.0001, 0), Formula (18) can be simplified as: When the tracing length is 2000 mm, the optical path analytical solution is: Using the algorithm proposed in this paper to trace the gradient refractive index medium from 0 to 2000 mm, it is possible to obtain the optical path length L1 = 2.3961413 × 10 3 .
The curve in Figure 5 shows the relationship between the geometric step length and the tracing distance of the two methods. The number of iterations for both of these two methods is 6821. For the traditional ray tracing method, the geometric step length is constant. For the method proposed in this paper, the geometric step length can be adjusted adaptively according to the refractive index distribution. In the first half of the tracing trajectory, the refractive index of the medium changes slower. Therefore, the geometric step size is larger. As the tracing trajectory gradually approaches 2000 mm, the refractive index of the medium changes intensify, and its geometric step size decreases rapidly to achieve higher accuracy.
When the incident point coordinate is (0, 0, 0) and the incident ray direction vector is (0, 0.0001, 0), Formula (18) can be simplified as: . (20) Using the algorithm proposed in this paper to trace the gradient refractive index medium from 0 to 2000 mm, it is possible to obtain the optical path length The curve in Figure 5 shows the relationship between the geometric step length and the tracing distance of the two methods. The number of iterations for both of these two methods is 6821. For the traditional ray tracing method, the geometric step length is constant. For the method proposed in this paper, the geometric step length can be adjusted adaptively according to the refractive index distribution. In the first half of the tracing trajectory, the refractive index of the medium changes slower. Therefore, the geometric step size is larger. As the tracing trajectory gradually approaches 2000 mm, the refractive index of the medium changes intensify, and its geometric step size decreases rapidly to achieve higher accuracy. In Figure 6, the RTRIS represents the error between the optical path obtained by the ray tracing method based on the refractive index step and the actual optical path from 0 to 2000 mm. The RTSGS represents the error between the optical path obtained by the ray tracing method based on the geometric step and the actual optical path. Comparing the tracing errors of the two tracing methods in Figure 5, it can be found that the tracing accuracy of the ray tracing method based on the refractive index step is significantly better In Figure 6, the RTRIS represents the error between the optical path obtained by the ray tracing method based on the refractive index step and the actual optical path from 0 to 2000 mm. The RTSGS represents the error between the optical path obtained by the ray tracing method based on the geometric step and the actual optical path. Comparing the tracing errors of the two tracing methods in Figure 5, it can be found that the tracing accuracy of the ray tracing method based on the refractive index step is significantly better than that of the traditional ray tracing method based on the geometric step. When the tracing distance is 2000 mm, the cumulative error of the RTRIS is 2.5 × 10 −4 mm, and the cumulative error of the RTSGS is up to 3.13 × 10 −4 mm.
Combine  in the initial stage of ray tracing, due to the slow change of the refractive index of the medium, the errors of the two methods are almost equal although the geometric step length of RTSGS is relatively large. As the tracing distance gradually increases, the advantages of the RTSGS gradually become apparent. Because the refractive index of the medium changes smoothly and the geometric step of RTSGS is larger, there are fewer iterations and smaller cumulative errors. When the tracing distance continues to increase, especially when it is close to 2000 mm, the refractive index of the medium decreases rapidly. The tracing step length for RTRIS cannot adaptively change, the error increases rapidly. RTSGS can adjust the geometric step length adaptively. Therefore, the geometric step is adjusted to be smaller in this section, which can better capture the refractive index change of the medium and improve the accuracy. It can be seen that the essence of RTSGS is to reasonably allocate the geometric step length in the ray tracing process.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 9 of 10 than that of the traditional ray tracing method based on the geometric step. When the tracing distance is 2000 mm, the cumulative error of the RTRIS is mm   4 2. 5 10 , and the cumulative error of the RTSGS is up to mm   4 3.13 10 . Combine  in the initial stage of ray tracing, due to the slow change of the refractive index of the medium, the errors of the two methods are almost equal although the geometric step length of RTSGS is relatively large. As the tracing distance gradually increases, the advantages of the RTSGS gradually become apparent. Because the refractive index of the medium changes smoothly and the geometric step of RTSGS is larger, there are fewer iterations and smaller cumulative errors. When the tracing distance continues to increase, especially when it is close to 2000 mm, the refractive index of the medium decreases rapidly. The tracing step length for RTRIS cannot adaptively change, the error increases rapidly. RTSGS can adjust the geometric step length adaptively. Therefore, the geometric step is adjusted to be smaller in this section, which can better capture the refractive index change of the medium and improve the accuracy. It can be seen that the essence of RTSGS is to reasonably allocate the geometric step length in the ray tracing process.

Conclusions
The ray tracing method based on the refractive index step proposed in this paper has extensive adaptability to the axial gradient refractive index, especially the inhomogeneous refractive index medium with a large refractive index gradient. It can efficiently solve the problem that the traditional ray tracing method based on geometric step cannot adapt to the inhomogeneous refractive index medium with large refractive index gradient or sudden change and improve ray tracing accuracy. At the same time, the proposed method can achieve high-efficiency ray tracing in a non-uniform refractive index medium with small refractive index change rate due to the refractive index step used in the recursion process. By using the proposed method, the efficiency of ray tracing can be achieved while the accuracy of ray tracing can be improved. Thus, the method realizes the adaptive process of the gradient index medium in the ray tracing process.
This method, however, has some problems which have to be solved. The most important one is that the constant refractive-index surfaces need to be fitted during the iterative process, so as to solve the iteration point through the fitting formula. The current method is to use high-order polynomials to fit the constant refractive-index surfaces, which might introduce a certain error. If the polynomial order is too high, it will cause Runge phenomenon. Therefore, different fitting models need to be considered in order to improve the ray tracing accuracy for different non-uniform refractive index media. Although there are some shortcomings in this method, we believe that this method will have a good application in ray tracing in gradient index media, especially in shock waves.

Conclusions
The ray tracing method based on the refractive index step proposed in this paper has extensive adaptability to the axial gradient refractive index, especially the inhomogeneous refractive index medium with a large refractive index gradient. It can efficiently solve the problem that the traditional ray tracing method based on geometric step cannot adapt to the inhomogeneous refractive index medium with large refractive index gradient or sudden change and improve ray tracing accuracy. At the same time, the proposed method can achieve high-efficiency ray tracing in a non-uniform refractive index medium with small refractive index change rate due to the refractive index step used in the recursion process. By using the proposed method, the efficiency of ray tracing can be achieved while the accuracy of ray tracing can be improved. Thus, the method realizes the adaptive process of the gradient index medium in the ray tracing process.
This method, however, has some problems which have to be solved. The most important one is that the constant refractive-index surfaces need to be fitted during the iterative process, so as to solve the iteration point through the fitting formula. The current method is to use high-order polynomials to fit the constant refractive-index surfaces, which might introduce a certain error. If the polynomial order is too high, it will cause Runge phenomenon. Therefore, different fitting models need to be considered in order to improve the ray tracing accuracy for different non-uniform refractive index media. Although there are some shortcomings in this method, we believe that this method will have a good application in ray tracing in gradient index media, especially in shock waves.