void powell_eyink(double* B, double divB) {
  double alpha = -divB / (2 * dot_product(B, B));
  for (int i = 0; i < 3; i++) {
    B[i] += alpha * B[i];
  }
}