| [x,alpha,kode,n,y,nz]=dbesi(x,alpha,kode,n,y,nz); |
function [x,alpha,kode,n,y,nz]=dbesi(x,alpha,kode,n,y,nz);
%***BEGIN PROLOGUE DBESI
%***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 nonnegative
% ALPHA and X.
%***LIBRARY SLATEC
%***CATEGORY C10B3
%***TYPE doubleprecision (BESI-S, DBESI-D)
%***KEYWORDS I BESSEL FUNCTION, SPECIAL FUNCTIONS
%***AUTHOR Amos, D. E., (SNLA)
% Daniel, S. L., (SNLA)
%***DESCRIPTION
%
% Abstract **** a doubleprecision routine ****
% DBESI 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 nonnegative 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.
%
% The maximum number of significant digits obtainable
% is the smaller of 14 and the number of digits carried in
% doubleprecision arithmetic.
%
% Description of Arguments
%
% Input X,ALPHA are doubleprecision
% X - X .GE. 0.0D0
% ALPHA - order of first member of the sequence,
% ALPHA .GE. 0.0D0
% 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 is doubleprecision
% 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.0D0, 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 D1MACH, DASYIK, DLNGAM, I1MACH, XERMSG
%***REVISION HISTORY (YYMMDD)
% 750101 DATE WRITTEN
% 890531 Changed all specific intrinsics to generic. (WRB)
% 890911 Removed unnecessary intrinsics. (WRB)
% 890911 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 DBESI
%
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.98942280401433d-01]; end;
if firstCall, inlim=[80]; end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT DBESI
nz = 0;
kt = 1;
% I1MACH(15) REPLACES I1MACH(12) IN A doubleprecision CODE
% I1MACH(14) REPLACES I1MACH(11) IN A doubleprecision CODE
[ra ]=d1mach(3);
tol = max(ra,1.0d-15);
i1 = fix(-i1mach(15));
[gln ]=d1mach(5);
elim = 2.303d0.*(i1.*gln-3.0d0);
% TOLLN = -LN(TOL)
i1 = fix(i1mach(14) + 1);
tolln = 2.303d0.*gln.*i1;
tolln = min(tolln,34.5388d0);
if( n<1 )
xermsg('SLATEC','DBESI','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','DBESI','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','DBESI','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','dBESI','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','dBESI','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','dBESI','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','dBESI','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','dBESI','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)]=dasyik(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)]=dasyik(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]=dlngam(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 dbesi
%!!ELSEIF ( X==0 ) THEN
%!! IF ( Alpha<0 ) GOTO 1400
%!! IF ( Alpha==0 ) THEN
%!! Y(1) = 1.0D0
%!! IF ( N==1 ) RETURN
%!! i1 = 2
%!! ELSE
%!! i1 = 1
%!! ENDIF
%!! DO i = i1 , N
%!! Y(i) = 0.0D0
%!! ENDDO
%!! RETURN
%!!ELSE
%!! IF ( Alpha<0.0D0 ) GOTO 1400
%!! !
%!! ialp = INT(Alpha)
%!! fni = ialp + N - 1
%!! fnf = Alpha - ialp
%!! dfn = fni + fnf
%!! fnu = dfn
%!! in = 0
%!! xo2 = X*0.5D0
%!! 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.0D0) ) THEN
%!! fn = fnu
%!! fnp1 = fn + 1.0D0
%!! xo2l = LOG(xo2)
%!! is = kt
%!! IF ( X<=0.5D0 ) GOTO 500
%!! ns = 0
%!! ELSEIF ( X<=12.0D0 ) THEN
%!! xo2l = LOG(xo2)
%!! ns = INT(sxo2-fnu)
%!! ELSE
%!! fn = 0.55D0*fnu*fnu
%!! fn = MAX(17.0D0,fn)
%!! IF ( X>=fn ) THEN
%!! !
%!! ! ASYMPTOTIC EXPANSION FOR X TO INFINITY
%!! !
%!! earg = rttpi/SQRT(X)
%!! IF ( Kode==2 ) GOTO 1000
%!! IF ( X>elim ) THEN
%!! CALL XERMSG('SLATEC','DBESI',%!!'OVERFLOW, X TOO LARGE FOR KODE = 1.',6,1)
%!! RETURN
%!! ELSE
%!! earg = earg*EXP(X)
%!! GOTO 1000
%!! ENDIF
%!! ELSE
%!! ansmlv = MAX(36.0D0-fnu,0.0D0)
%!! ns = INT(ansmlv)
%!! fni = fni + ns
%!! dfn = fni + fnf
%!! fn = dfn
%!! is = kt
%!! km = N - 1 + ns
%!! IF ( km>0 ) is = 3
%!! !
%!! ! OVERFLOW TEST ON UNIFORM ASYMPTOTIC EXPANSION
%!! !
%!! IF ( Kode==2 ) GOTO 100
%!! IF ( Alpha<1.0D0 ) THEN
%!! IF ( X<=elim ) GOTO 100
%!! CALL XERMSG('SLATEC','DBESI',%!!'OVERFLOW, X TOO LARGE FOR KODE = 1.',6,1)
%!! RETURN
%!! ELSE
%!! z = X/Alpha
%!! ra = SQRT(1.0D0+z*z)
%!! gln = LOG((1.0D0+ra)/z)
%!! t = ra*(1.0D0-etx) + etx/(z+ra)
%!! arg = Alpha*(t-gln)
%!! IF ( arg>elim ) THEN
%!! CALL XERMSG('SLATEC','DBESI',%!!'OVERFLOW, X TOO LARGE FOR KODE = 1.',6,1)
%!! RETURN
%!! ELSE
%!! IF ( km~=0 ) GOTO 100
%!! GOTO 200
%!! ENDIF
%!! ENDIF
%!! ENDIF
%!! ENDIF
%!! fni = fni + ns
%!! dfn = fni + fnf
%!! fn = dfn
%!! fnp1 = fn + 1.0D0
%!! is = kt
%!! IF ( N-1+ns>0 ) is = 3
%!! GOTO 500
%!!ENDIF
%!!!
%!!! UNDERFLOW TEST ON UNIFORM ASYMPTOTIC EXPANSION
%!!!
%!!100 z = X/fn
%!!ra = SQRT(1.0D0+z*z)
%!!gln = LOG((1.0D0+ra)/z)
%!!t = ra*(1.0D0-etx) + etx/(z+ra)
%!!arg = fn*(t-gln)
%!!200 IF ( arg>=(-elim) ) GOTO 400
%!!!
%!!! SET UNDERFLOW VALUE AND UPDATE PARAMETERS
%!!!
%!!Y(nn) = 0.0D0
%!!nn = nn - 1
%!!fni = fni - 1.0D0
%!!dfn = fni + fnf
%!!fn = dfn
%!!IF ( nn<1 ) GOTO 800
%!!IF ( nn==1 ) THEN
%!! kt = 2
%!! is = 2
%!!ENDIF
%!!GOTO 100
%!!300 is = 2
%!!fni = fni - 1.0D0
%!!dfn = fni + fnf
%!!fn = dfn
%!!IF ( i1==2 ) THEN
%!! !
%!! ! BACKWARD RECURSION SECTION
%!! !
%!! Nz = N - nn
%!! GOTO 900
%!!ELSE
%!! z = X/fn
%!! ra = SQRT(1.0D0+z*z)
%!! gln = LOG((1.0D0+ra)/z)
%!! t = ra*(1.0D0-etx) + etx/(z+ra)
%!! arg = fn*(t-gln)
%!!ENDIF
%!!400 i1 = ABS(3-is)
%!!i1 = MAX(i1,1)
%!!flgik = 1.0D0
%!!CALL DASYIK(X,fn,Kode,flgik,ra,arg,i1,temp(is))
%!!IF ( is==1 ) GOTO 300
%!!IF ( is==2 ) THEN
%!! Nz = N - nn
%!! GOTO 900
%!!ELSEIF ( is==3 ) THEN
%!! ! COMPUTATION OF LAST ORDER FOR ASYMPTOTIC EXPANSION NORMALIZATION
%!! t = 1.0D0/(fn*ra)
%!! ain = tolln/(gln+SQRT(gln*gln+t*tolln)) + 1.5D0
%!! in = INT(ain)
%!! IF ( in<=inlim ) GOTO 1200
%!! !
%!! ! UNIFORM ASYMPTOTIC EXPANSION FOR NU TO INFINITY
%!! !
%!! IF ( km~=0 ) THEN
%!! temp(1) = temp(3)
%!! in = ns
%!! kt = 1
%!! i1 = 0
%!! GOTO 300
%!! ELSE
%!! Y(1) = temp(3)
%!! RETURN
%!! ENDIF
%!!ENDIF
%!!!
%!!! SERIES FOR (X/2)**2.LE.NU+1
%!!!
%!!500 gln = DLNGAM(fnp1)
%!!arg = fn*xo2l - gln - sx
%!!IF ( arg<(-elim) ) GOTO 700
%!!earg = EXP(arg)
%!!600 s = 1.0D0
%!!IF ( X>=tol ) THEN
%!! ak = 3.0D0
%!! t2 = 1.0D0
%!! t = 1.0D0
%!! s1 = fn
%!! DO k = 1 , 17
%!! s2 = t2 + s1
%!! t = t*sxo2/s2
%!! s = s + t
%!! IF ( ABS(t)<tol ) EXIT
%!! t2 = t2 + ak
%!! ak = ak + 2.0D0
%!! s1 = s1 + fn
%!! ENDDO
%!!ENDIF
%!!temp(is) = s*earg
%!!IF ( is==2 ) THEN
%!! Nz = N - nn
%!! GOTO 900
%!!ELSEIF ( is==3 ) THEN
%!! !
%!! ! BACKWARD RECURSION WITH NORMALIZATION BY
%!! ! ASYMPTOTIC EXPANSION FOR NU TO INFINITY OR POWER SERIES.
%!! !
%!! ! COMPUTATION OF LAST ORDER FOR SERIES NORMALIZATION
%!! akm = MAX(3.0D0-fn,0.0D0)
%!! km = INT(akm)
%!! tfn = fn + km
%!! ta = (gln+tfn-0.9189385332D0-0.0833333333D0/tfn)/(tfn+0.5D0)
%!! ta = xo2l - ta
%!! tb = -(1.0D0-1.0D0/tfn)/tfn
%!! ain = tolln/(-ta+SQRT(ta*ta-tolln*tb)) + 1.5D0
%!! in = INT(ain)
%!! in = in + km
%!! GOTO 1200
%!!ELSE
%!! earg = earg*fn/xo2
%!! fni = fni - 1.0D0
%!! dfn = fni + fnf
%!! fn = dfn
%!! is = 2
%!! GOTO 600
%!!ENDIF
%!!700 Y(nn) = 0.0D0
%!!nn = nn - 1
%!!fnp1 = fn
%!!fni = fni - 1.0D0
%!!dfn = fni + fnf
%!!fn = dfn
%!!IF ( nn<1 ) GOTO 800
%!!IF ( nn==1 ) THEN
%!! kt = 2
%!! is = 2
%!!ENDIF
%!!IF ( sxo2>fnp1 ) GOTO 100
%!!arg = arg - xo2l + LOG(fnp1)
%!!IF ( arg>=(-elim) ) GOTO 500
%!!GOTO 700
%!!800 Nz = N - nn
%!!RETURN
%!!900 IF ( kt==2 ) THEN
%!! Y(1) = temp(2)
%!! RETURN
%!!ELSE
%!! s1 = temp(1)
%!! s2 = temp(2)
%!! trx = 2.0D0/X
%!! dtm = fni
%!! tm = (dtm+fnf)*trx
%!! IF ( in==0 ) THEN
%!! ! BACKWARD RECUR FROM INDEX ALPHA+NN-1 TO ALPHA
%!! Y(nn) = s1
%!! Y(nn-1) = s2
%!! IF ( nn==2 ) RETURN
%!! ELSE
%!! ! BACKWARD RECUR TO INDEX ALPHA+NN-1
%!! DO i = 1 , in
%!! s = s2
%!! s2 = tm*s2 + s1
%!! s1 = s
%!! dtm = dtm - 1.0D0
%!! tm = (dtm+fnf)*trx
%!! ENDDO
%!! Y(nn) = s1
%!! IF ( nn==1 ) RETURN
%!! Y(nn-1) = s2
%!! IF ( nn==2 ) RETURN
%!! ENDIF
%!! k = nn + 1
%!! DO i = 3 , nn
%!! k = k - 1
%!! Y(k-2) = tm*Y(k-1) + Y(k)
%!! dtm = dtm - 1.0D0
%!! tm = (dtm+fnf)*trx
%!! ENDDO
%!! RETURN
%!!ENDIF
%!!1000 etx = 8.0D0*X
%!!is = kt
%!!in = 0
%!!fn = fnu
%!!1100 dx = fni + fni
%!!tm = 0.0D0
%!!IF ( fni~=0.0D0 .OR. ABS(fnf)>=tol ) tm = 4.0D0*fnf*(fni+fni+fnf)
%!!dtm = dx*dx
%!!s1 = etx
%!!trx = dtm - 1.0D0
%!!dx = -(trx+tm)/etx
%!!t = dx
%!!s = 1.0D0 + dx
%!!atol = tol*ABS(s)
%!!s2 = 1.0D0
%!!ak = 8.0D0
%!!DO 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 ) EXIT
%!! ak = ak + 8.0D0
%!!ENDDO
%!!temp(is) = s*earg
%!!IF ( is==2 ) GOTO 900
%!!is = 2
%!!fni = fni - 1.0D0
%!!dfn = fni + fnf
%!!fn = dfn
%!!GOTO 1100
%!!1200 trx = 2.0D0/X
%!!dtm = fni + in
%!!tm = (dtm+fnf)*trx
%!!ta = 0.0D0
%!!tb = tol
%!!kk = 1
%!!!
%!!! BACKWARD RECUR UNINDEXED
%!!!
%!!1300 DO i = 1 , in
%!! s = tb
%!! tb = tm*tb + ta
%!! ta = s
%!! dtm = dtm - 1.0D0
%!! tm = (dtm+fnf)*trx
%!!ENDDO
%!!! NORMALIZATION
%!!IF ( kk==1 ) THEN
%!! ta = (ta/tb)*temp(3)
%!! tb = temp(3)
%!! kk = 2
%!! in = ns
%!! IF ( ns~=0 ) GOTO 1300
%!!ENDIF
%!!Y(nn) = tb
%!!Nz = N - nn
%!!IF ( nn==1 ) RETURN
%!!tb = tm*tb + ta
%!!k = nn - 1
%!!Y(k) = tb
%!!IF ( nn==2 ) RETURN
%!!dtm = dtm - 1.0D0
%!!tm = (dtm+fnf)*trx
%!!km = k - 1
%!!!
%!!! BACKWARD RECUR INDEXED
%!!!
%!!DO i = 1 , km
%!! Y(k-1) = tm*Y(k) + Y(k+1)
%!! dtm = dtm - 1.0D0
%!! tm = (dtm+fnf)*trx
%!! k = k - 1
%!!ENDDO
%!!RETURN
%!!1400 CALL XERMSG('SLATEC','DBESI','ORDER, ALPHA, LESS THAN ZERO.',2,1)
%!!RETURN
%!!99999 end subroutine DBESI
%DECK DBESJ0
|
|
