%Michael Battaglia
%Acoustics Project 2 Part 1b
%Modified from "tappered_onedimensionalpiston.m"

% You are required to set up a one-dimensional piston and plot its 
% beam pattern and particle displacement in the far field. 
% Select the piston parameters so that for frequencies between 
% 1 kHz and 5 kHz it produces a main lobe with maximum width equal 
% to 20o. Clearly state the selected piston parameters. 
% Provide beam plots for 1, 2, 3, 4 and 5 kHz. 
clear all, close all
b=-1/2:(1/200):1/2;
bpi=b*pi;
l=1;
c=354;
f=1500;
lambda=c./f;
i=1;
    for j=1:length(bpi)
    p0(i,j)=sin(pi*(l/lambda(i))*sin(bpi(j)))/(pi*(l/lambda(i))*sin(bpi(j)));
    p1(i,j)=0.5*sin(pi*((l/lambda(i))*sin(bpi(j))+1))/(pi*((l/lambda(i))*sin(bpi(j))+1));
    p2(i,j)=0.5*sin(pi*((l/lambda(i))*sin(bpi(j))-1))/(pi*((l/lambda(i))*sin(bpi(j))-1));
    p=p0;
    q(i,j)=p0(i,j)+p1(i,j)+p2(i,j);
    end
    F=c./lambda(i);
    %
    [X,Y] = meshgrid(0:0.01:1,-1.001:0.01:1); 
    k=2*pi*F/c;
    sin_t=abs(Y)./sqrt(X.^2+Y.^2);  %the sinusoid
    u=l*sin_t/lambda(i);
    %m0=pi*u;    m1=0.5*pi*(u+1));    m1=0.5*pi*((l*sin_t-1)/lambda(i));
    R=sqrt(X.^2+Y.^2);
    q=((sin(pi*u)./(pi*u))+0.5*(sin(pi*(u+1))./(pi*(u+1)))+0.5*(sin(pi*(u-1))./(pi*(u-1)))).*cos(k*R)./R;
    surf(X,Y,q+1), axis ([0 1 -1 1 -10 20])%q+1 is just a level adjustment for better color
    xlabel('X'), ylabel('Y'), zlabel('Z');
    caxis([-5 5])
    title(['Displacement pattern for piston size of size length/lambda= ', num2str(l/lambda(i))])
    view (30, 50)
    figure(gcf)