L=0.01;
neff=1.447;
c=2*1e-7/(2*neff);
deltaneff=0.0002;

M=2000;
lamda1=1546;
lamda2=1556;
lamda=linspace(lamda1*1e-9,lamda2*1e-9,M);
deltalamda=(lamda2*1e-9-lamda1*1e-9)/(M-1);

 global k;  


for k=1:M
deltaneff_z=deltaneff;
  %g=4;
   % deltaneff_z=deltaneff*exp(-g*g*((2*n-1-N)/(2*N))^2);    deltaneff_zh
  %  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*(n-1)*h;
   
    %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*h+h)/lamdaD^2;
    
    %kac=pi*deltaneff_z/lamda(k);
[t,x]=ode45('xprim',[0.5*L,-0.5*L],[1;0]);
R=x(:,1);
S=x(:,2);




    
  z(k)=S(end)/(R(end)+eps);
   r(k)=(abs(z))^2;
   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; 


subplot(1,2,1)

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





subplot(1,2,2)
plot(lamda1,tao,'k')
title('時延')


function xprim=xprim(t,x)
    neff=1.447;
    lamdaD=1550*1e+2*neff*c*t;
    c=2*1e-7/(2*neff);
lamda1=1546;
lamda2=1556;
   lamda=linspace(lamda1*1e-9,lamda2*1e-9,2000);
deltaneff=0.0002;
deltaneff_z=deltaneff;
    sigma=2*pi*neff*(1/lamda(k)-1/lamdaD)+2*pi*deltaneff_z/lamda(k)+(4*pi*neff)*c*(2*(neff+deltaneff))*t/lamdaD^2;
    
    kac=pi*deltaneff_z/lamda(k);

xprim=[1j*sigma*x(1)+1j*kac*x(2);-1j*sigma*x(2)-1j*kac*x(1)];