#!/usr/bin/env python3
"""horizon.py - terrain horizon seen from the tent (Copernicus GLO-30) and the resulting 'local sunset'
(sun disappearing behind Kholat Syakhl) on 1 Feb 1959 and local sunrise on 2 Feb.
Output: out/horizon.csv (horizon elevation angle per azimuth), out/horizon.md
"""
import csv, math, os
from datetime import datetime, timedelta, timezone
import numpy as np, tifffile, ephem
from scipy.interpolate import RegularGridInterpolator
from geo_common import RAW, OUT, load_coordinates, to_xy, MLAT, MLON, REF_LAT, REF_LON

LON0, LAT0, DLON, DLAT = 59.0, 62.0, 1 / 1800, 1 / 3600
with tifffile.TiffFile(os.path.join(RAW, "Copernicus_DSM_COG_10_N61_00_E059_00_DEM.tif")) as tf:
    full = tf.pages[0].asarray()
r0, r1 = int(round((LAT0 - 61.82) / DLAT)), int(round((LAT0 - 61.68) / DLAT)) + 1
c0, c1 = int(round((59.30 - LON0) / DLON)), int(round((59.60 - LON0) / DLON)) + 1
Z = full[r0:r1, c0:c1].astype(float); del full
lats = LAT0 - np.arange(r0, r1) * DLAT; lons = LON0 + np.arange(c0, c1) * DLON
X = (lons - REF_LON) * MLON; Y = (lats - REF_LAT) * MLAT
interp = RegularGridInterpolator((Y[::-1], X), Z[::-1, :], method="linear", bounds_error=False, fill_value=np.nan)
C = load_coordinates(); tent = C["tent_dyatlovpass_map"]
tx, ty = to_xy(tent["lat"], tent["lon"]); z0 = float(interp((ty, tx))) + 1.5   # eye height
hor = {}
for az in range(0, 360):
    ux, uy = math.sin(math.radians(az)), math.cos(math.radians(az))
    best = -90.0
    for d in np.arange(30, 6000, 10):
        z = float(interp((ty + uy * d, tx + ux * d)))
        if math.isnan(z): break
        ang = math.degrees(math.atan2(z - z0 - d * d / (2 * 6371000 * 1.17), d))  # with earth curvature/refraction (k=0.13 -> 1.17 factor)
        best = max(best, ang)
    hor[az] = best
with open(os.path.join(OUT, "horizon.csv"), "w", newline="") as f:
    w = csv.writer(f); w.writerow(["azimuth_deg", "horizon_elev_deg"]); [w.writerow([a, round(hor[a], 2)]) for a in hor]

UTC5 = timezone(timedelta(hours=5))
o = ephem.Observer(); o.lat, o.lon, o.elevation = str(tent["lat"]), str(tent["lon"]), 899; o.pressure = 1010; o.temperature = -15
sun = ephem.Sun()
def sun_altaz(t_local):
    o.date = ephem.Date(t_local.astimezone(timezone.utc).replace(tzinfo=None)); sun.compute(o)
    return math.degrees(sun.alt), math.degrees(sun.az)
events = {}
for day, label in [((1959, 2, 1), "1 Feb"), ((1959, 2, 2), "2 Feb"), ((1959, 2, 26), "26 Feb"), ((1959, 3, 5), "5 Mar")]:
    t = datetime(*day, 8, 0, tzinfo=UTC5); vis = None; first_vis = None; last_vis = None
    while t < datetime(*day, 19, 0, tzinfo=UTC5):
        alt, az = sun_altaz(t); h = hor[int(round(az)) % 360]
        up = alt + 0.27 > h  # upper limb
        if up and first_vis is None: first_vis = t
        if up: last_vis = t
        t += timedelta(minutes=1)
    events[label] = (first_vis, last_vis)
md = ["# Terrain horizon at the tent and local sun visibility (generated by horizon.py)", "",
      f"Eye point: tent TL 18.10, {z0-1.5:.0f} m + 1.5 m. Horizon from Copernicus GLO-30 rays to 6 km (10 m steps), curvature+refraction k=0.13.", "",
      "| azimuth sector | max horizon elevation (deg) | direction |", "|---|---|---|"]
for a0, a1, name in [(0, 45, "N-NE (Lozva valley / Otorten side)"), (45, 90, "NE-E (toward the cedar)"), (90, 135, "E-SE (pass, outlier rock)"), (135, 180, "SE-S"),
                     (180, 225, "S-SW (Kholat Syakhl massif)"), (225, 270, "SW-W (summit at 233 deg)"), (270, 315, "W-NW (spur crest)"), (315, 360, "NW-N")]:
    md.append(f"| {a0}-{a1} | {max(hor[a] for a in range(a0, a1)):.1f} | {name} |")
md += ["", "| date | sun first clears the terrain (upper limb) | sun last visible above terrain | astronomical sunset (flat horizon) |", "|---|---|---|---|"]
for label, (fv, lv) in events.items():
    md.append(f"| {label} | {fv.strftime('%H:%M') if fv else '-'} | {lv.strftime('%H:%M') if lv else '-'} | see astro_tables.md |")
alt17, az17 = sun_altaz(datetime(1959, 2, 1, 17, 0, tzinfo=UTC5)); alt15, az15 = sun_altaz(datetime(1959, 2, 1, 15, 0, tzinfo=UTC5))
md += ["", f"At 15:00 on 1 Feb the sun stood {alt15:.1f} deg high at azimuth {az15:.0f} deg where the terrain horizon from the tent is {hor[int(round(az15))]:.1f} deg; at 17:00 it was at {alt17:.1f} deg / {az17:.0f} deg (horizon {hor[int(round(az17))]:.1f} deg).",
       "So the tent site was in the shadow of Kholat Syakhl from mid-afternoon; direct sunlight ended there roughly 2 h before the astronomical sunset (17:03), which is relevant to any argument that dates the last photographs by their exposure."]
open(os.path.join(OUT, "horizon.md"), "w").write("\n".join(md)); print("\n".join(md))