clc;
clear;
l_o=0.453;  
l_r=0.475;  
l_1=0.446;  
l_2=0.446;  
n=9;  
w=3*10^8*2*pi;   
E=8.854*10^(-12);   
u=4*pi*10^(-7);  
a=5*10^(-3);   
d_r=0.3;  
d_o=0.3;  
d_l=0.3;  
div_o=l_o/(n+1);    
div_r=l_r/(n+1);     
div_1=l_1/(n+1);     
div_2=l_2/(n+1);     
k=2*pi;        
ita=120*pi;         
F1=@zuobiao;
zb=feval(F1,0);       %各段中點的坐標
zb1=feval(F1,1);      %各段+點的坐標
zb2=feval(F1,-1);     %各段-點的坐標
for p=1:4*n           %計算p,q兩點之間的距離
    for q=1:4*n
        R(p,q)=sqrt((zb(p,1)-zb(q,1))^2+(zb(p,2)-zb(q,2))^2+(zb(p,3)-zb(q,3))^2);
    end
end
for p=1:4*n          %計算p+,q-之間的距離
    for q=1:4*n
        R1(p,q)=sqrt((zb1(p,1)-zb2(q,1))^2+(zb1(p,2)-zb2(q,2))^2+(zb1(p,3)-zb2(q,3))^2);
    end
end
for p=1:4*n
    for q=1:4*n
         R2(p,q)=sqrt((zb2(p,1)-zb1(q,1))^2+(zb2(p,2)-zb1(q,2))^2+(zb2(p,3)-zb1(q,3))^2);
     end
 end
 F2=@zmn;
 S00=feval(F2,0,0,div_r,R,R1,R2);           %計算各個分塊阻抗矩陣
 S11=feval(F2,1,1,div_o,R,R1,R2); 
 S22=feval(F2,2,2,div_1,R,R1,R2); 
 S33=feval(F2,3,3,div_2,R,R1,R2);
 F3=@zmn1;
 S01=feval(F3,0,1,div_r,div_o,R,R1,R2);S02=feval(F3,0,2,div_r,div_1,R,R1,R2);
 S03=feval(F3,0,3,div_r,div_2,R,R1,R2);S12=feval(F3,1,2,div_o,div_1,R,R1,R2);
 S13=feval(F3,1,3,div_o,div_2,R,R1,R2); S23=feval(F3,2,3,div_1,div_2,R,R1,R2);
 S10=S01.'; S20=S02.';S30=S03.'; S21=S12.'; S31=S13.'; S32=S23.';
 S=[S00,S01,S02,S03;S10,S11,S12,S13;S20,S21,S22,S23;S30,S31,S32,S33];

 P=4*S;
 [C,lam]=eig(P);
 [D,mumu]=eig(conj(P));
 lamda=eig(P);
 mu=eig(conj(P));
 uu=C';
 vv=D';

 
 for p=1:4*n             %對電壓矩陣賦值
     if p==(n+(n+1)/2)
         V(p,1)=1;
     else
         V(p,1)=0;
     end
 end
 for p=1:4*n 
     AA(p,:)=dot(vv(p,:),V)/dot(vv(p,:),uu(p,:))/mu(p,:)*uu(p,:);
 end
A=zeros(4*n,4*n);

A(1,:)=AA(1,:);
for p=2:1:4*n
    A(p,:)=AA(p,:)+A(p-1,:);
end
     
I_n=inv(4*S)*V;%求電流矩陣
ss=1:4*n;
plot(ss,abs(I_n),'k',ss,abs(A(1,:)),'b*',ss,abs(A(2,:)),'g*',ss,abs(A(4*n,:)),'r*');
xlabel('分段數') ;
ylabel('電流大小 A');  
title('四元八木天線電流分布曲線');
legend('矩量法','一次模','前二次模加','前4*n次模加');


I1_n=A(4*n,:);              %把4*n修改
ii=n+(n+1)/2;
Z1_in=conj(1/I1_n(ii))          %正交模法輸入阻抗
Z_in=1/I_n(ii)           %矩量法輸入阻抗

syms Theta;
syms Fi;
 for q=1:n
     E_n(1,q)=i*w*u*sin(Theta)*exp(i*k*((-d_r)*sin(Theta)*cos(Fi)+zb(q,3)*cos(Theta)))*div_r/(4*pi);
 end
 for q=(n+1):2*n;
     E_n(1,q)=i*w*u*sin(Theta)*exp(i*k*zb(q,3)*cos(Theta))*div_o/(4*pi);
 end
 for q=(2*n+1):3*n
     E_n(1,q)=i*w*u*sin(Theta)*exp(i*k*((d_o)*sin(Theta)*cos(Fi)+zb(q,3)*cos(Theta)))*div_1/(4*pi);
 end
 for q=(3*n+1):4*n
      E_n(1,q)=i*w*u*sin(Theta)*exp(i*k*((d_o+d_l)*sin(Theta)*cos(Fi)+zb(q,3)*cos(Theta)))*div_2/(4*pi);
 end




  E=E_n*I_n;             %矩量法計算增益系數
  Emax=abs(subs(E,{Theta,Fi},{pi/2,0})); 
  G=4*pi*(Emax^2)/(ita*real(Z_in)*abs(I_n(14))^2)
  %G_db=10*log10(G)
  F=abs(E)/Emax;
  F1=subs(F,'Fi',0);
  F2=subs(F,'Theta',pi/2);
  h=0:2*pi/100:2*pi;
  F1=subs(F1,'Theta',h);
  F2=subs(F2,'Fi',h);
 
 
  E1=E_n*I1_n';             %正交模計算增益系數
  Emax=abs(subs(E1,{Theta,Fi},{pi/2,0})); 
  G=4*pi*(Emax^2)/(ita*real(Z1_in)*abs(I1_n(14))^2)
  H=abs(E1)/Emax;
  H1=subs(H,'Fi',0);
  H2=subs(H,'Theta',pi/2);
  H1=subs(H1,'Theta',h);
 H2=subs(H2,'Fi',h);
  
  figure
  polar(h,F1,'r*');
  hold on
  polar(h,H1,'b');
  grid on
  title('E面方向圖');
  legend('矩量法','雙正交模法');
  view(90,-90)
  figure
  polar(h,F2,'r*');
  hold on
  polar(h,H2,'b');
  grid on
  title('H面方向圖');
  legend('矩量法','雙正交模法');

 
% figure;
%x=F*sin(Theta)*cos(Fi);  
%y=F*sin(Theta)*sin(Fi);  
%z=F*cos(Theta);
%  ezsurf(x,y,z,120) ;
 % title('立體方向圖')