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