metapop2 = function(imi1,ext1,ini1,imi2,ext2,ini2,domcon,dim, iffig, animation, nostep) {
# This function calculates proportion of occupied islands for two species
# imi - immigration rate
# ext - extinction rate
# ini - proportion of habitats initially occupied  (0 - one individual in the centre)
# 1 and 2 refer to first (green) and second (blue) species
# domcon - 1: dominance control by the first species; 0: founder control
# dim - number of rows of cells (dim^2 is the number of habitats)
# iffig - whether to draw a spatial figure
# nostep - number of simulation steps

  if (ini1 > 0)
    a1=matrix(as.numeric(runif(dim^2)<ini1),nrow=dim)
  else {
    a1 = matrix(rep(0,dim^2),nrow=dim)
    a1[dim/2+1,dim/2+1]=1;
  }

  if (ini2 > 0)
    a2=matrix(as.numeric(runif(dim^2)<ini2),nrow=dim)
  else {
    a2 = matrix(rep(0,dim^2),nrow=dim)
    a2[dim/2+1,dim/2+1]=1;
  }

  anew1 = a1*(!(a1+a2==2));
  anew2 = a2*(!(a1+a2==2));

  a1=anew1;
  a2=anew2;

  n1 = NULL;
  n2 = NULL;

  for (t1 in 1:nostep) {
    nact1 = sum(a1>0)/ dim^2;
    n1 = c(n1, nact1)
    nact2 = sum(a2>0)/ dim^2;
    n2 = c(n2, nact2)
   
    if (domcon==1) {
      anew1 = a1 + (as.numeric(runif(dim^2) < (imi1 * nact1)))
      anew2 = a2 + (as.numeric(runif(dim^2) < (imi2 * nact2) & anew1 == 0))
      anew2 = anew2 * (!(anew1>0));
    }
    else {
      anew1 = a1 + (as.numeric(runif(dim^2) < (imi1 * nact1) & a2 == 0 & a1 == 0))
      anew2 = a2 + (as.numeric(runif(dim^2) < (imi2 * nact2) & a2 == 0 & a1 == 0))
      ax = as.numeric(runif(dim^2))
      ax1 = anew1+anew2==2;
      anew1 = anew1 * !(ax>0.5 & ax1);
      anew2 = anew2 * !(ax<0.5 & ax1);      
    }
      
    anew1 = anew1 * (as.numeric(runif(dim^2) > ext1))
    anew2 = anew2 * (as.numeric(runif(dim^2) > ext2))

    #a1 = as.numeric(anew1 > 0)
    #a2 = as.numeric(anew2 > 0)
    a1 = anew1 
    a2 = anew2
   
    #plotting
    if (iffig) {
      image(a1*2+a2,col=c("white","blue","green"),main=paste(t1),zlim=c(0,2))
      if (animation)
        Sys.sleep(0.3)
    }
  }
  cbind(n1,n2)
}