% Project 2 # 1: sound source dipole

%Beam Plot:

c=345;
h=.033;                      %Dipole Size
theta=-pi/2:pi/180:3*pi/2;
for f=1000         %Frequency
    k=2*pi*f/c;           %wave number
    p=1-exp(j*k*h*sin(theta)); %Gain as a function of the elevation angle, Q2/Q1=-Q1/Q1=-1, Beam Pattern
    p=abs(p);              %Take the Absolute Value to have positive numbers
    subplot(2,2,1)
    polar(theta,p,'m-')        %Plot
    title(['Beam Pattern for array with frequency = ',num2str(f),'Hz'])
    pause(.6)              %Amount of time before f is increased
end
for f=3000         %Frequency
    k=2*pi*f/c;           %wave number
    p=1-exp(j*k*h*sin(theta)); %Gain as a function of the elevation angle, Q2/Q1=-Q1/Q1=-1, Beam Pattern
    p=abs(p);              %Take the Absolute Value to have positive numbers
    subplot(2,2,2)
    polar(theta,p,'m-')        %Plot
    title(['Beam Pattern for array with frequency = ',num2str(f),'Hz'])
    pause(.6)              %Amount of time before f is increased
end
for f=5500        %Frequency
    k=2*pi*f/c;           %wave number
    p=1-exp(j*k*h*sin(theta)); %Gain as a function of the elevation angle, Q2/Q1=-Q1/Q1=-1, Beam Pattern
    p=abs(p);              %Take the Absolute Value to have positive numbers
    subplot(2,2,3)
    polar(theta,p,'m-')        %Plot
    title(['Beam Pattern for array with frequency = ',num2str(f),'Hz'])
    pause(.6)              %Amount of time before f is increased
end
for f=100000         %Frequency
    k=2*pi*f/c;           %wave number
    p=1-exp(j*k*h*sin(theta)); %Gain as a function of the elevation angle, Q2/Q1=-Q1/Q1=-1, Beam Pattern
    p=abs(p);              %Take the Absolute Value to have positive numbers
    subplot(2,2,4)
    polar(theta,p,'m-')        %Plot
    title(['Beam Pattern for array with frequency = ',num2str(f),'Hz'])
    pause(.6)              %Amount of time before f is increased
end

%Particle Displacement:

f=1000;
c=345;
k=2*pi*f/c;
h=.02;
[X,Y] = meshgrid(-10:.1:10);
r=sqrt(X.^2+(Y+(.5*h)).^2);
R=sqrt(X.^2+(Y-(.5*h)).^2);
p1=cos(k*r)./r;
p2=cos(k*R)./R;
figure(2);
meshc(X,Y,p2-p1);
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,p2-p1)
title(['Particle Displacement at Frequency = ',num2str(f),', Angled View using surf1'])
figure(4);
surfl(X,Y,p2-p1); 
shading interp;
colormap(pink);
title(['Particle Displacement at Frequency = ',num2str(f),', Side View using surf1'])