Essa função está em construção… ainda precisa de trabalho. Estamos no caminho!

meta.espacial=function(tf,c,l,fi,i,pe, pi)
{
paisag=array(0,dim=c(l,c,tf))
nmanchas=c*l
paisag[,,1]=matrix(sample(c(rep(1,fi*nmanchas), rep(0,round((1-fi)*nmanchas)))), ncol=c)
resultado=numeric()
	for(t in 2:tf)
	{
	vazio=which(paisag[,,(t-1)]==0, arr.ind=TRUE)
	nvazio=dim(vazio)[1]
	ocupa=which(paisag[,,(t-1)]==1, arr.ind=TRUE)
	nocupa=dim(ocupa)[1]
	dif=outer(ocupa,vazio, "-")
	difx=dif[,1,,1]
	dify=dif[,2,,2]
	dist=sqrt(difx^2+dify^2)
	dist.media=apply(dist,2,mean)
	dist.max=sqrt(c^2+l^2)
	dist.rel=dist.media/dist.max
	pi=dist.rel*(nocupa/nmanchas)
	paisag[,,t][paisag[,,(t-1)]==1]<-sample(c(0,1),nocupa,replace=T,prob=c(pe,1-pe))
	paisag[,,t][paisag[,,(t-1)]==0]<-rbinom(nvazio,1,prob=pi)
	resultado[t-1]=sum(paisag[,,t])/nmanchas
	}
F=1-(pe/i)
plot(1:tf,c(fi,resultado),type="l",xlab="Tempo",ylab="Fração de manchas ocupadas",
ylim=c(0,1),main="Dinâmica de ocupação de manchas",font.lab=2,lwd=2)
abline(h=F,col=2,lwd=2,lty=2)
return(paisag)
}

Os argumentos da função são:

• tf ⇒ número de iterações ou tempo final da simulação
• c = número de colunas do seu universo de manchas
• l = número de linhas
• fi= proporção de manchas ocupadas no tempo inicial
• pe= probabilidade de extinção

Um exemplo:

meta.espacial(tf=10,c=10,l=10,fi=0.2,pe=0.4)