Recreating code after origin /data crash.

git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/Cartoon2D/trunk@6 eec4d7dc-71c2-46d6-addf-10296150bf52

1 files changed, 73 insertions, 999 deletions

diff --git a/src/interpolate.c b/src/interpolate.c
index 346a1fa..fe8cbb0 100644
--- a/src/interpolate.c
+++ b/src/interpolate.c
@@ -1,247 +1,45 @@
 /* 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
  */
 
+/* To make porting easy, only the interpolate_local* functions 
+ * have been retained in this file. JT's original code (kept in 
+ * the archive directory) also has functions for derivatives and
+ * interpolation on non-uniform grids. ---Sai.
+*/
 
+#include "cctk.h"
 
-/* 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
+CCTK_REAL interpolate_local(CCTK_INT order, CCTK_REAL x0, CCTK_REAL dx, 
+			 CCTK_REAL y[], CCTK_REAL x);
 
 /* 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);
+PRIVATE CCTK_REAL interpolate_local_order1(CCTK_REAL x0, CCTK_REAL dx,
+					CCTK_REAL y[],
+					CCTK_REAL x);
+PRIVATE CCTK_REAL interpolate_local_order2(CCTK_REAL x0, CCTK_REAL dx,
+					CCTK_REAL y[],
+					CCTK_REAL x);
+PRIVATE CCTK_REAL interpolate_local_order3(CCTK_REAL x0, CCTK_REAL dx,
+					CCTK_REAL y[],
+					CCTK_REAL x);
+PRIVATE CCTK_REAL interpolate_local_order4(CCTK_REAL x0, CCTK_REAL dx,
+					CCTK_REAL y[],
+					CCTK_REAL x);
+PRIVATE CCTK_REAL interpolate_local_order5(CCTK_REAL x0, CCTK_REAL dx,
+					CCTK_REAL y[],
+					CCTK_REAL 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)
-      );
-}
 
 /*****************************************************************************/
 
@@ -269,10 +67,10 @@ return(
  * 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)
+PUBLIC CCTK_REAL interpolate_local(CCTK_INT order,
+				CCTK_REAL x0, CCTK_REAL dx,
+				CCTK_REAL y[],
+				CCTK_REAL x)
 {
 switch	(order)
 	{
@@ -295,11 +93,6 @@ case 4:
 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*/
 	}
 }
 
@@ -309,15 +102,15 @@ default:
  * 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)
+PRIVATE CCTK_REAL interpolate_local_order1(CCTK_REAL x0, CCTK_REAL dx,
+					CCTK_REAL y[],
+					CCTK_REAL x)
 {
-double c0 = (+1.0/2.0)*y[0] + (+1.0/2.0)*y[1];
-double c1 =          - y[0] +          + y[1];
+CCTK_REAL c0 = (+1.0/2.0)*y[0] + (+1.0/2.0)*y[1];
+CCTK_REAL c1 =          - y[0] +          + y[1];
 
-double xc = x0 + 0.5*dx;
-double xr = (x - xc) / dx;
+CCTK_REAL xc = x0 + 0.5*dx;
+CCTK_REAL xr = (x - xc) / dx;
 
 return(  c0 + xr*c1  );
 }
@@ -328,16 +121,16 @@ 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)
+PRIVATE CCTK_REAL interpolate_local_order2(CCTK_REAL x0, CCTK_REAL dx,
+					CCTK_REAL y[],
+					CCTK_REAL 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];
+CCTK_REAL c0 = y[1];
+CCTK_REAL c1 = (-1.0/2.0)*y[0]        + (+1.0/2.0)*y[2];
+CCTK_REAL 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;
+CCTK_REAL xc = x0 + 1.0*dx;
+CCTK_REAL xr = (x - xc) / dx;
 
 return(  c0 + xr*(c1 + xr*c2)  );
 }
@@ -348,21 +141,21 @@ 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)
+PRIVATE CCTK_REAL interpolate_local_order3(CCTK_REAL x0, CCTK_REAL dx,
+					CCTK_REAL y[],
+					CCTK_REAL x)
 {
-double c0 =   (-1.0/16.0)*y[0] + (+9.0/16.0)*y[1]
+CCTK_REAL 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]
+CCTK_REAL 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]
+CCTK_REAL 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]
+CCTK_REAL 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;
+CCTK_REAL xc = x0 + 1.5*dx;
+CCTK_REAL xr = (x - xc) / dx;
 
 return(  c0 + xr*(c1 + xr*(c2 + xr*c3))  );
 }
@@ -373,22 +166,22 @@ 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)
+PRIVATE CCTK_REAL interpolate_local_order4(CCTK_REAL x0, CCTK_REAL dx,
+					CCTK_REAL y[],
+					CCTK_REAL x)
 {
-double c0 = y[2];
-double c1 =   (+1.0/12.0)*y[0] + (-2.0/ 3.0)*y[1]
+CCTK_REAL c0 = y[2];
+CCTK_REAL 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]
+CCTK_REAL 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]
+CCTK_REAL 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]
+CCTK_REAL 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;
+CCTK_REAL xc = x0 + 2.0*dx;
+CCTK_REAL xr = (x - xc) / dx;
 
 return(  c0 + xr*(c1 + xr*(c2 + xr*(c3 + xr*c4)))  );
 }
@@ -399,751 +192,32 @@ 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)
+PRIVATE CCTK_REAL interpolate_local_order5(CCTK_REAL x0, CCTK_REAL dx,
+					CCTK_REAL y[],
+					CCTK_REAL x)
 {
-double c0 =   (+ 3.0/256.0)*y[0] + (-25.0/256.0)*y[1]
+CCTK_REAL 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]
+CCTK_REAL 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]
+CCTK_REAL 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]
+CCTK_REAL 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]
+CCTK_REAL 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]
+CCTK_REAL 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;
+CCTK_REAL xc = x0 + 2.5*dx;
+CCTK_REAL 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);
-}