[zr,zi,fnu,kode,n,yr,yi,nz,tol]=zmlri(zr,zi,fnu,kode,n,yr,yi,nz,tol);

function [zr,zi,fnu,kode,n,yr,yi,nz,tol]=zmlri(zr,zi,fnu,kode,n,yr,yi,nz,tol);
%***BEGIN PROLOGUE  ZMLRI
%***SUBSIDIARY
%***PURPOSE  Subsidiary to ZBESI and ZBESK
%***LIBRARY   SLATEC
%***TYPE      ALL (CMLRI-A, ZMLRI-A)
%***AUTHOR  Amos, D. E., (SNL)
%***DESCRIPTION
%
%     ZMLRI COMPUTES THE I BESSEL FUNCTION FOR RE(Z).GE.0.0 BY THE
%     MILLER ALGORITHM NORMALIZED BY A NEUMANN SERIES.
%
%***SEE ALSO  ZBESI, ZBESK
%***ROUTINES CALLED  D1MACH, DGAMLN, ZABS, ZEXP, ZLOG, ZMLT
%***REVISION HISTORY  (YYMMDD)
%   830501  DATE WRITTEN
%   910415  Prologue converted to Version 4.0 format.  (BAB)
%   930122  Added ZEXP and ZLOG to EXTERNAL statement.  (RWC)
%***end PROLOGUE  ZMLRI
%     COMPLEX CK,CNORM,CONE,CTWO,CZERO,PT,P1,P2,RZ,SUM,Y,Z
persistent ack ak ap at az bk cki ckr cnormi cnormr conei coner firstCall fkap fkk flam fnf i iaz idum ifnu igo inu itime k kk km m p1i p1r p2i p2r pti ptr raz rho rho2 rzi rzr scle sti str sumi sumr tfnf tst zeroi zeror ; if isempty(firstCall),firstCall=1;end; 

if isempty(ack), ack=0; end;
if isempty(ak), ak=0; end;
if isempty(ap), ap=0; end;
if isempty(at), at=0; end;
if isempty(az), az=0; end;
if isempty(bk), bk=0; end;
if isempty(cki), cki=0; end;
if isempty(ckr), ckr=0; end;
if isempty(cnormi), cnormi=0; end;
if isempty(cnormr), cnormr=0; end;
if isempty(conei), conei=0; end;
if isempty(coner), coner=0; end;
if isempty(fkap), fkap=0; end;
if isempty(fkk), fkk=0; end;
if isempty(flam), flam=0; end;
if isempty(fnf), fnf=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(raz), raz=0; end;
if isempty(rho), rho=0; end;
if isempty(rho2), rho2=0; end;
if isempty(rzi), rzi=0; end;
if isempty(rzr), rzr=0; end;
if isempty(scle), scle=0; end;
if isempty(sti), sti=0; end;
if isempty(str), str=0; end;
if isempty(sumi), sumi=0; end;
if isempty(sumr), sumr=0; end;
if isempty(tfnf), tfnf=0; end;
if isempty(tst), tst=0; end;
if isempty(zeroi), zeroi=0; end;
if isempty(zeror), zeror=0; end;
if isempty(i), i=0; end;
if isempty(iaz), iaz=0; end;
if isempty(idum), idum=0; end;
if isempty(ifnu), ifnu=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(km), km=0; end;
if isempty(m), m=0; end;
if isempty(igo), igo=0; end;
if firstCall,   zeror =[0.0d0];  end;
if firstCall,  zeroi =[0.0d0];  end;
if firstCall,  coner =[1.0d0];  end;
if firstCall,  conei=[0.0d0];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  ZMLRI
scle = d1mach(1)./tol;
nz = 0;
[az ,zr,zi]=zabs(zr,zi);
iaz = fix(az);
ifnu = fix(fnu);
inu = fix(ifnu + n - 1);
at = iaz + 1.0d0;
raz = 1.0d0./az;
str = zr.*raz;
sti = -zi.*raz;
ckr = str.*at.*raz;
cki = sti.*at.*raz;
rzr =(str+str).*raz;
rzi =(sti+sti).*raz;
p1r = zeror;
p1i = zeroi;
p2r = coner;
p2i = conei;
ack =(at+1.0d0).*raz;
rho = ack + sqrt(ack.*ack-1.0d0);
rho2 = rho.*rho;
tst =(rho2+rho2)./((rho2-1.0d0).*(rho-1.0d0));
tst = tst./tol;
%-----------------------------------------------------------------------
%     COMPUTE RELATIVE TRUNCATION ERROR INDEX FOR SERIES
%-----------------------------------------------------------------------
ak = at;
igo=0;
for i = 1 : 80;
ptr = p2r;
pti = p2i;
p2r = p1r -(ckr.*ptr-cki.*pti);
p2i = p1i -(cki.*ptr+ckr.*pti);
p1r = ptr;
p1i = pti;
ckr = ckr + rzr;
cki = cki + rzi;
[ap ,p2r,p2i]=zabs(p2r,p2i);
if( ap>tst.*ak.*ak )
k = 0;
if( inu<iaz )
igo=1;
break;
end;
%-----------------------------------------------------------------------
%     COMPUTE RELATIVE TRUNCATION ERROR FOR RATIOS
%-----------------------------------------------------------------------
p1r = zeror;
p1i = zeroi;
p2r = coner;
p2i = conei;
at = inu + 1.0d0;
str = zr.*raz;
sti = -zi.*raz;
ckr = str.*at.*raz;
cki = sti.*at.*raz;
ack = at.*raz;
tst = sqrt(ack./tol);
itime = 1;
for k = 1 : 80;
ptr = p2r;
pti = p2i;
p2r = p1r -(ckr.*ptr-cki.*pti);
p2i = p1i -(ckr.*pti+cki.*ptr);
p1r = ptr;
p1i = pti;
ckr = ckr + rzr;
cki = cki + rzi;
[ap ,p2r,p2i]=zabs(p2r,p2i);
if( ap>=tst )
if( itime==2 )
igo=1;
break;
end;
[ack ,ckr,cki]=zabs(ckr,cki);
flam = ack + sqrt(ack.*ack-1.0d0);
fkap = ap./zabs(p1r,p1i);
rho = min(flam,fkap);
tst = tst.*sqrt(rho./(rho.*rho-1.0d0));
itime = 2;
end;
end;
break;
else;
ak = ak + 1.0d0;
end;
end;
if(igo==0)
nz = -2;
return;
end;
%-----------------------------------------------------------------------
%     BACKWARD RECURRENCE AND SUM NORMALIZING RELATION
%-----------------------------------------------------------------------
k = fix(k + 1);
i = fix(i + 1);
kk = fix(max(i+iaz,k+inu));
fkk = kk;
p1r = zeror;
p1i = zeroi;
%-----------------------------------------------------------------------
%     SCALE P2 AND SUM BY SCLE
%-----------------------------------------------------------------------
p2r = scle;
p2i = zeroi;
fnf = fnu - ifnu;
tfnf = fnf + fnf;
bk = dgamln(fkk+tfnf+1.0d0,idum) - dgamln(fkk+1.0d0,idum)- dgamln(tfnf+1.0d0,idum);
bk = exp(bk);
sumr = zeror;
sumi = zeroi;
km = fix(kk - inu);
for i = 1 : km;
ptr = p2r;
pti = p2i;
p2r = p1r +(fkk+fnf).*(rzr.*ptr-rzi.*pti);
p2i = p1i +(fkk+fnf).*(rzi.*ptr+rzr.*pti);
p1r = ptr;
p1i = pti;
ak = 1.0d0 - tfnf./(fkk+tfnf);
ack = bk.*ak;
sumr = sumr +(ack+bk).*p1r;
sumi = sumi +(ack+bk).*p1i;
bk = ack;
fkk = fkk - 1.0d0;
end; i = fix(km+1);
yr(n) = p2r;
yi(n) = p2i;
if( n~=1 )
for i = 2 : n;
ptr = p2r;
pti = p2i;
p2r = p1r +(fkk+fnf).*(rzr.*ptr-rzi.*pti);
p2i = p1i +(fkk+fnf).*(rzi.*ptr+rzr.*pti);
p1r = ptr;
p1i = pti;
ak = 1.0d0 - tfnf./(fkk+tfnf);
ack = bk.*ak;
sumr = sumr +(ack+bk).*p1r;
sumi = sumi +(ack+bk).*p1i;
bk = ack;
fkk = fkk - 1.0d0;
m = fix(n - i + 1);
yr(m) = p2r;
yi(m) = p2i;
end; i = fix(n+1);
end;
if( ifnu>0 )
for i = 1 : ifnu;
ptr = p2r;
pti = p2i;
p2r = p1r +(fkk+fnf).*(rzr.*ptr-rzi.*pti);
p2i = p1i +(fkk+fnf).*(rzr.*pti+rzi.*ptr);
p1r = ptr;
p1i = pti;
ak = 1.0d0 - tfnf./(fkk+tfnf);
ack = bk.*ak;
sumr = sumr +(ack+bk).*p1r;
sumi = sumi +(ack+bk).*p1i;
bk = ack;
fkk = fkk - 1.0d0;
end; i = fix(ifnu+1);
end;
ptr = zr;
pti = zi;
if( kode==2 )
ptr = zeror;
end;
[rzr,rzi,str,sti,idum]=zlog(rzr,rzi,str,sti,idum);
p1r = -fnf.*str + ptr;
p1i = -fnf.*sti + pti;
[ap ,dumvar2,idum]=dgamln(1.0d0+fnf,idum);
ptr = p1r - ap;
pti = p1i;
%-----------------------------------------------------------------------
%     THE DIVISION CEXP(PT)/(SUM+P2) IS ALTERED TO AVOID OVERFLOW
%     IN THE DENOMINATOR BY SQUARING LARGE QUANTITIES
%-----------------------------------------------------------------------
p2r = p2r + sumr;
p2i = p2i + sumi;
[ap ,p2r,p2i]=zabs(p2r,p2i);
p1r = 1.0d0./ap;
[ptr,pti,str,sti]=zexp(ptr,pti,str,sti);
ckr = str.*p1r;
cki = sti.*p1r;
ptr = p2r.*p1r;
pti = -p2i.*p1r;
[ckr,cki,ptr,pti,cnormr,cnormi]=zmlt(ckr,cki,ptr,pti,cnormr,cnormi);
for i = 1 : n;
str = yr(i).*cnormr - yi(i).*cnormi;
yi(i) = yr(i).*cnormi + yi(i).*cnormr;
yr(i) = str;
end; i = fix(n+1);
return;
end %subroutine zmlri
%DECK ZMLT