Kometenseite erstmal fertig
This commit is contained in:
Eskimue
Binary file not shown.
+33
-2
@@ -13,6 +13,7 @@ if SCRIPT_DIR not in sys.path:
|
||||
|
||||
import astronomy
|
||||
import astronomical_conversions
|
||||
import comets
|
||||
|
||||
MOON_RADIUS_KM = 1737.4
|
||||
SYNODIC_MONTH = 29.530588853
|
||||
@@ -252,6 +253,31 @@ def action_astronomical_conversions(args: list[str]) -> dict:
|
||||
fail(str(exc), extra={"kind": kind})
|
||||
|
||||
|
||||
def action_comet_brightnesses(args: list[str]) -> dict:
|
||||
if len(args) != 1:
|
||||
fail(
|
||||
"Aktion comet_brightnesses erwartet 1 Argument: payloadBase64 oder @payloadDatei",
|
||||
extra={"argv": args},
|
||||
)
|
||||
|
||||
try:
|
||||
payload_arg = args[0]
|
||||
if payload_arg.startswith("@"):
|
||||
payload_path = payload_arg[1:]
|
||||
with open(payload_path, "r", encoding="utf-8") as handle:
|
||||
payload_text = handle.read()
|
||||
else:
|
||||
payload_text = base64.urlsafe_b64decode(payload_arg.encode("ascii")).decode("utf-8")
|
||||
payload = json.loads(payload_text)
|
||||
except Exception as exc:
|
||||
fail("Payload fuer comet_brightnesses ist ungueltig.", extra={"details": str(exc)})
|
||||
|
||||
try:
|
||||
return comets.handle_request(payload)
|
||||
except ValueError as exc:
|
||||
fail(str(exc), extra={"action": "comet_brightnesses"})
|
||||
|
||||
|
||||
def search_event(
|
||||
body: astronomy.Body,
|
||||
direction: astronomy.Direction,
|
||||
@@ -4173,7 +4199,7 @@ def action_satellite_passes(args: list[str]) -> dict:
|
||||
|
||||
def main() -> None:
|
||||
if len(sys.argv) < 2:
|
||||
fail("Es wurde keine Aktion uebergeben.", extra={"available_actions": ["sun_moon_rise_set", "solar_longitude_to_datetime", "current_solar_longitude", "moon_star_occultations", "moon_phase_details", "satellite_passes", "planet_rise_set", "month_sky_context", "moon_phases_for_month", "golden_handle_for_month", "season_changes_for_month", "time_changes_for_month", "moon_planet_approaches", "moon_deep_sky_approaches_for_month", "planet_bright_star_approaches_for_month", "planet_conjunctions_for_month", "eclipses_for_month", "moon_apsides_for_month", "sun_apsides_for_month", "inner_planet_elongations_for_month", "venus_peak_magnitude_for_month", "outer_planet_events_for_month", "outer_planet_stations_for_month", "jupiter_moons_one_side_for_month", "mercury_good_visibility_for_month", "planet_parades_for_month", "planet_constellation_changes_for_month", "moon_star_occultations_for_month"]})
|
||||
fail("Es wurde keine Aktion uebergeben.", extra={"available_actions": ["sun_moon_rise_set", "solar_longitude_to_datetime", "current_solar_longitude", "astronomical_conversions", "comet_brightnesses", "moon_star_occultations", "moon_phase_details", "satellite_passes", "planet_rise_set", "month_sky_context", "moon_phases_for_month", "golden_handle_for_month", "season_changes_for_month", "time_changes_for_month", "moon_planet_approaches", "moon_deep_sky_approaches_for_month", "planet_bright_star_approaches_for_month", "planet_conjunctions_for_month", "eclipses_for_month", "moon_apsides_for_month", "sun_apsides_for_month", "inner_planet_elongations_for_month", "venus_peak_magnitude_for_month", "outer_planet_events_for_month", "outer_planet_stations_for_month", "jupiter_moons_one_side_for_month", "mercury_good_visibility_for_month", "planet_parades_for_month", "planet_constellation_changes_for_month", "moon_star_occultations_for_month"]})
|
||||
|
||||
action = sys.argv[1]
|
||||
args = sys.argv[2:]
|
||||
@@ -4198,6 +4224,11 @@ def main() -> None:
|
||||
print(json.dumps(result, ensure_ascii=True))
|
||||
return
|
||||
|
||||
if action == "comet_brightnesses":
|
||||
result = action_comet_brightnesses(args)
|
||||
print(json.dumps(result, ensure_ascii=True))
|
||||
return
|
||||
|
||||
if action == "moon_star_occultations":
|
||||
result = action_moon_star_occultations(args)
|
||||
print(json.dumps(result, ensure_ascii=True))
|
||||
@@ -4323,7 +4354,7 @@ def main() -> None:
|
||||
print(json.dumps(result, ensure_ascii=True))
|
||||
return
|
||||
|
||||
fail("Unbekannte Aktion.", extra={"action": action, "available_actions": ["sun_moon_rise_set", "solar_longitude_to_datetime", "current_solar_longitude", "astronomical_conversions", "moon_star_occultations", "moon_phase_details", "satellite_passes", "planet_rise_set", "month_sky_context", "moon_phases_for_month", "golden_handle_for_month", "season_changes_for_month", "time_changes_for_month", "moon_planet_approaches", "moon_deep_sky_approaches_for_month", "planet_bright_star_approaches_for_month", "planet_conjunctions_for_month", "eclipses_for_month", "moon_apsides_for_month", "sun_apsides_for_month", "inner_planet_elongations_for_month", "venus_peak_magnitude_for_month", "outer_planet_events_for_month", "outer_planet_stations_for_month", "jupiter_moons_one_side_for_month", "mercury_good_visibility_for_month", "planet_parades_for_month", "planet_constellation_changes_for_month", "moon_star_occultations_for_month"]})
|
||||
fail("Unbekannte Aktion.", extra={"action": action, "available_actions": ["sun_moon_rise_set", "solar_longitude_to_datetime", "current_solar_longitude", "astronomical_conversions", "comet_brightnesses", "moon_star_occultations", "moon_phase_details", "satellite_passes", "planet_rise_set", "month_sky_context", "moon_phases_for_month", "golden_handle_for_month", "season_changes_for_month", "time_changes_for_month", "moon_planet_approaches", "moon_deep_sky_approaches_for_month", "planet_bright_star_approaches_for_month", "planet_conjunctions_for_month", "eclipses_for_month", "moon_apsides_for_month", "sun_apsides_for_month", "inner_planet_elongations_for_month", "venus_peak_magnitude_for_month", "outer_planet_events_for_month", "outer_planet_stations_for_month", "jupiter_moons_one_side_for_month", "mercury_good_visibility_for_month", "planet_parades_for_month", "planet_constellation_changes_for_month", "moon_star_occultations_for_month"]})
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
|
||||
@@ -0,0 +1,292 @@
|
||||
#!/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,
|
||||
}
|
||||
Reference in New Issue
Block a user