source_spectrum_temperaturechange.m

利用传输矩阵法计算光纤啁啾brag光纤光栅滤波器的温度滤波特性

%###############################################################
% In the following program, environmental temperature has changed,
% so, both of the original light source and the reflection spectrum
% of CFBG is also changed. Here, we obtain the final spectrum when
% considering the two foregoing changes.
%###############################################################

% hold on
close all
T0=300;          % initial temperature
delta_T=60;      % temperature variation
T=T0+delta_T    
alpha=-0.1e-6;    % CFBG(Ge-doped,SiO2)temperature sensitivity coefficent
delta_L=alpha*L*delta_T;
L=L+delta_L;
l=L/M;
kesi=8.9e-6;     % CFBG(Ge-doped,SiO2)thermo-optic coefficient 
Neff=Neff*(1+kesi*delta_T);
delta_N=delta_N*(1+kesi*delta_T);

% load sfs_all.mat;  
index=delta_T/5+1;
 
% for i=1:M
%    r(i)=obj(i)/(sfs_all(i,index)/(0.1250e-3));
%     rr=r(i);
%     lamda_B=lamda(i);
%     f=sym('tanh(pi*x/((1+x/Neff)*lamda_B)*l)^2=rr'); 
%     inter=solve(f);   
%     b=subs(inter);
%     delta_N(i)=b(1)/1.5;
%     r_object(i)=tanh(pi*delta_N(i)/((1+delta_N(i)/Neff)*lamda_B)*l)^2;
% end

count=0;
for i=1:M  % using transmit matrix method to obtain the reflecting spectrum                                                                                                                                                                                                                                                                                                                                                                                                                                                                                            
       F=[1,0;0,1];                                                                                                                                                                                                                                                                                                                                                                                                                                                                                               
   for k=1:M                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         
	    Period(k)=Period(k)*(1+alpha*delta_T);
        Beita(i,k)=Beita(i,k)*(1+kesi*delta_T);
        if k<2
           sigma=(Beita(i,k)-pi/Period(k))+2*pi*delta_N(k)/lamda(i); %(formula 17 from Erdogan's 'Fiber Grating Spectra' )
        else
           sigma=(Beita(i,k)-pi/Period(k))+2*pi*delta_N(k)/lamda(i)+8*pi*Neff^2*(l/2)/lamda(k)^2*(Period(k)-Period(k-1))/l;   
        end
        K=pi*v(k)*delta_N(k)/lamda(i);
        sub=sqrt(K^2-sigma^2);       
        F=F*[cosh(sub*L/M)-j*(sigma/sub)*sinh(sub*L/M),-j*(K/sub)*sinh(sub*L/M);j*(K/sub)*sinh(sub*L/M),cosh(sub*L/M)+j*(sigma/sub)*sinh(sub*L/M)];  %传输矩阵
        count=count+1
    end                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     
        r_design(i)=(abs(F(3)/F(1)))^2;                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  
end                                                                                                                                                                                                                                                                                                                                                                                                                                           
    
for i=1:M
final(i)=(sfs_all(i,index)/(0.1250e-3))*r_design(i);
end

b1=b/0.1250e-3;
plot(lamda,b1);
hold on
sfs_all1(:,13)=sfs_all(:,13)/0.1250e-3;
plot(lamda,sfs_all1(:,13),'m')

plot(lamda,final,'r');
plot(lamda,initial_final)