--- a/main/queueing/ChangeLog	Wed Jan 15 22:03:11 2014 +0000
+++ b/main/queueing/ChangeLog	Fri Jan 17 09:41:11 2014 +0000
@@ -5,6 +5,7 @@
 	* qnom accepts classes with zero arrival rate
 	* erlangb added
 	* erlangc added
+	* engset added
 
 2013-06-29 Moreno Marzolla <moreno.marzolla@unibo.it>
 

--- a/main/queueing/doc/singlestation.txi	Wed Jan 15 22:03:11 2014 +0000
+++ b/main/queueing/doc/singlestation.txi	Fri Jan 17 09:41:11 2014 +0000
@@ -30,6 +30,9 @@
 @menu
 * The M/M/1 System::    Single-server queueing station.
 * The M/M/m System::    Multiple-server queueing station.
+* The Erlang-B Formula::
+* The Erlang-C Formula::
+* The Engset Formula::
 * The M/M/inf System::  Infinite-server (delay center) station.
 * The M/M/1/K System::  Single-server, finite-capacity queueing station.
 * The M/M/m/K System::  Multiple-server, finite-capacity queueing station.
@@ -114,10 +117,6 @@
 
 @GETHELP{qsmmm}
 
-@GETHELP{erlangb}
-
-@GETHELP{erlangc}
-
 @noindent @strong{REFERENCES}
 
 @noindent G. Bolch, S. Greiner, H. de Meer and K. Trivedi, @cite{Queueing Networks
@@ -129,11 +128,37 @@
 @auindex de Meer, H.
 @auindex Trivedi, K.
 
+@c
+@c Erlang-B
+@c
+@node The Erlang-B Formula
+@section The Erlang-B Formula
+
+@GETHELP{erlangb}
+
+@noindent @strong{REFERENCES}
+
 @noindent G. Zeng, @cite{Two common properties of the erlang-B function, erlang-C function, and Engset blocking function}, Mathematical and Computer Modelling, Volume 37, Issues 12-13, June 2003, Pages 1287-1296
 
 @auindex Zeng, G.
 
 @c
+@c Erlang-c
+@c
+@node The Erlang-C Formula
+@section The Erlang-C Formula
+
+@GETHELP{erlangc}
+
+@c
+@c Engset
+@c
+@node The Engset Formula
+@section The Engset Formula
+
+@GETHELP{engset}
+
+@c
 @c M/M/inf
 @c
 @node The M/M/inf System

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/queueing/inst/engset.m	Fri Jan 17 09:41:11 2014 +0000
@@ -0,0 +1,126 @@
+## Copyright (C) 2014 Moreno Marzolla
+##
+## This file is part of the queueing toolbox.
+##
+## The queueing toolbox is free software: you can redistribute it and/or
+## modify it under the terms of the GNU General Public License as
+## published by the Free Software Foundation, either version 3 of the
+## License, or (at your option) any later version.
+##
+## The queueing toolbox is distributed in the hope that it will be
+## useful, but WITHOUT ANY WARRANTY; without even the implied warranty
+## of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with the queueing toolbox. If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+##
+## @deftypefn {Function File} {@var{B} =} engset (@var{A}, @var{m}, @var{n})
+##
+## @cindex Engset loss formula
+##
+## Compute the Engset blocking probability @math{P_b(A, m, n)} for a system
+## with a finite population of @math{n} users, @math{m} identical
+## servers, no queue, individual service rate @math{\mu}, individual
+## arrival rate @math{\lambda} (i.e., the time until a user tries to
+## request service is exponentially distributed with mean @math{1 /
+## \lambda}), and offered load @math{A = \lambda / \mu}.
+## 
+## @iftex
+##
+## @math{P_b(A, m, n)} is defined for @math{n > m} as:
+##
+## @tex
+## $$
+## P_b(A, m, n) = {\displaystyle{A^m {{n}\choose{m}}}} \over {\displaystyle{\sum_{k=0}^m {A^k {{n}\choose{k}}}}}
+## $$
+## @end tex
+##
+## and is 0 if @math{n @leq{} m}.
+## @end iftex
+##
+## @strong{INPUTS}
+##
+## @table @var
+##
+## @item A
+## Offered load, defined as @math{A = \lambda / \mu} where
+## @math{\lambda} is the mean arrival rate and @math{\mu} the mean
+## service rate of each individual server (real, @math{A > 0}).
+##
+## @item m
+## Number of identical servers (integer, @math{m @geq{} 1}). Default @math{m = 1}
+##
+## @item n
+## Number of requests (integer, @math{n @geq{} 1}). Default @math{n = 1}
+##
+## @end table
+##
+## @strong{OUTPUTS}
+##
+## @table @var
+##
+## @item B
+## The value @math{P_b(A, m, n)}
+##
+## @end table
+##
+## @var{A}, @var{m} or @math{n} can be vectors, and in this case, the
+## results will be vectors as well.
+##
+## @seealso{erlangb, erlangc}
+##
+## @end deftypefn
+
+## Author: Moreno Marzolla <moreno.marzolla(at)unibo.it>
+## Web: http://www.moreno.marzolla.name/
+function P = engset(A, m, n)
+  if ( nargin < 1 || nargin > 3 )
+    print_usage();
+  endif
+
+  ( isnumeric(A) && all( A(:) > 0 ) ) || error("A must be positive");
+  
+  if ( nargin < 2 )
+    m = 1;
+  else
+    ( isnumeric(m) && all( fix(m(:)) == m(:)) && all( m(:) > 0 ) ) || error("m must be a positive integer");
+  endif
+
+  if ( nargin < 3 )
+    n = 1;
+  else
+    ( isnumeric(n) && all( fix(n(:)) == n(:)) && all( n(:) > 0 ) ) || error("n must be a positive integer");
+  endif
+  
+  [err A m n] = common_size(A, m, n);
+  if ( err )
+    error("parameters are not of common size");
+  endif
+
+  P = arrayfun( @__engset_compute, A, m, n);
+endfunction
+
+## Compute P_b(A,m,n) recursively
+function P = __engset_compute(A, m, n)
+  if ( m >= n )
+    P = 0.0;
+  else
+    P = 1.0;
+    for i = 1:m
+      P=(A*(n-i)*P)/(i+A*(n-i)*P);
+    endfor
+  endif
+endfunction
+
+%!test
+%! fail("erlangb(1, -1)", "positive");
+%! fail("erlangb(-1, 1)", "positive");
+%! fail("erlangb(1, 0)", "positive");
+%! fail("erlangb(0, 1)", "positive");
+%! fail("erlangb('foo',1)", "positive");
+%! fail("erlangb(1,'bar')", "positive");
+%! fail("erlangb([1 1],[1 1 1])","common size");
+

--- a/main/queueing/inst/erlangb.m	Wed Jan 15 22:03:11 2014 +0000
+++ b/main/queueing/inst/erlangb.m	Fri Jan 17 09:41:11 2014 +0000
@@ -22,8 +22,9 @@
 ## @cindex Erlang-B formula
 ##
 ## Compute the value of the Erlang-B formula @math{E_B(A, m)} giving the
-## probability that a request arriving to a pool of @math{m} identical
-## servers is rejected. 
+## probability that an open system with @math{m} identical servers,
+## arrival rate @math{\lambda}, individual service rate @math{\mu}
+## and offered load @math{A = \lambda / \mu} has all servers busy.
 ## 
 ## @iftex
 ##
@@ -90,13 +91,15 @@
   B = arrayfun( @__erlangb_compute, A, m);
 endfunction
 
-## Compute E_B(A,m) recursively, as described in Guoping Zeng (June
-## 2003), Two common properties of the erlang-B function, erlang-C
-## function, and Engset blocking function, Elsevier Science, retrieved
-## 2011-02-03
+## Compute E_B(A,m) recursively, as described in:
+##
+## Guoping Zeng, Two common properties of the erlang-B function,
+## erlang-C function, and Engset blocking function, Mathematical and
+## Computer Modelling Volume 37, Issues 12–13, June 2003, Pages
+## 1287–1296 http://dx.doi.org/10.1016/S0895-7177(03)90040-9
 ##
 ## To improve numerical stability, the recursion is based on the inverse
-## of E_B rather than E_B itself.
+## 1 / E_B rather than E_B itself.
 function B = __erlangb_compute(A, m)
   Binv = 1.0;
   for k=1:m

--- a/main/queueing/inst/erlangc.m	Wed Jan 15 22:03:11 2014 +0000
+++ b/main/queueing/inst/erlangc.m	Fri Jan 17 09:41:11 2014 +0000
@@ -21,10 +21,11 @@
 ##
 ## @cindex Erlang-C formula
 ##
-## Compute the steady-state probability that a request finds all
-## @math{m} servers musy in a @math{M/M/m/\infty} queueing system with
-## @math{m} identical servers and infinite queue. This quantity is
-## also called Erlang-C formula @math{E_C(A, m)}.
+## Compute the steady-state probability that an open queueing system
+## with @math{m} identical servers, infinite wating space, arrival rate
+## @math{\lambda}, individual service rate @math{\mu} and offered load
+## @math{A = \lambda / \mu} are busy. This probability is also called
+## Erlang-C formula @math{E_C(A, m)}.
 ## 
 ## @iftex
 ##

--- a/main/queueing/inst/qsmmm.m	Wed Jan 15 22:03:11 2014 +0000
+++ b/main/queueing/inst/qsmmm.m	Fri Jan 17 09:41:11 2014 +0000
@@ -24,7 +24,7 @@
 ##
 ## Compute utilization, response time, average number of requests in
 ## service and throughput for a @math{M/M/m} queue, a queueing system
-## with @math{m} identical service centers connected to a single FCFS
+## with @math{m} identical servers connected to a single FCFS
 ## queue.
 ##
 ## @iftex