293 lines
10 KiB
Python
293 lines
10 KiB
Python
#!/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,
|
|
}
|