Longitude of lunar nodes using Skyfield

Question

I am trying to find out the longitude of ascending/descending moon nodes using Skyfield but unable to find any reference in documentation. Is it possible? Also do any of the JPL Files provide this data already?

Answer 1

Update:

Skyfield’s almanac module now supports this computation directly! See: Lunar Nodes

Original answer, for those wanting these details:

It is easy to at least find them relative to the J2000 ecliptic — which might be fine for dates far from the year 2000 as well, since I think that only the definition of ecliptic longitude changes with the passing years, but not latitude (which is what the nodes care about)?

In any case, you'd precede like this. Let's say you want the ascending node. It must happen within the next 30 days, because that's more than a full orbit of the Moon, so let's look for the day on which the latitude of the Moon passes from negative to positive:

from skyfield.api import load
ts = load.timescale()
eph = load('de421.bsp')
earth = eph['earth']
moon = eph['moon']

t = ts.utc(2018, 1, range(14, 14 + 30))
lat, lon, distance = earth.at(t).observe(moon).ecliptic_latlon()
angle = lat.radians

for i in range(len(angle)):
    if angle[i] < 0 and angle[i+1] > 0:
        break

print(t[i].utc_jpl(), angle[i])
print(t[i+1].utc_jpl(), angle[i+1])

The result is the discovery that the ascending node must happen sometime on January 31st:

A.D. 2018-Jan-31 00:00:00.0000 UT -0.0188679292421
A.D. 2018-Feb-01 00:00:00.0000 UT 0.00522392011676

To find the exact time, install the SciPy library, and ask one of its solvers to find the exact time at which the value reaches zero. You just have to create a little function that takes a number and returns a number, by converting the number to a Skyfield time and then the angle back to a plain number:

from scipy.optimize import brentq

def f(jd):
    t = ts.tt(jd=jd)
    angle, lon, distance = earth.at(t).observe(moon).ecliptic_latlon()
    return angle.radians

node_t = brentq(f, t[i].tt, t[i+1].tt)
print(ts.tt(jd=node_t).utc_jpl())

The result should be the exact moment of the node:

A.D. 2018-Jan-31 18:47:54.5856 UT