| [x,fnu,kode,n,y,nz]=dbesk(x,fnu,kode,n,y,nz); |
function [x,fnu,kode,n,y,nz]=dbesk(x,fnu,kode,n,y,nz);
%***BEGIN PROLOGUE DBESK
%***PURPOSE Implement forward recursion on the three term recursion
% relation for a sequence of non-negative order Bessel
% functions K/SUB(FNU+I-1)/(X), or scaled Bessel functions
% EXP(X)*K/SUB(FNU+I-1)/(X), I=1,...,N for real, positive
% X and non-negative orders FNU.
%***LIBRARY SLATEC
%***CATEGORY C10B3
%***TYPE doubleprecision (BESK-S, DBESK-D)
%***KEYWORDS K BESSEL FUNCTION, SPECIAL FUNCTIONS
%***AUTHOR Amos, D. E., (SNLA)
%***DESCRIPTION
%
% Abstract **** a doubleprecision routine ****
% DBESK implements forward recursion on the three term
% recursion relation for a sequence of non-negative order Bessel
% functions K/sub(FNU+I-1)/(X), or scaled Bessel functions
% EXP(X)*K/sub(FNU+I-1)/(X), I=1,..,N for real X .GT. 0.0D0 and
% non-negative orders FNU. If FNU .LT. NULIM, orders FNU and
% FNU+1 are obtained from DBSKNU to start the recursion. If
% FNU .GE. NULIM, the uniform asymptotic expansion is used for
% orders FNU and FNU+1 to start the recursion. NULIM is 35 or
% 70 depending on whether N=1 or N .GE. 2. Under and overflow
% tests are made on the leading term of the asymptotic expansion
% before any extensive computation is done.
%
% 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,FNU are doubleprecision
% X - X .GT. 0.0D0
% FNU - order of the initial K function, FNU .GE. 0.0D0
% KODE - a parameter to indicate the scaling option
% KODE=1 returns Y(I)= K/sub(FNU+I-1)/(X),
% I=1,...,N
% KODE=2 returns Y(I)=EXP(X)*K/sub(FNU+I-1)/(X),
% I=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 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 of Y set to zero due to
% underflow with KODE=1,
% NZ=0 , normal return, computation completed
% 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)
%
%***REFERENCES F. W. J. Olver, Tables of Bessel Functions of Moderate
% or Large Orders, NPL Mathematical Tables 6, Her
% Majesty's Stationery Office, London, 1962.
% 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, DASYIK, DBESK0, DBESK1, DBSK0E, DBSK1E,
% DBSKNU, I1MACH, XERMSG
%***REVISION HISTORY (YYMMDD)
% 790201 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)
% 920501 Reformatted the REFERENCES section. (WRB)
%***end PROLOGUE DBESK
%
persistent cn dnu do700 elim etx firstCall flgik fn fnn gln gnu goto200 goto300 goto400 goto500 goto600 goto78 i igo igo100a igo100b j k mz nb nd nn nud nulim rtz s s1 s2 t tm trx w xlim zn ; if isempty(firstCall),firstCall=1;end;
if isempty(i), i=0; end;
if isempty(j), j=0; end;
if isempty(k), k=0; end;
if isempty(mz), mz=0; end;
if isempty(nb), nb=0; end;
if isempty(nd), nd=0; end;
if isempty(nn), nn=0; end;
if isempty(nud), nud=0; end;
if isempty(nulim), nulim=zeros(1,2); end;
if isempty(igo), igo=0; end;
if isempty(igo100a), igo100a=0; end;
if isempty(igo100b), igo100b=0; end;
if isempty(do700), do700=0; end;
if isempty(goto78), goto78=0; end;
if isempty(goto200), goto200=0; end;
if isempty(goto300), goto300=0; end;
if isempty(goto400), goto400=0; end;
if isempty(goto500), goto500=0; end;
if isempty(goto600), goto600=0; end;
if isempty(cn), cn=0; end;
if isempty(dnu), dnu=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(fnn), fnn=0; end;
if isempty(gln), gln=0; end;
if isempty(gnu), gnu=0; end;
if isempty(rtz), rtz=0; end;
if isempty(s), s=0; end;
if isempty(s1), s1=0; end;
if isempty(s2), s2=0; end;
if isempty(t), t=0; end;
if isempty(tm), tm=0; end;
if isempty(trx), trx=0; end;
if isempty(w), w=zeros(1,2); end;
if isempty(xlim), xlim=0; end;
if isempty(zn), zn=0; end;
y_shape=size(y);y=reshape(y,1,[]);
if firstCall, nulim(1) =[35]; end;
if firstCall, nulim(2)=[70]; end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT DBESK
nn = fix(-i1mach(15));
elim = 2.303d0.*(nn.*d1mach(5)-3.0d0);
xlim = d1mach(1).*1.0d+3;
if( kode<1 || kode>2 )
%
%
%
xermsg('SLATEC','dBESK','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( fnu<0.0e0 ) ;
xermsg('SLATEC','dBESK','ORDER, FNU, LESS THAN ZERO',2,1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
elseif( x<=0.0e0 ) ;
xermsg('SLATEC','dBESK','X LESS THAN OR EQUAL TO ZERO',2,1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
elseif( x<xlim ) ;
xermsg('SLATEC','dBESK','OVERFLOW, FNU OR N TOO LARGE OR X TOO SMALL',6,1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
elseif( n<1 ) ;
xermsg('SLATEC','dBESK','N LESS THAN ONE',2,1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
else;
goto78=0;
goto200=0;
goto300=0;
goto600=0;
while (1);
etx = kode - 1;
%
% ND IS A DUMMY VARIABLE FOR N
% GNU IS A DUMMY VARIABLE FOR FNU
% NZ = NUMBER OF UNDERFLOWS ON KODE=1
%
nd = fix(n);
nz = 0;
nud = fix(fix(fnu));
dnu = fnu - nud;
gnu = fnu;
nn = fix(min(2,nd));
fn = fnu + n - 1;
fnn = fn;
if( fn<2.0e0 )
%
% UNDERFLOW TEST FOR KODE=1
if( kode==2 || x<=elim)
goto600=1;
break;
end;
goto78=1;
else;
%
% OVERFLOW TEST (LEADING EXPONENTIAL OF ASYMPTOTIC EXPANSION)
% FOR THE LAST ORDER, FNU+N-1.GE.NULIM
%
zn = x./fn;
if( zn==0.0e0 )
xermsg('SLATEC','dBESK','OVERFLOW, FNU OR N TOO LARGE OR X TOO SMALL',6,1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
else;
rtz = sqrt(1.0e0+zn.*zn);
gln = log((1.0e0+rtz)./zn);
t = rtz.*(1.0e0-etx) + etx./(zn+rtz);
cn = -fn.*(t-gln);
if( cn>elim )
xermsg('SLATEC','dBESK','OVERFLOW, FNU OR N TOO LARGE OR X TOO SMALL',6,1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
elseif( nud<nulim(nn) ) ;
%
if( kode==2 )
goto300=1;
break;
end;
%
% UNDERFLOW TEST (LEADING EXPONENTIAL OF ASYMPTOTIC EXPANSION IN X)
% FOR ORDER DNU
%
if( x<=elim )
goto300=1;
break;
end;
goto78=1;
elseif( nn==1 ) ;
goto200=1;
break;
end;
end;
end;
break;
end;
end;
%
% UNDERFLOW TEST (LEADING EXPONENTIAL OF ASYMPTOTIC EXPANSION)
% FOR THE FIRST ORDER, FNU.GE.NULIM
%
igo100a=1;
while(igo100a==1);
igo100a=0;
do700=1;
igo100b=1;
while(igo100b==1 & goto78==0);
igo100b=0;
goto500=0;
if(goto600==0)
if(goto300==0)
if(goto200==0)
fn = gnu;
zn = x./fn;
rtz = sqrt(1.0e0+zn.*zn);
gln = log((1.0e0+rtz)./zn);
t = rtz.*(1.0e0-etx) + etx./(zn+rtz);
cn = -fn.*(t-gln);
end;
goto200=0;
if( cn<-elim )
break;
end;
%
% ASYMPTOTIC EXPANSION FOR ORDERS FNU AND FNU+1.GE.NULIM
%
flgik = -1.0e0;
[x,gnu,kode,flgik,rtz,cn,nn,y]=dasyik(x,gnu,kode,flgik,rtz,cn,nn,y);
if( nn==1 )
do700=0;
break;
end;
trx = 2.0e0./x;
tm =(gnu+gnu+2.0e0)./x;
goto500=1;
%goto300
end;
goto300=0;
goto400=0;
if(goto500==0)
igo=1;
while(igo==1);
igo=0;
if( dnu~=0.0e0 )
nb = 2;
if( nud==0 && nd==1 )
nb = 1;
end;
[x,dnu,kode,nb,w,nz]=dbsknu(x,dnu,kode,nb,w,nz);
s1 = w(1);
if( nb==1 )
goto400=1;
break;
end;
s2 = w(2);
else;
if( kode==2 )
[s1 ,x]=dbsk0e(x);
else;
[s1 ,x]=dbesk0(x);
end;
%!! IF ( Kode==2 ) THEN
%!! s1 = BESK0E(X)
%!! ELSE
%!! s1 = BESK0(X)
%!! ENDIF
if( nud==0 && nd==1 )
goto400=1;
break;
end;
if( kode==2 )
[s2 ,x]=dbsk1e(x);
else;
[s2 ,x]=dbesk1(x);
end;
%!! IF ( Kode==2 ) THEN
%!! s2 = BESK1E(X)
%!! ELSE
%!! s2 = BESK1(X)
%!! ENDIF
end;
%igo
end;
if(goto400==0)
trx = 2.0e0./x;
tm =(dnu+dnu+2.0e0)./x;
% FORWARD RECUR FROM DNU TO FNU+1 TO GET Y(1) AND Y(2)
if( nd==1 )
nud = fix(nud - 1);
end;
if( nud>0 )
for i = 1 : nud;
s = s2;
s2 = tm.*s2 + s1;
s1 = s;
tm = tm + trx;
end; i = fix(nud+1);
if( nd==1 )
s1 = s2;
end;
elseif( nd<=1 ) ;
s1 = s2;
end;
end;
goto400=0;
y(1) = s1;
if( nd==1 )
do700=0;
break;
end;
y(2) = s2;
%goto500
end;
goto500=0;
if( nd~=2 )
% FORWARD RECUR FROM FNU+2 TO FNU+N-1
for i = 3 : nd;
y(i) = tm.*y(i-1) + y(i-2);
tm = tm + trx;
end; i = fix(nd+1);
end;
do700=0;
break;
%goto600
end;
goto600=0;
% OVERFLOW TEST
if( fn>1.0e0 )
if( -fn.*(log(x)-0.693e0)>elim )
xermsg('SLATEC','dBESK','OVERFLOW, FNU OR N TOO LARGE OR X TOO SMALL',6,1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
end;
if( dnu==0.0e0 )
j = fix(nud);
if( j~=1 )
j = fix(j + 1);
if( kode==2 )
[y(j) ,x]=dbsk0e(x);
else;
[y(j) ,x]=dbesk0(x);
end;
%!! IF ( Kode==2 ) THEN
%!! Y(j) = BESK0E(X)
%!! ELSE
%!! Y(j) = BESK0(X)
%!! ENDIF
if( nd==1 )
do700=0;
break;
end;
j = fix(j + 1);
end;
if( kode==2 )
[y(j) ,x]=dbsk1e(x);
else;
[y(j) ,x]=dbesk1(x);
end;
%!! IF ( Kode==2 ) THEN
%!! Y(j) = BESK1E(X)
%!! ELSE
%!! Y(j) = BESK1(X)
%!! ENDIF
else;
[x,fnu,kode,nd,y,mz]=dbsknu(x,fnu,kode,nd,y,mz);
end;
do700=0;
break;
%igo100b and goto78
end;
%
% UPDATE PARAMETERS ON UNDERFLOW
%
if(do700==1)
igo=1;
while(igo==1);
igo=0;
nud = fix(nud + 1);
nd = fix(nd - 1);
if( nd~=0 )
nn = fix(min(2,nd));
gnu = gnu + 1.0e0;
if( fnn<2.0e0 )
igo=1;
continue;
end;
if( nud>=nulim(nn) )
igo100a=1;
break;
end;
igo=1;
continue;
end;
%igo old 700
end;
if(igo100a==1)
continue;
end;
end;
nz = fix(n - nd);
if( nz==0 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
if( nd~=0 )
for i = 1 : nd;
j = fix(n - i + 1);
k = fix(nd - i + 1);
y(j) = y(k);
end; i = fix(nd+1);
end;
for i = 1 : nz;
y(i) = 0.0e0;
end; i = fix(nz+1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
%igo100a
end;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
end %subroutine dbesk
%!! INTEGER i , j , k , Kode , mz , N , nb , nd , nn , nud , nulim ,%!!Nz
%!! INTEGER I1MACH
%!! DOUBLEPRECISION cn , dnu , elim , etx , flgik , fn , fnn , Fnu ,%!!gln , gnu , rtz , s , s1 , s2 , t , tm , trx , w ,%!!X , xlim , Y , zn
%!! DOUBLEPRECISION DBESK0 , DBESK1 , DBSK1E , DBSK0E , D1MACH
%!! DIMENSION w(2) , nulim(2) , Y(*)
%!! SAVE nulim
%!! DATA nulim(1) , nulim(2)/35 , 70/
%!!!***FIRST EXECUTABLE STATEMENT DBESK
%!! nn = -I1MACH(15)
%!! elim = 2.303D0*(nn*D1MACH(5)-3.0D0)
%!! xlim = D1MACH(1)*1.0D+3
%!! IF ( Kode<1 .OR. Kode>2 ) THEN
%!!!
%!!!
%!!!
%!! CALL XERMSG('SLATEC','DBESK','SCALING OPTION, KODE, NOT 1 OR 2',%!!2,1)
%!! RETURN
%!! ELSEIF ( Fnu<0.0D0 ) THEN
%!! CALL XERMSG('SLATEC','DBESK','ORDER, FNU, LESS THAN ZERO',2,1)
%!! RETURN
%!! ELSEIF ( X<=0.0D0 ) THEN
%!! CALL XERMSG('SLATEC','DBESK','X LESS THAN OR EQUAL TO ZERO',2,1)
%!! RETURN
%!! ELSEIF ( X<xlim ) THEN
%!! CALL XERMSG('SLATEC','DBESK',%!!'OVERFLOW, FNU OR N TOO LARGE OR X TOO SMALL',6,1)
%!! RETURN
%!! ELSEIF ( N<1 ) THEN
%!! CALL XERMSG('SLATEC','DBESK','N LESS THAN ONE',2,1)
%!! RETURN
%!! ELSE
%!! etx = Kode - 1
%!!!
%!!! ND IS A DUMMY VARIABLE FOR N
%!!! GNU IS A DUMMY VARIABLE FOR FNU
%!!! NZ = NUMBER OF UNDERFLOWS ON KODE=1
%!!!
%!! nd = N
%!! Nz = 0
%!! nud = INT(Fnu)
%!! dnu = Fnu - nud
%!! gnu = Fnu
%!! nn = MIN(2,nd)
%!! fn = Fnu + N - 1
%!! fnn = fn
%!! IF ( fn<2.0D0 ) THEN
%!!!
%!!! UNDERFLOW TEST FOR KODE=1
%!! IF ( Kode==2 ) GOTO 600
%!! IF ( X<=elim ) GOTO 600
%!! GOTO 700
%!! ELSE
%!!!
%!!! OVERFLOW TEST (LEADING EXPONENTIAL OF ASYMPTOTIC EXPANSION)
%!!! FOR THE LAST ORDER, FNU+N-1.GE.NULIM
%!!!
%!! zn = X/fn
%!! IF ( zn==0.0D0 ) THEN
%!! CALL XERMSG('SLATEC','DBESK',%!!'OVERFLOW, FNU OR N TOO LARGE OR X TOO SMALL',6,1)
%!! RETURN
%!! ELSE
%!! rtz = SQRT(1.0D0+zn*zn)
%!! gln = LOG((1.0D0+rtz)/zn)
%!! t = rtz*(1.0D0-etx) + etx/(zn+rtz)
%!! cn = -fn*(t-gln)
%!! IF ( cn>elim ) THEN
%!! CALL XERMSG('SLATEC','DBESK',%!!'OVERFLOW, FNU OR N TOO LARGE OR X TOO SMALL',6,1)
%!! RETURN
%!! ELSEIF ( nud<nulim(nn) ) THEN
%!!!
%!! IF ( Kode==2 ) GOTO 300
%!!!
%!!! UNDERFLOW TEST (LEADING EXPONENTIAL OF ASYMPTOTIC EXPANSION IN X)
%!!! FOR ORDER DNU
%!!!
%!! IF ( X<=elim ) GOTO 300
%!! GOTO 700
%!! ELSEIF ( nn==1 ) THEN
%!! GOTO 200
%!! ENDIF
%!! ENDIF
%!! ENDIF
%!! ENDIF
%!!!
%!!! UNDERFLOW TEST (LEADING EXPONENTIAL OF ASYMPTOTIC EXPANSION)
%!!! FOR THE FIRST ORDER, FNU.GE.NULIM
%!!!
%!! 100 fn = gnu
%!! zn = X/fn
%!! rtz = SQRT(1.0D0+zn*zn)
%!! gln = LOG((1.0D0+rtz)/zn)
%!! t = rtz*(1.0D0-etx) + etx/(zn+rtz)
%!! cn = -fn*(t-gln)
%!! 200 IF ( cn<-elim ) GOTO 700
%!!!
%!!! ASYMPTOTIC EXPANSION FOR ORDERS FNU AND FNU+1.GE.NULIM
%!!!
%!! flgik = -1.0D0
%!! CALL DASYIK(X,gnu,Kode,flgik,rtz,cn,nn,Y)
%!! IF ( nn==1 ) GOTO 800
%!! trx = 2.0D0/X
%!! tm = (gnu+gnu+2.0D0)/X
%!! GOTO 500
%!! 300 IF ( dnu~=0.0D0 ) THEN
%!! nb = 2
%!! IF ( nud==0 .AND. nd==1 ) nb = 1
%!! CALL DBSKNU(X,dnu,Kode,nb,w,Nz)
%!! s1 = w(1)
%!! IF ( nb==1 ) GOTO 400
%!! s2 = w(2)
%!! ELSE
%!! IF ( Kode==2 ) THEN
%!! s1 = DBSK0E(X)
%!! ELSE
%!! s1 = DBESK0(X)
%!! ENDIF
%!! IF ( nud==0 .AND. nd==1 ) GOTO 400
%!! IF ( Kode==2 ) THEN
%!! s2 = DBSK1E(X)
%!! ELSE
%!! s2 = DBESK1(X)
%!! ENDIF
%!! ENDIF
%!! trx = 2.0D0/X
%!! tm = (dnu+dnu+2.0D0)/X
%!!! FORWARD RECUR FROM DNU TO FNU+1 TO GET Y(1) AND Y(2)
%!! IF ( nd==1 ) nud = nud - 1
%!! IF ( nud>0 ) THEN
%!! DO i = 1 , nud
%!! s = s2
%!! s2 = tm*s2 + s1
%!! s1 = s
%!! tm = tm + trx
%!! ENDDO
%!! IF ( nd==1 ) s1 = s2
%!! ELSEIF ( nd<=1 ) THEN
%!! s1 = s2
%!! ENDIF
%!! 400 Y(1) = s1
%!! IF ( nd==1 ) GOTO 800
%!! Y(2) = s2
%!! 500 IF ( nd~=2 ) THEN
%!!! FORWARD RECUR FROM FNU+2 TO FNU+N-1
%!! DO i = 3 , nd
%!! Y(i) = tm*Y(i-1) + Y(i-2)
%!! tm = tm + trx
%!! ENDDO
%!! ENDIF
%!! GOTO 800
%!!! OVERFLOW TEST
%!! 600 IF ( fn>1.0D0 ) THEN
%!! IF ( -fn*(LOG(X)-0.693D0)>elim ) THEN
%!! CALL XERMSG('SLATEC','DBESK',%!!'OVERFLOW, FNU OR N TOO LARGE OR X TOO SMALL',6,1)
%!! RETURN
%!! ENDIF
%!! ENDIF
%!! IF ( dnu==0.0D0 ) THEN
%!! j = nud
%!! IF ( j~=1 ) THEN
%!! j = j + 1
%!! IF ( Kode==2 ) THEN
%!! Y(j) = DBSK0E(X)
%!! ELSE
%!! Y(j) = DBESK0(X)
%!! ENDIF
%!! IF ( nd==1 ) GOTO 800
%!! j = j + 1
%!! ENDIF
%!! IF ( Kode==2 ) THEN
%!! Y(j) = DBSK1E(X)
%!! ELSE
%!! Y(j) = DBESK1(X)
%!! ENDIF
%!! ELSE
%!! CALL DBSKNU(X,Fnu,Kode,nd,Y,mz)
%!! ENDIF
%!! GOTO 800
%!!!
%!!! UPDATE PARAMETERS ON UNDERFLOW
%!!!
%!! 700 nud = nud + 1
%!! nd = nd - 1
%!! IF ( nd~=0 ) THEN
%!! nn = MIN(2,nd)
%!! gnu = gnu + 1.0D0
%!! IF ( fnn<2.0D0 ) GOTO 700
%!! IF ( nud>=nulim(nn) ) GOTO 100
%!! GOTO 700
%!! ENDIF
%!! 800 Nz = N - nd
%!! IF ( Nz==0 ) RETURN
%!! IF ( nd~=0 ) THEN
%!! DO i = 1 , nd
%!! j = N - i + 1
%!! k = nd - i + 1
%!! Y(j) = Y(k)
%!! ENDDO
%!! ENDIF
%!! DO i = 1 , Nz
%!! Y(i) = 0.0D0
%!! ENDDO
%!! RETURN
%!!99999 end
%DECK DBESKS
|
|
