Pi=3.142
L=0.04;
Neff=1.45;
c=10*1e-9;
deltaN=0.0001;
SectionN=50;
Lamda=[1549.5:0.001:1550.6]*1e-9;
%计算波长范围
for k=1:1101 %1100=(1549.5-1550.6)/0.001 为总计算点数
F=[1,0;0,1]; %初始化传输矩阵
for i=1:SectionN %计算相应波长下的传输矩阵
%deltaN=0.00005*exp(-(56*(-L/2+i*L/SectionN)^4)/(L)^4);%切趾
LamdaD=(1550-c*L/2+c*i*L/SectionN)*1e-9; %计算每段光栅本地中心波长
delta=2*Pi*Neff*(1/Lamda(k)-1/LamdaD)+2*Pi*deltaN/Lamda(k)+(4*3.142*Neff)*c*(-L/2+i*L/SectionN)/(LamdaD^2);
Kac=3.142*deltaN/Lamda(k);
RB=(Kac^2-delta^2)^0.5;
F=F*[cosh(RB*L/SectionN)-j*(delta/RB)*sinh(RB*L/SectionN),-j*(Kac/RB)*sinh(RB*L/SectionN);j*(Kac/RB)*sinh(RB*L/SectionN),cosh(RB*L/SectionN)+j*(delta/RB)*sinh(RB*L/SectionN)];
end
Y1(k)=(abs(-F(3)/F(1)))^2;
Q(k)=phase(-F(3)/F(1));
end
Y(1)=Q(1);Y(2)=Q(2);Y(3)=Q(3);
for i=4:1101
if(abs(Q(i-1)-Q(i))<=1) %对相位不连续处求导的简单处理
Y(i)=(1549.5+i*0.001)^2*1e-18/(2*Pi*3e-4)*((Q(i-1)-Q(i))/(0.001e-9));
else
Y(i)=(1549.5+i*0.001)^2*1e-18/(2*Pi*3e-4)*((Q(i-3)-Q(i-2))/(0.001e-9));
end
end
hold on;
figure(1)
plot(Lamda*1e9,Y,'b'); %绘制时延曲线
figure(2)
plot(Lamda*1e9,Y1*500,'b'); %绘制放大500倍的反射谱
