均匀fbg的反射谱.m

均匀布拉格光栅的仿真程序

% 本程序为：利用 Runger-Kutta 法求解均匀FBG的反射谱
% 单位统一为 纳米
%tongqing write in 2007.4-5
hold on
clear
tic  %估计运算时间
clf
neff=1.4508; % 光纤有效折射率
N=1878; %光纤光栅被分割成N个周期
L=1.002852e6; % 光栅总长度，单位统一为纳米，L=N*Lambda

fid1=fopen('lambdaCMT.txt','wt');
fid2=fopen('rhoCMT.txt','wt');

f1=inline('-i.*delta.*y1-k.*y2','t','lambda','wk','delta','k','y1','y2');    
f2=inline('-k.*y1+i.*delta.*y2','t','lambda','wk','delta','k','y1','y2');       

lambda=1547; %初始波长
while lambda<1552
t0=-501426;    % 从 -L/2 积分到 L/2      
h=2.0057e+003;   % 单位步长，这决定了本方法的精确程度
n=500;
y1=1;    %初始条件  
y2=0;   

wk=534/2; %每一段薄膜几何宽度
delta=2*pi*neff./lambda-pi/(2.*wk);
k=(1.4516-1.45)/(2*wk*1.4516);   % 文章中的式（13）

for j=0:n
    t=t0+j*h;
    
    k11=feval(f1,t,lambda,wk,delta,k,y1,y2);
    k12=feval(f2,t,lambda,wk,delta,k,y1,y2);
    
    k21=feval(f1,t+h/2,lambda,wk,delta,k,y1+h/2*k11,y2+h/2*k12);
    k22=feval(f2,t+h/2,lambda,wk,delta,k,y1+h/2*k11,y2+h/2*k12); 
    
    k31=feval(f1,t+h/2,lambda,wk,delta,k,y1+h/2*k21,y2+h/2*k22);
    k32=feval(f2,t+h/2,lambda,wk,delta,k,y1+h/2*k21,y2+h/2*k22);
    
    k41=feval(f1,t+h,lambda,wk,delta,k,y1+h*k31,y2+h*k32);
    k42=feval(f2,t+h,lambda,wk,delta,k,y1+h*k31,y2+h*k32);
    
    y1=y1+h*(k11/6+k21/3+k31/3+k41/6);    
    y2=y2+h*(k12/6+k22/3+k32/3+k42/6);     % 以上为R-K法的公式运算
end
    rho=(abs(y2./y1)).^2;   
    plot(lambda,rho,'r+');  
    hold on; 
    grid on;
   
    fprintf(fid1,'%1.06f\n',lambda);
    fprintf(fid2,'%1.06f\n',rho);
    lambda=lambda+0.05;
end

time=toc %显示运算时间