2023-01-05

# bessel function of the first kind

wikipedia: bessel function

there are multiple kinds of bessel functions that are grouped together just because they are similar. they are called first kind, second kind, and so on.

x values for which the bessel function of the first kind crosses zero:

• order 0: 2.404825557695773, 5.520078110286311, 8.653727912911013
• order 1: 3.8317059702075125, 7.015586669815619, 10.173468135062722
• order 2: 5.135622301840682556301, 8.417244140399864857784, 11.61984117214905942709
• order 3: 6.380161895923983506237, 9.761023129981669678545, 13.01520072169843441983

# formula

```bessel(order, x) = sum(n, 0, infinity,
((-1) ** n / (factorial(n) * gamma(n + order + 1)))
* (x / 2) ** (2 * n + order))```
• for integers, the gamma function can be replaced with the factorial
• the most direct translation to programming code quickly reaches the limits of numerical representation because of the exponential (2 * n + order)

# version for small values

this version is simple and compact, but in software it only gives accurate values for about x < 36 and then starts to rapidly move towards infinity.

## formula

```bessel(order, x) = (x / 2) ** order
* sum(n, 0, infinity, (-1 ** n * (x ** 2 / 4) ** n) / (n! * (n + order)!))```

the ratio between consecutively summed terms is

`-1 * (x ** 2 / 4) / (n * (order + n))`

source

## scheme

```(define (factorial n) (if (> 1 n) 1 (* n (factorial (- n 1)))))

(define (bessel order x term-count)
(* (expt (/ x 2) order)
(let loop ((n 1) (term (/ 1 (factorial order))))
(if (< n term-count)
(+ term (loop (+ 1 n) (* term -1 (/ (/ (* x x) 4) (* n (+ order n))))))
term))))```

## coffeescript

```factorial = (n) -> if 1 > n then 1 else n * factorial(n - 1)

bessel = (order, x, term_count) ->
term = 1 / factorial order
sum = term
for n in [1...term_count]
term = term * -1 * (x * x / 4) / (n * (order + n))
sum += term
sum * Math.pow(x / 2, order)```

# general version

extracted and translated from SheetJs/bessel

## coffeescript

```besselj = do ->
M = Math
W = 2 / M.PI

b0_a1a = [57568490574.0, -13362590354.0, 651619640.7, -11214424.18, 77392.33017, -184.9052456].reverse()
b0_a2a = [57568490411.0, 1029532985.0, 9494680.718, 59272.64853, 267.8532712, 1.0].reverse()
b0_a1b = [1.0, -0.1098628627e-2, 0.2734510407e-4, -0.2073370639e-5, 0.2093887211e-6].reverse()
b0_a2b = [-0.1562499995e-1, 0.1430488765e-3, -0.6911147651e-5, 0.7621095161e-6, -0.934935152e-7].reverse()
b1_a1a = [72362614232.0, -7895059235.0, 242396853.1, -2972611.439, 15704.48260, -30.16036606].reverse()
b1_a2a = [144725228442.0, 2300535178.0, 18583304.74, 99447.43394, 376.9991397, 1.0].reverse()
b1_a1b = [1.0, 0.183105e-2, -0.3516396496e-4, 0.2457520174e-5, -0.240337019e-6].reverse()
b1_a2b = [0.04687499995, -0.2002690873e-3, 0.8449199096e-5, -0.88228987e-6, 0.105787412e-6].reverse()

horner = (arr, v) ->
i = 0
z = 0
while i < arr.length
z = v * z + arr[i]
i += 1
z

bessel_iter = (x, n, f0, f1, sign) ->
return f0 if n == 0
return f1 if n == 1
tdx = 2 / x
f2 = f1
o = 1
while o < n
f2 = f1 * o * tdx + sign * f0
f0 = f1
f1 = f2
o += 1
f2

bessel0 = (x) ->
a = 0
a1 = 0
a2 = 0
y = x * x
if x < 8
a1 = horner(b0_a1a, y)
a2 = horner(b0_a2a, y)
a = a1 / a2
else
xx = x - 0.785398164
y = 64 / y
a1 = horner(b0_a1b, y)
a2 = horner(b0_a2b, y)
a = M.sqrt(W / x) * (M.cos(xx) * a1 - (M.sin(xx) * a2 * 8 / x))
a

bessel1 = (x) ->
a = 0
a1 = 0
a2 = 0
y = x * x
xx = M.abs(x) - 2.356194491
if M.abs(x) < 8
a1 = x * horner(b1_a1a, y)
a2 = horner(b1_a2a, y)
a = a1 / a2
else
y = 64 / y
a1 = horner(b1_a1b, y)
a2 = horner(b1_a2b, y)
a = M.sqrt(W / M.abs(x)) * (M.cos(xx) * a1 - (M.sin(xx) * a2 * 8 / M.abs(x)))
if x < 0
a = -a
a

(n, x) ->
n = M.round(n)
if !isFinite(x)
return if isNaN(x) then x else 0
if n < 0
return (if n % 2 then -1 else 1) * besselj(x, -n)
if x < 0
return (if n % 2 then -1 else 1) * besselj(-x, n)
return bessel0(x) if n == 0
return bessel1(x) if n == 1
return 0 if x == 0
ret = 0.0
if x > n
ret = bessel_iter(x, n, bessel0(x), bessel1(x), -1)
else
m = 2 * M.floor((n + M.floor(M.sqrt(40 * n))) / 2)
jsum = false
bjp = 0.0
sum = 0.0
bj = 1.0
bjm = 0.0
tox = 2 / x
j = m
while j > 0
bjm = j * tox * bj - bjp
bjp = bj
bj = bjm
if M.abs(bj) > 1e10
bj *= 1e-10
bjp *= 1e-10
ret *= 1e-10
sum *= 1e-10
if jsum
sum += bj
jsum = !jsum
if j == n
ret = bjp
j -= 1
sum = 2.0 * sum - bj
ret /= sum
ret```