## Abstract

Full characterization of optical systems, diffractive and geometric, is possible by using the Fresnel Gaussian Shape Invariant (FGSI) previously reported in the literature. The complex amplitude distribution in the object plane is represented by a linear superposition of complex Gaussians wavelets and then propagated through the optical system by means of the referred Gaussian invariant. This allows ray tracing through the optical system and at the same time allows calculating with high precision the complex wave-amplitude distribution at any plane of observation. This method is similar to conventional ray tracing additionally preserving the undulatory behavior of the field distribution. That is, we are propagating a linear combination of Gaussian shaped wavelets; keeping always track of both, the ray trajectory, and the wave phase of the whole complex optical field. This technique can be applied in a wide spectral range where the Fresnel diffraction integral applies including visible, X-rays, acoustic waves, etc. We describe the technique and we include one-dimensional illustrative examples.

©2011 Optical Society of America

## 1. Introduction

There are basically two approaches to characterize an optical system: geometric, based on ray tracing methods, and diffractive, based on calculations of diffraction integrals. Both techniques are necessary for a complete characterization. It is however difficult to conciliate both techniques in a given system. For example, the impulse response or the point spread function (PSF), which rigorously has to be calculated by diffractive means is commonly estimated by considering a high density of rays from which an approximation of the diffractive response is obtained [1–3]. A recent technique for improving the result obtained by ray density can be found in [4]. Furthermore, if one is interested in calculating the response of an optical system in regions other than the vicinity of the focal region, one may try to generate the so called spot diagrams [5]; this will give an estimation of the intensity distribution at a plane of interest; however, with geometric techniques, it will not be possible to calculate the corresponding complex diffraction patterns.

In this report we describe how the Fresnel Gaussian Shape Invariant (FGSI) recently reported [6] conciliates both techniques allowing calculating the complex amplitude distribution at any plane of observation and at the same time attaining ray tracing as in conventional geometrical methods. The name FGSI refers to the way the complex Gaussian is represented, preserving its mathematical shape under propagation using the Fresnel diffraction integral.

As described below, in order to obtain the diffraction pattern at an image plane by means of FGSI, it is necessary to represent the amplitude distribution in the object plane as a superposition of complex Gaussian wavelets as in [7,8]. That is, the object wave-front is sampled by a finite linear combination of complex, Gaussian shaped wavelets. Given that these Gaussian wavelets are spatially localized, they preserve both, the ray and the wave character of the propagated wave-front. The propagation of a Gaussian-shaped localized wave-front, always give a Gaussian-shaped localized response; this is the principle of the herein presented technique. The method consists in tracing rays using the laws of reflection and refraction at each interface, each ray being a FGSI, and then calculating the superposition of all the FGSI traced to obtain the resulting complex diffraction distribution.

We illustrate the use of FGSI by means of simple, one-dimensional, examples. Generalization to two dimensions is straightforward.

## 2. Analytical description

#### 2.1 Description of FGSI

FGSI is reported in [6] only for ray tracing purposes. In this report we extend its applicability by introducing new equations to calculate associated diffraction patterns. First, for convenience, a brief description on FGSI based in [6] follows.

Let us consider a one-dimensional (or cylindrical) Gaussian amplitude distribution located at a plane $\left(x,z=0\right)$ given by,

*n*being an integer introduced to be used for iterative purposes.In Eq. (1) ${P}_{n}$ is the amplitude of the Gaussian beam, in general a complex number. The values ${\alpha}_{n}$ stand for the tilt (direction) of the beam. The quadratic phase ${\beta}_{n}$ is introduced to allow an arbitrary defocus factor, ${A}_{n}$ corresponds to the spatial center and ${r}_{n}$ to the semi-width. Finally, the term ${B}_{n}$ gives the center of the Gaussian curvature${\gamma}_{n}$.By analyzing Eq. (1), one notices redundancy in the definition of some of the parameters in the way it is expressed; however, as we have described in [6], writing the amplitude distribution in this form Eq. (1) will allow us to obtain an iterative system of equations and the propagated field will preserve this mathematical shape. This is desirable for programming purposes. As usual,

*z*represents the optical axis.

To calculate the amplitude distribution, ${\Psi}_{n+1}\left(\xi ,z\right)$, at a coordinate plane $\left(\xi ,z\right)$, we will use the one-dimensional Fresnel diffraction integral,

*ξ*corresponds to a coordinate axis parallel to

*x*situated at a distance

*z*from the object plane;

*λ*is the wavelength of the field between the two planes.

By substituting Eq. (1) into Eq. (2) and performing the integral, one obtains the amplitude distribution (same shape) of the propagated field as,

The set of Eqs. (3)A-3G) represents our algebraic iterative system to calculate the propagated field from the n-plane to the (n + 1) plane. It will be noticed, in particular, that ${\alpha}_{n+1}=0$. This parameter will remain equal to zero until a change in the direction of the ray due to reflection or refraction has to be made, as described below.

First, conventional ray tracing is performed through the optical system. This can be done conventionally or by means of the iterative process of FGSI using the set of Eqs. (3)A)-(3G) as described and illustrated in [6] Then, each ${A}_{n}$ (the Gaussian spatial center), represents the local position of each ray. Additionally, the wave character of the whole wave-front is always available. This is because we are using a complex Gaussian wavelet where the wavy character is always present. Some examples on ray tracing by FGSI are presented in [6].

Now that FGSI for ray tracing has been described, we will introduce new equations to calculate the corresponding diffraction patterns. For this, let us begin the calculation of the complex wave amplitude of the Gaussian field, by considering the parameters${\alpha}_{n}$in Eq. (1). These parameters represent a tilt or a direction of the ray. As well-known, in the paraxial approximation the term $2\pi \mathrm{sin}\left(\theta \right)/\lambda $ will direct a wave-front with a wavelength *λ* at an angle equal to *θ*. However in our case, as described in [6], we are interested in non-paraxial rays, so we need to change the sine function by the tangent one as,

According to Eq. (3)C), the values ${\alpha}_{n}$ remain equal to zero until a change in the ray angle has to be made due to reflection or refraction after the rays have being traced. When the ray angle has to change, one has to update the value${\alpha}_{n+1}$. This is done by means of the following equation which is obtained from Eq. (3)E) as,

Obviously, if a change in the index of refraction occurs, the value${\lambda}_{n+1}$has also to be updated accordingly before applying Eq. (4)B).

At this point the value of the amplitude ${P}_{n+1}$ calculated by means of the iteration given in Eq. (3)B) will be incorrect due to the changes that have been introduced to the wave. Thus, by demanding continuity of the amplitude, one has to update its value as,

Once Eqs. (4)B-4C) have been applied, the iterative process continues by using (3A-3G). This procedure permits us to calculate the amplitude distribution at any desired plane of observation and at the same time to perform ray tracing. Equations (4)B-4C) are the required new set of equations that will be used to calculate the associated diffraction patterns.

In the following sub-section we describe how FGSI is applied to propagate the amplitude distribution of an object plane through an optical system to allow its whole characterization.

#### 2.2 FGSI applied to the propagation of complex fields

In order to propagate a complex field distribution located at an object plane through an optical system by means of FGSI, it is necessary first to represent the object field as a superposition of Gaussian wavelets. Our sampling proposal follows a Rayleigh-like criterion as described by us in [7,8], although other sampling methods had been proposed [9,10]. Sampling following a Ryleigh-like criterion greatly facilitates the calculations and the programming with the additional advantage that it represents accurately the object spatial spectrum as we have demonstrated in [8]. The accuracy of the superposition depends only on the number of Gaussian wavelets used. If additionally, each Gaussian wavelet in the superposition is expressed as a FGSI, then the propagation of the object field through an optical system will result to be a simple algebraic iterative process by using the equations introduced in the above sub-section.

Performing the propagation by means of FGSI will give additional advantages. To see this, let the one-dimensional complex amplitude distribution at the object plane be expressed as${\Psi}_{T}\left(x\right)$ = $R\left(x\right)\mathrm{exp}\left[i\phi \left(x\right)\right]$, being $R\left(x\right)$ and $\phi \left(x\right)$ real functions. We place the object plane at the zero-plane (initial conditions). Then, the object field is represented as,

In Eq. (5), the subscript zero refers to initial conditions at the object plane. The subscript *m* is used for the superposition process. As Eq. (5) refers to initial conditions, ${B}_{0,m}$ was made equal to${A}_{0,m}$, and $R\left({x}_{m}\right)={P}_{0,m}$ for each ${x}_{m}={A}_{0,m}$in the object plane.

If one rewrites Eq. (5) as,

From Eq. (8) the Gaussian local parameters are obtained as,

Equations (9)A) and (9B) reveal that the Gaussian parameters in each FGSI hold essential physical information of the object field, namely, the local slope, curvature and constant initial phase. These parameters will be available over the entire propagation. Additionally, as indicated above, the spatial frequency of the object plane will be properly represented. Thus, the herein technique turns out to be a powerful tool allowing a whole characterization of optical systems.

Now we have the set of equations required to completely characterize an optical system by means of FGSI. Before presenting our examples, a summary of the use of the technique follows:

First, the amplitude distribution at the object plane is represented by a superposition of Gaussian wavelets. For the one-dimensional case a rectangle function will typically represent the amplitude distribution at the object plane. This will be the case of the examples presented in the next section. To direct the rays to a desired direction, an initial tilting- phase (${\alpha}_{0}$) is assigned. At each interface the angle of the emerging ray has to be calculated using the laws of reflection or refraction and then the associated diffraction pattern is calculated by means of Eqs. (4)A)- (4C).

Second, for each Gaussian in the superposition, the values of their spatial centers will give the spatial local position for each ray. Initial conditions are assigned by means of Eqs. (9)A-9B). The propagation of each Gaussian is calculated using Eqs. (3)B-3G) and (4B-4C) each time a reflection or a refraction occurs. Finally, at the plane of observation, the amplitude of each Gaussian is superimposed (added coherently) to obtain the desired complex amplitude distribution. It is worth mentioning that a mutual coherence function can be assigned to the amplitude of each Gaussian, although this is not presented in this report.

As a final comment, we notice, that if instead of using the overall FGSI technique as describe above, one traces independent rays and then plans to apply the superposition of Gaussian wavelets only at the image plane, one would assign a complex Gaussian to each ray after propagation. However by doing this, the Gaussian representation (at the image plane) will contain incorrect phases and the diffraction pattern will only be noise. In contrast, in our FGSI technique the Gaussian sampling is performed at the object plane. At this plane the wave-front and its energy distribution are highly smooth. In a highly smooth light distribution, uniform spaced Gaussian wavelets will give the best representation [7,8]. Once the representation of the object wave-front is appropriate, and due to the linearity of the Fresnel integral, the subsequent propagated complex distributions will be equally well represented.

In the following section we present one-dimensional examples.

## 3. Examples

For our first set of examples a singlet with an index of refraction of 1.5163 (BK7) whose first surface has a radius of curvature of 2cm will be considered, (Fig. 1
). The illuminating source, at the left, is a plane wave (*λ* = 0.5 μm) emerging from a rectangular aperture centered at the origin and exhibiting a width of 2mm. These parameters were chosen arbitrarily for illustrative purposes.

For brevity we will only show normalized intensities of the diffraction patterns, although as indicated, FGSI allows obtaining the whole complex amplitude distributions.

For all the examples the rectangular aperture at the object plane are represented by 400 Gaussian wavelets (rays) equally spaced following a Rayleigh-like criterion as described in [7,8]. Additionally, 500 sampling points are used for representing each plot. These numbers can be chosen freely, according to the limitations of sampling (Nyquist) conditions. For clarity of the plots, only five rays are shown, the two outermost, two in the middle and, the central ray. The processing time in a double-core PC at 2 GHZ was approximately 15 seconds for each simulation.

In Fig. 1 the rays emerging from the object plane, at the left, are directed parallel to the optical axis.

By calculating the diffraction pattern by means of FGSI at different observation planes, we have found the back focal length (f) to be approximately at 3.74 cm. The corresponding normalized intensity distribution at the focal plane is shown in Fig. 2 .

As expected, a centered sinc function represents the PSF.

In Fig. 3 a similar calculation is performed but this time the incident beam is tilted 6°.

Figure 4 below shows the corresponding normalized intensity distribution at the back focal plane.

As expected the intensity distribution is shifted, defocused, and non-symmetric. It will be noticed that FGSI provides fine details of the distribution. It is important to remark that an additional advantage of the technique is that one can chose the field of view as desired maintaining the same sampling as described in [8].

In order to confirm the correctness of the shape of the intensity distribution of Fig. 4, it would be convenient to calculate numerically by direct use of the Fresnel diffraction integral the same pattern. This however is not an easy task as to sample the quadratic phase and at the same time due to such amount of tilt of the incoming beam it would require an extremely large number of sampling pixels. As an alternative approach we used the Fresnel diffraction integral to calculate a similar physical situation. First, we decreased the tilt in order to permit us attain adequate sampling; in this case we tilted the incoming beam to only 0.3°. Second, we heuristically propose a phase polynomial attaining defocus and spherical aberration terms which are symmetric around the optical axis to simulate the physical situation; the proposed polynomial reads:

The aberration values in Eq. (10) were chosen rather arbitrarily to allow adequate sampling and only for comparative purposes. By performing numerically the two propagations required using the Fresnel diffraction integral, one obtains the following intensity distribution.

Obviously Figs. 4 and 5 are not exactly the same as each one obeys to different numerical physical parameters. It can be appreciated that there is a clear similarity in their shape, making apparent the correctness and the advantages of using FGSI.

Figures 6 and 7 below, show normalized intensity distributions at an observation plane located a distance f/2 to the right of the back focal plane for two situations of the incident beam: i) directed parallel to the optical axis and ii) tilted 10°.

In our next examples a centered spherical mirror whose vertex is placed at the origin is considered. The radius of the mirror is 2 cm. The object rays emerge form a slit located at a distance of 2 cm at the right of the mirror. The aperture of the slit is 2mm. Figure 8 below shows the physical situation. The rays are traced parallel to the optical axis. By performing several calculations using FGSI, we found the focal length to be of approximately 1 cm. Figure 9 below shows the corresponding normalized intensity distribution at the focal plane.

In Fig. 10 below, the object rays are directed at an angle of 8°, Fig. 11 shows the corresponding normalized intensity distribution at the focal plane.

It will be noticed that, in contrast to other techniques, FGSI gives fine details of the diffraction patterns without requiring additional calculations.

Before finishing this report it is worth mentioning that we have applied the method in different regions of the spectrum and found that FGSI give accurate results using basically any wave-length, including acoustic waves and X-rays. FGSI may be also useful in considering polarization effects by adding Fresnel coefficients. As indicated above, the processing time has to be considered especially if one is interested in fast methods, as the processing time for each simulation took approximately 15 seconds in a double-core PC running at 2GHz.

## 4. Final remarks

As we have seen above, the linear superposition of Gaussian wavelets in combination with FGSI represent a powerful tool for a whole characterization of optical systems. It is now commercially available a program called ASAP developed by Breault Research Organization [11] which uses a superposition of Gaussian wavelets in a similar way as to our proposal. The inner-workings of this software are not of the public domain. However it seems by the information available that ASAP has to trace three rays linked to each Gaussian wavelet in the superposition; these rays are named by them the base, the waist and the divergence rays. Each one of these rays follows different trajectories and at the final step these rays are somehow combined to obtain the desired diffraction pattern. In contrast, each one of our FGSI wavelets carries the whole local optical information: amplitude, constant phase, tilt, and curvature. Thus, if the object plane is properly sampled by means of our Rayleigh-like criterion, due to the linearity of the Fresnel diffraction integral, the subsequent propagations will accurately represent the field under inspection.

Our proposed method can be applied to general optical systems including those using selfoc or gradient-index. In these cases the region where the index of refraction varies as a function of position is partitioned into several thin vertical strips. On each of these thin sub-regions a representative constant value for the index of refraction is assigned; this value may in turn be a function of the height of each ray. Then, conventional ray tracing is performed over the whole region by refracting the trajectory of the rays at each interface. If desired, Fresnel coefficients can be assigned under refraction at each interface. Finally FGSI equations are applied to each trajectory for obtaining the corresponding diffraction patterns.

It has to be mentioned that a limitation of these techniques which are based in the superposition of wavelets arises when an aperture in the optical system truncates part of the propagated wave-front. In this case the emerging diffraction pattern will consist of the product of the wave-front by the aperture function. Thus, the emerging diffraction pattern will in general be truncated and the superposition representation will be inaccurate. We might think of at least two ways to solve this problem: i) reassigning a new superposition to the truncated wave-front or ii) we may also represent the blocking aperture by another superposition of Gaussian wavelets. The resulting wave-front emerging from this blocking aperture will consist of the product of two wavelet representations, which in turn will be a new Gaussian-wavelet representation to be further propagated by FGSI. At the moment of this report we have not yet addressed this situation.

## 5. Conclusions

We have presented the analytical equations for extending the applicability of the Fresnel Gaussian shape invariant for calculating diffraction patterns of optical systems. With this extension it is now possible to calculate diffraction patterns with fine detail in basically any region of interest and at the same time attaining conventional ray tracing. In our Gaussian invariant technique the most important physical parameters of the wave under propagation are always available, the amplitude, direction and curvature; additionally, when properly sampled, also the spatial spectrum. The technique is simple as it requires only algebraic iterations. In view of all this, the present technique turns out to be simple and powerful for a whole characterization of optical systems.

We have illustrated the technique by means of simple one-dimensional examples. Although not shown here the technique can be applied for basically any wave length including acoustic waves and X-rays. The method can be also used for studies using partial or total coherence.

## References and links

**1. **W. T. Welford, *Aberrations of optical systems*, (Academic Press, 1974), pp. 220.

**2. **R. Kigslake, *Lens design fundamentals*, (Academic Press, 1978), pp. 8.

**3. **W. J. Smith, *Modern Optical Engineering*, 3rd ed., (Mc Graw-Hill, 2000) p. 372.

**4. **C.-S. Liu and P. D. Lin, “Computational method for deriving the geometric point spread function of an optical system,” Appl. Opt. **49**(1), 126–136 (2010). [CrossRef] [PubMed]

**5. **M. Herzberger, *Modern Geometrical Optics*, (Interscience Publishers, Inc. N.Y., 1958) pp. 383–400.

**6. **M. Cywiak, A. Morales, J. M. Flores, and M. Servín, “Fresnel-Gaussian shape invariant for optical ray tracing,” Opt. Express **17**(13), 10564–10572 (2009). [CrossRef] [PubMed]

**7. **M. Cywiak, M. Servín, and F. Mendoza-Santoyo, “Wave-front propagation by Gaussian superposition,” Opt. Commun. **195**(5-6), 351–359 (2001). [CrossRef]

**8. **M. Cywiak, A. Morales, M. Servín, and R. Gómez-Medina, “A technique for calculating the amplitude distribution of propagated fields by Gaussian sampling,” Opt. Express **18**(18), 19141–19155 (2010). [CrossRef] [PubMed]

**9. **A. W. Greynolds, “Propagation of generally astigmatic Gaussian beams along skew ray paths,” Proc. SPIE **560**, 33–50 (1985).

**10. **M. M. Popov, “A new method of computation of wave fields using Gaussian beams,” Wave Motion **4**(1), 85–97 (1982). [CrossRef]