| [x,alpha,n,y,nz]=dbesj(x,alpha,n,y,nz); |
function [x,alpha,n,y,nz]=dbesj(x,alpha,n,y,nz);
%***BEGIN PROLOGUE DBESJ
%***PURPOSE Compute an N member sequence of J Bessel functions
% J/SUB(ALPHA+K-1)/(X), K=1,...,N for non-negative ALPHA
% and X.
%***LIBRARY SLATEC
%***CATEGORY C10A3
%***TYPE doubleprecision (BESJ-S, DBESJ-D)
%***KEYWORDS J BESSEL FUNCTION, SPECIAL FUNCTIONS
%***AUTHOR Amos, D. E., (SNLA)
% Daniel, S. L., (SNLA)
% Weston, M. K., (SNLA)
%***DESCRIPTION
%
% Abstract **** a doubleprecision routine ****
% DBESJ computes an N member sequence of J Bessel functions
% J/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 of (NU,X) not covered by one of these
% formulae, the order is incremented or decremented by integer
% values into a region where one of the formulae apply. Backward
% recursion is applied to reduce orders by integer values except
% where the entire sequence lies in the oscillatory region. In
% this case forward recursion is stable and values from the
% asymptotic expansion for X to infinity start the recursion
% when it is efficient to do so. Leading terms of the series and
% uniform expansion are tested for underflow. 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 members are set to zero.
% Overflow cannot occur.
%
% 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
% N - number of members in the sequence, N .GE. 1
%
% Output Y is doubleprecision
% Y - a vector whose first N components contain
% values for J/sub(ALPHA+K-1)/(X), K=1,...,N
% 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
% 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, DASYJY, DJAIRY, 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 DBESJ
persistent ak akm ansmlv ap arg coef dalpha dfn dtm earg elim1 etx fidal firstCall flgjy fn fnf fni fnp1 fnu fnulim gln i i1 i2 ialp idalp iflw igo igo100a igo300 igo400 igo500 igo700 igo900 in inlim is k kk km kt nn ns pdf pidt pp rden relb rtol rttp rtwo rtx rzden s s1 s2 sa sb slim sxo2 t t1 t2 ta tau tb temp tfn tm tol tolln trx tx wk xo2 xo2l ; if isempty(firstCall),firstCall=1;end;
if isempty(i), i=0; end;
if isempty(ialp), ialp=0; end;
if isempty(idalp), idalp=0; end;
if isempty(iflw), iflw=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(i2), i2=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(igo300), igo300=0; end;
if isempty(igo900), igo900=0; end;
if isempty(igo100a), igo100a=0; end;
if isempty(igo400), igo400=0; end;
if isempty(igo500), igo500=0; end;
if isempty(igo700), igo700=0; end;
%!!INTEGER i , ialp , idalp , iflw , in , inlim , is , i1 , i2 , k ,%!!kk , km , kt , N , nn , ns , Nz
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(coef), coef=0; end;
if isempty(dalpha), dalpha=0; end;
if isempty(dfn), dfn=0; end;
if isempty(dtm), dtm=0; end;
if isempty(earg), earg=0; end;
if isempty(elim1), elim1=0; end;
if isempty(etx), etx=0; end;
if isempty(fidal), fidal=0; end;
if isempty(flgjy), flgjy=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(fnulim), fnulim=zeros(1,2); end;
if isempty(gln), gln=0; end;
if isempty(pdf), pdf=0; end;
if isempty(pidt), pidt=0; end;
if isempty(pp), pp=zeros(1,4); end;
if isempty(rden), rden=0; end;
if isempty(relb), relb=0; end;
if isempty(rttp), rttp=0; end;
if isempty(rtwo), rtwo=0; end;
if isempty(rtx), rtx=0; end;
if isempty(rzden), rzden=0; end;
if isempty(s), s=0; end;
if isempty(sa), sa=0; end;
if isempty(sb), sb=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(tau), tau=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(tolln), tolln=0; end;
if isempty(trx), trx=0; end;
if isempty(tx), tx=0; end;
if isempty(t1), t1=0; end;
if isempty(t2), t2=0; end;
if isempty(wk), wk=zeros(1,7); end;
if isempty(xo2), xo2=0; end;
if isempty(xo2l), xo2l=0; end;
if isempty(slim), slim=0; end;
if isempty(rtol), rtol=0; end;
y_shape=size(y);y=reshape(y,1,[]);
if firstCall, rtwo =[1.34839972492648d+00]; end;
if firstCall, pdf =[7.85398163397448d-01]; end;
if firstCall, rttp =[7.97884560802865d-01]; end;
if firstCall, pidt=[1.57079632679490d+00]; end;
if firstCall, pp(1) =[8.72909153935547d+00]; end;
if firstCall, pp(2) =[2.65693932265030d-01]; end;
if firstCall, pp(3) =[1.24578576865586d-01]; end;
if firstCall, pp(4)=[7.70133747430388d-04]; end;
if firstCall, inlim=[150]; end;
if firstCall, fnulim(1) =[100.0d0]; end;
if firstCall, fnulim(2)=[60.0d0]; end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT DBESJ
nz = 0;
kt = 1;
ns = 0;
% I1MACH(14) REPLACES I1MACH(11) IN A doubleprecision CODE
% I1MACH(15) REPLACES I1MACH(12) IN A doubleprecision CODE
[ta ]=d1mach(3);
tol = max(ta,1.0d-15);
i1 = fix(i1mach(14) + 1);
[i2 ]=i1mach(15);
[tb ]=d1mach(5);
elim1 = -2.303d0.*(i2.*tb+3.0d0);
rtol = 1.0d0./tol;
slim = d1mach(1).*rtol.*1.0d+3;
% TOLLN = -LN(TOL)
tolln = 2.303d0.*tb.*i1;
tolln = min(tolln,34.5388d0);
if( n<1 )
xermsg('SLATEC','DBESJ','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( x<0 )
xermsg('SLATEC','DBESJ','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','DBESJ','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;
if( alpha<0.0e0 )
xermsg('SLATEC','DBESJ','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;
xo2 = x.*0.5e0;
sxo2 = xo2.*xo2;
%
% DECISION TREE FOR REGION WHERE SERIES, ASYMPTOTIC EXPANSION FOR X
% TO INFINITY AND ASYMPTOTIC EXPANSION FOR NU TO INFINITY ARE
% APPLIED.
%
igo=1;
igo300=1;
igo900=1;
while(igo==1);
igo=0;
if( sxo2<=(fnu+1.0e0) )
fn = fnu;
fnp1 = fn + 1.0e0;
xo2l = log(xo2);
is = fix(kt);
if( x<=0.50e0 )
igo300=0;
break;
end;
ns = 0;
else;
ta = max(20.0e0,fnu);
if( x>ta )
rtx = sqrt(x);
tau = rtwo.*rtx;
ta = tau + fnulim(kt);
if( fnu<=ta )
%
% ASYMPTOTIC EXPANSION FOR X TO INFINITY WITH FORWARD RECURSION IN
% OSCILLATORY REGION X.GT.MAX(20, NU), PROVIDED THE LAST MEMBER
% OF THE SEQUENCE IS ALSO IN THE REGION.
%
in = fix(fix(alpha-tau+2.0e0));
if( in<=0 )
idalp = fix(ialp);
in = 0;
else;
idalp = fix(ialp - in - 1);
kt = 1;
end;
is = fix(kt);
fidal = idalp;
dalpha = fidal + fnf;
arg = x - pidt.*dalpha - pdf;
sa = sin(arg);
sb = cos(arg);
coef = rttp./rtx;
etx = 8.0e0.*x;
igo900=0;
break;
else;
fn = fnu;
is = fix(kt);
break;
end;
elseif( x>12.0e0 ) ;
ansmlv = max(36.0e0-fnu,0.0e0);
ns = fix(fix(ansmlv));
fni = fni + ns;
dfn = fni + fnf;
fn = dfn;
is = fix(kt);
if( n-1+ns>0 )
is = 3;
end;
break;
else;
xo2l = log(xo2);
ns = fix(fix(sxo2-fnu) + 1);
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;
igo300=0;
break;
%igo=1
end;
end;
%
% UNIFORM ASYMPTOTIC EXPANSION FOR NU TO INFINITY
%
igo100a=1;
while(igo100a==1);
igo100a=0;
if(igo900==1)
igo700=1;
while(igo700==1);
igo700=0;
if(igo300==1)
i1 = fix(abs(3-is));
i1 = fix(max(i1,1));
flgjy = 1.0e0;
[dumvar1,x,fn,flgjy,i1,temp(sub2ind(size(temp),max(is,1)):end),wk,iflw]=dasyjy(@djairy,x,fn,flgjy,i1,temp(sub2ind(size(temp),max(is,1)):end),wk,iflw);
if( iflw~=0 )
%
% SET UNDERFLOW VALUE AND UPDATE PARAMETERS
% UNDERFLOW CAN ONLY OCCUR FOR NS=0 SINCE THE ORDER MUST BE
% LARGER THAN 36. THEREFORE, NS NEED NOT BE CONSIDERED.
%
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;
else;
if( is~=1 )
if( is==2 )
%from 700
break;
end;
if( is==3 )
% COMPUTATION OF LAST ORDER FOR ASYMPTOTIC EXPANSION NORMALIZATION
gln = wk(3) + wk(2);
if( wk(6)>30.0e0 )
ta = 0.5e0.*tolln./wk(4);
ta =((0.0493827160e0.*ta-0.1111111111e0).*ta+0.6666666667e0).*ta.*wk(6);
if( wk(1)<0.10e0 )
tb =(1.259921049e0+(0.1679894730e0+0.0887944358e0.*wk(1)).*wk(1))./wk(7);
else;
tb = gln./wk(5);
end;
else;
rden =(pp(4).*wk(6)+pp(3)).*wk(6) + 1.0e0;
rzden = pp(1) + pp(2).*wk(6);
ta = rzden./rden;
if( wk(1)<0.10e0 )
tb =(1.259921049e0+(0.1679894730e0+0.0887944358e0.*wk(1)).*wk(1))./wk(7);
else;
tb = gln./wk(5);
end;
end;
in = fix(fix(ta./tb+1.5e0));
%old go to 1000
if( in<=inlim )
dtm = fni + in;
trx = 2.0e0./x;
tm =(dtm+fnf).*trx;
ta = 0.0e0;
tb = tol;
kk = 1;
ak = 1.0e0;
%
% BACKWARD RECUR UNINDEXED AND SCALE WHEN MAGNITUDES ARE CLOSE TO
% UNDERFLOW LIMITS (LESS THAN SLIM=R1MACH(1)*1.0E+3/TOL)
%
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 )
s = temp(3);
sa = ta./tb;
ta = s;
tb = s;
if( abs(s)<=slim )
ta = ta.*rtol;
tb = tb.*rtol;
ak = tol;
end;
ta = ta.*sa;
kk = 2;
in = fix(ns);
if( ns~=0 )
igo=1;
end;
end;
%igo
end;
y(nn) = tb.*ak;
nz = fix(n - nn);
if( nn==1 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
k = fix(nn - 1);
s = tb;
tb = tm.*tb - ta;
ta = s;
y(k) = tb.*ak;
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;
k = fix(nn - 2);
%
% BACKWARD RECUR INDEXED
%
for i = 3 : nn;
s = tb;
tb = tm.*tb - ta;
ta = s;
y(k) = tb.*ak;
dtm = dtm - 1.0e0;
tm =(dtm+fnf).*trx;
k = fix(k - 1);
end; i = fix(nn+1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
end;
temp(1) = temp(3);
kt = 1;
end;
end;
is = 2;
fni = fni - 1.0e0;
dfn = fni + fnf;
fn = dfn;
if( i1~=2 )
igo100a=1;
continue;
end;
%from 700
break;
%
% SERIES FOR (X/2)**2.LE.NU+1
%
igo300=1;
%igo300
end;
igo300=1;
while(igo300==1);
igo300=0;
[gln ,fnp1]=dlngam(fnp1);
arg = fn.*xo2l - gln;
if( arg>=(-elim1) )
earg = exp(arg);
igo400=1;
while(igo400==1);
igo400=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 )
igo700=1;
break;
end;
if( 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.5e0./tfn)./tfn;
akm = tolln./(-ta+sqrt(ta.*ta-tolln.*tb)) + 1.5e0;
in = fix(km + fix(akm));
%old go to 1000
dtm = fni + in;
trx = 2.0e0./x;
tm =(dtm+fnf).*trx;
ta = 0.0e0;
tb = tol;
kk = 1;
ak = 1.0e0;
%
% BACKWARD RECUR UNINDEXED AND SCALE WHEN MAGNITUDES ARE CLOSE TO
% UNDERFLOW LIMITS (LESS THAN SLIM=R1MACH(1)*1.0E+3/TOL)
%
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 )
s = temp(3);
sa = ta./tb;
ta = s;
tb = s;
if( abs(s)<=slim )
ta = ta.*rtol;
tb = tb.*rtol;
ak = tol;
end;
ta = ta.*sa;
kk = 2;
in = fix(ns);
if( ns~=0 )
igo=1;
end;
end;
%igo
end;
y(nn) = tb.*ak;
nz = fix(n - nn);
if( nn==1 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
k = fix(nn - 1);
s = tb;
tb = tm.*tb - ta;
ta = s;
y(k) = tb.*ak;
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;
k = fix(nn - 2);
%
% BACKWARD RECUR INDEXED
%
for i = 3 : nn;
s = tb;
tb = tm.*tb - ta;
ta = s;
y(k) = tb.*ak;
dtm = dtm - 1.0e0;
tm =(dtm+fnf).*trx;
k = fix(k - 1);
end; i = fix(nn+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;
igo400=1;
end;
%igo400
end;
if(igo700==1)
break;
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>=(-elim1) )
igo300=1;
break;
end;
%500
end;
if(igo100a==1)
break;
end;
%igo300
end;
if(igo700==1)
break;
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;
%
% BACKWARD RECURSION SECTION
%
%igo700
end;
igo700=1;
if( ns==0 )
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;
end;
% BACKWARD RECUR FROM INDEX ALPHA+NN-1 TO ALPHA
y(nn) = temp(1);
y(nn-1) = temp(2);
if( nn==2 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
end;
trx = 2.0e0./x;
dtm = fni;
tm =(dtm+fnf).*trx;
ak = 1.0e0;
ta = temp(1);
tb = temp(2);
if( abs(ta)<=slim )
ta = ta.*rtol;
tb = tb.*rtol;
ak = tol;
end;
kk = 2;
in = fix(ns - 1);
%1200 go to
if( in==0 )
y(nn) = tb.*ak;
nz = fix(n - nn);
if( nn==1 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
k = fix(nn - 1);
s = tb;
tb = tm.*tb - ta;
ta = s;
y(k) = tb.*ak;
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;
k = fix(nn - 2);
%
% BACKWARD RECUR INDEXED
%
for i = 3 : nn;
s = tb;
tb = tm.*tb - ta;
ta = s;
y(k) = tb.*ak;
dtm = dtm - 1.0e0;
tm =(dtm+fnf).*trx;
k = fix(k - 1);
end; i = fix(nn+1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
% old go to 1100
if( ns~=0 )
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 )
s = temp(3);
sa = ta./tb;
ta = s;
tb = s;
if( abs(s)<=slim )
ta = ta.*rtol;
tb = tb.*rtol;
ak = tol;
end;
ta = ta.*sa;
kk = 2;
in = fix(ns);
if( ns~=0 )
igo=1;
end;
end;
%igo
end;
y(nn) = tb.*ak;
nz = fix(n - nn);
if( nn==1 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
k = fix(nn - 1);
s = tb;
tb = tm.*tb - ta;
ta = s;
y(k) = tb.*ak;
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;
k = fix(nn - 2);
%
% BACKWARD RECUR INDEXED
%
for i = 3 : nn;
s = tb;
tb = tm.*tb - ta;
ta = s;
y(k) = tb.*ak;
dtm = dtm - 1.0e0;
tm =(dtm+fnf).*trx;
k = fix(k - 1);
end; i = fix(nn+1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
k = fix(nn - 2);
for i = 3 : nn;
s = tb;
tb = tm.*tb - ta;
ta = s;
y(k) = tb.*ak;
k = fix(k - 1);
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;
y(1) = temp(2);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
igo900=1;
%igo900
end;
igo=1;
while(igo==1);
igo=0;
%900 label
dtm = fidal + fidal;
dtm = dtm.*dtm;
tm = 0.0e0;
if( fidal~=0.0e0 || abs(fnf)>=tol )
tm = 4.0e0.*fnf.*(fidal+fidal+fnf);
end;
trx = dtm - 1.0e0;
t2 =(trx+tm)./etx;
s2 = t2;
relb = tol.*abs(t2);
t1 = etx;
s1 = 1.0e0;
fn = 1.0e0;
ak = 8.0e0;
for k = 1 : 13;
t1 = t1 + etx;
fn = fn + ak;
trx = dtm - fn;
ap = trx + tm;
t2 = -t2.*ap./t1;
s1 = s1 + t2;
t1 = t1 + etx;
ak = ak + 8.0e0;
fn = fn + ak;
trx = dtm - fn;
ap = trx + tm;
t2 = t2.*ap./t1;
s2 = s2 + t2;
if( abs(t2)<=relb )
break;
end;
ak = ak + 8.0e0;
end;
temp(is) = coef.*(s1.*sb-s2.*sa);
if( is==2 )
%
% FORWARD RECURSION SECTION
%
if( kt==2 )
y(1) = temp(2);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
s1 = temp(1);
s2 = temp(2);
tx = 2.0e0./x;
tm = dalpha.*tx;
if( in~=0 )
%
% FORWARD RECUR TO INDEX ALPHA
%
for i = 1 : in;
s = s2;
s2 = tm.*s2 - s1;
tm = tm + tx;
s1 = s;
end; i = fix(in+1);
if( nn==1 )
y(1) = s2;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
else;
s = s2;
s2 = tm.*s2 - s1;
tm = tm + tx;
s1 = s;
end;
end;
%
% FORWARD RECUR FROM INDEX ALPHA TO ALPHA+N-1
%
y(1) = s1;
y(2) = s2;
if( nn==2 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
for i = 3 : nn;
y(i) = tm.*y(i-1) - y(i-2);
tm = tm + tx;
end; i = fix(nn+1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
else;
fidal = fidal + 1.0e0;
dalpha = fidal + fnf;
is = 2;
tb = sa;
sa = -sb;
sb = tb;
igo=1;
end;
%igo==1
end;
dtm = fni + in;
trx = 2.0e0./x;
tm =(dtm+fnf).*trx;
ta = 0.0e0;
tb = tol;
kk = 1;
ak = 1.0e0;
%
% BACKWARD RECUR UNINDEXED AND SCALE WHEN MAGNITUDES ARE CLOSE TO
% UNDERFLOW LIMITS (LESS THAN SLIM=R1MACH(1)*1.0E+3/TOL)
%
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 )
s = temp(3);
sa = ta./tb;
ta = s;
tb = s;
if( abs(s)<=slim )
ta = ta.*rtol;
tb = tb.*rtol;
ak = tol;
end;
ta = ta.*sa;
kk = 2;
in = fix(ns);
if( ns~=0 )
igo=1;
end;
end;
%igo
end;
y(nn) = tb.*ak;
nz = fix(n - nn);
if( nn==1 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
k = fix(nn - 1);
s = tb;
tb = tm.*tb - ta;
ta = s;
y(k) = tb.*ak;
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;
k = fix(nn - 2);
%
% BACKWARD RECUR INDEXED
%
for i = 3 : nn;
s = tb;
tb = tm.*tb - ta;
ta = s;
y(k) = tb.*ak;
dtm = dtm - 1.0e0;
tm =(dtm+fnf).*trx;
k = fix(k - 1);
end; i = fix(nn+1);
%igo100a
end;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
%
%
%
end %subroutine dbesj
%!! IF ( N<1 ) THEN
%!! CALL XERMSG('SLATEC','DBESJ','N LESS THAN ONE.',2,1)
%!! RETURN
%!! ELSEIF ( N==1 ) THEN
%!! kt = 2
%!! ENDIF
%!! nn = N
%!! IF ( X<0 ) THEN
%!! CALL XERMSG('SLATEC','DBESJ','X LESS THAN ZERO.',2,1)
%!! RETURN
%!! ELSEIF ( X==0 ) THEN
%!! IF ( Alpha<0 ) GOTO 1300
%!! 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 1300
%!!!
%!! ialp = INT(Alpha)
%!! fni = ialp + N - 1
%!! fnf = Alpha - ialp
%!! dfn = fni + fnf
%!! fnu = dfn
%!! xo2 = X*0.5D0
%!! sxo2 = xo2*xo2
%!!!
%!!! 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.50D0 ) GOTO 300
%!! ns = 0
%!! ELSE
%!! ta = MAX(20.0D0,fnu)
%!! IF ( X>ta ) THEN
%!! rtx = SQRT(X)
%!! tau = rtwo*rtx
%!! ta = tau + fnulim(kt)
%!! IF ( fnu<=ta ) THEN
%!!!
%!!! ASYMPTOTIC EXPANSION FOR X TO INFINITY WITH FORWARD RECURSION IN
%!!! OSCILLATORY REGION X.GT.MAX(20, NU), PROVIDED THE LAST MEMBER
%!!! OF THE SEQUENCE IS ALSO IN THE REGION.
%!!!
%!! in = INT(Alpha-tau+2.0D0)
%!! IF ( in<=0 ) THEN
%!! idalp = ialp
%!! in = 0
%!! ELSE
%!! idalp = ialp - in - 1
%!! kt = 1
%!! ENDIF
%!! is = kt
%!! fidal = idalp
%!! dalpha = fidal + fnf
%!! arg = X - pidt*dalpha - pdf
%!! sa = SIN(arg)
%!! sb = COS(arg)
%!! coef = rttp/rtx
%!! etx = 8.0D0*X
%!! GOTO 900
%!! ELSE
%!! fn = fnu
%!! is = kt
%!! GOTO 100
%!! ENDIF
%!! ELSEIF ( X>12.0D0 ) THEN
%!! ansmlv = MAX(36.0D0-fnu,0.0D0)
%!! ns = INT(ansmlv)
%!! fni = fni + ns
%!! dfn = fni + fnf
%!! fn = dfn
%!! is = kt
%!! IF ( N-1+ns>0 ) is = 3
%!! GOTO 100
%!! ELSE
%!! xo2l = LOG(xo2)
%!! ns = INT(sxo2-fnu) + 1
%!! ENDIF
%!! ENDIF
%!! fni = fni + ns
%!! dfn = fni + fnf
%!! fn = dfn
%!! fnp1 = fn + 1.0D0
%!! is = kt
%!! IF ( N-1+ns>0 ) is = 3
%!! GOTO 300
%!! ENDIF
%!!!
%!!! UNIFORM ASYMPTOTIC EXPANSION FOR NU TO INFINITY
%!!!
%!! 100 i1 = ABS(3-is)
%!! i1 = MAX(i1,1)
%!! flgjy = 1.0D0
%!! CALL DASYJY(DJAIRY,X,fn,flgjy,i1,temp(is),wk,iflw)
%!! IF ( iflw~=0 ) THEN
%!!!
%!!! SET UNDERFLOW VALUE AND UPDATE PARAMETERS
%!!! UNDERFLOW CAN ONLY OCCUR FOR NS=0 SINCE THE ORDER MUST BE LARGER
%!!! THAN 36. THEREFORE, NS NEE NOT BE TESTED.
%!!!
%!! Y(nn) = 0.0D0
%!! nn = nn - 1
%!! fni = fni - 1.0D0
%!! dfn = fni + fnf
%!! fn = dfn
%!! IF ( nn<1 ) GOTO 600
%!! IF ( nn==1 ) THEN
%!! kt = 2
%!! is = 2
%!! ENDIF
%!! GOTO 100
%!! ELSE
%!! IF ( is==1 ) GOTO 200
%!! IF ( is==2 ) GOTO 700
%!! IF ( is==3 ) THEN
%!!! COMPUTATION OF LAST ORDER FOR ASYMPTOTIC EXPANSION NORMALIZATION
%!! gln = wk(3) + wk(2)
%!! IF ( wk(6)>30.0D0 ) THEN
%!! ta = 0.5D0*tolln/wk(4)
%!! ta = ((0.0493827160D0*ta-0.1111111111D0)*ta+0.6666666667D0)%!!*ta*wk(6)
%!! IF ( wk(1)<0.10D0 ) THEN
%!! tb = (1.259921049D0+(0.1679894730D0+0.0887944358D0*wk(1))%!!*wk(1))/wk(7)
%!! ELSE
%!! tb = gln/wk(5)
%!! ENDIF
%!! ELSE
%!! rden = (pp(4)*wk(6)+pp(3))*wk(6) + 1.0D0
%!! rzden = pp(1) + pp(2)*wk(6)
%!! ta = rzden/rden
%!! IF ( wk(1)<0.10D0 ) THEN
%!! tb = (1.259921049D0+(0.1679894730D0+0.0887944358D0*wk(1))%!!*wk(1))/wk(7)
%!! ELSE
%!! tb = gln/wk(5)
%!! ENDIF
%!! ENDIF
%!! in = INT(ta/tb+1.5D0)
%!! IF ( in<=inlim ) GOTO 1000
%!! ENDIF
%!! temp(1) = temp(3)
%!! kt = 1
%!! ENDIF
%!! 200 is = 2
%!! fni = fni - 1.0D0
%!! dfn = fni + fnf
%!! fn = dfn
%!! IF ( i1~=2 ) GOTO 100
%!! GOTO 700
%!!!
%!!! SERIES FOR (X/2)**2.LE.NU+1
%!!!
%!! 300 gln = DLNGAM(fnp1)
%!! arg = fn*xo2l - gln
%!! IF ( arg<(-elim1) ) GOTO 500
%!! earg = EXP(arg)
%!! 400 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 ) GOTO 700
%!! IF ( 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.5D0/tfn)/tfn
%!! akm = tolln/(-ta+SQRT(ta*ta-tolln*tb)) + 1.5D0
%!! in = km + INT(akm)
%!! GOTO 1000
%!! ELSE
%!! earg = earg*fn/xo2
%!! fni = fni - 1.0D0
%!! dfn = fni + fnf
%!! fn = dfn
%!! is = 2
%!! GOTO 400
%!! ENDIF
%!! 500 Y(nn) = 0.0D0
%!! nn = nn - 1
%!! fnp1 = fn
%!! fni = fni - 1.0D0
%!! dfn = fni + fnf
%!! fn = dfn
%!! IF ( nn<1 ) GOTO 600
%!! IF ( nn==1 ) THEN
%!! kt = 2
%!! is = 2
%!! ENDIF
%!! IF ( sxo2>fnp1 ) GOTO 100
%!! arg = arg - xo2l + LOG(fnp1)
%!! IF ( arg>=(-elim1) ) GOTO 300
%!! GOTO 500
%!! 600 Nz = N - nn
%!! RETURN
%!!!
%!!! BACKWARD RECURSION SECTION
%!!!
%!! 700 IF ( ns==0 ) THEN
%!! Nz = N - nn
%!! IF ( kt==2 ) GOTO 800
%!!! BACKWARD RECUR FROM INDEX ALPHA+NN-1 TO ALPHA
%!! Y(nn) = temp(1)
%!! Y(nn-1) = temp(2)
%!! IF ( nn==2 ) RETURN
%!! ENDIF
%!! trx = 2.0D0/X
%!! dtm = fni
%!! tm = (dtm+fnf)*trx
%!! ak = 1.0D0
%!! ta = temp(1)
%!! tb = temp(2)
%!! IF ( ABS(ta)<=slim ) THEN
%!! ta = ta*rtol
%!! tb = tb*rtol
%!! ak = tol
%!! ENDIF
%!! kk = 2
%!! in = ns - 1
%!! IF ( in==0 ) GOTO 1200
%!! IF ( ns~=0 ) GOTO 1100
%!! k = nn - 2
%!! DO i = 3 , nn
%!! s = tb
%!! tb = tm*tb - ta
%!! ta = s
%!! Y(k) = tb*ak
%!! dtm = dtm - 1.0D0
%!! tm = (dtm+fnf)*trx
%!! k = k - 1
%!! ENDDO
%!! RETURN
%!! 800 Y(1) = temp(2)
%!! RETURN
%!! 900 dtm = fidal + fidal
%!! dtm = dtm*dtm
%!! tm = 0.0D0
%!! IF ( fidal~=0.0D0 .OR. ABS(fnf)>=tol )%!!tm = 4.0D0*fnf*(fidal+fidal+fnf)
%!! trx = dtm - 1.0D0
%!! t2 = (trx+tm)/etx
%!! s2 = t2
%!! relb = tol*ABS(t2)
%!! t1 = etx
%!! s1 = 1.0D0
%!! fn = 1.0D0
%!! ak = 8.0D0
%!! DO k = 1 , 13
%!! t1 = t1 + etx
%!! fn = fn + ak
%!! trx = dtm - fn
%!! ap = trx + tm
%!! t2 = -t2*ap/t1
%!! s1 = s1 + t2
%!! t1 = t1 + etx
%!! ak = ak + 8.0D0
%!! fn = fn + ak
%!! trx = dtm - fn
%!! ap = trx + tm
%!! t2 = t2*ap/t1
%!! s2 = s2 + t2
%!! IF ( ABS(t2)<=relb ) EXIT
%!! ak = ak + 8.0D0
%!! ENDDO
%!! temp(is) = coef*(s1*sb-s2*sa)
%!! IF ( is==2 ) THEN
%!!!
%!!! FORWARD RECURSION SECTION
%!!!
%!! IF ( kt==2 ) GOTO 800
%!! s1 = temp(1)
%!! s2 = temp(2)
%!! tx = 2.0D0/X
%!! tm = dalpha*tx
%!! IF ( in~=0 ) THEN
%!!!
%!!! FORWARD RECUR TO INDEX ALPHA
%!!!
%!! DO i = 1 , in
%!! s = s2
%!! s2 = tm*s2 - s1
%!! tm = tm + tx
%!! s1 = s
%!! ENDDO
%!! IF ( nn==1 ) THEN
%!! Y(1) = s2
%!! RETURN
%!! ELSE
%!! s = s2
%!! s2 = tm*s2 - s1
%!! tm = tm + tx
%!! s1 = s
%!! ENDIF
%!! ENDIF
%!!!
%!!! FORWARD RECUR FROM INDEX ALPHA TO ALPHA+N-1
%!!!
%!! Y(1) = s1
%!! Y(2) = s2
%!! IF ( nn==2 ) RETURN
%!! DO i = 3 , nn
%!! Y(i) = tm*Y(i-1) - Y(i-2)
%!! tm = tm + tx
%!! ENDDO
%!! RETURN
%!! ELSE
%!! fidal = fidal + 1.0D0
%!! dalpha = fidal + fnf
%!! is = 2
%!! tb = sa
%!! sa = -sb
%!! sb = tb
%!! GOTO 900
%!! ENDIF
%!! 1000 dtm = fni + in
%!! trx = 2.0D0/X
%!! tm = (dtm+fnf)*trx
%!! ta = 0.0D0
%!! tb = tol
%!! kk = 1
%!! ak = 1.0D0
%!!!
%!!! BACKWARD RECUR UNINDEXED
%!!!
%!! 1100 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
%!! s = temp(3)
%!! sa = ta/tb
%!! ta = s
%!! tb = s
%!! IF ( ABS(s)<=slim ) THEN
%!! ta = ta*rtol
%!! tb = tb*rtol
%!! ak = tol
%!! ENDIF
%!! ta = ta*sa
%!! kk = 2
%!! in = ns
%!! IF ( ns~=0 ) GOTO 1100
%!! ENDIF
%!! 1200 Y(nn) = tb*ak
%!! Nz = N - nn
%!! IF ( nn==1 ) RETURN
%!! k = nn - 1
%!! s = tb
%!! tb = tm*tb - ta
%!! ta = s
%!! Y(k) = tb*ak
%!! IF ( nn==2 ) RETURN
%!! dtm = dtm - 1.0D0
%!! tm = (dtm+fnf)*trx
%!! k = nn - 2
%!!!
%!!! BACKWARD RECUR INDEXED
%!!!
%!! DO i = 3 , nn
%!! s = tb
%!! tb = tm*tb - ta
%!! ta = s
%!! Y(k) = tb*ak
%!! dtm = dtm - 1.0D0
%!! tm = (dtm+fnf)*trx
%!! k = k - 1
%!! ENDDO
%!! RETURN
%!!!
%!!!
%!!!
%!! 1300 CALL XERMSG('SLATEC','DBESJ','ORDER, ALPHA, LESS THAN ZERO.',2,1)
%!! RETURN
%!!99999 end
%DECK DBESK0
|
|
