1. Introduction

Numerical methods for effective calculation of diffracted fields have been intensively investigated and adopted in various areas of optics, typically in digital holography. Several types of computation methods, which have corresponding mathematical expressions, originate from the Rayleigh-Sommerfeld diffraction formula [1,2] and are selectively employed according to the specific purposes and physical conditions of the application. Two important factors that have to be considered in the field calculation are high accuracy and fast computation, which are inherently incompatible. Hence, most of the methods, such as Fresnel diffraction and the angular spectrum (AS) method, include fast Fourier transform (FFT) algorithms in their calculation procedure [3].

Among these, the AS method is rigorous, yielding results with no approximations. This method calculates field propagation by decomposing the wave field distribution on a plane into plane wave components. The weakness of the conventional AS propagation method is that it only guarantees precision in the near-field region that can be specified by Nyquist sampling theory. The band-limited AS method, however, can be exploited in far-field regions to a certain extent via band limitation of the propagation transfer function [4]. Utilization of a wide calculation window with linear convolution also broadens the applicability of the AS method in both far and near fields [5].

Our study is inspired from the fact that conventional numerical calculation methods are mainly limited to propagation between parallel planes, although there are several techniques for tilted or shifted planes [6–12]. Until recently, there have been a few studies on optical wave propagation from curved surfaces. To mention a few, Esmer et al. presented a field model based algorithm that is developed to achieve exact calculation results with high computational cost [13]. Shimobaba et al. proposed a calculation method for Fresnel diffraction from an arbitrary shaped surface to a flat plane by using non-uniform FFT [14]. The approach introduced by Sahin et al. uses local Gaussian beam decomposition of whole fields diffracted from curved geometries [15,16].

In this paper, we propose a novel diffraction calculation scheme, which can be employed to study the propagation between two arbitrary curved surfaces. To the best of our knowledge, this is the first report on a method, based on the AS method, that implements segmentation of the whole surface region. In Section 2, we present the analytical formulation of the proposed method. In Section 3, using the derived expressions, we show numerical results in order to validate our method. Finally, in Section 4, we provide concluding remarks, accompanied with a consideration of an efficient calculation issue.

2. Theory

In this section, we introduce a theoretical approach for diffraction analysis in which an arbitrary curved surface is involved. Throughout this paper, all the studies are conducted on the two-dimensional field calculation for the sake of simplicity.

2.1 Method 1: staircase approximation with uniform segmentation along the transverse direction

A basic outline of our method is illustrated in Fig. 1. As shown in the figure, the curved surface $z=g(x)$ is divided uniformly along the transverse (x-axis) direction and each subregion is approximated by a piecewise flat surface perpendicular to the longitudinal (z-axis) direction. Here and later, we restrict our analysis to the case that $g(x)$ is an arbitrary single-valued function.

 

Fig. 1 Schematics of the transversely uniform segmentation method for calculation of (a) propagation from an arbitrary curved surface to a flat plane and (b) the opposite case.

Download Full Size | PDF

For wave propagation from an arbitrary curved surface $z=g(x)$ to a plane at $z=z_0$ shown in Fig. 1(a), the angular spectrum of the field confined within the n-th subplane at $z=g(n\Delta x)$ can be expressed as

Note that by using the paraxial approximation in each subregion, the phase compensation term $g_p(n\Delta x,x)$ included in Eq. (1) reveals the effects of the local curvature. This compensation term plays an important role in improving accuracy when the number of segments is not sufficient.

Field propagation via the AS method yields the output field distribution at $z=z_0$ generated from the n-th subregion:

In this equation, we only consider the frequency component $f_x$ that satisfies $f_x^2\le 1/{\lambda }^2$ in order to exclude evanescent wave components. Finally, the field distribution at $z=z_0$ diffracted from the whole surface $z=g(x)$ can be obtained by summing all the field components as

In a similar fashion, as depicted in Fig. 1(b), the analytical expression for the opposite process can be written as follows:

Here, the rectangular window function can be replaced by other functions depending on the specific conditions. For example, we can use the Gaussian window function,

It is interesting to note that the approximate expressions represented by Eq. (5) and Eq. (6) become ideal forms as Δx approaches an infinitesimally small value (i.e., the number of subregions goes to infinity), and the corresponding two expressions are obtained as Eq. (10) and Eq. (12) respectively, as follows:

From Eq. (11) and Eq. (12), it is clear that the two ideal forms of expressions are similar to the conventional mathematical representation of angular spectrum propagation between two parallel planes, except that the curved surface is substituted instead of a flat plane. Consequently, due to this modification, Eq. (11) and Eq. (12) can no longer be considered as a form of a Fourier transform.

Unfortunately, depending on the geometry of some cases, this segmentation method can be inefficient requiring a large number of Fourier transforms. Therefore we can predict that the generation of non-equal sized subregions allows an improvement in precision, without any increase in computational cost.

2.2 Method 2: staircase approximation with uniform stratification along the longitudinal direction

Another possible way to generate subregions is shown in Fig. 2. Instead of creating equal-sized subdivisions along the transverse direction, we create subregions of non-uniform sizes through uniform longitudinal stratification. This way guarantees a certain degree of accuracy by equalizing the longitudinal distances between adjacent subplanes.

 

Fig. 2 Schematics of the longitudinally uniform stratification method for calculation of (a) propagation from an arbitrary curved surface to a flat plane and (b) the opposite case.

Download Full Size | PDF

Following the procedure of the previous subsection, we first establish an analytical formulation for propagation from an arbitrary curved surface to a plane. The angular spectrum of the n-th truncated field on the imaginary plane of $z=n\Delta x$ is

Using Eq. (14), we can find the distribution at $z=z_0$ by summing each field component that is obtained by taking the inverse Fourier transform of the propagated angular spectrum:

Similarly, we can easily derive an expression for the inverted case depicted in Fig. 2(b). The resulting field distribution on a curved surface $g(x)$ is given by

Here again, the phase compensation technique is used in both cases. As expected, when Δz approaches zero, Eq. (15) and Eq. (17) are also transformed into Eq. (10) and Eq. (12), respectively. Note that like before, we can use other window functions, such as the following:

3. Numerical results

To verify the proposed method, numerical simulations are performed with respect to a circularly curved geometry with a radius of curvature $R$, as shown in Fig. 3. In our calculations, the wavelength is $\lambda =500\text{nm}$ and the sampling interval and the number of samplings are $\Delta p=\lambda /2$ and $N_p=2048$, respectively. The size of the calculation window is $l=N_p\Delta p$ and the width of the input geometry is $w=l/2$. We assume that the thickness of the curved surface is infinitesimally thin and the initial amplitude and phase are constant over the entire surface.

 

Fig. 3 Curved geometry and calculation window for the numerical simulation.

Download Full Size | PDF

Figure 4 shows the amplitude of diffracted wave fields from the curved surface presented in Fig. 3. In these calculations, method 1 is used so that each subregion of equal width $w/N_1$ produces a local diffraction field. $N_1$ indicates the number of subregions generated by method 1 and the calculation results are plotted for three values of $N_1$: 5, 15, and 35. In Figs. 4(a)–4(c) and Figs. 4(d)–4(f), $R$ is $500\text{μm}$and $400\text{μm}$, respectively. The diffracted fields are focused at different longitudinal distances as $R$ changes, and these consequences are consistent with our theory.

 

Fig. 4 Calculated amplitude profiles of diffracted fields from the geometry shown in Fig. 3 using method 1 with R = 500μm and (a) N₁ = 5, (b) N₁ = 15, and (c) N₁ = 35. The corresponding results for R = 400μm are shown in (d)–(f).

Download Full Size | PDF

The effects caused by the size of the subregions and the radius of curvature can be found by comparing the figures shown in Fig. 4. Since light propagation becomes nonparaxial as $R$ decreases, the size of the subregions should be reduced to maintain a certain degree of accuracy. In other words, if the radius of curvature is shortened while the size of the subregions remains unchanged, the accuracy would be degraded. This can be observed by comparing the two columns of the results shown in Fig. 4. In addition, it can be clearly seen from the cases shown in Figs. 4(a) and 4(d) that an insufficient number of subregions causes unfavorable gaps between adjoining local field components. These gaps gradually disappear as the size of the subregions decreases.

From the viewpoint of computational complexity, method 1 requires $(N_1+1)$ FFTs to calculate the field distribution at a certain longitudinal distance, i.e. $N_1$ forward FFTs and one inverse FFT for the case shown in Fig. 1(a). The calculation times for the three cross-sectional field distributions of Figs. 4(a)–4(c) are 27.8ms $(N_1=5)$, 35.7ms $(N_1=15)$, and 50.9ms $(N_1=35)$ in our computing environment (Intel Core i5-2450M CPU, 8GB RAM, 64-bit Windows 7, and MATLAB R2013a).

In Fig. 5, we present numerical simulations using method 2. $N_2$ indicates the number of subregions generated by method 2. In this case, subregions have different widths along the x-axis due to the uniform stratification along the z-axis, so that the spacing between adjacent subplanes is constant. By adjusting the spacing between sampled subplanes, calculation errors arising from the staircase approximation can be circumscribed and controlled. As can be seen from the figures, subregions are more densely generated around the regions near both edges of the geometry, and noise field components are reduced compared to the results for method 1.

 

Fig. 5 Calculated amplitude profiles of diffracted fields from the geometry shown in Fig. 3 using method 2 with R = 500μm and (a) N₂ = 5, (b) N₂ = 15, and (c) N₂ = 35. The corresponding results for R = 400μm are shown in (d)–(f).

Download Full Size | PDF

In our numerical examples, the number of FFTs required by method 2 is less than that of method 1 for the same number of subregions since diffracted fields from multiple subregions can be simultaneously calculated using one FFT. The number of FFTs conducted by method 2 is less than or equal to $(N_2+1)$. The computation time of method 2 is not dependent on $N_2$ but depends on the number of stratified subplanes which are uniformly distributed along the z-axis.

To highlight and compare the accuracy of the two methods, we present normalized cross-sectional amplitude distributions at $z=R=400\text{μm}$. Figure 6 shows results with $N=35$. The result calculated from direct integration of Rayleigh-Sommerfeld diffraction formula is also plotted for comparison between the proposed methods and an accurate method. Although it seems that the two distributions are quite similar, signal-to-noise ratio (SNR) analysis [8] demonstrates that method 2 yields more accurate results: SNR values of method 1 and method 2 are 7.29dB and 7.75dB, respectively.

 

Fig. 6 Normalized amplitude distributions calculated from the Rayleigh-Sommerfeld diffraction formula (solid black line), method 1 (dashed blue line), and method 2 (dotted red line) at z = R = 400μm.

Download Full Size | PDF

Lastly, we demonstrate the effects brought about by using a different window function and the phase compensation technique with method 1. Figure 7 shows the amplitude distributions of diffracted fields in the region around the focused spot. Figures 7(a) and 7(b) correspond to Gaussian and rectangular window functions, respectively. In Fig. 7(a), the rectangular window is replaced by the Gaussian window presented in Eq. (9) and other conditions are the same as in Fig. 4(f), i.e., $N_1=35$ and $R=400\text{μm}$. Figure 7(b) is simply a magnified version of Fig. 4(f) and is inserted for the purpose of comparison. In Fig. 7(b), edge-diffracted noise components due to the abrupt change of field at both sides of the rectangular windows are evident. On the other hand, as demonstrated in Fig. 7(a), these unwanted edge diffraction phenomena can be alleviated by using a window function that varies less abruptly. In Fig. 7(c), the wave field without the phase compensation function $g_p(n\Delta x,x)$ is presented for comparison. From this comparison, the crucial role of the phase compensation function is manifested.

 

Fig. 7 Amplitude profiles around the focused spot when (a) Gaussian and (b) rectangular window functions are used, with method 1 and parameters of R = 400μm and N₁ = 35. (c) Wave field without the phase compensation function.

Download Full Size | PDF

4. Conclusion

We have proposed a novel method for diffraction calculation of an arbitrary curved surface based on the AS method. The main feature of our method is that the surface region is divided into multiple subregions and the AS method is applied with regard to each subregion. In addition, a phase compensation technique is utilized to reduce approximation errors that occur due to the finite number of subregions. The proposed method has been verified through numerical investigations on a circularly curved geometry. One can readily exploit our approach in three-dimensional diffraction calculations since the proposed method can be straightforwardly extended. Although the formulations presented are based on the AS method, our approach can be employed in other formulations such as Fresnel diffraction. A possible drawback is the requirement of a large number of FFT operations depending on the geometry. This disadvantage, however, can be overcome by adopting an optimized subregion-generating algorithm and parallel computation techniques [17]. We believe that our approach will provide a fundamental basis for field calculations with curved geometries.