metapop2_plot = function(imi1,ext1,ini1,imi2,ext2,ini2,domcon,dim, nostep, animation) {
#This function plots simulated metapopulation dynamics and isolines together

# 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

  a=metapop2(imi1,ext1,ini1,imi2,ext2,ini2,domcon,dim, 0, 0, nostep)



  nobin = 400

  incr1 = 1.1*max(c(1-ext1/imi1,1-ext2/imi2,(imi2-ext2)/(imi1+imi2)))/nobin
  x = NULL
  y1 = NULL
  y2 = NULL
  for (t1 in 0:nobin) {
    x1 = incr1*t1
    if (domcon)
      y1t = (imi2-ext2-(imi1+imi2)*x1)/imi2
    else {
      y1t = (1-x1-ext1/imi1)
      y2t = (1-x1-ext2/imi2)
    }
    x = c(x, x1)
    y1=c(y1, y1t)
    if (!domcon)
      y2=c(y2, y2t)       
  }

  if (!domcon)
    y=cbind(y1, y2)
  else
    y=y1

  if (!domcon)
    matplot(x,y, xlim=c(0, max(a[,1],x)), ylim = c(0, max(a[,2],y)),
            col=c("blue","green"),lty=1,lwd=2,type="l",
            xlab="Abundance of species 1",ylab="Abundance of species 2")   
    #matplot(x,y, col=c("blue","green"),lty=1,lwd=2,type="l")   
  else {
   plot(x,y, xlim=c(0, max(a[,1],x)), ylim = c(0, max(a[,2],y)),
        col="blue",lty=1,lwd=2,type="l",xlab="Abundance of species 1",ylab="Abundance of species 2")
   #plot(x,y, col="blue",lty=1,lwd=2,type="l")
   lines(c(1-ext1/imi1, 1-ext1/imi1),c(0, max(x,y)),col="green",lty=1,lwd=2)
  }
 
  if (animation > 0) {
    for (i in 1:(nostep-1)) {     
      title(main=paste(i),col.main="white")
      lines(c(a[i,1],a[i+1,1]),c(a[i,2],a[i+1,2]),col="red",lty=1,lwd=2,type="l")
      title(main=paste(i+1))
      Sys.sleep(0.1)
    }
  }
  else 
    lines(a[,1],a[,2],col="red",lty=1,lwd=2,type="l")
  
}