function n = metapop2(imi1,ext1,ini1,imi2,ext2,ini2,domcon,dim, iffig, 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 = rand(dim,dim) < ini1;
else
   a1=zeros(dim,dim);
   a1(floor(dim/3),floor(dim/2))=1
end

if ini2 > 0
   a2 = rand(dim,dim) < ini2;
else
   a2=zeros(dim,dim);
   a2(floor(2*dim/3),floor(dim/2))=1
end
anew1 = a1.*(~(a1+a2==2));
anew2 = a2.*(~(a1+a2==2));

a1=anew1;
a2=anew2;

n1 = [];
n2 = [];

for t=1:nostep
   nact1 = sum(sum(a1>0)) ./ dim^2;
   n1 = [n1 nact1];
   nact2 = sum(sum(a2>0)) ./ dim^2;
   n2 = [n2 nact2];
   
   if (domcon==1)
      anew1 = a1 + (rand(dim,dim) < (imi1 .* nact1));
      anew2 = a2 + (rand(dim,dim) < (imi2 .* nact2) & anew1 == 0);
      anew2 = anew2 .* (~(anew1>0));
   else
      anew1 = a1 + (rand(dim,dim) < (imi1 .* nact1) & a2 == 0 & a1 == 0);
      anew2 = a2 + (rand(dim,dim) < (imi2 .* nact2) & a2 == 0 & a1 == 0);
      ax = rand(dim,dim); 
      ax1 = anew1+anew2==2;
      anew1 = anew1 .* ~(ax>0.5 & ax1);
      anew2 = anew2 .* ~(ax<0.5 & ax1);
      
   end
      
   anew1 = anew1 .* (rand(dim,dim) > ext1);
   anew2 = anew2 .* (rand(dim,dim) > ext2);

   a1 = anew1 > 0;
   a2 = anew2 > 0;
   
   %plotting
   colormap([1 1 1;0 0 1; 0 1 0; 1 0 0]);
   if iffig
      clf;
      title(t);
      hold on;
      axis([0.5,dim+0.5 0.5 dim+0.5])
      image(a1+2*a2+1);
      pause(0.3);
   end
end
n = [n1;n2];