% Project 2 # 3: tapered one-dimensional piston
L=.5;
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)+sin(pi*(u+1))./(2*pi*(u+1))+sin(pi*(u-1))./(2*pi*(u-1));              
    G=abs(G);             %Take the Absolute Value to have positive numbers
    subplot(2,2,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)+sin(pi*(u+1))./(2*pi*(u+1))+sin(pi*(u-1))./(2*pi*(u-1));             
    G=abs(G);             %Take the Absolute Value to have positive numbers
    subplot(2,2,2)
    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)+sin(pi*(u+1))./(2*pi*(u+1))+sin(pi*(u-1))./(2*pi*(u-1));               
    G=abs(G);             %Take the Absolute Value to have positive numbers
    subplot(2,2,3)
    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)+sin(pi*(u+1))./(2*pi*(u+1))+sin(pi*(u-1))./(2*pi*(u-1));  
    G=abs(G);             %Take the Absolute Value to have positive numbers
    subplot(2,2,4)
    polar(theta,G,'m-')       %Plot
    title(['Frequency = ',num2str(f),'Hz'])
    pause(.6)
end

%Particle Displacement:

f=3000;
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),', Top View using surf1'])
figure(4);
surfl(X,Y,p); 
shading interp;
colormap(pink);
title(['Particle Displacement at Frequency = ',num2str(f),', Side View using surf1'])