finish splitting up (huge) LocalInterp_InterpLocalUniform()

function into smaller subfunctions, now it's "only" about 600 lines long
(before it was around 1200 lines :( :( :( )

also fix a few more bugs in error reporting
and add a few more comments


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

2 files changed, 1184 insertions, 844 deletions

diff --git a/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c b/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c
index 2c1582c..d0b4c78 100644
--- a/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c
+++ b/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c
@@ -1,3 +1,26 @@
+/* InterpLocalUniform.c -- driver for this interpolator */
+/* $Header: */
+/*
+** *****data structures and functions local to this file *****
+**
+ * LocalInterp_ILU_Lagrange - external entry point for Lagrange interpolator
+ * LocalInterp_ILU_Hermite - external entry point for Hermite interpolator
+**
+** LocalInterp_InterpLocalUniform - main driver routine
+**
+** check_boundary_tolerances - check boundary tolerances for consistency
+** get_and_decode_molecule_family - get & decode molecule_family from par table
+** get_molecule_positions - get molecule_positions from parameter table
+** get_Jacobian_query - get Jacobian-query info from parameter table
+** set_error_info - set error information in parameter table
+** set_molecule_structure - set molecule structure info in parameter table
+** set_molecule_min_max_m - set molecule size info in parameter table
+**
+** get_and_check_INT - get and range-check CCTK_INT from parameter table
+** get_INT_array - get CCTK_INT array from parameter table
+** get_REAL_array - get CCTK_REAL array from parameter table
+ */
+
 /*@@
   @file      InterpLocalUniform.c
   @date      22 Oct 2001
@@ -24,7 +47,6 @@
   @version   $Id$
   @@*/
 
-#include <float.h>
 #include <stdio.h>
 #include <stdlib.h>			/* malloc(), free() */
 #include <string.h>
@@ -47,6 +69,8 @@ CCTK_FILEVERSION(CactusBase_LocalInterp_src_GeneralizedPolynomialUniform_InterpL
  * *****data structures and functions local to this file *****
  */
 
+/**************************************/
+
 /*
  * data structures local to this file
  */
@@ -59,6 +83,7 @@ enum	interp_operator
 	N_INTERP_OPERATORS	/* this must be the last entry in the enum */
 	};
 
+/**************************************/
 
 /*
  * prototypes for functions local to this file
@@ -66,25 +91,233 @@ enum	interp_operator
 static
   int LocalInterp_InterpLocalUniform(enum interp_operator interp_operator,
 				     const char* interp_operator_string,
+				     CCTK_REAL boundary_off_cntr_tol_default,
+				     CCTK_REAL boundary_extrap_tol_default,
+				   /***** misc arguments *****/
 				     int N_dims,
 				     int param_table_handle,
-				     /***** coordinate system *****/
+				   /***** coordinate system *****/
 				     const CCTK_REAL coord_origin[],
 				     const CCTK_REAL coord_delta[],
-				     /***** interpolation points *****/
+				   /***** interpolation points *****/
 				     int N_interp_points,
 				     int interp_coords_type_code,
 				     const void *const interp_coords[],
-				     /***** input arrays *****/
+				   /***** input arrays *****/
 				     int N_input_arrays,
 				     const CCTK_INT input_array_dims[],
 				     const CCTK_INT input_array_type_codes[],
 				     const void *const input_arrays[],
-				     /***** output arrays *****/
+				   /***** output arrays *****/
 				     int N_output_arrays,
 				     const CCTK_INT output_array_type_codes[],
 				     void *const output_arrays[]);
 
+static
+  void check_boundary_tolerances
+	(int N_boundaries,
+	 const CCTK_REAL boundary_off_centering_tolerance[MAX_N_BOUNDARIES],
+	 const CCTK_REAL boundary_extrapolation_tolerance[MAX_N_BOUNDARIES]);
+static
+  int get_and_decode_molecule_family
+	(int param_table_handle,
+	 int buffer_size, char molecule_family_string_buffer[],
+	 enum molecule_family *p_molecule_family);
+static
+  int get_molecule_positions
+	(int param_table_handle,
+	 int N_dims, CCTK_INT* molecule_positions_array[MAX_N_DIMS]);
+static
+  int get_Jacobian_info(int param_table_handle,
+			int N_dims, int N_output_arrays,
+			struct Jacobian_info* p_Jacobian_info);
+static
+  int set_error_info(int param_table_handle,
+		     struct error_flags* p_error_flags);
+static
+  int set_molecule_structure
+	(int param_table_handle,
+	 const struct molecule_structure_flags* p_molecule_structure_flags);
+static
+  int set_molecule_min_max_m
+	(int param_table_handle,
+	 int N_dims,
+	 const struct molecule_min_max_m_info* p_molecule_min_max_m_info);
+
+static
+  int get_and_check_INT
+	(int param_table_handle, const char name[],
+	 bool mandatory_flag, int default_value,
+	 bool check_range_flag, int min_value, int max_value,
+				const char max_value_string[],
+	 int* p_value);
+static
+  int get_INT_array(int param_table_handle, const char name[],
+		    bool default_flag, int default_value,
+		    int N, CCTK_INT buffer[],
+		    bool* p_value_not_set);
+static
+  int get_REAL_array(int param_table_handle, const char name[],
+		      CCTK_REAL default_value,
+		      int N, CCTK_REAL buffer[]);
+
+/**************************************/
+
+/*
+ * table of pointers to actual interpolation functions
+ */
+
+/*
+ * typedef  p_interp_fn_t  as a function pointer pointing to an
+ * individual interpolator function of the sort defined by template.[ch]
+ */
+#undef  FUNCTION_NAME
+#define FUNCTION_NAME	(*p_interp_fn_t)
+typedef
+  #include "template.h"
+
+/* NULL (function) pointer of this type */
+#define NULL_INTERP_FN_PTR        ((p_interp_fn_t) NULL)
+
+/*
+ * For each axis of this array which is marked with a "see above" comment
+ * in the array declaration, the array size is declared as X+1, so the
+ * legal subscripts are of course 0...X.  But for these axes we actually
+ * only use 1...X.  (Having the array oversize like this slightly simplifies
+ * the array subscripting.)
+ */
+static const p_interp_fn_t p_interp_fn_table[N_INTERP_OPERATORS]
+					    [MAX_N_DIMS+1]	/* see above */
+					    [N_MOLECULE_FAMILIES]
+					    [MAX_ORDER+1]	/* see above */
+					    [MAX_SMOOTHING+1]
+    = {
+
+      /* in the comments on the initializers, "n.i." = "not implemented" */
+	{
+	/* interp_operator == interp_operator_Lagrange */
+	  {
+	  /* N_dims = 0 ==> unused */
+	    {
+	    /* molecule_family = molecule_family_cube */
+	      { NULL_INTERP_FN_PTR,                },	/* order=0 (unused) */
+	      { NULL_INTERP_FN_PTR,                },	/* order=1 (unused) */
+	      { NULL_INTERP_FN_PTR,                },	/* order=2 (unused) */
+	      { NULL_INTERP_FN_PTR,                },	/* order=3 (unused) */
+	      { NULL_INTERP_FN_PTR,                },	/* order=4 (unused) */
+	      { NULL_INTERP_FN_PTR,                },	/* order=5 (unused) */
+	      { NULL_INTERP_FN_PTR,                },	/* order=6 (unused) */
+	    },
+	  },
+
+	  {
+	  /* N_dims = 1 */
+	    {
+	    /* molecule_family = molecule_family_cube */
+	      { NULL_INTERP_FN_PTR,                },	/* order=0 (unused) */
+	      { LocalInterp_ULagrange_1d_cube10,   },	/* order=1 */
+	      { LocalInterp_ULagrange_1d_cube20,   },	/* order=2 */
+	      { LocalInterp_ULagrange_1d_cube30,   },	/* order=3 */
+	      { LocalInterp_ULagrange_1d_cube40,   },	/* order=4 */
+	      { LocalInterp_ULagrange_1d_cube50,   },	/* order=5 */
+	      { LocalInterp_ULagrange_1d_cube60,   },	/* order=6 */
+	    },
+	  },
+
+	  {
+	  /* N_dims = 2 */
+	    {
+	    /* molecule_family = molecule_family_cube */
+	      { NULL_INTERP_FN_PTR,                },	/* order=0 (unused) */
+	      { LocalInterp_ULagrange_2d_cube10,   },	/* order=1 */
+	      { LocalInterp_ULagrange_2d_cube20,   },	/* order=2 */
+	      { LocalInterp_ULagrange_2d_cube30,   },	/* order=3 */
+	      { LocalInterp_ULagrange_2d_cube40,   },	/* order=4 */
+	      { NULL_INTERP_FN_PTR,                },	/* order=5 (n.i.) */
+	      { NULL_INTERP_FN_PTR,                },	/* order=6 (n.i.) */
+	    },
+	  },
+
+	  {
+	  /* N_dims = 3 */
+	    {
+	    /* molecule_family = molecule_family_cube */
+	      { NULL_INTERP_FN_PTR,                },	/* order=0 (unused) */
+	      { LocalInterp_ULagrange_3d_cube10,   },	/* order=1 */
+	      { LocalInterp_ULagrange_3d_cube20,   },	/* order=2 */
+	      { LocalInterp_ULagrange_3d_cube30,   },	/* order=3 */
+	      { LocalInterp_ULagrange_3d_cube40,   },	/* order=4 */
+	      { NULL_INTERP_FN_PTR,                },	/* order=5 (n.i.) */
+	      { NULL_INTERP_FN_PTR,                },	/* order=6 (n.i.) */
+	    },
+	  },
+
+	/* end of interp_operator == interp_operator_Lagrange */
+	},
+
+	{
+	/* interp_operator == interp_operator_Hermite */
+	  {
+	  /* N_dims = 0 ==> unused */
+	    {
+	    /* molecule_family = molecule_family_cube */
+	      { NULL_INTERP_FN_PTR,                },	/* order=0 (unused) */
+	      { NULL_INTERP_FN_PTR,                },	/* order=1 (unused) */
+	      { NULL_INTERP_FN_PTR,                },	/* order=2 (unused) */
+	      { NULL_INTERP_FN_PTR,                },	/* order=3 (unused) */
+	      { NULL_INTERP_FN_PTR,                },	/* order=4 (unused) */
+	      { NULL_INTERP_FN_PTR,                },	/* order=5 (unused) */
+	      { NULL_INTERP_FN_PTR,                },	/* order=6 (unused) */
+	    },
+	  },
+
+	  {
+	  /* N_dims = 1 */
+	    {
+	    /* molecule_family = molecule_family_cube */
+	      { NULL_INTERP_FN_PTR,                },	/* order=0 (unused) */
+	      { NULL_INTERP_FN_PTR,                },	/* order=1 (unused) */
+	      { LocalInterp_UHermite_1d_cube2,     },	/* order=2 */
+	      { LocalInterp_UHermite_1d_cube3,     },	/* order=3 */
+	      { LocalInterp_UHermite_1d_cube4,     },	/* order=4 */
+	      { NULL_INTERP_FN_PTR,                },	/* order=5 (n.i.) */
+	      { NULL_INTERP_FN_PTR,                },	/* order=6 (n.i.) */
+	    },
+	  },
+
+	  {
+	  /* N_dims = 2 */
+	    {
+	    /* molecule_family = molecule_family_cube */
+	      { NULL_INTERP_FN_PTR,                },	/* order=0 (unused) */
+	      { NULL_INTERP_FN_PTR,                },	/* order=1 (unused) */
+	      { LocalInterp_UHermite_2d_cube2,     },	/* order=2 */
+	      { LocalInterp_UHermite_2d_cube3,     },	/* order=3 */
+	      { NULL_INTERP_FN_PTR,                },	/* order=4 (n.i.) */
+	      { NULL_INTERP_FN_PTR,                },	/* order=5 (n.i.) */
+	      { NULL_INTERP_FN_PTR,                },	/* order=6 (n.i.) */
+	    },
+	  },
+
+	  {
+	  /* N_dims = 3 */
+	    {
+	    /* molecule_family = molecule_family_cube */
+	      { NULL_INTERP_FN_PTR,                },	/* order=0 (unused) */
+	      { NULL_INTERP_FN_PTR,                },	/* order=1 (unused) */
+	      { LocalInterp_UHermite_3d_cube2,     },	/* order=2 */
+	      { LocalInterp_UHermite_3d_cube3,     },	/* order=3 */
+	      { NULL_INTERP_FN_PTR,                },	/* order=4 (n.i.) */
+	      { NULL_INTERP_FN_PTR,                },	/* order=5 (n.i.) */
+	      { NULL_INTERP_FN_PTR,                },	/* order=6 (n.i.) */
+	    },
+	  },
+
+	/* end of interp_operator == interp_operator_Hermite */
+	}
+
+      };
+
 /******************************************************************************/
 /******************************************************************************/
 /******************************************************************************/
@@ -131,6 +364,8 @@ int LocalInterp_ILU_Lagrange(int N_dims,
 {
 return LocalInterp_InterpLocalUniform(interp_operator_Lagrange,
 				      "Lagrange polynomial interpolation",
+				      LAGRANGE_BNDRY_OFF_CNTR_TOL_DEF,
+				      LAGRANGE_BNDRY_EXTRAP_TOL_DEF,
 				      N_dims,
 				      param_table_handle,
 				      coord_origin, coord_delta,
@@ -190,6 +425,8 @@ int LocalInterp_ILU_Hermite(int N_dims,
 {
 return LocalInterp_InterpLocalUniform(interp_operator_Hermite,
 				      "Hermite polynomial interpolation",
+				      HERMITE_BNDRY_OFF_CNTR_TOL_DEF,
+				      HERMITE_BNDRY_EXTRAP_TOL_DEF,
 				      N_dims,
 				      param_table_handle,
 				      coord_origin, coord_delta,
@@ -231,10 +468,6 @@ 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
@@ -245,7 +478,13 @@ return LocalInterp_InterpLocalUniform(interp_operator_Hermite,
             @var boundary_extrapolation_tolerance.
   @endhdesc
 
-  ***** misc arguments *****
+  @hdate    1.Feb.2003
+  @hauthor  Jonathan Thornburg <jthorn@aei.mpg.de>
+  @hdesc    Change defaults to "no off-centering" for Hermite interpolator,
+            split off most of this (huge) function into subfunctions
+  @endhdesc
+
+  ***** operator arguments *****
 
   @var		interp_operator
   @vdesc	describes the interpolation operator
@@ -258,6 +497,22 @@ return LocalInterp_InterpLocalUniform(interp_operator_Hermite,
   @vtype	const char* interp_operator_string
   @endvar
 
+  @var		boundary_off_cntr_tol_default
+  @vdesc	the default value for each element of
+		@var boundary_off_centering_default[]
+		if this key isn't present in the parameter table
+  @vtype	CCTK_REAL boundary_off_cntr_tol_default
+  @endvar
+
+  @var		boundary_extrap_tol_default
+  @vdesc	the default value for each element of
+		@var boundary_extrapolation_default[]
+		if this key isn't present in the parameter table
+  @vtype	CCTK_REAL boundary_extrap_tol_default
+  @endvar
+
+  ***** misc arguments *****
+
   @var		N_dims
   @vdesc	dimensionality of the interpolation
   @vtype	int N_dims			(must be >= 1)
@@ -358,7 +613,7 @@ return LocalInterp_InterpLocalUniform(interp_operator_Hermite,
   @vdesc	Sets the order of the interpolating polynomial
 		(1=linear, 2=quadratic, 3=cubic, ...).
 		This table entry is mandatory; all others are optional.
-  @vtype	CCTK_INT order
+  @vtype	int order
   @endvar
 
   @var		N_boundary_points_to_omit
@@ -611,7 +866,7 @@ return LocalInterp_InterpLocalUniform(interp_operator_Hermite,
 		UTIL_ERROR_BAD_INPUT	one of the input arguments is invalid
 		UTIL_ERROR_NO_MEMORY	unable to malloc() temporary memory
 		UTIL_ERROR_BAD_HANDLE	invalid parameter table handle
-		CCTK_ERROR_INTERP_POINT_X_RANGE
+		CCTK_ERROR_INTERP_POINT_OUTSIDE
 					interpolation point is out of range
 					(in this case additional information
 					is reported through the parameter table)
@@ -622,21 +877,24 @@ return LocalInterp_InterpLocalUniform(interp_operator_Hermite,
 static
   int LocalInterp_InterpLocalUniform(enum interp_operator interp_operator,
 				     const char* interp_operator_string,
+				     CCTK_REAL boundary_off_cntr_tol_default,
+				     CCTK_REAL boundary_extrap_tol_default,
+				   /***** misc arguments *****/
 				     int N_dims,
 				     int param_table_handle,
-				     /***** coordinate system *****/
+				   /***** coordinate system *****/
 				     const CCTK_REAL coord_origin[],
 				     const CCTK_REAL coord_delta[],
-				     /***** interpolation points *****/
+				   /***** interpolation points *****/
 				     int N_interp_points,
 				     int interp_coords_type_code,
 				     const void *const interp_coords[],
-				     /***** input arrays *****/
+				   /***** input arrays *****/
 				     int N_input_arrays,
 				     const CCTK_INT input_array_dims[],
 				     const CCTK_INT input_array_type_codes[],
 				     const void *const input_arrays[],
-				     /***** output arrays *****/
+				   /***** output arrays *****/
 				     int N_output_arrays,
 				     const CCTK_INT output_array_type_codes[],
 				     void *const output_arrays[])
@@ -662,7 +920,8 @@ static
 int N_input_arrays1  = N_input_arrays  + 1;
 int N_output_arrays1 = N_output_arrays + 1;
 
-int status, status1, status2, status3, status4;
+int status, status1, status2;
+bool value_not_set;
 
 /******************************************************************************/
 
@@ -707,285 +966,98 @@ if (N_dims > MAX_N_DIMS)
 /******************************************************************************/
 
 /*
- * interpolation order
+ * get the parameters from the parameter table, filling in defaults as needed
  */
+
   {
 CCTK_INT order;
-status = Util_TableGetInt(param_table_handle, &order, "order");
-if	(status == 1)
-   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 for\n"
-"                              \"order\" parameter!\n"
-"                              (this is a mandatory parameter)\n"
-"                              error status=%d",
-		   status);
-	return status;					/*** ERROR RETURN ***/
-	}
-if ((order < 1) || (order > MAX_ORDER))
-   then {
-	CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
-		   __LINE__, __FILE__, CCTK_THORNSTRING,
-"\n"
-"   CCTK_InterpLocalUniform(): implementation restriction: order=%d\n"
-"                              but we only support 1 <= order <= MAX_ORDER=%d!",
-		   order,
-		   MAX_ORDER);
-	return UTIL_ERROR_BAD_INPUT;			/*** ERROR RETURN ***/
-	}
-   
-/******************************************************************************/
+status = get_and_check_INT(param_table_handle, "order",
+			   true, 0,	/* this is a mandatory parameter */
+			   true, 1, MAX_ORDER, "MAX_ORDER",   /* range check */
+			   &order);
+if (status != 0)
+   then return status;					/*** ERROR RETURN ***/
 
-/*
- * out-of-range interpolation-point handling
- */
   {
-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 {
-	/* set the default value */
-	int ibndry;
-		for (ibndry = 0 ; ibndry < N_boundaries ; ++ibndry)
-		{
-		N_boundary_points_to_omit[ibndry] = 0;
-		}
-	}
-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 {
-	/* set the default value */
-	int ibndry;
-		for (ibndry = 0 ; ibndry < N_boundaries ; ++ibndry)
-		{
-		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(): bad table entry for\n"
-"                              \"boundary_off_centering_tolerance\" parameter!\n"
-"                              error status=%d",
-		   status);
-	return status;					/*** ERROR RETURN ***/
-	}
+CCTK_INT  N_boundary_points_to_omit[MAX_N_BOUNDARIES];
+status = get_INT_array(param_table_handle, "N_boundary_points_to_omit",
+		       true, 0,		/* default value */
+		       N_boundaries, N_boundary_points_to_omit,
+		       NULL);
+if (status != 0)
+   then 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 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():\n"
-"        warning: setting boundary_off_centering_tolerance[%d] == 0.0\n"
-"                 and     boundary_extrapolation_tolerance[%d] == %g > 0.0\n"
-"                 is almost certainly a mistake\n"
-"                 (the boundary_extrapolation_tolerance[%d] == %g\n"
-"                  setting will have no effect because of the\n"
-"                  boundary_off_centering_tolerance[%d] == 0.0 setting)"
-			   ,
-			   ibndry,
-			   ibndry,
-			      (double) boundary_extrapolation_tolerance[ibndry],
-			   ibndry,
-			      (double) boundary_extrapolation_tolerance[ibndry],
-			   ibndry);
-	}
-  }
+CCTK_REAL boundary_off_centering_tolerance[MAX_N_BOUNDARIES];
+CCTK_REAL boundary_extrapolation_tolerance[MAX_N_BOUNDARIES];
+status = get_REAL_array(param_table_handle, "boundary_off_centering_tolerance",
+			boundary_off_cntr_tol_default,
+			N_boundaries, boundary_off_centering_tolerance);
+if (status != 0)
+   then return status;					/*** ERROR RETURN ***/
+status = get_REAL_array(param_table_handle, "boundary_extrapolation_tolerance",
+			boundary_extrap_tol_default,
+			N_boundaries, boundary_extrapolation_tolerance);
+if (status != 0)
+   then return status;					/*** ERROR RETURN ***/
+check_boundary_tolerances(N_boundaries,
+			  boundary_off_centering_tolerance,
+			  boundary_extrapolation_tolerance);
 
-/******************************************************************************/
+/**************************************/
 
-/*
- * molecule family
- */
   {
 #define MOLECULE_FAMILY_BUFSIZ	(MAX_MOLECULE_FAMILY_STRLEN+1)
-char molecule_family_strbuf[MOLECULE_FAMILY_BUFSIZ];
-char* molecule_family_str = molecule_family_strbuf;
+char molecule_family_string[MOLECULE_FAMILY_BUFSIZ];
 enum molecule_family molecule_family;
-status = Util_TableGetString(param_table_handle,
-			     MOLECULE_FAMILY_BUFSIZ, molecule_family_strbuf,
-			     "molecule_family");
-if	(status == UTIL_ERROR_TABLE_NO_SUCH_KEY)
-   then {
-	/* set the default value for code here*/
-	/* (need enum for main code, string for error messages) */
-	molecule_family     = molecule_family_cube;
-	molecule_family_str = "cube";
+status = get_and_decode_molecule_family(param_table_handle,
+			     MOLECULE_FAMILY_BUFSIZ, molecule_family_string,
+			     &molecule_family);
+if (status != 0)
+   then return status;					/*** ERROR RETURN ***/
 
-	/* set this key in the parameter table */
-	/* to give the (default) molecule family we're going to use */
-	status = Util_TableSetString(param_table_handle,
-				     molecule_family_str,
-				     "molecule_family");
-	if (status < 0)
-	   then {
-		CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
-			   __LINE__, __FILE__, CCTK_THORNSTRING,
-"\n"
-"   CCTK_InterpLocalUniform(): error setting default molecule family\n"
-"                              in parameter table!\n"
-"                              error status=%d",
-			   status);
-		return status;				/*** ERROR RETURN ***/
-		}
-	}
-else if (status > 0)
-   then {
-	/* decode molecule family string */
-	if (strcmp(molecule_family_strbuf, "cube") == 0)
-	   then molecule_family = molecule_family_cube;
-	else	{
-		CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
-			   __LINE__, __FILE__, CCTK_THORNSTRING,
-"CCTK_InterpLocalUniform(): unknown molecule_family=\"%s\"!",
-			   molecule_family_strbuf);
-		return UTIL_ERROR_BAD_INPUT;		/*** ERROR RETURN ***/
-		}
-	}
-else	{
-	CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
-		   __LINE__, __FILE__, CCTK_THORNSTRING,
-"\n"
-"   CCTK_InterpLocalUniform(): bad table entry for\n"
-"                              \"molecule_family\" parameter!\n"
-"                              error status=%d",
-		   status);
-	return status;					/*** ERROR RETURN ***/
-	}
-
-/**************************************/
-
-/*
- * smoothing
- */
   {
 CCTK_INT smoothing;
-status = Util_TableGetInt(param_table_handle, &smoothing, "smoothing");
-if	(status == UTIL_ERROR_TABLE_NO_SUCH_KEY)
-   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\n"
-"                              \"smoothing\" parameter!\n"
-"                              error status=%d",
-		   status);
-	return status;					/*** ERROR RETURN ***/
-	}
-if ((smoothing < 0) || (smoothing > MAX_SMOOTHING))
-   then return UTIL_ERROR_BAD_INPUT;			/*** ERROR RETURN ***/
+status = get_and_check_INT(param_table_handle, "smoothing",
+			   false, 0,	/* optional, default value 0 */
+			   true, 0, MAX_SMOOTHING, "MAX_SMOOTHING",
+							/* range check */
+			   &smoothing);
+if (status != 0)
+   then return status;					/*** ERROR RETURN ***/
 
-/******************************************************************************/
+/**************************************/
 
-/*
- * input array offsets
- */
   {
 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,
-			       N_input_arrays, input_array_offsets,
-			       "input_array_offsets");
-if	 (status == UTIL_ERROR_TABLE_NO_SUCH_KEY)
+status = get_INT_array(param_table_handle, "input_array_offsets",
+		       true, 0,		/* default value */
+		       N_input_arrays, input_array_offsets,
+		       NULL);
+if (status != 0)
    then {
-	/* default offset = 0 along each axis */
-	LocalInterp_zero_int_array(N_input_arrays, input_array_offsets);
-	}
-else if (status == N_input_arrays)
-   then { }	/* ok ==> no-op here */
-else	{
 	free(input_array_offsets);
-	CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
-		   __LINE__, __FILE__, CCTK_THORNSTRING,
-"\n"
-"   CCTK_InterpLocalUniform(): bad table entry for\n"
-"                              \"input_array_offsets\" parameter!\n"
-"                              error status=%d",
-		   status);
 	return status;					/*** ERROR RETURN ***/
 	}
 
-/**************************************/
-
-/*
- * input array strides
- */
   {
 CCTK_INT input_array_strides[MAX_N_DIMS];
-status = Util_TableGetIntArray(param_table_handle,
-			       N_dims, input_array_strides,
-			       "input_array_strides");
-if	(status == UTIL_ERROR_TABLE_NO_SUCH_KEY)
+status = get_INT_array(param_table_handle, "input_array_strides",
+		       false, 0,	/* don't set default value */
+		       N_dims, input_array_strides,
+		       &value_not_set);
+if (status != 0)
+   then {
+	free(input_array_offsets);
+	return status;					/*** ERROR RETURN ***/
+	}
+if (value_not_set)
    then {
 	/* default stride = Fortran storage order, */
 	/* determined from input_array_dims[] */
@@ -1009,54 +1081,30 @@ if	(status == UTIL_ERROR_TABLE_NO_SUCH_KEY)
 		}
 	  }
 	}
-else if (status == N_dims)
-   then { }	/* ok ==> no-op here */
-else	{
-	free(input_array_offsets);
-	CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
-		   __LINE__, __FILE__, CCTK_THORNSTRING,
-"\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
- */
   {
 CCTK_INT input_array_min_subscripts[MAX_N_DIMS];
-status = Util_TableGetIntArray(param_table_handle,
-			       N_dims, input_array_min_subscripts,
-			       "input_array_min_subscripts");
-if	(status == UTIL_ERROR_TABLE_NO_SUCH_KEY)
+CCTK_INT input_array_max_subscripts[MAX_N_DIMS];
+status = get_INT_array(param_table_handle, "input_array_min_subscripts",
+		       true, 0,		/* default value */
+		       N_dims, input_array_min_subscripts,
+		       NULL);
+if (status != 0)
    then {
-	/* default min subscript = 0 along each axis */
-	LocalInterp_zero_int_array(N_dims, input_array_min_subscripts);
+	free(input_array_offsets);
+	return status;					/*** ERROR RETURN ***/
 	}
-else if (status == N_dims)
-   then { }	/* ok ==> no-op here */
-else	{
+
+status = get_INT_array(param_table_handle, "input_array_max_subscripts",
+		       false, 0,	/* don't set default value */
+		       N_dims, input_array_max_subscripts,
+		       &value_not_set);
+if (status != 0)
+   then {
 	free(input_array_offsets);
-	CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
-		   __LINE__, __FILE__, CCTK_THORNSTRING,
-"\n"
-"   CCTK_InterpLocalUniform(): bad table entry for\n"
-"                              \"input_array_min_subscripts\" parameter!\n"
-"                              error status=%d",
-		   status);
 	return status;					/*** ERROR RETURN ***/
 	}
-  {
-CCTK_INT input_array_max_subscripts[MAX_N_DIMS];
-status = Util_TableGetIntArray(param_table_handle,
-			       N_dims, input_array_max_subscripts,
-			       "input_array_max_subscripts");
-if	(status == UTIL_ERROR_TABLE_NO_SUCH_KEY)
+if (value_not_set)
    then {
 	/* default max subscript = input_array_dims[]-1 along each axis */
 	if (input_array_dims == NULL)
@@ -1077,25 +1125,9 @@ if	(status == UTIL_ERROR_TABLE_NO_SUCH_KEY)
 		}
 	  }
 	}
-else if (status == N_dims)
-   then { }	/* ok ==> no-op here */
-else	{
-	free(input_array_offsets);
-	CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
-		   __LINE__, __FILE__, CCTK_THORNSTRING,
-"\n"
-"   CCTK_InterpLocalUniform(): bad table entry for\n"
-"                              \"input_array_max_subscripts\" parameter!\n"
-"                              error status=%d",
-		   status);
-	return status;					/*** ERROR RETURN ***/
-	}
 
-/******************************************************************************/
+/**************************************/
 
-/*
- * operand indices
- */
   {
 CCTK_INT* const operand_indices
 	= (CCTK_INT*) malloc(N_output_arrays1 * sizeof(CCTK_INT));
@@ -1104,23 +1136,30 @@ if (operand_indices == NULL)
 	free(input_array_offsets);
 	return UTIL_ERROR_NO_MEMORY;			/*** ERROR RETURN ***/
 	}
-status = Util_TableGetIntArray(param_table_handle,
-			       N_output_arrays, operand_indices,
-			       "operand_indices");
-if	 (status == UTIL_ERROR_TABLE_NO_SUCH_KEY)
+status = get_INT_array(param_table_handle, "operand_indices",
+		       false, 0,	/* don't set default value */
+		       N_output_arrays, operand_indices,
+		       &value_not_set);
+if (status != 0)
+   then {
+	free(input_array_offsets);
+	free(operand_indices);
+	return UTIL_ERROR_NO_MEMORY;			/*** ERROR RETURN ***/
+	}
+if (value_not_set)
    then {
 	/* default operand will use each input exactly once, */
 	/* but this only makes since if N_input_arrays == N_output_arrays */
 	if (N_input_arrays != N_input_arrays)
 	   then {
-		free(operand_indices);
 		free(input_array_offsets);
+		free(operand_indices);
 		CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
 			   __LINE__, __FILE__, CCTK_THORNSTRING,
 "\n"
 "   CCTK_InterpLocalUniform(): can't default \"operand_indices\"\n"
 "                              with N_input_arrays=%d != N_output_arrays=%d!"
-,
+			   ,
 			   N_input_arrays, N_output_arrays);
 		return UTIL_ERROR_BAD_INPUT;		/*** ERROR RETURN ***/
 		}
@@ -1132,102 +1171,40 @@ if	 (status == UTIL_ERROR_TABLE_NO_SUCH_KEY)
 		}
 	  }
 	}
-else if (status == N_output_arrays)
-   then {
-	/* check that all operands are within range */
-	int out;
-		for (out = 0 ; out < N_output_arrays ; ++out)
-		{
-		if (    (operand_indices[out] < 0)
-		     || (operand_indices[out] >= N_input_arrays)    )
-		   then {
-			free(operand_indices);
-			free(input_array_offsets);
-			CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
-				   __LINE__, __FILE__, CCTK_THORNSTRING,
-"\n"
-"   CCTK_InterpLocalUniform(): operand_indices[out=%d]=%d specifies\n"
-"                              nonexistent input array!\n"
-"                              (valid range is 0 <= value < N_input_arrays=%d)"
-,
-				   out, operand_indices[out],
-				   N_input_arrays);
-			return UTIL_ERROR_BAD_INPUT;	/*** ERROR RETURN ***/
-			}
-		}
-	}
-else	{
-	free(operand_indices);
-	free(input_array_offsets);
-	CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
-		   __LINE__, __FILE__, CCTK_THORNSTRING,
-"\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
 	= (CCTK_INT*) malloc(N_output_arrays1 * sizeof(CCTK_INT));
 if (operation_codes == NULL)
    then {
-	free(operand_indices);
 	free(input_array_offsets);
+	free(operand_indices);
 	return UTIL_ERROR_NO_MEMORY;			/*** ERROR RETURN ***/
 	}
-status = Util_TableGetIntArray(param_table_handle,
-			       N_output_arrays, operation_codes,
-			       "operation_codes");
-if	 (status == UTIL_ERROR_TABLE_NO_SUCH_KEY)
+status = get_INT_array(param_table_handle, "operation_codes",
+		       true, 0,		/* default value */
+		       N_output_arrays, operation_codes,
+		       NULL);
+if (status != 0)
    then {
-	/* default operation_codes = all 0 */
-	/* but this only makes since if N_input_arrays == N_output_arrays */
-	if (N_input_arrays != N_input_arrays)
-	   then {
-		free(operation_codes);
-		free(operand_indices);
-		free(input_array_offsets);
-		CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
-			   __LINE__, __FILE__, CCTK_THORNSTRING,
-"\n"
-"   CCTK_InterpLocalUniform(): can't default \"operation_codes\"\n"
-"                              with N_input_arrays=%d != N_output_arrays=%d!"
-,
-			   N_input_arrays, N_output_arrays);
-		return UTIL_ERROR_BAD_INPUT;		/*** ERROR RETURN ***/
-		}
-	LocalInterp_zero_int_array(N_output_arrays, operation_codes);
-	}
-else if (status == N_output_arrays)
-   then { }	/* ok ==> no-op here */
-else	{
-	free(operation_codes);
-	free(operand_indices);
 	free(input_array_offsets);
-	CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
-		   __LINE__, __FILE__, CCTK_THORNSTRING,
-"\n"
-"   CCTK_InterpLocalUniform(): bad table entry for\n"
-"                              \"operation_codes\" parameter!\n"
-"                         error status=%d",
-		   status);
-	return status;					/*** ERROR RETURN ***/
+	free(operand_indices);
+	free(operation_codes);
+	return UTIL_ERROR_NO_MEMORY;			/*** ERROR RETURN ***/
 	}
 
 /******************************************************************************/
 
 /*
- * set up for any molecule min/max m and position queries
+ * set up for any molecule min/max m and/or position queries
  */
 
+  {
+struct molecule_min_max_m_info* p_molecule_min_max_m_info = NULL;
+struct molecule_min_max_m_info  molecule_min_max_m_info;
+
 /* molecule min/max m */
 status1 = Util_TableQueryValueInfo(param_table_handle,
 				   NULL, NULL,
@@ -1237,60 +1214,54 @@ status2 = Util_TableQueryValueInfo(param_table_handle,
 				   "molecule_max_m");
 if ((status1 < 0) || (status2 < 0))
    then {
-	free(operation_codes);
-	free(operand_indices);
 	free(input_array_offsets);
+	free(operand_indices);
+	free(operation_codes);
 	CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
 		   __LINE__, __FILE__, CCTK_THORNSTRING,
 "\n"
-"   CCTK_InterpLocalUniform(): bad \"interp_molecule_{min,max}_m\"\n"
-"                              table entry/entries!\n"
+"   CCTK_InterpLocalUniform(): bad \"molecule_{min,max}_m\"\n"
+"                              table entry/entries in query!\n"
 "                              error status1=%d status2=%d"
-,
+		   ,
 		   status1, status2);
 	return (status1 < 0) ? status1 : status2;	/*** ERROR RETURN ***/
 	}
-  {
-struct molecule_min_max_m_info molecule_min_max_m_info;
-struct molecule_min_max_m_info *p_molecule_min_max_m_info
-	= (status1 && status2) ? &molecule_min_max_m_info : NULL;
+if (status1 && status2)
+   then p_molecule_min_max_m_info = &molecule_min_max_m_info;
 
-/* molecule positions */
   {
-CCTK_INT* const* p_molecule_positions = NULL;
-CCTK_POINTER molecule_positions_temp [MAX_N_DIMS];
-CCTK_INT*    molecule_positions_array[MAX_N_DIMS];
-status = Util_TableGetPointerArray(param_table_handle,
-				   N_dims, molecule_positions_temp,
-				   "molecule_positions");
-if	(status == UTIL_ERROR_TABLE_NO_SUCH_KEY)
-   then {
-	/* no querying of interpolator molecule positions */
-	}
-else if (status == N_dims)
+CCTK_INT** p_molecule_positions = NULL;
+CCTK_INT*  molecule_positions_array[MAX_N_DIMS];
+
+/* are we doing a molecule-positions query? */
+status = Util_TableQueryValueInfo(param_table_handle,
+				  NULL, NULL,
+				  "molecule_positions");
+if (status < 0)
    then {
-	/* type-cast CCTK_POINTER pointers to CCTK_INT* pointers */
-	/* (which point to where the query results will be stored) */
-	int axis;
-		for (axis = 0 ; axis < N_dims ; ++axis)
-		{
-		molecule_positions_array[axis]
-			= (CCTK_INT *) molecule_positions_temp[axis];
-		}
-	p_molecule_positions = molecule_positions_array;
-	}
-else	{
-	free(operation_codes);
-	free(operand_indices);
 	free(input_array_offsets);
+	free(operand_indices);
+	free(operation_codes);
 	CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
 		   __LINE__, __FILE__, CCTK_THORNSTRING,
 "\n"
-"   CCTK_InterpLocalUniform(): bad \"molecule_positions\" table entry!\n"
-"                              error status=%d",
+"   CCTK_InterpLocalUniform(): bad \"molecule_positions\"\n"
+"                              table entry in query!\n"
+"                              error status=%d"
+		   ,
 		   status);
 	return status;					/*** ERROR RETURN ***/
 	}
+if (status)
+   then {
+	/* yes, we're doing a molecule-positions query */
+	status = get_molecule_positions(param_table_handle,
+				       N_dims, molecule_positions_array);
+	if (status != 0)
+	   then return status;				/*** ERROR RETURN ***/
+	p_molecule_positions = molecule_positions_array;
+	}
 
 /******************************************************************************/
 
@@ -1299,177 +1270,44 @@ else	{
  */
 
   {
-struct Jacobian_info *p_Jacobian_info = NULL;
-struct Jacobian_info Jacobian_info;
+struct Jacobian_info* p_Jacobian_info = NULL;
+struct Jacobian_info  Jacobian_info;
 Jacobian_info.Jacobian_pointer = NULL;
 Jacobian_info.Jacobian_offset = NULL;
 
-/* first find out if we are doing Jacobian queries at all */
+/* are we doing a Jacobian query? */
 status = Util_TableQueryValueInfo(param_table_handle,
 				  NULL, NULL,
 				  "Jacobian_pointer");
 if (status < 0)
    then {
-	free(operation_codes);
-	free(operand_indices);
 	free(input_array_offsets);
+	free(operand_indices);
+	free(operation_codes);
 	CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
 		   __LINE__, __FILE__, CCTK_THORNSTRING,
 "\n"
 "   CCTK_InterpLocalUniform(): bad \"Jacobian_pointer\" table entry!\n"
 "                              (preliminary check)\n"
 "                              error status=%d"
-,
+		   ,
 		   status);
 	 return status;					/*** ERROR RETURN ***/
 	}
 
 if (status)
    then {
-	/* yes, we're doing Jacobian queries */
-	Jacobian_info.Jacobian_pointer
-	    = (CCTK_REAL**) malloc(N_output_arrays1 * sizeof(CCTK_REAL *));
-	if (Jacobian_info.Jacobian_pointer == NULL)
-	   then {
-		free(operation_codes);
-		free(operand_indices);
-		free(input_array_offsets);
-		return UTIL_ERROR_NO_MEMORY;		/*** ERROR RETURN ***/
-		}
-	  {
-	CCTK_POINTER* Jacobian_pointer_temp
-	    = (CCTK_POINTER*) malloc(N_output_arrays1 * sizeof(CCTK_POINTER));
-	if (Jacobian_pointer_temp == NULL)
+	/* yes, we're doing a Jacobian query */
+	status = get_Jacobian_info(param_table_handle,
+				   N_dims, N_output_arrays,
+				   &Jacobian_info);
+	if (status != 0)
 	   then {
-		free(operation_codes);
-		free(operand_indices);
-		free(input_array_offsets);
-		free(Jacobian_info.Jacobian_pointer);
-		return UTIL_ERROR_NO_MEMORY;		/*** ERROR RETURN ***/
-		}
-	status = Util_TableGetPointerArray(param_table_handle,
-					   N_output_arrays,
-					   Jacobian_pointer_temp,
-					   "Jacobian_pointer");
-	if	(status == N_output_arrays)
-	   then {
-		/* type-cast CCTK_POINTER pointers to CCTK_REAL* pointers */
-		/* (which point to where the query results will be stored) */
-		int out;
-			for (out = 0 ; out < N_output_arrays ; ++out)
-			{
-			Jacobian_info.Jacobian_pointer[out]
-				= (CCTK_REAL *) Jacobian_pointer_temp[out];
-			}
-		free(Jacobian_pointer_temp);
-		}
-	else	{
-		free(operation_codes);
-		free(operand_indices);
 		free(input_array_offsets);
-		free(Jacobian_info.Jacobian_pointer);
-		free(Jacobian_pointer_temp);
-		CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
-			   __LINE__, __FILE__, CCTK_THORNSTRING,
-"\n"
-"   CCTK_InterpLocalUniform(): bad \"Jacobian_pointer\" table entry!\n"
-"                              error status=%d"
-,
-			   status);
-		return status;				/*** ERROR RETURN ***/
-		}
-	  }
-
-	Jacobian_info.Jacobian_offset
-		= (CCTK_INT*) malloc(N_output_arrays1 * sizeof(CCTK_INT));
-	if (Jacobian_info.Jacobian_offset == NULL)
-	   then {
-		free(operation_codes);
 		free(operand_indices);
-		free(input_array_offsets);
-		free(Jacobian_info.Jacobian_pointer);
-		return UTIL_ERROR_NO_MEMORY;		/*** ERROR RETURN ***/
-		}
-	status = Util_TableGetIntArray(param_table_handle,
-				       N_output_arrays,
-				       Jacobian_info.Jacobian_offset,
-				       "Jacobian_offset");
-	if	(status == N_output_arrays)
-	   then {
-		/* ok ==> no-op here */
-		}
-	else if (status == UTIL_ERROR_TABLE_NO_SUCH_KEY)
-	   then {
-		/* default offset = all 0 */
-		LocalInterp_zero_int_array(N_output_arrays,
-					   Jacobian_info.Jacobian_offset);
-		}
-	else	{
 		free(operation_codes);
-		free(operand_indices);
-		free(input_array_offsets);
-		free(Jacobian_info.Jacobian_pointer);
-		free(Jacobian_info.Jacobian_offset);
-		CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
-			   __LINE__, __FILE__, CCTK_THORNSTRING,
-"\n"
-"   CCTK_InterpLocalUniform(): bad \"Jacobian_offset\" table entry!\n"
-"                              error status=%d"
-,
-			   status);
 		return status;				/*** ERROR RETURN ***/
 		}
-
-	status1 = Util_TableGetInt(param_table_handle,
-				   & Jacobian_info.Jacobian_interp_point_stride,
-				   "Jacobian_interp_point_stride");
-	status2 = Util_TableGetIntArray(param_table_handle,
-					N_dims,
-					Jacobian_info.Jacobian_m_strides,
-					"Jacobian_m_strides");
-	if ((status1 < 0) || (status2 != N_dims))
-	   then {
-		free(operation_codes);
-		free(operand_indices);
-		free(input_array_offsets);
-		free(Jacobian_info.Jacobian_pointer);
-		free(Jacobian_info.Jacobian_offset);
-		CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
-			   __LINE__, __FILE__, CCTK_THORNSTRING,
-"\n"
-"   CCTK_InterpLocalUniform(): bad \"Jacobian_interp_point_stride\" and/or\n"
-"                              \"Jacobian_m_strides\" table entry/entries!\n"
-"                              error status1=%d status2=%d"
-,
-			   status1, status2);
-		return (status1 < 0) ? status1 : status2; /*** ERROR RETURN ***/
-		}
-
-	status = Util_TableGetInt(param_table_handle,
-				  & Jacobian_info.Jacobian_part_stride,
-				  "Jacobian_part_stride");
-	if	(status == 1)
-	   then {
-		/* ok ==> no-op here */
-		}
-	else if (status == UTIL_ERROR_TABLE_NO_SUCH_KEY)
-	   then Jacobian_info.Jacobian_part_stride = 0;	/* default */
-	else	{
-		free(operation_codes);
-		free(operand_indices);
-		free(input_array_offsets);
-		free(Jacobian_info.Jacobian_pointer);
-		free(Jacobian_info.Jacobian_offset);
-		CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
-			   __LINE__, __FILE__, CCTK_THORNSTRING,
-"\n"
-"   CCTK_InterpLocalUniform(): bad \"Jacobian_part_stride\" table entry!\n"
-"                              error status=%d"
-,
-			   status);
-		return status;				/*** ERROR RETURN ***/
-		}
-
 	p_Jacobian_info = & Jacobian_info;
 	}
 
@@ -1479,173 +1317,56 @@ if (status)
  * decode (interp_operator, N_dims, molecule_family, order, smoothing)
  * and call the appropriate subfunction to do the actual interpolation
  */
-  {
-/*
- * typedef  p_interp_fn_t  as a function pointer pointing to an
- * individual interpolator function of the sort defined by template.[ch]
- */
-#undef  FUNCTION_NAME
-#define FUNCTION_NAME	(*p_interp_fn_t)
-typedef
-  #include "template.h"
-
-/* NULL (function) pointer of this type */
-#define NULL_INTERP_FN_PTR        ((p_interp_fn_t) NULL)
-
-/*
- * For each axis of this array which is marked with a "see above" comment
- * in the array declaration, the array size is declared as X+1, so the
- * legal subscripts are of course 0...X.  But for these axes we actually
- * only use 1...X.  (Having the array oversize like this slightly simplifies
- * the array subscripting.)
- */
-static const p_interp_fn_t p_interp_fn_table[N_INTERP_OPERATORS]
-					    [MAX_N_DIMS+1]	/* see above */
-					    [N_MOLECULE_FAMILIES]
-					    [MAX_ORDER+1]	/* see above */
-					    [MAX_SMOOTHING+1]
-    = {
-
-      /* in the comments on the initializers, "n.i." = "not implemented" */
-	{
-	/* interp_operator == interp_operator_Lagrange */
-	  {
-	  /* N_dims = 0 ==> unused */
-	    {
-	    /* molecule_family = molecule_family_cube */
-	      { NULL_INTERP_FN_PTR,                },	/* order=0 (unused) */
-	      { NULL_INTERP_FN_PTR,                },	/* order=1 (unused) */
-	      { NULL_INTERP_FN_PTR,                },	/* order=2 (unused) */
-	      { NULL_INTERP_FN_PTR,                },	/* order=3 (unused) */
-	      { NULL_INTERP_FN_PTR,                },	/* order=4 (unused) */
-	      { NULL_INTERP_FN_PTR,                },	/* order=5 (unused) */
-	      { NULL_INTERP_FN_PTR,                },	/* order=6 (unused) */
-	    },
-	  },
-
-	  {
-	  /* N_dims = 1 */
-	    {
-	    /* molecule_family = molecule_family_cube */
-	      { NULL_INTERP_FN_PTR,                },	/* order=0 (unused) */
-	      { LocalInterp_ULagrange_1d_cube10,   },	/* order=1 */
-	      { LocalInterp_ULagrange_1d_cube20,   },	/* order=2 */
-	      { LocalInterp_ULagrange_1d_cube30,   },	/* order=3 */
-	      { LocalInterp_ULagrange_1d_cube40,   },	/* order=4 */
-	      { LocalInterp_ULagrange_1d_cube50,   },	/* order=5 */
-	      { LocalInterp_ULagrange_1d_cube60,   },	/* order=6 */
-	    },
-	  },
 
-	  {
-	  /* N_dims = 2 */
-	    {
-	    /* molecule_family = molecule_family_cube */
-	      { NULL_INTERP_FN_PTR,                },	/* order=0 (unused) */
-	      { LocalInterp_ULagrange_2d_cube10,   },	/* order=1 */
-	      { LocalInterp_ULagrange_2d_cube20,   },	/* order=2 */
-	      { LocalInterp_ULagrange_2d_cube30,   },	/* order=3 */
-	      { LocalInterp_ULagrange_2d_cube40,   },	/* order=4 */
-	      { NULL_INTERP_FN_PTR,                },	/* order=5 (n.i.) */
-	      { NULL_INTERP_FN_PTR,                },	/* order=6 (n.i.) */
-	    },
-	  },
 
-	  {
-	  /* N_dims = 3 */
-	    {
-	    /* molecule_family = molecule_family_cube */
-	      { NULL_INTERP_FN_PTR,                },	/* order=0 (unused) */
-	      { LocalInterp_ULagrange_3d_cube10,   },	/* order=1 */
-	      { LocalInterp_ULagrange_3d_cube20,   },	/* order=2 */
-	      { LocalInterp_ULagrange_3d_cube30,   },	/* order=3 */
-	      { LocalInterp_ULagrange_3d_cube40,   },	/* order=4 */
-	      { NULL_INTERP_FN_PTR,                },	/* order=5 (n.i.) */
-	      { NULL_INTERP_FN_PTR,                },	/* order=6 (n.i.) */
-	    },
-	  },
-
-	/* end of interp_operator == interp_operator_Lagrange */
-	},
-
-	{
-	/* interp_operator == interp_operator_Hermite */
-	  {
-	  /* N_dims = 0 ==> unused */
-	    {
-	    /* molecule_family = molecule_family_cube */
-	      { NULL_INTERP_FN_PTR,                },	/* order=0 (unused) */
-	      { NULL_INTERP_FN_PTR,                },	/* order=1 (unused) */
-	      { NULL_INTERP_FN_PTR,                },	/* order=2 (unused) */
-	      { NULL_INTERP_FN_PTR,                },	/* order=3 (unused) */
-	      { NULL_INTERP_FN_PTR,                },	/* order=4 (unused) */
-	      { NULL_INTERP_FN_PTR,                },	/* order=5 (unused) */
-	      { NULL_INTERP_FN_PTR,                },	/* order=6 (unused) */
-	    },
-	  },
-
-	  {
-	  /* N_dims = 1 */
-	    {
-	    /* molecule_family = molecule_family_cube */
-	      { NULL_INTERP_FN_PTR,                },	/* order=0 (unused) */
-	      { NULL_INTERP_FN_PTR,                },	/* order=1 (unused) */
-	      { LocalInterp_UHermite_1d_cube2,     },	/* order=2 */
-	      { LocalInterp_UHermite_1d_cube3,     },	/* order=3 */
-	      { LocalInterp_UHermite_1d_cube4,     },	/* order=4 */
-	      { NULL_INTERP_FN_PTR,                },	/* order=5 (n.i.) */
-	      { NULL_INTERP_FN_PTR,                },	/* order=6 (n.i.) */
-	    },
-	  },
-
-	  {
-	  /* N_dims = 2 */
-	    {
-	    /* molecule_family = molecule_family_cube */
-	      { NULL_INTERP_FN_PTR,                },	/* order=0 (unused) */
-	      { NULL_INTERP_FN_PTR,                },	/* order=1 (unused) */
-	      { LocalInterp_UHermite_2d_cube2,     },	/* order=2 */
-	      { LocalInterp_UHermite_2d_cube3,     },	/* order=3 */
-	      { NULL_INTERP_FN_PTR,                },	/* order=4 (n.i.) */
-	      { NULL_INTERP_FN_PTR,                },	/* order=5 (n.i.) */
-	      { NULL_INTERP_FN_PTR,                },	/* order=6 (n.i.) */
-	    },
-	  },
-
-	  {
-	  /* N_dims = 3 */
-	    {
-	    /* molecule_family = molecule_family_cube */
-	      { NULL_INTERP_FN_PTR,                },	/* order=0 (unused) */
-	      { NULL_INTERP_FN_PTR,                },	/* order=1 (unused) */
-	      { LocalInterp_UHermite_3d_cube2,     },	/* order=2 */
-	      { LocalInterp_UHermite_3d_cube3,     },	/* order=3 */
-	      { NULL_INTERP_FN_PTR,                },	/* order=4 (n.i.) */
-	      { NULL_INTERP_FN_PTR,                },	/* order=5 (n.i.) */
-	      { NULL_INTERP_FN_PTR,                },	/* order=6 (n.i.) */
-	    },
-	  },
-
-	/* end of interp_operator == interp_operator_Hermite */
+/* firewall: make sure all the lookup indices are in range */
+/* ... commented-out comparisons are always true since enums are unsigned */
+if (    /*(interp_operator >= 0) &&*/ (interp_operator <  N_INTERP_OPERATORS)
+     &&   (N_dims          >= 0) &&   (N_dims          <= MAX_N_DIMS)
+     && /*(molecule_family >= 0) &&*/ (molecule_family <  N_MOLECULE_FAMILIES)
+     &&   (order           >= 0) &&   (order           <= MAX_ORDER)
+     &&   (smoothing       >= 0) &&   (smoothing       <= MAX_SMOOTHING)    )
+   then {
+	/* ok ==> no-op here */
+	}
+   else {
+	CCTK_VWarn(BUG_MSG_SEVERITY_LEVEL,
+		   __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+"   CCTK_InterpLocalUniform():\n"
+"        internal error (interpolator bug):\n"
+"        one or more function-pointer-table indices is out of range!\n"
+"        interp_operator=(int)%d, should be in [0,%d)\n"
+"        N_dims         =     %d, should be in [0,%d]\n"
+"        molecule_family=(int)%d, should be in [0,%d)\n"
+"        order          =     %d, should be in [0,%d]\n"
+"        smoothing      =     %d, should be in [0,%d]\n"
+			   ,
+			   (int) interp_operator, (int) N_INTERP_OPERATORS,
+			   (int) N_dims         , (int) MAX_N_DIMS,
+			   (int) molecule_family, (int) N_MOLECULE_FAMILIES,
+			   (int) order          , (int) MAX_ORDER,
+			   (int) smoothing      , (int) MAX_SMOOTHING);
+								/*NOTREACHED*/
+	CCTK_Abort(NULL, BUG_ABORT_CODE);			/*NOTREACHED*/
 	}
-
-      };
 
 /* look up the subfunction to do the interpolation */
-p_interp_fn_t p_interp_fn = p_interp_fn_table[interp_operator]
-					     [N_dims]
-					     [molecule_family]
-					     [order]
-					     [smoothing];
+  {
+const p_interp_fn_t p_interp_fn = p_interp_fn_table[interp_operator]
+						   [N_dims]
+						   [molecule_family]
+						   [order]
+						   [smoothing];
 if (p_interp_fn == NULL_INTERP_FN_PTR)
    then {
-	free(operation_codes);
-	free(operand_indices);
 	free(input_array_offsets);
+	free(operand_indices);
+	free(operation_codes);
 	if (p_Jacobian_info != NULL)
 	   then {
-		free(Jacobian_info.Jacobian_pointer);
 		free(Jacobian_info.Jacobian_offset);
+		free(Jacobian_info.Jacobian_pointer);
 		}
 	CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
 		   __LINE__, __FILE__, CCTK_THORNSTRING,
@@ -1655,13 +1376,11 @@ 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_str, order, smoothing);
+		   molecule_family_string, order, smoothing);
 	return UTIL_ERROR_BAD_INPUT;			/*** ERROR RETURN ***/
 	}
 
-/*
- * call the subfunction to actually do the interpolation
- */
+/* call the subfunction to actually do the interpolation */
   {
 struct error_flags error_flags;
 struct molecule_structure_flags molecule_structure_flags;
@@ -1691,139 +1410,79 @@ const int return_code
 /*
  * is an interpolation point out of range?
  */
-if (return_code == CCTK_ERROR_INTERP_POINT_X_RANGE)
+if (return_code == CCTK_ERROR_INTERP_POINT_OUTSIDE)
    then {
-	/* store details about the error in the parameter table */
-
-	status1 = Util_TableSetInt(param_table_handle,
-				   error_flags.error_pt,
-				   "out_of_range_pt");
-	status2 = Util_TableSetInt(param_table_handle,
-				   error_flags.error_axis,
-				   "out_of_range_axis");
-	status3 = Util_TableSetInt(param_table_handle,
-				   error_flags.error_end,
-				   "out_of_range_end");
-	if ((status1 < 0) || (status2 < 0) || (status3 < 0))
+	status = set_error_info(param_table_handle,
+				&error_flags);
+	if (status != 0)
 	   then {
-		free(operation_codes);
-		free(operand_indices);
 		free(input_array_offsets);
+		free(operand_indices);
+		free(operation_codes);
 		if (p_Jacobian_info != NULL)
 		   then {
-			free(Jacobian_info.Jacobian_pointer);
 			free(Jacobian_info.Jacobian_offset);
+			free(Jacobian_info.Jacobian_pointer);
 			}
-		CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
-			   __LINE__, __FILE__, CCTK_THORNSTRING,
-"\n"
-"   CCTK_InterpLocalUniform(): error setting out-of-range table entry/entries!"
-     "                         error status1=%d status2=%d status3=%d",
-			   status1, status2, status3);
-		return   (status1 < 0) ? status1	/*** ERROR RETURN ***/
-		       : (status2 < 0) ? status2	/*** ERROR RETURN ***/
-		       :		 status3;	/*** ERROR RETURN ***/
+		return status;				/*** ERROR RETURN ***/
 		}
 	}
 
 /******************************************************************************/
 
 /*
- * store molecule structure info
+ * store query results
+ *
+ * ... only molecule structure and min/max m have to be stored explicitly;
+ *     the other queries are stored directly into their final destinatios
+ *     by the interpolation subfunction
  */
-status1 = Util_TableSetInt(param_table_handle,
-			   molecule_structure_flags.MSS_is_fn_of_interp_coords,
-			   "MSS_is_fn_of_interp_coords");
-status2 = Util_TableSetInt(param_table_handle,
-			   molecule_structure_flags
-			      .MSS_is_fn_of_which_operation,
-			   "MSS_is_fn_of_which_operation");
-status3 = Util_TableSetInt(param_table_handle,
-			   molecule_structure_flags
-			      .MSS_is_fn_of_input_array_values,
-			   "MSS_is_fn_of_input_array_values");
-status4 = Util_TableSetInt(param_table_handle,
-			   molecule_structure_flags
-			      .Jacobian_is_fn_of_input_array_values,
-			   "Jacobian_is_fn_of_input_array_values");
-if ((status1 < 0) || (status2 < 0) || (status3 < 0) || (status4 < 0))
+
+status = set_molecule_structure(param_table_handle,
+				&molecule_structure_flags);
+if (status != 0)
    then {
-	free(operation_codes);
-	free(operand_indices);
 	free(input_array_offsets);
+	free(operand_indices);
+	free(operation_codes);
 	if (p_Jacobian_info != NULL)
 	   then {
-		free(Jacobian_info.Jacobian_pointer);
 		free(Jacobian_info.Jacobian_offset);
+		free(Jacobian_info.Jacobian_pointer);
 		}
-	CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
-		   __LINE__, __FILE__, CCTK_THORNSTRING,
-"\n"
-"   CCTK_InterpLocalUniform(): error setting Jacobian structure\n"
-"                              table entry/entries!"
-"                              error status1=%d status2=%d\n"
-"                                    status3=%d status4=%d"
-			   ,
-			   status1, status2, status3, status4);
-	return   (status1 < 0) ? status1		/*** ERROR RETURN ***/
-	       : (status2 < 0) ? status2		/*** ERROR RETURN ***/
-	       : (status3 < 0) ? status3		/*** ERROR RETURN ***/
-	       :		 status4;		/*** ERROR RETURN ***/
+	return status;					/*** ERROR RETURN ***/
 	}
 
-/******************************************************************************/
-
-/*
- * store results from molecule min/max m queries
- *
- * ... molecule position queries are stored directly into
- *     their final destination by subfunction
- *     ==> no further store or cleanup needed here
- */
-
 if (p_molecule_min_max_m_info != NULL)
    then {
-	status1 = Util_TableSetIntArray(param_table_handle,
-					N_dims, p_molecule_min_max_m_info
-						   ->molecule_min_m,
-					"molecule_min_m");
-	status2 = Util_TableSetIntArray(param_table_handle,
-					N_dims, p_molecule_min_max_m_info
-						   ->molecule_max_m,
-					"molecule_max_m");
-	if ((status1 < 0) || (status2 < 0))
+	status = set_molecule_min_max_m(param_table_handle,
+					N_dims, p_molecule_min_max_m_info);
+	if (status != 0)
 	   then {
-		free(operation_codes);
-		free(operand_indices);
 		free(input_array_offsets);
+		free(operand_indices);
+		free(operation_codes);
 		if (p_Jacobian_info != NULL)
 		   then {
-			free(Jacobian_info.Jacobian_pointer);
 			free(Jacobian_info.Jacobian_offset);
+			free(Jacobian_info.Jacobian_pointer);
 			}
-		CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
-			   __LINE__, __FILE__, CCTK_THORNSTRING,
-"\n"
-"   CCTK_InterpLocalUniform(): error setting\n"
-"                              \"molecule_{min,max}_m\" table entry/entries!"
-"                              error status1=%d status2=%d",
-			   status1, status2);
-		return (status1 < 0) ? status1 : status2; /*** ERROR RETURN ***/
+		return status;				/*** ERROR RETURN ***/
 		}
 	}
 
 /******************************************************************************/
 
-free(operation_codes);
-free(operand_indices);
 free(input_array_offsets);
+free(operand_indices);
+free(operation_codes);
 if (p_Jacobian_info != NULL)
    then {
-	free(Jacobian_info.Jacobian_pointer);
 	free(Jacobian_info.Jacobian_offset);
+	free(Jacobian_info.Jacobian_pointer);
 	}
 
-return return_code;
+return return_code;					/*** NORMAL RETURN ***/
   }
   }
   }
@@ -1840,3 +1499,673 @@ return return_code;
   }
   }
 }
+
+/******************************************************************************/
+/******************************************************************************/
+/******************************************************************************/
+
+/*
+ * This function checks
+ *	CCTK_REAL boundary_off_centering_tolerance[N_boundaries]
+ *	CCTK_REAL boundary_extrapolation_tolerance[N_boundaries]
+ * for validity, warning about any dubious settings.
+ */
+static
+  void check_boundary_tolerances
+	(int N_boundaries,
+	 const CCTK_REAL boundary_off_centering_tolerance[MAX_N_BOUNDARIES],
+	 const CCTK_REAL boundary_extrapolation_tolerance[MAX_N_BOUNDARIES])
+{
+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():\n"
+"        warning: setting boundary_off_centering_tolerance[%d] == 0.0\n"
+"                 and     boundary_extrapolation_tolerance[%d] == %g > 0.0\n"
+"                 is almost certainly a mistake\n"
+"                 (the boundary_extrapolation_tolerance[%d] == %g\n"
+"                  setting will have no effect because of the\n"
+"                  boundary_off_centering_tolerance[%d] == 0.0 setting)"
+			   ,
+			   ibndry,
+			   ibndry,
+			      (double) boundary_extrapolation_tolerance[ibndry],
+			   ibndry,
+			      (double) boundary_extrapolation_tolerance[ibndry],
+			   ibndry);
+	}
+}
+
+/******************************************************************************/
+
+/*
+ * This function gets
+ *	const char molecule_family[]
+ * from the parameter table, or sets the default value "cube" if this
+ * key isn't present.  This function also decodes  molecule_family  into
+ * an  enum molecule_family .
+ *
+ * Results:
+ * This function returns 0 for ok, or the (nonzero) return code for error.
+ */
+static
+  int get_and_decode_molecule_family
+	(int param_table_handle,
+	 int buffer_size, char molecule_family_string_buffer[],
+	 enum molecule_family *p_molecule_family)
+{
+enum molecule_family molecule_family;
+int status;
+
+status = Util_TableGetString(param_table_handle,
+			     buffer_size, molecule_family_string_buffer,
+			     "molecule_family");
+
+if	(status == UTIL_ERROR_TABLE_NO_SUCH_KEY)
+   then {
+	/* set the default value */
+	LocalInterp_Strlcpy(molecule_family_string_buffer, "cube", buffer_size);
+	molecule_family = molecule_family_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,
+				     molecule_family_string_buffer,
+				     "molecule_family");
+	if (status < 0)
+	   then {
+		CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
+			   __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+"   CCTK_InterpLocalUniform(): error setting default molecule family\n"
+"                              in parameter table!\n"
+"                              error status=%d",
+			   status);
+		return status;				/*** ERROR RETURN ***/
+		}
+	}
+else if (status > 0)
+   then {
+	/* decode molecule family string */
+	if (strcmp(molecule_family_string_buffer, "cube") == 0)
+	   then molecule_family = molecule_family_cube;
+	else	{
+		CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
+			   __LINE__, __FILE__, CCTK_THORNSTRING,
+"CCTK_InterpLocalUniform(): unknown molecule_family string \"%s\"!",
+			   molecule_family_string_buffer);
+		return UTIL_ERROR_BAD_INPUT;		/*** ERROR RETURN ***/
+		}
+	}
+else	{
+	CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
+		   __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+"   CCTK_InterpLocalUniform(): bad table entry for\n"
+"                              \"molecule_family\" parameter!\n"
+"                              error status=%d",
+		   status);
+	return status;					/*** ERROR RETURN ***/
+	}
+
+*p_molecule_family = molecule_family;
+return 0;						/*** NORMAL RETURN ***/
+}
+
+/******************************************************************************/
+
+/*
+ * This function gets
+ *	CCTK_INT* molecule_positions[N_dims]
+ * from the parameter table.
+ *
+ * Results:
+ * This function returns 0 for ok, or the (nonzero) return code for error.
+ */
+static
+  int get_molecule_positions
+	(int param_table_handle,
+	 int N_dims, CCTK_INT* molecule_positions_array[MAX_N_DIMS])
+{
+CCTK_POINTER molecule_positions_temp[MAX_N_DIMS];
+int status;
+
+status = Util_TableGetPointerArray(param_table_handle,
+				   N_dims, molecule_positions_temp,
+				   "molecule_positions");
+
+if (status == N_dims)
+   then {
+	/* type-cast CCTK_POINTER pointers to CCTK_INT* pointers */
+	/* (which point to where the query results will be stored) */
+	int axis;
+		for (axis = 0 ; axis < N_dims ; ++axis)
+		{
+		molecule_positions_array[axis]
+			= (CCTK_INT *) molecule_positions_temp[axis];
+		}
+	}
+else	{
+	CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
+		   __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+"   CCTK_InterpLocalUniform(): bad \"molecule_positions\" table entry!\n"
+"                              error status=%d",
+		   status);
+	return status;					/*** ERROR RETURN ***/
+	}
+
+return 0;						/*** NORMAL RETURN ***/
+}
+
+/******************************************************************************/
+
+/*
+ * This function gets the Jacobian-query parameters
+ *	CCTK_REAL* Jacobian_pointer[N_output_arrays]
+ *	CCTK_INT   Jacobian_offset [N_output_arrays]   # optional, default=all 0
+ *	CCTK_INT   Jacobian_interp_point_stride
+ *	CCTK_INT   Jacobian_m_point_strides[N_dims]
+ *	CCTK_INT   Jacobian_part_stride                # optional, default=1
+ * and stores them in a  struct Jacobian_info .
+ *
+ * Results:
+ * This function returns 0 for ok, or the (nonzero) return code for error.
+ */
+static
+  int get_Jacobian_info(int param_table_handle,
+			int N_dims, int N_output_arrays,
+			struct Jacobian_info* p_Jacobian_info)
+{
+/* padded array size, cf. LocalInterp_InterpLocalUniform() header comments */
+const int N_output_arrays1 = N_output_arrays + 1;
+
+int status;
+
+p_Jacobian_info->Jacobian_pointer
+	= (CCTK_REAL**) malloc(N_output_arrays1 * sizeof(CCTK_REAL *));
+if (p_Jacobian_info->Jacobian_pointer == NULL)
+   then return UTIL_ERROR_NO_MEMORY;			/*** ERROR RETURN ***/
+
+  {
+CCTK_POINTER* Jacobian_pointer_temp
+	= (CCTK_POINTER*) malloc(N_output_arrays1 * sizeof(CCTK_POINTER));
+if (Jacobian_pointer_temp == NULL)
+   then {
+	free(p_Jacobian_info->Jacobian_pointer);
+	return UTIL_ERROR_NO_MEMORY;			/*** ERROR RETURN ***/
+	}
+status = Util_TableGetPointerArray(param_table_handle,
+				   N_output_arrays, Jacobian_pointer_temp,
+				   "Jacobian_pointer");
+if	(status == N_output_arrays)
+   then {
+	/* type-cast CCTK_POINTER pointers to CCTK_REAL* pointers */
+	/* (which point to where the query results will be stored) */
+	int out;
+		for (out = 0 ; out < N_output_arrays ; ++out)
+		{
+		p_Jacobian_info->Jacobian_pointer[out]
+			= (CCTK_REAL *) Jacobian_pointer_temp[out];
+		}
+	free(Jacobian_pointer_temp);
+	}
+else	{
+	free(p_Jacobian_info->Jacobian_pointer);
+	free(Jacobian_pointer_temp);
+	CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
+		   __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+"   CCTK_InterpLocalUniform(): bad \"Jacobian_pointer\" table entry!\n"
+"                              error status=%d"
+		   ,
+		   status);
+	return status;					/*** ERROR RETURN ***/
+	}
+  }
+
+p_Jacobian_info->Jacobian_offset
+	= (CCTK_INT*) malloc(N_output_arrays1 * sizeof(CCTK_INT));
+if (p_Jacobian_info->Jacobian_offset == NULL)
+   then {
+	free(p_Jacobian_info->Jacobian_pointer);
+	return UTIL_ERROR_NO_MEMORY;			/*** ERROR RETURN ***/
+	}
+status = get_INT_array(param_table_handle, "Jacobian_offset",
+		       true, 0,		/* default value */
+		       N_output_arrays, p_Jacobian_info->Jacobian_offset,
+		       NULL);
+if (status != 0)
+   then {
+	free(p_Jacobian_info->Jacobian_pointer);
+	free(p_Jacobian_info->Jacobian_offset);
+	return status;					/*** ERROR RETURN ***/
+	}
+
+status = get_and_check_INT(param_table_handle, "Jacobian_interp_point_stride",
+			   true, 0,	/* this is a mandatory parameter */
+					/* if we're executing this function */
+			   false, 0, 0, NULL,	/* no range check */
+			   & p_Jacobian_info->Jacobian_interp_point_stride);
+if (status != 0)
+   then {
+	free(p_Jacobian_info->Jacobian_offset);
+	free(p_Jacobian_info->Jacobian_pointer);
+	return status;					/*** ERROR RETURN ***/
+	}
+
+status = get_INT_array(param_table_handle, "Jacobian_m_strides",
+		       false, 0,	/* don't set default value */
+		       N_dims, p_Jacobian_info->Jacobian_m_strides,
+		       NULL);		/* this is a mandatory parameter */
+					/* if we're executing this function */
+if (status != 0)
+   then {
+	free(p_Jacobian_info->Jacobian_pointer);
+	free(p_Jacobian_info->Jacobian_offset);
+	return status;					 /*** ERROR RETURN ***/
+	}
+
+status = get_and_check_INT(param_table_handle, "Jacobian_part_stride",
+			   false, 1,	/* default value */
+			   false, 0, 0, NULL,	/* no range check */
+			   & p_Jacobian_info->Jacobian_part_stride);
+if (status != 0)
+   then {
+	free(p_Jacobian_info->Jacobian_pointer);
+	free(p_Jacobian_info->Jacobian_offset);
+	return status;					/*** ERROR RETURN ***/
+	}
+
+return 0;						/*** NORMAL RETURN ***/
+}
+
+/******************************************************************************/
+
+/*
+ * This function sets the error-info entries
+ *	CCTK_INT error_pt
+ *	CCTK_INT error_ibndry
+ *	CCTK_INT error_axis
+ *	CCTK_INT error_direction
+ * in the parameter table.
+ *
+ * Results:
+ * This function returns 0 for ok, or the (nonzero) return code for error.
+ */
+static
+  int set_error_info(int param_table_handle,
+		     struct error_flags* p_error_flags)
+{
+const int status1 = Util_TableSetInt(param_table_handle,
+				     p_error_flags->error_pt,
+				     "error_pt");
+const int status2 = Util_TableSetInt(param_table_handle,
+				     p_error_flags->error_ibndry,
+				     "error_ibndry");
+const int status3 = Util_TableSetInt(param_table_handle,
+				     p_error_flags->error_axis,
+				     "error_axis");
+const int status4 = Util_TableSetInt(param_table_handle,
+				     p_error_flags->error_direction,
+				     "error_direction");
+if ((status1 < 0) || (status2 < 0) || (status3 < 0) || (status4 < 0))
+   then {
+	CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
+		   __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+"   CCTK_InterpLocalUniform():\n"
+"        error setting error-info table entry/entries!\n"
+"        error status1=%d status2=%d status3=%d status4=%d"
+		   ,
+		   status1, status2, status3, status4);
+	return   (status1 < 0) ? status1		/*** ERROR RETURN ***/
+	       : (status2 < 0) ? status2		/*** ERROR RETURN ***/
+	       : (status3 < 0) ? status3		/*** ERROR RETURN ***/
+	       :		 status4;		/*** ERROR RETURN ***/
+	}
+
+return 0;						/*** NORMAL RETURN ***/
+}
+
+/******************************************************************************/
+
+/*
+ * This function sets the molecule structure flags
+ *	CCTK_INT MSS_is_fn_of_interp_coords
+ *	CCTK_INT MSS_is_fn_of_which_operation
+ *	CCTK_INT MSS_is_fn_of_input_array_values
+ *	CCTK_INT Jacobian_is_fn_of_input_array_values
+ * in the parameter table.
+ *
+ * Results:
+ * This function returns 0 for ok, or the (nonzero) return code for error.
+ */
+static
+  int set_molecule_structure
+	(int param_table_handle,
+	 const struct molecule_structure_flags* p_molecule_structure_flags)
+{
+const int status1
+   = Util_TableSetInt(param_table_handle,
+		      p_molecule_structure_flags->MSS_is_fn_of_interp_coords,
+		      "MSS_is_fn_of_interp_coords");
+const int status2
+   = Util_TableSetInt(param_table_handle,
+		      p_molecule_structure_flags->MSS_is_fn_of_which_operation,
+		      "MSS_is_fn_of_which_operation");
+const int status3
+   = Util_TableSetInt(param_table_handle,
+		      p_molecule_structure_flags->MSS_is_fn_of_input_array_values,
+		      "MSS_is_fn_of_input_array_values");
+const int status4
+   = Util_TableSetInt(param_table_handle,
+		      p_molecule_structure_flags->Jacobian_is_fn_of_input_array_values,
+		      "Jacobian_is_fn_of_input_array_values");
+if ((status1 < 0) || (status2 < 0) || (status3 < 0) || (status4 < 0))
+   then {
+	CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
+		   __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+"   CCTK_InterpLocalUniform():\n"
+"        error setting molecule structure flags table entry/entries!"
+"        error status1=%d status2=%d status3=%d status4=%d"
+			   ,
+			   status1, status2, status3, status4);
+	return   (status1 < 0) ? status1		/*** ERROR RETURN ***/
+	       : (status2 < 0) ? status2		/*** ERROR RETURN ***/
+	       : (status3 < 0) ? status3		/*** ERROR RETURN ***/
+	       :		 status4;		/*** ERROR RETURN ***/
+	}
+
+return 0;						/*** NORMAL RETURN ***/
+}
+
+/******************************************************************************/
+
+/*
+ * This function sets the molecule size parameters
+ *	CCTK_INT molecule_min_m[N_dims]
+ *	CCTK_INT molecule_max_m[N_dims]
+ * in the parameter table.
+ *
+ * Results:
+ * This function returns 0 for ok, or the (nonzero) return code for error.
+ */
+static
+  int set_molecule_min_max_m
+	(int param_table_handle,
+	 int N_dims,
+	 const struct molecule_min_max_m_info* p_molecule_min_max_m_info)
+{
+const int status1
+   = Util_TableSetIntArray(param_table_handle,
+			   N_dims, p_molecule_min_max_m_info->molecule_min_m,
+			   "molecule_min_m");
+const int status2
+   = Util_TableSetIntArray(param_table_handle,
+			   N_dims, p_molecule_min_max_m_info->molecule_max_m,
+			   "molecule_max_m");
+if ((status1 < 0) || (status2 < 0))
+   then {
+	CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
+		   __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+"   CCTK_InterpLocalUniform(): error setting\n"
+"                              \"molecule_{min,max}_m\" table entry/entries!"
+"                              error status1=%d status2=%d",
+		   status1, status2);
+	return (status1 < 0) ? status1 : status2;	/*** ERROR RETURN ***/
+	}
+
+return 0;						/*** NORMAL RETURN ***/
+}
+
+/******************************************************************************/
+/******************************************************************************/
+/******************************************************************************/
+
+/*
+ * This function gets and range-checks a CCTK_INT from the parameter table.
+ *
+ * Arguments:
+ * param_table_handle = handle to the Cactus key-value table
+ * name = the character-string key in the table
+ * mandatory_flag = true  ==> the value must be present in the parameter
+ *			      table (default_value is ignored)
+ *		    false ==> the value is optional
+ * default_value = the default value (for each array element) if the
+ *		   key isn't in the parameter table and mandatory_flag == false
+ * check_range_flag = true  ==> check that the value satisfies
+ *				min_value <= value <= max_value
+ *		      false ==> don't do a range check, and ignore
+ *				min_value, max_value, and max_value_string
+ * {min,max}_value = the inclusive range of valid values
+ * max_value_string = the character-string name of max_value (this is used
+ *		       only in formatting the out-of-range error message)
+ * p_value --> where we should store the value
+ *
+ * Results:
+ * This function returns 0 for ok, or the (nonzero) return code for error.
+ */
+static
+  int get_and_check_INT
+	(int param_table_handle, const char name[],
+	 bool mandatory_flag, int default_value,
+	 bool check_range_flag, int min_value, int max_value,
+				const char max_value_string[],
+	 CCTK_INT* p_value)
+{
+CCTK_INT value;
+
+const int status = Util_TableGetInt(param_table_handle, &value, name);
+
+if	(status == 1)
+   then { }	/* value set from table ==> no-op here */
+else if (status == UTIL_ERROR_TABLE_NO_SUCH_KEY)
+   then {
+	if (mandatory_flag)
+	   then {
+		CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
+			   __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+"   CCTK_InterpLocalUniform(): missing table entry for\n"
+"                              \"%s\" parameter!\n"
+"                              (this is a mandatory parameter)!\n"
+			   ,
+			   name);
+		return status;				/*** ERROR RETURN ***/
+		}
+	   else value = default_value;
+	}
+else	{
+	CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
+		   __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+"   CCTK_InterpLocalUniform(): bad table entry for\n"
+"                              \"%s\" parameter!\n"
+"                              error status=%d"
+		   ,
+		   name,
+		   status);
+	return status;					/*** ERROR RETURN ***/
+	}
+
+if (check_range_flag)
+   then {
+	if ((value < min_value) || (value > max_value))
+	   then {
+		CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
+			   __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+"   CCTK_InterpLocalUniform(): %s=%d is invalid!\n"
+"                              valid range is %d <= %s <= %s=%d",
+			   name, value,
+			   min_value, name, max_value_string, max_value);
+		return UTIL_ERROR_BAD_INPUT;		/*** ERROR RETURN ***/
+		}
+	}
+
+*p_value = value;
+return 0;						/*** NORMAL RETURN ***/
+}
+
+/******************************************************************************/
+
+/*
+ * This function gets an array of CCTK_INT values from the parameter table.
+ *
+ * Arguments:
+ * param_table_handle = handle to the Cactus key-value table
+ * name = the character-string key in the table
+ * default_flag = false ==> ignore  default_value  and set
+ *			     *p_value_not_set  to true if the
+ *			    key is not present in the parameter table,
+ *			    false if it is present
+ *		  true  ==> if the key is not present in the parameter table,
+ *			    set each array element to  default_value
+ * default_value = the default value (for each array element) if the
+ *		   key isn't in the parameter table
+ * N = the size of the array
+ * buffer --> a buffer in which to store the array
+ * p_value_not_set = NULL or a pointer to a flag to be set if and only if
+ *		     this function did *not* set values in the buffer.
+ *		     More precisely, if  default_flag  is set, this argument
+ *		     is ignored.  Otherwise, the logic is this:
+ *			if (this function set values in the buffer)
+ *			   then {
+ *				if (p_value_not_set is non-NULL)
+ *				   then *p_value_not_set = false
+ *				}
+ *			   else {
+ *				if (p_value_not_set is non-NULL)
+ *				   then *p_value_not_set = true
+ *				   else give an error message that
+ *					a mandatory value is missing
+ *					from the parameter table, then
+ *					return UTIL_ERROR_TABLE_NO_SUCH_KEY
+ *				}
+ *
+ * Results:
+ * This function returns 0 for ok, or the desired (nonzero) error return
+ * code if something is wrong.  Note that error return codes may be either
+ * positive or negative!
+ */
+static
+  int get_INT_array(int param_table_handle, const char name[],
+		    bool default_flag, int default_value,
+		    int N, CCTK_INT buffer[],
+		    bool* p_value_not_set)
+{
+const int status
+	= Util_TableGetIntArray(param_table_handle,
+				N, buffer,
+				name);
+
+if	 (status == UTIL_ERROR_TABLE_NO_SUCH_KEY)
+   then {
+	/* we didn't set values in the buffer */
+	/* ==> either set the default value ourself, */
+	/*     or print an error message, */
+	/*     or notify our caller that we didn't do so */
+	if (default_flag)
+	   then {
+		int i;
+			for (i = 0 ; i < N ; ++i)
+			{
+			buffer[i] = default_value;
+			}
+		}
+	   else {
+		if (p_value_not_set == NULL)
+		   then {
+			CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
+				   __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+"   CCTK_InterpLocalUniform(): missing table entry for\n"
+"                              \"%s\" parameter!\n"
+"                              (this is a mandatory parameter)!\n"
+				   ,
+				   name);
+			return UTIL_ERROR_TABLE_NO_SUCH_KEY;  /*** ERROR RETURN ***/
+			}
+		   else *p_value_not_set = true;
+		}
+	}
+else if (status == N)
+   then {
+	/* we did set values in the buffer */
+	if (! default_flag)
+	   then if (p_value_not_set != NULL)
+		   then *p_value_not_set = false;
+	}
+else	{
+	CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
+		   __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+"   CCTK_InterpLocalUniform(): bad table entry for\n"
+"                              \"%s\" parameter!\n"
+"                              error status=%d",
+		   name, status);
+	return status;					/*** ERROR RETURN ***/
+	}
+
+return 0;						/*** NORMAL RETURN ***/
+}
+
+/******************************************************************************/
+
+/*
+ * This function gets an array of CCTK_REAL values from the parameter table.
+ *
+ * Arguments:
+ * param_table_handle = handle to the Cactus key-value table
+ * name = the character-string key in the table
+ * default_value = the default value (for each array element) if the
+ *		   key isn't in the parameter table
+ * N = the size of the array
+ * buffer --> a buffer in which to store the array
+ *
+ * Results:
+ * This function returns 0 for ok, or the desired (nonzero) error return
+ * code if something is wrong.  Note that error return codes may be either
+ * positive or negative!
+ */
+static
+  int get_REAL_array(int param_table_handle, const char name[],
+		      CCTK_REAL default_value,
+		      int N, CCTK_REAL buffer[])
+{
+const int status
+	= Util_TableGetRealArray(param_table_handle,
+				 N, buffer,
+				 name);
+if	 (status == UTIL_ERROR_TABLE_NO_SUCH_KEY)
+   then {
+	/* set default value */
+	int i;
+		for (i = 0 ; i < N ; ++i)
+		{
+		buffer[i] = default_value;
+		}
+	}
+else if (status == N)
+   then { }	/* ok ==> no-op here */
+else	{
+	CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
+		   __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+"   CCTK_InterpLocalUniform(): bad table entry for\n"
+"                              \"%s\" parameter!\n"
+"                              error status=%d",
+		   name, status);
+	return status;					/*** ERROR RETURN ***/
+	}
+
+return 0;						/*** NORMAL RETURN ***/
+}
diff --git a/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.h b/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.h
index 851378a..8d9c551 100644
--- a/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.h
+++ b/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.h
@@ -109,15 +109,26 @@ enum	molecule_family
  * other compile-time settings
  */
 
-/* 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
+/* defaults for boundary_{off_centering,extrapolation}_tolerance[] */
+#ifdef NOT_YET
+  #define LAGRANGE_BNDRY_OFF_CNTR_TOL_DEF	999.0
+  #define LAGRANGE_BNDRY_EXTRAP_TOL_DEF		1.0e-10
+  #define HERMITE_BNDRY_OFF_CNTR_TOL_DEF	1.0e-10
+  #define HERMITE_BNDRY_EXTRAP_TOL_DEF		0.0
+#else
+  #define LAGRANGE_BNDRY_OFF_CNTR_TOL_DEF	999.0
+  #define LAGRANGE_BNDRY_EXTRAP_TOL_DEF		1.0e-10
+  #define HERMITE_BNDRY_OFF_CNTR_TOL_DEF	999.0
+  #define HERMITE_BNDRY_EXTRAP_TOL_DEF		1.0e-10
+#endif
 
 /* CCTK_VWarn() severity level for error/warning messages */
-#define ERROR_MSG_SEVERITY_LEVEL	0
-#define WARNING_MSG_SEVERITY_LEVEL	1
+#define BUG_MSG_SEVERITY_LEVEL			0
+#define ERROR_MSG_SEVERITY_LEVEL		0
+#define WARNING_MSG_SEVERITY_LEVEL		1
+
+/* CCTK_Abort() exit code for internal error (interpolator bug) aborts */
+#define BUG_ABORT_CODE				42
 
 /******************************************************************************/
 
@@ -133,7 +144,7 @@ struct	error_flags
 	int error_pt;
 	int error_ibndry;
 	int error_axis;
-	int error_end;
+	int error_direction;
 	};
 
 struct	molecule_structure_flags
@@ -230,5 +241,5 @@ int LocalInterp_molecule_posn(fp grid_origin, fp grid_delta,
 			      int* i_center, fp* x_rel);
 
 /* functions in util.c */
-void LocalInterp_zero_int_array(int N, CCTK_INT array[]);
 int LocalInterp_decode_N_parts(int type_code);
+size_t LocalInterp_Strlcpy(char* dst, const char* src, size_t dst_size);