%%%%矩量法解Hallen方程，求直線天線表面電流分布%%%%%
clear all
clc
ep=8.854e-12;%介電常數
c=3e8;%光速
ft=3e8;%頻率
wt=2*pi*ft;%角頻率
lamda=c/ft;%波長
a=0.0001*lamda;%天線半徑
L=1.7*lamda;%天線長度
k=2*pi/lamda;%波數
V=1;%激勵電壓
no=120*pi;%波阻抗


N=121;%分段數
pmn=zeros(N,N);
%qm=zeros(N,1);
%sm=zeros(N,1);
tm=zeros(N,1);
In=zeros(N,1);%電流
dL=L/N;%小段長

Q=linspace(-L/2,L/2,N+1);
for i=1:N
    W(i)=(Q(i)+Q(i+1))/2;
end

for m=1:N
  for n=1:N-1
     if m==n
        pmn(m,n)=log(dL/a)/(2*pi)-(j*k*dL)/(4*pi);
     else
        Rmn=sqrt((W(m)-W(n))^2+a^2);
        pmn(m,n)=exp(-j*k*Rmn)*dL/(4*pi*Rmn);
     end
  end
end

for i=1:N
pmn(i,N)=cos(k*W(i));%qm
end

for i=1:N
tm(i,1)=V*sin(k*abs(W(i)))/(j*2*no);
end

In=inv(pmn)*tm;
In(N,1)=0;
x=(W(1))/(lamda):L/(N*lamda):(W(N))/(lamda);
subplot(2,1,1);
plot(x,abs(In));
xlabel('L/lamda');
ylabel('電流幅值');
title('電流分布');
axis tight

subplot(2,1,2);
plot(x,180*angle(In)/pi);
xlabel('L/lamda');
ylabel('電流相位');
axis tight