#! /usr/bin/python3
from datetime import datetime, timedelta
import zoneinfo
from pathlib import Path

from astropy.coordinates import EarthLocation, SkyCoord, AltAz
from astropy.time import Time
import astropy.units as u
import numpy as np

# 1. Standort festlegen (München)
location = EarthLocation(lat=48.173333 * u.deg, lon=11.635 * u.deg, height=520 * u.m)
local_tz = zoneinfo.ZoneInfo("Europe/Berlin")

# 2. Himmelskoordinaten von Deneb (ICRS)
deneb = SkyCoord.from_name("Deneb")

def get_culmination_time(target_date: datetime.date) -> datetime:
    """Sucht den exakten Zeitpunkt der maximalen Höhe (Kulmination) für ein Datum."""
    # Startwert: Ungefähre Kulmination um 12:00 Uhr UTC als Ausgangspunkt
    start_utc = datetime(target_date.year, target_date.month, target_date.day, 12, 0, tzinfo=zoneinfo.ZoneInfo("UTC"))
    t_start = Time(start_utc)
    
    # 24-Stunden-Fenster abtasten (1-Minute-Schritte)
    times = t_start + np.linspace(0, 24, 1440) * u.hour
    altaz_frame = AltAz(obstime=times, location=location)
    altitudes = deneb.transform_to(altaz_frame).alt.deg
    
    # Ungefähren Höchststand ermitteln
    max_idx = np.argmax(altitudes)
    t_approx = times[max_idx]
    
    # Feinabstimmung: ±2 Minuten um das Maximum in Sekundenschritten
    fine_times = t_approx + np.linspace(-120, 120, 241) * u.second
    fine_altaz = AltAz(obstime=fine_times, location=location)
    fine_altitudes = deneb.transform_to(fine_altaz).alt.deg
    
    best_time_astropy = fine_times[np.argmax(fine_altitudes)]
    
    # In lokale Zeitzone (MEZ/MESZ) konvertieren
    dt_utc = best_time_astropy.to_datetime(timezone=zoneinfo.ZoneInfo("UTC"))
    return dt_utc.astimezone(local_tz)

# 3. Das Jahr 2026 durchlaufen
output_file = Path("deneb_kulmination_2027.txt")
start_day = datetime(2027, 1, 1).date()

print("Berechne Kulminationszeiten für Deneb (2027)...")

with output_file.open("w", encoding="utf-8") as f:
    f.write("# Kulminationszeiten für Deneb (Minga, 48.173333 N, 11.635 E)\n")
    f.write("# Format: YYYY-MM-DD | Lokale Uhrzeit (MEZ/MESZ) | Höhe (Grad)\n")
    f.write("# " + "="*60 + "\n")
    
    current_day = start_day
    while current_day.year == 2027:
        culm_dt = get_culmination_time(current_day)
        
        # Exakte Höhe berechnen
        t_final = Time(culm_dt)
        alt = deneb.transform_to(AltAz(obstime=t_final, location=location)).alt.deg
        
        line = f"{current_day.isoformat()} | {culm_dt.strftime('%H:%M:%S %Z')} | {alt:.2f} Grad\n"
        f.write(line)
        
        current_day += timedelta(days=1)

print(f"Fertig! Daten gespeichert in: {output_file.resolve()}")