aragog.core

def conducted_adiabatic_flow(self, r, t_cmb):
    """Heat flow [W] conducted along the adiabat through radius ``r`` (needs ``k_core``)."""
    if self.k_core is None:
        raise ValueError('conducted_adiabatic_flow needs k_core (a stratified budget)')
    return _q_ad(self.profiles, self.k_core, r, t_cmb)

def core_rates(self, t_cmb, t_shell, q_cmb, q_sources=0.0):
    """Rates of the convecting core [K/s] and of the shell cells [K/s] of a stratified core.

    The convecting core loses the flow into the shell base and keeps the mass share of the
    internal sources ``q_sources`` [W] below it; the shell gets the rest, cell by cell.
    """
    p, shell = self.profiles, self.shell
    heating = q_sources / p.enclosed_mass(p.r_cmb)
    d_shell, q_base = shell.rates(t_shell, t_cmb, q_cmb, heating)
    q_conv = heating * p.enclosed_mass(self.r_convecting)
    capacity = self.effective_capacity(
        t_cmb, gravitational_upper=shell.layer_base(t_shell, t_cmb)
    )
    return (q_conv - q_base) / capacity, d_shell

def dtcmb_dt(self, t_cmb, q_cmb, q_sources=0.0):
    """CMB cooling rate [K/s] for heat flow ``q_cmb`` [W] out of an unstratified core.

    ``dT_cmb/dt = (q_sources - q_cmb) / C_eff(T_cmb)``; positive
    ``q_cmb`` cools the core, and internal sources (radiogenic, tidal)
    offset it. A stratified core has shell state: see :meth:`core_rates`.
    """
    if self.stratification:
        raise ValueError('a stratified core evolves its shell too: use core_rates')
    return (q_sources - q_cmb) / self.effective_capacity(t_cmb)

def effective_capacity(self, t_cmb, *, gravitational_upper=None):
    """Total dQ/d(dT_cmb/dt) [J/K] of the convecting core: secular plus
    latent plus gravitational (profile mode).

    The secular integral runs to the convecting radius, the CMB or the shell base; the
    light elements mix up to ``gravitational_upper`` (default the same radius), the base of
    the stable layer when a shell is resolved.
    """
    if self.capacity_mode == 'legacy':
        return self.secular_capacity()
    upper = self.r_convecting
    top = upper if gravitational_upper is None else gravitational_upper
    return (
        self.secular_capacity(upper=upper)
        + self.latent_capacity(t_cmb)
        + self.gravitational_capacity(t_cmb, upper=top)
    )

def effective_tfac(self):
    """Mass-weighted mean of ``T_a/T_cmb`` over the core (profile mode)."""
    p = self.profiles

    def mass_integrand(r):
        return p.density(r) * 4.0 * jnp.pi * r**2

    mass = self._quad_0_upper(p.r_cmb, mass_integrand)
    return self.secular_capacity() / (p.c_p * mass)

def gravitational_capacity(self, t_cmb, upper=None):
    """Gravitational contribution to the effective heat capacity [J/K].

    Light elements rejected by the growing inner core mix through the
    outer core and release gravitational energy. Per unit cooling:
    ``[int_oc rho psi dV - M_oc psi(r_icb)] * alpha_c * Cc *
    |dr_icb/dT_cmb|`` with ``Cc = 4 pi r_icb^2 rho(r_icb) c_light /
    M_oc`` the enrichment per unit boundary advance in the
    complete-rejection limit; the structure of the Leeds
    ``thermal_history`` core model (Greenwood et al. 2021 lineage,
    Gubbins et al. 2003 formalism). Zero when either compositional
    parameter is zero, before onset, and after freeze-out. The
    outer-core region runs to ``upper`` (default the CMB radius),
    which the stratified budget reduces to the convecting radius.
    """
    p = self.profiles
    radius = self.r_icb(t_cmb)
    top = p.r_cmb if upper is None else upper
    # A layer below the ICB leaves no convecting shell: clamp the top at the ICB so the
    # integral closes to zero rather than flipping sign (a negative half-width would
    # negate both moments and pass the mass guard with a wrong-sign capacity).
    top = jnp.maximum(top, radius)

    def integrand(r):
        return p.density(r) * p.potential(r) * 4.0 * jnp.pi * r**2

    # Outer-core integrals on [r_icb, upper] with the shared GL panel.
    half = (top - radius) / 2.0
    centre = (top + radius) / 2.0
    s = centre[..., None] + half[..., None] * _GL_X
    rho_psi = half * jnp.sum(_GL_W * integrand(s), axis=-1)
    mass_oc = half * jnp.sum(_GL_W * p.density(s) * 4.0 * jnp.pi * s**2, axis=-1)
    potential_moment = rho_psi - mass_oc * p.potential(radius)

    enrichment = 4.0 * jnp.pi * radius**2 * p.density(radius) * self.c_light
    safe_mass = jnp.where(mass_oc > 0.0, mass_oc, 1.0)
    return (
        potential_moment
        * self.alpha_c
        * (enrichment / safe_mass)
        * self._boundary_sensitivity(t_cmb)
    )

def heat_content(self, t_cmb) -> float:
    """Core heat content [J] at CMB temperature ``t_cmb`` of the full, unstratified core.

    ``t_cmb * C_secular`` minus the latent and gravitational energy released
    between nucleation onset and ``t_cmb`` (zero above onset), so that
    ``heat_content(T2) - heat_content(T1)`` is the integral of the effective
    capacity from ``T1`` to ``T2``. The substitution ``T = T_onset - u^2``
    turns the square-root cusp of the latent and gravitational terms at onset,
    ``C ~ (T_onset - T)^(1/2)``, into a smooth integrand.
    Evaluated eagerly (``t_onset`` is a Python float).

    Raises
    ------
    ValueError
        For a stratified budget, whose capacity covers only the convecting
        volume and depends on the CMB heat flow.
    """
    if self.stratification:
        raise ValueError(
            'heat_content is the full-core content; a stratified budget evolves '
            'only its convecting volume'
        )
    secular = float(t_cmb) * float(self.secular_capacity())
    if self.capacity_mode == 'legacy' or float(t_cmb) >= self.t_onset:
        return secular
    u_max = float(_np.sqrt(self.t_onset - float(t_cmb)))
    boundary = jax.vmap(lambda t: self.latent_capacity(t) + self.gravitational_capacity(t))
    # One Gauss-Legendre panel per side of the freeze-out jump in the latent term.
    u_freeze = float(_np.sqrt(max(self.t_onset - self.t_freeze, 0.0)))
    edges = [0.0, u_freeze, u_max] if 0.0 < u_freeze < u_max else [0.0, u_max]
    released = 0.0
    for lo, hi in zip(edges[:-1], edges[1:]):
        u = lo + 0.5 * (hi - lo) * (_GL_X + 1.0)
        values = boundary(self.t_onset - u**2) * 2.0 * u
        released += 0.5 * (hi - lo) * float(jnp.sum(_GL_W * values))
    return secular - released

def latent_capacity(self, t_cmb):
    """Latent contribution to the effective heat capacity [J/K].

    ``L rho(r_icb) 4 pi r_icb^2 |dr_icb/dT_cmb|`` with the latent heat
    per unit mass ``L`` either the prescribed constant or
    ``T_icb * ds_fusion``. Zero before nucleation onset and after
    freeze-out completion when the boundary is pinned.
    """
    p = self.profiles
    radius = self.r_icb(t_cmb)
    if self.latent_heat is not None:
        heat = self.latent_heat
    else:
        heat = p.adiabat(radius, t_cmb) * self.ds_fusion
    area_mass = p.density(radius) * 4.0 * jnp.pi * radius**2
    return heat * area_mass * self._boundary_sensitivity(t_cmb)

def nucleation_factor(self, t_cmb):
    """Diagnostic activation in [0, 1]: sigmoid of centre subcooling."""
    return jax.nn.sigmoid(-self.superheat(0.0, t_cmb) / self.icn_width)

def r_icb(self, t_cmb):
    """Inner-core boundary radius [m] at ``t_cmb``.

    Fixed-iteration bisection of the superheat on [0, r_cmb]; with no
    crossing it converges to 0 (fully liquid) or r_cmb (fully frozen),
    so the value is always defined and trace-safe. A liquid centre
    returns zero regardless of frozen shells above it (the top-down
    and snow topologies have no inner core in the bottom-up sense);
    the boundary terms drop to zero there, and the solver refuses those
    states (``aragog.core.regime.refuse_unmodelled_regime``).

    The derivative is the implicit-function sensitivity attached as a
    custom JVP: a comparison-driven bisection carries no gradient of
    its own (autodiff sees only the converged constant trace), so
    without this rule anything differentiating through the boundary
    would silently read zero.
    """
    return self._r_icb_root(t_cmb)

def secular_capacity(self, upper=None):
    """Secular heat capacity dQ_s / d(dT_cmb/dt) [J/K].

    Profile mode: ``c_p int rho(r) f(r) 4 pi r^2 dr`` over the
    convecting volume (``upper`` defaults to the CMB radius), with
    ``f`` the adiabat shape ``T_a(r)/T_cmb``, i.e. ``c_p M
    tfac_eff``. Legacy mode: ``(4/3) pi r_cmb^3 rho c_p tfac``, the
    reservoir constant.
    """
    p = self.profiles
    if self.capacity_mode == 'legacy':
        volume = 4.0 / 3.0 * jnp.pi * p.r_cmb**3
        return volume * self.legacy_rho_core * p.c_p * self.legacy_tfac

    def integrand(r):
        shape = p.adiabat(r, 1.0)  # T_a / T_cmb: anchor-independent
        return p.density(r) * shape * 4.0 * jnp.pi * r**2

    top = p.r_cmb if upper is None else upper
    return p.c_p * self._quad_0_upper(top, integrand)

def superheat(self, r, t_cmb):
    """Adiabat minus melting curve [K] at radius ``r``; positive = liquid."""
    p = self.profiles
    return p.adiabat(r, t_cmb) - self.melting_curve.t_melt(p.pressure(r))

def adiabatic_heat_flow(self, t_cmb):
    """Heat conducted down the adiabat at the CMB, ``Qk`` [W]."""
    p = self.budget.profiles
    return -4.0 * jnp.pi * p.r_cmb**2 * self.k_core * p.adiabat_gradient(p.r_cmb, t_cmb)

def b_dipole_cmb(self, t_cmb, q_cmb, dipolarity: float = 1.0):
    """Dipole field at the CMB [T]: the core rms field times a
    dipolarity fraction; a potential-field continuation to the planet
    surface is the caller's ``(r_cmb / r_planet)^3`` factor."""
    return dipolarity * self.b_rms_core(t_cmb, q_cmb)

def b_rms_core(self, t_cmb, q_cmb):
    """Volume-averaged core field strength [T] from CHR09 Eq. 2.

    The reference flux is the superadiabatic part of the CMB heat flow
    over the CMB area; with nothing superadiabatic the field estimate
    is zero even where compositional convection keeps the entropy
    margin positive (a conservative floor, see the module docstring).
    """
    p = self.budget.profiles
    area = 4.0 * jnp.pi * p.r_cmb**2
    q_o = jnp.maximum(q_cmb - self.adiabatic_heat_flow(t_cmb), 0.0) / area
    mean_rho = p.enclosed_mass(p.r_cmb) / (4.0 / 3.0 * jnp.pi * p.r_cmb**3)
    energy_density = (
        _CHR09_C
        * self.f_ohm
        * mean_rho ** (1.0 / 3.0)
        * (self.chr09_efficiency_factor() * q_o) ** (2.0 / 3.0)
    )
    return jnp.sqrt(2.0 * MU0 * energy_density)

def chr09_efficiency_factor(self):
    """The printed Earth-core ``F``: ``0.88 (or 0.45) alpha g r / c_p``."""
    p = self.budget.profiles
    g_cmb = p.gravity(p.r_cmb)
    return _CHR09_F_GEOMETRY[self.flux_geometry] * p.alpha * g_cmb * p.r_cmb / p.c_p

def conduction_sink(self, upper=None):
    """Entropy sink of conduction along the adiabat, ``Ek`` [W/K].

    ``4 pi k int (dTa/dr / Ta)^2 r^2 dr`` over the convecting volume (``upper``
    defaults to the CMB radius), independent of ``T_cmb``: by quadrature, or in the
    small-radius mode, where the ratio is ``2 r / D^2``, as ``16 pi k upper^5 / (5 D^4)``.
    """
    p = self.budget.profiles
    top = p.r_cmb if upper is None else upper
    if p.adiabat_mode == 'small_radius':
        return 16.0 * jnp.pi * self.k_core * top**5 / (5.0 * p.d_scale**4)

    def integrand(r):  # d ln T / dr = -alpha g / c_p
        return (p.alpha * p.gravity(r) / p.c_p) ** 2 * r**2

    return 4.0 * jnp.pi * self.k_core * self.budget._quad_0_upper(top, integrand)

def entropy_margin(self, t_cmb, q_cmb, q_radio=0.0, *, t_shell=None):
    """Entropy production available to the dynamo, ``dE`` [W/K].

    The cooling rate follows from the energy budget for the given heat
    flow, the three cooling-proportional sources scale with it, and
    conduction subtracts. Positive margin sustains a dynamo. A stratified
    core needs its shell temperatures ``t_shell``: the convecting core and the shell
    each add their sources and sinks, all delivered at the temperature of the top shell
    cell (Greenwood et al. 2021, eq. 4).
    """
    b = self.budget
    if not b.stratification:
        cooling = -b.dtcmb_dt(t_cmb, q_cmb, q_sources=q_radio)
        capacity = (
            self.secular_entropy_capacity(t_cmb)
            + self.latent_entropy_capacity(t_cmb)
            + self.gravitational_entropy_capacity(t_cmb)
        )
        return (
            capacity * cooling
            + self.radiogenic_entropy(t_cmb, q_radio)
            - self.conduction_sink()
        )
    shell, upper = b.shell, b.r_convecting
    d_core, d_shell = b.core_rates(t_cmb, t_shell, q_cmb, q_radio)
    t_top = shell.top_temperature(t_shell)
    capacity = (
        self.secular_entropy_capacity(t_cmb, upper=upper, t_ref=t_top)
        + self.latent_entropy_capacity(t_cmb, t_ref=t_top)
        + self.gravitational_entropy_capacity(
            t_cmb, upper=shell.layer_base(t_shell, t_cmb), t_ref=t_top
        )
    )
    heating = q_radio / b.profiles.enclosed_mass(b.profiles.r_cmb)
    return (
        -capacity * d_core
        + self.radiogenic_entropy(t_cmb, q_radio, upper=upper, t_ref=t_top)
        - self.conduction_sink(upper=upper)
        + shell.entropy_rate(t_shell, d_shell, t_cmb, heating, t_top)
    )

def gravitational_entropy_capacity(self, t_cmb, upper=None, *, t_ref=None):
    """Entropy per unit cooling from gravitational energy [J/K^2].

    Gravitational energy dissipates in full within the convecting
    core: ``C_grav / T_ref`` (default ``T_cmb``).
    """
    ref = t_cmb if t_ref is None else t_ref
    return self.budget.gravitational_capacity(t_cmb, upper=upper) / ref

def latent_entropy_capacity(self, t_cmb, *, t_ref=None):
    """Entropy per unit cooling from latent heat [J/K^2].

    The latent capacity released at ``T_icb`` and delivered at
    ``T_ref`` (default ``T_cmb``): ``C_lat (T_icb - T_ref) / (T_icb T_ref)``.
    """
    b = self.budget
    radius = b.r_icb(t_cmb)
    t_icb = b.profiles.adiabat(radius, t_cmb)
    ref = t_cmb if t_ref is None else t_ref
    return b.latent_capacity(t_cmb) * (t_icb - ref) / (t_icb * ref)

def radiogenic_entropy(self, t_cmb, q_radio, upper=None, *, t_ref=None):
    """Entropy rate from internal heating [W/K].

    ``int h rho (1/T_ref - 1/Ta) dV`` with the heating rate per unit
    mass ``h = q_radio / M_core`` uniform, integrated over the
    convecting volume (``upper`` defaults to the CMB radius, ``T_ref`` to ``T_cmb``).
    """
    p = self.budget.profiles
    mass = p.enclosed_mass(p.r_cmb)
    ref = t_cmb if t_ref is None else t_ref

    def integrand(r):
        inv_gap = 1.0 / ref - 1.0 / p.adiabat(r, t_cmb)
        return p.density(r) * inv_gap * 4.0 * jnp.pi * r**2

    top = p.r_cmb if upper is None else upper
    return (q_radio / mass) * self.budget._quad_0_upper(top, integrand)

def secular_entropy_capacity(self, t_cmb, upper=None, *, t_ref=None):
    """Entropy per unit cooling from secular cooling [J/K^2].

    ``int rho c_p (Ta/T_ref - 1) dV / T_cmb``: heat extracted at
    temperature ``Ta`` and delivered at ``T_ref`` (default ``T_cmb``) produces entropy in
    proportion to the temperature drop. Integrated over the
    convecting volume (``upper`` defaults to the CMB radius).
    """
    p = self.budget.profiles
    scale = 1.0 if t_ref is None else t_cmb / t_ref

    def integrand(r):
        shape = p.adiabat(r, 1.0)
        return p.density(r) * p.c_p * (scale * shape - 1.0) * 4.0 * jnp.pi * r**2

    top = p.r_cmb if upper is None else upper
    return self.budget._quad_0_upper(top, integrand) / t_cmb

def diagnostics(self, q_cmb: float | None = None) -> dict:
    """Current core state as plain floats.

    Always contains ``t_cmb``, ``t_cen``, ``r_icb``,
    ``nucleation_factor``, ``effective_capacity``, and the phase flags
    ``inner_core_present`` / ``fully_frozen``. With an entropy budget
    attached and ``q_cmb`` given, adds ``entropy_margin``,
    ``dynamo_active``, and ``b_rms_core``.
    """
    b = self.budget
    t = self.t_cmb
    radius = float(b.r_icb(t))
    out = {
        't_cmb': t,
        't_cen': float(b.profiles.t_cen(t)),
        'r_icb': radius,
        'nucleation_factor': float(b.nucleation_factor(t)),
        'effective_capacity': float(b.effective_capacity(t)),
        'inner_core_present': bool(radius > 0.0),
        'fully_frozen': bool(~b._liquid_remains(t)),
        'regime': int(crystallization_regime(b, t)),
    }
    if self.entropy is not None and q_cmb is not None:
        margin = float(self.entropy.entropy_margin(t, q_cmb, self.q_radio))
        out['entropy_margin'] = margin
        out['dynamo_active'] = bool(margin > 0.0)
        out['b_rms_core'] = float(self.entropy.b_rms_core(t, q_cmb))
    return out

def step(self, q_cmb: float, dt: float, q_cmb_end: float | None = None) -> float:
    """Advance ``T_cmb`` over ``dt`` seconds of heat flow ``q_cmb`` [W].

    The heat flow ramps linearly from ``q_cmb`` to ``q_cmb_end`` (equal
    to ``q_cmb`` when omitted) across the interval, and the state
    advances with ``n_substeps`` classical RK4 sub-steps. The sub-step
    trajectory is stored on ``last_times`` / ``last_t_cmb`` /
    ``last_q_cmb`` so the caller can see through the coupling interval
    instead of only its endpoints. Returns the new ``T_cmb`` [K].

    Raises
    ------
    ValueError
        If ``dt`` is not positive.
    """
    if not float(dt) > 0.0:
        raise ValueError(f'dt must be positive, got {dt}')
    q_end = float(q_cmb) if q_cmb_end is None else float(q_cmb_end)
    h = float(dt) / self.n_substeps

    def q_at(t: float) -> float:
        return float(q_cmb) + (q_end - float(q_cmb)) * (t / float(dt))

    def rhs(t: float, temp):
        return self._rhs(temp, q_at(t), self.q_radio)

    times = [0.0]
    temps = [self.t_cmb]
    temp = jnp.asarray(self.t_cmb)
    t = 0.0
    for _ in range(self.n_substeps):
        k1 = rhs(t, temp)
        k2 = rhs(t + h / 2.0, temp + h / 2.0 * k1)
        k3 = rhs(t + h / 2.0, temp + h / 2.0 * k2)
        k4 = rhs(t + h, temp + h * k3)
        temp = temp + h / 6.0 * (k1 + 2.0 * k2 + 2.0 * k3 + k4)
        t += h
        times.append(t)
        temps.append(float(temp))
    self.last_times = jnp.asarray(times)
    self.last_t_cmb = jnp.asarray(temps)
    self.last_q_cmb = jnp.asarray([q_at(x) for x in times])
    self.t_cmb = float(temp)
    return self.t_cmb

def adiabat(self, r, t_cmb):
    """Adiabatic temperature [K] at radius ``r`` anchored at ``t_cmb``; see ``adiabat_mode``.

    Hotter inward, equal to ``t_cmb`` at the CMB by construction.
    """
    r = jnp.asarray(r)
    if self.adiabat_mode == 'exact':
        return t_cmb * jnp.exp(-self.alpha * self.potential(r) / self.c_p)
    return t_cmb * jnp.exp((self.r_cmb**2 - r**2) / self.d_scale**2)

def adiabat_gradient(self, r, t_cmb):
    """``dT/dr`` [K/m] of the adiabat, ``-alpha g T / c_p``, with the gravity of the mode."""
    r = jnp.asarray(r)
    g = (
        self.gravity(r)
        if self.adiabat_mode == 'exact'
        else 4.0 * jnp.pi / 3.0 * G * self.rho_cen * r
    )
    return -self.alpha * g * self.adiabat(r, t_cmb) / self.c_p

def density(self, r):
    """Density [kg m-3] at radius ``r`` [m]."""
    r = jnp.asarray(r)
    return self.rho_cen * jnp.exp(-((r / self.length_scale) ** 2))

def enclosed_mass(self, r):
    """Mass [kg] inside radius ``r``: the closed-form Gaussian integral.

    ``4 pi rho_cen [ (sqrt(pi) L^3 / 4) erf(r/L) - (L^2 r / 2)
    exp(-r^2/L^2) ]``, the exact antiderivative of
    ``4 pi rho(s) s^2``.
    """
    r = jnp.asarray(r)
    length = self.length_scale
    x = r / length
    integral = (_SQRT_PI * length**3 / 4.0) * erf(x) - (length**2 * r / 2.0) * jnp.exp(
        -(x**2)
    )
    return 4.0 * jnp.pi * self.rho_cen * integral

@classmethod
def from_structure(
    cls,
    *,
    m_core: float,
    p_cen: float,
    r_cmb: float,
    p_cmb: float,
    alpha: float,
    c_p: float,
    pressure_mode: str = 'quadrature',
    adiabat_mode: str = 'exact',
) -> GaussianCoreProfiles:
    """Fit central density and length scale to core mass and central pressure.

    Parameters
    ----------
    m_core : float
        Total core mass [kg], positive.
    p_cen : float
        Central pressure [Pa], must exceed ``p_cmb``.
    r_cmb : float
        Core-mantle boundary radius [m], positive.
    p_cmb : float
        Pressure at the core-mantle boundary [Pa], positive.
    alpha : float
        Thermal expansion coefficient [K-1], positive.
    c_p : float
        Isobaric specific heat capacity [J kg-1 K-1], positive.
    pressure_mode : str, optional
        Pressure mode: ``'quadrature'`` (default) or ``'labrosse'``.
    adiabat_mode : str, optional
        Adiabat mode: ``'exact'`` (default) or ``'small_radius'``.

    Returns
    -------
    GaussianCoreProfiles
        Profile instance with fitted ``rho_cen`` and ``length_scale``.

    Raises
    ------
    ValueError
        If parameters are non-positive, if ``p_cen <= p_cmb``, if the
        central pressure is below the incompressible sphere limit, or if the
        root solve does not converge within the Gaussian family regime.
    """
    params = {
        'm_core': m_core,
        'p_cen': p_cen,
        'r_cmb': r_cmb,
        'p_cmb': p_cmb,
        'alpha': alpha,
        'c_p': c_p,
    }
    for name, value in params.items():
        if not float(value) > 0.0:
            raise ValueError(f'{name} must be positive, got {value}')

    m_core = float(m_core)
    p_cen = float(p_cen)
    r_cmb = float(r_cmb)
    p_cmb = float(p_cmb)
    alpha = float(alpha)
    c_p = float(c_p)

    if p_cen <= p_cmb:
        raise ValueError(f'p_cen ({p_cen:.4e} Pa) must exceed p_cmb ({p_cmb:.4e} Pa)')

    # Uniform density sphere has lowest central pressure for given mass.
    vol = 4.0 / 3.0 * _np.pi * r_cmb**3
    rho_avg = m_core / vol
    p_incomp = p_cmb + (2.0 / 3.0) * _np.pi * G * rho_avg**2 * r_cmb**2
    if p_cen < p_incomp:
        raise ValueError(
            f'p_cen ({p_cen:.4e} Pa) is below the incompressible central pressure '
            f'({p_incomp:.4e} Pa)'
        )

    def build(rho_c: float, length: float) -> GaussianCoreProfiles:
        return cls(
            rho_cen=rho_c,
            length_scale=length,
            r_cmb=r_cmb,
            p_cmb=p_cmb,
            alpha=alpha,
            c_p=c_p,
            pressure_mode=pressure_mode,
            adiabat_mode=adiabat_mode,
        )

    def rho_for(length: float) -> float:
        return float(m_core / float(build(1.0, length).enclosed_mass(r_cmb)))

    def residual(length: float) -> float:
        return float(build(rho_for(length), length).pressure(0.0)) - p_cen

    # Enforce the Gaussian family validity regime (r_cmb < 3 * length_scale).
    l_min = r_cmb / 2.999
    l_max = 100.0 * r_cmb

    res_min = residual(l_min)
    res_max = residual(l_max)
    if res_min < 0.0:
        p_max = float(p_cen + res_min)
        raise ValueError(
            f'p_cen ({p_cen:.4e} Pa) exceeds maximum central pressure ({p_max:.4e} Pa) '
            'achievable within valid Gaussian regime (r_cmb < 3 L)'
        )
    if res_min * res_max > 0.0:
        raise ValueError(
            f'Gaussian core profile fit did not bracket a root for m_core={m_core:.4e} kg, '
            f'p_cen={p_cen:.4e} Pa (residual span [{res_min:.4e}, {res_max:.4e}])'
        )

    sol = root_scalar(
        residual,
        bracket=[l_min, l_max],
        method='brentq',
        xtol=1e-8,
        rtol=1e-10,
    )
    if not sol.converged:
        raise ValueError(
            'Gaussian core profile fit to M_core and P_cen did not converge: ' + str(sol)
        )

    l_fit = float(sol.root)
    return build(rho_for(l_fit), l_fit)

def gravity(self, r):
    """Gravitational acceleration magnitude [m s-2] at radius ``r``.

    Exact for the Gaussian density: ``G M(r) / r^2``, with the removable
    singularity at the centre replaced by its analytic limit
    ``(4 pi / 3) G rho_cen r``.
    """
    r = jnp.asarray(r)
    small = self.length_scale * 1e-6
    safe_r = jnp.where(r > small, r, small)
    exact = G * self.enclosed_mass(safe_r) / safe_r**2
    centre = 4.0 * jnp.pi / 3.0 * G * self.rho_cen * r
    return jnp.where(r > small, exact, centre)

def potential(self, r):
    """Gravitational potential [J/kg] at radius ``r``, zero at the CMB.

    ``psi(r) = -int_r^{r_cmb} g(s) ds`` on the same fixed 32-point
    Gauss-Legendre panel as the pressure; negative inside the core,
    which is the reference the gravitational-energy budget term uses.
    """
    r = jnp.asarray(r)
    half_span = (self.r_cmb - r) / 2.0
    centre = (self.r_cmb + r) / 2.0
    s = centre[..., None] + half_span[..., None] * _GL_X
    integral = half_span * jnp.sum(_GL_W * self.gravity(s), axis=-1)
    return -integral

def pressure(self, r):
    """Pressure [Pa] at radius ``r``, from hydrostatic balance.

    Quadrature mode: ``P(r) = p_cmb + int_r^{r_cmb} rho(s) g(s) ds`` on
    a fixed 32-point Gauss-Legendre panel, jit-safe with no adaptive
    control flow. Labrosse mode: the printed closed form.
    """
    if self.pressure_mode == 'labrosse':
        return self._pressure_labrosse(r)
    r = jnp.asarray(r)
    half_span = (self.r_cmb - r) / 2.0
    centre = (self.r_cmb + r) / 2.0
    # Broadcast quadrature nodes over any leading shape of r.
    s = centre[..., None] + half_span[..., None] * _GL_X
    integrand = self.density(s) * self.gravity(s)
    integral = half_span * jnp.sum(_GL_W * integrand, axis=-1)
    return self.p_cmb + integral

def t_cen(self, t_cmb):
    """Centre temperature [K] on the adiabat anchored at ``t_cmb``."""
    return self.adiabat(0.0, t_cmb)

def t_melt(self, pressure, light_element_fraction=None):
    """Alloy melting temperature [K] at ``pressure`` [Pa].

    ``light_element_fraction`` overrides the instance value when given,
    which allows callers to evolve the alloy composition without
    rebuilding the curve. A concrete override (a real scalar or array) is
    taken in float64 and checked like the instance value (in ``[0, 1)``, with
    the depressed curve kept positive); a traced override is the caller's
    responsibility, since a trace cannot raise on data.

    Raises
    ------
    ValueError
        If a concrete override is not a real number (a bool, string,
        bytes or complex value included), is outside ``[0, 1)``, or
        drives the depression factor to zero or below.
    """
    if light_element_fraction is None:
        x = self.light_element_fraction
    else:
        x = light_element_fraction
        if not isinstance(x, jax.core.Tracer):
            x = _checked_fraction(x, self.depression)
    return self.t_melt_pure(pressure) * (1.0 - self.depression * x)

@staticmethod
def t_melt_pure(pressure):
    """Pure-iron melting temperature [K] at ``pressure`` [Pa]: the two Anzellini
    branches joined by the smootherstep blend of the module docstring."""
    p = jnp.asarray(pressure)
    p_gpa = p / 1e9
    low = _T0 * ((p_gpa - _P0 / 1e9) / _DP_LOW + 1.0) ** _EXP_LOW
    high = _TT * ((p_gpa - _PT / 1e9) / _DP_HIGH + 1.0) ** _EXP_HIGH
    u = jnp.clip(
        0.5 * ((p - _PT) / IRON_MELTING_BRANCH_BLEND_HALF_WIDTH_PA + 1.0),
        0.0,
        1.0,
    )
    w = u * u * u * (u * (u * 6.0 - 15.0) + 10.0)
    return (1.0 - w) * low + w * high