%% 相移光柵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                     % 取1000個波長點 %
Matrix_g=[1,0;0,1];              % 初始化傳輸矩陣%
DeltaN_ac=v*DeltaN_dc;          
Kdc=(2*pi*DeltaN_ac)./lamda(j); 
                                        % Kac=pi*DeltaN_ac/lamda(j); 
Kac=4*1e+3;     %  kac為耦合系數  %
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; 
delta_phase=pi;               %            間隔產生的相移量              %  
T11=(cosh(alphaL)-i*(sigma_dot/gamma)*sinh(alphaL))*exp(i*sigma_dot*Lg);  %%%%bragg光柵的傳輸矩陣 
T12=-i*(Kac/gamma)*sinh(alphaL)*exp(i*sigma_dot*Lg)*exp(-i*delta_phase); 
T21=i*(Kac/gamma)*sinh(alphaL)*exp(-i*sigma_dot*Lg)*exp(i*delta_phase); 
T22=cosh(alphaL)+i*(sigma_dot/gamma)*sinh(alphaL)*exp(i*sigma_dot*Lg); 
Matrix_g=Matrix_g*[T11,T12;T21,T22]; 
 
    TT11=Matrix_g(1);TT12=Matrix_g(3);TT21=Matrix_g(2);TT22=Matrix_g(4); %Matlab中矩陣元素按列從左到右的順序標識 
 
    Reflect=TT21/TT11;                               
    Reflectivity=(abs(Reflect)).^2; 
    Permeance=1-Reflectivity; 
    rr(j)=Reflectivity; 
    pp(j)=Permeance;                 
end 
 
plot(lamda,rr); 
axis([1545e-9 1555e-9 0 1]); 
title('反射譜（等間距插入1個相移點的反射譜（pi））'); 
xlabel('wavelength'); 
ylabel('Reflectivity'); 
grid on;