clear all;
n=2;%彎折次數
a=0.0005;%天線半徑
l=0.06;%天線外部長度一半
h=0.015;%彎折高度
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;%天線總長度
        for k=1:201
            f=0.75+0.001*(k-1);
          lanbta=0.3/f;
          m=50;
     ol=L/m;%每段長度 
    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;%未彎折之前的坐標
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);
r1=sqrt((xq1-xp1)^2+(yq1-yp1)^2);
r2=sqrt((xq1-xp2)^2+(yq1-yp2)^2);
r3=sqrt((xq2-xp1)^2+(yq2-yp1)^2);
r4=sqrt((xq2-xp2)^2+(yq2-yp2)^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);
      end
  end
  y=inv(t);
 Zin(k)=1/y(m/2,m/2);
 S(k)=(Zin(k)-50)/(Zin(k)+50);
 zbb(k)=(1+abs(S(k)))/(1-abs(S(k)));
        end
        f=0.65:0.001:0.85;
         plot(f,zbb);
       
          grid minor
         else
     fanwei1=atan(2*h/sn)*180/pi  
     fanwei2=180-fanwei1
   disp('角度只能在fanwei1=argtan(2h/sn)與180-arg(tan2h/sn)之間變換');
 end