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

#include "cctk.h"

static const char *rcsid = "$Header$";

CCTK_FILEVERSION(CactusWave_WaveBinarySource_CoordinateStuff_c)

int IndexFloorC(cGH *GH, CCTK_REAL coord_value, int d);
int IndexCeilC(cGH *GH, CCTK_REAL coord_value, int d);


 /*@@
   @routine    IndexFloor
   @date       Fri Jan  7 10:34:29 2000
   @author     Gerd Lanfermann
   @desc 
     For a given physical, global coordinate, IndexFloor returns
     the closest lower grid index *locally* in the direction d,
     this includes ghostzones;

     Other return values: 
     -1 : gridindex not on this gridpatch
   @enddesc 
   @calls     
   @calledby   
   @history 
 
   @endhistory 

@@*/

int IndexFloorC(cGH *GH, CCTK_REAL coord_value, int d)
{ 
  int index_low,idummy;
  CCTK_REAL cmin,cmax;

  if (d==0)
  {
    CCTK_CoordRange(GH,&cmin,&cmax, 1, "x", "cart3d");
  }
  else if (d==1)
  {
    CCTK_CoordRange(GH,&cmin,&cmax, 2, "y", "cart3d");
  }
  else if (d==2)
  {
    CCTK_CoordRange(GH,&cmin,&cmax, 3, "z", "cart3d");
  }
  else 
  {
    CCTK_WARN(0,"IndexFloor: dimension is not valid. Fortran:1..3  C:0..2");
  }
 
  if ((coord_value<cmin)||(coord_value>cmax)) 
  {
    CCTK_VWarn(2,__LINE__,__FILE__,CCTK_THORNSTRING, 
	       "IndexFloor: coordinate value outside grid [%f,%f]:  %f \n", 
	       cmin, cmax, coord_value); 
    return(-3);
  }
  
  idummy    = floor((coord_value-cmin)/ GH->cctk_delta_space[d]);
  index_low = idummy-GH->cctk_lbnd[d];

  if (index_low<0)
  {
    index_low = -1;
  }
  if (index_low>=GH->cctk_lsh[d])
  {
    index_low = -1;
  } 
  return(index_low);
}


 /*@@
   @routine    IndexCeil
   @date       Fri Jan  7 10:34:29 2000
   @author     Gerd Lanfermann
   @desc 
     For a given physical, global coordinate, IndexCeil returns
     the closest upper grid index *locally* in the direction d
     including ghostzones.
     Other return values: 
     -1 : grid index not on this gridpatch
   @enddesc 
   @calls     
   @calledby   
   @history 
 
   @endhistory 

@@*/


int IndexCeilC(cGH *GH, CCTK_REAL coord_value, int d)
{  
  int index_up,idummy;
  CCTK_REAL cmin,cmax;

  if (d==0)
  {
    CCTK_CoordRange(GH,&cmin,&cmax, -1, "x", "cart3d");
  }
  else if (d==1)
  {
    CCTK_CoordRange(GH,&cmin,&cmax, -1, "y", "cart3d");
  }
  else if (d==2)
  {
    CCTK_CoordRange(GH,&cmin,&cmax, -1, "z", "cart3d");
  }
  else 
  {
    CCTK_WARN(0,"IndexCeilC: dimension is not valid. Fortran:1..3  C:0..2");
  }

  if ((coord_value<cmin)||(coord_value>cmax)) 
  { 
    CCTK_VWarn(2,__LINE__,__FILE__,CCTK_THORNSTRING,
	       "IndexCeil: coordinate value outside grid [%f,%f]:  %f \n", 
	       cmin,cmax, coord_value);
    return(-3);
  }
  
  idummy   = ceil((coord_value-cmin)/ GH->cctk_delta_space[d]);
  index_up = idummy-GH->cctk_lbnd[d];
  
  if (index_up<0)
  {
    index_up = -1;
  }
  if (index_up>=GH->cctk_lsh[d])
  {
    index_up = -1;
  } 
  return(index_up);
}