# Fraunhofer Diffraction - Multiple Square Apertures

I'm stuck for some days trying to get the diffraction pattern from a binary grating / square period grating / multi-aperture square screen... whatever you want to call it.

Imagine I have a 8 * 5 square apertures grating. Something like this:

The intensity profile from a single aperture is easy to calculate:

% Clean all
clc;
close all;
clear all;

% Light Source
lambda = 550e-9;                    % Wavelength [m]
k = (2*pi)/lambda;                  % Wavenumber [m^-1]
Io = 100;                           % Relative intensity

% Aperture
a = 1e-6;                           % Aperture size X-axis [m]
b = 1e-6;                           % Aperture size Y-axis [m]

% Screen
R = 1e-3;                           % Distance aperture-screen [m]
dsize = 0.5e-2;                     % Square size [m]
dpix = 1000;                        % Dimension [pixel]
X = -dsize:2*dsize/dpix:dsize;      % Screen X-axis
Y = X;                              % Screen Y-axis

% Intensity Profile
alpha = k*a*X/(2*pi*R);                     % Periodic parameter X-axis
beta = k*b*Y/(2*pi*R);                      % Periodic parameter Y-axis
I = Io*(sinc(alpha).^2)'*(sinc(beta).^2);   % Intensity profile

% Print
fig = figure();
imshow(I);


To calculate the same for a 8 * 5 grating, I have to do a convolution of my unit-square aperture by a comb function. I will get an infinite pattern. Then I just have to multiply it by rectangular (sinc) function to get a non-infinite pattern. But I'm struggling with it.

I would like to get something to input in Matlab.

So you have a square aperture convolved with a 2D comb. The Fourier transform of a convolution becomes a product; the Fourier transform of a comb is another comb; and the Fourier transform of a square aperture is a 2D sinc-function. So, you end up with the 2D sinc-function multiplied by a comb. For a finite array you them have to convolve the result with a small 2D sinc-function that comes from the overall square aperture. Where did you get stuck?

• I would like to do a numerical simulation in Matlab. But I don't know how to input my function on it. Matlab can interpret a sinc function but not a comb function. Perhaps I could do it in this way: SUM[rect[ (x-x_ij)/a , (y-y_ij)/a]. ("i" and "j" are the matrix indexs). But it's not working... Nov 2 '20 at 9:58
• Any simulation that you would do in Matlab would imply a finite domain. The implication is that you get everything convolved with since functions. Unless you understand what you see in Matlab, it would be better to do this analytically. Nov 2 '20 at 12:52

Could it be right?

% Clean all
clc;
close all;
clear all;

% Light Source
lambda = 550e-9;                    % Wavelength [m]
k = (2*pi)/lambda;                  % Wavenumber [m^-1]
Io = 100;                           % Relative intensity

% Aperture
a = 1e-6;                           % Aperture size X-axis [m]
b = 1e-6;                           % Aperture size Y-axis [m]

% Screen
R = 1e-3;                           % Distance aperture-screen [m]
dsize = 0.5e-2;                     % Square size [m]
dpix = 1000;                        % Dimension [pixel]
X = -dsize:2*dsize/dpix:dsize;      % Screen X-axis
Y = X;                              % Screen Y-axis

% Intensity Profile
alpha = k*a*X/(2*pi*R);                     % Periodic parameter X-axis
beta = k*b*Y/(2*pi*R);                      % Periodic parameter Y-axis
I = Io*(sinc(alpha).^2)'*(sinc(beta).^2);   % Intensity profile


% Square Aperture on Matrix
I2 = repmat(I,16);                          % Replicate square aperture on matrix


% Total Area
V = linspace(-8008,8008,16016);
M = sinc(V).^2'*sinc(V).^2;                 % Rectangular area

% Final Result
I3 = I2.*M*10^16;


Why do I have some bright central lines?

I know how to calculate de Fraunhofer diffraction for a single square aperture:

If I have a mask with multiple apertures, numerically my transmittance function will become:

I'm not sure if my fourier transform is a sum of sinc functions: