% succession3.m
%    Markov Chain model for landscape succession
%            from Ecobeaker - Math 151 
% This is to generate a simple landscape succession 
% model similar to Ecobeaker in which there are 
% 3 species (Bare soil, Grass and Shrub) with
% the transition matrix given by the matrix P 
%   The program shows a landscape of size 20 x 20 with
% initially all in Bare soil and then each location transitions 
% according to the matrix probabilities to a different type
% with the landscape displayed along with a bar chart of the 
% number of elements of the matrix of each type (e.g.
% total area in each vegetation type at each time ), 
% this version steps through each time step showing the
% landscape and bar chart rather than doing an animation
%
% the matrix which holds the landscape is L
%
%This version steps through rather than shows an animation
%
L=zeros(20,20);
%
% generate a random matrix for each location
%
%
%define the transition matrix such as
%
P=[.9 0 0; .1 .95 0; 0 .05 1 ]
%
% where the elements give transition probabilities from 
% the row to the column
%
% Show the first landscape
%
S=[400 0 0];
figure(1)
bar(S,1.5,'EraseMode','xor')
figure(2)
image(L,'EraseMode','xor'),pause
%
% k is the number of time steps
%
for k=1:20
%
% R holds the random numbers for each locations transition
% and S holds the current state vector
%
R=rand(20,20);
S=zeros(3);
%
% The landscape matrix element values are
%    0 for Bare Soil (state 1)
%    50 for Grass (state 2)
%    100 for Shrub (state 3)
% 
for i=1:20
    for j=1:20
%
% it is a check as to whether L(i,j) has been modified
%
      it=0;
      if L(i,j)==0 & it==0
           if R(i,j) <= P(1,1)
               L(i,j)=0;
           else 
                    if R(i,j) < P(1,1)+P(2,1)
                         L(i,j)=50;
                    else L(i,j)=100;
                    end;
           end;
      it=1;
      else;
      end;
      if L(i,j)==50 & it==0
           if R(i,j) < P(1,2)
               L(i,j)=0;
           else 
                     if R(i,j) < P(1,2)+P(2,2)
                         L(i,j)=50;
                     else L(i,j)=100;
                     end;
           end;
      it=1;
      else;
      end;
      if L(i,j)==100 & it==0
           if R(i,j) < P(1,3)
               L(i,j)=0;
           else 
                     if R(i,j) < P(1,3)+P(2,3)
                          L(i,j)=50;
                    else L(i,j)=100;
                    end;
          end;
       else;              
       end;
     if L(i,j)==0
        S(1)=S(1)+1;
        else
          if L(i,j)==50
           S(2)=S(2)+1;
           else
             S(3)=S(3)+1;
          end;
     end;   
   end;
end;
figure(1)
bar(S,2,'EraseMode','xor')
figure(2)
image(L,'EraseMode','xor'),pause
k
end;