[d9churesult,a,b,z]=d9chu(a,b,z);

function [d9churesult,a,b,z]=d9chu(a,b,z);
d9churesult=[];
persistent aa ab anbn bb bp c2 ct1 ct2 ct3 d1z eps first firstCall g1 g2 g3 i j sab sqeps x2i1 ; if isempty(firstCall),firstCall=1;end; 

;
if isempty(i), i=0; end;
if isempty(j), j=0; end;
%***BEGIN PROLOGUE  D9CHU
%***SUBSIDIARY
%***PURPOSE  Evaluate for large Z  Z**A * U(A,B,Z) where U is the
%            logarithmic confluent hypergeometric function.
%***LIBRARY   SLATEC (FNLIB)
%***CATEGORY  C11
%***TYPE      doubleprecision (R9CHU-S, D9CHU-D)
%***KEYWORDS  FNLIB, LOGARITHMIC CONFLUENT HYPERGEOMETRIC FUNCTION,
%             SPECIAL FUNCTIONS
%***AUTHOR  Fullerton, W., (LANL)
%***DESCRIPTION
%
% Evaluate for large Z  Z**A * U(A,B,Z)  where U is the logarithmic
% confluent hypergeometric function.  A rational approximation due to Y.
% L. Luke is used.  When U is not in the asymptotic region, i.e., when A
% or B is large compared with Z, considerable significance loss occurs.
% A warning is provided when the computed result is less than half
% precision.
%
%***REFERENCES  (NONE)
%***ROUTINES CALLED  D1MACH, 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)
%   900720  Routine changed from user-callable to subsidiary.  (WRB)
%***end PROLOGUE  D9CHU
if isempty(aa), aa=zeros(1,4); end;
if isempty(bb), bb=zeros(1,4); end;
if isempty(ab), ab=0; end;
if isempty(anbn), anbn=0; end;
if isempty(bp), bp=0; end;
if isempty(ct1), ct1=0; end;
if isempty(ct2), ct2=0; end;
if isempty(ct3), ct3=0; end;
if isempty(c2), c2=0; end;
if isempty(d1z), d1z=0; end;
if isempty(eps), eps=0; end;
if isempty(g1), g1=0; end;
if isempty(g2), g2=0; end;
if isempty(g3), g3=0; end;
if isempty(sab), sab=0; end;
if isempty(sqeps), sqeps=0; end;
if isempty(x2i1), x2i1=0; end;
if isempty(first), first=false; end;
if firstCall,   first=[true];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  D9CHU
if( first )
eps = 4.0d0.*d1mach(4);
sqeps = sqrt(d1mach(4));
end;
first = false;
%
bp = 1.0d0 + a - b;
ab = a.*bp;
ct2 = 2.0d0.*(z-ab);
sab = a + bp;
%
bb(1) = 1.0d0;
aa(1) = 1.0d0;
%
ct3 = sab + 1.0d0 + ab;
bb(2) = 1.0d0 + 2.0d0.*z./ct3;
aa(2) = 1.0d0 + ct2./ct3;
%
anbn = ct3 + sab + 3.0d0;
ct1 = 1.0d0 + 2.0d0.*z./anbn;
bb(3) = 1.0d0 + 6.0d0.*ct1.*z./ct3;
aa(3) = 1.0d0 + 6.0d0.*ab./anbn + 3.0d0.*ct1.*ct2./ct3;
%
for i = 4 : 300;
x2i1 = 2.*i - 3;
ct1 = x2i1./(x2i1-2.0d0);
anbn = anbn + x2i1 + sab;
ct2 =(x2i1-1.0d0)./anbn;
c2 = x2i1.*ct2 - 1.0d0;
d1z = x2i1.*2.0d0.*z./anbn;
%
ct3 = sab.*ct2;
g1 = d1z + ct1.*(c2+ct3);
g2 = d1z - c2;
g3 = ct1.*(1.0d0-ct3-2.0d0.*ct2);
%
bb(4) = g1.*bb(3) + g2.*bb(2) + g3.*bb(1);
aa(4) = g1.*aa(3) + g2.*aa(2) + g3.*aa(1);
if( abs(aa(4).*bb(1)-aa(1).*bb(4))<eps.*abs(bb(4).*bb(1)) )
d9churesult = aa(4)./bb(4);
if( d9churesult<sqeps || d9churesult>1.0D0./sqeps )
xermsg('SLATEC','D9CHU','ANSWER LT HALF PRECISION',2,1);
end;
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(3)), assignin('caller','FUntemp',z); 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;
%
% IF OVERFLOWS OR UNDERFLOWS PROVE TO BE A PROBLEM, THE STATEMENTS
% BELOW COULD BE ALTERED TO INCORPORATE A DYNAMICALLY ADJUSTED SCALE
% FACTOR.
%
for j = 1 : 3;
aa(j) = aa(j+1);
bb(j) = bb(j+1);
end; j = fix(3+1);
end; i = fix(300+1);
xermsg('SLATEC','D9CHU','NO CONVERGENCE IN 300 TERMS',2,2);
%
d9churesult = aa(4)./bb(4);
%
if( d9churesult<sqeps || d9churesult>1.0D0./sqeps )
xermsg('SLATEC','D9CHU','ANSWER LT HALF PRECISION',2,1);
end;
%
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(3)), assignin('caller','FUntemp',z); 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 D9GMIC