function I = PhasedArrayPlotMultiFrequency(N,d,a,lambda,td)
%N is the number of elements
%d is the distance between elements (center to center)
%a is the width of one element
%lambda is the wavelength
%td is the group time delay


%Initializations
theta = linspace(-90,90,10000);          %Create the angle vector to plot against
lambda_vect = [lambda;lambda/2;lambda/4;lambda/8;lambda/16];
phi = [td/lambda_vect(1)*343.3;td/lambda_vect(2)*343.3;td/lambda_vect(3)*343.3;td/lambda_vect(4)*343.3;td/lambda_vect(5)*343.3];     %Create the phi vector
I = zeros(5,length(theta));

%Start doing the calculations to measure the intensity of the wave relative
%to the angle
for k = 1:5
    A = sin(pi*a.*sind(theta)/lambda_vect(k));
    B = pi*a.*sind(theta)/lambda_vect(k);
    C = ((2*pi*d.*sind(theta)/lambda_vect(k))+phi(k));
    D = sin((N/2).*C);
    E = sin(C);
    I(k,:) = ((A./B).^2).*((D./E).^2);
end

%Plot I vs theta
plot(theta,I(1,:),'b',theta,I(2,:),'r',theta,I(3,:),'g',theta,I(4,:),'c',theta,I(5,:),'k');
xlabel('Viewing Angle Theta');
ylabel('Intensity of the wave');
title('Intensity vs Theta');
legend('200 Hz','400 Hz','800 Hz','1600 Hz','3200 Hz');

