function lifespan = GompertzRandomLifespan(p,badluck,treatment,RR,decay)

%returns a random lifespan duration from the Gompertz probability density function with parameters p1,p2 and p3.
%badluck=rand; %between 0 et 1; 0=infinite lifespan;1=no lifespan
%treatment:1=normal lifespan, 1.1=reduced by 10%
%RelativeRisk:1=normal death, 1.1=augmented by 10%

%GompertzRandomLifespan([0 0.004 5],rand)

%Sa0=exp(p(2)/treatment*p(3)*(1-exp(decay/p(3))));
%S0=exp(p(2)*p(3)*(1-exp(decay/p(3))));%pas de traitement jusqu'au début de l'étude de durée de vie

%maintenant on résoud badluck = exp(-treat[a t + b tau(exp(t/tau)-1)])
lifespan = decay + p(3)/treatment * log(1-log(badluck)/(RR*p(2)*p(3)*exp(decay/p(3)))) / (1+ p(1)/p(2));
%lifespan=fminsearch(inline(['abs(exp(' num2str(p(2)*p(3)/treatment) '*(1-exp( (t+' num2str(decay) ')/' num2str(p(3)) ')))/' num2str(Sa0) ' -' num2str(badluck) ')']), 100);

return;

Sa0=exp(p(2)/traitement*p(3)*(1-exp((1+decay)/p(3)))) / exp(p(2)*p(3)*(1-exp((1+decay)/p(3))));
t=1:700;Sa=exp(p(2)/traitement*p(3)*(1-exp((t+decay)/p(3))))/Sa0;plot(Sa,'k');