L = 100;  % size of lattice

% epidemiological parameters
beta = 6.2;
alpha = 1.5;

% what proportion of infection attempts are long-distance, as opposed to local
proportionLongDistance = 1;

% how many sites to make infected initially
numI0 = 50;

% how much elapsed time to let go by before recording the state of the system
recordDeltaT = 0.05;

% how many record-saving updates to do before drawing graphics
displayEvery = 1;

% we'll grow our vectors by this many elements at a time
chunksize = 500;

% black, red, blue; good for on-screen display
mycolormap = [0 0 0 ; 1 0 0 ; 0 0 1];
% black, white, gray; good for black&white printing
% mycolormap = [0 0 0 ; 1 1 1 ; .375 .375 .375];

% offsets used for local interactions
xoffsets = [0 0 -1 1];
yoffsets = [1 -1 0 0];

% vectorlengths keeps track of how long are state vectors currently are
vectorlengths = chunksize;
% Create vectors.
S = zeros(1,chunksize);
I = zeros(1,chunksize);
R = zeros(1,chunksize);
et = zeros(1,chunksize);

% stateArray will contain three values:
%   0=susceptible, 1=infectious, 2=recovered
stateArray = zeros(L);

% lists of x and y coordinates of infectious individuals
Ixv = zeros(1,L^2);
Iyv = zeros(1,L^2);

% how many individuals are currently infectious
currentI = 0;

% set up the initial population of infectious individuals
while (currentI < numI0)
    x = floor(rand*L)+1;
    y = floor(rand*L)+1;
    if (stateArray(x,y) == 0)  % if the site isn't already infectious
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % FILL IN CODE HERE TO DO THREE THINGS:
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % 1. mark that spot in the stateArray matrix
        % 2. add one to our total number currentI
        % 3. store the x and y coordinates in the Ixv and Iyv vectors
    end
end


% initialize our vectors
I(1) = currentI;
S(1) = L^2 - I(1);
R(1) = 0;
et(1) = 0;

% initialize the variables we'll actually use in the core simulation
currentTime = et(1);
currentS = S(1);
% would do "currentI = I(1);" too, but it's already set
currentR = R(1);

% keeps track of the last time we recorded the state of the system
% start out with a big enough negative value to force another update
% right away
lastTimeRecorded = -2*recordDeltaT;

i = 1;
while (currentI > 0)   % while there are still infectious individuals
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % FILL IN THESE THREE LINES
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    rIcontacts =  % total rate of infectious sites contacting others
    rIR =   % total rate of infectious sites recovering
    totalRate =   % total rate of all events

    % watch for problems, some people had errors in their code making rate=0
    if (totalRate == 0)
        error('Houston, we have a problem: totalRate is 0!')
    end

    currentTime = currentTime + exprnd(1/totalRate);  % time of next event

    % choose coordinates of the infectious individual producing the event
    % (in this simulation, all events originate with infectious sites)
    coordInd = floor(rand*currentI)+1;  % random index into coord vectors
    x = Ixv(coordInd);
    y = Iyv(coordInd);
    if (rand < rIcontacts/totalRate)   % if it's a contact / attempted infection
        % choose the target individual being contacted
        if (rand < proportionLongDistance)   % long-distance contact
            % choose coordinates at random from the entire lattice (from 1 to L)
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % FILL IN THESE TWO LINES
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            otherx = 
            othery = 
        else   % local contact
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % FILL IN CODE TO CHOOSE RANDOM NEIGHBOR
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % choose a random index from 1 to 4
            randInd = 
            % then use that as in index into the xoffsets and yoffsets vector
            % together with x and y to compute otherx and othery
            otherx =
            othery =
            % then do wraparound (or you can do it in the above 2 statements)
            otherx =
            othery =
        end
        if (stateArray(otherx,othery) == 0)   % if target is susceptible
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % FILL IN THIS SECTION
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % adjust numbers of individuals in the various states
            % (don't need to change currentR)

            % modify stateArray
            stateArray(otherx,othery) = 1;

            % add coordinates of newly-infected site to Ixv and Iyv
            ...adjust Ixv and Iyv here...


        end  % end "if".  If it wasn't true, we tried to infect a non-susceptible,
        % so we don't need to do anything
    else    % an infectious individual is recovering/dying
        % we've already chosen the x,y coords of the site recovering

        % modify stateArray
        stateArray(x,y) = 2;

        % remove coordinates of newly-recovered site from Ixv and Iyv,
        % by copying the coordinates from the ends of the vectors up to
        % fill in the "hole"
        Ixv(coordInd) = Ixv(currentI);
        Iyv(coordInd) = Iyv(currentI);

        % adjust numbers of individuals in the various states
        % (don't need to change currentS)
        currentI = currentI - 1;
        currentR = currentR + 1;
    end

    % if enough time has gone by, update record-keeping vectors
    % and possibly draw graphics
    if (currentTime - lastTimeRecorded > recordDeltaT)
        lastTimeRecorded = currentTime;
        disp(sprintf('I[%g] = %d', et(i), I(i)))
        et(i+1) = currentTime;
        S(i+1) = currentS;
        I(i+1) = currentI;
        R(i+1) = currentR;
        i = i + 1;

        % if displayEvery isn't 0, and it's time to draw stuff
        if (displayEvery && (mod(i, displayEvery) == 0))
            figure(1)
            % plot S,I,R as proportions of sites rather than numbers of sites
            % use this one for color printing
            plot(et(1:i), S(1:i)/L^2, 'k-', et(1:i), I(1:i)/L^2, 'r-', et(1:i), R(1:i)/L^2, 'b-')
            % use this one for black & white printing
            % plot(et(1:i), S(1:i)/L^2, 'k-', et(1:i), I(1:i)/L^2, 'k--', et(1:i), R(1:i)/L^2, 'k-.')

            grid on

            figure(2)
            pcolor(stateArray), shading flat, colormap(mycolormap)
            drawnow
        end

        % If our vectors have reached the length we've pre-allocated space for,
        % then "grow" them.  Setting the length this way pads with "NA" values.
        if (i == vectorlengths)
            vectorlengths = vectorlengths + chunksize;
            et = [et zeros(1,chunksize)];
            S = [S zeros(1,chunksize)];
            I = [I zeros(1,chunksize)];
            R = [R zeros(1,chunksize)];
        end
    end
end

% display final image showing the final state of the system
figure(2)
pcolor(stateArray), shading flat, colormap(mycolormap)


% when we're all done, strip off any extra elements remaining at the ends
% of the vectors, which we didn't use
et = et(1:i);
S = S(1:i);
I = I(1:i);
R = R(1:i);