# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: percent
#       format_version: '1.3'
#       jupytext_version: 1.19.4
#   kernelspec:
#     display_name: Python 3 (ipykernel)
#     language: python
#     name: python3
# ---

# %% [markdown]
# # TOI-791 b recovery and a publication-informed c-period ranking
#
# **Question:** Can a compact, transparent Box Least Squares (BLS) pipeline
# broadly recover TOI-791 b and show how the two c-like TESS event windows rank
# under a publication-informed period search?
#
# This notebook deliberately stops short of claiming an independent planet
# confirmation. The peer-reviewed study used detailed transit-timing models,
# ground-based follow-up, high-resolution imaging, and spectroscopy. Here, the
# narrower goal is to test whether the periods themselves are visible in the
# public space-based photometry.
#
# Primary references:
#
# - Peer-reviewed study: https://doi.org/10.1093/mnras/stag864
# - NASA announcement: https://science.nasa.gov/missions/tess/nasas-tess-mission-reveals-the-puffiest-planets-ever-found/
# - TESS archive/data products: https://archive.stsci.edu/missions-and-data/tess/data-products
# - ExoFOP target page: https://exofop.ipac.caltech.edu/tess/target.php?id=306472057

# %%
from __future__ import annotations

import csv
import hashlib
import json
import warnings
from datetime import datetime, timezone
from importlib.metadata import version
from pathlib import Path

warnings.filterwarnings("ignore", message="Warning: the tpfmodel")

import astropy.units as u
import lightkurve as lk
import matplotlib.pyplot as plt
import numpy as np
from astropy.timeseries import BoxLeastSquares

PROJECT_ROOT = Path(__file__).resolve().parent if "__file__" in globals() else Path.cwd()
CACHE_DIR = PROJECT_ROOT / ".cache" / "lightkurve"
OUTPUT_DIR = PROJECT_ROOT / "outputs"
CACHE_DIR.mkdir(parents=True, exist_ok=True)
OUTPUT_DIR.mkdir(parents=True, exist_ok=True)

TARGET = "TIC 306472057"  # TOI-791
QUERY_TIMESTAMP_UTC = datetime.now(timezone.utc).isoformat(timespec="seconds")
CADENCE_SECONDS = 120
BIN_SIZE = 30 * u.min
FLATTEN_WINDOW = 4001  # about 5.6 days at two-minute cadence
SEARCH_DURATIONS = np.linspace(0.30, 0.90, 9) * u.day
MASK_HALF_WIDTH_DAYS = 0.80

PUBLISHED = {
    "b": {
        "period_days": 139.29931,
        "duration_hours": 12.366,
        "depth_ppt": 5.428,
        "radius_jupiter": 0.993,
        "mass_earth": 9.5,
        "density_g_cm3": 0.037,
    },
    "c": {
        "period_days": 232.01570,
        "duration_hours": 11.926,
        "depth_ppt": 6.78,
        "radius_jupiter": 1.155,
        "mass_earth": 18.6,
        "density_g_cm3": 0.047,
    },
}

print(f"Lightkurve {lk.__version__}")
print(f"Target: {TARGET}")
print("Archive query details are recorded in the machine-readable manifest.")


# %% [markdown]
# ## 1. Discover and download the public TESS products
#
# The paper reports 44 observed sectors, including long-cadence products. For a
# uniform analysis, this notebook keeps only official SPOC light curves at
# two-minute cadence. `quality_bitmask="default"` removes cadences carrying the
# standard set of serious TESS quality flags. Lightkurve selects `PDCSAP_FLUX`
# for these products, which has already been corrected for common instrumental
# trends and crowding effects by the mission pipeline.

# %%
search_result = lk.search_lightcurve(TARGET, mission="TESS", author="SPOC")
two_minute_mask = np.isclose(np.asarray(search_result.table["exptime"], dtype=float), CADENCE_SECONDS)
two_minute_result = search_result[two_minute_mask]

if len(two_minute_result) == 0:
    raise RuntimeError("MAST returned no two-minute SPOC light curves for TOI-791.")

lightcurves = two_minute_result.download_all(
    download_dir=str(CACHE_DIR),
    quality_bitmask="default",
)
if lightcurves is None or len(lightcurves) != len(two_minute_result):
    raise RuntimeError("One or more TESS products did not download successfully.")


def sector_number(mission_name: str) -> int:
    """Extract the sector number from a string such as 'TESS Sector 04'."""

    return int(str(mission_name).split()[-1])


sectors = [sector_number(value) for value in two_minute_result.table["mission"]]
print(f"Downloaded {len(lightcurves)} two-minute SPOC light curves.")
print("Sectors:", sectors)

# Save a compact data manifest so the exact archive response is auditable.
manifest_path = OUTPUT_DIR / "toi791_tess_product_manifest.csv"
manifest_columns = [
    name
    for name in ("mission", "author", "exptime", "year", "productFilename", "dataURI")
    if name in two_minute_result.table.colnames
]
with manifest_path.open("w", newline="", encoding="utf-8") as handle:
    writer = csv.writer(handle)
    writer.writerow([*manifest_columns, "sha256", "query_timestamp_utc"])
    for row, lightcurve in zip(two_minute_result.table, lightcurves, strict=True):
        local_path = Path(lightcurve.filename)
        digest = hashlib.sha256(local_path.read_bytes()).hexdigest()
        writer.writerow(
            [*[row[name] for name in manifest_columns], digest, QUERY_TIMESTAMP_UTC]
        )


# %% [markdown]
# ## 2. Clean each sector without erasing the long transit
#
# Each sector is processed separately:
#
# 1. remove missing values;
# 2. normalize its median flux to one;
# 3. remove slow trends with a roughly 5.6-day Savitzky-Golay window;
# 4. clip strong upward outliers while using a loose lower threshold so real
#    transits are retained;
# 5. stitch the sectors and bin to 30 minutes.
#
# The planets' published transits last about 12 hours, so the detrending window
# is intentionally much longer than a transit. Even so, this processing can
# attenuate transit depths. The analysis therefore treats **period recovery**,
# not depth measurement, as its primary result.

# %%
def valid_odd_window(length: int, requested: int = FLATTEN_WINDOW) -> int:
    """Return a valid odd Savitzky-Golay window shorter than the light curve."""

    candidate = min(requested, length - 1 if length % 2 == 0 else length - 2)
    candidate = candidate if candidate % 2 == 1 else candidate - 1
    if candidate < 101:
        raise ValueError(f"Light curve is too short to flatten safely ({length} cadences).")
    return candidate


cleaned_sectors = []
for lightcurve in lightcurves:
    if lightcurve.meta.get("FLUX_ORIGIN") != "pdcsap_flux":
        raise RuntimeError("Expected PDCSAP flux, but Lightkurve selected another flux column.")

    cleaned = lightcurve.remove_nans().normalize()
    cleaned = cleaned.flatten(
        window_length=valid_odd_window(len(cleaned)),
        polyorder=2,
        break_tolerance=5,
        niters=3,
        sigma=5,
    )
    cleaned = cleaned.remove_outliers(sigma_lower=20, sigma_upper=5)
    cleaned_sectors.append(cleaned)

stitched = lk.LightCurveCollection(cleaned_sectors).stitch().remove_nans()
binned = stitched.bin(time_bin_size=BIN_SIZE).remove_nans()

time_days = np.asarray(binned.time.value, dtype=float)
relative_flux = np.asarray(binned.flux.value, dtype=float)
relative_flux_error = np.asarray(binned.flux_err.value, dtype=float)
finite = (
    np.isfinite(time_days)
    & np.isfinite(relative_flux)
    & np.isfinite(relative_flux_error)
    & (relative_flux_error > 0)
)
time_days = time_days[finite]
relative_flux = relative_flux[finite]
relative_flux_error = relative_flux_error[finite]

# Prevent a small number of unusually tiny pipeline uncertainties from
# dominating the BLS statistic.
relative_flux_error = np.maximum(relative_flux_error, np.percentile(relative_flux_error, 5))

print(f"Cleaned cadences: {len(stitched):,}")
print(f"Thirty-minute bins: {len(time_days):,}")
print(f"Time baseline: {np.ptp(time_days):,.1f} days")


# %% [markdown]
# ## 3. Search for the strongest long-period box-shaped signal
#
# Box Least Squares tests many candidate periods and transit durations for a
# repeating, box-shaped dip. The first search spans 80–280 days. Its power is a
# useful ranking statistic, but it is **not** a formal discovery significance:
# it does not include all look-elsewhere effects, instrumental systematics, or
# astrophysical false-positive tests. A dense local grid then refines the
# numerical maximum, but no formal period uncertainty is estimated.

# %%
def scalar(value) -> float:
    """Convert an Astropy scalar or NumPy scalar to a Python float."""

    return float(value.value if hasattr(value, "value") else value)


def run_bls(
    time: np.ndarray,
    flux: np.ndarray,
    flux_error: np.ndarray,
    minimum_period: float,
    maximum_period: float,
    grid_size: int,
):
    periods = np.linspace(minimum_period, maximum_period, grid_size) * u.day
    model = BoxLeastSquares(time * u.day, flux, flux_error)
    result = model.power(periods, SEARCH_DURATIONS, objective="snr")
    best_index = int(np.nanargmax(np.asarray(result.power)))
    return model, result, best_index


def refine_bls_peak(
    model: BoxLeastSquares,
    coarse_result,
    coarse_best_index: int,
    half_width_days: float = 0.10,
    grid_size: int = 4001,
):
    """Refine a coarse BLS maximum while keeping its preferred duration fixed.

    The dense grid locates the numerical peak more accurately, but it does not
    provide a formal period uncertainty. Results are therefore reported to only
    two decimal places in human-facing text.
    """

    center = scalar(coarse_result.period[coarse_best_index])
    duration = scalar(coarse_result.duration[coarse_best_index]) * u.day
    periods = np.linspace(center - half_width_days, center + half_width_days, grid_size) * u.day
    result = model.power(periods, duration, objective="snr")
    best_index = int(np.nanargmax(np.asarray(result.power)))
    return result, best_index


def recovered_signal(name: str, result, best_index: int, coarse_spacing_days: float) -> dict:
    """Create a labeled record for a refined BLS peak."""

    return {
        "planet": name,
        "period_days": scalar(result.period[best_index]),
        "duration_days": scalar(result.duration[best_index]),
        "transit_time_btjd": scalar(result.transit_time[best_index]),
        "depth_relative_flux": scalar(result.depth[best_index]),
        "bls_depth_snr": scalar(result.power[best_index]),
        "coarse_grid_spacing_days": coarse_spacing_days,
        "reporting_precision_days": 0.01,
        "formal_period_uncertainty_estimated": False,
    }


bls_b, periodogram_b, best_b_index = run_bls(
    time_days,
    relative_flux,
    relative_flux_error,
    minimum_period=80,
    maximum_period=280,
    grid_size=8000,
)
periodogram_b_refined, best_b_refined_index = refine_bls_peak(
    bls_b, periodogram_b, best_b_index
)

recovered_b = recovered_signal(
    "TOI-791 b",
    periodogram_b_refined,
    best_b_refined_index,
    coarse_spacing_days=200 / (8000 - 1),
)

print("First-pass strongest signal:", recovered_b)


# %% [markdown]
# ## 4. Mask TOI-791 b, then search for a second signal
#
# The first planet creates harmonics and aliases in an irregularly sampled
# survey. In a deliberately publication-informed step, remove a ±0.8-day window
# around each predicted TOI-791 b transit and repeat the BLS search from 170–280
# days. This ranks candidate periods; it does not independently resolve every
# integer-cycle alias for the sparsely observed c-like events.

# %%
def signed_phase_days(time: np.ndarray, period: float, transit_time: float) -> np.ndarray:
    """Return signed phase in days, centered on a transit at phase zero."""

    return ((time - transit_time + 0.5 * period) % period) - 0.5 * period


phase_b_days = signed_phase_days(
    time_days,
    recovered_b["period_days"],
    recovered_b["transit_time_btjd"],
)
keep_after_b_mask = np.abs(phase_b_days) > MASK_HALF_WIDTH_DAYS

bls_c, periodogram_c, best_c_index = run_bls(
    time_days[keep_after_b_mask],
    relative_flux[keep_after_b_mask],
    relative_flux_error[keep_after_b_mask],
    minimum_period=170,
    maximum_period=280,
    grid_size=7000,
)
periodogram_c_refined, best_c_refined_index = refine_bls_peak(
    bls_c, periodogram_c, best_c_index
)

recovered_c = recovered_signal(
    "TOI-791 c",
    periodogram_c_refined,
    best_c_refined_index,
    coarse_spacing_days=110 / (7000 - 1),
)

print("Second-pass strongest signal:", recovered_c)


def observed_event_windows(
    time: np.ndarray,
    flux: np.ndarray,
    signal: dict,
    minimum_depth_ppt: float = 2.0,
) -> list[dict]:
    """Identify observed transit-like windows aligned by a candidate period.

    This is an event inventory, not a planet-validation test. It makes the
    sparse support for the long-period c signal explicit.
    """

    phase = signed_phase_days(time, signal["period_days"], signal["transit_time_btjd"])
    cycles = np.rint(
        (time - signal["transit_time_btjd"]) / signal["period_days"]
    ).astype(int)
    events = []
    for cycle in np.unique(cycles[np.abs(phase) < 0.75]):
        in_transit = (cycles == cycle) & (np.abs(phase) < 0.25)
        local_baseline = (
            (cycles == cycle)
            & (np.abs(phase) >= 0.60)
            & (np.abs(phase) < 0.75)
        )
        if in_transit.sum() < 5 or local_baseline.sum() < 5:
            continue
        depth_ppt = (
            np.median(flux[local_baseline]) - np.median(flux[in_transit])
        ) * 1000
        if depth_ppt >= minimum_depth_ppt:
            events.append(
                {
                    "cycle_index": int(cycle),
                    "predicted_center_btjd": signal["transit_time_btjd"]
                    + cycle * signal["period_days"],
                    "in_transit_bins": int(in_transit.sum()),
                    "local_baseline_bins": int(local_baseline.sum()),
                    "local_median_depth_ppt": float(depth_ppt),
                }
            )
    return events


c_event_windows = observed_event_windows(time_days, relative_flux, recovered_c)
c_events_path = OUTPUT_DIR / "toi791_c_event_windows.csv"
c_event_fieldnames = [
    "cycle_index",
    "predicted_center_btjd",
    "in_transit_bins",
    "local_baseline_bins",
    "local_median_depth_ppt",
]
with c_events_path.open("w", newline="", encoding="utf-8") as handle:
    writer = csv.DictWriter(handle, fieldnames=c_event_fieldnames)
    writer.writeheader()
    writer.writerows(c_event_windows)

print(f"TOI-791 c event windows in selected data: {len(c_event_windows)}")
for event in c_event_windows:
    print(event)


# %% [markdown]
# ## 5. Compare the recovered periods with the peer-reviewed values
#
# A dense local grid refines each numerical peak after the broad search, but this
# notebook does not estimate formal period uncertainties. Human-facing results
# are therefore rounded to 0.01 day. The study also reports transit-timing
# variations of up to 50 minutes, and its physical model is much more
# sophisticated than the box model used here.

# %%
recovered = {"b": recovered_b, "c": recovered_c}
comparison_rows = []
for label in ("b", "c"):
    measured = recovered[label]
    reference = PUBLISHED[label]
    comparison_rows.append(
        {
            "planet": measured["planet"],
            "refined_bls_peak_period_days": measured["period_days"],
            "reported_recovered_period_days": round(measured["period_days"], 2),
            "published_period_days": reference["period_days"],
            "absolute_peak_difference_days": abs(
                measured["period_days"] - reference["period_days"]
            ),
            "coarse_grid_spacing_days": measured["coarse_grid_spacing_days"],
            "reporting_precision_days": measured["reporting_precision_days"],
            "formal_period_uncertainty_estimated": measured[
                "formal_period_uncertainty_estimated"
            ],
            "bls_box_depth_ppt": measured["depth_relative_flux"] * 1000,
            "published_transit_depth_ppt": reference["depth_ppt"],
            "bls_box_duration_hours": measured["duration_days"] * 24,
            "published_duration_hours": reference["duration_hours"],
            "bls_depth_snr_search_statistic": measured["bls_depth_snr"],
        }
    )

results_path = OUTPUT_DIR / "toi791_recovery_results.csv"
with results_path.open("w", newline="", encoding="utf-8") as handle:
    writer = csv.DictWriter(handle, fieldnames=list(comparison_rows[0].keys()))
    writer.writeheader()
    writer.writerows(comparison_rows)

summary_path = OUTPUT_DIR / "toi791_recovery_summary.json"
package_versions = {
    package: version(package)
    for package in ("astropy", "lightkurve", "matplotlib", "numpy")
}
summary = {
    "target": TARGET,
    "query_timestamp_utc": QUERY_TIMESTAMP_UTC,
    "package_versions": package_versions,
    "flux_column": "pdcsap_flux",
    "number_of_two_minute_spoc_products": len(lightcurves),
    "sectors": sectors,
    "cleaned_cadences": len(stitched),
    "thirty_minute_bins": len(time_days),
    "baseline_days": float(np.ptp(time_days)),
    "results": comparison_rows,
    "toi791_c_event_windows": c_event_windows,
    "scientific_scope": (
        "Approximate period-peak recovery from public TESS photometry. The c "
        "result is a publication-informed ranking supported by two TESS event "
        "windows, not independent alias resolution or planet confirmation."
    ),
}
summary_path.write_text(json.dumps(summary, indent=2), encoding="utf-8")

for row in comparison_rows:
    print(
        f"{row['planet']}: refined numerical peak "
        f"{row['refined_bls_peak_period_days']:.6f} d; report as "
        f"{row['reported_recovered_period_days']:.2f} d; published "
        f"{row['published_period_days']:.5f} d; no formal uncertainty estimated"
    )


# %% [markdown]
# ## 6. Publication-quality diagnostic chart
#
# The top panel shows the actual observation coverage. The middle panels show
# the two BLS searches. The bottom panels fold all observations onto each
# recovered period so repeated transit events line up at zero hours.

# %%
def robust_phase_bins(
    phase_hours: np.ndarray,
    flux_ppt: np.ndarray,
    limit_hours: float = 36,
    width_hours: float = 0.75,
):
    """Median-bin a folded light curve without assuming independent errors."""

    edges = np.arange(-limit_hours, limit_hours + width_hours, width_hours)
    centers = 0.5 * (edges[:-1] + edges[1:])
    bin_index = np.digitize(phase_hours, edges) - 1
    medians = np.full_like(centers, np.nan, dtype=float)
    counts = np.zeros_like(centers, dtype=int)
    for index in range(len(centers)):
        values = flux_ppt[bin_index == index]
        values = values[np.isfinite(values)]
        counts[index] = len(values)
        if len(values) >= 3:
            medians[index] = np.median(values)
    valid = np.isfinite(medians)
    return centers[valid], medians[valid], counts[valid]


NAVY = "#15253D"
TEAL = "#148A8A"
CORAL = "#DB5A42"
GOLD = "#D6A21E"
LIGHT = "#DDE6EC"

plt.rcParams.update(
    {
        "font.family": "DejaVu Sans",
        "font.size": 10,
        "axes.titlesize": 12,
        "axes.labelsize": 10,
        "axes.edgecolor": NAVY,
        "axes.labelcolor": NAVY,
        "xtick.color": NAVY,
        "ytick.color": NAVY,
        "text.color": NAVY,
    }
)

figure = plt.figure(figsize=(14, 13), constrained_layout=True)
grid = figure.add_gridspec(3, 2, height_ratios=[0.8, 1, 1])
ax_coverage = figure.add_subplot(grid[0, :])
ax_bls_b = figure.add_subplot(grid[1, 0])
ax_bls_c = figure.add_subplot(grid[1, 1])
ax_fold_b = figure.add_subplot(grid[2, 0])
ax_fold_c = figure.add_subplot(grid[2, 1])

flux_ppt = (relative_flux - 1.0) * 1000
coverage_low, coverage_high = np.percentile(flux_ppt, [0.5, 99.5])
ax_coverage.scatter(time_days, flux_ppt, s=2, color=NAVY, alpha=0.22, rasterized=True)
ax_coverage.set(
    title=f"Public TESS coverage: {len(lightcurves)} two-minute SPOC products",
    xlabel="Time (BTJD = BJD - 2,457,000)",
    ylabel="Detrended flux (ppt)",
    ylim=(coverage_low, coverage_high),
)
ax_coverage.grid(color=LIGHT, linewidth=0.7, alpha=0.8)

period_b_values = np.asarray(periodogram_b.period.value, dtype=float)
power_b_values = np.asarray(periodogram_b.power, dtype=float)
ax_bls_b.plot(period_b_values, power_b_values, color=TEAL, linewidth=1.2)
ax_bls_b.axvline(PUBLISHED["b"]["period_days"], color=NAVY, linestyle=":", linewidth=1.6)
ax_bls_b.scatter(
    [recovered_b["period_days"]],
    [recovered_b["bls_depth_snr"]],
    color=CORAL,
    s=50,
    zorder=5,
    label=f"Approx. peak: {recovered_b['period_days']:.2f} d",
)
ax_bls_b.set(
    title="First BLS search: strongest signal is TOI-791 b",
    xlabel="Candidate period (days)",
    ylabel="BLS depth S/N (search statistic)",
)
ax_bls_b.legend(frameon=False, loc="upper right")
ax_bls_b.grid(color=LIGHT, linewidth=0.7, alpha=0.8)

period_c_values = np.asarray(periodogram_c.period.value, dtype=float)
power_c_values = np.asarray(periodogram_c.power, dtype=float)
ax_bls_c.plot(period_c_values, power_c_values, color=GOLD, linewidth=1.2)
ax_bls_c.axvline(PUBLISHED["c"]["period_days"], color=NAVY, linestyle=":", linewidth=1.6)
ax_bls_c.scatter(
    [recovered_c["period_days"]],
    [recovered_c["bls_depth_snr"]],
    color=CORAL,
    s=50,
    zorder=5,
    label=f"Approx. peak: {recovered_c['period_days']:.2f} d",
)
ax_bls_c.set(
    title="Targeted search after masking b: c ranks highest",
    xlabel="Candidate period (days)",
    ylabel="BLS depth S/N (search statistic)",
)
ax_bls_c.legend(frameon=False, loc="upper right")
ax_bls_c.grid(color=LIGHT, linewidth=0.7, alpha=0.8)


def draw_folded_panel(
    axis,
    label: str,
    measured: dict,
    color: str,
    bls_model: BoxLeastSquares,
):
    phase_hours = (
        signed_phase_days(time_days, measured["period_days"], measured["transit_time_btjd"])
        * 24
    )
    near_transit = np.abs(phase_hours) <= 36
    axis.scatter(
        phase_hours[near_transit],
        flux_ppt[near_transit],
        s=8,
        color=NAVY,
        alpha=0.16,
        rasterized=True,
        label="30-min observations",
    )
    centers, medians, _ = robust_phase_bins(
        phase_hours[near_transit], flux_ppt[near_transit]
    )
    axis.plot(
        centers,
        medians,
        "o",
        markersize=4,
        color=color,
        label="45-min robust median",
    )
    model_phase_hours = np.linspace(-36, 36, 1000)
    model_times = measured["transit_time_btjd"] + model_phase_hours / 24
    model_flux = bls_model.model(
        model_times * u.day,
        measured["period_days"] * u.day,
        measured["duration_days"] * u.day,
        measured["transit_time_btjd"] * u.day,
    )
    axis.plot(
        model_phase_hours,
        (np.asarray(model_flux, dtype=float) - 1.0) * 1000,
        color=CORAL,
        linestyle="--",
        linewidth=1.7,
        label="Fitted BLS model",
    )
    axis.axvline(0, color=NAVY, linestyle=":", linewidth=1)
    axis.set(
        title=f"Folded signal: TOI-791 {label}",
        xlabel="Hours from recovered transit center",
        ylabel="Detrended flux (ppt)",
        xlim=(-36, 36),
    )
    axis.grid(color=LIGHT, linewidth=0.7, alpha=0.8)
    axis.legend(frameon=False, loc="lower left", fontsize=8)
    if label == "c":
        axis.text(
            0.98,
            0.96,
            f"{len(c_event_windows)} TESS event windows\npublication-informed period ranking",
            transform=axis.transAxes,
            ha="right",
            va="top",
            fontsize=8,
            color=NAVY,
        )


draw_folded_panel(ax_fold_b, "b", recovered_b, TEAL, bls_b)
draw_folded_panel(ax_fold_c, "c", recovered_c, GOLD, bls_c)

figure.suptitle(
    "TOI-791 b recovery and publication-informed c period ranking",
    fontsize=18,
    fontweight="bold",
    color=NAVY,
)
figure.text(
    0.5,
    -0.01,
    "Data: NASA TESS / MAST public 2-minute SPOC PDCSAP products. "
    "Dotted lines mark peer-reviewed periods; BLS models are diagnostic, not physical transit fits.",
    ha="center",
    va="bottom",
    fontsize=9,
    color=NAVY,
)

chart_path = OUTPUT_DIR / "toi791_transit_recovery.png"
figure.savefig(chart_path, dpi=220, bbox_inches="tight", facecolor="white")
if "agg" in plt.get_backend().lower():
    plt.close(figure)
else:
    plt.show()

print(f"Saved chart: {chart_path.relative_to(PROJECT_ROOT)}")
print(f"Saved results: {results_path.relative_to(PROJECT_ROOT)}")
print(f"Saved manifest: {manifest_path.relative_to(PROJECT_ROOT)}")
print(
    "Saved TOI-791 c event inventory: "
    f"{c_events_path.relative_to(PROJECT_ROOT)}"
)


# %% [markdown]
# ## Interpretation and limitations
#
# - **TOI-791 b:** the strongest broad-search BLS peak is approximately 139.31
#   days, consistent with the paper's 139.29931-day period at the reporting
#   precision used here. This is a numerical peak location, not a formal period
#   measurement with an uncertainty.
# - **TOI-791 c:** after b is masked, a publication-informed 170–280-day search
#   ranks approximately 232.01 days highest. The selected TESS data contain two
#   c-like event windows separated by about 928 days, or four 232-day cycles.
#   TESS-only folding does not independently eliminate every integer-cycle
#   alias; the paper's additional observations resolve that ambiguity.
# - **Folded panels:** folding stacks observations separated by one recovered
#   period. For b, multiple events align. For c, the two event windows align at
#   zero phase. The dashed curve is only the fitted BLS approximation; real
#   transits have curved ingress, egress, and limb-darkened bottoms.
# - **Do not overread the depths:** sector detrending and 30-minute binning
#   attenuate them. The peer-reviewed values, not this notebook's BLS boxes,
#   should be used for physical interpretation.
# - **Not a confirmation analysis:** BLS cannot rule out eclipsing binaries,
#   nearby contaminants, or every instrumental effect. The study's ground-based
#   and dynamical follow-up is what establishes the planetary interpretation.
# - **Not a blinded discovery search:** the period ranges and the second-stage
#   masking strategy were chosen with the published system in mind. This is a
#   transparent recovery exercise, not independent alias resolution.
#
# ## Possible next steps
#
# 1. Repeat the recovery across several defensible detrending windows and bin
#    sizes to measure how sensitive the periods, depths, and rankings are to
#    preprocessing choices.
# 2. Run injection/recovery tests and odd/even transit checks, then inspect TESS
#    pixel-level centroids to quantify signal completeness and contamination
#    risk.
# 3. Fit individual transit times from later sectors and compare them with a
#    linear ephemeris before attempting a full transit-timing-variation model.