function[Xtp1, dXt]=DoubleWeibull(Xt,t,pourcentbruit,  ...
		      nomparam, nomparamS, paramS)
%DoubleWeibull Compute DoubleWeibull model
% INPUT/OUTPUT
% See BaranyiRoberts.m
%% ========================================================================

% alpha etant une fonction Matlab, on rajoute 'v' devant le nom des
% parametres pour pas de confusion

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

Xtp1=Xt;
aux1= -(t./vdelta1).^ vp + valpha;
aux2=-(t./vdelta2).^vp;
aux3 = 10.^aux1;
aux4= 10.^aux2;
Xtp1(end,:)=vlog10N0 - log10(1+10 .^ valpha) + log10(aux3 + aux4);

%% Calcul de derivee de Nt par rapport au temps dans dXt
N1= aux3 .* (vp ./vdelta1) .* (t ./vdelta1) .^ (vp-1);
N2= aux4 .* (vp ./vdelta2) .* (t ./vdelta2) .^ (vp-1);
if t>0
  dXt = - (N1+N2) ./ (aux3+aux4);
else
  dXt = zeros(1, size(Xt,2));
end



%% Bruiter Nt
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;

