function [Curves,Lifespans] = GompertzRandomCurves(Ncurves,NperCurve,p,treatment,portion,duration)


Curves=zeros(Ncurves,duration);
Lifespans=zeros(Ncurves,NperCurve);

'Generating Lifespans'
decay=GompertzRandomLifespan(p,portion,1,0);
for i=1:Ncurves
    disp([num2str(i) '/' num2str(Ncurves)]);
    for j=1:NperCurve
        Lifespans(i,j)=GompertzRandomLifespan(p,rand*portion,treatment,decay);
    end
end

'Building survival curves'
for i=1:Ncurves
    maxLifespan=max(Lifespans(i,:));
    death=hist(round(Lifespans(i,:)),1:max(maxLifespan,duration));
    Curves(i,1)=NperCurve;
    for j=1:min(maxLifespan-1,duration-1)
        Curves(i,j+1)=Curves(i,j)-death(j);
    end
end

return;

figure;hold on;for i=1:100;plot(CurvesA21N50D700P15(i,:));end
t=1:700;Sa=50*exp(-p(1)*(1*t+decay)-p(2)*p(3)*exp((1*t+decay)/p(3))+p(2)*p(3))/portion;plot(Sa,'og');
StdA21N50D700P15=0;for j=1:700;StdA21N50D700P15(j)=std(CurvesA21N50D700P15(:,j));end
t=1:700;Sa=50*exp(-p(1)*(1*t+decay)-p(2)*p(3)*exp((1*t+decay)/p(3))+p(2)*p(3))/portion+2*StdA21N50D700(t);plot(Sa,'g');
t=1:700;Sa=50*exp(-p(1)*(1*t+decay)-p(2)*p(3)*exp((1*t+decay)/p(3))+p(2)*p(3))/portion-2*StdA21N50D700(t);plot(Sa,'g');