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


% Affectation des  paramètres à estimer 
lenp=length(nomparam);
for i=1:lenp
    var=Xt(i,:);
    eval([nomparam{i},'=var;'])
    clear var
end
clear lenp 
lens=length(nomparamS);
vect=ones(1,size(Xt,2));
for i=1:lens
    var=nomparamS{i};
    eval([var,'=paramS.(var).*vect;'])
end
clear vect lens 

Xtp1=Xt;

aux1=exp(-kmax1 .* t); 
aux2 = exp(kmax1 .* S1); 
grandf= f .* aux1;
aux3 = exp(-kmax2 .* t);
G =(1-f) .* aux3;
H= aux2 ./ (1+ (aux2 -1) .* aux1);
Xtp1(end,:)=log10N0 + log10((grandf +G) .* H);
% derivée
dFt = f .* kmax1 .* aux1;
dGt = - (1-f) .* kmax2 .* aux3;
num= aux2 .* (aux2-1) .* aux1 .* kmax1;
den = 1 + ( ( aux2 -1) .* aux1);
dHt= num ./ (den .* den);
Z= (dFt + dGt) .* H + (grandf + G) .* dHt;
log10e= log10(exp(1));
dXt = (log10e * Z) ./ ((grandf +G) .* H);


%% Bruiter
if pourcentbruit>0
    pbruit=pourcentbruit/100;
    Xtp1(end,:)= normrnd(Xtp1(end, :), pbruit);
    clear pbruit
end


%% Ramener les réponses non valides:
% log10(Nt) doit etre >0
Xtp1(end, Xtp1(end,:) < 0) = 0;

