rho=findrhoc(n,p)

%
% rho=findrhoc(n,p)
%
% Finds an intensity rho such that C(n,rho)=p.  
%
% Note: Must have 0<p<1.  Returns NaN if p is not in this range.
%
function rho=findrhoc(n,p)
%
% Sanity check- make sure that n is a positive integer.
%
  if ((floor(n) ~= n) || (n < 1))
    warning('n is not a positive integer');
    rho=NaN;
    return;
  end;
%
% Sanity check- make sure that p is a probability with 0<p<1.
%
  if ((p<0.0) || (p>1.0))
     warning('Invalid p value!');
     rho=NaN;
     return;
  end;
%
% We know that at rho=0, p=0, and at rho=+Inf, p=1.   We start by finding
% an interval [0,a] containing the root.
%
a=1.0;
testp=erlangc(n,a);
while (testp < p),
  a=a*2.0;
  testp=erlangc(n,a);
end;
%
% Now, the root is somewhere between 0 and a.  Use bisection to find it.
% 
%
left=0.0;
right=a;
mid=(left+right)/2;
midp=erlangc(n,mid);
while ((right-left) > 0.0001*max([1 left])),
  if (midp < p),
    left=mid;
    mid=(left+right)/2;
    midp=erlangc(n,mid);     
  else
    right=mid;
    mid=(left+right)/2;
    midp=erlangc(n,mid);     
  end;
end;
%
% Return the left end point of the current interval, which has prob < p.
%
rho=left;