| [x,fnu,kode,n,y,nz]=besk(x,fnu,kode,n,y,nz); |
function [x,fnu,kode,n,y,nz]=besk(x,fnu,kode,n,y,nz);
%***BEGIN PROLOGUE BESK
%***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 SINGLE PRECISION (BESK-S, DBESK-D)
%***KEYWORDS K BESSEL FUNCTION, SPECIAL FUNCTIONS
%***AUTHOR Amos, D. E., (SNLA)
%***DESCRIPTION
%
% Abstract
% BESK 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.0E0 and
% non-negative orders FNU. If FNU .LT. NULIM, orders FNU and
% FNU+1 are obtained from BESKNU 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.
%
% Description of Arguments
%
% Input
% X - X .GT. 0.0E0
% FNU - order of the initial K function, FNU .GE. 0.0E0
% 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 - 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.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)
%
%***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 ASYIK, BESK0, BESK0E, BESK1, BESK1E, BESKNU,
% I1MACH, R1MACH, XERMSG
%***REVISION HISTORY (YYMMDD)
% 790201 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 BESK
%
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 BESK
nn = fix(-i1mach(12));
elim = 2.303e0.*(nn.*r1mach(5)-3.0e0);
xlim = r1mach(1).*1.0e+3;
if( kode<1 || kode>2 )
%
%
%
xermsg('SLATEC','BESK','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','BESK','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','BESK','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','BESK','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','BESK','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','BESK','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','BESK','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]=asyik(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]=besknu(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]=besk0e(x);
else;
[s1 ,x]=besk0(x);
end;
if( nud==0 && nd==1 )
goto400=1;
break;
end;
if( kode==2 )
[s2 ,x]=besk1e(x);
else;
[s2 ,x]=besk1(x);
end;
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','BESK','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]=besk0e(x);
else;
[y(j) ,x]=besk0(x);
end;
if( nd==1 )
do700=0;
break;
end;
j = fix(j + 1);
end;
if( kode==2 )
[y(j) ,x]=besk1e(x);
else;
[y(j) ,x]=besk1(x);
end;
else;
[x,fnu,kode,nd,y,mz]=besknu(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 besk
%DECK BESKNU
|
|
