PropertyTable.h

Go to the documentation of this file.

// PropertyTable: a 1D temperature-indexed lookup with linear interpolation.
//
// This is the pluggable property backend used by ThermalPlasmaColumn1D for
// sigma(T), kappa(T), P_rad(T), and the optional seed profile T_seed(r).
// Today it is always fed rizer's tabulated LTE data (the same CSVs that
// ElenbaasHeller uses), so the numerical and analytical solutions share
// *identical* property models and any residual difference (~1–3 % in T_center
// and current for H2) is purely the analytical linearisation error.
//
// Design contract:
//   eval(T)  ->  linearly interpolated value, clamped to the table endpoints
//               outside [T_min, T_max] (no extrapolation, no NaN).
//   empty()  ->  true for a default-constructed table; eval returns 0.0.
//               Used by the P_rad slot when radiation is disabled.
//
// Future:  a Cantera-native backend (equilibrate("TP") + IonGasTransport for
// sigma/kappa; SRI/JANAF for ionization equilibrium) can satisfy the same
// eval(T) contract without touching ThermalPlasmaColumn1D.  Note that built-in
// Cantera plasma mechanisms do not cover H2 LTE thermal ionisation, so that
// path needs a suitable mechanism file.
 
#ifndef RIZER_PROPERTY_TABLE_H
#define RIZER_PROPERTY_TABLE_H
 
#include <vector>
#include <stdexcept>
#include <algorithm>
 
namespace rizer {
 
class PropertyTable {
public:
    PropertyTable() = default;

    PropertyTable(std::vector<double> T, std::vector<double> values)
        : m_T(std::move(T)), m_v(std::move(values))
    {
        if (m_T.size() != m_v.size() || m_T.size() < 2) {
            throw std::invalid_argument(
                "PropertyTable: T and values must have equal length >= 2");
        }
    }

    static PropertyTable constant(double value) {
        return PropertyTable({0.0, 1.0e9}, {value, value});
    }

    static PropertyTable fromTableOrScalar(const std::vector<double>& x,
                                            const std::vector<double>& v) {
        if (x.size() >= 2 && x.size() == v.size()) return PropertyTable(x, v);
        if (v.size() == 1) return PropertyTable::constant(v[0]);
        return PropertyTable();
    }

    bool empty() const { return m_T.empty(); }

    double eval(double T) const {
        if (m_T.empty()) {
            return 0.0;
        }
        if (T <= m_T.front()) {
            return m_v.front();
        }
        if (T >= m_T.back()) {
            return m_v.back();
        }
        // upper_bound gives the first node strictly greater than T; subtract 1
        // to get the left bracket of the containing interval.
        auto it = std::upper_bound(m_T.begin(), m_T.end(), T);
        size_t hi = static_cast<size_t>(it - m_T.begin());
        size_t lo = hi - 1;
        double w = (T - m_T[lo]) / (m_T[hi] - m_T[lo]);
        return m_v[lo] + w * (m_v[hi] - m_v[lo]);
    }
 
private:
    std::vector<double> m_T;  
    std::vector<double> m_v;  
};
 
} // namespace rizer
 
#endif