%矩量法：由Hallen方程求天線表面不同饋點電流分布；采用點匹配法，基函數為：脈沖函數
clear all
clc
Lxx=1;
ft=3e8;%頻率
c=3e8;%光速
lamda=c/ft;%波長
L=Lxx*lamda;%天線長度
a0=0.0025*lamda;%天線半徑
k0=2*pi/lamda;%波數
z0=120*pi;%波阻抗
E0=1;%激勵
N=41;
xs=0.5;%饋電點0.25處
for m=1:N
    zz1=-xs*L+(m-1)*L/N;
    zz2=-xs*L+m*L/N;
    zz=(zz1+zz2)/2;
    for n=2:(N-1)
           aa=-xs*L+(n-1)*L/N;%積分下限
           bb=-xs*L+n*L/N;%積分上限
           x=(aa+bb)/2;
           if(n==m)
               s1(m,n-1)=(log(L/(N*a0)))/(2*pi)-j*k0*(L/N)/(4*pi);
           else
               s1(m,n-1)=(L/N)*(exp(-j*k0*sqrt((zz-x)^2+a0^2)))/(4*pi*sqrt((zz-x)^2+a0^2));
           end
    end
       s1(m,N-1)=cos(k0*zz);
       s1(m,N)=sin(k0*zz);    
       s2(m,1)=-j*(E0*L/N)*sin(k0*abs(zz))/(2*z0);
end
%det_s1=det(s1)
s3=inv(s1)*s2;%系數Cn
%c1=s3(N-1,1)%系數c1
%c2=s3(N,1)%系數c2
xx=L/(2*N*lamda):L/(N*lamda):(2*N-1)*L/(2*N*lamda);
s(1,1:N)=0;
for n=1:(N-2)
    if(abs(s3(n,1))<1e-10)
        s(1,n+1)=0;
    else
        s(1,n+1)=s3(n,1);
    end
end  
figure(1)
subplot(2,2,1);
plot(xx,abs(s),'k');
legend('幅度');
xlabel('L(lamda)');
ylabel('Magnitude(mA)');
title('電流幅度與直線天線長度關系曲線');
axis tight

subplot(2,2,2);
plot(xx,angle(s)*180/pi,'k');
legend('相位');
xlabel('L(lamda)');
ylabel('Phase(degree)');
title('電流相位與直線天線長度關系曲線');
axis tight
%%%%%%%%%%%%%%%%%%%%%%畫E面方向圖%%%%%%%%%%%%%%%%%%%%%%%%
z=(-xs*L+L/(2*N)):L/N:(-xs*L+(2*N-1)*L/(2*N));
thita=0:0.01:2*pi;
for n=1:length(thita)
    u=0;
    for m=1:N
          u=u+s(1,m).*sin(thita(1,n)).*exp(j*k0*z(1,m).*cos(thita(1,n))).*L/N;       
    end
    f(1,n)=u;
end
subplot(2,2,3);
polar(thita,abs(f)./max(abs(f)),'k');
title('E面方向圖')

Zr=(E0*L/N)/s(1,round(N*xs))