"""Gravity-only envelope cost from two unchanged native stellar profiles.

Mass and radius are original MESA grid samples in solar units. The H-abundance
threshold defines an operational cut, not a unique future remnant. No internal
energy is reconstructed from temperature or inferred from a missing EOS column.
"""
from io import StringIO
from pathlib import Path
import shlex

import numpy as np
from scipy.integrate import quad

GRAVITATIONAL_CONSTANT_CGS = 6.67430e-8
SOLAR_MASS_CGS = 1.98847e33
SOLAR_RADIUS_CGS = 6.957e10
ENERGY_UNIT_ERG = GRAVITATIONAL_CONSTANT_CGS*SOLAR_MASS_CGS**2/SOLAR_RADIUS_CGS


def read_native(path: Path):
    lines = path.read_text().splitlines()
    header = dict(zip(shlex.split(lines[1]), shlex.split(lines[2])))
    rows = np.genfromtxt(StringIO("\n".join(lines[5:])), names=True)
    order = np.argsort(rows["mass"])
    rows = rows[order]
    radius = rows["radius"] if "radius" in rows.dtype.names else 10**rows["logR"]
    hydrogen = rows["h1"] if "h1" in rows.dtype.names else rows["x_mass_fraction_H"]
    return header, rows["mass"], radius, hydrogen, rows["zone"]


def envelope_cost(path: Path, threshold: float, stride=1):
    header, mass, radius, hydrogen, native_zone = read_native(path)
    candidates = np.flatnonzero(hydrogen < threshold)
    if not len(candidates) or candidates[-1] == len(mass)-1:
        raise ValueError("No interior H-depleted core at this operational threshold")
    cut = candidates[-1]
    indices = np.unique(np.r_[np.arange(cut, len(mass), stride), len(mass)-1])
    m, r = mass[indices], radius[indices]
    if np.any(np.diff(m) <= 0) or np.any(r <= 0):
        raise ValueError("Invalid native radial coordinate")
    # Adaptive integration for each interpolated native interval, independently
    # checked against a global composite trapezoid and deliberate grid coarsening.
    adaptive = sum(quad(lambda coordinate: coordinate/np.interp(coordinate, m, r),
                        float(left), float(right), epsabs=2e-13, epsrel=2e-12)[0]
                   for left, right in zip(m[:-1], m[1:]))
    trapezoid = np.trapezoid(m/r, m)
    donor = float(header["star_mass"])
    photosphere = float(header["photosphere_r"])
    core = float(mass[cut])
    return {"native_mass_solar": donor, "native_photospheric_radius_solar": photosphere,
            "hydrogen_threshold": threshold, "core_mass_solar": core,
            "cut_radius_solar": float(radius[cut]), "native_zones": len(mass),
            "native_cut_identity": {
                "rule": "Last original native zone with X_H below threshold; no interpolated cut used for integration",
                "below": {"zone": int(native_zone[cut]), "X_H": float(hydrogen[cut]),
                          "mass_solar": float(mass[cut]), "radius_solar": float(radius[cut])},
                "next_above": {"zone": int(native_zone[cut+1]), "X_H": float(hydrogen[cut+1]),
                               "mass_solar": float(mass[cut+1]), "radius_solar": float(radius[cut+1])},
                "mass_bracket_width_solar": float(mass[cut+1]-mass[cut]),
                "linear_threshold_mass_for_sensitivity_only": float(np.interp(threshold,
                    hydrogen[cut:cut+2], mass[cut:cut+2])),
                "precision_note": "Integration accuracy does not establish threshold location or a unique physical remnant boundary"},
            "envelope_nodes_used": len(m), "gravity_cost_units": adaptive,
            "gravity_cost_erg": adaptive*ENERGY_UNIT_ERG,
            "trapezoid_relative_difference": float(trapezoid/adaptive-1),
            "lambda_gravity": donor*(donor-core)/(photosphere*adaptive),
            "internal_energy_available": False,
            "structure_clock_note": "Native header age; a saved/restarted model's clock is not assumed to be total stellar age"}


def ce_final_separation(donor, core, companion, initial_separation, gravity_cost, alpha):
    if (not all(np.isfinite(value) for value in [donor, core, companion, initial_separation, gravity_cost, alpha])
            or not 0 < core < donor or companion <= 0 or alpha <= 0 or gravity_cost <= 0 or initial_separation <= 0):
        raise ValueError("Outside the positive gravity-only endpoint domain")
    return 1/(donor/(core*initial_separation)+2*gravity_cost/(alpha*core*companion))