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 \]

The six Seager materials are empirical fits valid for \(P < 10^{16}\) Pa, approximating the full merged Vinet/BME + TFD EOS to 2-12% accuracy across planetary pressures. The module also registers an exact n=1 polytrope (\(\rho_0 = 0\), \(n = 1/2\), so \(P = K \rho^2\)) as a verification material with a closed-form Lane-Emden solution, used to check the coupled structure solver end to end.

ANALYTIC_MATERIALS = {None: SEAGER2007_MATERIALS, None: VERIFICATION_MATERIALS} module-attribute

G = 6.67428e-11 module-attribute

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

USER_SELECTABLE_MATERIALS = set(SEAGER2007_MATERIALS) module-attribute

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

VERIFICATION_MATERIALS = {'polytrope_n1': (0.0, _K_POLYTROPE_N1 ** -0.5, 0.5)} module-attribute

_K_POLYTROPE_N1 = 2.0 * np.pi * G * (earth_radius / np.pi) ** 2 module-attribute

earth_radius = 6335439.0 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 ANALYTIC_MATERIALS (e.g., "iron", "MgSiO3", "H2O", "graphite", "SiC", "MgFeSiO3", or the verification material "polytrope_n1").

 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 
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
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 ANALYTIC_MATERIALS
        (e.g., "iron", "MgSiO3", "H2O", "graphite", "SiC", "MgFeSiO3",
        or the verification material "polytrope_n1").

    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 ANALYTIC_MATERIALS:
        raise ValueError(
            f"Unknown material key '{material_key}'. Valid keys: {sorted(VALID_MATERIAL_KEYS)}"
        )

    rho_0, c, n = ANALYTIC_MATERIALS[material_key]

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

    # The validity limit is an empirical bound on the Seager+2007 fits. The polytrope
    # verification material is an exact P = K rho^2 relation at every pressure, so it is
    # excluded from the warning to avoid mislabelling an exact EOS as inaccurate.
    if material_key in SEAGER2007_MATERIALS and 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)