function [pd,pb,mq,p] = erlangc(A,N,Q)
%ERLANGC Calculates probabilities using the Erlang-C formula
%Call [pd,pb,mq,p] = erlangc(A,N,Q)
%Where A = offered traffic in Erlang
%      N = number of serving trunks(Should be <~1000 to avoid numerical overflow)
%      Q = number of queue places 
%Returns
%      p = vector of Erlang-B probabilities
%      pd = probability of a call being delayed
%      pb = blocking probability 
%      mq = mean queue length for calls that are queued


u = exp((0:N)*log(A) - [0, cumsum(log(1:N))]);
u = [u, [exp((N + (1:Q))*log(A)-ones(1,Q).*sum(log(1:N)) - (1:Q).*log(N))]];
s =sum(u);
p =u/s;

pd = sum(p(N+2:N+1+Q));
mq = sum([1:Q].*p(N+2:N+1+Q));  %division bu pd removed???
pb = p(N+1+Q);



