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.
Could you advise me?

 A: 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?
A: 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?
A: I know how to calculate de Fraunhofer diffraction for a single square aperture:
$$\tau_A(\xi,\eta)=\mathrm{rect}\left(\frac{\xi}{l_x}\right)\mathrm{rect}\left(\frac{\eta}{l_y}\right)$$
$$I(x,y,z)=\frac{A^2}{\lambda^2z^2}\mathrm{sinc}^2\left(\frac{l_xx}{\lambda z}\right)\mathrm{sinc}^2\left(\frac{l_yy}{\lambda z}\right)$$
If I have a mask with multiple apertures, numerically my transmittance function will become:
$$\tau_A(\xi,\eta)=\sum_{i,j}\mathrm{rect}\left(\frac{\xi-\xi_{ij}}{a},\frac{\eta-\eta_{ij}}{b}\right)$$
I'm not sure if my fourier transform is a sum of sinc functions:
$$I(x,y,z)=\sum_{i,j}\mathrm{sinc}^2\left(\frac{l_xx-x_{ij}}{\lambda z}\right)\mathrm{sinc}^2\left(\frac{l_yy-y_{ij}}{\lambda z}\right)$$
