Code changes

* support the same off-centering and boundary options as the new global
  interpolator (which needs the support here to work):  this means changing
  the former parameter-table entry    off_centering_tolerance   to the
  new entries
	boundary_off_centering_tolerance
	boundary_extrapolation_tolerance
* implement the    N_boundary_points_to_omit   parameter-table entry,
  (previously this was documented but not implemented)
* completely rewrite the   LocalInterp_molecule_posn()   function
  and its test driver   to support the above changes
* change error return code from old CCTK_ERROR_INTERP_POINT_X_RANGE
  to new CCTK_ERROR_INTERP_POINT_OUTSIDE to avoid confusion with
  the old name suggesting that this was specific to the X coordinate
* fix documentation of interpolation operator names for older
  (uniform Cartesian) CCTK_InterpLocal() interpolator -- this
  closes Erik Schnetter's bug report CactusBase/1367

Documentation changes
* document the above code changes
* much expanded description of molecule centering and boundary handling,
  including 2 new diagrams (done in latex picture mode :)
* update README files to list all code files
* move caching parameters (not implemented, and probably won't be
  implemented any time soon) out of documentation.tex into a new file
  future-ideas.tex
* general cleanup of documentation.tex (not quite done, but 80% or so...)
  to more clearly describe what this thorn actually implements, rather
  than a generic spec of which this thorn only implements some pieces


git-svn-id: http://svn.cactuscode.org/arrangements/CactusBase/LocalInterp/trunk@125 df1f8a13-aa1d-4dd4-9681-27ded5b42416

7 files changed, 654 insertions, 539 deletions

diff --git a/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c b/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c
index 56f386a..1f63a03 100644
--- a/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c
+++ b/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c
@@ -38,7 +38,7 @@
 #include "Hermite/all_prototypes.h"
 
 /* the rcs ID and its dummy function to use it */
-static const char *rcsid = "$Header$";
+static const char* rcsid = "$Header$";
 CCTK_FILEVERSION(CactusBase_LocalInterp_src_GeneralizedPolynomialUniform_InterpLocalUniform_c)
 
 /******************************************************************************/
@@ -231,8 +231,20 @@ return LocalInterp_InterpLocalUniform(interp_operator_Hermite,
 	and calls the appropriate subfunction to do the actual
 	interpolation, then finally stores some results back in
 	the parameter table.
+
+	FIXME:
+	This function is huge, and should be split up into 
+	reasonable-sized subfunctions.
   @enddesc
 
+  @hdate    28.Jan.2003
+  @hauthor  Jonathan Thornburg <jthorn@aei.mpg.de>
+  @hdesc    Support parameter-table entries
+            @var N_boundary_points_to_omit,
+            @var boundary_off_centering_tolerance, and
+            @var boundary_extrapolation_tolerance.
+  @endhdesc
+
   ***** misc arguments *****
 
   @var		interp_operator
@@ -349,26 +361,43 @@ return LocalInterp_InterpLocalUniform(interp_operator_Hermite,
   @vtype	CCTK_INT order
   @endvar
 
-  @var		out_of_range_tolerance
-  @vdesc	Specifies how out-of-range interpolation points should
-		be handled.  The array elements are matched up with
-		the axes and minimum/maximum ends of the grid in the
+  @var		N_boundary_points_to_omit
+  @vdesc	If this key is present, then this array specifies how many
+		input grid points to omit (not use for the interpolation)
+		on each grid boundary.  The array elements are matched up
+		with the axes and minimum/maximum ends of the grid in the
 		order [x_min, x_max, y_min, y_max, z_min, z_max, ...].
-		An array value TOL is interpreted as follows:
-		If TOL >= 0.0,
-		   then an interpolation point is considered to be
-			"out of range" if and only if the interpolation
-			point is > TOL * coord_delta[axis]
-			outside the grid in this coordinate direction.
-		If TOL == -1.0,
-		   then an interpolation point is considered to be
-			"out of range" if and only if a centered molecule
-			(or more generally, a molecule whose centering
-			is chosen pretending that the grid is of infinite
-			extent), would require data from outside the grid
-			in this direction.
-		Other values of TOL are illegal.
-  @vtype	const CCTK_REAL out_of_range_tolerance[2*N_dims]
+		If this key isn't present, it defaults to all zeros, i.e.
+		to use all the input grid points in the interpolation.
+  @vtype	const CCTK_INT N_boundary_points_to_omit[2*N_dims]
+  @endvar
+
+  @var		boundary_off_centering_tolerance
+  @vdesc	If this key is present, then the interpolator allows
+		interpolation points to be up to (<=) this many grid
+		spacings outside the default-centering region on each
+		grid boundary, off-centering the interpolation molecules
+		as necessary.
+		The array elements are matched up with the axes and
+		minimum/maximum ends of the grid in the order
+		[x_min, x_max, y_min, y_max, z_min, z_max, ...].
+		If this key isn't present, it defaults to all infinities,
+		i.e. to place no restriction on the interpolation point.
+  @vtype	const CCTK_REAL boundary_off_centering_tolerance[2*N_dims]
+  @endvar
+
+  @var		boundary_extrapolation_tolerance
+  @vdesc	If this key is present, then the interpolator allows
+		interpolation points to be up to (<=) this many grid
+		spacings outside the valid-data region on each grid
+		boundary, doing extrapolation instead of interpolation
+		as necessary.
+		The array elements are matched up with the axes and
+		minimum/maximum ends of the grid in the order
+		[x_min, x_max, y_min, y_max, z_min, z_max, ...].
+		If this key isn't present, it defaults to all 1.0e-10,
+		i.e. to allow up to 1e-10 grid spacings of extrapolation.
+  @vtype	const CCTK_REAL boundary_extrapolation_tolerance[2*N_dims]
   @endvar
 
   @var		input_array_offsets
@@ -625,7 +654,7 @@ static
  * would falsely detect as an out-of-memory condition.
  *
  * As a workaround, we pad all our malloc request sizes, i.e.
- *	CCTK_INT *const p = malloc((N+1) * sizeof(CCTK_INT));
+ *	CCTK_INT* const p = (CCTK_INT*) malloc((N+1) * sizeof(CCTK_INT));
  *	if (p == NULL)
  *	   then return UTIL_ERROR_NO_MEMORY
  * This is a kludge, but so are the other plausible solutions. :( :(
@@ -657,7 +686,7 @@ if (    (N_dims <= 0)
    then {
 	CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
 		   __LINE__, __FILE__, CCTK_THORNSTRING,
-"InterpLocalUniform(): invalid arguments\n"
+"CCTK_InterpLocalUniform(): invalid arguments\n"
 "   (N_dims <= 0, param_table_handle < 0, N_input_arrays < 0,\n"
 "    N_output_arrays < 0, and/or NULL pointers-that-shouldn't-be-NULL)!");
 	return UTIL_ERROR_BAD_INPUT;			/*** ERROR RETURN ***/
@@ -684,14 +713,14 @@ if (N_dims > MAX_N_DIMS)
 CCTK_INT order;
 status = Util_TableGetInt(param_table_handle, &order, "order");
 if	(status == 1)
-   then { }	/* ok ==> no-op here */
+   then { }	/* value set from table ==> no-op here */
 else	{
 	/* n.b. order is a mandatory parameter (no default)!*/
 	CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
 		   __LINE__, __FILE__, CCTK_THORNSTRING,
 "\n"
-"   CCTK_InterpLocalUniform(): bad or missing table entry\n"
-"                              for \"order\" parameter!\n"
+"   CCTK_InterpLocalUniform(): bad or missing table entry for\n"
+"                              \"order\" parameter!\n"
 "                              (this is a mandatory parameter)\n"
 "                              error status=%d",
 		   status);
@@ -708,59 +737,122 @@ if ((order < 1) || (order > MAX_ORDER))
 		   MAX_ORDER);
 	return UTIL_ERROR_BAD_INPUT;			/*** ERROR RETURN ***/
 	}
-
+   
+/******************************************************************************/
 
 /*
  * out-of-range interpolation-point handling
  */
   {
-CCTK_REAL out_of_range_tolerance[2*MAX_N_DIMS];
-const int N_tolerances = 2*N_dims;
-status = Util_TableGetRealArray(param_table_handle,
-				N_tolerances, out_of_range_tolerance,
-				"out_of_range_tolerance");
+CCTK_INT  N_boundary_points_to_omit[MAX_N_BOUNDARIES];
+CCTK_REAL boundary_off_centering_tolerance[MAX_N_BOUNDARIES];
+CCTK_REAL boundary_extrapolation_tolerance[MAX_N_BOUNDARIES];
+const int N_boundaries = 2*N_dims;
+
+status = Util_TableGetIntArray(param_table_handle,
+			       N_boundaries, N_boundary_points_to_omit,
+			       "N_boundary_points_to_omit");
 if	(status == UTIL_ERROR_TABLE_NO_SUCH_KEY)
    then {
-	/* default */
-	int tol_index;
-		for (tol_index = 0 ; tol_index < N_tolerances ; ++tol_index)
+	/* set the default value */
+	int ibndry;
+		for (ibndry = 0 ; ibndry < N_boundaries ; ++ibndry)
 		{
-		out_of_range_tolerance[tol_index]
-			= OUT_OF_RANGE_TOLERANCE_DEFAULT;
+		N_boundary_points_to_omit[ibndry] = 0;
 		}
 	}
-else if (status == N_tolerances)
+else if (status == N_boundaries)
+   then { }	/* value set from table ==> no-op here */
+else	{
+	CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
+		   __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+"   CCTK_InterpLocalUniform(): bad table entry for\n"
+"                              \"N_boundary_points_to_omit\" parameter!\n"
+"                              error status=%d",
+		   status);
+	return status;					/*** ERROR RETURN ***/
+	}
+
+/**************************************/
+
+status = Util_TableGetRealArray(param_table_handle,
+				N_boundaries, boundary_off_centering_tolerance,
+				"boundary_off_centering_tolerance");
+if	(status == UTIL_ERROR_TABLE_NO_SUCH_KEY)
    then {
-	/* check that all values are valid */
-	int tol_index;
-		for (tol_index = 0 ; tol_index < N_tolerances ; ++tol_index)
+	/* set the default value */
+	int ibndry;
+		for (ibndry = 0 ; ibndry < N_boundaries ; ++ibndry)
 		{
-		if (! (    (out_of_range_tolerance[tol_index] >=  0.0)
-			|| (out_of_range_tolerance[tol_index] == -1.0)    )  )
-		   then {
-			CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
-				   __LINE__, __FILE__, CCTK_THORNSTRING,
+		boundary_off_centering_tolerance[ibndry]
+			= BOUNDARY_OFF_CENTER_TOL_DEFAULT;
+		}
+	}
+else if (status == N_boundaries)
+   then { }	/* value set from table ==> no-op here */
+else	{
+	CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
+		   __LINE__, __FILE__, CCTK_THORNSTRING,
 "\n"
-"   CCTK_InterpLocalUniform():\n"
-"        invalid out_of_range_tolerance[tol_index=%d] = %g!\n"
-"        (valid values are -1.0 or >= 0.0)",
-				   tol_index,
-				   out_of_range_tolerance[tol_index]);
-			return UTIL_ERROR_BAD_INPUT;	/*** ERROR RETURN ***/
-			}
+"   CCTK_InterpLocalUniform(): bad table entry for\n"
+"                              \"boundary_off_centering_tolerance\" parameter!\n"
+"                              error status=%d",
+		   status);
+	return status;					/*** ERROR RETURN ***/
+	}
+
+/**************************************/
+
+status = Util_TableGetRealArray(param_table_handle,
+				N_boundaries, boundary_extrapolation_tolerance,
+				"boundary_extrapolation_tolerance");
+if	(status == UTIL_ERROR_TABLE_NO_SUCH_KEY)
+   then {
+	/* set the default value */
+	int ibndry;
+		for (ibndry = 0 ; ibndry < N_boundaries ; ++ibndry)
+		{
+		boundary_extrapolation_tolerance[ibndry]
+			= BOUNDARY_EXTRAP_TOL_DEFAULT;
 		}
 	}
+else if (status == N_boundaries)
+   then { }	/* value set from table ==> no-op here */
 else	{
 	CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
 		   __LINE__, __FILE__, CCTK_THORNSTRING,
 "\n"
-"   CCTK_InterpLocalUniform(): bad table entry\n"
-"                              for \"out_of_range_tolerance\" parameter!\n"
+"   CCTK_InterpLocalUniform(): bad table entry for\n"
+"                              \"boundary_extrapolation_tolerance\" parameter!\n"
 "                              error status=%d",
 		   status);
 	return status;					/*** ERROR RETURN ***/
 	}
 
+/* check for almost-always-a-mistake settings, warn if found */
+  {
+int ibndry;
+	for (ibndry = 0 ; ibndry < N_boundaries ; ++ibndry)
+	{
+	if (    (boundary_off_centering_tolerance[ibndry] == 0.0)
+	     && (boundary_extrapolation_tolerance[ibndry] > 0.0)    )
+	   then CCTK_VWarn(WARNING_MSG_SEVERITY_LEVEL,
+			   __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+"   CCTK_InterpLocalUniform(): warning: setting\n"
+"                              boundary_off_centering_tolerance[%d] == 0.0\n"
+"                              and\n"
+"                              boundary_extrapolation_tolerance[%d] > 0.0\n"
+"                              is almost certainly a mistake\n"
+"                              (the boundary_extrapolation_tolerance[%d]\n"
+"                               setting will have no effect because of the\n"
+"                              boundary_off_centering_tolerance[%d] setting)"
+			   ,
+			   ibndry, ibndry, ibndry, ibndry);
+	}
+  }
+
 /******************************************************************************/
 
 /*
@@ -768,23 +860,23 @@ else	{
  */
   {
 #define MOLECULE_FAMILY_BUFSIZ	(MAX_MOLECULE_FAMILY_STRLEN+1)
-char molecule_family_string[MOLECULE_FAMILY_BUFSIZ];
+char molecule_family_strbuf[MOLECULE_FAMILY_BUFSIZ];
+char* molecule_family_str = molecule_family_strbuf;
 enum molecule_family molecule_family;
 status = Util_TableGetString(param_table_handle,
-			     MOLECULE_FAMILY_BUFSIZ, molecule_family_string,
+			     MOLECULE_FAMILY_BUFSIZ, molecule_family_strbuf,
 			     "molecule_family");
 if	(status == UTIL_ERROR_TABLE_NO_SUCH_KEY)
    then {
-	/* default */
+	/* set the default value for code here*/
 	/* (need enum for main code, string for error messages) */
-	strcpy(molecule_family_string, "cube");
-	molecule_family = molecule_family_cube;
-	}
-else if (status == 0)
-   then {
-	/* this is a query of our default molecule family */
+	molecule_family     = molecule_family_cube;
+	molecule_family_str = "cube";
+
+	/* set this key in the parameter table */
+	/* to give the (default) molecule family we're going to use */
 	status = Util_TableSetString(param_table_handle,
-				     "cube",
+				     molecule_family_str,
 				     "molecule_family");
 	if (status < 0)
 	   then {
@@ -801,13 +893,13 @@ else if (status == 0)
 else if (status > 0)
    then {
 	/* decode molecule family string */
-	if (strcmp(molecule_family_string, "cube") == 0)
+	if (strcmp(molecule_family_strbuf, "cube") == 0)
 	   then molecule_family = molecule_family_cube;
 	else	{
 		CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
 			   __LINE__, __FILE__, CCTK_THORNSTRING,
-"InterpLocalUniform(): unknown molecule_family=\"%s\"!",
-			   molecule_family_string);
+"CCTK_InterpLocalUniform(): unknown molecule_family=\"%s\"!",
+			   molecule_family_strbuf);
 		return UTIL_ERROR_BAD_INPUT;		/*** ERROR RETURN ***/
 		}
 	}
@@ -815,13 +907,15 @@ else	{
 	CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
 		   __LINE__, __FILE__, CCTK_THORNSTRING,
 "\n"
-"   CCTK_InterpLocalUniform(): bad table entry\n"
-"                              for \"molecule_family\" parameter!\n"
+"   CCTK_InterpLocalUniform(): bad table entry for\n"
+"                              \"molecule_family\" parameter!\n"
 "                              error status=%d",
 		   status);
 	return status;					/*** ERROR RETURN ***/
 	}
 
+/**************************************/
+
 /*
  * smoothing
  */
@@ -829,14 +923,15 @@ else	{
 CCTK_INT smoothing;
 status = Util_TableGetInt(param_table_handle, &smoothing, "smoothing");
 if	(status == UTIL_ERROR_TABLE_NO_SUCH_KEY)
-   then smoothing = 0;	/* default */
+   then smoothing = 0;	/* set default value */
 else if (status == 1)
    then { }		/* value set from table ==> no-op here */
 else	{
 	CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
 		   __LINE__, __FILE__, CCTK_THORNSTRING,
 "\n"
-"   CCTK_InterpLocalUniform(): bad table entry for \"smoothing\" parameter!\n"
+"   CCTK_InterpLocalUniform(): bad table entry for\n"
+"                              \"smoothing\" parameter!\n"
 "                              error status=%d",
 		   status);
 	return status;					/*** ERROR RETURN ***/
@@ -850,8 +945,8 @@ if ((smoothing < 0) || (smoothing > MAX_SMOOTHING))
  * input array offsets
  */
   {
-CCTK_INT *const input_array_offsets
-	= malloc(N_input_arrays1 * sizeof(CCTK_INT));
+CCTK_INT* const input_array_offsets
+	= (CCTK_INT*) malloc(N_input_arrays1 * sizeof(CCTK_INT));
 if (input_array_offsets == NULL)
    then return UTIL_ERROR_NO_MEMORY;			/*** ERROR RETURN ***/
 status = Util_TableGetIntArray(param_table_handle,
@@ -869,13 +964,15 @@ else	{
 	CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
 		   __LINE__, __FILE__, CCTK_THORNSTRING,
 "\n"
-"   CCTK_InterpLocalUniform(): bad table entry\n"
-"                              for \"input_array_offsets\" parameter!\n"
+"   CCTK_InterpLocalUniform(): bad table entry for\n"
+"                              \"input_array_offsets\" parameter!\n"
 "                              error status=%d",
 		   status);
 	return status;					/*** ERROR RETURN ***/
 	}
 
+/**************************************/
+
 /*
  * input array strides
  */
@@ -915,13 +1012,15 @@ else	{
 	CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
 		   __LINE__, __FILE__, CCTK_THORNSTRING,
 "\n"
-"   CCTK_InterpLocalUniform(): bad table entry\n"
-"                              for \"input_array_offsets\" parameter!\n"
+"   CCTK_InterpLocalUniform(): bad table entry for\n"
+"                              \"input_array_offsets\" parameter!\n"
 "                              error status=%d",
 		   status);
 	return status;					/*** ERROR RETURN ***/
 	}
 
+/**************************************/
+
 /*
  * input array min/max subscripts
  */
@@ -942,8 +1041,8 @@ else	{
 	CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
 		   __LINE__, __FILE__, CCTK_THORNSTRING,
 "\n"
-"   CCTK_InterpLocalUniform(): bad table entry\n"
-"                              for \"input_array_min_subscripts\" parameter!\n"
+"   CCTK_InterpLocalUniform(): bad table entry for\n"
+"                              \"input_array_min_subscripts\" parameter!\n"
 "                              error status=%d",
 		   status);
 	return status;					/*** ERROR RETURN ***/
@@ -981,8 +1080,8 @@ else	{
 	CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
 		   __LINE__, __FILE__, CCTK_THORNSTRING,
 "\n"
-"   CCTK_InterpLocalUniform(): bad table entry\n"
-"                              for \"input_array_max_subscripts\" parameter!\n"
+"   CCTK_InterpLocalUniform(): bad table entry for\n"
+"                              \"input_array_max_subscripts\" parameter!\n"
 "                              error status=%d",
 		   status);
 	return status;					/*** ERROR RETURN ***/
@@ -994,8 +1093,8 @@ else	{
  * operand indices
  */
   {
-CCTK_INT *const operand_indices
-	= malloc(N_output_arrays1 * sizeof(CCTK_INT));
+CCTK_INT* const operand_indices
+	= (CCTK_INT*) malloc(N_output_arrays1 * sizeof(CCTK_INT));
 if (operand_indices == NULL)
    then {
 	free(input_array_offsets);
@@ -1059,19 +1158,21 @@ else	{
 	CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
 		   __LINE__, __FILE__, CCTK_THORNSTRING,
 "\n"
-"   CCTK_InterpLocalUniform(): bad table entry\n"
-"                              for \"operand_indices\" parameter!\n"
+"   CCTK_InterpLocalUniform(): bad table entry for\n"
+"                              \"operand_indices\" parameter!\n"
 "                         error status=%d",
 		   status);
 	return status;					/*** ERROR RETURN ***/
 	}
 
+/**************************************/
+
 /*
  * operation_codes
  */
   {
-CCTK_INT *const operation_codes
-	= malloc(N_output_arrays1 * sizeof(CCTK_INT));
+CCTK_INT* const operation_codes
+	= (CCTK_INT*) malloc(N_output_arrays1 * sizeof(CCTK_INT));
 if (operation_codes == NULL)
    then {
 	free(operand_indices);
@@ -1110,8 +1211,8 @@ else	{
 	CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
 		   __LINE__, __FILE__, CCTK_THORNSTRING,
 "\n"
-"   CCTK_InterpLocalUniform(): bad table entry\n"
-"                              for \"operation_codes\" parameter!\n"
+"   CCTK_InterpLocalUniform(): bad table entry for\n"
+"                              \"operation_codes\" parameter!\n"
 "                         error status=%d",
 		   status);
 	return status;					/*** ERROR RETURN ***/
@@ -1223,7 +1324,7 @@ if (status)
    then {
 	/* yes, we're doing Jacobian queries */
 	Jacobian_info.Jacobian_pointer
-	    = (CCTK_REAL **) malloc(N_output_arrays1 * sizeof(CCTK_REAL *));
+	    = (CCTK_REAL**) malloc(N_output_arrays1 * sizeof(CCTK_REAL *));
 	if (Jacobian_info.Jacobian_pointer == NULL)
 	   then {
 		free(operation_codes);
@@ -1232,8 +1333,8 @@ if (status)
 		return UTIL_ERROR_NO_MEMORY;		/*** ERROR RETURN ***/
 		}
 	  {
-	CCTK_POINTER *Jacobian_pointer_temp
-	    = (CCTK_POINTER *) malloc(N_output_arrays1 * sizeof(CCTK_POINTER));
+	CCTK_POINTER* Jacobian_pointer_temp
+	    = (CCTK_POINTER*) malloc(N_output_arrays1 * sizeof(CCTK_POINTER));
 	if (Jacobian_pointer_temp == NULL)
 	   then {
 		free(operation_codes);
@@ -1276,7 +1377,7 @@ if (status)
 	  }
 
 	Jacobian_info.Jacobian_offset
-		= (CCTK_INT *) malloc(N_output_arrays1 * sizeof(CCTK_INT));
+		= (CCTK_INT*) malloc(N_output_arrays1 * sizeof(CCTK_INT));
 	if (Jacobian_info.Jacobian_offset == NULL)
 	   then {
 		free(operation_codes);
@@ -1550,7 +1651,7 @@ if (p_interp_fn == NULL_INTERP_FN_PTR)
 "        interp_operator=\"%s\", N_dims=%d\n"
 "        molecule_family=\"%s\", order=%d, smoothing=%d",
 		   interp_operator_string, N_dims,
-		   molecule_family_string, order, smoothing);
+		   molecule_family_str, order, smoothing);
 	return UTIL_ERROR_BAD_INPUT;			/*** ERROR RETURN ***/
 	}
 
@@ -1564,7 +1665,9 @@ const int return_code
 	= (*p_interp_fn)(coord_origin, coord_delta,
 			 N_interp_points,
 			    interp_coords_type_code, interp_coords,
-			    out_of_range_tolerance,
+			    N_boundary_points_to_omit,
+			    boundary_off_centering_tolerance,
+			    boundary_extrapolation_tolerance,
 			 N_input_arrays,
 			    input_array_offsets, input_array_strides,
 			    input_array_min_subscripts,
diff --git a/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.h b/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.h
index d2990df..851378a 100644
--- a/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.h
+++ b/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.h
@@ -70,11 +70,17 @@ typedef int bool;
  */
 #define MAX_N_DIMS	3
 
+/* a "boundary" is the combination of a dimension and a min/max "side" */
+#define MAX_N_BOUNDARIES	(2*MAX_N_DIMS)
+
 /*
  * if molecule_family_string is a C string specifying a molecule family
  * (i.e. value associated with the "molecule_family" key in the parameter
  * table), we must have
  *	strlen(molecule_family_string) <= MAX_MOLECULE_FAMILY_STRLEN
+ * n.b. exceeding this won't cause a buffer overflow, it will "just"
+ *      cause the string to be truncated (and probably not recognized
+ *      by the interpolator)
  */
 #define MAX_MOLECULE_FAMILY_STRLEN	20
 
@@ -103,12 +109,15 @@ enum	molecule_family
  * other compile-time settings
  */
 
-/* default for each element of out_of_range_tolerance[] */
-/* FIXME: this is based on C "double", which may not match CCTK_REAL */
-#define OUT_OF_RANGE_TOLERANCE_DEFAULT	(10000.0*DBL_EPSILON)
+/* default for each element of boundary_off_centering_tolerance[] */
+#define BOUNDARY_OFF_CENTER_TOL_DEFAULT		999.0
+
+/* default for each element of boundary_extrapolation_tolerance[] */
+#define BOUNDARY_EXTRAP_TOL_DEFAULT		1.0e-10
 
-/* CCTK_VWarn() severity level for error messages */
-#define ERROR_MSG_SEVERITY_LEVEL	1
+/* CCTK_VWarn() severity level for error/warning messages */
+#define ERROR_MSG_SEVERITY_LEVEL	0
+#define WARNING_MSG_SEVERITY_LEVEL	1
 
 /******************************************************************************/
 
@@ -122,6 +131,7 @@ enum	molecule_family
 struct	error_flags
 	{
 	int error_pt;
+	int error_ibndry;
 	int error_axis;
 	int error_end;
 	};
@@ -153,20 +163,16 @@ struct	Jacobian_info
 
 /*
  * error codes for LocalInterp_molecule_posn()
- * ... these are defined in terms of INT_MIN from <limits.h>
  */
 
 /* x < minimum allowable x in grid */
-#define MOLECULE_POSN_ERROR_X_LT_MIN	(INT_MIN + 0)
+#define MOLECULE_POSN_ERROR_X_LT_MIN	(-1)
 
 /* x > maximum allowable x in grid */
-#define MOLECULE_POSN_ERROR_X_GT_MAX	(INT_MIN + 1)
+#define MOLECULE_POSN_ERROR_X_GT_MAX	(-2)
 
 /* grid is smaller than molecule */
-#define MOLECULE_POSN_ERROR_GRID_TINY	(INT_MIN + 2)
-
-/* is a given integer an error code? */
-#define IS_MOLECULE_POSN_ERROR_CODE(x)	(x <= MOLECULE_POSN_ERROR_GRID_TINY)
+#define MOLECULE_POSN_ERROR_GRID_TINY	(-3)
 
 /******************************************************************************/
 
@@ -216,11 +222,12 @@ int LocalInterp_ILU_Hermite(int N_dims,
 int LocalInterp_molecule_posn(fp grid_origin, fp grid_delta,
 			      int grid_i_min, int grid_i_max,
 			      int molecule_size,
-			      fp out_of_range_tolerance_min,
-			      fp out_of_range_tolerance_max,
+			      fp boundary_off_centering_tolerance_min,
+			      fp boundary_off_centering_tolerance_max,
+			      fp boundary_extrapolation_tolerance_min,
+			      fp boundary_extrapolation_tolerance_max,
 			      fp x,
-			      fp *x_rel,
-			      int *molecule_m_min, int *molecule_m_max);
+			      int* i_center, fp* x_rel);
 
 /* functions in util.c */
 void LocalInterp_zero_int_array(int N, CCTK_INT array[]);
diff --git a/src/GeneralizedPolynomial-Uniform/README b/src/GeneralizedPolynomial-Uniform/README
index 0963854..b022489 100644
--- a/src/GeneralizedPolynomial-Uniform/README
+++ b/src/GeneralizedPolynomial-Uniform/README
@@ -7,7 +7,8 @@ under the name
 
 
 The source code files are as follows:
-* startup.c registers the interpolation operator
+* startup.c registers the interpolation operators
+* InterpLocalUniform.h is an overall header file for the whole interpolator
 * InterpLocalUniform.c is the top-level driver: it gets the various
   options from the parameter table, then decodes
 	(N_dims, molecule_family, order, smoothing)
@@ -16,6 +17,7 @@ The source code files are as follows:
 * [123]d.cube.order*.smooth*.c define the individual interpolation
   subfunctions.  Each of them just #defines a whole bunch of macros,
   then #includes template.c (which has the actual code).
+* template.h defines a prototype for the function defined by template.c
 * template.c is the actual interpolation code.  It is written in
   terms of a large number of macros, which should be #defined
   before #including template.c.  There's a long block comment
@@ -26,6 +28,10 @@ The source code files are as follows:
   Maple-generated code fragments in the [123]d.coeffs/directories;
   template.c uses various macros to tell it which fragments to
   #include.
+* molecule_posn.c contains the  LocalInterp_molecule_posn()  function
+  to compute where in the grid each (an) interpolation molecule should
+  be centered for each (a given) interpolation point.
+* util.c contains some low-level utility routines
 * [123]d.maple are the top-level Maple code files; they call
   various functions in interpolate.maple and util.maple to do
   the actual work.
@@ -33,7 +39,19 @@ The source code files are as follows:
   computing an interpolant and manipulating/printing it in various
   ways
 * util.maple contains low-level utility routines
+* makefile is a makefile for...
+* test_molecule_posn.c is a standalone test driver for the code
+  in  molecule_posn.c .  By default it runs a large set of (around 100)
+  test cases stored in a table in the code.
 
+The subdirectories are as follows:
+* [123]d.coeffs/ are empty directory trees left-over from previous
+  versions of this interpolator (CVS doesn't have any easy way to
+  remove directories)
+* Lagrange/ contains the code for the Lagrange polynomial interpolator.
+* Hermite/ contains the code for the Hermite polynomial interpolator.
+* common/ contains low-level code common to both the Lagrange and
+  the Hermite interpolators.
 
 
 To add a new combination of (N_dims, molecule_family, order, smoothing),
@@ -60,17 +78,10 @@ to this interpolator, you need to
 
 
 
-Other makefile targets:
-test_molecule_posn
-	This makes a standalone test driver for the
-	   LocalInterp_molecule_posn()
-	function defined in molecule_posn.c
-
-
 Copyright
 =========
 
-This interpolator is copyright (C) 2002-2003
+This interpolator is copyright (C) 2001-2003
 by Jonathan Thornburg <jthorn@aei.mpg.de>.
 
 This interpolator is free software; you can redistribute it and/or modify
diff --git a/src/GeneralizedPolynomial-Uniform/molecule_posn.c b/src/GeneralizedPolynomial-Uniform/molecule_posn.c
index e007396..6d42ef3 100644
--- a/src/GeneralizedPolynomial-Uniform/molecule_posn.c
+++ b/src/GeneralizedPolynomial-Uniform/molecule_posn.c
@@ -8,7 +8,6 @@
  @@*/
 
 #include <math.h>
-#include <limits.h>
 #include <stdlib.h>
 #include <stdio.h>
 
@@ -73,6 +72,18 @@ static const char *rcsid = "$Header$";
 	center, i.e. the above diagram has columns labelled with [i+m].
  @enddesc
 
+ @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,
+           also change to returning status code,
+           also drop returning @var min_m and @var max_m
+           since they were never used.
+ @endhdesc 
+
  @var	grid_origin
  @vdesc	The floating-point coordinate x of the grid point i=0.
  @vtype	fp grid_origin
@@ -98,50 +109,53 @@ static const char *rcsid = "$Header$";
  @vtype	int molecule_size
  @endvar
 
- @var	out_of_range_tolerance_min, out_of_range_tolerance_max
- @vdesc	Specifies how out-of-range interpolation points should
-	be handled for the {minimum,maximum} ends of the grid
-	respectively.
-	If out_of_range_tolerance >= 0.0,
-	   then an interpolation point is considered to be
-		"out of range" if and only if the interpolation
-		point is > out_of_range_tolerance * grid_delta
-		outside the grid in any coordinate.
-	If out_of_range_tolerance == -1.0,
-	   then an interpolation point is considered to be
-		"out of range" if and only if a centered molecule
-		would require data from outside the grid.
-	Other values of out_of_range_tolerance are illegal.
- @vtype fp out_of_range_tolerance_min, out_of_range_tolerance_max
+ @var	boundary_off_centering_tolerance_min,
+	boundary_off_centering_tolerance_max
+ @vdesc	Specifies how many grid spacings the interpolation point
+	may be beyond the default-centering region before being
+	declared "out of range" on the {minimum,maximum} boundaries
+	of the grid respectively.
+ @vtype fp boundary_off_centering_tolerance_min,
+	   boundary_off_centering_tolerance_max
+ @endvar
+
+ @var	boundary_extrapolation_tolerance_min,
+	boundary_extrapolation_tolerance_max
+ @vdesc	Specifies how many grid spacings the interpolation point
+	may be beyond the grid boundary before being declared
+	"out of range" on the {minimum,maximum} boundaries
+	of the grid respectively.
+ @vtype fp boundary_extrapolation_tolerance_min,
+	   boundary_extrapolation_tolerance_max
  @endvar
 
  @var	x
- @vdesc	The floating-point coordinate x at which we would like to center
-	the molecule.
+ @vdesc	The floating-point coordinate of the interpolation point.
  @vtype	fp x
  @endvar
 
- @var	x_rel
- @vdesc	A pointer to a floating-point value where this function should
-	store the x coordinate relative to the molecule center, measured
-	in units of the grid spacing; or NULL to skip storing this.
- @vtype	fp *x_rel
+ @var	i_center
+ @vdesc	A pointer to an value where this function should
+	store the integer coordinate of the molecule center,
+	or NULL to skip storing this
+ @vtype	int *i_center
  @vio	pointer to out
  @endvar
 
- @var	molecule_m_min, molecule_m_max
- @vdesc	A pointer to an int where this function should store the
-	{minimum,maximum} molecule coordinate m of the molecule;
+ @var	x_rel
+ @vdesc	A pointer to where this function should store the
+	interpolation point's x coordinate relative to the
+	molecule center, measured in units of the grid spacing;
 	or NULL to skip storing this.
- @vtype	int *molecule_m_min, *molecule_m_max
+ @vtype	fp *x_rel
  @vio	pointer to out
  @endvar
 
  @returntype	int
  @returndesc
-	This function returns the integer coordinate i_center at which
-	the molecule should be centered, or one of the error codes defined
-	in "InterpLocalUniform.h":
+	This function returns 0 if the interpolation point is "in range",
+	or one of the (negative) error codes defined in "InterpLocalUniform.h"
+	if the interpolation point is "out of range":
 	MOLECULE_POSN_ERROR_X_LT_MIN
 		if  x  is out-of-range on the min end of the grid
 		(i.e. x < the minimum allowable x)
@@ -155,73 +169,70 @@ static const char *rcsid = "$Header$";
 int LocalInterp_molecule_posn(fp grid_origin, fp grid_delta,
 			      int grid_i_min, int grid_i_max,
 			      int molecule_size,
-			      fp out_of_range_tolerance_min,
-			      fp out_of_range_tolerance_max,
+			      fp boundary_off_centering_tolerance_min,
+			      fp boundary_off_centering_tolerance_max,
+			      fp boundary_extrapolation_tolerance_min,
+			      fp boundary_extrapolation_tolerance_max,
 			      fp x,
-			      fp *x_rel,
-			      int *molecule_m_min, int *molecule_m_max)
+			      int* i_center, fp* x_rel)
 {
-/* min/max molecule coordinates */
-const int m_max = (molecule_size >> 1);		/* a positive number */
-const int m_min = m_max - molecule_size + 1;	/* a negative number */
+/* molecule radia (inherently positive) in +/- directions */
+const int mr_plus  = (molecule_size >> 1);
+const int mr_minus = molecule_size - mr_plus - 1;
 
-/* integer coordinate i of point x, as a floating-point number */
-const fp fp_i = (x - grid_origin) / grid_delta;
+/* 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;
 
-/* is point x out-of-range? */
-if (out_of_range_tolerance_min >= 0.0)
-   then {
-	const fp fp_effective_grid_i_min
-		= ((fp) grid_i_min) - out_of_range_tolerance_min;
-	if (fp_i < fp_effective_grid_i_min)
-	   then return INT_MIN;				/*** ERROR RETURN ***/
-	}
-if (out_of_range_tolerance_max >= 0.0)
-   then {
-	const fp fp_effective_grid_i_max
-		= ((fp) grid_i_max) + out_of_range_tolerance_max;
-	if (fp_i > fp_effective_grid_i_max)
-	   then return INT_MAX;				/*** ERROR RETURN ***/
-	}
+/* 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;
 
-/* integer coordinate i where we will center the molecule */
-/* (see diagram in header comment for explanation) */
-/* ... as a floating-point number */
-  {
-const fp fp_i_center = IS_EVEN(molecule_size) ? floor(fp_i) : JT_ROUND(fp_i);
-/* ... as an integer */
-int i_center = (int) fp_i_center;
+/* integer coordinate i of interpolation point, as a floating-point number */
+const fp fp_i = (x - grid_origin) / grid_delta;
 
-/* min/max integer coordinates in molecule assuming it's fully centered */
-const int centered_i_min = i_center + m_min;
-const int centered_i_max = i_center + m_max;
+/* 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 ***/
 
-/* check if off-centered molecules are forbidden? */
-if ((out_of_range_tolerance_min == -1.0) && (centered_i_min < grid_i_min))
+/* 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 ***/
-if ((out_of_range_tolerance_max == -1.0) && (centered_i_max > grid_i_max))
+if (fp_i > (fp)grid_i_max + boundary_extrapolation_tolerance_max)
    then return MOLECULE_POSN_ERROR_X_GT_MAX;		/*** ERROR RETURN ***/
 
-/* off-center as needed if we're close to the edge of the grid */
+/* is interpolation point x beyond the off-centering tolerance? */
+if (fp_i < fp_centered_min_possible_i - boundary_off_centering_tolerance_min)
+   then return MOLECULE_POSN_ERROR_X_LT_MIN;		/*** ERROR RETURN ***/
+if (fp_i > fp_centered_max_possible_i + boundary_off_centering_tolerance_max)
+   then return MOLECULE_POSN_ERROR_X_GT_MAX;		/*** ERROR RETURN ***/
+
+/* now choose the actual molecule position/centering: */
+/* first set up a centered molecule */
   {
-int shift = 0;
-if (molecule_size > HOW_MANY_IN_RANGE(grid_i_min,grid_i_max))
-   then return MOLECULE_POSN_ERROR_GRID_TINY;		/*** ERROR RETURN ***/
-if (centered_i_min < grid_i_min)
-   then shift = grid_i_min - centered_i_min;	/* a positive number */
-if (centered_i_max > grid_i_max)
-   then shift = grid_i_max - centered_i_max;	/* a negative number */
-i_center += shift;
+fp fp_i_center = IS_EVEN(molecule_size)	/* ... as a floating-point number */
+		 ? floor(fp_i)
+		 : JT_ROUND(fp_i);
+int int_i_center = (int) fp_i_center;	/* ... as an integer */
+
+/* then clamp the molecule at the grid boundaries */
+if (int_i_center - mr_minus < grid_i_min)
+   then {
+	int_i_center = grid_i_min + mr_minus;
+	fp_i_center = (fp) int_i_center;
+	}
+if (int_i_center + mr_plus > grid_i_max)
+   then {
+	int_i_center = grid_i_max - mr_plus;
+	fp_i_center = (fp) int_i_center;
+	}
 
 /* store the results */
+if (i_center != NULL)
+   then *i_center = int_i_center;
 if (x_rel != NULL)
-   then *x_rel = fp_i - ((fp) i_center);
-if (molecule_m_min != NULL)
-   then *molecule_m_min = m_min;
-if (molecule_m_max != NULL)
-   then *molecule_m_max = m_max;
+   then *x_rel = fp_i - fp_i_center;
 
-return i_center;
-  }
+return 0;						/*** NORMAL RETURN ***/
   }
 }
diff --git a/src/GeneralizedPolynomial-Uniform/template.c b/src/GeneralizedPolynomial-Uniform/template.c
index 0bb409b..a567119 100644
--- a/src/GeneralizedPolynomial-Uniform/template.c
+++ b/src/GeneralizedPolynomial-Uniform/template.c
@@ -20,6 +20,13 @@
 	    3-D molecules.
   @endhdesc
 
+  @hdate    28.Jan.2003
+  @hauthor  Jonathan Thornburg <jthorn@aei.mpg.de>
+  @hdesc    Support parameter-table entries
+            @var boundary_off_centering_tolerance and
+            @var boundary_extrapolation_tolerance.
+  @endhdesc 
+
   @desc
 	The following header files must be #included before
 	#including this file:
@@ -301,6 +308,9 @@
 	only briefly; see the  InterpLocalUniform()  documentation
 	and/or the thorn guide for this thorn for further details.
 
+	If you change the arguments for this function, note that you
+	must also change the prototype in "template.h".
+
   @var		error_flags
   @vdesc	If we encounter an out-of-range point and this pointer
 		is non-NULL, we store a description of the out-of-range
@@ -346,13 +356,15 @@
 
   @@*/
 int FUNCTION_NAME(/***** coordinate system *****/
-		  const CCTK_REAL coord_system_origin[],
-		  const CCTK_REAL grid_spacing[],
+		  const CCTK_REAL coord_origin[],
+		  const CCTK_REAL coord_delta[],
 		  /***** interpolation points *****/
 		  int N_interp_points,
 		  int interp_coords_type_code,
 		  const void* const interp_coords[],
-		  const CCTK_REAL out_of_range_tolerance[],
+		  const CCTK_INT N_boundary_points_to_omit[],
+		  const CCTK_REAL boundary_off_centering_tolerance[],
+		  const CCTK_REAL boundary_extrapolation_tolerance[],
 		  /***** input arrays *****/
 		  int N_input_arrays,
 		  const CCTK_INT input_array_offsets[],
@@ -691,8 +703,8 @@ int out;
  */
   {
 #if N_DIMS >= 1
-  const fp origin_x = coord_system_origin[X_AXIS];
-  const fp delta_x = grid_spacing[X_AXIS];
+  const fp origin_x = coord_origin[X_AXIS];
+  const fp delta_x = coord_delta[X_AXIS];
   #if    defined(HAVE_OP_DX) \
       || defined(HAVE_OP_DXY) || defined(HAVE_OP_DXZ) \
       || defined(HAVE_OP_DXX)
@@ -700,8 +712,8 @@ int out;
   #endif
 #endif
 #if N_DIMS >= 2
-  const fp origin_y = coord_system_origin[Y_AXIS];
-  const fp delta_y = grid_spacing[Y_AXIS];
+  const fp origin_y = coord_origin[Y_AXIS];
+  const fp delta_y = coord_delta[Y_AXIS];
   #if    defined(HAVE_OP_DY) \
       || defined(HAVE_OP_DXY) || defined(HAVE_OP_DYZ) \
       || defined(HAVE_OP_DYY)
@@ -709,8 +721,8 @@ int out;
   #endif
 #endif
 #if N_DIMS >= 3
-  const fp origin_z = coord_system_origin[Z_AXIS];
-  const fp delta_z = grid_spacing[Z_AXIS];
+  const fp origin_z = coord_origin[Z_AXIS];
+  const fp delta_z = coord_delta[Z_AXIS];
   #if    defined(HAVE_OP_DZ) \
       || defined(HAVE_OP_DXZ) || defined(HAVE_OP_DYZ) \
       || defined(HAVE_OP_DZZ)
@@ -764,6 +776,7 @@ int out;
 /*
  * interpolate at each point
  */
+int return_code = 0;
 int pt;
     for (pt = 0 ; pt < N_interp_points ; ++pt)
     {
@@ -890,124 +903,84 @@ int pt;
      *
      * n.b. we need the final answers in variables with the magic
      *      names (x,y,z) (machine-generated code uses these names),
-     *      but we use temp variables as intermediates for (likely)
-     *	better performance: the temp variables have their addresses
-     *	taken and so may not be register-allocated, whereas the
-     *	final (x,y,z) are "clean" and thus more likely to be
-     *      register-allocated
+     *      but we use temp variables as intermediates for these and
+     *      for  center_[ijk] for (likely) better performance:
+     *      the temp variables have their addresses taken and so
+     *      may not be register-allocated, whereas the final variables
+     *	    are "clean" and thus more likely to be register-allocated
      */
       {
-    #if (N_DIMS >= 1)
-      fp x_temp;
-      const int center_i
-	 = LocalInterp_molecule_posn(origin_x, delta_x,
-				     input_array_min_subscripts[X_AXIS],
-				     input_array_max_subscripts[X_AXIS],
-				     MOLECULE_SIZE,
-				     out_of_range_tolerance[X_AXIS_MIN],
-				     out_of_range_tolerance[X_AXIS_MAX],
-				     interp_coords_fp[X_AXIS],
-				     &x_temp,
-				     (int *) NULL, (int *) NULL);
-      const fp x = x_temp;
-    #endif
-    #if (N_DIMS >= 2)
-      fp y_temp;
-      const int center_j
-	 = LocalInterp_molecule_posn(origin_y, delta_y,
-    				     input_array_min_subscripts[Y_AXIS],
-    				     input_array_max_subscripts[Y_AXIS],
-    				     MOLECULE_SIZE,
-    				     out_of_range_tolerance[Y_AXIS_MIN],
-    				     out_of_range_tolerance[Y_AXIS_MAX],
-    				     interp_coords_fp[Y_AXIS],
-    				     &y_temp,
-    				     (int *) NULL, (int *) NULL);
-      const fp y = y_temp;
-    #endif
-    #if (N_DIMS >= 3)
-      fp z_temp;
-      const int center_k
-	 = LocalInterp_molecule_posn(origin_z, delta_z,
-    				     input_array_min_subscripts[Z_AXIS],
-    				     input_array_max_subscripts[Z_AXIS],
-    				     MOLECULE_SIZE,
-    				     out_of_range_tolerance[Z_AXIS_MIN],
-    				     out_of_range_tolerance[Z_AXIS_MAX],
-    				     interp_coords_fp[Z_AXIS],
-    				     &z_temp,
-    				     (int *) NULL, (int *) NULL);
-      const fp z = z_temp;
-    #endif
-    #if (N_DIMS > 3)
-      #error "N_DIMS may not be > 3!"
-    #endif
-
-    /* is the interpolation point out-of-range? */
-    #if (N_DIMS >= 1)
-      if (IS_MOLECULE_POSN_ERROR_CODE(center_i))
-	 then {
-	      if (center_i == MOLECULE_POSN_ERROR_GRID_TINY)
-		 then return CCTK_ERROR_INTERP_GRID_TOO_TINY;
-	      if (error_flags != NULL)
-		 then {
-		      error_flags->error_pt = pt;
-		      error_flags->error_axis = X_AXIS;
-		      error_flags->error_end
-			 = (center_i == MOLECULE_POSN_ERROR_X_GT_MAX) ? +1 : -1;
-		      }
-	      return CCTK_ERROR_INTERP_POINT_X_RANGE; /*** ERROR RETURN ***/
-	      }
-    #endif
-    #if (N_DIMS >= 2)
-      if (IS_MOLECULE_POSN_ERROR_CODE(center_j))
-	 then {
-	      if (center_j == MOLECULE_POSN_ERROR_GRID_TINY)
-		 then return CCTK_ERROR_INTERP_GRID_TOO_TINY;
-	      if (error_flags != NULL)
-		 then {
-		      error_flags->error_pt = pt;
-		      error_flags->error_axis = Y_AXIS;
-		      error_flags->error_end
-			 = (center_j == MOLECULE_POSN_ERROR_X_GT_MAX) ? +1 : -1;
-		      }
-	      return CCTK_ERROR_INTERP_POINT_X_RANGE; /*** ERROR RETURN ***/
-	      }
-    #endif
-    #if (N_DIMS >= 3)
-      if (IS_MOLECULE_POSN_ERROR_CODE(center_k))
-	 then {
-	      if (center_k == MOLECULE_POSN_ERROR_GRID_TINY)
-		 then return CCTK_ERROR_INTERP_GRID_TOO_TINY;
-	      if (error_flags != NULL)
-		 then {
-		      error_flags->error_pt = pt;
-		      error_flags->error_axis = Z_AXIS;
-		      error_flags->error_end
-			 = (center_k == MOLECULE_POSN_ERROR_X_GT_MAX) ? +1 : -1;
-		      }
-	      return CCTK_ERROR_INTERP_POINT_X_RANGE; /*** ERROR RETURN ***/
-	      }
-    #endif
+    int center_ijk_temp[MAX_N_DIMS];
+    fp  xyz_temp[MAX_N_DIMS];
+    int axis;
+	for (axis = 0 ; axis < N_DIMS ; ++axis)
+	{
+	const int ibndry_min = 2*axis;
+	const int ibndry_max = 2*axis + 1;
+	const int status = LocalInterp_molecule_posn
+				(coord_origin[axis], coord_delta[axis],
+				 input_array_min_subscripts[axis]
+				   + N_boundary_points_to_omit[ibndry_min],
+				 input_array_max_subscripts[axis]
+				   - N_boundary_points_to_omit[ibndry_max],
+				 MOLECULE_SIZE,
+				 boundary_off_centering_tolerance[ibndry_min],
+				 boundary_off_centering_tolerance[ibndry_max],
+				 boundary_extrapolation_tolerance[ibndry_min],
+				 boundary_extrapolation_tolerance[ibndry_max],
+				 interp_coords_fp[axis],
+				 &center_ijk_temp[axis], &xyz_temp[axis]);
+	if (status < 0)
+	   then {
+		if (status == MOLECULE_POSN_ERROR_GRID_TINY)
+		   then return CCTK_ERROR_INTERP_GRID_TOO_TINY;
+							/*** ERROR RETURN ***/
+		if (error_flags != NULL)
+		   then {
+			error_flags->error_pt = pt;
+			error_flags->error_axis = axis;
+			if (status == MOLECULE_POSN_ERROR_X_LT_MIN)
+			   then {
+				error_flags->error_ibndry = ibndry_min;
+				error_flags->error_end = -1;
+				}
+			   else {
+				error_flags->error_ibndry = ibndry_max;
+				error_flags->error_end = +1;
+				}
+			}
+		return CCTK_ERROR_INTERP_POINT_OUTSIDE;	/*** ERROR RETURN ***/
+		}
+	}
+      {
+  #if (N_DIMS >= 1)
+    const int center_i = center_ijk_temp[X_AXIS];
+    const fp  x        = xyz_temp       [X_AXIS];
+  #endif
+  #if (N_DIMS >= 2)
+    const int center_j = center_ijk_temp[Y_AXIS];
+    const fp  y        = xyz_temp       [Y_AXIS];
+  #endif
+  #if (N_DIMS >= 3)
+    const int center_k = center_ijk_temp[Z_AXIS];
+    const fp  z        = xyz_temp       [Z_AXIS];
+  #endif
 
 
     /*
      * compute 1-d position of molecule center in input arrays
      */
       {
-    #if (N_DIMS == 1)
+    #if   (N_DIMS == 1)
       const int molecule_center_posn =   stride_i*center_i;
-    #endif
-    #if (N_DIMS == 2)
+    #elif (N_DIMS == 2)
       const int molecule_center_posn =   stride_i*center_i
 				       + stride_j*center_j;
-    #endif
-    #if (N_DIMS == 3)
+    #elif (N_DIMS == 3)
       const int molecule_center_posn =   stride_i*center_i
 				       + stride_j*center_j
 				       + stride_k*center_k;
-    #endif
-    #if (N_DIMS > 3)
+    #else
       #error "N_DIMS may not be > 3!"
     #endif
 
@@ -1566,6 +1539,7 @@ CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
 	  }
       }
       }
+      }
 
     /* end of  for (pt = ...)  loop */
     }
diff --git a/src/GeneralizedPolynomial-Uniform/template.h b/src/GeneralizedPolynomial-Uniform/template.h
index 3bc2a13..b5eb5b4 100644
--- a/src/GeneralizedPolynomial-Uniform/template.h
+++ b/src/GeneralizedPolynomial-Uniform/template.h
@@ -8,13 +8,15 @@
   @version  $Header$
   @@*/
 int FUNCTION_NAME(/***** coordinate system *****/
-		  const CCTK_REAL coord_system_origin[],
-		  const CCTK_REAL grid_spacing[],
+		  const CCTK_REAL coord_origin[],
+		  const CCTK_REAL coord_delta[],
 		  /***** interpolation points *****/
 		  int N_interp_points,
 		  int interp_coords_type_code,
 		  const void* const interp_coords[],
-		  const CCTK_REAL out_of_range_tolerance[],
+		  const CCTK_INT N_boundary_points_to_omit[],
+		  const CCTK_REAL boundary_off_centering_tolerance[],
+		  const CCTK_REAL boundary_extrapolation_tolerance[],
 		  /***** input arrays *****/
 		  int N_input_arrays,
 		  const CCTK_INT input_array_offsets[],
diff --git a/src/GeneralizedPolynomial-Uniform/test_molecule_posn.c b/src/GeneralizedPolynomial-Uniform/test_molecule_posn.c
index bd3a1f9..4f6a5f5 100644
--- a/src/GeneralizedPolynomial-Uniform/test_molecule_posn.c
+++ b/src/GeneralizedPolynomial-Uniform/test_molecule_posn.c
@@ -7,11 +7,27 @@
  * Usage:
  *	test_molecule_posn			# run a preset set of tests
  *	test_molecule_posn molecule_size		\
- *	                   out_of_range_tolerance_min	\
- *	                   out_of_range_tolerance_max	\
+ *	                   boundary_off_centering_tolerance_min	\
+ *	                   boundary_off_centering_tolerance_max	\
+ *	                   boundary_extrapolation_tolerance_min	\
+ *	                   boundary_extrapolation_tolerance_max	\
+ *	                   boundary_off_centering_tolerance_max	\
  *	                   x			# do a single test as specified
  */
 
+const char* help_msg =
+"usage:\n"
+"# run a preset series of tests:\n"
+"   test_molecule_posn\n"
+"# run a single test as specified:\n"
+"   test_molecule_posn molecule_size \\\n"
+"                      boundary_off_centering_tolerance_min \\\n"
+"                      boundary_off_centering_tolerance_max \\\n"
+"                      boundary_extrapolation_tolerance_min \\\n"
+"                      boundary_extrapolation_tolerance_max \\\n"
+"                      x\n"
+;
+
 #include <math.h>
 #include <string.h>
 #include <limits.h>
@@ -39,8 +55,10 @@ const int grid_i_max = 105;	/* x_max = 13.6 */
 */
 static
   void run_interactive_test(int molecule_size,
-			    fp out_of_range_tolerance_min,
-			    fp out_of_range_tolerance_max,
+			    fp boundary_off_centering_tolerance_min,
+			    fp boundary_off_centering_tolerance_max,
+			    fp boundary_extrapolation_tolerance_min,
+			    fp boundary_extrapolation_tolerance_max,
 			    fp x);
 static
   int run_batch_tests(void);
@@ -56,152 +74,135 @@ struct	args_results
 	{
 	/* args */
 	int molecule_size;
-	fp out_of_range_tolerance_min;
-	fp out_of_range_tolerance_max;
+	fp boundary_off_centering_tolerance_min;
+	fp boundary_off_centering_tolerance_max;
+	fp boundary_extrapolation_tolerance_min;
+	fp boundary_extrapolation_tolerance_max;
 	fp x;
 	/* results */
+	int status;
 	int i_center;
 	fp x_rel;
-	int m_min, m_max;
 	};
 
 /* test data */
+#define OK		0
+#define X_LT_MIN	MOLECULE_POSN_ERROR_X_LT_MIN
+#define X_GT_MAX	MOLECULE_POSN_ERROR_X_GT_MAX
 static const struct args_results test_data[] =
   {
-    /* molecule size 2 */
-    { 2,  1.0,      1.0e-12,  7.19, INT_MIN, -1.1,       0,+1 },
-    { 2,  1.0,      1.0e-12,  7.21,      42, -0.9,       0,+1 },
-    { 2,  1.0,      1.0e-12,  7.24,      42, -0.6,       0,+1 },
-    { 2,  1.0,      1.0e-12,  7.26,      42, -0.4,       0,+1 },
-    { 2,  1.0,      1.0e-12,  7.29,      42, -0.1,       0,+1 },
-    { 2, -1.0,      1.0e-12,  7.29, INT_MIN, -0.1,       0,+1 },
-    { 2,  1.0e-12,  1.0e-12,  7.3,       42,  0.0,       0,+1 },/* grid_x_min */
-    { 2,  1.0e-12,  1.0e-12,  7.31,      42, +0.1,       0,+1 },
-    { 2, -1.0    ,  1.0e-12,  7.31,      42, +0.1,       0,+1 },
-    { 2,  1.0e-12,  1.0e-12,  7.35,      42, +0.5,       0,+1 },
-    { 2,  1.0e-12,  1.0e-12,  7.39,      42, +0.9,       0,+1 },
-    { 2,  1.0e-12,  1.0e-12,  7.41,      43,  0.1,       0,+1 },
-    { 2,  1.0e-12,  1.0e-12,  9.81,      67, +0.1,       0,+1 },
-    { 2,  1.0e-12,  1.0e-12,  9.85,      67, +0.5,       0,+1 },
-    { 2,  1.0e-12,  1.0e-12,  9.89,      67, +0.9,       0,+1 },
-    { 2,  1.0e-12,  1.0e-12, 13.45,     103, +0.5,       0,+1 },
-    { 2,  1.0e-12,  1.0e-12, 13.51,     104, +0.1,       0,+1 },
-    { 2,  1.0e-12,  1.0e-12, 13.55,     104, +0.5,       0,+1 },
-    { 2,  1.0e-12,  1.0e-12, 13.59,     104, +0.9,       0,+1 },
-    { 2,  1.0e-12, -1.0,     13.59,     104, +0.9,       0,+1 },
-    { 2,  1.0e-12,  1.0e-12, 13.6,      104, +1.0,       0,+1 },/* grid_x_max */
-    { 2,  1.0e-12,  1.0,     13.61,     104, +1.1,       0,+1 },
-    { 2,  1.0e-12, -1.0,     13.61, INT_MAX, +1.1,       0,+1 },
-    { 2,  1.0e-12,  1.0,     13.65,     104, +1.5,       0,+1 },
-    { 2,  1.0e-12,  1.0,     13.69,     104, +1.9,       0,+1 },
-    { 2,  1.0e-12,  1.0,     13.71, INT_MAX, +2.1,       0,+1 },
-
-    /* molecule size 3 */
-    { 3,  1.0,      1.0e-12,  7.19, INT_MIN, -2.1,      -1,+1 },
-    { 3,  1.0,      1.0e-12,  7.21,      43, -1.9,      -1,+1 },
-    { 3,  1.0,      1.0e-12,  7.24,      43, -1.6,      -1,+1 },
-    { 3,  1.0,      1.0e-12,  7.26,      43, -1.4,      -1,+1 },
-    { 3,  1.0,      1.0e-12,  7.29,      43, -1.1,      -1,+1 },
-    { 3, -1.0,      1.0e-12,  7.29, INT_MIN, -1.1,      -1,+1 },
-    { 3,  1.0e-12,  1.0e-12,  7.3,       43, -1.0,      -1,+1 },/* grid_x_min */
-    { 3,  1.0e-12,  1.0e-12,  7.31,      43, -0.9,      -1,+1 },
-    { 3,  1.0e-12,  1.0e-12,  7.34,      43, -0.6,      -1,+1 },
-    { 3, -1.0,      1.0e-12,  7.34, INT_MIN, -0.6,      -1,+1 },
-    { 3,  1.0e-12,  1.0e-12,  7.36,      43, -0.4,      -1,+1 },
-    { 3, -1.0,      1.0e-12,  7.36,      43, -0.4,      -1,+1 },
-    { 3,  1.0e-12,  1.0e-12,  7.39,      43, -0.1,      -1,+1 },
-    { 3,  1.0e-12,  1.0e-12,  7.4,       43,  0.0,      -1,+1 },
-    { 3,  1.0e-12,  1.0e-12,  7.44,      43, +0.4,      -1,+1 },
-    { 3,  1.0e-12,  1.0e-12,  7.46,      44, -0.4,      -1,+1 },
-    { 3,  1.0e-12,  1.0e-12,  9.8,       67,  0.0,      -1,+1 },
-    { 3,  1.0e-12,  1.0e-12,  9.81,      67, +0.1,      -1,+1 },
-    { 3,  1.0e-12,  1.0e-12,  9.84,      67, +0.4,      -1,+1 },
-    { 3,  1.0e-12,  1.0e-12,  9.86,      68, -0.4,      -1,+1 },
-    { 3,  1.0e-12,  1.0e-12,  9.89,      68, -0.1,      -1,+1 },
-    { 3,  1.0e-12,  1.0e-12,  9.9,       68,  0.0,      -1,+1 },
-    { 3,  1.0e-12,  1.0e-12, 13.44,     103, +0.4,      -1,+1 },
-    { 3,  1.0e-12,  1.0e-12, 13.46,     104, -0.4,      -1,+1 },
-    { 3,  1.0e-12,  1.0e-12, 13.5,      104,  0.0,      -1,+1 },
-    { 3,  1.0e-12,  1.0e-12, 13.51,     104, +0.1,      -1,+1 },
-    { 3,  1.0e-12,  1.0e-12, 13.54,     104, +0.4,      -1,+1 },
-    { 3,  1.0e-12, -1.0,     13.54,     104, +0.4,      -1,+1 },
-    { 3,  1.0e-12,  1.0e-12, 13.56,     104, +0.6,      -1,+1 },
-    { 3,  1.0e-12, -1.0,     13.56, INT_MAX, +0.6,      -1,+1 },
-    { 3,  1.0e-12,  1.0e-12, 13.59,     104, +0.9,      -1,+1 },
-    { 3,  1.0e-12,  1.0e-12, 13.6,      104, +1.0,      -1,+1 },/* grid_x_max */
-    { 3,  1.0e-12,  1.0,     13.61,     104, +1.1,      -1,+1 },
-    { 3,  1.0e-12,  1.0,     13.65,     104, +1.5,      -1,+1 },
-    { 3,  1.0e-12,  1.0,     13.69,     104, +1.9,      -1,+1 },
-    { 3,  1.0e-12,  1.0,     13.71, INT_MAX, +2.1,      -1,+1 },
-
-    /* molecule size 4 */
-    { 4,  0.2,      1.0e-12,  7.27, INT_MIN, -1.3,      -1,+2 },
-    { 4,  0.2,      1.0e-12,  7.29,      43, -1.1,      -1,+2 },
-    { 4,  1.0e-12,  1.0e-12,  7.3,       43, -1.0,      -1,+2 },/* grid_x_min */
-    { 4,  1.0e-12,  1.0e-12,  7.33,      43, -0.7,      -1,+2 },
-    { 4,  1.0e-12,  1.0e-12,  7.39,      43, -0.1,      -1,+2 },
-    { 4, -1.0,      1.0e-12,  7.39, INT_MIN, -0.1,      -1,+2 },
-    { 4,  1.0e-12,  1.0e-12,  7.4,       43,  0.0,      -1,+2 },
-    { 4, -1.0,      1.0e-12,  7.41,      43, +0.1,      -1,+2 },
-    { 4,  1.0e-12,  1.0e-12,  7.42,      43, +0.2,      -1,+2 },
-    { 4,  1.0e-12,  1.0e-12,  7.48,      43, +0.8,      -1,+2 },
-    { 4,  1.0e-12,  1.0e-12,  7.51,      44, +0.1,      -1,+2 },
-    { 4,  1.0e-12,  1.0e-12,  9.81,      67, +0.1,      -1,+2 },
-    { 4,  1.0e-12,  1.0e-12,  9.85,      67, +0.5,      -1,+2 },
-    { 4,  1.0e-12,  1.0e-12,  9.89,      67, +0.9,      -1,+2 },
-    { 4,  1.0e-12,  1.0e-12, 13.39,     102, +0.9,      -1,+2 },
-    { 4,  1.0e-12,  1.0e-12, 13.41,     103, +0.1,      -1,+2 },
-    { 4,  1.0e-12,  1.0e-12, 13.48,     103, +0.8,      -1,+2 },
-    { 4,  1.0e-12, -1.0,     13.48,     103, +0.8,      -1,+2 },
-    { 4,  1.0e-12,  1.0e-12, 13.5,      103, +1.0,      -1,+2 },
-    { 4,  1.0e-12,  1.0e-12, 13.51,     103, +1.1,      -1,+2 },
-    { 4,  1.0e-12, -1.0,     13.51, INT_MAX, +1.1,      -1,+2 },
-    { 4,  1.0e-12,  1.0e-12, 13.55,     103, +1.5,      -1,+2 },
-    { 4,  1.0e-12,  1.0e-12, 13.59,     103, +1.9,      -1,+2 },
-    { 4,  1.0e-12,  1.0e-12, 13.6,      103, +2.0,      -1,+2 },/* grid_x_max */
-    { 4,  1.0e-12,  2.0,     13.79,     103, +3.9,      -1,+2 },
-    { 4,  1.0e-12,  2.0,     13.81, INT_MAX, +4.1,      -1,+2 },
-
-    /* molecule size 5 */
-    { 5,  3.0,      1.0e-12,  6.99, INT_MIN, -5.1,      -2,+2 },
-    { 5,  3.0,      1.0e-12,  7.01,      44, -4.9,      -2,+2 },
-    { 5,  1.0e-12,  1.0e-12,  7.3,       44, -2.0,      -2,+2 },/* grid_x_min */
-    { 5,  1.0e-12,  1.0e-12,  7.4,       44, -1.0,      -2,+2 },
-    { 5,  1.0e-12,  1.0e-12,  7.44,      44, -0.6,      -2,+2 },
-    { 5, -1.0,      1.0e-12,  7.44, INT_MIN, -0.6,      -2,+2 },
-    { 5,  1.0e-12,  1.0e-12,  7.46,      44, -0.4,      -2,+2 },
-    { 5, -1.0,      1.0e-12,  7.46,      44, -0.4,      -2,+2 },
-    { 5,  1.0e-12,  1.0e-12,  7.49,      44, -0.1,      -2,+2 },
-    { 5,  1.0e-12,  1.0e-12,  7.5,       44,  0.0,      -2,+2 },
-    { 5,  1.0e-12,  1.0e-12,  7.54,      44, +0.4,      -2,+2 },
-    { 5,  1.0e-12,  1.0e-12,  7.56,      45, -0.4,      -2,+2 },
-    { 5,  1.0e-12,  1.0e-12,  7.6,       45,  0.0,      -2,+2 },
-    { 5,  1.0e-12,  1.0e-12,  7.64,      45, +0.4,      -2,+2 },
-    { 5,  1.0e-12,  1.0e-12,  9.8,       67,  0.0,      -2,+2 },
-    { 5,  1.0e-12,  1.0e-12,  9.81,      67, +0.1,      -2,+2 },
-    { 5,  1.0e-12,  1.0e-12,  9.84,      67, +0.4,      -2,+2 },
-    { 5,  1.0e-12,  1.0e-12,  9.86,      68, -0.4,      -2,+2 },
-    { 5,  1.0e-12,  1.0e-12,  9.89,      68, -0.1,      -2,+2 },
-    { 5,  1.0e-12,  1.0e-12,  9.9,       68,  0.0,      -2,+2 },
-    { 5,  1.0e-12,  1.0e-12, 13.34,     102, +0.4,      -2,+2 },
-    { 5,  1.0e-12,  1.0e-12, 13.36,     103, -0.4,      -2,+2 },
-    { 5,  1.0e-12,  1.0e-12, 13.39,     103, -0.1,      -2,+2 },
-    { 5,  1.0e-12,  1.0e-12, 13.41,     103, +0.1,      -2,+2 },
-    { 5,  1.0e-12,  1.0e-12, 13.44,     103, +0.4,      -2,+2 },
-    { 5,  1.0e-12, -1.0,     13.44,     103, +0.4,      -2,+2 },
-    { 5,  1.0e-12,  1.0e-12, 13.46,     103, +0.6,      -2,+2 },
-    { 5,  1.0e-12, -1.0,     13.46, INT_MAX, +0.6,      -2,+2 },
-    { 5,  1.0e-12,  1.0e-12, 13.48,     103, +0.8,      -2,+2 },
-    { 5,  1.0e-12,  1.0e-12, 13.5,      103, +1.0,      -2,+2 },
-    { 5,  1.0e-12,  1.0e-12, 13.51,     103, +1.1,      -2,+2 },
-    { 5,  1.0e-12,  1.0e-12, 13.54,     103, +1.4,      -2,+2 },
-    { 5,  1.0e-12, -1.0,     13.54, INT_MAX, +1.4,      -2,+2 },
-    { 5,  1.0e-12,  1.0e-12, 13.56,     103, +1.6,      -2,+2 },
-    { 5,  1.0e-12, -1.0,     13.56, INT_MAX, +1.6,      -2,+2 },
-    { 5,  1.0e-12,  1.0e-12, 13.59,     103, +1.9,      -2,+2 },
-    { 5,  1.0e-12,  1.0e-12, 13.6,      103, +2.0,      -2,+2 },/* grid_x_max */
-    { 5,  1.0e-12,  1.5,     13.74,     103, +3.4,      -2,+2 },
-    { 5,  1.0e-12,  1.5,     13.76, INT_MAX, +3.6,      -2,+2 },
+    /*** molecule size 2 ***/
+    /*                                        status           */
+    /*     off_centering     extrapolation        i_center     */
+    /*      min     max     min     max    x             x_rel */
+    { 2,  999.0,  999.0,    1.0,  999.0,  7.19, X_LT_MIN,   0,  0.0 },
+    { 2,  999.0,  999.0,    1.0,  999.0,  7.21,       OK,  42, -0.9 },
+    { 2,  999.0,  999.0,    1.0,  999.0,  7.24,       OK,  42, -0.6 },
+    { 2,  999.0,  999.0,    1.0,  999.0,  7.26,       OK,  42, -0.4 },
+    { 2,    0.0,    0.0,    1.0,  999.0,  7.26, X_LT_MIN,   0,  0.0 },
+    { 2, 1.0e-6,  999.0,    0.2,  999.0,  7.29, X_LT_MIN,   0,  0.0 },
+    { 2,    0.2,  999.0, 1.0e-6,  999.0,  7.29, X_LT_MIN,   0,  0.0 },
+    { 2,    0.2,  999.0,    0.2,  999.0,  7.29,       OK,  42, -0.1 },
+    { 2,    0.0,    0.0, 1.0e-6,  999.0,  7.31,       OK,  42, +0.1 },
+    { 2,    0.0,    0.0, 1.0e-6,  999.0,  7.34,       OK,  42, +0.4 },
+    { 2,    0.0,    0.0, 1.0e-6,  999.0,  7.36,       OK,  42, +0.6 },
+    { 2,    0.0,    0.0, 1.0e-6,  999.0,  7.39,       OK,  42, +0.9 },
+    { 2,    0.0,    0.0, 1.0e-6,  999.0,  7.41,       OK,  43, +0.1 },
+    { 2,    0.0,    0.0, 1.0e-6, 1.0e-6,  9.81,       OK,  67, +0.1 },
+    { 2,    0.0,    0.0, 1.0e-6, 1.0e-6,  9.84,       OK,  67, +0.4 },
+    { 2,    0.0,    0.0, 1.0e-6, 1.0e-6,  9.85,       OK,  67, +0.5 },
+    { 2,    0.0,    0.0, 1.0e-6, 1.0e-6,  9.86,       OK,  67, +0.6 },
+    { 2,    0.0,    0.0, 1.0e-6, 1.0e-6,  9.89,       OK,  67, +0.9 },
+    { 2,    0.0,    0.0, 1.0e-6, 1.0e-6, 13.45,       OK, 103, +0.5 },
+    { 2,    0.0,    0.0, 1.0e-6, 1.0e-6, 13.51,       OK, 104, +0.1 },
+    { 2,    0.0,    0.0, 1.0e-6, 1.0e-6, 13.59,       OK, 104, +0.9 },
+    { 2,    0.0,    0.2, 1.0e-6, 1.0e-6, 13.61, X_GT_MAX,   0,  0.0 },
+    { 2,    0.0, 1.0e-6, 1.0e-6, 1.0e-6, 13.61, X_GT_MAX,   0,  0.0 },
+    { 2,    0.0,    0.2,    0.0,    0.2, 13.61,       OK, 104, +1.1 },
+
+    /*** molecule size 3 ***/
+    /*                                        status           */
+    /*     off_centering     extrapolation        i_center     */
+    /*      min     max     min     max    x             x_rel */
+    { 3,    0.5,  999.0,  999.0,  999.0,  7.29, X_LT_MIN,   0,  0.0 },
+    { 3,  999.0,  999.0,    0.0,  999.0,  7.29, X_LT_MIN,   0,  0.0 },
+    { 3,    0.7,  999.0,    0.2,  999.0,  7.29,       OK,  43, -1.1 },
+    { 3,    0.3,  999.0,    0.0,  999.0,  7.31, X_LT_MIN,   0,  0.0 },
+    { 3,    0.5,  999.0,    0.0,  999.0,  7.31,       OK,  43, -0.9 },
+    { 3,    0.0,    0.0,    0.0,    0.0,  7.36,       OK,  43, -0.4 },
+    { 3,    0.0,    0.0,    0.0,    0.0,  7.44,       OK,  43, +0.4 },
+    { 3,    0.0,    0.0,    0.0,    0.0,  7.46,       OK,  44, -0.4 },
+    { 3,    0.0,    0.0,    0.0,    0.0,  9.81,       OK,  67, +0.1 },
+    { 3,    0.0,    0.0,    0.0,    0.0,  9.84,       OK,  67, +0.4 },
+    { 3,    0.0,    0.0,    0.0,    0.0,  9.86,       OK,  68, -0.4 },
+    { 3,    0.0,    0.0,    0.0,    0.0,  9.89,       OK,  68, -0.1 },
+    { 3,    0.0,    0.0,    0.0,    0.0, 13.44,       OK, 103, +0.4 },
+    { 3,    0.0,    0.0,    0.0,    0.0, 13.46,       OK, 104, -0.4 },
+    { 3,    0.0,    0.0,    0.0,    0.0, 13.54,       OK, 104, +0.4 },
+    { 3,    0.0,    0.2,    0.0,    0.0, 13.56,       OK, 104, +0.6 },
+    { 3,  999.0,    0.0,  999.0,  999.0, 13.56, X_GT_MAX,   0,  0.0 },
+    { 3,    0.0,    0.5,    0.0,    0.0, 13.59,       OK, 104, +0.9 },
+    { 3,  999.0,    0.3,  999.0,  999.0, 13.59, X_GT_MAX,   0,  0.0 },
+    { 3,  999.0,    0.5,  999.0,  999.0, 13.61, X_GT_MAX,   0,  0.0 },
+    { 3,  999.0,  999.0,  999.0,    0.0, 13.61, X_GT_MAX,   0,  0.0 },
+    { 3,  999.0,    0.7,  999.0,    0.2, 13.61,       OK, 104, +1.1 },
+    { 3,  999.0,  999.0,  999.0,    0.2, 13.63, X_GT_MAX,   0,  0.0 },
+    { 3,  999.0,    0.7,  999.0,  999.0, 13.63, X_GT_MAX,   0,  0.0 },
+    { 3,  999.0,    0.9,  999.0,    0.4, 13.63,       OK, 104, +1.3 },
+
+    /*** molecule size 4 ***/
+    /*                                        status           */
+    /*     off_centering     extrapolation        i_center     */
+    /*      min     max     min     max    x             x_rel */
+    { 4,    1.0,  999.0,  999.0,  999.0,  7.29, X_LT_MIN,   0,  0.0 },
+    { 4,  999.0,  999.0,    0.0,  999.0,  7.29, X_LT_MIN,   0,  0.0 },
+    { 4,    1.2,  999.0,    0.2,  999.0,  7.29,       OK,  43, -1.1 },
+    { 4,    0.5,  999.0,    0.0,  999.0,  7.31, X_LT_MIN,   0,  0.0 },
+    { 4,    1.0,  999.0,    0.0,  999.0,  7.31,       OK,  43, -0.9 },
+    { 4,    0.0,  999.0,    0.0,  999.0,  7.39, X_LT_MIN,   0,  0.0 },
+    { 4,    0.2,  999.0,    0.0,  999.0,  7.39,       OK,  43, -0.1 },
+    { 4,    0.0,    0.0,    0.0,  999.0,  7.41,       OK,  43, +0.1 },
+    { 4,    0.0,    0.0,    0.0,  999.0,  7.49,       OK,  43, +0.9 },
+    { 4,    0.0,    0.0,    0.0,  999.0,  7.51,       OK,  44, +0.1 },
+    { 4,    0.0,    0.0,    0.0,    0.0,  9.81,       OK,  67, +0.1 },
+    { 4,    0.0,    0.0,    0.0,    0.0,  9.84,       OK,  67, +0.4 },
+    { 4,    0.0,    0.0,    0.0,    0.0,  9.85,       OK,  67, +0.5 },
+    { 4,    0.0,    0.0,    0.0,    0.0,  9.86,       OK,  67, +0.6 },
+    { 4,    0.0,    0.0,    0.0,    0.0,  9.89,       OK,  67, +0.9 },
+    { 4,    0.0,    0.0,    0.0,    0.0, 13.44,       OK, 103, +0.4 },
+    { 4,    0.0,    0.0,    0.0,    0.0, 13.49,       OK, 103, +0.9 },
+    { 4,    0.0,    0.2,    0.0,    0.0, 13.51,       OK, 103, +1.1 },
+    { 4,  999.0,    0.0,  999.0,  999.0, 13.51, X_GT_MAX,   0,  0.0 },
+    { 4,  999.0,    0.8,  999.0,    0.0, 13.59, X_GT_MAX,   0,  0.0 },
+    { 4,    0.0,    1.0,    0.0,    0.0, 13.59,       OK, 103, +1.9 },
+    { 4,  999.0,  999.0,  999.0,    0.0, 13.61, X_GT_MAX,   0,  0.0 },
+    { 4,    0.0,    1.0,  999.0,  999.0, 13.61, X_GT_MAX,   0,  0.0 },
+    { 4,    0.0,    1.2,    0.0,    0.2, 13.61,       OK, 103, +2.1 },
+
+    /*** molecule size 5 ***/
+    /*                                        status           */
+    /*     off_centering     extrapolation        i_center     */
+    /*      min     max     min     max    x             x_rel */
+    { 5,    1.5,  999.0,  999.0,  999.0,  7.29, X_LT_MIN,   0,  0.0 },
+    { 5,  999.0,  999.0,    0.0,  999.0,  7.29, X_LT_MIN,   0,  0.0 },
+    { 5,    1.7,  999.0,    0.2,  999.0,  7.29,       OK,  44, -2.1 },
+    { 5,    0.9,  999.0,    0.0,  999.0,  7.35, X_LT_MIN,   0,  0.0 },
+    { 5,    1.1,  999.0,    0.0,  999.0,  7.35,       OK,  44, -1.5 },
+    { 5,    0.0,  999.0,  999.0,  999.0,  7.44, X_LT_MIN,   0,  0.0 },
+    { 5,    0.2,  999.0,    0.0,  999.0,  7.44,       OK,  44, -0.6 },
+    { 5,    0.0,    0.0,    0.0,    0.0,  7.46,       OK,  44, -0.4 },
+    { 5,    0.0,    0.0,    0.0,    0.0,  9.81,       OK,  67, +0.1 },
+    { 5,    0.0,    0.0,    0.0,    0.0,  9.84,       OK,  67, +0.4 },
+    { 5,    0.0,    0.0,    0.0,    0.0,  9.86,       OK,  68, -0.4 },
+    { 5,    0.0,    0.0,    0.0,    0.0,  9.89,       OK,  68, -0.1 },
+    { 5,    0.0,    0.0,    0.0,    0.0, 13.44,       OK, 103, +0.4 },
+    { 5,  999.0,    0.0,  999.0,  999.0, 13.46, X_GT_MAX,   0,  0.0 },
+    { 5,    0.0,    0.2,    0.0,    0.0, 13.46,       OK, 103, +0.6 },
+    { 5,    0.0,    1.1,    0.0,    0.0, 13.55,       OK, 103, +1.5 },
+    { 5,  999.0,    0.9,  999.0,  999.0, 13.55, X_GT_MAX,   0,  0.0 },
+    { 5,  999.0,  999.0,  999.0,    0.0, 13.61, X_GT_MAX,   0,  0.0 },
+    { 5,  999.0,    1.5,  999.0,  999.0, 13.61, X_GT_MAX,   0,  0.0 },
+    { 5,  999.0,    1.7,  999.0,    0.2, 13.61,       OK, 103, +2.1 },
   };
 #define N_TESTS ((int) (sizeof(test_data)/sizeof(test_data[0])))
 
@@ -212,7 +213,10 @@ int main(int argc, const char *const argv[])
 {
 bool N_fail;
 int molecule_size;
-double out_of_range_tolerance_min, out_of_range_tolerance_max;
+fp boundary_off_centering_tolerance_min;
+fp boundary_off_centering_tolerance_max;
+fp boundary_extrapolation_tolerance_min;
+fp boundary_extrapolation_tolerance_max;
 double x;
 
 switch	(argc)
@@ -222,44 +226,34 @@ case 1:
 	N_fail = run_batch_tests();
 	if (N_fail == 0)
 	   then {
-		printf("*** all tests successful ***\n");
+		printf("*** all %d tests successful ***\n", N_TESTS);
 		return 0;
 		}
 	   else {
-		printf("*** %d test(s) failed ***\n", N_fail);
+		printf("*** %d/%d test(s) failed ***\n", N_fail, N_TESTS);
 		return 1;
 		}
 
-case 5:
+case 7:
 	if (    (sscanf(argv[1], "%d",  &molecule_size) == 1)
-	     && (sscanf(argv[2], "%lf", &out_of_range_tolerance_min) == 1)
-	     && (sscanf(argv[3], "%lf", &out_of_range_tolerance_max) == 1)
-	     && (sscanf(argv[4], "%lf", &x) == 1)    )
+	     && (sscanf(argv[2], "%lf", &boundary_off_centering_tolerance_min) == 1)
+	     && (sscanf(argv[3], "%lf", &boundary_off_centering_tolerance_max) == 1)
+	     && (sscanf(argv[4], "%lf", &boundary_extrapolation_tolerance_min) == 1)
+	     && (sscanf(argv[5], "%lf", &boundary_extrapolation_tolerance_max) == 1)
+	     && (sscanf(argv[6], "%lf", &x) == 1)    )
 	   then {
 		run_interactive_test(molecule_size,
-				     out_of_range_tolerance_min,
-				     out_of_range_tolerance_max,
+				     boundary_off_centering_tolerance_min,
+				     boundary_off_centering_tolerance_max,
+				     boundary_extrapolation_tolerance_min,
+				     boundary_extrapolation_tolerance_max,
 				     x);
-		return 0;
+		return 0;					/*** NORMAL EXIT ***/
 		}
-	/* fall through */
+	/* error ==> fall through */
 
 default:
-	fprintf(stderr,
-		"usage:\n"
-		"# run a single test as specified:\n"
-		"  %s molecule_size			\\\n"
-		"  %*s out_of_range_tolerance_min	\\\n"
-		"  %*s out_of_range_tolerance_max	\\\n"
-		"  %*s x\n"
-		"# run a preset set of tests:\n"
-		"  %s\n"
-		,
-		argv[0],
-		(int) strlen(argv[0]), "",
-		(int) strlen(argv[0]), "",
-		(int) strlen(argv[0]), "",
-		argv[0]);
+	fprintf(stderr, help_msg);
 	return 1;
 	}
 }
@@ -271,33 +265,39 @@ default:
  */
 static
   void run_interactive_test(int molecule_size,
-			    fp out_of_range_tolerance_min,
-			    fp out_of_range_tolerance_max,
+			    fp boundary_off_centering_tolerance_min,
+			    fp boundary_off_centering_tolerance_max,
+			    fp boundary_extrapolation_tolerance_min,
+			    fp boundary_extrapolation_tolerance_max,
 			    fp x)
 {
+int i_center;
 fp x_rel;
-int m_min, m_max;
 
 printf("testing with molecule_size=%d\n", molecule_size);
-printf("             out_of_range_tolerance_[min,max]=[%g,%g]\n",
-       (double) out_of_range_tolerance_min,
-       (double) out_of_range_tolerance_max);
+printf("             boundary_off_centering_tolerance_[min,max]=[%g,%g]\n",
+       (double) boundary_off_centering_tolerance_min,
+       (double) boundary_off_centering_tolerance_max);
+printf("             boundary_extrapolation_tolerance_[min,max]=[%g,%g]\n",
+       (double) boundary_extrapolation_tolerance_min,
+       (double) boundary_extrapolation_tolerance_max);
 printf("             x=%g\n", (double) x);
 
   {
-const int i_center = LocalInterp_molecule_posn(grid_x0, grid_dx,
-					       grid_i_min, grid_i_max,
-					       molecule_size,
-					       out_of_range_tolerance_min,
-					       out_of_range_tolerance_max,
-					       x,
-					       &x_rel,
-					       &m_min, &m_max);
-
-if ((i_center == INT_MIN) || (i_center == INT_MAX))
-   then printf("i_center=%d\n", i_center);
-   else printf("i_center=%d x_rel=%g m_[min,max]=[%d,%d] i_[min,max]=[%d,%d]\n",
-	       i_center, x_rel, m_min, m_max, i_center+m_min, i_center+m_max);
+const int status
+   = LocalInterp_molecule_posn(grid_x0, grid_dx,
+			       grid_i_min, grid_i_max,
+			       molecule_size,
+			       boundary_off_centering_tolerance_min,
+			       boundary_off_centering_tolerance_max,
+			       boundary_extrapolation_tolerance_min,
+			       boundary_extrapolation_tolerance_max,
+			       x,
+			       &i_center, &x_rel);
+
+if (status < 0)
+   then printf("status=%d\n", status);
+   else printf("i_center=%d x_rel=%g\n", i_center, x_rel);
   }
 }
 
@@ -319,41 +319,48 @@ int failure_count = 0;
 	for (i = 0 ; i < N_TESTS ; ++i)
 	{
 	const struct args_results *p = & test_data[i];
-	fp x_rel = 0.0;
-	int m_min = 0, m_max = 0;
-	int i_center = LocalInterp_molecule_posn(grid_x0, grid_dx,
-						 grid_i_min, grid_i_max,
-						 p->molecule_size,
-						 p->out_of_range_tolerance_min,
-						 p->out_of_range_tolerance_max,
-						 p->x,
-						 &x_rel,
-						 &m_min, &m_max);
-	bool ok = ((i_center == INT_MIN) || (i_center == INT_MAX))
-		  ? (i_center == p->i_center)
-		  : (    (i_center == p->i_center)
-		      && fuzzy_EQ(x_rel, p->x_rel)
-		      && (m_min == p->m_min)
-		      && (m_max == p->m_max)    );
-
-	int msglen =
-	  printf("size=%d tol=[%g,%g] x=%g ==> ",
-		 p->molecule_size,
-		 (double) p->out_of_range_tolerance_min,
-		 (double) p->out_of_range_tolerance_max,
-		 (double) p->x);
-	printf("i_center=%d x_rel=%g m=[%d,%d]\n",
-	       i_center, (double) x_rel, m_min, m_max);
-
-	if (! ok)
-	   then {
+	int i_center = 999;
+	fp  x_rel = 999.999;
+	const int status = LocalInterp_molecule_posn
+				(grid_x0, grid_dx,
+				 grid_i_min, grid_i_max,
+				 p->molecule_size,
+				 p->boundary_off_centering_tolerance_min,
+				 p->boundary_off_centering_tolerance_max,
+				 p->boundary_extrapolation_tolerance_min,
+				 p->boundary_extrapolation_tolerance_max,
+				 p->x,
+				 &i_center, &x_rel);
+	bool ok = (status == 0)
+		  ? ( (status == p->status) && (i_center == p->i_center)
+					    && fuzzy_EQ(x_rel, p->x_rel) )
+		  :   (status == p->status);
+
+	printf("size=%d off_centering_tol=[%g,%g] extrapolation_tol=[%g,%g] x=%g",
+	       p->molecule_size,
+	       (double) p->boundary_off_centering_tolerance_min,
+	       (double) p->boundary_off_centering_tolerance_max,
+	       (double) p->boundary_extrapolation_tolerance_min,
+	       (double) p->boundary_extrapolation_tolerance_max,
+	       (double) p->x);
+	printf(": status=%d", status);
+	if (status == 0)
+	   then printf(" i_center=%d x_rel=%g: ",
+		       i_center, (double) x_rel);
+
+	if (ok)
+	   then printf(": ok\n");
+	   else {
 		++failure_count;
-		  {
-		int error_msglen = printf("***** FAIL: ");
-		printf("%-*s", msglen-error_msglen, "expected");
-		printf("i_center=%d x_rel=%g m_[min,max]=[%d,%d]\n",
-		       p->i_center, (double) p->x_rel, p->m_min, p->m_max);
-		  }
+		printf(": FAIL\n");
+		if (p->status == 0)
+		   then printf("*** expected status=%d i_center=%d x_rel=%g\n",
+			       p->status, p->i_center, (double) p->x_rel);
+		   else printf("*** expected status=%d\n", p->status);
+		if (status == 0)
+		   then printf("*** got      status=%d i_center=%d x_rel=%g\n",
+			       status, i_center, (double) x_rel);
+		   else printf("*** got      status=%d\n", status);
 		}
 	}