REQUIRES Arith

CarpetX::ncells_x = 64
CarpetX::ncells_y = 64
CarpetX::ncells_z = 64

CarpetX::ncells_x = 128   # 64
CarpetX::ncells_y = 128   # 64
CarpetX::ncells_z = 128   # 64

CarpetX::interpolation_order = 3

CarpetX::interpolation_order = 3   # 3

BrillLindquist::x0 = 0.1   # 0.0

AHFinder::lmax = 200

AHFinder::lmax_filter = 9999

AHFinder::z0 = 0.0
# AHFinder::z0 = -0.1
AHFinder::r0 = 0.5
# AHFinder::r0 = 0.6
AHFinder::r1z = 0.0
# AHFinder::r1z = 0.1
AHFinder::maxiters = 20 # 100

AHFinder::z0 = -0.1   # 0.0
AHFinder::r0 = 0.6    # 0.5
AHFinder::r1z = 0.1   # 0.0
AHFinder::maxiters = 100

CCTK_INT lmax "Number of modes for horizon" STEERABLE=always

CCTK_INT npoints "Number of sampling points (per direction) on the sphere" STEERABLE=always

} 16

} 81

CCTK_INT npoints "Number of sampling points (per direction) on the sphere" STEERABLE=always

CCTK_INT lmax_filter "Number of modes for horizon" STEERABLE=always

  1:* :: "must be at least nmodes"
} 76

  1:* :: ""
} 9999

#include <vect.hxx>

using namespace Arith;

#pragma omp simd

template <typename T>
metric_t<T> calculate_metric(const cGH *const cctkGH,
                             const coords_t<T> &coords) {
  const geom_t &geom = coords.geom;
  metric_t<T> metric(geom);
#pragma omp simd
  for (int n = 0; n < geom.npoints; ++n) {
    const T x0 = 0.0;
    const T y0 = 0.0;
    const T z0 = 0.0;
    const T M = 1.0;
    using U = vect<T, 3>;
    using DT = dual<T, U>;
    const DT x{coords.x.data()[n], {1, 0, 0}};
    const DT y{coords.y.data()[n], {0, 1, 0}};
    const DT z{coords.z.data()[n], {0, 0, 1}};
    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;
    metric.gxy.data()[n] = 0;
    metric.gxz.data()[n] = 0;
    metric.gyy.data()[n] = psi.val;
    metric.gyz.data()[n] = 0;
    metric.gzz.data()[n] = psi.val;
    metric.gxx_x.data()[n] = psi.eps[0];
    metric.gxy_x.data()[n] = 0;
    metric.gxz_x.data()[n] = 0;
    metric.gyy_x.data()[n] = psi.eps[0];
    metric.gyz_x.data()[n] = 0;
    metric.gzz_x.data()[n] = psi.eps[0];
    metric.gxx_y.data()[n] = psi.eps[1];
    metric.gxy_y.data()[n] = 0;
    metric.gxz_y.data()[n] = 0;
    metric.gyy_y.data()[n] = psi.eps[1];
    metric.gyz_y.data()[n] = 0;
    metric.gzz_y.data()[n] = psi.eps[1];
    metric.gxx_z.data()[n] = psi.eps[2];
    metric.gxy_z.data()[n] = 0;
    metric.gxz_z.data()[n] = 0;
    metric.gyy_z.data()[n] = psi.eps[2];
    metric.gyz_z.data()[n] = 0;
    metric.gzz_z.data()[n] = psi.eps[2];
    metric.kxx.data()[n] = 0;
    metric.kxy.data()[n] = 0;
    metric.kxz.data()[n] = 0;
    metric.kyy.data()[n] = 0;
    metric.kyz.data()[n] = 0;
    metric.kzz.data()[n] = 0;
  }

  return metric;
}

  T cx, cy, cz;

  T area = 0;

  aij_t<T> cxij(geom), cyij(geom), czij(geom);
  aij_t<T> Aij(geom);

#pragma omp simd

      // const dual<T> x = x0 + r * sin(theta) * cos(phi);
      // const dual<T> y = y0 + r * sin(theta) * sin(phi);
      // const dual<T> z = z0 + r * cos(theta);

      const dual<T> x = x0 + r * sin(theta) * cos(phi);
      const dual<T> y = y0 + r * sin(theta) * sin(phi);
      const dual<T> z = z0 + r * cos(theta);

      const T darea = sin(theta.val) * s_sqrt_detQ * geom.coord_dtheta(i, j) *
                      geom.coord_dphi(i, j);

      // const T darea = sin(theta.val) * s_sqrt_detQ * geom.coord_dtheta(i, j)
      // *
      //                 geom.coord_dphi(i, j);
      const T A = s_sqrt_detQ;

      area += darea;

      Aij(i, j) = A;
      cxij(i, j) = A * x.val;
      cyij(i, j) = A * y.val;
      czij(i, j) = A * z.val;

  const alm_t<T> Alm = expand(Aij);
  const alm_t<T> cxlm = expand(cxij);
  const alm_t<T> cylm = expand(cyij);
  const alm_t<T> czlm = expand(czij);

#pragma omp simd

  const T area = sqrt(4 * T(M_PI)) * real(Alm(0, 0));
  const T cx = real(cxlm(0, 0)) / real(Alm(0, 0));
  const T cy = real(cylm(0, 0)) / real(Alm(0, 0));
  const T cz = real(czlm(0, 0)) / real(Alm(0, 0));

#pragma omp simd

#pragma omp simd

  return {hlm, area, Thetalm, hlm_new};

  return {.hlm = hlm,
          .area = area,
          .cx = cx,
          .cy = cy,
          .cz = cz,
          .Thetalm = Thetalm,
          .hlm_new = hlm_new};

  const auto metric = interpolate_metric(cctkGH, coords);

#warning "TODO"
  // const auto metric = interpolate_metric(cctkGH, coords);
  const auto metric = calculate_metric(cctkGH, coords);

      make_unique<alm_t<T> >(filter(hlm_ini, lmax));

      make_unique<alm_t<T> >(filter(hlm_ini, lmax_filter));

    const auto hlm_new = filter(res.hlm_new, lmax);

    const auto hlm_new = filter(res.hlm_new, lmax_filter);

#pragma omp simd

    T h_maxabs{0}; // ignoring l=0
    for (int l = 1; l <= geom.lmax; ++l)
      for (int m = -l; m <= l; ++m)
        h_maxabs = fmax(h_maxabs, norm(hlm_new(l, m)));
    h_maxabs = sqrt(h_maxabs);

#pragma omp simd

    const T cx = res.cx;
    const T cy = res.cy;
    const T cz = res.cz;

    CCTK_VINFO("iter=%d, h=%g, R=%g |Θ|=%g, |∇h|=%g |Δh|=%g", iter, double(h),
               double(R), double(Theta_maxabs), double(h_maxabs),
               double(dh_maxabs));
    // for (int l = 0; l <= lmax; ++l) {
    //   T h_maxabs{0};
    //   for (int m = -l; m <= l; ++m)
    //     h_maxabs = fmax(h_maxabs, norm(hlm_new(l, m)));
    //   h_maxabs = sqrt(h_maxabs);
    //   CCTK_VINFO("l=%d |h|=%g", l, double(h_maxabs));
    // }
    // for (int l = 0; l <= lmax; ++l) {
    //   printf("l=%d", l);
    //   for (int m = -l; m <= l; ++m)
    //     printf(" %g", double(abs(hlm_new(l, m))));
    //   printf("\n");
    // }

    CCTK_VINFO("iter=%d h=%f c=[%f,%f,%f] R=%f", iter, double(h), double(cx),
               double(cy), double(cz), double(R));
    CCTK_VINFO("  |Θ|=%g |Δh|=%g", double(Theta_maxabs), double(dh_maxabs));

    // alm_t<T> hlm_next(geom);
    // for (int l = 0; l <= geom.lmax; ++l)
    //   for (int m = -l; m <= l; ++m)
    //     hlm_next(l, m) = hlm(l, m) - (hlm_new(l, m) - hlm(l, m));
    // hlm_ptr = make_unique<alm_t<T> >(move(hlm_next));