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\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}


\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.


  1. M. A. Heald, J. B. Marion. Classical Electromagnetic Radiation, 3rd edit.. Saunders, 1980.

  2. S. A. Teukolsky W. H. Press, B. P. Flannery. Numerical Recipes, The Art of Scientific Computing, 3rd edit.. Camb. Univ. Press, 2007.

[Someone else is editing this]

You are editing this file