[dchuresult,a,b,x]=dchu(a,b,x);

function [dchuresult,a,b,x]=dchu(a,b,x);
dchuresult=[];
persistent a0 aintb alnx b0 beps c0 eps factor firstCall gamri1 gamrni i istrt m n pch1ai pch1i pi pochai summlv t xeps1 xi xi1 xn xtoeps ; if isempty(firstCall),firstCall=1;end; 

;
if isempty(i), i=0; end;
if isempty(istrt), istrt=0; end;
if isempty(m), m=0; end;
if isempty(n), n=0; end;
%***BEGIN PROLOGUE  DCHU
%***PURPOSE  Compute the logarithmic confluent hypergeometric function.
%***LIBRARY   SLATEC (FNLIB)
%***CATEGORY  C11
%***TYPE      doubleprecision (CHU-S, DCHU-D)
%***KEYWORDS  FNLIB, LOGARITHMIC CONFLUENT HYPERGEOMETRIC FUNCTION,
%             SPECIAL FUNCTIONS
%***AUTHOR  Fullerton, W., (LANL)
%***DESCRIPTION
%
% DCHU(A,B,X) calculates the doubleprecision logarithmic confluent
% hypergeometric function U(A,B,X) for doubleprecision arguments
% A, B, and X.
%
% This routine is not valid when 1+A-B is close to zero if X is small.
%
%***REFERENCES  (NONE)
%***ROUTINES CALLED  D1MACH, D9CHU, DEXPRL, DGAMMA, DGAMR, DPOCH,
%                    DPOCH1, XERMSG
%***REVISION HISTORY  (YYMMDD)
%   770801  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890531  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
%   900727  Added EXTERNAL statement.  (WRB)
%***end PROLOGUE  DCHU
if isempty(aintb), aintb=0; end;
if isempty(alnx), alnx=0; end;
if isempty(a0), a0=0; end;
if isempty(beps), beps=0; end;
if isempty(b0), b0=0; end;
if isempty(c0), c0=0; end;
if isempty(eps), eps=0; end;
if isempty(factor), factor=0; end;
if isempty(gamri1), gamri1=0; end;
if isempty(gamrni), gamrni=0; end;
if isempty(pch1ai), pch1ai=0; end;
if isempty(pch1i), pch1i=0; end;
if isempty(pi), pi=0; end;
if isempty(pochai), pochai=0; end;
if isempty(summlv), summlv=0; end;
if isempty(t), t=0; end;
if isempty(xeps1), xeps1=0; end;
if isempty(xi), xi=0; end;
if isempty(xi1), xi1=0; end;
if isempty(xn), xn=0; end;
if isempty(xtoeps), xtoeps=0; end;
if firstCall,   pi=[3.141592653589793238462643383279503d0];  end;
if firstCall,   eps=[0.0d0];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  DCHU
if( eps==0.0d0 )
[ eps ]=d1mach(3);
end;
%
if( x==0.0D0 )
xermsg('SLATEC','DCHU','X IS ZERO SO DCHU IS INFINITE',1,2);
end;
if( x<0.0D0 )
xermsg('SLATEC','DCHU','X IS NEGATIVE, USE CCHU',2,2);
end;
%
if( max(abs(a),1.0d0).*max(abs(1.0d0+a-b),1.0d0)>=0.99d0.*abs(x) )
%
% THE ASCENDING SERIES WILL BE USED, BECAUSE THE DESCENDING RATIONAL
% APPROXIMATION (WHICH IS BASED ON THE ASYMPTOTIC SERIES) IS UNSTABLE.
%
if( abs(1.0D0+a-b)<sqrt(eps) )
xermsg('SLATEC','DCHU','ALGORITHMIS BAD WHEN 1+A-B IS NEAR ZERO FOR SMALL X',10,2);
end;
%
if( b>=0.0d0 )
aintb = fix(b+0.5d0);
end;
if( b<0.0d0 )
aintb = fix(b-0.5d0);
end;
beps = b - aintb;
n = fix(aintb);
%
alnx = log(x);
xtoeps = exp(-beps.*alnx);
%
% EVALUATE THE FINITE SUM.     -----------------------------------------
%
if( n>=1 )
%
% NOW CONSIDER THE CASE B .GE. 1.0.
%
summlv = 0.0d0;
m = fix(n - 2);
if( m>=0 )
t = 1.0d0;
summlv = 1.0d0;
if( m~=0 )
%
for i = 1 : m;
xi = i;
t = t.*(a-b+xi).*x./((1.0d0-b+xi).*xi);
summlv = summlv + t;
end; i = fix(m+1);
end;
%
summlv = dgamma(b-1.0d0).*dgamr(a).*x.^(1-n).*xtoeps.*summlv;
end;
else;
%
% CONSIDER THE CASE B .LT. 1.0 FIRST.
%
summlv = 1.0d0;
if( n~=0 )
%
t = 1.0d0;
m = fix(-n);
for i = 1 : m;
xi1 = i - 1;
t = t.*(a+xi1).*x./((b+xi1).*(xi1+1.0d0));
summlv = summlv + t;
end; i = fix(m+1);
end;
%
summlv = dpoch(1.0d0+a-b,-a).*summlv;
end;
%
% NEXT EVALUATE THE INFINITE SUM.     ----------------------------------
%
istrt = 0;
if( n<1 )
istrt = fix(1 - n);
end;
xi = istrt;
%
factor =(-1.0d0).^n.*dgamr(1.0d0+a-b).*x.^istrt;
if( beps~=0.0d0 )
factor = factor.*beps.*pi./sin(beps.*pi);
end;
%
[pochai ,a,xi]=dpoch(a,xi);
[gamri1 ]=dgamr(xi+1.0d0);
[gamrni ]=dgamr(aintb+xi);
b0 = factor.*dpoch(a,xi-beps).*gamrni.*dgamr(xi+1.0d0-beps);
%
if( abs(xtoeps-1.0d0)<=0.5d0 )
%
% X**(-BEPS) IS CLOSE TO 1.0D0, SO WE MUST BE CAREFUL IN EVALUATING THE
% DIFFERENCES.
%
[pch1ai ]=dpoch1(a+xi,-beps);
[pch1i ,dumvar2,beps]=dpoch1(xi+1.0d0-beps,beps);
c0 = factor.*pochai.*gamrni.*gamri1.*(-dpoch1(b+xi,-beps)+pch1ai-pch1i+beps.*pch1ai.*pch1i);
%
% XEPS1 = (1.0 - X**(-BEPS))/BEPS = (X**(-BEPS) - 1.0)/(-BEPS)
xeps1 = alnx.*dexprl(-beps.*alnx);
%
dchuresult = summlv + c0 + xeps1.*b0;
xn = n;
for i = 1 : 1000;
xi = istrt + i;
xi1 = istrt + i - 1;
b0 =(a+xi1-beps).*b0.*x./((xn+xi1).*(xi-beps));
c0 =(a+xi1).*c0.*x./((b+xi1).*xi)-((a-1.0d0).*(xn+2.0d0.*xi-1.0d0)+xi.*(xi-beps)).*b0./(xi.*(b+xi1).*(a+xi1-beps));
t = c0 + xeps1.*b0;
dchuresult = dchuresult + t;
if( abs(t)<eps.*abs(dchuresult) )
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(3)), assignin('caller','FUntemp',x); evalin('caller',[inputname(3),'=FUntemp;']); end
if csnil&&~isempty(inputname(2)), assignin('caller','FUntemp',b); evalin('caller',[inputname(2),'=FUntemp;']); end
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',a); evalin('caller',[inputname(1),'=FUntemp;']); end
return;
end;
end; i = fix(1000+1);
xermsg('SLATEC','DCHU','NO CONVERGENCE IN 1000 TERMS OF THE ASCENDING SERIES',3,2);
end;
%
% X**(-BEPS) IS VERY DIFFERENT FROM 1.0, SO THE STRAIGHTFORWARD
% FORMULATION IS STABLE.
%
a0 = factor.*pochai.*dgamr(b+xi).*gamri1./beps;
b0 = xtoeps.*b0./beps;
%
dchuresult = summlv + a0 - b0;
for i = 1 : 1000;
xi = istrt + i;
xi1 = istrt + i - 1;
a0 =(a+xi1).*a0.*x./((b+xi1).*xi);
b0 =(a+xi1-beps).*b0.*x./((aintb+xi1).*(xi-beps));
t = a0 - b0;
dchuresult = dchuresult + t;
if( abs(t)<eps.*abs(dchuresult) )
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(3)), assignin('caller','FUntemp',x); evalin('caller',[inputname(3),'=FUntemp;']); end
if csnil&&~isempty(inputname(2)), assignin('caller','FUntemp',b); evalin('caller',[inputname(2),'=FUntemp;']); end
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',a); evalin('caller',[inputname(1),'=FUntemp;']); end
return;
end;
end; i = fix(1000+1);
xermsg('SLATEC','DCHU','NO CONVERGENCE IN 1000 TERMS OF THE ASCENDING SERIES',3,2);
end;
%
% use LUKE-S RATIONAL APPROXIMATION IN THE ASYMPTOTIC REGION.
%
dchuresult = x.^(-a).*d9chu(a,b,x);
%
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(3)), assignin('caller','FUntemp',x); evalin('caller',[inputname(3),'=FUntemp;']); end
if csnil&&~isempty(inputname(2)), assignin('caller','FUntemp',b); evalin('caller',[inputname(2),'=FUntemp;']); end
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',a); evalin('caller',[inputname(1),'=FUntemp;']); end
end
%DECK DCKDER