ActiveThorns = "
    AMReX
    IOUtil
    WaveToyAMReX
"
 
$nlevels = 7
$ncells = 32
Cactus::cctk_show_schedule = no
# Cactus::terminate = "iteration"
# Cactus::cctk_itlast = 64   #TODO 0
Cactus::terminate = "time"
Cactus::cctk_final_time = 50.0
# AMReX::verbose = yes
AMReX::xmin = -100.0
AMReX::ymin = -100.0
AMReX::zmin = -100.0
AMReX::xmax = +100.0
AMReX::ymax = +100.0
AMReX::zmax = +100.0
AMReX::ncells_x = $ncells
AMReX::ncells_y = $ncells
AMReX::ncells_z = $ncells
# AMReX::periodic_x = no
# AMReX::periodic_y = no
# AMReX::periodic_z = no
AMReX::max_num_levels = $nlevels
AMReX::regrid_every = 16
AMReX::regrid_error_threshold = 0.1
AMReX::dtfac = 0.5
WaveToyAMReX::initial_condition = "central potential"
WaveToyAMReX::boundary_condition = "initial"
WaveToyAMReX::central_point_charge = 1.0
WaveToyAMReX::central_point_radius = 1.0
WaveToyAMReX::particle_charge = 1.0
IO::out_dir = $parfile
IO::out_every = $ncells * 2 ** ($nlevels - 1) / 32

CCTK_REAL particle_charge "Particle charge"
{
  *:* :: ""
} 0.0

SCHEDULE WaveToyAMReX_InitializeParticle AT initial AFTER WaveToyAMReX_Initialize
{
  OPTIONS: global
  LANG: C
} "Set up initial conditions for the scalar particle"

SCHEDULE WaveToyAMReX_EvolveParticle AT evol AFTER WaveToyAMReX_Evolve
{
  OPTIONS: global
  LANG: C
} "Evolve the scalar particle"

SCHEDULE WaveToyAMReX_OutputParticle AT analysis
{
  OPTIONS: global
  LANG: C
} "Output for the scalar particle"

  };
}
// Gradient
template <typename T> auto xderiv(T f(T t, T x, T y, T z), T dx) {
  return [=](T t, T x, T y, T z) {
    return (f(t, x + dx, y, z) - f(t, x - dx, y, z)) / (2 * dx);
  };
}
template <typename T> auto yderiv(T f(T t, T x, T y, T z), T dy) {
  return [=](T t, T x, T y, T z) {
    return (f(t, x, y + dy, z) - f(t, x, y - dy, z)) / (2 * dy);

template <typename T> auto zderiv(T f(T t, T x, T y, T z), T dz) {
  return [=](T t, T x, T y, T z) {
    return (f(t, x, y, z + dz) - f(t, x, y, z - dz)) / (2 * dz);
  };
}

  return central_point_charge * pow(central_point_radius, -dim) *

  return -central_point_charge * pow(central_point_radius, -dim) *

template <typename T> struct particle_t {
  T mass;
  T time;
  array<T, dim> pos;
  array<T, dim> vel;
};
particle_t<CCTK_REAL> particle;
////////////////////////////////////////////////////////////////////////////////

}
extern "C" void WaveToyAMReX_InitializeParticle(CCTK_ARGUMENTS) {
  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;
  if (particle_charge == 0)
    return;
  particle.mass = 1.0;
  particle.time = 1.0;
  for (int d = 0; d < dim; ++d)
    particle.pos[d] = 0.0;
  for (int d = 0; d < dim; ++d)
    particle.vel[d] = 0.0;
  particle.pos[0] = 10.0; // x axis
  particle.vel[1] = -0.1; // y direction, counter-clockwise

  const CCTK_REAL t = cctk_time;

    CCTK_REAL r = sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2));

               dt * (ddx_phi + ddy_phi + ddz_phi - pow(mass, 2) * phi_p[idx] -
                     4 * M_PI * central_point_charge *
                         pow(central_point_radius, -dim) *
                         spline(r / central_point_radius));

               dt * (ddx_phi + ddy_phi + ddz_phi - pow(mass, 2) * phi_p[idx] +
                     /*TODO*/ 0 * 4 * M_PI * central_potential(t, x, y, z));

}
extern "C" void WaveToyAMReX_EvolveParticle(CCTK_ARGUMENTS) {
  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;
  if (particle_charge == 0)
    return;
  int type;
  const void *maxlevel_par =
      CCTK_ParameterGet("max_num_levels", "AMReX", &type);
  assert(type == PARAMETER_INT);
  const int maxlevel = *static_cast<const CCTK_INT *>(maxlevel_par);
  const int maxlevfac = 1 << (maxlevel - 1);
  const CCTK_REAL dt = CCTK_DELTA_TIME;
  const CCTK_REAL dx[dim] = {cctk_delta_space[0] / maxlevfac,
                             cctk_delta_space[1] / maxlevfac,
                             cctk_delta_space[2] / maxlevfac};
  const auto q = particle_charge;
  const auto m = particle.mass;
  const auto t = particle.time;
  const auto p = particle.pos;
  const auto v = particle.vel;
  // http://www.livingreviews.org/lrr-2011-7
  const CCTK_REAL v2 = pow(v[0], 2) + pow(v[1], 2) + pow(v[2], 2);
  // vt = dt/dtau
  const CCTK_REAL vt = sqrt(1 - v2);
  const CCTK_REAL phit = timederiv(central_potential, dt)(t, p[0], p[1], p[2]);
  const array<CCTK_REAL, dim> dphi = {
      xderiv(central_potential, dx[0])(t, p[0], p[1], p[2]),
      yderiv(central_potential, dx[1])(t, p[0], p[1], p[2]),
      zderiv(central_potential, dx[2])(t, p[0], p[1], p[2])};
  // (17.8)
  // dm/dtau = - q phi,a u^a
  // alternatively: m = m0 - q phi
  particle.mass +=
      -dt / vt * q *
      (-phit * vt + dphi[0] * v[0] + dphi[1] * v[1] + dphi[2] * v[2]);
  particle.time += dt / vt;
  for (int a = 0; a < dim; ++a)
    particle.pos[a] += dt * v[a];
  // (17.7)
  // m du/dtau = q (g^ab + u^a u^b) phi,b
  for (int a = 0; a < dim; ++a) {
    CCTK_REAL dva = v[a] * vt * phit;
    for (int b = 0; b < dim; ++b)
      dva += ((a == b) + v[a] * v[b]) * dphi[b];
    particle.vel[a] += dt / vt * q / m * dva;
  }

extern "C" void WaveToyAMReX_OutputParticle(CCTK_ARGUMENTS) {
  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;

  if (particle_charge == 0)
    return;
  const auto q = particle_charge;
  const auto m = particle.mass;
  const auto t = particle.time;
  const auto p = particle.pos;
  const auto v = particle.vel;
  const CCTK_REAL r = sqrt(pow(p[0], 2) + pow(p[1], 2) + pow(p[2], 2));
  const CCTK_REAL v2 = pow(v[0], 2) + pow(v[1], 2) + pow(v[2], 2);
  const CCTK_REAL vt = sqrt(1 - v2);
  CCTK_VINFO("Particle:");
  cout << "  q=" << q << " m=" << m << "\n"
       << "  t=" << t << " x=[" << p[0] << "," << p[1] << "," << p[2]
       << "] r=" << r << "\n"
       << "  vt=" << vt << " v=[" << v[0] << "," << v[1] << "," << v[2]
       << "] v2=" << v2 << "\n";
}