% Project 2 # 2: one-dimensional piston
L=.6;
theta=-pi/2:pi/360:pi/2;      
c=345;
for f=500
    lambda=c/f;
    u=(L/lambda).*sin(theta);   %
    G=sin(pi*u)./(pi*u);          %sinc function  
    G=abs(G); %Take the Absolute Value to have positive numbers
    subplot(2,3,1)
    polar(theta,G,'m-')       %Plot
    title(['Beam Pattern for array of ',num2str(L),'m',' and Frequency = ',num2str(f),'Hz'])
    pause(.6)
end
for f=1000
    lambda=c/f;
    u=(L/lambda).*sin(theta);   %
    G=sin(pi*u)./(pi*u);          %sinc function  
    G=abs(G); %Take the Absolute Value to have positive numbers
    subplot(2,3,2)
    polar(theta,G,'m-')       %Plot
    title(['Frequency = ',num2str(f),'Hz'])
    pause(.6)
end
for f=2000
    lambda=c/f;
    u=(L/lambda).*sin(theta);   %
    G=sin(pi*u)./(pi*u);          %sinc function  
    G=abs(G); %Take the Absolute Value to have positive numbers
    subplot(2,3,3)
    polar(theta,G,'m-')       %Plot
       title(['Frequency = ',num2str(f),'Hz'])
    pause(.6)
end
for f=3000
    lambda=c/f;
    u=(L/lambda).*sin(theta);   %
    G=sin(pi*u)./(pi*u);          %sinc function  
    G=abs(G); %Take the Absolute Value to have positive numbers
    subplot(2,3,4)
    polar(theta,G,'m-')       %Plot
      title(['Frequency = ',num2str(f),'Hz'])
    pause(.6)
end
for f=4000
    lambda=c/f;
    u=(L/lambda).*sin(theta);   %
    G=sin(pi*u)./(pi*u);          %sinc function  
    G=abs(G); %Take the Absolute Value to have positive numbers
    subplot(2,3,5)
    polar(theta,G,'m-')       %Plot
       title(['Frequency = ',num2str(f),'Hz'])
    pause(.6)
end
for f=5000
    lambda=c/f;
    u=(L/lambda).*sin(theta);   %
    G=sin(pi*u)./(pi*u);          %sinc function  
    G=abs(G); %Take the Absolute Value to have positive numbers
    subplot(2,3,6)
    polar(theta,G,'m-')       %Plot
      title(['Frequency = ',num2str(f),'Hz'])
    pause(.6)
end

%The higher the frequency the more pointed the beam patter becomes.  This
%is why our ears have a hard time knowing the source of low frequencies.

%Particle Displacement:

f=2000;
c=345;
k=2*pi*f/c;
L=.0025;
lambda=c/f;
[X,Y] = meshgrid(0:.01:2,-3:.1:3);
R=sqrt((X.^2)+(Y.^2));
sin_thetaN=abs(X)./R;
m=((pi*L)*sin_thetaN)./lambda;
p=(cos(k*R)./R).*(sin(m)./m);
figure(2);
meshc(X,Y,p);
shading interp;
%colormap(pink);
COLORMAP('default')
title(['Particle Displacement at Frequency = ',num2str(f),', Angled View with a countour curtain using meshc'])
figure(3);
surfl(X,Y,p)
title(['Particle Displacement at Frequency = ',num2str(f),', Angled View using surf1'])
figure(4);
surfl(X,Y,p); 
shading interp;
colormap(pink);
title(['Particle Displacement at Frequency = ',num2str(f),', Side View using surf1'])