src/m_interp.h

   1 /*$Id: m_interp.h,v 1.1 2009-10-23 12:01:45 felix Exp $ -*- C++ -*-
   2  * interpolation on a sorted array
   3  *
   4  * Copyright (C) 2001 Albert Davis
   5  * Author: Albert Davis <aldavis@gnu.org>
   6  *
   7  * This file is part of "Gnucap", the Gnu Circuit Analysis Package
   8  *
   9  * This program is free software; you can redistribute it and/or modify
  10  * it under the terms of the GNU General Public License as published by
  11  * the Free Software Foundation; either version 3, or (at your option)
  12  * any later version.
  13  *
  14  * This program is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  17  * GNU General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU General Public License
  20  * along with this program; if not, write to the Free Software
  21  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
  22  * 02110-1301, USA.
  23  */
  24 //testing=script,sparse 2006.07.13
  25 #ifndef M_INTERP_H
  26 #define M_INTERP_H
  27 #include "m_cpoly.h"
  28 /*--------------------------------------------------------------------------*/
  29 /* interpolate:  linear interpolation on a table.
  30  * Keys must be sorted in increasing order.
  31  */
  32 template <class Iterator>
  33 FPOLY1 interpolate(Iterator begin, Iterator end, double x,
  34                    double below, double above)
  35 {
  36   double f1 = NOT_VALID;
  37   double f0 = NOT_VALID;
  38   if (begin == end) {
  39     untested();
  40     throw Exception("interpolate table is empty");
  41   }
  42   --end;
  43   if (begin == end) { // only 1 entry -- constant
  44     f1 = (x < (*begin).first)
  45       ? ((below != NOT_INPUT) ? below : 0.)
  46       : ((above != NOT_INPUT) ? above : 0.);
  47     f0 = (*begin).second + (x - (*begin).first) * f1;
  48   }else{
  49     ++begin;
  50
  51     // x is the search key
  52     // xx is the search key converted to a "pair" as required by upper_bound
  53     // upper might point to a match for the key, exact match
  54     // if not, it points just above it, no exact match
  55     // lower points just below where a match would go
  56     // so the real match is between upper and lower
  57     DPAIR    xx(x,BIGBIG);
  58     Iterator upper = upper_bound(begin, end, xx);
  59     Iterator lower = upper-1;
  60
  61     // set f1 (derivative)
  62     if ((upper == end) && (x > (*upper).first) && (above != NOT_INPUT)) {
  63       // x is out of bounds, above
  64       lower = upper;
  65       if (above != 0.) {
  66         untested();
  67       }
  68       f1 = above;
  69     }else if ((upper == begin) && (x < (*lower).first) && (below!=NOT_INPUT)) {
  70       // x is out of bounds, below
  71       untested();
  72       if (below != 0.) {
  73         untested();
  74       }else{
  75         untested();
  76       }
  77       f1 = below;
  78     }else if ((*upper).first <= (*lower).first) {untested();
  79       throw Exception("interpolate table is not sorted or has duplicate keys");
  80       f1 = 0.;
  81     }else{
  82       // ordinary interpolation
  83       assert((*upper).first != (*lower).first);
  84       f1 = ((*upper).second-(*lower).second) / ((*upper).first-(*lower).first);
  85     }
  86
  87     f0 = (*lower).second + (x - (*lower).first) * f1;
  88   }
  89   assert(f1 != NOT_VALID);
  90   assert(f0 != NOT_VALID);
  91   return FPOLY1(x, f0, f1);
  92 }
  93 /*--------------------------------------------------------------------------*/
  94 /*--------------------------------------------------------------------------*/
  95 #endif
  96 // vim:ts=8:sw=2:noet: