"""Independent native segment primitives, analytic limits and exact-volume failures."""
from pathlib import Path
import argparse
import json
import numpy as np
from scipy.integrate import quad
from scipy.optimize import brentq
from envelope_reference import read_native, envelope_cost, ce_final_separation
from orbit_reference import integrate_transfer, constant_gamma_closed
from stability_reference import run_stability_reference


def segment_primitive(mass, radius):
    interval = np.diff(mass); ratio = np.diff(radius)/radius[:-1]
    zeroth = np.empty(len(ratio)); first = np.empty(len(ratio)); small = abs(ratio) < 1e-4
    s = ratio[small]
    zeroth[small] = 1-s/2+s*s/3-s**3/4+s**4/5
    first[small] = .5-s/3+s*s/4-s**3/5+s**4/6
    s = ratio[~small]
    zeroth[~small] = np.log1p(s)/s; first[~small] = (s-np.log1p(s))/(s*s)
    return float(np.sum(interval/radius[:-1]*(mass[:-1]*zeroth+interval*first)))


def validate(inputs, output):
    output.mkdir(parents=True, exist_ok=True); native = []
    for member in json.loads(Path('native-members.json').read_text()):
        path = inputs/member['profile_file']; _, mass, radius, hydrogen, _ = read_native(path)
        for threshold in [.001, .01, .1]:
            original = envelope_cost(path, threshold); cut = np.flatnonzero(hydrogen < threshold)[-1]
            primitive = segment_primitive(mass[cut:], radius[cut:]); difference = abs(primitive/original['gravity_cost_units']-1)
            assert difference < 1e-10, difference
            donor, core = original['native_mass_solar'], original['core_mass_solar']; companion = 1.4; initial = 100.
            endpoint = ce_final_separation(donor, core, companion, initial, original['gravity_cost_units'], .5)
            root = brentq(lambda separation: .5*(core*companion/(2*separation)-donor*companion/(2*initial))-original['gravity_cost_units'], 1e-10, initial, xtol=1e-13)
            assert abs(endpoint/root-1) < 1e-10
            low_cost = ce_final_separation(donor, core, companion, initial, 1e-16, 1)
            assert abs(low_cost/((core/donor)*initial)-1) < 1e-12
            native.append({'member': member['name'], 'threshold': threshold, 'piecewise_analytic_vs_adaptive_relative_error': difference,
                           'scalar_root_vs_endpoint_relative_error': abs(endpoint/root-1), 'native_cut_identity': original['native_cut_identity']})
    wind = integrate_transfer(3, 2, 20, 1.5, 0, 'donor-wind')
    wind_error = float(np.max(abs(wind['separation']*(wind['donor']+wind['accretor'])/(20*5)-1))); assert wind_error < 1e-9
    removed = np.linspace(0,1.5,301); separation = constant_gamma_closed(3,2,1,removed,1,2)
    assert abs(removed[np.argmin(separation)]-.5) < 1e-12 and abs(np.min(separation)-.9216) < 1e-12
    sphere = quad(lambda mass: mass/mass**(1/3), .2, 1, epsabs=1e-13)[0]
    assert abs(sphere/(3/5*(1-.2**(5/3)))-1) < 1e-12
    stability = run_stability_reference()
    result = {'status': 'PASS independently recomputed native segment primitives, scalar roots, analytic limits and volume-response fit failures',
              'native_integrals': native, 'slow_wind_a_times_mass_maximum_relative_error': wind_error,
              'conservative_minimum_removed_solar_mass': .5, 'conservative_minimum_separation_ratio': .9216,
              'exact_volume_stability': stability,
              'boundary': 'Numerical integration does not locate a unique physical core; finite angular refinements are not a formal certificate. These are fixed mathematical/native-input checks, not MESA evolution, real CE calibration or fluid validation.'}
    (output/'independent-limits.json').write_text(json.dumps(result, indent=2)+'\n')
    print('PASS native integrals, closed limits,20volume responses and2retained classification failures')


if __name__ == '__main__':
    parser = argparse.ArgumentParser(); parser.add_argument('--inputs', type=Path, default=Path('inputs')); parser.add_argument('--output', type=Path, default=Path('results')); args = parser.parse_args(); validate(args.inputs, args.output)