diff options
Diffstat (limited to 'src/interpolate.c')
| -rw-r--r-- | src/interpolate.c | 186 |
1 files changed, 150 insertions, 36 deletions
diff --git a/src/interpolate.c b/src/interpolate.c index d990cf5..1d774a2 100644 --- a/src/interpolate.c +++ b/src/interpolate.c @@ -17,6 +17,7 @@ */ #include <math.h> +#include <stdio.h> #include "cctk.h" @@ -24,7 +25,7 @@ CCTK_REAL interpolate_local(int order, CCTK_REAL x0, CCTK_REAL dx, CCTK_REAL y[], CCTK_REAL x); -CCTK_REAL interpolate_eno(CCTK_REAL x0, CCTK_REAL dx, +CCTK_REAL interpolate_eno(CCTK_INT order, CCTK_REAL x0, CCTK_REAL dx, CCTK_REAL y[], CCTK_REAL x); @@ -214,58 +215,171 @@ return( c0 + xr*(c1 + xr*(c2 + xr*(c3 + xr*(c4 + xr*c5)))) ); /*****************************************************************************/ -CCTK_REAL interpolate_eno(CCTK_REAL x0, CCTK_REAL dx, +CCTK_REAL interpolate_eno(CCTK_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 diff[4]; */ +/* CCTK_INT seed = 0; */ +/* CCTK_INT j,r; */ - CCTK_REAL c0; - CCTK_REAL c1; - CCTK_REAL c2; +/* CCTK_REAL c0; */ +/* CCTK_REAL c1; */ +/* CCTK_REAL c2; */ - CCTK_REAL xc; - CCTK_REAL xr; +/* 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++; +/* while (x > x0 + dx * ((CCTK_REAL)seed+0.5)) seed++; */ - if (seed!=2) { +/* 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); - } +/* 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]; +/* 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; - } +/* 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]; +/* 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; +/* xc = x0 + dx * (CCTK_REAL)(r+1); */ +/* xr = (x - xc) / dx; */ +/* result = ( c0 + xr*(c1 + xr*c2) ); */ - return( 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]; + } + + 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; + } + + 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]; + + 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; + } + +/* #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]); + } + 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; } |