j0 function

double j0(double x)

Implementation

double j0(double x) {
  x = x.abs();
  if (x < 8.0) {
    double sum = 0, term = 1;
    for (int k = 0; k < 60; k++) {
      sum += term;
      term *= -(x * x) / (4.0 * (k + 1) * (k + 1));
      if (term.abs() < 1e-17 * sum.abs() && k > 0) break;
    }
    return sum;
  }
  double theta = x - M_PI_4;
  return math.sqrt(M_2_PI / x) * math.cos(theta);
}