"""Frozen-model Mercury check: local high-field tail and leading solar tide.

First-order Kepler orbit averaging, not an ephemeris fit or metric theory.
The orientation envelope covers every direction of the Galactic quadrupole.
No gravity parameter is adjusted using Mercury.
"""
from __future__ import annotations

import hashlib
import json
from pathlib import Path

import numpy as np
from scipy.special import roots_legendre

from dark_transition_tweak import REFERENCE

BASE = Path(__file__).resolve().parent
OUT = BASE / "outputs/dark_mercury"
PUBLICATION = BASE / "outputs/dark_transition_tweak/publication_summary.json"
AU = 149597870700.0
LIGHT_SPEED = 299792458.0
JULIAN_CENTURY = 36525 * 86400.0
RAD_S_TO_ARCSEC_CENTURY = (180 / np.pi) * 3600 * JULIAN_CENTURY
# Nominal solar GM for an illustrative test-particle orbit, not an ephemeris
# estimate. The existing solar Q2 was computed with a slightly different GM.
GM = 1.3271244e20
# JPL approximate mean elements, J2000 ecliptic, Table 1 (1800--2050).
SEMIMAJOR = .38709927 * AU
ECCENTRICITY = .20563593
INCLINATION = np.deg2rad(7.00497902)
ARG_PERIHELION = np.deg2rad(77.45779628 - 48.33076593)
MEAN_MOTION = np.sqrt(GM / SEMIMAJOR**3)


def sha256(path):
    return hashlib.sha256(Path(path).read_bytes()).hexdigest()


def orbit(nodes):
    """Perifocal coordinates: p toward perihelion, q in plane, h normal."""
    z, w = roots_legendre(nodes)
    anomaly = np.pi * (z + 1)
    co, si = np.cos(anomaly), np.sin(anomaly)
    s = np.sqrt(1 - ECCENTRICITY**2)
    factor = 1 - ECCENTRICITY * co
    r = SEMIMAJOR * np.column_stack((co - ECCENTRICITY, s * si, 0 * si))
    v = SEMIMAJOR * MEAN_MOTION / factor[:, None] * np.column_stack((-si, s * co, 0 * si))
    radius = SEMIMAJOR * factor
    radial = r / radius[:, None]
    tangent = np.column_stack((-radial[:, 1], radial[:, 0], 0 * si))
    return dict(r=r, v=v, radius=radius, radial=radial, tangent=tangent,
                # dt / period = (1-e*cos(E)) dE / (2*pi).
                weights=w * factor / 2)


def secular_rates(sample, acceleration):
    """Gauss equations, retaining the node term in longitude varpi=omega+Omega."""
    e = ECCENTRICITY
    s = np.sqrt(1 - e**2)
    radius = sample["radius"]
    cf, sf = sample["radial"][:, 0], sample["radial"][:, 1]
    radial = np.einsum("ij,ij->i", acceleration, sample["radial"])
    tangent = np.einsum("ij,ij->i", acceleration, sample["tangent"])
    normal = acceleration[:, 2]
    in_plane = s / (MEAN_MOTION * SEMIMAJOR * e) * (
        -radial * cf + (1 + radius / (SEMIMAJOR * s**2)) * tangent * sf)
    sin_u = np.sin(ARG_PERIHELION) * cf + np.cos(ARG_PERIHELION) * sf
    node = radius * sin_u * normal / (
        MEAN_MOTION * SEMIMAJOR**2 * s * np.sin(INCLINATION))
    average = lambda value: float(np.dot(sample["weights"], value))
    return {"in_plane_apsidal_rad_s": average(in_plane),
            "node_rad_s": average(node),
            "longitude_perihelion_rad_s": average(in_plane + (1 - np.cos(INCLINATION)) * node)}


def quadrupole_acceleration(sample, direction, q2):
    # delta Phi = -Q2/2 * [(k.r)^2-r^2/3]; acceleration = -grad(delta Phi).
    direction = np.asarray(direction, dtype=float)
    direction = direction / np.linalg.norm(direction)
    r = sample["r"]
    return q2 * ((r @ direction)[:, None] * direction - r / 3)


def orientation_matrix():
    """Exactly averaged longitude rate = (Q2/n) k.T @ B @ k.

    For an in-plane tide this reduces to the standard expression
    Q2*sqrt(1-e^2)/(4*n) * (1+5*cos(2*theta)). The off-diagonal
    terms retain the 3D nodal correction for the J2000 ecliptic longitude.
    """
    e = ECCENTRICITY
    s = np.sqrt(1 - e**2)
    b = np.diag([1.5*s, -s, -.5*s])
    prefactor = np.tan(INCLINATION / 2) / (4*s)
    b[0, 2] = b[2, 0] = prefactor * (1 + 4*e**2) * np.sin(ARG_PERIHELION)
    b[1, 2] = b[2, 1] = prefactor * (1 - e**2) * np.cos(ARG_PERIHELION)
    return b


def radial_tail(sample, candidate):
    """Leading perturbative spherical inversion of g*mu(g/gdark)=GM/r^2.

    Evaluate mu-1 directly: subtracting the implemented floating-point mu
    from one would incorrectly give zero at Mercury. Corrections from
    evaluating the tail at g rather than g_N are second order in mu-1.
    """
    g_newton = GM / sample["radius"]**2
    x = g_newton / candidate["galaxy_scale_m_s2"]
    inv_t = candidate["position"] / x
    original = -REFERENCE.c / x**3 + 2*REFERENCE.normalization*REFERENCE.d / x**4
    extension = (27/4)*candidate["amplitude"] * inv_t**4 / (1 + inv_t**2)**3
    # Positive radial acceleration points away from the Sun.
    def rate(epsilon):
        acceleration = (g_newton * epsilon)[:, None] * sample["radial"]
        return secular_rates(sample, acceleration)["longitude_perihelion_rad_s"] * RAD_S_TO_ARCSEC_CENTURY
    return {"x_min": float(x.min()), "x_max": float(x.max()),
            "max_abs_mu_minus_one": float(np.max(np.abs(original + extension))),
            "max_abs_radial_perturbation_m_s2": float(np.max(np.abs(g_newton * (original + extension)))),
            "original_tail_at_selected_scale_arcsec_century": rate(original),
            "added_term_arcsec_century": rate(extension),
            "total_arcsec_century": rate(original + extension)}


def main():
    publication = json.loads(PUBLICATION.read_text())
    candidate = publication["candidate"]
    assert (candidate["amplitude"], candidate["position"]) == (-.45, 4.25)
    # Follow the existing reports' provenance instead of rewriting them.
    for name, expected in publication["input_sha256"].items():
        assert sha256(BASE / name) == expected, name
    verification = json.loads((BASE / "outputs/dark_transition_tweak/Am0p45_T4p25_verification.json").read_text())
    bands = json.loads((BASE / "outputs/dark_summary.json").read_text())
    solar = json.loads((BASE / "outputs/dark_solar_summary.json").read_text())
    expected_sources = {"dark_transition_tweak.py": verification["input_sha256"]["dark_transition_tweak.py"],
                        "dark_band_report.py": bands["python_sources_sha256"]["dark_band_report.py"],
                        "dark_solar_report.py": solar["report_source_sha256"]}
    for name, expected in expected_sources.items():
        assert sha256(BASE / name) == expected, name
    b = orientation_matrix()
    eigenvalues, eigenvectors = np.linalg.eigh(b)
    directions = [np.array(v, dtype=float) for v in [(1, 0, 0), (0, 1, 0), (0, 0, 1), (2, 3, 5), (-1, 4, 2)]]
    directions += [eigenvectors[:, i] for i in range(3)]
    controls = []
    tails = []
    gr_formula = 3 * MEAN_MOTION * GM / (
        LIGHT_SPEED**2 * SEMIMAJOR * (1 - ECCENTRICITY**2))
    for nodes in [128, 256]:
        sample = orbit(nodes)
        tide_errors = []
        for direction in directions:
            direction = direction / np.linalg.norm(direction)
            numeric = secular_rates(sample, quadrupole_acceleration(sample, direction, 1e-27))["longitude_perihelion_rad_s"]
            expected = 1e-27 / MEAN_MOTION * float(direction @ b @ direction)
            tide_errors.append(abs(numeric - expected) / (1e-27 / MEAN_MOTION))
        r, v, radius = sample["r"], sample["v"], sample["radius"]
        gr_acceleration = GM / (LIGHT_SPEED**2 * radius[:, None]**3) * (
            (4*GM/radius - np.einsum("ij,ij->i", v, v))[:, None] * r
            + 4*np.einsum("ij,ij->i", r, v)[:, None] * v)
        gr_numeric = secular_rates(sample, gr_acceleration)["longitude_perihelion_rad_s"]
        zero = secular_rates(sample, np.zeros_like(r))
        assert max(map(abs, zero.values())) == 0
        assert max(tide_errors) < 1e-10
        assert abs(gr_numeric / gr_formula - 1) < 1e-10
        controls.append({"nodes": nodes, "quadrupole_max_error_in_units_Q2_over_n": max(tide_errors),
                         "GR_orbit_average_relative_error": abs(gr_numeric / gr_formula - 1),
                         "Newtonian_extra_precession_zero": True})
        tails.append(radial_tail(sample, candidate))
    tail_change = abs(tails[1]["total_arcsec_century"] - tails[0]["total_arcsec_century"])
    assert tail_change < 1e-25
    cases = []
    for row in publication["external_field_samples"]:
        interval = np.sort(row["Q2_s_minus2"] / MEAN_MOTION * eigenvalues * RAD_S_TO_ARCSEC_CENTURY)
        cases.append({**row,
                      "longitude_precession_orientation_interval_arcsec_century": [float(interval[0]), float(interval[-1])],
                      "max_abs_longitude_precession_arcsec_century": float(np.max(np.abs(interval)))})
    report = {
        "scope": "Frozen phenomenological LQ8 extension: first-order radial tail plus leading solar quadrupole; no metric prediction or ephemeris fit",
        "date": "2026-09-05",
        "model": {k: candidate[k] for k in ["name", "amplitude", "position", "galaxy_scale_m_s2"]},
        "orbit": {"semimajor_axis_au": SEMIMAJOR / AU, "eccentricity": ECCENTRICITY,
                  "inclination_deg": float(np.rad2deg(INCLINATION)),
                  "argument_perihelion_deg": float(np.rad2deg(ARG_PERIHELION)),
                  "frame": "J2000 ecliptic; longitude varpi=omega+Omega",
                  "elements_scope": "JPL approximate mean elements; illustrative frozen test-particle ellipse",
                  "GM_m3_s2": GM, "mean_motion_rad_s": float(MEAN_MOTION),
                  "solar_Q2_solver_GM_m3_s2": 6.67430e-11 * 1.98847e30,
                  "solar_Q2_GM_convention_relative_difference": (6.67430e-11*1.98847e30)/GM - 1},
        "GR_reference_arcsec_century": float(gr_formula * RAD_S_TO_ARCSEC_CENTURY),
        "PPN_reference_multiplier": "(2+2*gamma-beta)/3; gamma=beta=1 in GR; not derived for LQ8",
        "radial_tail": tails[-1],
        "orientation_matrix_B": b.tolist(),
        "orientation_eigenvalues": eigenvalues.tolist(),
        "quadrupole_cases": cases,
        "largest_quadrupole_envelope_arcsec_century": max(row["max_abs_longitude_precession_arcsec_century"] for row in cases),
        "controls": controls,
        "radial_tail_quadrature_change_arcsec_century": tail_change,
        "tail_constants": {"c3": REFERENCE.c, "normalization": REFERENCE.normalization, "D": REFERENCE.d},
        "source_sha256": sha256(__file__),
        "input_sha256": {str(path.relative_to(BASE)): sha256(path) for path in [PUBLICATION,
            BASE / "dark_band_report.py", BASE / "dark_solar_report.py", BASE / "dark_transition_tweak.py",
            BASE / "outputs/dark_lattice_bands.csv", BASE / "outputs/dark_lattice_controls.csv"]},
        "sources": {
            "orbital_elements": "https://ssd.jpl.nasa.gov/planets/approx_pos.html",
            "nominal_solar_GM": "https://arxiv.org/abs/1605.09788",
            "PPN_perihelion": "https://link.springer.com/article/10.12942/lrr-2014-4",
            "MESSENGER_context": "https://www.nature.com/articles/s41467-017-02558-1",
            "quadrupole_orbit_averaging": "https://arxiv.org/abs/1010.1349",
            "solar_Q2_convention": "https://arxiv.org/abs/2602.17884v2"},
        "limitations": [
            "The tiny calculated effects do not reproduce Mercury's 43 arcsec/century relativistic excess",
            "Adding GR's precession is an additional hypothesis, not a derivation or relativistic completion of LQ8",
            "Only the frozen high-field radial tail and the leading existing Q2 are propagated, not higher multipoles or a new full planetary solve",
            "Q2 comes from inner harmonic extrapolation of the existing solar calculation; Mercury was not directly resolved on that grid",
            "Orientation envelopes are geometric bounds at fixed Q2, not statistical confidence intervals",
            "Only five external-field magnitudes are sampled; the amplitudes were selected earlier using Cassini/Gaia",
            "No joint adjustment of planetary masses, solar oblateness, initial conditions, or relativistic parameters",
            "Approximate J2000 elements and nominal GM support scale estimates, not precision ephemeris predictions"]}
    OUT.mkdir(parents=True, exist_ok=True)
    (OUT / "summary.json").write_text(json.dumps(report, indent=2, allow_nan=False) + "\n")
    print(json.dumps({key: report[key] for key in ["GR_reference_arcsec_century", "radial_tail",
        "quadrupole_cases", "largest_quadrupole_envelope_arcsec_century", "controls"]}, indent=2))


if __name__ == "__main__":
    main()
