Merge branch 'main' of https://git.daheim.eskimue.de/Eskimue/skyview.astronomiemuseum.de

public/py/observation_suggestions_api.py

@@ -0,0 +1,296 @@
+#!/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)

public/py/test_observation_suggestions_api.py

@@ -0,0 +1,130 @@
+"""Run: python -m unittest discover -s public/py -p test_observation_suggestions_api.py"""
+import unittest
+
+import astronomy
+import observation_suggestions_api as api
+
+
+def payload():
+    return {'location': {'latitude': 50.3774, 'longitude': 11.1899, 'elevation': 640, 'timezone': 'Europe/Berlin'},
+            'filters': {'observation_date': '2026-09-15', 'object_types': ['moon', 'planets', 'galaxies'],
+                        'night_halves': ['first', 'second'], 'object_size_min': '', 'object_size_max': '',
+                        'object_magnitude': '12', 'object_altitude_min': '20', 'object_altitude_max': '90',
+                        'object_azimuth_start': '0', 'object_azimuth_end': '360', 'sun_altitude_max': '-12',
+                        'moon_separation_min': '', 'include_unknown_magnitude': ''},
+            'objects': [{'id': 'm31', 'name': 'M31', 'group': 'galaxies', 'ra_deg': 10.6847,
+                         'dec_deg': 41.269, 'magnitude': 3.44, 'size_arcsec': 10680}]}
+
+
+class ObservationTests(unittest.TestCase):
+    def test_azimuth_wrap_and_full_circle(self):
+        for value in (0, 20, 60, 300, 350):
+            self.assertTrue(api.azimuth_matches(value, 300, 60))
+        self.assertFalse(api.azimuth_matches(180, 300, 60))
+        for value in (0, 90, 180, 359.9):
+            self.assertTrue(api.azimuth_matches(value, 0, 360))
+        self.assertFalse(api.azimuth_matches(0, 90, 270))
+
+    def test_size_exempts_non_deep_sky(self):
+        f = api.validate(payload())[0]
+        f.update(object_size_min=60, object_size_max=600)
+        for group in ('moon', 'planets', 'double_stars', 'variable_stars', 'comets'):
+            self.assertTrue(api.size_matches(None, group, f))
+        self.assertFalse(api.size_matches(None, 'galaxies', f))
+        self.assertFalse(api.size_matches(601, 'galaxies', f))
+        self.assertTrue(api.size_matches(60, 'galaxies', f))
+
+    def test_magnitude_direction_unknown_and_zero(self):
+        f = api.validate(payload())[0]
+        self.assertTrue(api.magnitude_matches(-2, f))
+        self.assertTrue(api.magnitude_matches(12, f))
+        self.assertFalse(api.magnitude_matches(12.1, f))
+        self.assertFalse(api.magnitude_matches(None, f))
+        f['unknown'] = True
+        self.assertTrue(api.magnitude_matches(None, f))
+        f['object_magnitude'] = 0
+        self.assertFalse(api.magnitude_matches(1, f))
+
+    def test_invalid_input_rejected(self):
+        for key, value in [('object_altitude_min', '91'), ('sun_altitude_max', 'NaN'),
+                           ('object_size_min', '-1'), ('observation_date', '2026-02-30'),
+                           ('object_types', []), ('night_halves', [])]:
+            p = payload()
+            p['filters'][key] = value
+            with self.subTest(key=key), self.assertRaises(ValueError):
+                api.validate(p)
+        p = payload()
+        p['filters'].update(object_size_min='600', object_size_max='60')
+        with self.assertRaises(ValueError):
+            api.validate(p)
+
+    def test_frame_matches_library_horizon(self):
+        _, observer, _, start, _ = api.validate(payload())
+        t = api.astro_time(start)
+        rotation = astronomy.Rotation_EQJ_HOR(t, observer).rot
+        for body in (astronomy.Body.Moon, astronomy.Body.Mars, astronomy.Body.Sun):
+            j2000 = astronomy.Equator(body, t, observer, False, True)
+            ofdate = astronomy.Equator(body, t, observer, True, True)
+            expected = astronomy.Horizon(t, observer, ofdate.ra, ofdate.dec, astronomy.Refraction.Airless)
+            altitude, azimuth = api.horizontal(api.xyz(j2000.vec), rotation)
+            self.assertAlmostEqual(altitude, expected.altitude, places=7)
+            self.assertAlmostEqual(azimuth, expected.azimuth, places=7)
+
+    def test_night_bounds_and_half_split(self):
+        f, observer, _, start, end = api.validate(payload())
+        windows = api.night_windows(start, end, observer, f['sun_altitude_max'])
+        self.assertEqual(len(windows), 1)
+        a, b = windows[0]
+        for instant in (a, b):
+            self.assertLessEqual(api.sun_altitude(instant, observer), -12)
+            self.assertAlmostEqual(api.sun_altitude(instant, observer), -12, delta=0.01)
+        first = api.selected_windows(windows, ['first'])[0]
+        second = api.selected_windows(windows, ['second'])[0]
+        self.assertEqual(first[1], second[0])
+        self.assertEqual(first[1] - first[0], second[1] - second[0])
+
+    def test_dst_night_duration(self):
+        for date, hours in [('2026-03-28', 23), ('2026-10-24', 25)]:
+            p = payload()
+            p['filters']['observation_date'] = date
+            _, _, _, a, b = api.validate(p)
+            self.assertEqual((b - a).total_seconds() / 3600, hours)
+
+    def test_polar_day_empty_and_polar_night_explained(self):
+        p = payload()
+        p['location'].update(latitude=80, longitude=15, timezone='Arctic/Longyearbyen')
+        p['filters']['observation_date'] = '2026-06-21'
+        result = api.calculate(p)
+        self.assertEqual(result['rows'], [])
+        self.assertTrue(result['message'])
+        p['filters']['observation_date'] = '2026-12-21'
+        result = api.calculate(p)
+        self.assertTrue(any('12 Uhr' in warning for warning in result['warnings']))
+
+    def test_moon_separation_exemption_and_below_horizon(self):
+        f = api.validate(payload())[0]
+        f['moon_separation_min'] = 60
+        frame = {'rotation': [[1, 0, 0], [0, 1, 0], [0, 0, 1]], 'moon': (0, 0, 1), 'moon_alt': 90}
+        self.assertIsNone(api.sample_matches((0, 0, 1), 5, {'group': 'galaxies'}, frame, f))
+        self.assertIsNotNone(api.sample_matches((0, 0, 1), 5, {'group': 'moon'}, frame, f))
+        frame['moon_alt'] = -5
+        self.assertIsNotNone(api.sample_matches((0, 0, 1), 5, {'group': 'galaxies'}, frame, f))
+
+    def test_real_night_returns_m31_and_matching_positions(self):
+        p = payload()
+        result = api.calculate(p)
+        self.assertTrue(result['ok'])
+        self.assertIn('m31', [row['id'] for row in result['rows']])
+        for row in result['rows']:
+            self.assertGreaterEqual(row['best']['altitude'], 20)
+            self.assertLessEqual(row['best']['magnitude'], 12)
+            self.assertTrue(row['intervals'])
+        p['filters'].update(object_size_max='1')
+        small = api.calculate(p)
+        self.assertNotIn('m31', [row['id'] for row in small['rows']])
+        self.assertEqual({r['id'] for r in result['rows'] if r['group'] in ('moon', 'planets')},
+                         {r['id'] for r in small['rows'] if r['group'] in ('moon', 'planets')})
+
+
+if __name__ == '__main__':
+    unittest.main()