%ArrayPatGain - This program is an updated and combined version of the programs "ANTGAIN
% and "ARRAYPATGAIN" in the 4th edition of Kraus's Electromagnetics.  This program
% computes and plots the field pattern of a uniform linear array of sources.  
clear;
timesrun=0;
while timesrun<1000,
   if timesrun==0,
      SP=0.5;
      PH=0;
      N=5;
      MF=1;
   else,
      SP=input('Enter element spacing in wavelengths: ');
      PH=input('Enter phase difference between elements in degrees: ');
      N=input('Enter number of elements: ');
      MF=input('Enter pattern multiplication factor: ')';
   end;
   
% This version modified to include phase plots
% which can be either end referenced or center referenced.
% F. Merat  3/29/03

   
                                      % Compute fields for plotting
A=0.01:0.01:6.27;                     % Angle over 2pi radians

PH=-360*SP-180/N;                    %modified phase for maximum directivity

U=(2*pi*SP*cos(A)+(pi*PH/180))/2;     % compute argument
                                      % Note that /2 is included here rather
									  %   than in sine(NU)/sine(U) expression


FP=(sin(N*U)./sin(U));                % compute isotropic electric field
R=MF.*abs(FP);                        % multiply as appropriate

XI=(180/pi).*((N-1)*U-angle(FP));     % compute phase [end referenced]
                                      % the angle(FP) accounts for the sign
									  % of Sine(Nu)/Sine(u) in the phase

%XI=(180/pi).*(-angle(FP));           % compute phase [center referenced]
                                      % the angle(FP) accounts for the sign
									  % of Sine(Nu)/Sine(u) in the phase



pause(0.9);

                                      % compute the beam area in theta
B=0.01:0.01:3.14;                     % B is theta -- the angle from 0 degrees
                                      %    0 <theta (B) < 180 degrees
   W=(2*pi*SP*cos(B)+(pi*PH/180))/2;  % compute psi -- phase shift
   PP=(sin(N*W)./sin(W)).^2;          % compute unnormalized power
   Z=0.01*sin(B).*PP;                 % differential power is PP*sin(theta)*d(theta)
   SUM=sum(Z);                        % integrate the elements over 180 degrees

DR=(2*(N^2))/SUM;                     % numerical directivity is Pmax/Pavg.
                                      % N^2 is Pmax 
                                      % Not sure of 2; because I only did pi? 
                                      
DBI=10*log10(DR);                     % express in dBi


theta=(0:2*pi/626:2*pi);              % plot antenna pattern in polar coordinates


                                      % plot isotropic pattern

polar(theta,(round(R*100))/100,'r'); 

title(['Field pattern of ',num2str(N),' sources spaced by ',num2str(SP),' lambda,  ','phase delta= ',num2str(PH),' deg']);
text(max(R),max(R),['Dtheta = ',num2str((round(DBI*100))/100),' dBi']);

sep=max(R)/20;
xcent=sep*((-(N-1)/2)-1);
for elcount=1:N
   xcent=xcent+sep;
   ycent=0;
   radius=sep/5;
   xp=xcent+[-radius:radius/10:radius];
   yp=ycent+real(sqrt((radius^2)-(xp-xcent).^2));
   patch(xp,yp,'k');
   patch(xp,-yp,'k');
end;

rp=input('Enter 1 for rectangular plot: ');
if rp==1
   figure(2);
   plot(theta*(180/pi),R);
   grid on;
   zoom on;
   title(['Field Pattern of ',num2str(N),' sources spaced by ',num2str(SP),' lambda,  ','phase delta= ',num2str(PH),' deg']);
   xlabel('Angle (deg)');
   ylabel('Amplitude (linear)');
end;

logp=input('Enter 1 for dB plot: ');
if logp==1
   figure(3);
   plot(theta*(180/pi),10*log10(abs(R)));
   grid on;
   title(['Field Pattern of ',num2str(N),' sources spaced by ',num2str(SP),' lambda,  ','phase delta= ',num2str(PH),' deg']);
   xlabel('Angle (deg)');
   ylabel('Amplitude (dB)');
end;

powerp=input('Enter 1 for power pattern: ');
if powerp==1
   figure(4);
   plot(theta*(180/pi),20*log10(abs(R)));
   grid on;
   title(['Power Pattern of ',num2str(N),' sources spaced by ',num2str(SP),' lambda,  ','phase delta= ',num2str(PH),' deg']);
   xlabel('Angle (deg)');
   ylabel('Amplitude (dB)');
end;


phasep=input('Enter 1 for phase pattern: ');
if phasep==1
   figure(5);
   plot(theta*(180/pi),XI);
   grid on;
   title(['Phase Pattern of ',num2str(N),' sources spaced by ',num2str(SP),' lambda,  ','phase delta= ',num2str(PH),' deg']);
   xlabel('Angle (deg)');
   ylabel('Phase (deg)');
end;

timesrun=timesrun+1;
done=input('Enter 1 to modify parameters: ');
if done~=1,
   timesrun=1000;
end;
end;