Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Code of the Angular Spectrum Method

  1. Jul 16, 2015 #1
    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);
    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.
     
  2. jcsd
  3. Jul 20, 2015 #2

    jedishrfu

    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?
     
  4. Jul 21, 2015 #3
    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.
     
  5. Aug 6, 2015 #4

    kreil

    User Avatar
    Gold Member

    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.
     
  6. Aug 6, 2015 #5
    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!
     
  7. Aug 6, 2015 #6

    kreil

    User Avatar
    Gold Member

    Ok, glad you were able to resolve the problem! Consider posting the solution here for the benefit of others.
     
  8. Aug 6, 2015 #7
    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##.
     
  9. Sep 8, 2015 #8
    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;
    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');
     
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Code of the Angular Spectrum Method
  1. Matlab Code (Replies: 2)

Loading...