% http://numericalcomputation.blogspot.com/2012/08/matlab-double-slit-interference-and.html
% ********* N Slit Interference and Diffraction combined ***********
% ************************** By Mahesha MG ******************************
% Modified and Kenyon Page references added by Steven Sahyun 5-29-2014
clear;
clf();
thetamax=pi/50;
while(1)
    slits=input('Enter the number of slits: ');
% d is the slit width (p. 137)
% good valuses to use are 1 to 10 um.
    d=input('Enter slit width (in micro meters): ');
    d=d*1e-6;
% a is the slit separation (p. 139)
% good values to use are in the range a=4*d to 20*d
    a=input('Enter slit seperation (in  micro meters): ');
    a=a*1e-6;  
    %    l=input('Enter wavelength (in nm): ');
    lambda=632;
% using red light lambda = 632 nm
    lambda=lambda*1e-9; % convert to meters
    k = 2*pi / lambda; %  k = 2 pi / lambda; see p. 135.
%    s=input('Enter slit to screen distance (in m): ');
    s=0.2; % screen is 20 cm from slits.
    theta=-thetamax:1e-5:thetamax;

% Single slit intensity
% I(theta) = I(0)(sinc(kdsin(theta)/2))^2 [eq. 6.9, p. 137]
% intensity at point y is s*tan(theta)*I(theta)
   y=s*tan(theta);
 
% need to use beta2 rather than beta which is a reserved variable.
% Similar to Equations for alpha and beta are on p. 140. 
% There is a /2 since need to do that in 6.13 so it saves a processing step to do it here.
    alpha2=k*d*(sin(theta))/2;
    beta2=k*a*sin(theta)/2;


% Single slit intensity calculation
    I1=((sin(alpha2))./alpha2).^2;  % Single Slit Interference term
    
% N Slit diffraction term divided by slits^2 to normalize. Eq. 6.13 p. 140
    IN=(I1/slits^2).*(sin(slits.*beta2).^2)./(sin(beta2).^2); 
        
    plot(y,I1,'-.b',y,IN,'r');
    title([num2str(slits) ' slit diffraction']);
    xlabel('Distance in m');
    ylabel('Intensity');
    hold all;  %MatLab function
    ch= input('Press 1 to continue and 0 to exit: ');
if ch == 0
    break;
else
    clf();
end


end