ca = function(nospec, specpara, cc,dim,type, iffig, nostep, animation,inifile) {

# This function simulates spatially explicit community dynamics in a grid

# nospec: number of species (nospec x 4 matrix of data)

# specpara: parameters of individual species
#    imi - immigration rate
#    ext - extinction rate
#    long - long distance dispersal
#    ini - proportion of habitats initially occupied  (0 - one individual in the centre)

# cc: nospec x nospec matrix of interactions; used to define dominance control. Full founder control has all
# off-diagonal elements zero

# dim - number of rows of cells (dim^2 is the number of habitats)
# type - 1: local, von Neumann, 2: local Moore
# iffig - whether to draw a spatial figure
# nostep - number of simulation steps
# animation - animation step in seconds

  imi =specpara[,1];
  ext=specpara[,2];
  long=specpara[,3];
  ini=specpara[,4];
  for (i in 1:nospec)
    cc[i,i] = 1

  if (iffig) {

    a<-system("wmic desktopmonitor get screenwidth",intern=T) 
    scrw = as.numeric(a[2])
    #scrw=700
    graphics.off()

  }
moore = type - 1

if (nospec > 4) {
   print('Too many species')
   return();
} 

testm = cc*t(cc);
if (sum(testm) > nospec) {
   print('Conflicting off-diagonal elements in the interaction matrix')
   return();
} 

if (sum(long > 0.05)) {
   print('Long distance dispersal too high; sequential effects may be important')
   return()
}

a = matrix(rep(0,dim^2),nrow=dim);

if (!file.exists(as.character(inifile))) {
  if (prod(ini) == 0) {
    d3 = as.integer(dim/3)+1;
    d2 = as.integer(dim/2)+1;
    if (nospec == 1) {
      a[d2,d2]=1
    }
    if (nospec == 2) {
      a[d3,d2]=1;
      a[2*d3,d2]=2;
    }
    if (nospec == 3) {
      a[d3,d2]=1;
      a[2*d3,d3]=2;
      a[2*d3,2*d3]=3;
    }
    if (nospec == 4) {
      a[d3,d3]=1;
      a[2*d3,d3]=2;
      a[2*d3,2*d3]=3;
      a[d3,2*d3]=4;
    }
  }
   
  else {
    for (i in 1:dim) {
      for (j in 1:dim) {
         a1 = runif(nospec) < ini;
         if (sum(a1) == 1)
            a[i,j]=which(a1==1)
         if (sum(a1) > 1) {
            x=which(a1==1)
            s=as.integer(runif(1,max=length(x)))
            a[i,j]=x[1+s];
         }
      }
    }
  }
}
else {
  init = read.table("ini.txt",header=T)
  ll = length(init[,1])
  for (i in 1:ll) 
    a[init[i,2]:(init[i,2]+init[i,4]),init[i,3]:(init[i,3]+init[i,5])] = init[i,1]
}

if (iffig) {
  dev.new(2,xpos=10,ypos=10,width=5,height=5)
  plot(NULL,xlim=c(0,nostep+1),ylim=c(0,1), xlab='Time step',ylab='Proportion of cells occupied')
  dev.new(3,xpos=10+scrw/2,ypos=10,width=5,height=5)
  plot.new()
}

clr=c("blue","green","red","yellow")

if (iffig)
   plotca(a,0,animation,scrw);



  
n1 = NULL;
for (t1 in 1:nostep) {
   
   nact=NULL;
   for (s in 1:nospec) { 
      nact = c(nact, sum(a==s) / dim^2);
   }
   n1 = rbind(n1, nact)
   if (t1>1 & iffig) {
     dev.set(2)
     for (i in 1:nospec) {
       lines(c(t1,t1+1),c(n1[t1-1,i],n1[t1,i]),col=clr[i],lwd=2)
     }
   }
   else
     print(t1)

   asum = array(0,dim=c(dim,dim,8))
   
   #occupation of new cells by local process
      asum[,,1] = a[,c(dim,1:(dim-1))]
      asum[,,2] = a[,c(2:dim,1)]
      asum[,,3] = a[c(dim,1:(dim-1)),]
      asum[,,4] = a[c(2:dim,1),]
      if (moore) {
        asum[,,5] = a[c(dim,1:(dim-1)),c(dim,1:(dim-1))]
        asum[,,6] = a[c(dim,1:(dim-1)),c(2:dim,1)]
        asum[,,7] = a[c(2:dim,1),c(dim,1:(dim-1))]
        asum[,,8] = a[c(2:dim,1),c(2:dim,1)]
      }
      anew = matrix(rep(0,dim^2),nrow=dim)
      for (i in 1:dim) {
         for (j in  1:dim) {
            suma=rep(0,nospec);
            if (moore) 
              tbl = table(asum[i,j,])
            else
              tbl = table(asum[i,j,1:4])
            sp=as.numeric(names(tbl))
            qu = as.numeric(tbl)
            ind = as.numeric(sp[1]==0) #whether there is open space in the neighbourhood
            for (ij in (1+ind):length(tbl)) {
              suma[sp[ij]] = qu[ij]
            }
            nstemp = length(tbl) - ind
            if (a[i,j] > 0) { # occupied cell 
               if (nstemp == 1 ) { 
                  if (sp[ind+1] != a[i,j]) {  
                     if (cc[a[i,j],sp[ind+1]] == 1 ) {
                        if (runif(1) < (1 - (1 - imi[sp[ind+1]])^qu[ind+1]) ) 
                            anew[i,j]=sp[ind+1]
                     }
                     else {
                        anew[i,j]=a[i,j]
                     }
                  }
                  else {
                     anew[i,j]=a[i,j]
                  }
               }
               if (nstemp > 1) { 
                  state=(cc[a[i,j],1:nospec] * (a[i,j] != (1:nospec)) * (1 - (1-imi)^suma)); 
                  x=which(state>0)
                  if (length(x)>0) {
                     if (runif(1)<1-prod(1-state)) {
                        d=cumsum(state) / sum(state)
                        r=runif(1);
                        s=min(which(r<d & state>0))
                        anew[i,j]=s;
                     }
                     else {
                        anew[i,j]=a[i,j]
                     }   
                     
                  }   
                  else {
                     anew[i,j]=a[i,j]
                  }
               }      
               if (nstemp == 0) {
                  anew[i,j]=a[i,j]
               }
            }
            else {            # empty cells
               if (nstemp == 1 ) {
                  if (runif(1) < 1 - (1 - imi[sp[ind+1]]) ^ qu[ind+1]) {  
                     anew[i,j]=sp[ind+1];
                  }
               }   
               if (nstemp > 1) { 
                  state= 1-(1-imi)^suma ; 
                  if (runif(1)<1-prod(1-state)) {
                     d=cumsum(state)/sum(state); 
                     r=runif(1);
                     s=min(which(r<d & state>0));
                     anew[i,j]=s;
                  }
                  else {
                     anew[i,j]=a[i,j];
                        
                  }
                     
               }  
               if (nstemp == 0) { 
                  anew[i,j]=a[i,j]
               }
            }
         }
      }
   
   a = anew;
   
   #mortality and occupation of new cells by long-distance dispersal

   for (s in 1:nospec) {
      a = a * (!(matrix(runif(dim^2),nrow=dim) < ext[s] & a == s))
      a = a + s * (matrix(runif(dim^2),nrow=dim) < (long[s] * nact[s]) & a == 0); #sequential effects deemed unimportant for low long
   }
   
   if (iffig)
      plotca(a,t1,animation,scrw);
   
   

}
list(abund=n1,map=a)
}
   

plotca = function(a,t1,animation,scrw) {
   #plotting
   dev.set(3)
   #plot.new()
   title(main=paste(t1-1),col.main="white")
   image(a,add=T,col=c("white","blue","green","red","yellow"),zlim=c(0,4))
   title(main=paste(t1))
   axis(1,labels=F)
   axis(2,labels=F)
   axis(3,labels=F)
   axis(4,labels=F)
   if (animation)
      Sys.sleep(0.3)   
}