Code of the Angular Spectrum Method

Tags:
1. Jul 16, 2015

ecastro

I constructed my code of the Angular Spectrum Method. However, as the distance between the object and the plane of interest increases, the diffraction pattern never disappears; there is still some sort of a diffraction pattern, and I am expecting that it disappears as distance increases.

Here is the code:

Code (Matlab M):
clear; clc;

layer = zeros(499, 499);
diameter = 100*10^(-6);
wavelength = 500*10^(-9);
altitude = 2*10^(-3);
img_width = 1*10^(-3);
img_length = 1*10^(-3);

lambda = wavelength;                          % wavelength (meters)
z = altitude;                                 % altitude (meters)
k = 2*pi/lambda;                              % wavenumber
phy_x = img_width;                            % physical width (meters)
phy_y = img_length;                           % physical length (meters)
obj_size = size(layer);

% alpha and beta (wavenumber components)
dx = linspace(-phy_x/2, phy_x/2, obj_size(2));
dx = dx(1:length(dx));
dy = linspace(-phy_y/2, phy_y/2, obj_size(1));
dy = dy(1:length(dy));

for i = 1:obj_size(2);
for j = 1:obj_size(1)
if sqrt(dx(i)^2 + dy(j)^2) <= radius;
layer(j, i) = 1;
end;
end;
end;

U0 = fftshift(fft2(layer));

Fs_x = obj_size(2)/img_width;
Fs_y = obj_size(1)/img_length;

dx2 = Fs_x^(-1);
dy2 = Fs_y^(-1);

x2 = dx2*(0:(obj_size(2) - 1))';
y2 = dy2*(0:(obj_size(1) - 1))';

dFx = Fs_x/obj_size(2);
dFy = Fs_y/obj_size(1);

Fx = (-Fs_x/2:dFx:(Fs_x/2 - dFx));
Fy = (-Fs_y/2:dFy:(Fs_y/2 - dFy));

alpha = lambda.*Fx; beta = lambda.*Fy;

% gamma
gamma = zeros(length(beta), length(alpha));
for j = 1:length(beta);
for i = 1:length(alpha);
if (alpha(i)^2 + beta(j)^2) > 1;
gamma(j, i) = 0;
else
gamma(j, i) = sqrt(1 - alpha(i)^2 - beta(j)^2);
end;
end;
end;

U1 = ifft2(fftshift(fft2(layer)).*exp(1i*k.*gamma.*z));
I1 = (1/(16*pi)).*(U1.*conj(U1));
I think my error came from the calculation of the alpha and beta components. I just don't know where...

2. Jul 20, 2015

Staff: Mentor

I adjusted your post as it looked like MATLAB code and I split some lines where you did assignements like somex = some_expression; somey=some_expression;

Have you tried stepping through the computation to see where things go astray?

3. Jul 21, 2015

ecastro

Everything seems to be working fine at very small distances, where diffraction is expected to occur. However, as I increase the distance $z$, the diffraction never disappears.

4. Aug 6, 2015

kreil

Aside: gamma is a built-in MATLAB function, so you probably shouldn't name your variable with the same name as it will interfere with that functionality for you and anyone that runs your program.

The variable z only shows up in one place after it is defined: as an exponent to exp in the calculation of U1. I assume this U1 is the diffraction pattern that you want? (You might want to comment your code better so that it's clear what all the variables are.)

The only other quantity of exp that could change is gamma, which depends on alpha and beta. Upon looking at gamma, I realized that your section of code that sets elements of gamma to zero is "dead" (seemingly not used):

Code (Text):

if (alpha(i)^2 + beta(j)^2) > 1;
gamma(j, i) = 0;

When I ran your program, gamma has all nonzero values:

Code (Text):

nnz(gamma)/numel(gamma)

ans =

1

So there might be something up there, if you were expecting some elements of gamma to be zero. But ultimately, if gamma is unchanged and you vary z, then the exp term gets bigger and bigger. I would expect the exponent to be negative here if you wanted the exp term to disappear as z increases.

5. Aug 6, 2015

ecastro

I forgot that $\texttt{gamma}$ is a built-in function!

Yes, it is the diffraction pattern that I want. Sorry for not putting comments on my code.

The condition is to disregard evanescent waves as they do not contribute to the diffraction pattern. I guess these show-up at certain values of $z$.
Anyway, I just resolved the issue by not using spatial frequency coordinates, thank you!

6. Aug 6, 2015

kreil

Ok, glad you were able to resolve the problem! Consider posting the solution here for the benefit of others.

7. Aug 6, 2015

ecastro

Instead of guessing the spatial frequency values of the Fourier transformed electric field, which is done after $\texttt{U0}$ and before the definition of $\texttt{gamma}$, I consider this another form of the angular spectrum method:

$\mathbf{\hat{E}} \left(k_x, k_y; z\right) = \mathbf{\hat{E}} \left(k_x, k_y; 0\right) e^{\pm i k_z z}$

where $k_z = \frac{k z}{\sqrt{x^2 + y^2 + z^2}}$, and $k = \frac{2\pi}{\lambda}$. I then used the inverse Fourier transform to plot the electric field at $x, y, z$.

8. Sep 8, 2015

Congli

Hi ecastro,

I modified your code, delete some unnecessary parts and add a figure show to visualize the results. And note that instead of
Code (Text):
U1 = ifft2(fftshift(fft2(layer)).*exp(1i*k.*gamma_cust.*z));
you should use
Code (Text):
U1 = ifft2(ifftshift(fftshift(fft2(layer)).*exp(1i*k.*gamma_cust.*z)));
instead, since you fftshift the spectrem first, you need to first iffshift it back then can you impose the fft2 on the scalar field. Although the results seem to be the same because the captured intensity is equals abs(U1.*conj(U1))

The modified code is following:
Code (Text):
clc; clear; close all;

% parameters
layer = zeros(499, 499);
diameter = 100e-6;
lambda = 500e-9;            % wavelength (meters)
k = 2*pi/lambda;            % wavenumber
z = 2e-4;                      % altitude (meters)
phy_x = 1e-3;               % physical width (meters)
phy_y = 1e-3;               % physical length (meters)
obj_size = size(layer);

% generate meshgrid
dx = linspace(-phy_x/2, phy_x/2, obj_size(2));
dx = dx(1:length(dx));
dy = linspace(-phy_y/2, phy_y/2, obj_size(1));
dy = dy(1:length(dy));

% generate a circle aperture
for i = 1:obj_size(2)
for j = 1:obj_size(1)
if sqrt(dx(i)^2 + dy(j)^2) <= radius;
layer(j, i) = 1;
end
end
end

Fs_x = obj_size(2)/phy_x;
Fs_y = obj_size(1)/phy_y;
dx2 = Fs_x^(-1);
dy2 = Fs_y^(-1);
dFx = Fs_x/obj_size(2);
dFy = Fs_y/obj_size(1);
Fx = (-Fs_x/2:dFx:(Fs_x/2 - dFx));
Fy = (-Fs_y/2:dFy:(Fs_y/2 - dFy));

% alpha and beta (wavenumber components)
alpha = lambda.*Fx;
beta = lambda.*Fy;

% gamma_cust
gamma_cust = zeros(length(beta), length(alpha));
for j = 1:length(beta)
for i = 1:length(alpha)
if (alpha(i)^2 + beta(j)^2) > 1
gamma_cust(j, i) = 0;
else
gamma_cust(j, i) = sqrt(1 - alpha(i)^2 - beta(j)^2);
end
end
end

% angular spectrm based formula
U1 = ifft2(ifftshift(fftshift(fft2(layer)).*exp(1i*k.*gamma_cust.*z)));
I1 = (1/(16*pi)).*(U1.*conj(U1));

% show results
figure;     imshow(layer, []);      title('circular aperture');
figure;     imshow(I1, []);         title('angular spectrum based diffraction');