momantena.m

matlab求解半波天线,矩量法解双位方程

%%%%矩量法解双位方程，求直线天线表面电流分布 
clear all 
clc 
ft=3e8;%频率 
c=3e8;%光速 
ep=(10^(-9))/(36*pi); 
mu=4*pi*(10^(-7)); 
w=2*pi*ft; 
lamda=c/ft;%波长 
a0=0.0025*lamda;%天线半径 
L=1.77*lamda;%天线长度 
k0=2*pi/lamda;%波数 
z0=120*pi;%波阻抗 
N=51; 
xs0=1/2; 
v0=L/N; 
for m=1:(N-2) 
    zs1=-xs0*L+m*L/N;%m- 
    zs2=-xs0*L+(m+1)*L/N;%m+ 
    zs12=(zs1+zs2)/2;%m 
    for n=1:(N-2) 
        zs3=-xs0*L+n*L/N;%n- 
        zs4=-xs0*L+(n+1)*L/N;%n+ 
        zs34=(zs3+zs4)/2;%n 
        rn=sqrt((zs34-zs12)^2+a0^2);%m,n 
        rn11=sqrt((zs2-zs4)^2+a0^2);%m+,n+ 
        rn12=sqrt((zs1-zs4)^2+a0^2);%m-,n+ 
        rn21=sqrt((zs2-zs3)^2+a0^2);%m+,n- 
        rn22=sqrt((zs1-zs3)^2+a0^2);%m-,n- 
        %aa=m-n; 
        switch m-n 
            case 0  
            phi=(1/(2*pi*L/N))*log(L/(a0*N))-j*k0/(4*pi); 
            phi11=(1/(2*pi*L/N))*log(L/(a0*N))-j*k0/(4*pi); 
            phi12=exp(-j*k0*rn12)/(4*pi*rn12); 
            phi21=exp(-j*k0*rn21)/(4*pi*rn21); 
            phi22=phi11; 
            case 1 
            phi=exp(-j*k0*rn)/(4*pi*rn); 
            phi11=exp(-j*k0*rn11)/(4*pi*rn11); 
            phi12=(1/(2*pi*L/N))*log(L/(a0*N))-j*k0/(4*pi); 
            phi21=exp(-j*k0*rn21)/(4*pi*rn21); 
            phi22=exp(-j*k0*rn22)/(4*pi*rn22); 
            case -1 
            phi=exp(-j*k0*rn)/(4*pi*rn); 
            phi11=exp(-j*k0*rn11)/(4*pi*rn11); 
            phi12=exp(-j*k0*rn12)/(4*pi*rn12); 
            phi21=(1/(2*pi*L/N))*log(L/(a0*N))-j*k0/(4*pi); 
            phi22=exp(-j*k0*rn22)/(4*pi*rn22); 
            otherwise 
            phi=exp(-j*k0*rn)/(4*pi*rn); 
            phi11=exp(-j*k0*rn11)/(4*pi*rn11); 
            phi12=exp(-j*k0*rn12)/(4*pi*rn12); 
            phi21=exp(-j*k0*rn21)/(4*pi*rn21); 
            phi22=exp(-j*k0*rn22)/(4*pi*rn22); 
         end 
         z(m,n)=j*w*mu*phi*(L/N)^2+(phi11-phi12-phi21+phi22)/(j*w*ep); 
     end 
end 
%det_z=det(z) 
v=zeros(N-2,1); 
I1=zeros(1,N); 
v(round((N-2)*xs0),1)=v0; 
I2=inv(z)*v; 
for n=1:(N-2) 
    I1(1,n+1)=I2(n,1); 
end 
%I' 
figure(1) 
x=L/(2*N*lamda):L/(N*lamda):(2*N-1)*L/(2*N*lamda); 
subplot(2,2,1); 
plot(x,abs(I1),'k'); 
xlabel('L/lamda'); 
ylabel('Magnitude'); 
title('电流幅度'); 
axis tight 
 
subplot(2,2,2); 
plot(x,180*angle(I1)/pi,'k'); 
xlabel('L/lamda'); 
ylabel('Phase'); 
title('电流相位'); 
axis tight 
%%%%%%%%%%%%%%%%%%%%%%画E面方向图%%%%%%%%%%%%%%%%%%%%%%%% 
z=(1-N)*L/(2*N):L/N:(N-1)*L/(2*N); 
thita=0:0.01:2*pi; 
for n=1:length(thita) 
    u=0; 
    for m=1:N 
        u=u+I1(1,m).*sin(thita(1,n)).*exp(j*k0*z(1,m).*cos(thita(1,n))).*L/(N*lamda); 
    end 
    f(1,n)=u; 
end 
subplot(2,2,3); 
polar(thita,abs(f)./max(abs(f)),'k'); 
title('E面方向图') 
 
%s(1,round(N/2)) 
%max(s) 
Zr=v0./I1(1,round(N*xs0))