| [x,fnu,kode,n,y,nz]=besknu(x,fnu,kode,n,y,nz); |
function [x,fnu,kode,n,y,nz]=besknu(x,fnu,kode,n,y,nz);
%***BEGIN PROLOGUE BESKNU
%***SUBSIDIARY
%***PURPOSE Subsidiary to BESK
%***LIBRARY SLATEC
%***TYPE SINGLE PRECISION (BESKNU-S, DBSKNU-D)
%***AUTHOR Amos, D. E., (SNLA)
%***DESCRIPTION
%
% Abstract
% BESKNU computes N member sequences of K Bessel functions
% K/SUB(FNU+I-1)/(X), I=1,N for non-negative orders FNU and
% positive X. Equations of the references are implemented on
% small orders DNU for K/SUB(DNU)/(X) and K/SUB(DNU+1)/(X).
% Forward recursion with the three term recursion relation
% generates higher orders FNU+I-1, I=1,...,N. The parameter
% KODE permits K/SUB(FNU+I-1)/(X) values or scaled values
% EXP(X)*K/SUB(FNU+I-1)/(X), I=1,N to be returned.
%
% To start the recursion FNU is normalized to the interval
% -0.5.LE.DNU.LT.0.5. A special form of the power series is
% implemented on 0.LT.X.LE.X1 while the Miller algorithm for the
% K Bessel function in terms of the confluent hypergeometric
% function U(FNU+0.5,2*FNU+1,X) is implemented on X1.LT.X.LE.X2.
% For X.GT.X2, the asymptotic expansion for large X is used.
% When FNU is a half odd integer, a special formula for
% DNU=-0.5 and DNU+1.0=0.5 is used to start the recursion.
%
% BESKNU assumes that a significant digit SINH(X) function is
% available.
%
% Description of Arguments
%
% Input
% X - X.GT.0.0E0
% FNU - Order of initial K function, FNU.GE.0.0E0
% N - Number of members of the sequence, N.GE.1
% KODE - A parameter to indicate the scaling option
% KODE= 1 returns
% Y(I)= K/SUB(FNU+I-1)/(X)
% I=1,...,N
% = 2 returns
% Y(I)=EXP(X)*K/SUB(FNU+I-1)/(X)
% I=1,...,N
%
% Output
% Y - A vector whose first N components contain values
% for the sequence
% Y(I)= K/SUB(FNU+I-1)/(X), I=1,...,N or
% Y(I)=EXP(X)*K/SUB(FNU+I-1)/(X), I=1,...,N
% depending on KODE
% NZ - Number of components set to zero due to
% underflow,
% NZ= 0 , Normal return
% NZ.NE.0 , First NZ components of Y set to zero
% due to underflow, Y(I)=0.0E0,I=1,...,NZ
%
% Error Conditions
% Improper input arguments - a fatal error
% Overflow - a fatal error
% Underflow with KODE=1 - a non-fatal error (NZ.NE.0)
%
%***SEE ALSO BESK
%***REFERENCES N. M. Temme, On the numerical evaluation of the modified
% Bessel function of the third kind, Journal of
% Computational Physics 19, (1975), pp. 324-337.
%***ROUTINES CALLED GAMMA, I1MACH, R1MACH, XERMSG
%***REVISION HISTORY (YYMMDD)
% 790201 DATE WRITTEN
% 890531 Changed all specific intrinsics to generic. (WRB)
% 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)
% 900328 Added TYPE section. (WRB)
% 900727 Added EXTERNAL statement. (WRB)
% 910408 Updated the AUTHOR and REFERENCES sections. (WRB)
% 920501 Reformatted the REFERENCES section. (WRB)
%***end PROLOGUE BESKNU
%
persistent a a1 a2 ak b bk cc ck coef cx dk dnu dnu2 elim etest ex f fc fhs firstCall fk fks flrx fmu g1 g2 gt50 i iflag igo igo1 inu j k kk koded nn p p1 p2 pi pt q rthpi rx s s1 s2 smu sqk st t1 t2 tm tol x1 x2 ; if isempty(firstCall),firstCall=1;end;
if isempty(i), i=0; end;
if isempty(iflag), iflag=0; end;
if isempty(inu), inu=0; end;
if isempty(j), j=0; end;
if isempty(k), k=0; end;
if isempty(kk), kk=0; end;
if isempty(koded), koded=0; end;
if isempty(nn), nn=0; end;
if isempty(igo), igo=0; end;
if isempty(igo1), igo1=0; end;
if isempty(gt50), gt50=0; end;
if isempty(a), a=zeros(1,160); end;
if isempty(ak), ak=0; end;
if isempty(a1), a1=0; end;
if isempty(a2), a2=0; end;
if isempty(b), b=zeros(1,160); end;
if isempty(bk), bk=0; end;
if isempty(cc), cc=zeros(1,8); end;
if isempty(ck), ck=0; end;
if isempty(coef), coef=0; end;
if isempty(cx), cx=0; end;
if isempty(dk), dk=0; end;
if isempty(dnu), dnu=0; end;
if isempty(dnu2), dnu2=0; end;
if isempty(elim), elim=0; end;
if isempty(etest), etest=0; end;
if isempty(ex), ex=0; end;
if isempty(f), f=0; end;
if isempty(fc), fc=0; end;
if isempty(fhs), fhs=0; end;
if isempty(fk), fk=0; end;
if isempty(fks), fks=0; end;
if isempty(flrx), flrx=0; end;
if isempty(fmu), fmu=0; end;
if isempty(g1), g1=0; end;
if isempty(g2), g2=0; end;
if isempty(p), p=0; end;
if isempty(pi), pi=0; end;
if isempty(pt), pt=0; end;
if isempty(p1), p1=0; end;
if isempty(p2), p2=0; end;
if isempty(q), q=0; end;
if isempty(rthpi), rthpi=0; end;
if isempty(rx), rx=0; end;
if isempty(s), s=0; end;
if isempty(smu), smu=0; end;
if isempty(sqk), sqk=0; end;
if isempty(st), st=0; end;
if isempty(s1), s1=0; end;
if isempty(s2), s2=0; end;
if isempty(tm), tm=0; end;
if isempty(tol), tol=0; end;
if isempty(t1), t1=0; end;
if isempty(t2), t2=0; end;
if isempty(x1), x1=0; end;
if isempty(x2), x2=0; end;
y_shape=size(y);y=reshape(y,1,[]);
if firstCall, x1 =[2.0e0]; end;
if firstCall, x2=[17.0e0]; end;
if firstCall, pi =[3.14159265358979e+00]; end;
if firstCall, rthpi=[1.25331413731550e+00]; end;
if firstCall, cc(1) =[5.77215664901533e-01]; end;
if firstCall, cc(2) =[-4.20026350340952e-02]; end;
if firstCall, cc(3) =[-4.21977345555443e-02]; end;
if firstCall, cc(4) =[7.21894324666300e-03]; end;
if firstCall, cc(5) =[-2.15241674114900e-04]; end;
if firstCall, cc(6) =[-2.01348547807000e-05]; end;
if firstCall, cc(7) =[1.13302723200000e-06]; end;
if firstCall, cc(8)=[6.11609500000000e-09]; end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT BESKNU
kk = fix(-i1mach(12));
elim = 2.303e0.*(kk.*r1mach(5)-3.0e0);
[ak ]=r1mach(3);
tol = max(ak,1.0e-15);
if( x<=0.0e0 )
%
%
xermsg('SLATEC','BESKNU','X NOT GREATER THAN ZERO',2,1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
elseif( fnu<0.0e0 ) ;
xermsg('SLATEC','BESKNU','FNU NOT ZERO OR POSITIVE',2,1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
else;
if( kode<1 || kode>2 )
xermsg('SLATEC','BESKNU','KODE NOT 1 OR 2',2,1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
else;
if( n<1 )
xermsg('SLATEC','BESKNU','N NOT GREATER THAN 0',2,1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
else;
igo1=1;
while(igo1==1);
igo1=0;
gt50=0;
nz = 0;
iflag = 0;
koded = fix(kode);
rx = 2.0e0./x;
inu = fix(fix(fnu+0.5e0));
dnu = fnu - inu;
if( abs(dnu)~=0.5e0 )
dnu2 = 0.0e0;
if( abs(dnu)>=tol )
dnu2 = dnu.*dnu;
end;
if( x<=x1 )
%
% SERIES FOR X.LE.X1
%
a1 = 1.0e0 - dnu;
a2 = 1.0e0 + dnu;
t1 = 1.0e0./gamma(a1);
t2 = 1.0e0./gamma(a2);
if( abs(dnu)>0.1e0 )
g1 =(t1-t2)./(dnu+dnu);
else;
% SERIES FOR F0 TO RESOLVE INDETERMINACY FOR SMALL ABS(DNU)
s = cc(1);
ak = 1.0e0;
for k = 2 : 8;
ak = ak.*dnu2;
tm = cc(k).*ak;
s = s + tm;
if( abs(tm)<tol )
break;
end;
end;
g1 = -s;
end;
g2 =(t1+t2).*0.5e0;
smu = 1.0e0;
fc = 1.0e0;
flrx = log(rx);
fmu = dnu.*flrx;
if( dnu~=0.0e0 )
fc = dnu.*pi;
fc = fc./sin(fc);
if( fmu~=0.0e0 )
smu = sinh(fmu)./fmu;
end;
end;
f = fc.*(g1.*cosh(fmu)+g2.*flrx.*smu);
fc = exp(fmu);
p = 0.5e0.*fc./t2;
q = 0.5e0./(fc.*t1);
ak = 1.0e0;
ck = 1.0e0;
bk = 1.0e0;
s1 = f;
s2 = p;
if( inu>0 || n>1 )
if( x>=tol )
cx = x.*x.*0.25e0;
while (1);
f =(ak.*f+p+q)./(bk-dnu2);
p = p./(ak-dnu);
q = q./(ak+dnu);
ck = ck.*cx./ak;
t1 = ck.*f;
s1 = s1 + t1;
t2 = ck.*(p-ak.*f);
s2 = s2 + t2;
bk = bk + ak + ak + 1.0e0;
ak = ak + 1.0e0;
s = abs(t1)./(1.0e0+abs(s1)) + abs(t2)./(1.0e0+abs(s2));
if( s<=tol )
break;
end;
end;
end;
s2 = s2.*rx;
if( koded~=1 )
f = exp(x);
s1 = s1.*f;
s2 = s2.*f;
end;
break;
else;
if( x>=tol )
cx = x.*x.*0.25e0;
while (1);
f =(ak.*f+p+q)./(bk-dnu2);
p = p./(ak-dnu);
q = q./(ak+dnu);
ck = ck.*cx./ak;
t1 = ck.*f;
s1 = s1 + t1;
bk = bk + ak + ak + 1.0e0;
ak = ak + 1.0e0;
s = abs(t1)./(1.0e0+abs(s1));
if( s<=tol )
break;
end;
end;
end;
y(1) = s1;
if( koded==1 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
y(1) = s1.*exp(x);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
end;
end;
igo=1;
while(igo==1);
igo=0;
coef = rthpi./sqrt(x);
if( koded~=2 )
if( x>elim )
koded = 2;
iflag = 1;
igo=1;
continue;
else;
coef = coef.*exp(-x);
end;
end;
end;
if( abs(dnu)==0.5e0 )
%
% FNU=HALF ODD INTEGER CASE
%
s1 = coef;
s2 = coef;
elseif( x>x2 ) ;
%
% ASYMPTOTIC EXPANSION FOR LARGE X, X.GT.X2
%
% IFLAG=0 MEANS NO UNDERFLOW OCCURRED
% IFLAG=1 MEANS AN UNDERFLOW OCCURRED- COMPUTATION PROCEEDS WITH
% KODED=2 AND A TEST FOR ON SCALE VALUES IS MADE DURING FORWARD
% RECURSION
nn = 2;
if( inu==0 && n==1 )
nn = 1;
end;
dnu2 = dnu + dnu;
fmu = 0.0e0;
if( abs(dnu2)>=tol )
fmu = dnu2.*dnu2;
end;
ex = x.*8.0e0;
s2 = 0.0e0;
for k = 1 : nn;
s1 = s2;
s = 1.0e0;
ak = 0.0e0;
ck = 1.0e0;
sqk = 1.0e0;
dk = ex;
for j = 1 : 30;
ck = ck.*(fmu-sqk)./dk;
s = s + ck;
dk = dk + ex;
ak = ak + 8.0e0;
sqk = sqk + ak;
if( abs(ck)<tol )
break;
end;
end;
s2 = s.*coef;
fmu = fmu + 8.0e0.*dnu + 4.0e0;
end;
if( nn<=1 )
s1 = s2;
gt50=1;
break;
end;
else;
%
% MILLER ALGORITHM FOR X1.LT.X.LE.X2
%
etest = cos(pi.*dnu)./(pi.*x.*tol);
fks = 1.0e0;
fhs = 0.25e0;
fk = 0.0e0;
ck = x + x + 2.0e0;
p1 = 0.0e0;
p2 = 1.0e0;
k = 0;
while (1);
k = fix(k + 1);
fk = fk + 1.0e0;
ak =(fhs-dnu2)./(fks+fk);
bk = ck./(fk+1.0e0);
pt = p2;
p2 = bk.*p2 - ak.*p1;
p1 = pt;
a(k) = ak;
b(k) = bk;
ck = ck + 2.0e0;
fks = fks + fk + fk + 1.0e0;
fhs = fhs + fk + fk;
if( etest<=fk.*p1 )
break;
end;
end;
kk = fix(k);
s = 1.0e0;
p1 = 0.0e0;
p2 = 1.0e0;
for i = 1 : k;
pt = p2;
p2 =(b(kk).*p2-p1)./a(kk);
p1 = pt;
s = s + p2;
kk = fix(kk - 1);
end; i = fix(k+1);
s1 = coef.*(p2./s);
if( inu<=0 && n<=1 )
gt50=1;
break;
end;
s2 = s1.*(x+dnu+0.5e0-p1./p2)./x;
end;
%igo1
end;
end;
%
% FORWARD RECURSION ON THE THREE TERM RECURSION RELATION
%
if(gt50==0)
ck =(dnu+dnu+2.0e0)./x;
if( n==1 )
inu = fix(inu - 1);
end;
if( inu>0 )
for i = 1 : inu;
st = s2;
s2 = ck.*s2 + s1;
s1 = st;
ck = ck + rx;
end; i = fix(inu+1);
if( n==1 )
s1 = s2;
end;
elseif( n<=1 ) ;
s1 = s2;
end;
end;
end;
if( iflag==1 )
% IFLAG=1 CASES
s = -x + log(s1);
y(1) = 0.0e0;
nz = 1;
if( s>=-elim )
y(1) = exp(s);
nz = 0;
end;
if( n==1 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
s = -x + log(s2);
y(2) = 0.0e0;
nz = fix(nz + 1);
if( s>=-elim )
nz = fix(nz - 1);
y(2) = exp(s);
end;
if( n==2 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
kk = 2;
igo=0;
if( nz>=2 )
for i = 3 : n;
kk = fix(i);
st = s2;
s2 = ck.*s2 + s1;
s1 = st;
ck = ck + rx;
s = -x + log(s2);
nz = fix(nz + 1);
y(i) = 0.0e0;
if( s>=-elim )
y(i) = exp(s);
nz = fix(nz - 1);
igo=1;
break;
end;
end;
if(igo==0)
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
end;
if( kk==n )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
s2 = s2.*ck + s1;
ck = ck + rx;
kk = fix(kk + 1);
y(kk) = exp(-x+log(s2));
if( kk==n )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
kk = fix(kk + 1);
for i = kk : n;
y(i) = ck.*y(i-1) + y(i-2);
ck = ck + rx;
end; i = fix(n+1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
else;
y(1) = s1;
if( n==1 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
y(2) = s2;
if( n==2 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
for i = 3 : n;
y(i) = ck.*y(i-1) + y(i-2);
ck = ck + rx;
end; i = fix(n+1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
end;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
end %subroutine besknu
%DECK BESKS
|
|
