297 lines
13 KiB
Python
297 lines
13 KiB
Python
#!/usr/bin/env python3
|
|
"""Independent observation planner. JSON file in, JSON out; no dependency on api.py."""
|
|
import json
|
|
import math
|
|
import sys
|
|
from datetime import datetime, timedelta, timezone
|
|
from zoneinfo import ZoneInfo
|
|
|
|
import astronomy
|
|
import comets
|
|
|
|
UTC = timezone.utc
|
|
STEP = timedelta(minutes=5)
|
|
GROUPS = {
|
|
'open_clusters': ('OCl',), 'globular_clusters': ('GCl',),
|
|
'galaxies': ('G', 'GPair', 'GTrpl', 'GGroup'),
|
|
'nebulae': ('Neb', 'EmN', 'HII', 'RfN', 'SNR', 'Cl+N'),
|
|
'planetary_nebulae': ('PN',), 'double_stars': ('**',),
|
|
}
|
|
LABELS = {
|
|
'open_clusters': 'Offener Sternhaufen', 'globular_clusters': 'Kugelsternhaufen',
|
|
'galaxies': 'Galaxie / Galaxiengruppe', 'nebulae': 'Nebel',
|
|
'planetary_nebulae': 'Planetarischer Nebel', 'double_stars': 'Doppelstern',
|
|
'variable_stars': 'Veränderlicher Stern', 'comets': 'Komet',
|
|
'moon': 'Mond', 'planets': 'Planet',
|
|
}
|
|
BODIES = {'Moon': 'Mond', 'Mercury': 'Merkur', 'Venus': 'Venus', 'Mars': 'Mars',
|
|
'Jupiter': 'Jupiter', 'Saturn': 'Saturn', 'Uranus': 'Uranus', 'Neptune': 'Neptun'}
|
|
|
|
|
|
def number(value, label, minimum=None, maximum=None, optional=False):
|
|
if optional and (value is None or value == ''):
|
|
return None
|
|
try:
|
|
result = float(value)
|
|
except (TypeError, ValueError):
|
|
raise ValueError(f'{label}: Bitte eine Zahl eingeben.')
|
|
if not math.isfinite(result) or (minimum is not None and result < minimum) or (maximum is not None and result > maximum):
|
|
raise ValueError(f'{label}: Wert außerhalb des gültigen Bereichs.')
|
|
return result
|
|
|
|
|
|
def validate(payload):
|
|
raw = payload['filters']
|
|
result = {}
|
|
for key, low, high, optional in (
|
|
('object_size_min', 0, 1296000, True), ('object_size_max', 0, 1296000, True),
|
|
('object_magnitude', -40, 40, True), ('object_altitude_min', 0, 90, False),
|
|
('object_altitude_max', 0, 90, False), ('object_azimuth_start', 0, 360, False),
|
|
('object_azimuth_end', 0, 360, False), ('sun_altitude_max', -90, 0, False),
|
|
('moon_separation_min', 0, 180, True),
|
|
):
|
|
result[key] = number(raw.get(key), key, low, high, optional)
|
|
for prefix in ('object_size', 'object_altitude'):
|
|
lo, hi = result[prefix + '_min'], result[prefix + '_max']
|
|
if lo is not None and hi is not None and lo > hi:
|
|
raise ValueError('Das Minimum darf nicht größer als das Maximum sein.')
|
|
result['types'] = raw.get('object_types', [])
|
|
result['halves'] = raw.get('night_halves', [])
|
|
if not isinstance(result['types'], list) or not result['types'] or any(t not in LABELS for t in result['types']):
|
|
raise ValueError('Bitte mindestens einen gültigen Objekttyp auswählen.')
|
|
if not isinstance(result['halves'], list) or not result['halves'] or any(h not in ('first', 'second') for h in result['halves']):
|
|
raise ValueError('Bitte mindestens eine Nachthälfte auswählen.')
|
|
result['unknown'] = raw.get('include_unknown_magnitude') in (True, 1, '1')
|
|
date = datetime.strptime(raw.get('observation_date', ''), '%Y-%m-%d')
|
|
if not 1900 <= date.year <= 2100:
|
|
raise ValueError('Bitte ein Datum zwischen 1900 und 2100 wählen.')
|
|
loc = payload['location']
|
|
zone = ZoneInfo(loc['timezone'])
|
|
observer = astronomy.Observer(number(loc['latitude'], 'Breitengrad', -90, 90),
|
|
number(loc['longitude'], 'Längengrad', -180, 180),
|
|
number(loc.get('elevation') or 0, 'Standorthöhe', -500, 10000))
|
|
start = date.replace(hour=12, tzinfo=zone).astimezone(UTC)
|
|
end = (date + timedelta(days=1)).replace(hour=12, tzinfo=zone).astimezone(UTC)
|
|
return result, observer, zone, start, end
|
|
|
|
|
|
def astro_time(dt):
|
|
return comets.dt_to_time(dt)
|
|
|
|
|
|
def unit(vector):
|
|
length = math.sqrt(sum(v * v for v in vector))
|
|
return tuple(v / length for v in vector)
|
|
|
|
|
|
def xyz(vector):
|
|
return vector.x, vector.y, vector.z
|
|
|
|
|
|
def horizontal(vector, rotation):
|
|
x, y, z = (sum(rotation[j][i] * vector[j] for j in range(3)) for i in range(3))
|
|
return math.degrees(math.atan2(z, math.hypot(x, y))), math.degrees(math.atan2(-y, x)) % 360
|
|
|
|
|
|
def azimuth_matches(azimuth, start, end):
|
|
if abs(end - start) == 360:
|
|
return True
|
|
start, end = start % 360, end % 360
|
|
return start <= azimuth <= end if start <= end else azimuth >= start or azimuth <= end
|
|
|
|
|
|
def magnitude_matches(mag, filters):
|
|
limit = filters['object_magnitude']
|
|
return limit is None or (filters['unknown'] if mag is None else mag <= limit)
|
|
|
|
|
|
def size_matches(size, group, filters):
|
|
if group not in GROUPS or group == 'double_stars':
|
|
return True
|
|
lo, hi = filters['object_size_min'], filters['object_size_max']
|
|
if lo is None and hi is None:
|
|
return True
|
|
return size is not None and (lo is None or size >= lo) and (hi is None or size <= hi)
|
|
|
|
|
|
def sun_altitude(dt, observer):
|
|
t = astro_time(dt)
|
|
eq = astronomy.Equator(astronomy.Body.Sun, t, observer, True, True)
|
|
return astronomy.Horizon(t, observer, eq.ra, eq.dec, astronomy.Refraction.Airless).altitude
|
|
|
|
|
|
def night_windows(start, end, observer, limit):
|
|
"""Bracket solar threshold crossings in UTC, refine to one second."""
|
|
windows = []
|
|
previous = start
|
|
inside = sun_altitude(start, observer) <= limit
|
|
opened = start if inside else None
|
|
while previous < end:
|
|
current = min(previous + STEP, end)
|
|
next_inside = sun_altitude(current, observer) <= limit
|
|
if inside != next_inside:
|
|
lo, hi = previous, current
|
|
while (hi - lo).total_seconds() > 1:
|
|
mid = lo + (hi - lo) / 2
|
|
if (sun_altitude(mid, observer) <= limit) == inside:
|
|
lo = mid
|
|
else:
|
|
hi = mid
|
|
crossing = hi if next_inside else lo
|
|
if next_inside:
|
|
opened = crossing
|
|
else:
|
|
windows.append((opened, crossing))
|
|
opened = None
|
|
inside, previous = next_inside, current
|
|
if opened is not None:
|
|
windows.append((opened, end))
|
|
return windows
|
|
|
|
|
|
def selected_windows(windows, halves):
|
|
selected = []
|
|
for start, end in windows:
|
|
midpoint = start + (end - start) / 2
|
|
if 'first' in halves and 'second' in halves:
|
|
selected.append((start, end))
|
|
elif 'first' in halves:
|
|
selected.append((start, midpoint))
|
|
else:
|
|
selected.append((midpoint, end))
|
|
return selected
|
|
|
|
|
|
def frames_for(windows, observer):
|
|
frames = []
|
|
for segment, (start, end) in enumerate(windows):
|
|
dt = start
|
|
while True:
|
|
t = astro_time(dt)
|
|
rotation = astronomy.Rotation_EQJ_HOR(t, observer).rot
|
|
moon = astronomy.Equator(astronomy.Body.Moon, t, observer, False, True)
|
|
moon_vector = unit(xyz(moon.vec))
|
|
frames.append({'dt': dt, 'time': t, 'rotation': rotation, 'segment': segment,
|
|
'moon': moon_vector, 'moon_alt': horizontal(moon_vector, rotation)[0],
|
|
'earth': xyz(astronomy.HelioVector(astronomy.Body.Earth, t)),
|
|
'observer': xyz(astronomy.ObserverVector(t, observer, False))})
|
|
if dt == end:
|
|
break
|
|
dt = min(dt + STEP, end)
|
|
return frames
|
|
|
|
|
|
def fixed_vector(obj):
|
|
ra = math.radians(number(obj['ra_deg'], 'Rektaszension', 0, 360))
|
|
dec = math.radians(number(obj['dec_deg'], 'Deklination', -90, 90))
|
|
return math.cos(dec) * math.cos(ra), math.cos(dec) * math.sin(ra), math.sin(dec)
|
|
|
|
|
|
def comet_state(obj, frame):
|
|
orbit = obj['orbit']
|
|
# Iterated light time; shared two-body orbit helper, no api.py import.
|
|
emitted = frame['dt']
|
|
for _ in range(3):
|
|
helio = comets.comet_heliocentric_vector(orbit, emitted)
|
|
if helio is None:
|
|
raise ValueError('Unvollständige Kometenbahn')
|
|
geo = tuple(helio[i] - frame['earth'][i] for i in range(3))
|
|
distance = math.sqrt(sum(v * v for v in geo))
|
|
emitted = frame['dt'] - timedelta(days=distance / astronomy.C_AUDAY)
|
|
topocentric = tuple(geo[i] - frame['observer'][i] for i in range(3))
|
|
mag = comets.estimate_magnitude(comets.parse_float(orbit.get('absolute_magnitude_h')),
|
|
comets.parse_float(orbit.get('slope_parameter_g')),
|
|
math.sqrt(sum(v * v for v in helio)), distance)
|
|
return unit(topocentric), mag
|
|
|
|
|
|
def sample_matches(vector, mag, obj, frame, filters):
|
|
altitude, azimuth = horizontal(vector, frame['rotation'])
|
|
if not filters['object_altitude_min'] <= altitude <= filters['object_altitude_max']:
|
|
return None
|
|
if not azimuth_matches(azimuth, filters['object_azimuth_start'], filters['object_azimuth_end']):
|
|
return None
|
|
if not magnitude_matches(mag, filters):
|
|
return None
|
|
separation = math.degrees(math.acos(max(-1, min(1, sum(a * b for a, b in zip(vector, frame['moon']))))))
|
|
if obj['group'] != 'moon' and frame['moon_alt'] > 0 and filters['moon_separation_min'] is not None and separation < filters['moon_separation_min']:
|
|
return None
|
|
return {'altitude': altitude, 'azimuth': azimuth, 'magnitude': mag, 'moon_separation': separation}
|
|
|
|
|
|
def calculate(payload):
|
|
filters, observer, zone, start, end = validate(payload)
|
|
windows = night_windows(start, end, observer, filters['sun_altitude_max'])
|
|
selected = selected_windows(windows, filters['halves'])
|
|
frames = frames_for(selected, observer)
|
|
stamp = lambda dt: dt.astimezone(zone).isoformat(timespec='seconds')
|
|
warnings = []
|
|
if windows and (windows[0][0] == start or windows[-1][1] == end):
|
|
warnings.append('Die Dunkelphase reicht über die Nacht hinaus. Berechnet wird von 12 Uhr bis 12 Uhr am Folgetag in Standortzeit.')
|
|
objects = list(payload.get('objects', []))
|
|
if len(objects) > 30000:
|
|
raise ValueError('Zu viele Katalogobjekte für eine Anfrage.')
|
|
for body, name in BODIES.items():
|
|
group = 'moon' if body == 'Moon' else 'planets'
|
|
if group in filters['types']:
|
|
objects.append({'id': body, 'name': name, 'group': group, 'body': body, 'size_arcsec': None})
|
|
results, invalid = [], 0
|
|
for obj in objects if frames else []:
|
|
group = obj['group']
|
|
if group not in filters['types'] or not size_matches(obj.get('size_arcsec'), group, filters):
|
|
continue
|
|
moving = 'body' in obj or group == 'comets'
|
|
if not moving and not magnitude_matches(obj.get('magnitude'), filters):
|
|
continue
|
|
try:
|
|
vector = None if moving else fixed_vector(obj)
|
|
intervals, opened, last, best = [], None, None, None
|
|
for frame in frames:
|
|
mag = obj.get('magnitude')
|
|
if 'body' in obj:
|
|
body = getattr(astronomy.Body, obj['body'])
|
|
vector = unit(xyz(astronomy.Equator(body, frame['time'], observer, False, True).vec))
|
|
mag = astronomy.Illumination(body, frame['time']).mag
|
|
elif group == 'comets':
|
|
vector, mag = comet_state(obj, frame)
|
|
match = sample_matches(vector, mag, obj, frame, filters)
|
|
if last is not None and (match is None or last['segment'] != frame['segment']):
|
|
intervals.append({'start': stamp(opened['dt']), 'end': stamp(last['dt'])})
|
|
opened, last = None, None
|
|
if match is not None:
|
|
if opened is None:
|
|
opened = frame
|
|
last = frame
|
|
if best is None or match['altitude'] > best['altitude']:
|
|
best = dict(match, time=stamp(frame['dt']))
|
|
if last is not None:
|
|
intervals.append({'start': stamp(opened['dt']), 'end': stamp(last['dt'])})
|
|
if best:
|
|
results.append({'id': obj['id'], 'name': obj['name'], 'group': group, 'type': LABELS[group],
|
|
'size_arcsec': obj.get('size_arcsec'), 'best': best, 'intervals': intervals,
|
|
'magnitude_note': 'Modellschätzung' if moving else 'Katalogwert',
|
|
'note': obj.get('note', '')})
|
|
except (ValueError, TypeError, ZeroDivisionError, OverflowError):
|
|
invalid += 1
|
|
if 'variable_stars' in filters['types']:
|
|
warnings.append('Veränderliche: Sichtbarkeit und Kataloghelligkeit, keine aktuelle Helligkeit oder Minima-/Maxima-Vorhersage.')
|
|
if 'comets' in filters['types']:
|
|
warnings.append('Kometen: Näherung aus gespeicherten Bahnelementen; Helligkeiten sind unsichere Modellschätzungen.')
|
|
if invalid:
|
|
warnings.append(f'{invalid} Objekte wegen unbrauchbarer Koordinaten oder Bahnelemente ausgelassen.')
|
|
results.sort(key=lambda r: (r['best']['magnitude'] is None, -r['best']['altitude'], r['name']))
|
|
return {'ok': True, 'rows': results, 'count': len(results), 'timezone': str(zone),
|
|
'windows': [{'start': stamp(a), 'end': stamp(b)} for a, b in selected],
|
|
'step_minutes': 5, 'warnings': warnings,
|
|
'message': 'Die Sonne erreicht in dieser Nacht die gewählte Höhe nicht.' if not frames else ''}
|
|
|
|
|
|
if __name__ == '__main__':
|
|
try:
|
|
with open(sys.argv[1], encoding='utf-8-sig') as source:
|
|
response = calculate(json.load(source))
|
|
print(json.dumps(response, ensure_ascii=True, allow_nan=False))
|
|
except (ValueError, KeyError, TypeError) as exc:
|
|
print(json.dumps({'ok': False, 'error': str(exc)}, ensure_ascii=True))
|
|
sys.exit(1)
|