Jyotisha/scripts/birth_sky.py

#!/usr/bin/env python3
"""The real sky over a birth place at the birth instant (「那一刻的天空」).

Pure computation for the birth-sky cover (TASK-birth-sky-cover-20260928 T1).
Thinly registered as ``POST /api/birth_sky`` from ``jyotish_api_server.py``.

Conventions
-----------
* Planets: tropical apparent ecliptic positions from ``swe.calc_ut`` (default
  flags), converted with ``swe.azalt(..., swe.ECL2HOR, ...)``.
* Stars: Yale Bright Star Catalogue (J2000, V <= 4.0, ``data/bright_stars.json``),
  precessed to the birth epoch with IAU 1976 (Lieske) precession, then
  ``swe.azalt(..., swe.EQU2HOR, ...)``. Proper motion and nutation are not
  applied: at this drawing's scale they are far below one pixel.
* Altitude is the **geometric (true) altitude**: no atmospheric refraction
  (``atpress = 0``), so the golden values are reproducible anywhere.
* Azimuth is returned **from north, increasing eastward**. Swiss Ephemeris
  ``azalt`` measures azimuth from **south** (westward), so we add 180°.

The response never echoes the input date, time, or coordinates.
"""
from __future__ import annotations

import json
import math
from datetime import datetime, timedelta
from functools import lru_cache
from pathlib import Path
from typing import Any

try:
    import swisseph as swe
except ImportError:  # pragma: no cover
    swe = None

STAR_FILE = Path(__file__).resolve().parent / "data" / "bright_stars.json"
PLANETS = (
    ("sun", "SUN"),
    ("moon", "MOON"),
    ("mercury", "MERCURY"),
    ("venus", "VENUS"),
    ("mars", "MARS"),
    ("jupiter", "JUPITER"),
    ("saturn", "SATURN"),
)
TWILIGHT_LIMIT = -18.0
J2000 = 2451545.0


class BirthSkyError(ValueError):
    """Bad input. The HTTP layer maps it to 400."""


@lru_cache(maxsize=1)
def _catalog() -> dict[str, Any]:
    return json.loads(STAR_FILE.read_text(encoding="utf-8"))


def _number(body: dict[str, Any], key: str, low: float, high: float) -> float:
    value = body.get(key)
    if isinstance(value, bool) or value is None or value == "":
        raise BirthSkyError(f"{key} is required")
    try:
        number = float(value)
    except (TypeError, ValueError) as exc:
        raise BirthSkyError(f"{key} must be a number") from exc
    if not math.isfinite(number) or number < low or number > high:
        raise BirthSkyError(f"{key} out of range")
    return number


def _instant(date: str, time: str, tz_offset_hours: float) -> datetime:
    if not isinstance(date, str) or not isinstance(time, str):
        raise BirthSkyError("date and time are required")
    clock = time.strip()[:8]
    for pattern in ("%H:%M:%S", "%H:%M"):
        try:
            local = datetime.strptime(f"{date.strip()[:10]} {clock}", f"%Y-%m-%d {pattern}")
            break
        except ValueError:
            continue
    else:
        raise BirthSkyError("date must be YYYY-MM-DD and time HH:MM")
    return local - timedelta(hours=tz_offset_hours)


def _julian_day_ut(utc: datetime) -> float:
    hours = utc.hour + utc.minute / 60 + utc.second / 3600
    return swe.julday(utc.year, utc.month, utc.day, hours)


def precess_j2000(ra_deg: float, dec_deg: float, jd: float) -> tuple[float, float]:
    """IAU 1976 precession of a J2000 mean position to the mean equator of ``jd``."""
    t = (jd - J2000) / 36525.0
    arcsec = math.pi / (180 * 3600)
    zeta = (2306.2181 * t + 0.30188 * t * t + 0.017998 * t ** 3) * arcsec
    z = (2306.2181 * t + 1.09468 * t * t + 0.018203 * t ** 3) * arcsec
    theta = (2004.3109 * t - 0.42665 * t * t - 0.041833 * t ** 3) * arcsec
    ra0 = math.radians(ra_deg)
    dec0 = math.radians(dec_deg)
    a = math.cos(dec0) * math.sin(ra0 + zeta)
    b = math.cos(theta) * math.cos(dec0) * math.cos(ra0 + zeta) - math.sin(theta) * math.sin(dec0)
    c = math.sin(theta) * math.cos(dec0) * math.cos(ra0 + zeta) + math.cos(theta) * math.sin(dec0)
    ra = (math.degrees(math.atan2(a, b) + z)) % 360.0
    dec = math.degrees(math.asin(max(-1.0, min(1.0, c))))
    return ra, dec


def _horizontal(jd: float, flag: int, geopos: tuple[float, float, float], xin: tuple[float, float, float]) -> tuple[float, float]:
    # swe.azalt returns (azimuth from SOUTH, true altitude, apparent altitude).
    azimuth_from_south, true_altitude, _apparent = swe.azalt(jd, flag, geopos, 0, 10, xin)
    azimuth_from_north = (azimuth_from_south + 180.0) % 360.0
    return true_altitude, azimuth_from_north


def _phase(sun_altitude: float, sun_hour_angle: float) -> str:
    if sun_altitude > 0:
        return "day"
    if sun_altitude >= TWILIGHT_LIMIT:
        # Hour angle 0..180 = the Sun is west of the meridian: local apparent afternoon.
        return "twilight_evening" if 0 <= sun_hour_angle < 180 else "twilight_morning"
    return "night"


def compute_birth_sky(date: str, time: str, tz_offset_hours: float, lat: float, lon: float) -> dict[str, Any]:
    """Altitude / azimuth of the seven visible grahas and bright stars at one instant."""
    if swe is None:  # pragma: no cover
        raise RuntimeError("pyswisseph is not installed")
    jd = _julian_day_ut(_instant(date, time, tz_offset_hours))
    geopos = (float(lon), float(lat), 0.0)

    planets: list[dict[str, Any]] = []
    sun_altitude = 0.0
    sun_ra = 0.0
    for key, attr in PLANETS:
        body = getattr(swe, attr)
        ecliptic, _ = swe.calc_ut(jd, body)
        altitude, azimuth = _horizontal(jd, swe.ECL2HOR, geopos, (ecliptic[0], ecliptic[1], ecliptic[2]))
        planets.append({"id": key, "altitude": round(altitude, 2), "azimuth": round(azimuth, 2)})
        if key == "sun":
            sun_altitude = altitude
            equatorial, _ = swe.calc_ut(jd, body, swe.FLG_EQUATORIAL)
            sun_ra = equatorial[0]

    local_sidereal_deg = (swe.sidtime(jd) * 15.0 + lon) % 360.0
    sun_hour_angle = (local_sidereal_deg - sun_ra) % 360.0

    catalog = _catalog()
    stars: list[dict[str, Any]] = []
    for hr, ra, dec, vmag in catalog["stars"]:
        ra_date, dec_date = precess_j2000(ra, dec, jd)
        altitude, azimuth = _horizontal(jd, swe.EQU2HOR, geopos, (ra_date, dec_date, 1.0))
        if altitude <= 0:
            continue
        stars.append({"hr": hr, "altitude": round(altitude, 2), "azimuth": round(azimuth, 2), "magnitude": vmag})

    return {
        "phase": _phase(sun_altitude, sun_hour_angle),
        "sun_altitude": round(sun_altitude, 2),
        "planets": planets,
        "stars": stars,
        "lines": [{"id": line["id"], "strokes": line["strokes"]} for line in catalog["lines"]],
        "conventions": {
            "altitude": "geometric, no refraction",
            "azimuth": "degrees from north, eastward",
            "planets": "swe.calc_ut tropical apparent + swe.azalt ECL2HOR",
            "stars": "Yale BSC5 V<=4.0, IAU 1976 precession + swe.azalt EQU2HOR",
        },
    }


def compute_birth_sky_from_request(body: Any) -> dict[str, Any]:
    if not isinstance(body, dict):
        raise BirthSkyError("body must be an object")
    return compute_birth_sky(
        str(body.get("date") or ""),
        str(body.get("time") or ""),
        _number(body, "timezone_offset", -14, 14),
        _number(body, "latitude", -90, 90),
        _number(body, "longitude", -180, 180),
    )