L=0.01;
neff=1.447;
c=1*1e-7/(2*neff);
deltaneff=0.0004;
N=100;
M=2000;
lamda1=1548;
lamda2=1552;
lamda=linspace(lamda1*1e-9,lamda2*1e-9,M);
deltalamda=(lamda2*1e-9-lamda1*1e-9)/(M-1);
for k=1:M
    F=[1,0;0,1];
    for n=1:N

deltaneff_z=deltaneff;
  %g=4;
   % deltaneff_z=deltaneff*exp(-g*g*((2*n-1-N)/(2*N))^2);
  %  H=0.5;
  %deltaneff_z=deltaneff*(1+H*cos(pi*(2*n-1-N)/N))/(1+H);
    lamdaD=1550*1e-9-L*neff*c+2*neff*c*(L*n/N-0.5*L/N);
    sigma=2*pi*neff*(1/lamda(k)-1/lamdaD)+2*pi*deltaneff_z/lamda(k)+(4*pi*neff)*c*(2*(neff+deltaneff))*(-L/2+n*L/N)/lamdaD^2;
    kac=pi*deltaneff_z/lamda(k);
    RB=sqrt(kac^2-sigma^2)+eps;
    F=F*[cosh(RB*L/N)-(sigma/RB)*sinh(RB*L/N)*1i,-(kac/RB)*sinh(RB*L/N)*1i;(kac/RB)*sinh(RB*L/N)*1i,cosh(RB*L/N)+(sigma/RB)*sinh(RB*L/N)*1i];
    end
    R(k)=(abs(F(2,1)/(F(1,1)+eps)))^2;
   z=F(2,1)/(F(1,1)+eps);
   Q(k)=angle(z);

    if Q(k)<0                      %angle0-2pi
         Q(k)=2*pi+Q(k);
         end  
 
end
  for  a=2:M-1
     if   abs(Q(a)-Q(a+1))>=pi
       dQ(a)=(Q(a)-Q(a-1))/deltalamda;

      else
          dQ(a)=(Q(a+1)-Q(a))/deltalamda;
     end
end
dQ(1)=dQ(2);
lamda1=lamda(1:M-1);
tao=(-1/(2*pi*3*1e8))*(lamda1.*lamda1).*dQ; 

hold on;
%subplot(1,2,1)

plot(lamda,R,'k')
title('反射譜')

hold off

%{
hold on;
%subplot(1,2,2)
plot(lamda1,tao,'r')
title('時延')
hold off
%}