%==========================================================================
% Sngl_Slit_Diffn_Thy_1D.m
%
% Single Slit Diffraction - simplest theory - far-field/plane-wave approx!
% Sound waves assumed to be propagating in free air/great wide-open!
%
%==========================================================================
%
% Written by Prof. Steven Errede  Last Updated: Feb. 7, 2011 11:15 hr
%
%==========================================================================
close all;
clear all;

single  thtr(1800);
single  thtd(1800);
single Itot1(1800);
single  SIL1(1800);

single  yscr(2000);
single Itot2(2000);
single  SIL2(2000);

% Specify 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

% Specify slit width (m):
Wdth   =    0.1;    % 0.1; 1.0; slit width (m) n.b.   Wdth << Lobs

%===================================
% Calculate Itot vs. theta:
%===================================
Thetad =  -90.0;    % angle theta of observer in degrees
dTheta =    0.1;    % step angle in degrees
for i = 1:1800;
     thtd(i) = Thetad;            % angle theta of observer in degrees
     Thetar  = (pi/180.0)*Thetad; % angle theta of observer in radians
     thtr(i) = Thetar;
     
     delta   = ((pi*Wdth*sin(Thetar))/lambda); % phase (radians)
     
    Itot1(i) = Io*(sin(delta)/delta)^2; % total intensity (Watts/m^2)
     SIL1(i) = 10.0*log10(Itot1(i)/Ir); % Sound Intensity Level (dB)
     
     Thetad  = Thetad + dTheta;   % increment angle for next calculation
end

%===================================
% Calculate Itot vs. yscreen:
%===================================
 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
     
     delta   = ((pi*Wdth*sin(Thetar))/lambda); % phase (radians)
     
    Itot2(i) = Io*(sin(delta)/delta)^2; % total intensity (Watts/m^2)
     SIL2(i) = 10.0*log10(Itot2(i)/Ir); % Sound Intensity Level (dB)
     
     y  = y + dy;   % increment screen position for next calculation
end


figure(01);
plot(thtd,Itot1,'b');
grid on;
xlabel('theta (degrees)');
ylabel('Intensity (Watts/m^{2})');
title('Intensity vs. theta');

figure(02);
semilogy(thtd,Itot1,'b');
grid on;
xlabel('theta (degrees)');
ylabel('Intensity (Watts/m^{2})');
title('Log10 Intensity vs. theta');

figure(03);
plot(thtd,SIL1,'b');
grid on;
xlabel('theta (degrees)');
ylabel('SIL (dB)');
title('SIL vs. theta');

figure(04);
polar(thtr,SIL1,'b');
grid on;
xlabel('theta (degrees)');
ylabel('SIL (dB)');
title('SIL vs. theta');


figure(11);
plot(yscr,Itot2,'b');
grid on;
xlabel('Yscreen (m)');
ylabel('Intensity (Watts/m^{2})');
title('Intensity vs. Yscreen');

figure(12);
semilogy(yscr,Itot2,'b');
grid on;
xlabel('Yscreen (m)');
ylabel('Intensity (Watts/m^{2})');
title('Log10 Intensity vs. Yscreen');

figure(13);
plot(yscr,SIL2,'b');
grid on;
xlabel('Yscreen (m)');
ylabel('SIL (dB)');
title('SIL vs. Yscreen');

%==========================================================================
beep;
fprintf('\n Single Slit Diffraction Calculation Completed !!! \n')
%==========================================================================