jn function

double jn(int n, double x)

Implementation

double jn(int n, double x) {
  if (n < 0) return (n % 2 == 0 ? 1.0 : -1.0) * jn(-n, x);
  if (n == 0) return j0(x);
  if (n == 1) return j1(x);
  if (x == 0.0) return 0.0;

  if (x.abs() > n.toDouble()) {
    // Forward recurrence stable when |x| > n
    double jPrev = j0(x);
    double jCurr = j1(x);
    for (int k = 1; k < n; k++) {
      double jNext = (2 * k / x) * jCurr - jPrev;
      jPrev = jCurr;
      jCurr = jNext;
    }
    return jCurr;
  } else {
    // Miller's backward recurrence
    int m = n + 20 + (n ~/ 2);
    if (m < 50) m = 50;
    double jPrev = 0.0;
    double jCurr = 1.0;
    double scale = 0.0;
    double result = 0.0;
    for (int k = m; k >= 0; k--) {
      double jNext = (2 * (k + 1) / x) * jCurr - jPrev;
      jPrev = jCurr;
      jCurr = jNext;
      if (k == n) result = jCurr;
      if (k == 0) scale = jCurr;
    }
    return result * j0(x) / scale;
  }
}