Switch to using new version of Lagrange polynomial interpolation

which has separate tensor-product and maximum-degree interpolation
operators, with the former now the default.

This should fix the problem Yosef Zlowcher (sp?) found with
the interpolant not being continuous in multiple dimensions.


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

5 files changed, 276 insertions, 81 deletions

diff --git a/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c b/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c
index d0b4c78..9e0bf9f 100644
--- a/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c
+++ b/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c
@@ -3,8 +3,9 @@
 /*
 ** *****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_ILU_Lagrange_TP - driver for Lagrange interpolator (tensor prod)
+ * LocalInterp_ILU_Lagrange_MD - driver for Lagrange interpolator (max degree)
+ * LocalInterp_ILU_Hermite - driver point for Hermite interpolator
 **
 ** LocalInterp_InterpLocalUniform - main driver routine
 **
@@ -56,7 +57,8 @@
 #include "util_Table.h"
 #include "InterpLocalUniform.h"
 
-#include "Lagrange/all_prototypes.h"
+#include "Lagrange-tensor-product/all_prototypes.h"
+#include "Lagrange-maximum-degree/all_prototypes.h"
 #include "Hermite/all_prototypes.h"
 
 /* the rcs ID and its dummy function to use it */
@@ -78,8 +80,11 @@ CCTK_FILEVERSION(CactusBase_LocalInterp_src_GeneralizedPolynomialUniform_InterpL
 /* this enum describes which interpolation operator we're using */
 enum	interp_operator
 	{
-	interp_operator_Lagrange,     /* "Lagrange polynomial interpolation" */
-	interp_operator_Hermite,      /* "Hermite polynomial interpolation" */
+	interp_operator_Lagrange_TP,	/* "Lagrange polynomial interpolation */
+					/*  (tensor product)" */
+	interp_operator_Lagrange_MD,	/* "Lagrange polynomial interpolation
+					/*  (maximum degree)" */
+	interp_operator_Hermite,	/* "Hermite polynomial interpolation" */
 	N_INTERP_OPERATORS	/* this must be the last entry in the enum */
 	};
 
@@ -185,6 +190,8 @@ typedef
  * 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.)
+ *
+ * In the comments on the initializers, "n.i." = "not implemented".
  */
 static const p_interp_fn_t p_interp_fn_table[N_INTERP_OPERATORS]
 					    [MAX_N_DIMS+1]	/* see above */
@@ -193,9 +200,69 @@ static const p_interp_fn_t p_interp_fn_table[N_INTERP_OPERATORS]
 					    [MAX_SMOOTHING+1]
     = {
 
-      /* in the comments on the initializers, "n.i." = "not implemented" */
 	{
-	/* interp_operator == interp_operator_Lagrange */
+	/* interp_operator == interp_operator_Lagrange_TP */
+	  {
+	  /* 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_U_LagrTP_1d_cube10,   },	/* order=1 */
+	      { LocalInterp_U_LagrTP_1d_cube20,   },	/* order=2 */
+	      { LocalInterp_U_LagrTP_1d_cube30,   },	/* order=3 */
+	      { LocalInterp_U_LagrTP_1d_cube40,   },	/* order=4 */
+	      { LocalInterp_U_LagrTP_1d_cube50,   },	/* order=5 */
+	      { LocalInterp_U_LagrTP_1d_cube60,   },	/* order=6 */
+	    },
+	  },
+
+	  {
+	  /* N_dims = 2 */
+	    {
+	    /* molecule_family = molecule_family_cube */
+	      { NULL_INTERP_FN_PTR,                },	/* order=0 (unused) */
+	      { LocalInterp_U_LagrTP_2d_cube10,   },	/* order=1 */
+	      { LocalInterp_U_LagrTP_2d_cube20,   },	/* order=2 */
+	      { LocalInterp_U_LagrTP_2d_cube30,   },	/* order=3 */
+	      { LocalInterp_U_LagrTP_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_U_LagrTP_3d_cube10,   },	/* order=1 */
+	      { LocalInterp_U_LagrTP_3d_cube20,   },	/* order=2 */
+	      { LocalInterp_U_LagrTP_3d_cube30,   },	/* order=3 */
+	      { LocalInterp_U_LagrTP_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_TP */
+	},
+
+	{
+	/* interp_operator == interp_operator_Lagrange_MD */
 	  {
 	  /* N_dims = 0 ==> unused */
 	    {
@@ -215,12 +282,12 @@ static const p_interp_fn_t p_interp_fn_table[N_INTERP_OPERATORS]
 	    {
 	    /* 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 */
+	      { LocalInterp_U_LagrMD_1d_cube10,   },	/* order=1 */
+	      { LocalInterp_U_LagrMD_1d_cube20,   },	/* order=2 */
+	      { LocalInterp_U_LagrMD_1d_cube30,   },	/* order=3 */
+	      { LocalInterp_U_LagrMD_1d_cube40,   },	/* order=4 */
+	      { LocalInterp_U_LagrMD_1d_cube50,   },	/* order=5 */
+	      { LocalInterp_U_LagrMD_1d_cube60,   },	/* order=6 */
 	    },
 	  },
 
@@ -229,10 +296,10 @@ static const p_interp_fn_t p_interp_fn_table[N_INTERP_OPERATORS]
 	    {
 	    /* 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 */
+	      { LocalInterp_U_LagrMD_2d_cube10,   },	/* order=1 */
+	      { LocalInterp_U_LagrMD_2d_cube20,   },	/* order=2 */
+	      { LocalInterp_U_LagrMD_2d_cube30,   },	/* order=3 */
+	      { LocalInterp_U_LagrMD_2d_cube40,   },	/* order=4 */
 	      { NULL_INTERP_FN_PTR,                },	/* order=5 (n.i.) */
 	      { NULL_INTERP_FN_PTR,                },	/* order=6 (n.i.) */
 	    },
@@ -243,16 +310,16 @@ static const p_interp_fn_t p_interp_fn_table[N_INTERP_OPERATORS]
 	    {
 	    /* 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 */
+	      { LocalInterp_U_LagrMD_3d_cube10,   },	/* order=1 */
+	      { LocalInterp_U_LagrMD_3d_cube20,   },	/* order=2 */
+	      { LocalInterp_U_LagrMD_3d_cube30,   },	/* order=3 */
+	      { LocalInterp_U_LagrMD_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 */
+	/* end of interp_operator == interp_operator_Lagrange_MD */
 	},
 
 	{
@@ -314,7 +381,7 @@ static const p_interp_fn_t p_interp_fn_table[N_INTERP_OPERATORS]
 	  },
 
 	/* end of interp_operator == interp_operator_Hermite */
-	}
+	},
 
       };
 
@@ -323,7 +390,7 @@ static const p_interp_fn_t p_interp_fn_table[N_INTERP_OPERATORS]
 /******************************************************************************/
 
 /*@@
-  @routine    LocalInterp_ILU_Lagrange
+  @routine    LocalInterp_ILU_Lagrange_TP
   @date       22 Oct 2001
   @author     Jonathan Thornburg <jthorn@aei.mpg.de>
   @desc
@@ -332,10 +399,85 @@ static const p_interp_fn_t p_interp_fn_table[N_INTERP_OPERATORS]
 
 	It does Lagrange interpolation of a set of input arrays defined
 	on a uniform N-dimensional tensor-product grid, to a set of
-	interpolation points.  It can also do differentiation and/or
-	smoothing in the process of the interpolation.  It can also
-	return information about the Jacobian of the interpolated
-	values with respect to the input arrays.
+	interpolation points.  For multidimensional interpolation, it
+	uses the tensor-product basis for the interpolating function,
+	so the interpolant is guaranteed to be continuous and to pass
+	through the input data at the grid points (up to floating-point
+	roundoff errors).
+
+	This interpolator can also
+	* do differentiation and/or smoothing in the process of the
+	  interpolation
+	* return information about the Jacobian of the interpolated
+	  values with respect to the input arrays
+
+	This function does nothing except pass all its arguments
+	down to  LocalInterp_InterpLocalUniform()  with an extra
+	2 arguments added to indicate the operator handle.  See that
+	function for details on this function's arguments/results.
+  @enddesc
+  @@*/
+int LocalInterp_ILU_Lagrange_TP(int N_dims,
+				int param_table_handle,
+				/***** coordinate system *****/
+				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[],
+				/***** 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 *****/
+				int N_output_arrays,
+				const CCTK_INT output_array_type_codes[],
+				void *const output_arrays[])
+{
+return LocalInterp_InterpLocalUniform(interp_operator_Lagrange_TP,
+				      "Lagrange polynomial interpolation (tensor product)",
+				      LAGRANGE_BNDRY_OFF_CNTR_TOL_DEF,
+				      LAGRANGE_BNDRY_EXTRAP_TOL_DEF,
+				      N_dims,
+				      param_table_handle,
+				      coord_origin, coord_delta,
+				      N_interp_points,
+					interp_coords_type_code,
+					interp_coords,
+				      N_input_arrays,
+					input_array_dims,
+					input_array_type_codes,
+					input_arrays,
+				      N_output_arrays,
+					output_array_type_codes,
+					output_arrays);
+}
+
+/******************************************************************************/
+
+/*@@
+  @routine    LocalInterp_ILU_Lagrange_MD
+  @date       22 Oct 2001
+  @author     Jonathan Thornburg <jthorn@aei.mpg.de>
+  @desc
+	This function is a top-level driver for
+	LocalInterp::CCTK_InterpLocalUniform().
+
+	It does Lagrange interpolation of a set of input arrays defined
+	on a uniform N-dimensional tensor-product grid, to a set of
+	interpolation points.  For multidimensional interpolation, it
+	uses the maximum-degree basis for the interpolating function.
+	This means the interpolant is in general *not* continuous and
+	in general does *not* pass through the input data at the grid
+	points.
+
+	This interpolator can also
+	* do differentiation and/or smoothing in the process of the
+	  interpolation
+	* return information about the Jacobian of the interpolated
+	  values with respect to the input arrays
 
 	This function does nothing except pass all its arguments
 	down to  LocalInterp_InterpLocalUniform()  with an extra
@@ -343,27 +485,27 @@ static const p_interp_fn_t p_interp_fn_table[N_INTERP_OPERATORS]
 	function for details on this function's arguments/results.
   @enddesc
   @@*/
-int LocalInterp_ILU_Lagrange(int N_dims,
-			     int param_table_handle,
-			     /***** coordinate system *****/
-			     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[],
-			     /***** 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 *****/
-			     int N_output_arrays,
-			     const CCTK_INT output_array_type_codes[],
-			     void *const output_arrays[])
+int LocalInterp_ILU_Lagrange_MD(int N_dims,
+				int param_table_handle,
+				/***** coordinate system *****/
+				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[],
+				/***** 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 *****/
+				int N_output_arrays,
+				const CCTK_INT output_array_type_codes[],
+				void *const output_arrays[])
 {
-return LocalInterp_InterpLocalUniform(interp_operator_Lagrange,
-				      "Lagrange polynomial interpolation",
+return LocalInterp_InterpLocalUniform(interp_operator_Lagrange_MD,
+				      "Lagrange polynomial interpolation (maximum degree)",
 				      LAGRANGE_BNDRY_OFF_CNTR_TOL_DEF,
 				      LAGRANGE_BNDRY_EXTRAP_TOL_DEF,
 				      N_dims,
diff --git a/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.h b/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.h
index 8d9c551..d4821c4 100644
--- a/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.h
+++ b/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.h
@@ -192,24 +192,42 @@ struct	Jacobian_info
  */
 
 /* functions in InterpLocalUniform.c */
-int LocalInterp_ILU_Lagrange(int N_dims,
-			     int param_table_handle,
-			     /***** coordinate system *****/
-			     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[],
-			     /***** 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 *****/
-			     int N_output_arrays,
-			     const CCTK_INT output_array_type_codes[],
-			     void *const output_arrays[]);
+int LocalInterp_ILU_Lagrange_TP(int N_dims,
+				int param_table_handle,
+				/***** coordinate system *****/
+				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[],
+				/***** 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 *****/
+				int N_output_arrays,
+				const CCTK_INT output_array_type_codes[],
+				void *const output_arrays[]);
+int LocalInterp_ILU_Lagrange_MD(int N_dims,
+				int param_table_handle,
+				/***** coordinate system *****/
+				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[],
+				/***** 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 *****/
+				int N_output_arrays,
+				const CCTK_INT output_array_type_codes[],
+				void *const output_arrays[]);
 int LocalInterp_ILU_Hermite(int N_dims,
 			    int param_table_handle,
 			    /***** coordinate system *****/
diff --git a/src/GeneralizedPolynomial-Uniform/README b/src/GeneralizedPolynomial-Uniform/README
index b90cc26..15c6d27 100644
--- a/src/GeneralizedPolynomial-Uniform/README
+++ b/src/GeneralizedPolynomial-Uniform/README
@@ -3,8 +3,9 @@ $Header$
 This directory contains the interpolator registered for
    CCTK_InterpLocalUniform()
 under the names
-   "Lagrange polynomial interpolation"
-and
+   "Lagrange polynomial interpolation"	# a synonym for the next one...
+   "Lagrange polynomial interpolation (tensor product)"
+   "Lagrange polynomial interpolation (maximum degree)"
    "Hermite polynomial interpolation"
 
 
@@ -47,10 +48,13 @@ The source code files are as follows:
   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.
+* Lagrange/ and [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-tensor-product/ contains the code for the Lagrange polynomial
+  interpolator using tensor-product bases for multiple dimensions.
+* Lagrange-maximum-degree/ contains the code for the Lagrange polynomial
+  interpolator using maximum-degree bases for multiple dimensions.
 * Hermite/ contains the code for the Hermite polynomial interpolator.
 * common/ contains low-level code common to both the Lagrange and
   the Hermite interpolators.
@@ -63,9 +67,10 @@ to this interpolator, you need to
 * edit the makefile to create the appropriate directory or directories,
   and if necessary add any new targets
 * 'make' as appropriate to create the new coefficients etc
-  [note this can be rather computationally intensive -- the current
-  set of {1d,2d,3d} coefficients take {4 seconds, 15 seconds, 8 minutes}
-  to generate using Maple 7 on a xeon (1.7 GHz Pentium 4)]
+  [note this can be rather computationally intensive -- for example,
+  the set of {1d,2d,3d} Lagrange-tensor-product coefficients take about
+  {4 seconds, 40 seconds, 20 minutes} to generate using Maple 7 on
+  an AEI xeon (1.7 GHz Pentium 4 with 1GB memory)]
 * edit InterpLocalUniform.c to add the new case to the decoding
   switch statements
 * create an appropriate "script" file which defines the right macros,
diff --git a/src/GeneralizedPolynomial-Uniform/startup.c b/src/GeneralizedPolynomial-Uniform/startup.c
index 0c6be52..5b1ff01 100644
--- a/src/GeneralizedPolynomial-Uniform/startup.c
+++ b/src/GeneralizedPolynomial-Uniform/startup.c
@@ -33,16 +33,22 @@ void LocalInterp_GPU_Startup(void);
   @@*/
 void LocalInterp_GPU_Startup(void)
 {
-CCTK_InterpRegisterOpLocalUniform(LocalInterp_ILU_Lagrange,
-                                  "Lagrange polynomial interpolation",
+CCTK_InterpRegisterOpLocalUniform(LocalInterp_ILU_Lagrange_TP,
+				  "Lagrange polynomial interpolation (tensor product)",
 				  CCTK_THORNSTRING);
-
-/* old operator name for backwards compatability */
-CCTK_InterpRegisterOpLocalUniform(LocalInterp_ILU_Lagrange,
-                                  "generalized polynomial interpolation",
+CCTK_InterpRegisterOpLocalUniform(LocalInterp_ILU_Lagrange_MD,
+				  "Lagrange polynomial interpolation (maximum degree)",
 				  CCTK_THORNSTRING);
 
 CCTK_InterpRegisterOpLocalUniform(LocalInterp_ILU_Hermite,
                                   "Hermite polynomial interpolation",
 				  CCTK_THORNSTRING);
+
+/* synonym operator names for backwards compatability */
+CCTK_InterpRegisterOpLocalUniform(LocalInterp_ILU_Lagrange_TP,
+				  "Lagrange polynomial interpolation",
+				  CCTK_THORNSTRING);
+CCTK_InterpRegisterOpLocalUniform(LocalInterp_ILU_Lagrange_TP,
+                                  "generalized polynomial interpolation",
+				  CCTK_THORNSTRING);
 }
diff --git a/src/GeneralizedPolynomial-Uniform/template.c b/src/GeneralizedPolynomial-Uniform/template.c
index c12a308..454f733 100644
--- a/src/GeneralizedPolynomial-Uniform/template.c
+++ b/src/GeneralizedPolynomial-Uniform/template.c
@@ -863,6 +863,19 @@ int pt;
 	}
       }
 
+#ifdef PICKLE
+  #if (N_DIMS == 1)
+printf("pt=%d interp coord %.16g\n",
+	pt, (double) interp_coords_fp[0]);
+  #endif
+  #if (N_DIMS == 3)
+printf("pt=%d interp coords %.16g %.16g %.16g\n",
+	pt, (double) interp_coords_fp[0],
+	    (double) interp_coords_fp[1],
+	    (double) interp_coords_fp[2]);
+  #endif
+#endif
+
 
     /*
      * compute position of interpolation molecules with respect to
@@ -934,6 +947,17 @@ int pt;
     const fp  z        = xyz_temp       [Z_AXIS];
   #endif
 
+#ifdef PICKLE
+  #if (N_DIMS == 1)
+printf("interp center_i = %d\n", center_i);
+printf("interp grid-relative x = %.16g\n", (double) x);
+  #endif
+  #if (N_DIMS == 3)
+printf("interp center_ijk = %d %d %d\n", center_i, center_j, center_k);
+printf("interp grid-relative xyz = %.16g %.16g %.16g\n",
+	(double) x, (double) y, (double) z);
+  #endif
+#endif
 
     /*
      * compute 1-d position of molecule center in input arrays