%==========================================================================
% Intf_Diffn_2D_Circ_Aperture_Thy.m
%
% 2-D interference-diffraction associated with N circular apertures
% - simplest theory - far-field/plane-wave approx!
% Sound waves assumed to be propagating in free air/great wide-open!
%
% Ny circular apertures distributed symmetrically along horizontal y-axis
% Nx circular apertures distributed symmetrically along  vertical  x-axis
% We use a right-handed x-y-z coordinate system. 
%==========================================================================
%
% Written by Prof. Steven Errede  Last Updated: Feb. 11, 2011 09:10 hr
%
%==========================================================================
close all;
clear all;

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

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

single   Itotx1(1800);
single    SILx1(1800);

single   Itoty1(1800);
single    SILy1(1800);

single     xscr(2000);
single   Itotx2(2000);
single    SILx2(2000);

single     yscr(2000);
single   Itoty2(2000);
single    SILy2(2000);

single   Itotxy(2000,2000);
%single LgItotxy(2000,2000);
single    SILxy(2000,2000);
 
 
% Specify the # of y,x circular apertures:
Nyapr  = 2; % 1; 2; 1; 2; 2; 4; 4; horizontal
Nxapr  = 4; % 1; 1; 2; 2; 4; 2; 4;  vertical

% Specify the numerical values of 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;    % observer distance           (m) n.b. lambda << Lobs
Dx     =    1.0;    % x-distance between apertures(m) n.b.   Dx   << Lobs
Dy     =    3.0;    % y-distance between apertures(m) n.b.   Dy   << Lobs
Rapr   =    1.0;    % 0.6; 1.0; aperture radius   (m) n.b.   Rapr << Lobs

nu     =      1;    % order of bessel function of 1st kind (see below)

%===================================================================
% Calculate Itot, SIL vs. theta-y along horizontal y-axis @ x = 0:
%===================================================================
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*Dy)/lambda)*sin(Thetar);  % int'f phase (radians)
     rho     = ((2.0*pi*Rapr)/lambda)*sin(Thetar);  % diffn phase (radians)
     
   Itoty1(i) = Io*(Nxapr)^2*(sin(Nyapr*deltay/2.0)/sin(deltay/2.0))^2*((2.0*bessel(nu,rho))/rho)^2; % total intensity (Watts/m^2)
    SILy1(i) = 10.0*log10(Itoty1(i)/Ir); % Sound Intensity Level (dB)
     
     Thetad  = Thetad + dTheta;   % increment angle for next calculation
end

%===================================================================
% Calculate Itot, SIL vs. yscreen along horizontal y-axis @ x = 0:
%===================================================================
 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*Dy)/lambda)*sin(Thetar); % int'f phase (radians)
     rho     = ((2.0*pi*Rapr)/lambda)*sin(Thetar); % diffn phase (radians)
     
   Itoty2(i) = Io*(Nxapr)^2*(sin(Nyapr*deltay/2.0)/sin(deltay/2.0))^2*((2.0*bessel(nu,rho))/rho)^2; % total intensity (Watts/m^2)
    SILy2(i) = 10.0*log10(Itoty2(i)/Ir); % Sound Intensity Level (dB)
     
     y  = y + dy;   % increment screen position for next calculation
end

%======================================================================
% Calculate 1-D Itot, SIL vs. theta-x along horizontal x-axis @ y = 0:
%======================================================================
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*Dx)/lambda)*sin(Thetar);  % int'f phase (radians)
     rho     = ((2.0*pi*Rapr)/lambda)*sin(Thetar);  % diffn phase (radians)
     
   Itotx1(i) = Io*(Nyapr)^2*(sin(Nxapr*deltax/2.0)/sin(deltax/2.0))^2*((2.0*bessel(nu,rho))/rho)^2; % total intensity (Watts/m^2)
    SILx1(i) = 10.0*log10(Itotx1(i)/Ir); % Sound Intensity Level (dB)
     
     Thetad  = Thetad + dTheta;   % increment angle for next calculation
end

%======================================================================
% Calculate 1-D Itot, SIL vs. xscreen along horizontal x-axis @ y = 0:
%======================================================================
 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*Dx)/lambda)*sin(Thetar); % int'f phase (radians)
     rho     = ((2.0*pi*Rapr)/lambda)*sin(Thetar); % diffn phase (radians)
     
   Itotx2(i) = Io*(Nyapr)^2*(sin(Nxapr*deltax/2.0)/sin(deltax/2.0))^2*((2.0*bessel(nu,rho))/rho)^2; % total intensity (Watts/m^2)
    SILx2(i) = 10.0*log10(Itotx2(i)/Ir); % Sound Intensity Level (dB)
     
     x  = x + dx;   % increment screen position for next calculation
end


figure(01);
plot(thtyd,Itoty1,'b');
grid on;
xlabel('theta-y (degrees)');
ylabel('Intensity (Watts/m^{2})');
title('Intensity vs. theta-y at x = 0');

figure(02);
semilogy(thtyd,Itoty1,'b');
grid on;
xlabel('theta-y (degrees)');
ylabel('Intensity (Watts/m^{2})');
title('Log10 Intensity vs. theta-y at x = 0');

figure(03);
plot(thtyd,SILy1,'b');
grid on;
xlabel('theta-y (degrees)');
ylabel('SIL (dB)');
title('SIL vs. theta-y at x = 0');

figure(04);
polar(thtyr,SILy1,'b');
grid on;
xlabel('theta-y (degrees)');
ylabel('SIL (dB)');
title('Polar Plot of SIL vs. theta-y at x = 0');


figure(05);
plot(yscr,Itoty2,'b');
grid on;
xlabel('Yscreen (m)');
ylabel('Intensity (Watts/m^{2})');
title('Intensity vs. Yscreen at x = 0');

figure(06);
semilogy(yscr,Itoty2,'b');
grid on;
xlabel('Yscreen (m)');
ylabel('Intensity (Watts/m^{2})');
title('Log10 Intensity vs. Yscreen at x = 0');

figure(07);
plot(yscr,SILy2,'b');
grid on;
xlabel('Yscreen (m)');
ylabel('SIL (dB)');
title('SIL vs. Yscreen at x = 0');


figure(11);
plot(thtxd,Itotx1,'b');
grid on;
xlabel('theta-x (degrees)');
ylabel('Intensity (Watts/m^{2})');
title('Intensity vs. theta-x at y = 0');

figure(12);
semilogy(thtxd,Itotx1,'b');
grid on;
xlabel('theta-x (degrees)');
ylabel('Intensity (Watts/m^{2})');
title('Log10 Intensity vs. theta-x at y = 0');

figure(13);
plot(thtxd,SILx1,'b');
grid on;
xlabel('theta-x (degrees)');
ylabel('SIL (dB)');
title('SIL vs. theta-x at y = 0');

figure(14);
polar(thtxr,SILx1,'b');
grid on;
xlabel('theta-x (degrees)');
ylabel('SIL (dB)');
title('Polar Plot of SIL vs. theta-x at y = 0');


figure(15);
plot(xscr,Itotx2,'b');
grid on;
xlabel('Xscreen (m)');
ylabel('Intensity (Watts/m^{2})');
title('Intensity vs. Xscreen at y = 0');

figure(16);
semilogy(xscr,Itotx2,'b');
grid on;
xlabel('Xscreen (m)');
ylabel('Intensity (Watts/m^{2})');
title('Log10 Intensity vs. Xscreen at y = 0');

figure(17);
plot(xscr,SILx2,'b');
grid on;
xlabel('Xscreen (m)');
ylabel('SIL (dB)');
title('SIL vs. Xscreen at y = 0');

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

%=========================================
% Calculate 2-D Itot, SIL vs. x,y-screen:
%=========================================
 x = -50.00; % x-starting position on screen (m)
dx =   0.05; % x-step size on screen (m);
for j = 1:2000;
     xscr(j) = x; % x-position of observer on perp. screen (m)

     Thetax  = atan(x/Lobs);    % x-projected angle theta-x in radians
     deltax  = ((2.0*pi*Dx)/lambda)*sin(Thetax); % int'f phase (radians)
     
     y = -50.00; % y-starting position on screen (m)
    dy =   0.05; % y-step size on screen (m);
    for i = 1:2000;
         yscr(i) = y; % y-position of observer on perp. screen (m)
         rscr    = sqrt((xscr(j))^2 + (yscr(i))^2); % radial pos'n on perp. screen (m)
         Thetar  = atan(rscr/Lobs); % angle theta of observer in radians
         Thetay  = atan(y/Lobs);    % y-projected angle theta-y in radians
         deltay  =   ((2.0*pi*Dy)/lambda)*sin(Thetay); % int'f phase (radians)
         rho     = ((2.0*pi*Rapr)/lambda)*sin(Thetar); % diffn phase (radians)
     
     Itotxy(j,i) = Io*(sin(Nxapr*deltax/2.0)/sin(deltax/2.0))^2*(sin(Nyapr*deltay/2.0)/sin(deltay/2.0))^2*((2.0*bessel(nu,rho))/rho)^2; % total intensity (Watts/m^2)
%   LgItotxy(j,i) = log10(Itotxy(j,i));     % log10 of total intensity (Watts/m^2)
      SILxy(j,i) = 10.0*log10(Itotxy(j,i)/Ir);      % Sound Intensity Level (dB)
   
         y  = y + dy;   % increment y-screen position for next calculation
    end
     x  = x + dx;   % increment x-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,LgItotxy);
%shading interp;
%xlabel('Xscreen (m)');
%ylabel('Yscreen (m)');
%zlabel('Log10(Intensity) (Watts/m^{2})');
%title ('Log10(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 Calculation of intf-diffn for 2-D array of Nx, Ny circular apertures completed !!! \n')
%==========================================================================