Pimennyskaudet

Varsinkin paperilla laskemisen aikakaudella muitakin nopeampia mutta joskus epätarkempia ehtoja on käytetty pimennyksien kriteereinä.

Kuvan 1.5 tilanteessa voitaisiin määrittää kriteerit uudenkuun hetkellä sekä Kuun korkeudelle pystysuunnassa (ekliptinen leveysaste) ${\beta }_m$ että Auringon etäisyydelle solmusta sivusuunnassa (jana SN).

Koska radan inklinaatiossa $I$ sekä kappaleiden puolikkaissa kulmaläpimitoissa $s_s$ ja $s_m$ on ajan myötä pientä vaihtelua, tällaiset rajat eivät kuitenkaan ole aivan tarkkoja, toisin kuin yllä olevat $\sigma$-rajat.

Näille suureille voi kuitenkin antaa rajoja, jotka kertovat, milloin pimennys on varma, ja milloin se on mahdollinen. Vanhempi painos Explanatory Supplementistä (1961) antaa rajat pystysuunnassa: Auringonpimennys on mahdollinen, kun ${\beta }_m<1.58^{\circ }$ ja varma, kun ${\beta }_m<1.41^{\circ }$.

Espenak antaa rajat Auringon sivusuuntaiselle etäisyydelle auringonpimennyksille ($18.59^{\circ }$ ja $15.39^{\circ }$) ja kuunpimennyksille ($16.8^{\circ }$ ja $15.2^{\circ }$).

Nämä sivusuuntaiset etäisyydet ovat käteviä, jos halutaan merkitä alla esiteltävien pimennyskausien paikkoja kalenteriin.

from numpy import arctan, cos, degrees
from skyfield import almanac
from skyfield.api import load
from skyfield.framelib import ecliptic_frame

R_EARTH = 6378.1366           # km, IAU equatorial radius
R_SUN = 696000.0
R_MOON = 0.2725076 * R_EARTH  # Moon/Earth radius ratio, IAU

ts = load.timescale()
eph = load('de440s.bsp')
earth, sun, moon = eph['earth'], eph['sun'], eph['moon']

t0, t1 = ts.utc(2039, 1, 1), ts.utc(2040, 1, 1)
t, phase = almanac.find_discrete(t0, t1, almanac.moon_phases(eph))

def lonlat(body, ti):
    lat, lon, _ = earth.at(ti).observe(body).apparent().frame_latlon(ecliptic_frame)
    return lon.degrees, lat.degrees

def wrap(d):
    return (d + 180) % 360 - 180

for ti, ph in zip(t, phase):
    if ph not in (0, 2):                        # new and full moons only
        continue

    a, b = ts.tt_jd(ti.tt - 1 / 1440), ts.tt_jd(ti.tt + 1 / 1440)
    ml_a, mb_a = lonlat(moon, a)
    ml_b, mb_b = lonlat(moon, b)
    sl_a, _ = lonlat(sun, a)
    sl_b, _ = lonlat(sun, b)

    beta = lonlat(moon, ti)[1]
    lam = wrap(ml_b - ml_a) / wrap(sl_b - sl_a)
    tan_I = (mb_b - mb_a) / wrap(ml_b - ml_a)
    I_ = arctan(lam / (lam - 1) * tan_I)
    sigma = abs(beta) * cos(I_)

    e = earth.at(ti)
    d_sun = e.observe(sun).apparent().distance().km
    d_moon = e.observe(moon).apparent().distance().km

    s_s = degrees(R_SUN / d_sun)
    s_m = degrees(R_MOON / d_moon)
    pi_m = degrees(R_EARTH / d_moon)
    pi_s = degrees(R_EARTH / d_sun)

    if ph == 0:                                 # new moon: chapter 2's limits
        partial = s_s + s_m + pi_m - pi_s
        central = abs(s_s - s_m + pi_m - pi_s)
        if sigma > partial:
            verdict = '-'
        elif sigma > central:
            verdict = 'partial solar'
        else:
            verdict = 'total solar' if s_m > s_s else 'annular solar'
    else:                                       # full moon: chapter 3's limits
        pi_1 = 0.998340 * pi_m
        f_1 = 1.02 * (pi_1 + pi_s + s_s)
        f_2 = 1.02 * (pi_1 + pi_s - s_s)
        if sigma < f_2 - s_m:
            verdict = 'total lunar'
        elif sigma < f_2 + s_m:
            verdict = 'partial lunar'
        elif sigma < f_1 + s_m:
            verdict = 'penumbral lunar'
        else:
            verdict = '-'

    print(f"{ti.utc_strftime('%Y-%m-%d')}  {verdict}")

from numpy import cross
from skyfield import almanac
from skyfield.api import load
from skyfield.framelib import ecliptic_frame

ts = load.timescale()
eph = load('de440s.bsp')
earth, sun, moon = eph['earth'], eph['sun'], eph['moon']

def sun_south_of_orbit(t):
    """Is the Sun on the south side of the plane of the Moon's orbit?"""
    r, v = (moon - earth).at(t).frame_xyz_and_velocity(ecliptic_frame)
    h = cross(r.km, v.km_per_s, axis=0)   # normal to the plane of the orbit
    s = earth.at(t).observe(sun).apparent().frame_xyz(ecliptic_frame).km
    return (s * h).sum(axis=0) < 0

sun_south_of_orbit.step_days = 5     # the Sun needs half a year between nodes

t, south = almanac.find_discrete(ts.utc(2039, 1, 1), ts.utc(2040, 1, 1),
                                 sun_south_of_orbit)

for ti, s in zip(t, south):
    # the Sun passes to the south side at the ascending node
    print(f"{ti.utc_strftime('%Y-%m-%d %H:%M')}  "
          f"{'ascending' if s else 'descending'}")