[churesult,a,b,x]=chu(a,b,x);

function [churesult,a,b,x]=chu(a,b,x);
churesult=[];
persistent a0 aintb alnx b0 beps c0 chu 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(a0), a0=0; end;
if isempty(aintb), aintb=0; end;
if isempty(alnx), alnx=0; end;
if isempty(b0), b0=0; end;
if isempty(beps), beps=0; end;
if isempty(c0), c0=0; end;
if isempty(churesult), churesult=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 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  CHU
%***PURPOSE  Compute the logarithmic confluent hypergeometric function.
%***LIBRARY   SLATEC (FNLIB)
%***CATEGORY  C11
%***TYPE      SINGLE PRECISION (CHU-S, DCHU-D)
%***KEYWORDS  FNLIB, LOGARITHMIC CONFLUENT HYPERGEOMETRIC FUNCTION,
%             SPECIAL FUNCTIONS
%***AUTHOR  Fullerton, W., (LANL)
%***DESCRIPTION
%
% CHU computes the logarithmic confluent hypergeometric function,
% U(A,B,X).
%
% Input Parameters:
%       A   real
%       B   real
%       X   real and positive
%
% This routine is not valid when 1+A-B is close to zero if X is small.
%
%***REFERENCES  (NONE)
%***ROUTINES CALLED  EXPREL, GAMMA, GAMR, POCH, POCH1, R1MACH, R9CHU,
%                    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  CHU
if firstCall,   pi=[3.14159265358979324e0];  end;
if firstCall,   eps=[0.0];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  CHU
if( eps==0.0 )
[ eps ]=r1mach(3);
end;
%
if( x==0.0 )
xermsg('SLATEC','CHU','X IS ZERO SO CHU IS INFINITE',1,2);
end;
if( x<0.0 )
xermsg('SLATEC','CHU','X IS NEGATIVE, USE CCHU',2,2);
end;
%
if( max(abs(a),1.0).*max(abs(1.0+a-b),1.0)>=0.99.*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.0+a-b)<sqrt(eps) )
xermsg('SLATEC','CHU','ALGORITHM IS BAD WHEN 1+A-B IS NEAR ZERO FOR SMALL X',10,2);
end;
%
aintb = fix(b+0.5);
if( b<0.0 )
aintb = fix(b-0.5);
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.0;
m = fix(n - 2);
if( m>=0 )
t = 1.0;
summlv = 1.0;
if( m~=0 )
%
for i = 1 : m;
xi = i;
t = t.*(a-b+xi).*x./((1.0-b+xi).*xi);
summlv = summlv + t;
end; i = fix(m+1);
end;
%
summlv = gamma(b-1.0).*gamr(a).*x.^(1-n).*xtoeps.*summlv;
end;
else;
%
% CONSIDER THE CASE B .LT. 1.0 FIRST.
%
summlv = 1.0;
if( n~=0 )
%
t = 1.0;
m = fix(-n);
for i = 1 : m;
xi1 = i - 1;
t = t.*(a+xi1).*x./((b+xi1).*(xi1+1.0));
summlv = summlv + t;
end; i = fix(m+1);
end;
%
summlv = poch(1.0+a-b,-a).*summlv;
end;
%
% NOW EVALUATE THE INFINITE SUM.     -----------------------------------
%
istrt = 0;
if( n<1 )
istrt = fix(1 - n);
end;
xi = istrt;
%
factor =(-1.0).^n.*gamr(1.0+a-b).*x.^istrt;
if( beps~=0.0 )
factor = factor.*beps.*pi./sin(beps.*pi);
end;
%
[pochai ,a,xi]=poch(a,xi);
[gamri1 ]=gamr(xi+1.0);
[gamrni ]=gamr(aintb+xi);
b0 = factor.*poch(a,xi-beps).*gamrni.*gamr(xi+1.0-beps);
%
if( abs(xtoeps-1.0)<=0.5 )
%
% X**(-BEPS) IS CLOSE TO 1.0, SO WE MUST BE CAREFUL IN EVALUATING
% THE DIFFERENCES
%
[pch1ai ]=poch1(a+xi,-beps);
[pch1i ,dumvar2,beps]=poch1(xi+1.0-beps,beps);
c0 = factor.*pochai.*gamrni.*gamri1.*(-poch1(b+xi,-beps)+pch1ai-pch1i+beps.*pch1ai.*pch1i);
%
% XEPS1 = (1.0 - X**(-BEPS)) / BEPS
xeps1 = alnx.*exprel(-beps.*alnx);
%
churesult = 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.0).*(xn+2..*xi-1.0)+xi.*(xi-beps)).*b0./(xi.*(b+xi1).*(a+xi1-beps));
t = c0 + xeps1.*b0;
churesult = churesult + t;
if( abs(t)<eps.*abs(churesult) )
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','CHU','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.*gamr(b+xi).*gamri1./beps;
b0 = xtoeps.*b0./beps;
%
churesult = 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;
churesult = churesult + t;
if( abs(t)<eps.*abs(churesult) )
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','CHU','NO CONVERGENCE IN 1000 TERMS OF THE ASCENDING SERIES',3,2);
end;
%
% use LUKE-S RATIONAL APPROX IN THE ASYMPTOTIC REGION.
%
churesult = x.^(-a).*r9chu(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 CINVIT