clear all;
n=2;%彎折次數(shù)
lanbta=0.3/1.07;%波長(zhǎng)
a=0.0005;%天線半徑
l=0.06;%天線外部長(zhǎng)度一半
h=0.005;%彎折高度
sn=l/(2*n+1); 
fanwei1=atan(2*h/sn)*180/pi;
fanwei2=180-fanwei1;
zr=90;%彎折角度

z=(zr/180)*pi;
 if(atan(2*h/sn)<= z && z<=(pi-atan(2*h/sn)))
        L=(n*(2*sn+(2*h)/sin(z)-(2*h)/tan(z))+sn)*2;%天線總長(zhǎng)度
        pou=150;
        if(ceil((ceil(L/(lanbta/pou)))/2)==(ceil(L/(lanbta/pou)))/2)
           m=ceil(L/(lanbta/pou));
        else
           m=ceil(L/(lanbta/pou))+1;%分成的段數(shù),為偶數(shù)
        end
           
     ol=L/m;%每段長(zhǎng)度 
     for v=1:(m-1);
         if (v==m/2)
            u(v)=1;
         else
            u(v)=0;
         end
     end
  for  q=1:(m-1);
      for p=1:(m-1);
        Xq=(q-m/2)*ol;
        Xp=(p-m/2)*ol;
        Xq1=(q-m/2)*ol+ol/2;
        Xq2=(q-m/2)*ol-ol/2;
        Xp1=(p-m/2)*ol+ol/2;
        Xp2=(p-m/2)*ol-ol/2;%未彎折之前的坐標(biāo)
if(Xq>=0)
kq=floor(Xq/(2*sn+2*h/sin(z)-2*h/tan(z)));
else
    kq=ceil(Xq/(2*sn+2*h/sin(z)-2*h/tan(z)));
end
if(Xp>=0)
kp=floor(Xp/(2*sn+2*h/sin(z)-2*h/tan(z)));
else
 kp=ceil(Xp/(2*sn+2*h/sin(z)-2*h/tan(z)));
end
if(Xq1>=0)
kq1=floor(Xq1/(2*sn+2*h/sin(z)-2*h/tan(z)));
else
 kq1=ceil(Xq1/(2*sn+2*h/sin(z)-2*h/tan(z)));
end
if(Xq2>=0)
kq2=floor(Xq2/(2*sn+2*h/sin(z)-2*h/tan(z)));
else
kq2=ceil(Xq2/(2*sn+2*h/sin(z)-2*h/tan(z))); 
end
if(Xp1>=0)
kp1=floor(Xp1/(2*sn+2*h/sin(z)-2*h/tan(z)));
else
kp1=ceil(Xp1/(2*sn+2*h/sin(z)-2*h/tan(z)));
end
if(Xp2>=0)
kp2=floor(Xp2/(2*sn+2*h/sin(z)-2*h/tan(z)));
else
    kp2=ceil(Xp2/(2*sn+2*h/sin(z)-2*h/tan(z)));
end
location= abs(Xq-kq*(2*sn+2*h/sin(z)-2*h/tan(z)));
if(0<=location&&location<sn)
    xq=Xq-kq*(2*h/sin(z)-2*h/tan(z)) ;
    yq=0;
elseif(sn<=location&&location<sn+h/sin(z))
    if(q>=(m/2))
       xq=(2*kq+1)*sn+(Xq-sn-kq*(2*sn+2*h/sin(z)-2*h/tan(z)))*cos(z);
       yq=(Xq-sn-kq*(2*sn+2*h/sin(z)-2*h/tan(z)))* sin(z);
    else
       xq=(2*kq-1)*sn+(Xq-sn-kq*(2*sn+2*h/sin(z)-2*h/tan(z)))*cos(z);
       yq=(-Xq-sn+kq*(2*sn+2*h/sin(z)-2*h/tan(z)))*sin(z);
    end
elseif(sn+h/sin(z)<=location&&location<2*sn+h/sin(z)-2*h/tan(z))
    if(q>=(m/2))
       xq=Xq+(2*kq*h+h)*(1/tan(z)-1/sin(z));
       yq=h;
    else
        xq=Xq+(2*kq*h-h)*(1/tan(z)-1/sin(z));
         yq=h;
    end      
elseif(2*sn+h/sin(z)-2*h/tan(z)<=location&&location<2*sn+2*h/sin(z)-2*h/tan(z))
    if(q>=(m/2))
        xq=(2*kq+2)*sn-h/tan(z)+(Xq-kq*(2*sn+2*h/sin(z)-2*h/tan(z))-2*sn-h/sin(z)+2*h/tan(z))*cos(z);
        yq=h-(Xq-kq*(2*sn+2*h/sin(z)-2*h/tan(z))-2*sn-h/sin(z)+2*h/tan(z))*sin(z);
    else
        xq=(2*kq-2)*sn+h/tan(z)+(Xq-kq*(2*sn+2*h/sin(z)-2*h/tan(z))+2*sn+h/sin(z)-2*h/tan(z))*cos(z);
        yq=h+(Xq-kq*(2*sn+2*h/sin(z)-2*h/tan(z))+2*sn+h/sin(z)-2*h/tan(z))*sin(z);
    end
end

location=abs(Xp-kp*(2*sn+2*h/sin(z)-2*h/tan(z)));
if(0<=location&&location<sn)
    xp=Xp-kp*(2*h/sin(z)-2*h/tan(z)) ;
    yp=0;
elseif(sn<=location&&location<sn+h/sin(z))
    if(p>=m/2)
    xp=(2*kp+1)*sn+(Xp-sn-kp*(2*sn+2*h/sin(z)-2*h/tan(z)))*cos(z);
    yp=(Xp-sn-kp*(2*sn+2*h/sin(z)-2*h/tan(z)))* sin(z);
    else
     xp=(2*kp-1)*sn+(Xp-sn-kp*(2*sn+2*h/sin(z)-2*h/tan(z)))*cos(z);
       yp=(-Xp-sn+kp*(2*sn+2*h/sin(z)-2*h/tan(z)))*sin(z);
    end
        
elseif(sn+h/sin(z)<=location&&location<2*sn+h/sin(z)-2*h/tan(z))
    if(p>=m/2)
    xp=Xp+(2*kp*h+h)*(1/tan(z)-1/sin(z));
    yp=h;
    else
        xp=Xp+(2*kp*h-h)*(1/tan(z)-1/sin(z));
        yp=h;
    end
        
elseif(2*sn+h/sin(z)-2*h/tan(z)<=location&&location<2*sn+2*h/sin(z)-2*h/tan(z)) 
    if(p>=m/2)
    xp=(2*kp+2)*sn-h/tan(z)+(Xp-kp*(2*sn+2*h/sin(z)-2*h/tan(z))-2*sn-h/sin(z)+2*h/tan(z))*cos(z);
    yp=h-(Xp-kp*(2*sn+2*h/sin(z)-2*h/tan(z))-2*sn-h/sin(z)+2*h/tan(z))*sin(z);
    else
          xp=(2*kp-2)*sn+h/tan(z)+(Xp-kp*(2*sn+2*h/sin(z)-2*h/tan(z))+2*sn+h/sin(z)-2*h/tan(z))*cos(z);
        yp=h+(Xp-kp*(2*sn+2*h/sin(z)-2*h/tan(z))+2*sn+h/sin(z)-2*h/tan(z))*sin(z);
    end
        
end

location=abs(Xq1-kq1*(2*sn+2*h/sin(z)-2*h/tan(z)));
if(0<=location&&location<sn)
    xq1=Xq1-kq1*(2*h/sin(z)-2*h/tan(z)) ;
    yq1=0;
elseif(sn<=location&&location<sn+h/sin(z))
      if(q>=m/2)
    xq1=(2*kq1+1)*sn+(Xq1-sn-kq1*(2*sn+2*h/sin(z)-2*h/tan(z)))*cos(z);
    yq1=(Xq1-sn-kq1*(2*sn+2*h/sin(z)-2*h/tan(z)))* sin(z);
      else
           xq1=(2*kq1-1)*sn+(Xq1-sn-kq1*(2*sn+2*h/sin(z)-2*h/tan(z)))*cos(z);
       yq1=(-Xq1-sn+kq1*(2*sn+2*h/sin(z)-2*h/tan(z)))*sin(z);
      end
          
elseif(sn+h/sin(z)<=location&&location<2*sn+h/sin(z)-2*h/tan(z))
      if(q>=m/2)
    xq1=Xq1+(2*kq1*h+h)*(1/tan(z)-1/sin(z));
    yq1=h;
      else
          xq1=Xq1+(2*kq1*h-h)*(1/tan(z)-1/sin(z));
         yq1=h;
      end
elseif(2*sn+h/sin(z)-2*h/tan(z)<=location&&location<2*sn+2*h/sin(z)-2*h/tan(z)) 
    if(q>=m/2)
    xq1=(2*kq1+2)*sn-h/tan(z)+(Xq1-kq1*(2*sn+2*h/sin(z)-2*h/tan(z))-2*sn-h/sin(z)+2*h/tan(z))*cos(z);
    yq1=h-(Xq1-kq1*(2*sn+2*h/sin(z)-2*h/tan(z))-2*sn-h/sin(z)+2*h/tan(z))*sin(z);
    else
        xq1=(2*kq1-2)*sn+h/tan(z)+(Xq1-kq1*(2*sn+2*h/sin(z)-2*h/tan(z))+2*sn+h/sin(z)-2*h/tan(z))*cos(z);
        yq1=h+(Xq1-kq1*(2*sn+2*h/sin(z)-2*h/tan(z))+2*sn+h/sin(z)-2*h/tan(z))*sin(z);
    end
        
end

location=abs(Xq2-kq2*(2*sn+2*h/sin(z)-2*h/tan(z)));
if(0<=location&&location<sn)
    xq2=Xq2-kq2*(2*h/sin(z)-2*h/tan(z)) ;
    yq2=0;
elseif(sn<=location&&location<sn+h/sin(z))
    if(q>=m/2)
    xq2=(2*kq2+1)*sn+(Xq2-sn-kq2*(2*sn+2*h/sin(z)-2*h/tan(z)))*cos(z);
    yq2=(Xq2-sn-kq2*(2*sn+2*h/sin(z)-2*h/tan(z)))* sin(z);
    else
        xq2=(2*kq2-1)*sn+(Xq2-sn-kq2*(2*sn+2*h/sin(z)-2*h/tan(z)))*cos(z);
       yq2=(-Xq2-sn+kq2*(2*sn+2*h/sin(z)-2*h/tan(z)))*sin(z);
    end
        
elseif(sn+h/sin(z)<=location&&location<2*sn+h/sin(z)-2*h/tan(z))
       if(q>=m/2)
    xq2=Xq2+(2*kq2*h+h)*(1/tan(z)-1/sin(z));
    yq2=h;
       else
            xq2=Xq2+(2*kq2*h-h)*(1/tan(z)-1/sin(z));
         yq2=h;
       end
    
elseif(2*sn+h/sin(z)-2*h/tan(z)<=location&&location<2*sn+2*h/sin(z)-2*h/tan(z)) 
      if(q>=m/2)
    xq2=(2*kq2+2)*sn-h/tan(z)+(Xq2-kq2*(2*sn+2*h/sin(z)-2*h/tan(z))-2*sn-h/sin(z)+2*h/tan(z))*cos(z);
    yq2=h-(Xq2-kq2*(2*sn+2*h/sin(z)-2*h/tan(z))-2*sn-h/sin(z)+2*h/tan(z))*sin(z);
      else
          xq2=(2*kq2-2)*sn+h/tan(z)+(Xq2-kq2*(2*sn+2*h/sin(z)-2*h/tan(z))+2*sn+h/sin(z)-2*h/tan(z))*cos(z);
        yq2=h+(Xq2-kq2*(2*sn+2*h/sin(z)-2*h/tan(z))+2*sn+h/sin(z)-2*h/tan(z))*sin(z);
      end
          
end

location=abs(Xp1-kp1*(2*sn+2*h/sin(z)-2*h/tan(z)));
if(0<=location&&location<sn)
    xp1=Xp1-kp1*(2*h/sin(z)-2*h/tan(z)) ;
    yp1=0;
elseif(sn<=location&&location<sn+h/sin(z))
    if(p>=m/2)
    xp1=(2*kp1+1)*sn+(Xp1-sn-kp1*(2*sn+2*h/sin(z)-2*h/tan(z)))*cos(z);
    yp1=(Xp1-sn-kp1*(2*sn+2*h/sin(z)-2*h/tan(z)))* sin(z);
    else
        xp1=(2*kp1-1)*sn+(Xp1-sn-kp1*(2*sn+2*h/sin(z)-2*h/tan(z)))*cos(z);
       yp1=(-Xp1-sn+kp1*(2*sn+2*h/sin(z)-2*h/tan(z)))*sin(z);
    end
        
elseif(sn+h/sin(z)<=location&&location<2*sn+h/sin(z)-2*h/tan(z))
     if(p>=m/2)
    xp1=Xp1+(2*kp1*h+h)*(1/tan(z)-1/sin(z));
    yp1=h;
     else
          xp1=Xp1+(2*kp1*h-h)*(1/tan(z)-1/sin(z));
        yp1=h;
     end
         
elseif(2*sn+h/sin(z)-2*h/tan(z)<=location&&location<2*sn+2*h/sin(z)-2*h/tan(z)) 
    if(p>=m/2)
    xp1=(2*kp1+2)*sn-h/tan(z)+(Xp1-kp1*(2*sn+2*h/sin(z)-2*h/tan(z))-2*sn-h/sin(z)+2*h/tan(z))*cos(z);
    yp1=h-(Xp1-kp1*(2*sn+2*h/sin(z)-2*h/tan(z))-2*sn-h/sin(z)+2*h/tan(z))*sin(z);
    else
        xp1=(2*kp1-2)*sn+h/tan(z)+(Xp1-kp1*(2*sn+2*h/sin(z)-2*h/tan(z))+2*sn+h/sin(z)-2*h/tan(z))*cos(z);
        yp1=h+(Xp1-kp1*(2*sn+2*h/sin(z)-2*h/tan(z))+2*sn+h/sin(z)-2*h/tan(z))*sin(z);
    end
end

location=abs(Xp2-kp2*(2*sn+2*h/sin(z)-2*h/tan(z)));
 if(0<=location&&location<sn)
    xp2=Xp2-kp2*(2*h/sin(z)-2*h/tan(z)) ;
    yp2=0;
 elseif(sn<=location&&location<sn+h/sin(z))
     if(p>=m/2)
    xp2=(2*kp2+1)*sn+(Xp2-sn-kp2*(2*sn+2*h/sin(z)-2*h/tan(z)))*cos(z);
    yp2=(Xp2-sn-kp2*(2*sn+2*h/sin(z)-2*h/tan(z)))* sin(z);
     else
         xp2=(2*kp2-1)*sn+(Xp2-sn-kp2*(2*sn+2*h/sin(z)-2*h/tan(z)))*cos(z);
       yp2=(-Xp2-sn+kp2*(2*sn+2*h/sin(z)-2*h/tan(z)))*sin(z);
     end
 elseif(sn+h/sin(z)<=location&&location<2*sn+h/sin(z)-2*h/tan(z))
     if(p>=m/2)
    xp2=Xp2+(2*kp2*h+h)*(1/tan(z)-1/sin(z));
    yp2=h;
     else
          xp2=Xp2+(2*kp2*h-h)*(1/tan(z)-1/sin(z));
        yp2=h;
     end
         
 elseif(2*sn+h/sin(z)-2*h/tan(z)<=location&&location<2*sn+2*h/sin(z)-2*h/tan(z)) 
      if(p>=m/2)
    xp2=(2*kp2+2)*sn-h/tan(z)+(Xp2-kp2*(2*sn+2*h/sin(z)-2*h/tan(z))-2*sn-h/sin(z)+2*h/tan(z))*cos(z);
    yp2=h-(Xp2-kp2*(2*sn+2*h/sin(z)-2*h/tan(z))-2*sn-h/sin(z)+2*h/tan(z))*sin(z);
      else
           xp2=(2*kp2-2)*sn+h/tan(z)+(Xp2-kp2*(2*sn+2*h/sin(z)-2*h/tan(z))+2*sn+h/sin(z)-2*h/tan(z))*cos(z);
        yp2=h+(Xp2-kp2*(2*sn+2*h/sin(z)-2*h/tan(z))+2*sn+h/sin(z)-2*h/tan(z))*sin(z);
      end
          
 end 

dl=(xq2-xq1)*(xp2-xp1)+(yq2-yq1)*(yp2-yp1);
   r=sqrt((xq-xp)^2+(yq-yp)^2+a^2);
r1=sqrt((xq1-xp1)^2+(yq1-yp1)^2+a^2);
r2=sqrt((xq1-xp2)^2+(yq1-yp2)^2+a^2);
r3=sqrt((xq2-xp1)^2+(yq2-yp1)^2+a^2);
r4=sqrt((xq2-xp2)^2+(yq2-yp2)^2+a^2);

if(p==q)

    h0=1/(2*pi*ol)*log(ol/a)-i/(2*lanbta);
    h1=1/(2*pi*ol)*log(ol/a)-i/(2*lanbta);
    h4=1/(2*pi*ol)*log(ol/a)-i/(2*lanbta);
    h2=exp(-i*2*pi*r2/lanbta)/(4*pi*r2);    
h3=exp(-i*2*pi*r3/lanbta)/(4*pi*r3);

elseif(q-p==1)
    h3=1/(2*pi*ol)*log(ol/a)-i/(2*lanbta);
     h0=exp(-i*2*pi*r/lanbta)/(4*pi*r);
   h1=exp(-i*2*pi*r1/lanbta)/(4*pi*r1);
h2=exp(-i*2*pi*r2/lanbta)/(4*pi*r2);    
h4=exp(-i*2*pi*r4/lanbta)/(4*pi*r4); 
elseif(q-p==-1)
    h2=1/(2*pi*ol)*log(ol/a)-i/(2*lanbta);
     h0=exp(-i*2*pi*r/lanbta)/(4*pi*r);
   h1=exp(-i*2*pi*r1/lanbta)/(4*pi*r1); 
h3=exp(-i*2*pi*r3/lanbta)/(4*pi*r3);
h4=exp(-i*2*pi*r4/lanbta)/(4*pi*r4); 
else
   h0=exp(-i*2*pi*r/lanbta)/(4*pi*r);
   h1=exp(-i*2*pi*r1/lanbta)/(4*pi*r1);
h2=exp(-i*2*pi*r2/lanbta)/(4*pi*r2);    
h3=exp(-i*2*pi*r3/lanbta)/(4*pi*r3);
h4=exp(-i*2*pi*r4/lanbta)/(4*pi*r4); 
end

t(q,p)=i*240*pi*pi*dl*h0/lanbta-i*60*lanbta*(h1-h2-h3+h4);
sol(q)=(xq1-xq2)+i*(yq1-yq2);
srn(q)=xq+i*yq;

      end
  end
  y=inv(t);
 Zin=1/y(m/2,m/2)
 Pin=real(y(m/2,m/2))*u(m/2)*u(m/2)
 S=(Zin-50)/(Zin+50)
 zbb=(1+abs(S))/(1-abs(S))
c=u*y;
Emax1=0;
Emax2=0;
    for fan=0:(10/180)*pi:2*pi
        for sit=0:(10/180)*pi:2*pi
        Esit=0;
        Efan=0;
         for q=1:(m-1)
        Esit=Esit+dot(sol(q),cos(sit)*cos(fan)+i*cos(sit)*sin(fan))*exp(i*2*pi*dot(srn(q),sin(sit)*cos(fan)+i*sin(sit)*sin(fan))/lanbta)*c(q);
        Efan=Efan+dot(sol(q),-sin(fan)+i*cos(fan))*exp(i*2*pi*dot(srn(q),sin(sit)*cos(fan)+i*sin(sit)*sin(fan))/lanbta)*c(q);
         end
         E1=abs(Esit)^2;
         E2=abs(Efan)^2;
         if(E1>Emax1)
              Emax1=E1;
         end
         if(E2>Emax2)
             Emax2=E2;
         end
        end
    end
    G1=120*pi*pi*sqrt(Emax1^2+Emax2^2)/(lanbta*Pin);
    G=10*log(G1)
    
     b =1;
    for sit=0:(10/180)*pi:2*pi
        A=0;
         for q=1:(m-1)
          fan=0;  
        A=A+dot(sol(q),cos(sit)*cos(fan)+i*cos(sit)*sin(fan))*exp(i*2*pi*dot(srn(q),sin(sit)*cos(fan)+i*sin(sit)*sin(fan))/lanbta)*c(q);
         end
        E(b)=abs(A);
        b=b+1;
    end
     sit=0:(10/180)*pi:2*pi; 
     Smax=max(E);
     E=E/Smax;
   subplot(2,2,1); polar(sit,E);title('E面方向圖');
    b =1;
    for fan=0:(10/180)*pi:2*pi
        A=0;
         for q=1:(m-1)
           sit=0;  
        A=A+dot(sol(q),-sin(fan)+i*cos(fan))*exp(i*2*pi*dot(srn(q),sin(sit)*cos(fan)+i*sin(sit)*sin(fan))/lanbta)*c(q);
         end
        E(b)=abs(A);
        b=b+1;
    end
      fan=0:(10/180)*pi:2*pi; 
      Fmax=max(E);
      E=E/Fmax;
   subplot(2,2,2) ;polar(fan,E);title('H面方向圖');
cmax=max(abs(c));
        for k=1:(m-1);  %unit current matrix calculation
            c(k)=abs(c(k)/cmax);
        end

        k=1:(m-1);
      subplot(2,2,3);  plot(ol*(k-m/2),c(k),'k:p');title('電流分布');
        hold on
 else
     fanwei1=atan(2*h/sn)*180/pi  
     fanwei2=180-fanwei1
   disp('角度只能在fanwei1=argtan(2h/sn)與180-arg(tan2h/sn)之間變換');
 end