cea_cache.py

Go to the documentation of this file.

"""RocketCEA wrapper with MR/Pc caching and bilinear interpolation"""
 
import numpy as np
import os
import re
import json
import sys
import threading
import logging
from typing import Tuple, Optional, List, Dict, Any
from multiprocessing import Pool, cpu_count, Manager
from functools import partial
import time
from .config_schemas import CEAConfig
 
# ── In-process CEA-build gate ────────────────────────────────────────────────
# A cold cache build is ~n_points^3 (≈39k) rocketcea calls: hours of CPU and a
# steadily growing in-process cache that has taken the host to >12 GB when it
# runs inside a request handler. Building is an OFFLINE job (see
# .github/workflows/regenerate-cea-cache.yml); the committed output/cache/*.npz
# tables are what everyone loads.
#
# The gate is ENVIRONMENT-based, not code-path based: running the app locally
# over localhost is the same backend.main as production, just a different
# environment, so we can't distinguish them by which code is executing. Instead
# building defaults to ALLOWED, and the *production* environment opts out by
# setting ED_ALLOW_CEA_BUILD=0 (see docker-compose.yml). That way local dev, the
# regenerate workflow, tests, and direct runner/script use all keep building on
# demand, while the deployed server fails fast on a missing/stale cache instead
# of wedging the box. Set ED_ALLOW_CEA_BUILD=1 to force a build anywhere.
def _cea_build_allowed() -> bool:
    env = os.environ.get("ED_ALLOW_CEA_BUILD")
    if env is None:
        return True  # default: local/CLI/CI/tests may build on demand
    return env.strip().lower() in ("1", "true", "yes", "on")
 
 
# Version of the SEMANTICS of the cached tables (not the .npz container format).
# Bump whenever a table's meaning changes or a table is added, so stale caches on
# disk — including locally generated ones git never touches — get rebuilt instead
# of being silently served to code that expects the new meaning.
#   1 (implicit; pre-2026-07): "Cf" stored get_PambCf()[0]; no "Cf_vac" table.
#   2: "Cf" is the ambient-corrected coefficient get_PambCf()[1] (reproduces CEA's
#      own ambient Isp) and the "Cf_vac" table (RPA delivered-thrust basis) exists.
# Caches whose meta lacks "table_schema" are treated as version 1.
CEA_TABLE_SCHEMA_VERSION = 2
 
 
def _isentropic_cf_vac(gamma: float, eps: float) -> float:
    """Ideal vacuum thrust coefficient from gamma + area ratio (isentropic, frozen).
 
    Fallback for caches built before the Cf_vac table existed. Slightly low vs CEA's
    shifting-equilibrium Cf_vac (the same ~2% frozen gap documented in the thrust-model
    analysis); regenerate the cache to get the exact CEA value. Self-contained to avoid
    importing the nozzle solvers here.
    """
    g = float(gamma)
    e = float(eps)
    if not (g > 1.0) or not (e > 1.0):
        return float("nan")
    gm1, gp1 = g - 1.0, g + 1.0
    # Supersonic area-Mach solve (Newton), same relation as mach_solver.
    M = 2.5
    for _ in range(100):
        t = (2.0 / gp1) * (1.0 + 0.5 * gm1 * M * M)
        AR = (1.0 / M) * t ** (gp1 / (2.0 * gm1))
        dM = 1e-6
        t2 = (2.0 / gp1) * (1.0 + 0.5 * gm1 * (M + dM) ** 2)
        AR2 = (1.0 / (M + dM)) * t2 ** (gp1 / (2.0 * gm1))
        d = (AR2 - AR) / dM
        if d == 0.0:
            break
        Mn = M - (AR - e) / d
        if Mn <= 1.0:
            Mn = 0.5 * (M + 1.0)
        if abs(Mn - M) < 1e-10:
            M = Mn
            break
        M = Mn
    pe_pc = (1.0 + 0.5 * gm1 * M * M) ** (-g / gm1)
    cf = ((2.0 * g * g / gm1) * (2.0 / gp1) ** (gp1 / gm1) * (1.0 - pe_pc ** (gm1 / g))) ** 0.5 + pe_pc * e
    return float(cf)
 
 
def _scalar_clamp(x: float, lo: float, hi: float) -> float:
    """float(np.clip(x, lo, hi)) for a scalar, without numpy's per-call overhead.
    Matches np.clip exactly, including NaN pass-through (NaN < lo and NaN > hi are
    both False, so NaN is returned). Hot path: called per CEA lookup."""
    if x < lo:
        return lo
    if x > hi:
        return hi
    return x
 
 
# Lazy import of rocketcea - only import if we need to build cache
# This allows the module to work even if rocketcea is not installed/available
# as long as a cache file exists
_CEA_Obj = None
 
def _get_CEA_Obj():
    """Lazy import of CEA_Obj - only import when needed for cache building."""
    global _CEA_Obj
    if _CEA_Obj is None:
        try:
            from rocketcea.cea_obj import CEA_Obj
            _CEA_Obj = CEA_Obj
        except ImportError as e:
            raise ImportError(
                "rocketcea is required to build CEA cache. "
                "If you have an existing cache file, make sure it's in the correct location. "
                f"Original error: {e}"
            )
    return _CEA_Obj
 
# Create module-level logger
logger = logging.getLogger("evaluate")
 
# Fix console encoding issues (Windows, WSL, etc.)
_original_print = print
def safe_print(*args, **kwargs):
    """Log function that handles Unicode encoding errors gracefully"""
    try:
        # Build message string
        message = ' '.join(str(arg) for arg in args)
        logger.info(message)
    except UnicodeEncodeError:
        # Replace problematic characters with ASCII equivalents
        safe_args = []
        for arg in args:
            if isinstance(arg, str):
                # Replace Unicode characters that can't be encoded
                safe_str = arg.encode('ascii', errors='replace').decode('ascii')
                safe_args.append(safe_str)
            else:
                safe_args.append(arg)
        try:
            message = ' '.join(str(arg) for arg in safe_args)
            logger.info(message)
        except Exception:
            # Last resort: convert everything to string and sanitize
            safe_str_args = []
            for arg in safe_args:
                safe_str_args.append(str(arg).encode('ascii', errors='replace').decode('ascii'))
            message = ' '.join(safe_str_args)
            logger.info(message)
 
# Replace built-in print with safe version
print = safe_print
 
# Also try to reconfigure stdout if possible
try:
    if hasattr(sys.stdout, 'reconfigure'):
        sys.stdout.reconfigure(encoding='utf-8', errors='replace')
    elif hasattr(sys.stdout, 'buffer'):
        import io
        sys.stdout = io.TextIOWrapper(sys.stdout.buffer, encoding='utf-8', errors='replace')
except (AttributeError, ValueError, OSError):
    pass  # Safe print wrapper will handle it
 
 
def parse_cea_basic(out: str) -> Tuple[float, float, float, float, float, float]:
    """
    Extract Tc, gamma, R, c*, and Isp from NASA CEA text output (chamber values only).
    
    Returns:
    --------
    Tc : float [K]
    gamma : float
    R : float [J/(kg*K)]
    cstar : float [m/s]
    M : float [kg/kmol]
    Isp : float [s] (ambient Isp, or NaN if not found)
    """
    # Extract the main performance block
    block_match = re.search(r'THEORETICAL ROCKET PERFORMANCE[\s\S]+?MOLE FRACTIONS', out)
    block = block_match.group(0) if block_match else out
 
    # Extract chamber molecular weight M
    mw_match = re.search(r'M,\s*\‍(1/n\‍)\s+([\d.E+-]+)', block)
    if mw_match:
        M = float(mw_match.group(1))  # molecular weight [kg/kmol]
        R = 8314.462618 / M  # J/(kg*K)
    else:
        M = np.nan
        R = np.nan
 
    # Extract chamber temperature (Tc)
    Tc_match = re.search(r'T[, ]*K\s+([\d.E+-]+)', block)
    Tc = float(Tc_match.group(1)) if Tc_match else np.nan
 
    # Extract gamma (chamber)
    gamma_match = re.search(r'GAMMAs\s+([\d.E+-]+)', block)
    gamma = float(gamma_match.group(1)) if gamma_match else np.nan
 
    # Extract c*
    cstar_match = re.search(r'CSTAR[, ]*(?:M/SEC|FT/SEC)?\s+([\d.E+-]+)', block, re.IGNORECASE)
    cstar = float(cstar_match.group(1)) if cstar_match else np.nan
    # Convert ft/s -> m/s if needed
    if re.search(r'CSTAR.*FT/SEC', block, re.IGNORECASE):
        cstar *= 0.3048
 
    # Extract Isp (ambient Isp, typically appears as "Isp, M/SEC" or "Isp, SEC")
    # Look for Isp in the performance block - it may appear in different formats
    isp_match = re.search(r'Isp[, ]*(?:M/SEC|SEC|M/S)?\s+([\d.E+-]+)', block, re.IGNORECASE)
    if not isp_match:
        # Try alternative pattern: look for "Isp" followed by numbers
        isp_match = re.search(r'Isp[:\s]+([\d.E+-]+)', block, re.IGNORECASE)
    Isp = float(isp_match.group(1)) if isp_match else np.nan
 
    return Tc, gamma, R, cstar, M, Isp
 
 
def _compute_cea_point_chunk(
    chunk: List[Tuple[int, int, int, float, float, float]],
    ox_name: str,
    fuel_name: str,
    expansion_ratio: Optional[float] = None,
    lock: Optional[Any] = None,
) -> List[Tuple[int, int, int, float, float, float, float, float, float]]:
    """
    Worker function for parallel CEA cache building.
    
    Computes CEA properties for a chunk of grid points.
    This function must be at module level for multiprocessing.
    
    Parameters:
    -----------
    chunk : List[Tuple[int, int, int, float, float, float]]
        List of (i, j, k_idx, Pc_psia, MR, eps) tuples
    ox_name : str
        Oxidizer name
    fuel_name : str
        Fuel name
    expansion_ratio : float, optional
        Fixed expansion ratio (for 2D cache). If None, uses eps from chunk (3D cache)
    lock : multiprocessing.Lock, optional
        Lock to serialize RocketCEA calls (RocketCEA is not process-safe)
    
    Returns:
    --------
    results : List[Tuple[int, int, int, float, float, float, float, float, float]]
        List of (i, j, k_idx, cstar, Cf, Tc, gamma, R, M) tuples
    """
    # Lazy import CEA_Obj only when needed
    CEA_Obj = _get_CEA_Obj()
    # Create CEA object per worker (each process gets its own)
    chamber = CEA_Obj(oxName=ox_name, fuelName=fuel_name)
    results = []
    
    for i, j, k_idx, Pc_psia, MR, eps in chunk:
        # Use eps from chunk if provided (3D), otherwise use fixed expansion_ratio (2D)
        eps_to_use = eps if expansion_ratio is None else expansion_ratio
        
        try:
            # Serialize RocketCEA calls if lock is provided (RocketCEA uses shared temp files)
            if lock is not None:
                lock.acquire()
            try:
                out = chamber.get_full_cea_output(Pc=Pc_psia, MR=MR, eps=eps_to_use)
                Tc, gamma, R, cstar, M, Isp = parse_cea_basic(out)
 
                try:
                    # get_PambCf returns (CFcea, CFamb, mode); element [1] is the
                    # ambient-corrected Cf that reproduces CEA's own ambient Isp.
                    Cf_ideal = chamber.get_PambCf(Pc=Pc_psia, MR=MR, eps=eps_to_use)[1]
                except Exception:
                    # Fallback: estimate from Isp
                    try:
                        isp = chamber.estimate_Ambient_Isp(Pc=Pc_psia, MR=MR, eps=eps_to_use)[0]
                        Cf_ideal = isp * 9.80665 / cstar if cstar > 0 else np.nan
                    except Exception:
                        Cf_ideal = np.nan
                # Vacuum thrust coefficient (shifting equilibrium): Cf_vac = Isp_vac*g0/c*.
                # This is the RPA delivered-thrust basis: F = zeta_n*Cf_vac*Pc*At - Pa*Ae.
                try:
                    isp_vac = chamber.get_Isp(Pc=Pc_psia, MR=MR, eps=eps_to_use)
                    Cf_vac = isp_vac * 9.80665 / cstar if cstar > 0 else np.nan
                except Exception:
                    Cf_vac = np.nan
 
                results.append((i, j, k_idx, cstar, Cf_ideal, Tc, gamma, R, M, Cf_vac))
            finally:
                if lock is not None:
                    lock.release()
        except Exception as e:
            # On failure (including Fortran I/O errors), store NaN values
            # This handles "I/O past end of record" and "End of file" errors
            results.append((i, j, k_idx, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan))
 
    return results
 
 
# Decompressed CEA tables, keyed by (path, mtime, size). np.load returns a lazy
# NpzFile: every table access re-inflates that member from the zip, and a fresh
# CEACache is built for every PintleEngineRunner. Layer 1 builds a runner per parent
# objective() call, so one short run re-read this file 73 times -- 0.80 s of a 3.22 s
# parent, against 0.11 s for all the interpolation it actually fed. The tables are
# written only in _load_cache/_build_cache and never mutated afterwards (no caller
# assigns to a cache's attributes), so the arrays are safe to share by reference.
# Keying on mtime+size means a regenerated .npz is picked up rather than served stale.
_NPZ_TABLE_MEMO: Dict[Tuple[str, int, int], Dict[str, np.ndarray]] = {}
 
 
def _load_npz_tables(path: str) -> Dict[str, np.ndarray]:
    """Return the fully-materialised contents of a cache .npz, memoised per file version."""
    try:
        st = os.stat(path)
        key = (os.path.abspath(path), st.st_mtime_ns, st.st_size)
    except OSError:
        key = None
    if key is not None:
        hit = _NPZ_TABLE_MEMO.get(key)
        if hit is not None:
            return hit
    with np.load(path) as z:
        tables = {name: z[name] for name in z.files}
    if key is not None:
        # Bounded: one entry per propellant cache file, and they are a few hundred KB each.
        if len(_NPZ_TABLE_MEMO) > 8:
            _NPZ_TABLE_MEMO.clear()
        _NPZ_TABLE_MEMO[key] = tables
    return tables
 
 
class CEACache:
    """CEA cache with bilinear interpolation"""
    
    def __init__(self, config: CEAConfig):
        self.config = config
        # Resolve a deterministic cache path for this config.
        # Important: never silently swap to a different *.npz file, which can hide
        # fuel/range changes and make runs appear inconsistent.
        if os.path.isabs(config.cache_file):
            self.cache_file = config.cache_file
        else:
            project_root = os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
            candidates = [
                os.path.abspath(config.cache_file),
                os.path.join(project_root, config.cache_file),
            ]
            resolved = None
            for cand in candidates:
                if os.path.exists(cand):
                    resolved = cand
                    break
            self.cache_file = resolved if resolved is not None else candidates[1]
        
        # Determine if using 3D cache (Pc, MR, eps) or 2D cache (Pc, MR)
        self.use_3d = config.eps_range is not None
        
        # Grid parameters
        self.Pc_min, self.Pc_max = config.Pc_range
        self.MR_min, self.MR_max = config.MR_range
        self.n_points = config.n_points
        
        # Create grids
        self.Pc_grid = np.linspace(self.Pc_min, self.Pc_max, self.n_points)
        self.MR_grid = np.linspace(self.MR_min, self.MR_max, self.n_points)
        
        # 3D mode: add expansion ratio grid
        if self.use_3d:
            self.eps_min, self.eps_max = config.eps_range
            self.eps_grid = np.linspace(self.eps_min, self.eps_max, self.n_points)
            # print(f"[INFO] Using 3D CEA cache: {self.n_points}^3 = {self.n_points**3} points")
        else:
            self.eps_min = self.eps_max = config.expansion_ratio
            self.eps_grid = None
            # print(f"[INFO] Using 2D CEA cache: {self.n_points}^2 = {self.n_points**2} points")
        
        # Lookup tables (initialized as None, loaded from cache or computed)
        # Shape: (n, n, n) for 3D or (n, n) for 2D
        self.cstar_table = None
        self.Cf_table = None
        self.Tc_table = None
        self.gamma_table = None
        self.R_table = None
        self.M_table = None
        self.Cf_vac_table = None  # vacuum thrust coefficient (RPA delivered-thrust basis)
        
        # Load from cache or build
        if os.path.exists(self.cache_file):
            self._load_cache()
        else:
            self._build_cache()
    
    def _load_cache(self):
        """Load CEA data from cache file"""
        # print(f"[OK] Loading CEA cache from {self.cache_file}")
        data = _load_npz_tables(self.cache_file)
 
        meta_expected = {
            "table_schema": CEA_TABLE_SCHEMA_VERSION,
            "ox_name": self.config.ox_name,
            "fuel_name": self.config.fuel_name,
            "expansion_ratio": self.config.expansion_ratio,
            "Pc_range": list(self.config.Pc_range),
            "MR_range": list(self.config.MR_range),
            "n_points": self.n_points,
            "dimensions": 3 if self.use_3d else 2,
            "eps_range": list(self.config.eps_range) if self.use_3d else None,
        }
 
        meta_loaded = None
        if "meta" in data:
            try:
                meta_loaded = json.loads(data["meta"].tolist())
            except Exception:
                meta_loaded = None
 
        def _meta_matches(meta_loaded_dict: Optional[dict], meta_expected_dict: dict) -> bool:
            """Check if cache is usable - lenient matching: only check critical fields"""
            if meta_loaded_dict is None:
                return False
            try:
                # Critical: table SEMANTICS version must match. A stale-schema cache
                # (e.g. pre-Cf_vac, or "Cf" stored on the old get_PambCf[0] basis)
                # must be rebuilt, not silently served under the new meaning —
                # git updates code but never touches locally generated caches.
                loaded_schema = meta_loaded_dict.get("table_schema", 1)
                if loaded_schema != meta_expected_dict["table_schema"]:
                    logger.warning(
                        f"CEA cache table-schema mismatch: cache has v{loaded_schema}, "
                        f"code expects v{meta_expected_dict['table_schema']}; regenerating."
                    )
                    return False
                # Critical: propellant names must match
                if meta_loaded_dict.get("ox_name") != meta_expected_dict["ox_name"]:
                    return False
                if meta_loaded_dict.get("fuel_name") != meta_expected_dict["fuel_name"]:
                    return False
                # Critical: dimensions must match (2D vs 3D)
                loaded_dims = meta_loaded_dict.get("dimensions", 2)
                expected_dims = meta_expected_dict["dimensions"]
                if loaded_dims != expected_dims:
                    logger.warning(f"Cache dimension mismatch: {loaded_dims}D vs {expected_dims}D")
                    return False
                # Non-critical: ranges and n_points can differ - we can still interpolate
                # Just warn if they're very different
                loaded_pc_range = meta_loaded_dict.get("Pc_range", [])
                expected_pc_range = meta_expected_dict.get("Pc_range", [])
                if loaded_pc_range and expected_pc_range:
                    # Check if requested range is within cached range (with some margin)
                    if expected_pc_range[0] < loaded_pc_range[0] * 0.9 or expected_pc_range[1] > loaded_pc_range[1] * 1.1:
                        logger.info(f"Requested Pc range {expected_pc_range} extends beyond cached range {loaded_pc_range}")
                        logger.info(f"Will use trilinear interpolation - some extrapolation may occur")
                
                loaded_mr_range = meta_loaded_dict.get("MR_range", [])
                expected_mr_range = meta_expected_dict.get("MR_range", [])
                if loaded_mr_range and expected_mr_range:
                    if expected_mr_range[0] < loaded_mr_range[0] * 0.9 or expected_mr_range[1] > loaded_mr_range[1] * 1.1:
                        logger.info(f"Requested MR range {expected_mr_range} extends beyond cached range {loaded_mr_range}")
                        logger.info(f"Will use trilinear interpolation - some extrapolation may occur")
                
                # For 3D caches, check eps_range similarly
                if expected_dims == 3:
                    loaded_eps_range = meta_loaded_dict.get("eps_range", [])
                    expected_eps_range = meta_expected_dict.get("eps_range", [])
                    if loaded_eps_range and expected_eps_range:
                        if expected_eps_range[0] < loaded_eps_range[0] * 0.9 or expected_eps_range[1] > loaded_eps_range[1] * 1.1:
                            logger.info(f"Requested eps range {expected_eps_range} extends beyond cached range {loaded_eps_range}")
                            logger.info(f"Will use trilinear interpolation - some extrapolation may occur")
            except Exception as e:
                logger.warning(f"Error checking cache metadata: {e}")
                return False
            return True
 
        if not _meta_matches(meta_loaded, meta_expected):
            print("[WARNING] CEA cache metadata does not match configuration; regenerating...")
            try:
                os.remove(self.cache_file)
            except OSError:
                pass
            self._build_cache()
            return
        
        self.cstar_table = data["cstar"]
        self.Cf_table = data["Cf"]
        # Cf_vac added later; old caches lack it -> None triggers an isentropic
        # fallback in eval() until the cache is regenerated.
        self.Cf_vac_table = data["Cf_vac"] if "Cf_vac" in data else None
        self.Tc_table = data["Tc"]
        self.gamma_table = data["gamma"]
        self.R_table = data["R"]
        if "M" in data:
            self.M_table = data["M"]
        else:
            # Backwards compatibility: derive M from R
            self.M_table = 8314.462618 / self.R_table
        
        # Load grids from cache (use cached grids for interpolation)
        # Trilinear interpolation can handle different grid ranges - no need to rebuild
        Pc_loaded = data["Pc"]
        MR_loaded = data["MR"]
        
        # Use cached grids instead of regenerating - trilinear interpolation can handle range differences
        # Only rebuild if grids are completely incompatible (different sizes)
        if len(Pc_loaded) != len(self.Pc_grid) or len(MR_loaded) != len(self.MR_grid):
            print("[WARNING] Cache grid size doesn't match config, rebuilding...")
            self._build_cache()
            return
        
        # For 3D caches, check eps grid
        if self.use_3d:
            if "eps" in data:
                eps_loaded = data["eps"]
                if len(eps_loaded) != len(self.eps_grid):
                    print("[WARNING] Cache eps grid size doesn't match config, rebuilding...")
                    self._build_cache()
                    return
                # Use cached eps grid
                self.eps_grid = eps_loaded
            else:
                print("[WARNING] 3D cache missing eps grid, rebuilding...")
                self._build_cache()
                return
        
        # Use cached grids for interpolation (they may have different ranges, that's OK)
        # Update our grid references to use cached grids
        self.Pc_grid = Pc_loaded
        self.MR_grid = MR_loaded
        # Update min/max from loaded grids
        self.Pc_min = float(Pc_loaded.min())
        self.Pc_max = float(Pc_loaded.max())
        self.MR_min = float(MR_loaded.min())
        self.MR_max = float(MR_loaded.max())
        if self.use_3d:
            self.eps_min = float(self.eps_grid.min())
            self.eps_max = float(self.eps_grid.max())
        print(f"[INFO] Using cached grids - trilinear interpolation will handle any range differences")
    
    def _build_cache(self):
        """Build CEA lookup tables with parallel processing support"""
        # Fail fast where in-process builds are disabled (the web server). Every
        # build path -- missing file, metadata mismatch, schema bump -- funnels
        # through here, so this one check covers them all. A build in a request
        # handler is a multi-hour, memory-runaway hazard (host hit >12 GB); the
        # cure is a committed cache, not an on-demand build.
        if not _cea_build_allowed():
            raise RuntimeError(
                f"CEA cache missing or stale for {self.config.ox_name}/{self.config.fuel_name} "
                f"at {self.cache_file}, and in-process builds are disabled here. A cold build is "
                f"~{self.n_points}^3 rocketcea calls (hours) and can exhaust memory. Regenerate it "
                f"offline (.github/workflows/regenerate-cea-cache.yml) and commit the .npz, or set "
                f"ED_ALLOW_CEA_BUILD=1 to force a build in this process."
            )
        print(f"[BUILDING] Building CEA cache (this will take a while)...")
        print(f"   Grid: Pc in [{self.Pc_min/1e6:.1f}, {self.Pc_max/1e6:.1f}] MPa ({self.n_points} points)")
        print(f"         MR in [{self.MR_min:.2f}, {self.MR_max:.2f}] ({self.n_points} points)")
        if self.use_3d:
            print(f"         eps in [{self.eps_min:.2f}, {self.eps_max:.2f}] ({self.n_points} points)")
            total_points = self.n_points**3
            print(f"         Total: {self.n_points}^3 = {total_points} points")
            print(f"   [WARNING] Sequential build would take ~{total_points * 0.5 / 60:.0f} minutes!")
        else:
            total_points = self.n_points**2
            print(f"         Total: {self.n_points}^2 = {total_points} points")
            print(f"   [INFO] Sequential build would take ~{total_points * 0.5 / 60:.0f} minutes")
        
        # Check for parallel processing option
        # RocketCEA (Fortran) is NOT process-safe - multiple processes collide on temp files
        # Parallel processing causes Fortran I/O errors and hangs
        # Disable parallel by default - sequential is slower but reliable
        is_windows = sys.platform == 'win32'
        use_parallel = getattr(self.config, 'use_parallel_cea_build', False)  # Default to False
        n_workers = getattr(self.config, 'cea_parallel_workers', None)
        
        # Only parallelize if explicitly enabled and for large grids
        should_use_parallel = use_parallel and total_points > 100
 
        # macOS / Python: multiprocessing.Manager()/Pool started from a *non-main* thread
        # (e.g. FastAPI Layer‑1 optimizer via ThreadPoolExecutor) often raises
        # ``OSError: [Errno 1] Operation not permitted``. Sequential build avoids that.
        if (
            should_use_parallel
            and sys.platform == "darwin"
            and threading.current_thread() is not threading.main_thread()
        ):
            print(
                "   [WARNING] macOS + background thread: disabling parallel CEA build "
                "(avoids multiprocessing OSError errno 1); using sequential cache build."
            )
            should_use_parallel = False
        
        if should_use_parallel and not is_windows:
            # Linux/Mac: use parallel processing with lock (serialized RocketCEA calls)
            if n_workers is None:
                n_workers = min(cpu_count(), 4)  # Limit to 4 workers (RocketCEA is not safe)
            print(f"   [PARALLEL] Using {n_workers} workers with serialized RocketCEA calls")
            print(f"   [WARNING] RocketCEA is not process-safe - using lock to serialize calls")
            print(f"   [WARNING] Parallel mode may still hang - sequential is recommended")
            try:
                self._build_cache_parallel(n_workers, total_points)
            except (NameError, ImportError, MemoryError, Exception) as e:
                print(f"   [ERROR] Parallel processing failed: {e}")
                print(f"   [FALLBACK] Building cache sequentially")
                self._build_cache_sequential(total_points)
        elif should_use_parallel and is_windows:
            # Windows: Default to sequential due to multiprocessing spawn issues
            print(f"   [WARNING] Windows detected - using sequential mode to avoid multiprocessing issues")
            print(f"   [INFO] Sequential mode is slower but stable (ETA: ~{total_points * 0.5 / 60:.1f} min)")
            self._build_cache_sequential(total_points)
        else:
            print(f"   [SEQUENTIAL] Building cache sequentially (recommended for RocketCEA)")
            print(f"   [INFO] Sequential mode avoids Fortran I/O conflicts (ETA: ~{total_points * 0.5 / 60:.1f} min)")
            self._build_cache_sequential(total_points)
        
        # Save to cache
        self._save_cache()
        print(f"[SAVED] CEA cache saved to {self.cache_file}")
    
    def _build_cache_sequential(self, total_points: int):
        """Sequential CEA cache building (original method)"""
        # Lazy import CEA_Obj only when needed
        CEA_Obj = _get_CEA_Obj()
        chamber = CEA_Obj(oxName=self.config.ox_name, fuelName=self.config.fuel_name)
        
        # Initialize tables
        shape = (self.n_points, self.n_points, self.n_points) if self.use_3d else (self.n_points, self.n_points)
        self.cstar_table = np.zeros(shape)
        self.Cf_table = np.zeros(shape)
        self.Tc_table = np.zeros(shape)
        self.gamma_table = np.zeros(shape)
        self.R_table = np.zeros(shape)
        self.M_table = np.zeros(shape)
        self.Cf_vac_table = np.zeros(shape)
        
        Pc_psia_grid = self.Pc_grid / 6894.76
        point_count = 0
        start_time = time.time()
        
        for i, Pc_psia in enumerate(Pc_psia_grid):
            for j, MR in enumerate(self.MR_grid):
                eps_list = self.eps_grid if self.use_3d else [self.config.expansion_ratio]
                
                for k_idx, eps in enumerate(eps_list):
                    point_count += 1
                    elapsed = time.time() - start_time
                    if point_count % max(10, total_points // 100) == 0 or point_count == 1:
                        pct = 100 * point_count / total_points
                        rate = point_count / elapsed if elapsed > 0 else 0
                        eta = (total_points - point_count) / rate if rate > 0 else 0
                        print(f"   [{point_count}/{total_points}] {pct:.1f}% | Rate: {rate:.1f} pts/s | ETA: {eta/60:.1f} min | Pc={Pc_psia:.0f} psi, MR={MR:.2f}, eps={eps:.1f}")
                    
                    try:
                        out = chamber.get_full_cea_output(Pc=Pc_psia, MR=MR, eps=eps)
                        Tc, gamma, R, cstar, M, Isp = parse_cea_basic(out)
                        
                        try:
                            Cf_ideal = chamber.get_PambCf(Pc=Pc_psia, MR=MR, eps=eps)[1]
                        except:
                            isp = chamber.estimate_Ambient_Isp(Pc=Pc_psia, MR=MR, eps=eps)[0]
                            Cf_ideal = isp * 9.80665 / cstar if cstar > 0 else np.nan
                        try:
                            Cf_vac = chamber.get_Isp(Pc=Pc_psia, MR=MR, eps=eps) * 9.80665 / cstar if cstar > 0 else np.nan
                        except Exception:
                            Cf_vac = np.nan
 
                        if self.use_3d:
                            self.cstar_table[i, j, k_idx] = cstar
                            self.Cf_table[i, j, k_idx] = Cf_ideal
                            self.Cf_vac_table[i, j, k_idx] = Cf_vac
                            self.Tc_table[i, j, k_idx] = Tc
                            self.gamma_table[i, j, k_idx] = gamma
                            self.R_table[i, j, k_idx] = R
                            self.M_table[i, j, k_idx] = M
                        else:
                            self.cstar_table[i, j] = cstar
                            self.Cf_table[i, j] = Cf_ideal
                            self.Cf_vac_table[i, j] = Cf_vac
                            self.Tc_table[i, j] = Tc
                            self.gamma_table[i, j] = gamma
                            self.R_table[i, j] = R
                            self.M_table[i, j] = M
                        
                    except Exception as e:
                        print(f"   [WARNING] Error at Pc={Pc_psia:.1f} psia, MR={MR:.2f}, eps={eps:.1f}: {e}")
                        if self.use_3d:
                            self.cstar_table[i, j, k_idx] = np.nan
                            self.Cf_table[i, j, k_idx] = np.nan
                            self.Tc_table[i, j, k_idx] = np.nan
                            self.gamma_table[i, j, k_idx] = np.nan
                            self.R_table[i, j, k_idx] = np.nan
                            self.M_table[i, j, k_idx] = np.nan
                        else:
                            self.cstar_table[i, j] = np.nan
                            self.Cf_table[i, j] = np.nan
                            self.Tc_table[i, j] = np.nan
                            self.gamma_table[i, j] = np.nan
                            self.R_table[i, j] = np.nan
                            self.M_table[i, j] = np.nan
    
    def _build_cache_parallel(self, n_workers: int, total_points: int):
        """Parallel CEA cache building with chunking"""
        # Create all grid points as list of tuples (i, j, k, Pc_psia, MR, eps)
        grid_points = []
        Pc_psia_grid = self.Pc_grid / 6894.76
        eps_list = self.eps_grid if self.use_3d else [self.config.expansion_ratio]
        
        for i, Pc_psia in enumerate(Pc_psia_grid):
            for j, MR in enumerate(self.MR_grid):
                for k_idx, eps in enumerate(eps_list):
                    grid_points.append((i, j, k_idx, Pc_psia, MR, eps))
        
        # Initialize tables
        shape = (self.n_points, self.n_points, self.n_points) if self.use_3d else (self.n_points, self.n_points)
        self.cstar_table = np.zeros(shape)
        self.Cf_table = np.zeros(shape)
        self.Tc_table = np.zeros(shape)
        self.gamma_table = np.zeros(shape)
        self.R_table = np.zeros(shape)
        self.M_table = np.zeros(shape)
        self.Cf_vac_table = np.zeros(shape)
        
        # Chunk grid points for parallel processing
        chunk_size = max(1, len(grid_points) // (n_workers * 4))  # 4 chunks per worker
        chunks = [grid_points[i:i+chunk_size] for i in range(0, len(grid_points), chunk_size)]
        print(f"   [PARALLEL] Split into {len(chunks)} chunks of ~{chunk_size} points each")
        
        # Create a shared lock to serialize RocketCEA calls (RocketCEA uses shared temp files)
        # This prevents Fortran I/O errors when multiple processes access RocketCEA simultaneously
        manager = Manager()
        cea_lock = manager.Lock()
        
        # Prepare worker function
        # The function is defined at module level (above), so it's directly accessible
        # On Windows with 'spawn', workers will re-import this module and find the function
        worker_func = partial(
            _compute_cea_point_chunk,
            ox_name=self.config.ox_name,
            fuel_name=self.config.fuel_name,
            expansion_ratio=self.config.expansion_ratio if not self.use_3d else None,
            lock=cea_lock,  # Pass lock to serialize RocketCEA calls
        )
        
        # Process chunks in parallel
        start_time = time.time()
        completed = 0
        
        with Pool(processes=n_workers) as pool:
            results = pool.imap_unordered(worker_func, chunks)
            
            for chunk_results in results:
                completed += len(chunk_results)
                elapsed = time.time() - start_time
                pct = 100 * completed / total_points
                rate = completed / elapsed if elapsed > 0 else 0
                eta = (total_points - completed) / rate if rate > 0 else 0
                print(f"   [{completed}/{total_points}] {pct:.1f}% | Rate: {rate:.1f} pts/s | ETA: {eta/60:.1f} min")
                
                # Store results
                for i, j, k_idx, cstar, Cf, Tc, gamma, R, M, Cf_vac in chunk_results:
                    if self.use_3d:
                        self.cstar_table[i, j, k_idx] = cstar
                        self.Cf_table[i, j, k_idx] = Cf
                        self.Cf_vac_table[i, j, k_idx] = Cf_vac
                        self.Tc_table[i, j, k_idx] = Tc
                        self.gamma_table[i, j, k_idx] = gamma
                        self.R_table[i, j, k_idx] = R
                        self.M_table[i, j, k_idx] = M
                    else:
                        self.cstar_table[i, j] = cstar
                        self.Cf_table[i, j] = Cf
                        self.Cf_vac_table[i, j] = Cf_vac
                        self.Tc_table[i, j] = Tc
                        self.gamma_table[i, j] = gamma
                        self.R_table[i, j] = R
                        self.M_table[i, j] = M
 
    def _save_cache(self):
        """Save CEA data to cache file"""
        meta = {
            "table_schema": CEA_TABLE_SCHEMA_VERSION,
            "ox_name": self.config.ox_name,
            "fuel_name": self.config.fuel_name,
            "expansion_ratio": self.config.expansion_ratio,
            "Pc_range": list(self.config.Pc_range),
            "MR_range": list(self.config.MR_range),
            "n_points": self.n_points,
            "dimensions": 3 if self.use_3d else 2,
            "eps_range": list(self.config.eps_range) if self.use_3d else None,
        }
 
        save_dict = {
            "Pc": self.Pc_grid,
            "MR": self.MR_grid,
            "cstar": self.cstar_table,
            "Cf": self.Cf_table,
            "Cf_vac": self.Cf_vac_table,
            "Tc": self.Tc_table,
            "gamma": self.gamma_table,
            "R": self.R_table,
            "M": self.M_table,
            "meta": np.array(json.dumps(meta))
        }
        
        # Add eps grid for 3D caches
        if self.use_3d:
            save_dict["eps"] = self.eps_grid
        
        np.savez_compressed(self.cache_file, **save_dict)
    
    def _bilinear_interpolate(self, Pc: float, MR: float, table: np.ndarray) -> float:
        """
        Bilinear interpolation in (Pc, MR) space with robust error handling.
        
        Uses improved boundary handling and NaN management for numerical stability.
        """
        # Validate inputs
        if not (np.isfinite(Pc) and np.isfinite(MR)):
            raise ValueError(f"Non-finite interpolation point: Pc={Pc}, MR={MR}")
        
        # Find indices using searchsorted (returns insertion point)
        i_pc = np.searchsorted(self.Pc_grid, Pc)
        i_mr = np.searchsorted(self.MR_grid, MR)
        
        # Handle boundary cases robustly
        # If at or beyond upper bound, use last two points
        if i_pc >= len(self.Pc_grid):
            i_pc = len(self.Pc_grid) - 1
        if i_mr >= len(self.MR_grid):
            i_mr = len(self.MR_grid) - 1
        
        # If at or below lower bound, use first two points
        if i_pc == 0:
            i_pc = 1
        if i_mr == 0:
            i_mr = 1
        
        # Ensure we have valid indices for interpolation
        i_pc = np.clip(i_pc, 1, len(self.Pc_grid) - 1)
        i_mr = np.clip(i_mr, 1, len(self.MR_grid) - 1)
        
        # Get surrounding points
        Pc0, Pc1 = self.Pc_grid[i_pc - 1], self.Pc_grid[i_pc]
        MR0, MR1 = self.MR_grid[i_mr - 1], self.MR_grid[i_mr]
        
        # Validate grid spacing
        if Pc1 <= Pc0 or MR1 <= MR0:
            raise ValueError(f"Invalid grid spacing: Pc grid or MR grid not monotonic")
        
        # Get corner values
        f00 = table[i_pc - 1, i_mr - 1]
        f01 = table[i_pc - 1, i_mr]
        f10 = table[i_pc, i_mr - 1]
        f11 = table[i_pc, i_mr]
        
        # Count valid (non-NaN) values
        valid_values = [v for v in [f00, f01, f10, f11] if np.isfinite(v)]
        
        if len(valid_values) == 0:
            # All NaN - return NaN
            return np.nan
        elif len(valid_values) < 4:
            # Some NaN - use weighted average of valid values only
            # This is more robust than simple nearest neighbor
            weights = []
            values = []
            
            # Calculate weights for each corner
            for i, (f_val, pc_idx, mr_idx) in enumerate([
                (f00, i_pc - 1, i_mr - 1),
                (f01, i_pc - 1, i_mr),
                (f10, i_pc, i_mr - 1),
                (f11, i_pc, i_mr),
            ]):
                if np.isfinite(f_val):
                    # Distance-based weight (inverse distance)
                    pc_dist = abs(Pc - self.Pc_grid[pc_idx]) / (Pc1 - Pc0) if Pc1 != Pc0 else 1.0
                    mr_dist = abs(MR - self.MR_grid[mr_idx]) / (MR1 - MR0) if MR1 != MR0 else 1.0
                    weight = 1.0 / (1.0 + pc_dist + mr_dist)
                    weights.append(weight)
                    values.append(f_val)
            
            if len(weights) > 0:
                weights = np.array(weights)
                weights = weights / np.sum(weights)  # Normalize
                result = np.sum(np.array(values) * weights)
            else:
                result = np.nan
        else:
            # All values valid - standard bilinear interpolation
            # Interpolation weights (normalized to [0, 1])
            wx = (Pc - Pc0) / (Pc1 - Pc0) if Pc1 != Pc0 else 0.0
            wy = (MR - MR0) / (MR1 - MR0) if MR1 != MR0 else 0.0
            
            # Clamp weights to [0, 1] for numerical stability
            wx = np.clip(wx, 0.0, 1.0)
            wy = np.clip(wy, 0.0, 1.0)
            
            # Bilinear interpolation
            result = (f00 * (1 - wx) * (1 - wy) +
                     f10 * wx * (1 - wy) +
                     f01 * (1 - wx) * wy +
                     f11 * wx * wy)
        
        # Validate result
        if not np.isfinite(result):
            # Final fallback: use nearest valid neighbor
            distances = [
                (abs(Pc - Pc0) + abs(MR - MR0), f00),
                (abs(Pc - Pc0) + abs(MR - MR1), f01),
                (abs(Pc - Pc1) + abs(MR - MR0), f10),
                (abs(Pc - Pc1) + abs(MR - MR1), f11),
            ]
            valid_distances = [(d, v) for d, v in distances if np.isfinite(v)]
            if valid_distances:
                result = min(valid_distances, key=lambda x: x[0])[1]
            else:
                result = np.nan
        
        return float(result)
    
    def _trilinear_interpolate(self, Pc: float, MR: float, eps: float, table: np.ndarray) -> float:
        """Trilinear interpolation in (Pc, MR, eps) space"""
        # Find indices
        i_pc = np.searchsorted(self.Pc_grid, Pc)
        i_mr = np.searchsorted(self.MR_grid, MR)
        i_eps = np.searchsorted(self.eps_grid, eps)
        
        # Clamp to valid range. Scalar min/max on a Python int is far cheaper than
        # np.clip on a numpy scalar (which allocates + builds getlimits) and is
        # identical for an in-range index. This runs 6x per CEA lookup.
        i_pc = min(max(int(i_pc), 1), len(self.Pc_grid) - 1)
        i_mr = min(max(int(i_mr), 1), len(self.MR_grid) - 1)
        i_eps = min(max(int(i_eps), 1), len(self.eps_grid) - 1)
        
        # Get surrounding points
        Pc0, Pc1 = self.Pc_grid[i_pc - 1], self.Pc_grid[i_pc]
        MR0, MR1 = self.MR_grid[i_mr - 1], self.MR_grid[i_mr]
        eps0, eps1 = self.eps_grid[i_eps - 1], self.eps_grid[i_eps]
        
        # Get 8 corner values
        f000 = table[i_pc - 1, i_mr - 1, i_eps - 1]
        f001 = table[i_pc - 1, i_mr - 1, i_eps]
        f010 = table[i_pc - 1, i_mr, i_eps - 1]
        f011 = table[i_pc - 1, i_mr, i_eps]
        f100 = table[i_pc, i_mr - 1, i_eps - 1]
        f101 = table[i_pc, i_mr - 1, i_eps]
        f110 = table[i_pc, i_mr, i_eps - 1]
        f111 = table[i_pc, i_mr, i_eps]
        
        # Check for NaN values
        corners = [f000, f001, f010, f011, f100, f101, f110, f111]
        if any(np.isnan(c) for c in corners):
            # Fallback to nearest neighbor
            return table[i_pc - 1, i_mr - 1, i_eps - 1]
        
        # Interpolation weights
        wx = (Pc - Pc0) / (Pc1 - Pc0) if Pc1 != Pc0 else 0
        wy = (MR - MR0) / (MR1 - MR0) if MR1 != MR0 else 0
        wz = (eps - eps0) / (eps1 - eps0) if eps1 != eps0 else 0
        
        # Trilinear interpolation
        result = (
            f000 * (1 - wx) * (1 - wy) * (1 - wz) +
            f100 * wx * (1 - wy) * (1 - wz) +
            f010 * (1 - wx) * wy * (1 - wz) +
            f110 * wx * wy * (1 - wz) +
            f001 * (1 - wx) * (1 - wy) * wz +
            f101 * wx * (1 - wy) * wz +
            f011 * (1 - wx) * wy * wz +
            f111 * wx * wy * wz
        )
        
        return float(result)
    
    def eval(self, MR: float, Pc: float, Pa: float = 101325.0, eps: Optional[float] = None) -> dict:
        """
        Evaluate CEA properties at given conditions.
        
        Parameters:
        -----------
        MR : float
            Mixture ratio (O/F)
        Pc : float
            Chamber pressure [Pa]
        Pa : float
            Ambient pressure [Pa] (default: sea level)
        eps : float, optional
            Expansion ratio (uses config default if None)
        
        Returns:
        --------
        dict with keys:
            cstar_ideal : float [m/s]
            Cf_ideal : float, ambient thrust coefficient at THIS Pa: Cf_vac - Pa*eps/Pc
            Cf_sl : float, the stored table (CEA get_PambCf at its default 14.7 psia)
            Cf_vac : float
            Tc : float [K]
            gamma : float
            R : float [J/(kg*K)]
            M : float [kg/kmol] (molecular weight)
            extrapolated : bool, True when Pc, MR or eps lay outside the table and were clamped
            clamped : tuple of the clamped axis names ("Pc", "MR", "eps")
        """
        if eps is None:
            eps = self.config.expansion_ratio
 
        Pc_in, MR_in, eps_in = float(Pc), float(MR), float(eps)
 
        # Clamp (scalar; matches float(np.clip(...)) incl. NaN pass-through)
        Pc_clamped = _scalar_clamp(Pc_in, self.Pc_min, self.Pc_max)
        MR_clamped = _scalar_clamp(MR_in, self.MR_min, self.MR_max)
        eps_clamped = eps_in
        if self.use_3d:
            eps_clamped = _scalar_clamp(eps_in, self.eps_min, self.eps_max)
 
        # print(
        #     "[CEA DEBUG] in:",
        #     f"Pc={Pc_in:.6g} Pa, MR={MR_in:.6g}, Pa={Pa:.6g} Pa, eps={eps_in:.6g}, use_3d={self.use_3d}"
        # )
        # print(
        #     "[CEA DEBUG] bounds:",
        #     f"Pc=[{self.Pc_min:.6g}, {self.Pc_max:.6g}]",
        #     f"MR=[{self.MR_min:.6g}, {self.MR_max:.6g}]",
        #     f"eps=[{self.eps_min:.6g}, {self.eps_max:.6g}]"
        # )
        # if (Pc_in != Pc_clamped) or (MR_in != MR_clamped) or (eps_in != eps_clamped):
        #     print(
        #         "[CEA DEBUG][CLAMPED]:",
        #         f"Pc {Pc_in:.6g}->{Pc_clamped:.6g},",
        #         f"MR {MR_in:.6g}->{MR_clamped:.6g},",
        #         f"eps {eps_in:.6g}->{eps_clamped:.6g}"
        #     )
 
        # Interpolate
        if self.use_3d:
            cstar = float(self._trilinear_interpolate(Pc_clamped, MR_clamped, eps_clamped, self.cstar_table))
            Cf    = float(self._trilinear_interpolate(Pc_clamped, MR_clamped, eps_clamped, self.Cf_table))
            Tc    = float(self._trilinear_interpolate(Pc_clamped, MR_clamped, eps_clamped, self.Tc_table))
            gamma = float(self._trilinear_interpolate(Pc_clamped, MR_clamped, eps_clamped, self.gamma_table))
            R     = float(self._trilinear_interpolate(Pc_clamped, MR_clamped, eps_clamped, self.R_table))
            M     = float(self._trilinear_interpolate(Pc_clamped, MR_clamped, eps_clamped, self.M_table))
            Cf_vac = (float(self._trilinear_interpolate(Pc_clamped, MR_clamped, eps_clamped, self.Cf_vac_table))
                      if self.Cf_vac_table is not None else _isentropic_cf_vac(gamma, eps_clamped))
        else:
            cstar = float(self._bilinear_interpolate(Pc_clamped, MR_clamped, self.cstar_table))
            Cf    = float(self._bilinear_interpolate(Pc_clamped, MR_clamped, self.Cf_table))
            Tc    = float(self._bilinear_interpolate(Pc_clamped, MR_clamped, self.Tc_table))
            gamma = float(self._bilinear_interpolate(Pc_clamped, MR_clamped, self.gamma_table))
            R     = float(self._bilinear_interpolate(Pc_clamped, MR_clamped, self.R_table))
            M     = float(self._bilinear_interpolate(Pc_clamped, MR_clamped, self.M_table))
            Cf_vac = (float(self._bilinear_interpolate(Pc_clamped, MR_clamped, self.Cf_vac_table))
                      if self.Cf_vac_table is not None else _isentropic_cf_vac(gamma, eps_clamped))
 
        # The ambient coefficient belongs to the Pa asked about, not to the table's 14.7 psia.
        Cf_amb = (Cf_vac - float(Pa) * eps_in / Pc_in) if (Pc_in > 0 and np.isfinite(Cf_vac)) else Cf
        clamped = tuple(name for name, a, b in (("Pc", Pc_in, Pc_clamped), ("MR", MR_in, MR_clamped),
                                                ("eps", eps_in, eps_clamped)) if a != b)
        out = {
            "cstar_ideal": cstar,
            "Cf_ideal": float(Cf_amb),
            "Cf_sl": Cf,
            "Cf_vac": Cf_vac,   # vacuum thrust coefficient (RPA delivered-thrust basis)
            "Tc": Tc,
            "gamma": gamma,
            "R": R,
            "M": M,
            "extrapolated": bool(clamped),
            "clamped": clamped,
        }
 
        # print(
        #     "[CEA DEBUG] out:",
        #     f"c*={cstar:.3f} m/s, Cf={Cf:.4f}, Tc={Tc:.1f} K, gamma={gamma:.4f}, R={R:.3f}, M={M:.3f}"
        # )
 
        # Loud sanity flags (these are the ones that usually catch unit bugs instantly)
        if not (1200.0 <= cstar <= 2200.0):
            print("[CEA DEBUG][WARNING] c* is out of expected LOX/RP-1-ish range. Unit/table issue likely.")
        if not (1.05 <= gamma <= 1.40):
            print("[CEA DEBUG][WARNING] gamma is weird.")
        if not (2000.0 <= Tc <= 4200.0):
            print("[CEA DEBUG][WARNING] Tc is weird.")
        # Typical combustion-gas R is a few hundred J/(kg·K). If you see ~3000+, it's probably kJ/kmol-K or similar unit mismatch.
        if R > 2000.0 or R < 50.0:
            print("[CEA DEBUG][WARNING] R magnitude looks wrong. Possible units mismatch (J/kg-K vs something else).")
 
        return out
 
    def eval_cf_vac(self, MR: float, Pc: float, eps: Optional[float] = None) -> float:
        """Fast single-value vacuum thrust coefficient lookup (one interpolation).
 
        Hot path for the native Layer-1 inner loop: avoids the full ``eval`` dict
        build + sanity prints. Same value ``eval(...)["Cf_vac"]`` returns.
        """
        if eps is None:
            eps = self.config.expansion_ratio
        Pc_c = _scalar_clamp(float(Pc), self.Pc_min, self.Pc_max)
        MR_c = _scalar_clamp(float(MR), self.MR_min, self.MR_max)
        if self.use_3d:
            eps_c = _scalar_clamp(float(eps), self.eps_min, self.eps_max)
            if self.Cf_vac_table is not None:
                return float(self._trilinear_interpolate(Pc_c, MR_c, eps_c, self.Cf_vac_table))
            g = float(self._trilinear_interpolate(Pc_c, MR_c, eps_c, self.gamma_table))
            return _isentropic_cf_vac(g, eps_c)
        if self.Cf_vac_table is not None:
            return float(self._bilinear_interpolate(Pc_c, MR_c, self.Cf_vac_table))
        g = float(self._bilinear_interpolate(Pc_c, MR_c, self.gamma_table))
        return _isentropic_cf_vac(g, float(eps))
 
    @property
    def aux(self) -> "CEAAuxTables":
        """Transport, composition and exit-state tables for this propellant pair."""
        return get_aux_tables(self.config.ox_name, self.config.fuel_name)
 
 
# ── Auxiliary CEA tables ──────────────────────────────────────────────────────
# What the performance cache above does not carry, and the thermal and nozzle models need:
#   * frozen transport (mu, cp, Pr) at the chamber and throat stations -- the Bartz inputs;
#   * equilibrium mole fractions of the radiating and carbon-oxidising species, chamber and
#     throat, and the throat molecular weight -- gas emissivity and graphite B';
#   * the shifting-equilibrium exit state Pc/Pe, Te, Me over an area-ratio grid.
# About 50k rocketcea calls, ~40 s: built on demand, written next to the performance cache,
# and meant to be committed with it (CI has no rocketcea). The grid is fixed per propellant
# pair and wide enough for every preset, so one file serves every config of that pair.
CEA_AUX_SCHEMA_VERSION = 1
AUX_SPECIES = ("H2O", "CO2", "CO", "H2", "OH", "O2", "O", "H")
_AUX_MR = np.round(np.linspace(0.8, 4.0, 65), 6)
_AUX_PC = np.geomspace(1.0e5, 1.2e7, 16)        # Pa
_AUX_EPS = np.geomspace(1.1, 60.0, 48)
_PSI = 6894.757293168361
_BTU_LBM_R = 4186.8                              # J/(kg K)
_MILLIPOISE = 1.0e-4                             # Pa s
_MCAL_CM_S_K = 0.4184                            # W/(m K)
_AUX_MEMO: Dict[Tuple[str, str], "CEAAuxTables"] = {}
 
 
def _aux_axis(grid: np.ndarray, x: float, log: bool) -> Tuple[int, float, bool]:
    """Cell index, weight and clamp flag on a uniform (or log-uniform) grid."""
    lo, hi = float(grid[0]), float(grid[-1])
    xc = min(max(x, lo), hi)
    if log:
        f = (np.log(xc) - np.log(lo)) / (np.log(hi) - np.log(lo)) * (len(grid) - 1)
    else:
        f = (xc - lo) / (hi - lo) * (len(grid) - 1)
    i = min(int(f), len(grid) - 2)
    return i, f - i, xc != x
 
 
def _fill_nan_nearest(a: np.ndarray) -> np.ndarray:
    """Fill failed CEA points from the nearest good one along the first axis (MR), then Pc."""
    a = np.array(a, dtype=float)
    for axis in range(min(a.ndim, 2)):
        moved = np.moveaxis(a, axis, 0)
        for idx in np.ndindex(moved.shape[1:]):
            col = moved[(slice(None),) + idx]
            bad = ~np.isfinite(col)
            if bad.any() and (~bad).any():
                good = np.flatnonzero(~bad)
                for j in np.flatnonzero(bad):
                    col[j] = col[good[np.argmin(np.abs(good - j))]]
    return a
 
 
class CEAAuxTables:
    """Chamber/throat transport and composition on (MR, Pc); exit state on (MR, Pc, eps)."""
 
    def __init__(self, ox_name: str, fuel_name: str, path: Optional[str] = None):
        self.ox_name, self.fuel_name = ox_name, fuel_name
        if path is None:
            root = os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
            path = os.path.join(root, "output", "cache", f"cea_aux_{ox_name}_{fuel_name}.npz")
        self.path = path
        self.MR, self.Pc, self.eps = _AUX_MR, _AUX_PC, _AUX_EPS
        self.t: Dict[str, np.ndarray] = {}
        if not self._load():
            self._build()
 
    def _load(self) -> bool:
        if not os.path.exists(self.path):
            return False
        data = _load_npz_tables(self.path)
        try:
            meta = json.loads(data["meta"].tolist())
        except Exception:
            return False
        if (meta.get("schema") != CEA_AUX_SCHEMA_VERSION or meta.get("ox_name") != self.ox_name
                or meta.get("fuel_name") != self.fuel_name):
            return False
        if not (np.array_equal(data["MR"], self.MR) and np.allclose(data["Pc"], self.Pc)
                and np.allclose(data["eps"], self.eps)):
            return False
        self.t = {k: v for k, v in data.items() if k not in ("meta", "MR", "Pc", "eps")}
        return True
 
    def _build(self) -> None:
        C = _get_CEA_Obj()(oxName=self.ox_name, fuelName=self.fuel_name)
        nM, nP, nE = len(self.MR), len(self.Pc), len(self.eps)
        t = {k: np.full((nM, nP), np.nan) for k in
             ("mu_c", "cp_c", "Pr_c", "mu_t", "cp_t", "Pr_t", "MW_t")}
        for stn in ("c", "t"):
            for sp in AUX_SPECIES:
                t[f"x_{sp}_{stn}"] = np.zeros((nM, nP))
        for k in ("lnPcPe", "Te", "Me"):
            t[k] = np.full((nM, nP, nE), np.nan)
        for i, MR in enumerate(self.MR):
            for j, Pc in enumerate(self.Pc):
                pc = float(Pc) / _PSI
                try:
                    for stn, fn in (("c", C.get_Chamber_Transport), ("t", C.get_Throat_Transport)):
                        cp, mu, _k, pr = fn(Pc=pc, MR=float(MR), eps=2.0, frozen=1)
                        t[f"cp_{stn}"][i, j] = cp * _BTU_LBM_R
                        t[f"mu_{stn}"][i, j] = mu * _MILLIPOISE
                        t[f"Pr_{stn}"][i, j] = pr
                    t["MW_t"][i, j] = C.get_Throat_MolWt_gamma(Pc=pc, MR=float(MR), eps=2.0)[0]
                    _, xs = C.get_SpeciesMoleFractions(Pc=pc, MR=float(MR), eps=2.0, frozen=0,
                                                       frozenAtThroat=0, min_fraction=1e-8)
                    for key, vals in xs.items():
                        sp = key.lstrip("*")
                        if sp in AUX_SPECIES:
                            t[f"x_{sp}_c"][i, j] = vals[1]
                            t[f"x_{sp}_t"][i, j] = vals[2]
                except Exception:
                    pass
                for k, e in enumerate(self.eps):
                    try:
                        pcpe = C.get_PcOvPe(Pc=pc, MR=float(MR), eps=float(e), frozen=0)
                        te = C.get_Temperatures(Pc=pc, MR=float(MR), eps=float(e), frozen=0)[2] / 1.8
                        me = C.get_MachNumber(Pc=pc, MR=float(MR), eps=float(e), frozen=0)
                        if pcpe > 1.0 and te > 0 and me > 1.0:
                            t["lnPcPe"][i, j, k] = np.log(pcpe)
                            t["Te"][i, j, k] = te
                            t["Me"][i, j, k] = me
                    except Exception:
                        pass
        self.t = {k: _fill_nan_nearest(v) for k, v in t.items()}
        meta = json.dumps({"schema": CEA_AUX_SCHEMA_VERSION, "ox_name": self.ox_name,
                           "fuel_name": self.fuel_name,
                           "note": "frozen transport, equilibrium composition (index 1 chamber, "
                                   "2 throat), shifting-equilibrium exit state"})
        os.makedirs(os.path.dirname(self.path), exist_ok=True)
        tmp = self.path + f".tmp{os.getpid()}.npz"
        np.savez_compressed(tmp, meta=np.array(meta), MR=self.MR, Pc=self.Pc, eps=self.eps, **self.t)
        os.replace(tmp, self.path)
 
    def _bilinear(self, name: str, MR: float, Pc: float) -> Tuple[float, bool]:
        i, wi, ci = _aux_axis(self.MR, float(MR), False)
        j, wj, cj = _aux_axis(self.Pc, float(Pc), True)
        a = self.t[name]
        v = ((1 - wi) * (1 - wj) * a[i, j] + wi * (1 - wj) * a[i + 1, j]
             + (1 - wi) * wj * a[i, j + 1] + wi * wj * a[i + 1, j + 1])
        return float(v), (ci or cj)
 
    def transport(self, MR: float, Pc: float, station: str = "chamber") -> Dict[str, float]:
        """Frozen-composition mu [Pa s], cp [J/(kg K)], Pr, k = mu cp / Pr [W/(m K)]."""
        s = "c" if station == "chamber" else "t"
        mu, e1 = self._bilinear(f"mu_{s}", MR, Pc)
        cp, e2 = self._bilinear(f"cp_{s}", MR, Pc)
        pr, e3 = self._bilinear(f"Pr_{s}", MR, Pc)
        return {"mu": mu, "cp": cp, "Pr": pr, "k": mu * cp / pr, "extrapolated": e1 or e2 or e3}
 
    def composition(self, MR: float, Pc: float, station: str = "chamber") -> Dict[str, float]:
        """Equilibrium mole fractions of AUX_SPECIES; 'MW' is the throat molecular weight."""
        s = "c" if station == "chamber" else "t"
        out: Dict[str, float] = {}
        ext = False
        for sp in AUX_SPECIES:
            out[sp], e = self._bilinear(f"x_{sp}_{s}", MR, Pc)
            out[sp] = max(out[sp], 0.0)
            ext = ext or e
        out["MW"], _ = self._bilinear("MW_t", MR, Pc)
        out["extrapolated"] = ext
        return out
 
    def _exit(self, name: str, MR: float, Pc: float, eps: float) -> Tuple[float, bool]:
        i, wi, ci = _aux_axis(self.MR, float(MR), False)
        j, wj, cj = _aux_axis(self.Pc, float(Pc), True)
        k, wk, ck = _aux_axis(self.eps, float(eps), True)
        a = self.t[name]
        v = 0.0
        for di, fi in ((0, 1 - wi), (1, wi)):
            for dj, fj in ((0, 1 - wj), (1, wj)):
                for dk, fk in ((0, 1 - wk), (1, wk)):
                    v += fi * fj * fk * a[i + di, j + dj, k + dk]
        return float(v), (ci or cj or ck)
 
    def exit_state(self, MR: float, Pc: float, eps: float) -> Dict[str, float]:
        """Shifting-equilibrium exit pressure [Pa], temperature [K] and Mach number."""
        ln_r, e1 = self._exit("lnPcPe", MR, Pc, eps)
        Te, _ = self._exit("Te", MR, Pc, eps)
        Me, _ = self._exit("Me", MR, Pc, eps)
        return {"P_exit": float(Pc) / np.exp(ln_r), "T_exit": Te, "M_exit": Me,
                "PcOvPe": float(np.exp(ln_r)), "extrapolated": e1}
 
    def eps_for_exit_pressure(self, MR: float, Pc: float, Pe: float) -> float:
        """Area ratio whose shifting-equilibrium exit pressure is Pe (Sutton: optimum at Pe = Pa)."""
        target = np.log(float(Pc) / float(Pe))
        f = lambda e: self._exit("lnPcPe", MR, Pc, e)[0] - target
        lo, hi = float(self.eps[0]), float(self.eps[-1])
        flo, fhi = f(lo), f(hi)
        if flo >= 0:
            return lo
        if fhi <= 0:
            return hi
        from scipy.optimize import brentq
        return float(brentq(f, lo, hi, xtol=1e-10, rtol=1e-12))
 
 
def get_aux_tables(ox_name: str, fuel_name: str) -> CEAAuxTables:
    key = (ox_name, fuel_name)
    tab = _AUX_MEMO.get(key)
    if tab is None:
        tab = CEAAuxTables(ox_name, fuel_name)
        _AUX_MEMO[key] = tab
    return tab