% 1D diffraction; animate the sum of complex amplitudes as a function
% of diffraction angle

a = 1;             % lattice constant in mm
lambda = 0.5e-3;   % wavelength of illumination in mm
ylabel('imaginary') % labels for the axes
xlabel('real')

N=7;   %  use 7 lattice points for this animation
    
position = (-(N-1)/2:(N-1)/2)*a;        % positions of the lattice points
beta = 2*pi*position/lambda;            % scale by wavelength

%          start at zero angle and show firs and second order diffraction
%          make the divisions smaller as N increases
xlabel('real axis')
ylabel('imaginary axis')
for sin_angle = (0 : lambda/N/10/a : 2.5*lambda/a);

%    create a vector that gives the amplitude scattered by each lattice
%    point (rows) and each angle (columns)
amplitude = exp(i*beta'*sin_angle);
%    add the amplitudes and start vector at origin
sum_amplitude = cat(1,0,cumsum(amplitude));
hold off;    % must be a better way to do this but use hold off and on for now
%    calling plot with a complex number automatically creates a plot in the
%    complex plane
plot(sum_amplitude);
hold on;
% draw a thicker red line to show the resultant complex amplitude
resultantx = [0 real(sum_amplitude(N+1))];
resultanty = [0 imag(sum_amplitude(N+1))];
plot(resultantx,resultanty,'r','linewidth', 3);
axis([-N N -N N]);
xlabel('real axis');   % label the axis
ylabel('imaginary axis');
title(['sin\theta = ' num2str(1000*sin_angle,3)])    %  title the plot with the value of N
pause(0.3);
end;