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)];