sandbox/magfield/src
--- sandbox/magfield/src/spline.cpp 2014-04-04 19:02:38 UTC (rev 3084)
+++ sandbox/magfield/src/spline.cpp 2014-04-04 21:13:31 UTC (rev 3085)
@@ -1,987 +1,6827 @@
-// Magic Software, Inc.
-// http://www.magic-software.com
-// Copyright (c) 2000, 2001. All Rights Reserved
-//
-// Source code from Magic Software is supplied under the terms of a license
-// agreement and may not be copied or disclosed except in accordance with the
-// terms of that agreement. The various license agreements may be found at
-// the Magic Software web site. This file is subject to the license
-//
-// FREE SOURCE CODE
-// http://www.magic-software.com/License/free.pdf
-
-#include <fstream>
-#include <math.h>
-#include <float.h>
-#include "spline.h"
-//===========================================================================
-const double mgcSpline::INVALID_DATA = FLT_MAX;
-
-// error handling
-int mgcSpline::verbose = 0;
-unsigned mgcSpline::error = 0;
-const unsigned mgcSpline::invalid_dimensions = 0x00000001;
-const unsigned mgcSpline::invalid_degree = 0x00000002;
-const unsigned mgcSpline::invalid_dim = 0x00000004;
-const unsigned mgcSpline::null_data = 0x00000008;
-const unsigned mgcSpline::allocation_failed = 0x00000010;
-const char* mgcSpline::message[5] = {
- "invalid number of dimensions",
- "degree must be positive",
- "dimension must be at least DEGREE+1",
- "input data pointer is null",
- "failed to allocate memory"
-};
-//---------------------------------------------------------------------------
-void mgcSpline::
-InitNull ()
-{
- dimensions = 0;
- degree = 0;
- Dp1 = 0;
- Dp1_to_N = 0;
- Dp1_to_2N = 0;
- dim = 0;
- data = 0;
- dom_min = 0;
- dom_max = 0;
- grid_min = 0;
- grid_max = 0;
- base = 0;
- old_base = 0;
- cache = 0;
- inter = 0;
- poly = 0;
- coeff = 0;
- matrix = 0;
- product = 0;
- skip = 0;
-}
-//---------------------------------------------------------------------------
-void mgcSpline::
-FreeStorage ()
-{
- if ( dim == 0 ) // nothing to free
- return;
-
- delete[] dim;
- delete[] base;
- delete[] old_base;
- delete[] dom_min;
- delete[] dom_max;
- delete[] grid_min;
- delete[] grid_max;
- delete[] cache;
- delete[] inter;
- delete[] product;
- delete[] skip;
-
- int d;
- for (d = 0; d < dimensions; d++)
- delete[] poly[d];
- delete[] poly;
-
- for (d = 0; d <= degree; d++) {
- delete[] matrix[d];
- delete[] coeff[d];
- }
- delete[] matrix;
- delete[] coeff;
-}
-//---------------------------------------------------------------------------
-void mgcSpline::
-Create ()
-{
- int k, d;
-
- // setup degree constants
- Dp1 = degree+1;
- Dp1_to_N = 1;
- for (d = 0; d < dimensions; d++)
- Dp1_to_N *= Dp1;
- Dp1_to_2N = Dp1_to_N*Dp1_to_N;
-
- // initialize base indices
- base = new int[dimensions];
- if ( base == 0 ) {
- Report(allocation_failed);
- return;
- }
- old_base = new int[dimensions];
- if ( old_base == 0 ) {
- Report(allocation_failed);
- return;
- }
- for (d = 0; d < dimensions; d++)
- old_base[d] = -1;
-
- grid_min = new int[dimensions];
- if ( grid_min == 0 ) {
- Report(allocation_failed);
- return;
- }
- grid_max = new int[dimensions];
- if ( grid_max == 0 ) {
- Report(allocation_failed);
- return;
- }
- for (d = 0; d < dimensions; d++) {
- grid_min[d] = -1;
- grid_max[d] = -1;
- }
-
- // compute domain [min,max] for B-spline
- dom_min = new double[dimensions];
- if ( dom_min == 0 ) {
- Report(allocation_failed);
- return;
- }
-
- dom_max = new double[dimensions];
- if ( dom_max == 0 ) {
- Report(allocation_failed);
- return;
- }
-
- for (d = 0; d < dimensions; d++) {
- double tmpmin = 1.0f, tmpmax, dom_sup = double(dim[d]-degree+1);
- double next = (tmpmin+dom_sup)/2;
- do {
- tmpmax = next;
- next = (next+dom_sup)/2;
- } while ( next < dom_sup );
-
- dom_min[d] = tmpmin;
- dom_max[d] = tmpmax;
- }
-
- // cache for optimizing compute_intermediate()
- cache = new double[Dp1_to_N];
- if ( cache == 0 ) {
- Report(allocation_failed);
- return;
- }
-
- // storage for intermediate tensor product
- inter = new double[Dp1_to_N];
- if ( inter == 0 ) {
- Report(allocation_failed);
- return;
- }
-
- // polynomial allocations
- poly = new double*[dimensions];
- if ( poly == 0 ) {
- Report(allocation_failed);
- return;
- }
- for (d = 0; d < dimensions; d++) {
- poly[d] = new double[Dp1];
- if ( poly[d] == 0 ) {
- Report(allocation_failed);
- return;
- }
- }
-
- // coefficients for polynomial calculations
- coeff = new double*[Dp1];
- if ( coeff == 0 ) {
- Report(allocation_failed);
- return;
- }
- for (int row = 0; row <= degree; row++) {
- coeff[row] = new double[Dp1];
- if ( coeff[row] == 0 ) {
- Report(allocation_failed);
- return;
- }
- for (int col = row; col <= degree; col++) {
- coeff[row][col] = 1;
- for (d = 0; d <= row-1; d++)
- coeff[row][col] *= col-d;
- }
- }
-
- // generate spline blending matrix
- matrix = BlendMatrix(degree);
-
- // tensor product of m with itself N times
- product = new double[Dp1_to_2N];
- if ( product == 0 ) {
- Report(allocation_failed);
- return;
- }
- skip = new int[Dp1_to_2N];
- if ( skip == 0 ) {
- Report(allocation_failed);
- return;
- }
- int* coord = new int[2*dimensions]; // for address decoding
- if ( coord == 0 ) {
- Report(allocation_failed);
- return;
- }
- for (k = 0; k < Dp1_to_2N; k++) {
- int temp = k;
- for (d = 0; d < 2*dimensions; d++) {
- coord[d] = temp % Dp1;
- temp /= Dp1;
- }
- product[k] = 1;
- for (int s = 0; s < dimensions; s++)
- product[k] *= matrix[coord[s]][coord[s+dimensions]];
-
- skip[k] = 1;
- }
- delete[] coord;
-
- // compute increments to skip zero values of mtensor
- for (k = 0; k < Dp1_to_2N; ) {
- for (d = k+1; d < Dp1_to_2N && product[d] == 0; d++)
- skip[k]++;
- k = d;
- }
-
- evaluate_callback = 0;
-}
-//---------------------------------------------------------------------------
-mgcSpline& mgcSpline::
-Renew (int _dimensions, int _degree, const int* _dim, double* _data)
-{
- int d;
-
- // blow away old data
- FreeStorage();
-
- // get number of dimensions
- if ( (dimensions = _dimensions) <= 0 ) {
- Report(invalid_dimensions);
- return *this;
- }
-
- // get degree and check if valid
- if ( (degree = _degree) <= 0 ) {
- Report(invalid_degree);
- return *this;
- }
-
- // get dimension bounds and check if valid
- if ( _dim == 0 ) {
- Report(invalid_dim);
- return *this;
- }
-
- for (d = 0; d < dimensions; d++)
- if ( _dim[d] <= degree+1 ) {
- Report(invalid_dim);
- return *this;
- }
-
- dim = new int[dimensions];
- if ( dim == 0 ) {
- Report(allocation_failed);
- return *this;
- }
- for (d = 0; d < dimensions; d++)
- dim[d] = _dim[d];
-
- // get spline data
- if ( (data = _data) == 0 ) {
- Report(null_data);
- return *this;
- }
-
- Create();
-
- return *this;
-}
-//---------------------------------------------------------------------------
-mgcSpline::
-mgcSpline (int _dimensions, int _degree, const int* _dim, double* _data)
-{
- InitNull();
- Renew(_dimensions,_degree,_dim,_data);
-}
-//---------------------------------------------------------------------------
-void mgcSpline::
-SetPolynomial (int order, double diff, double* poly)
-{
- double diffpower = 1;
- for (int d = order; d <= degree; d++) {
- poly[d] = coeff[order][d]*diffpower;
- diffpower *= diff;
- }
-}
-//---------------------------------------------------------------------------
-int mgcSpline::
-Choose (int n, int k)
-{
- int i;
-
- // computes combination "n choose k"
-
- if ( n <= 1 || k >= n )
- return 1;
-
- int result = 1;
- for (i = 0; i < k; i++)
- result *= n-i;
- for (i = 1; i <= k; i++)
- result /= i;
-
- return result;
-}
-//---------------------------------------------------------------------------
-double** mgcSpline::
-BlendMatrix (int deg)
-{
- int degp1 = deg+1;
- int row, col, i, j, d;
-
- // allocate triple arrays
- int*** a = new int**[degp1];
- if ( a == 0 ) {
- Report(allocation_failed);
- return 0;
- }
- int*** b = new int**[degp1];
- if ( b == 0 ) {
- Report(allocation_failed);
- return 0;
- }
- for (d = 0; d <= deg; d++) {
- a[d] = new int*[degp1];
- if ( a[d] == 0 ) {
- Report(allocation_failed);
- return 0;
- }
- b[d] = new int*[degp1];
- if ( b[d] == 0 ) {
- Report(allocation_failed);
- return 0;
- }
- for (row = 0; row <= deg; row++) {
- a[d][row] = new int[degp1];
- if ( a[d][row] == 0 ) {
- Report(allocation_failed);
- return 0;
- }
- b[d][row] = new int[degp1];
- if ( b[d][row] == 0 ) {
- Report(allocation_failed);
- return 0;
- }
- for (col = 0; col <= deg; col++) {
- a[d][row][col] = 0;
- b[d][row][col] = 0;
- }
- }
- }
-
- a[0][0][0] = 1;
- b[0][0][0] = 1;
-
- for (d = 1; d <= deg; d++) {
- // compute a[]
- for (row = 0; row <= d; row++) {
- for (col = 0; col <= d; col++) {
- a[d][row][col] = 0;
- if ( col >= 1 ) {
- a[d][row][col] += a[d-1][row][col-1];
- if ( row >= 1 )
- a[d][row][col] -= b[d-1][row-1][col-1];
- }
- if ( row >= 1 )
- a[d][row][col] += (d+1)*b[d-1][row-1][col];
- }
- }
-
- // compute b[]
- for (row = 0; row <= d; row++) {
- for (col = 0; col <= d; col++) {
- b[d][row][col]= 0;
- for (i = col; i <= d; i++)
- if ( (i-col) % 2 )
- b[d][row][col] -= Choose(i,col)*a[d][row][i];
- else
- b[d][row][col] += Choose(i,col)*a[d][row][i];
- }
- }
- }
-
- double** imat = new double*[degp1];
- if ( imat == 0 ) {
- Report(allocation_failed);
- return 0;
- }
- for (row = 0; row <= deg; row++) {
- imat[row] = new double[degp1];
- if ( imat[row] == 0 ) {
- Report(allocation_failed);
- return 0;
- }
- for (col = 0; col <= deg; col++) {
- imat[row][col]= 0;
- for (i = col; i <= deg; i++) {
- int prod = 1;
- for (j = 1; j <= i-col; j++)
- prod *= deg-row;
- imat[row][col] += prod*Choose(i,col)*a[deg][deg-row][i];
- }
- }
- }
-
- double factorial = 1;
- for (d = 1; d <= deg; d++)
- factorial *= d;
- double** mat = new double*[degp1];
- if ( mat == 0 ) {
- Report(allocation_failed);
- return 0;
- }
- for (row = 0; row <= deg; row++) {
- mat[row] = new double[degp1];
- if ( mat[row] == 0 ) {
- Report(allocation_failed);
- return 0;
- }
- for (col = 0; col <= deg; col++)
- mat[row][col] = imat[row][col]/factorial;
- }
-
- // deallocate triple arrays
- for (d = 0; d <= deg; d++) {
- for (row = 0; row <= deg; row++) {
- delete[] b[d][row];
- delete[] a[d][row];
- }
- delete[] b[d];
- delete[] a[d];
- }
- delete[] b;
- delete[] a;
-
- // deallocate integer matrix
- for (d = 0; d <= deg; d++)
- delete[] imat[d];
- delete[] imat;
-
- return mat;
-}
-//---------------------------------------------------------------------------
-int mgcSpline::
-Number (unsigned single_error)
-{
- int result;
- for (result = -1; single_error; single_error >>= 1)
- result++;
- return result;
-}
-//---------------------------------------------------------------------------
-void mgcSpline::
-Report (unsigned single_error)
-{
- if ( mgcSpline::verbose )
- std::cout << "mgcSpline: " << message[Number(single_error)] << std::endl;
- else
- std::ofstream("spline.err",std::ios::out|std::ios::app)
- << "mgcSpline: " << message[Number(single_error)] << std::endl;
-
- error |= single_error;
-}
-//---------------------------------------------------------------------------
-void mgcSpline::
-Report (std::ostream& ostr)
-{
- for (unsigned single_error = 1; single_error; single_error <<= 1)
- if ( error & single_error )
- ostr << "mgcSpline: " << message[Number(single_error)] << std::endl;
-
- error = 0;
-}
-//===========================================================================
-const double mgcSplineD::INVALID_DATA = DBL_MAX;
-
-// error handling
-int mgcSplineD::verbose = 0;
-unsigned mgcSplineD::error = 0;
-const unsigned mgcSplineD::invalid_dimensions = 0x00000001;
-const unsigned mgcSplineD::invalid_degree = 0x00000002;
-const unsigned mgcSplineD::invalid_dim = 0x00000004;
-const unsigned mgcSplineD::null_data = 0x00000008;
-const unsigned mgcSplineD::allocation_failed = 0x00000010;
-const char* mgcSplineD::message[5] = {
- "invalid number of dimensions",
- "degree must be positive",
- "dimension must be at least DEGREE+1",
- "input data pointer is null",
- "failed to allocate memory"
-};
-//---------------------------------------------------------------------------
-void mgcSplineD::
-InitNull ()
-{
- dimensions = 0;
- degree = 0;
- Dp1 = 0;
- Dp1_to_N = 0;
- Dp1_to_2N = 0;
- dim = 0;
- data = 0;
- dom_min = 0;
- dom_max = 0;
- grid_min = 0;
- grid_max = 0;
- base = 0;
- old_base = 0;
- cache = 0;
- inter = 0;
- poly = 0;
- coeff = 0;
- matrix = 0;
- product = 0;
- skip = 0;
-}
-//---------------------------------------------------------------------------
-void mgcSplineD::
-FreeStorage ()
-{
- if ( dim == 0 ) // no storage to deallocate
- return;
-
- delete[] dim;
- delete[] base;
- delete[] old_base;
- delete[] dom_min;
- delete[] dom_max;
- delete[] grid_min;
- delete[] grid_max;
- delete[] cache;
- delete[] inter;
- delete[] product;
- delete[] skip;
-
- int d;
- for (d = 0; d < dimensions; d++)
- delete[] poly[d];
- delete[] poly;
-
- for (d = 0; d <= degree; d++) {
- delete[] matrix[d];
- delete[] coeff[d];
- }
- delete[] matrix;
- delete[] coeff;
-}
-//---------------------------------------------------------------------------
-void mgcSplineD::
-Create ()
-{
- int k, d;
-
- // setup degree constants
- Dp1 = degree+1;
- Dp1_to_N = 1;
- for (d = 0; d < dimensions; d++)
- Dp1_to_N *= Dp1;
- Dp1_to_2N = Dp1_to_N*Dp1_to_N;
-
- // initialize base indices
- base = new int[dimensions];
- if ( base == 0 ) {
- Report(allocation_failed);
- return;
- }
- old_base = new int[dimensions];
- if ( old_base == 0 ) {
- Report(allocation_failed);
- return;
- }
- for (d = 0; d < dimensions; d++)
- old_base[d] = -1;
-
- grid_min = new int[dimensions];
- if ( grid_min == 0 ) {
- Report(allocation_failed);
- return;
- }
- grid_max = new int[dimensions];
- if ( grid_max == 0 ) {
- Report(allocation_failed);
- return;
- }
- for (d = 0; d < dimensions; d++) {
- grid_min[d] = -1;
- grid_max[d] = -1;
- }
-
- // compute domain [min,max] for B-spline
- dom_min = new double[dimensions];
- if ( dom_min == 0 ) {
- Report(allocation_failed);
- return;
- }
-
- dom_max = new double[dimensions];
- if ( dom_max == 0 ) {
- Report(allocation_failed);
- return;
- }
-
- for (d = 0; d < dimensions; d++) {
- double tmpmin = 1, tmpmax, dom_sup = dim[d]-degree+1;
- double next = (tmpmin+dom_sup)/2;
- do {
- tmpmax = next;
- next = (next+dom_sup)/2;
- } while ( next < dom_sup );
-
- dom_min[d] = tmpmin;
- dom_max[d] = tmpmax;
- }
-
- // cache for optimizing compute_intermediate()
- cache = new double[Dp1_to_N];
- if ( cache == 0 ) {
- Report(allocation_failed);
- return;
- }
-
- // storage for intermediate tensor product
- inter = new double[Dp1_to_N];
- if ( inter == 0 ) {
- Report(allocation_failed);
- return;
- }
-
- // polynomial allocations
- poly = new double*[dimensions];
- if ( poly == 0 ) {
- Report(allocation_failed);
- return;
- }
- for (d = 0; d < dimensions; d++) {
- poly[d] = new double[Dp1];
- if ( poly[d] == 0 ) {
- Report(allocation_failed);
- return;
- }
- }
-
- // coefficients for polynomial calculations
- coeff = new double*[Dp1];
- if ( coeff == 0 ) {
- Report(allocation_failed);
- return;
- }
- for (int row = 0; row <= degree; row++) {
- coeff[row] = new double[Dp1];
- if ( coeff[row] == 0 ) {
- Report(allocation_failed);
- return;
- }
- for (int col = row; col <= degree; col++) {
- coeff[row][col] = 1;
- for (d = 0; d <= row-1; d++)
- coeff[row][col] *= col-d;
- }
- }
-
- // generate spline blending matrix
- matrix = BlendMatrix(degree);
-
- // tensor product of m with itself N times
- product = new double[Dp1_to_2N];
- if ( product == 0 ) {
- Report(allocation_failed);
- return;
- }
- skip = new int[Dp1_to_2N];
- if ( skip == 0 ) {
- Report(allocation_failed);
- return;
- }
- int* coord = new int[2*dimensions]; // for address decoding
- if ( coord == 0 ) {
- Report(allocation_failed);
- return;
- }
- for (k = 0; k < Dp1_to_2N; k++) {
- int temp = k;
- for (d = 0; d < 2*dimensions; d++) {
- coord[d] = temp % Dp1;
- temp /= Dp1;
- }
- product[k] = 1;
- for (int s = 0; s < dimensions; s++)
- product[k] *= matrix[coord[s]][coord[s+dimensions]];
-
- skip[k] = 1;
- }
- delete[] coord;
-
- // compute increments to skip zero values of mtensor
- for (k = 0; k < Dp1_to_2N; ) {
- for (d = k+1; d < Dp1_to_2N && product[d] == 0; d++)
- skip[k]++;
- k = d;
- }
-
- evaluate_callback = 0;
-}
-//---------------------------------------------------------------------------
-mgcSplineD& mgcSplineD::
-Renew (int _dimensions, int _degree, const int* _dim, double* _data)
-{
- // blow away old data
- FreeStorage();
-
- // get number of dimensions
- if ( (dimensions = _dimensions) <= 0 ) {
- Report(invalid_dimensions);
- return *this;
- }
-
- // get degree and check if valid
- if ( (degree = _degree) <= 0 ) {
- Report(invalid_degree);
- return *this;
- }
-
- // get dimension bounds and check if valid
- if ( _dim == 0 ) {
- Report(invalid_dim);
- return *this;
- }
-
- int d;
- for (d = 0; d < dimensions; d++)
- if ( _dim[d] <= degree+1 ) {
- Report(invalid_dim);
- return *this;
- }
-
- dim = new int[dimensions];
- if ( dim == 0 ) {
- Report(allocation_failed);
- return *this;
- }
- for (d = 0; d < dimensions; d++)
- dim[d] = _dim[d];
-
- // get spline data
- if ( (data = _data) == 0 ) {
- Report(null_data);
- return *this;
- }
-
- Create();
-
- return *this;
-}
-//---------------------------------------------------------------------------
-mgcSplineD::
-mgcSplineD (int _dimensions, int _degree, const int* _dim, double* _data)
-{
- InitNull();
- Renew(_dimensions,_degree,_dim,_data);
-}
-//---------------------------------------------------------------------------
-void mgcSplineD::
-SetPolynomial (int order, double diff, double* poly)
-{
- double diffpower = 1;
- for (int d = order; d <= degree; d++) {
- poly[d] = coeff[order][d]*diffpower;
- diffpower *= diff;
- }
-}
-//---------------------------------------------------------------------------
-int mgcSplineD::
-Choose (int n, int k)
-{
- int i;
-
- // computes combination "n choose k"
-
- if ( n <= 1 || k >= n )
- return 1;
-
- int result = 1;
- for (i = 0; i < k; i++)
- result *= n-i;
- for (i = 1; i <= k; i++)
- result /= i;
-
- return result;
-}
-//---------------------------------------------------------------------------
-double** mgcSplineD::
-BlendMatrix (int deg)
-{
- int degp1 = deg+1;
- int row, col, i, j, d;
-
- // allocate triple arrays
- int*** a = new int**[degp1];
- if ( a == 0 ) {
- Report(allocation_failed);
- return 0;
- }
- int*** b = new int**[degp1];
- if ( b == 0 ) {
- Report(allocation_failed);
- return 0;
- }
- for (d = 0; d <= deg; d++) {
- a[d] = new int*[degp1];
- if ( a[d] == 0 ) {
- Report(allocation_failed);
- return 0;
- }
- b[d] = new int*[degp1];
- if ( b[d] == 0 ) {
- Report(allocation_failed);
- return 0;
- }
- for (row = 0; row <= deg; row++) {
- a[d][row] = new int[degp1];
- if ( a[d][row] == 0 ) {
- Report(allocation_failed);
- return 0;
- }
- b[d][row] = new int[degp1];
- if ( b[d][row] == 0 ) {
- Report(allocation_failed);
- return 0;
- }
- for (col = 0; col <= deg; col++) {
- a[d][row][col] = 0;
- b[d][row][col] = 0;
- }
- }
- }
-
- a[0][0][0] = 1;
- b[0][0][0] = 1;
-
- for (d = 1; d <= deg; d++) {
- // compute a[]
- for (row = 0; row <= d; row++) {
- for (col = 0; col <= d; col++) {
- a[d][row][col] = 0;
- if ( col >= 1 ) {
- a[d][row][col] += a[d-1][row][col-1];
- if ( row >= 1 )
- a[d][row][col] -= b[d-1][row-1][col-1];
- }
- if ( row >= 1 )
- a[d][row][col] += (d+1)*b[d-1][row-1][col];
- }
- }
-
- // compute b[]
- for (row = 0; row <= d; row++) {
- for (col = 0; col <= d; col++) {
- b[d][row][col]= 0;
- for (i = col; i <= d; i++)
- if ( (i-col) % 2 )
- b[d][row][col] -= Choose(i,col)*a[d][row][i];
- else
- b[d][row][col] += Choose(i,col)*a[d][row][i];
- }
- }
- }
-
- double** imat = new double*[degp1];
- if ( imat == 0 ) {
- Report(allocation_failed);
- return 0;
- }
- for (row = 0; row <= deg; row++) {
- imat[row] = new double[degp1];
- if ( imat[row] == 0 ) {
- Report(allocation_failed);
- return 0;
- }
- for (col = 0; col <= deg; col++) {
- imat[row][col]= 0;
- for (i = col; i <= deg; i++) {
- int prod = 1;
- for (j = 1; j <= i-col; j++)
- prod *= deg-row;
- imat[row][col] += prod*Choose(i,col)*a[deg][deg-row][i];
- }
- }
- }
-
- double factorial = 1;
- for (d = 1; d <= deg; d++)
- factorial *= d;
- double** mat = new double*[degp1];
- if ( mat == 0 ) {
- Report(allocation_failed);
- return 0;
- }
- for (row = 0; row <= deg; row++) {
- mat[row] = new double[degp1];
- if ( mat[row] == 0 ) {
- Report(allocation_failed);
- return 0;
- }
- for (col = 0; col <= deg; col++)
- mat[row][col] = imat[row][col]/factorial;
- }
-
- // deallocate triple arrays
- for (d = 0; d <= deg; d++) {
- for (row = 0; row <= deg; row++) {
- delete[] b[d][row];
- delete[] a[d][row];
- }
- delete[] b[d];
- delete[] a[d];
- }
- delete[] b;
- delete[] a;
-
- // deallocate integer matrix
- for (d = 0; d <= deg; d++)
- delete[] imat[d];
- delete[] imat;
-
- return mat;
-}
-//---------------------------------------------------------------------------
-int mgcSplineD::
-Number (unsigned single_error)
-{
- int result;
- for (result = -1; single_error; single_error >>= 1)
- result++;
- return result;
-}
-//---------------------------------------------------------------------------
-void mgcSplineD::
-Report (unsigned single_error)
-{
- if ( mgcSplineD::verbose )
- std::cout << "mgcSplineD: " << message[Number(single_error)] << std::endl;
- else
- std::ofstream("spline.err",std::ios::out|std::ios::app)
- << "mgcSplineD: " << message[Number(single_error)] << std::endl;
-
- error |= single_error;
-}
-//---------------------------------------------------------------------------
-void mgcSplineD::
-Report (std::ostream& ostr)
-{
- for (unsigned single_error = 1; single_error; single_error <<= 1)
- if ( error & single_error )
- ostr << "mgcSplineD: " << message[Number(single_error)] << std::endl;
-
- error = 0;
-}
-//===========================================================================
-
+# include <cstdlib>
+# include <iostream>
+# include <iomanip>
+# include <cmath>
+# include <ctime>
+# include <cstring>
+
+using namespace std;
+
+# include "spline.h"
[truncated at 1000 lines; 6817 more skipped]
sandbox/magfield/src
--- sandbox/magfield/src/spline.h 2014-04-04 19:02:38 UTC (rev 3084)
+++ sandbox/magfield/src/spline.h 2014-04-04 21:13:31 UTC (rev 3085)
@@ -1,177 +1,99 @@
-// Magic Software, Inc.
-// http://www.magic-software.com
-// Copyright (c) 2000, 2001. All Rights Reserved
-//
-// Source code from Magic Software is supplied under the terms of a license
-// agreement and may not be copied or disclosed except in accordance with the
-// terms of that agreement. The various license agreements may be found at
-// the Magic Software web site. This file is subject to the license
-//
-// FREE SOURCE CODE
-// http://www.magic-software.com/License/free.pdf
-
-#ifndef SPLINE_H
-#define SPLINE_H
-
-#include <iostream>
-#include <math.h>
-
-class mgcSpline
-{
-public:
- mgcSpline () { InitNull(); }
- mgcSpline (int _dimensions, int _degree, const int* _dim, double* _data);
- virtual ~mgcSpline () { FreeStorage(); }
-
- int Dimensions () { return dimensions; }
- int Degree () { return degree; }
- const double** Matrix () { return (const double**)matrix; }
- double DomainMin (int d = 0) { return dom_min[d]; }
- double DomainMax (int d = 0) { return dom_max[d]; }
- int GridMin (int d = 0) { return grid_min[d]; }
- int GridMax (int d = 0) { return grid_max[d]; }
- double (*evaluate_callback)(int);
-
- static const double INVALID_DATA;
-
- // spline evaluation for function interpolation (no derivatives)
- virtual double operator() (double* x) = 0;
-
- // spline evaluation, derivative counts given in dx[]
- virtual double operator() (int* dx, double* x) = 0;
-
- // reallocate spline with new parameters
- mgcSpline& Renew ( int _dimensions, int _degree, const int* _dim,
- double* _data );
-
-protected:
- int dimensions; // N, number of dimensions
- int degree; // D, degree of polynomial spline
- int Dp1; // D+1
- int Dp1_to_N; // power(D+1,N)
- int Dp1_to_2N; // power(D+1,2N)
- int* dim; // dimension sizes dim[0] through dim[N-1]
- double* data; // N-dimensional array of data
- double* dom_min; // minimum allowed value of spline input vector
- double* dom_max; // maximum allowed value of spline input vector
- int* grid_min; // minimum allowed value for current local grid
- int* grid_max; // maximum allowed value for current local grid
- int* base; // base indices for grid of local control points
- int* old_base; // old base indices for grid of local control points
- double* cache; // cache for subblock of data
- double* inter; // intermediate product of data with blending matrix
- double** poly; // poly[N][D+1] for storing polynomials and derivatives
- double** coeff; // coefficients for polynomial construction
- double** matrix; // (D+1)x(D+1) blending matrix
- double* product; // outer tensor product of m with itself N times
- int* skip; // for skipping zero values of mtensor
-
- virtual void EvaluateUnknownData () = 0;
- virtual void ComputeIntermediate () = 0;
- void SetPolynomial (int order, double diff, double* poly);
-
-private:
- void InitNull ();
- void FreeStorage ();
- void Create ();
- static int Choose (int n, int k);
- static double** BlendMatrix (int deg);
-
-// error handling
-public:
- static int verbose;
- static unsigned error;
- static void Report (std::ostream& ostr);
-private:
- static const unsigned invalid_dimensions;
- static const unsigned invalid_degree;
- static const unsigned invalid_dim;
- static const unsigned null_data;
- static const unsigned allocation_failed;
- static const char* message[5];
- static int Number (unsigned single_error);
- static void Report (unsigned single_error);
-};
-
-
-
-class mgcSplineD
-{
-public:
- mgcSplineD () { InitNull(); }
- mgcSplineD (int _dimensions, int _degree, const int* _dim, double* _data);
- virtual ~mgcSplineD () { FreeStorage(); }
-
- int Dimensions () { return dimensions; }
- int Degree () { return degree; }
- const double** Matrix () { return (const double**)matrix; }
- double DomainMin (int d = 0) { return dom_min[d]; }
- double DomainMax (int d = 0) { return dom_max[d]; }
- int GridMin (int d = 0) { return grid_min[d]; }
- int GridMax (int d = 0) { return grid_max[d]; }
- double (*evaluate_callback)(int);
-
- static const double INVALID_DATA;
-
- // spline evaluation for function interpolation (no derivatives)
- virtual double operator() (double* x) = 0;
-
- // spline evaluation, derivative counts given in dx[]
- virtual double operator() (int* dx, double* x) = 0;
-
- // reallocate spline with new parameters
- mgcSplineD& Renew ( int _dimensions, int _degree, const int* _dim,
- double* _data );
-
-protected:
- int dimensions; // N, number of dimensions
- int degree; // D, degree of polynomial spline
- int Dp1; // D+1
- int Dp1_to_N; // power(D+1,N)
- int Dp1_to_2N; // power(D+1,2N)
- int* dim; // dimension sizes dim[0] through dim[N-1]
- double* data; // N-dimensional array of data
- double* dom_min; // minimum allowed value of spline input vector
- double* dom_max; // maximum allowed value of spline input vector
- int* grid_min; // minimum allowed value for current local grid
- int* grid_max; // maximum allowed value for current local grid
- int* base; // base indices for grid of local control points
- int* old_base; // old base indices for grid of local control points
- double* cache; // cache for subblock of data
- double* inter; // intermediate product of data with blending matrix
- double** poly; // poly[N][D+1] for storing polynomials and derivatives
- double** coeff; // coefficients for polynomial construction
- double** matrix; // (D+1)x(D+1) blending matrix
- double* product; // outer tensor product of m with itself N times
- int* skip; // for skipping zero values of mtensor
-
- virtual void EvaluateUnknownData () = 0;
- virtual void ComputeIntermediate () = 0;
- void SetPolynomial (int order, double diff, double* poly);
-
-private:
- void InitNull ();
- void FreeStorage ();
- void Create ();
- static int Choose (int n, int k);
- static double** BlendMatrix (int deg);
-
-// error handling
-public:
- static int verbose;
- static unsigned error;
- static void Report (std::ostream& ostr);
-private:
- static const unsigned invalid_dimensions;
- static const unsigned invalid_degree;
- static const unsigned invalid_dim;
- static const unsigned null_data;
- static const unsigned allocation_failed;
- static const char* message[5];
- static int Number (unsigned single_error);
- static void Report (unsigned single_error);
-};
-
-
-#endif
-
+#ifndef SPLINE_H
+#define SPLINE_H
+# include <string>
+double basis_function_b_val ( double tdata[], double tval );
+double basis_function_beta_val ( double beta1, double beta2, double tdata[],
+ double tval );
+double *basis_matrix_b_uni ( );
+double *basis_matrix_beta_uni ( double beta1, double beta2 );
+double *basis_matrix_bezier ( );
+double *basis_matrix_hermite ( );
+double *basis_matrix_overhauser_nonuni ( double alpha, double beta );
+double *basis_matrix_overhauser_nul ( double alpha );
+double *basis_matrix_overhauser_nur ( double beta );
+double *basis_matrix_overhauser_uni ( );
+double *basis_matrix_overhauser_uni_l ( );
+double *basis_matrix_overhauser_uni_r ( );
+double basis_matrix_tmp ( int left, int n, double mbasis[], int ndata,
+ double tdata[], double ydata[], double tval );
+void bc_val ( int n, double t, double xcon[], double ycon[], double *xval,
+ double *yval );
+double bez_val ( int n, double x, double a, double b, double y[] );
+double bpab_approx ( int n, double a, double b, double ydata[], double xval );
+double *bp01 ( int n, double x );
+double *bpab ( int n, double a, double b, double x );
+int chfev ( double x1, double x2, double f1, double f2, double d1, double d2,
+ int ne, double xe[], double fe[], int next[] );
+int d3_fs ( double a1[], double a2[], double a3[], int n, double b[], double x[] );
+double *d3_mxv ( int n, double a[], double x[] );
+double *d3_np_fs ( int n, double a[], double b[] );
+void d3_print ( int n, double a[], std::string title );
+void d3_print_some ( int n, double a[], int ilo, int jlo, int ihi, int jhi );
+double *d3_uniform ( int n, int *seed );
+void data_to_dif ( int ntab, double xtab[], double ytab[], double diftab[] );
+double dif_val ( int ntab, double xtab[], double diftab[], double xval );
+int i4_max ( int i1, int i2 );
+int i4_min ( int i1, int i2 );
+void least_set ( int point_num, double x[], double f[], double w[],
+ int nterms, double b[], double c[], double d[] );
+double least_val ( int nterms, double b[], double c[], double d[],
+ double x );
+void least_val2 ( int nterms, double b[], double c[], double d[], double x,
+ double *px, double *pxp );
+void least_set_old ( int ntab, double xtab[], double ytab[], int ndeg,
+ double ptab[], double b[], double c[], double d[], double *eps, int *ierror );
+double least_val_old ( double x, int ndeg, double b[], double c[], double d[] );
+void parabola_val2 ( int ndim, int ndata, double tdata[], double ydata[],
+ int left, double tval, double yval[] );
+double pchst ( double arg1, double arg2 );
+double *penta ( int n, double a1[], double a2[], double a3[], double a4[],
+ double a5[], double b[] );
+double r8_abs ( double x );
+double r8_max ( double x, double y );
+double r8_min ( double x, double y );
+double r8_uniform_01 ( int *seed );
+double *r8ge_fs_new ( int n, double a[], double b[] );
+void r8vec_bracket ( int n, double x[], double xval, int *left, int *right );
+void r8vec_bracket3 ( int n, double t[], double tval, int *left );
+double *r8vec_even_new ( int n, double alo, double ahi );
+double *r8vec_indicator_new ( int n );
+int r8vec_order_type ( int n, double x[] );
+void r8vec_print ( int n, double a[], std::string title );
+void r8vec_sort_bubble_a ( int n, double a[] );
+double *r8vec_uniform_new ( int n, double b, double c, int *seed );
+int r8vec_unique_count ( int n, double a[], double tol );
+void r8vec_zero ( int n, double a[] );
+double spline_b_val ( int ndata, double tdata[], double ydata[], double tval );
+double spline_beta_val ( double beta1, double beta2, int ndata, double tdata[],
+ double ydata[], double tval );
+double spline_constant_val ( int ndata, double tdata[], double ydata[],
+ double tval );
+double *spline_cubic_set ( int n, double t[], double y[], int ibcbeg,
+ double ybcbeg, int ibcend, double ybcend );
+double spline_cubic_val ( int n, double t[], double y[], double ypp[],
+ double tval, double *ypval, double *yppval );
+void spline_cubic_val2 ( int n, double t[], double tval, int *left, double y[],
+ double ypp[], double *yval, double *ypval, double *yppval );
+double *spline_hermite_set ( int ndata, double tdata[], double ydata[],
+ double ypdata[] );
+void spline_hermite_val ( int ndata, double tdata[], double c[], double tval,
+ double *sval, double *spval );
+double spline_linear_int ( int ndata, double tdata[], double ydata[], double a,
+ double b );
+void spline_linear_intset ( int int_n, double int_x[], double int_v[],
+ double data_x[], double data_y[] );
+void spline_linear_val ( int ndata, double tdata[], double ydata[],
+ double tval, double *yval, double *ypval );
+double spline_overhauser_nonuni_val ( int ndata, double tdata[],
+ double ydata[], double tval );
+double spline_overhauser_uni_val ( int ndata, double tdata[], double ydata[],
+ double tval );
+void spline_overhauser_val ( int ndim, int ndata, double tdata[], double ydata[],
+ double tval, double yval[] );
+void spline_pchip_set ( int n, double x[], double f[], double d[] );
+void spline_pchip_val ( int n, double x[], double f[], double d[], int ne,
+ double xe[], double fe[] );
+void spline_quadratic_val ( int ndata, double tdata[], double ydata[],
+ double tval, double *yval, double *ypval );
+void timestamp ( );
+#endif
\ No newline at end of file