Code of the Angular Spectrum Method

[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:

clear; clc;

layer = zeros(499, 499);
diameter = 100*10^(-6);
radius = diameter/2;
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...

Thank you in advance.

 

[jedishrfu]

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?

 

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

 

[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):

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

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

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.

 

[ecastro]

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.

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

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

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

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):

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

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

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.

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!

 

[kreil]

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

 

[ecastro]

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

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

 

[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

U1 = ifft2(fftshift(fft2(layer)).*exp(1i*k.*gamma_cust.*z));

you should use

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:

clc; clear; close all;

% parameters
layer = zeros(499, 499);
diameter = 100e-6;
radius = diameter/2;
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');