"""Per-window quality record for the committed sector 56 light curve.

Run from the repository root:

    python research/tess-timing/scripts/real_window_qc.py
    python research/tess-timing/scripts/real_window_qc.py --check

The command reads two committed files only, each of which is verified by
SHA-256 before it is used:

    research/tess-timing/data/real-run/lightcurve.fits
    research/tess-timing/data/real-run/tpf.fits

Both primary headers must report ``TICID = 420814525`` and ``SECTOR = 56``.
The eight window centres are read from the committed prediction record

    research/tess-timing/data/real-windows-run/windows.json

and are not re-derived here. Each window is the fixed interval

    centre_btjd - 0.5 d <= BTJD <= centre_btjd + 0.5 d

where the centre is the predicted mid-transit time already recorded in that
file. The 0.5 d half-width is a diagnostic choice, not a measured transit
duration, and the committed prediction record is not changed.

For each window the record reports cadence counts before and after the 175
quality mask, finite SAP and PDCSAP counts after the mask, gap count and the
longest gap at the existing 180 s threshold, SAP_BKG median and 5th to 95th
range, and MOM_CENTR1/MOM_CENTR2 median and peak-to-peak.

The committed outputs are:

    research/tess-timing/data/real-window-qc/windows.csv
    research/tess-timing/data/real-window-qc/windows.json
    docs/TESS-REAL-WINDOW-QC.md

This is a diagnostic record only. No contamination ratio or nearby-source
count is computed because no source catalogue is committed. No detrending,
transit fit, fold, or timing measurement is performed. No network request is
made and the input files are not modified.
"""

from __future__ import annotations

import argparse
import csv
import itertools
import json
import math
import shutil
import sys
from pathlib import Path
from typing import Any

PACKAGE_ROOT = Path(__file__).resolve().parents[1]
REPO_ROOT = PACKAGE_ROOT.parents[1]
SRC = PACKAGE_ROOT / "src"
if str(SRC) not in sys.path:
    sys.path.insert(0, str(SRC))

# isort: off
import numpy as np
from astropy.io import fits
from tess_timing.lightcurve import portable_path, sha256_file
# isort: on

TASK = "jdp-real-window-qc"
SCHEMA_VERSION = 1
TARGET_NAME = "HD 209458 b"
EXPECTED_TIC_ID = 420814525
EXPECTED_SECTOR = 56
EXPECTED_MASK_BITS = 175
EXPECTED_WINDOW_COUNT = 8
EXPECTED_CADENCE_SECONDS = 120.0
GAP_FACTOR = 1.5
GAP_THRESHOLD_SECONDS = 180.0
SECONDS_PER_DAY = 86400.0
WINDOW_HALF_WIDTH_DAYS = 0.5
HALF_WIDTH = WINDOW_HALF_WIDTH_DAYS

EXPECTED_LIGHTCURVE_SHA256 = (
    "1b76b4a4b73e24954fa6e29a7e28a174685113ecdbad1bf10ef0e7b766ce99e9"
)
EXPECTED_TPF_SHA256 = "2b5bd6ee4747cfe0c5405c869b871064fd74e4b5f0521c4e6be9fdc2b5c96030"

DEFAULT_LIGHTCURVE = PACKAGE_ROOT / "data" / "real-run" / "lightcurve.fits"
DEFAULT_TPF = PACKAGE_ROOT / "data" / "real-run" / "tpf.fits"
DEFAULT_WINDOWS = PACKAGE_ROOT / "data" / "real-windows-run" / "windows.json"
DEFAULT_CSV = PACKAGE_ROOT / "data" / "real-window-qc" / "windows.csv"
DEFAULT_JSON = PACKAGE_ROOT / "data" / "real-window-qc" / "windows.json"
DEFAULT_MARKDOWN = REPO_ROOT / "docs" / "TESS-REAL-WINDOW-QC.md"
CHECK_ROOT = REPO_ROOT / ".scratch" / "real-window-qc-check"

LIGHTCURVE_COLUMNS = (
    "TIME",
    "QUALITY",
    "SAP_FLUX",
    "PDCSAP_FLUX",
    "SAP_BKG",
    "MOM_CENTR1",
    "MOM_CENTR2",
)

CSV_COLUMNS = (
    "cycle",
    "btjd",
    "bjd",
    "sigma_days",
    "window_start_btjd",
    "window_end_btjd",
    "cadences_before_mask",
    "cadences_after_mask",
    "removed_by_quality_mask",
    "finite_sap_after_mask",
    "finite_pdcsap_after_mask",
    "gap_count",
    "longest_gap_seconds",
    "sap_bkg_cadences",
    "sap_bkg_median",
    "sap_bkg_percentile_5",
    "sap_bkg_percentile_95",
    "sap_bkg_range_5_95",
    "mom_centr1_cadences",
    "mom_centr1_median",
    "mom_centr1_peak_to_peak",
    "mom_centr2_cadences",
    "mom_centr2_median",
    "mom_centr2_peak_to_peak",
)

COUNT_COLUMNS = (
    "cadences_before_mask",
    "cadences_after_mask",
    "removed_by_quality_mask",
    "finite_sap_after_mask",
    "finite_pdcsap_after_mask",
    "gap_count",
    "sap_bkg_cadences",
    "mom_centr1_cadences",
    "mom_centr2_cadences",
)

LIMITS = [
    (
        "Diagnostic only. This record describes cadence availability and "
        "per-window column statistics for the committed sector 56 light "
        "curve. It is not a detection, a transit measurement, or a research "
        "result."
    ),
    (
        "No contamination ratio and no nearby-source count are computed "
        "because no source catalogue is committed in this repository."
    ),
    (
        "No detrending, transit fit, folding, or transit timing measurement "
        "is performed."
    ),
    (
        "Window centres are the eight predicted mid-transit times read from "
        "the committed windows record. They are not measured transit times."
    ),
    (
        "SAP_BKG and MOM_CENTR statistics are per-window summaries after the "
        "175 quality mask; they are not calibrated astrophysical quantities."
    ),
    (
        "Gap counts use finite TIME values inside each window before the "
        "quality mask, with a fixed 180 s threshold (1.5 times the 120 s "
        "cadence)."
    ),
    ("No target ranking, detection, or novelty claim is made."),
]


class WindowQCError(RuntimeError):
    """An input file, the windows record, or the generated record is unusable."""


def _display(path: Path) -> str:
    return portable_path(path)


def _number(value: object, context: str) -> float:
    if isinstance(value, bool) or not isinstance(value, (int, float)):
        raise WindowQCError(f"{context}: expected a number")
    number = float(value)
    if not math.isfinite(number):
        raise WindowQCError(f"{context}: expected a finite number")
    return number


def _integer(value: object, context: str) -> int:
    if isinstance(value, bool) or not isinstance(value, int):
        raise WindowQCError(f"{context}: expected an integer")
    return value


def _text(value: object, context: str) -> str:
    if not isinstance(value, str) or not value.strip():
        raise WindowQCError(f"{context}: expected a non-empty string")
    return value


def _header_int(header: fits.Header, key: str, path: Path) -> int:
    value = header.get(key)
    if value is None:
        raise WindowQCError(f"{path}: primary header has no {key} value")
    try:
        return int(value)
    except (TypeError, ValueError) as error:
        raise WindowQCError(
            f"{path}: primary header {key} value {value!r} is not an integer"
        ) from error


def _preflight_fits(path: Path, expected_sha256: str, label: str) -> dict[str, Any]:
    """Verify a committed FITS input before reading any table data."""

    path = Path(path)
    if not path.exists():
        raise WindowQCError(f"{label} input does not exist: {path}")
    digest = sha256_file(path)
    if digest != expected_sha256:
        raise WindowQCError(
            f"{label} input SHA-256 is {digest}; expected {expected_sha256}. "
            "The file is not the committed sector 56 product."
        )
    try:
        with fits.open(path, memmap=False) as hdul:
            header = hdul[0].header
            tic_id = _header_int(header, "TICID", path)
            sector = _header_int(header, "SECTOR", path)
            object_name = header.get("OBJECT")
            synthetic = bool(header.get("SYNTH", False))
            rows = None
            for hdu in hdul:
                if isinstance(hdu, fits.BinTableHDU) and hdu.data is not None:
                    rows = len(hdu.data)
                    break
    except OSError as error:
        raise WindowQCError(
            f"could not open {label} FITS file {path}: {error}"
        ) from error

    if tic_id != EXPECTED_TIC_ID:
        raise WindowQCError(
            f"{label} primary header TICID is {tic_id}; expected {EXPECTED_TIC_ID}"
        )
    if sector != EXPECTED_SECTOR:
        raise WindowQCError(
            f"{label} primary header SECTOR is {sector}; expected {EXPECTED_SECTOR}"
        )
    if synthetic:
        raise WindowQCError(f"{label} header marks the input as synthetic")

    return {
        "path": _display(path),
        "sha256": digest,
        "expected_sha256": expected_sha256,
        "size_bytes": int(path.stat().st_size),
        "tic_id": tic_id,
        "sector": sector,
        "object": str(object_name) if object_name is not None else None,
        "table_rows": rows,
        "verified": True,
    }


def _fits_column(table: fits.BinTableHDU, name: str) -> np.ndarray:
    if table.columns is None or name not in table.columns.names:
        raise WindowQCError(f"light-curve table has no {name} column")
    values = table.data[name]
    if np.ma.isMaskedArray(values):
        values = values.filled(np.nan)
    return np.asarray(values, dtype=float)


def _read_lightcurve(path: Path) -> dict[str, np.ndarray]:
    """Read the requested light-curve columns without normalizing them."""

    try:
        with fits.open(Path(path), memmap=False) as hdul:
            table = None
            for hdu in hdul:
                if (
                    isinstance(hdu, fits.BinTableHDU)
                    and hdu.columns is not None
                    and "TIME" in set(hdu.columns.names)
                ):
                    table = hdu
                    break
            if table is None:
                raise WindowQCError("the light curve has no binary table with TIME")
            columns = {name: _fits_column(table, name) for name in LIGHTCURVE_COLUMNS}
    except OSError as error:
        raise WindowQCError(f"could not read light curve {path}: {error}") from error

    quality_float = columns["QUALITY"]
    quality_finite = np.isfinite(quality_float)
    rounded = np.where(quality_finite, np.rint(quality_float), 0.0)
    if np.any(quality_finite & (quality_float != rounded)):
        raise WindowQCError("QUALITY column contains non-integer values")
    columns["QUALITY"] = np.where(quality_finite, rounded, 0.0).astype(np.int64)
    columns["_quality_finite"] = quality_finite
    return columns


def _read_tpf_time(path: Path) -> np.ndarray:
    try:
        with fits.open(Path(path), memmap=False) as hdul:
            for hdu in hdul:
                if (
                    isinstance(hdu, fits.BinTableHDU)
                    and hdu.columns is not None
                    and "TIME" in set(hdu.columns.names)
                ):
                    values = hdu.data["TIME"]
                    if np.ma.isMaskedArray(values):
                        values = values.filled(np.nan)
                    return np.asarray(values, dtype=float)
    except OSError as error:
        raise WindowQCError(
            f"could not read target pixel file {path}: {error}"
        ) from error
    raise WindowQCError("the target pixel file has no binary table with TIME")


def _read_windows(path: Path) -> dict[str, Any]:
    """Read and validate the committed prediction record without re-deriving it."""

    path = Path(path)
    try:
        record = json.loads(path.read_text(encoding="utf-8"))
    except (OSError, json.JSONDecodeError) as error:
        raise WindowQCError(f"could not read windows record {path}: {error}") from error
    if not isinstance(record, dict):
        raise WindowQCError("the windows record must be a JSON object")

    target = record.get("target")
    if not isinstance(target, dict):
        raise WindowQCError("the windows record has no target object")
    tic_id = _integer(target.get("tic_id"), "windows target tic_id")
    sector = _integer(target.get("sector"), "windows target sector")
    if tic_id != EXPECTED_TIC_ID:
        raise WindowQCError(
            f"windows record target TIC is {tic_id}; expected {EXPECTED_TIC_ID}"
        )
    if sector != EXPECTED_SECTOR:
        raise WindowQCError(
            f"windows record target sector is {sector}; expected {EXPECTED_SECTOR}"
        )

    predictions = record.get("predictions")
    if not isinstance(predictions, list):
        raise WindowQCError("the windows record predictions must be a list")
    if len(predictions) != EXPECTED_WINDOW_COUNT:
        raise WindowQCError(
            f"the windows record has {len(predictions)} predictions; expected "
            f"{EXPECTED_WINDOW_COUNT}"
        )

    windows: list[dict[str, Any]] = []
    for index, row in enumerate(predictions):
        context = f"windows prediction {index}"
        if not isinstance(row, dict):
            raise WindowQCError(f"{context}: expected an object")
        cycle = _integer(row.get("cycle"), f"{context} cycle")
        btjd = _number(row.get("btjd"), f"{context} btjd")
        bjd = _number(row.get("bjd"), f"{context} bjd")
        sigma_days = _number(row.get("sigma_days"), f"{context} sigma_days")
        if sigma_days <= 0.0:
            raise WindowQCError(f"{context} sigma_days: expected a positive number")
        if row.get("inside_sector") is not True:
            raise WindowQCError(f"{context}: inside_sector must be true")
        windows.append(
            {
                "cycle": cycle,
                "btjd": btjd,
                "bjd": bjd,
                "sigma_days": sigma_days,
            }
        )

    cycles = [row["cycle"] for row in windows]
    if len(set(cycles)) != len(cycles):
        raise WindowQCError("the windows record contains duplicate cycles")
    centres = [row["btjd"] for row in windows]
    if any(right <= left for left, right in itertools.pairwise(centres)):
        raise WindowQCError("window centres are not strictly increasing")
    if any(
        right - left <= 2.0 * HALF_WIDTH for left, right in itertools.pairwise(centres)
    ):
        raise WindowQCError(
            "adjacent windows overlap; per-window totals would double count"
        )

    return {
        "path": _display(path),
        "sha256": sha256_file(path),
        "task": record.get("task"),
        "prediction_count": len(windows),
        "target": {
            "name": target.get("name"),
            "tic_id": tic_id,
            "sector": sector,
        },
        "windows": windows,
    }


def _window_bounds(centre: float) -> tuple[float, float]:
    return centre - HALF_WIDTH, centre + HALF_WIDTH


def _gap_summary(time: np.ndarray, mask: np.ndarray) -> dict[str, Any]:
    """Count gaps at the fixed 180 s threshold inside one window."""

    selected = np.asarray(time, dtype=float)[mask]
    times = np.sort(selected[np.isfinite(selected)])
    if times.size < 2:
        return {"gap_count": 0, "longest_gap_seconds": 0.0}
    differences = np.diff(times)
    threshold_days = GAP_THRESHOLD_SECONDS / SECONDS_PER_DAY
    gaps = differences[differences > threshold_days]
    longest = float(gaps.max() * SECONDS_PER_DAY) if gaps.size else 0.0
    return {
        "gap_count": int(gaps.size),
        "longest_gap_seconds": round(longest, 6),
    }


def _column_stats(values: np.ndarray, mask: np.ndarray) -> dict[str, Any]:
    selected = np.asarray(values, dtype=float)[mask]
    finite = selected[np.isfinite(selected)]
    if finite.size == 0:
        return {
            "cadences": 0,
            "median": None,
            "percentile_5": None,
            "percentile_95": None,
            "range_5_95": None,
        }
    p5, median, p95 = np.percentile(finite, [5.0, 50.0, 95.0])
    return {
        "cadences": int(finite.size),
        "median": float(median),
        "percentile_5": float(p5),
        "percentile_95": float(p95),
        "range_5_95": float(p95 - p5),
    }


def _centroid_stats(values: np.ndarray, mask: np.ndarray) -> dict[str, Any]:
    selected = np.asarray(values, dtype=float)[mask]
    finite = selected[np.isfinite(selected)]
    if finite.size == 0:
        return {"cadences": 0, "median": None, "peak_to_peak": None}
    return {
        "cadences": int(finite.size),
        "median": float(np.median(finite)),
        "peak_to_peak": float(finite.max() - finite.min()),
    }


def _window_record(
    window: dict[str, Any],
    columns: dict[str, np.ndarray],
    quality_pass: np.ndarray,
    finite_time: np.ndarray,
) -> dict[str, Any]:
    start, end = _window_bounds(window["btjd"])
    in_window = finite_time & (columns["TIME"] >= start) & (columns["TIME"] <= end)
    after_mask = in_window & quality_pass
    before_count = int(np.count_nonzero(in_window))
    after_count = int(np.count_nonzero(after_mask))
    removed = before_count - after_count

    gaps = _gap_summary(columns["TIME"], in_window)
    background = _column_stats(columns["SAP_BKG"], after_mask)
    centr1 = _centroid_stats(columns["MOM_CENTR1"], after_mask)
    centr2 = _centroid_stats(columns["MOM_CENTR2"], after_mask)

    return {
        "cycle": window["cycle"],
        "btjd": window["btjd"],
        "bjd": window["bjd"],
        "sigma_days": window["sigma_days"],
        "window_start_btjd": start,
        "window_end_btjd": end,
        "cadences_before_mask": before_count,
        "cadences_after_mask": after_count,
        "removed_by_quality_mask": removed,
        "finite_sap_after_mask": int(
            np.count_nonzero(after_mask & np.isfinite(columns["SAP_FLUX"]))
        ),
        "finite_pdcsap_after_mask": int(
            np.count_nonzero(after_mask & np.isfinite(columns["PDCSAP_FLUX"]))
        ),
        "gap_count": gaps["gap_count"],
        "longest_gap_seconds": gaps["longest_gap_seconds"],
        "sap_bkg_cadences": background["cadences"],
        "sap_bkg_median": background["median"],
        "sap_bkg_percentile_5": background["percentile_5"],
        "sap_bkg_percentile_95": background["percentile_95"],
        "sap_bkg_range_5_95": background["range_5_95"],
        "mom_centr1_cadences": centr1["cadences"],
        "mom_centr1_median": centr1["median"],
        "mom_centr1_peak_to_peak": centr1["peak_to_peak"],
        "mom_centr2_cadences": centr2["cadences"],
        "mom_centr2_median": centr2["median"],
        "mom_centr2_peak_to_peak": centr2["peak_to_peak"],
    }


def _totals(windows: list[dict[str, Any]]) -> dict[str, int]:
    totals = {"windows": len(windows)}
    for key in COUNT_COLUMNS:
        totals[key] = int(sum(row[key] for row in windows))
    return totals


def _reconciliation(
    windows: list[dict[str, Any]], totals: dict[str, int]
) -> dict[str, bool]:
    per_window = all(
        row["cadences_before_mask"]
        == row["cadences_after_mask"] + row["removed_by_quality_mask"]
        and row["finite_sap_after_mask"] <= row["cadences_after_mask"]
        and row["finite_pdcsap_after_mask"] <= row["cadences_after_mask"]
        and row["removed_by_quality_mask"] >= 0
        for row in windows
    )
    totals_match = all(
        totals[key] == sum(row[key] for row in windows) for key in COUNT_COLUMNS
    )
    disjoint = all(
        right["window_start_btjd"] > left["window_end_btjd"]
        for left, right in itertools.pairwise(windows)
    )
    return {
        "per_window_counts_reconcile": per_window,
        "window_totals_match_sum": totals_match,
        "windows_disjoint": disjoint,
    }


def build_record(
    lightcurve_path: Path = DEFAULT_LIGHTCURVE,
    tpf_path: Path = DEFAULT_TPF,
    windows_path: Path = DEFAULT_WINDOWS,
) -> dict[str, Any]:
    """Build the per-window quality record after verifying both FITS inputs."""

    lightcurve_input = _preflight_fits(
        lightcurve_path, EXPECTED_LIGHTCURVE_SHA256, "light curve"
    )
    tpf_input = _preflight_fits(tpf_path, EXPECTED_TPF_SHA256, "target pixel file")
    windows_input = _read_windows(windows_path)

    columns = _read_lightcurve(lightcurve_path)
    tpf_time = _read_tpf_time(tpf_path)
    time_matches = bool(
        columns["TIME"].shape == tpf_time.shape
        and np.array_equal(columns["TIME"], tpf_time, equal_nan=True)
    )
    if not time_matches:
        raise WindowQCError(
            "the light curve and target pixel file TIME columns do not match"
        )
    lightcurve_input["time_matches_tpf"] = time_matches

    finite_time = np.isfinite(columns["TIME"])
    quality_pass = (
        finite_time
        & columns["_quality_finite"]
        & ((columns["QUALITY"] & EXPECTED_MASK_BITS) == 0)
    )

    windows = [
        _window_record(window, columns, quality_pass, finite_time)
        for window in windows_input["windows"]
    ]
    totals = _totals(windows)
    reconciliation = _reconciliation(windows, totals)
    if not all(reconciliation.values()):
        raise WindowQCError("per-window counts do not reconcile to the totals")

    return {
        "schema_version": SCHEMA_VERSION,
        "task": TASK,
        "status": "OK",
        "label": f"sector {EXPECTED_SECTOR} per-window quality record for {TARGET_NAME}",
        "target": {
            "name": TARGET_NAME,
            "tic_id": EXPECTED_TIC_ID,
            "sector": EXPECTED_SECTOR,
        },
        "inputs": {
            "lightcurve": lightcurve_input,
            "tpf": tpf_input,
            "windows": windows_input,
        },
        "mask_bits": EXPECTED_MASK_BITS,
        "mask_bit_values": [1, 2, 4, 8, 32, 128],
        "window_definition": {
            "source": windows_input["path"],
            "centre_field": "btjd from the committed predictions",
            "half_width_days": WINDOW_HALF_WIDTH_DAYS,
            "interval": ("centre_btjd - 0.5 <= BTJD <= centre_btjd + 0.5"),
            "inclusive": True,
            "note": (
                "The 0.5 d half-width is a fixed diagnostic interval, not a "
                "transit duration and not a fit. Window centres are read from "
                "the committed prediction record, not re-derived."
            ),
        },
        "cadence_definition": (
            "Cadence counts use rows with finite TIME. The quality cut removes "
            "a cadence when (QUALITY & 175) != 0. Flux counts require a finite "
            "value in that column after the quality cut."
        ),
        "nominal_cadence_seconds": EXPECTED_CADENCE_SECONDS,
        "gap_factor": GAP_FACTOR,
        "gap_threshold_seconds": GAP_THRESHOLD_SECONDS,
        "gap_scope": "finite TIME values before the quality mask",
        "window_totals": totals,
        "reconciliation": reconciliation,
        "windows": windows,
        "limits": list(LIMITS),
    }


def _format_float(value: float, digits: int) -> str:
    return f"{value:.{digits}f}"


def _optional(value: float | None, format_spec: str) -> str:
    return "n/a" if value is None else format_spec.format(value)


def render_markdown(record: dict[str, Any]) -> str:
    """Render the per-window quality record as Markdown."""

    target = record["target"]
    inputs = record["inputs"]
    definition = record["window_definition"]
    totals = record["window_totals"]
    windows = record["windows"]

    lines: list[str] = []
    lines.append(f"# TESS sector {target['sector']} per-window quality record")
    lines.append("")
    lines.append(
        "Status: diagnostic record for human verification. This document "
        "describes cadence availability and per-window column statistics in "
        "the committed sector 56 light curve. It is not a detection or a "
        "transit timing measurement."
    )
    lines.append("")
    lines.append(
        f"Target: {target['name']} (TIC {target['tic_id']}), sector {target['sector']}."
    )
    lines.append("")
    lines.append("## Inputs")
    lines.append("")
    lines.append("| Input | Repository path | SHA-256 |")
    lines.append("| --- | --- | --- |")
    lines.append(
        f"| Light curve | `{inputs['lightcurve']['path']}` | "
        f"`{inputs['lightcurve']['sha256']}` |"
    )
    lines.append(
        f"| Target pixel file | `{inputs['tpf']['path']}` | "
        f"`{inputs['tpf']['sha256']}` |"
    )
    lines.append(
        f"| Windows record | `{inputs['windows']['path']}` | "
        f"`{inputs['windows']['sha256']}` |"
    )
    lines.append("")
    lines.append(
        "Both FITS primary headers report TICID "
        f"{inputs['lightcurve']['tic_id']} and SECTOR "
        f"{inputs['lightcurve']['sector']}. The light curve and target pixel "
        "file TIME columns match row for row. The target pixel file is used "
        "for identity and time coverage only; no pixel photometry, "
        "contamination ratio, or nearby-source count is computed."
    )
    lines.append("")
    lines.append("## Window definition")
    lines.append("")
    lines.append(
        f"- Window source: `{definition['source']}` (the eight predicted "
        "mid-transit times, read without re-deriving them)."
    )
    lines.append(
        f"- Each window is "
        f"`centre_btjd - {definition['half_width_days']!r} d <= BTJD <= "
        f"centre_btjd + {definition['half_width_days']!r} d`, inclusive."
    )
    lines.append(
        "- The 0.5 d half-width is a fixed diagnostic interval. It is not a "
        "transit duration and it is not fitted to the data."
    )
    lines.append(
        f"- Quality mask: {record['mask_bits']} (bits "
        f"{', '.join(str(bit) for bit in record['mask_bit_values'])})."
    )
    lines.append(
        f"- Nominal cadence: {record['nominal_cadence_seconds']:.1f} s; gap "
        f"threshold: {record['gap_threshold_seconds']:.1f} s "
        f"({record['gap_factor']} times the cadence)."
    )
    lines.append(f"- Gaps are counted on {record['gap_scope']} inside each window.")
    lines.append("")
    lines.append("## Coverage per window")
    lines.append("")
    lines.append(
        "| Cycle | Centre BTJD | Window start | Window end | Before mask | "
        "After mask | Removed | SAP finite | PDCSAP finite | Gaps | "
        "Longest gap (s) |"
    )
    lines.append(
        "| ---: | ---: | ---: | ---: | ---: | ---: | ---: | ---: | ---: | ---: | ---: |"
    )
    for row in windows:
        lines.append(
            f"| {row['cycle']} | {_format_float(row['btjd'], 6)} | "
            f"{_format_float(row['window_start_btjd'], 6)} | "
            f"{_format_float(row['window_end_btjd'], 6)} | "
            f"{row['cadences_before_mask']:,} | "
            f"{row['cadences_after_mask']:,} | "
            f"{row['removed_by_quality_mask']:,} | "
            f"{row['finite_sap_after_mask']:,} | "
            f"{row['finite_pdcsap_after_mask']:,} | "
            f"{row['gap_count']} | "
            f"{row['longest_gap_seconds']:.6f} |"
        )
    lines.append(
        f"| All | - | - | - | {totals['cadences_before_mask']:,} | "
        f"{totals['cadences_after_mask']:,} | "
        f"{totals['removed_by_quality_mask']:,} | "
        f"{totals['finite_sap_after_mask']:,} | "
        f"{totals['finite_pdcsap_after_mask']:,} | "
        f"{totals['gap_count']} | - |"
    )
    lines.append("")
    lines.append(
        "The per-window counts reconcile: before mask equals after mask plus "
        "removed for every window, and each column of the totals equals the "
        "sum over the eight non-overlapping windows."
    )
    lines.append("")
    lines.append("## Background per window")
    lines.append("")
    lines.append(
        "SAP_BKG over cadences that pass the quality mask and have a finite "
        "background value."
    )
    lines.append("")
    lines.append(
        "| Cycle | Cadences | Median | 5th percentile | 95th percentile | "
        "5 to 95 range |"
    )
    lines.append("| ---: | ---: | ---: | ---: | ---: | ---: |")
    for row in windows:
        lines.append(
            f"| {row['cycle']} | {row['sap_bkg_cadences']:,} | "
            f"{_optional(row['sap_bkg_median'], '{:.6f}')} | "
            f"{_optional(row['sap_bkg_percentile_5'], '{:.6f}')} | "
            f"{_optional(row['sap_bkg_percentile_95'], '{:.6f}')} | "
            f"{_optional(row['sap_bkg_range_5_95'], '{:.6f}')} |"
        )
    lines.append("")
    lines.append("## Centroid per window")
    lines.append("")
    lines.append(
        "MOM_CENTR1 and MOM_CENTR2 over cadences that pass the quality mask "
        "and have a finite value. Peak-to-peak is the maximum minus the "
        "minimum in the window."
    )
    lines.append("")
    lines.append(
        "| Cycle | MOM_CENTR1 cadences | MOM_CENTR1 median | "
        "MOM_CENTR1 peak-to-peak | MOM_CENTR2 cadences | "
        "MOM_CENTR2 median | MOM_CENTR2 peak-to-peak |"
    )
    lines.append("| ---: | ---: | ---: | ---: | ---: | ---: | ---: |")
    for row in windows:
        lines.append(
            f"| {row['cycle']} | {row['mom_centr1_cadences']:,} | "
            f"{_optional(row['mom_centr1_median'], '{:.6f}')} | "
            f"{_optional(row['mom_centr1_peak_to_peak'], '{:.6f}')} | "
            f"{row['mom_centr2_cadences']:,} | "
            f"{_optional(row['mom_centr2_median'], '{:.6f}')} | "
            f"{_optional(row['mom_centr2_peak_to_peak'], '{:.6f}')} |"
        )
    lines.append("")
    lines.append("## Reproduction")
    lines.append("")
    lines.append("```")
    lines.append("python research/tess-timing/scripts/real_window_qc.py")
    lines.append("python research/tess-timing/scripts/real_window_qc.py --check")
    lines.append(
        "python -m pytest research/tess-timing/tests/test_real_window_qc.py -q"
    )
    lines.append("```")
    lines.append("")
    lines.append("## Limits")
    lines.append("")
    for limit in record["limits"]:
        lines.append(f"- {limit}")
    lines.append("")
    lines.append(
        "Generated by `research/tess-timing/scripts/real_window_qc.py` from "
        "the committed sector 56 products and prediction record. Offline; no "
        "archive, catalogue, or network request is made."
    )
    lines.append("")
    return "\n".join(lines)


def _csv_value(value: object) -> str:
    if value is None:
        return ""
    if isinstance(value, bool):
        return "true" if value else "false"
    if isinstance(value, int):
        return str(value)
    if isinstance(value, float):
        return repr(value)
    return str(value)


def write_csv(record: dict[str, Any], path: Path) -> Path:
    path = Path(path)
    path.parent.mkdir(parents=True, exist_ok=True)
    with path.open("w", encoding="utf-8", newline="") as handle:
        writer = csv.writer(handle, lineterminator="\n")
        writer.writerow(CSV_COLUMNS)
        for row in record["windows"]:
            writer.writerow([_csv_value(row[name]) for name in CSV_COLUMNS])
    return path


def write_json(record: dict[str, Any], path: Path) -> Path:
    path = Path(path)
    path.parent.mkdir(parents=True, exist_ok=True)
    path.write_text(
        json.dumps(record, indent=2, ensure_ascii=False) + "\n",
        encoding="utf-8",
    )
    return path


def write_markdown(record: dict[str, Any], path: Path) -> Path:
    path = Path(path)
    path.parent.mkdir(parents=True, exist_ok=True)
    path.write_text(render_markdown(record), encoding="utf-8")
    return path


def write_outputs(record: dict[str, Any], directory: Path) -> dict[str, Path]:
    directory = Path(directory)
    return {
        "csv": write_csv(record, directory / "windows.csv"),
        "json": write_json(record, directory / "windows.json"),
        "markdown": write_markdown(record, directory / "TESS-REAL-WINDOW-QC.md"),
    }


def check_outputs(
    csv_path: Path = DEFAULT_CSV,
    json_path: Path = DEFAULT_JSON,
    markdown_path: Path = DEFAULT_MARKDOWN,
) -> dict[str, Path]:
    """Regenerate all three files under .scratch and compare bytes."""

    committed = {
        "csv": Path(csv_path),
        "json": Path(json_path),
        "markdown": Path(markdown_path),
    }
    for name, path in committed.items():
        if not path.exists():
            raise WindowQCError(f"committed {name} output does not exist: {path}")

    if CHECK_ROOT.exists():
        shutil.rmtree(CHECK_ROOT)
    CHECK_ROOT.mkdir(parents=True, exist_ok=True)
    record = build_record()
    regenerated = write_outputs(record, CHECK_ROOT)

    issues: list[str] = []
    for name, path in committed.items():
        if path.read_bytes() != regenerated[name].read_bytes():
            issues.append(f"{name}: committed bytes differ from the regenerated bytes")
    if issues:
        raise WindowQCError(
            "committed outputs do not match the regenerated output: "
            + "; ".join(issues)
        )
    return regenerated


def _parser() -> argparse.ArgumentParser:
    parser = argparse.ArgumentParser(
        description=(
            "Write a per-window quality record for the committed sector 56 "
            "light curve. The eight window centres are read from the "
            "committed prediction record. The command verifies both FITS "
            "digests and headers, masks at 175, and reports per-window "
            "cadence, gap, background, and centroid statistics. It does not "
            "detrend, fit, fold, or measure a transit time."
        )
    )
    parser.add_argument(
        "--check",
        action="store_true",
        help=(
            "regenerate all three files under .scratch/real-window-qc-check "
            "and compare them with the committed outputs byte for byte"
        ),
    )
    return parser


def main(argv: list[str] | None = None) -> int:
    arguments = _parser().parse_args(argv)

    if arguments.check:
        try:
            check_outputs()
        except (WindowQCError, OSError) as error:
            print(f"error: {error}", file=sys.stderr)
            return 1
        print(
            f"status=OK check=byte_identical byte_identical=true "
            f"task={TASK} target={TARGET_NAME} tic_id={EXPECTED_TIC_ID} "
            f"sector={EXPECTED_SECTOR} windows={EXPECTED_WINDOW_COUNT} "
            f"mask_bits={EXPECTED_MASK_BITS} "
            f"wrote={_display(CHECK_ROOT / 'windows.csv')},"
            f"{_display(CHECK_ROOT / 'windows.json')},"
            f"{_display(CHECK_ROOT / 'TESS-REAL-WINDOW-QC.md')}"
        )
        return 0

    try:
        record = build_record()
        write_csv(record, DEFAULT_CSV)
        write_json(record, DEFAULT_JSON)
        write_markdown(record, DEFAULT_MARKDOWN)
    except (WindowQCError, OSError) as error:
        print(f"error: {error}", file=sys.stderr)
        return 1

    totals = record["window_totals"]
    print(
        f"status=OK task={TASK} target={TARGET_NAME} "
        f"tic_id={EXPECTED_TIC_ID} sector={EXPECTED_SECTOR} "
        f"windows={totals['windows']} mask_bits={EXPECTED_MASK_BITS} "
        f"cadences_before_mask={totals['cadences_before_mask']} "
        f"cadences_after_mask={totals['cadences_after_mask']} "
        f"finite_sap={totals['finite_sap_after_mask']} "
        f"finite_pdcsap={totals['finite_pdcsap_after_mask']} "
        f"gaps={totals['gap_count']} "
        f"reconciles={str(all(record['reconciliation'].values())).lower()} "
        f"wrote={_display(DEFAULT_CSV)},{_display(DEFAULT_JSON)},"
        f"{_display(DEFAULT_MARKDOWN)}"
    )
    return 0


if __name__ == "__main__":
    raise SystemExit(main())