function sir = sir_simulation ( m, n, a, k, tau, t_max )

%*****************************************************************************80
%
%% sir_simulation() simulates an SIR model of disease transmission.
%
%  Discussion:
%
%    We assume that a hospital ward comprises an array of M by N beds.
%
%    The status of each patient is recorded as an integer in an array A.
%
%    Susceptible patients, with a status of 0, have never had the disease.
%
%    Infected patients, with a positive status between 1 and K, have
%    had the disease for A(I,J) days.
%
%    Recovered patients, with a status of -1, have had the disease for K
%    days, are no longer infected, and cannot get the disease again.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    16 April 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Dianne OLeary,
%    Models of Infection: Person to Person,
%    Computing in Science and Engineering,
%    Volume 6, Number 1, January/February 2004.
%
%    Dianne OLeary,
%    Scientific Computing with Case Studies,
%    SIAM, 2008,
%    ISBN13: 978-0-898716-66-5,
%    LC: QA401.O44.
%
%  Input:
%
%    integer M, N, the number of rows and columns of beds.
%
%    integer A(M,N), status of each patient on day T = 1.
%
%    integer K, the maximum number of days of infection.
%
%    real TAU, the probability of transmission of the disease
%    over one day, due to one neighboring patient who is infected.
%
%    integer T_MAX, the number of days to model.
%
%  Output:
%
%    real SIR(3,T_MAX), the relative numbers of susceptible,
%    infected, and recovered patients.  The sum of these three values
%    will be 1 for each day.
%
% rng ( 'shuffle' );

  sir = zeros ( 3, t_max );

  for t = 1 : t_max
%
%  Plot the current data.
%
    timestep_display ( m, n, a, k, t );
%
%  Compute SIR.
%
    sir(1,t) = sum ( sum ( a == 0 ) );
    sir(2,t) = sum ( sum ( 1 <= a && a <= k ) );
    sir(3,t) = sum ( sum ( a == -1 ) );
%
%  Update the patient status for the next day.
%
    if ( t == t_max )
      break
    end

    a_new = zeros ( m, n );

    for i = 1 : m
      for j = 1 : n
%
%  Recovered patients never change.
%
        if ( a(i,j) == -1 )
          a_new(i,j) = -1;
%
%  Infected, and less than K days, increase by 1.
%
        elseif ( 1 <= a(i,j) && a(i,j) < k )
          a_new(i,j) = a(i,j) + 1;
%
%  Infected K days becomes recovered.
%
        elseif ( a(i,j) == k )
          a_new(i,j) = -1;
%
%  Susceptible, look for infected neighbors.
%
        else

          a_new(i,j) = 0;

          if ( 1 < i && 0 < a(i-1,j) && a(i-1,j) <= k )
            if ( rand ( 1, 1 ) <= tau )
              a_new(i,j) = 1;
            end
          end
          if ( i < m && 0 < a(i+1,j) && a(i+1,j) <= k )
            if ( rand ( 1, 1 ) <= tau )
              a_new(i,j) = 1;
            end
          end
          if ( 1 < j && 0 < a(i,j-1) && a(i,j-1) <= k )
            if ( rand ( 1, 1 ) <= tau )
              a_new(i,j) = 1;
            end
          end
          if ( j < n && 0 < a(i,j+1) && a(i,j+1) <= k )
            if ( rand ( 1, 1 ) <= tau )
              a_new(i,j) = 1;
            end
          end
     
        end

      end
    end 
%
%  Copy new data into A.
%
    a = a_new;

  end
%
%  Normalize SIR.
%
  sir(1:3,1:t_max) = sir(1:3,1:t_max) / m / n;

  return
end
