#!/usr/bin/env python3
"""Reproduce the article's hypothetical arithmetic; NOT a meteor forecast.

No sky/weather/ephemeris data are fetched. Assumptions: fixed activity and
population index, unobstructed view (F=1), standard perception (cp=1), and
the simple sine radiant correction (gamma=1). The domain checks below
are deliberate teaching restrictions, not a complete scientific model.
Reference: Dubietis & Arlt (2007), MNRAS 376, 890-894, equations 1 and 3.
https://articles.adsabs.harvard.edu/pdf/2007MNRAS.376..890D
"""

import json
import math
import platform


def teaching_rate(zhr, radiant_deg, limiting_mag, population_index):
    values = (zhr, radiant_deg, limiting_mag, population_index)
    if not all(math.isfinite(v) for v in values):
        raise ValueError("all parameters must be finite and known")
    if zhr < 0 or not (10 < radiant_deg <= 90):
        raise ValueError("require ZHR >= 0 and 10 < radiant altitude <= 90")
    if limiting_mag > 6.5 or population_index <= 1:
        raise ValueError("teaching domain: limiting magnitude <= 6.5 and r > 1")
    correction = population_index ** (6.5 - limiting_mag) / math.sin(
        math.radians(radiant_deg)
    )
    return zhr / correction, correction


def raw_hourly_rate(count, effective_minutes):
    if not math.isfinite(effective_minutes) or effective_minutes <= 0:
        raise ValueError("effective observing time must be positive and finite")
    if not isinstance(count, int) or count < 0:
        raise ValueError("count must be a nonnegative integer")
    return count * 60 / effective_minutes


def main():
    cases = [
        ("A", 90, 6.5),
        ("B", 30, 6.5),
        ("C", 30, 5.5),
        ("D", 30, 4.5),
        ("E", 60, 5.5),
    ]
    rows = []
    for name, altitude, lm in cases:
        rate, correction = teaching_rate(100, altitude, lm, 2.5)
        rows.append({
            "case": name,
            "assumed_zhr": 100,
            "assumed_radiant_altitude_deg": altitude,
            "assumed_limiting_magnitude": lm,
            "assumed_population_index": 2.5,
            "model_hourly_rate": round(rate, 6),
            "model_expected_count_in_30_min": round(rate / 2, 6),
            "inverse_correction_factor": round(correction, 6),
            "within_paper_correction_cutoff_5": correction <= 5 + 1e-12,
        })

    assert math.isclose(teaching_rate(100, 90, 6.5, 2.5)[0], 100)
    assert math.isclose(teaching_rate(100, 30, 6.5, 2.5)[0], 50)
    assert math.isclose(teaching_rate(100, 30, 5.5, 2.5)[0], 20)
    assert math.isclose(teaching_rate(100, 30, 4.5, 2.5)[0], 8)
    assert math.isclose(raw_hourly_rate(4, 30 - 5), 9.6)
    assert teaching_rate(0, 30, 5.5, 2.5)[0] == 0

    invalid = [
        ("radiant_at_cutoff_10_deg", lambda: teaching_rate(100, 10, 6.5, 2.5)),
        ("radiant_at_horizon", lambda: teaching_rate(100, 0, 6.5, 2.5)),
        ("unknown_population_index", lambda: teaching_rate(100, 30, 5.5, float("nan"))),
        ("zero_effective_observing_time", lambda: raw_hourly_rate(4, 0)),
    ]
    rejected = []
    for name, operation in invalid:
        try:
            operation()
        except ValueError:
            rejected.append(name)
        else:
            raise AssertionError("invalid input was accepted: " + name)

    report = {
        "kind": "hypothetical_teaching_calculations_not_observations_or_forecast",
        "environment": {
            "python": platform.python_version(),
            "system": platform.system(),
            "macos": platform.mac_ver()[0],
            "machine": platform.machine(),
            "dependencies": "Python standard library only",
        },
        "assumptions": {
            "field_obstruction_factor_F": 1,
            "perception_correction_cp": 1,
            "radiant_exponent_gamma": 1,
            "activity_and_conditions": "constant over each hypothetical interval",
            "sky_or_weather_data": "none",
        },
        "cases": rows,
        "hypothetical_log_example": {
            "elapsed_minutes": 30,
            "pause_minutes": 5,
            "effective_minutes": 25,
            "unclassified_meteor_count": 4,
            "all_meteor_raw_hourly_rate_not_zhr": raw_hourly_rate(4, 25),
        },
        "boundary_checks": {
            "zero_assumed_zhr_returns_zero_model_rate": True,
            "rejected_inputs": rejected,
        },
    }
    print(json.dumps(report, ensure_ascii=False, indent=2))


if __name__ == "__main__":
    main()