00001
00002
00003
00004 #include "cddefines.h"
00005 #include "monointerp.h"
00006
00010
00011 namespace {
00012
00013 inline double h00( double t )
00014 {
00015 return (1.0+2.0*t)*(1.0-t)*(1.0-t);
00016 }
00017 inline double h10( double t )
00018 {
00019 return t*(1.0-t)*(1.0-t);
00020 }
00021
00022 inline long bisect ( const std::vector<double> &x, double xval )
00023 {
00024 long n = x.size();
00025 long ilo = 0, ihi = n-1;
00026 if (xval <= x[0])
00027 return -1;
00028 if (xval >= x[n-1])
00029 return n-1;
00030 while (ihi-ilo!=1)
00031 {
00032 long imid = (ilo+ihi)/2;
00033 if (x[imid] > xval)
00034 ihi = imid;
00035 else
00036 ilo = imid;
00037 }
00038 return ilo;
00039 }
00040 }
00041
00042
00043 Monointerp::Monointerp ( const double x[], const double y[], long n )
00044 : m_x(x,x+n), m_y(y,y+n), m_g(n)
00045 {
00046 ASSERT(m_x.size() == m_y.size() && m_x.size() == m_g.size());
00047 std::vector<double> d(n-1),h(n-1);
00048 for (long k=0;k<n-1;++k)
00049 {
00050 h[k] = (x[k+1]-x[k]);
00051 d[k] = (y[k+1]-y[k])/h[k];
00052 }
00053 m_g[0] = d[0];
00054 for (long k=1;k<n-1;++k)
00055 {
00056 m_g[k] = d[k]*d[k-1];
00057 if (m_g[k] > 0.0)
00058 {
00059 double a = (h[k-1]+2.0*h[k])/(3.0*(h[k-1]+h[k]));
00060 m_g[k] /= (a*d[k]+(1-a)*d[k-1]);
00061 }
00062 else
00063 {
00064 m_g[k] = 0.0;
00065 }
00066 }
00067 m_g[n-1] = d[n-2];
00068 }
00069
00070 Monointerp::~Monointerp( void )
00071 {
00072 }
00073
00074
00075 double Monointerp::operator() ( double xval ) const
00076 {
00077 double yval;
00078
00079 if( xval <= m_x[0] )
00080 {
00081 yval = m_y[0];
00082 }
00083 else if( xval >= m_x[m_x.size()-1] )
00084 {
00085 yval = m_y[m_x.size()-1];
00086 }
00087 else
00088 {
00089 long k = bisect( m_x, xval );
00090 double h = m_x[k+1]-m_x[k], t = (xval-m_x[k])/h;
00091 yval = m_y[k]*h00(t) + h*m_g[k]*h10(t)
00092 + m_y[k+1]*h00(1.0-t) - h*m_g[k+1]*h10(1.0-t);
00093 }
00094 return yval;
00095 }