| [x,alpha,n,y,nz]=besj(x,alpha,n,y,nz); |
function [x,alpha,n,y,nz]=besj(x,alpha,n,y,nz);
%***BEGIN PROLOGUE BESJ
%***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 SINGLE PRECISION (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
% BESJ 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.
%
% Description of Arguments
%
% Input
% X - X .GE. 0.0E0
% ALPHA - order of first member of the sequence,
% ALPHA .GE. 0.0E0
% N - number of members in the sequence, N .GE. 1
%
% Output
% 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.0E0, 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 ALNGAM, ASYJY, I1MACH, JAIRY, 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 BESJ
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;
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(rtol), rtol=0; end;
if isempty(slim), slim=0; end;
y_shape=size(y);y=reshape(y,1,[]);
if firstCall, rtwo =[1.34839972492648e+00]; end;
if firstCall, pdf =[7.85398163397448e-01]; end;
if firstCall, rttp =[7.97884560802865e-01]; end;
if firstCall, pidt=[1.57079632679490e+00]; end;
if firstCall, pp(1) =[8.72909153935547e+00]; end;
if firstCall, pp(2) =[2.65693932265030e-01]; end;
if firstCall, pp(3) =[1.24578576865586e-01]; end;
if firstCall, pp(4)=[7.70133747430388e-04]; end;
if firstCall, inlim=[150]; end;
if firstCall, fnulim(1) =[100.0e0]; end;
if firstCall, fnulim(2)=[60.0e0]; end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT BESJ
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 ]=r1mach(3);
tol = max(ta,1.0e-15);
i1 = fix(i1mach(11) + 1);
[i2 ]=i1mach(12);
[tb ]=r1mach(5);
elim1 = -2.303e0.*(i2.*tb+3.0e0);
rtol = 1.0e0./tol;
slim = r1mach(1).*1.0e+3.*rtol;
% TOLLN = -LN(TOL)
tolln = 2.303e0.*tb.*i1;
tolln = min(tolln,34.5388e0);
if( n<1 )
xermsg('SLATEC','BESJ','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','BESJ','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','BESJ','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','BESJ','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]=asyjy(@jairy,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]=alngam(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 besj
%DECK BESK0E
|
|
