cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
monointerp.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2017 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 
10 // Internal utility functions
11 namespace {
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
43 Monointerp::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 
71 {
72 }
73 
74 // Evaluate interpolant
75 double 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 }
Monointerp(const double x[], const double y[], long n)
Definition: monointerp.cpp:43
std::vector< double > m_g
Definition: monointerp.h:13
double operator()(double xval) const
Definition: monointerp.cpp:75
const std::vector< double > m_x
Definition: monointerp.h:12
~Monointerp(void)
Definition: monointerp.cpp:70
#define ASSERT(exp)
Definition: cddefines.h:617
const std::vector< double > m_y
Definition: monointerp.h:12