#include <math.h>

#include <cctk.h>
#include <cctk_Arguments.h>
#include <cctk_Parameters.h>

CCTK_INT find_closest(const cGH *cctkGH, CCTK_INT *cctk_lsh,
                      CCTK_REAL *x, CCTK_REAL x_min, CCTK_INT dir)
{
    CCTK_INT i, ijk, min_i = -1;
    CCTK_REAL min = 1.e100;
    for(i=3; i<cctk_lsh[dir]-3; i++)
    {
        ijk = CCTK_GFINDEX3D(cctkGH, (dir==0)?i:0, (dir==1)?i:0, (dir==2)?i:0);
        if (fabs(x[ijk] - x_min) < min)
        {
            min = fabs(x[ijk] - x_min);
            min_i = i;
        }
    }
    return min_i;
}

void ADMMass_Surface(CCTK_ARGUMENTS)
{
    DECLARE_CCTK_ARGUMENTS
    DECLARE_CCTK_PARAMETERS

    CCTK_INT i,j,k, ijk;
    CCTK_REAL x_min, x_max, y_min, y_max, z_min, z_max;
    CCTK_INT  x_min_i, x_max_i, y_min_j, y_max_j, z_min_k, z_max_k, min;
    CCTK_REAL detg, idetg, ds[3];
    CCTK_REAL u[3][3], dg[3][3][3];
    CCTK_REAL local_sum = 0.0;

    /* grid-function strides for ADMMacros */
    const CCTK_INT di = 1;
    const CCTK_INT dj = cctk_lsh[0];
    const CCTK_INT dk = cctk_lsh[0]*cctk_lsh[1];

    if (ADMMass_boundary_distance[*ADMMass_LoopCounter] > 0.0)
    {
        x_min = ADMMass_x_pos[*ADMMass_LoopCounter] -
                ADMMass_boundary_distance[*ADMMass_LoopCounter];
        x_max = ADMMass_x_pos[*ADMMass_LoopCounter] +
                ADMMass_boundary_distance[*ADMMass_LoopCounter];
        y_min = ADMMass_y_pos[*ADMMass_LoopCounter] -
                ADMMass_boundary_distance[*ADMMass_LoopCounter];
        y_max = ADMMass_y_pos[*ADMMass_LoopCounter] +
                ADMMass_boundary_distance[*ADMMass_LoopCounter];
        z_min = ADMMass_z_pos[*ADMMass_LoopCounter] -
                ADMMass_boundary_distance[*ADMMass_LoopCounter];
        z_max = ADMMass_z_pos[*ADMMass_LoopCounter] +
                ADMMass_boundary_distance[*ADMMass_LoopCounter];
    }
    else
    {
        x_min = ADMMass_x_min[*ADMMass_LoopCounter];
        x_max = ADMMass_x_max[*ADMMass_LoopCounter];
        y_min = ADMMass_y_min[*ADMMass_LoopCounter];
        y_max = ADMMass_y_max[*ADMMass_LoopCounter];
        z_min = ADMMass_z_min[*ADMMass_LoopCounter];
        z_max = ADMMass_z_max[*ADMMass_LoopCounter];
    }
    x_min_i = find_closest(cctkGH, cctk_lsh, x, x_min, 0);
    x_max_i = find_closest(cctkGH, cctk_lsh, x, x_max, 0);
    y_min_j = find_closest(cctkGH, cctk_lsh, y, y_min, 1);
    y_max_j = find_closest(cctkGH, cctk_lsh, y, y_max, 1);
    z_min_k = find_closest(cctkGH, cctk_lsh, z, z_min, 2);
    z_max_k = find_closest(cctkGH, cctk_lsh, z, z_max, 2);

#include "CactusEinstein/ADMMacros/src/macro/UPPERMET_declare.h"
#include "CactusEinstein/ADMMacros/src/macro/DG_declare.h"
#include "CactusEinstein/ADMMacros/src/macro/ADM_Spacing_declare.h"

    for(i=3; i<cctk_lsh[0]-3; i++)
     for(j=3; j<cctk_lsh[1]-3; j++)
      for(k=3; k<cctk_lsh[2]-3; k++)
      {
          ijk = CCTK_GFINDEX3D(cctkGH, i, j, k);
          ADMMass_SurfaceMass_GF[ijk] = 0.0;
          if ((i >= x_min_i) &&
              (i <= x_max_i) &&
              (j >= y_min_j) &&
              (j <= y_max_j) &&
              (k >= z_min_k) &&
              (k <= z_max_k))
          {
              ds[0] = 0.0;
              ds[1] = 0.0;
              ds[2] = 0.0;
#include "CactusEinstein/ADMMacros/src/macro/ADM_Spacing.h"
              if ((i == x_min_i) && (j != y_max_j) && (k != z_max_k))
                  ds[0] = -dy*dz;
              if ((i == x_max_i) && (j != y_max_j) && (k != z_max_k))
                  ds[0] = dy*dz;
              if ((j == y_min_j) && (i != x_max_i) && (k != z_max_k))
                  ds[1] = -dx*dz;
              if ((j == y_max_j) && (i != x_max_i) && (k != z_max_k))
                  ds[1] = dx*dz;
              if ((k == z_min_k) && (j != y_max_j) && (i != x_max_i))
                  ds[2] = -dx*dy;
              if ((k == z_max_k) && (j != y_max_j) && (i != x_max_i))
                  ds[2] = dx*dy;
              if ((ds[0] != 0.0) || (ds[1] != 0.0) || (ds[2] != 0.0))
              {
#include "CactusEinstein/ADMMacros/src/macro/UPPERMET_guts.h"
#include "CactusEinstein/ADMMacros/src/macro/DG_guts.h"
                  u[0][0] = UPPERMET_UXX;
                  u[0][1] = UPPERMET_UXY;
                  u[0][2] = UPPERMET_UXZ;
                  u[1][0] = UPPERMET_UXY;
                  u[1][1] = UPPERMET_UYY;
                  u[1][2] = UPPERMET_UYZ;
                  u[2][0] = UPPERMET_UXZ;
                  u[2][1] = UPPERMET_UYZ;
                  u[2][2] = UPPERMET_UZZ;
                  dg[0][0][0] = DXDG_DXDGXX;
                  dg[0][0][1] = DYDG_DYDGXX;
                  dg[0][0][2] = DZDG_DZDGXX;
                  dg[0][1][0] = DXDG_DXDGXY;
                  dg[0][1][1] = DYDG_DYDGXY;
                  dg[0][1][2] = DZDG_DZDGXY;
                  dg[0][2][0] = DXDG_DXDGXZ;
                  dg[0][2][1] = DYDG_DYDGXZ;
                  dg[0][2][2] = DZDG_DZDGXZ;
                  dg[1][0][0] = DXDG_DXDGXY;
                  dg[1][0][1] = DYDG_DYDGXY;
                  dg[1][0][2] = DZDG_DZDGXY;
                  dg[1][1][0] = DXDG_DXDGYY;
                  dg[1][1][1] = DYDG_DYDGYY;
                  dg[1][1][2] = DZDG_DZDGYY;
                  dg[1][2][0] = DXDG_DXDGYZ;
                  dg[1][2][1] = DYDG_DYDGYZ;
                  dg[1][2][2] = DZDG_DZDGYZ;
                  dg[2][0][0] = DXDG_DXDGXZ;
                  dg[2][0][1] = DYDG_DYDGXZ;
                  dg[2][0][2] = DZDG_DZDGXZ;
                  dg[2][1][0] = DXDG_DXDGYZ;
                  dg[2][1][1] = DYDG_DYDGYZ;
                  dg[2][1][2] = DZDG_DZDGYZ;
                  dg[2][2][0] = DXDG_DXDGZZ;
                  dg[2][2][1] = DYDG_DYDGZZ;
                  dg[2][2][2] = DZDG_DZDGZZ;
                  for (int ti = 0; ti < 3; ti++)
                   for (int tj = 0; tj < 3; tj++)
                    for (int tk = 0; tk < 3; tk++)
                     for (int tl = 0; tl < 3; tl++)
                         ADMMass_SurfaceMass_GF[ijk] +=
                             u[ti][tj] * u[tk][tl] *
                             ( dg[ti][tk][tj] - dg[ti][tj][tk] ) * ds[tk];
                  ADMMass_SurfaceMass_GF[ijk] *= sqrt(DETG_DETG) / 16 / 3.14159265;
                  local_sum += ADMMass_SurfaceMass_GF[ijk];
              }
          }
      }
    #include "CactusEinstein/ADMMacros/src/macro/UPPERMET_undefine.h"
    #include "CactusEinstein/ADMMacros/src/macro/DG_undefine.h"
    #include "CactusEinstein/ADMMacros/src/macro/ADM_Spacing_undefine.h"
    /*CCTK_VInfo(CCTK_THORNSTRING,
     *           "ADM mass local to one processor (surface): %g\n", local_sum);
     */
}

void ADMMass_Surface_Global(CCTK_ARGUMENTS)
{
    DECLARE_CCTK_ARGUMENTS
    DECLARE_CCTK_PARAMETERS

    CCTK_INT reduction_handle;

    reduction_handle = CCTK_ReductionHandle("sum");
    if (reduction_handle < 0)
        CCTK_WARN(0, "Unable to ge reduction handle.");

    if (CCTK_Reduce(cctkGH, -1, reduction_handle, 1,
                    CCTK_VARIABLE_REAL,
                    &ADMMass_SurfaceMass[*ADMMass_LoopCounter], 1,
                    CCTK_VarIndex("ADMMass::ADMMass_SurfaceMass_GF")))
        CCTK_WARN(0, "Error while reducing ADMMass_SurfaceMass_GF");
    
    CCTK_VInfo(CCTK_THORNSTRING,
               "ADM mass, surface, detector %d: %g %g\n", *ADMMass_LoopCounter,
               ADMMass_boundary_distance[*ADMMass_LoopCounter],
               ADMMass_SurfaceMass[*ADMMass_LoopCounter]);
}

void ADMMass_Surface_Lapse(CCTK_ARGUMENTS)
{
    DECLARE_CCTK_ARGUMENTS
    DECLARE_CCTK_PARAMETERS

    CCTK_INT i,j,k, ijk;
    CCTK_REAL local_sum = 0.0;

    for(i=3; i<cctk_lsh[0]-3; i++)
     for(j=3; j<cctk_lsh[1]-3; j++)
      for(k=3; k<cctk_lsh[2]-3; k++)
      {
          ijk = CCTK_GFINDEX3D(cctkGH, i, j, k);
          ADMMass_SurfaceMass_GF[ijk] *= alp[ijk];
          local_sum += ADMMass_SurfaceMass_GF[ijk];
      }
}

void ADMMass_Surface_Lapse_Global(CCTK_ARGUMENTS)
{
    DECLARE_CCTK_ARGUMENTS
    DECLARE_CCTK_PARAMETERS

    CCTK_INT reduction_handle;

    reduction_handle = CCTK_ReductionHandle("sum");
    if (reduction_handle < 0)
        CCTK_WARN(0, "Unable to ge reduction handle.");

    if (CCTK_Reduce(cctkGH, -1, reduction_handle, 1,
                    CCTK_VARIABLE_REAL,
                    &ADMMass_SurfaceMass_Lapse[*ADMMass_LoopCounter], 1,
                    CCTK_VarIndex("ADMMass::ADMMass_SurfaceMass_GF")))
        CCTK_WARN(0, "Error while reducing ADMMass_SurfaceMass_GF");
    
    CCTK_VInfo(CCTK_THORNSTRING,
           "ADM mass, surface lapse, detector %d: %g %g\n",
           *ADMMass_LoopCounter,
           ADMMass_boundary_distance[*ADMMass_LoopCounter],
           ADMMass_SurfaceMass_Lapse[*ADMMass_LoopCounter]);
}