function [dbesy1result,x]=dbesy1(x);
dbesy1result=[];
persistent ampl by1cs first firstCall nty1 theta twodpi xmin xsml y ; if isempty(firstCall),firstCall=1;end;
;
if isempty(nty1), nty1=0; end;
%***BEGIN PROLOGUE DBESY1
%***PURPOSE Compute the Bessel function of the second kind of order
% one.
%***LIBRARY SLATEC (FNLIB)
%***CATEGORY C10A1
%***TYPE doubleprecision (BESY1-S, DBESY1-D)
%***KEYWORDS BESSEL FUNCTION, FNLIB, ORDER ONE, SECOND KIND,
% SPECIAL FUNCTIONS
%***AUTHOR Fullerton, W., (LANL)
%***DESCRIPTION
%
% DBESY1(X) calculates the doubleprecision Bessel function of the
% second kind of order for doubleprecision argument X.
%
% Series for BY1 on the interval 0. to 1.60000E+01
% with weighted error 8.65E-33
% log weighted error 32.06
% significant figures required 32.17
% decimal places required 32.71
%
%***REFERENCES (NONE)
%***ROUTINES CALLED D1MACH, D9B1MP, DBESJ1, 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 DBESY1
if isempty(by1cs), by1cs=zeros(1,20); end;
if isempty(ampl), ampl=0; end;
if isempty(theta), theta=0; end;
if isempty(twodpi), twodpi=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, by1cs(1)=[+.320804710061190862932352018628015d-1]; end;
if firstCall, by1cs(2)=[+.126270789743350044953431725999727d+1]; end;
if firstCall, by1cs(3)=[+.649996189992317500097490637314144d-2]; end;
if firstCall, by1cs(4)=[-.893616452886050411653144160009712d-1]; end;
if firstCall, by1cs(5)=[+.132508812217570954512375510370043d-1]; end;
if firstCall, by1cs(6)=[-.897905911964835237753039508298105d-3]; end;
if firstCall, by1cs(7)=[+.364736148795830678242287368165349d-4]; end;
if firstCall, by1cs(8)=[-.100137438166600055549075523845295d-5]; end;
if firstCall, by1cs(9)=[+.199453965739017397031159372421243d-7]; end;
if firstCall, by1cs(10)=[-.302306560180338167284799332520743d-9]; end;
if firstCall, by1cs(11)=[+.360987815694781196116252914242474d-11]; end;
if firstCall, by1cs(12)=[-.348748829728758242414552947409066d-13]; end;
if firstCall, by1cs(13)=[+.278387897155917665813507698517333d-15]; end;
if firstCall, by1cs(14)=[-.186787096861948768766825352533333d-17]; end;
if firstCall, by1cs(15)=[+.106853153391168259757070336000000d-19]; end;
if firstCall, by1cs(16)=[-.527472195668448228943872000000000d-22]; end;
if firstCall, by1cs(17)=[+.227019940315566414370133333333333d-24]; end;
if firstCall, by1cs(18)=[-.859539035394523108693333333333333d-27]; end;
if firstCall, by1cs(19)=[+.288540437983379456000000000000000d-29]; end;
if firstCall, by1cs(20)=[-.864754113893717333333333333333333d-32]; end;
if firstCall, twodpi=[0.636619772367581343075535053490057d0]; end;
if firstCall, first=[true]; end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT DBESY1
if( first )
[nty1 ,by1cs]=initds(by1cs,20,0.1.*real(d1mach(3)));
%
xmin = 1.571d0.*exp(max(log(d1mach(1)),-log(d1mach(2)))+0.01d0);
xsml = sqrt(4.0d0.*d1mach(3));
end;
first = false;
%
if( x<=0.0D0 )
xermsg('SLATEC','DBESY1','X IS ZERO OR NEGATIVE',1,2);
end;
if( x>4.0d0 )
%
[x,ampl,theta]=d9b1mp(x,ampl,theta);
dbesy1result = ampl.*sin(theta);
else;
%
if( x<xmin )
xermsg('SLATEC','DBESY1','X SO SMALL Y1 OVERFLOWS',3,2);
end;
y = 0.0d0;
if( x>xsml )
y = x.*x;
end;
dbesy1result = twodpi.*log(0.5d0.*x).*dbesj1(x)+(0.5d0+dcsevl(.125d0.*y-1.0d0,by1cs,nty1))./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 DBESY
