fix a nasty bug where, if the interpolation coordinate was close to

a grid boundary, the actual interpolation was sometimes done as if
the interpolation coordinate was displaced by one grid spacing :( :(

also add more debugging code (introduced to help track down this bug,
but there's no harm in leaving it in for possible future use)

also add explicit checks for a couple of error cases:
* interpolation coordinate is not a finite number (eg NaN or infinity)
* grid origin/delta not a finite number (eg NaN or infinity)
* grid spacing is 0 (we need to divide by this!)

also add explicit note about the test driver in this directory


git-svn-id: http://svn.aei.mpg.de/numrel/AEIThorns/AEILocalInterp/trunk@35 0f49ee68-0e4f-0410-9b9c-b2c123ded7ef

3 files changed, 53 insertions, 19 deletions

diff --git a/src/InterpLocalUniform.h b/src/InterpLocalUniform.h
index a111f02..bb412d0 100644
--- a/src/InterpLocalUniform.h
+++ b/src/InterpLocalUniform.h
@@ -35,7 +35,7 @@ typedef int bool;
 /* round floating-point value to nearest integer */
 /* ... result is expressed as floating point! */
 /* ... needs <math.h> for floor() */
-#define JT_ROUND(x)	floor((x) + 0.5)
+#define ROUND_TO_INTEGER__RESULT_IS_FP(x)	floor((x) + 0.5)
 
 /******************************************************************************/
 
@@ -175,6 +175,9 @@ struct	Jacobian_info
 /* grid is smaller than molecule */
 #define MOLECULE_POSN_ERROR_GRID_TINY	(-3)
 
+/* someone passed us a NaN or other non-finite floating-point number */
+#define MOLECULE_POSN_ERROR_NAN		(-4)
+
 /******************************************************************************/
 
 /*
diff --git a/src/molecule_posn.c b/src/molecule_posn.c
index 663d917..89debab 100644
--- a/src/molecule_posn.c
+++ b/src/molecule_posn.c
@@ -78,10 +78,10 @@ static const char *rcsid = "$Header$";
  @hdate    27.Jan.2003
  @hauthor  Jonathan Thornburg <jthorn@aei.mpg.de>
  @hdesc    Complete rewrite: now supports
-           @var boundary_off_centering_tolerance_min,
-           @var boundary_off_centering_tolerance_max,
-           @var boundary_extrapolation_tolerance_min,
-           @var boundary_extrapolation_tolerance_max,
+             @var boundary_off_centering_tolerance_min,
+             @var boundary_off_centering_tolerance_max,
+             @var boundary_extrapolation_tolerance_min,
+             @var boundary_extrapolation_tolerance_max,
            also change to returning status code,
            also drop returning @var min_m and @var max_m
            since they were never used.
@@ -185,6 +185,17 @@ int AEILocalInterp_molecule_posn(fp grid_origin, fp grid_delta,
 				 int debug,
 				 int* i_center, fp* x_rel)
 {
+/*
+ * ***** IMPORTANT *****
+ *
+ * This code is ++tricky.  There's a fairly rigorous test driver for
+ * it in the file  test_molecule_posn.c  in this directory.  This is a
+ * standalone test driver, which can be compiled via
+ *	gmake -f Makefile.standalone
+ * in this directory.  If you make any changes to this code, please
+ * verify that all the tests still pass!
+ */
+
 if (debug >= 8)
    then {
 	printf("AEILocalInterp_molecule_posn():\n");
@@ -201,17 +212,38 @@ if (debug >= 8)
 	printf("   x=%g\n", (double) x);
 	}
 
+/*
+ * basic sanity checks
+ */
+
+/* is the molecule larger than the grid? */
+if (molecule_size > HOW_MANY_IN_RANGE(grid_i_min,grid_i_max))
+   then return MOLECULE_POSN_ERROR_GRID_TINY;		/*** ERROR RETURN ***/
+
+/* has someone passed us NaNs for coordinates etc? */
+if (!finite(grid_origin) || !finite(grid_delta) || !finite(x))
+   then return MOLECULE_POSN_ERROR_NAN;
+
+/* ditto if the grid spacing is zero (we'll need to divide by it)! */
+if (grid_delta == 0.0)
+   then return MOLECULE_POSN_ERROR_NAN;
+
+
 /* molecule radia (inherently positive) in +/- directions */
 const int mr_plus  = (molecule_size >> 1);
 const int mr_minus = molecule_size - mr_plus - 1;
 
-/* range of x_rel for which this molecule is centered, cf. diagram in header comment */
+/* range of x_rel for which this molecule is centered, */
+/* cf. diagram in header comment */
 const fp centered_min_x_rel = IS_EVEN(molecule_size) ? 0.0 : -0.5;
 const fp centered_max_x_rel = IS_EVEN(molecule_size) ? 1.0 : +0.5;
 
-/* range of i where we *could* center the molecule, as floating-point numbers */
-const fp fp_centered_min_possible_i = grid_i_min + mr_minus + centered_min_x_rel;
-const fp fp_centered_max_possible_i = grid_i_max - mr_plus  + centered_max_x_rel;
+/* range of i where a centered molecule would fit within the grid, */
+/* as floating-point numbers */
+const fp fp_centered_min_possible_i = grid_i_min + mr_minus
+						 + centered_min_x_rel;
+const fp fp_centered_max_possible_i = grid_i_max - mr_plus
+						 + centered_max_x_rel;
 
 /* integer coordinate i of interpolation point, as a floating-point number */
 const fp fp_i = (x - grid_origin) / grid_delta;
@@ -227,10 +259,6 @@ if (debug > 8)
 	printf("   fp_i=%g\n", (double) fp_i);
 	}
 
-/* is the molecule larger than the grid? */
-if (molecule_size > HOW_MANY_IN_RANGE(grid_i_min,grid_i_max))
-   then return MOLECULE_POSN_ERROR_GRID_TINY;		/*** ERROR RETURN ***/
-
 /* is interpolation point x beyond the extrapolation tolerance? */
 if (fp_i < (fp)grid_i_min - boundary_extrapolation_tolerance_min)
    then return MOLECULE_POSN_ERROR_X_LT_MIN;		/*** ERROR RETURN ***/
@@ -246,22 +274,23 @@ if (fp_i > fp_centered_max_possible_i + boundary_off_centering_tolerance_max)
 /* now choose the actual molecule position/centering: */
 /* first set up a centered molecule */
   {
-fp fp_i_center = IS_EVEN(molecule_size)	/* ... as a floating-point number */
+/* compute as a floating point number */
+/* (because that's what the C math library provides), */
+/* but this value is guaranteed to be integral */
+fp fp_i_center = IS_EVEN(molecule_size)
 		 ? floor(fp_i)
-		 : JT_ROUND(fp_i);
-int int_i_center = (int) fp_i_center;	/* ... as an integer */
+		 : ROUND_TO_INTEGER__RESULT_IS_FP(fp_i);
+int int_i_center = (int) fp_i_center;	/* convert to integer */
 if (debug > 8)
    then printf("   initial: fp_i_center=%g int_i_center=%d\n",
 	       (double) fp_i_center, int_i_center);
 
-/* then clamp the molecule at the grid boundaries */
-if (int_i_center < grid_i_min) int_i_center = grid_i_min;
+/* clamp the molecule at the grid boundaries */
 if (int_i_center < grid_i_min + mr_minus)
    then {
 	int_i_center = grid_i_min + mr_minus;
 	fp_i_center = (fp) int_i_center;
 	}
-if (int_i_center > grid_i_max) int_i_center = grid_i_max;
 if (int_i_center > grid_i_max - mr_plus)
    then {
 	int_i_center = grid_i_max - mr_plus;
diff --git a/src/template.c b/src/template.c
index c1ad45f..09ed7b3 100644
--- a/src/template.c
+++ b/src/template.c
@@ -974,6 +974,8 @@ if (debug > 0)
 	case MOLECULE_POSN_ERROR_X_GT_MAX:
 		this_axis_status = CCTK_ERROR_INTERP_POINT_OUTSIDE;
 		break;
+	case MOLECULE_POSN_ERROR_NAN:
+		this_axis_status = CCTK_ERROR_INTERP_COORD_NAN;
 	default:
 		CCTK_VWarn(BUG_MSG_SEVERITY_LEVEL,
 			   __LINE__, __FILE__, CCTK_THORNSTRING,