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