# SPDX-License-Identifier: MIT

# astronomy — plan observations from sites, stars and bookings, and calculate
# time, the sky's rotation, the Sun and orbits.
#
# Two halves that share nothing but the subject.
#
# STRUCTURE. You supply sites, stars, nights and bookings as facts.
#
# Sites (`astronomy.Site`) by latitude, longitude and the faintest magnitude
# their instrument reaches, and stars (`astronomy.Star`) by right ascension,
# declination and magnitude. The rules conclude, for each site, the stars that
# never set, the stars that never rise, and the stars it can see:
#
#     import astronomy
#     astronomy.site("oslo", 59.91, 10.75, 6.0)
#     astronomy.star("polaris", 37.95, 89.26, 1.98)
#     astronomy.star("canopus", 95.99, -52.70, -0.74)
#     print(astronomy.circumpolar("oslo"), astronomy.hidden("oslo"))
#
# Nights (`astronomy.Night`), each the dark time at a site, and bookings
# (`astronomy.Booking`), each an instrument booked at a site for a star, are
# written with the windows they hold for. `up()` writes the window a star is
# above a site's horizon on a day. The rules conclude that a booking is
# observable for exactly the stretch where the three overlap, and report that
# stretch as the conclusion's own window:
#
#     astronomy.star("vega", 279.23, 38.78, 0.03)
#     astronomy.night("n1", "oslo", datetime("2026-11-08T16:00:00Z"), datetime("2026-11-09T06:00:00Z"))
#     astronomy.up("vega", "oslo", datetime("2026-11-08T16:00:00Z"))
#     astronomy.book("b1", "oslo", "vega", datetime("2026-11-08T18:00:00Z"), datetime("2026-11-08T22:00:00Z"))
#     print(astronomy.window("b1"))
#
# Each conclusion answers `why`, and withdraws itself when a fact it rests on
# is corrected. Where two sources disagree about a fact, both values are
# reasoned about and `unsettled()` names what they disagree about. No read
# refuses but `up()` and `window()`, which read one value and refuse a cell
# holding two.
#
# CALCULATION. Value scripts that return a value and write nothing: Julian
# dates, Greenwich and local sidereal time, a star's altitude and its rising
# and setting, the Sun's position and the length of the day, Kepler's
# equation by Newton's method, and the two-body problem by fourth-order
# Runge–Kutta. Angles are in degrees and times in UT. Rising and setting are
# for the geometric horizon unless an altitude is given: refraction lifts a
# body about half a degree, which `day_length` includes for the Sun.

import math

# ── Sites and stars as facts ─────────────────────────────────────────────────

class Site:
    # A site, by a name you choose, at latitude `lat` and longitude `lon` in
    # degrees, north and east positive, and `limit`, the faintest magnitude
    # its instrument reaches. Computed from the latitude: for a northern
    # site, `circumpolar_above`, 90° − φ, the declination above which a star
    # never sets, and `hidden_below`, φ − 90°, the declination below which
    # one never rises; for a southern site, `circumpolar_below`, −90° − φ,
    # and `hidden_above`, 90° + φ, the same limits for the other pole.
    lat: real
    lon: real
    limit: real
    circumpolar_above: real = 90 - lat
    hidden_below: real = lat - 90
    circumpolar_below: real = -90 - lat
    hidden_above: real = 90 + lat

class Star:
    # A star, by a name you choose: right ascension `ra` and declination
    # `dec` in degrees, and its magnitude `mag`.
    ra: real
    dec: real
    mag: real

class Circumpolar:
    # The relation `Circumpolar.star.site`, True for a star that never sets
    # at a site.
    pass

class Hidden:
    # The relation `Hidden.star.site`, True for a star that never rises at a
    # site.
    pass

class Visible:
    # The relation `Visible.star.site`, True for a star that rises at a site
    # and is no fainter than its limit.
    pass

# ── Nights, bookings and the sky as facts ────────────────────────────────────

class Night:
    # The dark time at the site named `site`, by a name you choose, written
    # with the window it lasts; `opens` and `closes` record that window as
    # facts of their own, valid at all times.
    site: string
    opens: datetime
    closes: datetime

class Up:
    # A star above a site's horizon, written by `up()` with the window from
    # its rising to its setting, and named `<star>_<site>_<yyyymmdd>`;
    # `rises` and `sets` record that window as facts valid at all times.
    star: string
    site: string
    rises: datetime
    sets: datetime

class Booking:
    # An instrument booked at the site named `site` for the star named
    # `star`, by a name you choose, written with the window it is booked
    # for, which `starts` and `ends` record as facts valid at all times. The
    # rules conclude `observable` where the booking, a night at its site and
    # its star's time above that site's horizon overlap, and the conclusion
    # holds for that overlap alone.
    site: string
    star: string
    starts: datetime
    ends: datetime
    observable: bool

# ── Standing rules: the sky at a site ────────────────────────────────────────

@rule
def circumpolar_north():
    for s in Star, t in Site:
        if t.lat >= 0 and s.dec > t.circumpolar_above:
            Circumpolar.s.t = True

@rule
def circumpolar_south():
    for s in Star, t in Site:
        if t.lat < 0 and s.dec < t.circumpolar_below:
            Circumpolar.s.t = True

@rule
def hidden_north():
    for s in Star, t in Site:
        if t.lat >= 0 and s.dec < t.hidden_below:
            Hidden.s.t = True

@rule
def hidden_south():
    for s in Star, t in Site:
        if t.lat < 0 and s.dec > t.hidden_above:
            Hidden.s.t = True

@rule
def visible_north():
    for s in Star, t in Site:
        if t.lat >= 0 and s.dec >= t.hidden_below and s.mag <= t.limit:
            Visible.s.t = True

@rule
def visible_south():
    for s in Star, t in Site:
        if t.lat < 0 and s.dec <= t.hidden_above and s.mag <= t.limit:
            Visible.s.t = True

# ── Standing rules: observing ────────────────────────────────────────────────

@rule
def observable():
    for b in Booking, n in Night, u in Up:
        if n.site == b.site and u.site == b.site and u.star == b.star:
            b.observable = True

# ── Writing facts ────────────────────────────────────────────────────────────

def site(name: string, latitude: real, longitude: real, faintest: real):
    update(Site(name).lat, latitude, "sites")
    update(Site(name).lon, longitude, "sites")
    update(Site(name).limit, faintest, "sites")

def star(name: string, right_ascension: real, declination: real, magnitude: real):
    update(Star(name).ra, right_ascension, "catalogue")
    update(Star(name).dec, declination, "catalogue")
    update(Star(name).mag, magnitude, "catalogue")

def night(name: string, where: string, start: datetime, end: datetime):
    # The dark time at a site, from start to end.
    update(Night(name).site, where, "nights", start, end)
    update(Night(name).opens, start, "nights")
    update(Night(name).closes, end, "nights")

def book(name: string, where: string, target: string, start: datetime, end: datetime):
    # An instrument at a site booked for a star, from start to end.
    update(Booking(name).site, where, "bookings", start, end)
    update(Booking(name).star, target, "bookings", start, end)
    update(Booking(name).starts, start, "bookings")
    update(Booking(name).ends, end, "bookings")

def number(c: active_claim) -> real:
    # The number a claim holds, an int read as a real; 0.0 for anything else.
    match c.value:
        case real as x:
            return x
        case int as k:
            return k * 1.0
        case default:
            return 0.0

def up(target: string, where: string, start: datetime):
    # Write the window a star is above a site's horizon on its first passage
    # that has not ended by `start`: from its rising to its setting, or the
    # day from `start` for a star that never sets. Nothing is written for a
    # star that never rises. Pass a night's start to get that night's
    # passage, which may have begun before it.
    r = current(Star.target.ra)
    d = current(Star.target.dec)
    p = current(Site.where.lat)
    q = current(Site.where.lon)
    h = semi_arc(number(d), number(p), 0.0)
    name = target + "_" + where + "_" + str(start)[0:10].replace("-", "")
    if h == 180:
        update(Up(name).star, target, "ephemeris", start, start + duration("P1D"))
        update(Up(name).site, where, "ephemeris", start, start + duration("P1D"))
        update(Up(name).rises, start, "ephemeris")
        update(Up(name).sets, start + duration("P1D"), "ephemeris")
    elif h > 0:
        w = duration("P1D") * (h / 360.98564736629)
        t0 = start - w
        m = ((math.fmod(number(r), 360.0) - lst(jd_of(t0), number(q))) % 360) / 360.98564736629
        transit = t0 + duration("P1D") * m
        update(Up(name).star, target, "ephemeris", transit - w, transit + w)
        update(Up(name).site, where, "ephemeris", transit - w, transit + w)
        update(Up(name).rises, transit - w, "ephemeris")
        update(Up(name).sets, transit + w, "ephemeris")

# ── Reading what the rules concluded ─────────────────────────────────────────
#
# The readers that return names or a range read every value in force, so a
# fact two sources disagree about is reported with both values. `up` reads a
# star's position and a site's, and `window` reads a booking's start and end,
# each with `current()`, which refuses a cell holding two different values
# rather than choose one. `observable_at` reads a truth, and a disputed booking
# is observable if any value of it is.

def stars_in(relation: string, where: string) -> list<string>:
    # The stars in "circumpolar", "hidden" or "visible" for a site, sorted.
    if relation == "circumpolar":
        vs = list_values(Circumpolar.where)
    elif relation == "hidden":
        vs = list_values(Hidden.where)
    else:
        vs = list_values(Visible.where)
    out = []
    for v in vs:
        parts = v.path.split(".")
        s = parts[len(parts) - 2]
        if s not in out:
            out = out + [s]
    return sorted(out)

def circumpolar(where: string) -> list<string>:
    # The stars that never set at a site, sorted.
    return stars_in("circumpolar", where)

def hidden(where: string) -> list<string>:
    # The stars that never rise at a site, sorted.
    return stars_in("hidden", where)

def visible(where: string) -> list<string>:
    # The stars that rise at a site and are no fainter than its limit,
    # sorted.
    return stars_in("visible", where)

def observable_at(t: datetime) -> list<string>:
    # The bookings observable at an instant, sorted.
    out = []
    vs = list_values(Booking.site, t)
    for v in vs:
        parts = v.path.split(".")
        b = parts[len(parts) - 2]
        if b not in out:
            c = current(Booking.b.observable, t)
            match c:
                case empty:
                    pass
                case default:
                    out = out + [b]
    return sorted(out)

def instant_of(c: active_claim) -> list<datetime>:
    # The instant a claim holds, as a list of one, or empty.
    match c.value:
        case datetime as t:
            return [t]
        case default:
            return []

def window(booking: string) -> list<string>:
    # The stretch a booking is observable for, as its first and last
    # instants, or empty where it never is. The stretch opens where the last
    # of its premises does: the booking itself, a night, or a rising. So the
    # conclusion is read at each of those instants inside the booking, in
    # order, and the first that holds reports its window.
    s0 = current(Booking.booking.starts)
    s1 = current(Booking.booking.ends)
    b0 = instant_of(s0)
    b1 = instant_of(s1)
    if len(b0) == 0 or len(b1) == 0:
        return []
    at = [b0[0]]
    vs = list_values(Night.opens)
    for v in vs:
        match v.value:
            case datetime as t:
                at = at + [t]
            case default:
                pass
    vs = list_values(Up.rises)
    for v in vs:
        match v.value:
            case datetime as t:
                at = at + [t]
            case default:
                pass
    for t in sorted(at):
        if t >= b0[0] and t <= b1[0]:
            c = current(Booking.booking.observable, t)
            match c:
                case empty:
                    pass
                case default:
                    out = []
                    match c.valid_from:
                        case datetime as f:
                            out = out + [str(f)]
                        case default:
                            out = out + [""]
                    match c.valid_to:
                        case datetime as e:
                            out = out + [str(e)]
                        case default:
                            out = out + [""]
                    return out
    return []

def unsettled() -> list<string>:
    # Each site and star that holds two different values for a fact you
    # write, as "site x lat", "star x dec" and so on, sorted.
    out = []
    for f in ["lat", "lon", "limit", "ra", "dec", "mag"]:
        if f == "lat":
            vs = list_values(Site.lat)
            label = "site "
        elif f == "lon":
            vs = list_values(Site.lon)
            label = "site "
        elif f == "limit":
            vs = list_values(Site.limit)
            label = "site "
        elif f == "ra":
            vs = list_values(Star.ra)
            label = "star "
        elif f == "dec":
            vs = list_values(Star.dec)
            label = "star "
        else:
            vs = list_values(Star.mag)
            label = "star "
        seen = {}
        for v in vs:
            parts = v.path.split(".")
            s = label + parts[len(parts) - 2] + " " + f
            x = str(v.value)
            if s in seen and seen[s] != x and s not in out:
                out = out + [s]
            seen[s] = x
    return sorted(out)

# ── Calculation: time ────────────────────────────────────────────────────────

def jd(year: int, month: int, day: real) -> real:
    # The Julian date of a calendar date, the day carrying its fraction, by
    # Meeus's formula 7.1: in the Julian calendar before 1582 October 15 and
    # the Gregorian from then, as astronomers date events.
    y = year
    m = month
    if m <= 2:
        y = y - 1
        m = m + 12
    a = y // 100
    b = 2 - a + a // 4
    if year < 1582 or (year == 1582 and (month < 10 or (month == 10 and day < 15))):
        b = 0
    return math.floor(365.25 * (y + 4716)) + math.floor(30.6001 * (m + 1)) + day + b - 1524.5

def jd_of(t: datetime) -> real:
    # The Julian date of an instant.
    return (t - datetime("2000-01-01T12:00:00Z")) / duration("P1D") + 2451545.0

def instant(j: real) -> datetime:
    # The instant of a Julian date, to the nanosecond. A Julian date near
    # 2.4 million held as a real steps in about 40 microseconds, so the
    # instant is that real's and not the decimal written for it.
    return datetime("2000-01-01T12:00:00Z") + duration("P1D") * (j - 2451545.0)

def gmst(j: real) -> real:
    # Greenwich mean sidereal time at a Julian date, in degrees from 0 to
    # 360, by Meeus's formula 12.4. A date so far from 2000 that the angle
    # passes the largest real has no sidereal time a real can give, and
    # refuses.
    t = (j - 2451545.0) / 36525
    th = 280.46061837 + 360.98564736629 * (j - 2451545.0) + 0.000387933 * t * t - t * t * t / 38710000
    if not math.isfinite(th) and math.isfinite(j):
        raise ValueError("a Julian date whose sidereal angle passes the largest real")
    return th % 360

def lst(j: real, longitude: real) -> real:
    # Local mean sidereal time at a longitude, east positive. The longitude is
    # taken modulo 360 degrees first, which math.fmod does exactly, so one past
    # a turn keeps the sidereal time beside it.
    return (gmst(j) + math.fmod(longitude, 360.0)) % 360

def hms(degrees: real) -> list<real>:
    # An angle of the sky's rotation as [hours, minutes, seconds], 15° an
    # hour.
    h = degrees / 15
    a = float(math.floor(h))
    b = float(math.floor((h - a) * 60))
    return [a, b, ((h - a) * 60 - b) * 60]

# ── Calculation: the sky's rotation ──────────────────────────────────────────

def altitude(declination: real, latitude: real, hour_angle: real) -> real:
    # A body's altitude above the horizon, from its declination, the
    # observer's latitude and its hour angle. Each angle is taken modulo 360
    # degrees first, which math.fmod does exactly, so one past a turn loses no
    # digits to radians.
    d = math.radians(math.fmod(declination, 360.0))
    p = math.radians(math.fmod(latitude, 360.0))
    return math.degrees(math.asin(math.sin(p) * math.sin(d) + math.cos(p) * math.cos(d) * math.cos(math.radians(math.fmod(hour_angle, 360.0)))))

def semi_arc(declination: real, latitude: real, horizon: real) -> real:
    # The hour angle at which a body crosses the altitude `horizon`, rising
    # and setting symmetrically about its transit: from 0, for one that never
    # reaches it, to 180, for one that never falls below it. Each angle is
    # taken modulo 360 degrees first, as in altitude.
    d = math.radians(math.fmod(declination, 360.0))
    p = math.radians(math.fmod(latitude, 360.0))
    c = (math.sin(math.radians(math.fmod(horizon, 360.0))) - math.sin(p) * math.sin(d)) / (math.cos(p) * math.cos(d))
    if c >= 1:
        return 0.0
    if c <= -1:
        return 180.0
    return math.degrees(math.acos(c))

def rise_set(ra: real, declination: real, latitude: real, longitude: real, day: datetime) -> list<datetime>:
    # A star's rising, transit and setting at a site around its transit on
    # the UT day starting at `day`, over the geometric horizon: [rising,
    # transit, setting], [transit] for a star that never sets, and [] for one
    # that never rises.
    # Each angle is taken modulo 360 degrees first, as in altitude.
    th = gmst(jd_of(day))
    m = ((math.fmod(ra, 360.0) - math.fmod(longitude, 360.0) - math.fmod(th, 360.0)) % 360) / 360.98564736629
    transit = day + duration("P1D") * m
    h = semi_arc(declination, latitude, 0.0)
    if h == 0:
        return []
    if h == 180:
        return [transit]
    w = duration("P1D") * (h / 360.98564736629)
    return [transit - w, transit, transit + w]

def sun(j: real) -> list<real>:
    # The Sun's right ascension and declination at a Julian date, by the
    # Astronomical Almanac's low-precision formulas, good to about 0.01°
    # from 1950 to 2050.
    n = j - 2451545.0
    l = (280.460 + 0.9856474 * n) % 360
    g = math.radians((357.528 + 0.9856003 * n) % 360)
    lam = math.radians(l + 1.915 * math.sin(g) + 0.020 * math.sin(2 * g))
    eps = math.radians(23.439 - 0.0000004 * n)
    ra = math.degrees(math.atan2(math.cos(eps) * math.sin(lam), math.cos(lam))) % 360
    return [ra, math.degrees(math.asin(math.sin(eps) * math.sin(lam)))]

def day_length(latitude: real, j: real) -> real:
    # Hours from sunrise to sunset at a latitude on the day of a Julian date:
    # the Sun's semi-arc to −0.833°, its upper limb at the horizon after
    # refraction, at 15° an hour. 24 where it never sets and 0 where it
    # never rises.
    return 2 * semi_arc(sun(j)[1], latitude, -0.833) / 15

# ── Calculation: orbits ──────────────────────────────────────────────────────

def kepler(mean_anomaly: real, e: real) -> list<real>:
    # The eccentric anomaly E with E − e·sin E = M, by Newton's method from M,
    # or from π where e is above 0.8 and M is small enough to stall it:
    # [E, steps], stopped when a step changes E by less than 1e-15. It is the
    # equation of an ellipse, and an eccentricity outside 0 to 1 refuses. A
    # mean anomaly past half a turn is taken to within half a turn of 0 first,
    # by math.fmod, which is exact, and the turns are added back to E; one
    # below 0 is solved for its size and E given its sign, as the equation is
    # odd. So E is had for any M, where Newton's method from π stalled for an
    # M below -π or past a few turns.
    if not (e >= 0.0 and e <= 1.0):
        raise ValueError("an eccentricity is from 0 up to 1")
    m = mean_anomaly
    base = 0.0
    if abs(m) > math.pi:
        r = math.fmod(m, 2 * math.pi)
        if r > math.pi:
            r = r - 2 * math.pi
        if r < -math.pi:
            r = r + 2 * math.pi
        base = m - r
        m = r
    sign = 1.0
    if m < 0:
        sign = -1.0
        m = -m
    x = m
    if e > 0.8:
        x = math.pi
    n = 0.0
    for k in range(64):
        step = (x - e * math.sin(x) - m) / (1 - e * math.cos(x))
        x = x - step
        n = k + 1.0
        if abs(step) < 1e-15:
            break
    big = sign * x
    if base != 0.0:
        big = base + big
    return [big, n]

def true_anomaly(eccentric_anomaly: real, e: real) -> real:
    # The true anomaly, the angle from perihelion seen from the focus, from
    # the eccentric anomaly.
    h = eccentric_anomaly / 2
    return 2 * math.atan2(math.sqrt(1 + e) * math.sin(h), math.sqrt(1 - e) * math.cos(h))

def two_body(e: real, orbits: real, per_orbit: int) -> list<real>:
    # A body on an orbit of eccentricity e about a mass with GM = 1 and a
    # semi-major axis of 1, from perihelion, by fourth-order Runge–Kutta at
    # `per_orbit` steps an orbit for `orbits` orbits: [the energy at the
    # start, the energy at the end, the distance from where Kepler's
    # equation puts it then]. The steps run in blocks of 64, since a run is
    # charged for each value a `range` holds.
    x = 1 - e
    y = 0.0
    vx = 0.0
    vy = math.sqrt((1 + e) / (1 - e))
    e0 = (vx * vx + vy * vy) / 2 - 1 / math.sqrt(x * x + y * y)
    dt = 2 * math.pi / per_orbit
    total = int(orbits * per_orbit + 0.5)
    done = 0
    for block in range(total // 64 + 1):
        for k in range(64):
            if done == total:
                break
            r1 = (x * x + y * y) ** 1.5
            ax1 = -x / r1
            ay1 = -y / r1
            x2 = x + dt / 2 * vx
            y2 = y + dt / 2 * vy
            u2 = vx + dt / 2 * ax1
            w2 = vy + dt / 2 * ay1
            r2 = (x2 * x2 + y2 * y2) ** 1.5
            ax2 = -x2 / r2
            ay2 = -y2 / r2
            x3 = x + dt / 2 * u2
            y3 = y + dt / 2 * w2
            u3 = vx + dt / 2 * ax2
            w3 = vy + dt / 2 * ay2
            r3 = (x3 * x3 + y3 * y3) ** 1.5
            ax3 = -x3 / r3
            ay3 = -y3 / r3
            x4 = x + dt * u3
            y4 = y + dt * w3
            u4 = vx + dt * ax3
            w4 = vy + dt * ay3
            r4 = (x4 * x4 + y4 * y4) ** 1.5
            ax4 = -x4 / r4
            ay4 = -y4 / r4
            x = x + dt / 6 * (vx + 2 * u2 + 2 * u3 + u4)
            y = y + dt / 6 * (vy + 2 * w2 + 2 * w3 + w4)
            vx = vx + dt / 6 * (ax1 + 2 * ax2 + 2 * ax3 + ax4)
            vy = vy + dt / 6 * (ay1 + 2 * ay2 + 2 * ay3 + ay4)
            done = done + 1
        if done == total:
            break
    e1 = (vx * vx + vy * vy) / 2 - 1 / math.sqrt(x * x + y * y)
    m = (total * dt) % (2 * math.pi)
    big = kepler(m, e)[0]
    kx = math.cos(big) - e
    ky = math.sqrt(1 - e * e) * math.sin(big)
    return [e0, e1, math.sqrt((x - kx) * (x - kx) + (y - ky) * (y - ky))]