From 8fe9c3a66c1dc44e137e176ff5018f4d8d65eb4d Mon Sep 17 00:00:00 2001 From: sai Date: Fri, 19 Nov 1999 10:25:13 +0000 Subject: Recreating code after origin /data crash. git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/Cartoon2D/trunk@6 eec4d7dc-71c2-46d6-addf-10296150bf52 --- src/interpolate.c | 1072 ++++------------------------------------------------- 1 file changed, 73 insertions(+), 999 deletions(-) (limited to 'src/interpolate.c') 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 */ /* -** <<>> - - * 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 -#include -#include -#include -#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); -} -- cgit v1.2.3