clear all close all N=20; %Number of rows and columns in the grid timesteps=700; %Number of time steps amax=60; %Maximum time steps F. cupressoides can live bmax=30; %Maximum time steps Nothofagus spp. can live disturbedratio=0.5; %Probability a cell is initially disturbed Fmature=7; %Age of maturity for F. cupressoides Nmature=2; %Age of maturity for Nothofagus spp. catastropheprob=0.002; %Probability a catastrophe happens in a cell at each time step %make blank "matrices" amatrix=zeros(N); %Creates a matrix of colonization times for F. cupressoides bmatrix=zeros(N); smatrix=zeros(N); %initialize: %Half of sites are disturbed smatrix(find(rand(N)0.5))=1 %Is the cell occupied by Nothofagus spp.? for i = 1:timesteps amatrixnew=zeros(N); bmatrixnew=zeros(N); smatrixnew=zeros(N); catastrophematrix=rand(N)0) ) ; amatrixnew(indices)=amatrix(indices)+1; %then, plants just get one year older (occupation time increases) indices=find(catastrophematrix==0 & smatrix==0 & amatrix==0) ; %candidates for occupation %OK, need to find nearest-neighbor indices, nnindices for j=1:length(indices) index=indices(j); nnindices=[]; if index==1 nnindices=[index+1,index+N]; elseif index<=N nnindices=[index-1,index+1,index+N]; elseif index<=(N-1)*N nnindices=[index-N,index-1,index+1,index+N]; elseif index<=N^2-1 nnindices=[index-N,index-1,index+1]; else %index=N^2 nnindices=[index-N,index-1]; end %compute occupation probability, based on those nnindices number_a=sum( amatrix(nnindices)>Fmature); proba_a=number_a/4; %with this proba, cell is colonized (i.e., its a value is set to 1) -- else it remains at 0! amatrixnew(index)=rand0 & bmatrixFmature); proba_a=number_a/4; %with this proba, cell is colonized (i.e., its a value is set to 1) -- else it remains at 0! amatrixnew(index)=rand0 & bmatrixNmature); proba_b=number_b/4; %with this proba, cell is colonized (i.e., its a value is set to 1) -- else it remains at 0! bmatrixnew(index)=rand15); smatrixnew(indices)=1; smatrix=smatrixnew amatrix=amatrixnew bmatrix=bmatrixnew meanFdensity(i)=length(find(amatrix>0))/N^2; meanNdensity(i)=length(find(bmatrix>0))/N^2; % figure(1) % bar3(smatrix); drawnow; % figure(2) % bar3(amatrix); drawnow; figure(1) subplot(121) set(gca,'FontSize',12) bar3(amatrix); title('Age of plant F at spatial site') subplot(122) set(gca,'FontSize',12) bar3(bmatrix); title('Age of plant N at spatial site') drawnow; end figure set(gca,'FontSize',20) plot(meanFdensity); hold on plot(meanNdensity,':'); xlabel('timesteps') ylabel('density') legend('F. cupressoides','Nothofagus spp.')