"""Historical ISS geometry/illumination example; not a live forecast or app test."""

import argparse
import hashlib
import importlib.metadata
import json
from pathlib import Path

import numpy as np
from skyfield.api import EarthSatellite, load, load_file, wgs84


def calculate(root):
    provenance = json.loads((root / 'ephemeris-source.json').read_text())
    ephemeris_file = root / 'de421.bsp'
    actual_hash = hashlib.sha256(ephemeris_file.read_bytes()).hexdigest()
    if actual_hash != provenance['sha256']:
        raise ValueError('Ephemeris checksum differs from the recorded input')
    ts = load.timescale(builtin=True)
    name, line1, line2 = (root / 'iss.tle').read_text().splitlines()
    satellite = EarthSatellite(line1, line2, name, ts)
    observer = wgs84.latlon(51.4769, -0.0005, elevation_m=46.0)
    eph = load_file(str(ephemeris_file))
    earth_observer = eph['earth'] + observer
    start, end = ts.utc(2018, 7, 4), ts.utc(2018, 7, 5)
    times, kinds = satellite.find_events(observer, start, end, altitude_degrees=10.0)
    events = []
    for t, kind in zip(times, kinds):
        altitude, azimuth, distance = (satellite - observer).at(t).altaz()
        sun_altitude = earth_observer.at(t).observe(eph['sun']).apparent().altaz()[0].degrees
        sunlit = bool(satellite.at(t).is_sunlit(eph))
        events.append({
            'utc': t.utc_iso(places=3),
            'event': ('rise_above_10_deg', 'culmination', 'set_below_10_deg')[int(kind)],
            'elevation_deg': round(float(altitude.degrees), 3),
            'azimuth_deg': round(float(azimuth.degrees), 3),
            'range_km': round(float(distance.km), 3),
            'sun_elevation_deg': round(float(sun_altitude), 3),
            'satellite_sunlit': sunlit,
            'passes_screen_at_this_instant': bool(sunlit and sun_altitude < -6.0 and altitude.degrees >= 9.95),
        })
    # A separate uniform grid catches transitions within a pass, rather than
    # inferring an entire pass's visibility from its single peak.
    seconds = np.arange(0, 86400, 10)
    samples = ts.utc(2018, 7, 4, 0, 0, seconds)
    elevations = (satellite - observer).at(samples).altaz()[0].degrees
    sun_elevations = earth_observer.at(samples).observe(eph['sun']).apparent().altaz()[0].degrees
    sunlit = satellite.at(samples).is_sunlit(eph)
    mask = (elevations >= 10.0) & (sun_elevations < -6.0) & sunlit
    windows = []
    indices = np.flatnonzero(mask)
    for group in np.split(indices, np.flatnonzero(np.diff(indices) > 1) + 1):
        if len(group):
            windows.append({
                'first_passing_sample_utc': samples[int(group[0])].utc_iso(),
                'last_passing_sample_utc': samples[int(group[-1])].utc_iso(),
                'passing_sample_count': int(len(group)),
            })
    peaks = [e for e in events if e['event'] == 'culmination']
    assert len(peaks) > 0 and all(e['elevation_deg'] >= 9.9 for e in peaks)
    assert all(abs(e['elevation_deg'] - 10.0) < 0.1 for e in events if e['event'] != 'culmination')
    assert list(times.tt) == sorted(times.tt)
    return {
        'scope': 'Historical independent calculation, not an observed sighting, live forecast or SatellitePlus validation',
        'versions': {p: importlib.metadata.version(p) for p in ('skyfield', 'sgp4', 'numpy', 'jplephem')},
        'observer': {'latitude_deg': 51.4769, 'longitude_deg': -0.0005, 'height_m': 46.0, 'label': 'Greenwich reference site'},
        'tle_epoch_utc': satellite.epoch.utc_iso(),
        'tle_sha256': hashlib.sha256((root / 'iss.tle').read_bytes()).hexdigest(),
        'ephemeris_sha256': actual_hash,
        'window_utc': ['2018-07-04T00:00:00Z', '2018-07-05T00:00:00Z'],
        'screen': {'minimum_elevation_deg': 10.0, 'sun_elevation_below_deg': -6.0, 'satellite_must_be_sunlit': True, 'sample_step_seconds': 10},
        'summary': {'geometric_culminations_above_10_deg': len(peaks), 'culminations_passing_screen': sum(e['passes_screen_at_this_instant'] for e in peaks), 'sampled_screen_windows': len(windows)},
        'events': events,
        'sampled_screen_windows': windows,
        'limits': ['Screening is not a brightness prediction or a guarantee of visibility.', 'Weather, phase angle, size, attitude, light pollution, eyesight and real obstructions are not modelled.', 'Sampled endpoints have 10-second resolution; events can have finer numerical timing without equivalent physical accuracy.', 'Use fresh orbital data and your real observing site for a current session.'],
    }


if __name__ == '__main__':
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument('--output', type=Path, default=Path('reproduced.json'))
    args = parser.parse_args()
    result = calculate(Path(__file__).resolve().parent)
    args.output.write_text(json.dumps(result, indent=2) + '\n')
    print(json.dumps(result['summary'], indent=2))