#include <stddef.h>
#include <math.h>

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

#define SQ(x) ((x)*(x))

/*
 * Initial data: gaussian pulse
 */
void WaveToy_InitialData(CCTK_ARGUMENTS)
{
  DECLARE_CCTK_ARGUMENTS
  DECLARE_CCTK_PARAMETERS

  CCTK_INT i,j,k,vindex;
  CCTK_REAL dt, X, Y, Z, R;

  dt = CCTK_DELTA_TIME;

  for(k=0; k<cctk_lsh[2]; k++) {
    for(j=0; j<cctk_lsh[1]; j++) {
      for(i=0; i<cctk_lsh[0]; i++) {
        vindex =  CCTK_GFINDEX3D(cctkGH,i,j,k);

        X = x[vindex];
        Y = y[vindex];
        Z = z[vindex];

        R = sqrt(X*X + Y*Y + Z*Z);

        phi[vindex] = amplitude*exp( - SQ( (R - radius) / sigma ) );

        if (R == 0.0) {
          phi_p[vindex] = amplitude*(1.0 - 2.0*dt*dt/sigma)*exp(-dt*dt/sigma);
        } else {
          phi_p[vindex] = amplitude/2.0*(R-dt)/R*
            exp( - SQ( (R - radius - dt)/ sigma ) )
            + amplitude/2.0*(R+dt)/R*
            exp( - SQ( (R - radius + dt)/ sigma ) );
        }
      }
    }
  }
}

/* 
 * Evolution routine
 */
void  WaveToy_Evolution(CCTK_ARGUMENTS)
{
  DECLARE_CCTK_ARGUMENTS;

  CCTK_INT i,j,k;
  CCTK_INT vindex;
  CCTK_REAL dx2i,dy2i,dz2i,dt2,factor;

  /* Set up shorthands */
  dx2i = 1.0/SQ(CCTK_DELTA_SPACE(0));
  dy2i = 1.0/SQ(CCTK_DELTA_SPACE(1));
  dz2i = 1.0/SQ(CCTK_DELTA_SPACE(2));
  dt2  = SQ(CCTK_DELTA_TIME);

  factor = 2*(1 - (dt2)*(dx2i + dy2i + dz2i));

  /* Main loop */
  for (k=1; k<cctk_lsh[2]-1; k++) {
    for (j=1; j<cctk_lsh[1]-1; j++) {
      for (i=1; i<cctk_lsh[0]-1; i++) {
        vindex = CCTK_GFINDEX3D(cctkGH,i,j,k);

        phi[vindex] = factor*
                     phi_p[vindex] - phi_p_p[vindex]
                     + (dt2) *
                   ( ( phi_p[CCTK_GFINDEX3D(cctkGH,i+1,j  ,k  )]
                      +phi_p[CCTK_GFINDEX3D(cctkGH,i-1,j  ,k  )] )*dx2i
                    +( phi_p[CCTK_GFINDEX3D(cctkGH,i  ,j+1,k  )]
                      +phi_p[CCTK_GFINDEX3D(cctkGH,i  ,j-1,k  )] )*dy2i
                    +( phi_p[CCTK_GFINDEX3D(cctkGH,i  ,j  ,k+1)]
                      +phi_p[CCTK_GFINDEX3D(cctkGH,i  ,j,  k-1)] )*dz2i);
      }
    }
  }
}


/*
 * Select boundary condition
 */
void WaveToy_Boundaries(CCTK_ARGUMENTS)
{
  DECLARE_CCTK_ARGUMENTS;

  const char *bctype = "scalar";

  /* Uses all default arguments, so invalid table handle -1 can be passed */
  if (Boundary_SelectVarForBC (cctkGH, CCTK_ALL_FACES, 1, -1,
                                         "wavetoy_simple::phi", bctype) < 0)  {
    CCTK_WARN (0, "WaveToy_Boundaries: Error selecting boundary condition");
  }
}