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

#include "cctk.h"

 /*@@
   @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, 
     which is owned by the grid path. (Ghostzone are not owned!)
     Other return values: 
     -1 : gridindex not on this gridpatch
     -3 : coordinate is outside global grid

     This routine and the Fortran wrapper are Fortran/C index aware.
   @enddesc 
   @calls     
   @calledby   
   @history 
 
   @endhistory 

@@*/


int IndexFloor(cGH *GH, CCTK_REAL coord_value, int d)
{ 
  int ierr,index_low,idummy,ugs,lgs;
  char *message;
  CCTK_REAL cmin,cmax;

  if (d==0)
    CCTK_CoordRange(GH,&cmin,&cmax,"x");
  else if (d==1)
    CCTK_CoordRange(GH,&cmin,&cmax,"y");
  else if (d==2)
    CCTK_CoordRange(GH,&cmin,&cmax,"z");
 else 
   CCTK_WARN(0,"IndexFloor: dimension is not valid. Fortran:1..3  C:0..2");
  
  if ((coord_value<cmin)||(coord_value>cmax)) {
    message = (char*)malloc(128*sizeof(char));
    sprintf(message,
	    "IndexFloor: coordinate value outside grid [%f,%f]:  %f \n", 
	    cmin, cmax, coord_value); 
    CCTK_WARN(2,message);
    free(message);
    return(-3);
  }
  
  idummy    = floor((coord_value-GH->cctk_origin_space[d])/ GH->cctk_delta_space[d]);
  index_low = idummy-GH->cctk_lbnd[d];

  if (GH->cctk_bbox[2*d]==1) 
    lgs = 0;
  else
    lgs = GH->cctk_nghostzones[d];

  if (GH->cctk_bbox[2*d+1]==1)
    ugs = 0;
  else
    ugs = GH->cctk_nghostzones[d];

  if (index_low<lgs)
    index_low = -1;

  if (index_low>=GH->cctk_lsh[d]-ugs)
    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, 
     which is owned by the grid path. (Ghostzone are not owned!)
     Other return values: 
     -1 : grid index not on this gridpatch
     -3 : coordinate is outside global grid

     This routine and the Fortran wrapper are Fortran/C index aware.
   @enddesc 
   @calls     
   @calledby   
   @history 
 
   @endhistory 

@@*/


int IndexCeil(cGH *GH, CCTK_REAL coord_value, int d)
{  
  int ierr,index_up,idummy,ugs,lgs;
  char *message;
  CCTK_REAL cmin,cmax;

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

  if ((coord_value<cmin)||(coord_value>cmax)) { 
    message = (char*) malloc(256*sizeof(char));
    sprintf(message,
	    "IndexCeil: coordinate value outside grid [%f,%f]:  %f \n", 
	    cmin,cmax, coord_value);
    CCTK_WARN(2,message);
    free(message);
    return(-3);
  }
  
  idummy   = ceil((coord_value-GH->cctk_origin_space[d])/ GH->cctk_delta_space[d]);
  index_up = idummy-GH->cctk_lbnd[d];

  if (GH->cctk_bbox[2*d]==1) 
    lgs = 0;
  else
    lgs = GH->cctk_nghostzones[d];

  if (GH->cctk_bbox[2*d+1]==1)
    ugs = 0;
  else
    ugs = GH->cctk_nghostzones[d];

  if (index_up<lgs)
    index_up = -1;

  if (index_up>=GH->cctk_lsh[d]-ugs)
    index_up = -1;
   
  return(index_up);
}

void FMODIFIER FORTRAN_NAME(IndexFloor)
  (cGH *GH, CCTK_REAL *coord_value, int *index_low, int *fdim) {
  
  /* note the increase/decrease of indeces for fortran compatability */
  *index_low = IndexFloor(GH, *coord_value,  (*fdim)-1);
  if ((*index_low)>=0) (*index_low)++; 
}

void FMODIFIER FORTRAN_NAME(IndexCeil)
  (cGH *GH, CCTK_REAL *coord_value, int *index_up, int *fdim) {
  /* note the increase/decrease of indeces for fortran compatability */
  *index_up = IndexCeil(GH, *coord_value, (*fdim)-1);
  if ((*index_up)>=0) (*index_up)++; 
}