[zr,zi,fnu,n,cyr,cyi,tol]=zrati(zr,zi,fnu,n,cyr,cyi,tol);

function [zr,zi,fnu,n,cyr,cyi,tol]=zrati(zr,zi,fnu,n,cyr,cyi,tol);
%***BEGIN PROLOGUE  ZRATI
%***SUBSIDIARY
%***PURPOSE  Subsidiary to ZBESH, ZBESI and ZBESK
%***LIBRARY   SLATEC
%***TYPE      ALL (CRATI-A, ZRATI-A)
%***AUTHOR  Amos, D. E., (SNL)
%***DESCRIPTION
%
%     ZRATI COMPUTES RATIOS OF I BESSEL FUNCTIONS BY BACKWARD
%     RECURRENCE.  THE STARTING INDEX IS DETERMINED BY FORWARD
%     RECURRENCE AS DESCRIBED IN J. RES. OF NAT. BUR. OF STANDARDS-B,
%     MATHEMATICAL SCIENCES, VOL 77B, P111-114, SEPTEMBER, 1973,
%     BESSEL FUNCTIONS I AND J OF COMPLEX ARGUMENT AND INTEGER ORDER,
%     BY D. J. SOOKNE.
%
%***SEE ALSO  ZBESH, ZBESI, ZBESK
%***ROUTINES CALLED  ZABS, ZDIV
%***REVISION HISTORY  (YYMMDD)
%   830501  DATE WRITTEN
%   910415  Prologue converted to Version 4.0 format.  (BAB)
%***end PROLOGUE  ZRATI
persistent ak amagz ap1 ap2 arg az cdfnui cdfnur conei coner czeroi czeror dfnu fdnu firstCall flam fnup i id idnu inu itime k kk magz p1i p1r p2i p2r pti ptr rak rap1 rho rt2 rzi rzr t1i t1r test test1 tti ttr ; if isempty(firstCall),firstCall=1;end; 

if isempty(ak), ak=0; end;
if isempty(amagz), amagz=0; end;
if isempty(ap1), ap1=0; end;
if isempty(ap2), ap2=0; end;
if isempty(arg), arg=0; end;
if isempty(az), az=0; end;
if isempty(cdfnui), cdfnui=0; end;
if isempty(cdfnur), cdfnur=0; end;
if isempty(conei), conei=0; end;
if isempty(coner), coner=0; end;
if isempty(czeroi), czeroi=0; end;
if isempty(czeror), czeror=0; end;
if isempty(dfnu), dfnu=0; end;
if isempty(fdnu), fdnu=0; end;
if isempty(flam), flam=0; end;
if isempty(fnup), fnup=0; end;
if isempty(pti), pti=0; end;
if isempty(ptr), ptr=0; end;
if isempty(p1i), p1i=0; end;
if isempty(p1r), p1r=0; end;
if isempty(p2i), p2i=0; end;
if isempty(p2r), p2r=0; end;
if isempty(rak), rak=0; end;
if isempty(rap1), rap1=0; end;
if isempty(rho), rho=0; end;
if isempty(rt2), rt2=0; end;
if isempty(rzi), rzi=0; end;
if isempty(rzr), rzr=0; end;
if isempty(test), test=0; end;
if isempty(test1), test1=0; end;
if isempty(tti), tti=0; end;
if isempty(ttr), ttr=0; end;
if isempty(t1i), t1i=0; end;
if isempty(t1r), t1r=0; end;
if isempty(i), i=0; end;
if isempty(id), id=0; end;
if isempty(idnu), idnu=0; end;
if isempty(inu), inu=0; end;
if isempty(itime), itime=0; end;
if isempty(k), k=0; end;
if isempty(kk), kk=0; end;
if isempty(magz), magz=0; end;
if firstCall,   czeror =[0.0d0];  end;
if firstCall,  czeroi =[0.0d0];  end;
if firstCall,  coner =[1.0d0];  end;
if firstCall,  conei =[0.0d0];  end;
if firstCall,  rt2=[1.41421356237309505d0];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  ZRATI
[az ,zr,zi]=zabs(zr,zi);
inu = fix(fnu);
idnu = fix(inu + n - 1);
magz = fix(az);
amagz = magz + 1;
fdnu = idnu;
fnup = max(amagz,fdnu);
id = fix(idnu - magz - 1);
itime = 1;
k = 1;
ptr = 1.0d0./az;
rzr = ptr.*(zr+zr).*ptr;
rzi = -ptr.*(zi+zi).*ptr;
t1r = rzr.*fnup;
t1i = rzi.*fnup;
p2r = -t1r;
p2i = -t1i;
p1r = coner;
p1i = conei;
t1r = t1r + rzr;
t1i = t1i + rzi;
if( id>0 )
id = 0;
end;
[ap2 ,p2r,p2i]=zabs(p2r,p2i);
[ap1 ,p1r,p1i]=zabs(p1r,p1i);
%-----------------------------------------------------------------------
%     THE OVERFLOW TEST ON K(FNU+I-1,Z) BEFORE THE CALL TO CBKNU
%     GUARANTEES THAT P2 IS ON SCALE. SCALE TEST1 AND ALL SUBSEQUENT
%     P2 VALUES BY AP1 TO ENSURE THAT AN OVERFLOW DOES NOT OCCUR
%     PREMATURELY.
%-----------------------------------------------------------------------
arg =(ap2+ap2)./(ap1.*tol);
test1 = sqrt(arg);
test = test1;
rap1 = 1.0d0./ap1;
p1r = p1r.*rap1;
p1i = p1i.*rap1;
p2r = p2r.*rap1;
p2i = p2i.*rap1;
ap2 = ap2.*rap1;
while( true );
k = fix(k + 1);
ap1 = ap2;
ptr = p2r;
pti = p2i;
p2r = p1r -(t1r.*ptr-t1i.*pti);
p2i = p1i -(t1r.*pti+t1i.*ptr);
p1r = ptr;
p1i = pti;
t1r = t1r + rzr;
t1i = t1i + rzi;
[ap2 ,p2r,p2i]=zabs(p2r,p2i);
if( ap1>test )
if( itime==2 )
break;
end;
ak = zabs(t1r,t1i).*0.5d0;
flam = ak + sqrt(ak.*ak-1.0d0);
rho = min(ap2./ap1,flam);
test = test1.*sqrt(rho./(rho.*rho-1.0d0));
itime = 2;
end;
end;
kk = fix(k + 1 - id);
ak = kk;
t1r = ak;
t1i = czeroi;
dfnu = fnu +(n-1);
p1r = 1.0d0./ap2;
p1i = czeroi;
p2r = czeror;
p2i = czeroi;
for i = 1 : kk;
ptr = p1r;
pti = p1i;
rap1 = dfnu + t1r;
ttr = rzr.*rap1;
tti = rzi.*rap1;
p1r =(ptr.*ttr-pti.*tti) + p2r;
p1i =(ptr.*tti+pti.*ttr) + p2i;
p2r = ptr;
p2i = pti;
t1r = t1r - coner;
end; i = fix(kk+1);
if( p1r==czeror && p1i==czeroi )
p1r = tol;
p1i = tol;
end;
[p2r,p2i,p1r,p1i,cyr(n),cyi(n)]=zdiv(p2r,p2i,p1r,p1i,cyr(n),cyi(n));
if( n==1 )
return;
end;
k = fix(n - 1);
ak = k;
t1r = ak;
t1i = czeroi;
cdfnur = fnu.*rzr;
cdfnui = fnu.*rzi;
for i = 2 : n;
ptr = cdfnur +(t1r.*rzr-t1i.*rzi) + cyr(k+1);
pti = cdfnui +(t1r.*rzi+t1i.*rzr) + cyi(k+1);
[ak ,ptr,pti]=zabs(ptr,pti);
if( ak==czeror )
ptr = tol;
pti = tol;
ak = tol.*rt2;
end;
rak = coner./ak;
cyr(k) = rak.*ptr.*rak;
cyi(k) = -rak.*pti.*rak;
t1r = t1r - coner;
k = fix(k - 1);
end; i = fix(n+1);
end %subroutine zrati
%DECK ZS1S2