Abstract
This article introduces two approximations that allow the evaluation of Fresnel integrals without the need for using numerical algorithms. These equations accomplish the characteristic of being continuous in the same interval as Fresnel. Both expressions have been determined applying the least squares method to suitable expressions. Accuracy of equations improves as increases; as for small values of , it is possible to achieve an absolute error less than . To probe the efficiency of the equations, two case studies are presented, both applied in the optics field. The first case is related to the semi-infinite opaque screen for Fresnel diffraction. In this case study Fresnel integrals are evaluated with the proposed equations to calculate the irradiance distribution and the Cornu spiral for diffraction computations of the Fresnel diffraction; obtained results show a good accuracy. The second case is related to the double aperture problem for Fresnel diffraction.
1. Introduction
Fresnel integrals are two transcendental functions named after Augustin-Jean Fresnel; these are widely used in physics, specially in optics and electromagnetic theory [1, 2]. For instance, in optics, solution of the Fresnel diffraction for a rectangular aperture is calculated, usually, in terms of integrals that cannot be obtained by analytical methods [3–6]. These integrals are known as Fresnel integrals and contain arguments in terms of the rectangular coordinate system (Cartesian coordinate system).
In [7] the iterative Fresnel integral method is introduced for numerical calculation of the Fresnel diffraction patterns from rectangular apertures tilted at an arbitrary angle to the optical axis. The work presented in [8] demonstrated that the propagation rate of truncated, nondiffracting, and accelerating beams could be accurately calculated using a method that represents the beams by a finite superposition of Gaussian [1] to be in turn propagated by means of the Fresnel diffraction integral. In [9], the iterative Fresnel integrals method (IFIM) has been applied, for the first time, to the simulation and generation of the complete near-field Fresnel diffraction images crated by -apertures. The article [2] provides a development for the double Fresnel integral based on adequate variable changes of the problem under study to predict electromagnetic reduction of the double razor edge with reflection on the earth surface. In [10], Fresnel integrals were applied to obtain the derivation of linear FM Pulse Compression Spectra.
For some years proposals to evaluate these integrals have been provided. The work from [11] introduces polynomial approximations to evaluate Fresnel integrals. These approximations are applied in diffraction analysis [12]. An algorithm to evaluate Fresnel integrals based on the continued fractions method is provided in [13]. On the other hand, to evaluate Fresnel integrals for high values of , the solution is resort to the expansion in terms of Bessel functions [14, 15]. Besides, in [16], one of the performed approximations was made by the sum of the first term of the Fresnel integrals asymptotic expansion and an exponential function which approaches infinity at the zero of Fresnel function argument.
Likewise, several different numerical alternatives have been proposed to evaluate Fresnel integrals. [17] presented the Chebyshev approximations for the Fresnel integrals which employ different values of the argument at different evaluation intervals, depending on its magnitude. In [18], a method for spreadsheet computations of Fresnel integrals to six significant figures is introduced; it is based on successive improvements of known relational approximations which are accurate to just three figures.
This work presents two expressions to calculate, analytically, the Fresnel sine integral and Fresnel cosine integral; with our approximations it is possible to evaluate Fresnel integrals in a simple way without the need for using numerical algorithms like Simpson’s rule, among others. These expressions can be evaluated using any calculus software or standard scientific calculator. In this article we assume that exact Fresnel functions are the ones provided by Maple .
The paper is organized as follows. Section 2 provides an introduction to the Fresnel integrals. In Section 3, the least squares method is presented. Deduction of equations that approximate Fresnel integrals is performed in Section 4. Study cases are introduced in Section 5. Section 6 presents the discussion about the obtained results and efficiency of the equations. Finally, conclusions of this work are provided in Section 7.
2. Fresnel Integrals
The Fresnel sine and cosine integrals are represented by and , respectively; these are two integral functions that originate by applying the analysis of Fresnel diffraction phenomena, which is defined by the following integrals:where is a real number and is a real variable. When tends to infinite, [19].
Fresnel integrals ( and ) are odd functions of . They can be extended to the complex numbers domain; thus, it is possible to obtain analytic functions for one variable [20].
Figure 1 shows the Fresnel integrals graph.
The Cornu spiral, also known as clothoid or Euler spiral, is the parametric curve generated by the Fresnel integrals and from (1) and (2). The spiral is a tangent curve to the abscissa axis at the origin and its radius of curvature decreases inversely in proportion to the distance travelled over it (see Figure 2).
In the optics field, the Cornu spiral can be employed to obtain a quantitative image from a diffraction pattern. In this way, when the diffraction is observed from a position close to the diffractor obstacle, it is referred to as near-field diffraction or Fresnel zone. Therefore, in order to perform the diffraction pattern analysis, the semi-periodic Fresnel zones are calculated from the observation point. In consequence, the contribution from each zone is a phasor that adds up. At the limit of the infinitesimal width zones the Cornu spiral is obtained.
The Cornu spiral [1] has the property that its curvature at any point is proportional to the distance along the curvature from the origin. This property allows using it as transition curve for tracing highways and railroads, given the fact that a mobile that has this displacement at constant speed will have a constant angular acceleration [21, 22]. The transition curve has infinite radius at the tangent point to the straight part of the path, radius at the tangency point to the uniform circular curve; therefore, the most common type of curves on highways is straight stretch-clothoid-circular-clothoid-straight stretch [21, 22].
3. Least Squares Approximations
Let (), (), (), , () be the coordinates of a data set, such that , and the adjustment curve be , where () are adjustment constants. The least squares approximations attempt to minimize the sum of the squares from the vertical distances of values to the ideal model and obtain the model function , which minimizes the square error defined by
Therefore, (5) is designed to minimize the error with respect to the sample set. To do this, partial derivatives are defined with respect to every adjustment constant, creating a system of nonlinear equations given as
Solving system (6) for adjustment constants , it is possible to obtain a model that provides the best fit to the data set. This method allows creating a continuous function in the space with a minimized error with respect to the samples [23]. In consequence, it is possible to test various models and select the one that provides the best fit. Figure 3 depicts the methodology. Notice that in Figure 3 the asterisk symbol represents the data set. For every data element a value is given, where . It is possible to obtain the best function that fits the data set, represented by the blue line.
4. Deduction of Approximations
This article proposes two equations which allow the numerical evaluation of the sine and cosine integrals ( and ) given by (1) and (2), respectively.
Notice that in (7) and (8) the first term decreases at the rate of [13]. In consequence, the first term in these equations dictates the behaviour for each one since it acts as the characteristic term of a solution just like the solution of a differential equation [5]. The following terms represent the stationary state, when approaching infinity tends to . Using as adjustment functions (7) and (8) and the NonlinearFit command in Maple , model parameters are adjusted to the data set reducing the error by least squares. In the procedure of obtaining the approximate models, the data set (, ) consists of 300 elements, obtained from the exact solutions and for the interval; therefore, the resultant equations are
In Figures 4 and 5 the absolute error is presented for (9) and (10), respectively. The maximum error for (9) is about when . It can be noticed in Figure 4 that accuracy improves as value increases.
For (10), the maximum absolute error occurs when and is . Likewise, in Figure 5 it can be noticed that the error for (10) tends to decrease as value increases.
In [18] two approximations are proposed to evaluate the sine and cosine integrals; there arewhere and .
As a way to compare the accuracy of our approximations, Figure 6 depicts the behaviour of the absolute error between (9) and (11), which approximate the cosine integral.
(a) Absolute error comparison in the interval
(b) Absolute error comparison in the interval
(c) Absolute error comparison in the interval
Initially the absolute error in (9) is higher than that in (11), but tends to decrease considerably. At , the error is lower than (11) (see Figure 6(a)). By , the absolute error for (11) is about , while that for (9) is lower. For instance, at , the absolute error for (11) has kept constant, but for (9) the absolute error has decreased to (see Figure 6(b)). For higher values the accuracy improves for (11) and (9); for example, at , the absolute error is and , respectively (see Figure 6(c)). From this analysis, it can be concluded that (9) provides better accuracy than (11). Figure 7 depicts the absolute error of approximations for the sine integral (10) and (12). When , absolute error for (10) is higher than that for (12), with a maximum error of (see Figure 7(a)).
(a) Absolute error comparison in the interval
(b) Absolute error comparison in the interval
(c) Absolute error comparison in the interval
When , the accuracy improves notoriously for (10). For instance, in Figure 7(b), notice that at the absolute error is about , while (12) keeps its absolute error at about . The accuracy of (10) continues to improve for larger values. At , the absolute error approaches , while (12) keeps its error magnitude close to (see Figure 7(c)). Nevertheless, although (9) and (10) exhibit an improved accuracy as increases, they just handle one term that describes the behaviour of the Fresnel integrals while (11) and (12) handle two terms.
Figure 8 shows a comparison between approximations, (9) and (10), and exact Fresnel integrals (1) and (2). On the other hand, Figure 9 shows the Cornu graphic generated by the exact equations (9) and (10). It can be noticed that inside the spiral, as value increases, accuracy improves notoriously.
(a) Comparison between (1) and (9)
(b) Comparison between (2) and (10)
5. Study Cases
This section provides two study cases applied to the optics field, specifically, the diffraction phenomenon. The first addresses the semi-infinite opaque screen subject; the second is about the double aperture problem for Fresnel diffraction.
5.1. The Semi-Infinite Opaque Screen for Fresnel Diffraction
For this study case we will consider the semi-infinite opaque screen defined by the plane (see Figure 10), where is the intensity at point in units or and is the intensity of these incident parallel wavefronts. From [1], intensity at is proportional to
For simplicity, we assume that the semi-infinite opaque screen for Fresnel diffraction is illuminated by a plane wave given bywhere is the source, the observation point, and the wavelength of the projected light beam. By substituting (14) in (9) and (10), equations are now in terms of ; thus, it is possible to calculate and . Now, substituting in (13), we obtain
Proposed numerical values are and , , , . To obtain the quantitative values of the intensities from the Cornu spiral, the length from point (,) to points , , , , and is measured, calculated using (9) and (10).
Likewise, Figures 11–14 depict the diagrams of the irradiance distribution and Cornu spiral for a semi-infinite opaque screen, obtained using (9) and (10), against the exact equations for different values of . The origin of the Cornu spiral diagram corresponds to the edge of the semi-infinite opaque screen, point (), where the intensity has decreased a quarter for negative values of ; in numerical terms, we have that , and because half of the waveform is obstructed. The amplitude of the perturbation has been reduced by half and the irradiance drops by a quarter. As already mentioned, this occurs at point of the plane (see Figure 10), which can be located at the corresponding calculated irradiance distribution and the Cornu spiral in Figures 11–14.
(a) The corresponding calculated irradiance distribution with and
(b) The Cornu spiral for a semi-infinite opaque screen with and
(a) The corresponding calculated irradiance distribution with and
(b) The Cornu spiral for a semi-infinite screen with and
(a) The corresponding calculated irradiance distribution with and
(b) The Cornu spiral for a semi-infinite screen with and
(a) The corresponding calculated irradiance distribution with and
(b) The Cornu spiral for a semi-infinite screen with and
Observation points and are on the shade of the plane. The arc cords decrease monotonically within the region; in consequence, the irradiance abruptly declines. In fact, as decreases, the bands move closer to the edge and become thinner [1]. Observation points were calculated using (9) and (10) and compared with the ones obtained by exact equations. These points are depicted in the magnitude diagrams for , , , and (see Figures 11(a), 12(a), 13(a), and 14(a), respectively). The points obtained by the approximate method in these diagrams keep a relationship with the approximated Cornu spiral in Figures 11(b), 12(b), 13(b), and 14(b), respectively.
5.2. The Double Aperture Problem for Fresnel Diffraction
Light diffraction coming through a rectangular aperture is analysed in [1, 25]. This study case considers the double aperture for Fresnel diffraction presented and studied in [24] (see Figure 15). We assume that the light source with wavelength is , at a distance , and the aperture is centred on the Cartesian system ; the origin of the coordinate system is located at the exact centre of the double aperture (see Figure 15). Upper and lower edges of the aperture are located at and , respectively. For this system, the width for each aperture is () and separation between centres is (). Besides, it is considered that diffracted light is observed on a screen located at away from it. Reference [24] can be consulted for further detailed analysis of this study case.
Simulations performed in Matlab for this study case (based on [24]) are presented in Figure 16(a). This work considers the same aperture proportions given in [24], these are , , , , image area , and step size of . Additionally, the illumination wavelength is .
![(a) Computer simulated Fresnel diffraction using approximate equations and iterative Fresnel integrals method [24]](https://static.hindawi.com/articles/mpe/volume-2018/4031793/figures/4031793.fig.0016a.jpg)
(a) Computer simulated Fresnel diffraction using approximate equations and iterative Fresnel integrals method [24]
(b) Computer simulated Fresnel diffraction using approximate equations (9) and (10)
To generate the simulation shown in Figure 16(b), built-in functions of Matlab were replaced by the proposed functions (9) and (10) (see the Appendix). The Pearson correlation coefficient, , is wildly employed in statistical analysis, pattern recognition, and image processing [26, 27]. For the monochromatic digital images, the Person correlation coefficient is defined aswhere and are intensity values of pixel in and image, respectively; and are mean intensity values of first and second image, respectively. The correlation coefficient is when both images are completely identical. For this study case, correlation value is very high between the images generated with a value of .
6. Discussion
This section presents numerical results of the previous study cases. For the first case, the comparison is performed using the Cornus spiral and observation points of the irradiance distribution. We calculated the absolute error difference from to every observation points of the diffraction pattern. Table 1 shows data for the case when and ; the error is and at observation points and , respectively.
Notice how the error decreases on observation points and (see Tables 1–4), because they are internally moving within the spiral since its geometry is contracting and are in function of the diffracted light wavelength. For observation points and , the error has increased to about because they have moved to the right, closer to the origin. This happens for two reasons: the wavelength is getting smaller and we approach the region where the error value is higher for (9) and (10). Nevertheless, in Table 4, the absolute error at has diminished due to the accuracy of our proposed expressions (see Figures 4 and 5). For observation point , the error is zero becuase the intersection with axis is always the same. This error is in good agreement with the exact spiral because (9) and (10) evaluated at are equal to zero.
For the second study case, the Pearson correlation coefficient value was , which indicates that properties of diffraction patterns from the reported simulation in [24] and the ones obtained from our approximations (see Figure 16) are identical. This result shows the accuracy that can be accomplished by using our proposal on applications for Fresnel diffraction.
7. Conclusions
Approximations proposed on this article allow evaluating Fresnel sine and Fresnel cosine integrals with high accuracy. In the first study case, we built the Cornu spiral and the irradiance distribution with good results since the lower wavelength of the diffraction, the higher accuracy on the shading zone, and the zero are kept at the origin of the Cornu diagram. The proposed equations in this work are continuous throughout the range of real numbers and are evaluated at the origin, with an exact value of zero. As for the second study case, simulations showed a correlation close to the unit. Therefore, the simulations of the double aperture problem for Fresnel diffraction produce near identical figures; thus, the simulation in Figure 16(a) presented in [24] keeps an almost perfect correlation against the simulation in Figure 16(b). With the obtained results, the proposed equations can be employed in diverse applications, specially for values of ; as the value of increases, accuracy improves considerably. Applying the experience of approximate Fresnel functions on this paper, as a future work, we will attempt to improve precision of the results. One possible path to follow would be the use of a different technique, besides least squares method, that allows improving precision closer to the origin region; among other techniques to be applied is the power series extender method [28, 29]. This method may be useful because it provides highly accurate approximations at the origin, while maintaining good prediction for large domains.
Appendix
Modified Matlab/Octave Program for the Double Aperture [24]
The original program of Matlab was
publicated on
Computer simulation of Fresnel diffraction
from double rectangular
apertures in one and two dimensions using
the iterative Fresnel
integrals method.
Abedin, Kazi Monowar and Rahman, SM Mujibur.
Optics & Laser Technology, vol.44(2),
pp.394--402.
2012, Elsevier
Program modified by Mario Sandoval H.
New name of funtion Matlab/Octave
function D=dbslitMprueba(amm,bmm,cmm,Wmm,
smm,q0mm,lnm,t)
Replace to-----------------------------
---------
function D=dbslit(amm,bmm,cmm,Wmm,smm,
q0mm,lnm,t)
----------------------------------
---------------
Date of modification
01/12/2017
Using the program dbslit, the values of the
inputs, i.e., aperture
dimension
Arguments
half-widths: a, b and c in [mm],
amm,
bmm,
cmm,
Image area half-width W in [mm],
Wmm,
Step size s in [mm],
smm,
The aperture-image plane distance
q0 in [mm],
q0mm,
The illuminating wavelength l in [nm]
lnm,
and the intensity factor t are input to the
program. The intensity
factor t is used to control the apparent
visual intensity in the
generated image. For most of the images,
a factor of unity is good
enough.
t,
format long
l=lnm*1e-6;
a=amm*sqrt(2/(l*q0mm));
b=bmm*sqrt(2/(l*q0mm));
c=cmm*sqrt(2/(l*q0mm));
W=Wmm*sqrt(2/(l*q0mm));
s=smm*sqrt(2/(l*q0mm));
dimenU=b:s:W+b;
longiU=length(dimenU);
Cu2=zeros(1,longiU);
Su2=zeros(1,longiU);
dimenU=a:s:W+a;
longiU=length(dimenU);
Cu1=zeros(1,longiU);
Su1=zeros(1,longiU);
dimenU=-a:s:W-a;
longiU=length(dimenU);
Cu3=zeros(1,longiU);
Su3=zeros(1,longiU);
dimenU=-b:s:W-b;
longiU=length(dimenU);
Cu4=zeros(1,longiU);
Su4=zeros(1,longiU);
dimenV=c:s:W+c;
longiV=length(dimenV);
Cv2=zeros(1,longiV);
Sv2=zeros(1,longiV);
dimenV=-c:s:W-c;
longiV=length(dimenV);
Cv1=zeros(1,longiV);
Sv1=zeros(1,longiV);
Functions modified
FFRESNELC and FFRESNELS replace to function
inter-build of matlab
FresnelC and FresnelS
You must be write the equations (9) and (10)
on file .M each one.
i=0;
for contador=b:s:W+b
i=i+1;
Cu2(1,i)=FFRESNELC(contador);
Su2(1,i)=FFRESNELS(contador);
end
i=0;
for contador=a:s:W+a
i=i+1;
Cu1(1,i)=FFRESNELC(contador);
Su1(1,i)=FFRESNELS(contador);
end
i=0;
for contador=-a:s:W-a
i=i+1;
Cu3(1,i)=FFRESNELC(contador);
Su3(1,i)=FFRESNELS(contador);
end
i=0;
for contador=-b:s:W-b
i=i+1;
Cu4(1,i)=FFRESNELC(contador);
Su4(1,i)=FFRESNELS(contador);
end
A=(Cu2+Cu3-Cu1-Cu4).∧2+(Su2+Su3-Su1-Su4)
. ∧2;
i=0;
for contador=c:s:W+c
i=i+1;
Cv2(1,i)=FFRESNELC(contador);
Sv2(1,i)=FFRESNELS(contador);
end
i=0;
for contador=-c:s:W-c
i=i+1;
Cv1(1,i)=FFRESNELC(contador);
Sv1(1,i)=FFRESNELS(contador);
end
B=(Cv1-Cv2).∧2+(Sv1-Sv2).∧2;
B=B′;
n=size(B);
B=repmat(B(:,1),1,n);
A=A′;
A=repmat(A(:,1),1,n);
A=(A)′;
D=B*A;
D=t*D/max(max(D));
m=2*fix(((2*W)/s)/2);
E=zeros(m+1,m+1);
for q=1:1:m/2+1
for p=1:1:m/2+1
E(m/2+2-p,q+m/2)=D(p,q);
end
end
for q=m+1:-1:m/2+2
for p=1:1:m/2+1
E(p,-q+2+m)=E(p,q);
end
end
for q=1:1:m+1
for p=1:1:m/2
E(-p+2+m,q)=E(p,q);
end
end
y=-W:s:W;
ymm=y*sqrt((l*q0mm)/2);
imagesc(ymm,ymm,E,[0 1]);colormap(gray);
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The first author is currently doing a postdoctoral stay at the Instituto Nacional de Astrofisica, Optica y Electronica in Puebla, Mexico, and wishes to thank M. Eng. Roberto Ruíz Gomez and M.TN. César Colocia García for their support and to acknowledge the financial support from the Secretary of Public Education of México (SEP) through Grant 2703 E476300.0275652.