/* $Header$ */

#include <assert.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>

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

#include "util_ErrorCodes.h"
#include "util_Table.h"

#define MAXDIM 3

void
setup_epsdis (CCTK_ARGUMENTS)
{
  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;

  int ni,nj,nk;
  int i,j,k,s;
  int index,indexP;
  int ierr;
  int npts;
  int *inds;
  CCTK_REAL *xa,*ya,*za,*rads;
  CCTK_REAL maxrad;
  CCTK_REAL xmin,xmax, ymin,ymax, zmin,zmax;
  CCTK_REAL odx,ody,odz;
  CCTK_REAL radp;
  const CCTK_INT  MAXSURFNUM=100; /* XXX hard limit */
  CCTK_INT doBC[2*MAXDIM],symbnd[2*MAXDIM];  
  CCTK_INT symtable;

  if (verbose) {
    CCTK_VInfo(CCTK_THORNSTRING,"Setting up spatially varying dissipation at T=%g",
               (double)cctk_time);
  }

  ni=cctk_lsh[0];
  nj=cctk_lsh[1];
  nk=cctk_lsh[2];

  for (i=0;i<ni*nj*nk;i++) {
    epsdisA[i]=epsdis;
  }



  if (extra_dissipation_at_outerbound) 
  {
    symtable = SymmetryTableHandleForGrid (cctkGH);
    if (symtable < 0) CCTK_WARN (1, "unable to get symmetry table");
    ierr=Util_TableGetIntArray(symtable, 6, symbnd, "symmetry_handle");
    if (ierr != 6) CCTK_WARN (1, "unable to get symmetry handle");
    for (i = 0; i < 6; i++) {
      doBC[i] = cctk_bbox[i]!=0 && symbnd[i] < 0;
    }
    if(doBC[0]) {
      for (k=0;k<nk;k++) {
        for (j=0;j<nj;j++) {
          for (i=0;i<outer_bound_npoints;i++) {
            index = CCTK_GFINDEX3D(cctkGH,i,j,k);
            indexP= CCTK_GFINDEX3D(cctkGH,outer_bound_npoints,j,k); 
            epsdisA[index]=epsdis+ob_slope*abs(x[index]-x[indexP]);
          }
        }
      }
    }
    if(doBC[1]) {
      for (k=0;k<nk;k++) {
        for (j=0;j<nj;j++) {
          for (i=ni;i>ni-outer_bound_npoints;i--) {
            index = CCTK_GFINDEX3D(cctkGH,i,j,k);
            indexP= CCTK_GFINDEX3D(cctkGH,ni-outer_bound_npoints,j,k); 
            epsdisA[index]=epsdis+ob_slope*abs(x[index]-x[indexP]);
          }
        }
      }
    }

    if(doBC[2]) {
      for (k=0;k<nk;k++) {
        for (j=0;j<outer_bound_npoints;j++) {
          for (i=0;i<ni;i++) {
            index = CCTK_GFINDEX3D(cctkGH,i,j,k);
            indexP= CCTK_GFINDEX3D(cctkGH,i,outer_bound_npoints,k); 
            epsdisA[index]=epsdis+ob_slope*abs(y[index]-y[indexP]);
          }
        }
      }
    }
    if(doBC[3]) {
      for (k=0;k<nk;k++) {
        for (j=nj;j>nj-outer_bound_npoints;j--) {
          for (i=0;i<ni;i++) {
            index = CCTK_GFINDEX3D(cctkGH,i,j,k);
            indexP= CCTK_GFINDEX3D(cctkGH,i,nj-outer_bound_npoints,k); 
            epsdisA[index]=epsdis+ob_slope*abs(y[index]-y[indexP]);
          }
        }
      }
    }

    if(doBC[4]) {
      for (k=0;k<outer_bound_npoints;k++) {
        for (j=0;j<nj;j++) {
          for (i=0;i<ni;i++) {
            index = CCTK_GFINDEX3D(cctkGH,i,j,k);
            indexP= CCTK_GFINDEX3D(cctkGH,i,j,outer_bound_npoints); 
            epsdisA[index]=epsdis+ob_slope*abs(z[index]-z[indexP]);
          }
        }
      }
    }
    if(doBC[5]) {
      for (k=nk;k>nk-outer_bound_npoints;k--) {
        for (j=0;j<nj;j++) {
          for (i=0;i<ni;i++) {
            index = CCTK_GFINDEX3D(cctkGH,i,j,k);
            indexP= CCTK_GFINDEX3D(cctkGH,i,j,nk-outer_bound_npoints); 
            epsdisA[index]=epsdis+ob_slope*abs(z[index]-z[indexP]);
          }
        }
      }
    }
  }



  if (extra_dissipation_in_horizons && cctk_iteration%update_ah_every == 0)
  {
    if (verbose) {
      CCTK_INFO("Linear Interpolation into AH surfaces");
    }

    for (s=0;s<MAXSURFNUM;s++)
    {
      if (surface_number[s]==-1 && horizon_number[s]==-1) {
        continue;
      }
      
      if (surface_number[s]<0 || horizon_number[s]<=0) {
        CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
                    "Invalid specification for horizon %d", s);
        continue;
      }
      
      if (! sf_valid[surface_number[s]]) {
        if (verbose) {
          CCTK_VInfo(CCTK_THORNSTRING,
                     "Invalid Surface: s=%d, sf_va=%d, surf_no=%d, ah_no=%d",
                     s,sf_valid[surface_number[s]],
                     surface_number[s],horizon_number[s]);
        }
        continue;
      }
      
      if (! CCTK_IsFunctionAliased ("HorizonRadiusInDirection"))
      {
        CCTK_WARN (1, "The aliased function \"HorizonRadiusInDirection\" must be defined when the parameter \"extra_dissipation_in_horizons\" is set and one of the sources is AHFinderDirect");
        continue;
      }

      maxrad=sf_max_radius[surface_number[s]];
      odx=sf_origin_x[surface_number[s]];
      ody=sf_origin_y[surface_number[s]];
      odz=sf_origin_z[surface_number[s]];
      xmin=odx-maxrad;
      xmax=odx+maxrad;
      ymin=ody-maxrad;
      ymax=ody+maxrad;
      zmin=odz-maxrad;
      zmax=odz+maxrad;

      npts=0;
      for (i=0;i<ni*nj*nk;i++) {
        if (  x[i]<=xmax&&x[i]>=xmin
            &&y[i]<=ymax&&y[i]>=ymin
            &&z[i]<=zmax&&z[i]>=zmin)
        {
          npts++;
        }
      }

      xa=(CCTK_REAL *)   malloc(npts*sizeof(CCTK_REAL));
      ya=(CCTK_REAL *)   malloc(npts*sizeof(CCTK_REAL));
      za=(CCTK_REAL *)   malloc(npts*sizeof(CCTK_REAL));
      rads=(CCTK_REAL *) malloc(npts*sizeof(CCTK_REAL));
      inds=(CCTK_INT *)  malloc(npts*sizeof(CCTK_INT));

      for (i=0;i<npts;i++) {
        rads[i]=0;
        inds[i]=0;
      }

      j=0;
      for (i=0;i<ni*nj*nk;i++) {
        if (  x[i]<=xmax&&x[i]>=xmin
            &&y[i]<=ymax&&y[i]>=ymin
            &&z[i]<=zmax&&z[i]>=zmin)
        {
          xa[j]=x[i];
          ya[j]=y[i];
          za[j]=z[i];
          inds[j]=i;
          j++;
        }
      }

      ierr=HorizonRadiusInDirection(horizon_number[s],
                                    npts,
                                    xa, ya, za,  rads);
      assert(!ierr);

      for (i=0;i<npts;i++) {
        radp=sqrt((xa[i]-odx)*(xa[i]-odx)+(ya[i]-ody)*(ya[i]-ody)+
                  (za[i]-odz)*(za[i]-odz));
        if (radp<=rads[i]) {
          epsdisA[inds[i]]=epsdis+ ah_slope*(rads[i]-radp);
          if (epsdisA[inds[i]] > ah_max_epsdis) {
            epsdisA[inds[i]] = ah_max_epsdis;
          }
        }
      }

      free(xa);
      free(ya);
      free(za);
      free(rads);
      free(inds);
    }
  }
}