function D = debye(x,k)%DEBYE Debye function.%   Y = DEBYE(X,K) returns the K-th order Debye function, evaluated at X.%   X and K are scalars.  For real X, this is defined as%%      (K/X^K) * integral from 0 to X of (t/(exp(t)-1)) dt%%   The Debye function satisfies D_K(-X) = D_K(X) + K*X/(K+1) however, DEBYE%   calculates D_K(X) directly for negative X.%%   Accurate to at least 6 digits (absolute) for orders 1 through 4.%   Copyright 2005 The MathWorks, Inc.%   $Revision: 1.1.6.1 $  $Date: 2005/11/18 14:28:43 $absx = abs(x);if absx > 1    D = k*quadl(@debye_integrand,0,x,1e-10,[],k)./(x.^k);elseif absx > 0    % see Abramowitz&Stegun, 27.1.1    j = 1:10;    logB2j = log([1/6 1/30 1/42 1/30 5/66 691/2730 7/6 3617/510 43867/798 174611/330]);    signB2j = [1 -1 1 -1 1 -1 1 -1 1 -1];    D = 1 - .5*x*k./(k+1) ...          + k*sum(signB2j.*exp(logB2j + 2*j.*log(absx) - log(2*j+k) - gammaln(2*j+1)));else    D = ones(class(x));endfunction y = debye_integrand(t,k)y = ones(size(t),class(t));nz = (abs(t) >= realmin(class(t)));y(nz) = (t(nz).^k) ./ expm1(t(nz));