 AEP 4830 HW2 Numerical Integration

•  Yen-Lin Chen

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\cite{Saunders1980}. 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\cite{Recipe}. 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). The $$I/I_{0}$$ ratio was plotted from $$x_{0}=-1\mu m$$ to $$x_{0}=4\mu m$$ by sampling 200 equally-spaced points in between.

Physical Intepretation

We can test if the diffraction pattern makes sense or not by checking the extreme values of $$x_{0}$$. First of all, at the point $$x_{0}=-\infty$$ where all the electromagnetic waves are blocked by the object, the light intensity will be zero. The plot tells us that just $$1\mu m$$ to the left of the edge, there is little diffracted light. On the other hand, for the point $$x_{0}=\infty$$, the presence of the object makes no difference on the light intensity so the ratio will be $$1$$. It is the interval near the edge of the object where Fresnel Diffraction comes in to play.