#!/usr/bin/env python3
import math
from datetime import datetime, timedelta, timezone

import astronomy

GAUSSIAN_GRAVITATIONAL_CONSTANT = 0.01720209895
J2000_OBLIQUITY_DEG = 23.439279444444445
PARABOLIC_ECCENTRICITY_TOLERANCE = 1.0e-6


def parse_float(value) -> float | None:
    if value is None:
        return None

    text = str(value).strip()
    if text == "":
        return None

    text = text.replace(",", ".")
    try:
        return float(text)
    except ValueError:
        return None


def parse_int(value) -> int | None:
    if value is None:
        return None

    text = str(value).strip()
    if text == "":
        return None

    try:
        return int(text)
    except ValueError:
        return None


def dt_to_time(dt_utc: datetime) -> astronomy.Time:
    dt_utc = dt_utc.astimezone(timezone.utc)
    return astronomy.Time.Make(
        dt_utc.year,
        dt_utc.month,
        dt_utc.day,
        dt_utc.hour,
        dt_utc.minute,
        dt_utc.second + (dt_utc.microsecond / 1_000_000.0),
    )


def build_perihelion_datetime(comet: dict) -> datetime | None:
    year = parse_int(comet.get("year_of_perihelion"))
    month = parse_int(comet.get("month_of_perihelion"))
    day_value = parse_float(comet.get("day_of_perihelion"))

    if year is None or month is None or day_value is None:
        return None

    if month < 1 or month > 12 or day_value <= 0.0:
        return None

    day = int(math.floor(day_value))
    if day < 1 or day > 31:
        return None

    fractional_day = day_value - day
    seconds = int(round(fractional_day * 86400.0))

    try:
        perihelion_dt = datetime(year, month, day, 0, 0, 0, tzinfo=timezone.utc)
    except ValueError:
        return None

    return perihelion_dt + timedelta(seconds=seconds)


def solve_elliptic_anomaly(mean_anomaly: float, eccentricity: float) -> float:
    anomaly = mean_anomaly if eccentricity < 0.8 else (math.pi if mean_anomaly >= 0.0 else -math.pi)
    for _ in range(30):
        delta = (anomaly - eccentricity * math.sin(anomaly) - mean_anomaly) / (1.0 - eccentricity * math.cos(anomaly))
        anomaly -= delta
        if abs(delta) < 1.0e-12:
            break
    return anomaly


def solve_hyperbolic_anomaly(mean_anomaly: float, eccentricity: float) -> float:
    if mean_anomaly == 0.0:
        anomaly = 0.0
    else:
        anomaly = math.asinh(mean_anomaly / eccentricity)

    for _ in range(40):
        sinh_value = math.sinh(anomaly)
        cosh_value = math.cosh(anomaly)
        delta = (eccentricity * sinh_value - anomaly - mean_anomaly) / (eccentricity * cosh_value - 1.0)
        anomaly -= delta
        if abs(delta) < 1.0e-12:
            break
    return anomaly


def solve_parabolic_parameter(delta_days: float, perihelion_distance_au: float) -> float:
    scale = GAUSSIAN_GRAVITATIONAL_CONSTANT * delta_days / math.sqrt(2.0 * perihelion_distance_au**3)
    parameter = scale

    for _ in range(40):
        numerator = parameter + (parameter**3) / 3.0 - scale
        denominator = 1.0 + parameter**2
        delta = numerator / denominator
        parameter -= delta
        if abs(delta) < 1.0e-12:
            break

    return parameter


def true_anomaly_and_radius(delta_days: float, perihelion_distance_au: float, eccentricity: float) -> tuple[float, float]:
    if perihelion_distance_au <= 0.0:
        raise ValueError("Periheldistanz muss positiv sein.")

    if eccentricity < 1.0 - PARABOLIC_ECCENTRICITY_TOLERANCE:
        semi_major_axis = perihelion_distance_au / (1.0 - eccentricity)
        mean_motion = GAUSSIAN_GRAVITATIONAL_CONSTANT / (semi_major_axis ** 1.5)
        mean_anomaly = math.fmod(mean_motion * delta_days, 2.0 * math.pi)
        eccentric_anomaly = solve_elliptic_anomaly(mean_anomaly, eccentricity)
        radius = semi_major_axis * (1.0 - eccentricity * math.cos(eccentric_anomaly))
        true_anomaly = 2.0 * math.atan2(
            math.sqrt(1.0 + eccentricity) * math.sin(eccentric_anomaly / 2.0),
            math.sqrt(1.0 - eccentricity) * math.cos(eccentric_anomaly / 2.0),
        )
        return true_anomaly, radius

    if eccentricity > 1.0 + PARABOLIC_ECCENTRICITY_TOLERANCE:
        semi_major_axis_abs = perihelion_distance_au / (eccentricity - 1.0)
        mean_anomaly = GAUSSIAN_GRAVITATIONAL_CONSTANT * delta_days / (semi_major_axis_abs ** 1.5)
        hyperbolic_anomaly = solve_hyperbolic_anomaly(mean_anomaly, eccentricity)
        radius = semi_major_axis_abs * (eccentricity * math.cosh(hyperbolic_anomaly) - 1.0)
        true_anomaly = 2.0 * math.atan2(
            math.sqrt(eccentricity + 1.0) * math.sinh(hyperbolic_anomaly / 2.0),
            math.sqrt(eccentricity - 1.0) * math.cosh(hyperbolic_anomaly / 2.0),
        )
        return true_anomaly, radius

    parabolic_parameter = solve_parabolic_parameter(delta_days, perihelion_distance_au)
    true_anomaly = 2.0 * math.atan(parabolic_parameter)
    radius = perihelion_distance_au * (1.0 + parabolic_parameter**2)
    return true_anomaly, radius


def ecliptic_to_equatorial(x_ecl: float, y_ecl: float, z_ecl: float) -> tuple[float, float, float]:
    epsilon = math.radians(J2000_OBLIQUITY_DEG)
    cos_epsilon = math.cos(epsilon)
    sin_epsilon = math.sin(epsilon)
    return (
        x_ecl,
        y_ecl * cos_epsilon - z_ecl * sin_epsilon,
        y_ecl * sin_epsilon + z_ecl * cos_epsilon,
    )


def comet_heliocentric_vector(comet: dict, dt_utc: datetime) -> tuple[float, float, float] | None:
    perihelion_distance_au = parse_float(comet.get("perihelion_dist_au"))
    eccentricity = parse_float(comet.get("eccentricity"))
    arg_perihelion_deg = parse_float(comet.get("arg_perihelion_deg"))
    ascending_node_deg = parse_float(comet.get("ascending_node_deg"))
    inclination_deg = parse_float(comet.get("inclination_deg"))
    perihelion_dt = build_perihelion_datetime(comet)

    if None in (
        perihelion_distance_au,
        eccentricity,
        arg_perihelion_deg,
        ascending_node_deg,
        inclination_deg,
        perihelion_dt,
    ):
        return None

    delta_days = (dt_utc - perihelion_dt).total_seconds() / 86400.0
    true_anomaly, radius = true_anomaly_and_radius(delta_days, perihelion_distance_au, eccentricity)

    arg_perihelion = math.radians(arg_perihelion_deg)
    ascending_node = math.radians(ascending_node_deg)
    inclination = math.radians(inclination_deg)
    argument_of_latitude = arg_perihelion + true_anomaly

    cos_node = math.cos(ascending_node)
    sin_node = math.sin(ascending_node)
    cos_inclination = math.cos(inclination)
    sin_inclination = math.sin(inclination)
    cos_argument = math.cos(argument_of_latitude)
    sin_argument = math.sin(argument_of_latitude)

    x_ecl = radius * (cos_node * cos_argument - sin_node * sin_argument * cos_inclination)
    y_ecl = radius * (sin_node * cos_argument + cos_node * sin_argument * cos_inclination)
    z_ecl = radius * (sin_argument * sin_inclination)

    return ecliptic_to_equatorial(x_ecl, y_ecl, z_ecl)


def estimate_magnitude(absolute_magnitude_h: float | None, slope_parameter_g: float | None, heliocentric_distance_au: float, geocentric_distance_au: float) -> float | None:
    if absolute_magnitude_h is None or slope_parameter_g is None:
        return None

    if heliocentric_distance_au <= 0.0 or geocentric_distance_au <= 0.0:
        return None

    return absolute_magnitude_h + (5.0 * math.log10(geocentric_distance_au)) + (2.5 * slope_parameter_g * math.log10(heliocentric_distance_au))


def equatorial_coordinates_from_vector(x: float, y: float, z: float) -> tuple[float, float]:
    distance = math.sqrt(x**2 + y**2 + z**2)
    if distance <= 0.0:
        raise ValueError("Geozentrischer Vektor ist null.")

    ra_hours = math.degrees(math.atan2(y, x)) / 15.0
    if ra_hours < 0.0:
        ra_hours += 24.0

    dec_deg = math.degrees(math.asin(z / distance))
    return ra_hours, dec_deg


def calculate_brightness(comet: dict, dt_utc: datetime) -> dict:
    comet_id = parse_int(comet.get("id"))
    result = {
        "id": comet_id,
        "heliocentric_distance_au": None,
        "geocentric_distance_au": None,
        "ra_hours": None,
        "dec_deg": None,
        "estimated_magnitude": None,
        "model": "stellarium_like",
    }

    heliocentric_position = comet_heliocentric_vector(comet, dt_utc)
    if heliocentric_position is None:
        result["error"] = "Bahnelemente unvollstaendig."
        return result

    comet_x, comet_y, comet_z = heliocentric_position
    heliocentric_distance_au = math.sqrt(comet_x**2 + comet_y**2 + comet_z**2)
    earth_vector = astronomy.HelioVector(astronomy.Body.Earth, dt_to_time(dt_utc))
    geo_x = comet_x - earth_vector.x
    geo_y = comet_y - earth_vector.y
    geo_z = comet_z - earth_vector.z
    geocentric_distance_au = math.sqrt(geo_x**2 + geo_y**2 + geo_z**2)
    ra_hours, dec_deg = equatorial_coordinates_from_vector(geo_x, geo_y, geo_z)

    absolute_magnitude_h = parse_float(comet.get("absolute_magnitude_h"))
    slope_parameter_g = parse_float(comet.get("slope_parameter_g"))
    estimated_magnitude = estimate_magnitude(
        absolute_magnitude_h,
        slope_parameter_g,
        heliocentric_distance_au,
        geocentric_distance_au,
    )

    result["heliocentric_distance_au"] = heliocentric_distance_au
    result["geocentric_distance_au"] = geocentric_distance_au
    result["ra_hours"] = ra_hours
    result["dec_deg"] = dec_deg
    result["estimated_magnitude"] = estimated_magnitude
    return result


def handle_request(payload: dict) -> dict:
    date_text = str(payload.get("date") or "").strip()
    if date_text == "":
        dt_utc = datetime.now(timezone.utc)
    else:
        normalized = date_text.replace("Z", "+00:00")
        dt_utc = datetime.fromisoformat(normalized)
        if dt_utc.tzinfo is None:
            dt_utc = dt_utc.replace(tzinfo=timezone.utc)
        else:
            dt_utc = dt_utc.astimezone(timezone.utc)

    comet_items = payload.get("comets")
    if not isinstance(comet_items, list):
        raise ValueError("Kometenliste fehlt oder ist ungueltig.")

    results = [calculate_brightness(comet, dt_utc) for comet in comet_items]
    return {
        "ok": True,
        "action": "comet_brightnesses",
        "date_utc": dt_utc.isoformat().replace("+00:00", "Z"),
        "results": results,
    }