function [nrecovered,infected,recovered,susceptible] = ...
             modelmc1d(n,k,tau)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function [nrecovered,infected,recovered,susceptible]   %
%                      = modelmc1d(n,k,tau)              %
%  for                                                   %
%  A solution to CSE Your Homework Assignment Project 5  %
%  Models of Infection: Person to Person                 %
%  Dianne P. O'Leary   09/03                             %
%                                                        %
% This model is used in Problem 6.  It is a simplified   %
% version of the model for Problems 1-3, with no         %
% mobility or vaccination, and only a 1-dimensional      %
% grid of beds arranged in a straight line.              %
%                                                        %
% The parameters for the model:                          %
%                                                        %
%   n:       dimension of the grid of hospital beds      %
%   k:       number of days of contagion                 %
%   tau:     transmission rate for the infection         %
%                                                        %
% We model the spread of an infection through a          %
% population, assuming that a person once infected       %
% cannot be reinfected.                                  %
%                                                        %
% Output:                                                %
%                                                        %
%                                                        %
%  nrecovered  = number of patients eventually infected  %
%  infected(t) = number of individuals infected on day t %
%  recovered(t) = number of individuals who have         %
%                 recovered from infection on day t      %
%                 and are immune                         %
%  susceptible(t) = number of individuals susceptible to %
%                   infection on day t                   %
%  Note that infected(t) + recovered(t) + susceptible(t) %
%            = 3.                                        %
%                                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Pop(i)   =  0  means not infected                      %
%          =  j  means in j-th day of infection          %
%          = -1  means recovered                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


Pop = zeros(1,n);
nPop = zeros(1,n);

Pop(1,2) = 1;
ninfected = 1;

t = 0;

while (ninfected > 0)

   t = t + 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% For the next timestep, calculate the infection status  %
% of each individual.                                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   for j=1:n,
      new = Pop(1,j);
      if (new > 0)
         new =  new + 1;
         if (new > k)  % if past the time of active infection 
            new = -1;  % then recovered
         end
      elseif (new == 0)
         new = infect(Pop,1,j,1,n,tau);
      end
      nPop(1,j) = new;
   end

   Pop = nPop;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Counts                                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   infected(t)    = sum(sum(Pop >   0));
   recovered(t)   = sum(sum(Pop == -1));
   susceptible(t) = sum(sum(Pop ==  0));
   ninfected = infected(t);

end

nrecovered = recovered(t);