function [besj1result,x]=besj1(x);
besj1result=[];
persistent ampl besj1 bj1cs bm1cs bth1cs first firstCall ntj1 ntm1 ntth1 pi4 theta xmax xmin xsml y z ; if isempty(firstCall),firstCall=1;end;
if isempty(ampl), ampl=0; end;
if isempty(besj1result), besj1result=0; end;
if isempty(bj1cs), bj1cs=zeros(1,12); end;
if isempty(bm1cs), bm1cs=zeros(1,21); end;
if isempty(bth1cs), bth1cs=zeros(1,24); end;
if isempty(pi4), pi4=0; end;
if isempty(theta), theta=0; end;
if isempty(xmax), xmax=0; end;
if isempty(xmin), xmin=0; end;
if isempty(xsml), xsml=0; end;
if isempty(y), y=0; end;
if isempty(z), z=0; end;
if isempty(ntj1), ntj1=0; end;
if isempty(ntm1), ntm1=0; end;
if isempty(ntth1), ntth1=0; end;
%***BEGIN PROLOGUE BESJ1
%***PURPOSE Compute the Bessel function of the first kind of order one.
%***LIBRARY SLATEC (FNLIB)
%***CATEGORY C10A1
%***TYPE SINGLE PRECISION (BESJ1-S, DBESJ1-D)
%***KEYWORDS BESSEL FUNCTION, FIRST KIND, FNLIB, ORDER ONE,
% SPECIAL FUNCTIONS
%***AUTHOR Fullerton, W., (LANL)
%***DESCRIPTION
%
% BESJ1(X) calculates the Bessel function of the first kind of
% order one for real argument X.
%
% Series for BJ1 on the interval 0. to 1.60000D+01
% with weighted error 4.48E-17
% log weighted error 16.35
% significant figures required 15.77
% decimal places required 16.89
%
% Series for BM1 on the interval 0. to 6.25000D-02
% with weighted error 5.61E-17
% log weighted error 16.25
% significant figures required 14.97
% decimal places required 16.91
%
% Series for BTH1 on the interval 0. to 6.25000D-02
% with weighted error 4.10E-17
% log weighted error 16.39
% significant figures required 15.96
% decimal places required 17.08
%
%***REFERENCES (NONE)
%***ROUTINES CALLED CSEVL, INITS, R1MACH, XERMSG
%***REVISION HISTORY (YYMMDD)
% 780601 DATE WRITTEN
% 890210 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)
%***end PROLOGUE BESJ1
if isempty(first), first=false; end;
if firstCall, bj1cs(1)=[-.11726141513332787e0]; end;
if firstCall, bj1cs(2)=[-.25361521830790640e0]; end;
if firstCall, bj1cs(3)=[.050127080984469569e0]; end;
if firstCall, bj1cs(4)=[-.004631514809625081e0]; end;
if firstCall, bj1cs(5)=[.000247996229415914e0]; end;
if firstCall, bj1cs(6)=[-.000008678948686278e0]; end;
if firstCall, bj1cs(7)=[.000000214293917143e0]; end;
if firstCall, bj1cs(8)=[-.000000003936093079e0]; end;
if firstCall, bj1cs(9)=[.000000000055911823e0]; end;
if firstCall, bj1cs(10)=[-.000000000000632761e0]; end;
if firstCall, bj1cs(11)=[.000000000000005840e0]; end;
if firstCall, bj1cs(12)=[-.000000000000000044e0]; end;
if firstCall, bm1cs(1)=[.1047362510931285e0]; end;
if firstCall, bm1cs(2)=[.00442443893702345e0]; end;
if firstCall, bm1cs(3)=[-.00005661639504035e0]; end;
if firstCall, bm1cs(4)=[.00000231349417339e0]; end;
if firstCall, bm1cs(5)=[-.00000017377182007e0]; end;
if firstCall, bm1cs(6)=[.00000001893209930e0]; end;
if firstCall, bm1cs(7)=[-.00000000265416023e0]; end;
if firstCall, bm1cs(8)=[.00000000044740209e0]; end;
if firstCall, bm1cs(9)=[-.00000000008691795e0]; end;
if firstCall, bm1cs(10)=[.00000000001891492e0]; end;
if firstCall, bm1cs(11)=[-.00000000000451884e0]; end;
if firstCall, bm1cs(12)=[.00000000000116765e0]; end;
if firstCall, bm1cs(13)=[-.00000000000032265e0]; end;
if firstCall, bm1cs(14)=[.00000000000009450e0]; end;
if firstCall, bm1cs(15)=[-.00000000000002913e0]; end;
if firstCall, bm1cs(16)=[.00000000000000939e0]; end;
if firstCall, bm1cs(17)=[-.00000000000000315e0]; end;
if firstCall, bm1cs(18)=[.00000000000000109e0]; end;
if firstCall, bm1cs(19)=[-.00000000000000039e0]; end;
if firstCall, bm1cs(20)=[.00000000000000014e0]; end;
if firstCall, bm1cs(21)=[-.00000000000000005e0]; end;
if firstCall, bth1cs(1)=[.74060141026313850e0]; end;
if firstCall, bth1cs(2)=[-.004571755659637690e0]; end;
if firstCall, bth1cs(3)=[.000119818510964326e0]; end;
if firstCall, bth1cs(4)=[-.000006964561891648e0]; end;
if firstCall, bth1cs(5)=[.000000655495621447e0]; end;
if firstCall, bth1cs(6)=[-.000000084066228945e0]; end;
if firstCall, bth1cs(7)=[.000000013376886564e0]; end;
if firstCall, bth1cs(8)=[-.000000002499565654e0]; end;
if firstCall, bth1cs(9)=[.000000000529495100e0]; end;
if firstCall, bth1cs(10)=[-.000000000124135944e0]; end;
if firstCall, bth1cs(11)=[.000000000031656485e0]; end;
if firstCall, bth1cs(12)=[-.000000000008668640e0]; end;
if firstCall, bth1cs(13)=[.000000000002523758e0]; end;
if firstCall, bth1cs(14)=[-.000000000000775085e0]; end;
if firstCall, bth1cs(15)=[.000000000000249527e0]; end;
if firstCall, bth1cs(16)=[-.000000000000083773e0]; end;
if firstCall, bth1cs(17)=[.000000000000029205e0]; end;
if firstCall, bth1cs(18)=[-.000000000000010534e0]; end;
if firstCall, bth1cs(19)=[.000000000000003919e0]; end;
if firstCall, bth1cs(20)=[-.000000000000001500e0]; end;
if firstCall, bth1cs(21)=[.000000000000000589e0]; end;
if firstCall, bth1cs(22)=[-.000000000000000237e0]; end;
if firstCall, bth1cs(23)=[.000000000000000097e0]; end;
if firstCall, bth1cs(24)=[-.000000000000000040e0]; end;
if firstCall, pi4=[0.78539816339744831e0]; end;
if firstCall, first=[true]; end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT BESJ1
if( first )
[ntj1 ,bj1cs]=inits(bj1cs,12,0.1.*r1mach(3));
[ntm1 ,bm1cs]=inits(bm1cs,21,0.1.*r1mach(3));
[ntth1 ,bth1cs]=inits(bth1cs,24,0.1.*r1mach(3));
%
xsml = sqrt(8.0.*r1mach(3));
xmin = 2.0.*r1mach(1);
xmax = 1.0./r1mach(4);
end;
first = false;
%
y = abs(x);
if( y>4.0 )
%
if( y>xmax )
xermsg('SLATEC','BESJ1','NO PRECISION BECAUSE ABS(X) IS TOO BIG',2,2);
end;
z = 32.0./y.^2 - 1.0;
ampl =(0.75+csevl(z,bm1cs,ntm1))./sqrt(y);
theta = y - 3.0.*pi4 + csevl(z,bth1cs,ntth1)./y;
besj1result = (abs(ampl).*sign(x)).*cos(theta);
else;
%
besj1result = 0.;
if( y==0.0 )
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','BESJ1','ABS(X) SO SMALL J1 UNDERFLOWS',1,1);
end;
if( y>xmin )
besj1result = 0.5.*x;
end;
if( y>xsml )
besj1result = x.*(.25+csevl(.125.*y.*y-1.,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 BESJ
