多种光纤光栅仿真程序

clc clear; np=0.5345e-6; %光栅周期% L=20000*np; %光栅长度2mm% n1=1.45; %纤芯原始折射率% v=1; %边缘可见度% DeltaN=1e-4; Neff=n1; wl_D=2*Neff*np z=L; % z2=L*2/3; F=zeros(2,2); X=[1;0] %初始化矩阵 a=pi/2; T=[exp(-i*a),0;0,exp(i*a)]; %相移参数 wavelength=[1540:0.01:1559.99]*1e-9; % 初始化输入波长区间% omega=2*pi*3e8./wavelength; NN=0; for j=1:2000 k=pi*v*DeltaN/wavelength(j); Beta=2*pi*Neff./wavelength(j); delta=2*pi*DeltaN/wavelength(j)+Beta-pi/np; g=sqrt(k^2-delta^2); A=cosh(g*z)-i*delta*sinh(g*z)/g; B=-i*k*sinh(g*z)/g; C=cosh(g*z)+i*delta*sinh(g*z)/g; F=[A,B;-B,C]; % A2=cosh(g*z2)-i*delta*sinh(g*z2)/g; % B2=-i*k*sinh(g*z2)/g; % C2=cosh(g*z2)+i*delta*sinh(g*z2)/g; % F2=[A2,B2;-B2,C2]; Y=F*X; % n(j)= angle(Y(2,1)/Y(1,1)); % 相位角theta在(-pi,pi)上的映射 % if j>1 % Dth(j)=n(j)-n(j-1); % if NN==1; % Dth(j-1)=(Dth(j)+Dth(j-2))/2; % NN=0; % end % if abs(n(j)-n(j-1))>pi/2 % NN=1; % end % end ta=Y(2,1)/Y(1,1); n(j)=atan(imag(ta)/real(ta)); if j>1 if n(j)>n(j-1) NN=NN+1; end end theta(j)=n(j)-NN*pi; % 相位角theta为连续量变化 end Dt=diff(theta)./diff(omega)*1e+12; plot(wavelength,n,'.',wavelength(2:2000),Dt); axis([1.5496e-6 1.5509e-6 -5,100]); title('反射率'); xlabel('wavelength'); ylabel('Reflectivity'); grid on;