Cantera  2.1.2
funcs.cpp
Go to the documentation of this file.
1 /**
2  * @file funcs.cpp file containing miscellaneous
3  * numerical functions.
4  */
5 /*
6  * Copyright 2001-2003 California Institute of Technology
7  * See file License.txt for licensing information
8  */
9
10 #include "cantera/numerics/funcs.h"
11
16
17 #include <algorithm>
18
19 using namespace std;
20
21 #ifndef FTN_TRAILING_UNDERSCORE
22 #define _DPOLFT_ dpolft
23 #define _DPCOEF_ dpcoef
24 #else
25 #define _DPOLFT_ dpolft_
26 #define _DPCOEF_ dpcoef_
27 #endif
28
29 extern "C" {
30  int _DPOLFT_(integer* n, doublereal* x, doublereal* y, doublereal* w,
31  integer* maxdeg, integer* ndeg, doublereal* eps, doublereal* r,
32  integer* ierr, doublereal* a);
33
34  int _DPCOEF_(integer* l, doublereal* c, doublereal* tc, doublereal* a);
35 }
36
37 namespace Cantera
38 {
39
40 doublereal linearInterp(doublereal x, const vector_fp& xpts,
41  const vector_fp& fpts)
42 {
43  if (x <= xpts[0]) {
44  return fpts[0];
45  }
46  if (x >= xpts.back()) {
47  return fpts.back();
48  }
49  vector_fp::const_iterator loc =
50  lower_bound(xpts.begin(), xpts.end(), x);
51  int iloc = int(loc - xpts.begin()) - 1;
52  doublereal ff = fpts[iloc] +
53  (x - xpts[iloc])*(fpts[iloc + 1]
54  - fpts[iloc])/(xpts[iloc + 1] - xpts[iloc]);
55  return ff;
56 }
57
58 //! Fits a polynomial function to a set of data points
59 /*!
60  * Given a collection of points X(I) and a set of values Y(I) which
61  * correspond to some function or measurement at each of the X(I),
62  * subroutine DPOLFT computes the weighted least-squares polynomial
63  * fits of all degrees up to some degree either specified by the user
64  * or determined by the routine. The fits thus obtained are in
65  * orthogonal polynomial form. Subroutine DP1VLU may then be
66  * called to evaluate the fitted polynomials and any of their
67  * derivatives at any point. The subroutine DPCOEF may be used to
68  * express the polynomial fits as powers of (X-C) for any specified
69  * point C.
70  *
71  * @param n The number of data points.
72  *
73  * @param x A set of grid points on which the data is specified.
74  * The array of values of the independent variable. These
75  * values may appear in any order and need not all be
76  * distinct. There are n of them.
77  *
78  * @param y array of corresponding function values. There are n of them
79  *
80  * @param w array of positive values to be used as weights. If
81  * W[0] is negative, DPOLFT will set all the weights
82  * to 1.0, which means unweighted least squares error
83  * will be minimized. To minimize relative error, the
84  * user should set the weights to: W(I) = 1.0/Y(I)**2,
85  * I = 1,...,N .
86  *
87  * @param maxdeg maximum degree to be allowed for polynomial fit.
88  * MAXDEG may be any non-negative integer less than N.
89  * Note -- MAXDEG cannot be equal to N-1 when a
90  * statistical test is to be used for degree selection,
91  * i.e., when input value of EPS is negative.
92  *
93  * @param ndeg output degree of the fit computed.
94  *
95  * @param eps Specifies the criterion to be used in determining
96  * the degree of fit to be computed.
97  * (1) If EPS is input negative, DPOLFT chooses the
98  * degree based on a statistical F test of
99  * significance. One of three possible
100  * significance levels will be used: .01, .05 or
101  * .10. If EPS=-1.0 , the routine will
102  * automatically select one of these levels based
103  * on the number of data points and the maximum
104  * degree to be considered. If EPS is input as
105  * -.01, -.05, or -.10, a significance level of
106  * .01, .05, or .10, respectively, will be used.
107  * (2) If EPS is set to 0., DPOLFT computes the
108  * polynomials of degrees 0 through MAXDEG .
109  * (3) If EPS is input positive, EPS is the RMS
110  * error tolerance which must be satisfied by the
111  * fitted polynomial. DPOLFT will increase the
112  * degree of fit until this criterion is met or
113  * until the maximum degree is reached.
114  *
115  * @param r Output vector containing the first LL+1 Taylor coefficients
116  * where LL=ABS(ndeg).
117  * P(X) = r[0] + r[1]*(X-C) + ... + r[ndeg] * (X-C)**ndeg
118  * ( here C = 0.0)
119  *
120  * @return returns the RMS error of the polynomial of degree ndeg .
121  */
122 doublereal polyfit(int n, doublereal* x, doublereal* y, doublereal* w, int maxdeg, int& ndeg, doublereal eps, doublereal* r)
123 {
124  integer nn = n;
125  integer mdeg = maxdeg;
126  integer ndg = ndeg;
127  doublereal epss = eps;
128  integer ierr;
129  int worksize = 3*n + 3*maxdeg + 3;
130  vector_fp awork(worksize,0.0);
131  vector_fp coeffs(n+1, 0.0);
132  doublereal zer = 0.0;
133
134  _DPOLFT_(&nn, x, y, w, &mdeg, &ndg, &epss, &coeffs[0],
135  &ierr, &awork[0]);
136  if (ierr != 1) throw CanteraError("polyfit",
137  "DPOLFT returned error code IERR = " + int2str(ierr) +
138  "while attempting to fit " + int2str(n) + " data points "
139  + "to a polynomial of degree " + int2str(maxdeg));
140  ndeg = ndg;
141  _DPCOEF_(&ndg, &zer, r, &awork[0]);
142  return epss;
143 }
144
145 }
std::string int2str(const int n, const std::string &fmt)
Convert an int to a string using a format converter.
Definition: stringUtils.cpp:40
C interface for Fortran DPOLFT subroutine.
doublereal polyfit(int n, doublereal *x, doublereal *y, doublereal *w, int maxdeg, int &ndeg, doublereal eps, doublereal *r)
Fits a polynomial function to a set of data points.
Definition: funcs.cpp:122
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:68
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
Definition: ct_defs.h:165
Contains declarations for string manipulation functions within Cantera.
Header for a file containing miscellaneous numerical functions.
Definitions for the classes that are thrown when Cantera experiences an error condition (also contain...
doublereal linearInterp(doublereal x, const vector_fp &xpts, const vector_fp &fpts)
Linearly interpolate a function defined on a discrete grid.
Definition: funcs.cpp:40