function [dbesk1result,x]=dbesk1(x);
dbesk1result=[];
persistent bk1cs first firstCall ntk1 xmax xmaxt xmin xsml y ; if isempty(firstCall),firstCall=1;end;
;
if isempty(ntk1), ntk1=0; end;
%***BEGIN PROLOGUE DBESK1
%***PURPOSE Compute the modified (hyperbolic) Bessel function of the
% third kind of order one.
%***LIBRARY SLATEC (FNLIB)
%***CATEGORY C10B1
%***TYPE doubleprecision (BESK1-S, DBESK1-D)
%***KEYWORDS FNLIB, HYPERBOLIC BESSEL FUNCTION,
% MODIFIED BESSEL FUNCTION, ORDER ONE, SPECIAL FUNCTIONS,
% THIRD KIND
%***AUTHOR Fullerton, W., (LANL)
%***DESCRIPTION
%
% DBESK1(X) calculates the doubleprecision modified (hyperbolic)
% Bessel function of the third kind of order one for doubleprecision
% argument X. The argument must be large enough that the result does
% not overflow and small enough that the result does not underflow.
%
% Series for BK1 on the interval 0. to 4.00000E+00
% with weighted error 9.16E-32
% log weighted error 31.04
% significant figures required 30.61
% decimal places required 31.64
%
%***REFERENCES (NONE)
%***ROUTINES CALLED D1MACH, DBESI1, DBSK1E, DCSEVL, INITDS, XERMSG
%***REVISION HISTORY (YYMMDD)
% 770701 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)
%***end PROLOGUE DBESK1
if isempty(bk1cs), bk1cs=zeros(1,16); end;
if isempty(xmax), xmax=0; end;
if isempty(xmaxt), xmaxt=0; end;
if isempty(xmin), xmin=0; end;
if isempty(xsml), xsml=0; end;
if isempty(y), y=0; end;
if isempty(first), first=false; end;
if firstCall, bk1cs(1)=[+.25300227338947770532531120868533d-1]; end;
if firstCall, bk1cs(2)=[-.35315596077654487566723831691801d+0]; end;
if firstCall, bk1cs(3)=[-.12261118082265714823479067930042d+0]; end;
if firstCall, bk1cs(4)=[-.69757238596398643501812920296083d-2]; end;
if firstCall, bk1cs(5)=[-.17302889575130520630176507368979d-3]; end;
if firstCall, bk1cs(6)=[-.24334061415659682349600735030164d-5]; end;
if firstCall, bk1cs(7)=[-.22133876307347258558315252545126d-7]; end;
if firstCall, bk1cs(8)=[-.14114883926335277610958330212608d-9]; end;
if firstCall, bk1cs(9)=[-.66669016941993290060853751264373d-12]; end;
if firstCall, bk1cs(10)=[-.24274498505193659339263196864853d-14]; end;
if firstCall, bk1cs(11)=[-.70238634793862875971783797120000d-17]; end;
if firstCall, bk1cs(12)=[-.16543275155100994675491029333333d-19]; end;
if firstCall, bk1cs(13)=[-.32338347459944491991893333333333d-22]; end;
if firstCall, bk1cs(14)=[-.53312750529265274999466666666666d-25]; end;
if firstCall, bk1cs(15)=[-.75130407162157226666666666666666d-28]; end;
if firstCall, bk1cs(16)=[-.91550857176541866666666666666666d-31]; end;
if firstCall, first=[true]; end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT DBESK1
if( first )
[ntk1 ,bk1cs]=initds(bk1cs,16,0.1.*real(d1mach(3)));
xmin = exp(max(log(d1mach(1)),-log(d1mach(2)))+0.01d0);
xsml = sqrt(4.0d0.*d1mach(3));
xmaxt = -log(d1mach(1));
xmax = xmaxt - 0.5d0.*xmaxt.*log(xmaxt)./(xmaxt+0.5d0);
end;
first = false;
%
if( x<=0.0D0 )
xermsg('SLATEC','DBESK1','X IS ZERO OR NEGATIVE',2,2);
end;
if( x>2.0d0 )
%
dbesk1result = 0.0d0;
if( x>xmax )
xermsg('SLATEC','DBESK1','X SO BIG K1 UNDERFLOWS',1,1);
end;
if( x>xmax )
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',x); evalin('caller',[inputname(1),'=FUntemp;']); end
return;
end;
%
dbesk1result = exp(-x).*dbsk1e(x);
else;
%
if( x<xmin )
xermsg('SLATEC','DBESK1','X SO SMALL K1 OVERFLOWS',3,2);
end;
y = 0.0d0;
if( x>xsml )
y = x.*x;
end;
dbesk1result = log(0.5d0.*x).*dbesi1(x)+(0.75d0+dcsevl(.5d0.*y-1.0d0,bk1cs,ntk1))./x;
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',x); evalin('caller',[inputname(1),'=FUntemp;']); end
return;
end;
%
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',x); evalin('caller',[inputname(1),'=FUntemp;']); end
end
%DECK DBESK
