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.

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 superior, 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 surfaces3,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.

Fig. 1.

(Color online) Inferior mirage over Lake Balaton, Hungary taken from a height of 2.7 m (top) and the scenery of the same area without a mirage (bottom) taken from a height of approximately 50 m, so that the surface level would not get cut off by Earth's curvature.

Fig. 1.

(Color online) Inferior mirage over Lake Balaton, Hungary taken from a height of 2.7 m (top) and the scenery of the same area without a mirage (bottom) taken from a height of approximately 50 m, so that the surface level would not get cut off by Earth's curvature.

Close modal

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 ° C and ambient air temperature 1 ° 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.

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 ° C, the water temperature is T w = 5 ° 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.

Fig. 2.

(Color online) Bending of light rays in air over a water surface with ambient temperature 1 ° C and water temperature 5 ° C, over a distance of 2000 m. The rays start at the observer at 2.7 m at angles ranging from 0.16 ° below horizontal to above horizontal + 0.04 ° in 0.02 ° steps. The curved surface of Earth is indicated by a solid black line that starts at h = 0 m, the object is indicated by a red arrow, and its appearances are shown by blue ones. We indicated the “mirrored” part of the object by blue on the object. The red-and-blue dashed line corresponds to the reference line (same as the blue dashed lines in Figs. 9–11).

Fig. 2.

(Color online) Bending of light rays in air over a water surface with ambient temperature 1 ° C and water temperature 5 ° C, over a distance of 2000 m. The rays start at the observer at 2.7 m at angles ranging from 0.16 ° below horizontal to above horizontal + 0.04 ° in 0.02 ° steps. The curved surface of Earth is indicated by a solid black line that starts at h = 0 m, the object is indicated by a red arrow, and its appearances are shown by blue ones. We indicated the “mirrored” part of the object by blue on the object. The red-and-blue dashed line corresponds to the reference line (same as the blue dashed lines in Figs. 9–11).

Close modal

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.

Fig. 3.

(Color online) Bending of light rays in air over a surface of water with ambient temperature 1 ° C and water temperature 5 ° C, on a distance of 2000 m. Rays start from heights 2.0, 2.7, and 3.5 m. Initial angles from horizontal are: 0.16 ° , 0.12 ° , 0.08 ° , 0.04 ° , 0 °, and 0.04 o. The surface of Earth is indicated with a solid black curve, and the arrow represents the object.

Fig. 3.

(Color online) Bending of light rays in air over a surface of water with ambient temperature 1 ° C and water temperature 5 ° C, on a distance of 2000 m. Rays start from heights 2.0, 2.7, and 3.5 m. Initial angles from horizontal are: 0.16 ° , 0.12 ° , 0.08 ° , 0.04 ° , 0 °, and 0.04 o. The surface of Earth is indicated with a solid black curve, and the arrow represents the object.

Close modal

Thick lines (with initial angle 0.16 °) 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 ° 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.

Fig. 4.

(Color online) Bending of light rays in air over a surface of water with ambient temperature 1 ° C and water temperature 5 ° C, from 2.7 m above the ground. Rays start from distances 1500, 2000, and 2500 m from the object. Initial angles from horizontal are: 0.16 ° , 0.12 ° , 0.08 ° , 0.04 ° C, 0 °, and 0.04 °. The surface of Earth is indicated with a solid black curve, and the arrow represents the object.

Fig. 4.

(Color online) Bending of light rays in air over a surface of water with ambient temperature 1 ° C and water temperature 5 ° C, from 2.7 m above the ground. Rays start from distances 1500, 2000, and 2500 m from the object. Initial angles from horizontal are: 0.16 ° , 0.12 ° , 0.08 ° , 0.04 ° C, 0 °, and 0.04 °. The surface of Earth is indicated with a solid black curve, and the arrow represents the object.

Close modal

As the temperature difference Δ T = T w 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 ° C) with fixed water temperature 5 ° 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).

Fig. 5.

(Color online) Bending of light rays in air over a surface of water with temperature differences 5, 4, 3, 2, and 1 ° C and the water temperature is 5 ° C. Rays start from 2.7 m above ground with a distance of 2000 m from the object. Initial angles from horizontal are: 0.16 ° , 0.12 ° , 0.08 ° , 0.04 ° C, 0 °, and 0.04 °. The figure is zoomed in to show where the paths differ, while the full path is shown in the inset. The surface of Earth is indicated with a black curve. Note that certain rays hit the ground, and naturally they do not continue after that.

Fig. 5.

(Color online) Bending of light rays in air over a surface of water with temperature differences 5, 4, 3, 2, and 1 ° C and the water temperature is 5 ° C. Rays start from 2.7 m above ground with a distance of 2000 m from the object. Initial angles from horizontal are: 0.16 ° , 0.12 ° , 0.08 ° , 0.04 ° C, 0 °, and 0.04 °. The figure is zoomed in to show where the paths differ, while the full path is shown in the inset. The surface of Earth is indicated with a black curve. Note that certain rays hit the ground, and naturally they do not continue after that.

Close modal

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.

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.

As a first step, we created a generalized parametrization for all the measured temperature profiles above water bodies for the simulations. This type of measured data should be fitted with an exponential function of the form
T = T a + Δ T e h / H ,
(1)
where T is temperature as a function of height h and Ta is the ambient temperature to which the exponential function relaxes. Δ T gives the change in temperature between water and ambient air, while 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.
In the second step, we calculated the refractive index n of the air, using the formula given in Ref. 14 
n = 1 + 7.86 × 10 4 ° C kPa × p 273.15 ° C + T ,
(2)
where pressure 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.

Fig. 6.

(Color online) Refractive index of air as a function of height at different ambient temperatures Ta written on the curves with error bands. The water temperature is T w = 5.0 ° C.

Fig. 6.

(Color online) Refractive index of air as a function of height at different ambient temperatures Ta written on the curves with error bands. The water temperature is T w = 5.0 ° C.

Close modal
Since we are dealing with a problem of refraction, we could use Snell's law to calculate ray paths. However, in a continuous medium, it is more natural to use the eikonal equation,15–18 which can be derived from the wave equation in the short wavelength limit. It describes the phase φ of the wave function, which is proportional to the propagation time of a wavefront
| φ | 2 = 1 v 2 ,
(3)
where 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 write19 the equations of motion as ordinary differential equations of the form
d r d s = v ( r ) u ( r ) , d u d s = 1 v ( r ) 2 v ( r ) ,
(4)
where r is the position vector, s is the length along the path, and u is the slowness vector defined by ( 1 / v ( r ) ) ( d r / d s ). For better understanding, see the definition of variables and vectors in Fig. 7.
Fig. 7.

(Color online) Variables of the eikonal equation.

Fig. 7.

(Color online) Variables of the eikonal equation.

Close modal

These equations give the d r difference of the position vector given an infinitesimal movement along the path ds. 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 du/ds, 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 dr/ds. Integrating Eq. (4), we can calculate rays. The integration is performed using the adaptive Runge–Kutta–Fehlberg method of the GNU Scientific Library20 (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, Tw and the Δ 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.

Fig. 8.

(a) Starting point and direction. (b) Coordinates of the ray's point of arrival. (c) and (d) Intensity (color) of object at the given point (black, if it is not on the object). The important C++ code objects are in italic.

Fig. 8.

(a) Starting point and direction. (b) Coordinates of the ray's point of arrival. (c) and (d) Intensity (color) of object at the given point (black, if it is not on the object). The important C++ code objects are in italic.

Close modal

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.

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 ° C and a water temperature equalling 5 ° 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 ° C, and the ambient changed with a Δ T difference of 0, 2, 4, and 5 ° 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: Δ T = 0 ° 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.

Fig. 9.

(Color online) Simulated mirages over Lake Balaton with T w = 5 ° C and different ambient temperatures ( Δ T = T w T a). Water shown underneath is not part of the simulation, and it is located between the object and the observer. Blue dashed line indicates the reference, red (same as the red and blue dashed line in Fig. 2) the “mirror” axis. Red arrow corresponds to that part of the upright image, that is “mirrored” upside down (blue arrow). They have the same meaning as in Fig. 2. Multimedia available online.

Fig. 9.

(Color online) Simulated mirages over Lake Balaton with T w = 5 ° C and different ambient temperatures ( Δ T = T w T a). Water shown underneath is not part of the simulation, and it is located between the object and the observer. Blue dashed line indicates the reference, red (same as the red and blue dashed line in Fig. 2) the “mirror” axis. Red arrow corresponds to that part of the upright image, that is “mirrored” upside down (blue arrow). They have the same meaning as in Fig. 2. Multimedia available online.

Close modal

Figure 10 was made using the same setup as Fig. 9, except for here we used the picture of a chess board25 as the base of the simulations for a better understanding. Also, temperature differences are shown for every ° 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 ° 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.

Fig. 10.

(Color online) Simulated mirages of a chess board with T w = 5 ° C and different ambient temperatures, such that Δ T = T w T a. The uniformly gray area shown underneath is not part of the simulation, and it is located between the object and the observer. Blue dashed line indicates the reference (same as the red-and-blue dashed line in Fig. 2), and red dashed line is the “mirror horizon.” Multimedia available online.

Fig. 10.

(Color online) Simulated mirages of a chess board with T w = 5 ° C and different ambient temperatures, such that Δ T = T w T a. The uniformly gray area shown underneath is not part of the simulation, and it is located between the object and the observer. Blue dashed line indicates the reference (same as the red-and-blue dashed line in Fig. 2), and red dashed line is the “mirror horizon.” Multimedia available online.

Close modal

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 ( Δ 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 ° C and T w = 5 ° 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 Ta = Tw, 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.

Fig. 11.

(Color online) Simulated mirages of a chess board viewed from different distances x. Height h of the base picture was adjusted to the distance. The uniformly gray area shown underneath is not part of the simulation, and it is located between the object and the observer. Blue dashed line (aligned for all cases) indicates the reference (where Ta = Tw, not shown here), and the red dashed one is the “mirror horizon” line.

Fig. 11.

(Color online) Simulated mirages of a chess board viewed from different distances x. Height h of the base picture was adjusted to the distance. The uniformly gray area shown underneath is not part of the simulation, and it is located between the object and the observer. Blue dashed line (aligned for all cases) indicates the reference (where Ta = Tw, not shown here), and the red dashed one is the “mirror horizon” line.

Close modal

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 Ta 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 Δ Ts. We got a value of T a = 4.3 ° C (with a range from 3.85 ° C to 4.65 ° C), with a temperature difference between air and water equalling Δ T = 0.7 ° C (with a range from 0.35 ° C to 1.15 ° C). This is based on the assumption that the temperature of water is 5 ° 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 ° 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.

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.

The author has no conflicts to disclose.

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 Ta, and the gradient here is small, the fitting of this parameter is more accurate than that of Δ 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 ° 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 Ta and dividing them by Δ T—both fit parameters. The final fit parameters for the collected data can be seen in the box in Fig. 12.

Fig. 12.

(Color online) Exponential fit (solid line) of the temperature as a function of height on normalized measurement data with error bands (dashed lines). Equation (A1) shows how we calculated actual temperature profiles using these fit parameters.

Fig. 12.

(Color online) Exponential fit (solid line) of the temperature as a function of height on normalized measurement data with error bands (dashed lines). Equation (A1) shows how we calculated actual temperature profiles using these fit parameters.

Close modal

Measurements were taken at different Ta and Tw values, the obtained parameters vary. Observing scaling of the fit parameters, normalization has been applied by shifting the temperature values with Tanorm and Δ T norm for each dataset. With this method, one single general form can be obtained for T at a given height, h; in agreement with any Ta and Tw values. The normalized data can also be seen in Fig. 12.

Knowing ambient temperature and the temperature of water, Ta and Tw, respectively, the temperature profiles can be calculated by multiplying the fit function by T w T a then adding Ta,
T = T a + ( T anorm + Δ T norm e h / H ) ( T w T a ) .
(A1)
A few examples have been calculated and presented in Fig. 13. Note that temperatures at the surface can get higher than that of water. This is due to uncertainties of the measurement and fit parameters.
Fig. 13.

(Color online) Modelled temperature of the air as a function of height above water with error bands for given ambient temperature values Ta written on the curves. The water temperature is T w = 5.0 ° C.

Fig. 13.

(Color online) Modelled temperature of the air as a function of height above water with error bands for given ambient temperature values Ta written on the curves. The water temperature is T w = 5.0 ° C.

Close modal

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).

Fig. 14.

(Color online) Refractive index of air as a function of the ambient temperature at different heights with error bands.

Fig. 14.

(Color online) Refractive index of air as a function of the ambient temperature at different heights with error bands.

Close modal

The obtained curves in Figs. 13 and 14 can be used to calculate the refractive index, n, using the program code setting the input parameter values for Ta and Tw. The n as a function of h is finally given in Fig. 6 in Sec. III.

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 ( log x log ϵ ) 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.

Fig. 15.

Iterations 1–4 for finding solution. s: Starting point and direction of iteration. a–e: Steps of adaptive Runge–Kutta. The one sufficiently close to the target condition (in our case x coordinate of the object) is the solution. If it goes past the target, in the next iteration, it restarts from the last point before the target.

Fig. 15.

Iterations 1–4 for finding solution. s: Starting point and direction of iteration. a–e: Steps of adaptive Runge–Kutta. The one sufficiently close to the target condition (in our case x coordinate of the object) is the solution. If it goes past the target, in the next iteration, it restarts from the last point before the target.

Close modal

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.

1.
J.
Blanco-García
,
V.
Dorrío
, and
F.
Ribas
, “
Photographing mirages above the sea surface
,” in
proceedings of 8th International Conference on Hands-on Science
, Ljubljana, Eslovenia (
2011
).
2.
Robert G.
Greenler
, “
Laboratory simulation of inferior and superior mirages
,”
J. Opt. Soc. Am. A
4
(
3
),
589
590
(
1987
).
3.
Nabil
Wakid
, “
Temperature profile and double images in the inferior mirage
,” e-print arXiv:1911.03507 (
2019
).
4.
Alistair B.
Fraser
and
William H.
Mach
, “
Mirages
,”
Sci. Am.
234
(
1
),
102
111
(
1976
).
5.
P. R.
Barker
,
P. R. M.
Crofts
, and
M.
Gal
, “
A superior ‘superior’ mirage
,”
Am. J. Phys.
57
(
10
),
953
954
(
1989
).
6.
Waldemar H.
Lehn
, “
Exact temperature profile for the hillingar mirage
,”
Am. J. Phys.
69
(
5
),
598
600
(
2001
).
7.
Himanshu
Kumar
,
Sumana
Gupta
, and
K. S.
Venkatesh
, “
A novel method for inferior mirage detection in video
,” in
proceedings of conference on Digital Image & Signal Processing (DISP)
(
ACM
,
2019
).
8.
Mario
Branca
, “
Simulation of the inferior mirage
,”
Phys. Teach.
48
(
6
),
372
373
(
2010
).
9.
T.
Kosa
and
P.
Palffy-Muhoray
, “
Mirage mirror on the wall
,”
Am. J. Phys.
68
(
12
),
1120
1122
(
2000
).
10.
E.
Fabri
,
G.
Fiorio
,
F.
Lazzeri
, and
P.
Violino
, “
Mirage in the laboratory
,”
Am. J. Phys.
50
(
6
),
517
520
(
1982
).
11.
J. S.
Selker
,
L.
Thevenaz
,
H.
Huwald
,
A.
Mallet
,
W.
Luxemburg
,
N.
van de Giesen
,
M.
Stejskal
,
J.
Zeman
,
M.
Westhoff
, and
M. B.
Parlange
, “
Distributed fiber-optic temperature sensing for hydrologic systems
,”
Water Resour. Res.
42
(
12
),
W12202
, https://doi.org/10.1029/2006WR005326 (
2006
).
12.
Song
Yang
,
Jing-Yuan
Zhang
,
Ying-Ying
Yang
,
Jun-Yuan
Huang
,
Yun-Rui
Bai
,
Yong
Zhang
, and
Xue-Chun
Lin
, “
Automatic compensation of thermal drift of laser beam through thermal balancing based on different linear expansions of metals
,”
Results Phys.
13
,
102201
(
2019
).
13.
J. F.
Croft
, “
The convective regime and temperature distribution above a horizontal heated surface
,”
Q. J. R. Met. Soc.
84
(
362
),
418
427
(
1958
).
14.
J.
Stone
and
J.
Zimmerman
, “
Index of refraction of air
” (
2001
), see https://www.nist.gov/publications/index-refraction-air.
15.
Qi
Mo
,
Hengchin
Yeh
, and
Dinesh
Manocha
, “Tracing analytic ray curves for light and sound propagation in non-linear media,”
IEEE Trans. Visual. Comput. Graph.
22
(11),
2493
2506
(
2014
).
16.
James Clerk
Maxwell
, “
Solutions of problems: (prob. 3, vol. VIII. p. 188)
,”
Cambridge Dublin Math. J.
9
,
9
11
. (
1854
);
reprinted by
The Scientific Papers of James Clerk Maxwell
, edited by
William Davidson
Nivin
(
Dover Publications
,
New York
,
1890
), pp.
76
79
.
17.
Rudolf Karl Luneberg
,
Mathematical Theory of Optics
(
University of California Press
,
Berkeley
,
1964
), p.
401
.
18.
Rudolf Karl Luneburg
,
Propagation of Electromagnetic Waves
(
New York University
,
New York
,
1948
).
19.
A.
Rüger
, “
Dynamic ray tracing and its application in triangulated media
,”
Dissertation
,
ProQuest
,
Ann Arbor
,
1993
.
20.
Free Software Foundation, Inc.
, “
GNU scientific library: Ordinary differential equations
,” <https://www.gnu.org/software/gsl/doc/html/ode-initval.html>, last accessed on April 22,
2022
.
21.
By the permission of the owner (John Wiley & Sons, Inc.), <https://www.dummies.com/article/home-auto-hobbies/games/board-games/chess/understanding-chess-notation-192295/>, last accessed on May 17, 2022;
James
Eade
,
Chess Openings For Dummies
(
Wiley
,
New York NY
,
2010
).
22.
We focused on inferior mirages; however, from this, data superior mirages could also be simulated, where the temperature of ambient air is higher than the temperature of water, which usually happens later in the day. For the latter case, the visibility is usually worse, and the direct reflection of sunlight from the surface could distinguish the mirage effect.
23.
David
Duncan
, “
Influence of pavement type on near surface air temperature
,” MS thesis,
Clemson University
,
2011
.
24.
Peter
Disimile
and
Norman
Toy
, “
The influence of engine pool fire dynamics due to external sources
,” in
First Thermal and Fluids Engineering Summer Conference
, August 9–12, 2015,
New York City, NY
(
ASTFE
,
2015
).
25.
See supplementary material online for the code for simulating mirages.

Supplementary Material