1428 lines
52 KiB
Python
1428 lines
52 KiB
Python
#!/usr/bin/env python3
|
|
import base64
|
|
import json
|
|
import math
|
|
import os
|
|
import sys
|
|
from datetime import datetime, timedelta, timezone
|
|
from zoneinfo import ZoneInfo
|
|
|
|
SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__))
|
|
if SCRIPT_DIR not in sys.path:
|
|
sys.path.insert(0, SCRIPT_DIR)
|
|
|
|
import astronomy
|
|
|
|
MOON_RADIUS_KM = 1737.4
|
|
SYNODIC_MONTH = 29.530588853
|
|
STAR_BODIES = [
|
|
astronomy.Body.Star1,
|
|
astronomy.Body.Star2,
|
|
astronomy.Body.Star3,
|
|
astronomy.Body.Star4,
|
|
astronomy.Body.Star5,
|
|
astronomy.Body.Star6,
|
|
astronomy.Body.Star7,
|
|
astronomy.Body.Star8,
|
|
]
|
|
PHASE_LABELS = [
|
|
"Neumond",
|
|
"Erstes Viertel / zunehmende Sichel",
|
|
"Zunehmender Halbmond",
|
|
"Zweites Viertel",
|
|
"Vollmond",
|
|
"Drittes Viertel",
|
|
"Abnehmender Halbmond",
|
|
"Letztes Viertel / abnehmende Sichel",
|
|
]
|
|
|
|
|
|
def fail(message: str, *, extra: dict | None = None, code: int = 1) -> None:
|
|
payload = {"ok": False, "error": message}
|
|
if extra:
|
|
payload.update(extra)
|
|
print(json.dumps(payload, ensure_ascii=False))
|
|
raise SystemExit(code)
|
|
|
|
|
|
def parse_float(value: str, label: str) -> float:
|
|
try:
|
|
return float(value)
|
|
except ValueError as exc:
|
|
fail(f"{label} ist ungueltig.", extra={"details": str(exc), "value": value})
|
|
|
|
|
|
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 serialize_event(label: str, event_time: astronomy.Time | None, tz: ZoneInfo) -> dict:
|
|
if event_time is None:
|
|
return {"label": label, "found": False}
|
|
|
|
utc_dt = time_to_datetime(event_time)
|
|
local_dt = utc_dt.astimezone(tz)
|
|
|
|
return {
|
|
"label": label,
|
|
"found": True,
|
|
"utc_iso": utc_dt.isoformat().replace("+00:00", "Z"),
|
|
"local_iso": local_dt.isoformat(),
|
|
"local_time": local_dt.strftime("%H:%M"),
|
|
"local_date": local_dt.strftime("%Y-%m-%d"),
|
|
}
|
|
|
|
|
|
def search_event(
|
|
body: astronomy.Body,
|
|
direction: astronomy.Direction,
|
|
observer: astronomy.Observer,
|
|
start_time: astronomy.Time,
|
|
end_utc: datetime,
|
|
) -> astronomy.Time | None:
|
|
search_limit_days = (end_utc - time_to_datetime(start_time)).total_seconds() / 86400.0
|
|
if search_limit_days <= 0:
|
|
return None
|
|
|
|
result = astronomy.SearchRiseSet(body, observer, direction, start_time, search_limit_days)
|
|
if result is None:
|
|
return None
|
|
|
|
result_dt = time_to_datetime(result)
|
|
if result_dt >= end_utc:
|
|
return None
|
|
|
|
return result
|
|
|
|
|
|
def action_sun_moon_rise_set(args: list[str]) -> dict:
|
|
if len(args) != 5:
|
|
fail(
|
|
"Aktion sun_moon_rise_set erwartet 5 Argumente: latitude longitude elevation date timezone",
|
|
extra={"argv": args},
|
|
)
|
|
|
|
latitude = parse_float(args[0], "Latitude")
|
|
longitude = parse_float(args[1], "Longitude")
|
|
elevation = parse_float(args[2], "Elevation")
|
|
date_text = args[3]
|
|
timezone_name = args[4]
|
|
|
|
try:
|
|
datetime.strptime(date_text, "%Y-%m-%d")
|
|
except ValueError:
|
|
fail("Datum muss im Format YYYY-MM-DD uebergeben werden.", extra={"date": date_text})
|
|
|
|
try:
|
|
tz = ZoneInfo(timezone_name)
|
|
except Exception as exc: # pragma: no cover
|
|
fail("Zeitzone ist ungueltig oder auf dem Server nicht verfuegbar.", extra={"details": str(exc)})
|
|
|
|
observer = astronomy.Observer(latitude, longitude, elevation)
|
|
|
|
local_start = datetime.strptime(date_text, "%Y-%m-%d").replace(tzinfo=tz)
|
|
local_end = local_start + timedelta(days=1)
|
|
utc_start = local_start.astimezone(timezone.utc)
|
|
utc_end = local_end.astimezone(timezone.utc)
|
|
|
|
start_time = dt_to_time(utc_start)
|
|
|
|
sun_rise = search_event(astronomy.Body.Sun, astronomy.Direction.Rise, observer, start_time, utc_end)
|
|
sun_set = search_event(astronomy.Body.Sun, astronomy.Direction.Set, observer, start_time, utc_end)
|
|
moon_rise = search_event(astronomy.Body.Moon, astronomy.Direction.Rise, observer, start_time, utc_end)
|
|
moon_set = search_event(astronomy.Body.Moon, astronomy.Direction.Set, observer, start_time, utc_end)
|
|
|
|
def search_altitude(direction: astronomy.Direction, altitude_deg: float) -> astronomy.Time | None:
|
|
limit_days = (utc_end - time_to_datetime(start_time)).total_seconds() / 86400.0
|
|
result = astronomy.SearchAltitude(astronomy.Body.Sun, observer, direction, start_time, limit_days, altitude_deg)
|
|
if result is None:
|
|
return None
|
|
if time_to_datetime(result) >= utc_end:
|
|
return None
|
|
return result
|
|
|
|
# Morgens: Sonne steigt auf → Direction.Rise (Höhe wird größer)
|
|
# Abends: Sonne sinkt ab → Direction.Set (Höhe wird kleiner)
|
|
astro_dawn = search_altitude(astronomy.Direction.Rise, -18.0)
|
|
nautical_dawn = search_altitude(astronomy.Direction.Rise, -12.0)
|
|
civil_dawn = search_altitude(astronomy.Direction.Rise, -6.0)
|
|
civil_dusk = search_altitude(astronomy.Direction.Set, -6.0)
|
|
nautical_dusk = search_altitude(astronomy.Direction.Set, -12.0)
|
|
astro_dusk = search_altitude(astronomy.Direction.Set, -18.0)
|
|
|
|
return {
|
|
"ok": True,
|
|
"action": "sun_moon_rise_set",
|
|
"date": date_text,
|
|
"timezone": timezone_name,
|
|
"observer": {
|
|
"latitude": latitude,
|
|
"longitude": longitude,
|
|
"elevation": elevation,
|
|
},
|
|
"window": {
|
|
"local_start": local_start.isoformat(),
|
|
"local_end": local_end.isoformat(),
|
|
"utc_start": utc_start.isoformat().replace("+00:00", "Z"),
|
|
"utc_end": utc_end.isoformat().replace("+00:00", "Z"),
|
|
},
|
|
"events": {
|
|
"astro_dawn": serialize_event("Astronomische Morgendämmerung", astro_dawn, tz),
|
|
"nautical_dawn": serialize_event("Nautische Morgendämmerung", nautical_dawn, tz),
|
|
"civil_dawn": serialize_event("Bürgerliche Morgendämmerung", civil_dawn, tz),
|
|
"sunrise": serialize_event("Sonnenaufgang", sun_rise, tz),
|
|
"sunset": serialize_event("Sonnenuntergang", sun_set, tz),
|
|
"civil_dusk": serialize_event("Bürgerliche Abenddämmerung", civil_dusk, tz),
|
|
"nautical_dusk": serialize_event("Nautische Abenddämmerung", nautical_dusk, tz),
|
|
"astro_dusk": serialize_event("Astronomische Abenddämmerung", astro_dusk, tz),
|
|
"moonrise": serialize_event("Mondaufgang", moon_rise, tz),
|
|
"moonset": serialize_event("Monduntergang", moon_set, tz),
|
|
},
|
|
}
|
|
|
|
|
|
PLANETS = [
|
|
("Merkur", astronomy.Body.Mercury),
|
|
("Venus", astronomy.Body.Venus),
|
|
("Mars", astronomy.Body.Mars),
|
|
("Jupiter", astronomy.Body.Jupiter),
|
|
("Saturn", astronomy.Body.Saturn),
|
|
("Uranus", astronomy.Body.Uranus),
|
|
("Neptun", astronomy.Body.Neptune),
|
|
]
|
|
|
|
|
|
def action_planet_rise_set(args: list[str]) -> dict:
|
|
if len(args) != 5:
|
|
fail(
|
|
"Aktion planet_rise_set erwartet 5 Argumente: latitude longitude elevation date timezone",
|
|
extra={"argv": args},
|
|
)
|
|
|
|
latitude = parse_float(args[0], "Latitude")
|
|
longitude = parse_float(args[1], "Longitude")
|
|
elevation = parse_float(args[2], "Elevation")
|
|
date_text = args[3]
|
|
timezone_name = args[4]
|
|
|
|
try:
|
|
datetime.strptime(date_text, "%Y-%m-%d")
|
|
except ValueError:
|
|
fail("Datum muss im Format YYYY-MM-DD uebergeben werden.", extra={"date": date_text})
|
|
|
|
try:
|
|
tz = ZoneInfo(timezone_name)
|
|
except Exception as exc:
|
|
fail("Zeitzone ist ungueltig.", extra={"details": str(exc)})
|
|
|
|
observer = astronomy.Observer(latitude, longitude, elevation)
|
|
|
|
local_start = datetime.strptime(date_text, "%Y-%m-%d").replace(tzinfo=tz)
|
|
local_end = local_start + timedelta(days=1)
|
|
utc_start = local_start.astimezone(timezone.utc)
|
|
utc_end = local_end.astimezone(timezone.utc)
|
|
start_time = dt_to_time(utc_start)
|
|
|
|
planets = []
|
|
for name, body in PLANETS:
|
|
rise = search_event(body, astronomy.Direction.Rise, observer, start_time, utc_end)
|
|
set_ = search_event(body, astronomy.Direction.Set, observer, start_time, utc_end)
|
|
planets.append({
|
|
"name": name,
|
|
"rise": serialize_event(name + "aufgang", rise, tz),
|
|
"set": serialize_event(name + "untergang", set_, tz),
|
|
})
|
|
|
|
return {
|
|
"ok": True,
|
|
"action": "planet_rise_set",
|
|
"date": date_text,
|
|
"timezone": timezone_name,
|
|
"planets": planets,
|
|
}
|
|
|
|
|
|
def serialize_body_position(
|
|
key: str,
|
|
label: str,
|
|
body: astronomy.Body,
|
|
time_value: astronomy.Time,
|
|
observer: astronomy.Observer,
|
|
) -> dict:
|
|
eq = astronomy.Equator(body, time_value, observer, True, True)
|
|
hor = astronomy.Horizon(time_value, observer, eq.ra, eq.dec, astronomy.Refraction.Normal)
|
|
|
|
payload = {
|
|
"key": key,
|
|
"label": label,
|
|
"ra_hours": float(eq.ra),
|
|
"dec_deg": float(eq.dec),
|
|
"azimuth_deg": float(hor.azimuth),
|
|
"altitude_deg": float(hor.altitude),
|
|
"visible": float(hor.altitude) >= 0.0,
|
|
}
|
|
|
|
try:
|
|
illum = astronomy.Illumination(body, time_value)
|
|
payload["distance_au"] = float(illum.geo_dist)
|
|
payload["phase_fraction"] = float(illum.phase_fraction)
|
|
payload["magnitude"] = float(illum.mag)
|
|
except Exception:
|
|
pass
|
|
|
|
return payload
|
|
|
|
|
|
def action_month_sky_context(args: list[str]) -> dict:
|
|
if len(args) != 8:
|
|
fail(
|
|
"Aktion month_sky_context erwartet 8 Argumente: latitude longitude elevation year month hour minute timezone",
|
|
extra={"argv": args},
|
|
)
|
|
|
|
latitude = parse_float(args[0], "Latitude")
|
|
longitude = parse_float(args[1], "Longitude")
|
|
elevation = parse_float(args[2], "Elevation")
|
|
|
|
try:
|
|
year = int(args[3])
|
|
month = int(args[4])
|
|
hour = int(args[5])
|
|
minute = int(args[6])
|
|
except ValueError as exc:
|
|
fail("Jahr, Monat oder Uhrzeit sind ungueltig.", extra={"details": str(exc), "argv": args})
|
|
|
|
if month < 1 or month > 12:
|
|
fail("Monat muss zwischen 1 und 12 liegen.", extra={"month": month})
|
|
if hour < 0 or hour > 23:
|
|
fail("Stunde muss zwischen 0 und 23 liegen.", extra={"hour": hour})
|
|
if minute < 0 or minute > 59:
|
|
fail("Minute muss zwischen 0 und 59 liegen.", extra={"minute": minute})
|
|
|
|
timezone_name = args[7]
|
|
try:
|
|
tz = ZoneInfo(timezone_name)
|
|
except Exception as exc:
|
|
fail("Zeitzone ist ungueltig oder auf dem Server nicht verfuegbar.", extra={"details": str(exc)})
|
|
|
|
local_dt = datetime(year, month, 15, hour, minute, 0, tzinfo=tz)
|
|
utc_dt = local_dt.astimezone(timezone.utc)
|
|
time_value = dt_to_time(utc_dt)
|
|
observer = astronomy.Observer(latitude, longitude, elevation)
|
|
|
|
gst_hours = float(astronomy.SiderealTime(time_value))
|
|
lst_hours = normalize_degrees(gst_hours * 15.0 + longitude) / 15.0
|
|
|
|
sun = serialize_body_position("sun", "Sonne", astronomy.Body.Sun, time_value, observer)
|
|
moon = serialize_body_position("moon", "Mond", astronomy.Body.Moon, time_value, observer)
|
|
|
|
moon_phase_angle = normalize_degrees(astronomy.MoonPhase(time_value))
|
|
moon_age_days = moon_phase_angle / 360.0 * SYNODIC_MONTH
|
|
moon["phase_angle_deg"] = float(moon_phase_angle)
|
|
moon["phase_label"] = get_phase_label(moon_age_days)
|
|
moon["age_days"] = float(moon_age_days)
|
|
moon["waxing"] = moon_phase_angle <= 180.0
|
|
|
|
planets = [
|
|
serialize_body_position("mercury", "Merkur", astronomy.Body.Mercury, time_value, observer),
|
|
serialize_body_position("venus", "Venus", astronomy.Body.Venus, time_value, observer),
|
|
serialize_body_position("mars", "Mars", astronomy.Body.Mars, time_value, observer),
|
|
serialize_body_position("jupiter", "Jupiter", astronomy.Body.Jupiter, time_value, observer),
|
|
serialize_body_position("saturn", "Saturn", astronomy.Body.Saturn, time_value, observer),
|
|
serialize_body_position("uranus", "Uranus", astronomy.Body.Uranus, time_value, observer),
|
|
serialize_body_position("neptune", "Neptun", astronomy.Body.Neptune, time_value, observer),
|
|
]
|
|
|
|
return {
|
|
"ok": True,
|
|
"action": "month_sky_context",
|
|
"selected": {
|
|
"year": year,
|
|
"month": month,
|
|
"day": 15,
|
|
"hour": hour,
|
|
"minute": minute,
|
|
"timezone": timezone_name,
|
|
"local_iso": local_dt.isoformat(),
|
|
"utc_iso": utc_dt.isoformat().replace("+00:00", "Z"),
|
|
"local_label": local_dt.strftime("%d.%m.%Y %H:%M") + f" {timezone_name}",
|
|
},
|
|
"observer": {
|
|
"latitude": latitude,
|
|
"longitude": longitude,
|
|
"elevation": elevation,
|
|
},
|
|
"sidereal": {
|
|
"greenwich_hours": gst_hours,
|
|
"local_hours": lst_hours,
|
|
},
|
|
"bodies": {
|
|
"sun": sun,
|
|
"moon": moon,
|
|
"planets": planets,
|
|
},
|
|
}
|
|
|
|
|
|
def classify_visibility(sun_alt_deg: float) -> dict:
|
|
if sun_alt_deg > -5:
|
|
return {"key": "day", "label": "Tag"}
|
|
if sun_alt_deg >= -10:
|
|
return {"key": "twilight", "label": "Daemmerung"}
|
|
return {"key": "night", "label": "Nacht"}
|
|
|
|
|
|
def spherical_separation_deg(ra_hours_a: float, dec_deg_a: float, ra_hours_b: float, dec_deg_b: float) -> float:
|
|
deg2rad = math.pi / 180.0
|
|
ra_a = ra_hours_a * 15.0 * deg2rad
|
|
dec_a = dec_deg_a * deg2rad
|
|
ra_b = ra_hours_b * 15.0 * deg2rad
|
|
dec_b = dec_deg_b * deg2rad
|
|
cos_sep = (
|
|
(math.sin(dec_a) * math.sin(dec_b))
|
|
+ (math.cos(dec_a) * math.cos(dec_b) * math.cos(ra_a - ra_b))
|
|
)
|
|
cos_sep = max(-1.0, min(1.0, cos_sep))
|
|
return math.degrees(math.acos(cos_sep))
|
|
|
|
|
|
def moon_angular_radius_deg(distance_au: float) -> float:
|
|
if not math.isfinite(distance_au) or distance_au <= 0:
|
|
return 0.0
|
|
ratio = MOON_RADIUS_KM / (distance_au * astronomy.KM_PER_AU)
|
|
ratio = max(-1.0, min(1.0, ratio))
|
|
return math.degrees(math.asin(ratio))
|
|
|
|
|
|
def normalize_delta_ra_hours(delta_ra_hours: float) -> float:
|
|
while delta_ra_hours > 12.0:
|
|
delta_ra_hours -= 24.0
|
|
while delta_ra_hours < -12.0:
|
|
delta_ra_hours += 24.0
|
|
return delta_ra_hours
|
|
|
|
|
|
def star_ecliptic_latitude_deg(ra_hours: float, dec_deg: float, reference_time: astronomy.Time) -> float:
|
|
sphere = astronomy.Spherical(dec_deg, ra_hours * 15.0, 1.0)
|
|
vector = astronomy.VectorFromSphere(sphere, reference_time)
|
|
ecliptic = astronomy.Ecliptic(vector)
|
|
return float(ecliptic.elat)
|
|
|
|
|
|
def build_moon_track_samples(
|
|
observer: astronomy.Observer,
|
|
utc_start: datetime,
|
|
utc_end: datetime,
|
|
step_minutes: int = 30,
|
|
) -> list[dict]:
|
|
samples: list[dict] = []
|
|
current = utc_start
|
|
step = timedelta(minutes=step_minutes)
|
|
|
|
while current <= utc_end:
|
|
time_value = dt_to_time(current)
|
|
moon_eq = astronomy.Equator(astronomy.Body.Moon, time_value, observer, True, True)
|
|
samples.append(
|
|
{
|
|
"dt_utc": current,
|
|
"ra": float(moon_eq.ra),
|
|
"dec": float(moon_eq.dec),
|
|
"radius_deg": moon_angular_radius_deg(float(moon_eq.dist)),
|
|
}
|
|
)
|
|
current += step
|
|
|
|
if not samples or samples[-1]["dt_utc"] < utc_end:
|
|
time_value = dt_to_time(utc_end)
|
|
moon_eq = astronomy.Equator(astronomy.Body.Moon, time_value, observer, True, True)
|
|
samples.append(
|
|
{
|
|
"dt_utc": utc_end,
|
|
"ra": float(moon_eq.ra),
|
|
"dec": float(moon_eq.dec),
|
|
"radius_deg": moon_angular_radius_deg(float(moon_eq.dist)),
|
|
}
|
|
)
|
|
|
|
return samples
|
|
|
|
|
|
def interpolate_track_value(left: float, right: float, fraction: float, *, wrap_hours: bool = False) -> float:
|
|
if wrap_hours:
|
|
delta = normalize_delta_ra_hours(right - left)
|
|
value = left + delta * fraction
|
|
while value < 0.0:
|
|
value += 24.0
|
|
while value >= 24.0:
|
|
value -= 24.0
|
|
return value
|
|
return left + (right - left) * fraction
|
|
|
|
|
|
def build_moon_track_segments(samples: list[dict]) -> list[dict]:
|
|
segments: list[dict] = []
|
|
|
|
for index in range(len(samples) - 1):
|
|
left = samples[index]
|
|
right = samples[index + 1]
|
|
duration_seconds = (right["dt_utc"] - left["dt_utc"]).total_seconds()
|
|
if duration_seconds <= 0:
|
|
continue
|
|
|
|
segments.append(
|
|
{
|
|
"start": left["dt_utc"],
|
|
"end": right["dt_utc"],
|
|
"duration_seconds": duration_seconds,
|
|
"start_ra": left["ra"],
|
|
"end_ra": right["ra"],
|
|
"start_dec": left["dec"],
|
|
"end_dec": right["dec"],
|
|
"start_radius_deg": left["radius_deg"],
|
|
"end_radius_deg": right["radius_deg"],
|
|
}
|
|
)
|
|
|
|
return segments
|
|
|
|
|
|
def minimum_interpolated_moon_distance_deg(
|
|
star_ra_hours: float,
|
|
star_dec_deg: float,
|
|
track_segments: list[dict],
|
|
candidate_padding_deg: float,
|
|
) -> float:
|
|
minimum_distance = float("inf")
|
|
|
|
for segment in track_segments:
|
|
for sample_index in range(7):
|
|
fraction = sample_index / 6.0
|
|
moon_ra = interpolate_track_value(segment["start_ra"], segment["end_ra"], fraction, wrap_hours=True)
|
|
moon_dec = interpolate_track_value(segment["start_dec"], segment["end_dec"], fraction)
|
|
moon_radius = interpolate_track_value(segment["start_radius_deg"], segment["end_radius_deg"], fraction)
|
|
|
|
if abs(star_dec_deg - moon_dec) > (candidate_padding_deg + moon_radius):
|
|
continue
|
|
|
|
distance_deg = spherical_separation_deg(star_ra_hours, star_dec_deg, moon_ra, moon_dec)
|
|
margin_deg = distance_deg - moon_radius
|
|
if margin_deg < minimum_distance:
|
|
minimum_distance = margin_deg
|
|
|
|
return minimum_distance
|
|
|
|
|
|
def evaluate_occultation_geometry(body: astronomy.Body, observer: astronomy.Observer, dt_utc: datetime) -> dict:
|
|
time_value = dt_to_time(dt_utc)
|
|
moon_eq = astronomy.Equator(astronomy.Body.Moon, time_value, observer, True, True)
|
|
star_eq = astronomy.Equator(body, time_value, observer, True, True)
|
|
sun_eq = astronomy.Equator(astronomy.Body.Sun, time_value, observer, True, True)
|
|
|
|
moon_hor = astronomy.Horizon(time_value, observer, moon_eq.ra, moon_eq.dec, astronomy.Refraction.Normal)
|
|
star_hor = astronomy.Horizon(time_value, observer, star_eq.ra, star_eq.dec, astronomy.Refraction.Normal)
|
|
sun_hor = astronomy.Horizon(time_value, observer, sun_eq.ra, sun_eq.dec, astronomy.Refraction.Normal)
|
|
|
|
separation_deg = spherical_separation_deg(moon_eq.ra, moon_eq.dec, star_eq.ra, star_eq.dec)
|
|
radius_deg = moon_angular_radius_deg(moon_eq.dist)
|
|
margin_deg = separation_deg - radius_deg
|
|
|
|
return {
|
|
"time": time_value,
|
|
"moon_alt_deg": float(moon_hor.altitude),
|
|
"star_alt_deg": float(star_hor.altitude),
|
|
"sun_alt_deg": float(sun_hor.altitude),
|
|
"margin_deg": float(margin_deg),
|
|
}
|
|
|
|
|
|
def find_root_time(body: astronomy.Body, observer: astronomy.Observer, start_utc: datetime, end_utc: datetime) -> datetime:
|
|
left = start_utc
|
|
right = end_utc
|
|
for _ in range(40):
|
|
mid = left + (right - left) / 2
|
|
geom = evaluate_occultation_geometry(body, observer, mid)
|
|
if geom["margin_deg"] <= 0:
|
|
right = mid
|
|
else:
|
|
left = mid
|
|
return left + (right - left) / 2
|
|
|
|
|
|
def find_minimum_time(body: astronomy.Body, observer: astronomy.Observer, start_utc: datetime, end_utc: datetime) -> datetime:
|
|
left = start_utc
|
|
right = end_utc
|
|
for _ in range(40):
|
|
span = (right - left) / 3
|
|
m1 = left + span
|
|
m2 = right - span
|
|
f1 = evaluate_occultation_geometry(body, observer, m1)["margin_deg"]
|
|
f2 = evaluate_occultation_geometry(body, observer, m2)["margin_deg"]
|
|
if f1 <= f2:
|
|
right = m2
|
|
else:
|
|
left = m1
|
|
return left + (right - left) / 2
|
|
|
|
|
|
def serialize_occultation_time(dt_utc: datetime | None, tz: ZoneInfo) -> dict | None:
|
|
if dt_utc is None:
|
|
return None
|
|
local_dt = dt_utc.astimezone(tz)
|
|
return {
|
|
"utc_iso": dt_utc.isoformat().replace("+00:00", "Z"),
|
|
"local_iso": local_dt.isoformat(),
|
|
"local_time": local_dt.strftime("%H:%M"),
|
|
"local_date": local_dt.strftime("%Y-%m-%d"),
|
|
}
|
|
|
|
|
|
def load_star_file(star_file_path: str) -> list[dict]:
|
|
try:
|
|
with open(star_file_path, "r", encoding="utf-8-sig") as handle:
|
|
payload = json.load(handle)
|
|
except OSError as exc:
|
|
fail("Sterndatei konnte nicht gelesen werden.", extra={"details": str(exc), "path": star_file_path})
|
|
except json.JSONDecodeError as exc:
|
|
fail("Sterndatei enthaelt kein gueltiges JSON.", extra={"details": str(exc), "path": star_file_path})
|
|
|
|
if not isinstance(payload, list):
|
|
fail("Sterndatei muss eine JSON-Liste sein.", extra={"path": star_file_path})
|
|
return payload
|
|
|
|
|
|
def normalize_degrees(value: float) -> float:
|
|
return ((value % 360.0) + 360.0) % 360.0
|
|
|
|
|
|
def normalize_signed_degrees(value: float) -> float:
|
|
return normalize_degrees(value + 180.0) - 180.0
|
|
|
|
|
|
def get_phase_label(age_days: float) -> str:
|
|
segment = SYNODIC_MONTH / 8.0
|
|
index = int(math.floor((age_days + segment / 2.0) / segment)) % 8
|
|
return PHASE_LABELS[index]
|
|
|
|
|
|
def vector_from_astronomy(vector: astronomy.Vector) -> tuple[float, float, float]:
|
|
return float(vector.x), float(vector.y), float(vector.z)
|
|
|
|
|
|
def vector_subtract(a: tuple[float, float, float], b: tuple[float, float, float]) -> tuple[float, float, float]:
|
|
return a[0] - b[0], a[1] - b[1], a[2] - b[2]
|
|
|
|
|
|
def vector_dot(a: tuple[float, float, float], b: tuple[float, float, float]) -> float:
|
|
return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]
|
|
|
|
|
|
def vector_normalize(v: tuple[float, float, float]) -> tuple[float, float, float]:
|
|
length = math.sqrt(vector_dot(v, v))
|
|
if length <= 0:
|
|
return 0.0, 0.0, 0.0
|
|
return v[0] / length, v[1] / length, v[2] / length
|
|
|
|
|
|
def calculate_moon_position_angle(time_value: astronomy.Time) -> float:
|
|
moon_vector = astronomy.GeoVector(astronomy.Body.Moon, time_value, True)
|
|
moon_equator = astronomy.EquatorFromVector(moon_vector)
|
|
axis = astronomy.RotationAxis(astronomy.Body.Moon, time_value)
|
|
pole_equator = astronomy.EquatorFromVector(axis.north)
|
|
|
|
delta_ra = math.radians((float(pole_equator.ra) - float(moon_equator.ra)) * 15.0)
|
|
moon_dec = math.radians(float(moon_equator.dec))
|
|
pole_dec = math.radians(float(pole_equator.dec))
|
|
|
|
return math.degrees(
|
|
math.atan2(
|
|
math.cos(pole_dec) * math.sin(delta_ra),
|
|
math.sin(pole_dec) * math.cos(moon_dec)
|
|
- math.cos(pole_dec) * math.sin(moon_dec) * math.cos(delta_ra),
|
|
)
|
|
)
|
|
|
|
|
|
def calculate_moon_axis_latitudes(time_value: astronomy.Time) -> dict:
|
|
axis = astronomy.RotationAxis(astronomy.Body.Moon, time_value)
|
|
moon_vector = astronomy.GeoVector(astronomy.Body.Moon, time_value, True)
|
|
sun_vector = astronomy.GeoVector(astronomy.Body.Sun, time_value, True)
|
|
|
|
north = vector_normalize(vector_from_astronomy(axis.north))
|
|
moon_xyz = vector_from_astronomy(moon_vector)
|
|
sun_xyz = vector_from_astronomy(sun_vector)
|
|
|
|
earth_from_moon = vector_normalize((-moon_xyz[0], -moon_xyz[1], -moon_xyz[2]))
|
|
sun_from_moon = vector_normalize(vector_subtract(sun_xyz, moon_xyz))
|
|
|
|
return {
|
|
"subearth_latitude": math.degrees(math.asin(max(-1.0, min(1.0, vector_dot(earth_from_moon, north))))),
|
|
"subsolar_latitude": math.degrees(math.asin(max(-1.0, min(1.0, vector_dot(sun_from_moon, north))))),
|
|
}
|
|
|
|
|
|
def action_moon_phase_details(args: list[str]) -> dict:
|
|
if len(args) != 2:
|
|
fail(
|
|
"Aktion moon_phase_details erwartet 2 Argumente: utc_iso timezone",
|
|
extra={"argv": args},
|
|
)
|
|
|
|
utc_iso = args[0]
|
|
timezone_name = args[1]
|
|
|
|
try:
|
|
tz = ZoneInfo(timezone_name)
|
|
except Exception as exc: # pragma: no cover
|
|
fail("Zeitzone ist ungueltig oder auf dem Server nicht verfuegbar.", extra={"details": str(exc)})
|
|
|
|
try:
|
|
dt_utc = datetime.fromisoformat(utc_iso.replace("Z", "+00:00")).astimezone(timezone.utc)
|
|
except ValueError as exc:
|
|
fail("utc_iso ist ungueltig.", extra={"details": str(exc), "utc_iso": utc_iso})
|
|
|
|
time_value = dt_to_time(dt_utc)
|
|
illumination = astronomy.Illumination(astronomy.Body.Moon, time_value)
|
|
elongation = normalize_degrees(astronomy.MoonPhase(time_value))
|
|
phase = elongation / 360.0
|
|
previous_new_moon = astronomy.SearchMoonPhase(0.0, time_value, -35.0)
|
|
age_days = (
|
|
(dt_utc - time_to_datetime(previous_new_moon)).total_seconds() / 86400.0
|
|
if previous_new_moon is not None
|
|
else phase * SYNODIC_MONTH
|
|
)
|
|
|
|
libration = astronomy.Libration(time_value)
|
|
axis_latitudes = calculate_moon_axis_latitudes(time_value)
|
|
position_angle = calculate_moon_position_angle(time_value)
|
|
local_dt = dt_utc.astimezone(tz)
|
|
|
|
return {
|
|
"ok": True,
|
|
"action": "moon_phase_details",
|
|
"selected": {
|
|
"utc_iso": dt_utc.isoformat().replace("+00:00", "Z"),
|
|
"local_iso": local_dt.isoformat(),
|
|
"timezone": timezone_name,
|
|
},
|
|
"phase": {
|
|
"age": age_days,
|
|
"phase": phase,
|
|
"illumination": float(illumination.phase_fraction),
|
|
"waxing": elongation <= 180.0,
|
|
"label": get_phase_label(age_days),
|
|
},
|
|
"libration": {
|
|
"longitude": -normalize_signed_degrees(float(libration.elon)),
|
|
"latitude": float(libration.elat),
|
|
"subsolarLatitude": float(axis_latitudes["subsolar_latitude"]),
|
|
"subearthLatitude": float(axis_latitudes["subearth_latitude"]),
|
|
"positionAngle": float(position_angle),
|
|
"distKm": float(libration.dist_km),
|
|
},
|
|
}
|
|
|
|
|
|
def compute_star_occultation(
|
|
star_rows: list[dict],
|
|
observer: astronomy.Observer,
|
|
local_start: datetime,
|
|
local_end: datetime,
|
|
tz: ZoneInfo,
|
|
) -> tuple[list[dict], dict]:
|
|
utc_start = local_start.astimezone(timezone.utc)
|
|
utc_end = local_end.astimezone(timezone.utc)
|
|
search_margin = timedelta(hours=8)
|
|
search_start = utc_start - search_margin
|
|
search_end = utc_end + search_margin
|
|
results: list[dict] = []
|
|
coarse_step = timedelta(minutes=2)
|
|
candidate_margin_deg = 0.32
|
|
ecliptic_band_deg = 8.0
|
|
reference_time = dt_to_time(utc_start)
|
|
moon_track_samples = build_moon_track_samples(observer, search_start, search_end, step_minutes=30)
|
|
moon_track_segments = build_moon_track_segments(moon_track_samples)
|
|
geometry_candidates: list[dict] = []
|
|
|
|
for star in star_rows:
|
|
star_ra = float(star["ra"])
|
|
star_dec = float(star["dec"])
|
|
ecliptic_lat_deg = abs(star_ecliptic_latitude_deg(star_ra, star_dec, reference_time))
|
|
|
|
if ecliptic_lat_deg > ecliptic_band_deg:
|
|
continue
|
|
|
|
estimated_min_margin_deg = minimum_interpolated_moon_distance_deg(
|
|
star_ra,
|
|
star_dec,
|
|
moon_track_segments,
|
|
candidate_margin_deg,
|
|
)
|
|
|
|
if estimated_min_margin_deg > candidate_margin_deg:
|
|
continue
|
|
|
|
geometry_candidates.append(star)
|
|
|
|
if not geometry_candidates:
|
|
debug = {
|
|
"input_star_count": len(star_rows),
|
|
"after_ecliptic_filter": 0,
|
|
"after_track_filter": 0,
|
|
"ecliptic_band_deg": ecliptic_band_deg,
|
|
"candidate_margin_deg": candidate_margin_deg,
|
|
}
|
|
return results, debug
|
|
|
|
ecliptic_count = 0
|
|
for star in star_rows:
|
|
star_ra = float(star["ra"])
|
|
star_dec = float(star["dec"])
|
|
if abs(star_ecliptic_latitude_deg(star_ra, star_dec, reference_time)) <= ecliptic_band_deg:
|
|
ecliptic_count += 1
|
|
|
|
for chunk_start in range(0, len(geometry_candidates), len(STAR_BODIES)):
|
|
chunk = geometry_candidates[chunk_start:chunk_start + len(STAR_BODIES)]
|
|
for body, star in zip(STAR_BODIES, chunk):
|
|
astronomy.DefineStar(
|
|
body,
|
|
float(star["ra"]),
|
|
float(star["dec"]),
|
|
max(1.0, float(star.get("distLy", 1000.0))),
|
|
)
|
|
|
|
for index, star in enumerate(chunk):
|
|
body = STAR_BODIES[index]
|
|
previous_margin = None
|
|
previous_time = None
|
|
ingress_bracket = None
|
|
egress_bracket = None
|
|
minimum_margin = float("inf")
|
|
minimum_time = search_start
|
|
start_margin = None
|
|
end_margin = None
|
|
ever_negative = False
|
|
near_hit = False
|
|
|
|
current = search_start
|
|
while current <= search_end:
|
|
geometry = evaluate_occultation_geometry(body, observer, current)
|
|
margin = geometry["margin_deg"]
|
|
|
|
if margin < minimum_margin:
|
|
minimum_margin = margin
|
|
minimum_time = current
|
|
|
|
if margin <= 0:
|
|
ever_negative = True
|
|
|
|
if margin <= candidate_margin_deg:
|
|
near_hit = True
|
|
|
|
if start_margin is None:
|
|
start_margin = margin
|
|
end_margin = margin
|
|
|
|
if previous_margin is not None and previous_time is not None:
|
|
if (previous_margin > 0 >= margin) and ingress_bracket is None:
|
|
ingress_bracket = (previous_time, current)
|
|
near_hit = True
|
|
if (previous_margin <= 0 < margin) and egress_bracket is None:
|
|
egress_bracket = (previous_time, current)
|
|
near_hit = True
|
|
|
|
previous_margin = margin
|
|
previous_time = current
|
|
current += coarse_step
|
|
|
|
if (not ever_negative) and (not near_hit):
|
|
continue
|
|
if (not ever_negative) or minimum_margin > 0:
|
|
continue
|
|
|
|
ingress_time = find_root_time(body, observer, ingress_bracket[0], ingress_bracket[1]) if ingress_bracket else None
|
|
egress_time = find_root_time(body, observer, egress_bracket[0], egress_bracket[1]) if egress_bracket else None
|
|
|
|
minimum_search_start = max(search_start, minimum_time - coarse_step)
|
|
minimum_search_end = min(search_end, minimum_time + coarse_step)
|
|
if ingress_time is not None and egress_time is not None:
|
|
minimum_search_start = ingress_time
|
|
minimum_search_end = egress_time
|
|
|
|
max_time = find_minimum_time(body, observer, minimum_search_start, minimum_search_end)
|
|
max_geometry = evaluate_occultation_geometry(body, observer, max_time)
|
|
|
|
if max_geometry["margin_deg"] > 0:
|
|
continue
|
|
if max_geometry["moon_alt_deg"] <= 0 or max_geometry["star_alt_deg"] <= 0:
|
|
continue
|
|
|
|
event_start = ingress_time or search_start
|
|
event_end = egress_time or search_end
|
|
if event_end < utc_start or event_start >= utc_end:
|
|
continue
|
|
|
|
duration_seconds = None
|
|
if ingress_time is not None and egress_time is not None:
|
|
duration_seconds = int(round((egress_time - ingress_time).total_seconds()))
|
|
|
|
results.append({
|
|
"star": {
|
|
"hip": int(star.get("hip", 0)),
|
|
"label": str(star.get("label", "Unbenannter Stern")),
|
|
"constellation": str(star.get("constellation", "")),
|
|
"mag": float(star.get("mag", 99.0)),
|
|
"ra": float(star.get("ra", 0.0)),
|
|
"dec": float(star.get("dec", 0.0)),
|
|
},
|
|
"ingress": serialize_occultation_time(ingress_time, tz),
|
|
"maximum": serialize_occultation_time(max_time, tz),
|
|
"egress": serialize_occultation_time(egress_time, tz),
|
|
"partial_start": ingress_time is None and start_margin is not None and start_margin <= 0,
|
|
"partial_end": egress_time is None and end_margin is not None and end_margin <= 0,
|
|
"duration_seconds": duration_seconds,
|
|
"limb_distance_arcmin": abs(max_geometry["margin_deg"]) * 60.0,
|
|
"moon_alt_deg": max_geometry["moon_alt_deg"],
|
|
"star_alt_deg": max_geometry["star_alt_deg"],
|
|
"sun_alt_deg": max_geometry["sun_alt_deg"],
|
|
"visibility": classify_visibility(max_geometry["sun_alt_deg"]),
|
|
})
|
|
|
|
results.sort(
|
|
key=lambda item: item["maximum"]["utc_iso"] if item["maximum"] else "9999-99-99T99:99:99Z"
|
|
)
|
|
debug = {
|
|
"input_star_count": len(star_rows),
|
|
"after_ecliptic_filter": ecliptic_count,
|
|
"after_track_filter": len(geometry_candidates),
|
|
"ecliptic_band_deg": ecliptic_band_deg,
|
|
"candidate_margin_deg": candidate_margin_deg,
|
|
"search_margin_hours": search_margin.total_seconds() / 3600.0,
|
|
}
|
|
return results, debug
|
|
|
|
|
|
def action_moon_star_occultations(args: list[str]) -> dict:
|
|
if len(args) != 6:
|
|
fail(
|
|
"Aktion moon_star_occultations erwartet 6 Argumente: latitude longitude elevation date timezone star_file",
|
|
extra={"argv": args},
|
|
)
|
|
|
|
latitude = parse_float(args[0], "Latitude")
|
|
longitude = parse_float(args[1], "Longitude")
|
|
elevation = parse_float(args[2], "Elevation")
|
|
date_text = args[3]
|
|
timezone_name = args[4]
|
|
star_file_path = args[5]
|
|
|
|
try:
|
|
datetime.strptime(date_text, "%Y-%m-%d")
|
|
except ValueError:
|
|
fail("Datum muss im Format YYYY-MM-DD uebergeben werden.", extra={"date": date_text})
|
|
|
|
try:
|
|
tz = ZoneInfo(timezone_name)
|
|
except Exception as exc: # pragma: no cover
|
|
fail("Zeitzone ist ungueltig oder auf dem Server nicht verfuegbar.", extra={"details": str(exc)})
|
|
|
|
star_rows = load_star_file(star_file_path)
|
|
observer = astronomy.Observer(latitude, longitude, elevation)
|
|
local_start = datetime.strptime(date_text, "%Y-%m-%d").replace(tzinfo=tz)
|
|
local_end = local_start + timedelta(days=1)
|
|
|
|
results, debug_profile = compute_star_occultation(star_rows, observer, local_start, local_end, tz)
|
|
|
|
return {
|
|
"ok": True,
|
|
"action": "moon_star_occultations",
|
|
"date": date_text,
|
|
"timezone": timezone_name,
|
|
"observer": {
|
|
"latitude": latitude,
|
|
"longitude": longitude,
|
|
"elevation": elevation,
|
|
},
|
|
"star_count": len(star_rows),
|
|
"window": {
|
|
"local_start": local_start.isoformat(),
|
|
"local_end": local_end.isoformat(),
|
|
"utc_start": local_start.astimezone(timezone.utc).isoformat().replace("+00:00", "Z"),
|
|
"utc_end": local_end.astimezone(timezone.utc).isoformat().replace("+00:00", "Z"),
|
|
},
|
|
"debug_profile": debug_profile,
|
|
"results": results,
|
|
}
|
|
|
|
|
|
# ── satellite_passes ──────────────────────────────────────────────────────────
|
|
|
|
import math as _math
|
|
|
|
_WGS84_A = 6378.137 # km
|
|
_WGS84_E2 = 0.00669437999014
|
|
_R_EARTH = 6371.0 # km (shadow check)
|
|
_KM_PER_AU_SAT = 149597870.7
|
|
_K_EXT = 0.18
|
|
_PHI90 = 1.0 / _math.pi
|
|
|
|
|
|
def _jday_from_dt(dt_utc: datetime) -> tuple[float, float]:
|
|
"""Return (jd_int, jd_fraction) from a UTC datetime."""
|
|
from sgp4.functions import jday as _jday_fn
|
|
return _jday_fn(dt_utc.year, dt_utc.month, dt_utc.day,
|
|
dt_utc.hour, dt_utc.minute,
|
|
dt_utc.second + dt_utc.microsecond / 1e6)
|
|
|
|
|
|
def _gmst_rad(jd_full: float) -> float:
|
|
T = (jd_full - 2451545.0) / 36525.0
|
|
gmst_deg = (280.46061837
|
|
+ 360.98564736629 * (jd_full - 2451545.0)
|
|
+ 0.000387933 * T * T
|
|
- T * T * T / 38710000.0)
|
|
return _math.radians(gmst_deg % 360.0)
|
|
|
|
|
|
def _eci_to_ecef(r_eci: tuple, gmst: float) -> tuple:
|
|
x, y, z = r_eci
|
|
c, s = _math.cos(gmst), _math.sin(gmst)
|
|
return (x * c + y * s, -x * s + y * c, z)
|
|
|
|
|
|
def _observer_ecef(lat_rad: float, lon_rad: float, elev_km: float) -> tuple:
|
|
N = _WGS84_A / _math.sqrt(1.0 - _WGS84_E2 * _math.sin(lat_rad) ** 2)
|
|
x = (N + elev_km) * _math.cos(lat_rad) * _math.cos(lon_rad)
|
|
y = (N + elev_km) * _math.cos(lat_rad) * _math.sin(lon_rad)
|
|
z = (N * (1.0 - _WGS84_E2) + elev_km) * _math.sin(lat_rad)
|
|
return (x, y, z)
|
|
|
|
|
|
def _ecef_to_altaz(r_ecef: tuple, obs_ecef: tuple, lat_rad: float, lon_rad: float) -> tuple[float, float, float]:
|
|
dx = r_ecef[0] - obs_ecef[0]
|
|
dy = r_ecef[1] - obs_ecef[1]
|
|
dz = r_ecef[2] - obs_ecef[2]
|
|
rng = _math.sqrt(dx*dx + dy*dy + dz*dz)
|
|
if rng < 1e-9:
|
|
return 0.0, 0.0, 0.0
|
|
|
|
# Topocentric South-East-Z (SEZ)
|
|
slat, clat = _math.sin(lat_rad), _math.cos(lat_rad)
|
|
slon, clon = _math.sin(lon_rad), _math.cos(lon_rad)
|
|
s = slat * clon * dx + slat * slon * dy - clat * dz
|
|
e = -slon * dx + clon * dy
|
|
z_top = clat * clon * dx + clat * slon * dy + slat * dz
|
|
|
|
el_rad = _math.asin(max(-1.0, min(1.0, z_top / rng)))
|
|
az_rad = _math.atan2(e, -s)
|
|
az_deg = (_math.degrees(az_rad) + 360.0) % 360.0
|
|
el_deg = _math.degrees(el_rad)
|
|
return el_deg, az_deg, rng
|
|
|
|
|
|
def _sun_eci_km(dt_utc: datetime) -> tuple:
|
|
t = dt_to_time(dt_utc)
|
|
vec = astronomy.GeoVector(astronomy.Body.Sun, t, True)
|
|
return (float(vec.x) * _KM_PER_AU_SAT,
|
|
float(vec.y) * _KM_PER_AU_SAT,
|
|
float(vec.z) * _KM_PER_AU_SAT)
|
|
|
|
|
|
def _in_shadow(r_eci: tuple, sun_eci_km: tuple) -> bool:
|
|
# Sun vector from Earth
|
|
sx, sy, sz = sun_eci_km
|
|
sun_len = _math.sqrt(sx*sx + sy*sy + sz*sz)
|
|
if sun_len < 1:
|
|
return False
|
|
# Satellite-to-sun vector direction
|
|
rx, ry, rz = r_eci
|
|
# Project satellite onto anti-sun direction
|
|
dot = -(rx*sx + ry*sy + rz*sz) / sun_len
|
|
if dot < 0:
|
|
return False # satellite on sun-side
|
|
perp2 = (rx*rx + ry*ry + rz*rz) - dot*dot
|
|
return perp2 < _R_EARTH * _R_EARTH
|
|
|
|
|
|
def _phi_lambert(alpha_rad: float) -> float:
|
|
if alpha_rad >= _math.pi:
|
|
return 1e-10
|
|
v = (_math.sin(alpha_rad) + (_math.pi - alpha_rad) * _math.cos(alpha_rad)) / _math.pi
|
|
return max(1e-10, v)
|
|
|
|
|
|
def _airmass(h_deg: float) -> float:
|
|
if h_deg <= 0:
|
|
return 40.0
|
|
sin_h = _math.sin(_math.radians(h_deg))
|
|
return 1.0 / (sin_h + 0.50572 * (h_deg + 6.07995) ** -1.6364)
|
|
|
|
|
|
def _apparent_magnitude(std_mag, rcs, rcs_size, r_eci, sun_eci_km, obs_ecef_km, lat_rad, lon_rad, gmst):
|
|
if std_mag is not None:
|
|
m1000 = float(std_mag)
|
|
estimated = False
|
|
elif rcs is not None and rcs > 0:
|
|
m1000 = 7.5 - 2.5 * _math.log10(rcs)
|
|
estimated = True
|
|
else:
|
|
size_map = {"LARGE": 3.5, "MEDIUM": 6.0, "SMALL": 8.5}
|
|
key = (rcs_size or "").upper()
|
|
if key not in size_map:
|
|
return None
|
|
m1000 = size_map[key]
|
|
estimated = True
|
|
|
|
rx, ry, rz = r_eci
|
|
r_ecef = _eci_to_ecef(r_eci, gmst)
|
|
dx = r_ecef[0] - obs_ecef_km[0]
|
|
dy = r_ecef[1] - obs_ecef_km[1]
|
|
dz = r_ecef[2] - obs_ecef_km[2]
|
|
range_km = _math.sqrt(dx*dx + dy*dy + dz*dz)
|
|
if range_km < 1:
|
|
return None
|
|
|
|
# Phase angle (sun → sat → observer)
|
|
sx, sy, sz = sun_eci_km
|
|
ox = -rx + obs_ecef_km[0] # obs ECI ≈ ECEF for this approximation
|
|
oy = -ry + obs_ecef_km[1]
|
|
oz = -rz + obs_ecef_km[2]
|
|
# to-sun from sat
|
|
ts_x, ts_y, ts_z = sx - rx, sy - ry, sz - rz
|
|
to_x, to_y, to_z = -rx, -ry, -rz # observer ≈ Earth center approx
|
|
ts_len = _math.sqrt(ts_x*ts_x + ts_y*ts_y + ts_z*ts_z)
|
|
to_len = _math.sqrt(to_x*to_x + to_y*to_y + to_z*to_z)
|
|
if ts_len < 1 or to_len < 1:
|
|
return None
|
|
cos_a = (ts_x*to_x + ts_y*to_y + ts_z*to_z) / (ts_len * to_len)
|
|
alpha_rad = _math.acos(max(-1.0, min(1.0, cos_a)))
|
|
|
|
el_deg, _, _ = _ecef_to_altaz(r_ecef, obs_ecef_km, lat_rad, lon_rad)
|
|
dist_corr = 5.0 * _math.log10(range_km / 1000.0)
|
|
phase_corr = -2.5 * _math.log10(_phi_lambert(alpha_rad) / _PHI90)
|
|
extinction = _K_EXT * _airmass(max(1.0, el_deg))
|
|
mag = m1000 + dist_corr + phase_corr + extinction
|
|
if not _math.isfinite(mag):
|
|
return None
|
|
return {"mag": round(mag, 2), "estimated": estimated}
|
|
|
|
|
|
def _bisect_crossing(satrec, obs_ecef, lat_rad, lon_rad, t1: datetime, t2: datetime, rising: bool) -> datetime:
|
|
for _ in range(30):
|
|
mid = t1 + (t2 - t1) / 2
|
|
jd, fr = _jday_from_dt(mid)
|
|
e, r, _ = satrec.sgp4(jd, fr)
|
|
if e != 0:
|
|
break
|
|
gmst = _gmst_rad(jd + fr)
|
|
r_ecef = _eci_to_ecef(r, gmst)
|
|
el, _, _ = _ecef_to_altaz(r_ecef, obs_ecef, lat_rad, lon_rad)
|
|
above = el >= 0
|
|
if rising:
|
|
if above:
|
|
t2 = mid
|
|
else:
|
|
t1 = mid
|
|
else:
|
|
if above:
|
|
t1 = mid
|
|
else:
|
|
t2 = mid
|
|
if (t2 - t1).total_seconds() < 0.5:
|
|
break
|
|
return t1 + (t2 - t1) / 2
|
|
|
|
|
|
def action_satellite_passes(args: list[str]) -> dict:
|
|
if len(args) != 7:
|
|
fail(
|
|
"Aktion satellite_passes erwartet 7 Argumente: latitude longitude elevation timezone window_start_iso window_end_iso satellites_payload",
|
|
extra={"argv": args},
|
|
)
|
|
|
|
latitude = parse_float(args[0], "Latitude")
|
|
longitude = parse_float(args[1], "Longitude")
|
|
elevation_m = parse_float(args[2], "Elevation")
|
|
timezone_name = args[3]
|
|
window_start_iso = args[4]
|
|
window_end_iso = args[5]
|
|
satellites_payload = args[6]
|
|
|
|
try:
|
|
tz = ZoneInfo(timezone_name)
|
|
except Exception as exc:
|
|
fail("Zeitzone ungueltig.", extra={"details": str(exc)})
|
|
|
|
try:
|
|
window_start = datetime.fromisoformat(window_start_iso.replace("Z", "+00:00")).astimezone(timezone.utc)
|
|
window_end = datetime.fromisoformat(window_end_iso.replace("Z", "+00:00")).astimezone(timezone.utc)
|
|
except ValueError as exc:
|
|
fail("Zeitfenster-ISO ungueltig.", extra={"details": str(exc)})
|
|
|
|
try:
|
|
satellites_json = base64.b64decode(satellites_payload.encode("ascii")).decode("utf-8")
|
|
satellites = json.loads(satellites_json)
|
|
except (ValueError, UnicodeDecodeError, json.JSONDecodeError) as exc:
|
|
fail("satellites_payload ist ungueltig.", extra={"details": str(exc)})
|
|
|
|
try:
|
|
import os as _os
|
|
_sgp4_dir = _os.path.join(SCRIPT_DIR, "sgp4")
|
|
if _sgp4_dir not in sys.path:
|
|
sys.path.insert(0, _sgp4_dir)
|
|
from sgp4.api import Satrec
|
|
except ImportError as exc:
|
|
fail("sgp4-Bibliothek konnte nicht importiert werden.", extra={"details": str(exc)})
|
|
|
|
lat_rad = _math.radians(latitude)
|
|
lon_rad = _math.radians(longitude)
|
|
elev_km = elevation_m / 1000.0
|
|
obs_ecef = _observer_ecef(lat_rad, lon_rad, elev_km)
|
|
|
|
SCAN_STEP = timedelta(seconds=30)
|
|
TRAJ_STEP = timedelta(seconds=5)
|
|
COLORS = ['#f0d990', '#4fc3d8', '#ff9f68', '#9ad77d', '#c7b0ff',
|
|
'#ffd166', '#7bdff2', '#ff7b9c', '#b8f2e6', '#f7a072']
|
|
|
|
results = []
|
|
|
|
for idx, sat_item in enumerate(satellites):
|
|
tle1 = str(sat_item.get("tle_line1", ""))
|
|
tle2 = str(sat_item.get("tle_line2", ""))
|
|
if not tle1 or not tle2:
|
|
continue
|
|
try:
|
|
satrec = Satrec.twoline2rv(tle1, tle2)
|
|
except Exception:
|
|
continue
|
|
|
|
color = COLORS[idx % len(COLORS)]
|
|
passes = []
|
|
|
|
# ── Phase 1: coarse scan to find pass windows ──
|
|
current = window_start
|
|
prev_el = None
|
|
in_pass = False
|
|
aos_time = None
|
|
aos_az = None
|
|
max_el = -999.0
|
|
max_el_time = None
|
|
prev_time = None
|
|
|
|
while current <= window_end:
|
|
jd, fr = _jday_from_dt(current)
|
|
e, r, _ = satrec.sgp4(jd, fr)
|
|
if e != 0:
|
|
current += SCAN_STEP
|
|
prev_el = None
|
|
continue
|
|
gmst = _gmst_rad(jd + fr)
|
|
r_ecef = _eci_to_ecef(r, gmst)
|
|
el, az, _ = _ecef_to_altaz(r_ecef, obs_ecef, lat_rad, lon_rad)
|
|
|
|
if prev_el is not None:
|
|
was_above = prev_el >= 0
|
|
is_above = el >= 0
|
|
|
|
if not in_pass and not was_above and is_above:
|
|
aos_time = _bisect_crossing(satrec, obs_ecef, lat_rad, lon_rad, prev_time, current, True)
|
|
jd2, fr2 = _jday_from_dt(aos_time)
|
|
e2, r2, _ = satrec.sgp4(jd2, fr2)
|
|
if e2 == 0:
|
|
gmst2 = _gmst_rad(jd2 + fr2)
|
|
r_ecef2 = _eci_to_ecef(r2, gmst2)
|
|
_, aos_az_v, _ = _ecef_to_altaz(r_ecef2, obs_ecef, lat_rad, lon_rad)
|
|
aos_az = aos_az_v
|
|
in_pass = True
|
|
max_el = el
|
|
max_el_time = current
|
|
|
|
if in_pass:
|
|
if el > max_el:
|
|
max_el = el
|
|
max_el_time = current
|
|
|
|
if was_above and not is_above:
|
|
los_time = _bisect_crossing(satrec, obs_ecef, lat_rad, lon_rad, prev_time, current, False)
|
|
jd3, fr3 = _jday_from_dt(los_time)
|
|
e3, r3, _ = satrec.sgp4(jd3, fr3)
|
|
los_az = None
|
|
if e3 == 0:
|
|
gmst3 = _gmst_rad(jd3 + fr3)
|
|
r_ecef3 = _eci_to_ecef(r3, gmst3)
|
|
_, los_az_v, _ = _ecef_to_altaz(r_ecef3, obs_ecef, lat_rad, lon_rad)
|
|
los_az = los_az_v
|
|
|
|
# ── Phase 2: fine trajectory (5s steps) ──
|
|
path = []
|
|
shadow_entry = None
|
|
t = aos_time
|
|
while t <= los_time + timedelta(seconds=5):
|
|
jd4, fr4 = _jday_from_dt(t)
|
|
e4, r4, _ = satrec.sgp4(jd4, fr4)
|
|
if e4 == 0:
|
|
gmst4 = _gmst_rad(jd4 + fr4)
|
|
r_ecef4 = _eci_to_ecef(r4, gmst4)
|
|
el4, az4, _ = _ecef_to_altaz(r_ecef4, obs_ecef, lat_rad, lon_rad)
|
|
if el4 >= -1:
|
|
sun_km = _sun_eci_km(t)
|
|
shad = _in_shadow(r4, sun_km)
|
|
path.append({
|
|
"t": t.isoformat().replace("+00:00", "Z"),
|
|
"az": round(az4, 2),
|
|
"el": round(el4, 2),
|
|
"shadow": shad,
|
|
})
|
|
t += TRAJ_STEP
|
|
|
|
# Find shadow entry in path
|
|
for pi in range(1, len(path)):
|
|
if not path[pi-1]["shadow"] and path[pi]["shadow"] and path[pi]["el"] >= 0:
|
|
shadow_entry = {
|
|
"t": path[pi]["t"],
|
|
"az": path[pi]["az"],
|
|
"el": path[pi]["el"],
|
|
}
|
|
break
|
|
|
|
# Refine peak
|
|
peak_time = max_el_time
|
|
peak_el = max_el
|
|
peak_az = None
|
|
for pt in path:
|
|
if pt["el"] > peak_el:
|
|
peak_el = pt["el"]
|
|
peak_time = datetime.fromisoformat(pt["t"].replace("Z", "+00:00"))
|
|
peak_az = pt["az"]
|
|
if peak_az is None and path:
|
|
mid_pt = path[len(path)//2]
|
|
peak_az = mid_pt["az"]
|
|
|
|
# Magnitude at peak
|
|
peak_mag = None
|
|
if peak_time:
|
|
jd5, fr5 = _jday_from_dt(peak_time)
|
|
e5, r5, _ = satrec.sgp4(jd5, fr5)
|
|
if e5 == 0:
|
|
sun5 = _sun_eci_km(peak_time)
|
|
peak_mag = _apparent_magnitude(
|
|
sat_item.get("std_mag"),
|
|
sat_item.get("rcs"),
|
|
sat_item.get("rcs_size"),
|
|
r5, sun5, obs_ecef, lat_rad, lon_rad,
|
|
_gmst_rad(jd5 + fr5)
|
|
)
|
|
|
|
duration_s = int(round((los_time - aos_time).total_seconds()))
|
|
local_tz = tz
|
|
|
|
def _iso(dt):
|
|
if dt is None:
|
|
return None
|
|
return dt.isoformat().replace("+00:00", "Z")
|
|
|
|
passes.append({
|
|
"aos_utc": _iso(aos_time),
|
|
"los_utc": _iso(los_time),
|
|
"peak_utc": _iso(peak_time),
|
|
"aos_az": round(aos_az, 2) if aos_az is not None else None,
|
|
"los_az": round(los_az, 2) if los_az is not None else None,
|
|
"peak_az": round(peak_az, 2) if peak_az is not None else None,
|
|
"peak_el": round(peak_el, 2),
|
|
"duration_s": duration_s,
|
|
"path": path,
|
|
"shadow_entry": shadow_entry,
|
|
"peak_mag": peak_mag,
|
|
})
|
|
|
|
in_pass = False
|
|
aos_time = None
|
|
max_el = -999.0
|
|
max_el_time = None
|
|
|
|
prev_el = el
|
|
prev_time = current
|
|
current += SCAN_STEP
|
|
|
|
results.append({
|
|
"id": sat_item.get("id"),
|
|
"norad_cat_id": sat_item.get("norad_cat_id"),
|
|
"name": sat_item.get("object_name", ""),
|
|
"is_favorite": bool(sat_item.get("is_favorite", False)),
|
|
"color": color,
|
|
"passes": passes,
|
|
})
|
|
|
|
return {
|
|
"ok": True,
|
|
"action": "satellite_passes",
|
|
"observer": {"latitude": latitude, "longitude": longitude, "elevation": elevation_m},
|
|
"timezone": timezone_name,
|
|
"window_start": window_start.isoformat().replace("+00:00", "Z"),
|
|
"window_end": window_end.isoformat().replace("+00:00", "Z"),
|
|
"results": results,
|
|
}
|
|
|
|
|
|
def main() -> None:
|
|
if len(sys.argv) < 2:
|
|
fail("Es wurde keine Aktion uebergeben.", extra={"available_actions": ["sun_moon_rise_set", "moon_star_occultations", "moon_phase_details", "satellite_passes", "planet_rise_set", "month_sky_context"]})
|
|
|
|
action = sys.argv[1]
|
|
args = sys.argv[2:]
|
|
|
|
if action == "sun_moon_rise_set":
|
|
result = action_sun_moon_rise_set(args)
|
|
print(json.dumps(result, ensure_ascii=False))
|
|
return
|
|
|
|
if action == "moon_star_occultations":
|
|
result = action_moon_star_occultations(args)
|
|
print(json.dumps(result, ensure_ascii=False))
|
|
return
|
|
|
|
if action == "moon_phase_details":
|
|
result = action_moon_phase_details(args)
|
|
print(json.dumps(result, ensure_ascii=False))
|
|
return
|
|
|
|
if action == "satellite_passes":
|
|
result = action_satellite_passes(args)
|
|
print(json.dumps(result, ensure_ascii=False))
|
|
return
|
|
|
|
if action == "planet_rise_set":
|
|
result = action_planet_rise_set(args)
|
|
print(json.dumps(result, ensure_ascii=False))
|
|
return
|
|
|
|
if action == "month_sky_context":
|
|
result = action_month_sky_context(args)
|
|
print(json.dumps(result, ensure_ascii=False))
|
|
return
|
|
|
|
fail("Unbekannte Aktion.", extra={"action": action, "available_actions": ["sun_moon_rise_set", "moon_star_occultations", "moon_phase_details", "satellite_passes", "planet_rise_set", "month_sky_context"]})
|
|
|
|
|
|
if __name__ == "__main__":
|
|
main()
|