function [Xtp1, dXt] =LnBaranyiRoberts(Xt,t,pourcentbruit, ...
		      nomparam, nomparamS, ...
                      paramS)
%LnBaranyiRoberts Compute BaranyiRoberts model in  version neperian log
% INPUT
% See BaranyiRoberts.m
% OUTPUT
% Xtp1(end, :) model output = ln(Nt)
% dXt: vector (number of particles)
%     derivatives of the model with respect to time.
%% ============================================================

%% Affectation des  parametres a estimer 

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

% Affectation des  parametres fixes
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;


%% Calcul de loge(Nt) dans Xtp1(end,:)
bexp=100; % borne admissible des exponentielles
e1= -t.*mumax - mumax.*lambda;
e1( e1 < -bexp)= -bexp; % protection des exponentielles
exp1 = exp(e1);
e2= -t.*mumax;
e2 (e2 < -bexp)= -bexp; % protection des exponentielles
exp2=exp(e2);
e3=-lambda .* mumax;
e3 (e3 < -bexp)= -bexp; % protection des exponentielles
exp3=exp(e3);
A0t= exp2 + exp3 - exp1;
A0t(A0t <=0) = 1e-10; % protection des log
At= t + (1./mumax .* log(A0t));
dAt= 1+ ((exp1- exp2) ./ A0t);
e4=mumax .* At;
e4 (e4 > bexp)= bexp; % protection des exponentielles
exp4= exp(e4);
z=exp(lnNmax - lnN0);
if any(z==0) || any(isinf(z))
  errordlg('Error in the model evaluation');
  return
end

Bt= 1 +  ( exp4 -1) ./ z;
if any(Bt <= 0)
  errordlg('Error in the model evaluation');
  return
end
Xtp1(end,:)=lnN0+  (mumax .* At) - log(Bt);


%% Calcul de derivee de loge(Nt) par rapport au temps dans dXt
dBt= exp(lnN0 - lnNmax) .* mumax .* dAt .* exp4;
dXt= mumax .* dAt - (dBt./Bt);

%% Bruiter loge(Nt)
if pourcentbruit>0 
    pbruit=pourcentbruit/100;
% Si on voulait calculer la reponse bruitée sous la forme exponentielle de celle-ci
%    lnX=exp(Xtp1(end,:));
%    lnX=(1-pbruit).*lnX+poissrnd(lnX*pbruit);
%      Xtp1(end,:)=log(lnX);

 Xtp1(end,:)= normrnd(Xtp1(end, :), pbruit); 
%Xtp1(end,:)+(normrnd(0,1) .* pbruit .* (Xtp1(end,:)/2.0));

    clear pbruit
end


clear mumax lambda LnN0 Nmax expmumaxA  


