cloudy trunk
Loading...
Searching...
No Matches
monointerp.cpp
Go to the documentation of this file.
1/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2 * others. For conditions of distribution and use see copyright notice in license.txt */
3
4#include "cddefines.h"
5#include "monointerp.h"
6
9
10// Internal utility functions
11namespace {
12 // Hermite polynomial basis functions
13 inline double h00( double t )
14 {
15 return (1.0+2.0*t)*(1.0-t)*(1.0-t);
16 }
17 inline double h10( double t )
18 {
19 return t*(1.0-t)*(1.0-t);
20 }
21 // Bisection search
22 inline long bisect ( const std::vector<double> &x, double xval )
23 {
24 long n = x.size();
25 long ilo = 0, ihi = n-1;
26 if (xval <= x[0])
27 return -1;
28 if (xval >= x[n-1])
29 return n-1;
30 while (ihi-ilo!=1)
31 {
32 long imid = (ilo+ihi)/2;
33 if (x[imid] > xval)
34 ihi = imid;
35 else
36 ilo = imid;
37 }
38 return ilo;
39 }
40}
41
42// Constructor for interpolation function object
43Monointerp::Monointerp ( const double x[], const double y[], long n )
44 : m_x(x,x+n), m_y(y,y+n), m_g(n)
45{
46 ASSERT(m_x.size() == m_y.size() && m_x.size() == m_g.size());
47 std::vector<double> d(n-1),h(n-1);
48 for (long k=0;k<n-1;++k)
49 {
50 h[k] = (x[k+1]-x[k]);
51 d[k] = (y[k+1]-y[k])/h[k];
52 }
53 m_g[0] = d[0];
54 for (long k=1;k<n-1;++k)
55 {
56 m_g[k] = d[k]*d[k-1];
57 if (m_g[k] > 0.0)
58 {
59 double a = (h[k-1]+2.0*h[k])/(3.0*(h[k-1]+h[k]));
60 m_g[k] /= (a*d[k]+(1-a)*d[k-1]);
61 }
62 else
63 {
64 m_g[k] = 0.0;
65 }
66 }
67 m_g[n-1] = d[n-2];
68}
69
73
74// Evaluate interpolant
75double Monointerp::operator() ( double xval ) const
76{
77 double yval;
78
79 if( xval <= m_x[0] )
80 {
81 yval = m_y[0];
82 }
83 else if( xval >= m_x[m_x.size()-1] )
84 {
85 yval = m_y[m_x.size()-1];
86 }
87 else
88 {
89 long k = bisect( m_x, xval );
90 double h = m_x[k+1]-m_x[k], t = (xval-m_x[k])/h;
91 yval = m_y[k]*h00(t) + h*m_g[k]*h10(t)
92 + m_y[k+1]*h00(1.0-t) - h*m_g[k+1]*h10(1.0-t);
93 }
94 return yval;
95}
#define ASSERT(exp)
Definition cddefines.h:578
std::vector< double > m_g
Definition monointerp.h:13
double operator()(double xval) const
~Monointerp(void)
const std::vector< double > m_y
Definition monointerp.h:12
Monointerp(const double x[], const double y[], long n)
const std::vector< double > m_x
Definition monointerp.h:12