"""Strict selection of measured times from the attributed CC BY 4.0 author table."""

from pathlib import Path
import csv
import hashlib
import io

SOURCE_HASH = "8f958e773fae5a6d17aa090c5d523e13829f6b69a2f3d51ceecbf5e2d0d56af9"


def read_observations(directory):
    raw = (Path(directory) / "ajac8553t1_mrt.txt").read_bytes()
    assert hashlib.sha256(raw).hexdigest() == SOURCE_HASH
    offsets = {"Kepler-289b": 0, "Kepler-289d": 259, "Kepler-289c": 395}
    rows = []
    all_rows = []
    for line in raw.decode().splitlines():
        if not line.startswith("Kepler-289"):
            continue
        planet = line[:11].strip()
        source_number = int(line[12:15])
        observed = line[16:25].strip()
        error = line[26:32].strip()
        all_rows.append((planet, source_number, observed, error))
        if not observed or not error:
            continue
        assert float(error) > 0
        rows.append(
            {
                "planet": planet,
                "source_num": source_number,
                "planet_epoch": source_number - offsets[planet],
                "observed_reduced_bjd": float(observed),
                "sigma_days": float(error),
            }
        )
    assert len(all_rows) == 466 and len(rows) == 51
    assert next(row for row in all_rows if row[1] == 430)[2:] == ("", "")
    selected = [
        row
        for row in rows
        if row["planet"] in {"Kepler-289d", "Kepler-289c"}
        and row["observed_reduced_bjd"] < 3000
    ]
    assert len(selected) == 25
    return rows, selected


def ephemeris(rows):
    weights = [1 / row["sigma_days"] ** 2 for row in rows]
    total = sum(weights)
    mean_epoch = sum(w * row["planet_epoch"] for w, row in zip(weights, rows)) / total
    mean_time = (
        sum(w * row["observed_reduced_bjd"] for w, row in zip(weights, rows)) / total
    )
    slope = sum(
        w
        * (row["planet_epoch"] - mean_epoch)
        * (row["observed_reduced_bjd"] - mean_time)
        for w, row in zip(weights, rows)
    ) / sum(
        w * (row["planet_epoch"] - mean_epoch) ** 2 for w, row in zip(weights, rows)
    )
    return mean_time - slope * mean_epoch, slope