#include "noise.h"

#include "cctk_Parameters.h"
#include "cctk_Arguments.h"
#include "cctk_FortranString.h"

#include "Symmetry.h"



static int ApplyBndNoise (const cGH *GH,
			  int stencil_dir,
			  const int *stencil_alldirs,
			  int dir,
			  int first_var,
			  int num_vars);

int
BndNoiseGI (const cGH *GH, const int *stencil, int gi)
{
  int first_vi, retval;

  first_vi = CCTK_FirstVarIndexI (gi);
  if (first_vi >= 0)
    {
      retval = ApplyBndNoise (GH, -1, stencil, 0, first_vi,
			      CCTK_NumVarsInGroupI (gi));
    }
  else
    {
      CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
		  "Invalid group index %d in BndNoiseGI", gi);
      retval = -1;
    }

  return (retval);
}

void
CCTK_FCALL CCTK_FNAME (BndNoiseGI) (int *ierr, const cGH *GH,
				    const int *stencil, const int *gi)
{
  *ierr = BndNoiseGI (GH, stencil, *gi);
}


int
BndNoiseGN (const cGH *GH, const int *stencil, const char *gn)
{
  int gi, retval;

  gi = CCTK_GroupIndex (gn);
  if (gi >= 0)
    {
      retval = BndNoiseGI (GH, stencil, gi);
    }
  else
    {
      CCTK_VWarn (2, __LINE__, __FILE__, CCTK_THORNSTRING,
		  "Invalid group name '%s' in BndNoiseGN", gn);
      retval = -1;
    }

  return (retval);
}

void CCTK_FCALL CCTK_FNAME (BndNoiseGN)
     (int *ierr,
      const cGH *GH,
      const int *stencil,
      ONE_FORTSTRING_ARG)
{
  ONE_FORTSTRING_CREATE (gn)
    *ierr = BndNoiseGN (GH, stencil, gn);
  free (gn);
}





static int ApplyBndNoise (const cGH *GH,
			  int stencil_dir,
			  const int *stencil_alldirs,
			  int dir,
			  int first_var,
			  int num_vars)
{
  int i, j, k;
  int var, vtypesize, gindex, gdim, timelvl;
  int doBC[2*MAXDIM], dstag[MAXDIM], lsh[MAXDIM], lssh[MAXDIM],
    stencil[MAXDIM];
  SymmetryGHex *sGHex;
  CCTK_REAL amplitude;
  CCTK_REAL offset;
  int type;

  stencil[0] = 1; stencil[1] = 1; stencil[2] = 1;

  /* get the group index of the variables */
  gindex = CCTK_GroupIndexFromVarI (first_var);

  amplitude = *((CCTK_REAL*) CCTK_ParameterGet("amplitude", "Noise", &type))
    / RAND_MAX;
  offset = 0.5 * amplitude;

  /* get the number of dimensions and the size of the variables' type */
  gdim      = CCTK_GroupDimI (gindex);
  vtypesize = CCTK_VarTypeSize (CCTK_VarTypeI (first_var));

  /* make sure we can deal with this number of dimensions */
  if (gdim > MAXDIM)
    {
      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
		  "ApplyBndNoise: Variable dimension of %d not supported", gdim);
      return (-1);
    }

  /* check the direction parameter */
  if (abs (dir) > gdim)
    {
      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
		  "ApplyBndNoise: direction %d greater than dimension %d",
		  dir, gdim);
      return (-2);
    }

  /* initialize arrays for variables with less dimensions than MAXDIM
     so that we can use the INDEX_3D macro later on */
  for (i = gdim; i < MAXDIM; i++)
    {
      lsh[i]  = 0;
      lssh[i] = 1;
    }

  /* get the directional staggering of the group */
  CCTK_GroupStaggerDirArrayGI (dstag, gdim, gindex);

  /* get the current timelevel */
  timelvl = 0;

  /* see if we have a symmetry array */
  sGHex = (SymmetryGHex *) CCTK_GHExtension (GH, "Symmetry");


  /* now loop over all variables */
  for (var = first_var; var < first_var + num_vars; var++)
    {
      /* Apply condition if:
	 + boundary is not a symmetry boundary (no symmetry or unset(=unsed))
	 + boundary is a physical boundary
	 + have enough grid points
      */
      memset (doBC, 1, sizeof (doBC));
      if (sGHex)
	{
	  for (i = 0; i < 2 * gdim; i++)
	    {
	      doBC[i] = sGHex->GFSym[var][i] == GFSYM_NOSYM ||
		sGHex->GFSym[var][i] == GFSYM_UNSET;
	    }
	}
      for (i = 0; i < gdim; i++)
	{
	  lsh[i]       = GH->cctk_lsh[i];
	  lssh[i]      = GH->cctk_lssh[CCTK_LSSH_IDX (dstag[i], i)];
	  doBC[i*2]   &= GH->cctk_lsh[i] > 1 && GH->cctk_bbox[i*2];
	  doBC[i*2+1] &= GH->cctk_lsh[i] > 1 && GH->cctk_bbox[i*2+1];
	  if (dir != 0)
	    {
	      doBC[i*2]   &= (dir < 0 && (i + 1 == abs (dir)));
	      doBC[i*2+1] &= (dir > 0 && (i + 1 == abs (dir)));
	    }
	}

      /* now apply the boundaries face by face */
      if (gdim > 0)
	{
#ifdef DEBUG_BOUNDARY
	  if (doBC[0])
	    {
	      printf("Boundary: Applying lower x noise to boundary\n");
	    }
	  if (doBC[1])
	    {
	      printf("Boundary: Applying upper x noise to boundary\n");
	    }
#endif /* DEBUG_BOUNDARY */
	  /* lower x */
	  BOUNDARY_NOISE (doBC[0], stencil[0], lssh[1], lssh[2],
			  i, j, k);
	  /* upper x */
	  BOUNDARY_NOISE (doBC[1], stencil[0], lssh[1], lssh[2],
			  lssh[0]-i-1, j, k);

	}
      if (gdim > 1)

	{
#ifdef DEBUG_BOUNDARY
	  if (doBC[2])
	    {
	      printf("Boundary: Applying lower y noise to boundary\n");
	    }
	  if (doBC[3])
	    {
	      printf("Boundary: Applying upper y noise to boundary\n");
	    }
#endif /* DEBUG_BOUNDARY */
	  /* lower y */
	  BOUNDARY_NOISE (doBC[2], lssh[0], stencil[1], lssh[2],
			  i, j, k);
	  /* upper y */
	  BOUNDARY_NOISE (doBC[3], lssh[0], stencil[1], lssh[2],
			  i, lssh[1]-j-1, k);
	}
      if (gdim > 2)
	{
#ifdef DEBUG_BOUNDARY
	  if (doBC[4])
	    {
	      printf("Boundary: Applying lower z noise to boundary\n");
	    }
	  if (doBC[5])
	    {
	      printf("Boundary: Applying upper z noise to boundary\n");
	    }
#endif /* DEBUG_BOUNDARY */
	  /* lower z */
	  BOUNDARY_NOISE (doBC[4], lssh[0], lssh[1], stencil[2],
			  i, j, k);
	  /* upper z */
	  BOUNDARY_NOISE (doBC[5], lssh[0], lssh[1], stencil[2],
			  i, j, lssh[2]-k-1);
	}
    }

  return(0);
}


void
add_bc_noise_to_group (int idx, const char* optstring, void* cctkGH)
{
  int i, j, k, ijk;
  CCTK_REAL* data;
  cGH* GH = cctkGH;
  int sw[3];

  sw[0] = 1;
  sw[1] = 1;
  sw[2] = 1;

  BndNoiseGI(GH, sw, idx);
}


void
bc_noise(CCTK_ARGUMENTS)
{
  DECLARE_CCTK_ARGUMENTS
    DECLARE_CCTK_PARAMETERS

    if (CCTK_TraverseString(noisy_bc_vars, add_bc_noise_to_group, cctkGH,
			    CCTK_GROUP_OR_VAR) < 0)
      {
	CCTK_WARN (1, "Failed to parse 'Noise::noisy_bc_vars' parameter");
      }
}