Commit in sandbox/magfield/src on MAIN | |||
HPS_b18d36_By_0_0_z.dat | +301 | added 3085 | |
OneDOneComponentMagneticFieldMap.cpp | +66 | added 3085 | |
OneDimOneComponentMagneticFieldMap.h | +41 | added 3085 | |
OneDimOneComponentMagneticFieldMap_t.cpp | +159 | added 3085 | |
spline.cpp | +6827 | -987 | 3084 -> 3085 |
spline.h | +99 | -177 | 3084 -> 3085 |
+7493 | -1164 |
Classes to handle one value of the BField measured along one coordinate.
--- sandbox/magfield/src/HPS_b18d36_By_0_0_z.dat (rev 0) +++ sandbox/magfield/src/HPS_b18d36_By_0_0_z.dat 2014-04-04 21:13:31 UTC (rev 3085) @@ -0,0 +1,301 @@
+0 -0.5006 +0.5 -0.5006 +1 -0.5006 +1.5 -0.5006 +2 -0.5006 +2.5 -0.5006 +3 -0.5006 +3.5 -0.5006 +4 -0.5006 +4.5 -0.5006 +5 -0.5006 +5.5 -0.5006 +6 -0.5006 +6.5 -0.5006 +7 -0.5006 +7.5 -0.5006 +8 -0.5006 +8.5 -0.5006 +9 -0.5006 +9.5 -0.5006 +10 -0.5006 +10.5 -0.5006 +11 -0.5006 +11.5 -0.5006 +12 -0.5006 +12.5 -0.5006 +13 -0.5006 +13.5 -0.5006 +14 -0.5006 +14.5 -0.5006 +15 -0.5006 +15.5 -0.5005 +16 -0.5005 +16.5 -0.5005 +17 -0.5005 +17.5 -0.5005 +18 -0.5005 +18.5 -0.5005 +19 -0.5005 +19.5 -0.5005 +20 -0.5005 +20.5 -0.5005 +21 -0.5005 +21.5 -0.5005 +22 -0.5004 +22.5 -0.5004 +23 -0.5004 +23.5 -0.5004 +24 -0.5004 +24.5 -0.5003 +25 -0.5003 +25.5 -0.5002 +26 -0.5002 +26.5 -0.5001 +27 -0.5001 +27.5 -0.5 +28 -0.4999 +28.5 -0.4998 +29 -0.4997 +29.5 -0.4995 +30 -0.4994 +30.5 -0.4992 +31 -0.4989 +31.5 -0.4987 +32 -0.4984 +32.5 -0.4981 +33 -0.4978 +33.5 -0.4973 +34 -0.4967 +34.5 -0.4961 +35 -0.4953 +35.5 -0.4945 +36 -0.4935 +36.5 -0.4925 +37 -0.4914 +37.5 -0.49 +38 -0.4883 +38.5 -0.4863 +39 -0.484 +39.5 -0.4814 +40 -0.4786 +40.5 -0.4754 +41 -0.4719 +41.5 -0.4679 +42 -0.4634 +42.5 -0.4584 +43 -0.4528 +43.5 -0.4467 +44 -0.44 +44.5 -0.4329 +45 -0.4252 +45.5 -0.417 +46 -0.4082 +46.5 -0.3991 +47 -0.3896 +47.5 -0.3798 +48 -0.3696 +48.5 -0.3593 +49 -0.3489 +49.5 -0.3385 +50 -0.3279 +50.5 -0.3175 +51 -0.3071 +51.5 -0.297 +52 -0.2869 +52.5 -0.2771 +53 -0.2675 +53.5 -0.2582 +54 -0.2492 +54.5 -0.2404 +55 -0.2319 +55.5 -0.2237 +56 -0.2158 +56.5 -0.2082 +57 -0.2008 +57.5 -0.1937 +58 -0.1869 +58.5 -0.1803 +59 -0.174 +59.5 -0.1679 +60 -0.162 +60.5 -0.1564 +61 -0.151 +61.5 -0.1458 +62 -0.1408 +62.5 -0.1359 +63 -0.1313 +63.5 -0.1268 +64 -0.1225 +64.5 -0.1184 +65 -0.1144 +65.5 -0.1106 +66 -0.1069 +66.5 -0.1033 +67 -0.0999 +67.5 -0.0966 +68 -0.0934 +68.5 -0.0904 +69 -0.0874 +69.5 -0.0846 +70 -0.0818 +70.5 -0.0792 +71 -0.0766 +71.5 -0.0742 +72 -0.0718 +72.5 -0.0695 +73 -0.0673 +73.5 -0.0652 +74 -0.0632 +74.5 -0.0612 +75 -0.0593 +75.5 -0.0575 +76 -0.0557 +76.5 -0.054 +77 -0.0523 +77.5 -0.0507 +78 -0.0492 +78.5 -0.0477 +79 -0.0463 +79.5 -0.0449 +80 -0.0435 +80.5 -0.0422 +81 -0.041 +81.5 -0.0398 +82 -0.0386 +82.5 -0.0375 +83 -0.0364 +83.5 -0.0353 +84 -0.0343 +84.5 -0.0333 +85 -0.0324 +85.5 -0.0315 +86 -0.0306 +86.5 -0.0297 +87 -0.0289 +87.5 -0.0282 +88 -0.0275 +88.5 -0.0268 +89 -0.0261 +89.5 -0.0254 +90 -0.0247 +90.5 -0.0241 +91 -0.0234 +91.5 -0.0228 +92 -0.0222 +92.5 -0.0216 +93 -0.021 +93.5 -0.0204 +94 -0.0198 +94.5 -0.0193 +95 -0.0188 +95.5 -0.0183 +96 -0.0178 +96.5 -0.0173 +97 -0.0168 +97.5 -0.0163 +98 -0.0159 +98.5 -0.0154 +99 -0.015 +99.5 -0.0146 +100 -0.0142 +100.5 -0.0139 +101 -0.0135 +101.5 -0.0131 +102 -0.0128 +102.5 -0.0125 +103 -0.0122 +103.5 -0.0119 +104 -0.0116 +104.5 -0.0113 +105 -0.0111 +105.5 -0.0109 +106 -0.0106 +106.5 -0.0104 +107 -0.0102 +107.5 -0.01 +108 -0.0099 +108.5 -0.0097 +109 -0.0095 +109.5 -0.0093 +110 -0.0091 +110.5 -0.0089 +111 -0.0087 +111.5 -0.0085 +112 -0.0083 +112.5 -0.0081 +113 -0.0079 +113.5 -0.0078 +114 -0.0076 +114.5 -0.0074 +115 -0.0073 +115.5 -0.0071 +116 -0.0069 +116.5 -0.0068 +117 -0.0066 +117.5 -0.0065 +118 -0.0064 +118.5 -0.0062 +119 -0.0061 +119.5 -0.006 +120 -0.0058 +120.5 -0.0057 +121 -0.0056 +121.5 -0.0055 +122 -0.0054 +122.5 -0.0052 +123 -0.0051 +123.5 -0.005 +124 -0.0049 +124.5 -0.0048 +125 -0.0047 +125.5 -0.0047 +126 -0.0046 +126.5 -0.0045 +127 -0.0044 +127.5 -0.0043 +128 -0.0043 +128.5 -0.0042 +129 -0.0041 +129.5 -0.0041 +130 -0.004 +130.5 -0.0039 +131 -0.0038 +131.5 -0.0038 +132 -0.0037 +132.5 -0.0037 +133 -0.0036 +133.5 -0.0035 +134 -0.0035 +134.5 -0.0034 +135 -0.0033 +135.5 -0.0033 +136 -0.0032 +136.5 -0.0032 +137 -0.0031 +137.5 -0.0031 +138 -0.003 +138.5 -0.003 +139 -0.0029 +139.5 -0.0029 +140 -0.0028 +140.5 -0.0028 +141 -0.0027 +141.5 -0.0027 +142 -0.0026 +142.5 -0.0026 +143 -0.0025 +143.5 -0.0025 +144 -0.0025 +144.5 -0.0024 +145 -0.0024 +145.5 -0.0023 +146 -0.0023 +146.5 -0.0023 +147 -0.0022 +147.5 -0.0022 +148 -0.0022 +148.5 -0.0021 +149 -0.0021 +149.5 -0.0021 +150 -0.0021
--- sandbox/magfield/src/OneDOneComponentMagneticFieldMap.cpp (rev 0) +++ sandbox/magfield/src/OneDOneComponentMagneticFieldMap.cpp 2014-04-04 21:13:31 UTC (rev 3085) @@ -0,0 +1,66 @@
+#include "OneDimOneComponentMagneticFieldMap.h" +#include <fstream> +#include <iostream> +#include <cmath> + +#include "spline.h" + +using namespace std; + +OneDimOneComponentMagneticFieldMap::OneDimOneComponentMagneticFieldMap( COORDINATE coord, BVAL bval, const char* filename, double offset, bool isSymmetric ) + :_coord(coord),_bval(bval),_offset(offset),_isSymmetric(false) +{ + cout << "\n-----------------------------------------------------------" + << "\n Magnetic field" + << "\n-----------------------------------------------------------"; + + cout << "\n ---> " "Reading the field grid from " << filename << " ... " << endl; + ifstream inFile( filename ); // Open the file for reading. + _dim = count(istreambuf_iterator<char>(inFile), + istreambuf_iterator<char>(), '\n'); + inFile.seekg(0); + + _ord = new double[_dim]; + _val = new double[_dim]; + + for(int i=0; i<_dim; ++i) + { + inFile >> _ord[i] >> _val[i]; + } + inFile.close(); + _min = _ord[0]; + _max = _ord[_dim-1]; + +} + +OneDimOneComponentMagneticFieldMap::~OneDimOneComponentMagneticFieldMap() +{ + delete[] _ord; + delete[] _val; +} + +void OneDimOneComponentMagneticFieldMap::GetFieldValue(const double point[4], double* Bfield ) const +{ + double val = point[_coord]; + Bfield[0] = 0.; + Bfield[1] = 0.; + Bfield[2] = 0.; + if(_isSymmetric == false) + { + if(val < _min || val >_max) return; + Bfield[_bval] = spline_b_val ( _dim, _ord, _val, val); + return; + } + else + { + if(val >= _min && val <= _max) + { + Bfield[_bval] = spline_b_val ( _dim, _ord, _val, val); + return; + } + if(abs(val) < abs(_min) || abs(val) > abs(_max)) return; + Bfield[_bval] = spline_b_val ( _dim, _ord, _val, -val); + return; + } +} +
--- sandbox/magfield/src/OneDimOneComponentMagneticFieldMap.h (rev 0) +++ sandbox/magfield/src/OneDimOneComponentMagneticFieldMap.h 2014-04-04 21:13:31 UTC (rev 3085) @@ -0,0 +1,41 @@
+#ifndef ONEDONECOMPONENTMAGNETICFIELDMAP_H +#define ONEDONECOMPONENTMAGNETICFIELDMAP_H + +class OneDimOneComponentMagneticFieldMap +{ +public: + enum COORDINATE {X, Y, Z}; + enum BVAL {BX,BY,BZ}; +private: + //ordinate axis + COORDINATE _coord; + // ordinate values + double* _ord; + // abscissa field component + BVAL _bval; + // abscissa values + double* _val; + //the dimension of the field + int _dim; + // The physical limits of the defined region + double _min, _max; + // Offset if field map is not in global coordinates + double _offset; + //is map symmetric in ordinate? + bool _isSymmetric; + + +public: + OneDimOneComponentMagneticFieldMap(COORDINATE coord, BVAL bval, const char* filename, double offset, bool isSymmetric ); + ~OneDimOneComponentMagneticFieldMap(); + double min() + { + return _min; + } + double max() + { + return _max; + } + void GetFieldValue( const double Point[4], double* Bfield ) const; +}; +#endif
\ No newline at end of file
--- sandbox/magfield/src/OneDimOneComponentMagneticFieldMap_t.cpp (rev 0) +++ sandbox/magfield/src/OneDimOneComponentMagneticFieldMap_t.cpp 2014-04-04 21:13:31 UTC (rev 3085) @@ -0,0 +1,159 @@
+#include <iostream> +#include <fstream> +#include <cassert> +#include <cmath> +#include "OneDimOneComponentMagneticFieldMap.h" +using namespace std; + +int main () +{ + + + // read in the field and find how many entries we have + ifstream inFile( "C:/work/svn/sandbox/magfield/src/HPS_b18d36_By_0_0_z.dat" ); // Open the file for reading. + const int dim = count(istreambuf_iterator<char>(inFile), + istreambuf_iterator<char>(), '\n'); + inFile.seekg(0); + + double* z = new double[dim]; + double* By = new double[dim]; + + for(int i=0; i<dim; ++i) + { + inFile >> z[i] >> By[i]; + } + double zmin = z[0]; + double zmax = z[dim-1]; + + + + OneDimOneComponentMagneticFieldMap::BVAL val = OneDimOneComponentMagneticFieldMap::BVAL::BX; + OneDimOneComponentMagneticFieldMap::COORDINATE coord = OneDimOneComponentMagneticFieldMap::COORDINATE::Z; + + OneDimOneComponentMagneticFieldMap map(coord,val, "C:/work/svn/sandbox/magfield/src/HPS_b18d36_By_0_0_z.dat", 0., true); + + double B[3]; + double p[4] = {0.,0.,0.}; + + double min = map.min(); + double max = map.max(); + double epsilon = .0005; + cout << "min "<<min << " " << max << endl; + for (int i=0; i<dim; ++i) + { + p[coord] = z[i]; + map.GetFieldValue( p, B ); + cout << p[0] << " " << p[1] << " "<< p[2] << " : "<< B[0] << " "<< B[1] << " "<< B[2] << " " << endl; + if(abs(B[val]-By[i]) >= epsilon ) cout << B[val] << " " << By[i] << endl; + assert( abs(B[val]-By[i]) < epsilon ); + } + // + // // read in the field and find how many entries we have + // ifstream inFile( "C:/work/svn/sandbox/magfield/src/HPS_b18d36_By_0_0_z.dat" ); // Open the file for reading. + // const int dim = count(istreambuf_iterator<char>(inFile), + // istreambuf_iterator<char>(), '\n'); + // inFile.seekg(0); + // + // double* z = new double[dim]; + // double* By = new double[dim]; + // double zmin = 0; + // double zmax = 0; + // for(int i=0; i<dim; ++i) + // { + // inFile >> z[i] >> By[i]; + // } + // + // int i; + // int j; + // int jhi; + // char mark; + // int nsample = 10; + // double pi = 3.141592653589793; + // double* tdata = new double[dim]; + // + // double thi; + // double tlo; + // double tval; + // double* ydata = new double[dim]; + // double yval; + // + // cout << "\n"; + // cout << "TEST13\n"; + // cout << " SPLINE_B_VAL evaluates the\n"; + // cout << " B spline.\n"; + // cout << "\n"; + // cout << " TDATA YDATA\n"; + // cout << "\n"; + // + // for ( i = 0; i < dim; i++ ) + // { + // tdata[i] = z[i]; + // ydata[i] = By[i]; + // cout << setw(12) << tdata[i] << " " + // << setw(12) << ydata[i] << "\n"; + // } + // + // cout << "\n" ; + // cout << " T, Spline(T)\n"; + // cout << "\n"; + // + // for ( i = 0; i <= dim; i++ ) + // { + // if ( i == 0 ) + // { + // tlo = tdata[1] - 0.5 * ( tdata[1] - tdata[0] ); + // thi = tdata[0]; + // } + // else if ( i < dim ) + // { + // tlo = tdata[i-1]; + // thi = tdata[i]; + // } + // else if ( dim <= i ) + // { + // tlo = tdata[dim-1]; + // thi = tdata[dim-1] + 0.5 * ( tdata[dim-1] - tdata[dim-2] ); + // } + // + // if ( i < dim ) + // { + // jhi = nsample - 1; + // } + // else + // { + // jhi = nsample; + // } + // + // for ( j = 0; j <= jhi; j++ ) + // { + // tval = ( ( double ) ( nsample - j ) * tlo + // + ( double ) ( j ) * thi ) + // / ( double ) ( nsample ); + // + // yval = spline_b_val ( dim, tdata, ydata, tval ); + // + // //if ( 0 < i && j == 0 ) + // //{ + // // mark = '*'; + // //} + // //else + // //{ + // mark = ' '; + // //} + // + // cout << " " + // << mark << " " + // << setw(12) << tval << " " + // << setw(12) << yval << "\n"; + // + // } + // + // } + // + // + ////%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + ////clean up + // delete[] z; + // delete[] By; + return 0; +}
--- 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/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
Use REPLY-ALL to reply to list
To unsubscribe from the LCDET-SVN list, click the following link:
https://listserv.slac.stanford.edu/cgi-bin/wa?SUBED1=LCDET-SVN&A=1