%==========================================================================
% Nslit_2D_Intf_Smpl_Thy.m
%
% 2D N-slit interference - simplest theory - far-field/plane-wave approx!
% sound waves assumed to be propagating in free air/great wide-open!
% each "slit" treated simplistically as point source - no spatial extent!
%
%==========================================================================
%
% Written by Prof. Steven Errede  Last Updated: Feb. 7, 2011 11:25 hr
%
%==========================================================================
close all;
clear all;

single     thtxr(1800);
single     thtxd(1800);

single     thtyr(1800);
single     thtyd(1800);

single    Itotx1(1800);
single    Itoty1(1800);
single    Itotxt(1800);
single    Itotyt(1800);
single     SILxt(1800);
single     SILyt(1800);

single      xscr(2000);
single      yscr(2000);

single    Itotx2(2000);
single    Itoty2(2000);
single    Itotxs(2000);
single    Itotys(2000);
single     SILxs(2000);
single     SILys(2000);

single    Itotxy(2000,2000);
%single LogItotxy(2000,2000);
single     SILxy(2000,2000);
 
% Specify # of slits in x and y (n.b. both must be >= 2 !):
Nxslits =  2; % 2; 5;  4;
Nyslits =  2; % 2; 3; 10;

% Specify other needed parameters:
Io     =    1.0;    % intensity from single slit (Watts/m^2)
Ir     = 1.0*10^-12;% reference sound intensity  (Watts/m^2)
Vair   =  343.0;    % speed of propagation of sound - free air (m/s)
freq   = 1000.0;    % frequency (Hz or cps)
lambda = Vair/freq; % wavelength (m)
Lobs   =   10.0;    % perp. distance observer from slits (m) n.b. lambda << Lobs
Dsrc   =    1.0;    % transv. distance between slits     (m) n.b.   Dsrc << Lobs

%=====================================
% Calculate Itot vs. theta-x,y:
%=====================================
Thetad =  -90.0;    % angle theta of observer in degrees
dTheta =    0.1;    %  step angle in degrees
for i = 1:1800;
    thtxd(i) = Thetad;            % angle theta of observer in degrees
    Thetar   = (pi/180.0)*Thetad; % angle theta of observer in radians
    thtxr(i) = Thetar;
     
     deltax  = (2.0*pi/lambda)*Dsrc*sin(Thetar); % resultant Nslit phase diff (radians)
     
   Itotx1(i) = Io*(sin(Nxslits*deltax/2.0)/sin(deltax/2.0))^2;
   Itoty1(i) = Io*(Nyslits)^2;
   Itotxt(i) =  Itotx1(i)*Itoty1(i)/Io;  % total intensity (Watts/m^2)
    SILxt(i) = 10.0*log10(Itotxt(i)/Ir); % Sound Intensity Level (dB)
   
     Thetad  = Thetad + dTheta;   % increment angle for next calculation
end

Thetad =  -90.0;    % angle theta of observer in degrees
dTheta =    0.1;    %  step angle in degrees
for i = 1:1800;
    thtyd(i) = Thetad;            % angle theta of observer in degrees
    Thetar   = (pi/180.0)*Thetad; % angle theta of observer in radians
    thtyr(i) = Thetar;
     
     deltay  = (2.0*pi/lambda)*Dsrc*sin(Thetar); % resultant Nslit phase diff (radians)
     
   Itotx1(i) = Io*(Nxslits)^2;
   Itoty1(i) = Io*(sin(Nyslits*deltay/2.0)/sin(deltay/2.0))^2; % total intensity (Watts/m^2)
   Itotyt(i) =  Itotx1(i)*Itoty1(i)/Io;  % total intensity (Watts/m^2)
    SILyt(i) = 10.0*log10(Itotyt(i)/Ir); % Sound Intensity Level (dB)
     
     Thetad  = Thetad + dTheta;   % increment angle for next calculation
end

figure(01);
plot(thtxd,Itotxt,'b');
grid on;
xlabel('theta-x (degrees)');
ylabel('I(theta-x) (Watts/m^{2})');
title('I(theta-x) vs. theta-x on Y = 0 axis');

figure(02);
semilogy(thtxd,Itotxt,'b');
grid on;
xlabel('theta-x (degrees)');
ylabel('I(theta-x) (Watts/m^{2})');
title('Log I(theta-x) vs. theta-x on Y = 0 axis');

figure(03);
plot(thtyd,Itotyt,'b');
grid on;
xlabel('theta-y (degrees)');
ylabel('I(theta-y) (Watts/m^{2})');
title('I(theta-y) vs. theta-y on X = 0 axis');

figure(04);
semilogy(thtyd,Itotyt,'b');
grid on;
xlabel('theta-y (degrees)');
ylabel('I(theta-y) (Watts/m^{2})');
title('Log I(theta-y) vs. theta-y on X = 0 axis');

figure(05);
plot(thtxd,SILxt,'b');
grid on;
xlabel('theta-x (degrees)');
ylabel('SIL(theta-x) (dB)');
title('SIL(theta-x) vs. theta-x on Y = 0 axis');

figure(06);
plot(thtyd,SILyt,'b');
grid on;
xlabel('theta-y (degrees)');
ylabel('SIL(theta-y) (dB)');
title('SIL(theta-y) vs. theta-y on X = 0 axis');

figure(07);
polar(thtxr,SILxt,'b');
grid on;
xlabel('theta-x (degrees)');
ylabel('SIL(theta-x) (dB)');
title('Polar plot of SIL(theta-x) vs. theta-x on Y = 0 axis');

figure(08);
polar(thtyr,SILyt,'b');
grid on;
xlabel('theta-y (degrees)');
ylabel('SIL(theta-y) (dB)');
title('Polar plot of SIL(theta-y) vs. theta-y on X = 0 axis');


%=====================================
% Calculate Itot vs. x,y-screen:
%=====================================
 x = -50.00; % starting position on screen (m)
dx =   0.05; % step size on screen (m);
for i = 1:2000;
     xscr(i) = x;            % position of observer on perp. screen (m)
     Thetar  = atan(x/Lobs); % angle theta of observer in radians
     
     deltax  = (2.0*pi/lambda)*Dsrc*sin(Thetar); % resultant Nslit phase diff (radians)
     
   Itotx2(i) = Io*(sin(Nxslits*deltax/2.0)/sin(deltax/2.0))^2; % total intensity (Watts/m^2)
   Itoty2(i) = Io*(Nyslits)^2;
   Itotxs(i) =  Itotx2(i)*Itoty2(i)/Io;  % total intensity (Watts/m^2)
    SILxs(i) = 10.0*log10(Itotxs(i)/Ir); % Sound Intensity Level (dB)
     
     x  = x + dx;   % increment screen position for next calculation
end

 y = -50.00; % starting position on screen (m)
dy =   0.05; % step size on screen (m);
for i = 1:2000;
     yscr(i) = y;            % position of observer on perp. screen (m)
     Thetar  = atan(y/Lobs); % angle theta of observer in radians
     
     deltay  = (2.0*pi/lambda)*Dsrc*sin(Thetar); % resultant Nslit phase diff (radians)
     
   Itotx2(i) = Io*(Nxslits)^2;
   Itoty2(i) = Io*(sin(Nyslits*deltay/2.0)/sin(deltay/2.0))^2; % total intensity (Watts/m^2)
   Itotys(i) =  Itotx2(i)*Itoty2(i)/Io;  % total intensity (Watts/m^2)
    SILys(i) = 10.0*log10(Itotys(i)/Ir); % Sound Intensity Level (dB)
     
     y  = y + dy;   % increment screen position for next calculation
end

figure(11);
plot(xscr,Itotxs,'b');
grid on;
xlabel('Xscreen (m)');
ylabel('Intensity (Watts/m^{2})');
title('Intensity vs. Xscreen on Y = 0 axis');

figure(12);
semilogy(xscr,Itotxs,'b');
grid on;
xlabel('Xscreen (m)');
ylabel('Intensity (Watts/m^{2})');
title('Log Intensity vs. Xscreen on Y = 0 axis');

figure(13);
plot(yscr,Itotys,'b');
grid on;
xlabel('Yscreen (m)');
ylabel('Intensity (Watts/m^{2})');
title('Intensity vs. Yscreen on X = 0 axis');

figure(14);
semilogy(yscr,Itotys,'b');
grid on;
xlabel('Yscreen (m)');
ylabel('Intensity (Watts/m^{2})');
title('Log Intensity vs. Yscreen on X = 0 axis');

figure(15);
plot(xscr,SILxs,'b');
grid on;
xlabel('Xscreen (m)');
ylabel('SIL(Xscreen) (dB)');
title('SIL(Xscreen) vs. Xscreen on Y = 0 axis');

figure(16);
plot(yscr,SILys,'b');
grid on;
xlabel('Yscreen (m)');
ylabel('SIL(Yscreen) (dB)');
title('SIL(Yscreen) vs. Yscreen on X = 0 axis');

%==========================================================================
beep;
fprintf('\n Very CPU-intensive I(x,y) vs. x,y calcs - please be patient!! \n')
%==========================================================================

%========================================
% Calculate 2D Itot vs. x,y-screen:
%========================================
 x = -50.00; % starting position on screen (m)
dx =   0.05; % step size on screen (m);
for j = 1:2000;
     xscr(j) = x;            % position of observer on perp. screen (m)
     Thetar  = atan(x/Lobs); % angle theta of observer in radians
     
     deltax  = (2.0*pi/lambda)*Dsrc*sin(Thetar); % resultant Nslit phase diff (radians)
     
   Itotx2(j) = Io*(sin(Nxslits*deltax/2.0)/sin(deltax/2.0))^2; % total intensity (Watts/m^2)
     
     y = -50.00; % starting position on screen (m)
    dy =   0.05; % step size on screen (m);
    for i = 1:2000;
         yscr(i) = y;            % position of observer on perp. screen (m)
         Thetar  = atan(y/Lobs); % angle theta of observer in radians
     
         deltay  = (2.0*pi/lambda)*Dsrc*sin(Thetar); % resultant Nslit phase diff (radians)
     
     Itoty2(i)   = Io*(sin(Nyslits*deltay/2.0)/sin(deltay/2.0))^2; % total intensity (Watts/m^2)
    
     Itotxy(j,i) =    Itotx2(j)*Itoty2(i)/Io;  % total intensity (Watts/m^2)
%  LogItotxy(j,i) = log10(Itotxy(j,i));
      SILxy(j,i) = 10.0*log10(Itotxy(j,i)/Ir); % Sound Intensity Level (dB)
      
         y  = y + dy;   % increment screen position for next calculation
    end
     x  = x + dx;   % increment screen position for next calculation
end


figure(21);
surf(xscr,yscr,Itotxy);
shading interp;
xlabel('Xscreen (m)');
ylabel('Yscreen (m)');
zlabel('Intensity (Watts/m^{2})');
title ('Intensity vs. Xscreen-Yscreen');

%figure(22);
%surf(xscr,yscr,LogItotxy);
%shading interp;
%xlabel('Xscreen (m)');
%ylabel('Yscreen (m)');
%zlabel('Log Intensity (Watts/m^{2})');
%title ('Log Intensity vs. Xscreen-Yscreen');

figure(23);
surf(xscr,yscr,SILxy);
shading interp;
xlabel('Xscreen (m)');
ylabel('Yscreen (m)');
zlabel('SIL (dB)');
title ('SIL vs. Xscreen-Yscreen');

%==========================================================================
beep;
fprintf('\n 2-D Nslit interference simple thy calculation completed !!! \n')
%==========================================================================