1 files changed, 55 insertions, 0 deletions
diff --git a/src/interpolate.c b/src/interpolate.c
index 8a39c44..d48da6e 100644
--- a/src/interpolate.c
+++ b/src/interpolate.c
@@ -16,6 +16,7 @@
 */
 
 #include "cctk.h"
+#include "math.h"
 
 #define PRIVATE 
 #define PUBLIC 
@@ -40,6 +41,9 @@ PRIVATE CCTK_REAL interpolate_local_order5(CCTK_REAL x0, CCTK_REAL dx,
 					CCTK_REAL y[],
 					CCTK_REAL x);
 
+PRIVATE CCTK_REAL interpolate_eno(CCTK_REAL x0, CCTK_REAL dx,
+                                  CCTK_REAL y[],
+                                  CCTK_REAL x);
 
 /*****************************************************************************/
 
@@ -225,3 +229,54 @@ CCTK_REAL xr = (x - xc) / dx;
 return(  c0 + xr*(c1 + xr*(c2 + xr*(c3 + xr*(c4 + xr*c5))))  );
 }
 
+
+/*****************************************************************************/
+
+CCTK_REAL interpolate_eno(CCTK_REAL x0, CCTK_REAL dx,
+					CCTK_REAL y[],
+					CCTK_REAL x)
+{
+
+  CCTK_REAL diff[4];
+  CCTK_INT seed = 0;
+  CCTK_INT j,r;  
+
+  /* Find seed index */
+  while (x > x0 + dx * ((float)seed+0.5)) seed++;
+
+  if (seed!=2) {
+    /* Not enough stencil, only perform linear interpolation */   
+    seed = 0;
+    while (x > x0 + dx * ((float)(seed+1))) seed++;
+    if (seed==4) seed=3;
+    return y[seed] + (y[seed+1]-y[seed])/dx * (x-x0-(float)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 */
+  double c0 = y[r+1];
+  double c1 = (-1.0/2.0)*y[r] + (+1.0/2.0)*y[r+2];
+  double c2 = (+1.0/2.0)*y[r] - y[r+1] + (+1.0/2.0)*y[r+2];
+
+  double xc = x0 + dx * (double)(r+1);
+  double xr = (x - xc) / dx;
+
+  return(  c0 + xr*(c1 + xr*c2)  );  
+
+}