Largely rewrite to work with mesh refinement.HEADmaster

As the actual work is now mostly done in the evolution code, this really
only handles the symmetry condition on the x-axis.

8 files changed, 435 insertions, 1002 deletions

diff --git a/interface.ccl b/interface.ccl
index f1df4ec..da70202 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -2,6 +2,7 @@
 # $Header$
 
 implements: cartoon2d
+inherits: grid
 
 INCLUDES HEADER: Cartoon2D_tensors.h in Cartoon2D_tensors.h
 INCLUDES HEADER: Cartoon2D.h in Cartoon2D.h
@@ -81,6 +82,15 @@ REQUIRES FUNCTION Boundary_SelectedGVs
 
 private:
 
+CCTK_REAL coord type=GF
+{
+    rho
+    zero
+    rxx, ryy, rxy, ryx
+    sxx, syy, sxy, syx
+}
+
+
 int excision_variables type=SCALAR
 {
   excision_active
diff --git a/schedule.ccl b/schedule.ccl
index 0ca570f..47a06eb 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -5,6 +5,20 @@ if (cartoon_active)
 {
   STORAGE: excision_variables
 
+  STORAGE: coord
+
+  schedule Cartoon2D_setup_coord at CCTK_POSTREGRID {
+    LANG: C
+  } ""
+
+  schedule Cartoon2D_setup_coord at CCTK_POSTREGRIDINITIAL {
+    LANG: C
+  } ""
+
+  schedule Cartoon2D_setup_coord at CCTK_BASEGRID {
+    LANG: C
+  } ""
+
   schedule Cartoon2D_CheckTensorTypes at CCTK_PARAMCHECK
   {
     LANG: C
diff --git a/src/ApplyCartoon.c b/src/ApplyCartoon.c
index b97c289..e5f3bc0 100644
--- a/src/ApplyCartoon.c
+++ b/src/ApplyCartoon.c
@@ -84,189 +84,86 @@ void Cartoon_ApplyBoundaries(CCTK_ARGUMENTS);
 
 void Cartoon_ApplyBoundaries(CCTK_ARGUMENTS) 
 {
+    DECLARE_CCTK_ARGUMENTS;
+    DECLARE_CCTK_PARAMETERS;
 
-  DECLARE_CCTK_ARGUMENTS;
-  DECLARE_CCTK_PARAMETERS;
+    int num_vars, err, i;
+    CCTK_INT * vars;
 
-  int num_vars, err, i, gi, group_tags_table;
-  CCTK_INT * vars;
-  char tensortype[TENSORTYPE_BUFF_SIZE];
-  char prolongtype[PROLONG_BUFF_SIZE];
-  int prolongmethod = -1;
+    //return;
 
-  int whiskycartoon, len;
+    /* Allocate memory to hold selected bcs */
+    num_vars = Boundary_SelectedGVs(cctkGH, 0, NULL, NULL, NULL, NULL, NULL);
+    vars = malloc(num_vars*sizeof *vars);
 
-  
-
-
-  /* Check grid size */
-  if(cctk_gsh[1] != 2*cctk_nghostzones[1]+1)
-    {
-      CCTK_WARN(0, "Grid size in y direction inappropriate.");
-    }
-
-  /* Allocate memory to hold selected bcs */
-  num_vars = Boundary_SelectedGVs(cctkGH, 0, NULL, NULL, NULL, NULL, NULL);
-#ifdef DEBUG
-  printf("Cartoon_ApplyBoundaries: num_vars is %d\n",num_vars);
-#endif
-  vars = malloc(num_vars*sizeof *vars);
-
-  /* get all selected vars */
-  err = Boundary_SelectedGVs(cctkGH, num_vars, vars, NULL, NULL, NULL, NULL);
-  if (err != num_vars)
-    {
-      CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
-		 "Boundary_SelectedGVs returned %d selected variables, but %d "
-		 "expected", err, num_vars);
+    /* get all selected vars */
+    err = Boundary_SelectedGVs(cctkGH, num_vars, vars, NULL, NULL, NULL, NULL);
+    if (err != num_vars) {
+        CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                   "Boundary_SelectedGVs returned %d selected variables, but %d "
+                   "expected", err, num_vars);
     }
 
-  /* Apply CartoonBC to each of them. */
-  /* One should probably check to see that the entire group has been
-   * selected, since Cartoon operates on an entire group at a time.
-   * For now I'll just skip a variable if it is not a 'group leader'.
-   */
-  for (i=0; i<num_vars; ++i) {
-#ifdef DEBUG
-    printf("Cartoon_ApplyBoundaries: i=%d applying cartoon to vi %d\n",i,
-	   vars[i]);
-#endif
-    gi = CCTK_GroupIndexFromVarI(vars[i]);
-    if (gi<0) {
-      CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
-		 "Invalid variable index %d selected for a boundary "
-		 "condition", gi);
-    }
-    if (vars[i] != CCTK_FirstVarIndexI(gi)) 
-      {
-	/* not a group leader -- skip to next var index */
-	continue;
-      }
-
-    /* Here one should check that the entire group is registered,
-     * using CCTK_NumVarsInGroupI.
+    /* Apply CartoonBC to each of them. */
+    /* One should probably check to see that the entire group has been
+     * selected, since Cartoon operates on an entire group at a time.
+     * For now I'll just skip a variable if it is not a 'group leader'.
      */
+    for (i = 0; i < num_vars; i++) {
+        char tensortype[TENSORTYPE_BUFF_SIZE];
+        char prolongtype[PROLONG_BUFF_SIZE];
+        int group_tags_table, ret;
 
-    /* Get table handle for group tags table */
-    group_tags_table = CCTK_GroupTagsTableI(gi);
-    if (group_tags_table<0) {
-      CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
-		 "Tags table for variable group %s is %d", CCTK_GroupName(gi),
-		 group_tags_table);
-    } 
-   
-
-    /* This lines are writen for Whisky2D to avoid to treat the hydrovariables with Cartoon */
- 
-    len = Util_TableGetString(group_tags_table,0,NULL,"whiskycartoon");
-
-    whiskycartoon = 1;          /* default is yes */
-    if (len >= 0)
-      {
-        char* value = malloc (len + 1);
-        Util_TableGetString (group_tags_table, len + 1, value, "whiskycartoon");
-        if (CCTK_Equals(value, "yes"))
-          {
-            whiskycartoon = 1;
-          }
-        else if (CCTK_Equals(value, "no"))
-          {
-            whiskycartoon = 0;
-          }
-        else
-          {
+        int gi = CCTK_GroupIndexFromVarI(vars[i]);
+        if (gi < 0) {
             CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
-                       "Tags table entry \"whiskycartoon\" for variable group %s should be either \"yes\" or \"no\", but is \"%s\"",
-                       CCTK_GroupName(gi), value);
-          }
-        free (value);
-      }
-
-    if (whiskycartoon)
-      {
-  
-	/*###################################################################*/
-
-
-	/* Get tensor type from group tags table */
-	err = Util_TableGetString(group_tags_table, TENSORTYPE_BUFF_SIZE,
-				  tensortype, "tensortypealias");
-	if (err<0) 
-	  {
-	    CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
-		       "Error in TAGS table for variable group %s",
-		       CCTK_GroupName(gi));
-	  }
-#ifdef DEBUG
-	printf("Cartoon_ApplyBoundaries: tensor type is %s\n",tensortype);
-#endif
-
-	/* Get prolongation type from group tags table */
-
-	err = Util_TableGetString(group_tags_table, PROLONG_BUFF_SIZE,
-				  prolongtype, "Prolongation");
-
-	if (err == UTIL_ERROR_TABLE_NO_SUCH_KEY) 
-	  {
-	    /* Use the default */
-
-	    prolongmethod = PROLONG_LAGRANGE;
-	  }
-	else if (err < 0)
-	  {
-	    CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,
-		       "Error (%d) in TAGS table for variable group %s", err,
-		       CCTK_GroupName(gi));
-	  }
-	else
-	  {
-	    if (CCTK_Equals(prolongtype, "None")) 
-	      {
-		prolongmethod = PROLONG_NONE; /* But why? */
-	      } 
-	    else if (CCTK_Equals(prolongtype, "Lagrange")) 
-	      {
-		prolongmethod = PROLONG_LAGRANGE;
-	      } 
-	    else if (CCTK_Equals(prolongtype, "TVD")) 
-	      {
-		prolongmethod = PROLONG_ENO;
-	      } 
-	    else if (CCTK_Equals(prolongtype, "ENO")) 
-	      {
-		prolongmethod = PROLONG_ENO;
-	      } 
-	    else 
-	      {
-		CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,
-			   "Error in TAGS table for variable group %s",
-			   CCTK_GroupName(gi));
-	      }
-	  }    
-    
-	/* Here one should check that the group's size is correct for the
-	 * specified tensortype.
-	 */
+                       "Invalid variable index %d selected for a boundary "
+                       "condition", gi);
+        }
+        if (vars[i] != CCTK_FirstVarIndexI(gi)) {
+            /* not a group leader -- skip to next var index */
+            continue;
+        }
+
+        /* Get table handle for group tags table */
+        group_tags_table = CCTK_GroupTagsTableI(gi);
+        if (group_tags_table < 0) {
+            CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                       "Tags table for variable group %s is %d", CCTK_GroupName(gi),
+                       group_tags_table);
+        }
+
+        /* Get tensor type from group tags table */
+        ret = Util_TableGetString(group_tags_table, TENSORTYPE_BUFF_SIZE,
+                                  tensortype, "tensortypealias");
+        if (ret < 0) {
+            CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                       "Error in TAGS table for variable group %s",
+                       CCTK_GroupName(gi));
+        }
+
+        /* Get prolongation type from group tags table */
+        ret = Util_TableGetString(group_tags_table, PROLONG_BUFF_SIZE,
+                                  prolongtype, "Prolongation");
+
+        if (!(ret == UTIL_ERROR_TABLE_NO_SUCH_KEY ||
+              (ret >= 0 && CCTK_Equals(prolongtype, "Lagrange"))))
+            continue;
+
+        /* Call BndCartoon2DVI, passing the appropriate tensor type integer
+           macro */
+        if (CCTK_Equals(tensortype, "scalar")) {
+            BndCartoon2DVI(cctkGH, TENSORTYPE_SCALAR, vars[i]);
+        } else if ((CCTK_Equals(tensortype, "u")) || (CCTK_Equals(tensortype, "d"))) {
+            BndCartoon2DVI(cctkGH, TENSORTYPE_U, vars[i]);
+        } else if (CCTK_Equals(tensortype, "dd_sym")) {
+            BndCartoon2DVI(cctkGH, TENSORTYPE_DDSYM, vars[i]);
+        } else {
+            CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,
+                       "invalid tensor type for group %s", CCTK_GroupName(gi));
+        }
+    }
 
-	/* Call BndCartoon2DVI, passing the appropriate tensor type integer 
-	   macro */
-	if (CCTK_Equals(tensortype, "scalar"))
-	  {
-	    BndCartoon2DVI(cctkGH, TENSORTYPE_SCALAR, prolongmethod, vars[i]);
-	  } else if ((CCTK_Equals(tensortype, "u")) ||
-		     (CCTK_Equals(tensortype, "d")))
-	  {
-	    BndCartoon2DVI(cctkGH, TENSORTYPE_U, prolongmethod, vars[i]);
-	  } else if (CCTK_Equals(tensortype, "dd_sym"))
-	  {
-	    BndCartoon2DVI(cctkGH, TENSORTYPE_DDSYM, prolongmethod, vars[i]);
-	  } else 
-	  {
-	    CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,
-		       "invalid tensor type for group %s", CCTK_GroupName(gi));
-	  }  
-      }
-  }
-  /* Free data */
-  free(vars);
+    /* Free data */
+    free(vars);
 }
diff --git a/src/Cartoon2D.h b/src/Cartoon2D.h
index b4fadc3..d371729 100644
--- a/src/Cartoon2D.h
+++ b/src/Cartoon2D.h
@@ -7,7 +7,7 @@
 extern "C" {
 #endif
 
-int BndCartoon2DVI(const cGH *GH, int tensortype, int prolongtype, int vi);
+int BndCartoon2DVI(const cGH *cctkGH, int tensortype, int vi);
 
 CCTK_INT
 Cartoon2D_SymmetryInterpolate (CCTK_POINTER_TO_CONST const cctkGH_,
diff --git a/src/Cartoon2DBC.c b/src/Cartoon2DBC.c
index 2dbca7f..40b3503 100644
--- a/src/Cartoon2DBC.c
+++ b/src/Cartoon2DBC.c
@@ -32,18 +32,9 @@ CCTK_FILEVERSION(BetaThorns_Cartoon2D_Cartoon2DBC_c);
 
 int BndCartoon2DGN(const cGH *GH, int tensortype, const char *group);
 int BndCartoon2DVN(const cGH *GH, int tensortype, const char *var);
-int BndCartoon2DVI(const cGH *GH, int tensortype, int prolongtype, int vi);
-static CCTK_REAL Cartoon2DInterp(const cGH *GH, CCTK_REAL *f, int i, int di,
-                             int ijk, CCTK_REAL x);
-static CCTK_REAL Cartoon2DInterp_ENO(const cGH *cctkGH, 
-                                    CCTK_REAL *f, int i, int di, 
-                                    int ijk, CCTK_REAL x);
+int BndCartoon2DVI(const cGH *GH, int tensortype, int vi);
 CCTK_REAL interpolate_local(int order, CCTK_REAL x0, CCTK_REAL dx, 
-                            CCTK_REAL y[], CCTK_REAL x, 
-                            CCTK_INT excised_points[]);
-CCTK_REAL interpolate_eno(CCTK_INT order, CCTK_REAL x0, CCTK_REAL dx,
-                          CCTK_REAL y[], CCTK_REAL x,
-                          CCTK_INT excised_points[]);
+                            CCTK_REAL y[], CCTK_REAL x);
 void CCTK_FCALL CCTK_FNAME(BndCartoon2DVI)(int *retval, const cGH **GH, 
                 int *tensortype, int *vi);
 void CCTK_FCALL CCTK_FNAME(BndCartoon2DVN)(int *retval, const cGH **GH,
@@ -55,6 +46,46 @@ void CCTK_FCALL CCTK_FNAME(BndCartoon2DGN)
 void Cartoon2D_InitExcisionVars(CCTK_ARGUMENTS);
 void SetupExcisionVars(CCTK_ARGUMENTS);
 
+/* interpolation on x-axis */
+static CCTK_REAL cartoon2d_interp(const cGH *cctkGH, int order, CCTK_REAL *f,
+                                  int i, int ijk, CCTK_REAL x_val)
+{
+    CCTK_REAL x0;
+    CCTK_REAL lx0, dx0;
+    int lnx;
+    int n, offset;
+
+    lnx = cctkGH->cctk_lsh[0];
+
+    dx0 = cctkGH->cctk_delta_space[0] / cctkGH->cctk_levfac[0];
+    lx0 = cctkGH->cctk_origin_space[0] + cctkGH->cctk_delta_space[0] / cctkGH->cctk_levfac[0] * cctkGH->cctk_levoff[0] / cctkGH->cctk_levoffdenom[0] + dx0 * cctkGH->cctk_lbnd[0];
+
+
+    /* 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_val > lx0 + dx0 * (i+1)) {i++; ijk++;}
+
+    /* 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;
+    while (lx0 + dx0 * (i - offset) < 0) offset--;
+
+    if (i-offset+order >= lnx) offset = i+order-lnx+1;
+
+    x0 = lx0 + dx0 * (i-offset);
+    if(x0 < 0)
+        fprintf(stderr, "x0 %g\n", x0);
+
+    /* call interpolator */
+    return interpolate_local(order, x0, dx0, f + ijk - offset, x_val);
+}
+
+
 /* set boundaries of a grid tensor assuming axisymmetry 
    - handles lower boundary in x
    - does not touch other boundaries
@@ -75,241 +106,175 @@ void SetupExcisionVars(CCTK_ARGUMENTS);
    group, using the tensortype argument to determine how many there
    are and how to treat them. */
 
-int BndCartoon2DVI(const cGH *cctkGH, int tensortype, int prolongtype, int vi)
+int BndCartoon2DVI(const cGH *cctkGH, int tensortype, int vi)
 
 {
-  DECLARE_CCTK_ARGUMENTS;
-  DECLARE_CCTK_PARAMETERS;
-
-  char * groupname;
-
-  int i, j, k, di, dj, ijk;
-  int jj, n, s;
-  int lnx, lnz, ny;
-  int gi, var0;
-  int n_vars;
-  int timelevel;
-  CCTK_REAL x, y, rho;
-  CCTK_REAL lx0, dx0, dy0;
-  CCTK_REAL rxx, rxy, ryx, ryy;
-  CCTK_REAL sxx, sxy, syx, syy;
-  CCTK_REAL f[100], fx, fy, fz, fxx, fxy, fxz, fyy, fyz, fzz;
-  CCTK_REAL *t, *tx, *ty, *tz, *txx, *txy, *txz, *tyy, *tyz, *tzz;
-
-  gi = CCTK_GroupIndexFromVarI (vi);
-  assert (gi >= 0);
-  n_vars = CCTK_NumVarsInGroupI (gi);
-  assert (n_vars > 0);
-  var0 = CCTK_FirstVarIndexI (gi);
-  assert (var0 >= 0);
-
-  if (n_vars > 100)
-  {
-    groupname = CCTK_GroupName (gi);
-    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
-                "Group \"%s\" has more than 100 variables -- cannot apply Cartoon boundary condition",
-                groupname);
-    free (groupname);
-    return -1;
-  }
-
-  switch (tensortype) {
-  case TENSORTYPE_SCALAR:
-    break;
-  case TENSORTYPE_U:
-    if (n_vars != 3)
-    {
-      groupname = CCTK_GroupName (gi);
-      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
-                  "Group \"%s\" is declared with the tensor type U, but does not contain 3 variables -- cannot apply Cartoon boundary condition",
-                  groupname);
-      free (groupname);
-      return -1;
-    }
-    break;
-  case TENSORTYPE_DDSYM:
-    if (n_vars != 6)
-    {
-      groupname = CCTK_GroupName (gi);
-      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
-                  "Group \"%s\" is declared with the tensor type DDSYM, but does not contain 6 variables -- cannot apply Cartoon boundary condition",
-                  groupname);
-      free (groupname);
-      return -1;
+    DECLARE_CCTK_ARGUMENTS;
+    DECLARE_CCTK_PARAMETERS;
+
+    int i, s;
+    int y_idx;
+    int gi, var0;
+    int n_vars;
+    int timelevel;
+
+    gi     = CCTK_GroupIndexFromVarI (vi);
+    n_vars = CCTK_NumVarsInGroupI (gi);
+    var0   = CCTK_FirstVarIndexI (gi);
+
+    if (n_vars > 100) {
+        char *groupname = CCTK_GroupName (gi);
+        CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                    "Group \"%s\" has more than 100 variables -- cannot apply Cartoon boundary condition",
+                    groupname);
+        free (groupname);
+        return -1;
     }
-    break;
-  default:
-    CCTK_WARN(0,"Tensor type not recognized by Cartoon2D.");
-  }
-    
-  s = cctkGH->cctk_nghostzones[0];
-  lnx = cctkGH->cctk_lsh[0];
-  lnz = cctkGH->cctk_lsh[2];
-  ny  = cctkGH->cctk_gsh[1];
-
-  dx0 = CCTK_DELTA_SPACE(0);
-  dy0 = CCTK_DELTA_SPACE(1);
-  lx0 = CCTK_ORIGIN_SPACE(0) + dx0 * cctk_lbnd[0];
-
-  timelevel = 0; 
-
-  /* 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 = lnx;
-
-  /* y = 0 */
-  assert (ny % 2 == 1);
-  j = ny/2;
-
-  /* make sure that the input data is synced */
-  CCTK_SyncGroupI(cctkGH, gi);
-
-  /* z-direction: include lower and upper boundary */
-  for (k = 0; k < lnz; k++)
-
-  /* y-direction: as many zombies as the evolution stencil needs */
-#if 0
-  for (jj = 1, dj = jj*lnx; jj <= ny/2; jj++, dj = jj*lnx) 
-#else
-  for (jj = 0, dj = jj*lnx; jj <= ny/2; jj++, dj = jj*lnx) 
-#endif
-
-  /* x-direction: zombie for x < 0, including upper boundary for extrapol */
-#if 0
-  for (i = -s; i < lnx; i++) {
-#else
-  for (i = 0; i < lnx; i++)
-    if (! (i >= s && jj == 0)) {
-#endif
-
-    /* index into linear memory */
-    ijk = CCTK_GFINDEX3D(cctkGH,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 = lx0 + dx0 * i;
-#if 0
-    y = (dj) ? dy0 * jj : 0;
-#else
-    y = dy0 * jj;
-#endif
-    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 */
-    switch (prolongtype)
-    {
-      case PROLONG_NONE:
-        /* no-op */
-        return (0);
+
+    switch (tensortype) {
+    case TENSORTYPE_SCALAR:
         break;
-      case PROLONG_LAGRANGE:
-        for (n = 0; n < n_vars; n++)
-          f[n] = Cartoon2DInterp(cctkGH, cctkGH->data[var0+n][timelevel], 
-                                 i, di, ijk, rho); 
+    case TENSORTYPE_U:
+        if (n_vars != 3) {
+            char *groupname = CCTK_GroupName (gi);
+            CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                        "Group \"%s\" is declared with the tensor type U, but does not contain 3 variables -- cannot apply Cartoon boundary condition",
+                        groupname);
+            free (groupname);
+            return -1;
+        }
         break;
-      case PROLONG_ENO:
-        for (n = 0; n < n_vars; n++)
-          f[n] = Cartoon2DInterp_ENO(cctkGH, cctkGH->data[var0+n][timelevel], i, 
-                                     di, ijk, rho); 
+    case TENSORTYPE_DDSYM:
+        if (n_vars != 6) {
+            char *groupname = CCTK_GroupName (gi);
+            CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                        "Group \"%s\" is declared with the tensor type DDSYM, but does not contain 6 variables -- cannot apply Cartoon boundary condition",
+                        groupname);
+            free (groupname);
+            return -1;
+        }
         break;
-      default:
-        CCTK_WARN(0, "Prolongation type not recognized by Cartoon2D.");
-    }
-
-    /* rotate grid tensor by matrix multiplication */
-    if (tensortype == TENSORTYPE_SCALAR) {
-      /* Scalar groups can have arbitrary many components; these
-         "components" are then all scalars.  */
-      for (n = 0; n < n_vars; n++)
-      {
-        t = cctkGH->data[var0 + n][timelevel];
-        t[ijk+dj] = f[n];
-        t[ijk-dj] = f[n];
-      }
-    }
-    else if (tensortype == TENSORTYPE_U) {
-      tx = cctkGH->data[vi][timelevel];
-      ty = cctkGH->data[vi+1][timelevel];
-      tz = cctkGH->data[vi+2][timelevel];
-      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 (tensortype == TENSORTYPE_DDSYM) {
-      txx = cctkGH->data[vi][timelevel];
-      txy = cctkGH->data[vi+1][timelevel];
-      txz = cctkGH->data[vi+2][timelevel];
-      tyy = cctkGH->data[vi+3][timelevel];
-      tyz = cctkGH->data[vi+4][timelevel];
-      tzz = cctkGH->data[vi+5][timelevel];
-      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 {
-      return(-1);
+    default:
+        CCTK_WARN(0,"Tensor type not recognized by Cartoon2D.");
     }
 
-#if 0
-    /* what a neat way to hack out the zombies for x < 0, y = 0 */
-    if (dj == 0) {i -= s; ijk -= s; dj = jj*lnx;}
-#endif
-  }
+    s = cctkGH->cctk_nghostzones[0];
 
-  /* syncs needed after interpolation (actually only for x direction) */
-  CCTK_SyncGroupI(cctkGH, gi);
+    timelevel = 0;
 
-  return(0);
+    /* find the index for which y == 0 */
+    for (i = 0; i < cctkGH->cctk_lsh[1]; i++) {
+        if (fabs(y[CCTK_GFINDEX3D(cctkGH, 0, i, 0)]) < 1e-12) {
+            y_idx = i;
+            break;
+        }
+    }
+    if (i == cctkGH->cctk_lsh[1])
+        CCTK_WARN(0, "y is not 0 anywhere");
+
+    /* make sure that the input data is synced */
+    CCTK_SyncGroupI(cctkGH, gi);
+
+#pragma omp parallel for
+  for (int k = 0; k < cctkGH->cctk_lsh[2]; k++)
+      for (int jj = 0; jj < cctkGH->cctk_lsh[1]; jj++)
+          for (int i = 0; i < cctkGH->cctk_lsh[0]; i++) {
+              CCTK_REAL f[100];
+              CCTK_REAL rho_val;
+              int idx_dst, ijk;
+
+              if (jj != y_idx || i >= s)
+                  continue;
+              //if (jj == y_idx && i >= s)
+              //    continue;
+
+              /* index into linear memory */
+              ijk     = CCTK_GFINDEX3D(cctkGH, i, y_idx, k);
+              idx_dst = CCTK_GFINDEX3D(cctkGH, i, jj, k);
+
+              rho_val = rho[idx_dst];
+
+              /* interpolate grid functions */
+              for (int n = 0; n < n_vars; n++) {
+                  f[n] = cartoon2d_interp(cctkGH, order, cctkGH->data[var0+n][timelevel],
+                                          i, ijk, rho_val);
+              }
+
+              /* rotate grid tensor by matrix multiplication */
+              if (tensortype == TENSORTYPE_SCALAR) {
+                  /* Scalar groups can have arbitrary many components; these
+                     "components" are then all scalars.  */
+                  for (int n = 0; n < n_vars; n++) {
+                      CCTK_REAL *t = cctkGH->data[var0 + n][timelevel];
+                      t[idx_dst] = f[n];
+                  }
+              } else if (tensortype == TENSORTYPE_U) {
+                  CCTK_REAL *tx = cctkGH->data[vi][timelevel];
+                  CCTK_REAL *ty = cctkGH->data[vi+1][timelevel];
+                  CCTK_REAL *tz = cctkGH->data[vi+2][timelevel];
+
+                  CCTK_REAL rxx_val = rxx[idx_dst];
+                  CCTK_REAL ryy_val = ryy[idx_dst];
+                  CCTK_REAL rxy_val = rxy[idx_dst];
+                  CCTK_REAL ryx_val = ryx[idx_dst];
+
+                  CCTK_REAL fx = f[0];
+                  CCTK_REAL fy = f[1];
+                  CCTK_REAL fz = f[2];
+
+                  tx[idx_dst] = rxx_val * fx + rxy_val * fy;
+                  ty[idx_dst] = ryx_val * fx + ryy_val * fy;
+                  tz[idx_dst] = fz;
+
+                  //tx[ijk-dj] = sxx * fx + sxy * fy;
+                  //ty[ijk-dj] = syx * fx + syy * fy;
+                  //tz[ijk-dj] = fz;
+              } else if (tensortype == TENSORTYPE_DDSYM) {
+                  CCTK_REAL *txx = cctkGH->data[vi][timelevel];
+                  CCTK_REAL *txy = cctkGH->data[vi+1][timelevel];
+                  CCTK_REAL *txz = cctkGH->data[vi+2][timelevel];
+                  CCTK_REAL *tyy = cctkGH->data[vi+3][timelevel];
+                  CCTK_REAL *tyz = cctkGH->data[vi+4][timelevel];
+                  CCTK_REAL *tzz = cctkGH->data[vi+5][timelevel];
+
+                  CCTK_REAL sxx_val = sxx[idx_dst];
+                  CCTK_REAL syy_val = syy[idx_dst];
+                  CCTK_REAL sxy_val = sxy[idx_dst];
+                  CCTK_REAL syx_val = syx[idx_dst];
+
+                  CCTK_REAL fxx = f[0];
+                  CCTK_REAL fxy = f[1];
+                  CCTK_REAL fxz = f[2];
+                  CCTK_REAL fyy = f[3];
+                  CCTK_REAL fyz = f[4];
+                  CCTK_REAL fzz = f[5];
+
+                  txx[idx_dst] = fxx * sxx_val * sxx_val + 2 * fxy * sxx_val * syx_val + fyy * syx_val * syx_val;
+                  tyy[idx_dst] = fxx * sxy_val * sxy_val + 2 * fxy * sxy_val * syy_val + fyy * syy_val * syy_val;
+                  txy[idx_dst] = fxx * sxx_val * sxy_val +     fxy * (sxy_val * syx_val + sxx_val * syy_val) + fyy * syx_val * syy_val;
+                  txz[idx_dst] = fxz * sxx_val + fyz * syx_val;
+                  tyz[idx_dst] = fxz * sxy_val + fyz * syy_val;
+                  tzz[idx_dst] = 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;
+              }
+          }
+
+    /* syncs needed after interpolation (actually only for x direction) */
+    CCTK_SyncGroupI(cctkGH, gi);
+
+    return 0;
 }
 
 void CCTK_FCALL CCTK_FNAME(BndCartoon2DVI)
                           (int *retval, const cGH **GH, int *tensortype, int *vi)
 {
-  *retval = BndCartoon2DVI(*GH, *tensortype, PROLONG_LAGRANGE, *vi);
+  *retval = BndCartoon2DVI(*GH, *tensortype, *vi);
 }
 
 int BndCartoon2DVN(const cGH *GH, int tensortype, const char *inpvarname)
@@ -321,7 +286,7 @@ int BndCartoon2DVN(const cGH *GH, int tensortype, const char *inpvarname)
   if (verbose) printf("Cartoon2D called for %s\n", inpvarname);
   vi = CCTK_VarIndex(inpvarname);
 
-  return(BndCartoon2DVI(GH, tensortype, PROLONG_LAGRANGE, vi));
+  return(BndCartoon2DVI(GH, tensortype, vi));
 }
 
 void CCTK_FCALL CCTK_FNAME(BndCartoon2DVN)
@@ -338,8 +303,7 @@ int BndCartoon2DGN(const cGH *GH, int tensortype, const char *group)
 
   if (verbose && GH->cctk_iteration==1) CCTK_VInfo(CCTK_THORNSTRING,"Cartoon2D called for %s (message appears only once)", group);
   
-  return(BndCartoon2DVI(GH, tensortype, PROLONG_LAGRANGE, 
-                        CCTK_FirstVarIndex(group)));
+  return(BndCartoon2DVI(GH, tensortype, CCTK_FirstVarIndex(group)));
 }
 
 void CCTK_FCALL CCTK_FNAME(BndCartoon2DGN)
@@ -350,246 +314,6 @@ void CCTK_FCALL CCTK_FNAME(BndCartoon2DGN)
   free(group);
 }
 
-/* interpolation on x-axis */
-static CCTK_REAL Cartoon2DInterp(const cGH *cctkGH, CCTK_REAL *f, 
-                                 int i, int di, int ijk, CCTK_REAL x)
-{
-  DECLARE_CCTK_ARGUMENTS;
-  DECLARE_CCTK_PARAMETERS;
-
-  CCTK_REAL x0;
-  CCTK_REAL lx0, dx0;
-  CCTK_REAL y[6], r;
-  int lnx;
-  int n, offset;
-
-  CCTK_REAL *old_mask;
-  CCTK_INT *new_mask;
-
-  CCTK_INT excised_points[6];
-
-  if (*excision_active < 0)
-  {
-    SetupExcisionVars((cGH*)cctkGH);
-  }
-
-  if (old_excision)
-  {
-    old_mask = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, *old_mask_vi);
-  }
-  else
-  {
-    old_mask = NULL;
-  }
-  
-  if (new_excision)
-  {
-    new_mask = (CCTK_INT*)CCTK_VarDataPtrI(cctkGH, 0, *new_mask_vi);
-  }
-  else
-  {
-    new_mask = NULL;
-  }
-
-  /* 
-     If this boundary point is excised then we return 
-     without doing anything 
-  */
-
-  if (*excision_active)
-  {
-    if ((*excision_active == 1)||(*excision_active == 3))
-    {
-      if (old_mask[ijk] < 0.75)
-      {
-        return f[ijk];
-      }
-    }
-    if ((*excision_active == 2)||(*excision_active == 3))
-    {
-      if (SpaceMask_CheckStateBits(new_mask, (ijk),
-                                   *new_excision_field,
-                                   *new_excision_descriptor))
-      {
-        return f[ijk];
-      }
-    }
-  }  
-
-  lnx = cctkGH->cctk_lsh[0];
-    
-  dx0 = CCTK_DELTA_SPACE(0);
-  lx0 = CCTK_ORIGIN_SPACE(0) + dx0 * cctk_lbnd[0];
-
-  /* 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 > lx0 + dx0 * (i+1)) {i++; ijk++;}
-
-  /* first attempt to interface to JT's interpolator */
-
-  /* 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 >= lnx) offset = i+order-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 = lx0 + dx0 * (i-offset);
-  for (n = 0; n <= order; n++) {
-    y[n] = f[ijk-offset+n];
-    excised_points[n] = 0;
-    if (*excision_active)
-    {
-      if ((*excision_active == 1)||(*excision_active == 3))
-      {
-        excised_points[n] = (old_mask[ijk-offset+n] > 0.75) ? 0 : 1;
-      }
-      if ((*excision_active == 2)||(*excision_active == 3))
-      {
-        excised_points[n] |= 
-          (SpaceMask_CheckStateBits(new_mask, (ijk-offset+n),
-                                    *new_excision_field,
-                                    *new_excision_descriptor));
-      }
-    }
-  }
-  
-  /* call interpolator */
-  r = interpolate_local(order, x0, dx0, y, x, excised_points);
-  
-  return r;
-}
-
-static CCTK_REAL Cartoon2DInterp_ENO(const cGH *cctkGH, 
-                                    CCTK_REAL *f, int i, int di, 
-                                    int ijk, CCTK_REAL x)
-{
-  
-  DECLARE_CCTK_PARAMETERS;
-  DECLARE_CCTK_ARGUMENTS;
-  
-  CCTK_REAL x0;
-  CCTK_REAL lx0, dx0;
-  CCTK_REAL y[11], r;
-  int lnx;
-  int n, offset;
-  /* int tmp_order = 2; */
-
-  CCTK_REAL *old_mask;
-  CCTK_INT *new_mask;
-
-  CCTK_INT excised_points[11];
-
-  if (*excision_active < 0)
-  {
-    SetupExcisionVars((cGH*)cctkGH);
-  }
-
-  if (old_excision)
-  {
-    old_mask = (CCTK_REAL*)CCTK_VarDataPtrI(cctkGH, 0, *old_mask_vi);
-  }
-  else
-  {
-    old_mask = NULL;
-  }
-  
-  if (new_excision)
-  {
-    new_mask = (CCTK_INT*)CCTK_VarDataPtrI(cctkGH, 0, *new_mask_vi);
-  }
-  else
-  {
-    new_mask = NULL;
-  }
-
-  /* 
-     If this boundary point is excised then we return 
-     without doing anything 
-  */
-
-  if (*excision_active)
-  {
-    if ((*excision_active == 1)||(*excision_active == 3))
-    {
-      if (old_mask[ijk] < 0.75)
-      {
-        return f[ijk];
-      }
-    }
-    if ((*excision_active == 2)||(*excision_active == 3))
-    {
-      if (SpaceMask_CheckStateBits(new_mask, (ijk),
-                                   *new_excision_field,
-                                   *new_excision_descriptor))
-      {
-        return f[ijk];
-      }
-    }
-  }
-  
-  lnx = cctkGH->cctk_lsh[0];
-  
-  dx0 = CCTK_DELTA_SPACE(0);
-  lx0 = CCTK_ORIGIN_SPACE(0) + dx0 * cctk_lbnd[0];
-  
-  /* 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 > lx0 + dx0 * (i+1)) {i++; ijk++;}
-  
-  /* first attempt to interface to JT's interpolator */
-  
-  /* offset of stencil, note that rotation leads to x close to x0 
-     for large x */
-  offset = eno_order;
-/*   offset = tmp_order; */
-  
-  /* 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+4 >= lnx) offset = i+5-lnx; */
-/*   if (i-offset+2*tmp_order >= lnx) offset = i+2*tmp_order+1-lnx; */
-  if (i-offset+2*eno_order >= lnx) offset = i+2*eno_order+1-lnx;
-  
-  /* fill in info */
-  /* fills in old data near axis, but as long as stencil at axis is 
-     centered, no interpolation happens anyway */  
-  x0 = lx0 + dx0 * (i-offset);
-/*   for (n = 0; n < 2 * tmp_order + 1; ++n)  */
-  for (n = 0; n < 2 * eno_order + 1; ++n) 
-  {
-    y[n] = f[ijk-offset+n];
-    excised_points[n] = 0;
-    if (*excision_active)
-    {
-      if ((*excision_active == 1)||(*excision_active == 3))
-      {
-        excised_points[n] = (old_mask[ijk-offset+n] > 0.75) ? 1 : 0;
-      }
-      if ((*excision_active == 2)||(*excision_active == 3))
-      {
-        excised_points[n] |= 
-          (SpaceMask_CheckStateBits(new_mask, (ijk-offset+n),
-                                    *new_excision_field,
-                                    *new_excision_descriptor));
-      }
-    }
-  }
-  
-  /* call interpolator */
-/*   r = interpolate_eno(tmp_order, x0, dx0, y, x); */
-  r = interpolate_eno(eno_order, x0, dx0, y, x, excised_points);
-  
-  return r;
-}
-
 void Cartoon2D_InitExcisionVars(CCTK_ARGUMENTS)
 {
   DECLARE_CCTK_ARGUMENTS;
diff --git a/src/RegisterSym.c b/src/RegisterSym.c
index 2145df0..dbed7be 100644
--- a/src/RegisterSym.c
+++ b/src/RegisterSym.c
@@ -60,10 +60,10 @@ void Cartoon2D_RegisterSymmetries (CCTK_ARGUMENTS)
     CCTK_WARN (0, "Could not register the Cartoon boundaries -- probably some other thorn has already registered the same boundary faces for a different symmetry");
   }
   
-  ierr = SymmetryRegisterGridInterpolator
-    (cctkGH, handle, Cartoon2D_SymmetryInterpolate);
-  if (ierr < 0) {
-    CCTK_WARN (0, "Could not register the symmetry interpolator");
-  }
+  //ierr = SymmetryRegisterGridInterpolator
+  //  (cctkGH, handle, Cartoon2D_SymmetryInterpolate);
+  //if (ierr < 0) {
+  //  CCTK_WARN (0, "Could not register the symmetry interpolator");
+  //}
   
 }
diff --git a/src/SetGrid.c b/src/SetGrid.c
index 954d301..28ff2f7 100644
--- a/src/SetGrid.c
+++ b/src/SetGrid.c
@@ -9,10 +9,12 @@
    @version   $Id$
  @@*/
 
+#include <math.h>
 #include <stdio.h>
 #include <string.h>
 
 #include "cctk.h"
+#include "cctk_Arguments.h"
 #include "cctk_Parameter.h"
 #include "cctk_Parameters.h"
 
@@ -267,3 +269,29 @@ int Cartoon2D_SetGrid(void)
 
   return 0;
 }
+
+void Cartoon2D_setup_coord(CCTK_ARGUMENTS)
+{
+    DECLARE_CCTK_ARGUMENTS;
+    DECLARE_CCTK_PARAMETERS;
+
+#pragma omp parallel for
+    for (int k = 0; k < cctkGH->cctk_lsh[2]; k++)
+        for (int j = 0; j < cctkGH->cctk_lsh[1]; j++)
+            for (int i = 0; i < cctkGH->cctk_lsh[0]; i++) {
+                int idx = CCTK_GFINDEX3D(cctkGH, i, j, k);
+                CCTK_REAL x_val = x[idx], y_val = y[idx];
+                rho[idx] = sqrt(x_val * x_val + y_val * y_val);
+                zero[idx] = 0.0;
+
+                rxx[idx] = x_val / rho[idx];
+                ryx[idx] = y_val / rho[idx];
+                rxy[idx] = -ryx[idx];
+                ryy[idx] = rxx[idx];
+
+                sxx[idx] =  ryy[idx];
+                syx[idx] = -ryx[idx];
+                sxy[idx] = -rxy[idx];
+                syy[idx] =  rxx[idx];
+            }
+}
diff --git a/src/interpolate.c b/src/interpolate.c
index 3ae6b3f..151f08c 100644
--- a/src/interpolate.c
+++ b/src/interpolate.c
@@ -23,128 +23,7 @@
 
 CCTK_REAL interpolate_local(int order, CCTK_REAL x0, CCTK_REAL dx, 
                             CCTK_REAL y[],
-                            CCTK_REAL x,
-                            CCTK_INT excised_points[]);
-
-CCTK_REAL interpolate_eno(CCTK_INT order, CCTK_REAL x0, CCTK_REAL dx,
-                          CCTK_REAL y[],
-                          CCTK_REAL x,
-                            CCTK_INT excised_points[]);
-
-/* prototypes for private functions defined in this file */
-static CCTK_REAL interpolate_local_order1(CCTK_REAL x0, CCTK_REAL dx,
-                                          CCTK_REAL y[],
-                                          CCTK_REAL x);
-static CCTK_REAL interpolate_local_order2(CCTK_REAL x0, CCTK_REAL dx,
-                                          CCTK_REAL y[],
-                                          CCTK_REAL x);
-static CCTK_REAL interpolate_local_order3(CCTK_REAL x0, CCTK_REAL dx,
-                                          CCTK_REAL y[],
-                                          CCTK_REAL x);
-static CCTK_REAL interpolate_local_order4(CCTK_REAL x0, CCTK_REAL dx,
-                                          CCTK_REAL y[],
-                                          CCTK_REAL x);
-static CCTK_REAL interpolate_local_order5(CCTK_REAL x0, CCTK_REAL dx,
-                                          CCTK_REAL y[],
-                                          CCTK_REAL 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.
- */
-CCTK_REAL interpolate_local(int order,
-                            CCTK_REAL x0, CCTK_REAL dx,
-                            CCTK_REAL y[],
-                            CCTK_REAL x,
-                            CCTK_INT excised_points[])
-{
-
-  /* 
-     We assume that the excised region has one of the following 4 forms:
-     
-     o o o o o         (No excision)
-     x o o o o         (left excision)
-     o o o o x         (right excision)
-     x o o o x         (left and right excision)
-
-     Anything else (e.g., x o x o x) will go horribly wrong.
-  */
-
-  int start, end, available_order;
-  
-  start = 0;
-  end   = order;
-  
-  while ((start < order + 1)&&(excised_points[start]))
-  {
-    ++start;
-  }
-  while ((end > -1)&&(excised_points[end]))
-  {
-    --end;
-  }
-  
-  if (end < start)
-  {
-    /* The whole block is excised */
-    return y[0];
-  }
-
-  /* The possible interpolation order that can be used */
-
-  available_order = end - start;
-
-  if (available_order == 0)
-  {
-    /* There is only one non-excised point */
-    return y[start];
-  }
-
-  if (dx * start + x0 > x)
-  {
-    /* We would have to extrapolate to the left */
-    return y[start];
-  }
-  
-  if (dx * end + x0 < x)
-  {
-    /* We would have to extrapolate to the right */
-    return y[end];
-  }
-
-  switch (available_order) 
-  {
-  case 1: return interpolate_local_order1(x0 + start * dx, dx, &y[start], x);
-  case 2: return interpolate_local_order2(x0 + start * dx, dx, &y[start], x);
-  case 3: return interpolate_local_order3(x0 + start * dx, dx, &y[start], x);
-  case 4: return interpolate_local_order4(x0 + start * dx, dx, &y[start], x);
-  case 5: return interpolate_local_order5(x0 + start * dx, dx, &y[start], x);
-  default: CCTK_WARN(0, "Interpolation order not supported");
-    return 0; /* Avoid warning */
-  }
-}
+                            CCTK_REAL x);
 
 /*****************************************************************************/
 
@@ -156,13 +35,13 @@ static CCTK_REAL interpolate_local_order1(CCTK_REAL x0, CCTK_REAL dx,
                                           CCTK_REAL y[],
                                           CCTK_REAL x)
 {
-CCTK_REAL c0 = (+1.0/2.0)*y[0] + (+1.0/2.0)*y[1];
-CCTK_REAL 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];
 
-CCTK_REAL xc = x0 + 0.5*dx;
-CCTK_REAL xr = (x - xc) / dx;
+    CCTK_REAL xc = x0 + 0.5*dx;
+    CCTK_REAL xr = (x - xc) / dx;
 
-return(  c0 + xr*c1  );
+    return  c0 + xr * c1;
 }
 
 /*****************************************************************************/
@@ -175,14 +54,14 @@ static CCTK_REAL interpolate_local_order2(CCTK_REAL x0, CCTK_REAL dx,
                                           CCTK_REAL y[],
                                           CCTK_REAL x)
 {
-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];
+    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];
 
-CCTK_REAL xc = x0 + 1.0*dx;
-CCTK_REAL xr = (x - xc) / dx;
+    CCTK_REAL xc = x0 + 1.0*dx;
+    CCTK_REAL xr = (x - xc) / dx;
 
-return(  c0 + xr*(c1 + xr*c2)  );
+    return c0 + xr*(c1 + xr*c2);
 }
 
 /*****************************************************************************/
@@ -195,19 +74,19 @@ static CCTK_REAL interpolate_local_order3(CCTK_REAL x0, CCTK_REAL dx,
                                           CCTK_REAL y[],
                                           CCTK_REAL x)
 {
-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];
-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];
-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];
-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];
-
-CCTK_REAL xc = x0 + 1.5*dx;
-CCTK_REAL xr = (x - xc) / dx;
-
-return(  c0 + xr*(c1 + xr*(c2 + xr*c3))  );
+    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];
+    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];
+    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];
+    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];
+
+    CCTK_REAL xc = x0 + 1.5*dx;
+    CCTK_REAL xr = (x - xc) / dx;
+
+    return c0 + xr*(c1 + xr*(c2 + xr*c3));
 }
 
 /*****************************************************************************/
@@ -220,20 +99,20 @@ static CCTK_REAL interpolate_local_order4(CCTK_REAL x0, CCTK_REAL dx,
                                           CCTK_REAL y[],
                                           CCTK_REAL x)
 {
-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];
-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];
-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];
-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];
-
-CCTK_REAL xc = x0 + 2.0*dx;
-CCTK_REAL xr = (x - xc) / dx;
-
-return(  c0 + xr*(c1 + xr*(c2 + xr*(c3 + xr*c4)))  );
+    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];
+    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];
+    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];
+    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];
+
+    CCTK_REAL xc = x0 + 2.0*dx;
+    CCTK_REAL xr = (x - xc) / dx;
+
+    return c0 + xr*(c1 + xr*(c2 + xr*(c3 + xr*c4)));
 }
 
 /*****************************************************************************/
@@ -246,200 +125,81 @@ static CCTK_REAL interpolate_local_order5(CCTK_REAL x0, CCTK_REAL dx,
                                           CCTK_REAL y[],
                                           CCTK_REAL x)
 {
-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];
-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];
-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];
-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];
-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];
-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];
-
-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))))  );
+    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];
+    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];
+    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];
+    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];
+    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];
+    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];
+
+    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))));
 }
 
-
-/*****************************************************************************/
-
-CCTK_REAL interpolate_eno(CCTK_INT order, 
-                          CCTK_REAL x0, CCTK_REAL dx,
-                          CCTK_REAL y[],
-                          CCTK_REAL x,
-                          CCTK_INT excised_points[])
+/*
+ * 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.
+ */
+CCTK_REAL interpolate_local(int order,
+                            CCTK_REAL x0, CCTK_REAL dx,
+                            CCTK_REAL y[],
+                            CCTK_REAL x)
 {
 
-/*   CCTK_REAL diff[4]; */
-/*   CCTK_INT seed = 0; */
-/*   CCTK_INT j,r;   */
-
-/*   CCTK_REAL c0; */
-/*   CCTK_REAL c1; */
-/*   CCTK_REAL c2; */
-
-/*   CCTK_REAL xc; */
-/*   CCTK_REAL xr; */
-
-  CCTK_REAL result;
-  CCTK_REAL undiv_diff [13][6];
-  CCTK_INT diff, i, ii, r;
-  CCTK_REAL fp_i, fp_ii, x_rel;
-
-  /* Find seed index */
-/*   while (x > x0 + dx * ((CCTK_REAL)seed+0.5)) seed++; */
-
-/*   if (seed!=2) { */
-    /* Not enough stencil, only perform linear interpolation */   
-/*     seed = 0; */
-/*     while (x > x0 + dx * ((CCTK_REAL)(seed+1))) seed++; */
-/*     if (seed==4) seed=3; */
-/*     return y[seed] + (y[seed+1]-y[seed])/dx * (x-x0-(CCTK_REAL)seed*dx); */
-/*   } */
-
-  /*  Calculate the undivided differences */
-/*   for (j=0; j<=3; j++) */
-/*     diff[j] = y[j+1] - y[j]; */
-
-  /*  Find the stencil */
-/*   if ( fabs(diff[1]) < fabs(diff[2]) ) { */
-/*     if ( fabs(diff[1]-diff[0]) < fabs(diff[2]-diff[1]) ) */
-/*       r = 0; */
-/*     else */
-/*       r = 1; */
-/*   } else { */
-/*     if ( fabs(diff[2]-diff[1]) < fabs(diff[3]-diff[2]) ) */
-/*       r = 1; */
-/*     else */
-/*       r = 2; */
-/*   } */
-
-  /* Interpolate second order */
-/*   c0 = y[r+1]; */
-/*   c1 = (-1.0/2.0)*y[r] + (+1.0/2.0)*y[r+2]; */
-/*   c2 = (+1.0/2.0)*y[r] - y[r+1] + (+1.0/2.0)*y[r+2]; */
-
-/*   xc = x0 + dx * (CCTK_REAL)(r+1); */
-/*   xr = (x - xc) / dx; */
-/*   result = (  c0 + xr*(c1 + xr*c2)  ); */
-
-  for (i = 1; i < 2 * order + 2; ++i)
-  {
-    undiv_diff[i][0] = y[i-1];
-  }
-  for (i = 0; i < 6; ++i)
-  {
-    undiv_diff[0][i] = 1e10;
-    undiv_diff[2 * order + 2][i] = 1e10;
-  }
-
-  for (diff = 1; diff < order + 1; ++diff)
-  {
-    for (i = 1; i < 2 * order + 2; ++i)
-    {
-      undiv_diff[i][diff] = undiv_diff[i][diff-1] - undiv_diff[i-1][diff-1];
-    }
-    undiv_diff[0][diff] = -10 * undiv_diff[1][diff];
-    undiv_diff[2 * order + 2][diff]= -10 * undiv_diff[2 * order + 1][diff];
-  }
+  int start, end;
   
-  fp_i = (x - x0) / dx;
-  fp_ii = floor(fp_i);
-  x_rel = fp_i - fp_ii;
-  
-  ii = (CCTK_INT)(fp_ii) + 1;
-  
-  if (ii < 1)
-  {
-    ii = 1;
-    x_rel = fp_i - (CCTK_REAL)(ii) + 1.0;
-  }
-  else if (ii > 2 * order + 1)
-  {
-    ii = 2 * order + 1;
-    x_rel = fp_i - (CCTK_REAL)(ii) + 1.0;
-  }
+  start = 0;
+  end   = order;
   
-  r = 0;
-  for (diff = 1; diff < order + 1; ++diff)
-  {
-    if (fabs(undiv_diff[ii-r][diff]) < fabs(undiv_diff[ii-r+1][diff]) )
-    {
-      ++r;
-    }
-    
-    if (ii - r < diff + 1)
-    {
-      --r;
-    }
-    else if (ii - r + diff > 2 * order + 1)
-    {
-      ++r;
-    }
-  }
-
-  result = undiv_diff[ii-r][0];
+  /* The possible interpolation order that can be used */
 
-  switch (order)
-  {
-    case 1:
-      result += (x_rel + r - 0.0) / 1.0 * undiv_diff[ii-r+1][1];
-      break;
-    case 2:
-      result += (x_rel + r - 0.0) / 1.0 * (undiv_diff[ii-r+1][1] +
-                (x_rel + r - 1.0) / 2.0 * (undiv_diff[ii-r+2][2]));
-      break;
-    case 3:
-      result += (x_rel + r - 0.0) / 1.0 * (undiv_diff[ii-r+1][1] +
-                (x_rel + r - 1.0) / 2.0 * (undiv_diff[ii-r+2][2] +
-                (x_rel + r - 2.0) / 3.0 * (undiv_diff[ii-r+3][3])));
-      break;
-    case 4:
-      result += (x_rel + r - 0.0) / 1.0 * (undiv_diff[ii-r+1][1] +
-                (x_rel + r - 1.0) / 2.0 * (undiv_diff[ii-r+2][2] +
-                (x_rel + r - 2.0) / 3.0 * (undiv_diff[ii-r+3][3] +
-                (x_rel + r - 3.0) / 4.0 * (undiv_diff[ii-r+4][4]))));
-      break;
-    case 5:
-      result += (x_rel + r - 0.0) / 1.0 * (undiv_diff[ii-r+1][1] +
-                (x_rel + r - 1.0) / 2.0 * (undiv_diff[ii-r+2][2] +
-                (x_rel + r - 2.0) / 3.0 * (undiv_diff[ii-r+3][3] +
-                (x_rel + r - 3.0) / 4.0 * (undiv_diff[ii-r+4][4] +
-                (x_rel + r - 4.0) / 5.0 * (undiv_diff[ii-r+5][5])))));
-      break;
-  }
+  if (x0 > x)
+      return y[start];
   
-/* #define CARTOON_ENO_DBG 1 */
-#ifdef CARTOON_ENO_DBG
-
-  printf("Cartoon ENO interpolation.\n"
-         "Input: order %d x0 %g dx %g x %g.\n",
-         order, x0, dx, x);
-  printf("Data is\n");
-  for (i = 0; i < 2 * order + 1; ++i)
-  {
-    printf(" %g", y[i]);
+  if (dx * end + x0 < x)
+      return y[end];
+
+  switch (order) {
+  case 1: return interpolate_local_order1(x0, dx, y, x);
+  case 2: return interpolate_local_order2(x0, dx, y, x);
+  case 3: return interpolate_local_order3(x0, dx, y, x);
+  case 4: return interpolate_local_order4(x0, dx, y, x);
+  case 5: return interpolate_local_order5(x0, dx, y, x);
+  default: CCTK_WARN(0, "Interpolation order not supported");
+           return 0; /* Avoid warning */
   }
-  printf("\nEntries used were\n");
-  printf("0: %g   1: %g   2: %g\n",
-         undiv_diff[ii-r][0], undiv_diff[ii-r+1][1], undiv_diff[ii-r+2][2]);
-  printf("Parameters are fp_i %g fp_ii %g x_rel %g ii %d r %d\n",
-         fp_i, fp_ii, x_rel, ii, r);
-  printf("Result is %g\n", result);
-  
-#endif
-
-  return result;  
-
 }