function [Xtp1, dXt]=BaranyiRoberts(Xt,t,pourcentbruit, ...
		      nomparam, nomparamS, ...
                      paramS)
%BaranyiRoberts Compute BaranyiRoberts model
% INPUT
% Xt: matrix (number of parameters+1 x number of particles)
%     Contains the values of the active parameters. In output, the last line 
%     contains the model value.
% t: current time in the computational time unit.
% pourcentbruit: percentage of noise (0 if none required)
% nomparam: names of the active parameters. 
% nomparamS: names of the fixed parameters.
% paramS: structure which contains the values of the fixed parameters
%    in components named nomparamS{i}, for i=1 to length(nomparamS).
% OUTPUT
% Xt(end, :):  the model value.
% dXt: vector (number of particles)
%     derivatives of the model with respect to time.
%% ============================================================

%% Set active parameters values
lenp=length(nomparam);
for i=1:lenp
    var=Xt(i,:);
    eval([nomparam{i},'=var;'])
    clear var
end
clear lenp

%% Set fixed parameters values
lens=length(nomparamS);
vect=ones(1,size(Xt,2));
for i=1:lens
    var=nomparamS{i};
    eval([var,'=paramS.(var).*vect;'])
    clear var
end
clear vect lens 
Xtp1=Xt;


%% Compute the model value in  Xtp1(end,:)
% Proteger les argu des exponentielles
e1 = mumax.*(t-lambda);
e1( e1 >100)=100;

expmumaxA=1.0 + exp(e1)-exp(-mumax.*lambda);

B=1+(expmumaxA-1).*N0./Nmax;
Xtp1(end,:)=log(N0.*expmumaxA./B);

%% Compute the derivative of the model in dXt
e1= -t.*mumax - mumax.*lambda;
e1( e1 < -100)= -100; % protection des exponentielles
exp1 = exp(e1);
e2= -t.*mumax;
e2 (e2 < -100)= -100; % protection des exponentielles
exp2=exp(e2);
e3=-lambda .* mumax;
e3 (e3 < -100)= -100; % protection des exponentielles
exp3=exp(e3);
A0t= exp2 + exp3 - exp1; 
A0t(A0t <=0) = 1e-10;
dAt= 1+ ((exp1- exp2) ./ A0t);
At= t + (1./mumax .* log(A0t));
e4=mumax .* At;
e4 (e4 > 100)= 100; % protection des exponentielles
exp4= exp(e4);
dBt= (N0./Nmax) .* mumax .* dAt .* exp4;
Bt= 1 + ( ( exp4 -1) ./ (N0./Nmax));
dXt= N0 .* exp4 .* 1./Bt .* (mumax .* dAt - 1./Bt .* dBt);
dXt=log(dXt);


%% Put noise on the model value
if pourcentbruit>0 
    pbruit=pourcentbruit/100;
    Xtp1(end,:)=(1-pbruit).*Xtp1(end,:)+poissrnd(Xtp1(end,:)*pbruit);
    clear pbruit
end


clear mumax lambda N0 Nmax expmumaxA B 


