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.cactuscode.org/arrangements/CactusNumerical/Cartoon2D/trunk@3 eec4d7dc-71c2-46d6-addf-10296150bf52

6 files changed, 1549 insertions, 0 deletions

diff --git a/src/Cartoon2DBC.c b/src/Cartoon2DBC.c
new file mode 100644
index 0000000..f1e1ebf
--- /dev/null
+++ b/src/Cartoon2DBC.c
@@ -0,0 +1,49 @@
+/*@@
+   @file      Cartoon2DBC.c
+   @date      Thu Nov  4 13:35:00 MET 1999
+   @author    Sai Iyer
+   @desc 
+              Apply Cartoon2D boundary conditions
+              Adapted from Bernd Bruegmann's cartoon_2d_tensorbound.c
+              and set_bound_axisym.c
+   @enddesc
+ @@*/
+
+static char *rcsid = "$Id$"
+
+#include "cctk.h"
+#include "cctk_arguments.h"
+#include "cctk_parameters.h"
+
+#include "cartoon2D.h"
+
+int Cartoon2DBCVarI(cGH *GH, int vi, int tensortype) { 
+  DECLARE_CCTK_PARAMETERS
+
+  static int firstcall = 1, pr = 0;
+
+  if (firstcall) {
+    firstcall = 0; 
+    pr = CCTK_Equals(verbose, "yes");
+  }
+  if (pr) printf("cartoon_2d_tensorbound called for %s\n", name);
+
+  return(SetBoundAxisym(GH, vi, tensortype));
+}
+
+void FMODIFIER FORTRAN_NAME(Cartoon2DBCVarI)(int *retval, cGH *GH, int *vi, int *tensortype) {
+  *retval = Cartoon2DBCVarI(GH, *vi, *tensortype);
+}
+
+int Cartoon2DBCVar(cGH *GH, const char *impvarname, int tensortype)
+{
+  int vi;
+  vi = CCTK_VarIndex(impvarname);
+  return(Cartoon2DBCVarI(GH, vi, tensortype));
+}
+
+void FMODIFIER FORTRAN_NAME(Cartoon2DBCVar)(int *retval, cGH *GH, ONE_FORTSTRING_ARG, int *tensortype) {
+  ONE_FORTSTRING_CREATE(impvarname)
+  *retval = Cartoon2DBCVar(GH, impvarname, *tensortype);
+  free(impvarname);
+}
diff --git a/src/CheckParameters.c b/src/CheckParameters.c
new file mode 100644
index 0000000..c7fcc79
--- /dev/null
+++ b/src/CheckParameters.c
@@ -0,0 +1,52 @@
+/*@@
+   @file      CheckParameters.c
+   @date      Wed Nov  3 10:17:46 MET 1999
+   @author    Sai Iyer
+   @desc 
+              Check Cartoon2D parameters
+   @enddesc 
+ @@*/
+
+#include "cctk.h" 
+#include "cctk_parameters.h"
+#include "cctk_arguments.h"
+
+static char *rcsid = "$Id$";
+
+ /*@@
+   @routine    Cartoon2D_CheckParameters
+   @date       Wed Nov  3 10:17:46 MET 1999
+   @author     Sai Iyer
+   @desc 
+               Check Cartoon2D parameters
+   @enddesc 
+   @calls      
+   @calledby   
+   @history 
+
+   @endhistory 
+
+@@*/
+
+void Cartoon2D_CheckParameters(CCTK_CARGUMENTS)
+{
+  DECLARE_CCTK_CARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+
+  if(cctk_gsh[1] != 2*cctk_nghostzones[1]+1)
+  {
+    CCTK_PARAMWARN("Grid size in y direction inappropriate.");
+  }
+
+  if(order < 2*cctk_nghostzones[0])
+  {
+    CCTK_PARAMWARN("Ghostzone width in x direction too small.");
+  }
+
+  if(order < 1 || order > 5)
+  {
+    CCTK_PARAMWARN("Interpolation order should be between 1 and 5.");
+  }
+  return;
+}
+
diff --git a/src/SetBoundAxisym.c b/src/SetBoundAxisym.c
new file mode 100644
index 0000000..37096e5
--- /dev/null
+++ b/src/SetBoundAxisym.c
@@ -0,0 +1,227 @@
+/* set_bound_axisym.c, BB */
+/*@@
+   @file      SetBoundAxisym.c
+   @date      Sat Nov  6 16:35:46 MET 1999
+   @author    Sai Iyer
+   @desc 
+              An implementation of Steve Brandt's idea for doing
+              axisymmetry with 3d stencils.
+              Cactus 4 version of Bernd Bruegmann's set_bound_axisym_GT.
+   @enddesc
+ @@*/
+
+static char *rcsid = "$Id$"
+
+#include "cctk.h"
+#include "cctk_arguments.h"
+#include "cctk_parameters.h"
+
+#include "cartoon2D.h"
+
+#include "math.h"
+
+/* Number of components */
+int tensor_type_n[N_TENSOR_TYPES] = {1, 3, 6};
+
+/* set boundaries of a grid tensor assuming axisymmetry 
+   - handles lower boundary in x
+   - does not touch other boundaries
+   - coordinates and rotation coefficients are independent of z and should
+     be precomputed
+   - this is also true for the constants in interpolator, but this may 
+     be too messy
+   - minimizes conceptual warpage, not computationally optimized
+*/
+/* uses rotation matrix and its inverse as linear transformation on
+   arbitrary tensor indices -- I consider this a good compromise
+   between doing index loops versus using explicit formulas in
+   cos(phi) etc with messy signs 
+*/
+int SetBoundAxiSym(cGH *GH, int vi, int tensortype) {
+
+{
+  int i, j, k, di, dj, dk, ijk;
+  int jj, n, s;
+  Double x, y, rho;
+  Double rxx, rxy, ryx, ryy;
+  Double sxx, sxy, syx, syy;
+  Double f[6], fx, fy, fz, fxx, fxy, fxz, fyy, fyz, fzz;
+  Double *t, *tx, *ty, *tz, *txx, *txy, *txz, *tyy, *tyz, *tzz;
+
+  if (0) printf("sba(%s)\n", GT->gf[0]->name);
+
+  /* stencil width: in x-direction use ghost zone width,
+                    in y-direction use ny/2, see below   */
+  s = GH->stencil_width;
+
+  /* Cactus loop in C:  grid points with y = 0 */
+  /* compare thorn_BAM_Elliptic */
+  /* strides used in stencils, should be provided by cactus */
+  di = 1;
+  dj = GH->lnx;
+  dk = GH->lnx*GH->lny;
+
+  /* y = 0 */
+  j = GH->ny/2;
+
+  /* z-direction: include lower and upper boundary */
+  for (k = 0; k < GH->lnz; k++)
+
+  /* y-direction: as many zombies as the evolution stencil needs */
+  for (jj = 1, dj = jj*GH->lnx; jj <= GH->ny/2; jj++, dj = jj*GH->lnx) 
+
+  /* x-direction: zombie for x < 0, including upper boundary for extrapol */
+  for (i = -s; i < GH->lnx; i++) {
+
+    /* index into linear memory */
+    ijk = DI(GH,i,j,k);
+   
+    /* what a neat way to hack in the zombie for x < 0, y = 0 */
+    if (i < 0) {i += s; ijk += s; dj = 0;}
+
+
+    /* compute coordinates (could also use Cactus grid function) */
+    x = GH->lx0 + GH->dx0 * i;
+    y = (dj) ? GH->dy0 * jj : 0;
+    rho = sqrt(x*x + y*y);
+
+    /* compute rotation matrix
+       note that this also works for x <= 0 (at lower boundary in x)
+       note that rho is nonzero by definition if cactus checks dy ... 
+    */
+    rxx = x/rho;
+    ryx = y/rho;
+    rxy = -ryx;
+    ryy = rxx;
+
+    /* inverse rotation matrix, assuming detr = 1 */
+    sxx = ryy;
+    syx = -ryx;
+    sxy = -rxy;
+    syy = rxx;
+
+    /* interpolate grid functions */
+    for (n = 0; n < GT->n; n++)
+      f[n] = axisym_interpolate(GH, GT->gf[n], i, di, ijk, rho); 
+
+    /* rotate grid tensor by matrix multiplication */
+    if (GT->tensor_type == TENSOR_TYPE_SCALAR) {
+      t = GT->gf[0]->data;
+      t[ijk+dj] = f[0];
+      t[ijk-dj] = f[0];
+    }
+    else if (GT->tensor_type == TENSOR_TYPE_U) {
+      tx = GT->gf[0]->data;
+      ty = GT->gf[1]->data;
+      tz = GT->gf[2]->data;
+      fx = f[0];
+      fy = f[1];
+      fz = f[2];
+
+      tx[ijk+dj] = rxx * fx + rxy * fy;
+      ty[ijk+dj] = ryx * fx + ryy * fy;
+      tz[ijk+dj] = fz;
+
+      tx[ijk-dj] = sxx * fx + sxy * fy;
+      ty[ijk-dj] = syx * fx + syy * fy;
+      tz[ijk-dj] = fz;
+    }
+    else if (GT->tensor_type == TENSOR_TYPE_DDSYM) {
+      txx = GT->gf[0]->data;
+      txy = GT->gf[1]->data;
+      txz = GT->gf[2]->data;
+      tyy = GT->gf[3]->data;
+      tyz = GT->gf[4]->data;
+      tzz = GT->gf[5]->data;
+      fxx = f[0];
+      fxy = f[1];
+      fxz = f[2];
+      fyy = f[3];
+      fyz = f[4];
+      fzz = f[5];
+      
+      txx[ijk+dj] = fxx*sxx*sxx + 2*fxy*sxx*syx + fyy*syx*syx;
+      tyy[ijk+dj] = fxx*sxy*sxy + 2*fxy*sxy*syy + fyy*syy*syy;
+      txy[ijk+dj] = fxx*sxx*sxy + fxy*(sxy*syx+sxx*syy) + fyy*syx*syy;
+      txz[ijk+dj] = fxz*sxx + fyz*syx;
+      tyz[ijk+dj] = fxz*sxy + fyz*syy;
+      tzz[ijk+dj] = fzz;
+
+      txx[ijk-dj] = fxx*rxx*rxx + 2*fxy*rxx*ryx + fyy*ryx*ryx;
+      tyy[ijk-dj] = fxx*rxy*rxy + 2*fxy*rxy*ryy + fyy*ryy*ryy;
+      txy[ijk-dj] = fxx*rxx*rxy + fxy*(rxy*ryx+rxx*ryy) + fyy*ryx*ryy;
+      txz[ijk-dj] = fxz*rxx + fyz*ryx;
+      tyz[ijk-dj] = fxz*rxy + fyz*ryy;
+      tzz[ijk-dj] = fzz;
+    } 
+    else {
+      /* oops */
+    }
+
+    /* what a neat way to hack out the zombies for x < 0, y = 0 */
+    if (dj == 0) {i -= s; ijk -= s; dj = jj*GH->lnx;}
+  }
+
+  /* syncs needed after interpolation (actually only for x direction) */
+  for (n = 0; n < GT->n; n++)
+    SyncGF(GH, GT->gf[n]); 
+}
+
+
+
+
+
+/* interpolation on x-axis */
+Double axisym_interpolate(pGH *GH, pGF *GF, int i, int di, int ijk, Double x)
+{
+  Double *f = GF->data;
+  Double c, d, x0;
+
+  /* find i such that x(i) < x <= x(i+1)
+     for rotation on entry always x > x(i), but sometimes also x > x(i+1) */
+  while (x > GH->lx0 + GH->dx0 * (i+1)) {i++; ijk++;}
+
+  /* the original demo */
+  if (0) {
+    c = (f[ijk+di] - f[ijk])/GH->dx0;
+    d = f[ijk];
+    x0 =  GH->lx0 + GH->dx0 * i; 
+    return c*(x-x0) + d;
+  }
+
+  /* first attempt to interface to JT's interpolator */
+  if (1) {
+    double interpolate_local(int order, double x0, double dx, 
+			     double y[], double x);
+    double y[6], r;
+    int n, offset;
+    int order = cartoon_order;
+
+    /* offset of stencil, note that rotation leads to x close to x0 
+       for large x */
+    offset = order/2;
+
+    /* shift stencil at boundaries */
+    /* note that for simplicity we don't distinguish between true boundaries
+       and decomposition boundaries: the sync fixes that! */
+    if (i-offset < 0) offset = i;
+    if (i-offset+order >= GH->lnx) offset = i+order-GH->lnx+1;
+
+    /* fill in info */
+    /* fills in old data near axis, but as long as stencil at axis is 
+       centered, no interpolation happens anyway */
+    x0 = GH->lx0 + GH->dx0 * (i-offset);
+    for (n = 0; n <= order; n++) {
+      y[n] = f[ijk-offset+n];
+      if (0) printf("y[%d]=%6.4f ", n, y[n]); 
+    }
+
+    /* call interpolator */
+    r = interpolate_local(order, x0, GH->dx0, y, x);
+
+    if (0) printf("i %d  x %7.4f  x0 %7.4f  r %7.4f\n", i, x, x0, r); 
+    return r;
+  }
+
+  return 0;
+}
diff --git a/src/error_exit.c b/src/error_exit.c
new file mode 100644
index 0000000..748de99
--- /dev/null
+++ b/src/error_exit.c
@@ -0,0 +1,63 @@
+/* error_exit.c -- print an error message and cleanly shutdown cactus */
+/* $Id$ */
+
+#include <stdio.h>
+#include <stdarg.h>
+#include <stdlib.h>
+
+#include "cactus.h"
+#include "thorn_cartoon_2d/src/error_exit.h"
+
+/*
+ * Jonathan Thornburg's favorite syntactic sugar to make the two branches
+ * of an if-else statement symmetrical:
+ *
+ * Also note (in)famous joke saying
+ *	"syntactic sugar causes cancer of the semicolon"
+ * :) :) :)
+ */
+#define then	/* empty */
+
+/*****************************************************************************/
+
+/*
+ * This function prints an error message (formatted via vfprintf(3S)),
+ * then cleanly shuts down cactus.  It does *not* return to its caller.
+ *
+ * Despite not actually returning, this function is declared as returning
+ * an  int , so it may be easily used in conditional expressions like
+ *      foo = bar_ok ? baz(bar) : error_exit(...);
+ *
+ * Usage:
+ *	error_exit(exit_code, message, args...)
+ * where  exit_code  is one of the codes defined in "error_exit.h"
+ *
+ * ***** THIS VERSION IS FOR ANSI/ISO C *****
+ *
+ * Arguments:
+ * exit_code = (in) One of the codes defined in "error_exit.h", or any
+ *		    other integer to pass back to the shell as our exit
+ *		    status.
+ * format = (in) vprintf(3S) format string for error message to be printed.
+ * args... = (in) Any additional arguments are (presumably) used in formatting
+ *		  the error message string.
+ */
+#ifdef __cplusplus
+extern "C"
+#endif
+/*VARARGS*/
+/* PUBLIC */
+int error_exit(const int exit_code, const char *const format, ...)
+{
+va_list ap;
+
+va_start(ap, format);
+vfprintf(stderr, format, ap);
+va_end(ap);
+
+ if	(exit_code == 0) {
+   then STOP;							/*NOTREACHED*/
+ } else if (exit_code > 0) {
+   then ERROR_STOP(exit_code);					/*NOTREACHED*/
+ } else	ERROR_ABORT(exit_code);					/*NOTREACHED*/
+}
diff --git a/src/interpolate.c b/src/interpolate.c
new file mode 100644
index 0000000..346a1fa
--- /dev/null
+++ b/src/interpolate.c
@@ -0,0 +1,1149 @@
+/* interpolate.c -- interpolation functions */
+/*
+** <<<design notes>>>
+
+ * interpolate_global - Lagrange polynomial interpolation (global)
+ * interpolate_local - ... (local)
+** interpolate_local_order1 - ... (order 1) (linear) (2-point)
+** interpolate_local_order2 - ... (order 2) (3-point)
+** interpolate_local_order3 - ... (order 3) (4-point)
+** interpolate_local_order4 - ... (order 4) (5-point)
+** interpolate_local_order5 - ... (order 5) (6-point)
+
+ * interpolate_deriv_global - Lagrange polynomial differentiation (global)
+ * interpolate_deriv_local - ... (local)
+** interpolate_deriv_local_order1 - ... (order 1) (linear) (2-point)
+** interpolate_deriv_local_order2 - ... (order 2) (3-point)
+** interpolate_deriv_local_order3 - ... (order 3) (4-point)
+** interpolate_deriv_local_order4 - ... (order 4) (5-point)
+** interpolate_deriv_local_order5 - ... (order 5) (6-point)
+
+ * interpolate_nu_global - Lagrange polynomial interp (nonuniform/global)
+ * interpolate_nu_local - ... (nonuniform/local)
+** interpolate_nu_local_order3 - ... (order 3) (4-point)
+** interpolate_nu_local_order4 - ... (order 4) (5-point)
+** interpolate_nu_local_order5 - ... (order 5) (6-point)
+
+ * subarray_posn - compute subarray positions for interpolation etc
+ */
+
+
+
+/* minimal thing to do so that this files compiles in cactus, BB */ 
+#define CACTUS_STANDALONE
+
+#ifdef CACTUS_STANDALONE
+
+#include "cactus.h"
+#include "thorn_cartoon_2d/src/error_exit.h"
+#define PRIVATE 
+#define PUBLIC 
+#define then
+#define boolean int
+#define FALSE 0
+#define TRUE  1
+#define HOW_MANY_IN_RANGE(a,b) ((b) - (a) + 1)
+#define floor_double_to_int(x) ((int)floor((double)x))
+
+double interpolate_local(int order, double x0, double dx, 
+			 double y[], double x);
+double interpolate_deriv_local(int order, double x0, double dx,
+			       double y[], double x);
+double interpolate_nu_local(int order, double xx[], double yy[], double x);
+
+#else
+
+#include <stdio.h>
+#include <jt/stdc.h>
+#include <jt/util.h>
+#include <jt/interpolate.h>
+
+#endif
+
+/* prototypes for private functions defined in this file */
+PRIVATE double interpolate_local_order1(double x0, double dx,
+					double y[],
+					double x);
+PRIVATE double interpolate_local_order2(double x0, double dx,
+					double y[],
+					double x);
+PRIVATE double interpolate_local_order3(double x0, double dx,
+					double y[],
+					double x);
+PRIVATE double interpolate_local_order4(double x0, double dx,
+					double y[],
+					double x);
+PRIVATE double interpolate_local_order5(double x0, double dx,
+					double y[],
+					double x);
+PRIVATE double interpolate_deriv_local_order1(double x0, double dx,
+					      double y[],
+					      double x);
+PRIVATE double interpolate_deriv_local_order2(double x0, double dx,
+					      double y[],
+					      double x);
+PRIVATE double interpolate_deriv_local_order3(double x0, double dx,
+					      double y[],
+					      double x);
+PRIVATE double interpolate_deriv_local_order4(double x0, double dx,
+					      double y[],
+					      double x);
+PRIVATE double interpolate_deriv_local_order5(double x0, double dx,
+					      double y[],
+					      double x);
+PRIVATE double interpolate_nu_local_order3(double xx[],
+					   double yy[],
+					   double x);
+PRIVATE double interpolate_nu_local_order4(double xx[],
+					   double yy[],
+					   double x);
+PRIVATE double interpolate_nu_local_order5(double xx[],
+					   double yy[],
+					   double x);
+
+/*****************************************************************************/
+/*****************************************************************************/
+/*****************************************************************************/
+
+/*
+ * ***** Design Notes *****
+ */
+
+/*
+ * The basic problem of interpolation is this: given a set of data points
+ *  (x[i],y[i])  sampled from a smooth function  y = y(x) , we wish to
+ * estimate  y  or one if its derivatives, at arbitrary  x  values.
+ *
+ * We classify interpolation functions along 3 dimensions:
+ * - local (the underlying routines, where all data points are used),
+ *   vs global (search for the right place in a large array, then do a
+ *   local interpolation),
+ * - uniform (the  x[i]  values are uniformly spaced, and may be implicitly
+ *   specified as a linear function x[i] = x0 + i*dx ), vs non-uniform
+ *   (the  x[i]  values are (potentially) non-uniformly spaced, and hence
+ *   must be supplied to the interpolation function as an explicit array),
+ * - function (we want to estimate  y ) vs derivative (we want to estimate
+ *    dy/dx ; higher derivatives are also possible, but typically don't
+ *   occur in an interpolation context).
+ */
+
+/*****************************************************************************/
+/*****************************************************************************/
+/*****************************************************************************/
+
+/*
+ * This function does global Lagrange polynomial interpolation on a
+ * uniform grid.  That is, given  N_points  uniformly spaced data points
+ *  (x0 + i*dx, y[i]) , i = 0...N_points-1, and given a specified  x
+ * coordinate, it locates a subarray of  order+1  data points centered
+ * around  x , interpolates a degree-order (Lagrange) polynomial through
+ * them, and evaluates this at  x .
+ *
+ * Except for possible end effects (see  end_action  (below)), the local
+ * interpolation is internally centered to improve its conditioning.  Even
+ * so, the interpolation is highly ill-conditioned for small  dx  and/or
+ * large  order  and/or  x  outside the domain of the data points.
+ *
+ * The interpolation formulas were (are) all derived via Maple, see
+ * "./interpolate.in" and "./interpolate.out".
+ *
+ * Arguments:
+ * N_points = (in) The number of data points.
+ * order = (in) The order of the local interpolation, i.e. the degree
+ *		of the interpolated polynomial.
+ * end_action = (in) Specifies what action to take if the (centered)
+ *		     local interpolant would need data from outside
+ *		     the domain of the data points.
+ *		     'e' ==> Do an  error_exit() .
+ *		     'o' ==> Off-center the local interpolant as
+ *			     necessary to keep it inside the domain
+ *			     of the data points.
+ * x0 = (in) The x coordinate corresponding to  y[0] .
+ * dx = (in) The x spacing between the data points.
+ * y[0...N_points) = (in) The y data array.
+ * x = (in) The x coordinate for the interpolation.
+ */
+PUBLIC double interpolate_global(int N_points, int order,
+				 char end_action,
+				 double x0, double dx,
+				 double y[],
+				 double x)
+{
+int N_interp, i, i_lo, i_hi, shift;
+
+/*
+ * locate subarray for local interpolation
+ */
+N_interp = order + 1;
+i = floor_double_to_int( (x - x0)/dx );
+subarray_posn(N_interp, i, &i_lo, &i_hi);
+if (! ((i_lo >= 0) && (i_hi < N_points)) )
+   then {
+	/* centered local interpolation would need data */
+	/* from outside  y[]  array */
+	switch	(end_action)
+		{
+	case 'e':
+		error_exit(ERROR_EXIT,
+"***** interpolate_global: local interpolation needs data\n"
+"                          from outside data array!\n"
+"                          N_points=%d   order=%d   N_interp=%d\n"
+"                          x0=%g   dx=%g\n"
+"                          x=%g\n"
+"                          ==> i=%d   i_[lo,hi]=[%d,%d]\n"
+,
+			   N_points, order, N_interp,
+			   x0, dx,
+			   x,
+			   i, i_lo, i_hi);			/*NOTREACHED*/
+
+	case 'o':
+		/* off-center the local interpolation as necessary */
+		/* to keep it within the  y[]  array */
+		shift = 0;
+		if (i_lo < 0) then shift = 0 - i_lo;
+		if (i_hi >= N_points) then shift = (N_points-1) - i_hi;
+
+		i_lo += shift;
+		i_hi += shift;
+
+		if (! ((i_lo >= 0) && (i_hi < N_points)) )
+		   then error_exit(ERROR_EXIT,
+"***** interpolate_global: shifted local interpolation needs data\n"
+"                          from outside data array on other side!\n"
+"                          (this should only happen if N_interp > N_points!)\n"
+"                          N_points=%d   order=%d   N_interp=%d\n"
+"                          x0=%g   dx=%g\n"
+"                          x=%g\n"
+"                          ==> i=%d\n"
+"                          ==> shift=%d   i_[lo,hi]=[%d,%d]\n"
+,
+				   N_points, order, N_interp,
+				   x0, dx,
+				   x,
+				   i,
+				   shift, i_lo, i_hi);		/*NOTREACHED*/
+		break;
+
+	default:
+		error_exit(PANIC_EXIT,
+"***** interpolate_global: unknown end_action=0x%02x='%c'\n",
+			   (int) end_action, end_action);	/*NOTREACHED*/
+		}
+	}
+
+/*
+ * do local interpolation
+ */
+return(
+   interpolate_local(order,
+		     x0 + i_lo*dx, dx,
+		     y + i_lo,
+		     x)
+      );
+}
+
+/*****************************************************************************/
+
+/*
+ * This function does local Lagrange polynomial interpolation on a
+ * uniform grid.  That is, given  order+1  uniformly spaced data points
+ *  (x0 + i*dx, y[i]) , i = 0...order, it interpolates a degree-(order)
+ * (Lagrange) polynomial through them and evaluates this polynomial at
+ * a specified  x  coordinate.  In general the error in this computed
+ * function value is  O(h^(order+1)) .
+ *
+ * Except for possible end effects (see  end_action  (below)), the local
+ * interpolation is internally centered to improve its conditioning.  Even
+ * so, the interpolation is highly ill-conditioned for small  dx  and/or
+ * large  order  and/or  x  outside the domain of the data points.
+ *
+ * The interpolation formulas were (are) all derived via Maple, see
+ * "./interpolate.in" and "./interpolate.out".
+ *
+ * Arguments:
+ * order = (in) The order of the local interpolation, i.e. the degree
+ *		of the interpolated polynomial.
+ * x0 = (in) The x coordinate corresponding to y[0].
+ * dx = (in) The x spacing between the data points.
+ * y[0...order] = (in) The y data array.
+ * x = (in) The x coordinate for the interpolation.
+ */
+PUBLIC double interpolate_local(int order,
+				double x0, double dx,
+				double y[],
+				double x)
+{
+switch	(order)
+	{
+case 1:
+	return(  interpolate_local_order1(x0, dx, y, x)  );
+	break;
+
+case 2:
+	return(  interpolate_local_order2(x0, dx, y, x)  );
+	break;
+
+case 3:
+	return(  interpolate_local_order3(x0, dx, y, x)  );
+	break;
+
+case 4:
+	return(  interpolate_local_order4(x0, dx, y, x)  );
+	break;
+
+case 5:
+	return(  interpolate_local_order5(x0, dx, y, x)  );
+	break;
+
+default:
+	error_exit(ERROR_EXIT,
+"***** interpolate_local: order=%d not implemented yet!\n",
+		   order);					/*NOTREACHED*/
+	}
+}
+
+/*****************************************************************************/
+
+/*
+ * This function interpolates a 1st-order (linear) Lagrange polynomial
+ * in a 2-point [0...1] data array centered about the [0.5] grid position.
+ */
+PRIVATE double interpolate_local_order1(double x0, double dx,
+					double y[],
+					double x)
+{
+double c0 = (+1.0/2.0)*y[0] + (+1.0/2.0)*y[1];
+double c1 =          - y[0] +          + y[1];
+
+double xc = x0 + 0.5*dx;
+double xr = (x - xc) / dx;
+
+return(  c0 + xr*c1  );
+}
+
+/*****************************************************************************/
+
+/*
+ * This function interpolates a 2nd-order (quadratic) Lagrange polynomial
+ * in a 3-point [0...2] data array centered about the [1] grid position.
+ */
+PRIVATE double interpolate_local_order2(double x0, double dx,
+					double y[],
+					double x)
+{
+double c0 = y[1];
+double c1 = (-1.0/2.0)*y[0]        + (+1.0/2.0)*y[2];
+double c2 = (+1.0/2.0)*y[0] - y[1] + (+1.0/2.0)*y[2];
+
+double xc = x0 + 1.0*dx;
+double xr = (x - xc) / dx;
+
+return(  c0 + xr*(c1 + xr*c2)  );
+}
+
+/*****************************************************************************/
+
+/*
+ * This function interpolates a 3rd-order (cubic) Lagrange polynomial
+ * in a 4-point [0...3] data array centered about the [1.5] grid position.
+ */
+PRIVATE double interpolate_local_order3(double x0, double dx,
+					double y[],
+					double x)
+{
+double c0 =   (-1.0/16.0)*y[0] + (+9.0/16.0)*y[1]
+	    + (+9.0/16.0)*y[2] + (-1.0/16.0)*y[3];
+double c1 =   (+1.0/24.0)*y[0] + (-9.0/ 8.0)*y[1]
+	    + (+9.0/ 8.0)*y[2] + (-1.0/24.0)*y[3];
+double c2 =   (+1.0/ 4.0)*y[0] + (-1.0/ 4.0)*y[1]
+	    + (-1.0/ 4.0)*y[2] + (+1.0/ 4.0)*y[3];
+double c3 =   (-1.0/ 6.0)*y[0] + (+1.0/ 2.0)*y[1]
+	    + (-1.0/ 2.0)*y[2] + ( 1.0/ 6.0)*y[3];
+
+double xc = x0 + 1.5*dx;
+double xr = (x - xc) / dx;
+
+return(  c0 + xr*(c1 + xr*(c2 + xr*c3))  );
+}
+
+/*****************************************************************************/
+
+/*
+ * This function interpolates a 4th-order (quartic) Lagrange polynomial
+ * in a 5-point [0...4] data array centered about the [2] grid position.
+ */
+PRIVATE double interpolate_local_order4(double x0, double dx,
+					double y[],
+					double x)
+{
+double c0 = y[2];
+double c1 =   (+1.0/12.0)*y[0] + (-2.0/ 3.0)*y[1]
+	    + (+2.0/ 3.0)*y[3] + (-1.0/12.0)*y[4];
+double c2 = + (-1.0/24.0)*y[0] + (+2.0/ 3.0)*y[1] + (-5.0/ 4.0)*y[2]
+	    + (+2.0/ 3.0)*y[3] + (-1.0/24.0)*y[4];
+double c3 = + (-1.0/12.0)*y[0] + (+1.0/ 6.0)*y[1]
+	    + (-1.0/ 6.0)*y[3] + (+1.0/12.0)*y[4];
+double c4 = + (+1.0/24.0)*y[0] + (-1.0/ 6.0)*y[1] + (+1.0/ 4.0)*y[2]
+	    + (-1.0/ 6.0)*y[3] + (+1.0/24.0)*y[4];
+
+double xc = x0 + 2.0*dx;
+double xr = (x - xc) / dx;
+
+return(  c0 + xr*(c1 + xr*(c2 + xr*(c3 + xr*c4)))  );
+}
+
+/*****************************************************************************/
+
+/*
+ * This function interpolates a 5th-order (quintic) Lagrange polynomial
+ * in a 6-point [0...5] data array centered about the [2.5] grid position.
+ */
+PRIVATE double interpolate_local_order5(double x0, double dx,
+					double y[],
+					double x)
+{
+double c0 =   (+ 3.0/256.0)*y[0] + (-25.0/256.0)*y[1]
+	    + (+75.0/128.0)*y[2] + (+75.0/128.0)*y[3]
+	    + (-25.0/256.0)*y[4] + (+ 3.0/256.0)*y[5];
+double c1 =   (- 3.0/640.0)*y[0] + (+25.0/384.0)*y[1]
+	    + (-75.0/ 64.0)*y[2] + (+75.0/ 64.0)*y[3]
+	    + (-25.0/384.0)*y[4] + (+ 3.0/640.0)*y[5];
+double c2 =   (- 5.0/ 96.0)*y[0] + (+13.0/ 32.0)*y[1]
+	    + (-17.0/ 48.0)*y[2] + (-17.0/ 48.0)*y[3]
+	    + (+13.0/ 32.0)*y[4] + (- 5.0/ 96.0)*y[5];
+double c3 =   (+ 1.0/ 48.0)*y[0] + (-13.0/ 48.0)*y[1]
+	    + (+17.0/ 24.0)*y[2] + (-17.0/ 24.0)*y[3]
+	    + (+13.0/ 48.0)*y[4] + (- 1.0/ 48.0)*y[5];
+double c4 =   (+ 1.0/ 48.0)*y[0] + (- 1.0/ 16.0)*y[1]
+	    + (+ 1.0/ 24.0)*y[2] + (+ 1.0/ 24.0)*y[3]
+	    + (- 1.0/ 16.0)*y[4] + (+ 1.0/ 48.0)*y[5];
+double c5 =   (- 1.0/120.0)*y[0] + (+ 1.0/ 24.0)*y[1]
+	    + (- 1.0/ 12.0)*y[2] + (+ 1.0/ 12.0)*y[3]
+	    + (- 1.0/ 24.0)*y[4] + (+ 1.0/120.0)*y[5];
+
+double xc = x0 + 2.5*dx;
+double xr = (x - xc) / dx;
+
+return(  c0 + xr*(c1 + xr*(c2 + xr*(c3 + xr*(c4 + xr*c5))))  );
+}
+
+/*****************************************************************************/
+/*****************************************************************************/
+/*****************************************************************************/
+
+/*
+ * This function does global Lagrange polynomial differentiation on a
+ * uniform grid.  That is, given  N_points  uniformly spaced data points
+ *  (x0 + i*dx, y[i]) , i = 0...N_points-1, and given a specified  x
+ * coordinate, it locates a subarray of  order+1  data points centered
+ * around  x , interpolates a degree-order (Lagrange) polynomial through
+ * them, differentiates this, and evaluates the derivative at  x .
+ *
+ * Except for possible end effects (see  end_action  (below)), the local
+ * interpolation is internally centered to improve its conditioning.  Even
+ * so, the interpolation is highly ill-conditioned for small  dx  and/or
+ * large  order  and/or  x  outside the domain of the data points.
+ *
+ * The interpolation formulas were (are) all derived via Maple, see
+ * "./interpolate.in" and "./interpolate.out".
+ *
+ * Arguments:
+ * N_points = (in) The number of data points.
+ * order = (in) The order of the local interpolation, i.e. the degree
+ *		of the interpolated polynomial.
+ * end_action = (in) Specifies what action to take if the (centered)
+ *		     local interpolant would need data from outside
+ *		     the domain of the data points.
+ *		     'e' ==> Do an  error_exit() .
+ *		     'o' ==> Off-center the local interpolant as
+ *			     necessary to keep it inside the domain
+ *			     of the data points.
+ * x0 = (in) The x coordinate corresponding to y[0].
+ * dx = (in) The x spacing between the data points.
+ * y[0...N_points) = (in) The y data array.
+ * x = (in) The x coordinate for the interpolation.
+ */
+PUBLIC double interpolate_deriv_global(int N_points, int order,
+				       char end_action,
+				       double x0, double dx,
+				       double y[],
+				       double x)
+{
+int N_interp, i, i_lo, i_hi, shift;
+
+/*
+ * locate subarray for local differentiation
+ */
+N_interp = order + 1;
+i = floor_double_to_int( (x - x0)/dx );
+subarray_posn(N_interp, i, &i_lo, &i_hi);
+if (! ((i_lo >= 0) && (i_hi < N_points)) )
+   then {
+	/* centered local differentiation would need data */
+	/* from outside  y[]  array */
+	switch	(end_action)
+		{
+	case 'e':
+		error_exit(ERROR_EXIT,
+"***** interpolate_deriv_global: local differentiation needs data\n"
+"                                from outside data array!\n"
+"                                N_points=%d   order=%d   N_interp=%d\n"
+"                                x0=%g   dx=%g\n"
+"                                x=%g\n"
+"                                ==> i=%d   i_[lo,hi]=[%d,%d]\n"
+,
+			   N_points, order, N_interp,
+			   x0, dx,
+			   x,
+			   i, i_lo, i_hi);			/*NOTREACHED*/
+
+	case 'o':
+		/* off-center the local differentiation as necessary */
+		/* to keep it within the  y[]  array */
+		shift = 0;
+		if (i_lo < 0) then shift = 0 - i_lo;
+		if (i_hi >= N_points) then shift = (N_points-1) - i_hi;
+
+		i_lo += shift;
+		i_hi += shift;
+
+		if (! ((i_lo >= 0) && (i_hi < N_points)) )
+		   then error_exit(ERROR_EXIT,
+"***** interpolate_deriv_global: shifted local differentiation needs data\n"
+"                                from outside data array on other side!\n"
+"                                (this should only happen\n"
+"                                 if N_interp > N_points!)\n"
+"                                N_points=%d   order=%d   N_interp=%d\n"
+"                                x0=%g   dx=%g\n"
+"                                x=%g\n"
+"                                ==> i=%d\n"
+"                                ==> shift=%d   i_[lo,hi]=[%d,%d]\n"
+,
+				   N_points, order, N_interp,
+				   x0, dx,
+				   x,
+				   i,
+				   shift, i_lo, i_hi);		/*NOTREACHED*/
+		break;
+
+	default:
+		error_exit(PANIC_EXIT,
+"***** interpolate_deriv_global: unknown end_action=0x%02x='%c'\n",
+			   (int) end_action, end_action);	/*NOTREACHED*/
+		}
+	}
+
+/*
+ * do local differentiation
+ */
+return(
+   interpolate_deriv_local(order,
+			   x0 + i_lo*dx, dx,
+			   y + i_lo,
+			   x)
+      );
+}
+
+/*****************************************************************************/
+
+/*
+ * This function does local Lagrange polynomial differentiation on a
+ * uniform grid.  That is, given  order+1  uniformly spaced data points
+ *  (x0 + i*dx, y[i]) , i = 0...order, it interpolates a degree-(order)
+ * (Lagrange) polynomial through them, differentiates it, and evaluates
+ * the derivative at a specified  x  coordinate.  In general the error
+ * in this computed function value is  O(h^(order)) .
+ *
+ * Except for possible end effects (see  end_action  (below)), the local
+ * interpolation is internally centered to improve its conditioning.  Even
+ * so, the interpolation is highly ill-conditioned for small  dx  and/or
+ * large  order  and/or  x  outside the domain of the data points.
+ *
+ * The interpolation formulas were (are) all derived via Maple.
+ *
+ * Arguments:
+ * order = (in) The order of the local interpolation, i.e. the degree
+ *		of the interpolated polynomial.
+ * x0 = (in) The x coordinate corresponding to y[0].
+ * dx = (in) The x spacing between the data points.
+ * y[0...order] = (in) The y data array.
+ * x = (in) The x coordinate for the interpolation.
+ */
+PUBLIC double interpolate_deriv_local(int order,
+				      double x0, double dx,
+				      double y[],
+				      double x)
+{
+switch	(order)
+	{
+case 1:
+	return(  interpolate_deriv_local_order1(x0, dx, y, x)  );
+	break;
+
+case 2:
+	return(  interpolate_deriv_local_order2(x0, dx, y, x)  );
+	break;
+
+case 3:
+	return(  interpolate_deriv_local_order3(x0, dx, y, x)  );
+	break;
+
+case 4:
+	return(  interpolate_deriv_local_order4(x0, dx, y, x)  );
+	break;
+
+case 5:
+	return(  interpolate_deriv_local_order5(x0, dx, y, x)  );
+	break;
+
+default:
+	error_exit(ERROR_EXIT,
+"***** interpolate_deriv_local: order=%d not implemented yet!\n",
+		   order);					/*NOTREACHED*/
+	}
+}
+
+/*****************************************************************************/
+
+/*
+ * This function interpolates a 1st-order (linear) Lagrange polynomial
+ * in a 2-point [0...1] data array centered about the [0.5] grid position,
+ * then evaluates the (0th-order) (constant) derivative of that polynomial.
+ */
+PRIVATE double interpolate_deriv_local_order1(double x0, double dx,
+					      double y[],
+					      double x)
+{
+double c0 = y[1] - y[0];
+
+return(  c0 / dx  );
+}
+
+/*****************************************************************************/
+
+/*
+ * This function interpolates a 2nd-order (quadratic) Lagrange polynomial
+ * in a 3-point [0...2] data array centered about the [1] grid position,
+ * then evaluates the (1st-order) (linear) derivative of that polynomial.
+ */
+PRIVATE double interpolate_deriv_local_order2(double x0, double dx,
+					      double y[],
+					      double x)
+{
+double c0 = (-1.0/2.0)*y[0]            + (+1.0/2.0)*y[2];
+double c1 =            y[0] - 2.0*y[1] +            y[2];
+
+double xc = x0 + 1.0*dx;
+double xr = (x - xc) / dx;
+
+return(  (c0 + xr*c1) / dx  );
+}
+
+/*****************************************************************************/
+
+/*
+ * This function interpolates a 3rd-order (cubic) Lagrange polynomial
+ * in a 4-point [0...3] data array centered about the [1.5] grid position,
+ * then evaluates the (2nd-order) (quadratic) derivative of that polynomial.
+ */
+PRIVATE double interpolate_deriv_local_order3(double x0, double dx,
+					      double y[],
+					      double x)
+{
+double c0 =   (+1.0/24.0)*y[0] + (-9.0/ 8.0)*y[1]
+	    + (+9.0/ 8.0)*y[2] + (-1.0/24.0)*y[3];
+double c1 =   (+1.0/ 2.0)*y[0] + (-1.0/ 2.0)*y[1]
+	    + (-1.0/ 2.0)*y[2] + (+1.0/ 2.0)*y[3];
+double c2 =   (-1.0/ 2.0)*y[0] + (+3.0/ 2.0)*y[1]
+	    + (-3.0/ 2.0)*y[2] + ( 1.0/ 2.0)*y[3];
+
+double xc = x0 + 1.5*dx;
+double xr = (x - xc) / dx;
+
+return(  (c0 + xr*(c1 + xr*c2)) / dx  );
+}
+
+/*****************************************************************************/
+
+/*
+ * This function interpolates a 4th-order (quartic) Lagrange polynomial
+ * in a 5-point [0...4] data array centered about the [2] grid position,
+ * then evaluates the (3rd-order) (cubic) derivative of that polynomial.
+ */
+PRIVATE double interpolate_deriv_local_order4(double x0, double dx,
+					      double y[],
+					      double x)
+{
+double c0 = + (+1.0/12.0)*y[0] + (-2.0/ 3.0)*y[1]
+	    + (+2.0/ 3.0)*y[3] + (-1.0/12.0)*y[4];
+double c1 = + (-1.0/12.0)*y[0] + (+4.0/ 3.0)*y[1] + (-5.0/ 2.0)*y[2]
+	    + (+4.0/ 3.0)*y[3] + (-1.0/12.0)*y[4];
+double c2 = + (-1.0/ 4.0)*y[0] + (+1.0/ 2.0)*y[1]
+	    + (-1.0/ 2.0)*y[3] + (+1.0/ 4.0)*y[4];
+double c3 = + (+1.0/ 6.0)*y[0] + (-2.0/ 3.0)*y[1] +             y[2]
+	    + (-2.0/ 3.0)*y[3] + (+1.0/ 6.0)*y[4];
+
+double xc = x0 + 2.0*dx;
+double xr = (x - xc) / dx;
+
+return(  (c0 + xr*(c1 + xr*(c2 + xr*c3))) / dx  );
+}
+
+/*****************************************************************************/
+
+/*
+ * This function interpolates a 5th-order (quartic) Lagrange polynomial
+ * in a 6-point [0...5] data array centered about the [2.5] grid position,
+ * then evaluates the (4th-order) (quartic) derivative of that polynomial.
+ */
+PRIVATE double interpolate_deriv_local_order5(double x0, double dx,
+					      double y[],
+					      double x)
+{
+double c0 =   (- 3.0/640.0)*y[0] + (+25.0/384.0)*y[1]
+	    + (-75.0/ 64.0)*y[2] + (+75.0/ 64.0)*y[3]
+	    + (-25.0/384.0)*y[4] + (+ 3.0/640.0)*y[5];
+double c1 =   (- 5.0/ 48.0)*y[0] + (+13.0/ 16.0)*y[1]
+	    + (-17.0/ 24.0)*y[2] + (-17.0/ 24.0)*y[3]
+	    + (+13.0/ 16.0)*y[4] + (- 5.0/ 48.0)*y[5];
+double c2 =   (+ 1.0/ 16.0)*y[0] + (-13.0/ 16.0)*y[1]
+	    + (+17.0/  8.0)*y[2] + (-17.0/  8.0)*y[3]
+	    + (+13.0/ 16.0)*y[4] + (- 1.0/ 16.0)*y[5];
+double c3 =   (+ 1.0/ 12.0)*y[0] + (- 1.0/  4.0)*y[1]
+	    + (+ 1.0/  6.0)*y[2] + (+ 1.0/  6.0)*y[3]
+	    + (- 1.0/  4.0)*y[4] + (+ 1.0/ 12.0)*y[5];
+double c4 =   (- 1.0/ 24.0)*y[0] + (+ 5.0/ 24.0)*y[1]
+	    + (- 5.0/ 12.0)*y[2] + (+ 5.0/ 12.0)*y[3]
+	    + (- 5.0/ 24.0)*y[4] + (+ 1.0/ 24.0)*y[5];
+
+double xc = x0 + 2.5*dx;
+double xr = (x - xc) / dx;
+
+return(  (c0 + xr*(c1 + xr*(c2 + xr*(c3 + xr*c4)))) / dx  );
+}
+
+/*****************************************************************************/
+/*****************************************************************************/
+/*****************************************************************************/
+
+/*
+ * This function does global Lagrange polynomial interpolation on a
+ * nonuniform grid.  That is, given  N_points  uniformly spaced data
+ * points  (x0 + i*dx, y[i]) , i = 0...N_points-1, and given a specified
+ *  x  coordinate, it locates a subarray of  order+1  data points centered
+ * around  x , interpolates a degree-order (Lagrange) polynomial through
+ * them, and evaluates this at  x .
+ *
+ * The local interpolation is done in the Lagrange form, which is fairly
+ * well-conditioned as ways of writing the interpolating polynomial go.
+ * Even so, the interpolation is still highly ill-conditioned for small
+ *  dx  and/or  large  order  and/or  x  outside the domain of the data
+ * points.
+ *
+ * Arguments:
+ * N_points = (in) The number of data points.
+ * order = (in) The order of the local interpolation, i.e. the degree
+ *		of the interpolated polynomial.
+ * end_action = (in) Specifies what action to take if the (centered)
+ *		     local interpolant would need data from outside
+ *		     the domain of the data points.
+ *		     'e' ==> Do an  error_exit() .
+ *		     'o' ==> Off-center the local interpolant as
+ *			     necessary to keep it inside the domain
+ *			     of the data points.
+ * xx[0...N_points) = (in) The x data array.  This is assumed to be in
+ *			   (increasing) order.
+ * yy[0...N_points) = (in) The y data array.
+ * x = (in) The x coordinate for the interpolation.
+ */
+PUBLIC double interpolate_nu_global(int N_points, int order,
+				    char end_action,
+				    double xx[],
+				    double yy[],
+				    double x)
+{
+/*
+ * Implementation note:
+ *
+ * We assume this function will typically be used to interpolate the
+ * same  xx[]  and  yy[]  values to a number of different  x  values,
+ * and that these  x  values will typically change monotonically.
+ * Therefore, we use a sequential search to locate  x  in the  xx[]
+ * array, starting from a cached position and going up or down as
+ * appropriate.
+ */
+static int cache_valid = FALSE;
+static int cache_N_points;	/* cache valid only if  N_points == this  */
+static double *cache_xx;	/*                      && xx == this */
+static int cache_i;
+
+int N_interp, i, i_lo, i_hi, shift;
+boolean ok;
+
+if (cache_valid && (N_points == cache_N_points) && (xx == cache_xx))
+   then i = cache_i;			/* use cache */
+   else i = 0;				/* default if cache isn't available */
+
+/*
+ * Find  i  such that (roughly speaking)  xx[i] <= x < xx[i+1] :
+ *
+ * More precisely, there are 3 cases to consider:
+ *	x < xx[0]			==> i = 0
+ *	xx[0] <= x <= xx[N_points-1]	==> i >= 0 && i+1 < N_points
+ *					    && xx[i] <= x < xx[i+1]
+ *	x > xx[N_points-1]		==> i = N_points-1
+ */
+if (x >= xx[i])
+   then {
+	/* search up from starting position */
+		for (++i ; ((i < N_points) && (x >= xx[i])) ; ++i)
+		{
+		}
+	--i;
+	}
+   else {
+	/* search down from starting position */
+		for (--i ; ((i >= 0) && (x < xx[i])) ; --i)
+		{
+		}
+	++i;
+	}
+
+/*
+ * firewall: sanity-check that we found  i  correctly
+ */
+if ((x >= xx[0]) && (x <= xx[N_points-1]))
+   then ok = (i >= 0) && (i+1 < N_points) && (x >= xx[i]) && (x < xx[i+1]);
+else if (x < xx[0])
+   then ok = (i == 0);
+else if (x > xx[N_points-1])
+   then ok = (i == N_points-1);
+else    ok = FALSE;
+if (! ok)
+   then error_exit(PANIC_EXIT,
+"***** interpolate_nu_global: search result failed consistency check!\n"
+"                             (this should never happen!)\n"
+"                             N_points=%d\n"
+"                             xx[0]=%g   xx[N_points-1]=%g\n"
+"                             x=%g\n"
+"                             ==> i=%d\n"
+"                                 xx[i]=%g   xx[i+1]=%g\n"
+,
+		   N_points,
+		   xx[0], xx[N_points-1],
+		   x,
+		   i,
+		   xx[i], xx[i+1]);				/*NOTREACHED*/
+
+/*
+ * update cache for future use
+ */
+cache_valid = TRUE;
+cache_N_points = N_points;
+cache_xx = xx;
+cache_i = i;
+
+/*
+ * locate subarray
+ */
+N_interp = order + 1;
+subarray_posn(N_interp, i, &i_lo, &i_hi);
+if (! ((i_lo >= 0) && (i_hi < N_points)) )
+   then {
+	/* centered local interpolation would need data */
+	/* from outside  y[]  array */
+	switch	(end_action)
+		{
+	case 'e':
+		error_exit(ERROR_EXIT,
+"***** interpolate_nu_global: local interpolation needs data\n"
+"                             from outside data array!\n"
+"                             N_points=%d   order=%d   N_interp=%d\n"
+"                             xx[0]=%g   xx[N_points-1]=%g\n"
+"                             x=%g\n"
+"                             ==> i=%d   i_[lo,hi]=[%d,%d]\n"
+,
+			   N_points, order, N_interp,
+			   xx[0], xx[N_points-1],
+			   x,
+			   i, i_lo, i_hi);			/*NOTREACHED*/
+
+	case 'o':
+		/* off-center the local interpolation as necessary */
+		/* to keep it within the  y[]  array */
+		shift = 0;
+		if (i_lo < 0) then shift = 0 - i_lo;
+		if (i_hi >= N_points) then shift = (N_points-1) - i_hi;
+
+		i_lo += shift;
+		i_hi += shift;
+
+		if (! ((i_lo >= 0) && (i_hi < N_points)) )
+		   then error_exit(ERROR_EXIT,
+"***** interpolate_nu_global: shifted local interpolation needs data\n"
+"                             from outside data array on other side!\n"
+"                             (this should only happen\n"
+"                              if N_interp > N_points!)\n"
+"                             N_points=%d   order=%d   N_interp=%d\n"
+"                             xx[0]=%g   xx[N_points-1]=%g\n"
+"                             x=%g\n"
+"                             ==> i=%d\n"
+"                             ==> shift=%d   i_[lo,hi]=[%d,%d]\n"
+,
+				   N_points, order, N_interp,
+				   xx[0], xx[N_points-1],
+				   x,
+				   i,
+				   shift, i_lo, i_hi);		/*NOTREACHED*/
+		break;
+
+	default:
+		error_exit(PANIC_EXIT,
+"***** interpolate_nu_global: unknown end_action=0x%02x='%c'\n",
+			   (int) end_action, end_action);	/*NOTREACHED*/
+		}
+	}
+
+/*
+ * do local interpolation
+ */
+return(
+   interpolate_nu_local(order,
+			xx + i_lo,
+			yy + i_lo,
+			x)
+      );
+}
+
+/*****************************************************************************/
+
+/*
+ * This function does local Lagrange polynomial interpolation on a
+ * nonuniform grid.  That is, given  order+1  data points  (xx[i],yy[i]) ,
+ *  i = 0...order, it interpolates a degree-(order) (Lagrange) polynomial
+ * through them and evaluates this polynomial at a specified  x  coordinate.
+ * In general the error in this computed function value is  O(h^order) .
+ *
+ * The interpolation is done in the Lagrange form, which is fairly
+ * well-conditioned as ways of writing the interpolating polynomial go.
+ * Even so, the interpolation is still highly ill-conditioned for small
+ *  dx  and/or  large  order  and/or  x  outside the domain of the data
+ * points.
+ *
+ * Arguments:
+ * N_points = (in) The number of data points.
+ * order = (in) The order of the local interpolation, i.e. the degree
+ *		of the interpolated polynomial.
+ * xx[0...N_points) = (in) The x data array.  This is assumed to be in
+ *			  (increasing) order.
+ * yy[0...N_points) = (in) The y data array.
+ * x = (in) The x coordinate for the interpolation.
+ */
+PUBLIC double interpolate_nu_local(int order,
+				   double xx[],
+				   double yy[],
+				   double x)
+{
+switch	(order)
+	{
+case 3:
+	return(  interpolate_nu_local_order3(xx, yy, x)  );
+	break;
+
+case 4:
+	return(  interpolate_nu_local_order4(xx, yy, x)  );
+	break;
+
+case 5:
+	return(  interpolate_nu_local_order5(xx, yy, x)  );
+	break;
+
+default:
+	error_exit(ERROR_EXIT,
+"***** interpolate_nu_local: order=%d not implemented yet!\n",
+		   order);					/*NOTREACHED*/
+	}
+}
+
+/*****************************************************************************/
+
+/*
+ * This function does 4-point local Lagrange polynomial interpolation
+ * on a nonuniform grid, interpolating a 3rd-order polynomial.
+ */
+PRIVATE double interpolate_nu_local_order3(double xx[],
+					   double yy[],
+					   double x)
+{
+double num0 = ( x    - xx[1]) * ( x    - xx[2]) * ( x    - xx[3]);
+double den0 = (xx[0] - xx[1]) * (xx[0] - xx[2]) * (xx[0] - xx[3]);
+double c0 = num0 / den0;
+
+double num1 = ( x    - xx[0]) * ( x    - xx[2]) * ( x    - xx[3]);
+double den1 = (xx[1] - xx[0]) * (xx[1] - xx[2]) * (xx[1] - xx[3]);
+double c1 = num1 / den1;
+
+double num2 = ( x    - xx[0]) * ( x    - xx[1]) * ( x    - xx[3]);
+double den2 = (xx[2] - xx[0]) * (xx[2] - xx[1]) * (xx[2] - xx[3]);
+double c2 = num2 / den2;
+
+double num3 = ( x    - xx[0]) * ( x    - xx[1]) * ( x    - xx[2]);
+double den3 = (xx[3] - xx[0]) * (xx[3] - xx[1]) * (xx[3] - xx[2]);
+double c3 = num3 / den3;
+
+return(  c0*yy[0] + c1*yy[1] + c2*yy[2] + c3*yy[3]  );
+}
+
+/*****************************************************************************/
+
+/*
+ * This function does 5-point local Lagrange polynomial interpolation
+ * on a nonuniform grid, interpolating a 4th-order polynomial.
+ */
+PRIVATE double interpolate_nu_local_order4(double xx[],
+					   double yy[],
+					   double x)
+{
+double num0 = ( x   -xx[1]) * ( x   -xx[2]) * ( x   -xx[3]) * ( x   -xx[4]);
+double den0 = (xx[0]-xx[1]) * (xx[0]-xx[2]) * (xx[0]-xx[3]) * (xx[0]-xx[4]);
+double c0 = num0 / den0;
+
+double num1 = ( x   -xx[0]) * ( x   -xx[2]) * ( x   -xx[3]) * ( x   -xx[4]);
+double den1 = (xx[1]-xx[0]) * (xx[1]-xx[2]) * (xx[1]-xx[3]) * (xx[1]-xx[4]);
+double c1 = num1 / den1;
+
+double num2 = ( x   -xx[0]) * ( x   -xx[1]) * ( x   -xx[3]) * ( x   -xx[4]);
+double den2 = (xx[2]-xx[0]) * (xx[2]-xx[1]) * (xx[2]-xx[3]) * (xx[2]-xx[4]);
+double c2 = num2 / den2;
+
+double num3 = ( x   -xx[0]) * ( x   -xx[1]) * ( x   -xx[2]) * ( x   -xx[4]);
+double den3 = (xx[3]-xx[0]) * (xx[3]-xx[1]) * (xx[3]-xx[2]) * (xx[3]-xx[4]);
+double c3 = num3 / den3;
+
+double num4 = ( x   -xx[0]) * ( x   -xx[1]) * ( x   -xx[2]) * ( x   -xx[3]);
+double den4 = (xx[4]-xx[0]) * (xx[4]-xx[1]) * (xx[4]-xx[2]) * (xx[4]-xx[3]);
+double c4 = num4 / den4;
+
+return(
+   c0*yy[0] + c1*yy[1] + c2*yy[2] + c3*yy[3] + c4*yy[4]
+      );
+}
+
+/*****************************************************************************/
+
+/*
+ * This function does 6-point local Lagrange polynomial interpolation
+ * on a nonuniform grid, interpolating a 5th-order polynomial.
+ */
+PRIVATE double interpolate_nu_local_order5(double xx[],
+					   double yy[],
+					   double x)
+{
+double num0
+   = ( x   -xx[1])*( x   -xx[2])*( x   -xx[3])*( x   -xx[4])*( x   -xx[5]);
+double den0
+   = (xx[0]-xx[1])*(xx[0]-xx[2])*(xx[0]-xx[3])*(xx[0]-xx[4])*(xx[0]-xx[5]);
+double c0 = num0 / den0;
+
+double num1
+   = ( x   -xx[0])*( x   -xx[2])*( x   -xx[3])*( x   -xx[4])*( x   -xx[5]);
+double den1
+   = (xx[1]-xx[0])*(xx[1]-xx[2])*(xx[1]-xx[3])*(xx[1]-xx[4])*(xx[1]-xx[5]);
+double c1 = num1 / den1;
+
+double num2
+   = ( x   -xx[0])*( x   -xx[1])*( x   -xx[3])*( x   -xx[4])*( x   -xx[5]);
+double den2
+   = (xx[2]-xx[0])*(xx[2]-xx[1])*(xx[2]-xx[3])*(xx[2]-xx[4])*(xx[2]-xx[5]);
+double c2 = num2 / den2;
+
+double num3
+   = ( x   -xx[0])*( x   -xx[1])*( x   -xx[2])*( x   -xx[4])*( x   -xx[5]);
+double den3
+   = (xx[3]-xx[0])*(xx[3]-xx[1])*(xx[3]-xx[2])*(xx[3]-xx[4])*(xx[3]-xx[5]);
+double c3 = num3 / den3;
+
+double num4
+   = ( x   -xx[0])*( x   -xx[1])*( x   -xx[2])*( x   -xx[3])*( x   -xx[5]);
+double den4
+   = (xx[4]-xx[0])*(xx[4]-xx[1])*(xx[4]-xx[2])*(xx[4]-xx[3])*(xx[4]-xx[5]);
+double c4 = num4 / den4;
+
+double num5
+   = ( x   -xx[0])*( x   -xx[1])*( x   -xx[2])*( x   -xx[3])*( x   -xx[4]);
+double den5
+   = (xx[5]-xx[0])*(xx[5]-xx[1])*(xx[5]-xx[2])*(xx[5]-xx[3])*(xx[5]-xx[4]);
+double c5 = num5 / den5;
+
+return(
+   c0*yy[0] + c1*yy[1] + c2*yy[2] + c3*yy[3] + c4*yy[4] + c5 * yy[5]
+      );
+}
+
+/*****************************************************************************/
+/*****************************************************************************/
+/*****************************************************************************/
+
+/*
+ * This function computes subarray positions for an interpolation or
+ * similar operation.  That is, it encapsulates the tricky odd/even
+ * computations needed to find the exact boundaries of the subarray in
+ * which an interpolation or similar operation should be done.
+ *
+ * Suppose we have a data array [i_min ... i_max] from which we wish to
+ * select an N point subarray for interpolation, and the subarray is to
+ * be centered between [i] and [i+1].  Then if N is even, the centering
+ * is unique, but if N is odd, we have to choose either the i-greater or
+ * i-less side on which to add the last point.  We choose the former here,
+ * as would be appropriate for functions varying most rapidly at smaller
+ * [i] and less rapidly at larger [i].  Thus, for example, we have
+ *
+ *					   center
+ *					    here
+ *   N   i_lo   i_hi   [i-2]   [i-1]   [i+0] | [i+1]   [i+2]   [i+3]   [i+4]
+ *   -   ----   ----   -----   -----   ----- | -----   -----   -----   -----
+ *   2   i-0    i+1                      x   |   x
+ *   3   i-0    i+2                      x   |   x       x
+ *   4   i-1    i+2              x       x   |   x       x
+ *   5   i-1    i+3              x       x   |   x       x       x
+ *   6   i-2    i+3      x       x       x   |   x       x       x
+ *   7   i-2    i+4      x       x       x   |   x       x       x       x
+ *
+ * The purpose of this function is to compute the bounds of the subarray,
+ *  i_lo  and  i_hi .  This function intentionally doesn't check for these
+ * overflowing the data array, because our caller can presumably provide
+ * a more useful error message if this occurs.
+ *
+ * Arguments:
+ * N = (in) The number of points in the subarray.
+ * i = (in) The subarray is to be centered between [i] and [i+1].
+ * pi_{lo,hi} --> (out) The inclusive {lower,upper} bound of the subarray.
+ *			Either or both of these arguments may be NULL, in
+ *			which case the corresponding result isn't stored.
+ *
+ * Results:
+ * The function returns i_lo, the inclusive lower bound of the subarray.
+ */
+PUBLIC int subarray_posn(int N, int i, int *pi_lo, int *pi_hi)
+{
+int i_lo = i + 1 - (N >> 1);
+int i_hi = i + ((N+1) >> 1);
+
+if (HOW_MANY_IN_RANGE(i_lo,i_hi) != N)
+   then error_exit(PANIC_EXIT,
+"***** interpolate_posn: subarray has wrong number of points!\n"
+"                        (this should never happen!)\n"
+"                        N=%d   i=%d   i_lo=%d   i_hi=%d\n"
+"                        HOW_MANY_IN_RANGE(i_lo,i_hi)=%d\n"
+,
+		   N, i, i_lo, i_hi,
+		   HOW_MANY_IN_RANGE(i_lo,i_hi));		/*NOTREACHED*/
+
+if (pi_lo != NULL)
+   then *pi_lo = i_lo;
+if (pi_hi != NULL)
+   then *pi_hi = i_hi;
+
+return(i_lo);
+}
diff --git a/src/make.code.defn b/src/make.code.defn
new file mode 100644
index 0000000..0483ce6
--- /dev/null
+++ b/src/make.code.defn
@@ -0,0 +1,9 @@
+# Main make.code.defn file for thorn Cartoon2D
+# $Header$
+
+# Source files in this directory
+SRCS = 
+
+# Subdirectories containing source files
+SUBDIRS = 
+