Branch data     Line data    Source code

       1                 :            : // *****************************************************************************
       2                 :            : /*!
       3                 :            :   \file      src/Base/Table.hpp
       4                 :            :   \copyright 2012-2015 J. Bakosi,
       5                 :            :              2016-2018 Los Alamos National Security, LLC.,
       6                 :            :              2019-2021 Triad National Security, LLC.
       7                 :            :              All rights reserved. See the LICENSE file for details.
       8                 :            :   \brief     Basic functionality for storing and sampling a discrete
       9                 :            :              (y1,y2,...,yN) = f(x) function
      10                 :            :   \details   Basic functionality for storing and sampling a discrete
      11                 :            :              (y1,y2,...,yN) = f(x) function.
      12                 :            : */
      13                 :            : // *****************************************************************************
      14                 :            : #ifndef Table_h
      15                 :            : #define Table_h
      16                 :            : 
      17                 :            : #include <array>
      18                 :            : #include <vector>
      19                 :            : #include <limits>
      20                 :            : 
      21                 :            : #include "Types.hpp"
      22                 :            : #include "Exception.hpp"
      23                 :            : 
      24                 :            : namespace tk {
      25                 :            : 
      26                 :            : //! Type alias for storing a discrete (y1,y2,...,yN) = f(x) function
      27                 :            : //! \tparam N Number of ordinates in the table
      28                 :            : template< std::size_t N >
      29                 :            : using Table = std::vector< std::array< real, N+1 > >;
      30                 :            : 
      31                 :            : //! Sample a discrete (y1,y2,...,yN) = f(x) function at x
      32                 :            : //! \tparam N Number of ordinates in the table
      33                 :            : //! \param[in] x Abscissa to sample at
      34                 :            : //! \param[in] table Table to sample
      35                 :            : //! \return Ordinates sampled
      36                 :            : template< std::size_t N >
      37                 :       4999 : std::array< real, N > sample( real x, const Table< N >& table ) {
      38                 :            : 
      39 [ -  + ][ -  - ]:       4999 :   Assert( not table.empty(), "Empty table to sample from" );
         [ -  - ][ -  - ]
      40                 :            : 
      41                 :            :   // Shortcut for the type of a single line
      42                 :            :   using Line = std::array< tk::real, N+1 >;
      43                 :            :   // Shortcut for the type of all ordinates
      44                 :            :   using Ord = std::array< tk::real, N >;
      45                 :            : 
      46                 :            :   // Lambda to return the abscissa of a Table (return the first value)
      47                 :      23751 :   auto abscissa = []( const Line& t ){ return t[0]; };
      48                 :            : 
      49                 :            :   // Lambda to return ordinates of a tk::Table
      50                 :       7711 :   auto ordinate = []( const Line& t ){
      51                 :            :     Ord o;
      52         [ +  + ]:      44858 :     for (std::size_t i=0; i<N; ++i) o[i] = t[i+1];
      53                 :       7711 :     return o;
      54                 :            :   };
      55                 :            : 
      56                 :       4999 :   auto eps = std::numeric_limits< real >::epsilon();
      57         [ +  + ]:       4999 :   if (x < abscissa(table.front())+eps) return ordinate( table.front() );
      58                 :            : 
      59         [ +  + ]:      10960 :   for (std::size_t i=0; i<table.size()-1; ++i) {
      60                 :       9376 :     auto t1 = abscissa( table[i] );
      61                 :       9376 :     auto t2 = abscissa( table[i+1] );
      62 [ +  - ][ +  + ]:       9376 :     if (t1 < x and x < t2) {
      63                 :       2712 :       auto d = (t2-t1)/(x-t1);
      64                 :       2712 :       auto p = ordinate( table[i] );
      65                 :       2712 :       auto q = ordinate( table[i+1] );
      66                 :            :       Ord r;
      67         [ +  + ]:      15600 :       for (std::size_t j=0; j<N; ++j) r[j] = p[j]+(q[j]-p[j])/d;
      68                 :       2712 :       return r;
      69                 :            :     }
      70                 :            :   }
      71                 :            : 
      72                 :       1584 :   return ordinate( table.back() );
      73                 :            : }
      74                 :            : 
      75                 :            : } // tk::
      76                 :            : 
      77                 :            : #endif // Table_h