A technique for calculating the amplitude distribution of propagated fields by Gaussian sampling

Moisés Cywiak; Arquímedes Morales; Manuel Servín; Rafael Gómez-Medina

doi:10.1364/OE.18.019141

1. Introduction

As is well known, the Fresnel diffraction integral is the key tool to calculate the amplitude distribution of propagated fields. However, due to the presence of a quadratic phase in the integrand, the integral is in general not analytically integrable and it is not possible to attain a closed calculation. Thus, several numerical methods have been developed as is the case of the Fast Fourier Transform (FFT), Lommel integrals, sampling expansions and techniques involving methods of stationary phase [1–7].

Sampling expansions consists basically in three steps: First, the diffraction integral is written in cylindrical coordinates giving rise to a Bessel function of the first kind and zero order in the integrand. Then, the product of the quadratic phase by the Bessel function is represented by a sampling series whose coefficients must be determined. Finally, the resulting integrals are evaluated. Unfortunately, sampling expansions as well as Lommel integrals are only applicable to situations with axial symmetry.

The numerical techniques mentioned above are limited when the quadratic term in the integrand oscillates rapidly imposing severe sampling (Nyquist) conditions; for example, for calculating the amplitude distribution of a converging lens, with a focal length of 5cm and an aperture of 2.5cm at a plane located 1cm away from the focal plane, would require sampling the quadratic term with approximately 17000 by 17000 pixels, resulting in a large computing time. If additionally, a system require multiple calculations the computing time may become prohibitive.

The above number of sampling pixels can be estimated as follows: Let the illuminating wavelength be 0.6μm. The real part of the quadratic phase corresponds to a cosine function; then, taking a line scan over one diameter results in 868 cycles with different widths, being the length of the narrower period 14.4 μm. If we assign ten pixels homogeneously distributed to this zone, then it will require 17000 pixels to homogeneously sample the whole diameter. Additionally, an equal number of pixels will be also necessary for the imaginary part. For interferometric and holographic applications where the complex phase has to be calculated with precision, it will be necessary to increase the number of pixels even more.

To overcome the above limitations we present a technique that is based in representing one or several terms in the integrand inside the Fresnel diffraction integral as a superposition of Gaussian functions. This allows us to calculate the integral in an analytical (closed) form. As in general the terms in the integral (not the quadratic phase) require a low sampling rate, the numerical calculations are performed faster. As indicated, this is especially important when considering multiple optical diffractive systems.

Although the possibility of representing a function by a superposition of Gaussians has previously been mentioned [8–10], it has not yet been used in the form presented in this report. The only restriction that will be imposed is that the functions to be represented by the superposition of Gaussians have to be band-limited.

2. Analytical description

In this section we describe the process of assigning a superposition of Gaussians to an arbitrary band-limited function. In order to illustrate how to assign the Gaussian functions to a particular case and also to find the limitations of the technique, it is convenient first to consider intuitively the consequences of applying the method to a non-band limited function. For this, let us consider a rectangle function $r e c t (x, A)$ defined as usual, unitary in the interval $- A / 2 \leq x \leq A / 2$ and zero otherwise.

For our intuitive approach an odd number of Gaussians will be assigned. Each Gaussian assigned will have the same width and will be placed following a Rayleigh-like criterion. Let us begin with three Gaussians, one Gaussian centered at the origin, one centered at $- A / 2$ and the other at $A / 2$ . The amplitude of each Gaussian is chosen equal to the value of the function to be represented at this point, in this case, to unity. Our superposition of Gaussians then reads Eq. (1),

f (x) = \exp (- \frac{{(x + A / 2)}^{2}}{σ^{2}}) + \exp (- \frac{{(x)}^{2}}{σ^{2}}) + \exp (- \frac{{(x - A / 2)}^{2}}{σ^{2}}) .

As a Rayleigh-like criterion is proposed the width achieves

2 σ = A

. Proceeding in this manner, we assign five Gaussians thus Eq. (2),

f (x) = \exp (- \frac{{(x + 2 A / 4)}^{2}}{σ^{2}}) + \exp (- \frac{{(x + A / 4)}^{2}}{σ^{2}}) + \exp (- \frac{{(x)}^{2}}{σ^{2}}) + \exp (- \frac{{(x - A / 4)}^{2}}{σ^{2}}) + \exp (- \frac{{(x - 2 A / 4)}^{2}}{σ^{2}}),

and now,

4 σ = A

.

For N Gaussians one gets the general equation [Eq. (3)],

f (x) = \sum_{n = - \frac{N - 1}{2}}^{\frac{N - 1}{2}} \exp (- \frac{{(x - n \frac{A}{N - 1})}^{2}}{{(\frac{A}{N - 1})}^{2}}),

N being an odd integer.

Now, let $F (u)$ be the Fourier transform of $f (x)$ as Eq. (4),

F (u) = \int_{- \infty}^{\infty} f (x) \exp (- i 2 π u x) d x,

then,

F (u) = \sqrt{π} \frac{A}{N - 1} \exp (- {(\frac{π A u}{N - 1})}^{2}) \sum_{n = - \frac{N - 1}{2}}^{\frac{N - 1}{2}} \exp (- i \frac{2 π n A u}{N - 1}) .

The summation in Eq. (5) is well known and can be written in closed form, see for example [11]. Thus, Eq. (5) can be written as,

F (u) = \sqrt{π} \frac{A}{N - 1} \exp (- {(\frac{π A u}{N - 1})}^{2}) \frac{\sin (\frac{N}{N - 1} π A u)}{\sin (\frac{π A u}{N - 1})} .

ForN large, and ubounded, Eq. (6) can be written in a good approximation as,

F (u) = \sqrt{π} A \frac{\sin (π A u)}{π A u} .

Aside of the term π, Eq. (7) is the exact Fourier transform of our rectangle function with arbitrary width A.

From the above description the following can be concluded:First, the superposition of equally spaced Gaussians following a Rayleigh-like criterion represents, with a reasonable accuracy, a band-limited function. Second, the amplitude of each Gaussian is taken equal to the value of the function under representation at its corresponding point (or pixel) multiplied by $1 / \sqrt{π}$ . Obviously when the function under representation is complex, actually two superposition of Gaussians are used, one for the real part and another for the imaginary part.

To generalize the above results let us consider an arbitrary band-limited function f exhibiting a maximum frequency $U_{M}$ , sampled with a function g as,

f_{s} (x) = \sum_{n = - \infty}^{\infty} f (n Δ) g (x - n Δ),

where Δ is the sampling interval.

Fourier transforming Eq. (8) one gets,

F_{s} (u) = G (u) \sum_{n = - \infty}^{\infty} f (n Δ) \exp (- i 2 π n Δ u),

where

F_{s} (u)

and

G (u)

are the Fourier transforms of

f_{s} (x)

and

g (x)

respectively.

By using Poisson’s formula, Eq. (9) can be rewritten as,

F_{s} (u) = \frac{1}{Δ} G (u) \sum_{n = - \infty}^{\infty} F (u - n / Δ) .

In our case the sampling function is a Gaussian Eq. (11),

g (x) = A \exp (- \frac{x^{2}}{σ^{2}}),

then its Fourier transform is given by Eq. (12),

G (u) = A \sqrt{π} σ \exp (- π^{2} σ^{2} u^{2}) .

Thus, Eq. (10) can be written as

F_{s} (u) = \frac{A}{Δ} \sqrt{π} σ \exp (- π^{2} σ^{2} u^{2}) \sum_{n = - \infty}^{\infty} F (u - n / Δ) .

The Gaussian function in Eq. (13) has a cut-off frequency equal to $1 / π σ$ . Then in order to limit aliasing, due to the properties of the Gaussian function, one can impose,

1 / π σ = K U_{M} and 1 / Δ - U_{M} = K (1 / π σ) with K > 1,

being K a real number.

To fit with the conditions given by Eq. (14) one suitable value is $K = \frac{π + \sqrt{π^{2} - 4}}{2} = 2.78$ , implying $Δ = σ$ . This choice stands for the Rayleigh-like criterion and additionally fixes $A = 1 / \sqrt{π}$ , as in the case of the rectangle function above. Obviously other suitable conditions or values of K can be chosen giving other characteristics to the superposition of Gaussians but the Rayleigh-like criterion is the most straightforward way to implement.

In the next section we present some applications of the superposition of Gaussians.

3. Amplitude distribution in a vicinity of a focal region

To illustrate the use of the superposition of Gaussians we have chosen some examples on the calculations of the amplitude distribution near the vicinity of the focal region of a focused wave emerging from a circular aperture. We have chosen this problem due to the complexity of the calculations. As indicated above, in this problem Lommel integrals and sampling expansions are commonly used.

Let us consider the physical situation depicted in Fig. 1 , where a spherical wave with amplitude distribution $g (x, y)$ is incident on a circular aperture of radius a, centered at the origin of a plane $(x, y)$ .

Fig. 1 Schematic diagram of the physical situation.

Download Full Size | PDF

In writing the incident amplitude distribution as an arbitrary smooth function $g (x, y)$ will allow us to consider different situations as is the case of apodization, tilted and/or defocused beams as well as different types of aberrations. Thus, at the incident plane with coordinates $(x, y)$ the amplitude distribution can be written as Eq. (15),

Ψ_{I N} (x, y) = c i r c (r, a) g (x, y) \exp [- i \frac{π}{λ f} (x^{2} + y^{2})],

where,

r = \sqrt{x^{2} + y^{2}}

, λ is the wavelength of the illuminating beam and

c i r c (r, a)

represents the usual circular aperture function of radius a centered at the origin.

As indicated above, we will use the Fresnel diffraction integral to calculate the amplitude distribution at the plane $(ξ, η)$ situated at a distance z of the incident plane. In order to compare our results with [1,2] we will limit ourselves to the case $| (f - z) / f | < < 1$ .Thus, the amplitude distribution in the coordinate plane $(ξ, η)$ can be written as Eq. (16),

ψ (ξ, μ) = \frac{\exp (i \frac{2 π}{λ} z)}{i λ z} {\int_{- \infty}^{\infty} \int_{- \infty}^{\infty} ψ}_{I N} (x, y) \exp [i \frac{π}{λ z} {{(x - ξ)}^{2} + {(y - η)}^{2}}] d x d y .

Due to the presence of the circular aperture, we will use cylindrical variables Eq. (17),

x = r \cos (θ), y = r \sin (θ); ξ = ρ \cos (φ), η = ρ \sin (φ),

where φ and θ are defined as usual (Fig. 1). Thus, Eq. (16) can be rewritten as,

\begin{matrix} Ψ (ρ, φ) = a^{2} \frac{\exp (i \frac{2 π}{λ} (f + Δ z))}{i λ (f + Δ z)} \exp (i \frac{π}{λ (f + Δ z)} ρ^{2}) \times \\ \int_{0}^{1} \int_{0}^{2 π} g (s, θ) \exp [- i \frac{π a^{2} Δ z}{λ f (f + Δ z)} s^{2}] \exp [- i \frac{2 π a}{λ (f + Δ z)} s ρ \cos (θ - φ)] s d s d θ . \end{matrix}

In Eq. (18) we have introduced $Δ z = z - f$ and the dimensionless variable $s = r / a$ .

We will approximate $f + Δ z \approx f$ in all the denominators of Eq. (18), this approximation will result in the well known Debye integral. This approximation is not a limitation of our technique and is introduced only for comparative purposes.

By introducing the following dimensionless variables Eq. (19),

u = \frac{2 π a^{2} Δ z}{λ f^{2}}; v = \frac{2 π a ρ}{λ f} .

Equation (18) can now be approximated as,

\begin{matrix} Ψ (ρ, φ) = a^{2} \frac{\exp (i \frac{2 π}{λ} (f + Δ z))}{i λ f} \exp (i \frac{π}{λ f} ρ^{2}) \times \\ \int_{0}^{1} \int_{0}^{2 π} g (s, θ) \exp [- i \frac{u s^{2}}{2}] \exp [- i v s \cos (θ - φ)] s d s d θ . \end{matrix}

Equation (20) represents the well known amplitude distribution previously reported in the literature [1,2] and, the methods above described (Lommel integrals and sampling expansions) are commonly used for its evaluation

Now the superposition of Gaussians will be applied to represent one of the integrands in Eq. (20). In this case, let

s \int_{0}^{2 π} g (s, θ) \exp [- i v s \cos (θ - φ)] d θ = \sum_{n = 0}^{N G} A_{n} \exp (- \frac{{(s - s_{n})}^{2}}{σ^{2}}),

where

N G + 1

is the total number of Gaussian functions needed to represent the terms at the left of Eq. (21). The number of Gaussians has to be chosen by trial and error until a desired accuracy is obtained.

In general, the terms in the Fresnel integral aside of the quadratic phase are slow varying functions. In our example, this can be confirmed by plotting the term at the left of the equal sign of Eq. (21) as a function of s. Figure 8 bellow shows a plot of one of these functions. We have found that 50 terms are sufficient to make the maximum error less than 2/1000. For an even better accuracy we are using 60 Gaussians in our calculations. Clearly, a much greater number of pixels would be required for sampling the quadratic phase in Eq. (20).

Fig. 8 Amplitude of Eq. (21) after performing the integral compared vs. its representation by superposition of 60 Gaussians for one of the line scans. The solid line corresponds to the integral, the circles to the superposition of Gaussians.

Download Full Size | PDF

Now the numerical part of our technique is performed. It is necessary to calculate numerically the integral involved in Eq. (21). This task is performed with precision by using the numerical integration technique of three point Simpson’s rule of parabolic segments [12].

Once the numerical integration is performed, we proceed by substituting the Gaussian superposition representation given by Eq. (21) into Eq. (20) to obtain,

\begin{matrix} Ψ (ρ, φ) = a^{2} \frac{\exp (i \frac{2 π}{λ} (f + Δ z))}{i λ f} \exp (i \frac{π}{λ f} ρ^{2}) \times \\ \sum_{n = 0}^{N G} A_{n} \int_{0}^{1} \exp [- i \frac{u s^{2}}{2}] \exp (- \frac{{(s - s_{n})}^{2}}{σ^{2}}) d s . \end{matrix}

As we are interested only in normalized intensity distributions, without loss of generality, the terms at the left of the summation symbol, for brevity, will not be considered; if desired these terms can be maintained to obtain the overall amplitude distribution. Thus, aside of the mentioned terms, Eq. (22) can be written as,

\begin{matrix} Ψ (ρ, φ) = \sum_{n = 0}^{N G} A_{n} \exp (- \frac{s_{n}^{2}}{σ^{2}}) \exp (\frac{2 s_{n}^{2}}{σ^{2} (i σ^{2} u + 2)}) \times \\ \int_{0}^{1} \exp [- (\frac{i σ^{2} u + 2}{2 σ^{2}}) {(s - \frac{2 s_{n}}{i σ^{2} u + 2})}^{2}] d s . \end{matrix}

To perform the integral in Eq. (23), we introduce the change of variable Eq. (24),

p = s - \frac{2 s_{n}}{i σ^{2} u + 2} .

Equation (23) can now be written as Eq. (25),

Ψ (ρ, φ) = \sum_{n = 0}^{N G} A_{n} \exp (- \frac{(i σ^{2} u + 2) s_{n}^{2} + 2 s_{n}^{2}}{σ^{2} (i σ^{2} u + 2)}) \int_{- (\frac{2 s_{n}}{i σ^{2} u + 2})}^{1 - \frac{2 s_{n}}{i σ^{2} u + 2}} \exp [- \frac{i σ^{2} u + 2}{2 σ^{2}} p^{2}] d p .

Next, by introducing a new change of variable Eq. (26),

q = \sqrt{\frac{i σ^{2} u + 2}{2 σ^{2}}} p,

Equation (25) becomes,

Ψ (ρ, φ) = \sum_{n = 0}^{N G} A_{n} \exp (- \frac{i u s_{n}^{2}}{i σ^{2} u + 2}) \sqrt{\frac{2 σ^{2}}{i σ^{2} u + 2}} \int_{- (\frac{2 s_{n}}{i σ^{2} u + 2}) \sqrt{\frac{i σ^{2} u + 2}{2 σ^{2}}}}^{1 - \frac{2 s_{n}}{i σ^{2} u + 2} \sqrt{\frac{i σ^{2} u + 2}{2 σ^{2}}}} \exp [- q^{2}] d q .

Finally, using well known properties of the error function, $e r f (x)$ , Eq. (27) is written as,

\begin{matrix} Ψ (ρ, φ) = \sqrt{\frac{π}{2}} \sqrt{\frac{2 σ^{2}}{i σ^{2} u + 2}} \sum_{n = 0}^{N G} A_{n} \exp (- \frac{i u s_{n}^{2}}{(i σ^{2} u + 2)}) \times \\ {e r f [(1 - \frac{2 s_{n}}{i σ^{2} u + 2}) \sqrt{\frac{i σ^{2} u + 2}{2 σ^{2}}}] + e r f [\frac{2 s_{n}}{i σ^{2} u + 2} \sqrt{\frac{i σ^{2} u + 2}{2 σ^{2}}}]} . \end{matrix}

Equation (28) represents the final step of the method. It will be noticed that the result is an analytical (closed) form. For our particular problem, Eq. (28) represents the main equation to calculate the amplitude distribution at any arbitrary meridional plane.

It will be noticed in Eq. (28) the existence of two error functions with complex arguments. As most programs do not accept complex arguments for calculating the error function, for completeness of our presentation, we will use the following Eqs. (29), (30a) and (30b) given in [13],

e r f (x + i y) = e r f (x) + \frac{\exp (- x^{2})}{2 π x} [\begin{array}{l} 1 - \cos (2 x y) + i \sin (2 x y) + \\ \frac{2}{π} \exp (- x^{2}) \sum_{n = 1}^{\infty} \frac{\exp (- \frac{n^{2}}{4})}{n^{2} + 4 x^{2}} [f (n, x, y) + i g (n, x, y)] \end{array}],

where,

f (n, x, y) = 2 x - 2 x \cosh (n y) \cos (2 x y) + n \sinh (n y) \sin (2 x y),

g (n, x, y) = 2 x \cosh (n y) \sin (2 x y) + n \sinh (n y) \cos (2 x y) .

The above set of equations to calculate the error function of a complex number is very convenient for the present application as it avoids using factorials which are difficult to calculate and diverges for large numbers. Additionally, the summation over n in the above equations can be truncated to moderate values due to the properties of the error function. To compare the results obtained with the above equations, and only for comparative purposes, we verified the accuracy of the results with those obtained with Wolfram’s Mathematica®.

In the following section, we present numerical calculations of some interesting cases. For brevity we will limit ourselves to only some descriptive examples.

3.1 Numerical examples

For our examples $g (s, θ)$ in Eq. (21) will be written as Eq. (31),

g (s, θ) = \exp [i a s \cos (θ) + i b s \sin (θ) + i c 1 s^{3} \cos (θ) + i c 2 s^{3} \sin (θ)],

where a and b stand for a tilted beam in the x and y axes respectively and

c 1

and

c 2

for coma in x and y respectively. Note however that our method may be applied to a far more general aberrated and/or apodizated wave-front. Our results shown below are normalized in intensity.

We start with the simpler case of a spherical aberration-free wave. The corresponding parameters are $a = b = c 1 = c 2 = 0$ . The intensity distribution obtained with the superposition of Gaussians is shown in Fig. 2(A) . The accuracy of the proposed technique may be verified by comparing Fig. 2(A) with the corresponding figures reported in [1,2,6] obtained by other techniques. We have calculated the same distribution using Lommel integrals as done in [1,2]. In Fig. 2(B) we present contours resulting from calculating the absolute value of the difference between both techniques. From Fig. 2(B), it can be appreciated the accuracy of the Gaussian sampling technique.

Fig. 2 (A). Isophotes for an aberration-free focusing wave, $φ = 0; a = b = c 1 = c 2 = 0.$ (B). Absolute value of the difference of the isophotes obtained by Gaussian sampling and by Lommel integrals.

Download Full Size | PDF

As Fig. 2(A) is plotted in adimensional units, its application to a particular case requires Eq. (19). As an example of the use of Fig. 2(A), we notice that the airy disk exhibits its first minimum at ( $u = 0, v = 3.77$ ). Thus, for a lens with an aperture of 2.5 cm, a focal length of 5 cm and an illuminating wavelength of 0.6μm, using Eq. (19), one obtains an Airy disk diameter of approximately 2.88 μm.

Figure 3 shows isophotes due to a tilted beam in the x axis, with corresponding values, $a = π; b = c 1 = c 2 = 0$ . The meridional plane is placed at $φ = 0$ .

Fig. 3 Isophotes for a tilted focusing wave, $φ = 0, a = π; b = c 1 = c 2 = 0.$

Download Full Size | PDF

The next simulation, Fig. 4 , corresponds to the above tilted wave, but now the meridional plane is placed at $φ = π / 2$ .

Fig. 4 Isophotes for a tilted focusing wave. $φ = π / 2, a = π; b = c 1 = c 2 = 0.$

Download Full Size | PDF

The following case, Fig. 5 , consists of a wave with coma in the x axis. The corresponding values are, $φ = 0, c 1 = π; a = b = c 2 = 0.$

Fig. 5 Isophotes for a focused wave with coma, $φ = 0, c 1 = π; a = b = c 2 = 0.$

Download Full Size | PDF

Finally, the two following simulations, Figs. 6 and 7 , correspond to the above case but now for meridional planes rotated at $φ = π / 4$ and $φ = π / 2$ respectively. These values were chosen arbitrarily only for illustrative purposes as basically any value can be chosen.

Fig. 6 Isophotes for a focused wave with coma, $φ = π / 4; c 1 = π; a = b = c 2 = 0$ .

Download Full Size | PDF

Fig. 7 Isophotes for a focused wave with coma, $φ = π / 2, c 1 = π; a = b = c 2 = 0$ .

Download Full Size | PDF

As for our knowledge, the plots given in Figs. 4–7 have not yet been reported; this confirms the merit of our technique.

The simulations were performed in a dual-core PC at 3.00 GHz. Each simulation consisted of a matrix of 60x60 pixels, and took approximately 4 minutes.

Figure 8 shows one of the slow-varying functions given by Eq. (21), for the aberration-free case and its corresponding superposition of Gaussians for one arbitrary line scan. Figure 9 shows the absolute value of the difference of both plots.

Fig. 9 Plot of the absolute value of the difference between the plots of Fig. 8.

Download Full Size | PDF

It will be noticed from Fig. 9 that the maximum error is less than 2/1000. For the simulations presented above, this error can be considered negligible and it can be lowered by increasing the number of Gaussians.

Finally, in order to illustrate the accuracy and versatility of the technique, we present an example on the calculation of the intensity distribution due to a one-dimensional sinusoidal transmittance grating with a period of 1.0 μm and width 20 μm. The grating is illuminated by a plane wave with wavelength equal to 0.5 μm. The plane of detection is placed at a distance of 80 μm. In Fig. 10 , the transmittance of the grating is plotted in a solid line. The red circles correspond to the plot of the transmittance represented by a superposition of 1000 Gaussians. For the simulation we used 1000 sampling pixels. Figure 11 depicts the normalized intensity distribution at the plane of detection. The processing time was 20 seconds.

Fig. 10 Sinusoidal transmittance grating (solid line), superposition of Gaussians (red circles).

Download Full Size | PDF

Fig. 11 Intensity distribution at the plane of detection due to the transmittance grating.

Download Full Size | PDF

To show an additional feature of our method, we performed again the calculation of the intensity distribution for the grating in the conditions described above, but this time with a field of view around the first diffraction order. The same number of pixels and Gaussians were used. Figure 12 shows the resulted intensity normalized with respect to Fig. 11.

Fig. 12 Intensity distribution at the plane of detection due to the transmittance grating in a region around the first diffraction order.

Download Full Size | PDF

It can be appreciated in Fig. 12 the structure of the diffraction pattern with high precision. This can be useful for applications where fine details of the diffraction pattern play an important role.

Before finishing this report it is worth mentioning that, as computer systems are becoming larger and faster, FFT techniques benefit by using larger sets of sampling pixels and speeding up the processing time [7]. However, by using FFT one is essentially calculating a convolution between pixels which may result in inaccuracies in the calculation of the complex phase. This has to be considered, especially in interferometry and holography, where the complex phase has to be calculated with high accuracy. The herein superposition of Gaussians avoids, in principle, this inherent convolution, as the integral that results can be calculated as a closed-form expression.

4. Conclusions

We have presented a technique for calculating the amplitude distribution of a propagated wave by means of the Fresnel diffraction integral by representing a given complex function as a finite superposition of complex Gaussians. We have shown that, once the superposition of Gaussians is obtained, it is possible to calculate analytically the amplitude distribution of the propagated field. We have illustrated the use and accuracy of the technique by calculating the intensity contours in the vicinity of the focal region of a focused, aberrated wave-front emerging from a circular aperture. Additionally, the technique has probed to be versatile and accurate and applicable to diffraction calculations when band-limited functions are present in the related integrals.

References and links

1. M. Born, and E. Wolf, Principles of Optics 7th ed. (Pergamon, 1980), Chap. 8.

2. E. H. Linfoot, Recent Advances in Optics (Clarendon Press, 1955), Chap. 3.

3. J. E. A. Landgrave and L. R. Berriel-Valdos, “Sampling expansions for three-dimensional light amplitude distribution in the vicinity of an axial image point,” J. Opt. Soc. Am. A 14(11), 2962–2976 (1997). [CrossRef]

4. J. J. Stamnes and H. Heier, “Scalar and electromagnetic diffraction point-spread functions,” Appl. Opt. 37(17), 3612–3622 (1998). [CrossRef]

5. Y. Li, “Expansions for irradiance distribution near the focus in systems of different Fresnel numbers,” J. Opt. Soc. Am. A 23(3), 730–740 (2006). [CrossRef]

6. J. J. Stamnes, “Waves in focal regions,” The Adam Hilger series on Optics and Optoelectronics, (1986).

7. M. Sypek, “Ligth propagation in the Fresnel region. new numerical approach,” Opt. Commun. 116(1-3), 43–48 (1995). [CrossRef]

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

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

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

11. H. P. Hsu, Fourier Analysis (Simon & Schuster, 1970), Chap. 9.

12. R. W. Southworth, and S. L. Deleeuw, Digital and Computation and Numerical Methods (McGraw-Hill, 1965) Chap. 9.

13. M. Abramovitz, and I. Stegun, Handbook of Mathematical Functions in Applied Mathematics series-55, 299 (1972).

A technique for calculating the amplitude distribution of propagated fields by Gaussian sampling

Abstract

1. Introduction

2. Analytical description

3. Amplitude distribution in a vicinity of a focal region

3.1 Numerical examples

4. Conclusions

References and links

Cited By

Figures (12)

Equations (32)

Optics Express