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)<disturbedratio))=1
%populate others with one or other type of plant, with age 1
R=rand(N);
amatrix(find(R<0.3))=1 %Is the cell occupied by F. cupressoides?
bmatrix(find(R>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)<catastropheprob;
    
    %Is there a catastrophe in the cell?
    %If there is a catastrophe...
    amatrixnew(find(catastrophematrix==1))=0; %Soil is disturbed, all vegetation dies
    bmatrixnew(find(catastrophematrix==1))=0;
    smatrixnew(find(catastrophematrix==1))=0;
    
    %If there is not a catastrophe...

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %AND soil disturbed
    

    %first, choose indices where this is true
    indices=find(catastrophematrix==0 & smatrix==0) ;
    bmatrixnew(indices)=0;  %set b to zero, Nothofagus spp. cannot grow in the cell
    
    %update rule for a -- further special cases
    indices=find(catastrophematrix==0 & smatrix==0 & amatrix==amax) ; 
    amatrixnew(indices)=0; %if occupation time has reach amax then F. cupressoides dies
    
    indices=find( (catastrophematrix==0 & smatrix==0 & amatrix<amax & amatrix>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)=rand<proba_a;
    end
        
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %next cases -- s=0, if no catastrophe and soil ***NOT*** disturbed
    
    
    %update rule for a -- further special cases
    indices=find(catastrophematrix==0 & smatrix==1 & (amatrix==amax|bmatrix==bmax) ) ; 
    amatrixnew(indices)=0; %if either occupation time has reached max then F. cupressoides dies
    
    indices=find( (catastrophematrix==0 & smatrix==1 & amatrix<amax & amatrix>0 & bmatrix<bmax) ) ; 
    amatrixnew(indices)=amatrix(indices)+1; %then, plants just get one year older (occupation time increases)
    

    indices=find( (catastrophematrix==0 & smatrix==1 & amatrix==0 & bmatrix==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)=rand<proba_a;
    end

    
    %%%%%%%%
    %update rule for b w/no catastrophe, and soil not disturbed -- further special cases
    indices=find(catastrophematrix==0 & smatrix==1 & bmatrix==bmax ) ; 
    bmatrixnew(indices)=0; %if occupation time has reached max then die
    
    indices=find( (catastrophematrix==0 & smatrix==1 &  bmatrix>0 & bmatrix<bmax) ) ; 
    bmatrixnew(indices)=bmatrix(indices)+1; %then, plants just get one year older (occupation time increases)
    

    indices=find( (catastrophematrix==0 & smatrix==1 & bmatrix==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_b=sum( bmatrix(nnindices)>Nmature);
        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)=rand<proba_b;
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    
    %finally, update soil disturbance 
    smatrixnew(find(catastrophematrix==0))=smatrix(find(catastrophematrix==0))    ;
    %If F. cupressoides occupation
    %time is 15 then soil in cell is no longer disturbed
    
    indices=find(amatrix>15);
    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.')