getStage<-function(object,stage)
{
it<-0;
x<-c();
for(i in 1:length(object))
{
if(object[i]==-1)
{

it=it+1;
}
else if(it==stage)
{

x<-c(x,object[i]);

}

}
return (x)
}








GaltonWatsonProcess<-function(start,lam,stage)
{
cont<-c(start,-1);
for(i in 1:stage)
{
family<-getStage(cont,i-1);
if(length(family)>0)
{
newFamily<-c();
	for(k in 1:length(family))
	{
		if(family[k]>0)
		{
		newFamily<-c(newFamily,rpois(family[k],lam));
		}
	}
	cont<-c(cont,newFamily,-1);


}
}
return (cont)
}

createPath<-function(max,population,lambda)
{

s<-c();
container<-GaltonWatsonProcess(population,lambda,max)
for(stage in 1:max)
{
k=0;
for(i in 1:1)
{

x<-getStage(container,stage)
if(length(x)!=0)
{
k=length(x)+k;
}
}
s<-c(s,k);
}
ax<-1:max
return (s)



}


createBench<-function(max,population,lambda,patches)
{
ax<-1:max
m<-matrix(ncol=max,nrow=patches,0);


for(i in 1:patches)
{
m[i,]=createPath(max,population,lambda);
}

plot(ax,rep(0,max),ylim=c(0,max(m)),col="white",type="l",ylab="Liczba potomków",xlab="Pokolenie",main="Proces Galtona Watsona")
for(i in 1:patches)
{
lines(ax,m[i,],col=rgb(g=runif(1),b=runif(1),r=runif(1)))
}

}


#Populacja wymierająca

createBench(50,20,0.9,20)


#Populacja rozwojowa

createBench(20,20,1.1,20)

#Populacja stabilna

createBench(20,20,1.0,20)