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

meta.espacial=function(tf,c,l,fi,i,pe, pi, canto=FALSE)
{
paisag=array(0,dim=c(l,c,tf))
nmanchas=c*l
  if(canto==TRUE)
  {
  paisag[1:round(l*fi),1:round(c*fi),1]=1
  }
  else
  {
  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.max-dist.media)/dist.max
	pi.aj=pi*dist.rel
	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.aj)
	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)
invisible(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
• canto=manchas ocupadas no tempo inicial agrupadas ou não (TRUE/FALSE)

Copie também a função abaixo para fazer a simulação animada.

anima2=function(dados)
{
	tf=dim(dados)[3]
	for(i in 1:tf){
	image(dados[,,i], main=("Ocupação de manchas"),col=c("white","red"),bty="n",xaxt='n',yaxt='n')
	grid(dim(dados)[1],dim(dados)[2])
	Sys.sleep(.2)
	}
	}

Pronto! Agora é só rodar o exemplo abaixo e depois mudar os parâmetros da primeira linha para rodar novas simulações:

simula=meta.espacial(tf=100,c=10,l=10,fi=0.1,pe=0.2,pi=0.3, canto=TRUE)
x11()
anima2(simula1)