This commit was generated by cvs2svn to compensate for changes in r2,

which included commits to RCS files with non-trunk default branches.


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

9 files changed, 2901 insertions, 0 deletions

diff --git a/archive/2d.cube.order2.smooth0.c b/archive/2d.cube.order2.smooth0.c
new file mode 100644
index 0000000..34e17ca
--- /dev/null
+++ b/archive/2d.cube.order2.smooth0.c
@@ -0,0 +1,534 @@
+ /*@@
+   @file      2d.cube.order2.smooth0.c
+   @date      23 Oct 2001
+   @author    Jonathan Thornburg <jthorn@aei.mpg.de>
+   @desc
+	Generalized interpolation for 2d, hypercube-shaped molecules,
+	order=2, smoothing=0.  For details, see the header comments
+	for "InterpLocalArrays.c" in this directory.
+   @enddesc
+
+   @version   $Id$
+ @@*/
+
+#include <math.h>
+#include <limits.h>
+#include <stdlib.h>
+#include <string.h>
+#include <stdio.h>
+
+#include "util_ErrorCodes.h"
+#include "cctk.h"
+#include "InterpLocalArrays.h"
+
+/* the rcs ID and its dummy function to use it */
+static const char *rcsid = "$Header$";
+CCTK_FILEVERSION(
+   CactusPUGH_LocalInterp_GeneralizedPolynomial_2d_cube_order2_smooth0_c
+		)
+
+/******************************************************************************/
+
+/*@@
+   @routine    LocalInterp_ILA_2d_cube_ord2_sm0
+   @date       23 Oct 2001
+   @author     Jonathan Thornburg <jthorn@aei.mpg.de>
+   @desc
+	This function does generalized interpolation of one or more
+	2d arrays to arbitrary points.  For details, see the header
+	comments for InterpLocalArrays() (in "InterpLocalArrays.c"
+	in this same directory).
+
+	This function's arguments are all a subset of those of
+	InterpLocalArrays() ; the only difference is that this function
+	takes all its arguments explicitly, whereas  InputLocalArrays()
+	takes some of them indirectly via a key/value parameter table.
+
+  @returntype   int
+  @returndesc	This function's return result is the same as that of
+		InterpLocalArrays():
+  		 0 ==> successful interpolation
+                -1 ==> in case of any errors
+  @endreturndesc
+
+  @@*/
+int LocalInterp_ILA_2d_cube_ord2_sm0
+	(int param_table_handle,
+	 const CCTK_REAL coord_system_origin[],
+	 const CCTK_REAL grid_spacing[],
+	 int N_interp_points,
+	 const CCTK_INT interp_coord_type_codes[],
+	 const void *const interp_coords[],
+	 int N_input_arrays,
+	 const CCTK_INT input_array_offsets[],
+	 const CCTK_INT input_array_strides[],
+	 const CCTK_INT input_array_min_subscripts[],
+	 const CCTK_INT input_array_max_subscripts[],
+	 const CCTK_INT input_array_type_codes[],
+	 const void *const input_arrays[],
+	 int N_output_arrays,
+	 const CCTK_INT output_array_type_codes[],
+	 void *const output_arrays[],
+	 const CCTK_INT operand_indices[], const CCTK_INT opcodes[])
+{
+/*
+ * Implementation notes:
+ * 
+ * The basic outline of this function is as follows:
+ *
+ *	compute "which derivatives are wanted" flags
+ *		for each interpolation point
+ *		{
+ *		declare all the coefficients
+ *		declare all the data-values variables
+ *		***fetch*** interpolation point coordinates
+ *		compute coefficients for all derivatives which are wanted
+ *			for each output array
+ *			{
+ *			int part;
+ *				for (part = 0 ; part <= 1 ; ++part)
+ *				{
+ *				if (this output array is computed
+ *				    using a different input array
+ *				    than the previous one || part != 0)
+ *				   then ***fetch*** the input array values
+ *					            into local variables
+ *				  {
+ *				fp result;
+ *				switch	(opcode)
+ *					{
+ *				case 0:
+ *					result = evaluate the interpolant
+ *					break;
+ *				case 1:
+ *					result = evaluate the interpolant
+ *					break;
+ *				case ...
+ *					}
+ *				***store*** result in output array
+ *				bool complex_flag = is datatype complex?
+ *				if (! complex_flag)
+ *				   then break;
+ *				  }
+ *				}
+ *			}
+ *		}
+ *
+ * Here "***fetch***" and "***store***" are all actually switches on
+ * the various array datatypes.  For complex datatypes they offset the
+ * 1D array position by  part  to handle real/imaginary components of
+ * the complex values.
+ *
+ * At present we do all floating-point computations in type "fp"
+ * (typically a typedef for CCTK_REAL), so arrays of higher precision
+ * than this will incur extra rounding errors.  In practice these should
+ * be negligible compared to the "truncation" interpolation errors.
+ */
+
+/*
+ * Naming conventions:
+ * input, output = 0-origin indices each selecting an input/output array
+ * point = 0-origin index selecting an interpolation point
+ */
+
+/*
+ * these are compile-time constants here; InterpLocalArrays() decoded
+ * them and called us (as opposed to another function) based in part
+ * on these values
+ */
+#define N_DIMS		2
+#define MOLECULE_SIZE	3
+
+/* layout of axes in N_dims-element arrays */
+#define X_AXIS	0
+#define Y_AXIS	1
+
+/* input array size, strides, and subscripting computation */
+const int stride_i = input_array_strides[X_AXIS];
+const int stride_j = input_array_strides[Y_AXIS];
+#define SUB2(i,j)	(i*stride_i + j*stride_j)
+
+/* macros used by machine-generated interpolation coefficient expressions */
+#define RATIONAL(num,den)	(num/den)
+
+/*
+ * compute flags specifying which derivatives are wanted
+ */
+bool want_I = false;
+bool want_dx = false, want_dy = false;
+  {
+int output;
+	for (output = 0 ; output < N_output_arrays ; ++output)
+	{
+	switch	(opcodes[output])
+		{
+	case 0:		want_I = true;		break;
+	case 1:		want_dx = true;		break;
+	case 2:		want_dy = true;		break;
+	default:
+		CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
+			   "Generalized interpolation opcode %d not supported",
+			   opcodes[output]);			/*NOTREACHED*/
+		return UTIL_ERROR_BAD_INPUT;		/*** ERROR RETURN ***/
+		}
+	}
+  }
+
+/*
+ * interpolate at each point
+ */
+  {
+int point;
+	for (point = 0 ; point < N_interp_points ; ++point)
+	{
+	/* declare all the interpolation coefficients */
+	#include "coeffs/2d.cube.order2.smooth0.I.dcl.c"
+	#include "coeffs/2d.cube.order2.smooth0.dx.dcl.c"
+	#include "coeffs/2d.cube.order2.smooth0.dy.dcl.c"
+
+	/* declare all the data-values variables */
+	#include "coeffs/2d.cube.size3.data-var.dcl.c"
+
+	/*
+	 * ***fetch*** interpolation point coordinates
+	 */
+	fp interp_coords_fp[N_DIMS];
+	int axis;
+		for (axis = 0 ; axis < N_DIMS ; ++axis)
+		{
+		/* pointer to array of interp coords for this axis */
+		const void *const interp_coords_ptr = interp_coords[axis];
+
+		switch	(interp_coord_type_codes[axis])
+			{
+#ifdef CCTK_REAL
+		case CCTK_VARIABLE_REAL:
+			  {
+			const CCTK_REAL *const interp_coords_ptr_real
+				= (const CCTK_REAL *) interp_coords_ptr;
+			interp_coords_fp[axis] = interp_coords_ptr_real[point];
+			break;
+			  }
+#endif
+#ifdef CCTK_REAL4
+		case CCTK_VARIABLE_REAL4:
+			  {
+			const CCTK_REAL4 *const interp_coords_ptr_real4
+				= (const CCTK_REAL4 *) interp_coords_ptr;
+			interp_coords_fp[axis] = interp_coords_ptr_real4[point];
+			break;
+			  }
+#endif
+#ifdef CCTK_REAL8
+		case CCTK_VARIABLE_REAL8:
+			  {
+			const CCTK_REAL8 *const interp_coords_ptr_real8
+				= (const CCTK_REAL8 *) interp_coords_ptr;
+			interp_coords_fp[axis] = interp_coords_ptr_real8[point];
+			break;
+			  }
+#endif
+		default:
+			CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
+				   "interp-coords datatype %d not supported",
+				   interp_coord_type_codes[axis]);
+								/*NOTREACHED*/
+			return UTIL_ERROR_BAD_INPUT;	/*** ERROR RETURN ***/
+			}
+		}
+
+	/*
+	 * locate interpolation molecules with respect to the grid,
+	 * i.e. compute interp_rel_(x,y)
+	 */
+	  {
+	fp interp_rel_x, interp_rel_y;	/* (x,y) coordinates of interpolation */
+					/* point relative to molecule center */
+					/* (in units of the grid spacing) */
+	const int center_i
+		= LocalInterp_molecule_posn(coord_system_origin[X_AXIS],
+					   grid_spacing[X_AXIS],
+					   input_array_min_subscripts[X_AXIS],
+					   input_array_max_subscripts[X_AXIS],
+					   MOLECULE_SIZE,
+					   interp_coords_fp[X_AXIS],
+					   &interp_rel_x,
+					   (int *) NULL, (int *) NULL);
+	const int center_j
+		= LocalInterp_molecule_posn(coord_system_origin[Y_AXIS],
+					   grid_spacing[Y_AXIS],
+					   input_array_min_subscripts[Y_AXIS],
+					   input_array_max_subscripts[Y_AXIS],
+					   MOLECULE_SIZE,
+					   interp_coords_fp[Y_AXIS],
+					   &interp_rel_y,
+					   (int *) NULL, (int *) NULL);
+	const int center_sub = SUB2(center_i, center_j);
+
+	/*
+	 * compute the coefficients for whichever derivatives are wanted
+	 * using machine-generated coefficient expressions
+	 * ... these expressions are polynomials in (x,y)
+	 *     ==> we need these names for the relative coordinates
+	 *	   (copying to fresh local variables will likely also
+	 *	    give better register allocation, since [xy]_rel
+	 *	    had their addresses taken and so probably won't be
+	 *	    register-allocated)
+	 */
+	const fp x = interp_rel_x;
+	const fp y = interp_rel_y;
+	if (want_I)
+	   then {
+		#include "coeffs/2d.cube.order2.smooth0.I.coeff.c"
+		}
+	if (want_dx)
+	   then {
+		#include "coeffs/2d.cube.order2.smooth0.dx.coeff.c"
+		}
+	if (want_dy)
+	   then {
+		#include "coeffs/2d.cube.order2.smooth0.dy.coeff.c"
+		}
+
+	/*
+	 * compute each output array at this point
+	 */
+	  {
+	int output;
+	const void *input_array_ptr = NULL;
+		for (output = 0 ; output < N_output_arrays ; ++output)
+		{
+		const int input = operand_indices[output];
+		const int input_offset = input_array_offsets[input];
+
+		/*
+		 * for each real/imag part of complex data values
+		 * ... for real we'll break out of this loop at the bottom
+		 *     after only a single iteration;
+		 *     for complex we'll do both iterations
+		 */
+		int part;
+			for (part = 0 ; part <= 1 ; ++part)
+			{
+			if ( (input_arrays[input] != input_array_ptr)
+			     || (part != 0) )
+			   then {
+				/*
+				 * ***fetch*** the input array values
+				 *             into local variables
+				 */
+				input_array_ptr = input_arrays[input];
+				switch	(input_array_type_codes[input])
+					{
+#ifdef CCTK_REAL
+case CCTK_VARIABLE_REAL:
+	  {
+	const CCTK_REAL *const input_array_ptr_real
+		= (const CCTK_REAL *) input_array_ptr;
+	#undef DATA
+	#define DATA(i,j)	\
+		input_array_ptr_real[input_offset + center_sub + SUB2(i,j)]
+	#include "coeffs/2d.cube.size3.data-var.assign.c"
+	break;
+	  }
+#endif
+#ifdef CCTK_REAL4
+case CCTK_VARIABLE_REAL4:
+	  {
+	const CCTK_REAL4 *const input_array_ptr_real4
+		= (const CCTK_REAL4 *) input_array_ptr;
+	#undef DATA
+	#define DATA(i,j)	\
+		input_array_ptr_real4[input_offset + center_sub + SUB2(i,j)]
+	#include "coeffs/2d.cube.size3.data-var.assign.c"
+	break;
+	  }
+#endif
+#ifdef CCTK_REAL8
+case CCTK_VARIABLE_REAL8:
+	  {
+	const CCTK_REAL8 *const input_array_ptr_real8
+		= (const CCTK_REAL8 *) input_array_ptr;
+	#undef DATA
+	#define DATA(i,j)	\
+		input_array_ptr_real8[input_offset + center_sub + SUB2(i,j)]
+	#include "coeffs/2d.cube.size3.data-var.assign.c"
+	break;
+	  }
+#endif
+#ifdef CCTK_COMPLEX
+case CCTK_VARIABLE_COMPLEX:
+	  {
+	const CCTK_COMPLEX *const input_array_ptr_complex
+		= (const CCTK_COMPLEX *) input_array_ptr;
+	#undef DATA
+	#define DATA(i,j)	\
+   ( (const CCTK_REAL *)	\
+     & input_array_ptr_complex[input_offset + center_sub + SUB2(i,j)] ) [part]
+	#include "coeffs/2d.cube.size3.data-var.assign.c"
+	break;
+	  }
+#endif
+#ifdef CCTK_COMPLEX8
+case CCTK_VARIABLE_COMPLEX8:
+	  {
+	const CCTK_COMPLEX8 *const input_array_ptr_complex8
+		= (const CCTK_COMPLEX8 *) input_array_ptr;
+	#undef DATA
+	#define DATA(i,j)	\
+   ( (const CCTK_REAL4 *)	\
+     & input_array_ptr_complex8[input_offset + center_sub + SUB2(i,j)] ) [part]
+	#include "coeffs/2d.cube.size3.data-var.assign.c"
+	break;
+	  }
+#endif
+#ifdef CCTK_COMPLEX16
+case CCTK_VARIABLE_COMPLEX16:
+	  {
+	const CCTK_COMPLEX16 *const input_array_ptr_complex16
+		= (const CCTK_COMPLEX16 *) input_array_ptr;
+	#undef DATA
+	#define DATA(i,j)	\
+   ( (const CCTK_REAL8 *)	\
+     & input_array_ptr_complex16[input_offset + center_sub + SUB2(i,j)] ) [part]
+	#include "coeffs/2d.cube.size3.data-var.assign.c"
+	break;
+	  }
+#endif
+default:
+	CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
+		   "input datatype %d not supported",
+		   input_array_type_codes[input]);		/*NOTREACHED*/
+	return UTIL_ERROR_BAD_INPUT;			/*** ERROR RETURN ***/
+					}
+				}
+
+			/*
+			 * evaluate the interpolant
+			 * as a linear combination of the data variables
+			 */
+			  {
+			fp result;
+			switch	(opcodes[output])
+				{
+			case 0:
+				result =
+	#include "coeffs/2d.cube.order2.smooth0.I.eval.c"
+				break;
+			case 1:
+				result =
+	#include "coeffs/2d.cube.order2.smooth0.dx.eval.c"
+				break;
+			case 2:
+				result =
+	#include "coeffs/2d.cube.order2.smooth0.dy.eval.c"
+				break;
+			default:
+				CCTK_VWarn(1, __LINE__, __FILE__,
+					   CCTK_THORNSTRING,
+					   "opcode %d not supported",
+					   opcodes[output]);	/*NOTREACHED*/
+				return UTIL_ERROR_BAD_INPUT;
+							/*** ERROR RETURN ***/
+				}
+
+			/*
+			 * ***store*** the result in the output array
+			 */
+			switch	(output_array_type_codes[output])
+				{
+case CCTK_VARIABLE_REAL:
+	  {
+	CCTK_REAL *const output_array_ptr_real
+		= (CCTK_REAL *) output_arrays[output];
+	output_array_ptr_real[point] = result;
+	break;
+	  }
+case CCTK_VARIABLE_REAL4:
+	  {
+	CCTK_REAL4 *const output_array_ptr_real4
+		= (CCTK_REAL4 *) output_arrays[output];
+	output_array_ptr_real4[point] = result;
+	break;
+	  }
+case CCTK_VARIABLE_REAL8:
+	  {
+	CCTK_REAL8 *const output_array_ptr_real8
+		= (CCTK_REAL8 *) output_arrays[output];
+	output_array_ptr_real8[point] = result;
+	break;
+	  }
+case CCTK_VARIABLE_COMPLEX:
+	  {
+	CCTK_COMPLEX *const output_array_ptr_complex
+		= (CCTK_COMPLEX *) output_arrays[output];
+	((CCTK_REAL *)  & output_array_ptr_complex[point]) [part]
+		= result;
+	break;
+	  }
+case CCTK_VARIABLE_COMPLEX8:
+	  {
+	CCTK_COMPLEX8 *const output_array_ptr_complex8
+		= (CCTK_COMPLEX8 *) output_arrays[output];
+	((CCTK_REAL4 *)  & output_array_ptr_complex8[point]) [part]
+		= result;
+	break;
+	  }
+case CCTK_VARIABLE_COMPLEX16:
+	  {
+	CCTK_COMPLEX16 *const output_array_ptr_complex16
+		= (CCTK_COMPLEX16 *) output_arrays[output];
+	((CCTK_REAL8 *)  & output_array_ptr_complex16[point]) [part]
+		= result;
+	break;
+	  }
+default:
+	CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
+		   "output datatype %d not supported",
+		   output_array_type_codes[output]);		/*NOTREACHED*/
+	return UTIL_ERROR_BAD_INPUT;			/*** ERROR RETURN ***/
+				}
+
+			/* decode datatype: is it real or complex? */
+			  {
+			bool complex_flag;
+			switch	(output_array_type_codes[output])
+				{
+			case CCTK_VARIABLE_REAL:
+			case CCTK_VARIABLE_REAL4:
+			case CCTK_VARIABLE_REAL8:
+				complex_flag = false;
+				break;
+			case CCTK_VARIABLE_COMPLEX:
+			case CCTK_VARIABLE_COMPLEX8:
+			case CCTK_VARIABLE_COMPLEX16:
+				complex_flag = true;
+				break;
+			default:
+				CCTK_VWarn(1, __LINE__, __FILE__,
+					   CCTK_THORNSTRING,
+					   "output datatype %d not supported",
+					   output_array_type_codes[output]);
+								/*NOTREACHED*/
+				return UTIL_ERROR_BAD_INPUT;
+							/*** ERROR RETURN ***/
+				}
+			
+			/* skip part=1 (imaginary part) for real datatypes */
+			if (! complex_flag)
+			   then break;
+			  }
+			  }
+			/* end of  for (part = ...)  loop */
+			}
+		/* end of  for (output = ...)  loop */
+		}
+	  }
+	  }
+
+	/* end of  for (point = ...)  loop */
+	}
+
+return 0;						/*** NORMAL RETURN ***/
+  }
+}
diff --git a/archive/C_str.maple b/archive/C_str.maple
new file mode 100644
index 0000000..e0335a6
--- /dev/null
+++ b/archive/C_str.maple
@@ -0,0 +1,42 @@
+#
+# This function is a wrapper around  codegen[C]()  which returns the
+# genenerated code explictly as a Maple string.
+#
+# Arguments:
+# expr = (in) The expression for which code is to be generated.
+# ... = (in) Any further arguments are taken as options to be passed
+#	     to codegen[C]()
+#
+# Results:
+# The function returns a maple string of C code.
+#
+C_str :=
+proc(expr::algebraic)
+local tempname, str, temp;
+
+# name of temp file
+# FIXME: should use process number to ensure uniqueness
+tempname := "/tmp/C_str.tmp.c";
+
+# truncate temp file to zero length
+fopen(tempname, WRITE);
+fclose(tempname);
+
+# generate the code
+codegen[C](args, filename=tempname);
+
+# read the code back in
+str := "";
+	while true
+	do
+	temp := readline(tempname);
+	if (temp = 0)
+	   then break;
+	end if;
+	str := cat(str, temp);
+	end do;
+fclose(tempname);
+
+# strip off the leading "   t0 = "
+return op(2,sscanf(str, "%s = %[^;];"));
+end proc;
diff --git a/archive/README b/archive/README
new file mode 100644
index 0000000..6601932
--- /dev/null
+++ b/archive/README
@@ -0,0 +1,6 @@
+This directory contains stuff that's no longer used, but that I
+thought might be useful for future reference, so I didn't want to
+throw it away.
+
+Compared to the CVS attic, having a separate directory like this
+seems more convenient -- you can 'cat, awk, and grep' more easily.
diff --git a/archive/api1 b/archive/api1
new file mode 100644
index 0000000..5bc5f22
--- /dev/null
+++ b/archive/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/archive/api2 b/archive/api2
new file mode 100644
index 0000000..422be0d
--- /dev/null
+++ b/archive/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/archive/api2.1 b/archive/api2.1
new file mode 100644
index 0000000..e316aba
--- /dev/null
+++ b/archive/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/archive/api2.2 b/archive/api2.2
new file mode 100644
index 0000000..35eca3c
--- /dev/null
+++ b/archive/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/archive/api3 b/archive/api3
new file mode 100644
index 0000000..9aeedd5
--- /dev/null
+++ b/archive/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/archive/caching.tex b/archive/caching.tex
new file mode 100644
index 0000000..9149de6
--- /dev/null
+++ b/archive/caching.tex
@@ -0,0 +1,61 @@
+% $Header$
+
+% This file contains documentation for features which aren't implemented
+% yet, and which are far enough in the future that I (Jonathan) decided
+% to remove them from documentation.tex for the time being.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\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}