ou_tangshucheng_xiaoyongliang.m

用于模拟光纤光栅(FBG)的反射特性和时延特性,特别是用于模拟各种切趾包络下光纤光栅的反射和时延特性。

%布拉格光栅谱特性模拟
clc
clear all
format long
%定义光栅长度,有效折射率,周秋系数,光栅所分段数,波长选点数
%L=0.04;neff=1.45;C=10.*1e-9;N=50;M=1501;
%L=0.01;neff=1.45;C=0;N=50;M=1501;         %T.ERDOGAN光栅光谱中图(7)a数据
L=0.100;neff=1.45;C=2.2*1e-9;N=100;M=3201;         %刘勇图1数据，其中C=delta_lamda_D/delta_z,为啁啾系数
F=0.48;%T1=100;  %QINZIXIONG切趾包络函数（1）数据
lamda1=1549.5;lamda2=1550.5;%计算波长范围
lamda=linspace(lamda1,lamda2,M)*1e-9;
deltalamda=(lamda2-lamda1)/M*1e-9;
for k=1:M        %M为总计算点数，即分别求出对应M个波长反射率和时延，使这M个波长组成一段波长范围[lamda1,lamda2]
    F=[1,0;0,1]; 
    for i=1:N    %计算相应波长lamda(k)下的传输矩阵
        Z(i)=-L/2+i*L/N;
        %(1)均匀p        %deltaneff=10*1e-5;%折射率调制深度         汤树成
        %deltaneff=2*1e-5;
        %(2)超高斯包络f(z)=exp(-alpha.Z^(2*m)/L^(2*m)) m=1为高斯包络
        %deltaneff=2*1e-5*exp(-(4*log(2)*(-L/2+i*L/N)^2/(0.8*L)^2)^6);
        %deltaneff=2*1e-5*exp(-(56*(-L/2+i*L/N)^4)/(L)^4);
        %(3)升余弦
        %FWHM=L/3;deltaneff_raiandcos=10*1e-5;deltaneff=deltaneff_raiandcos*1/2*(1+cos(pi*Z(i)/FWHM)); 
        %(4)余弦
        %deltaneff_cos=5*1e-5;deltaneff=deltaneff_cos*(1+cos(400*(-L/2+i*L/N)))/2;
        %(5)Hamming
        %FWHM=L/3;H=0.75;deltaneff_Hamming=0.0005;deltaneff=deltaneff_Hamming*(1+H*cos(2*pi*Z(i)/FWHM))/(1+H);
        %(6)Blackman
        %FWHM=L/3;B=0.19;deltaneff_Blackman=0.0005;deltaneff=deltaneff_Blackman*(1+(1+B)*cos(2*pi*Z(i)/FWHM)+B*cos(4*pi*Z(i)/FWHM))/(2+B);
        %(7)tanh
        %alpha=2;beta=0.7;deltaneff_tanh=10*1e-5;deltaneff=deltaneff_tanh*(1+tanh(beta*(1-2*(abs(2*Z(i)/L)^alpha))));
        %(8)Sinc
        %FWHM=L/3;A=4;B=2;deltaneff_Sinc=0.0005;deltaneff=deltaneff_Sinc*(Sinc(1/2*(2*Z(i)/FWHM)^B))^2;
        %(9)Cauchy
        %Cauchy=2;FWHM=L/3;deltaneff_Cauchy=10*1e-5;deltaneff=deltaneff_Cauchy*(1-(2*(-L/2+i*L/N)/FWHM)^2)/(1-(2*Cauchy*Z(i)/FWHM)^2);
        %(10)QINZIXIONG
        % if  Z(i)<(L*(1-0.48)/2)
         if  Z(i)<(-(L*0.26))
               %deltaneff=2*1e-5*((1+T1/4-T1*((-L/2+i*L/N)+(1-F)*L/2)/(L*F))^2)/(1+T1/4);
               %deltaneff=(2*1e-5)*(1+0.25*T1-T1*(((-L/2+i*L/N)+(1-F)*L*0.5)/(L*F))^2)/(1+T1*0.25);
               deltaneff=(2*1e-5)*(26-100*((Z(i)+(0.52)*L*0.5)/(L*0.48))^2)/(26);
          elseif abs(Z(i))<=(L*0.26)
               deltaneff=2*1e-5;
          else Z(i)>(L*0.26)
               %deltaneff=(2*1e-5)*(1+0.25*100-100*((Z(i)-(1-0.48)*L*0.5)/(L*0.48))^2)/(1+100*0.25);
               deltaneff=(2*1e-5)*(26-100*((Z(i)-(0.52)*L*0.5)/(L*0.48))^2)/(26);
          end 
        %高斯型包络
        %deltaneff=40*1e-5*exp(-16*((-L/2+Z(i))/L)^2);
        
          
        lamda_D=(1550-C*L/2+C*i*L/N)*1e-9;  %计算每段光栅本地中心波长lamda_B
        sigma=2*pi*neff*(1/lamda(k)-1/lamda_D)+2*pi*deltaneff/lamda(k)+(4*pi*neff)*C*Z(i)/lamda_D^2;
        kac=pi*deltaneff/lamda(k);
        RB=sqrt(kac^2-sigma^2);
        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)]*F;
    end
    R(k)=(abs(-F(3)/F(1)))^2;
    Q(k)=angle(-F(3)/F(1)); 
end
%%%%%%%%%%%%%%%%计算time delay(1)%%%%%%%%%%%%%%%%%%%%%%%%%%%
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.002)^2*1e-18/(2*pi*3e-4)*(Q(i-1)-Q(i))/deltalamda;    % /* 右极限 */
    else
       tao(i)=(lamda1+i*0.002)^2*1e-18/(2*pi*3e-4)*(Q(i-3)-Q(i-2))/deltalamda;  % /* 左极限 */ 
    end
end
tao(1)=tao(4);tao(2)=tao(4);tao(3)=tao(4);

%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%

% tao(1)=Q(1);tao(2)=Q(2);tao(3)=Q(3);%初始化传输矩阵
% for i=4:M
% %   double  h=1.0;
% %   int i=1;
%   
% %   ds1=Q(1),ds2=Q(2);
%   if(abs(Q(i-1)-Q(i))<=1)
%     tao(i)=(lamda1+i*0.002)^2*1e-18/(2*pi*3e-4)*(Q(i-1)-Q(i))/deltalamda;     % /* 右极限 */
%   else
%     tao(i)=(lamda1+i*0.002)^2*1e-18/(2*pi*3e-4)*(Q(i-3)-Q(i-2))/deltalamda;       % /* 左极限 */
%     
% %     i*=2;
% %     h=h/i;
% %     ds1=((*fun)(x)-(*fun)(x-h))/h;  % /* 左极限 */
% %     ds2=((*fun)(x+h)-(*fun)(x))/h;  % /* 右极限 */
% %     return(ds2);
%   end
% end  
% tao(1)=tao(4);tao(2)=tao(4);tao(3)=tao(4);
  
%%%%%%%%%%%%%%%%计算time delay(2)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%产生双坐标图形
%%%%%%%%%反射率图形%%%%%%%%%%%
clf reset,h_ap=axes('position',[0.13,0.13,0.7,0.75]);  %(0.13 0.13)为图形的左下角坐标
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')  %如果没有这一步则坐标上将没有2.5 7.5等
h_linet=line(lamda*1e9,R,'color','b');
set(get(h_ap,'xlabel'),'string','Wavelength/\mum')
set(get(h_ap,'ylabel'),'string','反射率')
%%%%%%%%%相位图形%%%%%%%%%%%%%  !!!!!!严格来说，Q(k)=angle(-F(3)/F(1)); 才是相位，而tao是光栅时延。
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,'ylabel'),'string','相位   /picosecend')
a=floor(min(tao));b=ceil(max(tao));%floor意思是向负方向舍入而ceil刚好相反
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) 


%%%%%%%%%%%%  在给定的范围内将x轴和y周分别分成十格  %%%%%%%%%%%%%%%%%%%%
% clf reset,h_ap=axes('position',[0.13,0.13,0.7,0.75]);  %(0.13 0.13)为图形的左下角坐标
% 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;  %在给定的范围内将x轴和y周分别分成十格
% set(h_ap,'xtick',pxtick,'ytick',pytick,'xgrid','on','ygrid','on') %如果没有这一步则坐标上将没有2.5 7.5等
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 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,'ylabel'),'string','相位')
% a=floor(min(tao));b=ceil(max(tao));%floor意思是向负方向舍入而ceil刚好相反
% set(h_at,'xlim',[lamda1,lamda2],'ylim',[a,b])
% xpm=get(h_at,'xlim');
% txtick=lamda1:((lamda2-lamda1)/nx):lamda2;
% tytick=a:((b-a)/ny):b;
% set(h_at,'xtick',txtick,'ytick',tytick) 

% [haxes,hline1,hline2]=plotyy(lamda*1e9,R,lamda*1e9,tao,'plot','plot');grid on;
% set(get(haxes(1),'Ylabel'),'String','反射率');
% set(get(haxes(2),'Ylabel'),'String','相位');
% set(hline1,'LineStyle','-','color','b')
% set(hline2,'LineStyle','-','color','r')

% axis([1549 1552]);
% xlabel('Wavelength/\mum');
title('在双坐标中的反射谱和相位.')

% axes(haxes(1));
% ylabel('反射率');
% axes(haxes(2));
% ylabel('相位');
% xlabel('lamda(1549--1552)','FontWeight','bold');