function [p_n,rho_n] = machine_repair_model(lambda,mu,M,s)
% MACHINE_REPAIR_MODEL - Returns the rho_n, p_n, and p_0 for a machine repair model
%   
% Input: 
%  lambda: arrival rate
%  mu: departure rate
%  M: the number of machines 
%  s: the number of servers
% 
% Output: 
%   p_0: is the steady-state probabity we are in the state X(t)=0
%   p_n: is the steady-state probabity we are in the state X(t)=n
%   rho_n: the proportionality factor for p_n = rho_n p_0 
% 
% References: 
% 
% Introduction to Probability Models: Chapter 10
% by Roe Goodman 
% 
% Written by:
% -- 
% John L. Weatherwax                2007-09-10
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

rho = lambda/mu;

% initialize space (rho_n(1) is \rho_0=1): 
rho_n = ones(1,M+1);

% explicitly compute the first value of \rho_n: \rho_1 = (M/1!) rho
n=1; 
num        = M;
den        = 1; 
rho_n(n+1) = (num/den)*rho;
% fill the elements for n=2:s (\rho_{2:s}):
for n=2:s,
  num        = num*(M-n+1);
  den        = n*den; 
  rho_n(n+1) = (num/den)*(rho^n);
end
% fill the elements for n=s+1:M
for n=s+1:M,
  num        = num*(M-n+1); 
  den        = factorial(s)*(s^(n-s));
  rho_n(n+1) = (num/den)*(rho^n);
end

% compute the normilization probability: 
p_0 = 1/sum(rho_n); 

% compute the steady state probabilities: 
p_n = p_0*rho_n;