hermite_cubic.cpp

Go to the documentation of this file.

 
/* ===============================================================================================================================
 
   This file is part of CoastalME, the Coastal Modelling Environment.
 
   CoastalME is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version.
 
   This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
 
   You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 
===============================================================================================================================*/
#include <iostream>
using std::cout;
using std::cerr;
using std::endl;
using std::ios;
 
#include "cme.h"
#include "hermite_cubic.h"
 
//===============================================================================================================================
//
// Purpose:
//
// R8VEC_BRACKET3 finds the interval containing or nearest a given value.
//
// Discussion:
//
// An R8VEC is a vector of R8's.
//
// The routine always returns the index LEFT of the sorted array
// T with the property that either
// *  T is contained in the interval [ T[LEFT], T[LEFT+1] ], or
// *  T < T[LEFT] = T[0], or
// *  T > T[LEFT+1] = T[N-1].
//
// The routine is useful for interpolation problems, where
// the abscissa must be located within an interval of data
// abscissas for interpolation, or the "nearest" interval
// to the (extreme) abscissa must be found so that extrapolation
// can be carried out.
//
// This version of the function has been revised so that the value of
// LEFT that is returned uses the 0-based indexing natural to C++.
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 30 April 2009
//
// Author:
//
// John Burkardt
//
// Parameters:
//
// Input, int N, length of the input array.
//
// Input, double T[N], an array that has been sorted into ascending order.
//
// Input, double TVAL, a value to be bracketed by entries of T.
//
// Input/output, int *LEFT.
// On input, if 0 <= LEFT <= N-2, LEFT is taken as a suggestion for the
// interval [ T[LEFT-1] T[LEFT] ] in which TVAL lies.  This interval
// is searched first, followed by the appropriate interval to the left
// or right.  After that, a binary search is used.
// On output, LEFT is set so that the interval [ T[LEFT], T[LEFT+1] ]
// is the closest to TVAL; it either contains TVAL, or else TVAL
// lies outside the interval [ T[0], T[N-1] ].
//
//===============================================================================================================================
//===============================================================================================================================
void r8vec_bracket3(int const n, double const* t, double const tval, int* left)
{
   int high;
   int low;
   int mid;
 
   // Check the input data
   if (n < 2)
   {
      cerr << endl << "R8VEC_BRACKET3 - Fatal error! N must be at least 2." << endl;
      return;
   }
 
   // If *left is not between 0 and n-2, set it to the middle value
   if (( * left < 0) || (n - 2 < * left))
   {
      * left = (n - 1) / 2;
   }
 
   // CASE 1: TVAL < T[*LEFT]: search for TVAL in (T[I],T[I+1]), for I = 0 to *LEFT-1.
   if ( tval < t[ * left] )
   {
      if ( * left == 0 )
      {
         return;
      }
 
      else if ( * left == 1 )
      {
         * left = 0;
         return;
      }
 
      else if ( t[ * left - 1] <= tval )
      {
         * left = * left - 1;
         return;
      }
 
      else if ( tval <= t[1] )
      {
         * left = 0;
         return;
      }
 
      //
      // ...Binary search for TVAL in (T[I],T[I+1]), for I = 1 to *LEFT-2.
      //
      low = 1;
      high = * left - 2;
 
      for ( ; ; )
      {
         if ( low == high )
         {
            * left = low;
            return;
         }
 
         mid = ( low + high + 1 ) / 2;
 
         if ( t[mid] <= tval )
         {
            low = mid;
         }
 
         else
         {
            high = mid - 1;
         }
      }
   }
 
   //
   // CASE 2: T[*LEFT+1] < TVAL:
   // Search for TVAL in (T[I],T[I+1]) for intervals I = *LEFT+1 to N-2.
   //
   else if ( t[ * left + 1] < tval )
   {
      if ( * left == n - 2 )
      {
         return;
      }
 
      else if ( * left == n - 3 )
      {
         * left = * left + 1;
         return;
      }
 
      else if ( tval <= t[ * left + 2] )
      {
         * left = * left + 1;
         return;
      }
 
      else if ( t[n - 2] <= tval )
      {
         * left = n - 2;
         return;
      }
 
      // ...Binary search for TVAL in (T[I],T[I+1]) for intervals I = *LEFT+2 to N-3.
      low = * left + 2;
      high = n - 3;
 
      for ( ; ; )
      {
         if ( low == high )
         {
            * left = low;
            return;
         }
 
         mid = ( low + high + 1 ) / 2;
 
         if ( t[mid] <= tval )
         {
            low = mid;
         }
 
         else
         {
            high = mid - 1;
         }
      }
   }
 
   // CASE 3: T[*LEFT] <= TVAL <= T[*LEFT+1]: T is just where the user said it might be.
   else
   {
   }
 
   return;
}
 
//===============================================================================================================================
//
// Purpose:
//
// HERMITE_CUBIC_VALUE evaluates a Hermite cubic polynomial.
//
// Discussion:
//
// The input arguments can be vectors.
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 13 February 2011
//
// Author:
//
// John Burkardt
//
// Reference:
//
// Fred Fritsch, Ralph Carlson,
// Monotone Piecewise Cubic Interpolation,
// SIAM Journal on Numerical Analysis,
// Volume 17, Number 2, April 1980, pages 238-246.
//
// Parameters:
//
// Input, double X1, F1, D1, the left endpoint, function value
// and derivative.
//
// Input, double X2, F2, D2, the right endpoint, function value
// and derivative.
//
// Input, int N, the number of evaluation points.
//
// Input, double X[N], the points at which the Hermite cubic
// is to be evaluated.
//
// Output, double F[N], D[N], S[N], T[N], the value and first
// three derivatives of the Hermite cubic at X.
//
//===============================================================================================================================
void hermite_cubic_value(double const x1, double const f1, double const d1, double const x2, double const f2, double const d2, int const n, double const* x, double* f, double* d, double* s, double* t)
{
   double c2;
   double c3;
   double df;
   double h;
   int i;
 
   h = x2 - x1;
   df = (f2 - f1) / h;
 
   c2 = - ( 2.0 * d1 - 3.0 * df + d2 ) / h;
   c3 = ( d1 - 2.0 * df + d2 ) / h / h;
 
   for (i = 0; i < n; i++)
   {
      f[i] = f1 + ( x[i] - x1 ) * ( d1
                                    + ( x[i] - x1 ) * ( c2
                                       + ( x[i] - x1 ) * c3 ) );
      d[i] = d1 + ( x[i] - x1 ) * ( 2.0 * c2
                                    + ( x[i] - x1 ) * 3.0 * c3 );
      s[i] = 2.0 * c2 + ( x[i] - x1 ) * 6.0 * c3;
      t[i] = 6.0 * c3;
   }
 
   return;
}
 
//===============================================================================================================================
//
// Purpose:
//
// HERMITE_CUBIC_SPLINE_VALUE evaluates a Hermite cubic spline.
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 13 February 2011
//
// Author:
//
// John Burkardt
//
// Reference:
//
// Fred Fritsch, Ralph Carlson,
// Monotone Piecewise Cubic Interpolation,
// SIAM Journal on Numerical Analysis,
// Volume 17, Number 2, April 1980, pages 238-246.
//
// Parameters:
//
// Input, int NN, the number of data points.
//
// Input, double XN[NN], the coordinates of the data points.
// The entries in XN must be in strictly ascending order.
//
// Input, double FN[NN], the function values.
//
// Input, double DN[NN], the derivative values.
//
// Input, int N, the number of sample points.
//
// Input, double X[N], the coordinates of the sample points.
//
// Output, double F[N], the function value at the sample points.
//
// Output, double D[N], the derivative value at the sample points.
//
// Output, double S[N], the second derivative value at the
// sample points.
//
// Output, double T[N], the third derivative value at the
// sample points.
//
//===============================================================================================================================
//===============================================================================================================================
void hermite_cubic_spline_value(int const nn, double* const xn, double* const fn, double* const dn, int const n, double* const x, double* f, double* d, double* s, double* t)
{
   int left = n / 2;
 
   for (int i = 0; i < n; i++)
   {
      r8vec_bracket3(nn, xn, x[i], &left);
 
      hermite_cubic_value(xn[left], fn[left], dn[left], xn[left + 1], fn[left + 1], dn[left + 1], 1, x + i, f + i, d + i, s + i, t + i);
   }
 
   return;
}