% Simulate sharks and fish
% Inputs (run sfinit to get a sample set of inputs)
%   d = size of square, toroidally connected ocean
%   n = number of time steps to simulate
%   ofish(d,d) = initial distribution of fish  (1 for fish, 0 for no fish)
%   oshark(d,d) = initial distribution of sharks  (1 for shark, 0 for no shark)
%
% Output
%   plot of shark ('+') and fish ('o') at the end of every time step
%   (computation pauses after every plot; hit return to continue)
%
% Other routines used:
%      UpdateBoundaries
%      GetSharkDirection
%      StarveAgeSharks
%      UpdateBoundaries
%      MoveSharkRight
%      MoveSharkLeft
%      MoveSharkUp
%      MoveSharkDown
%      GetFishDirection
%      AgeFish
%      UpdateBoundaries
%      MoveFishRight
%      MoveFishLeft
%      MoveFishUp
%      MoveFishDown
%      PlotSharkFish
%
% Main data structures:
%
%   shark(i,j) = 1 if shark present at (i,j), 0 otherwise
%   sharkage(i,j) = age of shark at (i,j), 0 if none present
%   sharkeat(i,j) = time since shark's last meal, 0 if none present
%   fish(i,j) = 1 if fish present at (i,j), 0 otherwise
%   fishage(i,j) = age of fish at (i,j), 0 if none present
%
%   All these are of dimension d+2 by d+2, with an extra row and column
%   on all edges to accomodate toroidal connection
%
% Rules of the system
%
%   At each time step, first the sharks are updated and then the fish.
%
%   For each shark, a direction is computed, which includes
%      a random vector
%      an "electrostatic" attraction (proportional to SHARKATTRACT) to all 
%           the fish
%      an attraction to any fish who are nearest neighbors (proportional 
%           to EATNOW)
%   If the shark can't move in the desired direction, it stays put.
%   If the age of the shark exceeds SBREED, and it moves, it leaves a new shark behind.
%   If a shark has not eaten for STARVE time steps, it dies.
%
%   For each fish, a direction is computed, which includes
%      a random vector
%      an "electrostatic" repulsion (proportional to FISHREPEL) away from all 
%            the shark
%   If the fish can't move in the desired direction, it stays put.
%   If the age of the fish exceeds FBREED, and it moves, it leaves a new fish behind.
%
%   Initialize fish and shark data
    fish=[zeros(1,d+2);[zeros(d,1),ofish,zeros(d,1)];zeros(1,d+2)];
    shark=[zeros(1,d+2);[zeros(d,1),oshark,zeros(d,1)];zeros(1,d+2)];
    fishage=zeros(d+2,d+2);
    sharkage=zeros(d+2,d+2);
    sharkeat=zeros(d+2,d+2);
    sharkdir=zeros(d+2,d+2);
    fishdir=zeros(d+2,d+2);
    UpdateBoundaries
%
%   Assign constants
    rand('normal')
    II = sqrt(-1);
%   Age at which Sharks breed (sharks breed if they are old enough and can move)
    SBREED = 7;
%   Age at which Fish breed (fish breed if they are old enough and can move)
    FBREED = 0;
%   Degree to which fish are ``electrostatically'' repelled by sharks
    FISHREPEL = 1;
%   Degree to which sharks are ``electrostatically'' attracted by fish
    SHARKATTRACT = 2;
%   Degree to which a shark wants to each a neighboring fish
    EATNOW = 100;
%   Length of time a shark can go without eating and without starving
    STARVE = 1;
%   ZZ stores locations of each cell in ocean as complex number
    [XX,YY]=meshgrid(1:d,1:d);
    ZZ=XX+II*YY;
%   Turn debugging output off; set dbug=1 to turn it on
    dbug=0;
%
%   Loop over time steps
    for i=1:n,
       disp(['Computing time step i= ',num2str(i)])
%
       UpdateBoundaries
       GetSharkDirection
       StarveAgeSharks
       UpdateBoundaries
       MoveSharkRight
       MoveSharkLeft
       MoveSharkUp
       MoveSharkDown
%
       GetFishDirection
       AgeFish
       UpdateBoundaries
       MoveFishRight
       MoveFishLeft
       MoveFishUp
       MoveFishDown
%
%      Plot output
       PlotSharkFish
    end