#!/usr/bin/env python3
import math
from datetime import datetime, timedelta, timezone
from decimal import Decimal, InvalidOperation, ROUND_HALF_UP, localcontext

import astronomy
try:
    from astropy.time import Time as AstropyTime
    ASTROPY_AVAILABLE = True
except ImportError:
    AstropyTime = None
    ASTROPY_AVAILABLE = False

SPEED_OF_LIGHT_KM_S = 299792.458
H0 = 70.0
OMEGA_M = 0.3
OMEGA_LAMBDA = 0.7
MPC_IN_KM = 3.0856775814913673e19
SECONDS_PER_GYR = 365.25 * 24 * 3600 * 1e9

DECIMAL_PI = Decimal("3.14159265358979323846264338327950288419716939937510")
DECIMAL_180 = Decimal("180")
DECIMAL_36525 = Decimal("36525")
DECIMAL_2400000_5 = Decimal("2400000.5")
DECIMAL_2440587_5 = Decimal("2440587.5")
DECIMAL_86400 = Decimal("86400")
DECIMAL_1000 = Decimal("1000")
DECIMAL_ONE = Decimal("1")
DECIMAL_HUNDRED = Decimal("100")

DISTANCE_TO_KM = {
    "km": Decimal("1"),
    "au": Decimal("149597870.7"),
    "lm": Decimal("17987547.48"),
    "lh": Decimal("1079252848.8"),
    "ly": Decimal("9460730472580.8"),
    "pc": Decimal("30856775814913.67"),
    "mpc": Decimal("30856775814913670000"),
}

DISTANCE_ORDER = ["km", "au", "lm", "lh", "ly", "pc", "mpc"]
REDSHIFT_ORDER = [
    "z",
    "restNm",
    "factor",
    "scaleFactor",
    "observedNm",
    "frequencyFactor",
    "velocityApprox",
    "velocityRel",
    "velocityC",
    "lookbackGyr",
    "comovingMpc",
    "luminosityMpc",
    "angularDiameterMpc",
]
TIME_ORDER = ["utc", "ut1", "tai", "jd", "mjd", "unix", "trueSolar", "sidereal"]
COORDINATE_ORDER = ["ra", "dec", "galLon", "galLat"]
ANGLE_ORDER = ["deg", "arcmin", "arcsec", "rad"]
ANGLE_TO_DEGREES = {
    "deg": Decimal("1"),
    "arcmin": Decimal("0.01666666666666666666666666667"),
    "arcsec": Decimal("0.0002777777777777777777777777778"),
    "rad": DECIMAL_180 / DECIMAL_PI,
}
LEAP_SECONDS = [
    ("1972-01-01T00:00:00Z", 10),
    ("1972-07-01T00:00:00Z", 11),
    ("1973-01-01T00:00:00Z", 12),
    ("1974-01-01T00:00:00Z", 13),
    ("1975-01-01T00:00:00Z", 14),
    ("1976-01-01T00:00:00Z", 15),
    ("1977-01-01T00:00:00Z", 16),
    ("1978-01-01T00:00:00Z", 17),
    ("1979-01-01T00:00:00Z", 18),
    ("1980-01-01T00:00:00Z", 19),
    ("1981-07-01T00:00:00Z", 20),
    ("1982-07-01T00:00:00Z", 21),
    ("1983-07-01T00:00:00Z", 22),
    ("1985-07-01T00:00:00Z", 23),
    ("1988-01-01T00:00:00Z", 24),
    ("1990-01-01T00:00:00Z", 25),
    ("1991-01-01T00:00:00Z", 26),
    ("1992-07-01T00:00:00Z", 27),
    ("1993-07-01T00:00:00Z", 28),
    ("1994-07-01T00:00:00Z", 29),
    ("1996-01-01T00:00:00Z", 30),
    ("1997-07-01T00:00:00Z", 31),
    ("1999-01-01T00:00:00Z", 32),
    ("2006-01-01T00:00:00Z", 33),
    ("2009-01-01T00:00:00Z", 34),
    ("2012-07-01T00:00:00Z", 35),
    ("2015-07-01T00:00:00Z", 36),
    ("2017-01-01T00:00:00Z", 37),
]

EQ_TO_GAL_MATRIX = (
    (-0.0548755604162154, -0.8734370902348850, -0.4838350155487132),
    (0.4941094278755837, -0.4448296299600112, 0.7469822444972189),
    (-0.8676661490190047, -0.1980763734312015, 0.4559837761750669),
)


def parse_decimal_input(value: str) -> Decimal | None:
    text = str(value or "").strip()
    if text == "":
        return None

    normalized = text.replace(" ", "")
    if "," in normalized and "." in normalized:
        if normalized.rfind(",") > normalized.rfind("."):
            normalized = normalized.replace(".", "").replace(",", ".")
        else:
            normalized = normalized.replace(",", "")
    elif "," in normalized:
        normalized = normalized.replace(".", "").replace(",", ".")
    if normalized.count(".") > 1:
        return None

    try:
        return Decimal(normalized)
    except InvalidOperation:
        return None


def format_decimal_german(value: Decimal | float | int | None, *, max_fraction_digits: int = 10, min_fraction_digits: int = 0) -> str:
    if value is None:
        return ""

    decimal_value = value if isinstance(value, Decimal) else Decimal(str(value))
    quantize_digits = max(0, int(max_fraction_digits))

    with localcontext() as context:
        context.prec = max(50, quantize_digits + 20)
        if quantize_digits > 0:
            quantum = Decimal("1").scaleb(-quantize_digits)
            decimal_value = decimal_value.quantize(quantum, rounding=ROUND_HALF_UP)

    sign = "-" if decimal_value < 0 else ""
    decimal_value = abs(decimal_value)
    text = format(decimal_value, "f")
    integer_part, _, fractional_part = text.partition(".")

    fractional_part = fractional_part.rstrip("0")
    if len(fractional_part) < min_fraction_digits:
        fractional_part = fractional_part + ("0" * (min_fraction_digits - len(fractional_part)))

    integer_part = format(int(integer_part or "0"), ",").replace(",", ".")
    return f"{sign}{integer_part},{fractional_part}" if fractional_part else f"{sign}{integer_part}"


def format_utc_datetime(dt_utc: datetime) -> str:
    dt_utc = dt_utc.astimezone(timezone.utc)
    return dt_utc.strftime("%Y-%m-%d %H:%M:%S") + f".{dt_utc.microsecond // 1000:03d}"


def format_german_datetime(dt_utc: datetime) -> str:
    dt_utc = dt_utc.astimezone(timezone.utc)
    return dt_utc.strftime("%d.%m.%Y %H:%M:%S") + f",{dt_utc.microsecond // 1000:03d}"


def format_clock(hours_value: float) -> str:
    normalized = (hours_value % 24.0 + 24.0) % 24.0
    total_seconds = int(round(normalized * 3600.0))
    hours = (total_seconds // 3600) % 24
    minutes = (total_seconds % 3600) // 60
    seconds = total_seconds % 60
    return f"{hours:02d}:{minutes:02d}:{seconds:02d}"


def normalize_angle_degrees(value: float) -> float:
    normalized = math.fmod(value, 360.0)
    if normalized < 0.0:
        normalized += 360.0
    return normalized


def transpose_matrix(matrix: tuple[tuple[float, float, float], ...]) -> tuple[tuple[float, float, float], ...]:
    return tuple(tuple(matrix[row][column] for row in range(3)) for column in range(3))


def multiply_matrix_vector(
    matrix: tuple[tuple[float, float, float], ...], vector: tuple[float, float, float]
) -> tuple[float, float, float]:
    return tuple(sum(matrix_row[index] * vector[index] for index in range(3)) for matrix_row in matrix)


GAL_TO_EQ_MATRIX = transpose_matrix(EQ_TO_GAL_MATRIX)


def leap_seconds_for_utc(dt_utc: datetime) -> int:
    count = 0
    for effective_iso, total_offset in LEAP_SECONDS:
        effective_dt = datetime.strptime(effective_iso, "%Y-%m-%dT%H:%M:%SZ").replace(tzinfo=timezone.utc)
        if dt_utc >= effective_dt:
            count = total_offset
        else:
            break
    return count


def parse_utc_datetime(value: str) -> datetime | None:
    text = str(value or "").strip()
    if text == "":
        return None

    normalized = text.replace(",", ".")
    for fmt in ("%Y-%m-%d %H:%M:%S.%f", "%Y-%m-%dT%H:%M:%S.%f", "%Y-%m-%d %H:%M:%S", "%Y-%m-%dT%H:%M:%S", "%Y-%m-%d %H:%M", "%Y-%m-%d"):
        try:
            parsed = datetime.strptime(normalized, fmt)
            return parsed.replace(tzinfo=timezone.utc)
        except ValueError:
            continue

    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 time_to_datetime(time_value: astronomy.Time) -> datetime:
    year, month, day, hour, minute, second = time_value.Calendar()
    second_int = int(second)
    microsecond = int(round((second - second_int) * 1_000_000))

    if microsecond >= 1_000_000:
        second_int += 1
        microsecond -= 1_000_000

    return datetime(year, month, day, hour, minute, second_int, microsecond, tzinfo=timezone.utc)


def julian_decimal_to_datetime(jd_decimal: Decimal) -> datetime:
    ut = float(jd_decimal - Decimal("2451545.0"))
    return time_to_datetime(astronomy.Time(ut))


def datetime_to_julian_decimal(dt_utc: datetime) -> Decimal:
    time_value = dt_to_time(dt_utc)
    return Decimal(str(time_value.ut)) + Decimal("2451545.0")


def cosmology_e(current_z: float) -> float:
    return math.sqrt((OMEGA_M * ((1 + current_z) ** 3)) + OMEGA_LAMBDA)


def adaptive_simpson(fn, start: float, end: float, epsilon: float = 1e-10, max_depth: int = 20) -> float:
    def simpson(left: float, right: float, left_value: float, middle_value: float, right_value: float) -> float:
        return (right - left) * (left_value + 4 * middle_value + right_value) / 6

    def recurse(left: float, right: float, eps: float, whole: float, left_value: float, middle_value: float, right_value: float, depth: int) -> float:
        middle = (left + right) / 2
        left_middle = (left + middle) / 2
        right_middle = (middle + right) / 2
        left_middle_value = fn(left_middle)
        right_middle_value = fn(right_middle)
        left_area = simpson(left, middle, left_value, left_middle_value, middle_value)
        right_area = simpson(middle, right, middle_value, right_middle_value, right_value)
        delta = left_area + right_area - whole

        if depth <= 0 or abs(delta) <= 15 * eps:
            return left_area + right_area + (delta / 15)

        return recurse(left, middle, eps / 2, left_area, left_value, left_middle_value, middle_value, depth - 1) + recurse(
            middle, right, eps / 2, right_area, middle_value, right_middle_value, right_value, depth - 1
        )

    left_value = fn(start)
    right_value = fn(end)
    middle = (start + end) / 2
    middle_value = fn(middle)
    whole = simpson(start, end, left_value, middle_value, right_value)
    if not math.isfinite(whole):
        return math.nan

    return recurse(start, end, epsilon, whole, left_value, middle_value, right_value, max_depth)


def spherical_to_cartesian(longitude_deg: float, latitude_deg: float) -> tuple[float, float, float]:
    longitude_rad = math.radians(longitude_deg)
    latitude_rad = math.radians(latitude_deg)
    cos_latitude = math.cos(latitude_rad)
    return (
        cos_latitude * math.cos(longitude_rad),
        cos_latitude * math.sin(longitude_rad),
        math.sin(latitude_rad),
    )


def cartesian_to_spherical(x: float, y: float, z: float) -> tuple[float, float]:
    radius = math.sqrt((x * x) + (y * y) + (z * z))
    if radius == 0.0:
        raise ValueError("Koordinatenvektor darf nicht null sein.")

    longitude_deg = normalize_angle_degrees(math.degrees(math.atan2(y, x)))
    latitude_deg = math.degrees(math.asin(max(-1.0, min(1.0, z / radius))))
    return longitude_deg, latitude_deg


def convert_distance(payload: dict) -> dict:
    source_unit = str(payload.get("sourceUnit") or "").strip().lower()
    source_value = parse_decimal_input(payload.get("sourceValue", ""))

    if source_unit not in DISTANCE_TO_KM or source_value is None:
        return {"valid": False, "fields": {key: "" for key in DISTANCE_ORDER}}

    value_in_km = source_value * DISTANCE_TO_KM[source_unit]
    fields: dict[str, str] = {}
    for unit in DISTANCE_ORDER:
        converted = value_in_km / DISTANCE_TO_KM[unit]
        fields[unit] = format_decimal_german(converted)

    return {"valid": True, "fields": fields}


def convert_redshift(payload: dict) -> dict:
    z_decimal = parse_decimal_input(payload.get("z", ""))
    rest_nm_decimal = parse_decimal_input(payload.get("restNm", ""))
    if z_decimal is None or rest_nm_decimal is None:
        return {"valid": False, "fields": {key: "" for key in REDSHIFT_ORDER}}

    z = float(z_decimal)
    rest_nm = float(rest_nm_decimal)
    if not math.isfinite(z) or not math.isfinite(rest_nm) or z <= -1.0:
        return {"valid": False, "fields": {key: "" for key in REDSHIFT_ORDER}}

    factor_decimal = DECIMAL_ONE + z_decimal
    scale_factor_decimal = DECIMAL_ONE / factor_decimal
    observed_nm_decimal = rest_nm_decimal * factor_decimal
    frequency_factor_decimal = DECIMAL_ONE / factor_decimal
    velocity_approx_decimal = z_decimal * Decimal(str(SPEED_OF_LIGHT_KM_S))
    factor_squared_decimal = factor_decimal * factor_decimal
    beta_rel_decimal = (factor_squared_decimal - DECIMAL_ONE) / (factor_squared_decimal + DECIMAL_ONE)
    velocity_rel_decimal = beta_rel_decimal * Decimal(str(SPEED_OF_LIGHT_KM_S))
    velocity_c_decimal = beta_rel_decimal * DECIMAL_HUNDRED

    factor = float(factor_decimal)
    comoving_integral = adaptive_simpson(lambda current_z: 1 / cosmology_e(current_z), 0, z, 1e-10, 20)
    lookback_integral = adaptive_simpson(lambda current_z: 1 / ((1 + current_z) * cosmology_e(current_z)), 0, z, 1e-10, 20)
    hubble_time_seconds = MPC_IN_KM / H0
    comoving_mpc = (SPEED_OF_LIGHT_KM_S / H0) * comoving_integral
    luminosity_mpc = comoving_mpc * factor
    angular_diameter_mpc = comoving_mpc / factor
    lookback_gyr = (hubble_time_seconds * lookback_integral) / SECONDS_PER_GYR

    fields = {
        "z": format_decimal_german(z_decimal),
        "restNm": format_decimal_german(rest_nm_decimal),
        "factor": format_decimal_german(factor_decimal),
        "scaleFactor": format_decimal_german(scale_factor_decimal),
        "observedNm": format_decimal_german(observed_nm_decimal),
        "frequencyFactor": format_decimal_german(frequency_factor_decimal),
        "velocityApprox": format_decimal_german(velocity_approx_decimal),
        "velocityRel": format_decimal_german(velocity_rel_decimal),
        "velocityC": format_decimal_german(velocity_c_decimal),
        "lookbackGyr": format_decimal_german(lookback_gyr),
        "comovingMpc": format_decimal_german(comoving_mpc),
        "luminosityMpc": format_decimal_german(luminosity_mpc),
        "angularDiameterMpc": format_decimal_german(angular_diameter_mpc),
    }
    return {"valid": True, "fields": fields}


def convert_time(payload: dict) -> dict:
    source_unit = str(payload.get("sourceUnit") or "").strip()
    longitude = payload.get("longitude")
    longitude_value = None
    try:
        longitude_value = float(longitude) if longitude not in (None, "") else None
    except (TypeError, ValueError):
        longitude_value = None

    dt_utc: datetime | None = None
    if source_unit == "utc":
        dt_utc = parse_utc_datetime(payload.get("sourceValue", ""))
    elif source_unit == "jd":
        jd_decimal = parse_decimal_input(payload.get("sourceValue", ""))
        if jd_decimal is not None:
            dt_utc = julian_decimal_to_datetime(jd_decimal)
    elif source_unit == "mjd":
        mjd_decimal = parse_decimal_input(payload.get("sourceValue", ""))
        if mjd_decimal is not None:
            dt_utc = julian_decimal_to_datetime(mjd_decimal + DECIMAL_2400000_5)
    elif source_unit == "unix":
        unix_decimal = parse_decimal_input(payload.get("sourceValue", ""))
        if unix_decimal is not None:
            unix_ms = int((unix_decimal * DECIMAL_1000).to_integral_value(rounding=ROUND_HALF_UP))
            dt_utc = datetime.fromtimestamp(unix_ms / 1000, tz=timezone.utc)

    if dt_utc is None:
        return {"valid": False, "fields": {key: "" for key in TIME_ORDER}}

    jd_decimal = datetime_to_julian_decimal(dt_utc)
    mjd_decimal = jd_decimal - DECIMAL_2400000_5
    unix_decimal = Decimal(str(dt_utc.timestamp()))

    fields = {
        "utc": format_utc_datetime(dt_utc),
        "ut1": "",
        "tai": "",
        "jd": format_decimal_german(jd_decimal, min_fraction_digits=1),
        "mjd": format_decimal_german(mjd_decimal),
        "unix": format_decimal_german(unix_decimal),
        "trueSolar": "",
        "sidereal": "",
    }

    leap_seconds = leap_seconds_for_utc(dt_utc)
    ut1_fallback = format_german_datetime(dt_utc) + " (nahezu UTC)"
    tai_fallback = format_german_datetime(dt_utc + timedelta(seconds=leap_seconds))

    if ASTROPY_AVAILABLE:
        astropy_time = AstropyTime(dt_utc, scale="utc")
        try:
            ut1_datetime = astropy_time.ut1.to_datetime(timezone=timezone.utc)
            fields["ut1"] = format_german_datetime(ut1_datetime)
        except Exception:
            fields["ut1"] = ut1_fallback

        try:
            tai_datetime = astropy_time.tai.to_datetime(timezone=timezone.utc)
            fields["tai"] = format_german_datetime(tai_datetime)
        except Exception:
            fields["tai"] = tai_fallback
    else:
        fields["ut1"] = ut1_fallback
        fields["tai"] = tai_fallback

    if longitude_value is not None and math.isfinite(longitude_value):
        observer = astronomy.Observer(0.0, longitude_value, 0.0)
        time_value = dt_to_time(dt_utc)
        sidereal_hours = (astronomy.SiderealTime(time_value) + longitude_value / 15.0) % 24.0
        hour_angle = astronomy.HourAngle(astronomy.Body.Sun, time_value, observer)
        true_solar_hours = (hour_angle + 12.0) % 24.0
        utc_hours = dt_utc.hour + (dt_utc.minute / 60.0) + (dt_utc.second / 3600.0) + (dt_utc.microsecond / 3_600_000_000.0)
        delta_hours = true_solar_hours - utc_hours
        if delta_hours <= -12.0:
            delta_hours += 24.0
        elif delta_hours > 12.0:
            delta_hours -= 24.0

        true_solar_dt = dt_utc + timedelta(hours=delta_hours)
        fields["trueSolar"] = format_german_datetime(true_solar_dt)
        fields["sidereal"] = format_clock(sidereal_hours)
    else:
        fields["trueSolar"] = "Kein Standard-Standort verfuegbar"
        fields["sidereal"] = "Kein Standard-Standort verfuegbar"

    return {"valid": True, "fields": fields}


def convert_coordinates(payload: dict) -> dict:
    source_system = str(payload.get("sourceSystem") or "").strip().lower()
    first_value = parse_decimal_input(payload.get("value1", ""))
    second_value = parse_decimal_input(payload.get("value2", ""))

    if first_value is None or second_value is None:
        return {"valid": False, "fields": {key: "" for key in COORDINATE_ORDER}}

    longitude = float(first_value)
    latitude = float(second_value)
    if not math.isfinite(longitude) or not math.isfinite(latitude) or latitude < -90.0 or latitude > 90.0:
        return {"valid": False, "fields": {key: "" for key in COORDINATE_ORDER}}

    if source_system == "equatorial":
        equatorial_vector = spherical_to_cartesian(normalize_angle_degrees(longitude), latitude)
        galactic_vector = multiply_matrix_vector(EQ_TO_GAL_MATRIX, equatorial_vector)
        ra_deg, dec_deg = cartesian_to_spherical(*equatorial_vector)
        gal_lon_deg, gal_lat_deg = cartesian_to_spherical(*galactic_vector)
    elif source_system == "galactic":
        galactic_vector = spherical_to_cartesian(normalize_angle_degrees(longitude), latitude)
        equatorial_vector = multiply_matrix_vector(GAL_TO_EQ_MATRIX, galactic_vector)
        gal_lon_deg, gal_lat_deg = cartesian_to_spherical(*galactic_vector)
        ra_deg, dec_deg = cartesian_to_spherical(*equatorial_vector)
    else:
        raise ValueError("Unbekanntes Koordinatensystem.")

    fields = {
        "ra": format_decimal_german(ra_deg, max_fraction_digits=6),
        "dec": format_decimal_german(dec_deg, max_fraction_digits=6),
        "galLon": format_decimal_german(gal_lon_deg, max_fraction_digits=6),
        "galLat": format_decimal_german(gal_lat_deg, max_fraction_digits=6),
    }
    return {"valid": True, "fields": fields}


def convert_angle(payload: dict) -> dict:
    source_unit = str(payload.get("sourceUnit") or "").strip().lower()
    source_value = parse_decimal_input(payload.get("sourceValue", ""))

    if source_unit not in ANGLE_TO_DEGREES or source_value is None:
        return {"valid": False, "fields": {key: "" for key in ANGLE_ORDER}}

    value_in_degrees = source_value * ANGLE_TO_DEGREES[source_unit]
    fields: dict[str, str] = {}
    for unit in ANGLE_ORDER:
        converted = value_in_degrees / ANGLE_TO_DEGREES[unit]
        fields[unit] = format_decimal_german(converted, max_fraction_digits=12)

    return {"valid": True, "fields": fields}


def handle_request(payload: dict) -> dict:
    kind = str(payload.get("kind") or "").strip().lower()
    if kind == "distance":
        result = convert_distance(payload)
    elif kind == "redshift":
        result = convert_redshift(payload)
    elif kind == "time":
        result = convert_time(payload)
    elif kind == "coordinates":
        result = convert_coordinates(payload)
    elif kind == "angle":
        result = convert_angle(payload)
    else:
        raise ValueError("Unbekannter Umrechnungsbereich.")

    return {
        "ok": True,
        "action": "astronomical_conversions",
        "kind": kind,
        **result,
    }