function y = nctpdf(x,v,delta)
%NCTPDF Noncentral T probability density function (pdf).
%   Y = NCTPDF(X,V,DELTA) Returns the noncentral T pdf with V degrees 
%   of freedom and noncentrality parameter, DELTA, at the values in X. 
%
%   The size of Y is the common size of the input arguments. A scalar input  
%   functions as a constant matrix of the same size as the other inputs.     

%   Reference:
%      [1]  Johnson, Norman, and Kotz, Samuel, "Distributions in
%      Statistics: Continuous Univariate Distributions-2", Wiley
%      1970 pp. 204-205 equation 8.
%      [2]  Evans, Merran, Hastings, Nicholas and Peacock, Brian,
%      "Statistical Distributions, Second Edition", Wiley
%      1993 pp. 147-148.

%   B.A. Jones 4-21-95
%   Copyright 1993-2000 The MathWorks, Inc. 
%   $Revision: 2.14 $  $Date: 2000/05/26 18:53:08 $


if nargin < 1, 
    error('Requires at least one input argument.');
end

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

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

% if some delta==0, call tcdf for those entries, call nctcdf for other entries.
f = (delta == 0);
if any(f(:))
   y = zeros(size(delta));
   y(f) = tpdf(x(f),v(f));
   f1 = ~f;
   if any(f1)
      y(f1) = nctpdf(x(f1),v(f1),delta(f1));
   end
   return
end

%   Initialize Y to zero.
y = zeros(size(x));

k = find(v > 0);
if any(k)
   nu = v(k);
   d = delta(k);
   t = x(k);
   lnu = log(nu+t.^2);
   lnu1 = gammaln((nu+1)/2);
   lnu2 = gammaln((nu+2)/2);
   term = exp(-0.5.*(log(pi)+log(nu)) -(d.^2)/2  + lnu1 - gammaln(nu/2) ...
      + ((nu+1)/2).*(log(nu)-lnu));
   
   infsum = zeros(size(k));
   qeps = eps^(3/4);
   tsq = t.^2;
   ratio = tsq .* (d.^2) ./ (nu + tsq);

   % Set up for infinite sum by computing first two terms.
   j = 0;
   t0 = 1;
   t1 = t .* d .* sqrt(2) .* exp(lnu2 - lnu1) ./ sqrt(nu + tsq);
   infsum = infsum + term .* (t0 + t1);
   
   % Sum the series by computing term j from term j-2
   while 1
      t0 = t0 .* ratio .* (nu+j+1) / ((j+1) * (j+2));
      j = j + 1;
      t1 = t1 .* ratio .* (nu+j+1) / ((j+1) * (j+2));
      j = j + 1;
      twoterms = term .* (t0 + t1);
      infsum = infsum + twoterms;

      % Convergence test.
      if all(abs(twoterms)./infsum < qeps), break; end
   end
   y(k) = infsum;
end

p(v <= 0) = NaN;
