function [besy1result,x]=besy1(x);
besy1result=[];
persistent ampl besy1 bm1cs bth1cs by1cs first firstCall ntm1 ntth1 nty1 pi4 theta twodpi xmax xmin xsml y z ; if isempty(firstCall),firstCall=1;end;
if isempty(ampl), ampl=0; end;
if isempty(besy1result), besy1result=0; end;
if isempty(bm1cs), bm1cs=zeros(1,21); end;
if isempty(bth1cs), bth1cs=zeros(1,24); end;
if isempty(by1cs), by1cs=zeros(1,14); end;
if isempty(pi4), pi4=0; end;
if isempty(theta), theta=0; end;
if isempty(twodpi), twodpi=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(ntm1), ntm1=0; end;
if isempty(ntth1), ntth1=0; end;
if isempty(nty1), nty1=0; end;
%***BEGIN PROLOGUE BESY1
%***PURPOSE Compute the Bessel function of the second kind of order
% one.
%***LIBRARY SLATEC (FNLIB)
%***CATEGORY C10A1
%***TYPE SINGLE PRECISION (BESY1-S, DBESY1-D)
%***KEYWORDS BESSEL FUNCTION, FNLIB, ORDER ONE, SECOND KIND,
% SPECIAL FUNCTIONS
%***AUTHOR Fullerton, W., (LANL)
%***DESCRIPTION
%
% BESY1(X) calculates the Bessel function of the second kind of
% order one for real argument X.
%
% Series for BY1 on the interval 0. to 1.60000D+01
% with weighted error 1.87E-18
% log weighted error 17.73
% significant figures required 17.83
% decimal places required 18.30
%
% 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 BESJ1, CSEVL, INITS, R1MACH, XERMSG
%***REVISION HISTORY (YYMMDD)
% 770401 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)
%***end PROLOGUE BESY1
if isempty(first), first=false; end;
if firstCall, by1cs(1)=[.03208047100611908629e0]; end;
if firstCall, by1cs(2)=[1.262707897433500450e0]; end;
if firstCall, by1cs(3)=[.00649996189992317500e0]; end;
if firstCall, by1cs(4)=[-.08936164528860504117e0]; end;
if firstCall, by1cs(5)=[.01325088122175709545e0]; end;
if firstCall, by1cs(6)=[-.00089790591196483523e0]; end;
if firstCall, by1cs(7)=[.00003647361487958306e0]; end;
if firstCall, by1cs(8)=[-.00000100137438166600e0]; end;
if firstCall, by1cs(9)=[.00000001994539657390e0]; end;
if firstCall, by1cs(10)=[-.00000000030230656018e0]; end;
if firstCall, by1cs(11)=[.00000000000360987815e0]; end;
if firstCall, by1cs(12)=[-.00000000000003487488e0]; end;
if firstCall, by1cs(13)=[.00000000000000027838e0]; end;
if firstCall, by1cs(14)=[-.00000000000000000186e0]; 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, twodpi=[0.63661977236758134e0]; end;
if firstCall, pi4=[0.78539816339744831e0]; end;
if firstCall, first=[true]; end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT BESY1
if( first )
[nty1 ,by1cs]=inits(by1cs,14,0.1.*r1mach(3));
[ntm1 ,bm1cs]=inits(bm1cs,21,0.1.*r1mach(3));
[ntth1 ,bth1cs]=inits(bth1cs,24,0.1.*r1mach(3));
%
xmin = 1.571.*exp(max(log(r1mach(1)),-log(r1mach(2)))+.01);
xsml = sqrt(4.0.*r1mach(3));
xmax = 1.0./r1mach(4);
end;
first = false;
%
if( x<=0. )
xermsg('SLATEC','BESY1','X IS ZERO OR NEGATIVE',1,2);
end;
if( x>4.0 )
%
if( x>xmax )
xermsg('SLATEC','BESY1','NO PRECISION BECAUSE X IS BIG',2,2);
end;
%
z = 32.0./x.^2 - 1.0;
ampl =(0.75+csevl(z,bm1cs,ntm1))./sqrt(x);
theta = x - 3.0.*pi4 + csevl(z,bth1cs,ntth1)./x;
besy1result = ampl.*sin(theta);
else;
%
if( x<xmin )
xermsg('SLATEC','BESY1','X SO SMALL Y1 OVERFLOWS',3,2);
end;
y = 0.;
if( x>xsml )
y = x.*x;
end;
besy1result = twodpi.*log(0.5.*x).*besj1(x)+(0.5+csevl(.125.*y-1.,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 BESY
