#!/usr/bin/env python3
"""Fact-Checked Science, /white-dwarf-double-disc/: our numbers for WD J1328-1450.

Recomputes everything the page marks "Our analysis" from public inputs:

1. ESO's 11 public, pipeline-reduced X-shooter VIS spectra of WDJ132800.85-145019.95 (ESO Phase 3,
   "XSHOOTER ECHELLE pipeline-reduced spectroscopic data", https://doi.org/10.18727/archive/71; programmes
   105.20G8.001, 107.22RX.001 and 111.24T9.001). They are fetched from ESO's data portal, which serves
   public files without a login, and checked against the SHA-256 values we saw on 2026-10-31.
   From them: the four peaks of the Ca II lines, their separations, the Keplerian radii those imply,
   the 2021-against-2023 comparison and the Mg I 8806 test.
2. Arithmetic on numbers printed in Gentile Fusillo et al. (2026), Nature Astronomy,
   https://doi.org/10.1038/s41550-026-03004-3 (CC BY 4.0): ring edges in km, the strengthless (fluid)
   Roche densities, and the tensile strength its size limits imply under eq. S16 of Manser et al. (2019).

Usage:
    python3 j1328_numbers.py            downloads about 15 MB into ./eso-xshooter-j1328/ and runs (a few seconds)
    python3 j1328_numbers.py <folder>   uses copies already in <folder> (file names as ESO serves them,
                                        with or without ':' replaced by '_')
Writes j1328-profiles.csv and j1328-numbers.json next to this script. Standard library only.

Based on data obtained from the ESO Science Archive Facility with DOI: https://doi.org/10.18727/archive/71.
ESO's processed data may be redistributed under CC BY 4.0 with credit; this script hosts none, it fetches them.
"""
import hashlib
import json
import math
import os
import random
import statistics
import struct
import sys
import urllib.request

HERE = os.path.dirname(os.path.abspath(__file__))
PORTAL = "https://dataportal.eso.org/dataPortal/file/"
FILES = {  # ESO product ID: SHA-256 seen on 2026-10-31
    "ADP.2021-06-01T10:29:43.801": "545088422e0b18cd2d1683b74db1d85af73b605408e9c057a9eecc383f7bc04c",
    "ADP.2021-06-01T10:29:43.589": "821953cdbff711a781da83e7dfe3e35340d50206280b64fcfc9abe321d90b54c",
    "ADP.2021-06-01T10:29:43.624": "9d26ba5d4073b484beda65c89fe1683477a9dc38de604fa528626ad8558ffc08",
    "ADP.2021-07-06T09:48:03.697": "cd25d9303eb9bbe3f7e1570cf84ab851d9d42a30f4e4b44c0f56272fb2a4bdf9",
    "ADP.2021-07-06T09:48:03.745": "477c94f0334483232ce5bd50c54f18b1d0ec7ba881e5a14e94df561d7879c085",
    "ADP.2021-07-16T11:53:45.234": "939a1102c468a5b2b1fa72ea2972175d6074b0db2dddcb442cffa8c496a2c9bd",
    "ADP.2021-07-16T11:53:45.265": "401d00ac14292cd289f3a2bd98e369e6d331ed7461372a8245ee04c0297bd534",
    "ADP.2021-07-16T11:53:45.260": "cc9d9beea15e23eff035444d16f6a7cc18439be9eb1e3ef4fea89da1823ecc1e",
    "ADP.2023-06-01T07:51:00.515": "87463186b275e4fdcfa3e3fc0fa4ee4810142aaaca3cc97c6efac9506bd39d1d",
    "ADP.2023-06-01T07:51:00.520": "5666a28a669c8bb09fb243de088ec07ac2e36d2fb7b58676d49a5f4929afbb96",
    "ADP.2023-06-01T07:51:00.553": "4f473eaa8fa358e2b72a2ceb0ab715a605e98e57832edb7299b645f53079776d",
}

C_KMS = 299792.458
G = 6.67430e-11
MSUN = 1.98847e30
RSUN_M = 6.957e8
KB = 1.380649e-23
MH = 1.6735575e-27
LINES_AIR = {"CaII8498": 8498.02, "CaII8542": 8542.09, "CaII8662": 8662.14, "MgI8806": 8806.757}  # Angstrom, air
M_WD, INC = 0.47, 77.0          # the paper's fixed mass (Msun) and best-fit inclination (deg)
CONT_WINDOWS = [(8440, 8470), (8570, 8635), (8695, 8770)]
MG_CONT_WINDOWS = [(8745, 8782), (8850, 8870)]  # clear of Mg I itself and of Mg II 8824/8835


# ---------- FITS (one-row binary table, as ESO's Phase 3 spectra are) ----------

def _header(f):
    cards = []
    while True:
        block = f.read(2880)
        if len(block) < 2880:
            raise ValueError("truncated FITS header")
        for i in range(0, 2880, 80):
            c = block[i:i + 80].decode("ascii", "replace")
            cards.append(c)
            if c.startswith("END") and not c[3:].strip():
                return cards


def _cards(cards):
    h = {}
    for c in cards:
        if c.startswith("HIERARCH"):
            k, _, v = c[9:].partition("=")
        elif c[8:10] == "= ":
            k, v = c[:8], c[10:]
        else:
            continue
        k, v = k.strip(), v.strip()
        if v.startswith("'"):
            v = v[1:v.find("'", 1)].strip()
        else:
            v = v.split("/")[0].strip()
            for cast in (int, float):
                try:
                    v = cast(v)
                    break
                except ValueError:
                    pass
        h[k] = v
    return h


def read_spectrum(path):
    with open(path, "rb") as f:
        h0 = _cards(_header(f))
        h1 = _cards(_header(f))
        raw = f.read(h1["NAXIS1"] * h1["NAXIS2"])
    fmt = {"D": ("d", 8), "E": ("f", 4), "J": ("i", 4), "I": ("h", 2), "K": ("q", 8)}
    cols, off = {}, 0
    for j in range(1, h1["TFIELDS"] + 1):
        form = h1["TFORM%d" % j]
        n = int("".join(ch for ch in form if ch.isdigit()) or 1)
        code, size = fmt[[ch for ch in form if ch.isalpha()][0]]
        cols[h1["TTYPE%d" % j]] = struct.unpack(">%d%s" % (n, code), raw[off:off + n * size])
        off += n * size
    bary = float(h0.get("ESO QC VRAD BARYCOR", 0.0))
    return {"date": str(h0["DATE-OBS"])[:10], "exptime": h0["EXPTIME"], "bary": bary,
            "wave": [10 * w * (1 + bary / C_KMS) for w in cols["WAVE"]],  # nm -> Angstrom, barycentric
            "flux": cols["FLUX"], "err": cols["ERR"], "qual": cols["QUAL"]}


# ---------- fetching ----------

def fetch(folder):
    os.makedirs(folder, exist_ok=True)
    paths, notes = [], []
    for pid, sha in FILES.items():
        cands = [os.path.join(folder, pid + ".fits"), os.path.join(folder, pid.replace(":", "_") + ".fits")]
        path = next((p for p in cands if os.path.isfile(p)), None)
        if path is None:
            path = cands[1]
            req = urllib.request.Request(PORTAL + pid, headers={"User-Agent": "factchecked.science j1328_numbers.py"})
            with urllib.request.urlopen(req, timeout=120) as r, open(path, "wb") as out:
                out.write(r.read())
        got = hashlib.sha256(open(path, "rb").read()).hexdigest()
        if got != sha:
            notes.append(f"{pid}: SHA-256 {got} differs from the {sha} we saw on 2026-10-31")
        paths.append(path)
    return paths, notes


# ---------- spectra ----------

def resample(s, grid):
    w, f, e, q = s["wave"], s["flux"], s["err"], s["qual"]
    out_f, out_e, j = [], [], 0
    for x in grid:
        while j < len(w) - 2 and w[j + 1] < x:
            j += 1
        if not (w[j] <= x <= w[j + 1]) or q[j] or q[j + 1] or e[j] <= 0 or e[j + 1] <= 0:
            out_f.append(None)
            out_e.append(None)
            continue
        t = (x - w[j]) / (w[j + 1] - w[j])
        out_f.append(f[j] * (1 - t) + f[j + 1] * t)
        out_e.append(math.hypot(e[j] * (1 - t), e[j + 1] * t))
    return out_f, out_e


def coadd(specs, grid):
    num, den = [0.0] * len(grid), [0.0] * len(grid)
    for s in specs:
        f, e = resample(s, grid)
        for i, (fi, ei) in enumerate(zip(f, e)):
            if fi is not None:
                num[i] += fi / ei ** 2
                den[i] += 1 / ei ** 2
    return ([n / d if d else None for n, d in zip(num, den)], [1 / math.sqrt(d) if d else None for d in den])


def continuum(grid, flux, windows):
    xs = [x for x, y in zip(grid, flux) if y is not None and any(a <= x <= b for a, b in windows)]
    ys = [y for x, y in zip(grid, flux) if y is not None and any(a <= x <= b for a, b in windows)]
    mx, my = statistics.fmean(xs), statistics.fmean(ys)
    slope = sum((x - mx) * (y - my) for x, y in zip(xs, ys)) / sum((x - mx) ** 2 for x in xs)
    a = my - slope * mx
    rms = statistics.pstdev([y - (a + slope * x) for x, y in zip(xs, ys)])
    return a, slope, rms


def smooth(vals, k=5):
    h, out = k // 2, []
    for i in range(len(vals)):
        seg = [v for v in vals[max(0, i - h):i + h + 1] if v is not None]
        out.append(sum(seg) / len(seg) if seg else None)
    return out


def peak(vel, prof, lo, hi):
    best = None
    for v, y in zip(vel, prof):
        if y is not None and lo <= v <= hi and (best is None or y > best[1]):
            best = (v, y)
    return best


WINDOWS = {"outer_blue": (-330, -110), "outer_red": (110, 330), "inner_blue": (-720, -380), "inner_red": (380, 720)}


def peaks_of(grid, flux, err, cont, line, rng, nmc):
    a, b = cont[0], cont[1]
    vel = [C_KMS * (x / LINES_AIR[line] - 1) for x in grid]
    norm = [y / (a + b * x) if y is not None else None for x, y in zip(grid, flux)]
    nerr = [e / (a + b * x) if e is not None else None for x, e in zip(grid, err)]
    sm = smooth(norm)
    res = {k: peak(vel, sm, lo, hi) for k, (lo, hi) in WINDOWS.items()}
    mc = {k: [] for k in WINDOWS}
    for _ in range(nmc):
        pert = smooth([y + rng.gauss(0, s) if y is not None and s is not None else None for y, s in zip(norm, nerr)])
        for k, (lo, hi) in WINDOWS.items():
            mc[k].append(peak(vel, pert, lo, hi)[0])
    sd = {k: statistics.pstdev(v) for k, v in mc.items()}
    dip = statistics.fmean([y for v, y in zip(vel, sm) if y is not None and 300 <= abs(v) <= 380])
    return {
        "inner_kms": [round(res["inner_blue"][0]), round(res["inner_red"][0])],
        "outer_kms": [round(res["outer_blue"][0]), round(res["outer_red"][0])],
        "heights_inner": [round(res["inner_blue"][1], 3), round(res["inner_red"][1], 3)],
        "heights_outer": [round(res["outer_blue"][1], 3), round(res["outer_red"][1], 3)],
        "dip_300_380": round(dip, 3),
        "sep_inner_kms": round(res["inner_red"][0] - res["inner_blue"][0]),
        "sep_inner_sd": round(math.hypot(sd["inner_blue"], sd["inner_red"])),
        "sep_outer_kms": round(res["outer_red"][0] - res["outer_blue"][0]),
        "sep_outer_sd": round(math.hypot(sd["outer_blue"], sd["outer_red"])),
        "centre_outer_kms": round((res["outer_red"][0] + res["outer_blue"][0]) / 2),
    }, vel, norm, nerr


def kepler_r(sep_kms, m=M_WD, inc=INC):
    v = sep_kms / 2 * 1e3
    return G * m * MSUN * math.sin(math.radians(inc)) ** 2 / v ** 2 / RSUN_M


def main():
    folder = sys.argv[1] if len(sys.argv) > 1 else os.path.join(os.getcwd(), "eso-xshooter-j1328")
    paths, notes = fetch(folder)
    specs = sorted((read_spectrum(p) for p in paths), key=lambda s: s["date"])
    rng = random.Random(20261031)
    grid = [8400.0 + 0.2 * i for i in range(2501)]  # 8400-8900 A, the pipeline's 0.2 A step
    groups = {"2021": [s for s in specs if s["date"] < "2022"], "2023": [s for s in specs if s["date"] > "2022"], "all": specs}
    out = {"inputs": {"files": len(specs), "dates": sorted({s["date"] for s in specs}),
                      "per_date": {d: sum(1 for s in specs if s["date"] == d) for d in sorted({s["date"] for s in specs})},
                      "checksum_notes": notes},
           "acknowledgement": "Based on data obtained from the ESO Science Archive Facility with DOI: https://doi.org/10.18727/archive/71"}
    coadds, rows = {}, []
    for g, ss in groups.items():
        flux, err = coadd(ss, grid)
        cont = continuum(grid, flux, CONT_WINDOWS)
        coadds[g] = (flux, err, cont)
        out.setdefault("peaks", {})[g] = {}
        for line in ("CaII8498", "CaII8542", "CaII8662"):
            res, vel, norm, nerr = peaks_of(grid, flux, err, cont, line, rng, 200 if g == "all" else 60)
            out["peaks"][g][line] = res
            if line == "CaII8662" and g != "all":
                rows += [("CaII8662", g, round(v, 1), round(y, 4), round(s, 4)) for v, y, s in zip(vel, norm, nerr)
                         if y is not None and abs(v) <= 900]
    # Keplerian radii at the paper's mass and inclination, and at the ends of their 1-sigma ranges
    out["radii_rsun"] = {}
    for line in ("CaII8542", "CaII8662"):
        p = out["peaks"]["all"][line]
        out["radii_rsun"][line] = {
            "inner_ring_outer_edge": [round(kepler_r(p["sep_inner_kms"], m, i), 3) for m, i in ((M_WD, INC), (0.38, 64.0), (0.56, 86.0))],
            "outer_ring_outer_edge": [round(kepler_r(p["sep_outer_kms"], m, i), 2) for m, i in ((M_WD, INC), (0.38, 64.0), (0.56, 86.0))],
            "note": "values at (M, i) = (0.47, 77), (0.38, 64), (0.56, 86)"}
    # 2021 against 2023: chi-square per pixel of the difference, against the same in line-free windows
    f21, e21, c21 = coadds["2021"]
    f23, e23, c23 = coadds["2023"]

    def chi2(windows):
        tot, n = 0.0, 0
        for x, a1, s1, a2, s2 in zip(grid, f21, e21, f23, e23):
            if None in (a1, s1, a2, s2) or not any(lo <= x <= hi for lo, hi in windows):
                continue
            k1, k2 = c21[0] + c21[1] * x, c23[0] + c23[1] * x
            tot += ((a1 / k1 - a2 / k2) / math.hypot(s1 / k1, s2 / k2)) ** 2
            n += 1
        return tot / n
    ref = chi2(CONT_WINDOWS)
    out["stability"] = {"chi2_per_pixel_continuum": round(ref, 2)}
    for line in ("CaII8542", "CaII8662"):
        lam = LINES_AIR[line]
        val = chi2([(lam * (1 - 800 / C_KMS), lam * (1 + 800 / C_KMS))])
        out["stability"][line] = {"chi2_per_pixel_within_800kms": round(val, 2), "ratio_to_continuum": round(val / ref, 2)}
    # Mg I 8806: flux where the outer ring's peaks would fall, against the paper's own inner-peak level
    flux, err, _ = coadds["all"]
    cm = continuum(grid, flux, MG_CONT_WINDOWS)
    vel = [C_KMS * (x / LINES_AIR["MgI8806"] - 1) for x in grid]
    norm = [y / (cm[0] + cm[1] * x) if y is not None else None for x, y in zip(grid, flux)]
    rms = cm[2] / (cm[0] + cm[1] * LINES_AIR["MgI8806"])
    p8542 = out["peaks"]["all"]["CaII8542"]
    h_in, h_out = statistics.fmean(p8542["heights_inner"]), statistics.fmean(p8542["heights_outer"])
    paper_inner = 1.245  # the paper's inner Mg I peak, normalised (1.94e-17 over its continuum), read off its Extended Data Fig. 2b by eye
    predicted = {"paper_method": round(1.22 * paper_inner, 3), "line_scaling": round(1 + (h_out - 1) / (h_in - 1) * (paper_inner - 1), 3)}
    mg = {"continuum_rms_per_pixel": round(rms, 3), "predicted_outer_peak": predicted,
          "ca8542_outer_over_inner_peak_height": {"of_the_totals": round(h_out / h_in, 2), "above_continuum": round((h_out - 1) / (h_in - 1), 2)},
          "windows": {}}
    for lo, hi in ((-235, -180), (180, 235)):
        vals = [y for v, y in zip(vel, norm) if y is not None and lo <= v <= hi]
        m, se = statistics.fmean(vals), rms / math.sqrt(len(vals))
        mg["windows"][f"{lo}..{hi}"] = {"mean": round(m, 3), "n": len(vals), "se": round(se, 3),
                                        "shortfall_se": [round((predicted["paper_method"] - m) / se, 1), round((predicted["line_scaling"] - m) / se, 1)]}
    out["mg8806"] = mg
    rows += [("MgI8806", "all", round(v, 1), round(y, 4), None) for v, y in zip(vel, norm) if y is not None and abs(v) <= 900]
    # Arithmetic on the paper's own numbers
    m = M_WD * MSUN
    out["arithmetic"] = {
        "ring_edges_km": {str(r): round(r * RSUN_M / 1e3, -2) for r in (0.14, 0.34, 0.54, 1.93)},
        "two_vk_sini_kms": {str(r): round(2 * math.sqrt(G * m / (r * RSUN_M)) / 1e3 * math.sin(math.radians(INC))) for r in (0.34, 1.93)},
        "fluid_roche_density_gcc": {str(r): round(2.44 ** 3 * 3 * m * 1e3 / (4 * math.pi * (r * RSUN_M * 100) ** 3), 1) for r in (0.14, 0.34, 0.54, 1.93)},
        "implied_strength_mpa_eqS16": {f"{s} km at {a} Rsun, rho {rho}": round((s * 1e3) ** 2 * G * m * rho * 1e3 / (2 * (a * RSUN_M) ** 3) / 1e6)
                                       for s, a, rho in ((45, 0.14, 7.8), (172, 0.34, 7.8), (346, 0.54, 7.8), (2340, 1.93, 7.8), (3500, 1.93, 3.5))},
        "si_case_eqS16_km": round(math.sqrt(2 * 5e8 * (0.62 * RSUN_M) ** 3 / (G * m * 7900)) / 1e3),
        # eq. S16's s is a "size" of about 1.65 radii (Manser et al. 2019, supplementary section 3): if the paper's figures
        # are true diameters, s = 0.825 D and the implied strengths fall by 0.825 squared
        "implied_strength_mpa_if_true_diameters": {f"{s} km at {a} Rsun, rho {rho}": round((0.825 * s * 1e3) ** 2 * G * m * rho * 1e3 / (2 * (a * RSUN_M) ** 3) / 1e6)
                                                   for s, a, rho in ((45, 0.14, 7.8), (172, 0.34, 7.8), (346, 0.54, 7.8), (2340, 1.93, 7.8), (3500, 1.93, 3.5))},
        # scaled to the supplement's own case (440 km at 0.62 Rsun, 7.9 g/cm3, 500 MPa), which needs no reading of s
        "implied_strength_mpa_scaled_to_si_case": {f"{s} km at {a} Rsun, rho {rho}": round(500 * (s / 440) ** 2 * (0.62 / a) ** 3 * (rho / 7.9))
                                                   for s, a, rho in ((45, 0.14, 7.8), (172, 0.34, 7.8), (346, 0.54, 7.8), (2340, 1.93, 7.8), (3500, 1.93, 3.5))},
        # the supplement's iron body of 440 km at 0.62 Rsun: escape speed and Hill radius (the Hill radius over the body's radius
        # does not depend on its size)
        "si_case_escape_speed_ms": {"440 km as diameter": round(math.sqrt(8 / 3 * math.pi * G * 7900) * 220e3),
                                    "440 km as radius": round(math.sqrt(8 / 3 * math.pi * G * 7900) * 440e3)},
        "si_case_hill_radius_over_body_radius": round(0.62 * RSUN_M * (4 * math.pi * 7900 / (9 * m)) ** (1 / 3), 3),
        "max_size_km_at_10MPa_iron": {str(a): round(math.sqrt(2 * 1e7 * (a * RSUN_M) ** 3 / (G * m * 7800)) / 1e3) for a in (0.14, 0.34, 0.54, 1.93)},
        "scale_height_mu40": {f"{t} K at {r} Rsun": round(math.sqrt(KB * t / (40 * MH)) / math.sqrt(G * m / (r * RSUN_M)), 4)
                              for t, r in ((3000, 0.14), (3000, 0.34), (1300, 0.54), (1300, 1.93))},
        "distance_pc_from_parallax": [round(1000 / 2.014), round(1000 / (2.014 + 0.292)), round(1000 / (2.014 - 0.292))],
        "abundance_ratios_from_ED_table_1": {"log_Mg_Ca": round(-7.00 + 7.50, 2), "log_Si_Ca": round(-6.30 + 7.50, 2), "Mg_Si": round(10 ** (-7.00 + 6.30), 2)},
    }
    with open(os.path.join(HERE, "j1328-profiles.csv"), "w", encoding="utf-8", newline="\n") as f:
        f.write("# Co-added, continuum-normalised ESO X-shooter VIS spectra of WD J1328-1450 (our analysis).\n")
        f.write("# Based on data obtained from the ESO Science Archive Facility with DOI: https://doi.org/10.18727/archive/71 (CC BY 4.0).\n")
        f.write("# Velocity in km/s from the line's air rest wavelength, barycentric; 2021 = 8 spectra, 2023 = 3, all = 11.\n")
        f.write("line,epochs,velocity_kms,flux_normalised,error_normalised\n")
        for r in rows:
            f.write(",".join("" if x is None else str(x) for x in r) + "\n")
    with open(os.path.join(HERE, "j1328-numbers.json"), "w", encoding="utf-8", newline="\n") as f:
        json.dump(out, f, indent=1, ensure_ascii=False)
        f.write("\n")
    print(json.dumps(out, indent=1, ensure_ascii=False))


if __name__ == "__main__":
    main()