Fresnel-Gaussian shape invariant for optical ray tracing

Moisés Cywiak; A. Morales; J. M. Flores; Manuel Servín

doi:10.1364/OE.17.010564

1. Introduction

The general problem of ray tracing, consisting of a series of refracting or reflecting surfaces, and/or homogeneous/inhomogeneous media (or a combination of them), can be determined by repeated application of common ray-tracing equations [1-4]. When light is strongly focussed well-known caustic theories are applied [5].

Evolved optical devices, as graded-index media, used for example in the design of inhomogeneous lenses, require new methods for ray tracing [6]. Moreover, today’s technological demands, as for example, the use of ultra-short laser pulses strongly focussed, require more versatile techniques [7, 8]. On these conditions, common ray tracing programs need to evolve to accurately describe the complete behavior of light through new optical devices and complex optical media. Examples on this issue can be found in Ref. [9], where two commercially available programs, one for ray tracing and another that includes diffraction effects were combined for ray tracing of ultra-short focussed beams through optical devices.

Common techniques for ray tracing relay in only one driving term that directs the ray properly, calculating it at each point of the optical path and then, updating this angle at each point of interest using mainly Snell’s law. More elaborate techniques has been proposed, based on finding the parameters for reflected and transmitted Gaussian beams (the waist, the center, the direction of propagation and the complex angle of rotation), using phase matching by approximating the refractive or reflective surface with a quadratic function [10].

The tracing technique that is proposed here provides more versatility from a computational point of view, in the sense that it presents two driving terms which give the direction of the ray at each point. Each of the two driving terms can be used independently or can be combined. The first driving term is related to the angle of propagation, but differently to common techniques, this term is included in the phase of our proposed Gaussian invariant. The second driving term consists in updating the value of the wavelength as a function of the index of refraction at each point. We describe bellow the use of the two driving terms and, in the sections that follow we describe our Gaussian shape invariant proposal and we give some simple examples of its use.

2. Analytical description

To simplify our analytical description, we will refer to one-dimensional waves (cylindrical-like waves) generalization to more dimensions is straightforward. For the description of our Gaussian invariant, we will begin with some preliminary ideas.

2.1 Motivation of changing sin(θ) by tan(θ)

Before establishing our main results, let us motivate the reason why we have to change sin(θ) by tan(θ)in our model for a tilted Gaussian beam.

Let us consider a one-dimensional Gaussian amplitude transmittance, centered at z=0(z being the optical axis), and placed along an axis x, in a coordinate system (x, z), illuminated by a plane wave with an angle of incidence θ. In an equivalent manner, it can also be viewed as a tilted real Gaussian amplitude distribution.

In z=0, this physical situation would be typically expressed as,

$Ψ (x, z = 0) = A \exp (i \frac{2 π \sin (θ)}{λ} x) \exp (- \frac{χ^{2}}{r_{0}^{2}}),$ ,

λ being the wavelength of the illuminating beam, and r ₀ the semi-width of the Gaussian, A being a real amplitude. However, if we introduce this tilted Gaussian into the Fresnel integral, we would not obtain the expected geometrical projection for large θ. However, we have found that changing sin(θ) by tan(θ) allow us to obtain the expected geometrical projection. With this change, the above equation reads,

Ψ (x, y = 0) = A \exp (i \frac{2 π \tan (θ)}{λ} x) \exp (- \frac{x^{2}}{r_{0}^{2}}),

In order to calculate the amplitude distribution at a plane with a coordinate system (ξ, z) located at a distance z, the well-known one-dimensional Fresnel diffraction integral will be used. This integral is expressed as,

Ψ_{F} (ξ, z) = \frac{\exp (i 2 π z ⁄ λ)}{\sqrt{i λz}} \int_{- \infty}^{\infty} Ψ (x, 0) \exp (ì \frac{π}{λ z} {(x - ξ)}^{2}) d x .

By substituting Eq. (1) into Eq. (2), and calculating the integral involved, one obtains,

Ψ_{F} (ξ, z) = \frac{A}{\sqrt{i λ z}} \sqrt{\frac{π r_{0}^{2} λ z}{λ z - i π r_{0}^{2}}} \exp (- i \frac{π r_{o}^{4} α^{2}}{4 λ z}) \exp (- \frac{{(ξ - ξ_{A})}^{2}}{R_{0}^{2}}) \exp (i \frac{π}{λ R_{c}} {(ξ - ξ_{B})}^{2}),

As it can be seen from Eq. (3), after propagation, the tilted beam in Eq. (1) preserves its Gaussian distribution; however it presents a complex amplitude factor and, additionally, a quadratic phase. This result is well known.

The main point to notice in Eq. (3), is that the spatial center of the propagated Gaussian beam, ξ_A=ztan(θ), for any arbitrary angle (|θ|≠π/2), is located precisely at the geometrical projection of the incoming beam center. This result would fail if one maintains sin(θ) for large angles, this change of sin(θ) by tan(θ) removes the paraxial limitation of the Fresnel diffraction integral in this particular case.

Other interesting aspect to notice is that the above result is also true for any arbitrary distance of propagation even for very small ones. It also holds for any value of λ. Additionally it is important the fact that the center of the quadratic phase of the propagated beam, ξ_B, does not coincide in general with the spatial center of the beam, ξ_A. These are the key elements of our proposal that will be generalized in our description of our Gaussian shape invariant under Fresnel propagation.

2.2 Gaussian shape invariant

In order to generalize the above results, we propose a general (one-dimensional) Gaussian amplitude distribution at (z=0) to be represented as

Ψ_{n} (x) = A_{n} \exp (i α_{n} x) \exp (i β_{n} x^{2}) \exp (- \frac{{(x - x A_{n})}^{2}}{r_{n}^{2}}) \exp (i γ_{n} {[x - x B_{n}]}^{2}),

where A_n is the amplitude of the Gaussian beam (which can be complex) and α_n stands for the tilt or the direction of the proposed incoming beam, as described above. The quadratic phase involving β_n is introduced to allow an arbitrary defocus factor in an independent manner. xA_n and r_n are parameters corresponding to the spatial center and to the semi-width of the Gaussian respectively. Finally, the term γ_n is a factor that includes the Gaussian curvature whose quadratic phase is centered at xB_n. As shown above, the center of curvature does not coincide in general with the spatial center of the Gaussian.

Equation (4), represents our proposed Gaussian (shape invariant) amplitude distribution at the origin (z=0), that will be propagated a distance z, towards a coordinate system (ξ, z), by means of the Fresnel diffraction integral. One obtains for the propagated beam, after a straightforward but lengthy calculation,

Ψ_{n + 1} (ξ, z) = A_{n + 1} \exp (i α_{n + 1} ξ) \exp (i β_{n + 1} ξ^{2}) \exp (- \frac{{(ξ - ξ A_{n + 1})}^{2}}{r_{n + 1}^{2}}) \exp (i γ_{n + 1} {[ξ - ξ B_{n + 1}]}^{2}),

where,

A_{n + 1} = A_{n} \frac{\exp (i 2 π z ⁄ λ)}{\sqrt{i λ z}} \sqrt{\frac{π r_{n}^{2} λ z}{λ z - i r_{n}^{2} (β_{n} λ z + γ_{n} λ z + π)}} \times

is the amplitude of the propagated beam, and,

α_{n + 1} = 0, β_{n + 1} = \frac{π}{λ z}, D_{n} = λ^{2} z^{2} + r_{n}^{4} {(β_{n} λ z + γ_{n} λ z + π)}^{2},

The above equations can be seen to represent a set of recursive equations. It can be noticed that as $α_{n} = \frac{2 π}{λ} \tan (θ)$ implies that α_n+1=0. Thus α_n can be identified as a driving term that permits changing the angle or the ray, or the direction of propagation when needed; this term must remain equal to zero all the way that the direction of the beam is maintained unchanged. When the direction of the beam needs to be changed, α _n+1 is updated by means of a law of reflection or refraction; this has to be done for example at each interface, where the beam is reflected or refracted by a surface. In the case that the light passes from one medium to another, α_n+1 is updated by applying Snell’s law to give the beam the appropriate direction and then it is propagated again by means of the recursive equations above. Near field calculations can be used instead for more precision [11, 12].

Additionally,

r_{n + 1} = \frac{\sqrt{D_{n}}}{π r_{n}}, ξ A_{n + 1} = x A_{n} + \frac{α_{n} λ z}{2 π} - \frac{γ_{n} λ z (x B_{n})}{π} + \frac{(β_{n} + γ_{n}) λ z (x A_{n})}{π}

and

ξ B_{n + 1} = \frac{α_{n} λ z}{2 π} - \frac{γ_{n} λ z (x B_{n})}{π} - \frac{λ^{2} z^{2}}{β_{n} λ z + γ_{n} λ z + π} \frac{x A_{n}}{π r_{n}^{4}} .

The set of Eqs. (5A)-(5D) corresponds to exact calculations and can be applied in general. These equations represent the complete set of recursive equations needed in our proposal. It has to be mentioned however, that the Fresnel diffraction integral imposes restrictions because it corresponds to a paraxial theory and it is not considered accurate in the very-near field theory. Thus, in principle, the Fresnel integral demands that two restrictions have to be overcome: i) It is only valid for moderate angles of incidence, due to the paraxial restriction and ii) it is invalid for very-small distances of propagation.

As shown above, in our proposal we are using two facts to overcome the mentioned limitations. First $α = \frac{2 π \tan (θ)}{λ}$ substitutes $α = \frac{2 π \sin (θ)}{λ}$ , allowing arbitrary angles of propagation, to any value because, as we have shown, an accurate geometrical projection of the spatial center of the Gaussian is obtained at any arbitrary angle of incidence. Second, the moderate distance restriction can be seen to be overcome as we have shown that the geometrical projection of the ray is obtained independently of the values of z, additionally, also of λ. Thus the ray is properly projected even if both, z and λ, are very small.

It is immediate to recognize that the second driving term to direct the ray is the wavelength, λ. For its use, at the beginning of the calculations an initial value α ₀ and λ ₀ must be provided. The smaller the wavelength the more a ray-like behavior because the width of the beam will remain basically unchanged under propagation. Additionally, a very small semi-width can be chosen, thus actually propagating rays.

In the cases where the index of refraction of the media varies continuously, it is sufficient to update the value of the wavelength as λ/n (at each iteration), n being the index of refraction. This automatically causes the ray to change its direction accordingly, making the calculations very easy.

Thus, two driving direction terms are available, α and λ. Which of them or which combination of them are to be used as the driving parameters in a specific ray tracing situation, depends on which parameter makes the calculations easier. For example, if a ray has to be traced when there is only a limited number of interfaces, it is easier to provide the new value of α after refraction or reflection at an interface by applying Snell’s law. On the other hand, if a ray travels through an optical media where the index of refraction changes continuously, updating the value of λ as a function of position results easier.

Based on the above description, we proceed to show some simple simulations of ray tracing with the proposed technique. As indicated, the technique is applied in an iterative manner; the overall path of propagation is divided in small sub-paths. Each sub-path corresponds to a fraction of the total path length. The calculations are as follows: we calculate the path of the ray for a very small propagation (sub-path), and then to the next, using the former result, and so on until the ray arrives to the end of the path. When very high accuracy is needed, the length z in each sub-path can be made as small as desired but of course the limit is the computer precision; additionally, very small values will require more iterations, thus, more computing time. Finally, for the easy of the calculations, the length of each sub-path is considered equal, although this value can be made different as needed. In the next section we show some simple examples.

3. Computer examples

In the examples that follow the iterative Eqs. (5A)-(5D) are used, the spatial center of the Gaussian invariant, which represents the ray path, is plotted as described above. We have chosen arbitrarily very small values, r ₀=10^-7 m and λ=10^-15 m for a ray-like behavior, although this is not a necessary condition. We repeated the calculations with more realistic values, and we obtained the same ray tracing. 4000 iterations were performed for each ray trace. The total length of the optical axis was divided in 4000 equal sub-paths, thus, for each iteration, z=z_T/4000, being z_T the length of the optical axis. Due to the very small wavelength, the semi-width shows to remain basically constant in all the calculations. The simulations were performed in a personal computer running at 1.2 GHz, using commercial mathematical software, the total processing time for each ray being approximately 1/2 second. The results were compared with the one obtained by means of a commercially available program and found to coincide.

Our first example, shown in Fig. 1, consists in rays in a gap between two reflective planes. The distance between planes can be chosen arbitrarily, very close or very far as desired, there are no limitations. In our example, the upper plane has a slope of 0.1 and is placed at a height of 1 cm. The lower plane is placed at a height of −2.0 cm with a slope of −0.1. The index of refraction in the gap is taken equal to one. For illustrative purposes only two rays are traced. The dash-dotted ray is launched at the origin at the left, with an initial inclination of -0.7 rad. The solid ray is also launched at the left at a height -1.5 cm, with an angle of 0.8 rad. The length of the optical axis is, z_T=0.5 m.

Fig. 1. Ray tracing in a gap between two reflective planes (units in m).

Download Full Size | PDF

Our second example consists in ray tracing through a selfoc media whose equation is given by $n (x) = 1.5 \sqrt{1 - (8) 10^{4} x^{2}}$ (x in units of meters). According to the Eikonal equation a ray launched at the origin, with an inclination ϑ=0.056 rad, will follow a sinusoidal path with a period $T = 2 π \cos ϑ ⁄ \sqrt{(8) 10^{4}}$ and with amplitude $A = \sin ϑ ⁄ \sqrt{(8)} 10^{4}$ , Ref. [13]. Figure 2 shows the path traced by means of the Gaussian invariant. The resulting path agrees well with the Eikonal theory, in this case, resulting in a period of 0.022m, and an amplitude (1.9)10^-4 m.

Fig. 2. Ray tracing in a selfoc media as described in the text.

Download Full Size | PDF

Our next example, shown in Fig. 3, consists in a ray tracing through a plano-convex lens. The convex surface is given as $z (x) = \sqrt{{(64) 10}^{- 4} {- (x - 0.15)}^{2}}$ (units in meters). The indexes of refraction are 1.5 inside the lens and 1.0 for the optical media.

Fig. 3. Ray tracing through a convex lens as described in the text.

Download Full Size | PDF

In Fig. 4, the convex surface above was replaced by an aspheric in order to lower defocus. The equation of the aspheric surface is given as $z (x) = \sqrt{(64) 10^{- 4} - {(x - 0.15)}^{2}} + (5) 10^{3} x^{6}$ (units in meters).

Fig. 4. Ray tracing through an aspheric lens as described in the text.

Download Full Size | PDF

In all of the above examples the driving factor α remains equal to zero and it is updated by using Snell’s law only when the ray intersects a refractive surface.

Our final example, Fig. 5, uses a combination of both driving terms, α and λ. In this example, an inhomogeneous media is placed near the back focal plane of the above aspheric lens; the index of refraction as a function of z being n(z)=1.8-cos(2πz/0.2). In this case we are using 8000 iterations for a better resolution. The driving term used in the inhomogeneous media is λ. When the ray enters the media an initial value of λ=10^-15 m is assigned and this value is then updated at each iteration as λ/n. It is only necessary to apply Snell’s law at the point where the beam enters the inhomogeneous media to provide an initial angle and then, α remains zero in the entire region. This example shows the simplicity of the calculations in the inhomogeneous region as it is only necessary to update the value of the wavelength in each iteration.

Fig. 5. Ray tracing through an aspheric lens and an inhomogeneous media placed near the back focal plane of the lens as described in the text.

Download Full Size | PDF

4. Comparison with traditional ray tracing

In this section we establish the advantages of using our Gaussian shape invariant method over traditional ray tracing.

In traditional ray tracing only two mathematical terms are tracked after the ray has traversed any given complex optical system; one is, the overall phase due to the whole optical path traveled; the other one is the output tilt of the ray at the output of the system. Additionally in standard ray tracing the amplitude of the ray remains unchanged as it travels throughout the system. In contrast in our proposed technique all the inherent (four) properties of the (physical) ray as it propagates are preserved. These are:

1. The now (overall) complex amplitude which includes the amplitude of the entrance ray as well as its optical path traversed (let, A_n=|A_n|exp(iϕ_n)) where ϕ_n corresponds to the accumulated optical path.

2. The output ray’s tilt (as in traditional ray tracing) is also included in our (physical) ray model in the term $\exp (i \frac{2 π \tan (θ)}{λ} x)$ .

3. The propagated ray curvature is tracked throughout the whole optical system. This output value may be useful in interferometric and holographic applications where the overall ray phase is important.

4. Finally, as one can select real values for the wavelength’s (physical) ray, the amplitude spread may also be tracked. In contrast, to obtain an estimation of the energy spread in traditional ray tracing it would be necessary to estimate the relative density of the rays over the image surface.

All of the above comments are mainly valid in the paraxial approximation that includes most of the optical systems in use.

5. Conclusions

We have described a technique for ray tracing in homogeneous and/or inhomogeneous media, based in the use of two driving terms to direct the ray in an appropriate direction. We have shown by means of simple examples how both driving terms can be used in an independent way or combined according to the optical situation, making apparent the versatility of the technique from a computational point of view. Although the technique is based in the propagation of a Gaussian shape invariant under the Fresnel diffraction integral, valid only in a paraxial region and not considered as a very-near field propagator, these two limitations were overcome by replacing the commonly used sine function by a tangent in the phase of the amplitude distribution of the Gaussian, and by showing that it possible to obtain an accurate geometrical projection at any arbitrary distance. With these considerations, we have shown the feasibility of the technique by means of simple examples of ray tracing through some imaging devices and optical media, where very small values of wavelengths and very small beam widths were used, so that very narrow like-rays beams were attained. Although not shown here, as the technique is based on the theory of diffraction, it can result useful in the study of ray tracing on systems where diffractive effects have to be included, this can be taken into consideration in future work.

References and links

1. V. N. Mahajan, Optical Imaging and Aberrations Part I Ray Geometrical Optics, (SPIE Press USA, 1998).

2. M. Herzberger, Modern Geometrical Optics, (Interscience Publishers, Inc. N.Y., 1958).

3. O. N. Stavroudis, “Simpler derivation of the formulas for generalized ray tracing,” J. Opt. Soc. Am. 66, 1330–1333 (1976). [CrossRef]

4. G. H. Spencer and M. V. R. K. Murty, “General ray-tracing procedure,” J. Opt. Soc. Am. 52, 672–678 (1962). [CrossRef]

5. O. N. Stavroudis, The Optics of Rays, Wavefronts, and Caustics, (Academic Press, Inc. N.Y., 1972).

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

7. S. Nolte, M. Will, J. Burghoff, and A. Tünerman, “Ultrafast laser processing: new options for three-dimensional photonic structures,” J. Mod. Opt. 51, 2533–2542 (2004). [CrossRef]

8. S. Szatmári and G. Kühnle, “Pulse front and pulse duration distortion in refractive optics, and its compensation,” Opt. Commun. 69, 60–65 (1988). [CrossRef]

9. U. Fuchs, U. D. Zeitner, and A. Tőnnermann, “Ultra-short pulse propagation in complex optical systems,” Opt. Express 13, 3852–3861 (2005). [CrossRef] [PubMed]

10. A. Rohani, A. A. Shishegar, and S. Safavi-Naeini, “A fast Gaussian beam tracing method for reflection and refraction of general vectorial astigmatic Gaussian beams from general curved surfaces,” Opt. Commun. 232, 1–10 (2004). [CrossRef]

11. J. J. Stamnes, Waves in Focal Regions, (IOP Publishing Limited, 1986).

12. L. Novotny and B. Hecht, Principles of Nano-Optics, (Cambridge University Press, 2006).

13. K. Iizuka, Engineering Optics, (Springer Verlag, 1986).

Fresnel-Gaussian shape invariant for optical ray tracing

Abstract

1. Introduction

2. Analytical description

2.1 Motivation of changing sin(θ) by tan(θ)

2.2 Gaussian shape invariant

3. Computer examples

4. Comparison with traditional ray tracing

5. Conclusions

References and links

Cited By

Figures (5)

Equations (12)

Optics Express