Fitting NASA-7 polynomials from data

This notebook shows how to generate NASA-7 coefficients from tabulated \(C_p(T)\) data using pyglenn — the reverse of the usual read-only workflow.

Two functions are covered:

  • fit_nasa7 — fit a single temperature interval;

  • fit_nasa7_multi — fit several intervals with continuity of \(C_p\) and its first derivative at each breakpoint.

Both honor the NASA-7 convention H°(298.15) = h_ref (the enthalpy integration constant is b1, not a6).

import numpy as np

from pyglenn import (
    NASACoefficients,
    R_GLENN,
    ThermoDBQuery,
    fit_nasa7,
    fit_nasa7_multi,
)

print('pyglenn imported successfully')

pyglenn imported successfully

Synthetic reference data

We generate “measurements” from a known NASA-7 polynomial (the real O₂ low-temperature coefficients) so we can check how well the fit recovers the source.

# Real O2 NASA-7 coefficients (200-1000 K interval, from the bundled database).
TRUE = NASACoefficients(
    a1=-34255.6342,
    a2=484.700097,
    a3=1.119010961,
    a4=0.00429388924,
    a5=-6.83630052e-07,
    a6=-2.0233727e-09,
    a7=1.039040018e-12,
    b1=-3391.45487,
    b2=18.4969947,
)


def cp(coeffs, T):
    return ThermoDBQuery.calculate_cp(coeffs, T) * R_GLENN


def h(coeffs, T):
    return ThermoDBQuery.calculate_h(coeffs, T) * R_GLENN * T


T = np.linspace(200.0, 1000.0, 50)
cp_meas = np.array([cp(TRUE, t) for t in T])

print(f'{len(T)} points from {T[0]:.1f} to {T[-1]:.1f} K')
print(f'Cp(298.15 K) = {cp(TRUE, 298.15):.3f} J/(mol·K)')

50 points from 200.0 to 1000.0 K
Cp(298.15 K) = 29.378 J/(mol·K)

Single-interval fit

fit_nasa7 recovers the seven heat-capacity coefficients and anchors the enthalpy so that H°(t_ref) = h_ref exactly.

h_ref = h(TRUE, 298.15)

fitted = fit_nasa7(T, cp_meas, t_ref=298.15, h_ref=h_ref)

print('fitted a1..a7 vs true:')
for name in ('a1', 'a2', 'a3', 'a4', 'a5', 'a6', 'a7'):
    true = getattr(TRUE, name)
    got = getattr(fitted, name)
    print(f'  {name:>3}: {got: .8e}   (true {true: .8e})')

fitted a1..a7 vs true:
   a1: -3.42556342e+04   (true -3.42556342e+04)
   a2:  4.84700097e+02   (true  4.84700097e+02)
   a3:  1.11901096e+00   (true  1.11901096e+00)
   a4:  4.29388924e-03   (true  4.29388924e-03)
   a5: -6.83630052e-07   (true -6.83630052e-07)
   a6: -2.02337270e-09   (true -2.02337270e-09)
   a7:  1.03904002e-12   (true  1.03904002e-12)

cp_fit = np.array([cp(fitted, t) for t in T])
max_rel_err = float(np.max(np.abs(cp_fit - cp_meas) / np.abs(cp_meas)))
print(f'max relative Cp error  : {max_rel_err:.2e}')
print(f'H°(298.15) fit vs true : {h(fitted, 298.15):.6f} vs {h_ref:.6f} J/mol')

max relative Cp error  : 6.40e-16
H°(298.15) fit vs true : -0.000013 vs -0.000013 J/mol

Multi-interval fit with continuity

fit_nasa7_multi fits one coefficient set per interval and enforces continuity of \(C_p\) and its first derivative at each breakpoint (e.g. 1000 K), with \(H°\) continuous across intervals.

T_wide = np.linspace(200.0, 6000.0, 300)
cp_wide = np.array([cp(TRUE, t) for t in T_wide])

low, high = fit_nasa7_multi(
    T_wide, cp_wide, breakpoints=[1000.0], t_ref=298.15, h_ref=h_ref
)

print('interval 0 (200-1000 K):', low)
print('interval 1 (1000-6000 K):', high)

interval 0 (200-1000 K): NASACoefficients(a1=-34255.634198868036, a2=484.7000969860279, a3=1.119010961067088, a4=0.004293889239839271, a5=-6.836300517947506e-07, a6=-2.0233727001341847e-09, a7=1.039040018035353e-12, b1=-3391.4548699310194, b2=0.0)
interval 1 (1000-6000 K): NASACoefficients(a1=-34255.634201091336, a2=484.7000970025349, a3=1.1190109609975611, a4=0.004293889240001214, a5=-6.836300520003201e-07, a6=-2.023372699999958e-09, a7=1.0390400179999982e-12, b1=-3391.454870016678, b2=0.0)

Tb = 1000.0


def dcp_dt(coeffs, T):
    a = coeffs
    return (
        -2 * a.a1 / T**3
        - a.a2 / T**2
        + a.a4
        + 2 * a.a5 * T
        + 3 * a.a6 * T**2
        + 4 * a.a7 * T**3
    ) * R_GLENN


print(f'Cp  left vs right : {cp(low, Tb):.6f} vs {cp(high, Tb):.6f} J/(mol·K)')
print(f'dCp left vs right : {dcp_dt(low, Tb):.6f} vs {dcp_dt(high, Tb):.6f} J/(mol·K²)')
print(f'H°  left vs right : {h(low, Tb):.3f} vs {h(high, Tb):.3f} J/mol')

Cp  left vs right : 34.882545 vs 34.882545 J/(mol·K)
dCp left vs right : 0.004959 vs 0.004959 J/(mol·K²)
H°  left vs right : 22707.211 vs 22707.211 J/mol

Summary

  • fit_nasa7 reproduces the reference \(C_p(T)\) to near machine precision and anchors \(H°(298.15)\) exactly.

  • fit_nasa7_multi extends the same to several intervals, with \(C_p\) and its derivative continuous at every breakpoint and \(H°\) continuous across intervals.

  • NASA datasets store one coefficient set per interval — this is exactly how the multiple intervals in thermo.inp are produced.