/*@@
   @file      CoordinateStuff.c
   @date      
   @author    Cactus Maintainers
   @desc 
   Extra functions for coordinate information
   @enddesc 
   @version $Header$
 @@*/

#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)

/********************************************************************
 *********************     Local Data Types   ***********************
 ********************************************************************/


/********************************************************************
 ********************* Local Routine Prototypes *********************
 ********************************************************************/


/********************************************************************
 *********************     Local Data   *****************************
 ********************************************************************/


/********************************************************************
 *****************   External Routines Prototype ********************
 ********************************************************************/

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


/********************************************************************
 **********************   External Routines  ************************
 ********************************************************************/


 /*@@
   @routine    IndexFloorC
   @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
     -2 : illegal dimension
     -3 : coordinate value outside grid
   @enddesc 
   @calls     
   @calledby   
   @history 
 
   @endhistory 

@@*/

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

  if (d>=0 || d<3)
  {
    CCTK_CoordRange(GH,&cmin,&cmax, d+1, NULL, "cart3d");

    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); 
      index_low = -3;
    }
    else
    {
      index_low    = floor((coord_value-cmin)/ GH->cctk_delta_space[d])
	  - GH->cctk_lbnd[d];

      if (index_low<0 || index_low>=GH->cctk_lsh[d])
      {
	index_low = -1;
      }
    }
  }
  else 
  {
    CCTK_WARN(1,"IndexFloorC: dimension is not valid");
    index_low = -2;
  }
 
  return(index_low);
}


 /*@@
   @routine    IndexCeilC
   @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
     -2 : illegal dimension
     -3 : coordinate value outside grid
   @enddesc 
   @calls     
   @calledby   
   @history 
 
   @endhistory 

@@*/


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

  if (d>=0 || d<3)
  {
    CCTK_CoordRange(GH,&cmin,&cmax, d+1, NULL, "cart3d");

    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); 
      index_up = -3;
    }
    else
    {
      index_up   = ceil((coord_value-cmin)/ GH->cctk_delta_space[d])
	  -GH->cctk_lbnd[d];
  
      if (index_up<0 || index_up>=GH->cctk_lsh[d])
      {
	index_up = -1;
      }
    }
  }
  else 
  {
    CCTK_WARN(1,"IndexCeilC: dimension is not valid");
    index_up = -2;
  }

  return(index_up);
}