% MatLab script to calculate diffraction by a finite-sized
% one-dimensional lattice
a = 1;             % lattice constant in mm
lambda = 0.5e-3;   % wavelength of illumination in mm
ylabel('intensity') % labels for the axes
xlabel('sin\theta')

for N=3:2:15     %  do a loop starting from 3 lattice points
    
position = (-(N-1)/2:(N-1)/2)*a;        % positions of the lattice points
beta = 2*pi*position/lambda;            % scale by wavelength

%          include 5 peaks in the range of diffraction angles
%          make the divisions smaller as N increases
sin_angle = (-2.5*lambda/a : lambda/N/10/a : 2.5*lambda/a);

%    create a matrix that gives the amplitude scattered by each lattice
%    point (rows) and each angle (columns)
amplitude = exp(i*beta'*sin_angle);

%  add up the contribution from all lattice points to get total
%  amplitude at each angle and take the modulus squared to get the intensity
intensity = abs(sum(amplitude)).^2;

% plot this on a linear x-axis (angles) and log y-axis (intensity)
semilogy(sin_angle, intensity);
% automatic y-axis limits doesn't work here because some values are zero
% set my own limits for the plot
axis([-3*lambda/a 3*lambda/a 0.5 N^2]);
title(['N=' int2str(N)])    %  title the plot with the value of N
pause;
end;
