bamu.m

计算八木天线的电流分布

clear
lambda=0.6263;
k=2*pi/lambda;
u=4*pi*10^(-7);
e=8.854*10^(-12);
a=0.0026*lambda;
LR=0.5*lambda;
L=0.47*lambda;
LD=0.43*lambda;
SR=0.25*lambda;
SD=0.3*lambda;
w=k/sqrt(u*e);
y=120*pi;
n=6;
N=5;
dlr=LR/(N+1);
dl=L/(N+1);
dld=LD/(N+1);
point=zeros(n*(2*N+1),4);
mid=zeros(n*N,3);
for ii=1:2*N+1
    point(ii,1:3)=[-SR LR/2-ii*LR/(2*(N+1)) dlr];
    if rem(ii+point(ii,4),2)==0
        mid((ii+point(ii,4))/2,:)=point(ii,1:3);
    end
end
for ii=2*N+1+1:2*(2*N+1)
    point(ii,2:4)=[L/2-(ii-(2*N+1))*L/(2*(N+1)) dl 1];
    if rem(ii+point(ii,4),2)==0
        mid((ii-point(ii,4))/2,2:3)=point(ii,2:3);
    end
end
for ii=2*(2*N+1)+1:3*(2*N+1)
    point(ii,:)=[SD LD/2-(ii-2*(2*N+1))*LD/(2*(N+1)) dld 2];
    if rem(ii+point(ii,4),2)==0
        mid((ii-point(ii,4))/2,:)=point(ii,1:3);
    end
end
for ii=3*(2*N+1)+1:4*(2*N+1)
    point(ii,:)=[2*SD LD/2-(ii-3*(2*N+1))*LD/(2*(N+1)) dld 3];
    if rem(ii+point(ii,4),2)==0
        mid((ii-point(ii,4))/2,:)=point(ii,1:3);
    end
end
for ii=4*(2*N+1)+1:5*(2*N+1)
    point(ii,:)=[3*SD LD/2-(ii-4*(2*N+1))*LD/(2*(N+1)) dld 4];
    if rem(ii+point(ii,4),2)==0
        mid((ii-point(ii,4))/2,:)=point(ii,1:3);
    end
end
for ii=5*(2*N+1)+1:6*(2*N+1)
    point(ii,:)=[4*SD LD/2-(ii-5*(2*N+1))*LD/(2*(N+1)) dld 5];
    if rem(ii+point(ii,4),2)==0
        mid((ii-point(ii,4))/2,:)=point(ii,1:3);
    end
end
V=zeros(n*N,1);
V(N+(N+1)/2)=1;
U=ones(n*N,1);
psi=zeros(n*(2*N+1));
for jj=1:n*(2*N+1)
    for kk=1:n*(2*N+1)
        if jj==kk
            psi(jj,kk)=log(point(jj,3)/a)/(2*pi*point(jj,3))-(j*k)/(4*pi);
        else
            psi(jj,kk)=exp(-j*k*sqrt((point(kk,1)-point(jj,1))^2+(point(kk,2)-point(jj,2))^2))/(4*pi*sqrt((point(kk,1)-point(jj,1))^2+(point(kk,2)-point(jj,2))^2));
        end
    end
end
Z=zeros(n*N);
for pp=1:n*N
    for qq=1:n*N
        Z(pp,qq)=j*w*u*point(pp,3)*point(qq,3)*psi(2*pp+point(pp,4),2*qq+point(qq,4))+(psi(2*pp+point(pp,4)+1,2*qq+point(qq,4)+1)-psi(2*pp+point(pp,4)+1,2*qq+point(qq,4)-1)-psi(2*pp+point(pp,4)-1,2*qq+point(qq,4)+1)+psi(2*pp+point(pp,4)-1,2*qq+point(qq,4)-1))/(j*w*e);
    end
end
si=Z\V;    %In
t=1:n*N;
figure(1);
plot(t,abs(si)),ylabel('I'),title('电流分布')
in=U'*(Z\V);
i=V'*si;
Zin=1/i
theta=(-pi:pi/100:pi)+eps;
for m=1:length(theta)
    E1=-j*w*u*exp(-j*k).*exp(j*k.*sqrt(mid(:,1).^2+mid(:,2).^2).*cos(abs(atan(mid(:,1)./(mid(:,2)+eps))-theta(m)))).*mid(:,3).*sin(theta(m))/(4*pi);
    Etheta(m)=E1'*si;
end
Etheta=Etheta./max(Etheta);
figure(2);
polar(theta,abs(Etheta)/max(abs(Etheta))),title('E平面方向图 (\Phi = 0)');
Lo=find((abs(Etheta-1/sqrt(2))<0.05)==1);
G=abs(4*pi.*Etheta.*conj(Etheta)/(y*real(Zin).*si((N+1)/2).*conj(si((N+1)/2))));
Gmax=max(G)
phi=(0:pi/100:2*pi)+eps;
for m=1:length(theta)
    E2=-j*w*u*exp(-j*k).*exp(j*k.*sqrt(mid(:,1).^2+mid(:,2).^2).*cos(abs(atan(mid(:,1)./(mid(:,2)+eps))-theta(m)))).*mid(:,3).*sin(theta(m))/(4*pi);
    Ephi(m)=E2'*abs(si);
end
figure(3);
polar(theta-pi/2,(abs(Etheta)/120/pi)/max(abs(Etheta)/120/pi)),title('H平面方向图 (\theta = \pi / 2)');