// MAUS WARNING: THIS IS LEGACY CODE. // @(#) $Id: Spline1D.hh,v 1.1 2005-10-28 08:22:42 rogers Exp $ $Name: $ // // ******************************************************************** // * DISCLAIMER * // * * // * Neither the authors of this software system, nor their employing * // * institutes,nor the agencies providing financial support for this * // * work make any representation or warranty, express or implied, * // * regarding this software system or assume any liability for its * // * use. * // * * // * This code implementation is the intellectual property of * // * FERMILAB. * // * By copying, distributing or modifying the Program (or any work * // * based on the Program) you indicate your acceptance of this * // * statement, and all its terms. * // ******************************************************************** // // // Spline1D.hh // // Created: Mark Fishler (5/00) // #ifndef SPLINE1D_H #define SPLINE1D_H // The Spline1D class represents a table-driven cubic spline in 1 dimension. // The user supplies the table or the function to the constructor; thereafter, // the spline is invoked as if it were a simple function of the input value x, // as in // // double answer = mySpline ( x ); // // The spline requires (or pre-computes) values of gradient at each point in // the table. For cases where the gradient is known analytically, these may be // supplied directly. Normally, the spline would be fed just the function or // values at the grid points; the gradients are pre-computed from the condition // that the function has a continuous second derivative. // A related class is Spline, for an N-dimensional cubic spline taking // an array instead of a value for x. The general Spline class is correct // for 1 dimension, but involves a bit of extra overhead, and the inconvenience // of requiring a pointer to the value x (rather than the value itself) because // it will treat x[] as an array of dimension 1. // See Spline.h. #ifdef Spline_DESIGN The Spline1D class is an object representing the transformation between some value x and a value phi, based on tables of values and gradients (derivatives). The values and gradients are computed (or supplied) at a set of points (we call these node points) upon construction. These node points may be irregularly spaced. From then on, the user can invoke the () operator supplying the value of x; the operator returns the value of the function phi, based on an interpolating function which exactly matches the values and slopes at every node point. A given instance of a Spline is constructed from: (1) An integer extent of the array defining the node points. (2) An array of N numbers representing the coordinates of each node points. (3) The values and optionally the gradients of phi at those node points. For convenience, the value and gradient information can be supplied in any of several ways: * A pointer to a function returning a double which takes a double x. If this is supplied, then the constructor will call the function however many times to fill in the tables of values at the node points. It will then compute the gradients needed for continuity in the second derivative, using "natural" boundary conditions (second derivative zero at both ends). * A single table containing the value for each point. Again, derivatives are comuted using "natural" boundary conditions. * Either a pointer to a function or a table, as above, followed by two additional arguments which represent the values of slope at the beginnng and end boundaries; * A table of values and a separate table of one gradients (derivative) per point. * A pointer to a function returning a void, which takes a double x and two pointers to value and gradient. If this is supplied, then the constructor will call the function however many times to fill in the tables of values and gradients at the node points. The tables are captured in the class instance data, and may be discarded after construction. The only user-relevant method other than the constructor is double operator() ( double x ) const; that is, supplying a x for the desired sample point and receiving a double value of phi. The value returned is a cubic function in each interval which will match the values and the gradients at every node. The class does NOT verify that x lies within the domain defined by the node points; if it is outside, then the value obtained is still a valid approximation, but one which will get worse as the distance from the nearest node grows large. #endif class Spline1D { // The usage of the spline, and any functions provided, are based on doubles. // However, for space purposes, the class might be reconfigured to store (and // accept) tables in the float format. typedef double Data_t; // In the future we MAY want to template off // the type of table data; using Data_t for // this purpose will make this easy. public: // Invoke the spline to get a value: double operator() ( double x ) const; // (Deep) copy constructor Spline1D ( const Spline1D & s ); Spline1D ( ); // Constructors A -- the class accepts value info in two ways, and optionally // accepts boundary slope info (if that is not provided, // then the second derivative is set to 0 at boundaries). #ifdef FUTURE Spline1D ( const std::vector nodes, // Locations in x of nodes const std::vector values // Values being approximated ); Spline1D ( const std::vector nodes, // Locations in x of nodes const std::vector values,// Values being approximated Data_t slope0, Data_t slopeN ); #endif // FUTURE Spline1D ( int extent, // Number of node points const Data_t* nodes, // Locations in x of nodes const Data_t* values // Values being approximated ); Spline1D ( int extent, // Number of node points const Data_t* nodes, // Locations in x of nodes double values_function (double x)// Function being approximated ); Spline1D ( int extent, // Number of node points const Data_t* nodes, // Locations in x of nodes const Data_t* values, // Values being approximated Data_t slope0, Data_t slopeN ); Spline1D ( int extent, // Number of node points const Data_t* nodes, // Locations in x of nodes double values_function (double x), // Function being approximated Data_t slope0, Data_t slopeN ); // Constructors B -- the class accepts value and gradient info in three ways: #ifdef FUTURE Spline1D ( const std::vector nodes, // Locations in x of nodes const std::vector values, // Values at each node point const std::vector gradients // derivatives at each node point ); // And some appropriate form for passing a values/grads function. #endif // FUTURE Spline1D ( int extent, // Number of node points const Data_t* nodes, // Locations in x of nodes const Data_t* values, // values at each node point const Data_t* gradients // derivatives at each node point // nodes, values, gradients are arrays // of extent Data_ts. ); Spline1D ( int extent, // Number of node points const Data_t* nodes, // Locations in x of nodes void values_grads_function ( double x, double* val, double* grads) // Function filling value and // gradient for input x ); ~Spline1D(); protected: // member data: int extent_; Data_t* nodePoints_; Data_t* distances_; // distances separating nodes Data_t* values_; Data_t* grads_; Data_t* secondDerivs_; bool naturalBoundaryConditions; Data_t slope0_; Data_t slopeN_; // internal methods: void captureGrid( int extent, const Data_t* nodes ); void captureValues( int extent, const Data_t* values ); void fillValues( int extent, double values_function(double x)); void computeSecondDerivs(); }; // Spline #endif // SPLINE1D_H