clc;clear all;close all;
w=2*pi*3e8; % 頻率
E=8.854e-12;% 真空介電常數(shù)
u=pi*4e-7;% 磁導(dǎo)率
eta=120*pi;% 自由空間的波阻抗
lm=1;% 波長
k=2*pi/lm;% 波數(shù)
l=0.5*lm;% 天線長度為半波長
a=0.0000005;% 天線半徑
N=31;% 分段數(shù)
div=l/(N+1);% 單位長度
alpha=pi/2;%半張角
for mm=1:((N+1)/2)          
    for nn=1:((N+1)/2)     
        R(mm,nn)=abs(div*(mm-nn));
    end
end
for mm=1:((N+1)/2)
    for nn=((N+3)/2):N 
        R(mm,nn)=div*sqrt((N/2+1/2-mm)^2+(nn-(N+1)/2)^2-2*((N+1)/2-mm)*(nn-(N+1)/2)*cos(2*alpha));
    end
end
for mm=((N+3)/2):N       
    for nn=((N+3)/2):N 
        R(mm,nn)=abs(div*(mm-nn));
    end
end
for mm=((N+3)/2):N 
    for nn=1:((N+1)/2)
        R(mm,nn)=div*sqrt((mm-(N+1)/2)^2+((N+1)/2-nn)^2-2*(mm-(N+1)/2)*((N+1)/2-nn)*cos(2*alpha));
    end
end
for mm=1:((N+1)/2)
    for nn=1:((N+1)/2)
        R1(mm,nn)=abs(div*(mm-nn+1));
    end
end
for mm=1:((N-1)/2)
    for nn=((N+3)/2):N
        R1(mm,nn)=div*sqrt((N/2-mm)^2+(nn-N/2-1)^2-2*(N/2-mm)*(nn-N/2-1)*cos(2*alpha));
    end
end
for mm=((N+1)/2):N
    for nn=((N+3)/2):N
        R1(mm,nn)=abs(div*(mm-nn+1));
    end
end
for mm=((N+3)/2):N
    for nn=1:((N+1)/2)
        R1(mm,nn)=div*sqrt((mm-N/2)^2+(N/2+1-nn)^2-2*(mm-N/2)*(N/2+1-nn)*cos(2*alpha));
    end
end
for mm=1:((N+1)/2)
    for nn=1:((N-1)/2)
        R2(mm,nn)=abs(div*(mm-nn-1));
    end
end
for mm=1:((N+1)/2)
    for nn=((N+1)/2):N
        R2(mm,nn)=div*sqrt((N/2+1-mm)^2+(nn-N/2)^2-2*(N/2+1-mm)*(nn-N/2)*cos(2*alpha));
    end
end
for mm=((N+3)/2):N
    for nn=1:((N-1)/2)
        R2(mm,nn)=div*sqrt((mm-N/2-1)^2+(N/2-nn)^2-2*(mm-N/2-1)*(N/2-nn)*cos(2*alpha));
    end
end
for mm=((N+3)/2):N
    for nn=((N+1)/2):N
         R2(mm,nn)=abs(div*(mm-nn-1));
     end
end
for mm=1:N
    for nn=1:N
        if mm==nn
            Fi(mm,nn)=log(div/a)/(2*pi*div)-j*k/(4*pi);
        else
            Fi(mm,nn)=exp(-j*k*R(mm,nn))/(4*pi*R(mm,nn));
        end
    end
end
for mm=1:N
    for nn=1:N
        if mm==nn-1
            Fi1(mm,nn)=log(div/a)/(2*pi*div)-j*k/(4*pi);
        else
            Fi1(mm,nn)=exp(-j*k*R1(mm,nn))/(4*pi*R1(mm,nn));
        end
    end
end
for mm=1:N
    for nn=1:N
        if mm==nn+1
            Fi2(mm,nn)=log(div/a)/(2*pi*div)-j*k/(4*pi);
        else
            Fi2(mm,nn)=exp(-j*k*R2(mm,nn))/(4*pi*R2(mm,nn));
        end
    end
end
for mm=1:N
    for nn=1:N
        Z(mm,nn)=j*w*u*div*div*Fi(mm,nn)+(2*Fi(mm,nn)-Fi1(mm,nn)-Fi2(mm,nn))/(j*w*E);
    end
end
V=zeros(N,1);
V((N+1)/2,1)=1;
I_n=inv(Z)*V;%計算電流矩陣
I_temp=zeros(N+2,1);
for mm=1:N
    I_temp(mm+1,1)=I_n(mm,1);
end
mm=0:N+1;
plot(mm,abs(I_temp),'r');
grid off;
title('夾角天線的電流分布圖');
xlabel('分段數(shù)');
ylabel('電流')
Z_in=1/I_n((N+1)/2,1)
for mm=1:N   %計算各點到原點的距離
    L_n(mm,1)=((N+1)/2-mm)*div;
end
syms thet fai;
for mm=0:pi/10:pi
   fai=mm;
   for nn=1:N
       Enn(1,nn)=-j*w*u*exp(-j*k)*div*(cos(alpha)*cos(thet)*sin(fai)-sin(alpha)*sin(thet)+cos(alpha)*cos(fai))*exp(j*k*L_n(nn,1)*(cos(alpha)*sin(thet)*sin(fai)+sin(alpha)*sin(thet)))/(4*pi);
   end
   for nn=1:N
      E_n(1,nn)=subs(Enn(1,nn),{thet},pi/2);%計算電場中的極值點
  end
E_r=Enn*I_n;
E_max=E_n*I_n;%計算增益系數(shù)
G=4*pi*(abs(E_max))^2/(eta*real(Z_in)*(abs(I_n((N+1)/2,1)))^2);
F=abs(E_r)/abs(E_max);
end
figure;
ezpolar(F) %畫垂直方向性圖
title('夾角天線的E面方向性圖');