sfun_fbg.m

fiber bragg grating file

function[sys,x0,str,ts]=sfun_fbg(t,x,u,flag,length,neff,lamda_wg,lamda_cr,N,kapa)


% SFUN_FBG simulates fiber bragg gratings.

switch flag
    %%%%%%%%%%%%%%%%
    % Initialization
    %%%%%%%%%%%%%%%%
    case 0,
        [sys,x0,str,ts]= mdlInitializeSizes(length,neff,lamda_wg,lamda_cr,N,kapa);
        
    %%%%%%%%
    %output
    %%%%%%%%
    case 3,
        sys=mdlOutputs(length,neff,lamda_wg,lamda_cr,N,kapa);
        
    %%%%%%%%%%%%%
    %TErminate
    %%%%%%%%%%%%
    case {1 2 4 9},
        sys=[];     % do  nothing
        
    %%%%%%%%%%%%%%%%%%
    %unexpected Flags
    %%%%%%%%%%%%%%%%%%
    otherwise
        error(['unhandled flag=',num2str(flag)]);
end

%end dsfunc
%mdlInitialSizes
% Return the sizes,initial conditions and sample times for the S-function.

function[sys,x0,str,ts]=mdlInitializeSizes(length,neff,lamda_wg,lamda_cr,N,kapa)

sizes=simsizes;
sizes.NumContStates = 0;
sizes.NumDiscStates = 0;
sizes.NumOutputs    = 4;
sizes.NumInputs     = 2;
sizes.DirFeedthrough= 1;
sizes.NumSampleTimes= 1;

sys = simsizes(sizes);

x0 = ones(sizes.NumDiscStates,1);
str=[];
sample_time=length/30*neff/N;
ts = [sample_time 0];

set_param(gcb,'UserData',[zeros(6,N+1)]);
%Initialization values

temp_f=[];
temp_b=[];
Y_f_estimate=[];
Y_b_estimate=[];
tempvalue=[];

%end mdlInitializeSizes
%mdlOutputs
% Return the output vector for the S-function

function sys=mdlOutputs(t,x,u,length,neff,lamda_wg,lamda_cr,N,kapa)
sys=[x];
% coefficients

a21=5/24;
a22=1/3;
a23=-1/24;
a31=1/6;
a32=2/3;
a33=1/6;

b21=3/8;
b22=0;
b23=9/8;
b31=4/3;
b32=-8/3;
b33=10/3;

k0=2*pi/lamda_wg*neff;      %propagation constant
L=k0*length;                %nominal length

N=100;

h=L/N;                      %step size

% lamda_cr=1550;
% kapa=1.5;
%lamda_cr=lamda_cr;
delta=lamda_wg/lamda_cr-1;  %detuning parameter

%input values

tempvalue=get_param(gcb,'UserData');

temp_f=tempvalue(1,:);
temp_b=tempvalue(2,:);

Y_f_estimate(1:2,:)=tempvalue(3:4,:);
Y_f_estimate(1:2,:)=tempvalue(5:6,:);
Y_b_estimate(1:2,:)=tempvalue(3:4,:);
Y_b_estimate(1:2,:)=tempvalue(5:6,:);


%forward signal
Y_f(1,1:N+1)=temp_f(1:N+1);
Y_f(2,1)=Y_f(1,1)+0.5*h*(i)*(delta*temp_f(1)+kapa*temp_b(1));

%by Euler's method
Y_f(2,2:N)=Y_f_estimate(1,1:N-1);       %obtain by extrapolation

%backward signal
Y_b(1,1:N)=temp_b(1:N);
Y_b(1,N+1)=0;

Y_b(2,1:N-1)=Y_b_estimate(1,1:N-1);     %obtain by extrapolation
Y_b(2,N)=Y_b(1,N+1)+0.5*h*(i)*kapa*temp_f(N+1);  %obtain by Euler;s method

%forward signal
Y_f(3,1)=real(u)+i*imag(u);       %boundary condition, u(1):real part; u(2):imag part.
Y_f(3,2)=Y_f(1,1)+h*(i)*(delta*temp_f(1)+kapa*temp_b(1));   %by Euler's method

Y_f(3,3:N+1)=Y_f_estimate(2,1:N-1); % by extrapolation

%backward signal
Y_b(3,1:N-1)=Y_b_estimate(2,1:N-1); %obtain by extrapolation
Y_b(3,N)=Y_b(1,N+1)+h*(i)*kapa*temp_f(N+1); % by euler's method
Y_b(3,N+1)=0;

function y=g(x1,delta,x2,kapa)
y=i*(delta*x1+kapa*x2);

%calculation
NI=3;
for ii=1:NI
    G1_f=g(Y_b(1,1:N+1));
    G2_f=g(Y_b(2,1:N));
    G3_f=g(Y_b(3,1:N+1));
    
    G1_b=g(Y_f(1,1:N+1));
    G2_b=g(Y_f(2,1:N));
    G3_b=g(Y_f(3,1:N+1));
    
    Y_f(2,1:N)=Y_f(1,1:N)+h*(a21*G1_f(1,1:N)+a22*G2_f(1,1:N)+a23*G3_f(1,2:N+1));
    Y_f(3,2:N+1)=Y_f(1,1:N)+h*(a31*G1_f(1,1:N)+a32*G2_f(1,1:N)+a33*G3_f(1,2:N+1));
    
    Y_b(2,1:N)=Y_b(1,2:N+1)+h*(a21*G1_b(1,2:N+1)+a22*G2_b(1,1:N)+a23*G3_b(1,1:N));
    Y_b(3,1:N)=Y_b(1,2:N+1)+h*(a31*G1_b(1,2:N+1)+a32*G2_b(1,1:N)+a33*G3_b(1,1:N));
end

%calculation of extrapolating values

Y_f_estimate(1,1:N-1)=Y_f(1,1:N-1)+h*g((b21*Y_f(1,1:N-1)+b22*Y_f(2,1:N-1)+b23*Y_f(3,2:N)),delta,(b21*Y_b(1,1:N-1)+b22*Y_b(2,1:N-1)+b23*Y_b(3,2:N)),kapa);
Y_f_estimate(2,1:N-1)=Y_f(1,1:N-1)+h*g((b31*Y_f(1,1:N-1)+b32*Y_f(2,1:N-1)+b33*Y_f(3,2:N)),delta,(b31*Y_b(1,1:N-1)+b32*Y_b(2,1:N-1)+b33*Y_b(3,2:N)),kapa);

Y_b_estimate(1,1:N-1)=Y_b(1,3:N+1)+h*g((b21*Y_b(1,3:N+1)+b22*Y_b(2,2:N)+b23*Y_b(3,2:N)),delta,(b21*Y_f(1,3:N+1)+b22*Y_f(2,2:N)+b23*Y_f(3,2:N)),kapa);
Y_b_estimate(2,1:N-1)=Y_b(1,3:N+1)+h*g((b31*Y_b(1,3:N+1)+b32*Y_b(2,2:N)+b33*Y_b(3,2:N)),delta,(b31*Y_f(1,3:N+1)+b32*Y_f(2,2:N)+b33*Y_f(3,2:N)),kapa);

%output the values

set_param(gcb,'UserData',[Y_f(3,1:N+1);Y_b(3,1:N+1);Y_f_estimate(1,1:N-1);Y_f_estimate(2,1:N-1);Y_b_estimate(1,1:N-1);Y_b_estimate(2,1:N-1)]);

sys=[real(Y_b(3,1)) imag(Y_b(3,1)) real(Y_f(3,N+1)) imag(Y_f(3,N+1))];

%endmdlupdate