# bessel function of the first kind

[wikipedia: bessel function](https://en.wikipedia.org/wiki/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

# computation
* [libc](https://www.gnu.org/software/libc/manual/html_node/Special-Functions.html) ([standard](https://www.ibm.com/docs/en/i/7.5?topic=extensions-standard-c-library-functions-table-by-name)) has jn(n, x) and more
* [gsl](https://www.gnu.org/software/gsl/doc/html/specfunc.html#bessel-functions) has gsl_sf_bessel_Jn(n, x)
* [implementation in gsl](https://www.gnu.org/software/gsl/#downloading) see specfunc/bessel_Jn.c
* [implementation of jn in glibc](https://sourceware.org/git/?p=glibc.git;a=tree;f=sysdeps/ieee754) see e_jnl.c in the directories sysdeps/ieee754/dbl-64 and similar (128 bit implementation in ldbl-128)
* [j. bremer, 2017, an algorithm for the rapid numerical evaluation of bessel functions of real orders and arguments](https://arxiv.org/abs/1705.07820) ([code](https://github.com/JamesCBremerJr/BesselEval))
* [leung, po kin; gammie, charles f.; noble, scott c. - bessel: fast bessel function jn(z) routine for large n,z](https://ascl.net/1306.013) ([code](https://rainman.astro.illinois.edu/codelib/))
* [j. harrison, 2009, fast and accurate bessel function computation](https://www.cl.cam.ac.uk/~jrh13/papers/bessel.pdf)
* [digital library of mathematical functions: methods of computation](https://dlmf.nist.gov/10.74)
* [bessel programs in c/c++](http://jean-pierre.moreau.pagesperso-orange.fr/c_bessel.html)
* [bessel.c](https://www.astro.rug.nl/~gipsy/sub/bessel.c) (1998) approximations based on ideas from "Numerical Recipes, Press et. al."

# 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](https://www.bragitoff.com/2017/08/bessel-function-series-c-program/)

## 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](https://github.com/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
~~~