%%%%%Pred/Prey cellar automata model

%%%%Initial conditions
clear all
close all


NL =  10;  %%Spatial dimension of grid box
WI =  0.1;   %%Initial Number of spatial boxes with wolves
RI =  0.5;  %%Initial Number of spatial boxes with rabbits
G  =  rand(NL);  %%The initial grid
TS = 1000;

time = 1:TS;

NW = zeros(1,TS);  %%% Number of wolves
NR = zeros(1,TS);  %%% Number of Rabbits
NE = zeros(1,TS);  %%% Number of Empties
count = 0;

% for animation:
z=zeros(NL,NL); imh = image(cat(3,z,z,z));
set(imh, 'erasemode', 'none'); axis equal; axis tight


%%%SETTING UP INITIAL GRID
for i = 1:NL;
  for j = 1:NL;
    count = count +1;
    if G(i,j) > (WI + RI)  %%%EMPTY GRID SPOT
      G(i,j) = 0;
      NE(1,1)= NE(1,1) + 1;
    elseif G(i,j) > (WI)   %%%RABBITS IN GRID SPOT
      G(i,j) = 1;
      NR(1,1)= NR(1,1) + 1;
    else
      G(i,j) = 2;        %%%WOLVES IN GRID SPOT
      NW(1,1)= NW(1,1) + 1;
    end
  end
end

Ginitial = G;
%%% INTEGRATE THE MODEL FORWARD IN TIME (three for loops: 1st time, then i,
%%%then j

%%% Define tansition probabilities:
%%% For an empty cell:
Prborn  = 0.05;  %%%The probability of an empty cell next to at least 1 rabbit becoming a rabbit
%%%Rabbit
Peaten = 0.03;  %%%The probability of an rabbit next to a wolf gets eaten
Pwborn = 0.25;  %%%If a rabbit was just eaten by a wolf, then the probability that the wolves will produce babies
               %%%NOTE:  If you are not eaten, then you stay a rabbit
Pwdeath = 0.005; %%%The probability that a wolf just dies

for t = 2:TS;
  for i = 1:NL;
    for j = 1:NL;
      %%%Assign check indicies (checkL = check left, checkR =
      %%%check right, checkU = check up, checkD = checkD
      checkL = i - 1;
      checkR = i + 1;
      checkU = j - 1;
      checkD = j + 1;
      
      %%%Re-assign check indices for looping boundary conditions
      if checkU == 0;
        checkU = NL;
      end
      if  checkD == NL+1;
        checkD = 1;
      end
      if checkR == NL + 1;
        checkR = 1;
      end
      if checkL == 0;
        checkL = NL;
      end
      
      if G(i,j) == 0  %%%Empty rules apply-->check for rabbits
        if G(checkL,j) == 2;
          if rand(1) > 1 - Prborn;  %%%Rabbit birth
            G(i,j) = 2;
          end
        elseif G(checkR,j) == 2;
          if rand(1) > 1 - Prborn;
            G(i,j) = 2;
          end
        elseif G(i,checkU) == 2
          if rand(1) > 1 - Prborn;
            G(i,j) = 2;
          end
        elseif G(i,checkD) == 2
          if rand(1) > 1 - Prborn;
            G(i,j) = 2;
          end
        end
                
      elseif G(i,j) == 2 %%%Rabbit rules apply-->check for wolves
        if G(checkL,j) == 1;  %%%Checking for a wolf
          if rand(1) > 1 - Peaten;
            G(i,j) = 0;
            %% If the rabbits are eaten, then wolves are born with the probabilty
            %% = Pwborn
            if rand(1) > 1 - Pwborn
              G(i,j) = 1;
            end
          end
        end
                    
        if G(checkR,j) == 1;  %%%Checking for a wolf
          if rand(1) > 1 - Peaten;
            G(i,j) = 0;
            %% If the rabbits are eaten, then wolves are born with the probabilty
            %% = Pwborn
            if rand(1) > 1 - Pwborn
              G(i,j) = 1;
            end
          end
        end
        if G(i,checkU) == 1;  %%%Checking for a wolf
          if rand(1) > 1 - Peaten;
            G(i,j) = 0;
            %% If the rabbits are eaten, then wolves are born with the probabilty
            %% = Pwborn
            if rand(1) > 1 - Pwborn
              G(i,j) = 1;
            end
          end
        end
        if G(i,checkD) == 1;  %%%Checking for a wolf
          if rand(1) > 1 - Peaten;
            G(i,j) = 0;
            %% If the rabbits are eaten, then wolves are born with the probabilty
            %% = Pwborn
            if rand(1) > 1 - Pwborn
              G(i,j) = 1;
            end
          end 
        end
        
      elseif G(i,j) == 1 %%%Wolf rules apply-->constant death rate
        if rand(1) > 1 - Pwdeath
          G(i,j) = 0;
        end
      end
    
      if 0
        %% implement random movement of wolves and rabbits:  
        velocity=ceil(4*rand(1));
        %% velocity: 1: right; 2: left; 3: up; 4: down;
        if     velocity==1; imove=i-1;jmove=j;
        elseif velocity==2; imove=i+1;jmove=j;
        elseif velocity==3; imove=i;  jmove=j+1;
        elseif velocity==4; imove=i;  jmove=j-1;
        end
        %% periodic b.c. for motion:
        if imove==NL+1; imove=1; end
        if imove==0; imove=NL; end
        if jmove==NL+1; jmove=1; end
        if jmove==0; jmove=NL; end
      
        %% move if rabbit or wolf, but only if cell to which we need
        %to move is empty:
        if (G(i,j)==1||G(i,j)==2)&&G(imove,jmove)==0
          G(imove,jmove)=G(i,j);
          G(i,j)=0;
        end
      end  %% end of random motion

    end  %%%j loop end
  end      %%%i loop end
  
  %%%Count the number of wolves and rabbits at the end of this time step
  numberwolves = 0;
  numberrabbits = 0;
  for r = 1:NL;
    for c = 1:NL;
      if G(r,c) == 1;
        numberwolves = numberwolves + 1;
      elseif G(r,c) == 2;
        numberrabbits = numberrabbits + 1;
      end
    end
  end
  
  NW(1,t) = numberwolves;
  NR(1,t) = numberrabbits;

  %% G = 0 - empty; 1 - wolf; 2 - rabbit:
  set(imh, 'cdata', cat(3,(G==1),(G==2),z))
  drawnow
  %% contour(G);
  pause(0.01)
end   %%%t loop end

%Ginitial
%G
figure
plot(time,NW,'r',time,NR,'g','linewidth',3)
xlabel('Time')
ylabel('Number of Wolves and Rabbits')
legend('Wolves','Rabbits')
