A mirage or *fata* *morgana* is typically an upside-down “mirror” image of a scenery in deserts, over sun-heated roads, or above bodies of water. When the temperature gradient of air is large, as can happen near a surface, it results in a large gradient of the density and the refractive index as a function of height. Mirages appear when light travels through a medium with a gradient in its refractive index and, therefore, get bent towards the higher values, generating reflected images. A computer program that simulates mirages above water using the method of ray tracing has been developed and is presented here in detail for educational purposes. Results on the effect are shown by simulated images for various water-air temperature-difference cases with values ranging from 5 to 0 °C. Comparison of the simulations to a real-life scenario at Lake Balaton, Hungary has also been provided.

## I. INTRODUCTION

When light travels through a medium with a non-uniform refractive index, it bends towards the higher values. This effect can produce upside down “mirror” images of a scenery in deserts, over sun-heated roads, or above bodies of water. The so-called *mirage* or *fata* *morgana* appears in stable weather conditions if the temperature difference between two layers is high-enough to create a large refractive index gradient.

Due to the high heat capacity of water relative to the ground or air, in the morning lakes and sea surfaces have a higher temperature than ambient air. If one goes to the shore and looks at an object across the water, a mirror-like downward reflection appears, which in reality is caused by a mirage effect.^{1} This is called an *inferior mirage*, in contrast with the *s*uperior, where the “mirror” image appears above.^{1–5} A special case of the superior mirage is the so-called *hillingar*, in which objects normally not visible due to Earth's curvature rise into sight.^{6}^{,} *Inferior mirages* can also be experienced on hot summer days over heated road surfaces^{3,7,8} or next to heated walls.^{9} Double “reflections” are also possible.^{10}

We implemented a computer program that simulates mirages above water using the method of ray tracing. We built a model for the temperature profile based on temperature measurements, from which we calculated the refractive index of air. Rays of light are traced by solving the eikonal equation using different Runge–Kutta methods. Given a picture, a physical setup, and a value for both the temperature of the water body and that of the ambient air sufficiently far from the surface (where it can be considered constant), our program can realistically reproduce pictures of photographed mirages. The numerical code is capable of simulating the “mirrored” images of objects at large distances, from practically zero to tens of kilometres away. Because of this, the curvature of Earth was also taken into account during the simulations, which affects the visibility of distant objects depending on the observer's height as in Fig. 1.

We made a comparison between a real photographed mirage and a simulated one with the same physical parameters. In our example, the distance between the observer and the scenery was about 17 km with water temperature equalling $ 5 \xb0$ C and ambient air temperature $ 1 \xb0$ C.

The structure of this paper is as follows: We present the optical setup and the effect for different cases in Sec. II. The model for the temperature-profile above the water is presented in Sec. III. We explain the ray tracing and the numerical solution of the eikonal equation in Sec. IV. Finally, results by numerical simulations are shown in Sec. V, which are compared to scenarios in nature.

In this paper, we will refer to visible physical phenomena (virtual or real) as an *image. Object* will mean a physical object whose image is formed by some physical transformation of light. By *picture*, we mean a photograph or a digital picture.

## II. REFRACTION OF LIGHT

Here, we model the optical conditions of mirage-formation above water surfaces. The air temperature well above the water surface is the ambient temperature $ T a = 1 \xb0$ C, the water temperature is $ T w = 5 \xb0$ C, and the air pressure is *p* = 101 kPa. All these quantities are assumed to be uniform.

In the following figures, the method of ray tracing was used (explained later in detail). Rays are traced in reverse: starting at the observer's eyes, going towards the object (red arrow in Fig. 2) that the observer sees. This has been a standard method in computer graphics for decades. It is important because otherwise we would need to trace every single ray that is emitted from a given point on the object, going in all possible directions. Most of these rays never reach the observer, since they are travelling in other directions, so there is no point in calculating those. However, from optics, we know that ray paths are reversible, so we can simply trace back the rays that land in our eyes, as if they were emitted from there, which is a computationally much more efficient method.

Figure 2 shows how rays travel with various viewing angle settings in a medium with a non-uniform refractive index. The object we see is located on the right side (red arrow that starts at the Earth's surface). From each point on the object, rays go in many directions, but the figure only shows those rays that reach the observer, which is located on the left side of Fig. 2 at a height of 2.7 m. In normal conditions, there is only one ray from each point reaching the observer; however, in a medium with a non-uniform refractive index, two paths are possible. Rays shown in red produce a real image of the object that is in an upright position (red arrow), while the ones with blue color give us an inverted virtual image (dashed blue line) that appears to be below the Earth's surface. Only the blue part of the original object is seen in the virtual image. Red-and-blue dashed line represents two important *extrema*, which are distinct, but non-resolvable in this figure: The reference line and “mirror horizon.” The reference line (which is the same as the blue dashed lines in Figs. 9–11) is defined by the case with no temperature difference between air and water, so light rays are straight lines. The “mirror horizont” is defined as the ray that contacts the lowest point on the object.

The black dashed-line path corresponds to the ray starting horizontally (with 0° angle). This can also get bent slightly upwards; however, in this case where the ray starts 2.7 m above the ground, the amount of bending is negligible. As seen in the figure, small portion of the object's bottom will not be visible at all when there is a mirage. Also, the bottom part of the upright image is not the same as the “bottom” of the inverted one. The latter is more distorted, as indicated by the rays on each side of the reference line (see the vanishing line in Ref. 3). Note that the curvature of Earth (surface indicated by solid black curve) might also crop the bottom of the images.

We investigated the dependency on the observer's elevation. Figure 3 shows light rays starting from different heights. Physical parameters (temperature distribution and pressure) are the same as in the case of Fig. 2. The black arrow represents the object.

Thick lines (with initial angle $ \u2212 0.16 \xb0$) represent the steepest rays that can be seen “reflected” (bent in reality) from each height under the current conditions. The lower the point of view, the taller the objects one can see in the upside-down image. Notice that, as the height of the observer increases, the thick lines intercept a lower point on the right side of the figure, telling us that less of the height of the object can be seen in the inverted image. This is nicely illustrated by pictures taken in Ref. 1.

Distance between the object and the observer is also important; it affects the size of the mirage. Figure 4 shows observational distances of 1500, 2000, and 2500 m. The further we are from the object, the higher the point on the object that is visible in the inverted image. In the case of the largest distance (2500 m), the last ray (indicated with thick line) crosses the $ 0 \xb0$ one, while for the other scenarios, it can only reach lower heights. Note that this change in relative sizes is not just a proportionality, and it is not simply a magnification of the object caused by moving closer to it.

As the temperature difference $ \Delta T = T w \u2212 T a$ decreases (while keeping the water temperature fixed), the gradient of the refractive index decreases as well, thus rays bend less and less. In Fig. 5, we can see this effect for five different temperature difference values (1, 2, 3, 4, and $ 5 \xb0$ C) with fixed water temperature $ 5 \xb0$ C. Rays start from 2.7 m above ground and travel a horizontal distance of 2000 m; however, since the temperature gradient is confined to a region less than 0.5 m above the surface, Fig. 5 is zoomed in to show where the paths differ, while the full path is shown in the inset. We can see that rays with a smaller temperature difference ray-traced from the observer hit the ground (solid black line) rather than reflecting, so they let the observer see the surface itself. Also, since the bending of rays with a smaller temperature difference is smaller as well, these can get closer to the ground, thus making the bottom of objects more visible (less cutoff).

Note that the angular size of the mirage is relatively small compared to the object: it is only a little above a quarter of a degree, as specified in the range of initial angles in Figs. 2–5. The relative size of the vertical and horizontal axes is not realistic: the horizontal axis was shortened for a better visibility. Also in Fig. 1, we are zooming in to show the mirage, because otherwise the effect would be too tiny to be properly presented.

## III. THE TEMPERATURE PROFILE

Since the temperature gradient is the driving force of the mirage generation, it is required to make a realistic model for the temperature as a function of height above the surface. Such measurement-based models can be formulated using high precision measurements of the air temperature, especially near the surface. Here, we provide an exponential temperature function based on fits to measurements from Ref. 11. In this work, both water and air temperature were sampled with a 1500-point time series close to the surface, within a 2 m height region. Note that other types of fit functions can also be used, see Refs. 12 and 13.

*T*is temperature as a function of height

*h*and

*T*is the ambient temperature to which the exponential function relaxes. $ \Delta T$ gives the change in temperature between water and ambient air, while

_{a}*H*is the characteristic distance, over which these variations become relevant. Obtaining the general formula for the air temperature required the analysis of data from Ref. 11. Our method is explained in detail in Appendix A.

*n*of the air, using the formula given in Ref. 14

*p*is given in kPa and air temperature

*T*is given in Eq. (1) in units of ° C. A term with relative humidity could be added, but it is negligible in our case, since its effect is less than 0.5% of that of the corresponding temperature and pressure. Latter was taken to be

*p*= 101 kPa throughout the simulations, deviations from this value were not taken into account.

The dependence of the refractive index on the height above the water surface is presented in Fig. 6 with fixed ambient temperature values. It can be clearly seen that most of the variation in the refractive index happens on the first 10–20 cm above ground; thus, the relevant bending of light traces occurs here as well. At higher elevations, we assume saturation to the ambient temperature. It is interesting to note that the source of the whole mirage effect is almost exclusively in this 10–20 cm thick region right above the surface.

## IV. RAY TRACING

^{15–18}which can be derived from the wave equation in the short wavelength limit. It describes the phase $\phi $ of the wave function, which is proportional to the propagation time of a wavefront

*v*is the speed of light in the medium given by

*c*/

*n*with

*n*being the refractive index. From this partial differential equation, one can write

^{19}the equations of motion as ordinary differential equations of the form

**r**is the position vector,

*s*is the length along the path, and

**u**is the slowness vector defined by $ ( 1 / v ( r ) ) \u2009 ( d \u2009 r / d \u2009 s )$. For better understanding, see the definition of variables and vectors in Fig. 7.

These equations give the d **r** difference of the position vector given an infinitesimal movement along the path d*s*. The input of the equations is the value of the refractive index at all points in space, which can be used to calculate speed, *v* and its gradient. From here, we have the difference of the slowness vector d**u**/d*s*, which we can use to calculate **u**. Now we know every term on the right hand side of the first equation, so we can get d**r**/d*s*. Integrating Eq. (4), we can calculate rays. The integration is performed using the adaptive Runge–Kutta–Fehlberg method of the GNU Scientific Library^{20} (GSL).

Figure 8 shows the structure of the program. After initialization, the code gets an input, which is on the one hand: a black and white PNG picture as the object, placed at fixed distance from the observer. On the other hand, one has to set the water temperature, *T _{w}* and the $ \Delta T$ inputs as well. The picture should be a scenery without a mirage: A photo of a distant shore, for example. Rays travel towards this object from the observer, who is represented by a pinhole camera. The output picture is created on the backplane of the camera, the optical screen, which is then turned to the normal view in this work. From every “pixel” on the screen, rays are cast towards the pinhole of the camera, which sets the initial direction of each ray. Then, rays are traced through the medium by numerically solving the ODEs derived from the eikonal equation in the main OdeSolverGsl part. Note that this figure shows the main C++ objects and the data flow between them. That's why there are distinct arrows marked by

*c*and

*d*, although the data carried by them are the same.

When rays reach the plane of the object PNG picture, it is evaluated if they hit it or not. If no, the pixel corresponding to the given ray in the back of the camera will be assigned black. If a ray hits the object, it will be assigned the intensity (color) of the object at that point. In the code, more sophisticated averaging methods are used for the determination of the intensity of each pixel. Unfortunately, we could not use GSL directly, because the independent variable cannot connect directly to object positions in space. To solve this problem, we have implemented an indirect method similar to binary search on top of GSL's integration solution. For more details, please refer to Fig. 15 and the algorithm description in Appendix B. Nevertheless, the output picture is built up pixel by pixel by sending out rays from every single one of them, tracing their path, and then checking what they hit.

Our program can be downloaded from the supplementary material online and can be freely used for educational purposes. A more detailed description of the simulation code is added in Appendix B.

## V. NUMERICAL RESULTS

Using the program described above, it is possible to reproduce real life mirages. In the top of Fig. 1, we can see a picture taken at a distance of about 17 km at Lake Balaton, Hungary. The viewpoint was around 2.7 m above water, with an ambient temperature of about $ 1 \xb0$ C and a water temperature equalling $ 5 \xb0$ C. In the bottom part of Fig. 1, one can see a picture that was taken 50 m above ground. The height is needed so that lower parts of the scenery would also be visible. Otherwise (even with no mirage) the opposite shore and small buildings near the shore could not be seen because of Earth's curvature. Thus, the bottom picture was used as a base for simulations, placed at sea level (its lower edge is cut at the shore). The scenery is considered two-dimensional; 3D variations (depth) were not taken into account. Black and white pictures were used for simplicity. This model does not take into account the negligible effect of dispersion: the difference between bending of light rays with different colors (wavelengths). Here, we used monochromatic assumption; therefore, the wavelength is fixed; so each ray, regardless of energy, bends the same way. Contrast in the photo was increased to make details more distinctive.

In Fig. 9, we can see simulated mirages, where the temperature of water was taken to be $ 5 \xb0$ C, and the ambient changed with a $ \Delta T$ difference of 0, 2, 4, and $ 5 \xb0$ C. Water shown underneath is not part of the simulation, and it is located between the object and the observer. The top image is a reference (indicated by the blue dashed line in all the figures), meaning that the ambient temperature equals that of water: $ \Delta T = 0 \xb0$ C. Here, light rays do not get bent; however, we still cannot see the very shore at *h* = 0 because of the curvature of Earth, only higher structures. Basically, the blue line indicates the boundary between water and land in the absence of a mirage, so it is *exactly* at the same location in all the pictures. From the other three images, it is clear that the greater the difference in temperature, the “taller” (extends further towards the bottom) the “mirror” image gets. Basically, the blue arrow gets longer. (This is equivalent to a greater temperature gradient, thus also a greater refractive index gradient.) Red arrows on the left correspond to the part of the original, upright image, that gets “mirrored” upside down (blue arrows). They have the same meaning as in Fig. 2. The “mirror” axis is shown with a red dashed line (same as the red and blue dashed line in Fig. 2). One can see that for the cases with a temperature gradient, it does not coincide with the blue one, and it does change, although to a quite moderate extent. A narrow band at the bottom of the scene gets cut out solely by the mirage effect—not due to the curvature of Earth.

Figure 10 was made using the same setup as Fig. 9, except for here we used the picture of a chess board^{25} as the base of the simulations for a better understanding. Also, temperature differences are shown for every $ \xb0$ C between 0 and 5, where the last one is the reference again (indicated by the blue dashed line in every figure). One can see that both the upright and the upside down images get cropped by the mirage (red dashed “mirror horizon” is higher up than the dashed blue reference line); however, the bottom of the “mirrored” one is much more compressed. For example, where $ T a = 0 \xb0$ C, row 4 is clearly less compressed, than row 3 or 2, and row 1 is barely visible. As the ambient temperature rises, this relative compression seems to stay the same, but the upside down image shrinks gradually, until it eventually disappears completely in the absence of a temperature gradient.

We created two videos, which show how the scenery at Lake Balaton and the view of the chess board changes during a linear rise in ambient temperature ( $ \Delta T$ decreases with time), as the air is getting heated up in the morning.

In Fig. 11, we can see variations due to viewing from different distances. In Fig. 11, for every image $ T a = 1 \xb0$ C and $ T w = 5 \xb0$ C. Size (height) of the base image was adjusted to the distance with a simple proportionality. Note that there is also an autozoom applied to fill the screen with the images. We can see that the dependence on distance is not a straightforward one: The curvature of Earth and the mirage both play a role. Former is represented by the blue dashed (reference) line, which would be the bottom at *T _{a}* =

*T*, skyline, the boundary between surface and the object. As the distance decreases, the upside down images get less compressed, and they slightly shrink in vertical size. The greatest change in size happens on the last few kilometers. It is interesting to note that for the last image (the 1 km one), the “mirror horizon” (red dashed line) is located underneath the blue one, unlike for the other cases. This suggests that between the upright image and the “mirrored” one, we would see objects from the foreground of the scenery.

_{w}## VI. CONCLUSIONS

For a better understanding of mirages in nature, we developed a program code that can be used for educational purposes.

As it can be seen in Fig. 9, the most apparent feature that changes as a function of air temperature is the vertical size of the mirage. We compared the ratio *size of tree*/*size of mirage* in the simulated and the photographed (Fig. 1) case, thus determined an air temperature *T _{a}* at the time of taking the photo. Basically, we checked which picture matches the real one the best from an ensemble of 50 simulated images with different $ \Delta T$s. We got a value of $ T a = 4.3 \xb0$ C (with a range from $ 3.85 \xb0$ C to $ 4.65 \xb0$ C), with a temperature difference between air and water equalling $ \Delta T = 0.7 \xb0$ C (with a range from $ 0.35 \xb0$ C to $ 1.15 \xb0$ C). This is based on the assumption that the temperature of water is $ 5 \xb0$ C, measured exactly. (We did not take into account its uncertainty.) Indicated errors correspond to distance measurements (for determining ratios) in the pictures. The temperature of air at the time was in fact $ 1 \xb0$ C, which falls outside of the margin of error. However, this value was measured only at the observer's point and most probably varied a significant amount along the 17 km path, where the light rays had to travel. Finally, our simple model of the temperature distribution does not take into account horizontal changes as well, which would also be difficult to determine experimentally. The temperature of water (although more stable over larger distances than that of air) of course also varied.

Understanding natural phenomena is easier with visualization and with the possibility of changing the visualization interactively. Some physical setups would be difficult—or even impossible—to study based solely on calculations. Our program serves the purpose of being a tool for those, who would like to have a better understanding of the phenomenon of *fata morgana*.

## ACKNOWLEDGMENTS

The research was supported by the Hungarian National Research, Development and Innovation Office (NKFIH) under Contract Nos. OTKA K123815, K135515, and 2020-2.1.1-ED-2021-00179. Computation resources were provided by the Wigner Scientific Computational Laboratory (WSCLAB) (the former Wigner GPU Laboratory). A.H. acknowledges the discussion with Dániel Berényi.

## AUTHOR DECLARATIONS

### Conflict of Interest

The author has no conflicts to disclose.

### APPENDIX A: REFRACTIVE INDEX OF AIR

The refractive index of air can be calculated via Eqs. (1) and (2). These are phenomenological inputs; therefore it is required to measure and fit the model experimentally. In this work, we analysed the data from Ref. 11 to obtain the air temperature, which was used in Eq. (2) providing the refractive index. Here, we present the data analysis for the air temperature first, and then how the refractive index can be given.

Since relatively more data points represent the temperature of ambient air *T _{a}*, and the gradient here is small, the fitting of this parameter is more accurate than that of $ \Delta T$ and

*H*. However, temperatures higher above the surface can vary more due to air turbulence and other effects.

^{3}Here, measurements were taken utilizing Raman scattering in an optical cable and might have been subject to unwanted heating due to direct sunshine. In this case, altogether 12 measurements were taken at different times of the day at the surface of Lake Geneva. For some of these, the ambient temperature was lower; for others, it was higher than that of the water.

^{22}

The temperature of water in the data was about $ 10 \xb0$ C, and the ambient temperature varied approximately between 8°C and 15°C. Fortunately, typical temperature profiles over heated surfaces have similar characteristics over a wide range of temperatures and circumstances; therefore, normalized fits can be done.^{12,23,24} Temperature changes most in the first, approximately 10–20 cm from ground, and so does the refractive index, since the two have a connection close to linear. The data in Fig. 12 were created by fitting exponential functions to each of the 12 measurements. Then we normalized them by subtracting the ambient temperature *T _{a}* and dividing them by $ \Delta T$—both fit parameters. The final fit parameters for the collected data can be seen in the box in Fig. 12.

Measurements were taken at different *T _{a}* and

*T*values, the obtained parameters vary. Observing scaling of the fit parameters, normalization has been applied by shifting the temperature values with

_{w}*T*and $ \Delta T norm$ for each dataset. With this method, one single general form can be obtained for

_{anorm}*T*at a given height,

*h*; in agreement with any

*T*and

_{a}*T*values. The normalized data can also be seen in Fig. 12.

_{w}*T*and

_{a}*T*, respectively, the temperature profiles can be calculated by multiplying the fit function by $ T w \u2212 T a$ then adding

_{w}*T*,

_{a}Next, we plotted the calculated refractive index *n* in Fig. 14 by substituting Eq. (1) into Eq. (2). The result is a monotonically decreasing function of the air temperature (*T*), and the relationship is approximately linear on a wide range. Here, we plotted various refractive index profiles as a function of the ambient temperature at fixed heights and water temperature based on Eq. (A1).

### APPENDIX B: THE STRUCTURE OF THE PROGRAM

Here, we provide the structure of the C++ simulation program code,^{25} which can be downloaded and freely used for educational purposes. The code traces optical rays in a medium with a nonuniform refractive index. Using a landscape photo as input, it generates a new picture with a mirage. Although the GitLab page contains instructions for compiling, it may not be straightforward, especially for other Linux distributions than Ubuntu. Feel free to contact the authors.

The structure of the program is shown in Fig. 8. Simulated rays start from a pixel on the optical screen of a pinhole camera. The initial angle of the rays is determined by the relative position of each pixel and the hole. To let the program render better-looking pictures (actually reduce the aliasing effect of the program code near high contrast edges), pixel subsampling is possible. This means each camera pixel can be divided into *k* × *k* cells with each cell casting a ray. The pixel value will be the average of intensities (colors) these rays hit.

After the rays leave the camera (a), they travel through the medium with a predetermined temperature distribution and pressure. Each one of them is traced until it reaches the water surface or the plane of the target (object) (b). Here, if it hits the object, it will be assigned an intensity characteristic of the target at that point (c). If it does not hit anything, it will be assigned black. This intensity will be subject to the subpixel average calculation, thus forming the image that the camera sees (d).

For input and output picture management, we use *png++*, which is a C++ wrapper around *libpng*. We have chosen the PNG format, because it is efficient, widely used, and allows us to examine high-contrast pictures without JPEG compression artifacts.

For low-level differential equation manipulation, we use the *Odeiv2* module from GNU Scientific Library (GSL),^{20} a proven and robust mathematical library written in C. *Odeiv2* contains a set of low-level to high-level methods to numerically solve *n*-dimensional first-order differential equations (ODE). To describe Eq. (4) in *Odeiv2*, we use a system of 6 ODEs, one for each axis of each vectorial differential equation.

As *Odeiv2* can only solve the system of ODEs at specific values of the independent variable (here the track travelled along the ray), the *Odeiv2* library cannot be used directly to solve the system of ODEs for an intersection with some geometry (like the object's plane). This is probably a common problem, so we have developed a C++ class around *Odeiv2* capable of solving general problems similar to this.

From an algorithmic point of view, the *OdeSolverGsl* class is the most important among the ones drawn in Fig. 8. It is a general solution, a template in C++ terminology. It is templated by the class describing the actual physical problem, the eikonal. This class is capable of describing flat and round Earth physics. It is parameterized by the Earth form, radius (for round Earth), surface model, and ambient and base temperature. There is an important simplification for the round Earth case: The object plane is not perpendicular to the surface at its location, but to that at the camera's position. However, since the average visibility distances are not in the magnitude of Earth's radius, its effect cannot be noticed.

The *OdeSolverGsl* class receives a configured physics description class (here *Eikonal*) as template parameter. Further parameters are the *Odeiv2* stepper method, integration limits (here distance along the ray), and relative and absolute tolerances to meet during approximation and stepping parameters. These are required to control the behaviour of the adaptive Runge–Kutta stepper.

Figure 15 shows how the algorithm finds the target along a ray. This figure depicts four iterations on the same ray to perform the integration until a planar target (dashed line). To check where the intersection of the ray and the object (here plane) is, we can only inspect sample points calculated along the ray. We use consecutive steps calculated by the adaptive Runge–Kutta method. Since an adaptive Runge–Kutta method uses larger-and-larger steps if the integrated function does not vary rapidly, the object most probably gets stepped over. Every time this happens, the integration is restarted (using the tiny initial step size) from the last sampled point before the target, until it gets into a predefined vicinity of the target. In Fig. 15, the last step point before the target on the previous iteration becomes the start point for stepping on the next iteration. As the adaptive step size grows approximately exponentially, the cost of this solution in term of steps is similar to binary search, so around $ O ( \u2212 log \u2009 x \u2009 log \u2009 \u03f5 )$ where *x* is the distance along the *x*-axis, which the ray should travel and *ϵ* is the required accuracy to approximate the *x* coordinate of the object.

The above algorithm is only the core idea of the solution in *OdeSolverGsl—*it needed several enhancements. First, steps must be checked to see if the next calculated ray point is below the surface, which would be nonphysical, of course. If it is, the *Eikonal* ODE calculation returns an *Odeiv2* error to indicate that the step was too big. It will return with a half step, since it is possible that the ray does not truly hit the ground: the adaptive solver just got inaccurate due to the growing step size.

If the step size gets below a configurable minimal value due to the above retries, the *OdeSolverGsl* concludes that the direction is definitely wrong, because the ray will hit the surface (cannot pass over it) before reaching the target. In this case, it returns indicating that the result is invalid.

If the step is larger than a configurable maximum value, and the step along the ray changes direction to a greater extent than a configured value, we must restart the steps from the initial small value. This is necessary because the *Odeiv2* step control algorithm sometimes fails to detect the inaccuracy against the given tolerance requirements.

This check and the hit check against the target geometry are performed by *λ* functions, so *OdeSolverGsl* remains general purpose, and can be reused for other problems as well.

## REFERENCES

*Cambridge Dublin Math. J.*

*The Scientific Papers of James Clerk Maxwell*

*Mathematical Theory of Optics*

*Propagation of Electromagnetic Waves*

*Chess Openings For Dummies*

*First Thermal and Fluids Engineering Summer Conference*