% linear array of dipoles with a ground plane
% z is array axis; ground plane is the xz plane
% y is normal to ground plane
% option to select dipoles along x or z
% fixed to include Etheta and Ephi
% two dimensional array SCANS IN BOTH X AND Z DIMENSION
% modified to do half wave dipoles
clear
warning off
rad=pi/180;
f=2.4e9;
wave=3e8/f;
beta=2*pi/wave;   
dz=0.065; 
Nz=8;
dx=0.1;
Nx=2;
h=0.25*wave;
% select the dipole direction: 'x' or 'z'
dipdir='x';  % dipoles x directed (parallel case)
%dipdir='z';  % dipoles z directed (collinear case)
% select the dipole type: 'h'=halfwave, 's'=short (constant current) 
diptype='s';
%diptype='h';

%%%%%%%%%%%%%%%%%%%%%%% phi=0 cut %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% phi=0 vary theta (x-z plane)
phi=90*rad;
theta=[0:1:180]*rad; w=cos(theta);
thetas=90; phis=90;
cts=cos(thetas*rad); u=sin(theta)*cos(phi); v=sin(theta)*sin(phi);
cps=cos(phis*rad);
alphaz=-beta*dz*cts;
alphax=-beta*dx*cps;
psiz=beta*dz*w+alphaz;
psix=beta*dx*u+alphax;
if dipdir=='z'   % colinear dipoles (z directed)
    if diptype=='s'
    EFtheta=2*sqrt(1-w.^2).*sin(beta*h*v);
    EFphi=0;
    end
    if diftype=='h'
    EFtheta=2*sqrt(1-w.^2).*sin(beta*h*v).*cos(pi*w/2)./sqrt(1-w.^2+1e-5);
    EFphi=0;
    end
end
if dipdir=='x' % parallel case (x directed)
    if diptype=='s'
    EFtheta=2*cos(phi).*w.*sin(beta*h*v); 
    EFphi=-2*sin(phi).*sin(beta*h*v);
    end
    if diptype=='h'
    EFtheta=2*cos(phi).*w.*sin(beta*h*v).*cos(pi*u/2)./(1-u.^2+1e-5); 
    EFphi=-2*sin(phi).*sin(beta*h*v).*cos(pi*u/2);
    end
end
AF=sin(Nx*psix/2)./sin(psix/2).*sin(Nz*psiz/2)./sin(psiz/2);
caftheta=abs(AF.*EFtheta)/2/Nx/Nz;
cafphi=abs(AF.*EFphi)/2/Nx/Nz;
cafthdb=20*log10(caftheta+1e-5);
cafphdb=20*log10(cafphi+1e-5);
caf=sqrt(cafphi.^2+caftheta.^2);
cafdb=10*log10(caf.^2);
figure(1),clf
subplot(121)
polardb(theta/rad,cafthdb,'-')
title(['Polar dB plot E_\theta, \phi=',num2str(phi/rad),'^o'])
subplot(122)
polardb(theta/rad,cafphdb,'-')
title(['Polar dB plot E_\phi, \phi=',num2str(phi/rad),'^o'])
hold off
%%%%%%%%%%%%%%%%%%%%%%%%%% theta=90 cut %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear theta phi AF caf cafthdb cafphdb EFth EFph psi u v w cts ct
% theta=90 vary theta (x-y plane)
phi=[0:1:180]*rad;
theta=89.99*rad; ct=cos(theta);
u=sin(theta)*cos(phi);
cts=cos(thetas*rad); v=sin(theta).*sin(phi); w=ones(1,length(v))*ct;
alphaz=-beta*dz*cts;
psiz=(beta*dz*cos(theta)+alphaz)*ones(length(phi),1)';
psix=(beta*dx*u+alphax);
if dipdir=='z'  % colinear dipoles (z directed)
    if diptype=='s'
    EFth=2*sqrt(1-w.^2).*sin(beta*h*v); 
    EFph=0;
    end
    if diftype=='h'
    EFth=2*sqrt(1-w.^2).*sin(beta*h*v).*cos(pi*w/2)./sqrt(1-w.^2+1e-5);
    EFph=0;
    end
end
if dipdir=='x'  % parallel case (x directed)
    if diptype=='s'
    EFth=2*cos(phi).*w.*sin(beta*h*v);
    EFph=-2*sin(phi).*sin(beta*h*v);
    end
    if diptype=='h'
    EFth=2*cos(phi).*w.*sin(beta*h*v).*cos(pi*u/2)./(1-u.^2+1e-5); 
    EFph=-2*sin(phi).*sin(beta*h*v).*cos(pi*u/2);
    end 
end
AF=sin(Nx*psix/2)./sin(psix/2).*sin(Nz*psiz/2)./sin(psiz/2);
caftheta=abs(AF.*EFth)/2/Nx/Nz;
cafphi=abs(AF.*EFph)/2/Nx/Nz;
cafthdb=20*log10(caftheta+1e-5);
cafphdb=20*log10(cafphi+1e-5);
caf=sqrt(cafphi.^2+caftheta.^2);
cafdb=10*log10(caf.^2);
figure(2),clf
subplot(121)
polardb(phi/rad,cafthdb,'-')
title(['Polar dB plot E_\theta, \theta=',num2str(theta/rad),'^o'])
subplot(122)
polardb(phi/rad,cafphdb,'-')
title(['Polar dB plot E_\phi, \theta=',num2str(theta/rad),'^o'])
hold off
% clear theta phi AF caf cafdb EFth EFph psi
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3-D plot %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
iplt3d=0; % skip 3d plots if =1
if iplt3d==0
 it=0; 
 for th=0:3:180
	 jt=0; it=it+1;
   for ph=0:3:360
       jt=jt+1;
       thr=th*rad;
	   phr=ph*rad;
	   phd(it,jt)=ph;
	   thd(it,jt)=th; 
       ct=cos(thr);
       st=sin(thr);
       u=sin(thr)*cos(phr);
	   v=sin(thr)*sin(phr);
	   w=cos(thr);
       xsiz=beta*dz*cos(thr)+alphaz;
       xsix=beta*dx*u+alphax;
       E=0;
       if v>=0;  % upper half space (lower half space zero)
            Ex=Nx;
            Ez=Nz;
            if abs(xsiz)>1e-6
                Ez=sin(Nz*xsiz/2)/sin(xsiz/2);
            end
            if abs(xsix)>1e-6
                Ex=sin(Nx*xsix/2)/sin(xsix/2);
            end
            E=Ex*Ez;
       end
       if dipdir=='z'
           EFth=2*sqrt(1-w^2)*sin(beta*h*v); % colinear dipoles (z directed)
           EFph=0;
       end
       if dipdir=='x'
           EFth=2*cos(phr).*w*sin(beta*h*v);  % parallel case (x directed)
           EFph=-2*sin(phr)*sin(beta*h*v);
       end
	   ffth=abs(EFth*E);
       ffph=abs(EFph*E);
	   fth(it,jt)=ffth;
       fph(it,jt)=ffph;
	   Xth(it,jt)=fth(it,jt)*u;
	   Yth(it,jt)=fth(it,jt)*v;
	   Zth(it,jt)=fth(it,jt)*w;
       Xph(it,jt)=fph(it,jt)*u;
	   Yph(it,jt)=fph(it,jt)*v;
	   Zph(it,jt)=fph(it,jt)*w;
   end
 end
%plot theta component
MAX=max([max(max(Xth)),max(max(Yth)),max(max(Zth)),max(max(Xph)),max(max(Yph)),max(max(Zph))]);
figure(3),clf
mesh(Xth/MAX,Yth/MAX,Zth/MAX)
axis([-1,1,-1,1,-1,1])
axis square
hold on
% spread N elements along the z axis between [-1,1]
z0=[-(Nz-1)/2:(Nz-1)/2]*.75/Nz;
x0=[-(Nx-1)/2:(Nx-1)/2]*.25/Nx;
[X,Z]=meshgrid(x0,z0); Y=ones(size(X))*.05;
plot3(X,Y,Z,'ko')
xlabel('X=|E|sin(t)cos(p)')
ylabel('Y=|E|sin(t)sin(p)')
zlabel('Z=|E|cos(t)');
if dipdir=='x', title('E_\theta for x directed dipoles (parallel case)'), end
if dipdir=='z', title('E_\theta for z directed dipoles (collinear case)'), end

%plot phi component
figure(4),clf
mesh(Xph/MAX,Yph/MAX,Zph/MAX)
axis([-1,1,-1,1,-1,1])
axis square
hold on
% spread N elements along the z axis between [-1,1]
z0=[-(Nz-1)/2:(Nz-1)/2]*.75/Nz;
x0=[-(Nx-1)/2:(Nx-1)/2]*.25/Nx;
[X,Z]=meshgrid(x0,z0); Y=ones(size(X))*.05;
plot3(X,Y,Z,'ko')
xlabel('X=|E|sin(t)cos(p)')
ylabel('Y=|E|sin(t)sin(p)')
zlabel('Z=|E|cos(t)');
if dipdir=='x', title('E_\phi for x directed dipoles (parallel case)'), end
if dipdir=='z', title('E_\phi for z directed dipoles (collinear case)'), end
end


