%  本程序利用傳輸矩陣法計算啁啾光纖光柵光譜特性
%  采用了科西函數進行切趾（變跡）；
%  定義常數及變量的物理意義
%   gating_length             % 光纖光柵的長度
%   grating_period            % 光纖光柵的周期
%   index_0                   % 光纖的折射率
%   index_modu                % 光纖光柵的折射率調制
%   index_modu_am             % 光纖光柵的折射率調制幅度
%   index_aver                % 光纖光柵的折射率平均增加值
%   index_eff                   光纖光柵的有效折射率
%   wavelength_Bragg          % 光纖光柵Bragg波長
%   v                         % 光纖光柵的折射率調制可見度
%   wave_l1                   % 起始光譜波長
%   wave_l2                   % 終止光譜波長
%   wave_incr                 % 光譜波長掃描的步長
%   wavelength                % 光譜波長
%   detune                    % 失諧Δβ
%   couple_coe                % 模間耦合系數
%   wave_num                  % 光纖光柵的波數失諧
%   phase                     % 相位
%   time_delay                % 時延
%   dispersion                % 色散
%   disp_index                % 色散系數
%  主程序
%  初始賦值過程
clear all;
grating_length=20e6;              % 單位納米nm     
length_incr=0.2e6;
%grating_period=0.5e-3;           
index_0=1.45;                 
index_modu_am=0.25e-4;
%index_modu=index_modu_am
wavelength_center=1400.0;         
grating_period_0=0.5*wavelength_center/index_0;
index_aver=0.0;
index_eff=index_0+index_aver;
%wavelength_Bragg=1.550e-3;         
wave_l1=1399.9;                   
wave_l2=1400.1;                   
wave_incr=0.0005;
v=1.0; 
F=0;
C=F*grating_period_0/(pi*grating_length.^2);
C1=0;
optic_vel=3.0e5;
pai_2=pi+pi;
para=1/(optic_vel*pai_2);
n=round((wave_l2-wave_l1)/wave_incr)+1;
m=grating_length/length_incr+1;
wave_length=[wave_l1:wave_incr:wave_l2];
distance=[0.0:length_incr:grating_length];
%velocity_group=optic_vel/
wavelength_Bragg=2*index_eff*grating_period_0*(1-(distance-0.5*grating_length)*C);
%   光纖光柵折射率調制的變跡函數
% y=(1-4*(distance-0.5*grating_length).^2/grating_length.^2)./(1-4*C1.^2*(distance-0.5*grating_length).^2/grating_length.^2);
index_modu=index_modu_am;
couple_coe=pi*v*index_modu./wavelength_Bragg;
%  波長掃描
 for i=1:n;
 detune= 2.0*index_eff*pi*(1.0./wavelength_Bragg-1.0/wave_length(i));      
 wave_num=sqrt(couple_coe.*couple_coe-detune.*detune)+eps;
 s_distance=wave_num*length_incr;
 s_inh=sinh(s_distance);
 tr_mat_0=[1,0;0,1];
 for k=1:m;
  tr_mat11=cosh(s_distance(k))+(detune(k)./wave_num(k)).*s_inh(k)*j;
  tr_mat12=couple_coe(k)./wave_num(k).*s_inh(k)*j;
  tr_mat21=-tr_mat12;
  tr_mat22=conj(tr_mat11);
  tr_mat_0=tr_mat_0*[tr_mat11,tr_mat12;tr_mat21,tr_mat22];
 end
  tr_mat11=1/tr_mat_0(1,1);
  tr_mat_tr(i)=tr_mat11;
  tr_mat_re(i)=-tr_mat_0(2,1)*tr_mat11;
  refl_index(i)=abs(tr_mat_re(i))^2;
  tran_index(i)=abs(tr_mat11)^2;
  end;
  phase_tr=phase(tr_mat_tr);
  phase_re=phase(tr_mat_re);
  para_2=-wave_length.^2*para;
  [time_delay_tr]=para_2.*gradient(phase_tr,wave_incr);
  [time_delay_re]=para_2.*gradient(phase_re,wave_incr);
  [dispersion_tr]=gradient(time_delay_tr,wave_incr);
  [dispersion_re]=gradient(time_delay_re,wave_incr);
  disp_index_tr=para_2.*dispersion_tr;
  disp_index_re=para_2.*dispersion_re;
  time_delay_re_0=time_delay_re((n+1)/2);
  refl_index_max=max(refl_index);
%顯示光譜
 figure;
plot(wave_length,refl_index),title('反射譜'); 
figure;
plot(wave_length,time_delay_re),axis tight,title('反射時延');
figure;
plot(wave_length,phase_re),axis tight,title('相位');