Different output from running Python script and cea binary directly · nasa/cea · Discussion #63

minh-dng
Mar 21, 2026

I compared the same problem running through the binary and the Python distribution in bin/cea. You could see the numerical differences between the two program (despite the rounding of Python output), e.g. (C*).

Just want to double check that I haven't made any mistake in my Python program.

The input file is

reac
 fuel Acrylonitrile C 3 H 3 N 1 wt%=39.92
 h,kj/mol=140. t(k)=298.15
 fuel Butadiene C 4 H 6 wt%=47.37
 h,kj/mol=90.50 t(k)=298.15
 fuel Styrene C 8 H 8 wt%=12.71
 h,kj/mol=103.4 t(k)=298.15
 oxid N2O wt%=100 t(k)=298.15
prob rocket eql o/f=5.098591, p,bar=1.01325, supar=4.000000 ,pi/p=1
end

The input code to cea is

import cea
import numpy as np
acrylonitrile_args = {
 "name": "Acrylonitrile",
 "formula": {"C": 3, "H": 3, "N": 1},
 "enthalpy": 140,
 "temperature": 298.15,
 "molecular_weight": 53.0626,
}
"""https://webbook.nist.gov/cgi/cbook.cgi?ID=C107131&Mask=2"""
butadiene_args = {
 "name": "Butadiene",
 "formula": {"C": 4, "H": 6},
 "enthalpy": 90.50,
 "temperature": 298.15,
 "molecular_weight": 54.0904,
}
"""https://webbook.nist.gov/cgi/cbook.cgi?ID=C106990&Mask=7"""
styrene_args = {
 "name": "Styrene",
 "formula": {"C": 8, "H": 8},
 "enthalpy": 103.4,
 "temperature": 298.15,
 "molecular_weight": 104.1491,
}
"""https://webbook.nist.gov/cgi/cbook.cgi?ID=C100425&Mask=F"""
reac_names = [
 cea.Reactant(**acrylonitrile_args),
 cea.Reactant(**butadiene_args),
 cea.Reactant(**styrene_args),
 "N2O",
]
T_reactant = np.array([298.15] * 4) # Reactant temperatures (K)
fuel_weights = np.array([0.3992, 0.4737, 0.1271, 0.0])
oxidant_weights = np.array([0.0] * 3 + [1.0])
of_ratio = 5.098591
pi_p = 1.0 # Pressure ratio
supar = 4.0 # Supersonic area ratio
of = 5.098591 # O/F ratio
pc = 1.01325 # Chamber pressure (bar)
reac = cea.Mixture(reac_names)
prod = cea.Mixture(reac_names, products_from_reactants=True)
solver = cea.RocketSolver(prod, reactants=reac)
solution = cea.RocketSolution(solver)
weights = reac.of_ratio_to_weights(oxidant_weights, fuel_weights, of_ratio)
hc = reac.calc_property(cea.ENTHALPY, weights, T_reactant) / cea.R
solver.solve(solution, weights, pc, pi_p=pi_p, supar=supar, hc=hc, iac=True)
# see below for the printing stuff

The output of cea binary

THEORETICAL ROCKET PERFORMANCE ASSUMING EQUILIBRIUM
COMPOSITION DURING EXPANSION FROM INFINITE AREA COMBUSTOR
 
Pc = 1.01 bar
 
pi/p = 1.0000
supar = 4.0000
 
 
 REACTANT MOLES ENERGY TEMP
 KJ/MOL K
 Acrylonitrile 39.920 140000.000 298.15
 Butadiene 47.370 90500.000 298.15
 Styrene 12.710 103400.000 298.15
 N2O 100.000 81600.000 298.15
 
O/F = 5.09859 % Fuel = 16.39723 r, Eq. Ratio = 0.00000 phi, Eq. Ratio = 0.00000
 
 CHAMBER THROAT EXIT EXIT
 Pinf/P 1.0000 1.7327 1.000 21.448
 P, bar 1.0132 0.5848 1.0132 0.0472
 T, K 3054.12 2899.49 3054.12 2138.78
 Density, kg/m^3 0.9993E-1 0.6161E-1 0.9993E-1 0.7041E-2
 H, kJ/kg 1873.352 1334.037 1873.354 -707.805
 U, kJ/kg 859.344 384.915 859.345 -1378.776
 G, kJ/kg -32116.28 -30934.65 -32116.29 -24510.47
 S, kJ/kg-K 11.129 11.129 11.129 11.129
 
 M, (1/n) 25.0427 25.4001 25.0427 26.5032
 (dln(V)/dln(P))t -1.0299 -1.0241 -1.0299 -1.0016
 (dln(V)/dln(T))p 1.6003 1.5139 1.6003 1.0431
 Cp, kJ/(kg-K) 5.6400 5.2020 5.6400 1.9123
 Gamma_s 1.1375 1.1365 1.1375 1.2150
 Son. Vel., m/s 1073.99 1038.57 1073.99 902.89
 Mach 0.0000 1.0000 NaN 2.5164
 
 PERFORMANCE PARAMETERS
 
 Ae/At 1.0000 NaN 4.0000
 C*, m/s 1583.45 1583.45 1583.45
 Cf 0.6559 NaN 1.4349
 Ivac, m/s 1952.44 NaN 2567.38
 Isp, m/s 1038.57 NaN 2272.07
 
MOLE FRACTIONS
 
 CO 0.226498 0.222595 0.226498 0.206825
 CO2 0.050113 0.057965 0.050113 0.085919
 H 0.042405 0.034232 0.042405 0.004395
 H2 0.052549 0.051875 0.052549 0.058441
 H2O 0.087324 0.097099 0.087324 0.123134
 N 1.359E-5 6.549E-6 1.359E-5 1.947E-8
 N2 0.487266 0.495331 0.487266 0.519707
 NO 0.007722 0.005618 0.007722 1.445E-4
 O 0.011919 0.007993 0.011919 6.210E-5
 O2 0.007806 0.005955 0.007806 5.550E-5
 OH 0.026376 0.021326 0.026376 0.001316

whereas the output of cea Python package

P, bar 1.013 0.584 1.013 0.046
T, K 2994.246 2836.530 2994.246 2000.870
Density, kg/m^3 0.103 0.064 0.103 0.007
H, kJ/kg 1550.02 1026.02 1550.02 -953.04
U, kJ/kg 566.87 106.96 566.87 -1579.59
G, kJ/kg -31453.1 -30238.8 -31453.1 -23007.0
S, kJ/kg-K 11.022 11.022 11.022 11.022
M, (1/n) 25.322 25.661 25.322 26.552
MW 25.322 25.661 25.322 26.552
Cp_eq, kJ/kg-K 5.156 4.604 5.156 1.684
Cp_fr, kJ/kg-K 1.566 1.560 1.566 1.506
Cv_eq, kJ/kg-K 4.415 3.961 4.415 1.360
Cv_eq, kJ/kg-K 1.237 1.236 1.237 1.192
Gamma_s 1.139 1.140 1.139 1.237
Son. vel., m/s 1058.25 1023.72 1058.25 880.53
Mach 0.000 1.000 nan 2.541
PERFORMANCE PARAMETERS
Ae/At 0.000 1.000 nan 4.000
C*, m/s 1557.81 1557.81 1557.81 1557.81
Cf 0.0000 0.6572 nan 1.4363
Isp, vac., m/s 0.000 1921.485 nan 2517.464
Isp, m/s 0.000 1023.721 nan 2237.435
MOLE FRACTIONS
CO 0.22434 0.22018 0.22434 0.20474
CO2 0.055364 0.063269 0.055364 0.088546
H 0.03523 0.027642 0.03523 0.0019052
H2 0.052091 0.051538 0.052091 0.061562
H2O 0.094997 0.10437 0.094997 0.12204
N 9.3531e-06 4.2282e-06 9.3531e-06 3.0869e-09
N2 0.49336 0.50104 0.49336 0.52072
NO 0.0065091 0.0044611 0.0065091 3.67e-05
O 0.0087627 0.0054229 0.0087627 8.5008e-06
O2 0.0063217 0.0043898 0.0063217 7.2603e-06
OH 0.023015 0.01769 0.023015 0.00042668

# printing stuff
num_pts = solution.num_pts
T = solution.T
P = solution.P
rho = solution.density
enthalpy = solution.enthalpy
energy = solution.energy
gibbs = solution.gibbs_energy
entropy = solution.entropy
M_1n = solution.M
MW = solution.MW
cp_eq = solution.cp_eq
cp_fr = solution.cp_fr
cv_eq = solution.cv_eq
cv_fr = solution.cv_fr
Mach = solution.Mach
gamma_s = solution.gamma_s
v_sonic = solution.sonic_velocity
ae_at = solution.ae_at
c_star = solution.c_star
Cf = solution.coefficient_of_thrust
Isp = solution.Isp
Isp_vac = solution.Isp_vacuum
print("P, bar ", end=" ")
for i in range(num_pts):
 if i < num_pts - 1:
 print("{0:10.3f}".format(P[i]), end=" ")
 else:
 print("{0:10.3f}".format(P[i]))
print("T, K ", end=" ")
for i in range(num_pts):
 if i < num_pts - 1:
 print("{0:10.3f}".format(T[i]), end=" ")
 else:
 print("{0:10.3f}".format(T[i]))
print("Density, kg/m^3", end=" ")
for i in range(num_pts):
 if i < num_pts - 1:
 print("{0:10.3f}".format(rho[i]), end=" ")
 else:
 print("{0:10.3f}".format(rho[i]))
print("H, kJ/kg ", end=" ")
for i in range(num_pts):
 if i < num_pts - 1:
 print("{0:10.2f}".format(enthalpy[i]), end=" ")
 else:
 print("{0:10.2f}".format(enthalpy[i]))
print("U, kJ/kg ", end=" ")
for i in range(num_pts):
 if i < num_pts - 1:
 print("{0:10.2f}".format(energy[i]), end=" ")
 else:
 print("{0:10.2f}".format(energy[i]))
print("G, kJ/kg ", end=" ")
for i in range(num_pts):
 if i < num_pts - 1:
 print("{0:10.1f}".format(gibbs[i]), end=" ")
 else:
 print("{0:10.1f}".format(gibbs[i]))
print("S, kJ/kg-K ", end=" ")
for i in range(num_pts):
 if i < num_pts - 1:
 print("{0:10.3f}".format(entropy[i]), end=" ")
 else:
 print("{0:10.3f}".format(entropy[i]))
print("M, (1/n) ", end=" ")
for i in range(num_pts):
 if i < num_pts - 1:
 print("{0:10.3f}".format(M_1n[i]), end=" ")
 else:
 print("{0:10.3f}".format(M_1n[i]))
print("MW ", end=" ")
for i in range(num_pts):
 if i < num_pts - 1:
 print("{0:10.3f}".format(MW[i]), end=" ")
 else:
 print("{0:10.3f}".format(MW[i]))
print("Cp_eq, kJ/kg-K ", end=" ")
for i in range(num_pts):
 if i < num_pts - 1:
 print("{0:10.3f}".format(cp_eq[i]), end=" ")
 else:
 print("{0:10.3f}".format(cp_eq[i]))
print("Cp_fr, kJ/kg-K ", end=" ")
for i in range(num_pts):
 if i < num_pts - 1:
 print("{0:10.3f}".format(cp_fr[i]), end=" ")
 else:
 print("{0:10.3f}".format(cp_fr[i]))
print("Cv_eq, kJ/kg-K ", end=" ")
for i in range(num_pts):
 if i < num_pts - 1:
 print("{0:10.3f}".format(cv_eq[i]), end=" ")
 else:
 print("{0:10.3f}".format(cv_eq[i]))
print("Cv_eq, kJ/kg-K ", end=" ")
for i in range(num_pts):
 if i < num_pts - 1:
 print("{0:10.3f}".format(cv_fr[i]), end=" ")
 else:
 print("{0:10.3f}".format(cv_fr[i]))
print("Gamma_s ", end=" ")
for i in range(num_pts):
 if i < num_pts - 1:
 print("{0:10.3f}".format(gamma_s[i]), end=" ")
 else:
 print("{0:10.3f}".format(gamma_s[i]))
print("Son. vel., m/s ", end=" ")
for i in range(num_pts):
 if i < num_pts - 1:
 print("{0:10.2f}".format(v_sonic[i]), end=" ")
 else:
 print("{0:10.2f}".format(v_sonic[i]))
print("Mach ", end=" ")
for i in range(num_pts):
 if i < num_pts - 1:
 print("{0:10.3f}".format(Mach[i]), end=" ")
 else:
 print("{0:10.3f}".format(Mach[i]))
print()
print("PERFORMANCE PARAMETERS")
print()
print("Ae/At ", end=" ")
for i in range(num_pts):
 if i < num_pts - 1:
 print("{0:10.3f}".format(ae_at[i]), end=" ")
 else:
 print("{0:10.3f}".format(ae_at[i]))
print("C*, m/s ", end=" ")
for i in range(num_pts):
 if i < num_pts - 1:
 print("{0:10.2f}".format(c_star[i]), end=" ")
 else:
 print("{0:10.2f}".format(c_star[i]))
print("Cf ", end=" ")
for i in range(num_pts):
 if i < num_pts - 1:
 print("{0:10.4f}".format(Cf[i]), end=" ")
 else:
 print("{0:10.4f}".format(Cf[i]))
print("Isp, vac., m/s ", end=" ")
for i in range(num_pts):
 if i < num_pts - 1:
 print("{0:10.3f}".format(Isp_vac[i]), end=" ")
 else:
 print("{0:10.3f}".format(Isp_vac[i]))
print("Isp, m/s ", end=" ")
for i in range(num_pts):
 if i < num_pts - 1:
 print("{0:10.3f}".format(Isp[i]), end=" ")
 else:
 print("{0:10.3f}".format(Isp[i]))
print()
print("MOLE FRACTIONS")
print("")
trace_species = []
for prod in solution.mole_fractions:
 if np.any(solution.mole_fractions[prod] > 5e-6):
 print("{0:15s}".format(prod), end=" ")
 for j in range(len(solution.mole_fractions[prod])):
 if j < len(solution.mole_fractions[prod]) - 1:
 print("{0:10.5g}".format(solution.mole_fractions[prod][j]), end=" ")
 else:
 print("{0:10.5g}".format(solution.mole_fractions[prod][j]))
 else:
 trace_species.append(prod)

Answered by markleader

Mar 23, 2026

Thanks for reporting this and for comparing the binary and Python results so carefully.

I took a closer look, and I do not think this is a solver divergence. The difference comes from the enthalpy units used in the Python API. Right now, the Python Reactant API interprets enthalpy using J/kg, while the .inp file used kJ/mol.
Changing the values in your Python input to match these default units should make the Python and binary results consistent today.

View full answer

Replies: 1 comment

markleader
Mar 23, 2026
Maintainer

Thanks for reporting this and for comparing the binary and Python results so carefully.

I am also working on improving this in the Python API so the units are more explicit. The plan is to add an enthalpy_units argument going forward. Initially this will be optional, but recommended, so existing code keeps working. After that transition, I expect to require users to pass enthalpy units explicitly rather than relying on defaults, which should avoid this kind of confusion entirely.

Thanks again for flagging it. This was a helpful catch.

0 replies

Answer selected by markleader

Navigation Menu

Search code, repositories, users, issues, pull requests...

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Different output from running Python script and cea binary directly #63

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{editor}}'s edit

{{editor}}'s edit

Uh oh!

minh-dng
Mar 21, 2026

Replies: 1 comment

Uh oh!

{{title}}

Uh oh!

markleader
Mar 23, 2026
Maintainer

Select a reply

Uh oh!

Uh oh!

Different output from running Python script and cea binary directly #63

Uh oh!

Uh oh!

minh-dng Mar 21, 2026

Replies: 1 comment

Uh oh!

markleader Mar 23, 2026 Maintainer

minh-dng
Mar 21, 2026

markleader
Mar 23, 2026
Maintainer