simgear/math/interpolater.cxx

   1 //
   2 // interpolater.cxx -- routines to handle linear interpolation from a table of
   3 //                     x,y   The table must be sorted by "x" in ascending order
   4 //
   5 // Written by Curtis Olson, started April 1998.
   6 //
   7 // Copyright (C) 1998  Curtis L. Olson  - http://www.flightgear.org/~curt
   8 //
   9 // This library is free software; you can redistribute it and/or
  10 // modify it under the terms of the GNU Library General Public
  11 // License as published by the Free Software Foundation; either
  12 // version 2 of the License, or (at your option) any later version.
  13 //
  14 // This library 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 GNU
  17 // Library 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  02110-1301, USA.
  22 //
  23 // $Id$
  24
  25 #ifdef HAVE_CONFIG_H
  26 #  include <simgear_config.h>
  27 #endif
  28
  29 #include <simgear/compiler.h>
  30
  31 #include <string>
  32
  33 #include <simgear/debug/logstream.hxx>
  34 #include <simgear/misc/sgstream.hxx>
  35 #include <simgear/misc/sg_path.hxx>
  36 #include <simgear/props/props.hxx>
  37
  38 #include "interpolater.hxx"
  39
  40 // Constructor -- starts with an empty table.
  41 SGInterpTable::SGInterpTable()
  42 {
  43 }
  44
  45 SGInterpTable::SGInterpTable(const SGPropertyNode* interpolation)
  46 {
  47     if (!interpolation)
  48         return;
  49     std::vector<SGPropertyNode_ptr> entries = interpolation->getChildren("entry");
  50     for (unsigned i = 0; i < entries.size(); ++i)
  51         addEntry(entries[i]->getDoubleValue("ind", 0.0),
  52                  entries[i]->getDoubleValue("dep", 0.0));
  53 }
  54
  55 // Constructor -- loads the interpolation table from the specified
  56 // file
  57 SGInterpTable::SGInterpTable( const std::string& file )
  58 {
  59     SG_LOG( SG_MATH, SG_INFO, "Initializing Interpolator for " << file );
  60
  61     sg_gzifstream in( SGPath::fromUtf8(file) );
  62     if ( !in.is_open() ) {
  63         SG_LOG( SG_GENERAL, SG_ALERT, "Cannot open file: " << file );
  64         return;
  65     }
  66
  67     in >> skipcomment;
  68     while ( in ) {
  69         double ind, dep;
  70         in >> ind >> dep;
  71         in >> std::skipws;
  72         _table[ind] = dep;
  73     }
  74 }
  75
  76
  77 // Constructor -- loads the interpolation table from the specified
  78 // file
  79 SGInterpTable::SGInterpTable( const SGPath& file )
  80 {
  81     SG_LOG( SG_MATH, SG_INFO, "Initializing Interpolator for " << file );
  82
  83     sg_gzifstream in( file );
  84     if ( !in.is_open() ) {
  85         SG_LOG( SG_GENERAL, SG_ALERT, "Cannot open file: " << file );
  86         return;
  87     }
  88
  89     in >> skipcomment;
  90     while ( in ) {
  91       double ind, dep;
  92       in >> ind >> dep;
  93       in >> std::skipws;
  94       _table[ind] = dep;
  95     }
  96 }
  97
  98
  99 // Add an entry to the table.
 100 void SGInterpTable::addEntry (double ind, double dep)
 101 {
 102   _table[ind] = dep;
 103 }
 104
 105 // Given an x value, linearly interpolate the y value from the table
 106 double SGInterpTable::interpolate(double x) const
 107 {
 108   // Empty table??
 109   if (_table.empty())
 110     return 0;
 111
 112   // Find the table bounds for the requested input.
 113   Table::const_iterator upBoundIt = _table.upper_bound(x);
 114   // points to a value outside the map. That is we are out of range.
 115   // use the last entry
 116   if (upBoundIt == _table.end())
 117     return _table.rbegin()->second;
 118
 119   // points to the first key must be lower
 120   // use the first entry
 121   if (upBoundIt == _table.begin())
 122     return upBoundIt->second;
 123
 124   // we know that we do not stand at the beginning, so it is safe to do so
 125   Table::const_iterator loBoundIt = upBoundIt;
 126   --loBoundIt;
 127
 128   // Just do linear interpolation.
 129   double loBound = loBoundIt->first;
 130   double upBound = upBoundIt->first;
 131   double loVal = loBoundIt->second;
 132   double upVal = upBoundIt->second;
 133
 134   // division by zero should not happen since the std::map
 135   // has sorted out duplicate entries before. Also since we have a
 136   // map, we know that we get different first values for different iterators
 137   return loVal + (upVal - loVal)*(x - loBound)/(upBound - loBound);
 138 }
 139
 140
 141 // Destructor
 142 SGInterpTable::~SGInterpTable() {
 143 }
 144
 145