%cas=1  Circular aperture with square obstruction
%cas=2  Separated square apertures 
%cas=3  Periodic slit

cas=3;
lf=0.5e-6*10e-3;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%cas=1
if cas==1
% This file calculates the diffraction pattern of a 
% circular aperture with a square obstruction

figure(1) % The 2 diffraction patterns

R=2.5e-3;   %radius of circle aperture
%W=5e-3/sqrt(2);     %width of square obscuration
W=2e-3;

x=linspace(-6e-6,6e-6,101)+1e-20;
y=linspace(-6e-6,6e-6,101)+1e-20;

[X,Y]=meshgrid(x,y);
r=sqrt(X.^2+Y.^2);

E1=pi*R^2*(2*besselj(1,2*pi*R/lf*r)./(2*pi*R/lf*r));
I1=abs(E1).^2;
I1=I1/max(max(I1));

I1_db=10*log10(I1);
val=-30;
II=find(I1_db<val);
I1_db(II)=val*ones(size(II));

E1x=pi*R^2*(2*besselj(1,2*pi*R/lf*x)./(2*pi*R/lf*x));
I1x=abs(E1x).^2;



subplot(3,2,1)
surf(X*1e6,Y*1e6,I1)
view(2)
shading interp
axis image
title('linear units')

subplot(3,2,2)
surf(X*1e6,Y*1e6,I1_db)
view(2)
shading interp
axis image
title('dB units')


E2=W^2*sinc(W/lf*X).*sinc(W/lf*Y);
I2=abs(E2).^2;
I2=I2/max(max(I2));

I2_db=10*log10(I2);
val=-30;
II=find(I2_db<val);
I2_db(II)=val*ones(size(II));


E2x=W^2*sinc(W/lf*x);
I2x=abs(E2x).^2;


subplot(3,2,3)
surf(X*1e6,Y*1e6,I2)
view(2)
shading interp
axis image

subplot(3,2,4)
surf(X*1e6,Y*1e6,I2_db)
view(2)
shading interp
axis image


%The complete diffraction pattern
Itot=abs(E1-E2).^2;
Itot=Itot/max(max(Itot));

Itot_db=10*log10(Itot);
val=-30;
II=find(Itot_db<val);
Itot_db(II)=val*ones(size(II));

subplot(3,2,5)
surf(X,Y,Itot)
view(2)
shading interp
axis image


subplot(3,2,6)
surf(X,Y,Itot_db)
view(2)
shading interp
axis image

figure(2)  % a cut through the diffraction pattern
subplot(2,1,1)
plot(x,I1x,'b')
hold
plot(x,I2x,'r')

Ix=abs(E1x-E2x).^2;
plot(x,Ix,'g')
hold

subplot(2,1,2)
plot(x,I1x/max(I1x),'b')
hold
plot(x,Ix/max(Ix),'r')

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%case 2
if cas==2
%cas=2
% This file calculates the diffraction pattern of  
% 2 separated square apertures

figure(1) % The diffraction pattern

W=5e-3;     %width of square aperture
D=1*W;        %separation of apertures

del=3e-6;
x=linspace(-del,del,101)+1e-20;
y=linspace(-del,del,101)+1e-20;
[X,Y]=meshgrid(x,y);


E=sinc(W*X/lf).*sinc(W*Y/lf).*cos(pi*D*X./lf);
I=abs(E).^2;
I=I/max(max(I));

I_db=10*log10(I);
val=-30;
II=find(I_db<val);
I_db(II)=val*ones(size(II));

xx=linspace(-del,del,1001)+1e-20;
Ex=sinc(W*xx/lf).*sinc(0).*cos(pi*D*xx./lf);
Ix=abs(Ex).^2;

E1x=sinc(W*xx/lf).*sinc(0);
I1x=abs(E1x).^2;

figure(1)
clf
subplot(1,2,1)
surf(X,Y,I)
view(2)
shading interp
axis image

subplot(1,2,2)
surf(X,Y,I_db)
view(2)
shading interp
axis image


figure(2)
clf
plot(xx,Ix/max(Ix))
hold
plot(xx,I1x/max(I1x),'r')
hold

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Case 3
%Periodic slits
if cas==3
    
    del=0.3e-3; % width of the calculated spectrum
    
    x=linspace(-del,del,501)+1e-20;
    y=linspace(-0.3*del,0.3*del,501)+1e-20;
    [X,Y]=meshgrid(x,y);
    
    xx=linspace(-del,del,1001)+1e-20;
    yy=0;
    
    %N=100;
    period=200e-6;  %separation of slits
    W=0.2*period;   %width of slits
    Wy=1000e-6;    %keep the aperture square
    N=round(Wy/period)  % number of slits
    
    E1=sinc(Wy*Y/lf);
    E1x=sinc(Wy*yy/lf);
    
    Etmp=0;
    Etmpx=0;
    ntmp=round(del*period/lf);
    for nn=-ntmp:ntmp
        E2=sinc(W*nn/period);
        E3=sinc(N*period/lf*(X-nn*lf/period));
        E3x=sinc(N*period/lf*(xx-nn*lf/period));
        
        Etmp=Etmp+E2*E3;
        Etmpx=Etmpx+E2.*E3x;
    end
    Etot=E1.*Etmp;
    Itot=abs(Etot).^2;
    Itot=Itot/max(max(Itot));
    
    Itot_db=10*log10(Itot);
    val=-60
    II=find(Itot_db<val);
    Itot_db(II)=val*ones(size(II));
    
    Etotx=E1x.*Etmpx;
    Itotx=abs(Etotx).^2;
    
    figure(1)
    subplot(2,1,1)
    surf(X,Y,Itot)
    axis([min(x) max(x) min(y) max(y) 0 1])
    view(2)
    shading interp
    axis image
    
    subplot(2,1,2)
    surf(X,Y,Itot_db)
    axis([min(x) max(x) min(y) max(y) min(min(Itot_db)) max(max(Itot_db))])
    view(2)
    shading interp

    
    
    figure(2)
    plot(xx,Itotx./max(Itotx))
    hold
    plot(xx,(sinc(W/lf*xx)).^2,'r:')
    hold
   
    
end