initial checkin of new LocalInterp thorn

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

10 files changed, 3390 insertions, 0 deletions

diff --git a/doc/api1 b/doc/api1
new file mode 100644
index 0000000..5bc5f22
--- /dev/null
+++ b/doc/api1
@@ -0,0 +1,323 @@
+As I've explained in a recent message to the cactususers mailing list,
+I'm working on a new "generalized interpolator" for Cactus, which will
+be able to interpolate derivatives of grid functions as well as the
+grid functions themselves.  Corresponding to this, I (in consultation
+with Thomas Radke and Tom Goodale) have also come up with a new design
+for the Cactus interpolation API.
+
+The purpose of this message is to solicit comments and/or suggestions
+about the new API.
+
+
+Problems with the Current Interpolator API
+==========================================
+
+I see two major problems with the current interpolator API
+(CCTK_InterpGV() and CCTK_InterpLocal()):
+
+The current API uses variable argument lists to allow interpolating
+many arrays in a single call.  This has turned out to be somewhat
+error-prone.  (There is *no* compiler checking of variable argument lists!)
+To allow better compiler error checking, we propose replacing all the
+variable argument lists with explicit arrays of pointers.
+
+The current interpolator interface encodes all parameters for the
+interpolator (eg things like the order) into the character-string
+interpolator name.  This is a bit awkward, but managable if there are
+only a few parameters, all integers.  But it's completely unworkable for
+a generalized interpolator, which needs a whole bunch of extra parameters
+to tell it what derivatives to take.  Instead, we propose that a
+generalized interpolator get a handle to one of the new key-value
+tables, and that any extra parameters (such as "which derivatives to
+take" info) be passed in the table.  This makes it easy to add new
+interpolation operators which take additional parameters (eg a rational
+interpolator might have two order parameters, one for numerator and
+one for denominator), without having to change the API.
+
+
+The Generalized Interpolator API (version 1)
+============================================
+
+In more detail, here is my proposal for the new API:
+(don't be scared by the length, it's not that complicated!)
+(there are some examples below)
+
+  int status = CCTK_GInterpGV(arguments described below)
+  int status = CCTK_GInterpLocal(arguments described below)
+
+arguments (both CCTK_GInterpGV() and CCTK_GInterpLocal())
+  const cGH *GH;
+  int N_dims;
+  int operator_handle, int coord_system_handle;
+  int param_table_handle;	/* handle to key-value table, see below */
+
+  int N_interp_points;
+  /* array of CCTK_VARIABLE_* codes giving data types */
+  /* of arrays pointed to by  interp_coords[]  arrays */
+  const int interp_coord_types[N_dims];
+  /* array of pointers to arrays giving x,y,... coordinates of interp points */
+  const void *const interp_coords[N_dims];
+
+  int N_input_arrays, N_output_arrays;
+
+for CCTK_GInterpGV():
+  /* array of CCTK variable indices of input arrays */
+  const int input_array_indices[N_input_arrays];
+
+for CCTK_GInterpLocal():
+  /* array of CCTK_VARIABLE_* codes giving data types of input arrays */
+  const int input_array_types[N_input_arrays];
+  /* array of input array dimensions (common to all input arrays) */
+  const int input_array_dims[N_dims];
+  /* array of pointers to input arrays */
+  const void *const input_arrays[N_input_arrays];
+
+for both functions again:
+  /* array of CCTK_VARIABLE_* codes giving data types of output arrays */
+  const int output_array_types[N_output_arrays];
+  /* array of pointers to output arrays */
+  void *const output_arrays[N_output_arrays];
+
+The information passed in the table would depend on the interpolation
+operator.  For what I am currently implementing, I will register only
+a single operator, "generalized polynomial interpolation", and take the
+following parameters in the table:
+
+  int order;			/* mandatory */
+
+Thus the simplest call can just use Util_TableCreateFromString("order=3")
+for a 3rd-order interpolation, for example.
+
+All the remaining parameters in the table are optional; if they're
+omitted defaults will be supplied:
+
+  const char *molecule_type;	/* this selects one of a family of related */
+				/* operators; the default (and the only one */
+				/* I'm implementing right now) is to use the */
+				/* usual hypercube-shaped molecules */
+
+  int smoothing;		/* the way I'm implementing the interpolation */
+				/* it's easy to also do Savitzky-Golay type */
+				/* smoothing; this parameter gives how much */
+				/* to enlarge the interpolation molecule for */
+				/* this; the default is 0 (no smoothing) */
+
+Optionally, the caller can specify a mask gridfn
+for CCTK_GInterpGV():
+  int mask_gridfn_index;
+for CCTK_GInterpLocal();
+  int mask_type;		/* one of the CCTK_VARIABLE_* codes */
+  const void *const mask_array;
+and the range of values for the mask which correspond to valid points:
+  CCTK_REAL mask_min, mask_max;	/* valid <--> closed interval [min,max] */
+
+And last but not least, optionally, the caller can specify taking
+derivatives as part of the interpolation:
+  const int operand_indices[N_output_arrays];
+  const int opcodes[N_output_arrays];
+The semantics here are that
+   output array[i] = op(input array[ operand_indices[i] ])
+where  op  is specified as an integer operation code
+(we'll put #defines for these in one of the cctk header files).
+Note that the array operand_indices[] doesn't directly name the inputs,
+but rather gives the indices in the list of inputs.  This allows for a
+more efficient implementation in the case where some of the input arrays
+have many different operations applied to them.
+
+
+
+Pointers in Fortran
+===================
+
+One possible problem area with this API is that it requires creating
+arrays of pointers pointing to other arrays.  In C this is no problem,
+but in Fortran 77 this is difficult.  So, I propose adding a new Cactus
+function  pointer_to()  to make this easier for Fortran users:
+
+   CCTK_POINTER pointer_to(any array in your code)
+
+This would be #defined to %loc on those compilers which have that
+extension to standard Fortran, or would be a Cactus-provided utility
+routine for other cases.  It's trivial to write the latter case in C
+so long as the Fortran compiler actually uses call by reference, as
+almost all Fortran compilers do anyway for arrays.
+
+
+
+A Simple Example
+================
+
+Here's a simple example, written in Fortran 77, interpolating a real
+and a complex grid function in 3D:
+
+c
+c input grid functions:
+c    real_gridfn     (CCTK_REAL)
+c    complex_gridfn  (CCTK_COMPLEX)
+c
+c interpolation coordinates
+c    xcoord, ycoord, zcoord   (CCTK_REAL arrays of size N)
+c
+c output arrays:
+c    real_at_xyz     (CCTK_REAL array of size N)
+c    complex_at_xyz  (CCTK_COMPLEX array of size N)
+c
+	integer status
+	CCTK_POINTER input_gridfn_indices(2)
+	integer interp_coord_types(3)
+	CCTK_POINTER interp_coords(3)
+	integer output_array_types(2)
+	CCTK_POINTER output_arrays(2)
+
+	call CCTK_VarIndex(input_gridfn_indices(1), "mythorn::real_gridfn")
+	call CCTK_VarIndex(input_gridfn_indices(2), "mythorn::complex_gridfn")
+
+c could also initialize this with a DATA statement
+	interp_coord_types(1) = CCTK_VARIABLE_REAL
+	interp_coord_types(2) = CCTK_VARIABLE_REAL
+	interp_coord_types(3) = CCTK_VARIABLE_REAL
+	interp_coords(1) = pointer_to(xcoord)
+	interp_coords(2) = pointer_to(ycoord)
+	interp_coords(3) = pointer_to(zcoord)
+c could also initialize this with a DATA statement
+	output_array_types(1) = CCTK_VARIABLE_REAL
+	output_array_types(2) = CCTK_VARIABLE_COMPLEX
+	output_arrays(1) = pointer_to(real_at_xyz)
+	output_arrays(2) = pointer_to(complex_at_xyz)
+	call CCTK_InterpGV(status,
+			   cctk_GH,
+			   3,			# number of dimensions
+			   operator_handle, coord_system_handle,
+			   Util_TableCreateFromString("order=3"),
+			   N, interp_coord_types, interp_coords,
+			   2, 2,		# number of input/output arrays
+			   input_gridfn_indices,
+			   output_array_types, output_arrays);
+	if (status .lt. 0) then
+c		help, an error occured!
+	end if
+
+
+
+A More Complicated Example
+==========================
+
+Here's a more complicated example, written in C++.  (I'm really only using
+C++ to get cleaner initialization of the various arrays, this is still
+"almost C".)  This example is a simplified form of what I will be doing
+in my new apparent horizon finder:
+
+//
+// input grid functions (12 of them, all 3D CCTK_REAL):
+//	gxx, gxy, gxz, gyy, gyz, gzz,
+//	Kxx, Kxy, Kxz, Kyy, Kyz, Kzz
+//
+// interpolation coordinates:
+//	xcoord, ycoord, zcoord   (CCTK_REAL arrays of size N)
+//
+// output arrays (30 of them, all 3D CCTK_REAL)
+// (interpolate the gij and Kij, also interpolate all first derivs of the gij)
+//	I_gxx, dx_gxx, dy_gxx, dz_gxx,
+//	I_gxy, dx_gxy, dy_gxy, dz_gxy,
+//	I_gxz, dx_gxz, dy_gxz, dz_gxz,
+//	I_gyy, dx_gyy, dy_gyy, dz_gyy,
+//	I_gyz, dx_gyz, dy_gyz, dz_gyz,
+//	I_gzz, dx_gzz, dy_gzz, dz_gzz,
+//	I_Kxx, I_Kxy, I_Kxz, I_Kyy, I_Kyz, I_Kzz,
+//
+
+const int gxx_index = CCTK_VarIndex("somethorn::gxx");
+const int gxy_index = CCTK_VarIndex("somethorn::gxy");
+const int gxz_index = CCTK_VarIndex("somethorn::gxz");
+const int gyy_index = CCTK_VarIndex("somethorn::gyy");
+const int gyz_index = CCTK_VarIndex("somethorn::gyz");
+const int gzz_index = CCTK_VarIndex("somethorn::gzz");
+const int Kxx_index = CCTK_VarIndex("somethorn::Kxx");
+const int Kxy_index = CCTK_VarIndex("somethorn::Kxy");
+const int Kxz_index = CCTK_VarIndex("somethorn::Kxz");
+const int Kyy_index = CCTK_VarIndex("somethorn::Kyy");
+const int Kyz_index = CCTK_VarIndex("somethorn::Kyz");
+const int Kzz_index = CCTK_VarIndex("somethorn::Kzz");
+
+const int N_dims = 3;
+int interp_coord_types[N_dims]
+	= { CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL };
+const void *const interp_coords[N_dims]
+	= { (const void *) &xcoord,
+	    (const void *) &ycoord,
+	    (const void *) &zcoord };
+
+const int N_input_arrays = 12;
+int input_array_types[N_input_arrays]
+	= { CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+	    CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+	    CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+	    CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL };
+
+const int N_output_arrays = 30;
+int output_array_types[N_output_arrays];
+	for (int oi = 0 ; oi < N_output_arrays ; ++oi)
+	{
+	output_array_types[oi] = CCTK_VARIABLE_REAL;
+	}
+
+#define VP(x) static_cast<void *>(x)
+void *const output_arrays[N_output_arrays]
+	= {
+	  VP(I_gxx), VP(dx_gxx), VP(dy_gxx), VP(dz_gxx),
+	  VP(I_gxy), VP(dx_gxy), VP(dy_gxy), VP(dz_gxy),
+	  VP(I_gxz), VP(dx_gxz), VP(dy_gxz), VP(dz_gxz),
+	  VP(I_gyy), VP(dx_gyy), VP(dy_gyy), VP(dz_gyy),
+	  VP(I_gyz), VP(dx_gyz), VP(dy_gyz), VP(dz_gyz),
+	  VP(I_gzz), VP(dx_gzz), VP(dy_gzz), VP(dz_gzz),
+	  VP(I_Kxx), VP(I_Kxy), VP(I_Kxz), VP(I_Kyy), VP(I_Kyz), VP(I_Kzz),
+	  };
+
+int operand_indices[N_output_arrays];
+	= {
+	  gxx_index, gxx_index, gxx_index, gxx_index,
+	  gxy_index, gxy_index, gxy_index, gxy_index,
+	  gxz_index, gxz_index, gxz_index, gxz_index,
+	  gyy_index, gyy_index, gyy_index, gyy_index,
+	  gyz_index, gyz_index, gyz_index, gyz_index,
+	  gzz_index, gzz_index, gzz_index, gzz_index,
+	  Kxx_index, Kxy_index, Kxz_index, Kyy_index, Kyz_index, Kzz_index,
+	  };
+
+const int op_I  = CCTK_INTERP_OPCODE_I;
+const int op_dx = CCTK_INTERP_OPCODE_DX;
+const int op_dy = CCTK_INTERP_OPCODE_DY;
+const int op_dz = CCTK_INTERP_OPCODE_DZ;
+int opcodes[N_output_arrays]
+	= {
+	  op_I, op_dx, op_dy, op_dz,
+	  op_I, op_dx, op_dy, op_dz,
+	  op_I, op_dx, op_dy, op_dz,
+	  op_I, op_dx, op_dy, op_dz,
+	  op_I, op_dx, op_dy, op_dz,
+	  op_I, op_dx, op_dy, op_dz,
+	  op_I, op_I, op_I, op_I, op_I, op_I,
+	  };
+
+int param_table_handle = Util_TableCreate(UTIL_TABLE_DEFAULT);
+Util_TableSet1Int(param_table_handle,
+		  3, "order");
+Util_TableSetInt(param_table_handle,
+		 N_output_arrays, operand_indices, "operand_indices");
+Util_TableSetInt(param_table_handle,
+		 N_output_arrays, opcodes, "opcodes");
+
+int status = CCTK_GInterpGV(GH,
+			    N_dims,
+			    operator_handle, coord_system_handle,
+			    param_table_handle,
+			    N_interp_points, interp_coord_types, interp_coords,
+			    N_input_arrays, N_output_arrays,
+			    input_array_indices,
+			    output_array_types, output_arrays);
+if (status < 0)
+	{ /* something bad happened */ }
+Util_TableDestroy(param_table_handle);
+
+
+
diff --git a/doc/api2 b/doc/api2
new file mode 100644
index 0000000..422be0d
--- /dev/null
+++ b/doc/api2
@@ -0,0 +1,460 @@
+This is version 2 of my proposal for the new interpolator API.  The
+main changes since version 1 are
+* function names have changed (longer, more descriptive)
+* output arrays may now be N-dimensional (before they were only 1-D)
+* a bunch more optional stuff in the table to specify array strides etc
+
+Don't be scared by the length, for most uses it's not that complicated!
+There are some examples below...
+
+
+
+Synopsis
+========
+
+  int status = CCTK_InterpLocalArrays(arguments described below)
+  int status = CCTK_InterpGridArrays(arguments described below)
+
+return is 0 for ok, various -ve values for error codes
+
+
+
+Function Arguments
+==================
+
+for both CCTK_InterpLocalArrays() and CCTK_InterpGridArrays():
+  /***** misc arguments *****/
+  const cGH *GH;
+  /* note N_dims is the number of dimensions in the *interpolation*; */
+  /* this may be smaller than the number of dimensions of the input arrays */
+  /* if the storage indexing is set up appropriately (eg to interpolate */
+  /* in 1-D lines or 2-D planes of 3-D grid arrays) */
+  int N_dims;
+  int operator_handle, int coord_system_handle;
+  int param_table_handle;       /* handle to "parameter table", a key-value */
+                                /* table, see below for table entries */
+
+  /***** arguments specifying the interpolation points *****/
+  int N_interp_points;
+  /* array of CCTK_VARIABLE_* codes giving the data types */
+  /* of the arrays pointed to by  interp_coords[] */
+  const int interp_coord_type_codes[N_dims];
+  /* array[N_dims] of pointers to arrays[N_interp_points] */
+  /* giving x,y,z,... coordinates of interpolation points */
+  const void *const interp_coords[N_dims];
+
+  /***** arguments specifying the inputs (the data to be interpolated) *****/
+  int N_input_arrays;
+for CCTK_InterpLocalArrays():
+  /* array of CCTK_VARIABLE_* codes giving data types of input arrays */
+  const int input_array_type_codes[N_input_arrays];
+  /* array of input array dimensions (common to all input arrays) */
+  const int input_array_dims[N_dims];
+  /* array of pointers to input arrays */
+  const void *const input_arrays[N_input_arrays];
+for CCTK_InterpGridArrays():
+  /* array of CCTK variable indices of input arrays */
+  const int input_array_variable_indices[N_input_arrays];
+
+for both CCTK_InterpLocalArrays() and CCTK_InterpGridArrays() again:
+  /***** arguments specifying the outputs (the interpolation results) *****/
+  int N_output_arrays;
+  /* array of CCTK_VARIABLE_* codes giving data types of output arrays */
+  const int output_array_type_codes[N_output_arrays];
+  /* array[N_output_arrays] of pointers to output arrays[N_interp_points] */
+  void *const output_arrays[N_output_arrays];
+
+
+
+Information Passed in the Parameter Table
+=========================================
+
+The "parameter table" may be used to specify non-default storage indexing
+for input or output arrays, and/or various options for the interpolation
+itself.  Some interpolators may not implement all of these options.
+
+
+Storage Indexing Options
+------------------------
+
+Sometimes one of the "arrays" used by the interpolator isn't contiguous
+in memory.  For example, we might want to do 2-D interpolation within a
+plane of a 3-D grid array, and/or the grid array might be a member of a
+compact group.  To support this, we use several optional table entries
+(these should be supported by all interpolation operators):
+
+For the input arrays, we use
+
+  const int input_array_offsets[N_input_arrays];
+  /* next 3 table entries are shared by all input arrays */
+  const int input_array_strides       [N_dims];
+  const int input_array_min_subscripts[N_dims];
+  const int input_array_max_subscripts[N_dims];
+
+Then for input array number a, the generic subscripting expression for
+the 3-D case is
+  data_pointer[offset + i*istride + j*jstride + k*kstride]
+where
+  data_pointer = input_arrays[a]
+  offset = input_array_offsets[a]
+  (istride,jstride,kstride) = input_array_stride[]
+and where (i,j,k) run from input_array_min_subscripts[] to
+input_array_max_subscripts[] inclusive.
+
+The defaults are offset=0, stride=determined from input_array_dims[]
+in the usual Fortran manner, input_array_min_subscripts[] = 0,
+input_array_max_subscripts[] = input_array_dims[]-1.  If the stride
+and max subscript are both specified explicitly, then the
+input_array_dims[] function argument is ignored.
+
+For CCTK_InterpGridArrays() operating on a member of a compact group
+the offset and strides are interpreted in units of _grid_points_.  This
+has the advantage that interpolator calls need not be changed if a grid
+array is changed from being simple to/from compact.  In terms of actual
+memory addressing, then, the internal subscripting expression for this
+case would be
+  group_data_pointer[offset + member_number + i*istride*N_members
+                                            + j*jstride*N_members
+                                            + k*kstride*N_members]
+
+By default the interpolation-point coordinates and the output arrays
+are all contiguous 1-D arrays.  This may be changed with the table
+entries
+
+  const int interp_coords_offsets[N_dims];
+  const int output_array_offsets[N_output_arrays];
+  /* next 4 table entries are shared by all interp coords and output arrays */
+  const int interp_point_N_dims;
+  const int interp_point_strides       [interp_point_N_dims];
+  const int interp_point_min_subscripts[interp_point_N_dims];
+  const int interp_point_max_subscripts[interp_point_N_dims];
+
+For example, if we wanted to do 3-D interpolation, interpolating a value
+at each non-ghost-zone point of a 2-D grid of points, with the grid point
+coordinates stored as 2-D arrays, we would use
+  N_dims = 3
+  interp_point_N_dims = 2
+  interp_point_strides[] = set up from the full size of the 2-D grid
+  interp_point_{min,max}_subscripts[] = specify the non-ghost-zone points
+                                        of the 2-D grid
+
+Excision Options
+----------------
+
+Some interpolators may specifically support excision, where a mask array
+(same dimensionality and indexing as the input arrays) is used to mark
+some grid points as valid (ok to use data there) and other grid points
+as invalid (the interpolator isn't allowed to use data there).
+
+If an interpolator supports this, it should use the following optional
+parameters:
+
+for CCTK_InterpLocalArrays();
+  const int mask_type_code;     /* one of the CCTK_VARIABLE_* codes */
+  const void *const mask_array;
+for CCTK_InterpGridArrays():
+  const int mask_variable_index;
+
+for both CCTK_InterpLocalArrays() and CCTK_InterpGridArrays():
+  /* we consider a grid point to be valid if and only if the mask */
+  /* has a value in the closed interval [mask_valid_min,mask_valid_max] */
+  /* n.b. the caller should beware of possible rounding problems here; */
+  /*      it may be appropriate to widen the valid interval slightly */
+  /*      if the endpoints aren't exactly-representable floating-point */
+  /*      values */
+  const mask_type mask_valid_min, mask_valid_max;
+
+The same type of storage options supported for the input and/or output
+arrays, are also supported for the mask; the mask may have its own offset,
+but shares any input-array stride and/or min/max subscript specification:
+
+  const int mask_offset;
+
+
+The remaining parameter-table options are specific to the new interpolator
+I'm currently implementing for PUGHInterp.  This registers (only) a single
+operator, "generalized polynomial interpolation".
+
+
+Interpolation Order and Molecule Family
+---------------------------------------
+
+The mandatory parameter
+
+  const int order;
+
+sets the order of the interpolating polynomial (1=linear, 2=quadratic,
+3=cubic, etc).  Thus the simplest call can just use (eg)
+  Util_TableCreateFromString("order=3")
+for cubic interpolation.
+
+All the remaining parameters in the table are optional; if they're
+omitted defaults will be supplied.
+
+  /* this selects one of a family of related operators */
+  /* the default (and the only one I'm implementing right now) */
+  /* is "cube" to use the usual hypercube-shaped molecules */
+  const char *const molecule_family;
+
+Smoothing
+---------
+
+The way I'm implementing the interpolation it's easy to also do
+Savitzky-Golay type smoothing (= moving least-square fitting, cf
+Numerical Recipes 2nd edition section 14.8).  This is controlled by
+the parameter
+
+  const int smoothing;
+
+which says how much (how many points) to enlarge the interpolation
+molecule for this.  The default is 0 (no smoothing).  1 would mean to
+enlarge the molecule by 1 point, e.g. to use a 5-point molecule instead
+of the usual 4-point one for cubic interpolation.  2 would mean to
+enlarge by 2 points, e.g. to use a 6-point molecule for cubic
+interpolation.  Etc etc.
+
+This type of smoothing is basically free apart from the increase in
+the molecule size, e.g. a smoothing=2 cubic interpolation has exactly
+the same cost as any other 6-point-molecule interpolation.
+
+Derivatives
+-----------
+
+This interpolator can optionally (and again at no extra cost) take
+partial derivatives as part of the interpolation:
+  const int operand_indices[N_output_arrays];
+  const int opcodes        [N_output_arrays];
+The semantics here are that
+   output array[i] = op(input array[ operand_indices[i] ])
+where  op  is specified as an integer operation code as described below.
+
+Note that the array operand_indices[] doesn't directly name the inputs,
+but rather gives the indices (0-origin) in the list of inputs.  This
+allows for a more efficient implementation in the case where some of
+the input arrays have many different operations applied to them.
+
+The operations are coded based on the decimal digits of the integer:
+each decimal digit means to take the derivative in that direction;
+the order of the digits in a number is ignored.  For example:
+  0 = no derivative, "just" interpolate
+  1 = interpolate d/dx1 (derivative along x1 coordinate)
+  2 = interpolate d/dx2 (derivative along x2 coordinate)
+  11 = interpolate d^2/dx1^2 (2nd derivative along x1 coordinate)
+  22 = interpolate d^2/dx2^2 (2nd derivative along x2 coordinate)
+  12 = 21 = interpolate d^2/dx1 dx2 (mixed 2nd partial derivative in x1 and x2)
+  122 = 212 = 221 = interpolate d^3/dx1 dx2^2 (mixed 3rd partial derivative)
+  222 = interpolate d^3/dx2^3 (3rd derivative along x2 coordinate)
+  123 = 132 = 213 = 231 = 312 = 321
+      = interpolate d^3/dx1 dx2 dx3 (mixed 3rd partial derivative)
+
+
+
+Pointers in Fortran
+===================
+
+One possible problem area with this API is that it requires creating
+arrays of pointers pointing to other arrays.  In C this is no problem,
+but in Fortran 77 this is difficult.  So, I propose adding two new Cactus
+functions to make this easier for Fortran users:
+
+   CCTK_POINTER Util_PointerTo(any Fortran variable or array)
+   CCTK_POINTER Util_NullPointer()
+
+Util_PointerTo would be #defined to %loc on those compilers which have
+that extension to standard Fortran, or would be a Cactus-provided utility
+routine for other cases.  It's trivial to write the latter case in C so
+long as the Fortran compiler actually uses call by reference; I've never
+heard of a Fortran compiler doing otherwise for arrays.  (And even for
+Fortran scalar variables it would be very hard for a compiler to do otherwise
+in light of separate compilation and 1-element arrays being allowed to be
+passed to/from scalar variables.)
+
+
+
+A Simple Example
+================
+
+Here's a simple example, written in Fortran 77, to do quadratic interpolation
+of a real and a complex local array in 3-D:
+
+c input arrays:
+        integer ni, nj, nk
+        parameter (ni=..., nj=..., nk=...)
+        CCTK_REAL    real_gridfn   (ni,nj,nk)
+        CCTK_COMPLEX complex_gridfn(ni,nj,nk)
+
+c interpolation coordinates
+        integer N_interp
+        parameter (N_interp = ...)
+        CCTK_REAL xcoord(N_interp), y_coord(N_interp), z_coord(N_interp)
+
+c output arrays:
+        CCTK_REAL    real_at_xyz   (N_interp)
+        CCTK_COMPLEX complex_at_xyz(N_interp)
+
+        integer status, dummy
+        integer input_array_type_codes(2)
+        data input_array_type_codes /CCTK_VARIABLE_REAL,
+     $                               CCTK_VARIABLE_COMPLEX/
+        integer input_array_dims(3)
+        CCTK_POINTER input_arrays(2)
+        integer interp_coord_type_codes(3)
+        data interp_coord_type_codes /CCTK_VARIABLE_REAL,
+     $                                CCTK_VARIABLE_REAL,
+     $                                CCTK_VARIABLE_REAL/
+        CCTK_POINTER interp_coords(3)
+        integer output_array_type_codes(2)
+        data output_array_type_codes /CCTK_VARIABLE_REAL,
+     $                                CCTK_VARIABLE_COMPLEX/
+        CCTK_POINTER output_arrays(2)
+
+        input_array_dims(1) = ni
+        input_array_dims(2) = nj
+        input_array_dims(3) = nk
+        interp_coords(1) = Util_PointerTo(xcoord)
+        interp_coords(2) = Util_PointerTo(ycoord)
+        interp_coords(3) = Util_PointerTo(zcoord)
+        output_arrays(1) = Util_PointerTo(real_at_xyz)
+        output_arrays(2) = Util_PointerTo(complex_at_xyz)
+
+        call CCTK_InterpLocalArrays
+     $          (status,                ! return code
+                 cctk_GH,
+                 3,                     ! number of dimensions
+                 operator_handle, coord_system_handle,
+                 Util_TableCreateFromString("order=2"),
+                 N_interp,
+                 interp_coord_type_codes, interp_coords,
+                 2,                     ! number of input arrays
+                 input_array_type_codes, input_array_dims, input_arrays,
+                 2,                     ! number of output arrays
+                 output_array_type_codes, output_arrays)
+
+        if (status .lt. 0) then
+                call CCTK_WARN(status, "Error return from interpolator!")
+                call CCTK_Exit(dummy, GH, status)
+        end if
+
+
+
+A More Complicated Example
+==========================
+
+Here's a more complicated example, written in C++.  (I'm really only using
+C++ to get cleaner initialization of the various arrays, this is still
+"almost C".)  This example is a simplified form of what I will be doing
+in my new apparent horizon finder:
+
+//
+// input grid functions (12 of them, all 3-D CCTK_REAL):
+//      gxx, gxy, gxz, gyy, gyz, gzz,
+//      Kxx, Kxy, Kxz, Kyy, Kyz, Kzz
+//
+// interpolation coordinates:
+//      xcoord, ycoord, zcoord   (all CCTK_REAL[N_interp_points])
+//
+// we want to interpolate the gij and Kij, and also interpolate all the
+// first derivatives of the gij, so the output arrays are
+// (30 of them, all CCTK_REAL[N_interp_points])
+//      I_gxx, dx_gxx, dy_gxx, dz_gxx,
+//      I_gxy, dx_gxy, dy_gxy, dz_gxy,
+//      I_gxz, dx_gxz, dy_gxz, dz_gxz,
+//      I_gyy, dx_gyy, dy_gyy, dz_gyy,
+//      I_gyz, dx_gyz, dy_gyz, dz_gyz,
+//      I_gzz, dx_gzz, dy_gzz, dz_gzz,
+//      I_Kxx, I_Kxy, I_Kxz, I_Kyy, I_Kyz, I_Kzz
+//
+
+#define VP(x) static_cast<void *>(x)
+
+const int N_dims = 3;
+const int interp_coord_type_codes[N_dims]
+        = { CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL };
+const void *const interp_coords[N_dims]
+        = { VP(xcoord), VP(ycoord), VP(zcoord) };
+
+const int N_input_arrays = 12;
+const int input_array_types[N_input_arrays]
+        = { CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+            CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+            CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+            CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL };
+
+const int input_array_variable_indices[N_input_arrays]
+        = { CCTK_VarIndex("somethorn::gxx"),
+            CCTK_VarIndex("somethorn::gxy"),
+            CCTK_VarIndex("somethorn::gxz"),
+            CCTK_VarIndex("somethorn::gyy"),
+            CCTK_VarIndex("somethorn::gyz"),
+            CCTK_VarIndex("somethorn::gzz"),
+            CCTK_VarIndex("somethorn::Kxx"),
+            CCTK_VarIndex("somethorn::Kxy"),
+            CCTK_VarIndex("somethorn::Kxz"),
+            CCTK_VarIndex("somethorn::Kyy"),
+            CCTK_VarIndex("somethorn::Kyz"),
+            CCTK_VarIndex("somethorn::Kzz") };
+
+const int N_output_arrays = 30;
+int output_array_type_codes[N_output_arrays];
+        for (int oi = 0 ; oi < N_output_arrays ; ++oi)
+        {
+        output_array_type_codes[oi] = CCTK_VARIABLE_REAL;
+        }
+
+void *const output_arrays[N_output_arrays]
+        = {
+          VP(I_gxx), VP(dx_gxx), VP(dy_gxx), VP(dz_gxx),
+          VP(I_gxy), VP(dx_gxy), VP(dy_gxy), VP(dz_gxy),
+          VP(I_gxz), VP(dx_gxz), VP(dy_gxz), VP(dz_gxz),
+          VP(I_gyy), VP(dx_gyy), VP(dy_gyy), VP(dz_gyy),
+          VP(I_gyz), VP(dx_gyz), VP(dy_gyz), VP(dz_gyz),
+          VP(I_gzz), VP(dx_gzz), VP(dy_gzz), VP(dz_gzz),
+          VP(I_Kxx), VP(I_Kxy), VP(I_Kxz), VP(I_Kyy), VP(I_Kyz), VP(I_Kzz)
+          };
+
+const int operand_indices[N_output_arrays];
+        = {
+          0, 0, 0, 0,           // gxx
+          1, 1, 1, 1,           // gxy
+          2, 2, 2, 2,           // gxz
+          3, 3, 3, 3,           // gyy
+          4, 4, 4, 4,           // gyz
+          5, 5, 5, 5,           // gzz
+          6, 7, 8, 9, 10, 11    // Kxx-Kzz
+          };
+
+const int opcodes[N_output_arrays]
+        = {
+          0, 1, 2, 3,           // I, dx, dy, dz
+          0, 1, 2, 3,           // I, dx, dy, dz
+          0, 1, 2, 3,           // I, dx, dy, dz
+          0, 1, 2, 3,           // I, dx, dy, dz
+          0, 1, 2, 3,           // I, dx, dy, dz
+          0, 1, 2, 3,           // I, dx, dy, dz
+          0, 0, 0, 0, 0, 0      // all I
+          };
+
+int param_table_handle = Util_TableCreate(UTIL_TABLE_DEFAULT);
+Util_TableSet1Int(param_table_handle,
+                  3, "order");
+Util_TableSetInt(param_table_handle,
+                 N_output_arrays, operand_indices, "operand_indices");
+Util_TableSetInt(param_table_handle,
+                 N_output_arrays, opcodes, "opcodes");
+
+int status = CCTK_InterpGridArrays(GH,
+                                   N_dims,
+                                   operator_handle, coord_system_handle,
+                                   param_table_handle,
+                                   N_interp_points,
+                                   interp_coord_type_codes, interp_coords,
+                                   N_input_arrays,
+                                   input_array_variable_indices,
+                                   N_output_arrays,
+                                   output_array_type_codes, output_arrays);
+if (status < 0)
+        {
+        CCTK_WARN(status, "error return from CCTK_InterpGridArrays()!");
+        CCTK_Exit(GH, status);                                  /*NOTREACHED*/
+        }
+Util_TableDestroy(param_table_handle);
diff --git a/doc/api2.1 b/doc/api2.1
new file mode 100644
index 0000000..e316aba
--- /dev/null
+++ b/doc/api2.1
@@ -0,0 +1,488 @@
+Author: Jonathan Thornburg <jthorn@aei.mpg.de>
+Date: December 2001
+
+This is version 2.1 of my proposal for the new interpolator API.  The
+main changes since version 1 are
+* function names have changed (longer, more descriptive)
+* output arrays may now be N-dimensional (before they were only 1-D)
+* a bunch more optional stuff in the table to specify array strides etc
+* no #defines in the header file for derivative operation codes
+* int arrays changed to CCTK_INT
+
+Don't be scared by the length, for most uses it's not that complicated!
+There are some examples below...
+
+
+
+Synopsis
+========
+
+  int status = CCTK_InterpLocalArrays(arguments described below)
+  int status = CCTK_InterpGridArrays(arguments described below)
+
+return is 0 for ok, various -ve values for error codes
+
+(N.b. the flesh APIs to register interpolation operators will also need
+their C function prototypes changed to match the changes here.)
+
+
+
+Function Arguments
+==================
+
+  /***** misc arguments *****/
+for CCTK_InterpGridArrays() only:
+  const cGH *GH;
+for both CCTK_InterpLocalArrays() and CCTK_InterpGridArrays():
+  /* note N_dims is the number of dimensions in the *interpolation*; */
+  /* this may be smaller than the number of dimensions of the input arrays */
+  /* if the storage indexing is set up appropriately (eg to interpolate */
+  /* in 1-D lines or 2-D planes of 3-D grid arrays) */
+  int N_dims;
+  int operator_handle;
+  int param_table_handle;       /* handle to "parameter table", a key-value */
+                                /* table, see below for table entries */
+
+  /***** arguments specifying the coordinate system *****/
+for CCTK_InterpLocalArrays():
+  /* the local coordinate system is specified as a generic linear mapping */
+  /* from (integer) local input array subscripts --> (global) coordinates: */
+  /*   coordinate = coord_system_origin[axis] + subscript*grid_spacing[axis] */
+  const CCTK_REAL coord_system_origin[N_dims];	/* coords of subscript 0 */
+  const CCTK_REAL grid_spacing[N_dims];
+for CCTK_InterpGridArrays():
+  int coord_system_handle;	/* handle to Cactus coordinate system */
+				/* specifying mapping of */
+				/* (integer) input array subscripts */
+				/* <--> (floating point) coordinates */
+
+  /***** arguments specifying the interpolation points *****/
+  int N_interp_points;
+  /* (pointer to) array of CCTK_VARIABLE_* codes giving the */
+  /* data types of the arrays pointed to by  interp_coords[] */
+  const CCTK_INT interp_coord_type_codes[N_dims];
+  /* (pointer to) array[N_dims] of pointers to arrays[N_interp_points] */
+  /* giving x,y,z,... coordinates of interpolation points */
+  const void *const interp_coords[N_dims];
+
+  /***** arguments specifying the inputs (the data to be interpolated) *****/
+  int N_input_arrays;
+for CCTK_InterpLocalArrays():
+  /* array of input array dimensions (common to all input arrays) */
+  const CCTK_INT input_array_dims[N_dims];
+  /* array of CCTK_VARIABLE_* codes giving data types of input arrays */
+  const CCTK_INT input_array_type_codes[N_input_arrays];
+  /* array of pointers to input arrays */
+  const void *const input_arrays[N_input_arrays];
+for CCTK_InterpGridArrays():
+  /* array of CCTK variable indices of input arrays */
+  const CCTK_INT input_array_variable_indices[N_input_arrays];
+
+for both CCTK_InterpLocalArrays() and CCTK_InterpGridArrays() again:
+  /***** arguments specifying the outputs (the interpolation results) *****/
+  int N_output_arrays;
+  /* array of CCTK_VARIABLE_* codes giving data types of output arrays */
+  const CCTK_INT output_array_type_codes[N_output_arrays];
+  /* array[N_output_arrays] of pointers to output arrays[N_interp_points] */
+  void *const output_arrays[N_output_arrays];
+
+
+
+Information Passed in the Parameter Table
+=========================================
+
+The "parameter table" may be used to specify non-default storage indexing
+for input or output arrays, and/or various options for the interpolation
+itself.  Some interpolators may not implement all of these options.
+
+
+Subscripting Options
+--------------------
+
+Sometimes one of the "arrays" used by the interpolator isn't contiguous
+in memory.  For example, we might want to do 2-D interpolation within a
+plane of a 3-D grid array, and/or the grid array might be a member of a
+compact group.  To support this, we use several optional table entries
+(these should be supported by all interpolation operators):
+
+For the input arrays, we use
+
+  const CCTK_INT input_array_offsets[N_input_arrays];
+  /* next 3 table entries are shared by all input arrays */
+  const CCTK_INT input_array_strides       [N_dims];
+  const CCTK_INT input_array_min_subscripts[N_dims];
+  const CCTK_INT input_array_max_subscripts[N_dims];
+
+Then for input array number a, the generic subscripting expression for
+the 3-D case is
+  data_pointer[offset + i*istride + j*jstride + k*kstride]
+where
+  data_pointer = input_arrays[a]
+  offset = input_array_offsets[a]
+  (istride,jstride,kstride) = input_array_stride[]
+and where (i,j,k) run from input_array_min_subscripts[] to
+input_array_max_subscripts[] inclusive.
+
+The defaults are offset=0, stride=determined from input_array_dims[]
+in the usual Fortran manner, input_array_min_subscripts[] = 0,
+input_array_max_subscripts[] = input_array_dims[]-1.  If the stride
+and max subscript are both specified explicitly, then the
+input_array_dims[] function argument is ignored.
+
+For CCTK_InterpGridArrays() operating on a member of a compact group
+the offset and strides are interpreted in units of _grid_points_.  This
+has the advantage that interpolator calls need not be changed if a grid
+array is changed from being simple to/from compact.  In terms of actual
+memory addressing, then, the internal subscripting expression for this
+case would be
+  group_data_pointer[offset + member_number + i*istride*N_members
+                                            + j*jstride*N_members
+                                            + k*kstride*N_members]
+
+By default the interpolation-point coordinates and the output arrays
+are all contiguous 1-D arrays.  This may be changed with the table
+entries
+
+  const CCTK_INT interp_coords_offsets[N_dims];
+  const CCTK_INT output_array_offsets[N_output_arrays];
+  /* next 4 table entries are shared by all interp coords and output arrays */
+  const CCTK_INT interp_point_N_dims;
+  const CCTK_INT interp_point_strides       [interp_point_N_dims];
+  const CCTK_INT interp_point_min_subscripts[interp_point_N_dims];
+  const CCTK_INT interp_point_max_subscripts[interp_point_N_dims];
+
+For example, if we wanted to do 3-D interpolation, interpolating a value
+at each non-ghost-zone point of a 2-D grid of points, with the grid point
+coordinates stored as 2-D arrays, we would use
+  N_dims = 3
+  interp_point_N_dims = 2
+  interp_point_strides[] = set up from the full size of the 2-D grid
+  interp_point_{min,max}_subscripts[] = specify the non-ghost-zone points
+                                        of the 2-D grid
+
+Excision Options
+----------------
+
+Some interpolators may specifically support excision, where a mask array
+(same dimensionality and indexing as the input arrays) is used to mark
+some grid points as valid (ok to use data there) and other grid points
+as invalid (the interpolator isn't allowed to use data there).
+
+If an interpolator supports this, it should use the following optional
+parameters:
+
+for CCTK_InterpLocalArrays();
+  const CCTK_INT mask_type_code;  /* one of the CCTK_VARIABLE_* codes */
+  const void *const mask_array;   /* same dimensions/indexing as input arrays */
+for CCTK_InterpGridArrays():
+  const CCTK_INT mask_variable_index;
+
+for both CCTK_InterpLocalArrays() and CCTK_InterpGridArrays():
+  /* we consider a grid point to be valid if and only if the mask */
+  /* has a value in the closed interval [mask_valid_min,mask_valid_max] */
+  /* n.b. the caller should beware of possible rounding problems here; */
+  /*      it may be appropriate to widen the valid interval slightly */
+  /*      if the endpoints aren't exactly-representable floating-point */
+  /*      values */
+  const mask_type mask_valid_min, mask_valid_max;
+
+The same type of storage options supported for the input and/or output
+arrays, are also supported for the mask; the mask may have its own offset,
+but shares any input-array stride and/or min/max subscript specification:
+
+  const CCTK_INT mask_offset;
+
+
+The remaining parameter-table options are specific to the new interpolator
+I'm currently implementing for PUGHInterp.  This registers (only) a single
+operator, "generalized polynomial interpolation".
+
+
+Interpolation Order and Molecule Family
+---------------------------------------
+
+The mandatory parameter
+
+  const CCTK_INT order;
+
+sets the order of the interpolating polynomial (1=linear, 2=quadratic,
+3=cubic, etc).  Thus the simplest call can just use (eg)
+  Util_TableCreateFromString("order=3")
+for cubic interpolation.
+
+All the remaining parameters in the table are optional; if they're
+omitted defaults will be supplied.
+
+  /* this selects one of a family of related operators */
+  /* the default (and the only one I'm implementing right now) */
+  /* is "cube" to use the usual hypercube-shaped molecules */
+  const char *const molecule_family;
+
+Smoothing
+---------
+
+The way I'm implementing the interpolation it's easy to also do
+Savitzky-Golay type smoothing (= moving least-square fitting, cf
+Numerical Recipes 2nd edition section 14.8).  This is controlled by
+the parameter
+
+  const CCTK_INT smoothing;
+
+which says how much (how many points) to enlarge the interpolation
+molecule for this.  The default is 0 (no smoothing).  1 would mean to
+enlarge the molecule by 1 point, e.g. to use a 5-point molecule instead
+of the usual 4-point one for cubic interpolation.  2 would mean to
+enlarge by 2 points, e.g. to use a 6-point molecule for cubic
+interpolation.  Etc etc.
+
+This type of smoothing is basically free apart from the increase in
+the molecule size, e.g. a smoothing=2 cubic interpolation has exactly
+the same cost as any other 6-point-molecule interpolation.
+
+Derivatives
+-----------
+
+This interpolator can optionally (and again at no extra cost) take
+partial derivatives as part of the interpolation:
+  const CCTK_INT operand_indices[N_output_arrays];
+  const CCTK_INT opcodes        [N_output_arrays];
+The semantics here are that
+   output array[i] = op(input array[ operand_indices[i] ])
+where  op  is specified as an integer operation code as described below.
+
+Note that the array operand_indices[] doesn't directly name the inputs,
+but rather gives the indices (0-origin) in the list of inputs.  This
+allows for a more efficient implementation in the case where some of
+the input arrays have many different operations applied to them.
+
+The operations are coded based on the decimal digits of the integer:
+each decimal digit means to take the derivative in that direction;
+the order of the digits in a number is ignored.  For example:
+  0 = no derivative, "just" interpolate
+  1 = interpolate d/dx1 (derivative along x1 coordinate)
+  2 = interpolate d/dx2 (derivative along x2 coordinate)
+  11 = interpolate d^2/dx1^2 (2nd derivative along x1 coordinate)
+  22 = interpolate d^2/dx2^2 (2nd derivative along x2 coordinate)
+  12 = 21 = interpolate d^2/dx1 dx2 (mixed 2nd partial derivative in x1 and x2)
+  122 = 212 = 221 = interpolate d^3/dx1 dx2^2 (mixed 3rd partial derivative)
+  222 = interpolate d^3/dx2^3 (3rd derivative along x2 coordinate)
+  123 = 132 = 213 = 231 = 312 = 321
+      = interpolate d^3/dx1 dx2 dx3 (mixed 3rd partial derivative)
+
+After discussion with Tom Goodale, we have decided *not* to put #defines
+for the operation codes in any of the interpolator header files -- the
+operation codes are specific to this particular interpolation operator,
+not common to all operators, so they don't belong in the overall
+common-to-all header files.
+
+
+
+Pointers in Fortran
+===================
+
+One possible problem area with this API is that it requires creating
+arrays of pointers pointing to other arrays.  In C this is no problem,
+but in Fortran 77 this is difficult.  So, I propose adding two new Cactus
+functions to make this easier for Fortran users:
+
+   CCTK_POINTER Util_PointerTo(any Fortran variable or array)
+   CCTK_POINTER Util_NullPointer()
+
+Util_PointerTo would be #defined to %loc on those compilers which have
+that extension to standard Fortran, or would be a Cactus-provided utility
+routine for other cases.  It's trivial to write the latter case in C so
+long as the Fortran compiler actually uses call by reference; I've never
+heard of a Fortran compiler doing otherwise for arrays.  (And even for
+Fortran scalar variables it would be very hard for a compiler to do otherwise
+in light of separate compilation and 1-element arrays being allowed to be
+passed to/from scalar variables.)
+
+
+
+A Simple Example
+================
+
+Here's a simple example, written in Fortran 77, to do quadratic interpolation
+of a real and a complex local array in 3-D:
+
+c input arrays:
+        integer ni, nj, nk
+        parameter (ni=..., nj=..., nk=...)
+        CCTK_REAL    real_gridfn   (ni,nj,nk)
+        CCTK_COMPLEX complex_gridfn(ni,nj,nk)
+
+c interpolation coordinates
+        integer N_interp
+        parameter (N_interp = ...)
+        CCTK_REAL xcoord(N_interp), y_coord(N_interp), z_coord(N_interp)
+
+c output arrays:
+        CCTK_REAL    real_at_xyz   (N_interp)
+        CCTK_COMPLEX complex_at_xyz(N_interp)
+
+        integer status, dummy
+        CCTK_INT input_array_type_codes(2)
+        data input_array_type_codes /CCTK_VARIABLE_REAL,
+     $                               CCTK_VARIABLE_COMPLEX/
+        CCTK_INT input_array_dims(3)
+        CCTK_POINTER input_arrays(2)
+        CCTK_INT interp_coord_type_codes(3)
+        data interp_coord_type_codes /CCTK_VARIABLE_REAL,
+     $                                CCTK_VARIABLE_REAL,
+     $                                CCTK_VARIABLE_REAL/
+        CCTK_POINTER interp_coords(3)
+        CCTK_INT output_array_type_codes(2)
+        data output_array_type_codes /CCTK_VARIABLE_REAL,
+     $                                CCTK_VARIABLE_COMPLEX/
+        CCTK_POINTER output_arrays(2)
+
+        input_array_dims(1) = ni
+        input_array_dims(2) = nj
+        input_array_dims(3) = nk
+        interp_coords(1) = Util_PointerTo(xcoord)
+        interp_coords(2) = Util_PointerTo(ycoord)
+        interp_coords(3) = Util_PointerTo(zcoord)
+        output_arrays(1) = Util_PointerTo(real_at_xyz)
+        output_arrays(2) = Util_PointerTo(complex_at_xyz)
+
+        call CCTK_InterpLocalArrays
+     $          (status,                ! return code
+                 3,                     ! number of dimensions
+                 operator_handle, coord_system_handle,
+                 Util_TableCreateFromString("order=2"),
+                 N_interp,
+                 interp_coord_type_codes, interp_coords,
+                 2,                     ! number of input arrays
+                 input_array_type_codes, input_array_dims, input_arrays,
+                 2,                     ! number of output arrays
+                 output_array_type_codes, output_arrays)
+
+        if (status .lt. 0) then
+                call CCTK_WARN(status, "Error return from interpolator!")
+                call CCTK_Exit(dummy, Util_NullPointer(), status)
+        end if
+
+
+
+A More Complicated Example
+==========================
+
+Here's a more complicated example, written in C++.  (I'm really only using
+C++ to get cleaner initialization of the various arrays, this is still
+"almost C".)  This example is a simplified form of what I will be doing
+in my new apparent horizon finder:
+
+//
+// input grid functions (12 of them, all 3-D CCTK_REAL):
+//      gxx, gxy, gxz, gyy, gyz, gzz,
+//      Kxx, Kxy, Kxz, Kyy, Kyz, Kzz
+//
+// interpolation coordinates:
+//      xcoord, ycoord, zcoord   (all CCTK_REAL[N_interp_points])
+//
+// we want to interpolate the gij and Kij, and also interpolate all the
+// first derivatives of the gij, so the output arrays are
+// (30 of them, all CCTK_REAL[N_interp_points])
+//      I_gxx, dx_gxx, dy_gxx, dz_gxx,
+//      I_gxy, dx_gxy, dy_gxy, dz_gxy,
+//      I_gxz, dx_gxz, dy_gxz, dz_gxz,
+//      I_gyy, dx_gyy, dy_gyy, dz_gyy,
+//      I_gyz, dx_gyz, dy_gyz, dz_gyz,
+//      I_gzz, dx_gzz, dy_gzz, dz_gzz,
+//      I_Kxx, I_Kxy, I_Kxz, I_Kyy, I_Kyz, I_Kzz
+//
+
+#define VP(x) static_cast<void *>(x)
+
+const int N_dims = 3;
+const CCTK_INT interp_coord_type_codes[N_dims]
+        = { CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL };
+const void *const interp_coords[N_dims]
+        = { VP(xcoord), VP(ycoord), VP(zcoord) };
+
+const int N_input_arrays = 12;
+const CCTK_INT input_array_types[N_input_arrays]
+        = { CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+            CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+            CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+            CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL };
+
+const CCTK_INT input_array_variable_indices[N_input_arrays]
+        = { CCTK_VarIndex("somethorn::gxx"),
+            CCTK_VarIndex("somethorn::gxy"),
+            CCTK_VarIndex("somethorn::gxz"),
+            CCTK_VarIndex("somethorn::gyy"),
+            CCTK_VarIndex("somethorn::gyz"),
+            CCTK_VarIndex("somethorn::gzz"),
+            CCTK_VarIndex("somethorn::Kxx"),
+            CCTK_VarIndex("somethorn::Kxy"),
+            CCTK_VarIndex("somethorn::Kxz"),
+            CCTK_VarIndex("somethorn::Kyy"),
+            CCTK_VarIndex("somethorn::Kyz"),
+            CCTK_VarIndex("somethorn::Kzz") };
+
+const int N_output_arrays = 30;
+CCTK_INT output_array_type_codes[N_output_arrays];
+        for (int oi = 0 ; oi < N_output_arrays ; ++oi)
+        {
+        output_array_type_codes[oi] = CCTK_VARIABLE_REAL;
+        }
+
+void *const output_arrays[N_output_arrays]
+        = {
+          VP(I_gxx), VP(dx_gxx), VP(dy_gxx), VP(dz_gxx),
+          VP(I_gxy), VP(dx_gxy), VP(dy_gxy), VP(dz_gxy),
+          VP(I_gxz), VP(dx_gxz), VP(dy_gxz), VP(dz_gxz),
+          VP(I_gyy), VP(dx_gyy), VP(dy_gyy), VP(dz_gyy),
+          VP(I_gyz), VP(dx_gyz), VP(dy_gyz), VP(dz_gyz),
+          VP(I_gzz), VP(dx_gzz), VP(dy_gzz), VP(dz_gzz),
+          VP(I_Kxx), VP(I_Kxy), VP(I_Kxz), VP(I_Kyy), VP(I_Kyz), VP(I_Kzz)
+          };
+
+const CCTK_INT operand_indices[N_output_arrays];
+        = {
+          0, 0, 0, 0,           // gxx
+          1, 1, 1, 1,           // gxy
+          2, 2, 2, 2,           // gxz
+          3, 3, 3, 3,           // gyy
+          4, 4, 4, 4,           // gyz
+          5, 5, 5, 5,           // gzz
+          6, 7, 8, 9, 10, 11    // Kxx-Kzz
+          };
+
+const CCTK_INT opcodes[N_output_arrays]
+        = {
+          0, 1, 2, 3,           // I, dx, dy, dz
+          0, 1, 2, 3,           // I, dx, dy, dz
+          0, 1, 2, 3,           // I, dx, dy, dz
+          0, 1, 2, 3,           // I, dx, dy, dz
+          0, 1, 2, 3,           // I, dx, dy, dz
+          0, 1, 2, 3,           // I, dx, dy, dz
+          0, 0, 0, 0, 0, 0      // all I
+          };
+
+int param_table_handle = Util_TableCreate(UTIL_TABLE_DEFAULT);
+Util_TableSetInt(param_table_handle, 3, "order");
+Util_TableSetIntArray(param_table_handle,
+		      N_output_arrays, operand_indices,
+		      "operand_indices");
+Util_TableSetIntArray(param_table_handle,
+		      N_output_arrays, opcodes,
+		      "opcodes");
+
+int status = CCTK_InterpGridArrays(GH,
+                                   N_dims,
+                                   operator_handle, coord_system_handle,
+                                   param_table_handle,
+                                   N_interp_points,
+                                   interp_coord_type_codes, interp_coords,
+                                   N_input_arrays,
+                                   input_array_variable_indices,
+                                   N_output_arrays,
+                                   output_array_type_codes, output_arrays);
+if (status < 0)
+        {
+        CCTK_WARN(status, "error return from CCTK_InterpGridArrays()!");
+        CCTK_Exit(GH, status);                                  /*NOTREACHED*/
+        }
+Util_TableDestroy(param_table_handle);
diff --git a/doc/api2.2 b/doc/api2.2
new file mode 100644
index 0000000..35eca3c
--- /dev/null
+++ b/doc/api2.2
@@ -0,0 +1,496 @@
+Author: Jonathan Thornburg <jthorn@aei.mpg.de>
+Date: 30 December 2001
+
+This is version 2.2 of my proposal for the new interpolator API.  The
+main changes since version 1 are
+* function names have changed (longer, more descriptive)
+* output arrays may now be N-dimensional (before they were only 1-D)
+* a bunch more optional stuff in the table to specify array strides,
+  time levels, etc
+* no #defines in the header file for derivative operation codes
+* int arrays changed to CCTK_INT
+
+Don't be scared by the length, for most uses it's not that complicated!
+There are some examples below...
+
+
+
+Synopsis
+========
+
+  int status = CCTK_InterpLocalArrays(arguments described below)
+  int status = CCTK_InterpGridArrays(arguments described below)
+
+return is 0 for ok, various -ve values for error codes
+
+(N.b. the flesh APIs to register interpolation operators will also need
+their C function prototypes changed to match the changes here.)
+
+
+
+Function Arguments
+==================
+
+  /***** misc arguments *****/
+for CCTK_InterpGridArrays() only:
+  const cGH *GH;
+for both CCTK_InterpLocalArrays() and CCTK_InterpGridArrays():
+  /* note N_dims is the number of dimensions in the *interpolation*; */
+  /* this may be smaller than the number of dimensions of the input arrays */
+  /* if the storage indexing is set up appropriately (eg to interpolate */
+  /* in 1-D lines or 2-D planes of 3-D grid arrays) */
+  int N_dims;
+  int operator_handle;
+  int param_table_handle;       /* handle to "parameter table", a key-value */
+                                /* table, see below for table entries */
+
+  /***** arguments specifying the coordinate system *****/
+for CCTK_InterpLocalArrays():
+  /* the local coordinate system is specified as a generic linear mapping */
+  /* from (integer) local input array subscripts --> (global) coordinates: */
+  /*   coordinate = coord_system_origin[axis] + subscript*grid_spacing[axis] */
+  const CCTK_REAL coord_system_origin[N_dims];	/* coords of subscript 0 */
+  const CCTK_REAL grid_spacing[N_dims];
+for CCTK_InterpGridArrays():
+  int coord_system_handle;	/* handle to Cactus coordinate system */
+				/* specifying mapping of */
+				/* (integer) input array subscripts */
+				/* <--> (floating point) coordinates */
+
+  /***** arguments specifying the interpolation points *****/
+  int N_interp_points;
+  /* (pointer to) array of CCTK_VARIABLE_* codes giving the */
+  /* data types of the arrays pointed to by  interp_coords[] */
+  const CCTK_INT interp_coord_type_codes[N_dims];
+  /* (pointer to) array[N_dims] of pointers to arrays[N_interp_points] */
+  /* giving x,y,z,... coordinates of interpolation points */
+  const void *const interp_coords[N_dims];
+
+  /***** arguments specifying the inputs (the data to be interpolated) *****/
+  int N_input_arrays;
+for CCTK_InterpLocalArrays():
+  /* array of input array dimensions (common to all input arrays) */
+  const CCTK_INT input_array_dims[N_dims];
+  /* array of CCTK_VARIABLE_* codes giving data types of input arrays */
+  const CCTK_INT input_array_type_codes[N_input_arrays];
+  /* array of pointers to input arrays */
+  const void *const input_arrays[N_input_arrays];
+for CCTK_InterpGridArrays():
+  /* array of CCTK variable indices of input arrays */
+  const CCTK_INT input_array_variable_indices[N_input_arrays];
+
+for both CCTK_InterpLocalArrays() and CCTK_InterpGridArrays() again:
+  /***** arguments specifying the outputs (the interpolation results) *****/
+  int N_output_arrays;
+  /* array of CCTK_VARIABLE_* codes giving data types of output arrays */
+  const CCTK_INT output_array_type_codes[N_output_arrays];
+  /* array[N_output_arrays] of pointers to output arrays[N_interp_points] */
+  void *const output_arrays[N_output_arrays];
+
+
+
+Information Passed in the Parameter Table
+=========================================
+
+The "parameter table" may be used to specify non-default storage indexing
+for input or output arrays, and/or various options for the interpolation
+itself.  Some interpolators may not implement all of these options.
+
+
+Array Addressing/Subscripting Options
+-------------------------------------
+
+Sometimes one of the "arrays" used by the interpolator isn't contiguous
+in memory.  For example, we might want to do 2-D interpolation within a
+plane of a 3-D grid array, and/or the grid array might be a member of a
+compact group.  To support this, we use several optional table entries
+(these should be supported by all interpolation operators):
+
+For the input arrays, we use
+
+  const CCTK_INT input_array_offsets[N_input_arrays];
+  /* next 3 table entries are shared by all input arrays */
+  const CCTK_INT input_array_strides       [N_dims];
+  const CCTK_INT input_array_min_subscripts[N_dims];
+  const CCTK_INT input_array_max_subscripts[N_dims];
+
+Then for input array number a, the generic subscripting expression for
+the 3-D case is
+  data_pointer[offset + i*istride + j*jstride + k*kstride]
+where
+  data_pointer = input_arrays[a]
+  offset = input_array_offsets[a]
+  (istride,jstride,kstride) = input_array_stride[]
+and where (i,j,k) run from input_array_min_subscripts[] to
+input_array_max_subscripts[] inclusive.
+
+The defaults are offset=0, stride=determined from input_array_dims[]
+in the usual Fortran manner, input_array_min_subscripts[] = 0,
+input_array_max_subscripts[] = input_array_dims[]-1.  If the stride
+and max subscript are both specified explicitly, then the
+input_array_dims[] function argument is ignored.
+
+For CCTK_InterpGridArrays() operating on a member of a compact group
+the offset and strides are interpreted in units of _grid_points_.  This
+has the advantage that interpolator calls need not be changed if a grid
+array is changed from being simple to/from compact.  In terms of actual
+memory addressing, then, the internal subscripting expression for this
+case would be
+  group_data_pointer[offset + member_number + i*istride*N_members
+                                            + j*jstride*N_members
+                                            + k*kstride*N_members]
+
+For CCTK_InterpGridArrays(), by default the input (grid) arrays are at
+the "current" Cactus time level (level 0).  This may be changed with the
+table entry
+  const CCTK_INT input_array_time_levels[N_input_arrays];
+
+By default the interpolation-point coordinates and the output arrays
+are all contiguous 1-D arrays.  This may be changed with the table
+entries
+
+  const CCTK_INT interp_coords_offsets[N_dims];
+  const CCTK_INT output_array_offsets[N_output_arrays];
+  /* next 4 table entries are shared by all interp coords and output arrays */
+  const CCTK_INT interp_point_N_dims;
+  const CCTK_INT interp_point_strides       [interp_point_N_dims];
+  const CCTK_INT interp_point_min_subscripts[interp_point_N_dims];
+  const CCTK_INT interp_point_max_subscripts[interp_point_N_dims];
+
+For example, if we wanted to do 3-D interpolation, interpolating a value
+at each non-ghost-zone point of a 2-D grid of points, with the grid point
+coordinates stored as 2-D arrays, we would use
+  N_dims = 3
+  interp_point_N_dims = 2
+  interp_point_strides[] = set up from the full size of the 2-D grid
+  interp_point_{min,max}_subscripts[] = specify the non-ghost-zone points
+                                        of the 2-D grid
+
+Excision Options
+----------------
+
+Some interpolators may specifically support excision, where a mask array
+(same dimensionality and indexing as the input arrays) is used to mark
+some grid points as valid (ok to use data there) and other grid points
+as invalid (the interpolator isn't allowed to use data there).
+
+If an interpolator supports this, it should use the following optional
+parameters:
+
+for CCTK_InterpLocalArrays();
+  const CCTK_INT mask_type_code;  /* one of the CCTK_VARIABLE_* codes */
+  const void *const mask_array;   /* same dimensions/indexing as input arrays */
+for CCTK_InterpGridArrays():
+  const CCTK_INT mask_variable_index;
+
+for both CCTK_InterpLocalArrays() and CCTK_InterpGridArrays():
+  /* we consider a grid point to be valid if and only if the mask */
+  /* has a value in the closed interval [mask_valid_min,mask_valid_max] */
+  /* n.b. the caller should beware of possible rounding problems here; */
+  /*      it may be appropriate to widen the valid interval slightly */
+  /*      if the endpoints aren't exactly-representable floating-point */
+  /*      values */
+  const mask_type mask_valid_min, mask_valid_max;
+
+The same type of storage options supported for the input and/or output
+arrays, are also supported for the mask; the mask may have its own offset
+and/or time level, but shares any input-array stride and/or min/max
+subscript specification:
+
+  const CCTK_INT mask_offset;
+  const CCTK_INT mask_time_level;
+
+
+The remaining parameter-table options are specific to the new interpolator
+I'm currently implementing for PUGHInterp.  This registers (only) a single
+operator, "generalized polynomial interpolation".
+
+
+Interpolation Order and Molecule Family
+---------------------------------------
+
+The mandatory parameter
+
+  const CCTK_INT order;
+
+sets the order of the interpolating polynomial (1=linear, 2=quadratic,
+3=cubic, etc).  Thus the simplest call can just use (eg)
+  Util_TableCreateFromString("order=3")
+for cubic interpolation.
+
+All the remaining parameters in the table are optional; if they're
+omitted defaults will be supplied.
+
+  /* this selects one of a family of related operators */
+  /* the default (and the only one I'm implementing right now) */
+  /* is "cube" to use the usual hypercube-shaped molecules */
+  const char *const molecule_family;
+
+Smoothing
+---------
+
+The way I'm implementing the interpolation it's easy to also do
+Savitzky-Golay type smoothing (= moving least-square fitting, cf
+Numerical Recipes 2nd edition section 14.8).  This is controlled by
+the parameter
+
+  const CCTK_INT smoothing;
+
+which says how much (how many points) to enlarge the interpolation
+molecule for this.  The default is 0 (no smoothing).  1 would mean to
+enlarge the molecule by 1 point, e.g. to use a 5-point molecule instead
+of the usual 4-point one for cubic interpolation.  2 would mean to
+enlarge by 2 points, e.g. to use a 6-point molecule for cubic
+interpolation.  Etc etc.
+
+This type of smoothing is basically free apart from the increase in
+the molecule size, e.g. a smoothing=2 cubic interpolation has exactly
+the same cost as any other 6-point-molecule interpolation.
+
+Derivatives
+-----------
+
+This interpolator can optionally (and again at no extra cost) take
+partial derivatives as part of the interpolation:
+  const CCTK_INT operand_indices[N_output_arrays];
+  const CCTK_INT opcodes        [N_output_arrays];
+The semantics here are that
+   output array[i] = op(input array[ operand_indices[i] ])
+where  op  is specified as an integer operation code as described below.
+
+Note that the array operand_indices[] doesn't directly name the inputs,
+but rather gives the indices (0-origin) in the list of inputs.  This
+allows for a more efficient implementation in the case where some of
+the input arrays have many different operations applied to them.
+
+The operations are coded based on the decimal digits of the integer:
+each decimal digit means to take the derivative in that direction;
+the order of the digits in a number is ignored.  For example:
+  0 = no derivative, "just" interpolate
+  1 = interpolate d/dx1 (derivative along x1 coordinate)
+  2 = interpolate d/dx2 (derivative along x2 coordinate)
+  11 = interpolate d^2/dx1^2 (2nd derivative along x1 coordinate)
+  22 = interpolate d^2/dx2^2 (2nd derivative along x2 coordinate)
+  12 = 21 = interpolate d^2/dx1 dx2 (mixed 2nd partial derivative in x1 and x2)
+  122 = 212 = 221 = interpolate d^3/dx1 dx2^2 (mixed 3rd partial derivative)
+  222 = interpolate d^3/dx2^3 (3rd derivative along x2 coordinate)
+  123 = 132 = 213 = 231 = 312 = 321
+      = interpolate d^3/dx1 dx2 dx3 (mixed 3rd partial derivative)
+
+After discussion with Tom Goodale, we have decided *not* to put #defines
+for the operation codes in any of the interpolator header files -- the
+operation codes are specific to this particular interpolation operator,
+not common to all operators, so they don't belong in the overall
+common-to-all header files.
+
+
+
+Pointers in Fortran
+===================
+
+One possible problem area with this API is that it requires creating
+arrays of pointers pointing to other arrays.  In C this is no problem,
+but in Fortran 77 this is difficult.  So, I propose adding two new Cactus
+functions to make this easier for Fortran users:
+
+   CCTK_POINTER Util_PointerTo(any Fortran variable or array)
+   CCTK_POINTER Util_NullPointer()
+
+Util_PointerTo would be #defined to %loc on those compilers which have
+that extension to standard Fortran, or would be a Cactus-provided utility
+routine for other cases.  It's trivial to write the latter case in C so
+long as the Fortran compiler actually uses call by reference; I've never
+heard of a Fortran compiler doing otherwise for arrays.  (And even for
+Fortran scalar variables it would be very hard for a compiler to do otherwise
+in light of separate compilation and 1-element arrays being allowed to be
+passed to/from scalar variables.)
+
+
+
+A Simple Example
+================
+
+Here's a simple example, written in Fortran 77, to do quadratic interpolation
+of a real and a complex local array in 3-D:
+
+c input arrays:
+        integer ni, nj, nk
+        parameter (ni=..., nj=..., nk=...)
+        CCTK_REAL    real_gridfn   (ni,nj,nk)
+        CCTK_COMPLEX complex_gridfn(ni,nj,nk)
+
+c interpolation coordinates
+        integer N_interp
+        parameter (N_interp = ...)
+        CCTK_REAL xcoord(N_interp), y_coord(N_interp), z_coord(N_interp)
+
+c output arrays:
+        CCTK_REAL    real_at_xyz   (N_interp)
+        CCTK_COMPLEX complex_at_xyz(N_interp)
+
+        integer status, dummy
+        CCTK_INT input_array_type_codes(2)
+        data input_array_type_codes /CCTK_VARIABLE_REAL,
+     $                               CCTK_VARIABLE_COMPLEX/
+        CCTK_INT input_array_dims(3)
+        CCTK_POINTER input_arrays(2)
+        CCTK_INT interp_coord_type_codes(3)
+        data interp_coord_type_codes /CCTK_VARIABLE_REAL,
+     $                                CCTK_VARIABLE_REAL,
+     $                                CCTK_VARIABLE_REAL/
+        CCTK_POINTER interp_coords(3)
+        CCTK_INT output_array_type_codes(2)
+        data output_array_type_codes /CCTK_VARIABLE_REAL,
+     $                                CCTK_VARIABLE_COMPLEX/
+        CCTK_POINTER output_arrays(2)
+
+        input_array_dims(1) = ni
+        input_array_dims(2) = nj
+        input_array_dims(3) = nk
+        interp_coords(1) = Util_PointerTo(xcoord)
+        interp_coords(2) = Util_PointerTo(ycoord)
+        interp_coords(3) = Util_PointerTo(zcoord)
+        output_arrays(1) = Util_PointerTo(real_at_xyz)
+        output_arrays(2) = Util_PointerTo(complex_at_xyz)
+
+        call CCTK_InterpLocalArrays
+     $          (status,                ! return code
+                 3,                     ! number of dimensions
+                 operator_handle, coord_system_handle,
+                 Util_TableCreateFromString("order=2"),
+                 N_interp,
+                 interp_coord_type_codes, interp_coords,
+                 2,                     ! number of input arrays
+                 input_array_type_codes, input_array_dims, input_arrays,
+                 2,                     ! number of output arrays
+                 output_array_type_codes, output_arrays)
+
+        if (status .lt. 0) then
+                call CCTK_WARN(status, "Error return from interpolator!")
+                call CCTK_Exit(dummy, Util_NullPointer(), status)
+        end if
+
+
+
+A More Complicated Example
+==========================
+
+Here's a more complicated example, written in C++.  (I'm really only using
+C++ to get cleaner initialization of the various arrays, this is still
+"almost C".)  This example is a simplified form of what I will be doing
+in my new apparent horizon finder:
+
+//
+// input grid functions (12 of them, all 3-D CCTK_REAL):
+//      gxx, gxy, gxz, gyy, gyz, gzz,
+//      Kxx, Kxy, Kxz, Kyy, Kyz, Kzz
+//
+// interpolation coordinates:
+//      xcoord, ycoord, zcoord   (all CCTK_REAL[N_interp_points])
+//
+// we want to interpolate the gij and Kij, and also interpolate all the
+// first derivatives of the gij, so the output arrays are
+// (30 of them, all CCTK_REAL[N_interp_points])
+//      I_gxx, dx_gxx, dy_gxx, dz_gxx,
+//      I_gxy, dx_gxy, dy_gxy, dz_gxy,
+//      I_gxz, dx_gxz, dy_gxz, dz_gxz,
+//      I_gyy, dx_gyy, dy_gyy, dz_gyy,
+//      I_gyz, dx_gyz, dy_gyz, dz_gyz,
+//      I_gzz, dx_gzz, dy_gzz, dz_gzz,
+//      I_Kxx, I_Kxy, I_Kxz, I_Kyy, I_Kyz, I_Kzz
+//
+
+#define VP(x) static_cast<void *>(x)
+
+const int N_dims = 3;
+const CCTK_INT interp_coord_type_codes[N_dims]
+        = { CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL };
+const void *const interp_coords[N_dims]
+        = { VP(xcoord), VP(ycoord), VP(zcoord) };
+
+const int N_input_arrays = 12;
+const CCTK_INT input_array_types[N_input_arrays]
+        = { CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+            CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+            CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+            CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL };
+
+const CCTK_INT input_array_variable_indices[N_input_arrays]
+        = { CCTK_VarIndex("somethorn::gxx"),
+            CCTK_VarIndex("somethorn::gxy"),
+            CCTK_VarIndex("somethorn::gxz"),
+            CCTK_VarIndex("somethorn::gyy"),
+            CCTK_VarIndex("somethorn::gyz"),
+            CCTK_VarIndex("somethorn::gzz"),
+            CCTK_VarIndex("somethorn::Kxx"),
+            CCTK_VarIndex("somethorn::Kxy"),
+            CCTK_VarIndex("somethorn::Kxz"),
+            CCTK_VarIndex("somethorn::Kyy"),
+            CCTK_VarIndex("somethorn::Kyz"),
+            CCTK_VarIndex("somethorn::Kzz") };
+
+const int N_output_arrays = 30;
+CCTK_INT output_array_type_codes[N_output_arrays];
+        for (int oi = 0 ; oi < N_output_arrays ; ++oi)
+        {
+        output_array_type_codes[oi] = CCTK_VARIABLE_REAL;
+        }
+
+void *const output_arrays[N_output_arrays]
+        = {
+          VP(I_gxx), VP(dx_gxx), VP(dy_gxx), VP(dz_gxx),
+          VP(I_gxy), VP(dx_gxy), VP(dy_gxy), VP(dz_gxy),
+          VP(I_gxz), VP(dx_gxz), VP(dy_gxz), VP(dz_gxz),
+          VP(I_gyy), VP(dx_gyy), VP(dy_gyy), VP(dz_gyy),
+          VP(I_gyz), VP(dx_gyz), VP(dy_gyz), VP(dz_gyz),
+          VP(I_gzz), VP(dx_gzz), VP(dy_gzz), VP(dz_gzz),
+          VP(I_Kxx), VP(I_Kxy), VP(I_Kxz), VP(I_Kyy), VP(I_Kyz), VP(I_Kzz)
+          };
+
+const CCTK_INT operand_indices[N_output_arrays];
+        = {
+          0, 0, 0, 0,           // gxx
+          1, 1, 1, 1,           // gxy
+          2, 2, 2, 2,           // gxz
+          3, 3, 3, 3,           // gyy
+          4, 4, 4, 4,           // gyz
+          5, 5, 5, 5,           // gzz
+          6, 7, 8, 9, 10, 11    // Kxx-Kzz
+          };
+
+const CCTK_INT opcodes[N_output_arrays]
+        = {
+          0, 1, 2, 3,           // I, dx, dy, dz
+          0, 1, 2, 3,           // I, dx, dy, dz
+          0, 1, 2, 3,           // I, dx, dy, dz
+          0, 1, 2, 3,           // I, dx, dy, dz
+          0, 1, 2, 3,           // I, dx, dy, dz
+          0, 1, 2, 3,           // I, dx, dy, dz
+          0, 0, 0, 0, 0, 0      // all I
+          };
+
+int param_table_handle = Util_TableCreate(UTIL_TABLE_DEFAULT);
+Util_TableSetInt(param_table_handle, 3, "order");
+Util_TableSetIntArray(param_table_handle,
+		      N_output_arrays, operand_indices,
+		      "operand_indices");
+Util_TableSetIntArray(param_table_handle,
+		      N_output_arrays, opcodes,
+		      "opcodes");
+
+int status = CCTK_InterpGridArrays(GH,
+                                   N_dims,
+                                   operator_handle, coord_system_handle,
+                                   param_table_handle,
+                                   N_interp_points,
+                                   interp_coord_type_codes, interp_coords,
+                                   N_input_arrays,
+                                   input_array_variable_indices,
+                                   N_output_arrays,
+                                   output_array_type_codes, output_arrays);
+if (status < 0)
+        {
+        CCTK_WARN(status, "error return from CCTK_InterpGridArrays()!");
+        CCTK_Exit(GH, status);                                  /*NOTREACHED*/
+        }
+Util_TableDestroy(param_table_handle);
diff --git a/doc/api3 b/doc/api3
new file mode 100644
index 0000000..9aeedd5
--- /dev/null
+++ b/doc/api3
@@ -0,0 +1,491 @@
+Author: Jonathan Thornburg <jthorn@aei.mpg.de>
+Date: 21 January 2002
+
+This is version 3.0 of my proposal for the new interpolator API.  The
+main changes since previous versions are
+* This only describes processor-local interpolation -- global (or grid)
+  interpolation is trickier, and we're deferring it till after the
+  processor-local code is finished
+* The API has split again to separate out {uniform, nonuniform, irregular}
+  grids.  This proposal only addresses the first two; irregular is hard
+  and needs lots of work all over Cactus before we can support it.
+* Since noone seems to care about it, I've dropped support for different
+  coordinates being of different types (i.e. x is CCTK_REAL4 but y is
+  CCTK_REAL16); now all the coordinates must be the same type.
+* Ditto for the interpolation coordinates.
+
+Don't be scared by the length, for most uses it's not that complicated!
+There are some examples below...
+
+
+
+Synopsis
+========
+
+  int status = CCTK_InterpLocalUniform(arguments described below)
+  int status = CCTK_InterpLocalNonUniform(arguments described below)
+
+return is 0 for ok, various -ve values for error codes
+
+(N.b. the flesh APIs to register interpolation operators will also need
+their C function prototypes changed to match the changes here.)
+
+
+
+Function Arguments
+==================
+
+  /***** misc arguments *****/
+  /* note N_dims is the number of dimensions in the *interpolation*; */
+  /* this may be smaller than the number of dimensions of the input arrays */
+  /* if the storage indexing is set up appropriately (eg to interpolate */
+  /* in 1-D lines or 2-D planes of 3-D grid arrays) */
+  int N_dims;
+  int operator_handle;
+  int param_table_handle;       /* handle to a key-value table used to pass */
+                                /* additional parameters to the interpolator */
+
+  /***** arguments specifying the local coordinate system *****/
+for CCTK_InterpLocalUniform()
+  /* the local coordinate system is specified as a generic linear mapping */
+  /* from (integer) local input array subscripts --> (global) coordinates: */
+  /*   coordinate = origin[axis] + subscript*delta[axis] */
+  const CCTK_REAL origin[N_dims];	/* coords of subscript 0 */
+  const CCTK_REAL delta[N_dims];
+for CCTK_InterpLocalNonuniform()
+  /* the local coordinate system is specified by user-supplied arrays: */
+  /*   coordinate = array[subscript]  (separate array for each axis) */
+  /* FIXME: what if subscript=0 is out-of-range??? */
+  const void *const coords[N_dims];	/* coords of subscript 0 */
+
+  /***** arguments specifying the interpolation points *****/
+  int N_interp_points;
+  int interp_coords_type_code;		/* CCTK_VARIABLE_* type code for */
+					/* interp_coords arrays */
+  /* (pointer to) array[N_dims] of pointers to arrays[N_interp_points] */
+  /* giving x,y,z,... coordinates of interpolation points */
+  const void *const interp_coords[N_dims];
+
+  /***** arguments specifying the inputs (the data to be interpolated) *****/
+  int N_input_arrays;
+  /* array of input array dimensions (common to all input arrays) */
+  const CCTK_INT input_array_dims[N_dims];
+  /* array of CCTK_VARIABLE_* codes giving data types of input arrays */
+  const CCTK_INT input_array_type_codes[N_input_arrays];
+  /* array of pointers to input arrays */
+  const void *const input_arrays[N_input_arrays];
+
+  /***** arguments specifying the outputs (the interpolation results) *****/
+  int N_output_arrays;
+  /* array of CCTK_VARIABLE_* codes giving data types of output arrays */
+  const CCTK_INT output_array_type_codes[N_output_arrays];
+  /* array[N_output_arrays] of pointers to output arrays[N_interp_points] */
+  void *const output_arrays[N_output_arrays];
+
+
+
+Information Passed in the Parameter Table
+=========================================
+
+The "parameter table" may be used to specify non-default storage indexing
+for input or output arrays, and/or various options for the interpolation
+itself.  Some interpolators may not implement all of these options.
+
+
+Array Addressing/Subscripting Options
+-------------------------------------
+
+Sometimes one of the "arrays" used by the interpolator isn't contiguous
+in memory.  For example, we might want to do 2-D interpolation within a
+plane of a 3-D grid array, and/or the grid array might be a member of a
+compact group.  To support this, we use several optional table entries
+(these should be supported by all interpolation operators):
+
+For the input arrays, we use
+
+  const CCTK_INT input_array_offsets[N_input_arrays];
+  /* next 3 table entries are shared by all input arrays */
+  const CCTK_INT input_array_strides       [N_dims];
+  const CCTK_INT input_array_min_subscripts[N_dims];
+  const CCTK_INT input_array_max_subscripts[N_dims];
+
+Then for input array number a, the generic subscripting expression for
+the 3-D case is
+  data_pointer[offset + i*istride + j*jstride + k*kstride]
+where
+  data_pointer = input_arrays[a]
+  offset = input_array_offsets[a]
+  (istride,jstride,kstride) = input_array_stride[]
+and where (i,j,k) run from input_array_min_subscripts[] to
+input_array_max_subscripts[] inclusive.
+
+The defaults are offset=0, stride=determined from input_array_dims[]
+in the usual Fortran manner, input_array_min_subscripts[] = 0,
+input_array_max_subscripts[] = input_array_dims[]-1.  If the stride
+and max subscript are both specified explicitly, then the
+input_array_dims[] function argument is ignored.
+
+For CCTK_InterpGridArrays() operating on a member of a compact group
+the offset and strides are interpreted in units of _grid_points_.  This
+has the advantage that interpolator calls need not be changed if a grid
+array is changed from being simple to/from compact.  In terms of actual
+memory addressing, then, the internal subscripting expression for this
+case would be
+  group_data_pointer[offset + member_number + i*istride*N_members
+                                            + j*jstride*N_members
+                                            + k*kstride*N_members]
+
+For CCTK_InterpGridArrays(), by default the input (grid) arrays are at
+the "current" Cactus time level (level 0).  This may be changed with the
+table entry
+  const CCTK_INT input_array_time_levels[N_input_arrays];
+
+By default the interpolation-point coordinates and the output arrays
+are all contiguous 1-D arrays.  This may be changed with the table
+entries
+
+  const CCTK_INT interp_coords_offsets[N_dims];
+  const CCTK_INT output_array_offsets[N_output_arrays];
+  /* next 4 table entries are shared by all interp coords and output arrays */
+  const CCTK_INT interp_point_N_dims;
+  const CCTK_INT interp_point_strides       [interp_point_N_dims];
+  const CCTK_INT interp_point_min_subscripts[interp_point_N_dims];
+  const CCTK_INT interp_point_max_subscripts[interp_point_N_dims];
+
+For example, if we wanted to do 3-D interpolation, interpolating a value
+at each non-ghost-zone point of a 2-D grid of points, with the grid point
+coordinates stored as 2-D arrays, we would use
+  N_dims = 3
+  interp_point_N_dims = 2
+  interp_point_strides[] = set up from the full size of the 2-D grid
+  interp_point_{min,max}_subscripts[] = specify the non-ghost-zone points
+                                        of the 2-D grid
+
+Excision Options
+----------------
+
+Some interpolators may specifically support excision, where a mask array
+(same dimensionality and indexing as the input arrays) is used to mark
+some grid points as valid (ok to use data there) and other grid points
+as invalid (the interpolator isn't allowed to use data there).
+
+If an interpolator supports this, it should use the following optional
+parameters:
+
+for CCTK_InterpLocalArrays();
+  const CCTK_INT mask_type_code;  /* one of the CCTK_VARIABLE_* codes */
+  const void *const mask_array;   /* same dimensions/indexing as input arrays */
+for CCTK_InterpGridArrays():
+  const CCTK_INT mask_variable_index;
+
+for both CCTK_InterpLocalArrays() and CCTK_InterpGridArrays():
+  /* we consider a grid point to be valid if and only if the mask */
+  /* has a value in the closed interval [mask_valid_min,mask_valid_max] */
+  /* n.b. the caller should beware of possible rounding problems here; */
+  /*      it may be appropriate to widen the valid interval slightly */
+  /*      if the endpoints aren't exactly-representable floating-point */
+  /*      values */
+  const mask_type mask_valid_min, mask_valid_max;
+
+The same type of storage options supported for the input and/or output
+arrays, are also supported for the mask; the mask may have its own offset
+and/or time level, but shares any input-array stride and/or min/max
+subscript specification:
+
+  const CCTK_INT mask_offset;
+  const CCTK_INT mask_time_level;
+
+
+The remaining parameter-table options are specific to the new interpolator
+I'm currently implementing for PUGHInterp.  This registers (only) a single
+operator, "generalized polynomial interpolation".
+
+
+Interpolation Order and Molecule Family
+---------------------------------------
+
+The mandatory parameter
+
+  const CCTK_INT order;
+
+sets the order of the interpolating polynomial (1=linear, 2=quadratic,
+3=cubic, etc).  Thus the simplest call can just use (eg)
+  Util_TableCreateFromString("order=3")
+for cubic interpolation.
+
+All the remaining parameters in the table are optional; if they're
+omitted defaults will be supplied.
+
+  /* this selects one of a family of related operators */
+  /* the default (and the only one I'm implementing right now) */
+  /* is "cube" to use the usual hypercube-shaped molecules */
+  const char *const molecule_family;
+
+Smoothing
+---------
+
+The way I'm implementing the interpolation it's easy to also do
+Savitzky-Golay type smoothing (= moving least-square fitting, cf
+Numerical Recipes 2nd edition section 14.8).  This is controlled by
+the parameter
+
+  const CCTK_INT smoothing;
+
+which says how much (how many points) to enlarge the interpolation
+molecule for this.  The default is 0 (no smoothing).  1 would mean to
+enlarge the molecule by 1 point, e.g. to use a 5-point molecule instead
+of the usual 4-point one for cubic interpolation.  2 would mean to
+enlarge by 2 points, e.g. to use a 6-point molecule for cubic
+interpolation.  Etc etc.
+
+This type of smoothing is basically free apart from the increase in
+the molecule size, e.g. a smoothing=2 cubic interpolation has exactly
+the same cost as any other 6-point-molecule interpolation.
+
+Derivatives
+-----------
+
+This interpolator can optionally (and again at no extra cost) take
+partial derivatives as part of the interpolation:
+  const CCTK_INT operand_indices[N_output_arrays];
+  const CCTK_INT operation_codes[N_output_arrays];
+The semantics here are that
+   output array[i] = op(input array[ operand_indices[i] ])
+where  op  is specified as an integer operation code as described below.
+
+Note that the array operand_indices[] doesn't directly name the inputs,
+but rather gives the indices (0-origin) in the list of inputs.  This
+allows for a more efficient implementation in the case where some of
+the input arrays have many different operations applied to them.
+
+The operations are coded based on the decimal digits of the integer:
+each decimal digit means to take the derivative in that direction;
+the order of the digits in a number is ignored.  For example:
+  0 = no derivative, "just" interpolate
+  1 = interpolate d/dx1 (derivative along x1 coordinate)
+  2 = interpolate d/dx2 (derivative along x2 coordinate)
+  11 = interpolate d^2/dx1^2 (2nd derivative along x1 coordinate)
+  22 = interpolate d^2/dx2^2 (2nd derivative along x2 coordinate)
+  12 = 21 = interpolate d^2/dx1 dx2 (mixed 2nd partial derivative in x1 and x2)
+  122 = 212 = 221 = interpolate d^3/dx1 dx2^2 (mixed 3rd partial derivative)
+  222 = interpolate d^3/dx2^3 (3rd derivative along x2 coordinate)
+  123 = 132 = 213 = 231 = 312 = 321
+      = interpolate d^3/dx1 dx2 dx3 (mixed 3rd partial derivative)
+
+After discussion with Tom Goodale, we have decided *not* to put #defines
+for the operation codes in any of the interpolator header files -- the
+operation codes are specific to this particular interpolation operator,
+not common to all operators, so they don't belong in the overall
+common-to-all header files.
+
+
+
+Pointers in Fortran
+===================
+
+One possible problem area with this API is that it requires creating
+arrays of pointers pointing to other arrays.  In C this is no problem,
+but in Fortran 77 this is difficult.  So, I propose adding two new Cactus
+functions to make this easier for Fortran users:
+
+   CCTK_POINTER Util_PointerTo(any Fortran variable or array)
+   CCTK_POINTER Util_NullPointer()
+
+Util_PointerTo would be #defined to %loc on those compilers which have
+that extension to standard Fortran, or would be a Cactus-provided utility
+routine for other cases.  It's trivial to write the latter case in C so
+long as the Fortran compiler actually uses call by reference; I've never
+heard of a Fortran compiler doing otherwise for arrays.  (And even for
+Fortran scalar variables it would be very hard for a compiler to do otherwise
+in light of separate compilation and 1-element arrays being allowed to be
+passed to/from scalar variables.)
+
+
+
+A Simple Example
+================
+
+Here's a simple example, written in Fortran 77, to do quadratic interpolation
+of a real and a complex local array in 3-D:
+
+c input arrays:
+        integer ni, nj, nk
+        parameter (ni=..., nj=..., nk=...)
+        CCTK_REAL    real_gridfn   (ni,nj,nk)
+        CCTK_COMPLEX complex_gridfn(ni,nj,nk)
+
+c interpolation coordinates
+        integer N_interp
+        parameter (N_interp = ...)
+        CCTK_REAL xcoord(N_interp), y_coord(N_interp), z_coord(N_interp)
+
+c output arrays:
+        CCTK_REAL    real_at_xyz   (N_interp)
+        CCTK_COMPLEX complex_at_xyz(N_interp)
+
+        integer status, dummy
+        CCTK_INT input_array_type_codes(2)
+        data input_array_type_codes /CCTK_VARIABLE_REAL,
+     $                               CCTK_VARIABLE_COMPLEX/
+        CCTK_INT input_array_dims(3)
+        CCTK_POINTER input_arrays(2)
+        CCTK_INT interp_coord_type_codes(3)
+        data interp_coord_type_codes /CCTK_VARIABLE_REAL,
+     $                                CCTK_VARIABLE_REAL,
+     $                                CCTK_VARIABLE_REAL/
+        CCTK_POINTER interp_coords(3)
+        CCTK_INT output_array_type_codes(2)
+        data output_array_type_codes /CCTK_VARIABLE_REAL,
+     $                                CCTK_VARIABLE_COMPLEX/
+        CCTK_POINTER output_arrays(2)
+
+        input_array_dims(1) = ni
+        input_array_dims(2) = nj
+        input_array_dims(3) = nk
+        interp_coords(1) = Util_PointerTo(xcoord)
+        interp_coords(2) = Util_PointerTo(ycoord)
+        interp_coords(3) = Util_PointerTo(zcoord)
+        output_arrays(1) = Util_PointerTo(real_at_xyz)
+        output_arrays(2) = Util_PointerTo(complex_at_xyz)
+
+        call CCTK_InterpLocalArrays
+     $          (status,                ! return code
+                 3,                     ! number of dimensions
+                 operator_handle, coord_system_handle,
+                 Util_TableCreateFromString("order=2"),
+                 N_interp,
+                 interp_coord_type_codes, interp_coords,
+                 2,                     ! number of input arrays
+                 input_array_type_codes, input_array_dims, input_arrays,
+                 2,                     ! number of output arrays
+                 output_array_type_codes, output_arrays)
+
+        if (status .lt. 0) then
+                call CCTK_WARN(status, "Error return from interpolator!")
+                call CCTK_Exit(dummy, Util_NullPointer(), status)
+        end if
+
+
+
+A More Complicated Example
+==========================
+
+Here's a more complicated example, written in C++.  (I'm really only using
+C++ to get cleaner initialization of the various arrays, this is still
+"almost C".)  This example is a simplified form of what I will be doing
+in my new apparent horizon finder:
+
+//
+// input grid functions (12 of them, all 3-D CCTK_REAL):
+//      gxx, gxy, gxz, gyy, gyz, gzz,
+//      Kxx, Kxy, Kxz, Kyy, Kyz, Kzz
+//
+// interpolation coordinates:
+//      xcoord, ycoord, zcoord   (all CCTK_REAL[N_interp_points])
+//
+// we want to interpolate the gij and Kij, and also interpolate all the
+// first derivatives of the gij, so the output arrays are
+// (30 of them, all CCTK_REAL[N_interp_points])
+//      I_gxx, dx_gxx, dy_gxx, dz_gxx,
+//      I_gxy, dx_gxy, dy_gxy, dz_gxy,
+//      I_gxz, dx_gxz, dy_gxz, dz_gxz,
+//      I_gyy, dx_gyy, dy_gyy, dz_gyy,
+//      I_gyz, dx_gyz, dy_gyz, dz_gyz,
+//      I_gzz, dx_gzz, dy_gzz, dz_gzz,
+//      I_Kxx, I_Kxy, I_Kxz, I_Kyy, I_Kyz, I_Kzz
+//
+
+#define VP(x) static_cast<void *>(x)
+
+const int N_dims = 3;
+const CCTK_INT interp_coord_type_codes[N_dims]
+        = { CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL };
+const void *const interp_coords[N_dims]
+        = { VP(xcoord), VP(ycoord), VP(zcoord) };
+
+const int N_input_arrays = 12;
+const CCTK_INT input_array_types[N_input_arrays]
+        = { CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+            CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+            CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+            CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL };
+
+const CCTK_INT input_array_variable_indices[N_input_arrays]
+        = { CCTK_VarIndex("somethorn::gxx"),
+            CCTK_VarIndex("somethorn::gxy"),
+            CCTK_VarIndex("somethorn::gxz"),
+            CCTK_VarIndex("somethorn::gyy"),
+            CCTK_VarIndex("somethorn::gyz"),
+            CCTK_VarIndex("somethorn::gzz"),
+            CCTK_VarIndex("somethorn::Kxx"),
+            CCTK_VarIndex("somethorn::Kxy"),
+            CCTK_VarIndex("somethorn::Kxz"),
+            CCTK_VarIndex("somethorn::Kyy"),
+            CCTK_VarIndex("somethorn::Kyz"),
+            CCTK_VarIndex("somethorn::Kzz") };
+
+const int N_output_arrays = 30;
+CCTK_INT output_array_type_codes[N_output_arrays];
+        for (int oi = 0 ; oi < N_output_arrays ; ++oi)
+        {
+        output_array_type_codes[oi] = CCTK_VARIABLE_REAL;
+        }
+
+void *const output_arrays[N_output_arrays]
+        = {
+          VP(I_gxx), VP(dx_gxx), VP(dy_gxx), VP(dz_gxx),
+          VP(I_gxy), VP(dx_gxy), VP(dy_gxy), VP(dz_gxy),
+          VP(I_gxz), VP(dx_gxz), VP(dy_gxz), VP(dz_gxz),
+          VP(I_gyy), VP(dx_gyy), VP(dy_gyy), VP(dz_gyy),
+          VP(I_gyz), VP(dx_gyz), VP(dy_gyz), VP(dz_gyz),
+          VP(I_gzz), VP(dx_gzz), VP(dy_gzz), VP(dz_gzz),
+          VP(I_Kxx), VP(I_Kxy), VP(I_Kxz), VP(I_Kyy), VP(I_Kyz), VP(I_Kzz)
+          };
+
+const CCTK_INT operand_indices[N_output_arrays];
+        = {
+          0, 0, 0, 0,           // gxx
+          1, 1, 1, 1,           // gxy
+          2, 2, 2, 2,           // gxz
+          3, 3, 3, 3,           // gyy
+          4, 4, 4, 4,           // gyz
+          5, 5, 5, 5,           // gzz
+          6, 7, 8, 9, 10, 11    // Kxx-Kzz
+          };
+
+const CCTK_INT operation_codes[N_output_arrays]
+        = {
+          0, 1, 2, 3,           // I, dx, dy, dz
+          0, 1, 2, 3,           // I, dx, dy, dz
+          0, 1, 2, 3,           // I, dx, dy, dz
+          0, 1, 2, 3,           // I, dx, dy, dz
+          0, 1, 2, 3,           // I, dx, dy, dz
+          0, 1, 2, 3,           // I, dx, dy, dz
+          0, 0, 0, 0, 0, 0      // all I
+          };
+
+int param_table_handle = Util_TableCreate(UTIL_TABLE_DEFAULT);
+Util_TableSetInt(param_table_handle, 3, "order");
+Util_TableSetIntArray(param_table_handle,
+		      N_output_arrays, operand_indices,
+		      "operand_indices");
+Util_TableSetIntArray(param_table_handle,
+		      N_output_arrays, operation_codes,
+		      "operation_codes");
+
+int status = CCTK_InterpGridArrays(GH,
+                                   N_dims,
+                                   operator_handle, coord_system_handle,
+                                   param_table_handle,
+                                   N_interp_points,
+                                   interp_coord_type_codes, interp_coords,
+                                   N_input_arrays,
+                                   input_array_variable_indices,
+                                   N_output_arrays,
+                                   output_array_type_codes, output_arrays);
+if (status < 0)
+        {
+        CCTK_WARN(status, "error return from CCTK_InterpGridArrays()!");
+        CCTK_Exit(GH, status);                                  /*NOTREACHED*/
+        }
+Util_TableDestroy(param_table_handle);
diff --git a/doc/documentation.tex b/doc/documentation.tex
new file mode 100644
index 0000000..5ac372c
--- /dev/null
+++ b/doc/documentation.tex
@@ -0,0 +1,1063 @@
+%version   $Header$
+\documentclass{article}
+
+\def\eqref#1{$(\ref{#1})$}
+\def\cf{cf.\hbox{}}
+\def\ie{i.e.\hbox{}}
+\def\eg{e.g.\hbox{}}
+\def\etal{{\it et~al.\/}}
+
+\def\defn#1{{\bf #1}}
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\begin{document}
+\title{LocalInterp}
+\author{Jonathan Thornburg, incorporating code from Thomas Radke}
+%
+% We want CVS to expand the Id keyword on the next line, but we don't
+% want TeX to go into math mode to typeset the expansion (because that
+% wouldn't look as nice in the output), so we use the "$ $" construct
+% to get TeX out of math mode again when typesetting the expansion.
+%
+\date{$ $Id$ $}
+\maketitle
+
+\abstract{
+	 This thorn provides processor-local interpolation of
+	 1-D, 2-D, and 3-D data arrays.  It provides several
+	 different interpolators, and supports both the older
+	 \verb|CCTK_InterpLocal()| API, and the newer
+	 \verb|CCTK_InterpLocalUniform()| and
+	 \verb|CCTK_InterpLocalNonUniform()| APIs.
+	 }  
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\section{Introduction}
+
+This thorn provides processor-local interpolation of 1, 2, and 3-D
+arrays.  At present there are 2~interpolators provided (we hope to
+add other interpolators soon):
+\begin{description}
+\item[Uniform Cartesian]
+	This is the local interpolator which used to live
+	in the PUGHInterp thorn.  It was written in C by
+	Thomas Radke in early 2001, drawing on earlier Fortran code
+	by Paul Walker.  It supports the \verb|CCTK_InterpLocal()| API.
+	It provides Lagrange polynomial interpolation for all
+	combinations of~1, 2, and 3~dimensions, and orders~1,
+	2, and~3.  (It would be easy to add additional dimensions
+	and/or orders if desired.)
+\item[Generalized Polynomial]
+	This interpolator was written in C (plus Maple to generate
+	the coefficients) by Jonathan Thornburg in winter 2001-2002.
+	It supports the \verb|CCTK_InterpLocalUniform()| API.
+	It provides Lagrange polynomial interpolation in
+	1~dimension (orders~1-6) and~2 and 3~dimensions (orders~1-4).
+	(Again, it would be easy to add additional dimensions and/or
+	orders.)  It offers a number of options, described below.
+\end{description}
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\section{Terminology}
+
+Within Cactus, each interpolator has a character-string name;
+this is mapped to a Cactus \defn{interpolator handle} by
+\verb|CCTK_InterpHandle()|.  For any given interpolator handle,
+there may be a separate interpolator defined for each of the
+interpolation APIs (both the processor-local ones provided by this
+thorn, and the global ones provided by driver-specific interpolator
+thorns such as PUGHInterp).
+
+Terminology for interpolation seems to differ a bit from one author
+to another.  Here we refer to the \defn{point-centered} interpolation
+of a set of \defn{input arrays} (defining data on a uniformly or
+nonuniformly spaced \defn{grid} of \defn{data points}) to a set of
+\defn{interpolation points} (specified by a corresponding set of
+\defn{interpolation coordinates}), with the results being stored
+in a corresponding set of \defn{output arrays}.  Alternatively,
+we may refer to \defn{cell-centered} interpolation, using a grid
+of \defn{data cells} and a set of \defn{interpolation cells}.
+(This latter terminology is common in hydrodynamics interpolators.)
+
+At present all the interpolators do polynomial interpolation, and
+allow the interpolation of multiple input arrays (to the same set of
+interpolation points) in a single interpolator call, using the basic
+algorithm:
+\begin{verbatim}
+for each interpolation point
+{
+choose a finite difference molecule position
+   somewhere near the interpolation point
+        for each input array
+        {
+        compute an interpolating polynomial which approximates
+           the input data at the molecule points
+        output = polynomial(interpolation point)
+        }
+}
+\end{verbatim}
+In the future, we may add other interpolators where the choice of
+molecule is data-dependent (and thus may vary from one input array to
+the next), and/or where the entire input grid is used in each interpolation.
+
+We define the \defn{order} of the interpolation to be the order of
+the fitted polynomial.  That is, in our terminology, linear interpolation
+is order~1, quadratic is order~2, cubic is order~3, etc.  An order~$n$
+interpolator thus has $O(\Delta x^{n+1})$ errors for smooth input data.
+
+However, because the interpolating polynomial generally changes if
+the interpolation point moves from one grid cell to another, unless
+something special is done the interpolating function isn't smooth,
+\ie{} its The 1st~derivative is generically {\em discontinuous\/},
+with $O(\Delta x^n)$ jump discontinuities each time the interpolating
+polynomial changes.  Correspondingly, the interpolation error is
+generically a ``bump function'' which is zero at each grid point
+and rises to a local maximum in each grid cell.  There are interpolation
+algorithms (\eg{} Hermite polynomial and spline interpolation) which
+give better smoothness, but none of the present interpolators implement
+them. :(
+
+As given in the function reference section of the Cactus User's Guide,
+\verb|interp_coords|, \verb|input_arrays|, and \verb|output_arrays| are
+actually all pointers to arrays of \verb|void *| pointers, since we
+support a number of different Cactus data types.  Internally, the
+interpolator casts these \verb|void *| pointers to \verb|CCTK_REAL *|
+or whatever the correct Cactus data types are.  But when describing
+how the interpolator accesses the various arrays, for simplicity we
+usually gloss over this casting, \ie{}~we pretend that \verb|interp_coords|,
+\verb|input_arrays|, and \verb|output_arrays| are pointers to arrays
+of \verb|CCTK_REAL *| pointers.  (This may become clearer once you
+read the next example.)
+
+We use \verb|pt|, \verb|in|, and \verb|out| as generic 0-origin integer
+subscripts into the arrays of interpolation points, input arrays, and
+output arrays respectively.  We use \verb|(i,j,k)| as a generic
+\verb|N_dims|-vector of integer subscripts into the input array
+\verb|input_arrays[in]|.  (Thus \defn{{\tt |(i,j,k)|} space} refers to
+the grid of data points.)  We usually only write array subscripting
+expressions for the 3-D case; the restrictions/generalizations to
+other dimensions should be obvious.
+
+For example, for 3-D interpolation, the (x,y,z) coordinates of a typical
+interpolation point are given by
+\begin{verbatim}
+x = interp_coords[0][pt]
+y = interp_coords[1][pt]
+z = interp_coords[2][pt]
+\end{verbatim}
+(Notice here that as described above, we've implicitly taken
+\verb|interp_coords| to have the C~type
+\verb|const CCTK_REAL *interp_coords[]|, glossing over the casting
+from its actual C~type of \verb|const void *interp_coords[]|.)
+
+We take \verb|axis| to be an integer specifying a dimension,
+\ie{}~0~for~$x$, 1~for~$y$, 2~for~$z$, \dots.
+
+When describing Jacobians and domains of dependence, it's useful to
+introduce the notion of \defn{molecule position}, a nominal reference
+point for the interpolation molecule in \verb|(i,j,k)| coordinate
+space.  (For example, the molecule position might just be the \verb|(i,j,k)|
+coordinates of the molecule's center.)  We also introduce
+\defn{molecule coordinates} \verb|(mi,mj,jk)|, which are just
+\verb|(i,j,k)| coordinates relative to the molecule position.
+We use \verb|m| or as a generic molecule coordinate.  Thus
+(in notation which should be obvious) a generic molecule operation
+can be written
+\begin{equation}
+\verb|output| = \sum_{\tt m} \verb|coeff[posn+m]| \times \verb|input[posn+m]|
+\end{equation}
+Note that there is no requirement that the output be semantically
+located at the position \verb|posn|!  (This may become clearer once
+you look at the example in
+section~\ref{sect-generic-options/Jacobian/fixed-sized-hypercube}.)
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\section{Generic Interpolator Options}
+\label{sect-generic-options}
+
+The newer \verb|CCTK_InterpLocalUniform()| and
+\verb|CCTK_InterpLocalNonUniform()| APIs specify a \defn{parameter table}
+(a Cactus key-value table, manipulated by the \verb|Util_Table*()| APIs)
+as a generic mechanism for passing optional input to, and/or returning
+optional results from, the interpolator.  Different interpolators support
+different options; in this section we describe some options which are,
+or will plausibly, be common to multiple interpolators.
+
+Note that except as described in section~\ref{sect-generic-options/caching}
+(``Caching''), all interpolator arguments and parameters apply only to
+the current interpolator call: there is no visible state kept inside
+the interpolator from one call to another.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\subsection{Interpolation Order}
+
+The uniform Cartesian interpolator encodes the order in the operator
+name, but other interpolators should use a (mandatory) parameter-table
+parameter
+\begin{verbatim}
+const CCTK_INT order;
+\end{verbatim}
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\subsection{Handling of Out-of-Range Interpolation Points}
+
+By default, interpolators should consider it an error if any interpolation
+point lies outside the grid, \ie{} if the ``interpolation'' is actually
+an extrapolation.  (Some interpolators may apply a small ``fuzz'' to
+this test to help avoid floating-point rounding problems.)
+
+Some interpolators may allow this behavior to be changed by the
+optional parameter
+\begin{verbatim}
+const CCTK_REAL out_of_range_tolerance[N_dims];
+\end{verbatim}
+The semantics of this are as follows:  For each \verb|axis|,
+\begin{description}
+\item[\rm If $\hbox{\tt out\_of\_range\_tolerance[axis]} \ge 0.0$,]
+	then an interpolation point is considered to be ``out of range''
+	if and only if the interpolation point is
+	$> \verb|out_of_range_tolerance[axis]|
+	   \times \verb|coord_delta[axis]|$
+	outside the grid in any coordinate.
+\item[\rm If $\hbox{\tt out\_of\_range\_tolerance[axis]} = -1.0$,]
+	then an interpolation point is considered to be ``out of range''
+	if and only if a centered molecule (or more generally,
+	a molecule whose centering is chosen pretending that
+	the grid is of infinite extent), would require data
+	from outside the grid.
+\end{description}
+Other values of \verb|out_of_range_tolerance[axis]| are illegal.
+
+To provide the ``fuzz'' noted above, \verb|out_of_range_tolerance[]|
+should default to having all elements set to a very small positive
+value, say $10^{-14}$.  (This value should really be scaled with the
+floating-point precision used,%%%
+\footnote{%%%
+	 If this scaling is done, a reasonable value for
+	 {\tt out\_of\_range\_tolerance[]} would be around
+	 $100 \epsilon$, where the \defn{machine epsilon}
+	 $\epsilon$ is defined to be the smallest positive
+	 floating-point number such that $1.0 + \epsilon$
+	 compares ``not equal to'' 1.0 in floating-point
+	 arithmetic.  Machine epsilon values can be found
+	 in the standard header file {\tt <float.h>};
+	 for IEEE single and double precision they are
+	 about $1.19{\times}10^{-7}$ and $2.22{\times}10^{-16}$
+	 respectively.
+	 }%%%
+{} but none of the interpolators do this at present.)
+
+If any interpolation points are out of range, interpolators should
+return the error code \verb|CCTK_ERROR_INTERP_POINT_X_RANGE|, and
+report further details about the error by setting the following
+parameter-table entries:
+\begin{verbatim}
+/* which interpolation point is out of range? */
+/* (value is  pt  for out-of-range point) */
+CCTK_INT out_of_range_pt;
+
+/* in which coordinate axis is the point out of range? */
+/* (value is  axis  for out-of-range point) */
+CCTK_INT out_of_range_axis;
+
+/* on which end of the grid is the point out of range? */
+/* (value is -1 for min, +1 for max) */
+CCTK_INT out_of_range_end;
+\end{verbatim}
+
+Note that if multiple points and/or axes are out of range, different
+interpolators may vary in which out-of-range error they report.  That
+is, it is {\em not\/} necessarily the case that the ``first'' such
+error will be the one reported.%%%
+\footnote{%%%
+	 This is to make life simpler for interpolators
+	 which work in parallel over multiple processors.
+	 }%%%
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\subsection{Molecule Family}
+\label{sect-generic-options/molecule-family}
+
+An interpolator may support different families/shapes of interpolation
+molecules.  Hypercube-shaped molecules are the simplest and most common
+case, but one could also imagine (say) octagon-shaped molecules in 2-D,
+or some generalization of this in higher numbers of dimensions.
+
+The parameter
+\begin{verbatim}
+const char molecule_family[];   /* set with Util_TableSetString() */
+\end{verbatim}
+may be used to query or change the strategy for choosing the molecules.
+
+The semantics are that if this key is present with the value \verb|""|
+(a 0-length ``empty'' string), this queries the default molecule family:
+the interpolator will (re)set the value to a string describing the default
+molecule family (\verb|"cube"| for hypercube-shaped molecules%%%
+\footnote{%%%
+	 Strictly speaking, these should be called
+	 ``parallelipiped-shaped'' molecules, but this
+	 term is so clumsy (and hard to spell!) that
+	 we just call them hypercube-shaped instead.
+	 }%%%
+).
+If this key is present with a value which is a non-empty string, this
+sets the molecule family/shape based on that string.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\subsection{Non-Contiguous Input Arrays}
+\label{sect-generic-options/non-contiguous-inputs}
+
+Sometimes one of the input ``arrays'' used by the interpolator might not
+be contiguous in memory.  For example, we might want to do 2-D interpolation
+within a plane of a 3-D grid array, but the plane might or might not be
+contiguous in memory.  (Another example would be that the input array
+might be a member of a compact group.)
+
+One way to do this would be to use the hyperslab API to explicitly compute
+the input data on an appropriate hyperslab, then interpolate within that.
+However, some interpolators may support accessing the appropriate hyperslab
+of the input grid directly inside the interpolator.  If this is supported,
+it's probably considerably cheaper than explicitly computing the hyperslab.
+
+If an interpolator supports this, it should use the following (optional)
+parameter-table entries to specify non-contiguous inputs:
+\begin{verbatim}
+const CCTK_INT input_array_offsets[N_input_arrays];
+
+/* the next 3 table entries are shared by all input arrays */
+const CCTK_INT input_array_strides       [N_dims];
+const CCTK_INT input_array_min_subscripts[N_dims];
+const CCTK_INT input_array_max_subscripts[N_dims];
+\end{verbatim}
+
+Then the interpolator would access the input using the generic
+subscripting expression
+\begin{verbatim}
+input_array[in][offset + i*stride_i + j*stride_j + k*stride_k]
+\end{verbatim}
+where
+\begin{verbatim}
+offset = input_array_offsets[in]
+(stride_i,stride_j,stride_k) = input_array_strides[]
+\end{verbatim}
+\verb|(i,j,k)| run from \verb|input_array_min_subscripts[]|
+to \verb|input_array_max_subscripts[]| inclusive
+(n.b.~this is an {\em inclusive\/} range, \ie{}
+$\verb|min| \le \verb|(i,j,k)| \le \verb|max|$).
+
+The defaults are that each input array is contiguous in memory,
+\ie{} \verb|input_array_offsets[]| = 0,
+\verb|stride| determined from \verb|input_array_dims[]|
+              in the usual Fortran manner,
+\verb|input_array_min_subscripts[]| = all 0s, and
+\verb|input_array_max_subscripts[]| = \verb|input_array_dims[]|-1.
+If the stride and max subscript are both specified explicitly, then the
+explicit \verb|input_array_dims[]| argument to
+\verb|CCTK_InterpLocalUniform()| is ignored.
+
+At present all the interpolators take the output arrays to be
+contiguous 1-D arrays.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\subsection{Derivatives}
+
+Some interpolators can optionally take derivatives as part of the
+interpolation, \ie{} if we view the input data as being samples of
+a smooth function, then instead of estimating values of that function
+at the interpolation points, the interpolator can instead or additionally
+estimate values of various derivatives of that function at the
+interpolation points.  We don't currently implement it, but one could
+also imagine interpolating more general molecule operations such as
+Laplacians, running integrals, etc.
+
+To specify such operations, an interpolator should use the parameter-table
+entries
+\begin{verbatim}
+const CCTK_INT operand_indices[N_output_arrays];
+const CCTK_INT operation_codes[N_output_arrays];
+\end{verbatim}
+The semantics here are that
+\begin{verbatim}
+output_array[out] = op(input_array[operand_indices[out]])
+\end{verbatim}
+where \verb|op| is an operator specified by the \verb|operation_codes[out]|
+value as described below.
+
+Note that \verb|operand_indices[]| doesn't directly name the inputs,
+but rather gives their (0-origin) subscripts in the list of inputs.
+This allows for a more efficient implementation in the (common) case
+where some of the input arrays have many different operations applied
+to them.
+
+Negative \verb|operation_codes[out]| values  are reserved for future
+expansion.  An \verb|operation_codes[out]| value which is $\ge 0$ is
+taken as a decimal integer encoding a coordinate partial derivative:
+each decimal digit means to take the coordinate partial derivative along
+that (1-origin) axis; the order of the digits in a number is ignored.
+For example:
+\begin{flushleft}
+\begin{tabular}{cl}
+\verb|operation_codes[out]|
+	& What it means								\\
+\hline %-----------------------------------------------------------------------
+0	& no derivative, ``just'' interpolate the input array			\\
+1	& interpolate $\partial \big/ \partial x^1$ of the input array		\\
+2	& interpolate $\partial \big/ \partial x^2$ of the input array		\\
+12 or 21
+	& interpolate $\partial^2 \big/ \partial x^1 \partial x^2$
+	  of the input array							\\
+33	& interpolate $\partial^2 \big/ \partial (x^3)^2$ of the input array	%%%\\
+\end{tabular}
+\end{flushleft}
+
+At present we do {\em not\/} have any \verb|#define|s for the
+operation codes in the Cactus header files.  However, you can avoid
+most of the software-engineering problems of having ``magic numbers''
+for the operation codes, by using the macro
+\begin{verbatim}
+#define DERIV(op)   op
+\end{verbatim}
+to mark all such \verb|operation_codes[]| values in your code.
+There's an example of this in section~\ref{sect-example-derivatives}.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\subsection{Jacobian and Domain of Dependence}
+\label{sect-generic-options/Jacobian}
+
+Sometimes we want to know the Jacobian of the interpolator,
+\begin{equation}
+\frac{\partial \hbox{\tt output\_array[out][pt]}}
+     {\partial \hbox{\tt input\_array[in][(i,j,k)]}}
+							\label{eqn-Jacobian}
+\end{equation}
+and/or ``just'' the set of \verb|(i,j,k)| for which this is nonzero
+(\ie{} the sparsity structure of the Jacobian, equivalently the domain
+of dependence of the output result, or equivalently the interpolation
+molecule size and shape) for a given \verb|out|, \verb|in|, and \verb|pt|.
+
+The complexity of doing this depends (strongly!) on the structure
+of the Jacobian, and in particular on the answers to the following
+questions:
+\begin{itemize}
+\item	What molecule family is being used?%%%
+\footnote{%%%
+	 For something like a spline interpolator the Jacobian
+	 is generally nonzero for all {\tt |(i,j,k)|}, but
+	 exponentially small for most of them; in this case
+	 the Jacobian-query API would probably specify a cutoff
+	 to return an approximate Jacobian with reasonable sparsity.
+	 }%%%
+\item	Does the interpolator always use the same molecule size
+	and shape (possibly depending on the interpolation order,
+	molecule family, or other such parameters), independent
+	of where the interpolation points are in the grid?%%%
+\footnote{%%%
+	 We can always make the interpolation molecules be
+	 the ``same size and shape'' by padding them with
+	 zero entries to bring them all up to the size of
+	 the largest molecule used anywhere in the grid,
+	 but this would be very inefficient if many molecules
+	 were padded in this way.
+	 }%%%
+\item	Does the interpolator use the same interpolation molecule
+	size and shape for each output array?  (The answer to this
+	question may depend on whether and/or what derivatives are
+	being computed.)
+\end{itemize}
+Because the different cases differ so much in their complexity,
+we define several distinct APIs for obtaining the interpolator's
+Jacobian and/or domain of dependence.
+
+%%%%%%%%%%%%%%%%%%%%
+
+\subsubsection{Querying the Interpolator about the Jacobian's structure}
+
+To allow generic code to determine which of the different Jacobian-structure
+cases applies, (and thus which APIs to use), an interpolator which
+supports Jacobian operations should support using the following
+parameter-table entries to query the interpolator about the Jacobian
+structure:
+
+The parameter \verb|molecule_family| may may be used to query what
+molecule family is being used.  This is described in detail in
+section~\ref{sect-generic-options/molecule-family}.
+
+This parameter may be used to query whether the interpolation molecule's
+size and/or shape varies with the interpolation coordinates:
+\begin{verbatim}
+CCTK_INT MSS_is_fn_of_interp_coords;
+\end{verbatim}
+The semantics of this are that if this key is present%%%
+\footnote{%%%
+	 ``MSS'' abbreviates ``{\bf m}olecule {\bf s}ize
+	 and/or {\bf s}hape''.
+	 }%%%
+{} (the value doesn't matter), then the interpolator will (re)set
+the value to 0~if the molecule size and shape do {\em not\/} vary
+with the interpolation coordinates, or 1~if the molecule size and/or
+shape {\em do\/} vary with the interpolation coordinates.
+
+This parameter may be used to query whether the interpolation molecule's
+size and/or shape varies from one output array to another:
+\begin{verbatim}
+CCTK_INT MSS_is_fn_of_which_output;
+\end{verbatim}
+The semantics of this are that if this key is present (the value
+doesn't matter), then the interpolator will (re)set the value to
+0~if the molecule size and shape do {\em not\/} vary from one output
+array to another (this includes the case where there is only 1~output
+array!), or 1~if the molecule size and/or shape {\em do\/} vary from
+one output array to another.  Note that since the answer to this may
+depend on whether and/or what derivatives are being computed, setting
+this key may force the interpolator to scan completely through the
+derivative specifications.  This probably isn't {\em very\/} costly,
+but you may want to avoid unnecessarily paying that cost by deleting
+the query key (\verb|MSS_is_fn_of_which_output|) from the table
+on future interpolator calls.
+
+%%%%%%%%%%%%%%%%%%%%
+
+\subsubsection{Fixed-Size Hypercube-Shaped Molecules}
+\label{sect-generic-options/Jacobian/fixed-sized-hypercube}
+
+The simplest case (and the only one for which we have defined an API
+at present) is when the molecules are hypercube-shaped and of (small)
+fixed size, independent of the interpolation coordinates (though likely
+depending on the interpolation order or other such parameters), and
+not varying from one output array to another.
+
+The following parameters may be used to query the molecule size:
+\begin{verbatim}
+CCTK_INT const interp_molecule_min_m[N_dims];
+CCTK_INT const interp_molecule_max_m[N_dims];
+\end{verbatim}
+The semantics of these are that if these keys are present (the values
+don't matter), then the interpolator will (re)set the values to give
+the (inclusive) minimum and maximum \verb|m|~molecule coordinates.
+
+The following parameter may be used to query the molecule positions:
+\begin{verbatim}
+CCTK_INT *const interp_molecule_positions[N_dims];
+\end{verbatim}
+The semantics of this is that the caller should set
+\verb|interp_molecule_positions[]| to an array of
+\verb|N_dims| pointers to (caller-supplied) arrays of
+\verb|N_interp_points| \verb|CCTK_INT|s each.  If this key exists,
+then the interpolator will store the molecule positions in the
+pointed-to arrays.
+
+The following parameters may be used to query the
+Jacobian~\eqref{eqn-Jacobian} itself:
+\begin{verbatim}
+CCTK_REAL *const Jacobian_pointer[N_output_arrays];
+const CCTK_INT   Jacobian_offset [N_output_arrays];
+
+/* the next 2 table entries are shared by all Jacobians */
+const CCTK_INT Jacobian_interp_point_stride;
+const CCTK_INT Jacobian_m_strides[N_dims];
+\end{verbatim}
+Then the interpolator would store the Jacobian~\eqref{eqn-Jacobian} in
+\begin{verbatim}
+Jacobian_pointer[out][offset + pt*Jacobian_interp_point_stride
+                             + mi*stride_i + mj*stride_j + mk*stride_k]
+\end{verbatim}
+where
+\begin{verbatim}
+offset = Jacobian_offset[out]
+(stride_i,stride_j,stride_k) = Jacobian_m_strides[]
+\end{verbatim}
+
+By appropriately setting the various stride parameters, this allows
+a fairly wide variety of possible storage layouts for the Jacobian.
+
+An example may help to clarify this:  Suppose we have a 1-D grid
+with 11~grid points, with integer subscripts~0 through~10 inclusive,
+and interpolation coordinates given by \verb|coord_origin = 0.0|
+and \verb|coord_delta = 0.1|.
+Suppose further that we're doing quadratic interpolation, using
+3-point (vertex-centered) molecules, and that the interpolator
+centers the interpolation molecules as close to the interpolation
+point as possible subject to the constraint that the molecules
+never require data from outside the grid.
+Finally, suppose that we query the Jacobian molecule positions for
+the \verb|N_interp_points=14| interpolation points 0.0, 0.04, 0.06,
+0.10, 0.14, 0.16, 0.20, 0.80, 0.84, 0.86, 0.90, 0.94, 0.96, 1.00.
+Then the queries might return
+\begin{verbatim}
+interp_molecule_min_m = -1
+interp_molecule_max_m = +1
+                                         /* interp_x      molecule     */
+interp_molecule_positions[0][ 0] =  1    /*   0.00     [0.0, 0.1, 0.2] */
+interp_molecule_positions[0][ 1] =  1    /*   0.04     [0.0, 0.1, 0.2] */
+interp_molecule_positions[0][ 2] =  1    /*   0.06     [0.0, 0.1, 0.2] */
+interp_molecule_positions[0][ 3] =  1    /*   0.10     [0.0, 0.1, 0.2] */
+interp_molecule_positions[0][ 4] =  1    /*   0.14     [0.0, 0.1, 0.2] */
+interp_molecule_positions[0][ 5] =  2    /*   0.16     [0.1, 0.2, 0.3] */
+interp_molecule_positions[0][ 6] =  2    /*   0.20     [0.1, 0.2, 0.3] */
+interp_molecule_positions[0][ 7] =  8    /*   0.80     [0.7, 0.8, 0.9] */
+interp_molecule_positions[0][ 8] =  8    /*   0.84     [0.7, 0.8, 0.9] */
+interp_molecule_positions[0][ 9] =  9    /*   0.86     [0.8, 0.9, 1.0] */
+interp_molecule_positions[0][10] =  9    /*   0.90     [0.8, 0.9, 1.0] */
+interp_molecule_positions[0][11] =  9    /*   0.94     [0.8, 0.9, 1.0] */
+interp_molecule_positions[0][12] =  9    /*   0.96     [0.8, 0.9, 1.0] */
+interp_molecule_positions[0][13] =  9    /*   1.00     [0.8, 0.9, 1.0] */
+\end{verbatim}
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\subsection{Caching}
+\label{sect-generic-options/caching}
+
+Some interpolators may support special ``caching'' optimizations to
+speed repeated interpolations where some or all of the interpolator
+arguments and/or parameters are the same.  For example, when interpolating
+a tabulated equation of state the number of dimensions, the coordinate
+origin and grid spacing, and the input arrays (the tabulated equation
+of state data), will probably all be the same from one interpolator
+call to another.
+
+If an interpolator supports caching, the following parameters should
+be used to control this:
+
+\begin{verbatim}
+const char cache_type[];        /* set with Util_TableSetString() */
+const char cache_control[];     /* set with Util_TableSetString() */
+CCTK_INT cache_handle;
+\end{verbatim}
+
+There are three basic operations supported:
+\begin{description}
+\item[Create a Cache]
+	To set up caching, call the interpolator with \verb|cache_type|
+	set to describe what arguments and/or parameters will remain
+	the same in future interpolator calls, \verb|cache_control|
+	set to the string \verb|"create"|, and \verb|cache_handle|
+	{\em not\/} in the parameter table.  The interpolator will
+	then do extra (possibly quite time-consuming) work to set
+	up cached information.  The interpolator will delete the
+	key \verb|cache_control|, and return a handle to the cached
+	information in \verb|cache_handle|; this allows multiple caches
+	to be active concurrently.
+\item[Use a Cache]
+	To use a cache (\ie{}~to make an interpolation with the
+	hoped-for speedup), just call the interpolator with
+	\verb|cache_handle| set to the value returned when the cache
+	was created.  Note that you still have to provide all the
+	``will be the same'' interpolator arguments and/or parameters;
+	providing a cache handle is essentially just a promise that
+	these will be the same as in the cache-create interpolator
+	call.  The details of what information is cached, and if/how
+	the ``will be the same'' arguments are still used, are up to
+	the interpolator.
+\item[Destroy a Cache]
+	To destroy a cache (\ie{} free any memory allocated when
+	the cache was created), call the interpolator with
+	\verb|cache_handle| set to the value returned when the cache
+	was created, and  \verb|cache_control| set to the string
+	\verb|"destroy"|.  The interpolator will delete the keys
+	\verb|cache_handle| and \verb|cache_control|, and destroy
+	the cache.
+\end{description}
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\section{The Uniform Cartesian Interpolator}
+
+The uniform Cartesian interpolator is the one which used to live
+in the PUGHInterp thorn.  It was written in C by Thomas Radke in
+early 2001, drawing on earlier Fortran code by Paul Walker.  It
+supports the \verb|CCTK_InterpLocal()| API, and thus doesn't use
+any parameter table at all.  It provides Lagrange polynomial
+interpolation of orders 1, 2, and 3, registered under the operator
+names \verb|"first-order uniform Cartesian"|,
+\verb|"second-order uniform Cartesian"|, and
+\verb|"third-order uniform Cartesian"| respectively.  Each of these
+supports 1, 2, and 3-D arrays.  The code is hard-wired for
+hypercube-shaped interpolation molecules, but it would be easy
+to add additional orders and/or dimensions if desired.
+
+Although the \verb|CCTK_InterpLocal()| API supports both uniform
+and nonuniform grids for the input data, the present implementation
+assumes a uniform grid (and silently gives wrong results for a
+nonuniform grid).
+
+See the Cactus User's Guide ``Full Description of Functions''
+appendix for a full description of the \verb|CCTK_InterpLocal()|
+API, and some examples of how to use it.
+
+%%%%%%%%%%%%%%%%%%%%
+
+\subsection{Implementation}
+
+Internally, this interpolator always does 3-D interpolation, inserting
+zero coordinates as appropriate for lower dimensionalities.  The
+interpolation is done by successive 1-D interpolations along each
+axis.  See the \verb|README| file in the source code directory
+\verb|LocalInterp/src/UniformCartesian/| for further details.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\section{The Generalized Polynomial Interpolator}
+
+The generalized polynomial interpolator was written in C
+(plus Maple to generate the coefficients) by Jonathan Thornburg in
+winter 2001-2002.  It provides Lagrange polynomial interpolation of
+orders 1-6 for 1-D arrays, and 1-4 for 2- and 3-D
+arrays, all registered under the operator name
+\verb|"generalized polynomial"|.  (Again, it would be easy to add
+additional orders and/or dimensions if desired.)  The code allows
+arbitrarily-shaped interpolation molecules, but at the moment only
+hypercube-shaped molecules are implemented.
+
+This interpolator supports a number of the options described in
+section~\ref{sect-generic-options}:
+\begin{itemize}
+\item	interpolation order
+\item	handling of out-of-range interpolation points
+	(if there are multiple out-of-range points/axes, the
+	one reported will be the first, \ie{}~the out-of-range
+	point with the smallest \verb|pt|, and of that points
+	out-of-range axes, the one with the smallest \verb|axis|)
+	
+\item	non-contiguous input arrays
+\item	derivatives
+	(when taking derivatives with this interpolator,
+	it's most efficient to group all operations on
+	a given input array together in the \verb|operand_indices|
+	and \verb|operation_codes| arrays, as in the example in
+	section~\ref{sect-example-derivatives})
+\end{itemize}
+At present this interpolator does not support computing the Jacobian,
+but we hope to add this in the near future.  This interpolator also
+does not support any caching features; at present (Feb 2002) we have
+no plans to add these.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\subsection{Smoothing}
+
+The way the generalized polynomial interpolator is implemented it's
+easy to also do Savitzky-Golay smoothing.%%%
+\footnote{%%%
+	 See section 14.8 of Numerical Recipes
+	 2nd edition for a general discussion of
+	 Savitzky-Golay smoothing.
+	 }%%%
+{}  This is best described by way of an example:  Suppose we're doing
+1-D cubic interpolation, \ie{} at each interpolation point we're using
+a cubic polynomial fitted to 4~surrounding data points.  For
+Savitzky-Golay smoothing, we would instead {\em least-squares fit\/}
+a cubic polynomial to some {\em larger\/} number of surrounding data
+points.  This combines interpolation with smoothing, so there's less
+amplification of noise in the input data in the interpolation outputs.
+
+The optional input
+\begin{verbatim}
+const CCTK_INT smoothing;
+\end{verbatim}
+specifies how much (how many points) to enlarge the interpolation
+molecule for this.  The default is 0 (no smoothing).  1 would mean to
+enlarge the molecule by 1 point, \eg{} to use a 5-point molecule instead
+of the usual 4-point one for cubic interpolation.  2 would mean to
+enlarge by 2 points, \eg{} to use a 6-point molecule for cubic
+interpolation.  Etc etc.
+
+Note that in $>1$~dimension, the default hypercube-shaped molecules
+already use more data points than the number of free parameters in
+the interpolating polynomials, \ie{} they already do some Savitzky-Golay
+smoothing.  For example, in 2-D a generic cubic polynomial
+$f(x,y) = \sum_{i+j \leq 3} c_{ij} x^i y^j$ has 10 free parameters
+$c_{ij}$, which we least-squares fit to the 16 data points in the
+$4 \times 4$ molecule.
+
+Savitzky-Golay smoothing is basically free apart from the increase in
+the molecule size, e.g. a \verb|smoothing|=2 cubic interpolation has
+exactly the same cost as any other 6-point--molecule interpolation.
+
+Alas, at the moment only the trivial case \verb|smoothing|=0 is
+implemented, but the framework is all there for more general cases.
+
+%%%%%%%%%%%%%%%%%%%%
+
+\subsection{Implementation}
+
+This interpolator's basic design is to use separate specialized code
+for each combination of
+$(\verb|N_dims|, \verb|molecule_family|, \verb|order|, \verb|smoothing|)$,
+\ie{}~in practice for each distinct choice of interpolation molecule.
+Maple is used to generate all the interpolation coefficients.
+The C preprocessor is then used to generate all the specialized code
+from a single master ``template'' file.  The template code uses
+\verb|#ifdef|s to handle lower dimensions with no extra overhead,
+\eg{}~1-D/2-D interpolation is basically just as efficient as in
+a purpose-built 1-D/2-D interpolator, with no extra overhead imposed
+by the interpolator also supporting higher-dimensional interpolation.
+
+The Maple code which generates the interpolation coefficients is quite
+general-purpose, and can handle an arbitrary dimensionality and molecule
+size/shape.  Generating new coefficients can be rather time-consuming,
+though, \eg{}~the current coefficients for 3-D for orders~1-4 take about
+8~cpu minutes to generate using Maple~7 on a 1.7~GHz~P4.
+
+See the \verb|README| file in the source code directory
+\verb|LocalInterp/src/UniformCartesian/| for further details on the
+implementation.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\subsection{A Simple Example of {\tt CCTK\_InterpLocalUniform} Usage}
+
+Here's a simple example of interpolating a \verb|CCTK_REAL| and a
+\verb|CCTK_COMPLEX| $10 \times 20$ 2-D array, at 5 interpolation points,
+using cubic interpolation.
+
+\begin{verbatim}
+#define N_DIMS   2
+#define N_INTERP_POINTS   5
+#define N_INPUT_ARRAYS    2
+#define N_OUTPUT_ARRAYS   2
+
+/* (x,y,z) coordinates of data grid points */
+#define X_ORIGIN   ...
+#define X_DELTA    ...
+#define Y_ORIGIN   ...
+#define Y_DELTA    ...
+const CCTK_REAL origin[N_DIMS] = { X_ORIGIN, Y_ORIGIN };
+const CCTK_REAL delta [N_DIMS] = { X_DELTA,  Y_DELTA  };
+
+/* (x,y,z) coordinates of interpolation points */
+const CCTK_REAL interp_x[N_INTERP_POINTS];
+const CCTK_REAL interp_y[N_INTERP_POINTS];
+const void *const interp_coords[N_DIMS]
+        = { (const void *) interp_x, (const void *) interp_y };
+
+/* input arrays */
+/* ... note Cactus uses Fortran storage ordering, i.e. X is contiguous */
+#define N_X   10
+#define N_Y   20
+const CCTK_REAL    input_real   [N_Y][N_X];
+const CCTK_COMPLEX input_complex[N_Y][N_X];
+const CCTK_INT input_array_dims[N_DIMS] = { N_X, N_Y };
+const CCTK_INT input_array_type_codes[N_INPUT_ARRAYS]
+        = { CCTK_VARIABLE_REAL, CCTK_VARIABLE_COMPLEX };
+const void *const input_arrays[N_INPUT_ARRAYS]
+        = { (const void *) input_real, (const void *) input_complex };
+
+/* output arrays */
+CCTK_REAL    output_real   [N_INTERP_POINTS];
+CCTK_COMPLEX output_complex[N_INTERP_POINTS];
+const CCTK_INT output_array_type_codes[N_OUTPUT_ARRAYS]
+        = { CCTK_VARIABLE_REAL, CCTK_VARIABLE_COMPLEX };
+void *const output_arrays[N_OUTPUT_ARRAYS]
+        = { (void *) output_real, (void *) output_complex };
+
+int operator_handle, param_table_handle;
+operator_handle = CCTK_InterpHandle("my interpolation operator");
+if (operator_handle < 0)
+        CCTK_WARN(-1, "can't get interpolation handle!");
+param_table_handle = Util_TableCreateFromString("order=3");
+if (param_table_handle < 0)
+        CCTK_WARN(-1, "can't create parameter table!");
+
+if (CCTK_InterpLocalUniform(N_DIMS,
+                            operator_handle, param_table_handle,
+                            origin, delta,
+                            N_INTERP_POINTS,
+                               CCTK_VARIABLE_REAL,
+                               interp_coords,
+                            N_INPUT_ARRAYS,
+                               input_array_dims,
+                               input_array_type_codes,
+                               input_arrays,
+                            N_OUTPUT_ARRAYS,
+                               output_array_type_codes,
+                               output_arrays) < 0)
+        CCTK_WARN(-1, "error return from interpolator!");
+\end{verbatim}
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\subsection{An Example of Interplating Derivatives}
+\label{sect-example-derivatives}
+
+Consider the problem described earlier: the computation of all the
+2nd~derivatives of the 3-metric at a set of interpolation points on
+a 2-sphere.  Here's how we could do this in C:
+\begin{verbatim}
+#define N_DIMS        3
+
+/* interpolation points */
+#define N_INTERP_POINTS 1000
+const CCTK_REAL interp_x[N_INTERP_POINTS],
+                interp_y[N_INTERP_POINTS],
+                interp_z[N_INTERP_POINTS];
+const void *const interp_coords[N_DIMS]
+        = { (const void *) interp_x,
+            (const void *) interp_y,
+            (const void *) interp_z };
+
+/* dimensions of the data grid */
+#define NX   30
+#define NY   40
+#define NZ   50
+
+/* input arrays */
+#define N_INPUT_ARRAYS  6
+const CCTK_REAL gxx[NX][NY][NZ], gxy[NX][NY][NZ], gxz[NX][NY][NZ],
+                                 gyy[NX][NY][NZ], gyz[NX][NY][NZ],
+                                                  gzz[NX][NY][NZ];
+
+const CCTK_INT input_array_dims[N_DIMS] = {NX, NY, NZ};
+const CCTK_INT input_array_type_codes[N_INPUT_ARRAYS]
+        = { CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+                                CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+                                                    CCTK_VARIABLE_REAL };
+const void *const interp_coords[N_INPUT_ARRAYS]
+        = { (const void *) gxx, (const void *) gxy, (const void *) gxz,
+                                (const void *) gyy, (const void *) gyz,
+                                                    (const void *) gzz };
+
+/* output arrays */
+#define N_OUTPUT_ARRAYS 36
+CCTK_REAL
+ dxx_gxx[N_INTERP_POINTS], dxy_gxx[N_INTERP_POINTS], dxz_gxx[N_INTERP_POINTS],
+                           dyy_gxx[N_INTERP_POINTS], dyz_gxx[N_INTERP_POINTS],
+                                                     dzz_gxx[N_INTERP_POINTS],
+ dxx_gxy[N_INTERP_POINTS], dxy_gxy[N_INTERP_POINTS], dxz_gxy[N_INTERP_POINTS],
+                           dyy_gxy[N_INTERP_POINTS], dyz_gxy[N_INTERP_POINTS],
+                                                     dzz_gxy[N_INTERP_POINTS],
+ dxx_gxz[N_INTERP_POINTS], dxy_gxz[N_INTERP_POINTS], dxz_gxz[N_INTERP_POINTS],
+                           dyy_gxz[N_INTERP_POINTS], dyz_gxz[N_INTERP_POINTS],
+                                                     dzz_gxz[N_INTERP_POINTS],
+ dxx_gyy[N_INTERP_POINTS], dxy_gyy[N_INTERP_POINTS], dxz_gyy[N_INTERP_POINTS],
+                           dyy_gyy[N_INTERP_POINTS], dyz_gyy[N_INTERP_POINTS],
+                                                     dzz_gyy[N_INTERP_POINTS],
+ dxx_gyz[N_INTERP_POINTS], dxy_gyz[N_INTERP_POINTS], dxz_gyz[N_INTERP_POINTS],
+                           dyy_gyz[N_INTERP_POINTS], dyz_gyz[N_INTERP_POINTS],
+                                                     dzz_gyz[N_INTERP_POINTS],
+ dxx_gzz[N_INTERP_POINTS], dxy_gzz[N_INTERP_POINTS], dxz_gzz[N_INTERP_POINTS],
+                           dyy_gzz[N_INTERP_POINTS], dyz_gzz[N_INTERP_POINTS],
+                                                     dzz_gzz[N_INTERP_POINTS];
+const CCTK_INT output_array_type_codes[N_OUTPUT_ARRAYS]
+   = { CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+                           CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+                                               CCTK_VARIABLE_REAL,
+       CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+                           CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+                                               CCTK_VARIABLE_REAL,
+       CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+                           CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+                                               CCTK_VARIABLE_REAL,
+       CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+                           CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+                                               CCTK_VARIABLE_REAL,
+       CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+                           CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+                                               CCTK_VARIABLE_REAL,
+       CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+                           CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+                                               CCTK_VARIABLE_REAL };
+void *const output_arrays[N_OUTPUT_ARRAYS]
+   = { (void *) dxx_gxx, (void *) dxy_gxx, (void *) dxz_gxx,
+                         (void *) dyy_gxx, (void *) dyz_gxx,
+                                           (void *) dzz_gxx,
+       (void *) dxx_gxy, (void *) dxy_gxy, (void *) dxz_gxy,
+                         (void *) dyy_gxy, (void *) dyz_gxy,
+                                           (void *) dzz_gxy,
+       (void *) dxx_gxz, (void *) dxy_gxz, (void *) dxz_gxz,
+                         (void *) dyy_gxz, (void *) dyz_gxz,
+                                           (void *) dzz_gxz,
+       (void *) dxx_gyy, (void *) dxy_gyy, (void *) dxz_gyy,
+                         (void *) dyy_gyy, (void *) dyz_gyy,
+                                           (void *) dzz_gyy,
+       (void *) dxx_gyz, (void *) dxy_gyz, (void *) dxz_gyz,
+                         (void *) dyy_gyz, (void *) dyz_gyz,
+                                           (void *) dzz_gyz,
+       (void *) dxx_gzz, (void *) dxy_gzz, (void *) dxz_gzz,
+                         (void *) dyy_gzz, (void *) dyz_gzz,
+                                           (void *) dzz_gzz };
+
+/* integer codes to specify the derivatives */
+/* (for best efficiency we group all operations on a given input together) */
+const CCTK_INT operand_indices[N_OUTPUT_ARRAYS]
+   = { 0, 0, 0, 0, 0, 0,
+       1, 1, 1, 1, 1, 1,
+       2, 2, 2, 2, 2, 2,
+       3, 3, 3, 3, 3, 3,
+       4, 4, 4, 4, 4, 4,
+       5, 5, 5, 5, 5, 5 };
+#define DERIV(x)   x
+const CCTK_INT operation_codes[N_OUTPUT_ARRAYS]
+   = { DERIV(11), DERIV(12), DERIV(13), DERIV(22), DERIV(23), DERIV(33),
+       DERIV(11), DERIV(12), DERIV(13), DERIV(22), DERIV(23), DERIV(33),
+       DERIV(11), DERIV(12), DERIV(13), DERIV(22), DERIV(23), DERIV(33),
+       DERIV(11), DERIV(12), DERIV(13), DERIV(22), DERIV(23), DERIV(33),
+       DERIV(11), DERIV(12), DERIV(13), DERIV(22), DERIV(23), DERIV(33),
+       DERIV(11), DERIV(12), DERIV(13), DERIV(22), DERIV(23), DERIV(33) };
+
+int operator_handle, param_table_handle;
+operator_handle = CCTK_InterpHandle("my interpolation operator");
+if (operator_handle < 0)
+        CCTK_WARN(-1, "can't get interpolation handle!");
+
+param_table_handle = Util_TableCreate(UTIL_TABLE_FLAGS_DEFAULT);
+if (param_table_handle < 0)
+        CCTK_WARN(-1, "can't create parameter table!");
+if (Util_TableSetInt(param_table_handle, 3, "order") < 0)
+        CCTK_WARN(-1, "can't set order in parameter table!");
+if (Util_TableSetIntArray(param_table_handle,
+                          N_OUTPUT_ARRAYS, operand_indices,
+                          "operand_indices") < 0)
+        CCTK_WARN(-1, "can't set operand_indices array in parameter table!");
+if (Util_TableSetIntArray(param_table_handle,
+                          N_OUTPUT_ARRAYS, operation_codes,
+                          "operation_codes") < 0)
+        CCTK_WARN(-1, "can't set operation_codes array in parameter table!");
+
+if (CCTK_InterpLocalUniform(N_DIMS,
+                            operator_handle, param_table_handle,
+                            origin, delta,
+                            N_INTERP_POINTS,
+                               CCTK_VARIABLE_REAL,
+                               interp_coords,
+                            N_INPUT_ARRAYS,
+                               input_array_dims,
+                               input_array_type_codes,
+                               input_arrays,
+                            N_OUTPUT_ARRAYS,
+                               output_array_type_codes,
+                               output_arrays) < 0)
+        CCTK_WARN(-1, "error return from interpolator!");
+\end{verbatim}
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%
+% Automatically created from the ccl files 
+% Do not worry for now.
+%
+\include{interface}
+\include{param}
+\include{schedule}
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\section{Acknowledgments}
+
+Thanks to Ian Hawke for helpful comments on the documentation,
+to Tom Goodale and Thomas Radke for many useful design discussions,
+and to all the Cactus crew for a great infrastructure!
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\end{document}
diff --git a/doc/interface.aux b/doc/interface.aux
new file mode 100644
index 0000000..76caf46
--- /dev/null
+++ b/doc/interface.aux
@@ -0,0 +1,19 @@
+\relax 
+\@setckpt{interface}{
+\setcounter{page}{18}
+\setcounter{equation}{2}
+\setcounter{enumi}{0}
+\setcounter{enumii}{0}
+\setcounter{enumiii}{0}
+\setcounter{enumiv}{0}
+\setcounter{footnote}{7}
+\setcounter{mpfootnote}{0}
+\setcounter{part}{0}
+\setcounter{section}{5}
+\setcounter{subsection}{4}
+\setcounter{subsubsection}{0}
+\setcounter{paragraph}{0}
+\setcounter{subparagraph}{0}
+\setcounter{figure}{0}
+\setcounter{table}{0}
+}
diff --git a/doc/param.aux b/doc/param.aux
new file mode 100644
index 0000000..4164907
--- /dev/null
+++ b/doc/param.aux
@@ -0,0 +1,19 @@
+\relax 
+\@setckpt{param}{
+\setcounter{page}{18}
+\setcounter{equation}{2}
+\setcounter{enumi}{0}
+\setcounter{enumii}{0}
+\setcounter{enumiii}{0}
+\setcounter{enumiv}{0}
+\setcounter{footnote}{7}
+\setcounter{mpfootnote}{0}
+\setcounter{part}{0}
+\setcounter{section}{5}
+\setcounter{subsection}{4}
+\setcounter{subsubsection}{0}
+\setcounter{paragraph}{0}
+\setcounter{subparagraph}{0}
+\setcounter{figure}{0}
+\setcounter{table}{0}
+}
diff --git a/doc/references b/doc/references
new file mode 100644
index 0000000..53a38c0
--- /dev/null
+++ b/doc/references
@@ -0,0 +1,12 @@
+Thermodynamic consistency in EOS interpolation
+@ARTICLE{2000ApJS..126..501T,
+    author = {{Timmes}, F.~X. and {Swesty}, F.~D.},
+    title = "{The Accuracy, Consistency, and Speed of an Electron-Positron Equation of State Based on Table Interpolation of the Helmholtz Free Energy}",
+    journal = {Astrophysical Journal Supplement Series},
+    year = 2000,
+    month = feb,
+    volume = 126,
+    pages = {501--516},
+    url = {http://esoads.eso.org/cgi-bin/nph-bib_query?bibcode=2000ApJS..126..501T&db_key=AST},
+    adsnote = {Provided by the NASA Astrophysics Data System}
+}
diff --git a/doc/schedule.aux b/doc/schedule.aux
new file mode 100644
index 0000000..9a76364
--- /dev/null
+++ b/doc/schedule.aux
@@ -0,0 +1,19 @@
+\relax 
+\@setckpt{schedule}{
+\setcounter{page}{18}
+\setcounter{equation}{2}
+\setcounter{enumi}{0}
+\setcounter{enumii}{0}
+\setcounter{enumiii}{0}
+\setcounter{enumiv}{0}
+\setcounter{footnote}{7}
+\setcounter{mpfootnote}{0}
+\setcounter{part}{0}
+\setcounter{section}{5}
+\setcounter{subsection}{4}
+\setcounter{subsubsection}{0}
+\setcounter{paragraph}{0}
+\setcounter{subparagraph}{0}
+\setcounter{figure}{0}
+\setcounter{table}{0}
+}