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

% 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

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

if type == 2
   moore = 1;
elseif type == 1
   moore = 0;
else
   'Illegal type'
   return;
end

if (nospec > 4)
   'Too many species'
   return;
end 

testm = cc.*cc';
if sum(sum(testm ~= eye(nospec,nospec)))
   'Illegal dominance type'
   return;
end

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

a = zeros(dim,dim);
if prod(ini) == 0
   d3 = round(dim/3);
   d2 = round(dim/2);
   switch nospec
   case 1
      a(d2,d2)=1;
   case 2
      a(d3,d2)=1;
      a(2*d3,d2)=2;
   case 3
      a(d3,d2)=1;
      a(2*d3,d3)=2;
      a(2*d3,2*d3)=3;
   case 4
      a(d3,d3)=1;
      a(2*d3,d3)=2;
      a(2*d3,2*d3)=3;
      a(d3,2*d3)=4;
   otherwise
      return;
   end
   
else
   for i=1:dim
      for j=1:dim
         a1 = rand(nospec,1) < ini;
         if sum(a1) == 1
            a(i,j)=find(a1==1);
         elseif sum(a1) > 1
            x=find(a1==1);
            s=floor(rand*length(x));
            a(i,j)=x(1+s);
         end
      end
   end
end

if iffig
   plotca(a,0,dim,animation);
end

  
n1 = [];
for t=1:nostep
   
   nact=[];
   for s=1:nospec
      nact = [nact sum(sum(a==s)) ./ dim^2];
   end
   n1 = [n1; nact];
   
   %occupation of new cells by local process
   
      anew = zeros(dim,dim);
      for i=1:dim
         for j=1:dim
            suma=zeros(nospec,1);
            for ii=-1:1
               for jj=-1:1
                  indx=i+ii;
                  indy=j+jj;
                  if i+ii>dim
                     indx = 1;
                  end
                  if i+ii<1
                     indx = dim;
                  end
                  if j+jj>dim
                     indy = 1;
                  end
                  if j+jj<1
                     indy = dim;
                  end
                  if (moore == 1 | ii .* jj == 0) & a(indx,indy)>0 
                     suma(a(indx,indy)) = suma(a(indx,indy))+1;
                  end
                  
               end
            end
                       
            if (a(i,j) > 0)  % occupied cell 
               if length(find(suma>0)) == 1 
                  if find(suma>0) ~= a(i,j) 
                     if cc(a(i,j),suma(find(suma>0))) == 1
                        if rand(1,1)<imi(find(suma>0)).*suma(find(suma>0))  
                            anew(i,j)=find(suma>0);
                        end
                     else
                        anew(i,j)=a(i,j);
                     end
                  else
                     anew(i,j)=a(i,j);
                  end
               elseif length(find(suma>0)) > 1 
                  state=(cc(a(i,j),1:nospec).*suma'.*(a(i,j)~=(1:nospec)).*imi');
                  x=find(state>0);
                  if (length(x)>0)
                     if (rand<sum(state))
                        d=cumsum(state)./sum(state);
                        r=rand;
                        s=max(find(r<d & state>0));
                        anew(i,j)=s;
                     else
                        anew(i,j)=a(i,j);
                        
                     end
                     
                  else
                     anew(i,j)=a(i,j);
                  end
                     
               else
                  anew(i,j)=a(i,j);
               end
               
            else            % empty cells
               if length(find(suma>0)) == 1 
                  if rand(1,1)<imi(find(suma>0)).*suma(find(suma>0))  
                     anew(i,j)=find(suma>0);
                  end
                  
               elseif length(find(suma>0)) > 1 
                  state=(suma.*imi);
                  if (rand<sum(state))
                     d=cumsum(state)./sum(state);
                     r=rand;
                     s=max(find(r<d & state>0));
                     anew(i,j)=s;
                  else
                     anew(i,j)=a(i,j);
                        
                  end
                     
                end  
            end
         end
      end
   
   a = anew;
   
   %mortality and occupation of new cells by long-distance dispersal

   for s=1:nospec
      a = a .* (~(rand(dim,dim) < ext(s) & a == s));
      a = a + s.*(rand(dim,dim) < (long(s) .* nact(s)) & a == 0); %sequential effects deemed unimportant for low long
   end
   
   if iffig
      plotca(a,t,dim,animation);
   end
   

end
return;
   

function plotca(a,t,dim,animation)
   %plotting
   colour=['w','b','g','r','m'];
   clf;
   hold on;
   title(t);
   axis([1,dim+1 1 dim+1]);
   for (i=1:dim)
      for (j=1:dim)
         fill([i,i,i+1,i+1],[j,j+1,j+1,j],colour(a(i,j)+1))
      end
   end
   pause(animation);
return;