clc;clear all;close all;
e=8.854e-12;%介電常數
u=pi*4e-7;%磁導率
c=3e+8;%
N=45;%段數
lambda=1;%波長
L=lambda/2;%天線長度
f=c/lambda;%頻率
w=2*pi*f;%
k=2*pi/lambda;%波數
eta=120*pi;%波阻抗
a=5*10^(-3);%天線半徑   
delta=L/(N+1);
for mm=1:N
    Zb(mm)=((N+1)/2-mm)*delta;
end
psi=zeros(N,N);
for mm=1:N
    for nn=2:N-1
        if mm==nn
            psi(mm,nn)=log(delta/a)/(2*pi)-(j*k*delta)/(4*pi);
        else
            R1=sqrt((Zb(mm)-Zb(nn))^2+a^2);
            psi(mm,nn)=exp(-j*k*R1)*delta/(4*pi*R1);
        end;
    end;
end;
for mm=1:N
    for nn=1:N-2
        PSI(mm,nn)=psi(mm,nn+1);
    end;
end;
for mm=1:N
    PSI(mm,N-1)=cos(k*Zb(mm));
end;
for mm=1:N
    PSI(mm,N)=sin(k*Zb(mm));
end;
for mm=1:N
    temp(mm,1)=sin(k*abs(Zb(mm)))/(j*2*eta);
end;
I=inv(PSI)*temp;
current=zeros(N,1);
for mm=1:(N-2)
    current(mm+1)=I(mm);
end;
mm=-(N-1)/2:1:(N-1)/2;
plot(mm,abs(current),'b'),xlabel('z/delta（z表示天線在ｚ軸上的坐標，delta為天線分割長度）'),ylabel('電流強度'),title('電流分布圖');
Z=1/current((N+1)/2,1)
for mm=1:N
    L_n(mm,1)=((N+1)/2-mm)*delta;
end;
syms thet phi;
for nn=1:N;
    Enn(1,nn)=j*w*u*sin(thet)*exp(j*k*L_n(nn,1)*cos(thet))*delta/(4*pi);
end;
E_n=subs(Enn,{thet},pi/2);
E_max=E_n*current;
G=4*pi*(abs(E_max))^2/(eta*real(Z)*(abs(current((N+1)/2,1)))^2)
E=Enn*current;
F=abs(E)/abs(E_max);
figure;
ezpolar(F);
title('E面方向圖');
figure;
x=F*sin(thet)*cos(phi);
y=F*sin(thet)*sin(phi);
z=F*cos(thet);
ezsurf(x,y,z,120);
title('三維方向圖');