From 593f30e3475617fe87b90d2e79c52075a300af89 Mon Sep 17 00:00:00 2001 From: sai Date: Sat, 6 Nov 1999 16:19:25 +0000 Subject: This commit was generated by cvs2svn to compensate for changes in r2, which included commits to RCS files with non-trunk default branches. git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/Cartoon2D/trunk@3 eec4d7dc-71c2-46d6-addf-10296150bf52 --- README | 7 + interface.ccl | 2 + param.ccl | 22 + schedule.ccl | 2 + src/Cartoon2DBC.c | 49 +++ src/CheckParameters.c | 52 +++ src/SetBoundAxisym.c | 227 ++++++++++ src/error_exit.c | 63 +++ src/interpolate.c | 1149 +++++++++++++++++++++++++++++++++++++++++++++++++ src/make.code.defn | 9 + 10 files changed, 1582 insertions(+) create mode 100644 README create mode 100644 interface.ccl create mode 100644 param.ccl create mode 100644 schedule.ccl create mode 100644 src/Cartoon2DBC.c create mode 100644 src/CheckParameters.c create mode 100644 src/SetBoundAxisym.c create mode 100644 src/error_exit.c create mode 100644 src/interpolate.c create mode 100644 src/make.code.defn diff --git a/README b/README new file mode 100644 index 0000000..ff6a86f --- /dev/null +++ b/README @@ -0,0 +1,7 @@ +Cactus Code Thorn Cartoon2D +Authors : ... +CVS info : $Header$ +-------------------------------------------------------------------------- + +Purpose of the thorn: + diff --git a/interface.ccl b/interface.ccl new file mode 100644 index 0000000..52065c5 --- /dev/null +++ b/interface.ccl @@ -0,0 +1,2 @@ +# Interface definition for thorn Cartoon2D +# $Header$ diff --git a/param.ccl b/param.ccl new file mode 100644 index 0000000..52937e3 --- /dev/null +++ b/param.ccl @@ -0,0 +1,22 @@ +# Parameter definitions for thorn Cartoon2D +# $Header$ + +private: + +BOOLEAN active "Use cartoon to implement axisymmetry" +{ +} "no" + +BOOLEAN verbose "Verbose information" +{ +} "no" + + +BOOLEAN stencil "Use custom 2D stencil if available" +{ +} "yes" + +INT interp_order "Cartoon's interpolation order" +{ + 0:* :: +} 4 diff --git a/schedule.ccl b/schedule.ccl new file mode 100644 index 0000000..730af80 --- /dev/null +++ b/schedule.ccl @@ -0,0 +1,2 @@ +# Schedule definitions for thorn Cartoon2D +# $Header$ diff --git a/src/Cartoon2DBC.c b/src/Cartoon2DBC.c new file mode 100644 index 0000000..f1e1ebf --- /dev/null +++ b/src/Cartoon2DBC.c @@ -0,0 +1,49 @@ +/*@@ + @file Cartoon2DBC.c + @date Thu Nov 4 13:35:00 MET 1999 + @author Sai Iyer + @desc + Apply Cartoon2D boundary conditions + Adapted from Bernd Bruegmann's cartoon_2d_tensorbound.c + and set_bound_axisym.c + @enddesc + @@*/ + +static char *rcsid = "$Id$" + +#include "cctk.h" +#include "cctk_arguments.h" +#include "cctk_parameters.h" + +#include "cartoon2D.h" + +int Cartoon2DBCVarI(cGH *GH, int vi, int tensortype) { + DECLARE_CCTK_PARAMETERS + + static int firstcall = 1, pr = 0; + + if (firstcall) { + firstcall = 0; + pr = CCTK_Equals(verbose, "yes"); + } + if (pr) printf("cartoon_2d_tensorbound called for %s\n", name); + + return(SetBoundAxisym(GH, vi, tensortype)); +} + +void FMODIFIER FORTRAN_NAME(Cartoon2DBCVarI)(int *retval, cGH *GH, int *vi, int *tensortype) { + *retval = Cartoon2DBCVarI(GH, *vi, *tensortype); +} + +int Cartoon2DBCVar(cGH *GH, const char *impvarname, int tensortype) +{ + int vi; + vi = CCTK_VarIndex(impvarname); + return(Cartoon2DBCVarI(GH, vi, tensortype)); +} + +void FMODIFIER FORTRAN_NAME(Cartoon2DBCVar)(int *retval, cGH *GH, ONE_FORTSTRING_ARG, int *tensortype) { + ONE_FORTSTRING_CREATE(impvarname) + *retval = Cartoon2DBCVar(GH, impvarname, *tensortype); + free(impvarname); +} diff --git a/src/CheckParameters.c b/src/CheckParameters.c new file mode 100644 index 0000000..c7fcc79 --- /dev/null +++ b/src/CheckParameters.c @@ -0,0 +1,52 @@ +/*@@ + @file CheckParameters.c + @date Wed Nov 3 10:17:46 MET 1999 + @author Sai Iyer + @desc + Check Cartoon2D parameters + @enddesc + @@*/ + +#include "cctk.h" +#include "cctk_parameters.h" +#include "cctk_arguments.h" + +static char *rcsid = "$Id$"; + + /*@@ + @routine Cartoon2D_CheckParameters + @date Wed Nov 3 10:17:46 MET 1999 + @author Sai Iyer + @desc + Check Cartoon2D parameters + @enddesc + @calls + @calledby + @history + + @endhistory + +@@*/ + +void Cartoon2D_CheckParameters(CCTK_CARGUMENTS) +{ + DECLARE_CCTK_CARGUMENTS + DECLARE_CCTK_PARAMETERS + + if(cctk_gsh[1] != 2*cctk_nghostzones[1]+1) + { + CCTK_PARAMWARN("Grid size in y direction inappropriate."); + } + + if(order < 2*cctk_nghostzones[0]) + { + CCTK_PARAMWARN("Ghostzone width in x direction too small."); + } + + if(order < 1 || order > 5) + { + CCTK_PARAMWARN("Interpolation order should be between 1 and 5."); + } + return; +} + diff --git a/src/SetBoundAxisym.c b/src/SetBoundAxisym.c new file mode 100644 index 0000000..37096e5 --- /dev/null +++ b/src/SetBoundAxisym.c @@ -0,0 +1,227 @@ +/* set_bound_axisym.c, BB */ +/*@@ + @file SetBoundAxisym.c + @date Sat Nov 6 16:35:46 MET 1999 + @author Sai Iyer + @desc + An implementation of Steve Brandt's idea for doing + axisymmetry with 3d stencils. + Cactus 4 version of Bernd Bruegmann's set_bound_axisym_GT. + @enddesc + @@*/ + +static char *rcsid = "$Id$" + +#include "cctk.h" +#include "cctk_arguments.h" +#include "cctk_parameters.h" + +#include "cartoon2D.h" + +#include "math.h" + +/* Number of components */ +int tensor_type_n[N_TENSOR_TYPES] = {1, 3, 6}; + +/* set boundaries of a grid tensor assuming axisymmetry + - handles lower boundary in x + - does not touch other boundaries + - coordinates and rotation coefficients are independent of z and should + be precomputed + - this is also true for the constants in interpolator, but this may + be too messy + - minimizes conceptual warpage, not computationally optimized +*/ +/* uses rotation matrix and its inverse as linear transformation on + arbitrary tensor indices -- I consider this a good compromise + between doing index loops versus using explicit formulas in + cos(phi) etc with messy signs +*/ +int SetBoundAxiSym(cGH *GH, int vi, int tensortype) { + +{ + int i, j, k, di, dj, dk, ijk; + int jj, n, s; + Double x, y, rho; + Double rxx, rxy, ryx, ryy; + Double sxx, sxy, syx, syy; + Double f[6], fx, fy, fz, fxx, fxy, fxz, fyy, fyz, fzz; + Double *t, *tx, *ty, *tz, *txx, *txy, *txz, *tyy, *tyz, *tzz; + + if (0) printf("sba(%s)\n", GT->gf[0]->name); + + /* stencil width: in x-direction use ghost zone width, + in y-direction use ny/2, see below */ + s = GH->stencil_width; + + /* 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 = GH->lnx; + dk = GH->lnx*GH->lny; + + /* y = 0 */ + j = GH->ny/2; + + /* z-direction: include lower and upper boundary */ + for (k = 0; k < GH->lnz; k++) + + /* y-direction: as many zombies as the evolution stencil needs */ + for (jj = 1, dj = jj*GH->lnx; jj <= GH->ny/2; jj++, dj = jj*GH->lnx) + + /* x-direction: zombie for x < 0, including upper boundary for extrapol */ + for (i = -s; i < GH->lnx; i++) { + + /* index into linear memory */ + ijk = DI(GH,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 = GH->lx0 + GH->dx0 * i; + y = (dj) ? GH->dy0 * jj : 0; + 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 */ + for (n = 0; n < GT->n; n++) + f[n] = axisym_interpolate(GH, GT->gf[n], i, di, ijk, rho); + + /* rotate grid tensor by matrix multiplication */ + if (GT->tensor_type == TENSOR_TYPE_SCALAR) { + t = GT->gf[0]->data; + t[ijk+dj] = f[0]; + t[ijk-dj] = f[0]; + } + else if (GT->tensor_type == TENSOR_TYPE_U) { + tx = GT->gf[0]->data; + ty = GT->gf[1]->data; + tz = GT->gf[2]->data; + 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 (GT->tensor_type == TENSOR_TYPE_DDSYM) { + txx = GT->gf[0]->data; + txy = GT->gf[1]->data; + txz = GT->gf[2]->data; + tyy = GT->gf[3]->data; + tyz = GT->gf[4]->data; + tzz = GT->gf[5]->data; + 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 { + /* oops */ + } + + /* what a neat way to hack out the zombies for x < 0, y = 0 */ + if (dj == 0) {i -= s; ijk -= s; dj = jj*GH->lnx;} + } + + /* syncs needed after interpolation (actually only for x direction) */ + for (n = 0; n < GT->n; n++) + SyncGF(GH, GT->gf[n]); +} + + + + + +/* interpolation on x-axis */ +Double axisym_interpolate(pGH *GH, pGF *GF, int i, int di, int ijk, Double x) +{ + Double *f = GF->data; + Double c, d, x0; + + /* 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 > GH->lx0 + GH->dx0 * (i+1)) {i++; ijk++;} + + /* the original demo */ + if (0) { + c = (f[ijk+di] - f[ijk])/GH->dx0; + d = f[ijk]; + x0 = GH->lx0 + GH->dx0 * i; + return c*(x-x0) + d; + } + + /* first attempt to interface to JT's interpolator */ + if (1) { + double interpolate_local(int order, double x0, double dx, + double y[], double x); + double y[6], r; + int n, offset; + int order = cartoon_order; + + /* 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 >= GH->lnx) offset = i+order-GH->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 = GH->lx0 + GH->dx0 * (i-offset); + for (n = 0; n <= order; n++) { + y[n] = f[ijk-offset+n]; + if (0) printf("y[%d]=%6.4f ", n, y[n]); + } + + /* call interpolator */ + r = interpolate_local(order, x0, GH->dx0, y, x); + + if (0) printf("i %d x %7.4f x0 %7.4f r %7.4f\n", i, x, x0, r); + return r; + } + + return 0; +} diff --git a/src/error_exit.c b/src/error_exit.c new file mode 100644 index 0000000..748de99 --- /dev/null +++ b/src/error_exit.c @@ -0,0 +1,63 @@ +/* error_exit.c -- print an error message and cleanly shutdown cactus */ +/* $Id$ */ + +#include +#include +#include + +#include "cactus.h" +#include "thorn_cartoon_2d/src/error_exit.h" + +/* + * Jonathan Thornburg's favorite syntactic sugar to make the two branches + * of an if-else statement symmetrical: + * + * Also note (in)famous joke saying + * "syntactic sugar causes cancer of the semicolon" + * :) :) :) + */ +#define then /* empty */ + +/*****************************************************************************/ + +/* + * This function prints an error message (formatted via vfprintf(3S)), + * then cleanly shuts down cactus. It does *not* return to its caller. + * + * Despite not actually returning, this function is declared as returning + * an int , so it may be easily used in conditional expressions like + * foo = bar_ok ? baz(bar) : error_exit(...); + * + * Usage: + * error_exit(exit_code, message, args...) + * where exit_code is one of the codes defined in "error_exit.h" + * + * ***** THIS VERSION IS FOR ANSI/ISO C ***** + * + * Arguments: + * exit_code = (in) One of the codes defined in "error_exit.h", or any + * other integer to pass back to the shell as our exit + * status. + * format = (in) vprintf(3S) format string for error message to be printed. + * args... = (in) Any additional arguments are (presumably) used in formatting + * the error message string. + */ +#ifdef __cplusplus +extern "C" +#endif +/*VARARGS*/ +/* PUBLIC */ +int error_exit(const int exit_code, const char *const format, ...) +{ +va_list ap; + +va_start(ap, format); +vfprintf(stderr, format, ap); +va_end(ap); + + if (exit_code == 0) { + then STOP; /*NOTREACHED*/ + } else if (exit_code > 0) { + then ERROR_STOP(exit_code); /*NOTREACHED*/ + } else ERROR_ABORT(exit_code); /*NOTREACHED*/ +} diff --git a/src/interpolate.c b/src/interpolate.c new file mode 100644 index 0000000..346a1fa --- /dev/null +++ b/src/interpolate.c @@ -0,0 +1,1149 @@ +/* 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 + */ + + + +/* 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 + +/* 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); + +/*****************************************************************************/ +/*****************************************************************************/ +/*****************************************************************************/ + +/* + * ***** 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) + ); +} + +/*****************************************************************************/ + +/* + * 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. + */ +PUBLIC double interpolate_local(int order, + double x0, double dx, + double y[], + double x) +{ +switch (order) + { +case 1: + return( interpolate_local_order1(x0, dx, y, x) ); + break; + +case 2: + return( interpolate_local_order2(x0, dx, y, x) ); + break; + +case 3: + return( interpolate_local_order3(x0, dx, y, x) ); + break; + +case 4: + return( interpolate_local_order4(x0, dx, y, x) ); + break; + +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*/ + } +} + +/*****************************************************************************/ + +/* + * 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) +{ +double c0 = (+1.0/2.0)*y[0] + (+1.0/2.0)*y[1]; +double c1 = - y[0] + + y[1]; + +double xc = x0 + 0.5*dx; +double xr = (x - xc) / dx; + +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) +{ +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]; + +double xc = x0 + 1.0*dx; +double xr = (x - xc) / dx; + +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) +{ +double 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] + + (+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] + + (-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] + + (-1.0/ 2.0)*y[2] + ( 1.0/ 6.0)*y[3]; + +double xc = x0 + 1.5*dx; +double xr = (x - xc) / dx; + +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) +{ +double c0 = y[2]; +double 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] + + (+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] + + (-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] + + (-1.0/ 6.0)*y[3] + (+1.0/24.0)*y[4]; + +double xc = x0 + 2.0*dx; +double xr = (x - xc) / dx; + +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) +{ +double 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] + + (-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] + + (-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] + + (+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] + + (+ 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] + + (- 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; + +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); +} diff --git a/src/make.code.defn b/src/make.code.defn new file mode 100644 index 0000000..0483ce6 --- /dev/null +++ b/src/make.code.defn @@ -0,0 +1,9 @@ +# Main make.code.defn file for thorn Cartoon2D +# $Header$ + +# Source files in this directory +SRCS = + +# Subdirectories containing source files +SUBDIRS = + -- cgit v1.2.3