L=0.04;neff=1.45;C=10*1e-9;N=50;M=1501;Fai=pi;
deltaneff = 0.0001;
lamda1=1549;lamda2=1551;
lamda=linspace(lamda1,lamda2,M) * 1e-9;
deltalamda=(lamda2-lamda1)/M * 1e-9 ;
tic
for k=1:M
F=[1,0;0,1];
for i =1:N/2
% (1)均匀
deltaneff = 0.00005;
lamda_D =(1550 -C*L/2+C*i * L/N) * 1e-9;
sigma= 2 * pi * neff * (1/lamda(k) -1/lamda_D) + 2 * pi * deltaneff/lamda (k)+ (4 * pi * neff) * C * ( -L/2+i *L/N)/lamda_D^ 2;
kac = pi * deltaneff/lamda (k);
RB= sqrt(kac^2 -sigma^2);
F=F*[cosh(RB* L/N)-j*(sigma/RB)*sinh(RB*L/N),-j*(kac/RB)*sinh(RB*L/N);j*(kac/RB)*sinh(RB*L/N),cosh(RB*L/N)+j*(sigma/RB)*sinh(RB*L/N)];
end
F=F*[exp(-i*Fai/2) 0;0 exp(i*Fai/2)];
for i =N/2:N
% (1)均匀
deltaneff = 0.00005;
lamda_D =(1550 -C*L/2+C*i * L/N) * 1e-9;
sigma= 2 * pi * neff * (1/lamda(k) -1/lamda_D) + 2 * pi * deltaneff/lamda (k)+ (4 * pi * neff) * C * ( -L/2+i *L/N)/lamda_D^ 2;
kac = pi * deltaneff/lamda (k);
RB= sqrt(kac^2 -sigma^2);
F=F*[cosh(RB* L/N)-j*(sigma/RB)*sinh(RB*L/N),-j*(kac/RB)*sinh(RB*L/N);j*(kac/RB)*sinh(RB*L/N),cosh(RB*L/N)+j*(sigma/RB)*sinh(RB*L/N)];
end
R(k) = (abs( -F(3)/F(1)))^2;
Q(k)=phase(( -F(3)/F(1)));
end
toc
tao(1)=Q(1) ;tao(2) =Q(2) ;tao(3) =Q(3);
for i=4:M
if (abs(Q(i -1)-Q(i))<=1)
tao(i) = ((lamda1+i * 0.001)^2 * 1e-18/(2 *pi * 3e-4) * (Q(i-1)-Q(i))/deltalamda);
else
tao(i) = ((lamda1+i * 0.001)^ 2 * 1e-18/(2 * pi * 3e-4) * (Q(i-3) -Q(i-2))/deltalamda);
end
end
%%%%%%%%反射率图形%%%%%%%%%%
clf reset,h_ap = axes ('position',[0.13,0.13,0.7,0.75]);
set (h_ap,'xcolor','b','ycolor', 'b','xlim',[lamda1,lamda2],'ylim',[0,1]);
nx = 10;ny = 10;
pxtick = lamda1: ( (lamda2- lamda1) /nx) :lamda2 ;pytick = 0:((1 -0)/ny):1;
set (h_ap,'xtick',pxtick,'ytick',pytick ,'xgrid' ,'on','ygrid','on')
h_linet = line (lamda * 1e9,R,'color', 'b' );
set(get(h_ap,'xlabel'),'string','Wavelength/nm')
set(get(h_ap,'ylabel' ),'string','反射率')
%%%%%%% 时延图形%%%%%%%%%
h_at=axes ('position', get(h_ap ,'position' ));
set (h_at,'color', 'none', 'xcolor', 'b', 'ycolor','r' );
set (h_at ,'yaxislocation','right')
set(get(h_at,'ylabe' ) ,'string',' 时延' )
a = floor (min (tao ) ) ; b = ceil (max (tao));
set (h_at,'xlim',[lamda1,lamda2], 'ylim',[a,b])
line(lamda*1e9 ,tao ,'color','r','parent' ,h_at)
xpm=get(h_at,'xlim');
txtick=lamda1:((lamda2-lamda1)/nx) :lamda2 ;
tytick=a:((b-a)/ny):b;set(h_at,'xtick', txtick ,'ytick', tytick)