Dear r-help members,

Could any of you help me with this model, please?

This model gives error when some value touch whatever border and I do not know how to correct it. The 80% of the seeds produced by a plant will fell into the parent cell, the 15% in the first ring according to the king movement (in chess), and a 5% in the second ring defined by the queen2 matrix. Someone said me the functions that order the position of the seed dispersed, they are:

  neighb <- aperm(apply(aperm(outer(king,IJ,'+'),c(2,4,1,3)),3:4,diag),c(2,1,3)) #Seed going back into source cell???
  for(i in 1:nrow(IJ)) SB[matrix(neighb[,,i],nc=2)]<- SR1[matrix(IJ[i,],1,2)]/8 + SB[neighb[,,i]] #Seed going to anuli 1 neighbor cells from all possible sources
  neighb <- aperm(apply(aperm(outer(queen2,IJ,'+'),c(2,4,1,3)),3:4,diag),c(2,1,3))
  for(i in 1:nrow(IJ)) SB[matrix(neighb[,,i],nc=2)]<- SR2[matrix(IJ[i,],1,2)]/8 + SB[neighb[,,i]] #Seed going to anuli 2 neighbor cells from all possible sources

But they produce that when a dispersed seed goes into a border, it stops the model and gives an error. I guess I should put a buffer in the border to get the value 0 each time some seed reach that buffer but I do not how to manage the model with those complicated dispersal functions.
I do not understand how these functions disperse the seeds, Do they it in a particular order? The queen2 ring should be 16 cells, why in the 'SR2[]' expression divides by 8 and not by 16?
Do anyone know if the dispersed seed are properly included in the SB for the next iteration?

Here I attach you the full model to be this message understandable:
width = 10
height= 10
endt = 6
nPop = 1
SBLL = 1
SBUL = 10
LrD <- matrix(0, width, height)
x <- as.integer(runif(nPop,min=1,max=10))  #uniform random x ccordinates for occupied cells
y <- as.integer(runif(nPop,min=1,max=10))  #uniform random y ccordinates for occupied cells
z <- cbind(x,y)
LrD[z] = runif(nPop,min=SBLL,max=SBUL)      #uniform random SB densities assigned to random locations

par(mfrow=c(1,1))
image(x=1:width,y=1:height,z=LrD, col=terrain.colors(6))

#Parameter values for population simulation
sl1 =       matrix(0.93, ncol=width, nrow=height)
gl1 =       matrix(0.64, ncol=width, nrow=height)
vB = matrix(0.76, ncol=width, nrow=height)
f = matrix(100, ncol=width, nrow=height)
a = matrix(0.34, ncol=width, nrow=height)

#Initial population values at time = 0
SB <- matrix(LrD, ncol=width, nrow=height)
SDL <- matrix(gl1*LrD, ncol=width, nrow=height)
AP <- matrix(vB*SDL, ncol=width, nrow=height)
SP <- matrix((f*AP)/(1+a*AP), ncol=width, nrow=height)
gN <- matrix(0, ncol=5, nrow=endt)

#Define dispersal anuli at radius 1 (king) and 2 (queen2).
  king<- matrix(c(1,1,1,-1,-1,1,-1,-1,0,1,0,-1,1,0,-1,0),nc=2,byrow=T) #first anuli
  queen2<- as.matrix(rbind(expand.grid(seq(-2,2,1),c(-2,2)),expand.grid(c(-1,0,1),c(-2,2))))  #second anuli

# Main population simulation loop
for (t in 1:endt){
  SB<- round(SB)
  SDL = SB*gl1
  AP = SDL*vB
  SP = (f*AP)/(1+a*AP)
  SR1 = 0.15 * SP    #Total Seed dispersed to radius 1
  SR2 = 0.05 * SP    #Total Seed dispersed to radius 2

  IJ = which(SP>0,arr.ind=T)
  SB[IJ] = 0.8 * SP[IJ]
  neighb <- aperm(apply(aperm(outer(king,IJ,'+'),c(2,4,1,3)),3:4,diag),c(2,1,3)) #Seed going back into source cell???
  for(i in 1:nrow(IJ)) SB[matrix(neighb[,,i],nc=2)]<- SR1[matrix(IJ[i,],1,2)]/8 + SB[neighb[,,i]] #Seed going to anuli 1 neighbor cells from all possible sources
  neighb <- aperm(apply(aperm(outer(queen2,IJ,'+'),c(2,4,1,3)),3:4,diag),c(2,1,3))
  for(i in 1:nrow(IJ)) SB[matrix(neighb[,,i],nc=2)]<- SR2[matrix(IJ[i,],1,2)]/8 + SB[neighb[,,i]] #Seed going to anuli 2 neighbor cells from all possible sources

  image(1:width,1:height, SB, col=terrain.colors(6), xlab="", ylab="")
  gN[t,1:5] <- cbind(t, mean(SB), mean(SDL), mean(AP), mean(SP))
}
gN

Thank you very much in advance, and congratulations for this magnificent help that all of you offer,
Judit Barroso

	[[alternative HTML version deleted]]

