AEP 4830 HW2 Numerical Integration

Calculation of \(I/I_{0}\) with varying \(n\)

The analytical Fresnel Diffraction Theory gives the light intensity at \(x_{0}\) on a viewing plane with distance \(z\) behind the edge the following form.

\begin{equation} I(x_{0})=\frac{I_{0}}{\lambda z}\left|\int^{x_{0}}_{-\infty}\exp\left(\frac{i\pi x^{2}}{\lambda z}\right)dx\right|^{2}\\ \end{equation}

where \(\lambda\) and \(I_{0}\) are the light wavelength and intensity respectively(Heald 1980). It can be transformed into the following form by substituting \(u_{0}=x_{0}\sqrt{\frac{2}{\lambda z}}\):

\begin{equation} \frac{I(u_{0})}{I_{0}}=\frac{1}{2}\left[(C(u_{0})-C(-\infty))^{2}+(S(u_{0})-S(-\infty))^{2}\right]\\ \end{equation}

where

\begin{equation} C(u_{0})=\int^{u_{0}}_{0}\cos\left(\frac{\pi u^{2}}{2}\right)du\\ S(u_{0})=\int^{u_{0}}_{0}\sin\left(\frac{\pi u^{2}}{2}\right)du\\ \\ \end{equation}

As a result, instead of calculating the complex integral in Eq. (1), we simplify the task by evaluating \(C(u_{0})\) and \(S(u_{0})\).

First, I’ll start with calculating both \(C(u_{0})\) and \(S(u_{0})\) by the trapezoid rule(W. H. Press 2007). The inputs of the program are the number of iterations \(n\) and the integral upper limit \(u_{0}\). For \(m^{th}\) iteration, the interval will be cut into \(2^{m}\) equally spaced trapezoids and the intensity ratio in Eq. (2) as well as the errors in \(C(u_{0})\) and \(S(u_{0})\) will be printed out. The output is shown in the following tables. As \(n\) increases, the errors in both functions are reduced dramatically, giving the convergence of the intensity ratio, shown as bold face numbers.

The output of the program \(u_{0}=0.5\) with the calculation result in bold face. \(I/I_{0}=0.65183492\).
\(n=\) \(I(u_{0}=0.5)\) \(C(u_{0})\) Error \(S(u_{0})\) Error
2 0.653124945 0.0083112401 -0.023331144
4 0.652133288 0.00228455698 -0.00571454848
8 0.651907994 0.000583128139 -0.00142013267
16 0.651853074 0.000146514596 -0.000354486758
32 0.651839431 3.66741661e-005 -8.85873175e-005
64 0.651836026 9.17138215e-006 -2.21446777e-005
128 0.651835175 2.29302301e-006 -5.5360349e-006
256 0.651834963 5.73266843e-007 -1.38400032e-006
512 0.651834909 1.43317404e-007 -3.45999553e-007
1024 0.651834896 3.58293944e-008 -8.64998555e-008
2048 0.651834893 8.95735192e-009 -2.16249618e-008
4096 0.651834892 2.23933799e-009 -5.4062403e-009
8192 0.651834892 5.5983429e-010 -1.35156006e-009
16384 0.651834892 1.39958101e-010 -3.37890119e-010
The output of the program \(u_{0}=3.0\) with the calculation result in bold face. \(I/I_{0}=1.107629\).
\(n=\) \(I(u_{0}=0.5)\) \(C(u_{0})\) Error \(S(u_{0})\) Error
2 0.237694441 -2.1358193 -1.32402515
4 2.00269971 0.720191757 1.23815896
8 0.968263507 0.379218525 -0.910132102
16 1.07600671 0.112968633 -0.00749138798
32 1.09992974 0.0221667746 -0.000187176326
64 1.10571672 0.00526303121 -9.54043379e-006
128 1.10715173 0.0012995844 -5.67577483e-007
256 1.10750975 0.000323903392 -3.50423354e-008
512 1.10759921 8.09140821e-005 -2.18347262e-009
1024 1.10762157 2.02246645e-005 -1.36363032e-010
2048 1.10762716 5.05592519e-006 -8.5204066e-012
4096 1.10762856 1.26396624e-006 -5.32296429e-013
8192 1.10762891 3.1599062e-007 -3.09197112e-014
16384 1.107629 7.89975954e-008 -1.66533454e-015

Plot Intensity Ratio for a Given Interval

It is also interesting to calculate the diffraction intensity distribution along the observation plane. We choose the interval from \(x_{0}=-1\mu m\) to \(x_{0}=+4\mu m\) and \(\lambda=0.5\mu m\). This gives the relationship \(u_{0}=2x_{0}\) and we can use the program in the first section to evaluate various points on the plane. It is also convenient to write a c++ function (header file, .h) for both \(C(u_{0})\) and \(S(u_{0})\) and include them in the main function. I chose the number of iteration \(2^{14}=8192\) since it gave enough significant digits in section 1. The program calculates the intensity ratio \(I/I_{0}\) for 200 equally spaced points in the interval and output a .dat file containing two columns of the data: \(x_{0}\) and \(I/I_{0}\). The data are then plotted with MatLab, shown in Fig. (1).