ecoloop2.m

一个利用元胞自动机进行生态仿真的matlab程序

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  ECOLOOP2
%  Simulation Tool for Ecological Phenomena
%  in 2 dimensional lattice
%
%    First Coded   26   February  1997
%    Last Modified 3    June     1998
%
%       Copyright (c) T. Suzudo
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Whether cells stay alive, die, or generate new cells depends
% upon how many of their eight possible neighbors are alive.
% Here we generate index vectors for four of the eight neighbors.
% We use periodic (torus) boundary conditions at the edges of the universe.
s = [nl 1:nl-1]; % this means [nl 1 2 ... nl-1].
e = [2:nl 1];
n = [2:nl 1];
w = [nl 1:nl-1];

% Loop starts
while (t < tfinal) & (cp(t+1)>0)
	if (t>0) 
		X=Y;
	end

	t=t+1;

% If number of neighbors is 5
	Y = zeros(nl,nl);
	if lnei==5
		if tot(2)==1 %Totalistic rule
			N=X+X(e,:)+X(w,:)+X(:,n)+X(:,s);
			for ii =1:lnei+1 Y = Y + (N==(ii-1))*tr_code225(ii); end
		elseif tot(2)==2 % Outer-totalistic rule
			N=X(e,:)+X(w,:)+X(:,n)+X(:,s);
			for ii=1:lnei Y=Y+(N==(ii-1)&~X)*otr_code225(ii); end
			for ii=lnei+1:2*lnei Y=Y+((N==(ii-lnei-1))&X)*otr_code225(ii); end
		elseif tot(2)==3 % Symmetric rule
			N=X(e,:)+X(w,:)+X(:,n)+X(:,s);
			ND=mod(X(e,:)+X(w,:),2); % If not diagonal, ND=1 
			Y=Y+((N==0)&~X)		*sym_code225(1);
			Y=Y+((N==1)&~X)		*sym_code225(2);
			Y=Y+((N==2)&~X&ND)	*sym_code225(3);
			Y=Y+((N==2)&~X&~ND)	*sym_code225(4);
			Y=Y+((N==3)&~X)		*sym_code225(5);
			Y=Y+((N==4)&~X)		*sym_code225(6);
			Y=Y+((N==0)&X)			*sym_code225(7);
			Y=Y+((N==1)&X)			*sym_code225(8);
			Y=Y+((N==2)&X&ND)		*sym_code225(9);
			Y=Y+((N==2)&X&~ND)	*sym_code225(10);
			Y=Y+((N==3)&X)			*sym_code225(11);
			Y=Y+((N==4)&X)			*sym_code225(12);
		end	
% If number of neighbors is 9
	else
		if tot(2)==1
			N=X+X(e,:)+X(w,:)+X(:,n)+X(:,s)+X(e,s)+X(e,n)+X(w,s)+X(w,n);
			for i =1:lnei+1 Y = Y + (N==(i-1))*tr_code229(i); end
		elseif tot(2)==2
			N=X(e,:)+X(w,:)+X(:,n)+X(:,s)+X(e,s)+X(e,n)+X(w,s)+X(w,n);
			for ii=1:lnei Y=Y+(N==(ii-1)&~X)*otr_code229(ii); end
			for ii=lnei+1:2*lnei Y=Y+((N==(ii-lnei-1))&X)*otr_code229(ii); end
		end	
	end

%%  This is special for time-entropy calculation (Disable when not necessary)
%	if ii==10 Y3(:,:,size(Y3,3)+1)=X; end

% animation
	if (ani(2)>1 & gra)
	      [x,y]=find(Y);
		set(plothandle,'xdata',x,'ydata',y);
   	end

% *** The following two lines shall be disabled for cataloged procedure ***
	set(timehandle,'String',num2str(t));
	drawnow; 
%	if ani(2)==1 drawnow; end
	cp(t+1,1)=sum(sum(Y));
end
% graphics
if (ani(2)==1 & gra)
	[x y]=find(Y);
	set(plothandle,'xdata',x,'ydata',y);
	drawnow	;
end