clear all
%   all length units are mm
lambda = 0.000635;
L = 110; 
for i = 1:801
    yp(i) = (i-401)*0.001;
    A(i) = 0;
end
%   a = slit width, d = slit separation, N = number of slits
a = 0.04;
d = 0.125;
N = 2;
%   generate aperature function (must fit into 0.8 mm)
Na = round(1000*a/2);
Nd = round(1000*d/2);
for i = 1:N
    NN = -(N-1) + 2*(i-1);
    for j = 1:2*Na
        jj = 401+NN*Nd-Na+(j-1);
        A(jj) = 1;
    end    
end
%   compute diffraction pattern (Fourier Transform)
for j = 1:833
    y(j) = (j-417)*0.009;
    for i = 1:801
        Acos(i) = (1/(lambda*L))*A(i)*cos(2*pi*yp(i)*y(j)/(lambda*L));
        Asin(i) = (1/(lambda*L))*A(i)*sin(2*pi*yp(i)*y(j)/(lambda*L));
    end
    Er(j) = trapz(yp,Acos);
    Ei(j) = trapz(yp,Asin);
    TheorI(j) = Er(j)^2 + Ei(j)^2;
end
plot(y, TheorI)
figure
%   slits2.dat is intensity data (a=0.04, d=0.125, N=2)
load slits2.dat
yexp = slits2(:,1);
ExpI = slits2(:,2);
plot(yexp,ExpI,'r',y,230*TheorI,'b')
xlim([-4 +4])
hold off