Skip to content

Zalmoxis.eos analytic

eos_analytic

Analytic EOS from Seager et al. (2007), Table 3, Eq. 11.

Modified polytropic fit:

\[ \rho(P) = \rho_0 + c \cdot P^n \]

Valid for \(P < 10^{16}\) Pa. Approximates the full merged Vinet/BME + TFD EOS to 2-12% accuracy for all planetary pressures.

P_MAX = 1e+16 module-attribute

SEAGER2007_MATERIALS = {'iron': (8300.0, 0.00349, 0.528), 'MgSiO3': (4100.0, 0.00161, 0.541), 'MgFeSiO3': (4260.0, 0.00127, 0.549), 'H2O': (1460.0, 0.00311, 0.513), 'graphite': (2250.0, 0.0035, 0.514), 'SiC': (3220.0, 0.00172, 0.537)} module-attribute

VALID_MATERIAL_KEYS = set(SEAGER2007_MATERIALS.keys()) module-attribute

logger = logging.getLogger(__name__) module-attribute

get_analytic_density(pressure, material_key)

Compute density from the Seager et al. (2007) analytic modified polytrope.

Parameters:

Name Type Description Default
pressure float

Pressure in Pa.

 required
material_key str

One of the keys in SEAGER2007_MATERIALS (e.g., "iron", "MgSiO3", "H2O", "graphite", "SiC", "MgFeSiO3").

 required

Returns:

Type Description
float

Density in kg/m^3. Returns rho_0 for non-positive pressure (allows the ODE solver to remain stable during pressure adjustment iterations). Returns None only for NaN pressure.

Raises:

Type Description
ValueError

If material_key is not recognized.
Source code in src/zalmoxis/eos_analytic.py 
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
def get_analytic_density(pressure: float, material_key: str) -> float | None:
    """
    Compute density from the Seager et al. (2007) analytic modified polytrope.

    Parameters
    ----------
    pressure : float
        Pressure in Pa.
    material_key : str
        One of the keys in SEAGER2007_MATERIALS
        (e.g., "iron", "MgSiO3", "H2O", "graphite", "SiC", "MgFeSiO3").

    Returns
    -------
    float
        Density in kg/m^3. Returns rho_0 for non-positive pressure (allows
        the ODE solver to remain stable during pressure adjustment iterations).
        Returns None only for NaN pressure.

    Raises
    ------
    ValueError
        If material_key is not recognized.
    """
    if material_key not in SEAGER2007_MATERIALS:
        raise ValueError(
            f"Unknown material key '{material_key}'. Valid keys: {sorted(VALID_MATERIAL_KEYS)}"
        )

    rho_0, c, n = SEAGER2007_MATERIALS[material_key]

    # Guard against nonphysical pressure
    if np.isnan(pressure):
        return None
    if pressure <= 0:
        return float(rho_0)

    if pressure > P_MAX:
        logger.warning(
            f'Pressure {pressure:.2e} Pa exceeds validity limit of {P_MAX:.0e} Pa '
            f'for Seager+2007 analytic EOS. Results may be inaccurate.'
        )

    density = rho_0 + c * pressure**n

    return float(density)