410 lines
15 KiB
Python
410 lines
15 KiB
Python
#!/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"]
|
|
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),
|
|
]
|
|
|
|
|
|
def parse_decimal_input(value: str) -> Decimal | None:
|
|
text = str(value or "").strip()
|
|
if text == "":
|
|
return None
|
|
|
|
normalized = text.replace(" ", "").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:
|
|
return dt_utc.astimezone(timezone.utc).strftime("%Y-%m-%d %H:%M:%S")
|
|
|
|
|
|
def format_german_datetime(dt_utc: datetime) -> str:
|
|
return dt_utc.astimezone(timezone.utc).strftime("%d.%m.%Y %H:%M:%S")
|
|
|
|
|
|
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 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
|
|
|
|
for fmt in ("%Y-%m-%d %H:%M:%S", "%Y-%m-%dT%H:%M:%S", "%Y-%m-%d %H:%M", "%Y-%m-%d"):
|
|
try:
|
|
parsed = datetime.strptime(text, 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 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": "",
|
|
}
|
|
|
|
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"] = "Nicht verfuegbar"
|
|
|
|
try:
|
|
tai_datetime = astropy_time.tai.to_datetime(timezone=timezone.utc)
|
|
fields["tai"] = format_german_datetime(tai_datetime)
|
|
except Exception:
|
|
fields["tai"] = "Nicht verfuegbar"
|
|
else:
|
|
leap_seconds = leap_seconds_for_utc(dt_utc)
|
|
fields["ut1"] = format_german_datetime(dt_utc) + " (nahezu UTC)"
|
|
fields["tai"] = format_german_datetime(dt_utc + timedelta(seconds=leap_seconds))
|
|
|
|
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 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)
|
|
else:
|
|
raise ValueError("Unbekannter Umrechnungsbereich.")
|
|
|
|
return {
|
|
"ok": True,
|
|
"action": "astronomical_conversions",
|
|
"kind": kind,
|
|
**result,
|
|
}
|