1. Introduction

Wavefront reconstruction of an optical wave has found wide applications in analyzing and measuring topography, aberrations of optical systems [1,2], adaptive optics for astronomy [3,4]. Because of high frequency of optical wave there is no direct method for detecting the phase of CW monochromatic light. Hence, different techniques have been proposed and implemented to reconstruct wavefront indirectly. All wavefront sensing methods fall into two general categories: interferometric methods and noninterferometric methods [4].

The noninterferometric methods have been developed on the basis of light intensity measurement e.g. the Foucault knife-edge test [4], Shack–Hartmann wavefront sensor [4,5], intensity transport equation [6,7], pyramid wavefront sensors [8], and Gerchberg–Saxton algorithm [9,10] to name a few. In practice, these methods are applied in situations where high accuracy is not required or high speed data processing is crucial and harsh environment conditions are exist.

Interferometry with laser has become one of the most commonly used methods due to availability of solid-state detectors and fast processors. These interferometric techniques require highly coherent reference wave sources, high spatial resolution detectors, and a stable and clean setup. Generally speaking, the precision of interferometric methods is limited by the quality of the reference wavefront. There are also interferometric methods without reference wave for instance; shearing interferometry [11], phase conjugate technique [12], three flats method [13,14] and diffracted wave interferometry [15]. But these methods also have limitations, they add to complexity of the setup, or need extra measurements, and furthermore they require very good optics.

We have already shown that by analyzing the intensity distribution in an interference field in three dimensions one can reconstruct the wave vectors of the two interfering wavefronts [16–18]. Authors have introduced and implemented a new method to reconstruct the wavefront based on nonlocalized two-wave interferometry in which the overall angle of two wavefronts are small and satisfy the paraxial approximation. The algorithm of the method is summarized as following; A) the interference patterns are recorded at two parallel planes at $z_1$ and $z_2$ separated by $\Delta z=z_2-z_1$, and perpendicular to the average propagation direction of two waves, represented by z. B) The phase difference distributions (PDDs) $\Delta {\varphi }_1$and $\Delta {\varphi }_2$on planes $z_1$ and $z_2$ are reconstructed from the fringe patterns by using available fringe analyzing methods. C) Then, the phase difference gradient $\nabla (\Delta \varphi )$and propagation directions of the interfering wavefronts are determined at enough number of points on the interferograms. D) Finally, the propagation directions and consequently the phase distributions of two interfering wavefronts are reconstructed on each element of the fringe pattern.

In present paper, we introduce a new numerical method for reconstruction of two interfering wavefronts. It is a robust and simpler than previous method [16], but takes more computational time. It is assumed that the angle between the two interfering wavefronts is small enough to satisfy the paraxial approximation. The first two steps of the new algorithm are similar to the algorithm used in previous method in which the PDD is utilized instead of intensity distribution on fringe pattern. After these two steps, we assume one of the wavefronts is known, say it’s a plane wave, and construct the other wave by considering the PDD on plane $z=z_{\text{1}}$. Then, we propagate the two wavefronts to the plane $z=z_{\text{2}}$, and calculate the PDD over the plane. We obtain the mean absolute values of the differences between experimental phase distribution and the calculated one at each point of the plane, the COST function. Finally, using an optimization algorithm, such as simulated annealing, the COST function changes the parameters of the assumed wave in an iterative process that minimize the COST function. The wavefront corresponding to the latter value is the sought wavefront.

The article is organized as follows. The algorithm of the method is described in section two. The simulation results for different wavefronts are presented in section three, and conclusion remarks are given in section four.

2. Two wavefronts reconstruction

2.1. Theory and algorithm

In paraxial approximation a wave that is propagating along z direction can be express as:

 

Fig. 1 Schematic representation of the two interfering waves at planes $z=z_1$ and $z=z_2$.

Download Full Size | PDF

The intensity distribution interference patterns of these two waves at two planes $z=z_{\text{1}}$ and $z=z_{\text{2}}$ can be described as,

The PDD can be obtained from fringe pattern by well known fringe pattern analysis methods [1,2], for instance the Fourier or the phase shift fringe pattern analysis methods. To reduce the calculation time, the resolution of experimental phase distributions ($\Delta {\varphi }_{z_1(\exp )}$ and$\Delta {\varphi }_{z_2(\exp )}$) are decreased to N×N meshes, where simulation results have shown that N should be about more than15 and its maximum should be correspond to the lowest required resolution. To start with, we assume one of the wavefronts as a plane wave (${\varphi }_{1z_1}=cte$). Therefore, the second wavefront can be given by,

The two constructed waves propagate distance $\Delta z=z_{\text{2}}-z_{\text{1}}$ in the z direction. So the phase difference of the two wavefronts in the second plane is given by,

We can define a COST function by comparing the measured and calculated distribution of phase differences ($\Delta {\varphi }_{z_2(\text{exp})}$ and$\Delta {\varphi }_{z_2(\text{cal})}$). The value of the function is a parameter related to the accuracy of the estimated wavefronts in the first plane. By Using an optimization algorithm and iterative forward propagation of the waves, we can arrive at a solution that best describes the initial interfering wavefronts. We can summarize the algorithm as follows:

(a) We reduce the resolution of experimental data ($\Delta {\varphi }_{z_1(\exp )}$ and$\Delta {\varphi }_{z_2(\exp )}$) to N×N meshes. (b) We choose the initial conditions for ${\varphi }_{1z_1}$ as a plane wave, ${\varphi }_{1z_1}=cte$. (c) We use all N×N elements of ${\varphi }_{1z_1}$ as optimization parameters (input of COST function). (d) We construct the second wavefront (${\varphi }_{2z_1}$) by adding experimental $\Delta {\varphi }_{z_2}$ to ${\varphi }_{1z_1}$. (e) We propagate the both wavefronts distance$\Delta z$ in the z direction and construct ${\varphi }_{1z_2}$, ${\varphi }_{2z_2}$and calculate the PDD on the second plane. f) We compare the experimental PDD with the calculated one at the second plane, and calculate COST function (Eq. (9). (g) We use an optimization algorithm and the value of the COST function to update the optimization parameters and we repeat the iteration between step (d) and step (g) until the cost function is minimized. (h) End.

2.2. Wavefront propagator

In this method, a wavefront propagator is required to generate the wavefronts in the second plane from the wavefronts in the first plane. There are different free space optical wavefront propagators that are valid under different regimes of operation. The Fourier method that is based on Fresnel-Kirchhoff diffraction integral is used extensively for simulation of wavefront propagation [19–21]. If wave propagation is in geometrical optics regime and the effects of diffraction optics can be neglected, the ray tracing method for wavefront propagation can be used that is simple and need short calculation times [19,22].

Because the physics of our problem satisfies geometrical optics requirements, we use the ray tracing method to propagate wavefront. It is assumed that the wavefront propagates in free space and its phase distribution in the first plane ($z=z_{\text{1}}$) is ${\varphi }_z{}_{{}_{\text{1}}}$. Since a wavefront propagates along its gradient ($\overrightarrow{\nabla {\text{S}}_{{\text{z}}_{\text{1}}}}$) which is proportional to negative of gradient of phase distribution ($\overrightarrow{\nabla {\varphi }_{{\text{z}}_{\text{1}}}}$) [19], if the gradient of phase distribution in the first plane ($z=z_{\text{1}}$) is known, the phase distribution in second observation plane ($z=z_{\text{2}}$) can be calculated. We can obtain the components of $\overrightarrow{\nabla {\varphi }_{{\text{z}}_{\text{1}}}}$ by using spatial derivatives of ${\varphi }_z{}_{{}_{\text{1}}}$ in x and y directions, and eikonal equation [19] that gives the third component of $\overrightarrow{\nabla {\varphi }_{{\text{z}}_{\text{1}}}}$,

The new wave field is known on the primed coordinates, which will generally no longer coincide with the initial N×N grid. Therefore, we transformed the rays’ data to the data grid points by Spline 2D interpolation. This algorithm was implemented in MATLAB, and its accuracy was tested by aberrated waves which were generated by ZEMAX [23] and was propagated the same distance. The waves were generated by reflecting a plane wave from a mirror surface that consisted of 20th order Zernike polynomial aberration, Fig. 2(a) and the propagation distance was 15cm, Fig. 2(b). The size of detector is taken as 4mm×4mm with 400×400 pixels. As shown in Fig. 2(c), peak-to-valley of the phase error between ZEMAX propagation and the method described in section 2.2 is about 0.01rad (optical path difference OPD error is about$\lambda /\text{1}{0}^3$) which is quite accurate.

 

Fig. 2 (a) Phase distribution of a wavefront that is generated by ZEMAX. (b) Propagated wavefront of Fig. 2(a) by 15cm distance. (c) The difference in phase distribution of the propagated wavefront by our algorithm and by ZEMAX.

Download Full Size | PDF

2.3. COST unction

As it was mentioned previously, if the calculated PDD, $\Delta {\varphi }_{z_2(\text{cal})}$, matches the measured PDD, $\Delta {\varphi }_{z_2(\text{exp})}$, then the estimated answer (${\varphi }_{1z_1}$) is the sought wavefront. To compare the results quantitatively a COST functions should be defined (E) which measures the degree of phase matching. For this purpose we compare $\Delta {\varphi }_{z_2\left(\text{cal}\right)}$ to $\Delta {\varphi }_{z_2\left(\text{exp}\right)}$ as follows.

Here e is a 2D array and represents difference between measured and calculated phase distribution at the coordinate $z_2$ (second plane).

3. Simulation results

The spherical and cylindrical wave propagation can be done analytically. A combination of these waves with different radii was used to simulate the wavefronts and their interference patterns in space. The PDDs at both sides of the reference plane were used to find out the individual phase distributions. In this section, the simulation results under two different conditions are presented (Interference of spherical-cylindrical and two cylindrical waves). For both of these simulations, it is assumed that the observation plane is perpendicular to z axis (mean propagation direction) and its size is 5mm×5mm with 500×500 resolution. In order to lowering the calculation time, the resolution of phase distribution is decreased to 50×50 pixels by averaging over 10×10 window. Simulations were performed on a PC with 2.6 GHz clock rate and triple cores processor (AMD Phenom(tm) II X3 710). The optimization time and convergence behavior of the COST function strongly depends on complexity of the wavefronts, grid size of the reconstructed waves and optimization algorithm. Normally with simple optimization algorithm the optimization loop is trapped on local minima and results has less accuracy. We have used the simulated annealing for optimization algorithm which reduce the computation time and also improve the accuracy of the reconstruction.

3.1 Interference of spherical and cylindrical waves

In this case we investigate reconstructing a spherical and a cylindrical wave, whose radii are −1m and 4m, respectively (Figs. 3(a) and 3(b)). Their interference patterns are inscribed on two planes separated by 70cm from each other.

 

Fig. 3 (a) phase distribution of a −1m-radius spherical wave. (b) Phase distribution of a 4m-radius cylindrical wave. (c) Interference patterns of spherical and cylindrical waves at first observation plane (Media 1). (d) Difference of estimated and simulated phase distribution.

Download Full Size | PDF

The difference between the estimated and simulated wavefronts is demonstrated in Fig. 3(d). A 0.7rad peak-to-valley error in PDD and a 6% relative error in the estimated wavefront are achieved in a 52min processing timeframe.

3.2 Interference of two cylindrical waves

In this case we investigate reconstructing two cylindrical waves with −2m and −1m radii (Figs. 4(a) and 4(b)). Their axis is rotated 60° respect to each other in plane of observation. Their interference is inscribed on two planes separated 70cm from each other.

 

Fig. 4 (a) phase distribution of a −2m-radius cylindrical wave. (b) Phase distribution of a −1m-radius cylindrical wave. (c) Interference patterns of two cylindrical waves at first observation plane (Media 2). (d) Difference of estimated and simulated phase distribution.

Download Full Size | PDF

The difference between the estimated and simulated wavefronts is demonstrated in Fig. 4(d). A 0.35rad peak-to-valley error in PDD and a 3% relative error in the estimated wavefront are achieved in a 39min processing timeframe.

4. Conclusion

A technique for reconstruction of two interfering wavefronts has been introduced. The method is based on an iterative algorithm and numerical analysis. The algorithm has been evaluated with some simulated wavefronts such as spherical and cylindrical or two cylindrical wavefronts. The simulation results confirm this method has potentiality to be used for reconstruction of two interfering wavefronts from interference field. Although simulation results show some systematic errors, the attainable accuracy and precision is about $\lambda /\text{1}0$in noiseless condition. The error sources are due to numerical errors in wavefront propagator mostly at corner of the wavefronts and also optimization algorithm which normally trapped in local minima.

For comparison; the algorithm of the reference [16] is an analytical method which for simple wavefronts is much faster than method of current paper. In order to has accurate calculation of PDD with respect to $\Delta z$, there is a limitation for maximum value of $\Delta z$which will affect the attainable accuracy for complex wavefronts. Also, it is difficult to write a automated program for whole process and in the points of fringe patterns at which gradient of PDD in (x, y) plane become zero, it fails to reconstruct the wavefronts. Current method is a numerical method. It is possible to reconstruct more detail of the wavefronts far from boarders by increasing$\Delta z$. At boarders the divergence and convergence of the wave front will limit the accessible information and fails to reconstruct wavefronts accurately. The accuracy depends on grid size and $\Delta z$. It can be easily written in a form of automated program. The method has application to find out the two interfering unknown wavefronts or survive for accuracy of reference wavefront, for instance: two reflected waves from front and back surfaces of a transparent plate, finding the position of two coherent optical point sources for range finding and reflection from lens interfaces for finding the lens parameters [17]. Authors believe the presented method open new possibilities for removing the necessity of good reference wavefront in optical testing and other filed of metrology.