#include <limits>

  const T fx = 0.8;
  const T fy = 1.0;
  const T fz = 1.0;

    const DT r = sqrt(pow(x - x0, 2) + pow(y - y0, 2) + pow(z - z0, 2));
    const DT psi = pow(1 + M / (2 * r), 4);
    metric.gxx.data()[n] = psi.val;

    const DT X = (x - x0) / fx;
    const DT Y = (y - y0) / fy;
    const DT Z = (z - z0) / fz;
    const DT R = sqrt(pow(X, 2) + pow(Y, 2) + pow(Z, 2));
    const DT Psi = pow(1 + M / (2 * R), 4);
    const vect<T, 3> eta{pow(X.eps[0], 2), pow(Y.eps[1], 2), pow(Z.eps[2], 2)};
    metric.gxx.data()[n] = Psi.val * eta[0];

    metric.gyy.data()[n] = psi.val;

    metric.gyy.data()[n] = Psi.val * eta[1];

    metric.gzz.data()[n] = psi.val;
    metric.gxx_x.data()[n] = psi.eps[0];

    metric.gzz.data()[n] = Psi.val * eta[2];
    metric.gxx_x.data()[n] = Psi.eps[0] * eta[0];

    metric.gyy_x.data()[n] = psi.eps[0];

    metric.gyy_x.data()[n] = Psi.eps[0] * eta[1];

    metric.gzz_x.data()[n] = psi.eps[0];
    metric.gxx_y.data()[n] = psi.eps[1];

    metric.gzz_x.data()[n] = Psi.eps[0] * eta[2];
    metric.gxx_y.data()[n] = Psi.eps[1] * eta[0];

    metric.gyy_y.data()[n] = psi.eps[1];

    metric.gyy_y.data()[n] = Psi.eps[1] * eta[1];

    metric.gzz_y.data()[n] = psi.eps[1];
    metric.gxx_z.data()[n] = psi.eps[2];

    metric.gzz_y.data()[n] = Psi.eps[1] * eta[2];
    metric.gxx_z.data()[n] = Psi.eps[2] * eta[0];

    metric.gyy_z.data()[n] = psi.eps[2];

    metric.gyy_z.data()[n] = Psi.eps[2] * eta[1];

    metric.gzz_z.data()[n] = psi.eps[2];

    metric.gzz_z.data()[n] = Psi.eps[2] * eta[2];

    alm_t<T> dhlm(geom);
    for (int l = 0; l <= geom.lmax; ++l)
#pragma omp simd
      for (int m = -l; m <= l; ++m)
        dhlm(l, m) = hlm_new(l, m) - hlm(l, m);
    auto dhij = evaluate(dhlm);
    aij_t<T> dh2ij(geom);
    for (int i = 0; i < geom.ntheta; ++i)
#pragma omp simd
      for (int j = 0; j < geom.nphi; ++j)
        dh2ij(i, j) = pow(dhij(i, j), 2);
    const auto dh2lm = expand(dh2ij);
    const T dh_norm2 = sqrt(sqrt(4 * M_PI) * real(dh2lm(0, 0)));

    aij_t<T> Theta2ij(geom);
    for (int i = 0; i < geom.ntheta; ++i)
#pragma omp simd
      for (int j = 0; j < geom.nphi; ++j)
        Theta2ij(i, j) = pow(Thetaij(i, j), 2);
    const auto Theta2lm = expand(Theta2ij);
    const T Theta_norm2 = sqrt(sqrt(4 * M_PI) * real(Theta2lm(0, 0)));

    CCTK_VINFO("  |Θ|=%g |Δh|=%g", double(Theta_maxabs), double(dh_maxabs));
    if (iter >= maxiters || dh_maxabs <= 1.0e-12)

    CCTK_VINFO("  h_0m=%f", double(real(hlm(0, 0))));
    CCTK_VINFO("  h_1m=%f (%f,%f)", double(real(hlm(1, 0))),
               double(real(hlm(1, 1))), double(imag(hlm(1, 1))));
    CCTK_VINFO("  h_2m=%f (%f,%f) (%f,%f)", double(real(hlm(2, 0))),
               double(real(hlm(2, 1))), double(imag(hlm(2, 1))),
               double(real(hlm(2, 2))), double(imag(hlm(2, 2))));
    CCTK_VINFO("  |Θ|∞=%g |Θ|2=%g |Δh|∞=%g |Δh|2=%g", double(Theta_maxabs),
               double(Theta_norm2), double(dh_maxabs), double(dh_norm2));
    const T eps = pow(numeric_limits<T>::epsilon(), T(3) / 4);
    if (iter >= maxiters || dh_maxabs <= eps)