function [dbesj1result,x]=dbesj1(x);
dbesj1result=[];
persistent ampl bj1cs first firstCall ntj1 theta xmin xsml y ; if isempty(firstCall),firstCall=1;end;
;
if isempty(ntj1), ntj1=0; end;
%***BEGIN PROLOGUE DBESJ1
%***PURPOSE Compute the Bessel function of the first kind of order one.
%***LIBRARY SLATEC (FNLIB)
%***CATEGORY C10A1
%***TYPE doubleprecision (BESJ1-S, DBESJ1-D)
%***KEYWORDS BESSEL FUNCTION, FIRST KIND, FNLIB, ORDER ONE,
% SPECIAL FUNCTIONS
%***AUTHOR Fullerton, W., (LANL)
%***DESCRIPTION
%
% DBESJ1(X) calculates the doubleprecision Bessel function of the
% first kind of order one for doubleprecision argument X.
%
% Series for BJ1 on the interval 0. to 1.60000E+01
% with weighted error 1.16E-33
% log weighted error 32.93
% significant figures required 32.36
% decimal places required 33.57
%
%***REFERENCES (NONE)
%***ROUTINES CALLED D1MACH, D9B1MP, DCSEVL, INITDS, XERMSG
%***REVISION HISTORY (YYMMDD)
% 780601 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)
% 910401 Corrected error in code which caused values to have the
% wrong sign for arguments less than 4.0. (WRB)
%***end PROLOGUE DBESJ1
if isempty(bj1cs), bj1cs=zeros(1,19); end;
if isempty(ampl), ampl=0; end;
if isempty(theta), theta=0; end;
if isempty(xsml), xsml=0; end;
if isempty(xmin), xmin=0; end;
if isempty(y), y=0; end;
if isempty(first), first=false; end;
if firstCall, bj1cs(1)=[-.117261415133327865606240574524003d+0]; end;
if firstCall, bj1cs(2)=[-.253615218307906395623030884554698d+0]; end;
if firstCall, bj1cs(3)=[+.501270809844695685053656363203743d-1]; end;
if firstCall, bj1cs(4)=[-.463151480962508191842619728789772d-2]; end;
if firstCall, bj1cs(5)=[+.247996229415914024539124064592364d-3]; end;
if firstCall, bj1cs(6)=[-.867894868627882584521246435176416d-5]; end;
if firstCall, bj1cs(7)=[+.214293917143793691502766250991292d-6]; end;
if firstCall, bj1cs(8)=[-.393609307918317979229322764073061d-8]; end;
if firstCall, bj1cs(9)=[+.559118231794688004018248059864032d-10]; end;
if firstCall, bj1cs(10)=[-.632761640466139302477695274014880d-12]; end;
if firstCall, bj1cs(11)=[+.584099161085724700326945563268266d-14]; end;
if firstCall, bj1cs(12)=[-.448253381870125819039135059199999d-16]; end;
if firstCall, bj1cs(13)=[+.290538449262502466306018688000000d-18]; end;
if firstCall, bj1cs(14)=[-.161173219784144165412118186666666d-20]; end;
if firstCall, bj1cs(15)=[+.773947881939274637298346666666666d-23]; end;
if firstCall, bj1cs(16)=[-.324869378211199841143466666666666d-25]; end;
if firstCall, bj1cs(17)=[+.120223767722741022720000000000000d-27]; end;
if firstCall, bj1cs(18)=[-.395201221265134933333333333333333d-30]; end;
if firstCall, bj1cs(19)=[+.116167808226645333333333333333333d-32]; end;
if firstCall, first=[true]; end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT DBESJ1
if( first )
[ntj1 ,bj1cs]=initds(bj1cs,19,0.1.*real(d1mach(3)));
%
xsml = sqrt(8.0d0.*d1mach(3));
xmin = 2.0d0.*d1mach(1);
end;
first = false;
%
y = abs(x);
if( y>4.0d0 )
%
[y,ampl,theta]=d9b1mp(y,ampl,theta);
dbesj1result = (abs(ampl).*sign(x)).*cos(theta);
else;
%
dbesj1result = 0.0d0;
if( y==0.0d0 )
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;
if( y<=xmin )
xermsg('SLATEC','DBESJ1','ABS(X) SO SMALL J1 UNDERFLOWS',1,1);
end;
if( y>xmin )
dbesj1result = 0.5d0.*x;
end;
if( y>xsml )
dbesj1result = x.*(.25d0+dcsevl(.125d0.*y.*y-1.0d0,bj1cs,ntj1));
end;
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 DBESJ
