pocklington.m

Pocklington方法的ＭＡＴＬＡＢ环境下半波振子仿真

clear all;close all;
e=8.854e-12;%介电常数
u=pi*4e-7;%磁导率
c=3e+8;
N=41;
a=0.0000001;%天线半径
lambda=1;
L=lambda/2;
f=c/lambda;
k=2*pi/lambda;%波数
z=zeros(1,N);
delta=L/(N+1);
z(1)=-L/2+delta;
for n=2:N
    z(n)=z(n-1)+delta;
end;
for mm=1:N+2
    for nn=1:N+2
        if mm==nn
            psi(mm,nn)=log(delta/a)/(2*pi*delta)-j*k/4/pi;
        else
            R(mm,nn)=abs(mm-nn)*delta;
            psi(mm,nn)=exp(-j*k*R(mm,nn))/(4*pi*R(mm,nn));
        end;
    end;
end;
for mm=1:N
    for nn=1:N
        Z(mm,nn)=((2-k*k*delta*delta)*psi(mm+1,nn+1)-psi(mm+2,nn+1)-psi(mm,nn+1))/(j*2*pi*f*e);
    end;
end;
U=zeros(N,1);
U((N+1)/2,1)=1;
I=eye(N,N)*inv(Z)*U;
Zin=1/I((N+1)/2);%输入阻抗
current=zeros(N+2,1);%加入两端零电流
for nn=2:N+1
    current(nn)=I(nn-1);
end;
figure(1);%电流分布图
nn=-(N+1)/2:1:(N+1)/2
plot(nn,abs(current),'r'),xlabel('z/delta（z表示天线在ｚ轴上的坐标，delta为天线分割长度）'),ylabel('电流强度'),title('电流分布图');
mm=1,r=1;
for thet=0:0.01*pi:pi
    for nn=1:N
        temp=2*pi*f*u*delta/(4*pi*r);
        En(1,nn)=temp*(exp(-j*k*r)*exp(j*k*delta*nn*cos(thet))*sin(thet));
    end;
    E=En*inv(Z)*U;
    G(mm)=(4*pi*r^2)/sqrt(u/e)/abs(I((N+1)/2))^2/real(Zin)*abs(E)^2;
    mm=mm+1;
end;
Gmax=max(G)
figure(2);%E方向性图
mm=1;
for thet=0:pi/100:2*pi;
    for nn=1:nn
        F_E(nn)=current(nn+1)*delta*exp(j*k*delta*nn*cos(thet))*sin(thet);
    end;
    FF_E(mm)=sum(F_E);
    mm=mm+1;
end;
FF_E_max=max(abs(FF_E));
for nn=1:mm-1
    FF_E(nn)=FF_E(nn)/FF_E_max;
end;
thet=0:pi/100:2*pi;
polar(thet,abs(FF_E),'r');title('E面方向图');
figure(3);%H方向性图
mm=1,thet=pi/2;
for psi=0:pi/100:2*pi
    for nn=1:N
        F_H(nn)=current(nn+1)*delta*exp(j*k*delta*nn*cos(thet))*sin(thet);
    end;
    FF_H(mm)=sum(F_H);
    mm=mm+1;
end;
FF_H_max=max(abs(FF_H));
for nn=1:mm-1
    FF_H(nn)=FF_H(nn)/FF_H_max;
end;
psi=0:pi/100:2*pi;
polar(psi,abs(FF_H),'r');title('H面方向图');
f0=c/lambda;%对于既定的天线长度,其对应的半波频率.
counter=1;
for alpha=0.5:0.01:1.5
    f_temp=alpha*f0;
    lambda_temp=c/f_temp;
    k_temp=2*pi/lambda_temp;
    for mm=1:N+2
        for nn=1:N+2
            if mm==nn
                psi(mm,nn)=log(delta/a)/(2*pi*delta)-j*k_temp/4/pi;
            else R(mm,nn)=abs(mm-nn)*delta;
                psi(mm,nn)=exp(-j*k_temp*R(mm,nn))/(4*pi*R(mm,nn));
            end;
        end;
    end;
    for mm=1:N
        for nn=1:N
            Z(mm,nn)=((2-k_temp*k_temp*delta*delta)*psi(mm+1,nn+1)-psi(mm+2,nn+1)-psi(mm,nn+1))/(j*2*pi*f_temp*e);
        end;
    end;
    U((N+1)/2,1)=1;
    current=eye(N,N)*inv(Z)*U;
    Zin(counter)=1/current((N+1)/2);%输入阻抗
    Real_In(counter)=real(Zin(counter));
    Imaginary_In(counter)=imag(Zin(counter));
    counter=counter+1;
end;
figure(4);%输入阻抗频率响应图
mm=0.5:0.01:1.5;
subplot(1,2,1),plot(mm,Real_In,'r'),xlabel('f/f0');ylabel('电阻');title('天线长度固定,');
subplot(1,2,2),plot(mm,Imaginary_In,'r'),xlabel('f/f0');ylabel('电抗');title('f0等于当其为半波振子时对应的频率');