function [survival,maxLifespan,Lifespans] = GompertzRandomSurvival(p,n,treatment,RR,startportion,duration)
%portion=1-selection bias. Ex: 0.9354
%returns a random lifespan duration from the Gompertz probability density function with parameters p1,p2 and p3.
%n=number of people
%treatment:1=normal lifespan, 1.1=accelerated by 10%

if n~=round(n)
    %figure
    hold on;MonStyle='-';%o
    n=round(n);
    %p=[0.1 0.005 5];
    [S,L]=GompertzRandomSurvival(p,n,treatment,RR,startportion,duration);for i=1:L;line([i-1 i i],[S(i) S(i) S(i+1)],'LineWidth',1,'Color',[0;1;1],'LineStyle',MonStyle);end
    [S,L]=GompertzRandomSurvival(p,n,treatment,RR,startportion,duration);for i=1:L;line([i-1 i i],[S(i) S(i) S(i+1)],'LineWidth',1,'Color',[1;0;1],'LineStyle',MonStyle);end
    [S,L]=GompertzRandomSurvival(p,n,treatment,RR,startportion,duration);for i=1:L;line([i-1 i i],[S(i) S(i) S(i+1)],'LineWidth',1,'Color',[1;1;0],'LineStyle',MonStyle);end
    [S,L]=GompertzRandomSurvival(p,n,treatment,RR,startportion,duration);for i=1:L;line([i-1 i i],[S(i) S(i) S(i+1)],'LineWidth',1,'Color',[1;0;0],'LineStyle',MonStyle);end
    [S,L]=GompertzRandomSurvival(p,n,treatment,RR,startportion,duration);for i=1:L;line([i-1 i i],[S(i) S(i) S(i+1)],'LineWidth',1,'Color',[0;1;0],'LineStyle',MonStyle);end
    [S,L]=GompertzRandomSurvival(p,n,treatment,RR,startportion,duration);for i=1:L;line([i-1 i i],[S(i) S(i) S(i+1)],'LineWidth',1,'Color',[0;0;1],'LineStyle',MonStyle);end
    %title(['p=[' num2str(p) '], n=', num2str(n)]);
    survival=S;maxLifespan=L;Lifespans=0;
    return;
end

%GompertzRandomLifespan([0 0.004 5],rand)
decay=GompertzRandomLifespan(p,startportion,1,1,0);

Lifespans=zeros(n,1);
for i=1:n
    Lifespans(i)=min(decay+duration,GompertzRandomLifespan(p,rand,treatment,RR,decay));%0.9354
end

maxLifespan=round(max(Lifespans)+1);
death=hist(round(Lifespans),1:maxLifespan);

survival = 100*ones(1,maxLifespan-1);
Nsurvivors=n;
for i=1:maxLifespan
    if Nsurvivors == 0
        survival(i+1)=0;
    else
        %survival(i+1)=survival(i)*(1-death(i)/Nsurvivors);
        Nsurvivors=Nsurvivors-death(i);
        survival(i)=Nsurvivors;
    end
end

return;