#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"
#include "cctk_Functions.h"

#include "loopcontrol.h"
#include "HydroBase.h"

void SetMask_SphericalSurface (CCTK_ARGUMENTS)
{
  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;

  CCTK_INT *mask = (CCTK_INT*) CCTK_VarDataPtr(cctkGH, 0, SetMask_MaskName);
  if (!mask)
    CCTK_WARN(0, "No such variable, or no storage enabled");

  /* Delete mask first!  */
  #pragma omp parallel
      {
        LC_LOOP3(mask_zero, i,j,k, 0,0,0,
                 cctk_lsh[0], cctk_lsh[1], cctk_lsh[2],
                 cctk_lsh[0], cctk_lsh[1], cctk_lsh[2])
        {
          CCTK_INT i3D = CCTK_GFINDEX3D(cctkGH, i, j, k);
          mask[i3D] = HYDRO_EXCISION_NORMAL;
        }
        LC_ENDLOOP3(mask_zero);
      }

  /* Now set excision! */
  for (int smi = 0; smi < 10; smi++)
  {
    CCTK_INT sfi = sf_IdFromName(SetMask_SurfaceIndex[smi], SetMask_SurfaceName[smi]);
    if (sfi >= 0 && sf_active[sfi] && sf_valid[sfi] >= 0)
    {
      #pragma omp parallel
      {
        LC_LOOP3(setsurface, i,j,k, 0,0,0,
                 cctk_lsh[0], cctk_lsh[1], cctk_lsh[2],
                 cctk_lsh[0], cctk_lsh[1], cctk_lsh[2])
        {
          CCTK_INT i3D = CCTK_GFINDEX3D(cctkGH, i, j, k);
          CCTK_REAL dist2 = (sf_centroid_x[sfi]-x[i3D]) * (sf_centroid_x[sfi]-x[i3D]) +
                            (sf_centroid_y[sfi]-y[i3D]) * (sf_centroid_y[sfi]-y[i3D]) +
                            (sf_centroid_z[sfi]-z[i3D]) * (sf_centroid_z[sfi]-z[i3D]);
          if (dist2 < SetMask_RadiusFactor[smi] * sf_min_radius[sfi] *
                      SetMask_RadiusFactor[smi] * sf_min_radius[sfi])
          {
            mask[i3D] = HYDRO_EXCISION_EXCISED;
          }
        }
        LC_ENDLOOP3(setsurface);
      }
    }
  }
}