function  outcomes = fire(numRows, numCols, startRow, startCol, trials)
    % Simulate the spread of a fire.
    % USAGE: outcomes = fire(numRows, numCols, startRow, startCol, trials)
    %    The simulation will be performed on a grid with
    %    specified number of rows and colums, assumed to be numbered
    %    starting at (1,1) at the top-left.  The fire begins at location
    %    (startRow, startCol) within that grid.
    %    
    %    An individual trial will have one of five outcomes:
    %    1) it reaches the right-side of the grid
    %    2) it reaches the bottom-side of the grid
    %    3) it reaches the left-side of the grid
    %    4) it reaches the top-side of the grid
    %    5) it burns out within the grid
    %    The returned outcome will be a row-vector of length five
    %    designating the percentage of trials that resulted in
    %    such outcomes. For example, the value outcomes(5) is the
    %    percentage of trials that burned out.
    %
    %    input argument trials specifies the number of independent
    %    trials.  If not specified, one trial will be performed.
    %
    %    When only one trial is to be performed, the function animates
    %    the spread of the fire.
    if nargin < 5
        trials = 1;
    end
    
    % set to 4 for regular motion, 8 for extra credit
    numDirections = 4;
    
    outcomes = zeros(1,5);
    for t = 1:trials
        mark = ones(numRows, numCols);
        r = startRow; c = startCol;
        if trials == 1
            mark(r, c) = 2;
            plotGrid(mark);
            pause(.1);
        end
        mark(r, c) = 3;
        
        done = false;
        while ~done
            % determine legal directions of motion
            possible = [];
            for dir = 1:numDirections
                [tr tc] = move(dir,r,c);
                if mark(tr,tc) == 1  % legal move
                possible = [possible dir];
                end
            end
    
            % pick from the legal possibilities, if any
            if length(possible) > 0
                ran = ceil(rand()*length(possible));
                choice = possible(ran);
                [r c] = move(choice, r, c);
            end   
    
            % see if we have reached a final outcome
            location = checkPosition(numRows, numCols, r, c);
            if mark(r,c) ~= 1 || location(5) ~= 1
                done = true;   % either stuck or on boundary
                outcomes = outcomes + location;
            end
    
            % create frame of animation
            if trials == 1
                mark(r, c) = 2; % momentarily orange
                plotGrid(mark);
                pause(.1);
            end
            mark(r, c) = 3;     % mark as gray
        end
    end
    outcomes = outcomes ./ trials;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function outcome = checkPosition(nr, nc, r, c)
    outcome = zeros(1,5);
    if c == nc
        outcome(1) = 1;   % right edge
    end
    if r == nr
        outcome(2) = 1;   % bottom edge
    end
    if c == 1
        outcome(3) = 1;   % left edge
    end
    if r == 1
        outcome(4) = 1;   % top edge
    end
  
    total = sum(outcome);
    if total == 0         % within interior
        outcome(5) = 1;   
    elseif total > 1      % at a corner
        outcome = outcome ./ total;
    end
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [newR newC] = move(dir, r, c)
    delta = [  0  +1;     % right one column
              +1   0;     % down one row
               0  -1;     % left one column
              -1   0;     % up one row
                          % And, for extra credit:
              +1  +1;     % down/right
              +1  -1;     % down/left
              -1  -1;     % up/left
              -1  +1;     % up/right
            ];
    newR = r + delta(dir,1);
    newC = c + delta(dir,2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function plotGrid(mark)
    colormap([.035 .200 .153; 1.00 .367 .063; .561 .510 .592])
    image(mark);
    axis equal;
    s = size(mark);
    axis([0.5 s(2)+0.5 0.5 0.5+s(1)]);
