Homework 7

Problem 1: Goodman 4-4
you may use Mathematica to help with analytic calculations if you wish.

Problem 2: Goodman 4-16

Problem 3

Write a program to do non-paraxial free-space propagation for one transverse dimension, using the concept developed in section 3.10. The idea is very similar to what you did in Homework 5 in t-ω space. Starting from the input field U(x,y), the first step is to calculate the angular spectrum of the input field A((Subscript[f, x]),(Subscript[f, y]))=(Subscript[FT, xy]){U(x,y)}. Multiplication of the result by the propagation phase as shown in eqn 3-73 in Goodman yields the new angular spectrum. Finally, inverse FT gives the field in position space at the target screen. This code will enable you to calculate diffracted fields from any structure, such as diffraction gratings. As with any code, it is best to test it by comparing it with a situation where the solution is known analytically. To make things less complicated, we will do the calculation in one dimension with a test function which is a simple aperture with uniform illumination: rect[x/x0].

Most of the work is in properly handling the sampling, especially since the point of doing this in a nonparaxial way is to be able to calculate wide-angle propagation.

a. sampling rate

Recall that the function of the circ[ ] in eqn 3-73 is to limit propagation to those angular plane waves that are not evanescent. This sets an upper limit to what frequencies we care about, which therefore determines the sampling rate. Given that our input wavelength is 1μm, calculate the minimum period in x-space for sampling.

b. set up grids, sample the input function

In principle, we could Fourier-transform the input, then multiply by the circ[ ] function (which in one-dimension is just a rect[ ] ). Since we'll be using the FFT, we have to avoid aliasing in the spectrum. Before we sample, then, we can smooth the input function in x-space. There are two ways to accomplish this: The most general way is to convolve the input with the FT of the bandwidth limiting rect[ ] function. This is somewhat cumbersome to program, so as a shortcut, we will simply use a high-order super-Gaussian function as a model of a smoothed rectangular aperture function. 

Set up your grid of points in x space, using a power of 2 for the number of points. Also generate the grid for the spatial frequencies. 
Create a function for the input, where the aperture function is exp[-(((x/w))^40)]. 
Create a list for the aperture function and do a list plot. Let w = 200μm.

c. calculate the angular spectrum of the input function

Now use the FFTx[ ] function for the spatial Fourier transform defined in the FFTdemo program to calculate the angular spectrum. Make a plot of the magnitude of the angular spectrum, and ensure that the power of the result goes close to zero near the boundaries.

d. calculate the propagation phase and the angular spectrum in the target plane

Multiply the spectrum by the propagation phase function as shown in Eqn 3-73. A version of the propagation phase function you can use is:

propPhase[fx_,z_]:=Exp[I ((2 π)/λ)((Sqrt[1-(fx^2)(λ^2)])-1)z]UnitStep[1-(fx^2) (λ^2)]

In this version, the -1 in parentheses represents dropping the Exp[ i k z] term from the propagation, which is unimportant in this calculation.

e. calculate field in target plane

Calculate the field for a distance that corresponds to a Fresnel number of 4 (see ch4) and compare it to the analytic Fresnel integral shown below. If you have done it correctly, the overall shape should be very close, though the fine structure is slightly different because the approximations made are different. Change the propagation phase by expanding the square root to first order, then do the comparison again. This time, you should get an exact match to the Fresnel integral (zoom in on your plot to see the details of the oscillations as shown in the book).

NFresnel[a_,z_]:=((((a/2))^2)/(λ z));

α1[w_,x_,z_]:=-(Sqrt[(2/(λ z))])(w+x);
α2[w_,x_,z_]:=(Sqrt[(2/(λ z))])(w-x);
intenFresnelAper[a_,x_,z_]:=(1/2)(((FresnelC[α2[a/2,x,z]]-FresnelC[α1[a/2,x,z]]))^2)+(1/2)(((FresnelS[α2[a/2,x,z]]-FresnelS[α1[a/2,x,z]]))^2);

aTst=20;λ=1;zTst=25;
Print["Fresnel number: ",NFresnel[aTst,zTst]//N];
fresnelPlt=Plot[intenFresnelAper[aTst,x,zTst],{x,-20,20},PlotStyle->Red];

Show[fresnelPlt]