[x,alpha,kode,n,y,nz]=besi(x,alpha,kode,n,y,nz);

function [x,alpha,kode,n,y,nz]=besi(x,alpha,kode,n,y,nz);
%***BEGIN PROLOGUE  BESI
%***PURPOSE  Compute an N member sequence of I Bessel functions
%            I/SUB(ALPHA+K-1)/(X), K=1,...,N or scaled Bessel functions
%            EXP(-X)*I/SUB(ALPHA+K-1)/(X), K=1,...,N for non-negative
%            ALPHA and X.
%***LIBRARY   SLATEC
%***CATEGORY  C10B3
%***TYPE      SINGLE PRECISION (BESI-S, DBESI-D)
%***KEYWORDS  I BESSEL FUNCTION, SPECIAL FUNCTIONS
%***AUTHOR  Amos, D. E., (SNLA)
%           Daniel, S. L., (SNLA)
%***DESCRIPTION
%
%     Abstract
%         BESI computes an N member sequence of I Bessel functions
%         I/sub(ALPHA+K-1)/(X), K=1,...,N or scaled Bessel functions
%         EXP(-X)*I/sub(ALPHA+K-1)/(X), K=1,...,N for non-negative ALPHA
%         and X.  A combination of the power series, the asymptotic
%         expansion for X to infinity, and the uniform asymptotic
%         expansion for NU to infinity are applied over subdivisions of
%         the (NU,X) plane.  For values not covered by one of these
%         formulae, the order is incremented by an integer so that one
%         of these formulae apply.  Backward recursion is used to reduce
%         orders by integer values.  The asymptotic expansion for X to
%         infinity is used only when the entire sequence (specifically
%         the last member) lies within the region covered by the
%         expansion.  Leading terms of these expansions are used to test
%         for over or underflow where appropriate.  If a sequence is
%         requested and the last member would underflow, the result is
%         set to zero and the next lower order tried, etc., until a
%         member comes on scale or all are set to zero.  An overflow
%         cannot occur with scaling.
%
%     Description of Arguments
%
%         Input
%           X      - X .GE. 0.0E0
%           ALPHA  - order of first member of the sequence,
%                    ALPHA .GE. 0.0E0
%           KODE   - a parameter to indicate the scaling option
%                    KODE=1 returns
%                           Y(K)=        I/sub(ALPHA+K-1)/(X),
%                                K=1,...,N
%                    KODE=2 returns
%                           Y(K)=EXP(-X)*I/sub(ALPHA+K-1)/(X),
%                                K=1,...,N
%           N      - number of members in the sequence, N .GE. 1
%
%         Output
%           Y      - a vector whose first N components contain
%                    values for I/sub(ALPHA+K-1)/(X) or scaled
%                    values for EXP(-X)*I/sub(ALPHA+K-1)/(X),
%                    K=1,...,N depending on KODE
%           NZ     - number of components of Y set to zero due to
%                    underflow,
%                    NZ=0   , normal return, computation completed
%                    NZ .NE. 0, last NZ components of Y set to zero,
%                             Y(K)=0.0E0, K=N-NZ+1,...,N.
%
%     Error Conditions
%         Improper input arguments - a fatal error
%         Overflow with KODE=1 - a fatal error
%         Underflow - a non-fatal error (NZ .NE. 0)
%
%***REFERENCES  D. E. Amos, S. L. Daniel and M. K. Weston, CDC 6600
%                 subroutines IBESS and JBESS for Bessel functions
%                 I(NU,X) and J(NU,X), X .GE. 0, NU .GE. 0, ACM
%                 Transactions on Mathematical Software 3, (1977),
%                 pp. 76-92.
%               F. W. J. Olver, Tables of Bessel Functions of Moderate
%                 or Large Orders, NPL Mathematical Tables 6, Her
%                 Majesty's Stationery Office, London, 1962.
%***ROUTINES CALLED  ALNGAM, ASYIK, I1MACH, R1MACH, XERMSG
%***REVISION HISTORY  (YYMMDD)
%   750101  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)
%   900326  Removed duplicate information from DESCRIPTION section.
%           (WRB)
%   920501  Reformatted the REFERENCES section.  (WRB)
%***end PROLOGUE  BESI
%
persistent ain ak akm ansmlv ap arg atol dfn dtm dx earg elim etx firstCall flgik fn fnf fni fnp1 fnu gln i i1 ialp igo igo1 igo100 igo1000 igo100a igo300 igo500 igo500a in inlim is k kk km kt nn ns ra rttpi s s1 s2 sx sxo2 t t2 ta tb temp tfn tm tol tolln trx xo2 xo2l z ; if isempty(firstCall),firstCall=1;end; 

if isempty(i), i=0; end;
if isempty(ialp), ialp=0; end;
if isempty(in), in=0; end;
if isempty(inlim), inlim=0; end;
if isempty(is), is=0; end;
if isempty(i1), i1=0; end;
if isempty(k), k=0; end;
if isempty(kk), kk=0; end;
if isempty(km), km=0; end;
if isempty(kt), kt=0; end;
if isempty(nn), nn=0; end;
if isempty(ns), ns=0; end;
if isempty(igo), igo=0; end;
if isempty(igo1), igo1=0; end;
if isempty(igo100), igo100=0; end;
if isempty(igo100a), igo100a=0; end;
if isempty(igo1000), igo1000=0; end;
if isempty(igo500), igo500=0; end;
if isempty(igo500a), igo500a=0; end;
if isempty(igo300), igo300=0; end;
if isempty(ain), ain=0; end;
if isempty(ak), ak=0; end;
if isempty(akm), akm=0; end;
if isempty(ansmlv), ansmlv=0; end;
if isempty(ap), ap=0; end;
if isempty(arg), arg=0; end;
if isempty(atol), atol=0; end;
if isempty(tolln), tolln=0; end;
if isempty(dfn), dfn=0; end;
if isempty(dtm), dtm=0; end;
if isempty(dx), dx=0; end;
if isempty(earg), earg=0; end;
if isempty(elim), elim=0; end;
if isempty(etx), etx=0; end;
if isempty(flgik), flgik=0; end;
if isempty(fn), fn=0; end;
if isempty(fnf), fnf=0; end;
if isempty(fni), fni=0; end;
if isempty(fnp1), fnp1=0; end;
if isempty(fnu), fnu=0; end;
if isempty(gln), gln=0; end;
if isempty(ra), ra=0; end;
if isempty(rttpi), rttpi=0; end;
if isempty(s), s=0; end;
if isempty(sx), sx=0; end;
if isempty(sxo2), sxo2=0; end;
if isempty(s1), s1=0; end;
if isempty(s2), s2=0; end;
if isempty(t), t=0; end;
if isempty(ta), ta=0; end;
if isempty(tb), tb=0; end;
if isempty(temp), temp=zeros(1,3); end;
if isempty(tfn), tfn=0; end;
if isempty(tm), tm=0; end;
if isempty(tol), tol=0; end;
if isempty(trx), trx=0; end;
if isempty(t2), t2=0; end;
if isempty(xo2), xo2=0; end;
if isempty(xo2l), xo2l=0; end;
if isempty(z), z=0; end;
y_shape=size(y);y=reshape(y,1,[]);
if firstCall,   rttpi=[3.98942280401433e-01];  end;
if firstCall,   inlim=[80];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  BESI
nz = 0;
kt = 1;
%     I1MACH(15) REPLACES I1MACH(12) IN A doubleprecision CODE
%     I1MACH(14) REPLACES I1MACH(11) IN A doubleprecision CODE
[ra ]=r1mach(3);
tol = max(ra,1.0e-15);
i1 = fix(-i1mach(12));
[gln ]=r1mach(5);
elim = 2.303e0.*(i1.*gln-3.0e0);
%     TOLLN = -LN(TOL)
i1 = fix(i1mach(11) + 1);
tolln = 2.303e0.*gln.*i1;
tolln = min(tolln,34.5388e0);
if( n<1 )
xermsg('SLATEC','BESI','N LESS THAN ONE.',2,1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
elseif( n==1 ) ;
kt = 2;
end;
nn = fix(n);
if( kode<1 || kode>2 )
%
%
%
xermsg('SLATEC','BESI','SCALING OPTION, KODE, NOT 1 OR 2.',2,1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
elseif( x<0 ) ;
xermsg('SLATEC','BESI','X LESS THAN ZERO.',2,1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
elseif( x==0 ) ;
if( alpha<0 )
xermsg('SLATEC','BESI','ORDER, ALPHA, LESS THAN ZERO.',2,1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
if( alpha==0 )
y(1) = 1.0e0;
if( n==1 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
i1 = 2;
else;
i1 = 1;
end;
for i = i1 : n;
y(i) = 0.0e0;
end; i = fix(n+1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
else;
igo=1;
igo100=1;
igo500=1;
igo1000=1;
while(igo==1);
igo=0;
if( alpha<0.0e0 )
xermsg('SLATEC','BESI','ORDER, ALPHA, LESS THAN ZERO.',2,1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
%
ialp = fix(fix(alpha));
fni = ialp + n - 1;
fnf = alpha - ialp;
dfn = fni + fnf;
fnu = dfn;
in = 0;
xo2 = x.*0.5e0;
sxo2 = xo2.*xo2;
etx = kode - 1;
sx = etx.*x;
%
%     DECISION TREE FOR REGION WHERE SERIES, ASYMPTOTIC EXPANSION FOR X
%     TO INFINITY AND ASYMPTOTIC EXPANSION FOR NU TO INFINITY ARE
%     APPLIED.
%
if( sxo2<=(fnu+1.0e0) )
fn = fnu;
fnp1 = fn + 1.0e0;
xo2l = log(xo2);
is = fix(kt);
if( x<=0.5e0 )
igo500=0;
break;
end;
ns = 0;
elseif( x<=12.0e0 ) ;
xo2l = log(xo2);
ns = fix(fix(sxo2-fnu));
else;
fn = 0.55e0.*fnu.*fnu;
fn = max(17.0e0,fn);
if( x>=fn )
%
%     ASYMPTOTIC EXPANSION FOR X TO INFINITY
%
earg = rttpi./sqrt(x);
if( kode==2 )
igo1000=0;
break;
end;
if( x>elim )
xermsg('SLATEC','BESI','OVERFLOW, X TOO LARGE FOR KODE = 1.',6,1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
else;
earg = earg.*exp(x);
igo1000=0;
break;
end;
else;
ansmlv = max(36.0e0-fnu,0.0e0);
ns = fix(fix(ansmlv));
fni = fni + ns;
dfn = fni + fnf;
fn = dfn;
is = fix(kt);
km = fix(n - 1 + ns);
if( km>0 )
is = 3;
end;
%
%     OVERFLOW TEST ON UNIFORM ASYMPTOTIC EXPANSION
%
if( kode==2 )
break;
end;
if( alpha<1.0e0 )
if( x<=elim )
break;
end;
xermsg('SLATEC','BESI','OVERFLOW, X TOO LARGE FOR KODE = 1.',6,1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
else;
z = x./alpha;
ra = sqrt(1.0e0+z.*z);
gln = log((1.0e0+ra)./z);
t = ra.*(1.0e0-etx) + etx./(z+ra);
arg = alpha.*(t-gln);
if( arg>elim )
xermsg('SLATEC','BESI','OVERFLOW, X TOO LARGE FOR KODE = 1.',6,1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
else;
if( km~=0 )
break;
end;
igo100=0;
break;
end;
end;
end;
end;
fni = fni + ns;
dfn = fni + fnf;
fn = dfn;
fnp1 = fn + 1.0e0;
is = fix(kt);
if( n-1+ns>0 )
is = 3;
end;
igo500=0;
break;
end;
end;
%
%     UNDERFLOW TEST ON UNIFORM ASYMPTOTIC EXPANSION
%
igo100a=1;
while(igo100a==1);
igo100a=0;
if(igo1000==1)
if(igo500==1)
if(igo100==1)
z = x./fn;
ra = sqrt(1.0e0+z.*z);
gln = log((1.0e0+ra)./z);
t = ra.*(1.0e0-etx) + etx./(z+ra);
arg = fn.*(t-gln);
end;
igo100=1;
if( arg<(-elim) )
%
%     SET UNDERFLOW VALUE AND UPDATE PARAMETERS
%
y(nn) = 0.0e0;
nn = fix(nn - 1);
fni = fni - 1.0e0;
dfn = fni + fnf;
fn = dfn;
if( nn<1 )
nz = fix(n - nn);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
if( nn==1 )
kt = 2;
is = 2;
end;
igo100a=1;
continue;
is = 2;
fni = fni - 1.0e0;
dfn = fni + fnf;
fn = dfn;
if( i1==2 )
%
%     BACKWARD RECURSION SECTION
%
nz = fix(n - nn);
if( kt==2 )
y(1) = temp(2);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
else;
s1 = temp(1);
s2 = temp(2);
trx = 2.0e0./x;
dtm = fni;
tm =(dtm+fnf).*trx;
if( in==0 )
%     BACKWARD RECUR FROM INDEX ALPHA+NN-1 TO ALPHA
y(nn) = s1;
y(nn-1) = s2;
if( nn==2 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
else;
%     BACKWARD RECUR TO INDEX ALPHA+NN-1
for i = 1 : in;
s = s2;
s2 = tm.*s2 + s1;
s1 = s;
dtm = dtm - 1.0e0;
tm =(dtm+fnf).*trx;
end; i = fix(in+1);
y(nn) = s1;
if( nn==1 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
y(nn-1) = s2;
if( nn==2 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
end;
k = fix(nn + 1);
for i = 3 : nn;
k = fix(k - 1);
y(k-2) = tm.*y(k-1) + y(k);
dtm = dtm - 1.0e0;
tm =(dtm+fnf).*trx;
end; i = fix(nn+1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
else;
z = x./fn;
ra = sqrt(1.0e0+z.*z);
gln = log((1.0e0+ra)./z);
t = ra.*(1.0e0-etx) + etx./(z+ra);
arg = fn.*(t-gln);
end;
end;
i1 = fix(abs(3-is));
i1 = fix(max(i1,1));
flgik = 1.0e0;
[x,fn,kode,flgik,ra,arg,i1,temp(sub2ind(size(temp),max(is,1)):end)]=asyik(x,fn,kode,flgik,ra,arg,i1,temp(sub2ind(size(temp),max(is,1)):end));
igo300=1;
while(igo300==1);
igo300=0;
%old go to 300 branch
if( is==1 )
is = 2;
fni = fni - 1.0e0;
dfn = fni + fnf;
fn = dfn;
if( i1==2 )
%
%     BACKWARD RECURSION SECTION
%
nz = fix(n - nn);
if( kt==2 )
y(1) = temp(2);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
else;
s1 = temp(1);
s2 = temp(2);
trx = 2.0e0./x;
dtm = fni;
tm =(dtm+fnf).*trx;
if( in==0 )
%     BACKWARD RECUR FROM INDEX ALPHA+NN-1 TO ALPHA
y(nn) = s1;
y(nn-1) = s2;
if( nn==2 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
else;
%     BACKWARD RECUR TO INDEX ALPHA+NN-1
for i = 1 : in;
s = s2;
s2 = tm.*s2 + s1;
s1 = s;
dtm = dtm - 1.0e0;
tm =(dtm+fnf).*trx;
end; i = fix(in+1);
y(nn) = s1;
if( nn==1 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
y(nn-1) = s2;
if( nn==2 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
end;
k = fix(nn + 1);
for i = 3 : nn;
k = fix(k - 1);
y(k-2) = tm.*y(k-1) + y(k);
dtm = dtm - 1.0e0;
tm =(dtm+fnf).*trx;
end; i = fix(nn+1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
else;
z = x./fn;
ra = sqrt(1.0e0+z.*z);
gln = log((1.0e0+ra)./z);
t = ra.*(1.0e0-etx) + etx./(z+ra);
arg = fn.*(t-gln);
end;
i1 = fix(abs(3-is));
i1 = fix(max(i1,1));
flgik = 1.0e0;
[x,fn,kode,flgik,ra,arg,i1,temp(sub2ind(size(temp),max(is,1)):end)]=asyik(x,fn,kode,flgik,ra,arg,i1,temp(sub2ind(size(temp),max(is,1)):end));
end;
if( is==2 )
nz = fix(n - nn);
if( kt==2 )
y(1) = temp(2);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
else;
s1 = temp(1);
s2 = temp(2);
trx = 2.0e0./x;
dtm = fni;
tm =(dtm+fnf).*trx;
if( in==0 )
%     BACKWARD RECUR FROM INDEX ALPHA+NN-1 TO ALPHA
y(nn) = s1;
y(nn-1) = s2;
if( nn==2 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
else;
%     BACKWARD RECUR TO INDEX ALPHA+NN-1
for i = 1 : in;
s = s2;
s2 = tm.*s2 + s1;
s1 = s;
dtm = dtm - 1.0e0;
tm =(dtm+fnf).*trx;
end; i = fix(in+1);
y(nn) = s1;
if( nn==1 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
y(nn-1) = s2;
if( nn==2 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
end;
k = fix(nn + 1);
for i = 3 : nn;
k = fix(k - 1);
y(k-2) = tm.*y(k-1) + y(k);
dtm = dtm - 1.0e0;
tm =(dtm+fnf).*trx;
end; i = fix(nn+1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
elseif( is==3 ) ;
%     COMPUTATION OF LAST ORDER FOR ASYMPTOTIC EXPANSION NORMALIZATION
t = 1.0e0./(fn.*ra);
ain = tolln./(gln+sqrt(gln.*gln+t.*tolln)) + 1.5e0;
in = fix(fix(ain));
if( in<=inlim )
trx = 2.0e0./x;
dtm = fni + in;
tm =(dtm+fnf).*trx;
ta = 0.0e0;
tb = tol;
kk = 1;
%
%     BACKWARD RECUR UNINDEXED
igo=1;
while(igo==1);
igo=0;
for i = 1 : in;
s = tb;
tb = tm.*tb + ta;
ta = s;
dtm = dtm - 1.0e0;
tm =(dtm+fnf).*trx;
end; i = fix(in+1);
%     NORMALIZATION
if( kk==1 )
ta =(ta./tb).*temp(3);
tb = temp(3);
kk = 2;
in = fix(ns);
if( ns~=0 )
igo=1;
end;
end;
end;
y(nn) = tb;
nz = fix(n - nn);
if( nn==1 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
tb = tm.*tb + ta;
k = fix(nn - 1);
y(k) = tb;
if( nn==2 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
dtm = dtm - 1.0e0;
tm =(dtm+fnf).*trx;
km = fix(k - 1);
%
%     BACKWARD RECUR INDEXED
%
for i = 1 : km;
y(k-1) = tm.*y(k) + y(k+1);
dtm = dtm - 1.0e0;
tm =(dtm+fnf).*trx;
k = fix(k - 1);
end; i = fix(km+1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
%
%     UNIFORM ASYMPTOTIC EXPANSION FOR NU TO INFINITY
%
if( km~=0 )
temp(1) = temp(3);
in = fix(ns);
kt = 1;
i1 = 0;
igo300=1;
continue;
else;
y(1) = temp(3);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
end;
%igo300
end;
%
%     SERIES FOR (X/2)**2.LE.NU+1
%
igo500=1;
end;
igo500a=1;
while(igo500a==1);
igo500a=0;
[gln ,fnp1]=alngam(fnp1);
arg = fn.*xo2l - gln - sx;
if( arg>=(-elim) )
earg = exp(arg);
igo1=1;
while(igo1==1);
igo1=0;
s = 1.0e0;
if( x>=tol )
ak = 3.0e0;
t2 = 1.0e0;
t = 1.0e0;
s1 = fn;
for k = 1 : 17;
s2 = t2 + s1;
t = t.*sxo2./s2;
s = s + t;
if( abs(t)<tol )
break;
end;
t2 = t2 + ak;
ak = ak + 2.0e0;
s1 = s1 + fn;
end;
end;
temp(is) = s.*earg;
if( is==2 )
nz = fix(n - nn);
if( kt==2 )
y(1) = temp(2);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
else;
s1 = temp(1);
s2 = temp(2);
trx = 2.0e0./x;
dtm = fni;
tm =(dtm+fnf).*trx;
if( in==0 )
%     BACKWARD RECUR FROM INDEX ALPHA+NN-1 TO ALPHA
y(nn) = s1;
y(nn-1) = s2;
if( nn==2 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
else;
%     BACKWARD RECUR TO INDEX ALPHA+NN-1
for i = 1 : in;
s = s2;
s2 = tm.*s2 + s1;
s1 = s;
dtm = dtm - 1.0e0;
tm =(dtm+fnf).*trx;
end; i = fix(in+1);
y(nn) = s1;
if( nn==1 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
y(nn-1) = s2;
if( nn==2 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
end;
k = fix(nn + 1);
for i = 3 : nn;
k = fix(k - 1);
y(k-2) = tm.*y(k-1) + y(k);
dtm = dtm - 1.0e0;
tm =(dtm+fnf).*trx;
end; i = fix(nn+1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
elseif( is==3 ) ;
%
%     BACKWARD RECURSION WITH NORMALIZATION BY
%     ASYMPTOTIC EXPANSION FOR NU TO INFINITY OR POWER SERIES.
%
%     COMPUTATION OF LAST ORDER FOR SERIES NORMALIZATION
akm = max(3.0e0-fn,0.0e0);
km = fix(fix(akm));
tfn = fn + km;
ta =(gln+tfn-0.9189385332e0-0.0833333333e0./tfn)./(tfn+0.5e0);
ta = xo2l - ta;
tb = -(1.0e0-1.0e0./tfn)./tfn;
ain = tolln./(-ta+sqrt(ta.*ta-tolln.*tb)) + 1.5e0;
in = fix(fix(ain));
in = fix(in + km);
%1200
trx = 2.0e0./x;
dtm = fni + in;
tm =(dtm+fnf).*trx;
ta = 0.0e0;
tb = tol;
kk = 1;
%
%     BACKWARD RECUR UNINDEXED
igo=1;
while(igo==1);
igo=0;
for i = 1 : in;
s = tb;
tb = tm.*tb + ta;
ta = s;
dtm = dtm - 1.0e0;
tm =(dtm+fnf).*trx;
end; i = fix(in+1);
%     NORMALIZATION
if( kk==1 )
ta =(ta./tb).*temp(3);
tb = temp(3);
kk = 2;
in = fix(ns);
if( ns~=0 )
igo=1;
end;
end;
end;
y(nn) = tb;
nz = fix(n - nn);
if( nn==1 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
tb = tm.*tb + ta;
k = fix(nn - 1);
y(k) = tb;
if( nn==2 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
dtm = dtm - 1.0e0;
tm =(dtm+fnf).*trx;
km = fix(k - 1);
%
%     BACKWARD RECUR INDEXED
%
for i = 1 : km;
y(k-1) = tm.*y(k) + y(k+1);
dtm = dtm - 1.0e0;
tm =(dtm+fnf).*trx;
k = fix(k - 1);
end; i = fix(km+1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
else;
earg = earg.*fn./xo2;
fni = fni - 1.0e0;
dfn = fni + fnf;
fn = dfn;
is = 2;
igo1=1;
end;
%igo1==1
end;
end;
while (1);
y(nn) = 0.0e0;
nn = fix(nn - 1);
fnp1 = fn;
fni = fni - 1.0e0;
dfn = fni + fnf;
fn = dfn;
if( nn<1 )
nz = fix(n - nn);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
if( nn==1 )
kt = 2;
is = 2;
end;
if( sxo2>fnp1 )
igo100a=1;
break;
end;
arg = arg - xo2l + log(fnp1);
if( arg>=(-elim) )
igo500a=1;
break;
end;
end;
if(igo100a==1)
break;
end;
%igo500a
end;
if(igo100a==1)
continue;
end;
nz = fix(n - nn);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
if( kt==2 )
y(1) = temp(2);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
else;
s1 = temp(1);
s2 = temp(2);
trx = 2.0e0./x;
dtm = fni;
tm =(dtm+fnf).*trx;
if( in==0 )
%     BACKWARD RECUR FROM INDEX ALPHA+NN-1 TO ALPHA
y(nn) = s1;
y(nn-1) = s2;
if( nn==2 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
else;
%     BACKWARD RECUR TO INDEX ALPHA+NN-1
for i = 1 : in;
s = s2;
s2 = tm.*s2 + s1;
s1 = s;
dtm = dtm - 1.0e0;
tm =(dtm+fnf).*trx;
end; i = fix(in+1);
y(nn) = s1;
if( nn==1 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
y(nn-1) = s2;
if( nn==2 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
end;
k = fix(nn + 1);
for i = 3 : nn;
k = fix(k - 1);
y(k-2) = tm.*y(k-1) + y(k);
dtm = dtm - 1.0e0;
tm =(dtm+fnf).*trx;
end; i = fix(nn+1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
igo1000=1;
%igo1000
end;
etx = 8.0e0.*x;
is = fix(kt);
in = 0;
fn = fnu;
while (1);
dx = fni + fni;
tm = 0.0e0;
if( fni~=0.0e0 || abs(fnf)>=tol )
tm = 4.0e0.*fnf.*(fni+fni+fnf);
end;
dtm = dx.*dx;
s1 = etx;
trx = dtm - 1.0e0;
dx = -(trx+tm)./etx;
t = dx;
s = 1.0e0 + dx;
atol = tol.*abs(s);
s2 = 1.0e0;
ak = 8.0e0;
for k = 1 : 25;
s1 = s1 + etx;
s2 = s2 + ak;
dx = dtm - s2;
ap = dx + tm;
t = -t.*ap./s1;
s = s + t;
if( abs(t)<=atol )
break;
end;
ak = ak + 8.0e0;
end;
temp(is) = s.*earg;
if( is==2 )
if( kt==2 )
y(1) = temp(2);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
else;
s1 = temp(1);
s2 = temp(2);
trx = 2.0e0./x;
dtm = fni;
tm =(dtm+fnf).*trx;
if( in==0 )
%     BACKWARD RECUR FROM INDEX ALPHA+NN-1 TO ALPHA
y(nn) = s1;
y(nn-1) = s2;
if( nn==2 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
else;
%     BACKWARD RECUR TO INDEX ALPHA+NN-1
for i = 1 : in;
s = s2;
s2 = tm.*s2 + s1;
s1 = s;
dtm = dtm - 1.0e0;
tm =(dtm+fnf).*trx;
end; i = fix(in+1);
y(nn) = s1;
if( nn==1 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
y(nn-1) = s2;
if( nn==2 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
end;
k = fix(nn + 1);
for i = 3 : nn;
k = fix(k - 1);
y(k-2) = tm.*y(k-1) + y(k);
dtm = dtm - 1.0e0;
tm =(dtm+fnf).*trx;
end; i = fix(nn+1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
end;
is = 2;
fni = fni - 1.0e0;
dfn = fni + fnf;
fn = dfn;
end;
trx = 2.0e0./x;
dtm = fni + in;
tm =(dtm+fnf).*trx;
ta = 0.0e0;
tb = tol;
kk = 1;
%
%     BACKWARD RECUR UNINDEXED
igo=1;
while(igo==1);
igo=0;
for i = 1 : in;
s = tb;
tb = tm.*tb + ta;
ta = s;
dtm = dtm - 1.0e0;
tm =(dtm+fnf).*trx;
end; i = fix(in+1);
%     NORMALIZATION
if( kk==1 )
ta =(ta./tb).*temp(3);
tb = temp(3);
kk = 2;
in = fix(ns);
if( ns~=0 )
igo=1;
end;
end;
end;
y(nn) = tb;
nz = fix(n - nn);
if( nn==1 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
tb = tm.*tb + ta;
k = fix(nn - 1);
y(k) = tb;
if( nn==2 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
dtm = dtm - 1.0e0;
tm =(dtm+fnf).*trx;
km = fix(k - 1);
%
%     BACKWARD RECUR INDEXED
%
for i = 1 : km;
y(k-1) = tm.*y(k) + y(k+1);
dtm = dtm - 1.0e0;
tm =(dtm+fnf).*trx;
k = fix(k - 1);
end; i = fix(km+1);
%igo100a loop ends
end;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end %subroutine besi
%DECK BESJ0