add new files which will be (are) needed by Hermite interpolator

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

26 files changed, 1671 insertions, 46 deletions

diff --git a/src/GeneralizedPolynomial-Uniform/common/1d.log b/src/GeneralizedPolynomial-Uniform/common/1d.log
index c3e9168..bf3e9bf 100644
--- a/src/GeneralizedPolynomial-Uniform/common/1d.log
+++ b/src/GeneralizedPolynomial-Uniform/common/1d.log
@@ -444,12 +444,13 @@ ftruncate :=
 
     proc(file_name::string) fopen(file_name, 'WRITE'); fclose(%); NULL end proc
 
-# interpolate.maple -- compute generalized interpolation formulas/coefficients
-# $Id: interpolate.maple,v 1.4 2002/05/14 15:52:50 jthorn Exp $
+# interpolate.maple -- compute interpolation formulas/coefficients
+# $Id: interpolate.maple,v 1.7 2002/08/19 14:30:43 jthorn Exp $
 > 
 #
 # <<<representation of numbers, data values, etc>>>
-# polynomial_interpolant - compute polynomial interpolant
+# Lagrange_polynomial_interpolant - compute Lagrange polynomial interpolant
+# Hermite_polynomial_interpolant - compute Hermite polynomial interpolant
 # coeff_as_lc_of_data - coefficients of ... (linear combination of data)
 #
 # print_coeff__lc_of_data - print C code to compute coefficients
@@ -487,7 +488,8 @@ ftruncate :=
 ################################################################################
 > 
 #
-# This function computes a polynomial interpolant in any number of dimensions.
+# This function computes a Lagrange polynomial interpolant in any
+# number of dimensions.
 #
 # Arguments:
 # fn = The interpolation function.  This should be a procedure in the
@@ -507,7 +509,7 @@ ftruncate :=
 # This function returns the interpolating polynomial, in the form of
 # an algebraic expression in the coordinates and the data values.
 #
-> polynomial_interpolant :=
+> Lagrange_polynomial_interpolant :=
 > proc(
 >       fn::procedure, coeff_list::list(name),
 >       coord_list::list(name), posn_list::list(list(numeric))
@@ -524,7 +526,7 @@ ftruncate :=
 # interpolant as a polynomial in the coordinates
 > return subs(coeff_eqns, eval(fn))(op(coord_list));
 > end proc;
-polynomial_interpolant := proc(fn::procedure, coeff_list::list(name),
+Lagrange_polynomial_interpolant := proc(fn::procedure, coeff_list::list(name),
 coord_list::list(name), posn_list::list(list(numeric)))
 local posn, data_eqns, coeff_eqns;
     data_eqns := {seq(fn(op(posn)) = 'DATA'(op(posn)), posn = posn_list)};
@@ -539,6 +541,171 @@ end proc
 ################################################################################
 > 
 #
+# This function computes a Hermite polynomial interpolant in any
+# number of dimensions.  This is a polynomial which
+#  has values which match the given data DATA() at a specified set of
+#   points, and
+# * has derivatives which match the specified finite-difference derivatives
+#   of the given data DATA() at a specified set of points
+#
+# For the derivative matching, we actually match all possible products
+# of 1st derivatives, i.e. in 2-D we match dx, dy, and dxy, in 3-D we
+# match dx, dy, dz, dxy, dxz, dyz, and dxyz, etc etc.
+#
+# Arguments:
+# fn = The interpolation function.  This should be a procedure in the
+#      coordinates, having the coefficients as global variables.  For
+#      example,
+#		proc(x,y)
+#		+ c03*y^3 + c13*x*y^3 + c23*x^2*y^3 + c33*x^3*y^3
+#		+ c02*y^2 + c12*x*y^2 + c22*x^2*y^2 + c32*x^3*y^2
+#		+ c01*y   + c11*x*y   + c21*x^2*y   + c31*x^3*y
+#		+ c00     + c10*x     + c20*x^2     + c30*x^3
+#		end proc;
+# coeff_list = A set of the interpolation coefficients (coefficients in
+#	       the interpolation function), for example
+#			[
+#			c03, c13, c23, c33,
+#			c02, c12, c22, c32,
+#			c01, c11, c21, c31,
+#			c00, c10, c20, c30
+#			]
+# coord_list = A list of the coordinates (independent variables in the
+#	       interpolation function), for example [x,y].
+# deriv_coord_list = A list of lists of coordinates specifying which
+#		     derivatives are computed by the deriv_proc_list[]
+#		     procedures, for example
+#			[[x], [y], [x,y]]
+# deriv_proc_list = A list of procedures for computing finite-difference
+#		    derivatives.  Each procedure should take N_dims integer
+#		    arguments specifying an evaluation point, and return
+#		    a suitable linear combination of the DATA() for the
+#		    derivative at that point.  For example
+#		    example,
+#			[
+#			  proc(i::integer, j::integer)
+#			  - 1/2*DATA(i-1,j) + 1/2*DATA(i+1,j)
+#			  end proc
+#			,
+#			  proc(i::integer, j::integer)
+#			  - 1/2*DATA(i,j-1) + 1/2*DATA(i,j+1)
+#			  end proc
+#			,
+#			  proc(i::integer, j::integer)
+#			  - 1/4*DATA(i-1,j+1) + 1/4*DATA(i+1,j+1)
+#			  + 1/4*DATA(i-1,j-1) - 1/4*DATA(i+1,j-1)
+#			  end proc
+#			]
+# fn_posn_list = A list of positions (each a list of numeric values)
+#		 where the interpolant is to match the given data DATA(),
+#		 for example
+#			[[0,0], [0,1], [1,0], [1,1]]
+# deriv_posn_list = A list of positions (each a list of numeric values)
+#		    where the interpolant is to match each possible product
+#		    of 1st derivatives"1st"
+#		    derivatives (actually all derivatives the given data DATA(),
+#		    for example
+#			[[0,0], [0,1], [1,0], [1,1]]
+#
+# Results:
+# This function returns the interpolating polynomial, in the form of
+# an algebraic expression in the coordinates and the data values.
+#
+> Hermite_polynomial_interpolant :=
+> proc(
+>       fn::procedure, coeff_list::list(name), coord_list::list(name),
+>       deriv_coord_list::list(list(name)), fd_deriv_proc_list::list(procedure),
+>       fn_posn_list::list(list(numeric)), deriv_posn_list::list(list(numeric))
+>     )
+> local fn_eset,
+>       fn_expr, deriv_eset;
+> 
+# set of equations {fn(posn) = DATA(posn)}
+> fn_eset := map(
+> 		  # return equation that fn(this point) = DATA(this point)
+> 		  proc(posn::list(integer))
+> 		  fn(op(posn)) = 'DATA'(op(posn));
+> 		  end proc
+> 		,
+> 		  {op(fn_posn_list)}
+> 	      );
+> 
+# set of sets of equations
+#	{
+#	  { deriv1(posn1) = fd_deriv1(posn1),
+#	    deriv2(posn1) = fd_deriv2(posn1),
+#	    ... },
+#	  { deriv1(posn2) = fd_deriv1(posn2),
+#	    deriv2(posn2) = fd_deriv2(posn2),
+#	    ... },
+#	  ...
+#	}
+> fn_expr := fn(op(coord_list));
+> map(
+>        # return set of equations
+>        #	{deriv(posn) = fd_deriv(posn)}
+>        # for this point
+>        proc(posn::list(integer))
+>        {
+> 	op(
+> 	 zip(
+> 		# return equation that
+> 		#	deriv(posn) = fd_deriv(posn)
+> 		# for this deriv and this point
+> 		proc(deriv_coords::list(name), fd_deriv_proc::procedure)
+> 		local fn_deriv_proc;
+> 		fn_deriv_proc := unapply( diff(fn_expr,deriv_coords),
+> 					  op(coord_list) );
+> 		fn_deriv_proc(op(posn)) = fd_deriv_proc(op(posn));
+> 		end proc
+> 	      ,
+> 		[op(deriv_coord_list)]
+> 	      ,
+> 		[op(fd_deriv_proc_list)]
+> 	    )
+> 	  )
+>        }
+>        end proc
+>      ,
+>        {op(deriv_posn_list)}
+>    );
+> 
+# set of equations {deriv-i(posn-j) = fd_deriv-i(posn-j)}
+> deriv_eset := `union`(op(%));
+> 
+# solve equations for coeffs
+> solve(fn_eset union deriv_eset, {op(coeff_list)});
+> 
+# interpolant as a polynomial in the coordinates
+> return subs(%, eval(fn))(op(coord_list));
+> end proc;
+Hermite_polynomial_interpolant := proc(fn::procedure, coeff_list::list(name),
+coord_list::list(name), deriv_coord_list::list(list(name)),
+fd_deriv_proc_list::list(procedure), fn_posn_list::list(list(numeric)),
+deriv_posn_list::list(list(numeric)))
+local fn_eset, fn_expr, deriv_eset;
+    fn_eset := map(
+        proc(posn::list(integer)) fn(op(posn)) = 'DATA'(op(posn)) end proc,
+        {op(fn_posn_list)});
+    fn_expr := fn(op(coord_list));
+    map(proc(posn::list(integer))
+            {op(zip(proc(deriv_coords::list(name), fd_deriv_proc::procedure
+            )
+            local fn_deriv_proc;
+                fn_deriv_proc :=
+                    unapply(diff(fn_expr, deriv_coords), op(coord_list));
+                fn_deriv_proc(op(posn)) = fd_deriv_proc(op(posn))
+            end proc, [op(deriv_coord_list)], [op(fd_deriv_proc_list)]))}
+        end proc, {op(deriv_posn_list)});
+    deriv_eset := `union`(op(%));
+    solve(fn_eset union deriv_eset, {op(coeff_list)});
+    return subs(%, eval(fn))(op(coord_list))
+end proc
+
+> 
+################################################################################
+> 
+#
 # This function takes as input an interpolating polynomial, expresses
 # it as a linear combination of the data values, and returns the coefficeints
 # of that form.
@@ -546,8 +713,7 @@ end proc
 # Arguments:
 # interpolant = The interpolating polynomial (an algebraic expression
 #		in the coordinates and the data values).
-# posn_list = The same list of positions as was used to compute the
-#	      interpolating polynomial.
+# posn_list = The same list of data positions used in the interpolant.
 #
 # Results:
 # This function returns the coefficients, as a list of equations of the
@@ -915,7 +1081,7 @@ data_var_name := proc(posn::list(numeric), name_prefix::string)
 end proc
 
 # Maple code to compute lists of point positions in hypercube-shaped molecules
-# $Id: $
+# $Id: cube_posns.maple,v 1.1.1.1 2002/08/18 15:05:38 jthorn Exp $
 > 
 ################################################################################
 > 
@@ -974,6 +1140,18 @@ posn_list_2d_size5 := [[-2, -2], [-1, -2], [0, -2], [1, -2], [2, -2], [-2, -1],
 
     [1, 2], [2, 2]]
 
+> posn_list_2d_size6 := map(ListTools[Reverse],
+> 			  hypercube_points([-2,-2], [+3,+3]));
+posn_list_2d_size6 := [[-2, -2], [-1, -2], [0, -2], [1, -2], [2, -2], [3, -2],
+
+    [-2, -1], [-1, -1], [0, -1], [1, -1], [2, -1], [3, -1], [-2, 0], [-1, 0],
+
+    [0, 0], [1, 0], [2, 0], [3, 0], [-2, 1], [-1, 1], [0, 1], [1, 1], [2, 1],
+
+    [3, 1], [-2, 2], [-1, 2], [0, 2], [1, 2], [2, 2], [3, 2], [-2, 3], [-1, 3],
+
+    [0, 3], [1, 3], [2, 3], [3, 3]]
+
 > 
 ################################################################################
 > 
@@ -1068,8 +1246,84 @@ posn_list_3d_size5 := [[-2, -2, -2], [-1, -2, -2], [0, -2, -2], [1, -2, -2],
 
     [0, 2, 2], [1, 2, 2], [2, 2, 2]]
 
+> posn_list_3d_size6 := map(ListTools[Reverse],
+> 			  hypercube_points([-2,-2,-2], [+3,+3,+3]));
+posn_list_3d_size6 := [[-2, -2, -2], [-1, -2, -2], [0, -2, -2], [1, -2, -2],
+
+    [2, -2, -2], [3, -2, -2], [-2, -1, -2], [-1, -1, -2], [0, -1, -2],
+
+    [1, -1, -2], [2, -1, -2], [3, -1, -2], [-2, 0, -2], [-1, 0, -2], [0, 0, -2],
+
+    [1, 0, -2], [2, 0, -2], [3, 0, -2], [-2, 1, -2], [-1, 1, -2], [0, 1, -2],
+
+    [1, 1, -2], [2, 1, -2], [3, 1, -2], [-2, 2, -2], [-1, 2, -2], [0, 2, -2],
+
+    [1, 2, -2], [2, 2, -2], [3, 2, -2], [-2, 3, -2], [-1, 3, -2], [0, 3, -2],
+
+    [1, 3, -2], [2, 3, -2], [3, 3, -2], [-2, -2, -1], [-1, -2, -1], [0, -2, -1],
+
+    [1, -2, -1], [2, -2, -1], [3, -2, -1], [-2, -1, -1], [-1, -1, -1],
+
+    [0, -1, -1], [1, -1, -1], [2, -1, -1], [3, -1, -1], [-2, 0, -1],
+
+    [-1, 0, -1], [0, 0, -1], [1, 0, -1], [2, 0, -1], [3, 0, -1], [-2, 1, -1],
+
+    [-1, 1, -1], [0, 1, -1], [1, 1, -1], [2, 1, -1], [3, 1, -1], [-2, 2, -1],
+
+    [-1, 2, -1], [0, 2, -1], [1, 2, -1], [2, 2, -1], [3, 2, -1], [-2, 3, -1],
+
+    [-1, 3, -1], [0, 3, -1], [1, 3, -1], [2, 3, -1], [3, 3, -1], [-2, -2, 0],
+
+    [-1, -2, 0], [0, -2, 0], [1, -2, 0], [2, -2, 0], [3, -2, 0], [-2, -1, 0],
+
+    [-1, -1, 0], [0, -1, 0], [1, -1, 0], [2, -1, 0], [3, -1, 0], [-2, 0, 0],
+
+    [-1, 0, 0], [0, 0, 0], [1, 0, 0], [2, 0, 0], [3, 0, 0], [-2, 1, 0],
+
+    [-1, 1, 0], [0, 1, 0], [1, 1, 0], [2, 1, 0], [3, 1, 0], [-2, 2, 0],
+
+    [-1, 2, 0], [0, 2, 0], [1, 2, 0], [2, 2, 0], [3, 2, 0], [-2, 3, 0],
+
+    [-1, 3, 0], [0, 3, 0], [1, 3, 0], [2, 3, 0], [3, 3, 0], [-2, -2, 1],
+
+    [-1, -2, 1], [0, -2, 1], [1, -2, 1], [2, -2, 1], [3, -2, 1], [-2, -1, 1],
+
+    [-1, -1, 1], [0, -1, 1], [1, -1, 1], [2, -1, 1], [3, -1, 1], [-2, 0, 1],
+
+    [-1, 0, 1], [0, 0, 1], [1, 0, 1], [2, 0, 1], [3, 0, 1], [-2, 1, 1],
+
+    [-1, 1, 1], [0, 1, 1], [1, 1, 1], [2, 1, 1], [3, 1, 1], [-2, 2, 1],
+
+    [-1, 2, 1], [0, 2, 1], [1, 2, 1], [2, 2, 1], [3, 2, 1], [-2, 3, 1],
+
+    [-1, 3, 1], [0, 3, 1], [1, 3, 1], [2, 3, 1], [3, 3, 1], [-2, -2, 2],
+
+    [-1, -2, 2], [0, -2, 2], [1, -2, 2], [2, -2, 2], [3, -2, 2], [-2, -1, 2],
+
+    [-1, -1, 2], [0, -1, 2], [1, -1, 2], [2, -1, 2], [3, -1, 2], [-2, 0, 2],
+
+    [-1, 0, 2], [0, 0, 2], [1, 0, 2], [2, 0, 2], [3, 0, 2], [-2, 1, 2],
+
+    [-1, 1, 2], [0, 1, 2], [1, 1, 2], [2, 1, 2], [3, 1, 2], [-2, 2, 2],
+
+    [-1, 2, 2], [0, 2, 2], [1, 2, 2], [2, 2, 2], [3, 2, 2], [-2, 3, 2],
+
+    [-1, 3, 2], [0, 3, 2], [1, 3, 2], [2, 3, 2], [3, 3, 2], [-2, -2, 3],
+
+    [-1, -2, 3], [0, -2, 3], [1, -2, 3], [2, -2, 3], [3, -2, 3], [-2, -1, 3],
+
+    [-1, -1, 3], [0, -1, 3], [1, -1, 3], [2, -1, 3], [3, -1, 3], [-2, 0, 3],
+
+    [-1, 0, 3], [0, 0, 3], [1, 0, 3], [2, 0, 3], [3, 0, 3], [-2, 1, 3],
+
+    [-1, 1, 3], [0, 1, 3], [1, 1, 3], [2, 1, 3], [3, 1, 3], [-2, 2, 3],
+
+    [-1, 2, 3], [0, 2, 3], [1, 2, 3], [2, 2, 3], [3, 2, 3], [-2, 3, 3],
+
+    [-1, 3, 3], [0, 3, 3], [1, 3, 3], [2, 3, 3], [3, 3, 3]]
+
 # Maple code to compute common coefficients for all 1d interpolation schemes
-# $Id: 1d.maple,v 1.4 2002/05/14 15:54:01 jthorn Exp $
+# $Id: 1d.maple,v 1.1.1.1 2002/08/18 15:05:38 jthorn Exp $
 > 
 ################################################################################
 > 
@@ -1306,4 +1560,4 @@ data_var_list_1d_size7 := [
 > 
 ################################################################################
 > quit
-bytes used=716820, alloc=655240, time=0.16
+bytes used=797948, alloc=720764, time=0.12
diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-I.dcl.c b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-I.dcl.c
new file mode 100644
index 0000000..0f93f17
--- /dev/null
+++ b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-I.dcl.c
@@ -0,0 +1,36 @@
+fp coeff_I_m2_m2,
+   coeff_I_m1_m2,
+   coeff_I_0_m2,
+   coeff_I_p1_m2,
+   coeff_I_p2_m2,
+   coeff_I_p3_m2,
+   coeff_I_m2_m1,
+   coeff_I_m1_m1,
+   coeff_I_0_m1,
+   coeff_I_p1_m1,
+   coeff_I_p2_m1,
+   coeff_I_p3_m1,
+   coeff_I_m2_0,
+   coeff_I_m1_0,
+   coeff_I_0_0,
+   coeff_I_p1_0,
+   coeff_I_p2_0,
+   coeff_I_p3_0,
+   coeff_I_m2_p1,
+   coeff_I_m1_p1,
+   coeff_I_0_p1,
+   coeff_I_p1_p1,
+   coeff_I_p2_p1,
+   coeff_I_p3_p1,
+   coeff_I_m2_p2,
+   coeff_I_m1_p2,
+   coeff_I_0_p2,
+   coeff_I_p1_p2,
+   coeff_I_p2_p2,
+   coeff_I_p3_p2,
+   coeff_I_m2_p3,
+   coeff_I_m1_p3,
+   coeff_I_0_p3,
+   coeff_I_p1_p3,
+   coeff_I_p2_p3,
+   coeff_I_p3_p3;
diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-I.store.c b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-I.store.c
new file mode 100644
index 0000000..dfdbfda
--- /dev/null
+++ b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-I.store.c
@@ -0,0 +1,36 @@
+COEFF(-2,-2) = coeff_I_m2_m2;
+COEFF(-1,-2) = coeff_I_m1_m2;
+COEFF(0,-2) = coeff_I_0_m2;
+COEFF(1,-2) = coeff_I_p1_m2;
+COEFF(2,-2) = coeff_I_p2_m2;
+COEFF(3,-2) = coeff_I_p3_m2;
+COEFF(-2,-1) = coeff_I_m2_m1;
+COEFF(-1,-1) = coeff_I_m1_m1;
+COEFF(0,-1) = coeff_I_0_m1;
+COEFF(1,-1) = coeff_I_p1_m1;
+COEFF(2,-1) = coeff_I_p2_m1;
+COEFF(3,-1) = coeff_I_p3_m1;
+COEFF(-2,0) = coeff_I_m2_0;
+COEFF(-1,0) = coeff_I_m1_0;
+COEFF(0,0) = coeff_I_0_0;
+COEFF(1,0) = coeff_I_p1_0;
+COEFF(2,0) = coeff_I_p2_0;
+COEFF(3,0) = coeff_I_p3_0;
+COEFF(-2,1) = coeff_I_m2_p1;
+COEFF(-1,1) = coeff_I_m1_p1;
+COEFF(0,1) = coeff_I_0_p1;
+COEFF(1,1) = coeff_I_p1_p1;
+COEFF(2,1) = coeff_I_p2_p1;
+COEFF(3,1) = coeff_I_p3_p1;
+COEFF(-2,2) = coeff_I_m2_p2;
+COEFF(-1,2) = coeff_I_m1_p2;
+COEFF(0,2) = coeff_I_0_p2;
+COEFF(1,2) = coeff_I_p1_p2;
+COEFF(2,2) = coeff_I_p2_p2;
+COEFF(3,2) = coeff_I_p3_p2;
+COEFF(-2,3) = coeff_I_m2_p3;
+COEFF(-1,3) = coeff_I_m1_p3;
+COEFF(0,3) = coeff_I_0_p3;
+COEFF(1,3) = coeff_I_p1_p3;
+COEFF(2,3) = coeff_I_p2_p3;
+COEFF(3,3) = coeff_I_p3_p3;
diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dx.dcl.c b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dx.dcl.c
new file mode 100644
index 0000000..bf2d3ef
--- /dev/null
+++ b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dx.dcl.c
@@ -0,0 +1,36 @@
+fp coeff_dx_m2_m2,
+   coeff_dx_m1_m2,
+   coeff_dx_0_m2,
+   coeff_dx_p1_m2,
+   coeff_dx_p2_m2,
+   coeff_dx_p3_m2,
+   coeff_dx_m2_m1,
+   coeff_dx_m1_m1,
+   coeff_dx_0_m1,
+   coeff_dx_p1_m1,
+   coeff_dx_p2_m1,
+   coeff_dx_p3_m1,
+   coeff_dx_m2_0,
+   coeff_dx_m1_0,
+   coeff_dx_0_0,
+   coeff_dx_p1_0,
+   coeff_dx_p2_0,
+   coeff_dx_p3_0,
+   coeff_dx_m2_p1,
+   coeff_dx_m1_p1,
+   coeff_dx_0_p1,
+   coeff_dx_p1_p1,
+   coeff_dx_p2_p1,
+   coeff_dx_p3_p1,
+   coeff_dx_m2_p2,
+   coeff_dx_m1_p2,
+   coeff_dx_0_p2,
+   coeff_dx_p1_p2,
+   coeff_dx_p2_p2,
+   coeff_dx_p3_p2,
+   coeff_dx_m2_p3,
+   coeff_dx_m1_p3,
+   coeff_dx_0_p3,
+   coeff_dx_p1_p3,
+   coeff_dx_p2_p3,
+   coeff_dx_p3_p3;
diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dx.store.c b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dx.store.c
new file mode 100644
index 0000000..233b5f3
--- /dev/null
+++ b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dx.store.c
@@ -0,0 +1,36 @@
+COEFF(-2,-2) = factor * coeff_dx_m2_m2;
+COEFF(-1,-2) = factor * coeff_dx_m1_m2;
+COEFF(0,-2) = factor * coeff_dx_0_m2;
+COEFF(1,-2) = factor * coeff_dx_p1_m2;
+COEFF(2,-2) = factor * coeff_dx_p2_m2;
+COEFF(3,-2) = factor * coeff_dx_p3_m2;
+COEFF(-2,-1) = factor * coeff_dx_m2_m1;
+COEFF(-1,-1) = factor * coeff_dx_m1_m1;
+COEFF(0,-1) = factor * coeff_dx_0_m1;
+COEFF(1,-1) = factor * coeff_dx_p1_m1;
+COEFF(2,-1) = factor * coeff_dx_p2_m1;
+COEFF(3,-1) = factor * coeff_dx_p3_m1;
+COEFF(-2,0) = factor * coeff_dx_m2_0;
+COEFF(-1,0) = factor * coeff_dx_m1_0;
+COEFF(0,0) = factor * coeff_dx_0_0;
+COEFF(1,0) = factor * coeff_dx_p1_0;
+COEFF(2,0) = factor * coeff_dx_p2_0;
+COEFF(3,0) = factor * coeff_dx_p3_0;
+COEFF(-2,1) = factor * coeff_dx_m2_p1;
+COEFF(-1,1) = factor * coeff_dx_m1_p1;
+COEFF(0,1) = factor * coeff_dx_0_p1;
+COEFF(1,1) = factor * coeff_dx_p1_p1;
+COEFF(2,1) = factor * coeff_dx_p2_p1;
+COEFF(3,1) = factor * coeff_dx_p3_p1;
+COEFF(-2,2) = factor * coeff_dx_m2_p2;
+COEFF(-1,2) = factor * coeff_dx_m1_p2;
+COEFF(0,2) = factor * coeff_dx_0_p2;
+COEFF(1,2) = factor * coeff_dx_p1_p2;
+COEFF(2,2) = factor * coeff_dx_p2_p2;
+COEFF(3,2) = factor * coeff_dx_p3_p2;
+COEFF(-2,3) = factor * coeff_dx_m2_p3;
+COEFF(-1,3) = factor * coeff_dx_m1_p3;
+COEFF(0,3) = factor * coeff_dx_0_p3;
+COEFF(1,3) = factor * coeff_dx_p1_p3;
+COEFF(2,3) = factor * coeff_dx_p2_p3;
+COEFF(3,3) = factor * coeff_dx_p3_p3;
diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dxx.dcl.c b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dxx.dcl.c
new file mode 100644
index 0000000..be3aa1c
--- /dev/null
+++ b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dxx.dcl.c
@@ -0,0 +1,36 @@
+fp coeff_dxx_m2_m2,
+   coeff_dxx_m1_m2,
+   coeff_dxx_0_m2,
+   coeff_dxx_p1_m2,
+   coeff_dxx_p2_m2,
+   coeff_dxx_p3_m2,
+   coeff_dxx_m2_m1,
+   coeff_dxx_m1_m1,
+   coeff_dxx_0_m1,
+   coeff_dxx_p1_m1,
+   coeff_dxx_p2_m1,
+   coeff_dxx_p3_m1,
+   coeff_dxx_m2_0,
+   coeff_dxx_m1_0,
+   coeff_dxx_0_0,
+   coeff_dxx_p1_0,
+   coeff_dxx_p2_0,
+   coeff_dxx_p3_0,
+   coeff_dxx_m2_p1,
+   coeff_dxx_m1_p1,
+   coeff_dxx_0_p1,
+   coeff_dxx_p1_p1,
+   coeff_dxx_p2_p1,
+   coeff_dxx_p3_p1,
+   coeff_dxx_m2_p2,
+   coeff_dxx_m1_p2,
+   coeff_dxx_0_p2,
+   coeff_dxx_p1_p2,
+   coeff_dxx_p2_p2,
+   coeff_dxx_p3_p2,
+   coeff_dxx_m2_p3,
+   coeff_dxx_m1_p3,
+   coeff_dxx_0_p3,
+   coeff_dxx_p1_p3,
+   coeff_dxx_p2_p3,
+   coeff_dxx_p3_p3;
diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dxx.store.c b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dxx.store.c
new file mode 100644
index 0000000..08c3da8
--- /dev/null
+++ b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dxx.store.c
@@ -0,0 +1,36 @@
+COEFF(-2,-2) = factor * coeff_dxx_m2_m2;
+COEFF(-1,-2) = factor * coeff_dxx_m1_m2;
+COEFF(0,-2) = factor * coeff_dxx_0_m2;
+COEFF(1,-2) = factor * coeff_dxx_p1_m2;
+COEFF(2,-2) = factor * coeff_dxx_p2_m2;
+COEFF(3,-2) = factor * coeff_dxx_p3_m2;
+COEFF(-2,-1) = factor * coeff_dxx_m2_m1;
+COEFF(-1,-1) = factor * coeff_dxx_m1_m1;
+COEFF(0,-1) = factor * coeff_dxx_0_m1;
+COEFF(1,-1) = factor * coeff_dxx_p1_m1;
+COEFF(2,-1) = factor * coeff_dxx_p2_m1;
+COEFF(3,-1) = factor * coeff_dxx_p3_m1;
+COEFF(-2,0) = factor * coeff_dxx_m2_0;
+COEFF(-1,0) = factor * coeff_dxx_m1_0;
+COEFF(0,0) = factor * coeff_dxx_0_0;
+COEFF(1,0) = factor * coeff_dxx_p1_0;
+COEFF(2,0) = factor * coeff_dxx_p2_0;
+COEFF(3,0) = factor * coeff_dxx_p3_0;
+COEFF(-2,1) = factor * coeff_dxx_m2_p1;
+COEFF(-1,1) = factor * coeff_dxx_m1_p1;
+COEFF(0,1) = factor * coeff_dxx_0_p1;
+COEFF(1,1) = factor * coeff_dxx_p1_p1;
+COEFF(2,1) = factor * coeff_dxx_p2_p1;
+COEFF(3,1) = factor * coeff_dxx_p3_p1;
+COEFF(-2,2) = factor * coeff_dxx_m2_p2;
+COEFF(-1,2) = factor * coeff_dxx_m1_p2;
+COEFF(0,2) = factor * coeff_dxx_0_p2;
+COEFF(1,2) = factor * coeff_dxx_p1_p2;
+COEFF(2,2) = factor * coeff_dxx_p2_p2;
+COEFF(3,2) = factor * coeff_dxx_p3_p2;
+COEFF(-2,3) = factor * coeff_dxx_m2_p3;
+COEFF(-1,3) = factor * coeff_dxx_m1_p3;
+COEFF(0,3) = factor * coeff_dxx_0_p3;
+COEFF(1,3) = factor * coeff_dxx_p1_p3;
+COEFF(2,3) = factor * coeff_dxx_p2_p3;
+COEFF(3,3) = factor * coeff_dxx_p3_p3;
diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dxy.dcl.c b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dxy.dcl.c
new file mode 100644
index 0000000..6b7f784
--- /dev/null
+++ b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dxy.dcl.c
@@ -0,0 +1,36 @@
+fp coeff_dxy_m2_m2,
+   coeff_dxy_m1_m2,
+   coeff_dxy_0_m2,
+   coeff_dxy_p1_m2,
+   coeff_dxy_p2_m2,
+   coeff_dxy_p3_m2,
+   coeff_dxy_m2_m1,
+   coeff_dxy_m1_m1,
+   coeff_dxy_0_m1,
+   coeff_dxy_p1_m1,
+   coeff_dxy_p2_m1,
+   coeff_dxy_p3_m1,
+   coeff_dxy_m2_0,
+   coeff_dxy_m1_0,
+   coeff_dxy_0_0,
+   coeff_dxy_p1_0,
+   coeff_dxy_p2_0,
+   coeff_dxy_p3_0,
+   coeff_dxy_m2_p1,
+   coeff_dxy_m1_p1,
+   coeff_dxy_0_p1,
+   coeff_dxy_p1_p1,
+   coeff_dxy_p2_p1,
+   coeff_dxy_p3_p1,
+   coeff_dxy_m2_p2,
+   coeff_dxy_m1_p2,
+   coeff_dxy_0_p2,
+   coeff_dxy_p1_p2,
+   coeff_dxy_p2_p2,
+   coeff_dxy_p3_p2,
+   coeff_dxy_m2_p3,
+   coeff_dxy_m1_p3,
+   coeff_dxy_0_p3,
+   coeff_dxy_p1_p3,
+   coeff_dxy_p2_p3,
+   coeff_dxy_p3_p3;
diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dxy.store.c b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dxy.store.c
new file mode 100644
index 0000000..32313f6
--- /dev/null
+++ b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dxy.store.c
@@ -0,0 +1,36 @@
+COEFF(-2,-2) = factor * coeff_dxy_m2_m2;
+COEFF(-1,-2) = factor * coeff_dxy_m1_m2;
+COEFF(0,-2) = factor * coeff_dxy_0_m2;
+COEFF(1,-2) = factor * coeff_dxy_p1_m2;
+COEFF(2,-2) = factor * coeff_dxy_p2_m2;
+COEFF(3,-2) = factor * coeff_dxy_p3_m2;
+COEFF(-2,-1) = factor * coeff_dxy_m2_m1;
+COEFF(-1,-1) = factor * coeff_dxy_m1_m1;
+COEFF(0,-1) = factor * coeff_dxy_0_m1;
+COEFF(1,-1) = factor * coeff_dxy_p1_m1;
+COEFF(2,-1) = factor * coeff_dxy_p2_m1;
+COEFF(3,-1) = factor * coeff_dxy_p3_m1;
+COEFF(-2,0) = factor * coeff_dxy_m2_0;
+COEFF(-1,0) = factor * coeff_dxy_m1_0;
+COEFF(0,0) = factor * coeff_dxy_0_0;
+COEFF(1,0) = factor * coeff_dxy_p1_0;
+COEFF(2,0) = factor * coeff_dxy_p2_0;
+COEFF(3,0) = factor * coeff_dxy_p3_0;
+COEFF(-2,1) = factor * coeff_dxy_m2_p1;
+COEFF(-1,1) = factor * coeff_dxy_m1_p1;
+COEFF(0,1) = factor * coeff_dxy_0_p1;
+COEFF(1,1) = factor * coeff_dxy_p1_p1;
+COEFF(2,1) = factor * coeff_dxy_p2_p1;
+COEFF(3,1) = factor * coeff_dxy_p3_p1;
+COEFF(-2,2) = factor * coeff_dxy_m2_p2;
+COEFF(-1,2) = factor * coeff_dxy_m1_p2;
+COEFF(0,2) = factor * coeff_dxy_0_p2;
+COEFF(1,2) = factor * coeff_dxy_p1_p2;
+COEFF(2,2) = factor * coeff_dxy_p2_p2;
+COEFF(3,2) = factor * coeff_dxy_p3_p2;
+COEFF(-2,3) = factor * coeff_dxy_m2_p3;
+COEFF(-1,3) = factor * coeff_dxy_m1_p3;
+COEFF(0,3) = factor * coeff_dxy_0_p3;
+COEFF(1,3) = factor * coeff_dxy_p1_p3;
+COEFF(2,3) = factor * coeff_dxy_p2_p3;
+COEFF(3,3) = factor * coeff_dxy_p3_p3;
diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dy.dcl.c b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dy.dcl.c
new file mode 100644
index 0000000..489d18d
--- /dev/null
+++ b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dy.dcl.c
@@ -0,0 +1,36 @@
+fp coeff_dy_m2_m2,
+   coeff_dy_m1_m2,
+   coeff_dy_0_m2,
+   coeff_dy_p1_m2,
+   coeff_dy_p2_m2,
+   coeff_dy_p3_m2,
+   coeff_dy_m2_m1,
+   coeff_dy_m1_m1,
+   coeff_dy_0_m1,
+   coeff_dy_p1_m1,
+   coeff_dy_p2_m1,
+   coeff_dy_p3_m1,
+   coeff_dy_m2_0,
+   coeff_dy_m1_0,
+   coeff_dy_0_0,
+   coeff_dy_p1_0,
+   coeff_dy_p2_0,
+   coeff_dy_p3_0,
+   coeff_dy_m2_p1,
+   coeff_dy_m1_p1,
+   coeff_dy_0_p1,
+   coeff_dy_p1_p1,
+   coeff_dy_p2_p1,
+   coeff_dy_p3_p1,
+   coeff_dy_m2_p2,
+   coeff_dy_m1_p2,
+   coeff_dy_0_p2,
+   coeff_dy_p1_p2,
+   coeff_dy_p2_p2,
+   coeff_dy_p3_p2,
+   coeff_dy_m2_p3,
+   coeff_dy_m1_p3,
+   coeff_dy_0_p3,
+   coeff_dy_p1_p3,
+   coeff_dy_p2_p3,
+   coeff_dy_p3_p3;
diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dy.store.c b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dy.store.c
new file mode 100644
index 0000000..4ab90c2
--- /dev/null
+++ b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dy.store.c
@@ -0,0 +1,36 @@
+COEFF(-2,-2) = factor * coeff_dy_m2_m2;
+COEFF(-1,-2) = factor * coeff_dy_m1_m2;
+COEFF(0,-2) = factor * coeff_dy_0_m2;
+COEFF(1,-2) = factor * coeff_dy_p1_m2;
+COEFF(2,-2) = factor * coeff_dy_p2_m2;
+COEFF(3,-2) = factor * coeff_dy_p3_m2;
+COEFF(-2,-1) = factor * coeff_dy_m2_m1;
+COEFF(-1,-1) = factor * coeff_dy_m1_m1;
+COEFF(0,-1) = factor * coeff_dy_0_m1;
+COEFF(1,-1) = factor * coeff_dy_p1_m1;
+COEFF(2,-1) = factor * coeff_dy_p2_m1;
+COEFF(3,-1) = factor * coeff_dy_p3_m1;
+COEFF(-2,0) = factor * coeff_dy_m2_0;
+COEFF(-1,0) = factor * coeff_dy_m1_0;
+COEFF(0,0) = factor * coeff_dy_0_0;
+COEFF(1,0) = factor * coeff_dy_p1_0;
+COEFF(2,0) = factor * coeff_dy_p2_0;
+COEFF(3,0) = factor * coeff_dy_p3_0;
+COEFF(-2,1) = factor * coeff_dy_m2_p1;
+COEFF(-1,1) = factor * coeff_dy_m1_p1;
+COEFF(0,1) = factor * coeff_dy_0_p1;
+COEFF(1,1) = factor * coeff_dy_p1_p1;
+COEFF(2,1) = factor * coeff_dy_p2_p1;
+COEFF(3,1) = factor * coeff_dy_p3_p1;
+COEFF(-2,2) = factor * coeff_dy_m2_p2;
+COEFF(-1,2) = factor * coeff_dy_m1_p2;
+COEFF(0,2) = factor * coeff_dy_0_p2;
+COEFF(1,2) = factor * coeff_dy_p1_p2;
+COEFF(2,2) = factor * coeff_dy_p2_p2;
+COEFF(3,2) = factor * coeff_dy_p3_p2;
+COEFF(-2,3) = factor * coeff_dy_m2_p3;
+COEFF(-1,3) = factor * coeff_dy_m1_p3;
+COEFF(0,3) = factor * coeff_dy_0_p3;
+COEFF(1,3) = factor * coeff_dy_p1_p3;
+COEFF(2,3) = factor * coeff_dy_p2_p3;
+COEFF(3,3) = factor * coeff_dy_p3_p3;
diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dyy.dcl.c b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dyy.dcl.c
new file mode 100644
index 0000000..056404a
--- /dev/null
+++ b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dyy.dcl.c
@@ -0,0 +1,36 @@
+fp coeff_dyy_m2_m2,
+   coeff_dyy_m1_m2,
+   coeff_dyy_0_m2,
+   coeff_dyy_p1_m2,
+   coeff_dyy_p2_m2,
+   coeff_dyy_p3_m2,
+   coeff_dyy_m2_m1,
+   coeff_dyy_m1_m1,
+   coeff_dyy_0_m1,
+   coeff_dyy_p1_m1,
+   coeff_dyy_p2_m1,
+   coeff_dyy_p3_m1,
+   coeff_dyy_m2_0,
+   coeff_dyy_m1_0,
+   coeff_dyy_0_0,
+   coeff_dyy_p1_0,
+   coeff_dyy_p2_0,
+   coeff_dyy_p3_0,
+   coeff_dyy_m2_p1,
+   coeff_dyy_m1_p1,
+   coeff_dyy_0_p1,
+   coeff_dyy_p1_p1,
+   coeff_dyy_p2_p1,
+   coeff_dyy_p3_p1,
+   coeff_dyy_m2_p2,
+   coeff_dyy_m1_p2,
+   coeff_dyy_0_p2,
+   coeff_dyy_p1_p2,
+   coeff_dyy_p2_p2,
+   coeff_dyy_p3_p2,
+   coeff_dyy_m2_p3,
+   coeff_dyy_m1_p3,
+   coeff_dyy_0_p3,
+   coeff_dyy_p1_p3,
+   coeff_dyy_p2_p3,
+   coeff_dyy_p3_p3;
diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dyy.store.c b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dyy.store.c
new file mode 100644
index 0000000..1a77da9
--- /dev/null
+++ b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/coeff-dyy.store.c
@@ -0,0 +1,36 @@
+COEFF(-2,-2) = factor * coeff_dyy_m2_m2;
+COEFF(-1,-2) = factor * coeff_dyy_m1_m2;
+COEFF(0,-2) = factor * coeff_dyy_0_m2;
+COEFF(1,-2) = factor * coeff_dyy_p1_m2;
+COEFF(2,-2) = factor * coeff_dyy_p2_m2;
+COEFF(3,-2) = factor * coeff_dyy_p3_m2;
+COEFF(-2,-1) = factor * coeff_dyy_m2_m1;
+COEFF(-1,-1) = factor * coeff_dyy_m1_m1;
+COEFF(0,-1) = factor * coeff_dyy_0_m1;
+COEFF(1,-1) = factor * coeff_dyy_p1_m1;
+COEFF(2,-1) = factor * coeff_dyy_p2_m1;
+COEFF(3,-1) = factor * coeff_dyy_p3_m1;
+COEFF(-2,0) = factor * coeff_dyy_m2_0;
+COEFF(-1,0) = factor * coeff_dyy_m1_0;
+COEFF(0,0) = factor * coeff_dyy_0_0;
+COEFF(1,0) = factor * coeff_dyy_p1_0;
+COEFF(2,0) = factor * coeff_dyy_p2_0;
+COEFF(3,0) = factor * coeff_dyy_p3_0;
+COEFF(-2,1) = factor * coeff_dyy_m2_p1;
+COEFF(-1,1) = factor * coeff_dyy_m1_p1;
+COEFF(0,1) = factor * coeff_dyy_0_p1;
+COEFF(1,1) = factor * coeff_dyy_p1_p1;
+COEFF(2,1) = factor * coeff_dyy_p2_p1;
+COEFF(3,1) = factor * coeff_dyy_p3_p1;
+COEFF(-2,2) = factor * coeff_dyy_m2_p2;
+COEFF(-1,2) = factor * coeff_dyy_m1_p2;
+COEFF(0,2) = factor * coeff_dyy_0_p2;
+COEFF(1,2) = factor * coeff_dyy_p1_p2;
+COEFF(2,2) = factor * coeff_dyy_p2_p2;
+COEFF(3,2) = factor * coeff_dyy_p3_p2;
+COEFF(-2,3) = factor * coeff_dyy_m2_p3;
+COEFF(-1,3) = factor * coeff_dyy_m1_p3;
+COEFF(0,3) = factor * coeff_dyy_0_p3;
+COEFF(1,3) = factor * coeff_dyy_p1_p3;
+COEFF(2,3) = factor * coeff_dyy_p2_p3;
+COEFF(3,3) = factor * coeff_dyy_p3_p3;
diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/data-var.assign.c b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/data-var.assign.c
new file mode 100644
index 0000000..0172f4a
--- /dev/null
+++ b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/data-var.assign.c
@@ -0,0 +1,36 @@
+data_m2_m2 = DATA(-2,-2);
+data_m1_m2 = DATA(-1,-2);
+data_0_m2 = DATA(0,-2);
+data_p1_m2 = DATA(1,-2);
+data_p2_m2 = DATA(2,-2);
+data_p3_m2 = DATA(3,-2);
+data_m2_m1 = DATA(-2,-1);
+data_m1_m1 = DATA(-1,-1);
+data_0_m1 = DATA(0,-1);
+data_p1_m1 = DATA(1,-1);
+data_p2_m1 = DATA(2,-1);
+data_p3_m1 = DATA(3,-1);
+data_m2_0 = DATA(-2,0);
+data_m1_0 = DATA(-1,0);
+data_0_0 = DATA(0,0);
+data_p1_0 = DATA(1,0);
+data_p2_0 = DATA(2,0);
+data_p3_0 = DATA(3,0);
+data_m2_p1 = DATA(-2,1);
+data_m1_p1 = DATA(-1,1);
+data_0_p1 = DATA(0,1);
+data_p1_p1 = DATA(1,1);
+data_p2_p1 = DATA(2,1);
+data_p3_p1 = DATA(3,1);
+data_m2_p2 = DATA(-2,2);
+data_m1_p2 = DATA(-1,2);
+data_0_p2 = DATA(0,2);
+data_p1_p2 = DATA(1,2);
+data_p2_p2 = DATA(2,2);
+data_p3_p2 = DATA(3,2);
+data_m2_p3 = DATA(-2,3);
+data_m1_p3 = DATA(-1,3);
+data_0_p3 = DATA(0,3);
+data_p1_p3 = DATA(1,3);
+data_p2_p3 = DATA(2,3);
+data_p3_p3 = DATA(3,3);
diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/data-var.dcl.c b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/data-var.dcl.c
new file mode 100644
index 0000000..342b513
--- /dev/null
+++ b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/data-var.dcl.c
@@ -0,0 +1,36 @@
+fp data_m2_m2,
+   data_m1_m2,
+   data_0_m2,
+   data_p1_m2,
+   data_p2_m2,
+   data_p3_m2,
+   data_m2_m1,
+   data_m1_m1,
+   data_0_m1,
+   data_p1_m1,
+   data_p2_m1,
+   data_p3_m1,
+   data_m2_0,
+   data_m1_0,
+   data_0_0,
+   data_p1_0,
+   data_p2_0,
+   data_p3_0,
+   data_m2_p1,
+   data_m1_p1,
+   data_0_p1,
+   data_p1_p1,
+   data_p2_p1,
+   data_p3_p1,
+   data_m2_p2,
+   data_m1_p2,
+   data_0_p2,
+   data_p1_p2,
+   data_p2_p2,
+   data_p3_p2,
+   data_m2_p3,
+   data_m1_p3,
+   data_0_p3,
+   data_p1_p3,
+   data_p2_p3,
+   data_p3_p3;
diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/interp-I.compute.c b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/interp-I.compute.c
new file mode 100644
index 0000000..90eaa22
--- /dev/null
+++ b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/interp-I.compute.c
@@ -0,0 +1,37 @@
+result =
+	coeff_I_m2_m2*data_m2_m2
+	+ coeff_I_m1_m2*data_m1_m2
+	+ coeff_I_0_m2*data_0_m2
+	+ coeff_I_p1_m2*data_p1_m2
+	+ coeff_I_p2_m2*data_p2_m2
+	+ coeff_I_p3_m2*data_p3_m2
+	+ coeff_I_m2_m1*data_m2_m1
+	+ coeff_I_m1_m1*data_m1_m1
+	+ coeff_I_0_m1*data_0_m1
+	+ coeff_I_p1_m1*data_p1_m1
+	+ coeff_I_p2_m1*data_p2_m1
+	+ coeff_I_p3_m1*data_p3_m1
+	+ coeff_I_m2_0*data_m2_0
+	+ coeff_I_m1_0*data_m1_0
+	+ coeff_I_0_0*data_0_0
+	+ coeff_I_p1_0*data_p1_0
+	+ coeff_I_p2_0*data_p2_0
+	+ coeff_I_p3_0*data_p3_0
+	+ coeff_I_m2_p1*data_m2_p1
+	+ coeff_I_m1_p1*data_m1_p1
+	+ coeff_I_0_p1*data_0_p1
+	+ coeff_I_p1_p1*data_p1_p1
+	+ coeff_I_p2_p1*data_p2_p1
+	+ coeff_I_p3_p1*data_p3_p1
+	+ coeff_I_m2_p2*data_m2_p2
+	+ coeff_I_m1_p2*data_m1_p2
+	+ coeff_I_0_p2*data_0_p2
+	+ coeff_I_p1_p2*data_p1_p2
+	+ coeff_I_p2_p2*data_p2_p2
+	+ coeff_I_p3_p2*data_p3_p2
+	+ coeff_I_m2_p3*data_m2_p3
+	+ coeff_I_m1_p3*data_m1_p3
+	+ coeff_I_0_p3*data_0_p3
+	+ coeff_I_p1_p3*data_p1_p3
+	+ coeff_I_p2_p3*data_p2_p3
+	+ coeff_I_p3_p3*data_p3_p3;
diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/interp-dx.compute.c b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/interp-dx.compute.c
new file mode 100644
index 0000000..b67994e
--- /dev/null
+++ b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/interp-dx.compute.c
@@ -0,0 +1,37 @@
+result =
+	coeff_dx_m2_m2*data_m2_m2
+	+ coeff_dx_m1_m2*data_m1_m2
+	+ coeff_dx_0_m2*data_0_m2
+	+ coeff_dx_p1_m2*data_p1_m2
+	+ coeff_dx_p2_m2*data_p2_m2
+	+ coeff_dx_p3_m2*data_p3_m2
+	+ coeff_dx_m2_m1*data_m2_m1
+	+ coeff_dx_m1_m1*data_m1_m1
+	+ coeff_dx_0_m1*data_0_m1
+	+ coeff_dx_p1_m1*data_p1_m1
+	+ coeff_dx_p2_m1*data_p2_m1
+	+ coeff_dx_p3_m1*data_p3_m1
+	+ coeff_dx_m2_0*data_m2_0
+	+ coeff_dx_m1_0*data_m1_0
+	+ coeff_dx_0_0*data_0_0
+	+ coeff_dx_p1_0*data_p1_0
+	+ coeff_dx_p2_0*data_p2_0
+	+ coeff_dx_p3_0*data_p3_0
+	+ coeff_dx_m2_p1*data_m2_p1
+	+ coeff_dx_m1_p1*data_m1_p1
+	+ coeff_dx_0_p1*data_0_p1
+	+ coeff_dx_p1_p1*data_p1_p1
+	+ coeff_dx_p2_p1*data_p2_p1
+	+ coeff_dx_p3_p1*data_p3_p1
+	+ coeff_dx_m2_p2*data_m2_p2
+	+ coeff_dx_m1_p2*data_m1_p2
+	+ coeff_dx_0_p2*data_0_p2
+	+ coeff_dx_p1_p2*data_p1_p2
+	+ coeff_dx_p2_p2*data_p2_p2
+	+ coeff_dx_p3_p2*data_p3_p2
+	+ coeff_dx_m2_p3*data_m2_p3
+	+ coeff_dx_m1_p3*data_m1_p3
+	+ coeff_dx_0_p3*data_0_p3
+	+ coeff_dx_p1_p3*data_p1_p3
+	+ coeff_dx_p2_p3*data_p2_p3
+	+ coeff_dx_p3_p3*data_p3_p3;
diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/interp-dxx.compute.c b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/interp-dxx.compute.c
new file mode 100644
index 0000000..cef82ff
--- /dev/null
+++ b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/interp-dxx.compute.c
@@ -0,0 +1,37 @@
+result =
+	coeff_dxx_m2_m2*data_m2_m2
+	+ coeff_dxx_m1_m2*data_m1_m2
+	+ coeff_dxx_0_m2*data_0_m2
+	+ coeff_dxx_p1_m2*data_p1_m2
+	+ coeff_dxx_p2_m2*data_p2_m2
+	+ coeff_dxx_p3_m2*data_p3_m2
+	+ coeff_dxx_m2_m1*data_m2_m1
+	+ coeff_dxx_m1_m1*data_m1_m1
+	+ coeff_dxx_0_m1*data_0_m1
+	+ coeff_dxx_p1_m1*data_p1_m1
+	+ coeff_dxx_p2_m1*data_p2_m1
+	+ coeff_dxx_p3_m1*data_p3_m1
+	+ coeff_dxx_m2_0*data_m2_0
+	+ coeff_dxx_m1_0*data_m1_0
+	+ coeff_dxx_0_0*data_0_0
+	+ coeff_dxx_p1_0*data_p1_0
+	+ coeff_dxx_p2_0*data_p2_0
+	+ coeff_dxx_p3_0*data_p3_0
+	+ coeff_dxx_m2_p1*data_m2_p1
+	+ coeff_dxx_m1_p1*data_m1_p1
+	+ coeff_dxx_0_p1*data_0_p1
+	+ coeff_dxx_p1_p1*data_p1_p1
+	+ coeff_dxx_p2_p1*data_p2_p1
+	+ coeff_dxx_p3_p1*data_p3_p1
+	+ coeff_dxx_m2_p2*data_m2_p2
+	+ coeff_dxx_m1_p2*data_m1_p2
+	+ coeff_dxx_0_p2*data_0_p2
+	+ coeff_dxx_p1_p2*data_p1_p2
+	+ coeff_dxx_p2_p2*data_p2_p2
+	+ coeff_dxx_p3_p2*data_p3_p2
+	+ coeff_dxx_m2_p3*data_m2_p3
+	+ coeff_dxx_m1_p3*data_m1_p3
+	+ coeff_dxx_0_p3*data_0_p3
+	+ coeff_dxx_p1_p3*data_p1_p3
+	+ coeff_dxx_p2_p3*data_p2_p3
+	+ coeff_dxx_p3_p3*data_p3_p3;
diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/interp-dxy.compute.c b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/interp-dxy.compute.c
new file mode 100644
index 0000000..7c66027
--- /dev/null
+++ b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/interp-dxy.compute.c
@@ -0,0 +1,37 @@
+result =
+	coeff_dxy_m2_m2*data_m2_m2
+	+ coeff_dxy_m1_m2*data_m1_m2
+	+ coeff_dxy_0_m2*data_0_m2
+	+ coeff_dxy_p1_m2*data_p1_m2
+	+ coeff_dxy_p2_m2*data_p2_m2
+	+ coeff_dxy_p3_m2*data_p3_m2
+	+ coeff_dxy_m2_m1*data_m2_m1
+	+ coeff_dxy_m1_m1*data_m1_m1
+	+ coeff_dxy_0_m1*data_0_m1
+	+ coeff_dxy_p1_m1*data_p1_m1
+	+ coeff_dxy_p2_m1*data_p2_m1
+	+ coeff_dxy_p3_m1*data_p3_m1
+	+ coeff_dxy_m2_0*data_m2_0
+	+ coeff_dxy_m1_0*data_m1_0
+	+ coeff_dxy_0_0*data_0_0
+	+ coeff_dxy_p1_0*data_p1_0
+	+ coeff_dxy_p2_0*data_p2_0
+	+ coeff_dxy_p3_0*data_p3_0
+	+ coeff_dxy_m2_p1*data_m2_p1
+	+ coeff_dxy_m1_p1*data_m1_p1
+	+ coeff_dxy_0_p1*data_0_p1
+	+ coeff_dxy_p1_p1*data_p1_p1
+	+ coeff_dxy_p2_p1*data_p2_p1
+	+ coeff_dxy_p3_p1*data_p3_p1
+	+ coeff_dxy_m2_p2*data_m2_p2
+	+ coeff_dxy_m1_p2*data_m1_p2
+	+ coeff_dxy_0_p2*data_0_p2
+	+ coeff_dxy_p1_p2*data_p1_p2
+	+ coeff_dxy_p2_p2*data_p2_p2
+	+ coeff_dxy_p3_p2*data_p3_p2
+	+ coeff_dxy_m2_p3*data_m2_p3
+	+ coeff_dxy_m1_p3*data_m1_p3
+	+ coeff_dxy_0_p3*data_0_p3
+	+ coeff_dxy_p1_p3*data_p1_p3
+	+ coeff_dxy_p2_p3*data_p2_p3
+	+ coeff_dxy_p3_p3*data_p3_p3;
diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/interp-dy.compute.c b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/interp-dy.compute.c
new file mode 100644
index 0000000..ae07fcb
--- /dev/null
+++ b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/interp-dy.compute.c
@@ -0,0 +1,37 @@
+result =
+	coeff_dy_m2_m2*data_m2_m2
+	+ coeff_dy_m1_m2*data_m1_m2
+	+ coeff_dy_0_m2*data_0_m2
+	+ coeff_dy_p1_m2*data_p1_m2
+	+ coeff_dy_p2_m2*data_p2_m2
+	+ coeff_dy_p3_m2*data_p3_m2
+	+ coeff_dy_m2_m1*data_m2_m1
+	+ coeff_dy_m1_m1*data_m1_m1
+	+ coeff_dy_0_m1*data_0_m1
+	+ coeff_dy_p1_m1*data_p1_m1
+	+ coeff_dy_p2_m1*data_p2_m1
+	+ coeff_dy_p3_m1*data_p3_m1
+	+ coeff_dy_m2_0*data_m2_0
+	+ coeff_dy_m1_0*data_m1_0
+	+ coeff_dy_0_0*data_0_0
+	+ coeff_dy_p1_0*data_p1_0
+	+ coeff_dy_p2_0*data_p2_0
+	+ coeff_dy_p3_0*data_p3_0
+	+ coeff_dy_m2_p1*data_m2_p1
+	+ coeff_dy_m1_p1*data_m1_p1
+	+ coeff_dy_0_p1*data_0_p1
+	+ coeff_dy_p1_p1*data_p1_p1
+	+ coeff_dy_p2_p1*data_p2_p1
+	+ coeff_dy_p3_p1*data_p3_p1
+	+ coeff_dy_m2_p2*data_m2_p2
+	+ coeff_dy_m1_p2*data_m1_p2
+	+ coeff_dy_0_p2*data_0_p2
+	+ coeff_dy_p1_p2*data_p1_p2
+	+ coeff_dy_p2_p2*data_p2_p2
+	+ coeff_dy_p3_p2*data_p3_p2
+	+ coeff_dy_m2_p3*data_m2_p3
+	+ coeff_dy_m1_p3*data_m1_p3
+	+ coeff_dy_0_p3*data_0_p3
+	+ coeff_dy_p1_p3*data_p1_p3
+	+ coeff_dy_p2_p3*data_p2_p3
+	+ coeff_dy_p3_p3*data_p3_p3;
diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/interp-dyy.compute.c b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/interp-dyy.compute.c
new file mode 100644
index 0000000..dea1cd3
--- /dev/null
+++ b/src/GeneralizedPolynomial-Uniform/common/2d.cube.size6/interp-dyy.compute.c
@@ -0,0 +1,37 @@
+result =
+	coeff_dyy_m2_m2*data_m2_m2
+	+ coeff_dyy_m1_m2*data_m1_m2
+	+ coeff_dyy_0_m2*data_0_m2
+	+ coeff_dyy_p1_m2*data_p1_m2
+	+ coeff_dyy_p2_m2*data_p2_m2
+	+ coeff_dyy_p3_m2*data_p3_m2
+	+ coeff_dyy_m2_m1*data_m2_m1
+	+ coeff_dyy_m1_m1*data_m1_m1
+	+ coeff_dyy_0_m1*data_0_m1
+	+ coeff_dyy_p1_m1*data_p1_m1
+	+ coeff_dyy_p2_m1*data_p2_m1
+	+ coeff_dyy_p3_m1*data_p3_m1
+	+ coeff_dyy_m2_0*data_m2_0
+	+ coeff_dyy_m1_0*data_m1_0
+	+ coeff_dyy_0_0*data_0_0
+	+ coeff_dyy_p1_0*data_p1_0
+	+ coeff_dyy_p2_0*data_p2_0
+	+ coeff_dyy_p3_0*data_p3_0
+	+ coeff_dyy_m2_p1*data_m2_p1
+	+ coeff_dyy_m1_p1*data_m1_p1
+	+ coeff_dyy_0_p1*data_0_p1
+	+ coeff_dyy_p1_p1*data_p1_p1
+	+ coeff_dyy_p2_p1*data_p2_p1
+	+ coeff_dyy_p3_p1*data_p3_p1
+	+ coeff_dyy_m2_p2*data_m2_p2
+	+ coeff_dyy_m1_p2*data_m1_p2
+	+ coeff_dyy_0_p2*data_0_p2
+	+ coeff_dyy_p1_p2*data_p1_p2
+	+ coeff_dyy_p2_p2*data_p2_p2
+	+ coeff_dyy_p3_p2*data_p3_p2
+	+ coeff_dyy_m2_p3*data_m2_p3
+	+ coeff_dyy_m1_p3*data_m1_p3
+	+ coeff_dyy_0_p3*data_0_p3
+	+ coeff_dyy_p1_p3*data_p1_p3
+	+ coeff_dyy_p2_p3*data_p2_p3
+	+ coeff_dyy_p3_p3*data_p3_p3;
diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.log b/src/GeneralizedPolynomial-Uniform/common/2d.log
index 835c96e..909de57 100644
--- a/src/GeneralizedPolynomial-Uniform/common/2d.log
+++ b/src/GeneralizedPolynomial-Uniform/common/2d.log
@@ -4,7 +4,7 @@
  <____ ____>  Waterloo Maple Inc.
       |       Type ? for help.
 # util.maple -- misc utility routines
-# $Id: util.maple,v 1.3 2002/05/19 13:12:18 jthorn Exp $
+# $Header: /cactusdevcvs/CactusBase/LocalInterp/src/GeneralizedPolynomial-Uniform/util.maple,v 1.4 2002/08/20 16:46:06 jthorn Exp $
 > 
 #
 # fix_rationals - convert numbers to RATIONAL() calls
@@ -444,12 +444,13 @@ ftruncate :=
 
     proc(file_name::string) fopen(file_name, 'WRITE'); fclose(%); NULL end proc
 
-# interpolate.maple -- compute generalized interpolation formulas/coefficients
-# $Id: interpolate.maple,v 1.4 2002/05/14 15:52:50 jthorn Exp $
+# interpolate.maple -- compute interpolation formulas/coefficients
+# $Header: /cactusdevcvs/CactusBase/LocalInterp/src/GeneralizedPolynomial-Uniform/interpolate.maple,v 1.8 2002/08/20 16:46:05 jthorn Exp $
 > 
 #
 # <<<representation of numbers, data values, etc>>>
-# polynomial_interpolant - compute polynomial interpolant
+# Lagrange_polynomial_interpolant - compute Lagrange polynomial interpolant
+# Hermite_polynomial_interpolant - compute Hermite polynomial interpolant
 # coeff_as_lc_of_data - coefficients of ... (linear combination of data)
 #
 # print_coeff__lc_of_data - print C code to compute coefficients
@@ -487,7 +488,8 @@ ftruncate :=
 ################################################################################
 > 
 #
-# This function computes a polynomial interpolant in any number of dimensions.
+# This function computes a Lagrange polynomial interpolant in any
+# number of dimensions.
 #
 # Arguments:
 # fn = The interpolation function.  This should be a procedure in the
@@ -507,7 +509,7 @@ ftruncate :=
 # This function returns the interpolating polynomial, in the form of
 # an algebraic expression in the coordinates and the data values.
 #
-> polynomial_interpolant :=
+> Lagrange_polynomial_interpolant :=
 > proc(
 >       fn::procedure, coeff_list::list(name),
 >       coord_list::list(name), posn_list::list(list(numeric))
@@ -524,7 +526,7 @@ ftruncate :=
 # interpolant as a polynomial in the coordinates
 > return subs(coeff_eqns, eval(fn))(op(coord_list));
 > end proc;
-polynomial_interpolant := proc(fn::procedure, coeff_list::list(name),
+Lagrange_polynomial_interpolant := proc(fn::procedure, coeff_list::list(name),
 coord_list::list(name), posn_list::list(list(numeric)))
 local posn, data_eqns, coeff_eqns;
     data_eqns := {seq(fn(op(posn)) = 'DATA'(op(posn)), posn = posn_list)};
@@ -539,6 +541,171 @@ end proc
 ################################################################################
 > 
 #
+# This function computes a Hermite polynomial interpolant in any
+# number of dimensions.  This is a polynomial which
+#  has values which match the given data DATA() at a specified set of
+#   points, and
+# * has derivatives which match the specified finite-difference derivatives
+#   of the given data DATA() at a specified set of points
+#
+# For the derivative matching, we actually match all possible products
+# of 1st derivatives, i.e. in 2-D we match dx, dy, and dxy, in 3-D we
+# match dx, dy, dz, dxy, dxz, dyz, and dxyz, etc etc.
+#
+# Arguments:
+# fn = The interpolation function.  This should be a procedure in the
+#      coordinates, having the coefficients as global variables.  For
+#      example,
+#		proc(x,y)
+#		+ c03*y^3 + c13*x*y^3 + c23*x^2*y^3 + c33*x^3*y^3
+#		+ c02*y^2 + c12*x*y^2 + c22*x^2*y^2 + c32*x^3*y^2
+#		+ c01*y   + c11*x*y   + c21*x^2*y   + c31*x^3*y
+#		+ c00     + c10*x     + c20*x^2     + c30*x^3
+#		end proc;
+# coeff_list = A set of the interpolation coefficients (coefficients in
+#	       the interpolation function), for example
+#			[
+#			c03, c13, c23, c33,
+#			c02, c12, c22, c32,
+#			c01, c11, c21, c31,
+#			c00, c10, c20, c30
+#			]
+# coord_list = A list of the coordinates (independent variables in the
+#	       interpolation function), for example [x,y].
+# deriv_coord_list = A list of lists of coordinates specifying which
+#		     derivatives are computed by the deriv_proc_list[]
+#		     procedures, for example
+#			[[x], [y], [x,y]]
+# deriv_proc_list = A list of procedures for computing finite-difference
+#		    derivatives.  Each procedure should take N_dims integer
+#		    arguments specifying an evaluation point, and return
+#		    a suitable linear combination of the DATA() for the
+#		    derivative at that point.  For example
+#		    example,
+#			[
+#			  proc(i::integer, j::integer)
+#			  - 1/2*DATA(i-1,j) + 1/2*DATA(i+1,j)
+#			  end proc
+#			,
+#			  proc(i::integer, j::integer)
+#			  - 1/2*DATA(i,j-1) + 1/2*DATA(i,j+1)
+#			  end proc
+#			,
+#			  proc(i::integer, j::integer)
+#			  - 1/4*DATA(i-1,j+1) + 1/4*DATA(i+1,j+1)
+#			  + 1/4*DATA(i-1,j-1) - 1/4*DATA(i+1,j-1)
+#			  end proc
+#			]
+# fn_posn_list = A list of positions (each a list of numeric values)
+#		 where the interpolant is to match the given data DATA(),
+#		 for example
+#			[[0,0], [0,1], [1,0], [1,1]]
+# deriv_posn_list = A list of positions (each a list of numeric values)
+#		    where the interpolant is to match each possible product
+#		    of 1st derivatives"1st"
+#		    derivatives (actually all derivatives the given data DATA(),
+#		    for example
+#			[[0,0], [0,1], [1,0], [1,1]]
+#
+# Results:
+# This function returns the interpolating polynomial, in the form of
+# an algebraic expression in the coordinates and the data values.
+#
+> Hermite_polynomial_interpolant :=
+> proc(
+>       fn::procedure, coeff_list::list(name), coord_list::list(name),
+>       deriv_coord_list::list(list(name)), fd_deriv_proc_list::list(procedure),
+>       fn_posn_list::list(list(numeric)), deriv_posn_list::list(list(numeric))
+>     )
+> local fn_eset,
+>       fn_expr, deriv_eset;
+> 
+# set of equations {fn(posn) = DATA(posn)}
+> fn_eset := map(
+> 		  # return equation that fn(this point) = DATA(this point)
+> 		  proc(posn::list(integer))
+> 		  fn(op(posn)) = 'DATA'(op(posn));
+> 		  end proc
+> 		,
+> 		  {op(fn_posn_list)}
+> 	      );
+> 
+# set of sets of equations
+#	{
+#	  { deriv1(posn1) = fd_deriv1(posn1),
+#	    deriv2(posn1) = fd_deriv2(posn1),
+#	    ... },
+#	  { deriv1(posn2) = fd_deriv1(posn2),
+#	    deriv2(posn2) = fd_deriv2(posn2),
+#	    ... },
+#	  ...
+#	}
+> fn_expr := fn(op(coord_list));
+> map(
+>        # return set of equations
+>        #	{deriv(posn) = fd_deriv(posn)}
+>        # for this point
+>        proc(posn::list(integer))
+>        {
+> 	op(
+> 	 zip(
+> 		# return equation that
+> 		#	deriv(posn) = fd_deriv(posn)
+> 		# for this deriv and this point
+> 		proc(deriv_coords::list(name), fd_deriv_proc::procedure)
+> 		local fn_deriv_proc;
+> 		fn_deriv_proc := unapply( diff(fn_expr,deriv_coords),
+> 					  op(coord_list) );
+> 		fn_deriv_proc(op(posn)) = fd_deriv_proc(op(posn));
+> 		end proc
+> 	      ,
+> 		[op(deriv_coord_list)]
+> 	      ,
+> 		[op(fd_deriv_proc_list)]
+> 	    )
+> 	  )
+>        }
+>        end proc
+>      ,
+>        {op(deriv_posn_list)}
+>    );
+> 
+# set of equations {deriv-i(posn-j) = fd_deriv-i(posn-j)}
+> deriv_eset := `union`(op(%));
+> 
+# solve equations for coeffs
+> solve(fn_eset union deriv_eset, {op(coeff_list)});
+> 
+# interpolant as a polynomial in the coordinates
+> return subs(%, eval(fn))(op(coord_list));
+> end proc;
+Hermite_polynomial_interpolant := proc(fn::procedure, coeff_list::list(name),
+coord_list::list(name), deriv_coord_list::list(list(name)),
+fd_deriv_proc_list::list(procedure), fn_posn_list::list(list(numeric)),
+deriv_posn_list::list(list(numeric)))
+local fn_eset, fn_expr, deriv_eset;
+    fn_eset := map(
+        proc(posn::list(integer)) fn(op(posn)) = 'DATA'(op(posn)) end proc,
+        {op(fn_posn_list)});
+    fn_expr := fn(op(coord_list));
+    map(proc(posn::list(integer))
+            {op(zip(proc(deriv_coords::list(name), fd_deriv_proc::procedure
+            )
+            local fn_deriv_proc;
+                fn_deriv_proc :=
+                    unapply(diff(fn_expr, deriv_coords), op(coord_list));
+                fn_deriv_proc(op(posn)) = fd_deriv_proc(op(posn))
+            end proc, [op(deriv_coord_list)], [op(fd_deriv_proc_list)]))}
+        end proc, {op(deriv_posn_list)});
+    deriv_eset := `union`(op(%));
+    solve(fn_eset union deriv_eset, {op(coeff_list)});
+    return subs(%, eval(fn))(op(coord_list))
+end proc
+
+> 
+################################################################################
+> 
+#
 # This function takes as input an interpolating polynomial, expresses
 # it as a linear combination of the data values, and returns the coefficeints
 # of that form.
@@ -546,8 +713,7 @@ end proc
 # Arguments:
 # interpolant = The interpolating polynomial (an algebraic expression
 #		in the coordinates and the data values).
-# posn_list = The same list of positions as was used to compute the
-#	      interpolating polynomial.
+# posn_list = The same list of data positions used in the interpolant.
 #
 # Results:
 # This function returns the coefficients, as a list of equations of the
@@ -915,7 +1081,7 @@ data_var_name := proc(posn::list(numeric), name_prefix::string)
 end proc
 
 # Maple code to compute lists of point positions in hypercube-shaped molecules
-# $Id: $
+# $Id: cube_posns.maple,v 1.1.1.1 2002/08/18 15:05:38 jthorn Exp $
 > 
 ################################################################################
 > 
@@ -974,6 +1140,18 @@ posn_list_2d_size5 := [[-2, -2], [-1, -2], [0, -2], [1, -2], [2, -2], [-2, -1],
 
     [1, 2], [2, 2]]
 
+> posn_list_2d_size6 := map(ListTools[Reverse],
+> 			  hypercube_points([-2,-2], [+3,+3]));
+posn_list_2d_size6 := [[-2, -2], [-1, -2], [0, -2], [1, -2], [2, -2], [3, -2],
+
+    [-2, -1], [-1, -1], [0, -1], [1, -1], [2, -1], [3, -1], [-2, 0], [-1, 0],
+
+    [0, 0], [1, 0], [2, 0], [3, 0], [-2, 1], [-1, 1], [0, 1], [1, 1], [2, 1],
+
+    [3, 1], [-2, 2], [-1, 2], [0, 2], [1, 2], [2, 2], [3, 2], [-2, 3], [-1, 3],
+
+    [0, 3], [1, 3], [2, 3], [3, 3]]
+
 > 
 ################################################################################
 > 
@@ -1068,8 +1246,84 @@ posn_list_3d_size5 := [[-2, -2, -2], [-1, -2, -2], [0, -2, -2], [1, -2, -2],
 
     [0, 2, 2], [1, 2, 2], [2, 2, 2]]
 
+> posn_list_3d_size6 := map(ListTools[Reverse],
+> 			  hypercube_points([-2,-2,-2], [+3,+3,+3]));
+posn_list_3d_size6 := [[-2, -2, -2], [-1, -2, -2], [0, -2, -2], [1, -2, -2],
+
+    [2, -2, -2], [3, -2, -2], [-2, -1, -2], [-1, -1, -2], [0, -1, -2],
+
+    [1, -1, -2], [2, -1, -2], [3, -1, -2], [-2, 0, -2], [-1, 0, -2], [0, 0, -2],
+
+    [1, 0, -2], [2, 0, -2], [3, 0, -2], [-2, 1, -2], [-1, 1, -2], [0, 1, -2],
+
+    [1, 1, -2], [2, 1, -2], [3, 1, -2], [-2, 2, -2], [-1, 2, -2], [0, 2, -2],
+
+    [1, 2, -2], [2, 2, -2], [3, 2, -2], [-2, 3, -2], [-1, 3, -2], [0, 3, -2],
+
+    [1, 3, -2], [2, 3, -2], [3, 3, -2], [-2, -2, -1], [-1, -2, -1], [0, -2, -1],
+
+    [1, -2, -1], [2, -2, -1], [3, -2, -1], [-2, -1, -1], [-1, -1, -1],
+
+    [0, -1, -1], [1, -1, -1], [2, -1, -1], [3, -1, -1], [-2, 0, -1],
+
+    [-1, 0, -1], [0, 0, -1], [1, 0, -1], [2, 0, -1], [3, 0, -1], [-2, 1, -1],
+
+    [-1, 1, -1], [0, 1, -1], [1, 1, -1], [2, 1, -1], [3, 1, -1], [-2, 2, -1],
+
+    [-1, 2, -1], [0, 2, -1], [1, 2, -1], [2, 2, -1], [3, 2, -1], [-2, 3, -1],
+
+    [-1, 3, -1], [0, 3, -1], [1, 3, -1], [2, 3, -1], [3, 3, -1], [-2, -2, 0],
+
+    [-1, -2, 0], [0, -2, 0], [1, -2, 0], [2, -2, 0], [3, -2, 0], [-2, -1, 0],
+
+    [-1, -1, 0], [0, -1, 0], [1, -1, 0], [2, -1, 0], [3, -1, 0], [-2, 0, 0],
+
+    [-1, 0, 0], [0, 0, 0], [1, 0, 0], [2, 0, 0], [3, 0, 0], [-2, 1, 0],
+
+    [-1, 1, 0], [0, 1, 0], [1, 1, 0], [2, 1, 0], [3, 1, 0], [-2, 2, 0],
+
+    [-1, 2, 0], [0, 2, 0], [1, 2, 0], [2, 2, 0], [3, 2, 0], [-2, 3, 0],
+
+    [-1, 3, 0], [0, 3, 0], [1, 3, 0], [2, 3, 0], [3, 3, 0], [-2, -2, 1],
+
+    [-1, -2, 1], [0, -2, 1], [1, -2, 1], [2, -2, 1], [3, -2, 1], [-2, -1, 1],
+
+    [-1, -1, 1], [0, -1, 1], [1, -1, 1], [2, -1, 1], [3, -1, 1], [-2, 0, 1],
+
+    [-1, 0, 1], [0, 0, 1], [1, 0, 1], [2, 0, 1], [3, 0, 1], [-2, 1, 1],
+
+    [-1, 1, 1], [0, 1, 1], [1, 1, 1], [2, 1, 1], [3, 1, 1], [-2, 2, 1],
+
+    [-1, 2, 1], [0, 2, 1], [1, 2, 1], [2, 2, 1], [3, 2, 1], [-2, 3, 1],
+
+    [-1, 3, 1], [0, 3, 1], [1, 3, 1], [2, 3, 1], [3, 3, 1], [-2, -2, 2],
+
+    [-1, -2, 2], [0, -2, 2], [1, -2, 2], [2, -2, 2], [3, -2, 2], [-2, -1, 2],
+
+    [-1, -1, 2], [0, -1, 2], [1, -1, 2], [2, -1, 2], [3, -1, 2], [-2, 0, 2],
+
+    [-1, 0, 2], [0, 0, 2], [1, 0, 2], [2, 0, 2], [3, 0, 2], [-2, 1, 2],
+
+    [-1, 1, 2], [0, 1, 2], [1, 1, 2], [2, 1, 2], [3, 1, 2], [-2, 2, 2],
+
+    [-1, 2, 2], [0, 2, 2], [1, 2, 2], [2, 2, 2], [3, 2, 2], [-2, 3, 2],
+
+    [-1, 3, 2], [0, 3, 2], [1, 3, 2], [2, 3, 2], [3, 3, 2], [-2, -2, 3],
+
+    [-1, -2, 3], [0, -2, 3], [1, -2, 3], [2, -2, 3], [3, -2, 3], [-2, -1, 3],
+
+    [-1, -1, 3], [0, -1, 3], [1, -1, 3], [2, -1, 3], [3, -1, 3], [-2, 0, 3],
+
+    [-1, 0, 3], [0, 0, 3], [1, 0, 3], [2, 0, 3], [3, 0, 3], [-2, 1, 3],
+
+    [-1, 1, 3], [0, 1, 3], [1, 1, 3], [2, 1, 3], [3, 1, 3], [-2, 2, 3],
+
+    [-1, 2, 3], [0, 2, 3], [1, 2, 3], [2, 2, 3], [3, 2, 3], [-2, 3, 3],
+
+    [-1, 3, 3], [0, 3, 3], [1, 3, 3], [2, 3, 3], [3, 3, 3]]
+
 # Maple code to compute common coefficients for all 2d interpolation schemes
-# $Id: 2d.maple,v 1.3 2002/05/14 15:54:01 jthorn Exp $
+# $Id: 2d.maple,v 1.1.1.1 2002/08/18 15:05:38 jthorn Exp $
 > 
 ################################################################################
 > 
@@ -1218,6 +1472,7 @@ data_var_list_2d_size4 := ["data_m1_m1", "data_0_m1", "data_p1_m1",
 > 		    "2d.cube.size4/coeff-dxy.dcl.c");
 > print_name_list_dcl(map(coeff_name, posn_list_2d_size4, "coeff_dyy_"), "fp",
 > 		    "2d.cube.size4/coeff-dyy.dcl.c");
+bytes used=1000212, alloc=917336, time=0.13
 > 
 > print_interp_cmpt__lc_of_data(posn_list_2d_size4,
 > 			      "result", "coeff_I_", "data_",
@@ -1228,7 +1483,6 @@ data_var_list_2d_size4 := ["data_m1_m1", "data_0_m1", "data_p1_m1",
 > print_interp_cmpt__lc_of_data(posn_list_2d_size4,
 > 			      "result", "coeff_dy_", "data_",
 > 			      "2d.cube.size4/interp-dy.compute.c");
-bytes used=1000024, alloc=917336, time=0.17
 > print_interp_cmpt__lc_of_data(posn_list_2d_size4,
 > 			      "result", "coeff_dxx_", "data_",
 > 			      "2d.cube.size4/interp-dxx.compute.c");
@@ -1310,5 +1564,79 @@ data_var_list_2d_size5 := ["data_m2_m2", "data_m1_m2", "data_0_m2",
 > 			      "2d.cube.size5/interp-dyy.compute.c");
 > 
 ################################################################################
+#
+# generic stuff for 2d, cube, size=5
+#
+> 
+> data_var_list_2d_size6 := map(data_var_name, posn_list_2d_size6, "data_");
+data_var_list_2d_size6 := ["data_m2_m2", "data_m1_m2", "data_0_m2",
+
+    "data_p1_m2", "data_p2_m2", "data_p3_m2", "data_m2_m1", "data_m1_m1",
+
+    "data_0_m1", "data_p1_m1", "data_p2_m1", "data_p3_m1", "data_m2_0",
+
+    "data_m1_0", "data_0_0", "data_p1_0", "data_p2_0", "data_p3_0",
+
+    "data_m2_p1", "data_m1_p1", "data_0_p1", "data_p1_p1", "data_p2_p1",
+
+    "data_p3_p1", "data_m2_p2", "data_m1_p2", "data_0_p2", "data_p1_p2",
+
+    "data_p2_p2", "data_p3_p2", "data_m2_p3", "data_m1_p3", "data_0_p3",
+
+    "data_p1_p3", "data_p2_p3", "data_p3_p3"]
+
+> 
+> print_name_list_dcl(data_var_list_2d_size6, "fp",
+> 		    "2d.cube.size6/data-var.dcl.c");
+> print_data_var_assign(posn_list_2d_size6, "data_",
+> 		      "2d.cube.size6/data-var.assign.c");
+> 
+> print_interp_coeff_var_store(posn_list_2d_size6, "", "coeff_I_",
+> 			     "2d.cube.size6/coeff-I.store.c");
+> print_interp_coeff_var_store(posn_list_2d_size6, "factor", "coeff_dx_",
+> 			     "2d.cube.size6/coeff-dx.store.c");
+> print_interp_coeff_var_store(posn_list_2d_size6, "factor", "coeff_dy_",
+> 			     "2d.cube.size6/coeff-dy.store.c");
+> print_interp_coeff_var_store(posn_list_2d_size6, "factor", "coeff_dxx_",
+> 			     "2d.cube.size6/coeff-dxx.store.c");
+bytes used=2000440, alloc=1179432, time=0.27
+> print_interp_coeff_var_store(posn_list_2d_size6, "factor", "coeff_dxy_",
+> 			     "2d.cube.size6/coeff-dxy.store.c");
+> print_interp_coeff_var_store(posn_list_2d_size6, "factor", "coeff_dyy_",
+> 			     "2d.cube.size6/coeff-dyy.store.c");
+> 
+> print_name_list_dcl(map(coeff_name, posn_list_2d_size6, "coeff_I_"), "fp",
+> 		    "2d.cube.size6/coeff-I.dcl.c");
+> print_name_list_dcl(map(coeff_name, posn_list_2d_size6, "coeff_dx_"), "fp",
+> 		    "2d.cube.size6/coeff-dx.dcl.c");
+> print_name_list_dcl(map(coeff_name, posn_list_2d_size6, "coeff_dy_"), "fp",
+> 		    "2d.cube.size6/coeff-dy.dcl.c");
+> print_name_list_dcl(map(coeff_name, posn_list_2d_size6, "coeff_dxx_"), "fp",
+> 		    "2d.cube.size6/coeff-dxx.dcl.c");
+> print_name_list_dcl(map(coeff_name, posn_list_2d_size6, "coeff_dxy_"), "fp",
+> 		    "2d.cube.size6/coeff-dxy.dcl.c");
+> print_name_list_dcl(map(coeff_name, posn_list_2d_size6, "coeff_dyy_"), "fp",
+> 		    "2d.cube.size6/coeff-dyy.dcl.c");
+> 
+> print_interp_cmpt__lc_of_data(posn_list_2d_size6,
+> 			      "result", "coeff_I_", "data_",
+> 			      "2d.cube.size6/interp-I.compute.c");
+> print_interp_cmpt__lc_of_data(posn_list_2d_size6,
+> 			      "result", "coeff_dx_", "data_",
+> 			      "2d.cube.size6/interp-dx.compute.c");
+> print_interp_cmpt__lc_of_data(posn_list_2d_size6,
+> 			      "result", "coeff_dy_", "data_",
+> 			      "2d.cube.size6/interp-dy.compute.c");
+> print_interp_cmpt__lc_of_data(posn_list_2d_size6,
+> 			      "result", "coeff_dxx_", "data_",
+> 			      "2d.cube.size6/interp-dxx.compute.c");
+> print_interp_cmpt__lc_of_data(posn_list_2d_size6,
+> 			      "result", "coeff_dxy_", "data_",
+> 			      "2d.cube.size6/interp-dxy.compute.c");
+> print_interp_cmpt__lc_of_data(posn_list_2d_size6,
+> 			      "result", "coeff_dyy_", "data_",
+> 			      "2d.cube.size6/interp-dyy.compute.c");
+> 
+################################################################################
 > quit
-bytes used=1720848, alloc=917336, time=0.31
+bytes used=2688044, alloc=1179432, time=0.39
diff --git a/src/GeneralizedPolynomial-Uniform/common/2d.maple b/src/GeneralizedPolynomial-Uniform/common/2d.maple
index 4a7b939..6efcc67 100644
--- a/src/GeneralizedPolynomial-Uniform/common/2d.maple
+++ b/src/GeneralizedPolynomial-Uniform/common/2d.maple
@@ -213,3 +213,60 @@ print_interp_cmpt__lc_of_data(posn_list_2d_size5,
 			      "2d.cube.size5/interp-dyy.compute.c");
 
 ################################################################################
+#
+# generic stuff for 2d, cube, size=5
+#
+
+data_var_list_2d_size6 := map(data_var_name, posn_list_2d_size6, "data_");
+
+print_name_list_dcl(data_var_list_2d_size6, "fp",
+		    "2d.cube.size6/data-var.dcl.c");
+print_data_var_assign(posn_list_2d_size6, "data_",
+		      "2d.cube.size6/data-var.assign.c");
+
+print_interp_coeff_var_store(posn_list_2d_size6, "", "coeff_I_",
+			     "2d.cube.size6/coeff-I.store.c");
+print_interp_coeff_var_store(posn_list_2d_size6, "factor", "coeff_dx_",
+			     "2d.cube.size6/coeff-dx.store.c");
+print_interp_coeff_var_store(posn_list_2d_size6, "factor", "coeff_dy_",
+			     "2d.cube.size6/coeff-dy.store.c");
+print_interp_coeff_var_store(posn_list_2d_size6, "factor", "coeff_dxx_",
+			     "2d.cube.size6/coeff-dxx.store.c");
+print_interp_coeff_var_store(posn_list_2d_size6, "factor", "coeff_dxy_",
+			     "2d.cube.size6/coeff-dxy.store.c");
+print_interp_coeff_var_store(posn_list_2d_size6, "factor", "coeff_dyy_",
+			     "2d.cube.size6/coeff-dyy.store.c");
+
+print_name_list_dcl(map(coeff_name, posn_list_2d_size6, "coeff_I_"), "fp",
+		    "2d.cube.size6/coeff-I.dcl.c");
+print_name_list_dcl(map(coeff_name, posn_list_2d_size6, "coeff_dx_"), "fp",
+		    "2d.cube.size6/coeff-dx.dcl.c");
+print_name_list_dcl(map(coeff_name, posn_list_2d_size6, "coeff_dy_"), "fp",
+		    "2d.cube.size6/coeff-dy.dcl.c");
+print_name_list_dcl(map(coeff_name, posn_list_2d_size6, "coeff_dxx_"), "fp",
+		    "2d.cube.size6/coeff-dxx.dcl.c");
+print_name_list_dcl(map(coeff_name, posn_list_2d_size6, "coeff_dxy_"), "fp",
+		    "2d.cube.size6/coeff-dxy.dcl.c");
+print_name_list_dcl(map(coeff_name, posn_list_2d_size6, "coeff_dyy_"), "fp",
+		    "2d.cube.size6/coeff-dyy.dcl.c");
+
+print_interp_cmpt__lc_of_data(posn_list_2d_size6,
+			      "result", "coeff_I_", "data_",
+			      "2d.cube.size6/interp-I.compute.c");
+print_interp_cmpt__lc_of_data(posn_list_2d_size6,
+			      "result", "coeff_dx_", "data_",
+			      "2d.cube.size6/interp-dx.compute.c");
+print_interp_cmpt__lc_of_data(posn_list_2d_size6,
+			      "result", "coeff_dy_", "data_",
+			      "2d.cube.size6/interp-dy.compute.c");
+print_interp_cmpt__lc_of_data(posn_list_2d_size6,
+			      "result", "coeff_dxx_", "data_",
+			      "2d.cube.size6/interp-dxx.compute.c");
+print_interp_cmpt__lc_of_data(posn_list_2d_size6,
+			      "result", "coeff_dxy_", "data_",
+			      "2d.cube.size6/interp-dxy.compute.c");
+print_interp_cmpt__lc_of_data(posn_list_2d_size6,
+			      "result", "coeff_dyy_", "data_",
+			      "2d.cube.size6/interp-dyy.compute.c");
+
+################################################################################
diff --git a/src/GeneralizedPolynomial-Uniform/common/3d.log b/src/GeneralizedPolynomial-Uniform/common/3d.log
index 2fed690..12f76ca 100644
--- a/src/GeneralizedPolynomial-Uniform/common/3d.log
+++ b/src/GeneralizedPolynomial-Uniform/common/3d.log
@@ -4,7 +4,7 @@
  <____ ____>  Waterloo Maple Inc.
       |       Type ? for help.
 # util.maple -- misc utility routines
-# $Id: util.maple,v 1.3 2002/05/19 13:12:18 jthorn Exp $
+# $Header: /cactusdevcvs/CactusBase/LocalInterp/src/GeneralizedPolynomial-Uniform/util.maple,v 1.4 2002/08/20 16:46:06 jthorn Exp $
 > 
 #
 # fix_rationals - convert numbers to RATIONAL() calls
@@ -444,12 +444,13 @@ ftruncate :=
 
     proc(file_name::string) fopen(file_name, 'WRITE'); fclose(%); NULL end proc
 
-# interpolate.maple -- compute generalized interpolation formulas/coefficients
-# $Id: interpolate.maple,v 1.4 2002/05/14 15:52:50 jthorn Exp $
+# interpolate.maple -- compute interpolation formulas/coefficients
+# $Header: /cactusdevcvs/CactusBase/LocalInterp/src/GeneralizedPolynomial-Uniform/interpolate.maple,v 1.8 2002/08/20 16:46:05 jthorn Exp $
 > 
 #
 # <<<representation of numbers, data values, etc>>>
-# polynomial_interpolant - compute polynomial interpolant
+# Lagrange_polynomial_interpolant - compute Lagrange polynomial interpolant
+# Hermite_polynomial_interpolant - compute Hermite polynomial interpolant
 # coeff_as_lc_of_data - coefficients of ... (linear combination of data)
 #
 # print_coeff__lc_of_data - print C code to compute coefficients
@@ -487,7 +488,8 @@ ftruncate :=
 ################################################################################
 > 
 #
-# This function computes a polynomial interpolant in any number of dimensions.
+# This function computes a Lagrange polynomial interpolant in any
+# number of dimensions.
 #
 # Arguments:
 # fn = The interpolation function.  This should be a procedure in the
@@ -507,7 +509,7 @@ ftruncate :=
 # This function returns the interpolating polynomial, in the form of
 # an algebraic expression in the coordinates and the data values.
 #
-> polynomial_interpolant :=
+> Lagrange_polynomial_interpolant :=
 > proc(
 >       fn::procedure, coeff_list::list(name),
 >       coord_list::list(name), posn_list::list(list(numeric))
@@ -524,7 +526,7 @@ ftruncate :=
 # interpolant as a polynomial in the coordinates
 > return subs(coeff_eqns, eval(fn))(op(coord_list));
 > end proc;
-polynomial_interpolant := proc(fn::procedure, coeff_list::list(name),
+Lagrange_polynomial_interpolant := proc(fn::procedure, coeff_list::list(name),
 coord_list::list(name), posn_list::list(list(numeric)))
 local posn, data_eqns, coeff_eqns;
     data_eqns := {seq(fn(op(posn)) = 'DATA'(op(posn)), posn = posn_list)};
@@ -539,6 +541,171 @@ end proc
 ################################################################################
 > 
 #
+# This function computes a Hermite polynomial interpolant in any
+# number of dimensions.  This is a polynomial which
+#  has values which match the given data DATA() at a specified set of
+#   points, and
+# * has derivatives which match the specified finite-difference derivatives
+#   of the given data DATA() at a specified set of points
+#
+# For the derivative matching, we actually match all possible products
+# of 1st derivatives, i.e. in 2-D we match dx, dy, and dxy, in 3-D we
+# match dx, dy, dz, dxy, dxz, dyz, and dxyz, etc etc.
+#
+# Arguments:
+# fn = The interpolation function.  This should be a procedure in the
+#      coordinates, having the coefficients as global variables.  For
+#      example,
+#		proc(x,y)
+#		+ c03*y^3 + c13*x*y^3 + c23*x^2*y^3 + c33*x^3*y^3
+#		+ c02*y^2 + c12*x*y^2 + c22*x^2*y^2 + c32*x^3*y^2
+#		+ c01*y   + c11*x*y   + c21*x^2*y   + c31*x^3*y
+#		+ c00     + c10*x     + c20*x^2     + c30*x^3
+#		end proc;
+# coeff_list = A set of the interpolation coefficients (coefficients in
+#	       the interpolation function), for example
+#			[
+#			c03, c13, c23, c33,
+#			c02, c12, c22, c32,
+#			c01, c11, c21, c31,
+#			c00, c10, c20, c30
+#			]
+# coord_list = A list of the coordinates (independent variables in the
+#	       interpolation function), for example [x,y].
+# deriv_coord_list = A list of lists of coordinates specifying which
+#		     derivatives are computed by the deriv_proc_list[]
+#		     procedures, for example
+#			[[x], [y], [x,y]]
+# deriv_proc_list = A list of procedures for computing finite-difference
+#		    derivatives.  Each procedure should take N_dims integer
+#		    arguments specifying an evaluation point, and return
+#		    a suitable linear combination of the DATA() for the
+#		    derivative at that point.  For example
+#		    example,
+#			[
+#			  proc(i::integer, j::integer)
+#			  - 1/2*DATA(i-1,j) + 1/2*DATA(i+1,j)
+#			  end proc
+#			,
+#			  proc(i::integer, j::integer)
+#			  - 1/2*DATA(i,j-1) + 1/2*DATA(i,j+1)
+#			  end proc
+#			,
+#			  proc(i::integer, j::integer)
+#			  - 1/4*DATA(i-1,j+1) + 1/4*DATA(i+1,j+1)
+#			  + 1/4*DATA(i-1,j-1) - 1/4*DATA(i+1,j-1)
+#			  end proc
+#			]
+# fn_posn_list = A list of positions (each a list of numeric values)
+#		 where the interpolant is to match the given data DATA(),
+#		 for example
+#			[[0,0], [0,1], [1,0], [1,1]]
+# deriv_posn_list = A list of positions (each a list of numeric values)
+#		    where the interpolant is to match each possible product
+#		    of 1st derivatives"1st"
+#		    derivatives (actually all derivatives the given data DATA(),
+#		    for example
+#			[[0,0], [0,1], [1,0], [1,1]]
+#
+# Results:
+# This function returns the interpolating polynomial, in the form of
+# an algebraic expression in the coordinates and the data values.
+#
+> Hermite_polynomial_interpolant :=
+> proc(
+>       fn::procedure, coeff_list::list(name), coord_list::list(name),
+>       deriv_coord_list::list(list(name)), fd_deriv_proc_list::list(procedure),
+>       fn_posn_list::list(list(numeric)), deriv_posn_list::list(list(numeric))
+>     )
+> local fn_eset,
+>       fn_expr, deriv_eset;
+> 
+# set of equations {fn(posn) = DATA(posn)}
+> fn_eset := map(
+> 		  # return equation that fn(this point) = DATA(this point)
+> 		  proc(posn::list(integer))
+> 		  fn(op(posn)) = 'DATA'(op(posn));
+> 		  end proc
+> 		,
+> 		  {op(fn_posn_list)}
+> 	      );
+> 
+# set of sets of equations
+#	{
+#	  { deriv1(posn1) = fd_deriv1(posn1),
+#	    deriv2(posn1) = fd_deriv2(posn1),
+#	    ... },
+#	  { deriv1(posn2) = fd_deriv1(posn2),
+#	    deriv2(posn2) = fd_deriv2(posn2),
+#	    ... },
+#	  ...
+#	}
+> fn_expr := fn(op(coord_list));
+> map(
+>        # return set of equations
+>        #	{deriv(posn) = fd_deriv(posn)}
+>        # for this point
+>        proc(posn::list(integer))
+>        {
+> 	op(
+> 	 zip(
+> 		# return equation that
+> 		#	deriv(posn) = fd_deriv(posn)
+> 		# for this deriv and this point
+> 		proc(deriv_coords::list(name), fd_deriv_proc::procedure)
+> 		local fn_deriv_proc;
+> 		fn_deriv_proc := unapply( diff(fn_expr,deriv_coords),
+> 					  op(coord_list) );
+> 		fn_deriv_proc(op(posn)) = fd_deriv_proc(op(posn));
+> 		end proc
+> 	      ,
+> 		[op(deriv_coord_list)]
+> 	      ,
+> 		[op(fd_deriv_proc_list)]
+> 	    )
+> 	  )
+>        }
+>        end proc
+>      ,
+>        {op(deriv_posn_list)}
+>    );
+> 
+# set of equations {deriv-i(posn-j) = fd_deriv-i(posn-j)}
+> deriv_eset := `union`(op(%));
+> 
+# solve equations for coeffs
+> solve(fn_eset union deriv_eset, {op(coeff_list)});
+> 
+# interpolant as a polynomial in the coordinates
+> return subs(%, eval(fn))(op(coord_list));
+> end proc;
+Hermite_polynomial_interpolant := proc(fn::procedure, coeff_list::list(name),
+coord_list::list(name), deriv_coord_list::list(list(name)),
+fd_deriv_proc_list::list(procedure), fn_posn_list::list(list(numeric)),
+deriv_posn_list::list(list(numeric)))
+local fn_eset, fn_expr, deriv_eset;
+    fn_eset := map(
+        proc(posn::list(integer)) fn(op(posn)) = 'DATA'(op(posn)) end proc,
+        {op(fn_posn_list)});
+    fn_expr := fn(op(coord_list));
+    map(proc(posn::list(integer))
+            {op(zip(proc(deriv_coords::list(name), fd_deriv_proc::procedure
+            )
+            local fn_deriv_proc;
+                fn_deriv_proc :=
+                    unapply(diff(fn_expr, deriv_coords), op(coord_list));
+                fn_deriv_proc(op(posn)) = fd_deriv_proc(op(posn))
+            end proc, [op(deriv_coord_list)], [op(fd_deriv_proc_list)]))}
+        end proc, {op(deriv_posn_list)});
+    deriv_eset := `union`(op(%));
+    solve(fn_eset union deriv_eset, {op(coeff_list)});
+    return subs(%, eval(fn))(op(coord_list))
+end proc
+
+> 
+################################################################################
+> 
+#
 # This function takes as input an interpolating polynomial, expresses
 # it as a linear combination of the data values, and returns the coefficeints
 # of that form.
@@ -546,8 +713,7 @@ end proc
 # Arguments:
 # interpolant = The interpolating polynomial (an algebraic expression
 #		in the coordinates and the data values).
-# posn_list = The same list of positions as was used to compute the
-#	      interpolating polynomial.
+# posn_list = The same list of data positions used in the interpolant.
 #
 # Results:
 # This function returns the coefficients, as a list of equations of the
@@ -915,7 +1081,7 @@ data_var_name := proc(posn::list(numeric), name_prefix::string)
 end proc
 
 # Maple code to compute lists of point positions in hypercube-shaped molecules
-# $Id: $
+# $Id: cube_posns.maple,v 1.1.1.1 2002/08/18 15:05:38 jthorn Exp $
 > 
 ################################################################################
 > 
@@ -974,6 +1140,18 @@ posn_list_2d_size5 := [[-2, -2], [-1, -2], [0, -2], [1, -2], [2, -2], [-2, -1],
 
     [1, 2], [2, 2]]
 
+> posn_list_2d_size6 := map(ListTools[Reverse],
+> 			  hypercube_points([-2,-2], [+3,+3]));
+posn_list_2d_size6 := [[-2, -2], [-1, -2], [0, -2], [1, -2], [2, -2], [3, -2],
+
+    [-2, -1], [-1, -1], [0, -1], [1, -1], [2, -1], [3, -1], [-2, 0], [-1, 0],
+
+    [0, 0], [1, 0], [2, 0], [3, 0], [-2, 1], [-1, 1], [0, 1], [1, 1], [2, 1],
+
+    [3, 1], [-2, 2], [-1, 2], [0, 2], [1, 2], [2, 2], [3, 2], [-2, 3], [-1, 3],
+
+    [0, 3], [1, 3], [2, 3], [3, 3]]
+
 > 
 ################################################################################
 > 
@@ -1068,8 +1246,84 @@ posn_list_3d_size5 := [[-2, -2, -2], [-1, -2, -2], [0, -2, -2], [1, -2, -2],
 
     [0, 2, 2], [1, 2, 2], [2, 2, 2]]
 
+> posn_list_3d_size6 := map(ListTools[Reverse],
+> 			  hypercube_points([-2,-2,-2], [+3,+3,+3]));
+posn_list_3d_size6 := [[-2, -2, -2], [-1, -2, -2], [0, -2, -2], [1, -2, -2],
+
+    [2, -2, -2], [3, -2, -2], [-2, -1, -2], [-1, -1, -2], [0, -1, -2],
+
+    [1, -1, -2], [2, -1, -2], [3, -1, -2], [-2, 0, -2], [-1, 0, -2], [0, 0, -2],
+
+    [1, 0, -2], [2, 0, -2], [3, 0, -2], [-2, 1, -2], [-1, 1, -2], [0, 1, -2],
+
+    [1, 1, -2], [2, 1, -2], [3, 1, -2], [-2, 2, -2], [-1, 2, -2], [0, 2, -2],
+
+    [1, 2, -2], [2, 2, -2], [3, 2, -2], [-2, 3, -2], [-1, 3, -2], [0, 3, -2],
+
+    [1, 3, -2], [2, 3, -2], [3, 3, -2], [-2, -2, -1], [-1, -2, -1], [0, -2, -1],
+
+    [1, -2, -1], [2, -2, -1], [3, -2, -1], [-2, -1, -1], [-1, -1, -1],
+
+    [0, -1, -1], [1, -1, -1], [2, -1, -1], [3, -1, -1], [-2, 0, -1],
+
+    [-1, 0, -1], [0, 0, -1], [1, 0, -1], [2, 0, -1], [3, 0, -1], [-2, 1, -1],
+
+    [-1, 1, -1], [0, 1, -1], [1, 1, -1], [2, 1, -1], [3, 1, -1], [-2, 2, -1],
+
+    [-1, 2, -1], [0, 2, -1], [1, 2, -1], [2, 2, -1], [3, 2, -1], [-2, 3, -1],
+
+    [-1, 3, -1], [0, 3, -1], [1, 3, -1], [2, 3, -1], [3, 3, -1], [-2, -2, 0],
+
+    [-1, -2, 0], [0, -2, 0], [1, -2, 0], [2, -2, 0], [3, -2, 0], [-2, -1, 0],
+
+    [-1, -1, 0], [0, -1, 0], [1, -1, 0], [2, -1, 0], [3, -1, 0], [-2, 0, 0],
+
+    [-1, 0, 0], [0, 0, 0], [1, 0, 0], [2, 0, 0], [3, 0, 0], [-2, 1, 0],
+
+    [-1, 1, 0], [0, 1, 0], [1, 1, 0], [2, 1, 0], [3, 1, 0], [-2, 2, 0],
+
+    [-1, 2, 0], [0, 2, 0], [1, 2, 0], [2, 2, 0], [3, 2, 0], [-2, 3, 0],
+
+    [-1, 3, 0], [0, 3, 0], [1, 3, 0], [2, 3, 0], [3, 3, 0], [-2, -2, 1],
+
+    [-1, -2, 1], [0, -2, 1], [1, -2, 1], [2, -2, 1], [3, -2, 1], [-2, -1, 1],
+
+    [-1, -1, 1], [0, -1, 1], [1, -1, 1], [2, -1, 1], [3, -1, 1], [-2, 0, 1],
+
+    [-1, 0, 1], [0, 0, 1], [1, 0, 1], [2, 0, 1], [3, 0, 1], [-2, 1, 1],
+
+    [-1, 1, 1], [0, 1, 1], [1, 1, 1], [2, 1, 1], [3, 1, 1], [-2, 2, 1],
+
+    [-1, 2, 1], [0, 2, 1], [1, 2, 1], [2, 2, 1], [3, 2, 1], [-2, 3, 1],
+
+    [-1, 3, 1], [0, 3, 1], [1, 3, 1], [2, 3, 1], [3, 3, 1], [-2, -2, 2],
+
+    [-1, -2, 2], [0, -2, 2], [1, -2, 2], [2, -2, 2], [3, -2, 2], [-2, -1, 2],
+
+    [-1, -1, 2], [0, -1, 2], [1, -1, 2], [2, -1, 2], [3, -1, 2], [-2, 0, 2],
+
+    [-1, 0, 2], [0, 0, 2], [1, 0, 2], [2, 0, 2], [3, 0, 2], [-2, 1, 2],
+
+    [-1, 1, 2], [0, 1, 2], [1, 1, 2], [2, 1, 2], [3, 1, 2], [-2, 2, 2],
+
+    [-1, 2, 2], [0, 2, 2], [1, 2, 2], [2, 2, 2], [3, 2, 2], [-2, 3, 2],
+
+    [-1, 3, 2], [0, 3, 2], [1, 3, 2], [2, 3, 2], [3, 3, 2], [-2, -2, 3],
+
+    [-1, -2, 3], [0, -2, 3], [1, -2, 3], [2, -2, 3], [3, -2, 3], [-2, -1, 3],
+
+    [-1, -1, 3], [0, -1, 3], [1, -1, 3], [2, -1, 3], [3, -1, 3], [-2, 0, 3],
+
+    [-1, 0, 3], [0, 0, 3], [1, 0, 3], [2, 0, 3], [3, 0, 3], [-2, 1, 3],
+
+    [-1, 1, 3], [0, 1, 3], [1, 1, 3], [2, 1, 3], [3, 1, 3], [-2, 2, 3],
+
+    [-1, 2, 3], [0, 2, 3], [1, 2, 3], [2, 2, 3], [3, 2, 3], [-2, 3, 3],
+
+    [-1, 3, 3], [0, 3, 3], [1, 3, 3], [2, 3, 3], [3, 3, 3]]
+
 # Maple code to compute common coefficients for all 3d interpolation schemes
-# $Id: 3d.maple,v 1.3 2002/05/14 15:54:01 jthorn Exp $
+# $Id: 3d.maple,v 1.1.1.1 2002/08/18 15:05:38 jthorn Exp $
 > 
 ################################################################################
 > 
@@ -1168,12 +1422,12 @@ data_var_list_3d_size3 := ["data_m1_m1_m1", "data_0_m1_m1", "data_p1_m1_m1",
 > 			     "3d.cube.size3/coeff-dyz.store.c");
 > print_interp_coeff_var_store(posn_list_3d_size3, "factor", "coeff_dzz_",
 > 			     "3d.cube.size3/coeff-dzz.store.c");
+bytes used=1000020, alloc=917336, time=0.16
 > 
 > print_name_list_dcl(map(coeff_name, posn_list_3d_size3, "coeff_I_"), "fp",
 > 		    "3d.cube.size3/coeff-I.dcl.c");
 > print_name_list_dcl(map(coeff_name, posn_list_3d_size3, "coeff_dx_"), "fp",
 > 		    "3d.cube.size3/coeff-dx.dcl.c");
-bytes used=1000120, alloc=917336, time=0.18
 > print_name_list_dcl(map(coeff_name, posn_list_3d_size3, "coeff_dy_"), "fp",
 > 		    "3d.cube.size3/coeff-dy.dcl.c");
 > print_name_list_dcl(map(coeff_name, posn_list_3d_size3, "coeff_dz_"), "fp",
@@ -1229,6 +1483,7 @@ bytes used=1000120, alloc=917336, time=0.18
 #
 > 
 > data_var_list_3d_size4 := map(data_var_name, posn_list_3d_size4, "data_");
+bytes used=2000252, alloc=1179432, time=0.25
 data_var_list_3d_size4 := ["data_m1_m1_m1", "data_0_m1_m1", "data_p1_m1_m1",
 
     "data_p2_m1_m1", "data_m1_0_m1", "data_0_0_m1", "data_p1_0_m1",
@@ -1266,7 +1521,6 @@ data_var_list_3d_size4 := ["data_m1_m1_m1", "data_0_m1_m1", "data_p1_m1_m1",
 > 		    "3d.cube.size4/data-var.dcl.c");
 > print_data_var_assign(posn_list_3d_size4, "data_",
 > 		      "3d.cube.size4/data-var.assign.c");
-bytes used=2000380, alloc=1179432, time=0.28
 > 
 > print_interp_coeff_var_store(posn_list_3d_size4, "", "coeff_I_",
 > 			     "3d.cube.size4/coeff-I.store.c");
@@ -1293,9 +1547,9 @@ bytes used=2000380, alloc=1179432, time=0.28
 > 		    "3d.cube.size4/coeff-I.dcl.c");
 > print_name_list_dcl(map(coeff_name, posn_list_3d_size4, "coeff_dx_"), "fp",
 > 		    "3d.cube.size4/coeff-dx.dcl.c");
+bytes used=3000424, alloc=1244956, time=0.36
 > print_name_list_dcl(map(coeff_name, posn_list_3d_size4, "coeff_dy_"), "fp",
 > 		    "3d.cube.size4/coeff-dy.dcl.c");
-bytes used=3000632, alloc=1244956, time=0.40
 > print_name_list_dcl(map(coeff_name, posn_list_3d_size4, "coeff_dz_"), "fp",
 > 		    "3d.cube.size4/coeff-dz.dcl.c");
 > print_name_list_dcl(map(coeff_name, posn_list_3d_size4, "coeff_dxx_"), "fp",
@@ -1320,7 +1574,7 @@ bytes used=3000632, alloc=1244956, time=0.40
 > print_interp_cmpt__lc_of_data(posn_list_3d_size4,
 > 			      "result", "coeff_dy_", "data_",
 > 			      "3d.cube.size4/interp-dy.compute.c");
-bytes used=4000828, alloc=1310480, time=0.49
+bytes used=4000728, alloc=1310480, time=0.46
 > print_interp_cmpt__lc_of_data(posn_list_3d_size4,
 > 			      "result", "coeff_dz_", "data_",
 > 			      "3d.cube.size4/interp-dz.compute.c");
@@ -1342,6 +1596,7 @@ bytes used=4000828, alloc=1310480, time=0.49
 > print_interp_cmpt__lc_of_data(posn_list_3d_size4,
 > 			      "result", "coeff_dzz_", "data_",
 > 			      "3d.cube.size4/interp-dzz.compute.c");
+bytes used=5000896, alloc=1376004, time=0.60
 > 
 ########################################
 > 
@@ -1350,7 +1605,6 @@ bytes used=4000828, alloc=1310480, time=0.49
 #
 > 
 > data_var_list_3d_size5 := map(data_var_name, posn_list_3d_size5, "data_");
-bytes used=5001008, alloc=1310480, time=0.63
 data_var_list_3d_size5 := ["data_m2_m2_m2", "data_m1_m2_m2", "data_0_m2_m2",
 
     "data_p1_m2_m2", "data_p2_m2_m2", "data_m2_m1_m2", "data_m1_m1_m2",
@@ -1431,9 +1685,9 @@ data_var_list_3d_size5 := ["data_m2_m2_m2", "data_m1_m2_m2", "data_0_m2_m2",
 > 			     "3d.cube.size5/coeff-dz.store.c");
 > print_interp_coeff_var_store(posn_list_3d_size5, "factor", "coeff_dxx_",
 > 			     "3d.cube.size5/coeff-dxx.store.c");
+bytes used=6001052, alloc=1376004, time=0.74
 > print_interp_coeff_var_store(posn_list_3d_size5, "factor", "coeff_dxy_",
 > 			     "3d.cube.size5/coeff-dxy.store.c");
-bytes used=6001252, alloc=1310480, time=0.77
 > print_interp_coeff_var_store(posn_list_3d_size5, "factor", "coeff_dxz_",
 > 			     "3d.cube.size5/coeff-dxz.store.c");
 > print_interp_coeff_var_store(posn_list_3d_size5, "factor", "coeff_dyy_",
@@ -1447,9 +1701,9 @@ bytes used=6001252, alloc=1310480, time=0.77
 > 		    "3d.cube.size5/coeff-I.dcl.c");
 > print_name_list_dcl(map(coeff_name, posn_list_3d_size5, "coeff_dx_"), "fp",
 > 		    "3d.cube.size5/coeff-dx.dcl.c");
+bytes used=7001272, alloc=1376004, time=0.86
 > print_name_list_dcl(map(coeff_name, posn_list_3d_size5, "coeff_dy_"), "fp",
 > 		    "3d.cube.size5/coeff-dy.dcl.c");
-bytes used=7001440, alloc=1376004, time=0.89
 > print_name_list_dcl(map(coeff_name, posn_list_3d_size5, "coeff_dz_"), "fp",
 > 		    "3d.cube.size5/coeff-dz.dcl.c");
 > print_name_list_dcl(map(coeff_name, posn_list_3d_size5, "coeff_dxx_"), "fp",
@@ -1462,7 +1716,7 @@ bytes used=7001440, alloc=1376004, time=0.89
 > 		    "3d.cube.size5/coeff-dyy.dcl.c");
 > print_name_list_dcl(map(coeff_name, posn_list_3d_size5, "coeff_dyz_"), "fp",
 > 		    "3d.cube.size5/coeff-dyz.dcl.c");
-bytes used=8001776, alloc=1441528, time=0.98
+bytes used=8001444, alloc=1441528, time=0.98
 > print_name_list_dcl(map(coeff_name, posn_list_3d_size5, "coeff_dzz_"), "fp",
 > 		    "3d.cube.size5/coeff-dzz.dcl.c");
 > 
@@ -1475,10 +1729,10 @@ bytes used=8001776, alloc=1441528, time=0.98
 > print_interp_cmpt__lc_of_data(posn_list_3d_size5,
 > 			      "result", "coeff_dy_", "data_",
 > 			      "3d.cube.size5/interp-dy.compute.c");
+bytes used=9001600, alloc=1441528, time=1.08
 > print_interp_cmpt__lc_of_data(posn_list_3d_size5,
 > 			      "result", "coeff_dz_", "data_",
 > 			      "3d.cube.size5/interp-dz.compute.c");
-bytes used=9001948, alloc=1441528, time=1.09
 > print_interp_cmpt__lc_of_data(posn_list_3d_size5,
 > 			      "result", "coeff_dxx_", "data_",
 > 			      "3d.cube.size5/interp-dxx.compute.c");
@@ -1488,10 +1742,10 @@ bytes used=9001948, alloc=1441528, time=1.09
 > print_interp_cmpt__lc_of_data(posn_list_3d_size5,
 > 			      "result", "coeff_dxz_", "data_",
 > 			      "3d.cube.size5/interp-dxz.compute.c");
+bytes used=10001804, alloc=1507052, time=1.20
 > print_interp_cmpt__lc_of_data(posn_list_3d_size5,
 > 			      "result", "coeff_dyy_", "data_",
 > 			      "3d.cube.size5/interp-dyy.compute.c");
-bytes used=10002200, alloc=1441528, time=1.21
 > print_interp_cmpt__lc_of_data(posn_list_3d_size5,
 > 			      "result", "coeff_dyz_", "data_",
 > 			      "3d.cube.size5/interp-dyz.compute.c");
@@ -1501,4 +1755,4 @@ bytes used=10002200, alloc=1441528, time=1.21
 > 
 ################################################################################
 > quit
-bytes used=10722536, alloc=1441528, time=1.29
+bytes used=10794072, alloc=1507052, time=1.28
diff --git a/src/GeneralizedPolynomial-Uniform/common/cube_posns.maple b/src/GeneralizedPolynomial-Uniform/common/cube_posns.maple
index 0d23e1c..13c134c 100644
--- a/src/GeneralizedPolynomial-Uniform/common/cube_posns.maple
+++ b/src/GeneralizedPolynomial-Uniform/common/cube_posns.maple
@@ -26,6 +26,8 @@ posn_list_2d_size4 := map(ListTools[Reverse],
 			  hypercube_points([-1,-1], [+2,+2]));
 posn_list_2d_size5 := map(ListTools[Reverse],
 			  hypercube_points([-2,-2], [+2,+2]));
+posn_list_2d_size6 := map(ListTools[Reverse],
+			  hypercube_points([-2,-2], [+3,+3]));
 
 ################################################################################
 
@@ -40,3 +42,5 @@ posn_list_3d_size4 := map(ListTools[Reverse],
 			  hypercube_points([-1,-1,-1], [+2,+2,+2]));
 posn_list_3d_size5 := map(ListTools[Reverse],
 			  hypercube_points([-2,-2,-2], [+2,+2,+2]));
+posn_list_3d_size6 := map(ListTools[Reverse],
+			  hypercube_points([-2,-2,-2], [+3,+3,+3]));
diff --git a/src/GeneralizedPolynomial-Uniform/common/makefile b/src/GeneralizedPolynomial-Uniform/common/makefile
index c34120b..77f2ead 100644
--- a/src/GeneralizedPolynomial-Uniform/common/makefile
+++ b/src/GeneralizedPolynomial-Uniform/common/makefile
@@ -52,6 +52,7 @@ dirs	: 1d.dirs 2d.dirs 3d.dirs
 	mkdir 2d.cube.size3
 	mkdir 2d.cube.size4
 	mkdir 2d.cube.size5
+	mkdir 2d.cube.size6
 
 .PHONY	: 3d.dirs
 3d.dirs:
@@ -60,3 +61,4 @@ dirs	: 1d.dirs 2d.dirs 3d.dirs
 	mkdir 3d.cube.size3
 	mkdir 3d.cube.size4
 	mkdir 3d.cube.size5
+	mkdir 3d.cube.size6