ecoloop1.m

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  ECOLOOP1
%  Simulation Tool for Complex System with CA
%
%    First Coded		26	  February	1997
%    Last Modified	28    December  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.

e = [2:nl 1];
w = [nl 1:nl-1];
ee = [3:nl 1 2];
ww = [nl-1 nl 1:nl-2];
eee = [4:nl 1 2 3];
www = [nl-2 nl-1 nl 1:nl-3];
eeee = [5:nl 1 2 3 4];
wwww = [nl-3 nl-2 nl-1 nl 1:nl-4];
e5 = [6:nl 1:5];
w5 = [nl-4:nl 1:nl-5];
e6 = [7:nl 1:6];
w6 = [nl-5:nl 1:nl-6];
e7 = [8:nl 1:7];
w7 = [nl-6:nl 1:nl-7];
e8 = [9:nl 1:8];
w8 = [nl-7:nl 1:nl-8];
e9 = [10:nl 1:9];
w9 = [nl-8:nl 1:nl-9];
e10 = [11:nl 1:10];
w10 = [nl-9:nl 1:nl-10];
e11 = [12:nl 1:11];
w11 = [nl-10:nl 1:nl-11];
e12 = [13:nl 1:12];
w12 = [nl-11:nl 1:nl-12];
e13 = [14:nl 1:13];
w13 = [nl-12:nl 1:nl-13];

lnei=nei(1); % Number of neighbor
lns=ns(1); % Number of cell states

while (t < tfinal) & (sum(cp(t+1,:))>0)
	t = t+1;

% Copy the previous time step
	X(:,t+1)=X(:,t);

% Neighbor = 3
	Y=zeros(nl,1);
	Y0 = ~X(:,t+1);
	if lnei==3
		if tot(1)==1 % Totalistic
			N=X(w,t+1)+X(:,t+1)+X(e,t+1);
			for ii=1:lnei+1
				Y=Y+(N==(ii-1))*tr_code123(ii);
			end
		elseif tot(1)==2 % Outer-Totalistic
			N=X(e,t+1)+X(w,t+1);
			for ii=1:lnei Y=Y+(N==(ii-1)&~X(:,t+1))*otr_code123(ii); end
			for ii=lnei+1:2*lnei Y=Y+((N==(ii-lnei-1))&X(:,t+1))*otr_code123(ii); end
		elseif tot(1)==3 % General
			if lns==2
				Y1=X(:,t+1);
				Y = Y + (Y0(w) & Y0(:) & Y0(e))*r_code123(1);
				Y = Y + (Y0(w) & Y0(:) & Y1(e))*r_code123(2);
				Y = Y + (Y0(w) & Y1(:) & Y0(e))*r_code123(3);
				Y = Y + (Y0(w) & Y1(:) & Y1(e))*r_code123(4);
				Y = Y + (Y1(w) & Y0(:) & Y0(e))*r_code123(5);
				Y = Y + (Y1(w) & Y0(:) & Y1(e))*r_code123(6);
				Y = Y + (Y1(w) & Y1(:) & Y0(e))*r_code123(7);
				Y = Y + (Y1(w) & Y1(:) & Y1(e))*r_code123(8);
			elseif lns==3
				Y1=X(:,t+1)-2 & X(:,t+1);
				Y2=X(:,t+1)-1 & X(:,t+1);
				Y = Y + (Y0(w) & Y0(:) & Y0(e))*r_code133(1);
				Y = Y + (Y0(w) & Y0(:) & Y1(e))*r_code133(2);
				Y = Y + (Y0(w) & Y0(:) & Y2(e))*r_code133(3);
				Y = Y + (Y0(w) & Y1(:) & Y0(e))*r_code133(4);
				Y = Y + (Y0(w) & Y1(:) & Y1(e))*r_code133(5);
				Y = Y + (Y0(w) & Y1(:) & Y2(e))*r_code133(6);
				Y = Y + (Y0(w) & Y2(:) & Y0(e))*r_code133(7);
				Y = Y + (Y0(w) & Y2(:) & Y1(e))*r_code133(8);
				Y = Y + (Y0(w) & Y2(:) & Y2(e))*r_code133(9);
	
				Y = Y + (Y1(w) & Y0(:) & Y0(e))*r_code133(10);
				Y = Y + (Y1(w) & Y0(:) & Y1(e))*r_code133(11);
				Y = Y + (Y1(w) & Y0(:) & Y2(e))*r_code133(12);
				Y = Y + (Y1(w) & Y1(:) & Y0(e))*r_code133(13);
				Y = Y + (Y1(w) & Y1(:) & Y1(e))*r_code133(14);
				Y = Y + (Y1(w) & Y1(:) & Y2(e))*r_code133(15);
				Y = Y + (Y1(w) & Y2(:) & Y0(e))*r_code133(16);
				Y = Y + (Y1(w) & Y2(:) & Y1(e))*r_code133(17);
				Y = Y + (Y1(w) & Y2(:) & Y2(e))*r_code133(18);
	
				Y = Y + (Y2(w) & Y0(:) & Y0(e))*r_code133(19);
				Y = Y + (Y2(w) & Y0(:) & Y1(e))*r_code133(20);
				Y = Y + (Y2(w) & Y0(:) & Y2(e))*r_code133(21);
				Y = Y + (Y2(w) & Y1(:) & Y0(e))*r_code133(22);
				Y = Y + (Y2(w) & Y1(:) & Y1(e))*r_code133(23);
				Y = Y + (Y2(w) & Y1(:) & Y2(e))*r_code133(24);
				Y = Y + (Y2(w) & Y2(:) & Y0(e))*r_code133(25);
				Y = Y + (Y2(w) & Y2(:) & Y1(e))*r_code133(26);
				Y = Y + (Y2(w) & Y2(:) & Y2(e))*r_code133(27);
			end
		elseif tot(1)==4 %Symmetric 
			if lns==3
				N=X(e,t+1)+X(w,t+1); % Outer-Total
				E=(X(e,t+1)==X(w,t+1)); % If 1-1, E=1, if 0-2 E=0. 
				Y1=X(:,t+1)-2 & X(:,t+1);
				Y2=X(:,t+1)-1 & X(:,t+1);
				Y = Y + (Y0 & (N==0))            *sym_code133(1);
				Y = Y + (Y0 & (N==1))            *sym_code133(2);
				Y = Y + (Y0 & (N==2) & (~E)) *sym_code133(3);
				Y = Y + (Y0 & (N==2) &  E)  *sym_code133(4);
				Y = Y + (Y0 & (N==3))            *sym_code133(5);
				Y = Y + (Y0 & (N==4))            *sym_code133(6);

				Y = Y + (Y1 & (N==0))            *sym_code133(7);
				Y = Y + (Y1 & (N==1))            *sym_code133(8);
				Y = Y + (Y1 & (N==2) & (~E)) *sym_code133(9);
				Y = Y + (Y1 & (N==2) &  E)  *sym_code133(10);
				Y = Y + (Y1 & (N==3))            *sym_code133(11);
				Y = Y + (Y1 & (N==4))            *sym_code133(12);

				Y = Y + (Y2 & (N==0))            *sym_code133(13);
				Y = Y + (Y2 & (N==1))            *sym_code133(14);
				Y = Y + (Y2 & (N==2) & (~E)) *sym_code133(15);
				Y = Y + (Y2 & (N==2) &  E)  *sym_code133(16);
				Y = Y + (Y2 & (N==3))            *sym_code133(17);
				Y = Y + (Y2 & (N==4))            *sym_code133(18);
			else
				error('Invalid Rule Type');
			end
		end
	elseif lnei==5
		if tot(1)==1
			N=X(ww,t+1)+X(w,t+1)+X(:,t+1)+X(e,t+1)+X(ee,t+1);
			for ii=1:lnei+1
				Y=Y+(N==(ii-1))*tr_code125(ii);
			end
		elseif tot(1)==2
			N=X(ww,t+1)+X(w,t+1)+X(e,t+1)+X(ee,t+1);
			for ii=1:lnei Y=Y+(N==(ii-1)&~X(:,t+1))*otr_code125(ii); end
			for ii=lnei+1:2*lnei Y=Y+((N==(ii-lnei-1))&X(:,t+1))*otr_code125(ii); end
		elseif tot(1)==3
			Y1=X(:,t+1);
			Y = Y + (Y0(ww) & Y0(w) & Y0(:) & Y0(e) & Y0(ee))*r_code125(1);
			Y = Y + (Y0(ww) & Y0(w) & Y0(:) & Y0(e) & Y1(ee))*r_code125(2);
			Y = Y + (Y0(ww) & Y0(w) & Y0(:) & Y1(e) & Y0(ee))*r_code125(3);
			Y = Y + (Y0(ww) & Y0(w) & Y0(:) & Y1(e) & Y1(ee))*r_code125(4);
			Y = Y + (Y0(ww) & Y0(w) & Y1(:) & Y0(e) & Y0(ee))*r_code125(5);
			Y = Y + (Y0(ww) & Y0(w) & Y1(:) & Y0(e) & Y1(ee))*r_code125(6);
			Y = Y + (Y0(ww) & Y0(w) & Y1(:) & Y1(e) & Y0(ee))*r_code125(7);
			Y = Y + (Y0(ww) & Y0(w) & Y1(:) & Y1(e) & Y1(ee))*r_code125(8);
			Y = Y + (Y0(ww) & Y1(w) & Y0(:) & Y0(e) & Y0(ee))*r_code125(9);
			Y = Y + (Y0(ww) & Y1(w) & Y0(:) & Y0(e) & Y1(ee))*r_code125(10);
			Y = Y + (Y0(ww) & Y1(w) & Y0(:) & Y1(e) & Y0(ee))*r_code125(11);
			Y = Y + (Y0(ww) & Y1(w) & Y0(:) & Y1(e) & Y1(ee))*r_code125(12);
			Y = Y + (Y0(ww) & Y1(w) & Y1(:) & Y0(e) & Y0(ee))*r_code125(13);
			Y = Y + (Y0(ww) & Y1(w) & Y1(:) & Y0(e) & Y1(ee))*r_code125(14);
			Y = Y + (Y0(ww) & Y1(w) & Y1(:) & Y1(e) & Y0(ee))*r_code125(15);
			Y = Y + (Y0(ww) & Y1(w) & Y1(:) & Y1(e) & Y1(ee))*r_code125(16);
			Y = Y + (Y1(ww) & Y0(w) & Y0(:) & Y0(e) & Y0(ee))*r_code125(17);

			Y = Y + (Y1(ww) & Y0(w) & Y0(:) & Y0(e) & Y1(ee))*r_code125(18);
			Y = Y + (Y1(ww) & Y0(w) & Y0(:) & Y1(e) & Y0(ee))*r_code125(19);
			Y = Y + (Y1(ww) & Y0(w) & Y0(:) & Y1(e) & Y1(ee))*r_code125(20);
			Y = Y + (Y1(ww) & Y0(w) & Y1(:) & Y0(e) & Y0(ee))*r_code125(21);
			Y = Y + (Y1(ww) & Y0(w) & Y1(:) & Y0(e) & Y1(ee))*r_code125(22);
			Y = Y + (Y1(ww) & Y0(w) & Y1(:) & Y1(e) & Y0(ee))*r_code125(23);
			Y = Y + (Y1(ww) & Y0(w) & Y1(:) & Y1(e) & Y1(ee))*r_code125(24);
			Y = Y + (Y1(ww) & Y1(w) & Y0(:) & Y0(e) & Y0(ee))*r_code125(25);
			Y = Y + (Y1(ww) & Y1(w) & Y0(:) & Y0(e) & Y1(ee))*r_code125(26);
			Y = Y + (Y1(ww) & Y1(w) & Y0(:) & Y1(e) & Y0(ee))*r_code125(27);
			Y = Y + (Y1(ww) & Y1(w) & Y0(:) & Y1(e) & Y1(ee))*r_code125(28);
			Y = Y + (Y1(ww) & Y1(w) & Y1(:) & Y0(e) & Y0(ee))*r_code125(29);
			Y = Y + (Y1(ww) & Y1(w) & Y1(:) & Y0(e) & Y1(ee))*r_code125(30);
			Y = Y + (Y1(ww) & Y1(w) & Y1(:) & Y1(e) & Y0(ee))*r_code125(31);
			Y = Y + (Y1(ww) & Y1(w) & Y1(:) & Y1(e) & Y1(ee))*r_code125(32);
		end
	elseif lnei==7
		if tot(1)==1
			N=X(www,t+1)+X(ww,t+1)+X(w,t+1)+X(:,t+1)+X(e,t+1)+X(ee,t+1)+X(eee,t+1);
			for ii=1:lnei+1
				Y=Y+(N==(ii-1))*tr_code127(ii);
			end
		elseif tot(1)==2
			N=X(www,t+1)+X(ww,t+1)+X(w,t+1)+X(e,t+1)+X(ee,t+1)+X(eee,t+1);
			for ii=1:lnei Y=Y+(N==(ii-1)&~X(:,t+1))*otr_code127(ii); end
			for ii=lnei+1:2*lnei Y=Y+((N==(ii-lnei-1))&X(:,t+1))*otr_code127(ii); end
		elseif tot(1)==3
			Y1=X(:,t+1);
			Y = Y + (Y0(www) & Y0(ww) & Y0(w) & Y0(:) & Y0(e) & Y0(ee) & Y0(eee))*r_code127(1);