[z,fnu,kode,n,y,nz,tol,elim,alim]=cseri(z,fnu,kode,n,y,nz,tol,elim,alim);

function [z,fnu,kode,n,y,nz,tol,elim,alim]=cseri(z,fnu,kode,n,y,nz,tol,elim,alim);
%***BEGIN PROLOGUE  CSERI
%***SUBSIDIARY
%***PURPOSE  Subsidiary to CBESI and CBESK
%***LIBRARY   SLATEC
%***TYPE      ALL (CSERI-A, ZSERI-A)
%***AUTHOR  Amos, D. E., (SNL)
%***DESCRIPTION
%
%     CSERI COMPUTES THE I BESSEL FUNCTION FOR REAL(Z).GE.0.0 BY
%     MEANS OF THE POWER SERIES FOR LARGE ABS(Z) IN THE
%     REGION ABS(Z).LE.2*SQRT(FNU+1). NZ=0 IS A NORMAL RETURN.
%     NZ.GT.0 MEANS THAT THE LAST NZ COMPONENTS WERE SET TO ZERO
%     DUE TO UNDERFLOW. NZ.LT.0 MEANS UNDERFLOW OCCURRED, BUT THE
%     CONDITION ABS(Z).LE.2*SQRT(FNU+1) WAS VIOLATED AND THE
%     COMPUTATION MUST BE COMPLETED IN ANOTHER ROUTINE WITH N=N-ABS(NZ).
%
%***SEE ALSO  CBESI, CBESK
%***ROUTINES CALLED  CUCHK, GAMLN, R1MACH
%***REVISION HISTORY  (YYMMDD)
%   830501  DATE WRITTEN
%   910415  Prologue converted to Version 4.0 format.  (BAB)
%***end PROLOGUE  CSERI
persistent aa acz ak ak1 arm ascle atol az ck coef cone crsc cz czero dfnu firstCall fnup hz i ib idum iflag igo il k l m nn nw rak1 rs rtr1 rz s s1 s2 ss w x ; if isempty(firstCall),firstCall=1;end; 

if isempty(ak1), ak1=0; end;
if isempty(ck), ck=0; end;
if isempty(coef), coef=0; end;
if isempty(cone), cone=0; end;
if isempty(crsc), crsc=0; end;
if isempty(cz), cz=0; end;
if isempty(czero), czero=0; end;
if isempty(hz), hz=0; end;
if isempty(rz), rz=0; end;
if isempty(s1), s1=0; end;
if isempty(s2), s2=0; end;
if isempty(w), w=zeros(1,2); end;
if isempty(aa), aa=0; end;
if isempty(acz), acz=0; end;
if isempty(ak), ak=0; end;
if isempty(arm), arm=0; end;
if isempty(ascle), ascle=0; end;
if isempty(atol), atol=0; end;
if isempty(az), az=0; end;
if isempty(dfnu), dfnu=0; end;
if isempty(fnup), fnup=0; end;
if isempty(rak1), rak1=0; end;
if isempty(rs), rs=0; end;
if isempty(rtr1), rtr1=0; end;
if isempty(s), s=0; end;
if isempty(ss), ss=0; end;
if isempty(x), x=0; end;
if isempty(i), i=0; end;
if isempty(ib), ib=0; end;
if isempty(idum), idum=0; end;
if isempty(iflag), iflag=0; end;
if isempty(il), il=0; end;
if isempty(k), k=0; end;
if isempty(l), l=0; end;
if isempty(m), m=0; end;
if isempty(nn), nn=0; end;
if isempty(nw), nw=0; end;
if isempty(igo), igo=0; end;
if firstCall,   czero =[complex(0.0e0,0.0e0)];  end;
if firstCall,  cone=[complex(1.0e0,0.0e0)];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  CSERI
nz = 0;
az = abs(z);
igo=0;
if( az~=0.0e0 )
x = real(z);
arm = 1.0e3.*r1mach(1);
rtr1 = sqrt(arm);
crsc = complex(1.0e0,0.0e0);
iflag = 0;
if( az<arm )
nz = fix(n);
if( fnu==0.0e0 )
nz = fix(nz - 1);
end;
else;
hz = z.*complex(0.5e0,0.0e0);
cz = czero;
if( az>rtr1 )
cz = hz.*hz;
end;
acz = abs(cz);
nn = fix(n);
ck = log(hz);
while( true );
igo=0;
dfnu = fnu +(nn-1);
fnup = dfnu + 1.0e0;
%-----------------------------------------------------------------------
%     UNDERFLOW TEST
%-----------------------------------------------------------------------
ak1 = ck.*complex(dfnu,0.0e0);
[ak ,fnup,idum]=gamln(fnup,idum);
ak1 = ak1 - complex(ak,0.0e0);
if( kode==2 )
ak1 = ak1 - complex(x,0.0e0);
end;
rak1 = real(ak1);
if( rak1>(-elim) )
if( rak1<=(-alim) )
iflag = 1;
ss = 1.0e0./tol;
crsc = complex(tol,0.0e0);
ascle = arm.*ss;
end;
ak = imag(ak1);
aa = exp(rak1);
if( iflag==1 )
aa = aa.*ss;
end;
coef = complex(aa,0.0e0).*complex(cos(ak),sin(ak));
atol = tol.*acz./fnup;
il = fix(min(2,nn));
for i = 1 : il;
dfnu = fnu +(nn-i);
fnup = dfnu + 1.0e0;
s1 = cone;
if( acz>=tol.*fnup )
ak1 = cone;
ak = fnup + 2.0e0;
s = fnup;
aa = 2.0e0;
while( true );
rs = 1.0e0./s;
ak1 = ak1.*cz.*complex(rs,0.0e0);
s1 = s1 + ak1;
s = s + ak;
ak = ak + 2.0e0;
aa = aa.*acz.*rs;
if( aa<=atol )
break;
end;
end;
end;
m = fix(nn - i + 1);
s2 = s1.*coef;
w(i) = s2;
if( iflag~=0 )
[s2,nw,ascle,tol]=cuchk(s2,nw,ascle,tol);
if( nw~=0 )
igo=1;
break;
end;
end;
y(m) = s2.*crsc;
if( i~=il )
coef = coef.*complex(dfnu,0.0e0)./hz;
end;
end;
if(igo==0)
break;
end;
end;
nz = fix(nz + 1);
y(nn) = czero;
if( acz>dfnu )
%-----------------------------------------------------------------------
%     RETURN WITH NZ.LT.0 IF ABS(Z*Z/4).GT.FNU+N-NZ-1 COMPLETE
%     THE CALCULATION IN CBINU WITH N=N-ABS(NZ)
%-----------------------------------------------------------------------
nz = fix(-nz);
return;
else;
nn = fix(nn - 1);
if( nn==0 )
return;
end;
end;
end;
if( nn>2 )
k = fix(nn - 2);
ak = k;
rz =(cone+cone)./z;
igo=0;
if( iflag==1 )
%-----------------------------------------------------------------------
%     RECUR BACKWARD WITH SCALED VALUES
%-----------------------------------------------------------------------
%-----------------------------------------------------------------------
%     EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION ABOVE THE
%     UNDERFLOW LIMIT = ASCLE = R1MACH(1)*CSCL*1.0E+3
%-----------------------------------------------------------------------
s1 = w(1);
s2 = w(2);
for l = 3 : nn;
ck = s2;
s2 = s1 + complex(ak+fnu,0.0e0).*rz.*s2;
s1 = ck;
ck = s2.*crsc;
y(k) = ck;
ak = ak - 1.0e0;
k = fix(k - 1);
if( abs(ck)>ascle )
ib = fix(l + 1);
if( ib>nn )
return;
end;
igo=1;
break;
end;
end;
if(igo==0)
return;
end;
else;
ib = 3;
end;
for i = ib : nn;
y(k) = complex(ak+fnu,0.0e0).*rz.*y(k+1) + y(k+2);
ak = ak - 1.0e0;
k = fix(k - 1);
end; i = fix(nn+1);
end;
return;
end;
end;
y(1) = czero;
if( fnu==0.0e0 )
y(1) = cone;
end;
if( n~=1 )
for i = 2 : n;
y(i) = czero;
end; i = fix(n+1);
end;
end %subroutine cseri
%DECK CSEVL