aboutsummaryrefslogtreecommitdiff
path: root/src/interpolate.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/interpolate.c')
-rw-r--r--src/interpolate.c454
1 files changed, 107 insertions, 347 deletions
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;
-
}
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-15 18:34:34 +0000