shife_phase_grating.m

用于模拟相移光纤光栅(FBG)的传输特性,包括多相移光栅

%% shife_phase_grating,用矩阵法据求bragg的反射谱
clear all
clc
Speed=3*1e+8;           % 光速 %
Pg=0.53e-006;           % 光栅周期 %
Lg=0.001;             % 光栅长度 %
n1=1.45;             % 纤芯原始折射率 %             
v=0.878;                    % 调制指数（边缘可见度） %

% n(z)=n1+n1*sigma(z)*(1+m*cos(2*pi*z/Pg+Fai(z))) %

DeltaN_dc=3.2*1e-5;                                                          % DeltaN_dc=n1*sigma(z) %
DeltaN_ac=v*DeltaN_dc;  % 这是折射率改变的交流分量，见Erdogan式（28） %    % DeltaN_dc=n1*sigma(z)*m %
Neff=n1;
lamda_b=1550*1e-9;

S_N=50;
% Dz=Lg/S_N;
% FWHM=0.8*Lg;

lamda=[1545:0.01:1555]*1e-9;      % 初始化输出向量 %

%deltalamda=(lamda2-lamda1)/M*1e-9;
M=5;
Delta_Lg=Lg/M;
% for m=1:M
for j=1:1001
Matrix_g=[1,0;0,1];

DeltaN_ac=v*DeltaN_dc;%*exp(-4*log(2)*(z(k)-(Lg/2))^2/(FWHM^2)),引入高斯切趾会造成左右不平衡
Kdc=(2*pi*DeltaN_ac)./lamda(j);
% Kac=pi*DeltaN_ac/lamda(j);
Kac=4*1e+3;
sigma_dot=Kdc+2*pi*Neff.*(1./lamda(j)-1./lamda_b);%啁啾量为零，所以不用对光栅轴向求z=-L/2+i*L/N.
gamma=sqrt(Kac.^2-sigma_dot.^2);
alphaL=gamma.*Lg;


T11=cosh(alphaL)-i*(sigma_dot/gamma)*sinh(alphaL);%%%%bragg光栅的传输矩阵
T12=-i*(Kac/gamma)*sinh(alphaL);
T21=i*(Kac/gamma)*sinh(alphaL);
T22=cosh(alphaL)+i*(sigma_dot/gamma)*sinh(alphaL);
Matrix_g=Matrix_g*[T11,T12;T21,T22];

    Matrix_g1=Matrix_g;
    Matrix_g2=Matrix_g;

    delta_phase=pi;               %            间隔产生的相移量              %              
    pp_1=exp(-i*delta_phase/2);                        %          FBG光栅对应的相位矩阵           %
    pp_2=exp(i*delta_phase/2);                         % 光栅是LPFG时的相位矩阵是该相位矩阵的共轭  %
    Matrix_fb=[pp_1,0;0,pp_2];

    % Matrix_total=Matrix_g1*Matrix_fb*Matrix_g2;   
    % Matrix_total=Matrix_g1;

    Matrix_total=Matrix_g2*Matrix_fb*Matrix_g1;

    TT11=Matrix_total(1);TT12=Matrix_total(3);TT21=Matrix_total(2);TT22=Matrix_total(4);%Matlab中矩阵元素按列从左到右的顺序标识

    Reflect=TT21/TT11;                              
    Reflectivity=(abs(Reflect)).^2;
    Permeance=1-Reflectivity;
    rr(j)=Reflectivity;
    pp(j)=Permeance;                
end

%产生双坐标图形
%%%%%%%%%反射率图形%%%%%%%%%%%
% 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,Permeance,'color','b');
% set(get(h_ap,'xlabel'),'string','Wavelength/\mum')
% set(get(h_ap,'ylabel'),'string','透射率')
plot(lamda,pp);
axis([1545e-9 1555e-9 0 1]);
title('透射谱（相移光栅，等间距插入3个pi/2相移）');
xlabel('wavelength');
ylabel('Transmission');
grid on;

% subplot(2,1,2);
% plot(temp_lamda,Delay);
% axis([1.5400e-6 1.5520e-6 0 20]);
% title('群延时（相移光栅）');
% xlabel('wavelength');
% ylabel('Group Delay(ps)');
% grid on;