In fact, floating point values of x
close to 0
also return Nan
. You actually have three cases to possibly worry about. x = ±∞ results in NaN
as well. Here's a vectorized function that handles those:
function js = sphbesselj(nu,x)
if isscalar(nu) && isscalar(x)
js = 0;
elseif isscalar(nu)
js = zeros(size(x));
nu = js+nu;
else
js = zeros(size(nu));
x = js+x;
end
x0 = (abs(x) < realmin);
x0nult0 = (x0 & nu < 0);
x0nueq0 = (x0 & nu == 0);
js(x0nult0) = Inf; % Re(Nu) < 0, X == 0
js(x0nueq0) = 1; % Re(Nu) == 0, X == 0
i = ~x0nult0 & ~x0nueq0 & ~(x0 & nu > 0) & (abs(x) < realmax);
js(i) = sign(x(i)).*sqrt(pi./(2*x(i))).*besselj(nu(i)+0.5,x(i));
end
A useful resource when developing such functions is http://functions.wolfram.com. The page on the spherical Bessel function of the first kind has many useful relationships.