| [x,fnu,kode,n,y,nz]=dbsknu(x,fnu,kode,n,y,nz); |
function [x,fnu,kode,n,y,nz]=dbsknu(x,fnu,kode,n,y,nz);
%***BEGIN PROLOGUE DBSKNU
%***SUBSIDIARY
%***PURPOSE Subsidiary to DBESK
%***LIBRARY SLATEC
%***TYPE doubleprecision (BESKNU-S, DBSKNU-D)
%***AUTHOR Amos, D. E., (SNLA)
%***DESCRIPTION
%
% Abstract **** A doubleprecision routine ****
% DBSKNU 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.
%
% The maximum number of significant digits obtainable
% is the smaller of 14 and the number of digits carried in
% doubleprecision arithmetic.
%
% DBSKNU assumes that a significant digit SINH function is
% available.
%
% Description of Arguments
%
% INPUT X,FNU are doubleprecision
% X - X.GT.0.0D0
% FNU - Order of initial K function, FNU.GE.0.0D0
% 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 is doubleprecision
% 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.0D0,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 DBESK
%***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 D1MACH, DGAMMA, I1MACH, XERMSG
%***REVISION HISTORY (YYMMDD)
% 790201 DATE WRITTEN
% 890531 Changed all specific intrinsics to generic. (WRB)
% 890911 Removed unnecessary intrinsics. (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 DBSKNU
%
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.0d0]; end;
if firstCall, x2=[17.0d0]; end;
if firstCall, pi =[3.14159265358979d+00]; end;
if firstCall, rthpi=[1.25331413731550d+00]; end;
if firstCall, cc(1) =[5.77215664901533d-01]; end;
if firstCall, cc(2) =[-4.20026350340952d-02]; end;
if firstCall, cc(3) =[-4.21977345555443d-02]; end;
if firstCall, cc(4) =[7.21894324666300d-03]; end;
if firstCall, cc(5) =[-2.15241674114900d-04]; end;
if firstCall, cc(6) =[-2.01348547807000d-05]; end;
if firstCall, cc(7) =[1.13302723200000d-06]; end;
if firstCall, cc(8)=[6.11609500000000d-09]; end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT DBSKNU
kk = fix(-i1mach(15));
elim = 2.303d0.*(kk.*d1mach(5)-3.0d0);
[ak ]=d1mach(3);
tol = max(ak,1.0d-15);
if( x<=0.0e0 )
%
%
xermsg('SLATEC','DBSKNU','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','DBSKNU','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','DBSKNU','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','DBSKNU','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./dgamma(a1);
t2 = 1.0e0./dgamma(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 dbsknu
%!!
%!!IF ( X<=0.0D0 ) THEN
%!! !
%!! !
%!! CALL XERMSG('SLATEC','DBSKNU','X NOT GREATER THAN ZERO',2,1)
%!!ELSEIF ( Fnu<0.0D0 ) THEN
%!! CALL XERMSG('SLATEC','DBSKNU','FNU NOT ZERO OR POSITIVE',2,1)
%!!ELSEIF ( Kode<1 .OR. Kode>2 ) THEN
%!! CALL XERMSG('SLATEC','DBSKNU','KODE NOT 1 OR 2',2,1)
%!!ELSEIF ( N<1 ) THEN
%!! CALL XERMSG('SLATEC','DBSKNU','N NOT GREATER THAN 0',2,1)
%!!ELSE
%!! Nz = 0
%!! iflag = 0
%!! koded = Kode
%!! rx = 2.0D0/X
%!! inu = INT(Fnu+0.5D0)
%!! dnu = Fnu - inu
%!! IF ( ABS(dnu)~=0.5D0 ) THEN
%!! dnu2 = 0.0D0
%!! IF ( ABS(dnu)>=tol ) dnu2 = dnu*dnu
%!! IF ( X<=x1 ) THEN
%!! !
%!! ! SERIES FOR X.LE.X1
%!! !
%!! a1 = 1.0D0 - dnu
%!! a2 = 1.0D0 + dnu
%!! t1 = 1.0D0/DGAMMA(a1)
%!! t2 = 1.0D0/DGAMMA(a2)
%!! IF ( ABS(dnu)>0.1D0 ) THEN
%!! g1 = (t1-t2)/(dnu+dnu)
%!! ELSE
%!! ! SERIES FOR F0 TO RESOLVE INDETERMINACY FOR SMALL ABS(DNU)
%!! s = cc(1)
%!! ak = 1.0D0
%!! DO k = 2 , 8
%!! ak = ak*dnu2
%!! tm = cc(k)*ak
%!! s = s + tm
%!! IF ( ABS(tm)<tol ) EXIT
%!! ENDDO
%!! g1 = -s
%!! ENDIF
%!! g2 = (t1+t2)*0.5D0
%!! smu = 1.0D0
%!! fc = 1.0D0
%!! flrx = LOG(rx)
%!! fmu = dnu*flrx
%!! IF ( dnu~=0.0D0 ) THEN
%!! fc = dnu*pi
%!! fc = fc/SIN(fc)
%!! IF ( fmu~=0.0D0 ) smu = SINH(fmu)/fmu
%!! ENDIF
%!! f = fc*(g1*COSH(fmu)+g2*flrx*smu)
%!! fc = EXP(fmu)
%!! p = 0.5D0*fc/t2
%!! q = 0.5D0/(fc*t1)
%!! ak = 1.0D0
%!! ck = 1.0D0
%!! bk = 1.0D0
%!! s1 = f
%!! s2 = p
%!! IF ( inu>0 .OR. N>1 ) THEN
%!! IF ( X>=tol ) THEN
%!! cx = X*X*0.25D0
%!! DO WHILE ( true )
%!! 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.0D0
%!! ak = ak + 1.0D0
%!! s = ABS(t1)/(1.0D0+ABS(s1)) + ABS(t2)/(1.0D0+ABS(s2))
%!! IF ( s<=tol ) EXIT
%!! ENDDO
%!! ENDIF
%!! s2 = s2*rx
%!! IF ( koded~=1 ) THEN
%!! f = EXP(X)
%!! s1 = s1*f
%!! s2 = s2*f
%!! ENDIF
%!! GOTO 100
%!! ELSE
%!! IF ( X>=tol ) THEN
%!! cx = X*X*0.25D0
%!! DO WHILE ( true )
%!! 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.0D0
%!! ak = ak + 1.0D0
%!! s = ABS(t1)/(1.0D0+ABS(s1))
%!! IF ( s<=tol ) EXIT
%!! ENDDO
%!! ENDIF
%!! Y(1) = s1
%!! IF ( koded~=1 ) Y(1) = s1*EXP(X)
%!! RETURN
%!! ENDIF
%!! ENDIF
%!! ENDIF
%!! DO WHILE ( true )
%!! coef = rthpi/SQRT(X)
%!! IF ( koded==2 ) GOTO 50
%!! IF ( X<=elim ) EXIT
%!! koded = 2
%!! iflag = 1
%!! ENDDO
%!! coef = coef*EXP(-X)
%!!50 IF ( ABS(dnu)==0.5D0 ) THEN
%!! !
%!! ! FNU=HALF ODD INTEGER CASE
%!! !
%!! s1 = coef
%!! s2 = coef
%!! ELSEIF ( X>x2 ) THEN
%!! !
%!! ! 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 .AND. N==1 ) nn = 1
%!! dnu2 = dnu + dnu
%!! fmu = 0.0D0
%!! IF ( ABS(dnu2)>=tol ) fmu = dnu2*dnu2
%!! ex = X*8.0D0
%!! s2 = 0.0D0
%!! DO k = 1 , nn
%!! s1 = s2
%!! s = 1.0D0
%!! ak = 0.0D0
%!! ck = 1.0D0
%!! sqk = 1.0D0
%!! dk = ex
%!! DO j = 1 , 30
%!! ck = ck*(fmu-sqk)/dk
%!! s = s + ck
%!! dk = dk + ex
%!! ak = ak + 8.0D0
%!! sqk = sqk + ak
%!! IF ( ABS(ck)<tol ) EXIT
%!! ENDDO
%!! s2 = s*coef
%!! fmu = fmu + 8.0D0*dnu + 4.0D0
%!! ENDDO
%!! IF ( nn<=1 ) THEN
%!! s1 = s2
%!! GOTO 150
%!! ENDIF
%!! ELSE
%!! !
%!! ! MILLER ALGORITHM FOR X1.LT.X.LE.X2
%!! !
%!! etest = COS(pi*dnu)/(pi*X*tol)
%!! fks = 1.0D0
%!! fhs = 0.25D0
%!! fk = 0.0D0
%!! ck = X + X + 2.0D0
%!! p1 = 0.0D0
%!! p2 = 1.0D0
%!! k = 0
%!! DO WHILE ( true )
%!! k = k + 1
%!! fk = fk + 1.0D0
%!! ak = (fhs-dnu2)/(fks+fk)
%!! bk = ck/(fk+1.0D0)
%!! pt = p2
%!! p2 = bk*p2 - ak*p1
%!! p1 = pt
%!! a(k) = ak
%!! b(k) = bk
%!! ck = ck + 2.0D0
%!! fks = fks + fk + fk + 1.0D0
%!! fhs = fhs + fk + fk
%!! IF ( etest<=fk*p1 ) EXIT
%!! ENDDO
%!! kk = k
%!! s = 1.0D0
%!! p1 = 0.0D0
%!! p2 = 1.0D0
%!! DO i = 1 , k
%!! pt = p2
%!! p2 = (b(kk)*p2-p1)/a(kk)
%!! p1 = pt
%!! s = s + p2
%!! kk = kk - 1
%!! ENDDO
%!! s1 = coef*(p2/s)
%!! IF ( inu<=0 .AND. N<=1 ) GOTO 150
%!! s2 = s1*(X+dnu+0.5D0-p1/p2)/X
%!! ENDIF
%!! !
%!! ! FORWARD RECURSION ON THE THREE TERM RECURSION RELATION
%!! !
%!!100 ck = (dnu+dnu+2.0D0)/X
%!! IF ( N==1 ) inu = inu - 1
%!! IF ( inu>0 ) THEN
%!! DO i = 1 , inu
%!! st = s2
%!! s2 = ck*s2 + s1
%!! s1 = st
%!! ck = ck + rx
%!! ENDDO
%!! IF ( N==1 ) s1 = s2
%!! ELSEIF ( N<=1 ) THEN
%!! s1 = s2
%!! ENDIF
%!!150 IF ( iflag==1 ) THEN
%!! ! IFLAG=1 CASES
%!! s = -X + LOG(s1)
%!! Y(1) = 0.0D0
%!! Nz = 1
%!! IF ( s>=-elim ) THEN
%!! Y(1) = EXP(s)
%!! Nz = 0
%!! ENDIF
%!! IF ( N~=1 ) THEN
%!! s = -X + LOG(s2)
%!! Y(2) = 0.0D0
%!! Nz = Nz + 1
%!! IF ( s>=-elim ) THEN
%!! Nz = Nz - 1
%!! Y(2) = EXP(s)
%!! ENDIF
%!! IF ( N~=2 ) THEN
%!! kk = 2
%!! IF ( Nz>=2 ) THEN
%!! DO i = 3 , N
%!! kk = i
%!! st = s2
%!! s2 = ck*s2 + s1
%!! s1 = st
%!! ck = ck + rx
%!! s = -X + LOG(s2)
%!! Nz = Nz + 1
%!! Y(i) = 0.0D0
%!! IF ( s>=-elim ) THEN
%!! Y(i) = EXP(s)
%!! Nz = Nz - 1
%!! GOTO 155
%!! ENDIF
%!! ENDDO
%!! RETURN
%!! ENDIF
%!!155 IF ( kk~=N ) THEN
%!! s2 = s2*ck + s1
%!! ck = ck + rx
%!! kk = kk + 1
%!! Y(kk) = EXP(-X+LOG(s2))
%!! IF ( kk~=N ) THEN
%!! kk = kk + 1
%!! DO i = kk , N
%!! Y(i) = ck*Y(i-1) + Y(i-2)
%!! ck = ck + rx
%!! ENDDO
%!! ENDIF
%!! ENDIF
%!! ENDIF
%!! ENDIF
%!! ELSE
%!! Y(1) = s1
%!! IF ( N~=1 ) THEN
%!! Y(2) = s2
%!! IF ( N~=2 ) THEN
%!! DO i = 3 , N
%!! Y(i) = ck*Y(i-1) + Y(i-2)
%!! ck = ck + rx
%!! ENDDO
%!! ENDIF
%!! ENDIF
%!! ENDIF
%!!ENDIF
%!!end subroutine DBSKNU
%DECK DBSPDR
|
|
