function p = ncx2cdf(x,v,delta)
%NCX2CDF Noncentral chi-square cumulative distribution function (cdf).
%   P = NCX2CDF(X,V,DELTA) Returns the noncentral chi-square cdf with V 
%   degrees of freedom and noncentrality parameter, DELTA, at the values 
%   in X.
%
%   The size of P is the common size of the input arguments. A scalar input  
%   functions as a constant matrix of the same size as the other inputs.     
%
%   Some texts refer to this distribution as the generalized Rayleigh,
%   Rayleigh-Rice, or Rice distribution.

%   Reference:
%      [1]  Evans, Merran, Hastings, Nicholas and Peacock, Brian,
%      "Statistical Distributions, Second Edition", Wiley
%      1993 p. 50-52.

%   B.A. Jones 4-29-94, ZP. You, 1-29-99
%   Copyright 1993-2000 The MathWorks, Inc. 
%   $Revision: 2.11 $  $Date: 2000/05/26 18:53:17 $

if nargin <  3, 
    error('Requires three input arguments.'); 
end

[errorcode x v delta] = distchck(3,x,v,delta);

if errorcode > 0
    error('Requires non-scalar arguments to match in size.');
end

% Initialize P to zero.
p = zeros(size(x));

p(x <= 0) = 0;  
p(delta < 0) = NaN; % can't have negative non-centrality parameter.
k = find((x > 0) & (delta >=0));
k1 = k;
delta = delta(k);
hdelta = delta/2;  % this is used in the loop, per-compute.
v = v(k);
x = x(k);
x1 = x;  % make a copy to be used when reverse direction.
v1 = v;

% when non-centrality parameter is very large, the initial values of
% the poisson numbers used in the approximation are very small,
% smaller than epsilon. This would cause premature convergence. To
% avoid that, we start from counter=hdelta, which is the peak of the
% poisson numbers, and go in both directions.

counter = floor(hdelta);
direction = 1;

% Sum the series.
while any(counter >= 0)
   pplus  = poisspdf(counter,hdelta).*chi2cdf(x,v+2*counter);
   j = isnan(pplus);
   if any(j) % if there's nan in pplus, we need to skip them and update indeces.
      if all(j), break; end
      j = ~j;
      x = x(j);
      v = v(j);
      hdelta = hdelta(j);
      counter = counter(j);
      k = k(j);
      pplus = pplus(k);
   end
   p(k) = p(k) + pplus; % accumulate p for k indices
   j = find((pplus > p(k)*eps) & ((counter > 0) | (direction == 1))); % check for exit
   if isempty(j) % no more computation needed in this direction
      if direction == 1 % if we haven't done the negative direction, reverse it.
         direction = -1;
         counter = floor(delta/2); % set counter back to the peak;
         j = find(counter>0);  % get rid of all the peaks that happens to be 0.
         if isempty(j), break; end % if nothing left, exit.
         x = x1(j);  % only keep the ones that doesn't have a peak at 0.
         v = v1(j);
         delta = delta(j);
         hdelta = delta/2;
         k = k1(j);
         counter = counter(j);
      else
         break;
      end
   else
      x = x(j);
      v = v(j);
      hdelta = hdelta(j);
      counter = counter(j);
      k = k(j);
   end
   counter = counter + direction;
end
