clc,format compact
AI=0.1; H=0.03;

FEI=[];   GN_C=[];   GN_A=[];   GN_R=[];   GN_CR=[];
AA=[];
Win=[9.255,12.71,12.21,1,1.074,0.8505,3.251,9.262,1,1.073];
Win=[0.1,6,12,20];

for W=Win
%for W=[0.05]
   Fei=[];   GN_c=[];   GN_a=[];   GN_r=[];   GN_cr=[];
   k=0;aa=[];
   while k<100
      A=W*AI/2*(1+0.0001)+k*W*AI/100;
      k=k+2;dt=2*pi/W/100;
      tt=[];t0=0.0; time=0;
      while (time<=4*pi/W)
         time=time+dt;
         f=A/W*(cos(W*t0)-cos(W*time));
         if abs(f)<0.9*AI
            dt=2*pi/W/100;
         elseif abs(f)<0.99*AI
            dt=2*pi/W/1000;
         elseif abs(f)<0.999*AI
            dt=2*pi/W/10000;
         elseif abs(f)<0.9999*AI
            dt=2*pi/W/100000;
         elseif abs(f)<AI
            dt=2*pi/W/1000000;
        end
         if f>=AI
            t0=time;
            if time>2*pi/W+2*pi/W/10000
               tt=[tt time];               
            end
         elseif f<=-AI
            t0=time;
            if time>2*pi/W+2*pi/W/10000               
               tt=[tt -time];
            end
         end
      end
      b1=sum(sign(tt).*(cos(W*tt)-cos(W*(tt+sign(tt).*H)))/pi);
      a1=sum((-sin(W*tt)+sin(W*(tt+sign(tt)*H)))/pi);
      N_c=sqrt(b1^2+a1^2)/A;
      N_a=W/2/A*H*length(tt)/2;
      N_r=1/A*sqrt(W*H*length(tt)/pi);
      N_cr=1/A*sqrt(((sin(W*H/2)*4/pi)^2+...
         (sin(3*W*H/2)*4/pi/3)^2)*length(tt));
      fei=atan2(a1,b1);
      Fei=[Fei;fei];
      GN_c=[GN_c;N_c];   GN_a=[GN_a;N_a];
      GN_r=[GN_r;N_r];   GN_cr=[GN_cr;N_cr];
      aa=[aa;A];
   disp(['** working.. Waiting.. W=', num2str(W),'  A=',num2str(A),' tt=',mat2str(tt)])
   end
   FEI=[FEI,Fei];
   GN_C=[GN_C,GN_c];   GN_A=[GN_A,GN_a];
   GN_R=[GN_R,GN_r];   GN_CR=[GN_CR,GN_cr];
   AA=[AA,aa];
end
%plot(FEI,GN_c)
%plot(Win,FEI)
plot(-FEI*57.3-180,-20*log10(GN_C)),hold
plot(-FEI(1,:)*57.3-180,-20*log10(GN_C(1,:)))
plot(-FEI*57.3-180,-20*log10(GN_A))
plot(-FEI(1,:)*57.3-180,-20*log10(GN_A(1,:)))
plot(-FEI*57.3-180,-20*log10(GN_CR))
plot(-FEI(1,:)*57.3-180,-20*log10(GN_CR(1,:)))
for nn=1:length(Win),text(-FEI(1,nn)*57.3-180,-20*log10(GN_C(1,nn)),num2str(Win(nn))),end
hold off
save ipfdf1 Win FEI GN_C GN_A GN_R GN_CR AA