PyEphem: how to calculate position of the Moon nodes?

Question

Is there a way to get the longitude of the Moon's ascending/descending nodes (North/South nodes) in pyEphem?

Thanks

Answer 1

It does not directly provide a method for calculating the positions of the Moon's nodes. But you could calculate the positions of the ascending and descending nodes by using the positions of the Moon and the Earth.

Calculate the ecliptic longitude of the Moon and the Sun. Find the points where the Moon's orbit intersects the ecliptic plane (nodes).

import ephem
import math

def moon_ascending_node(date):
    observer = ephem.Observer()
    observer.date = date

    # Sun's position
    sun = ephem.Sun(observer)
    sun_ecl = ephem.Ecliptic(sun)

    # Moon's position
    moon = ephem.Moon(observer)
    moon_ecl = ephem.Ecliptic(moon)

    # Longitude of the ascending node
    delta_lambda = moon_ecl.lon - sun_ecl.lon
    while delta_lambda < 0:
        delta_lambda += 2 * math.pi
    while delta_lambda >= 2 * math.pi:
        delta_lambda -= 2 * math.pi

    if delta_lambda > math.pi:
        delta_lambda -= 2 * math.pi

    ascending_node_lon = moon_ecl.lon - delta_lambda
    while ascending_node_lon < 0:
        ascending_node_lon += 2 * math.pi
    while ascending_node_lon >= 2 * math.pi:
        ascending_node_lon -= 2 * math.pi

    return ascending_node_lon

date = '2024/07/01'
ascending_node_lon = moon_ascending_node(date)
print(f"Longitude of the Moon's ascending node on {date}: {math.degrees(ascending_node_lon)} degr

For both Rahu (North Lunar Node) and Ketu (South Lunar Node). The positions of these nodes can be calculated using the ecliptic longitudes of the Moon and the Sun.

import ephem
import math

def moon_nodes(date):
    observer = ephem.Observer()
    observer.date = date

    # Sun's position
    sun = ephem.Sun(observer)
    sun_ecl = ephem.Ecliptic(sun)

    # Moon's position
    moon = ephem.Moon(observer)
    moon_ecl = ephem.Ecliptic(moon)

    # Longitude of the ascending node (Rahu)
    delta_lambda = moon_ecl.lon - sun_ecl.lon
    while delta_lambda < 0:
        delta_lambda += 2 * math.pi
    while delta_lambda >= 2 * math.pi:
        delta_lambda -= 2 * math.pi

    if delta_lambda > math.pi:
        delta_lambda -= 2 * math.pi

    ascending_node_lon = moon_ecl.lon - delta_lambda
    while ascending_node_lon < 0:
        ascending_node_lon += 2 * math.pi
    while ascending_node_lon >= 2 * math.pi:
        ascending_node_lon -= 2 * math.pi

    # Longitude of the descending node (Ketu) is 180 degrees away from Rahu
    descending_node_lon = ascending_node_lon + math.pi
    while descending_node_lon >= 2 * math.pi:
        descending_node_lon -= 2 * math.pi

    return math.degrees(ascending_node_lon), math.degrees(descending_node_lon)

date = '2024/07/01'
rahu_lon, ketu_lon = moon_nodes(date)
print(f"Longitude of Rahu (North Lunar Node) on {date}: {rahu_lon} degrees")
print(f"Longitude of Ketu (South Lunar Node) on {date}: {ketu_lon} degrees")

Answer 2

Pyephem can compute positions of planets on certain dates, but not vice versa, so I usually store planetary positions calculated over a period of time in a data frame and work from there. This is how I search for Moon's nodes:

import ephem
import pandas as pd
import numpy as np

# Let's create a data frame with 2 columns: lon and lat to store
# the daily ecliptic coordinates of Earth's Moon for the year 2001.
# I like to work in degrees of the circle, so I always convert from radians.
Mo = ephem.Moon()
df = pd.DataFrame(index=pd.date_range(start='2001-01-01', end='2001-12-31', freq='D'))
df['lon'] = df.index.map(lambda d: [Mo.compute(d, epoch=d), np.rad2deg(ephem.Ecliptic(Mo).lon)][1])
df['lat'] = df.index.map(lambda d: [Mo.compute(d, epoch=d), np.rad2deg(ephem.Ecliptic(Mo).lat)][1])

# To find the Moon's nodes you need to first find the dates when Moon
# is crossing the ecliptic (so its latitude changes sign), then find
# the longitudes for those dates.
xdates = df.index[np.where(np.diff(np.sign(df['lat'].values)) != 0)]
nodes = df['lon'][xdates]

To verify the results first print the df stats and you should see something like this:

In [180]: df.describe()
Out[180]:
              lon         lat
count  365.000000  365.000000
mean   181.984888   -0.285830
std    107.690034    3.642805
min      0.910882   -5.298351
25%     85.118205   -3.918240
50%    184.577992   -0.595535
75%    277.629225    3.290252
max    359.851182    5.285288

Note the max latitude is 5.285288 which equals to the angle between Moon's and Sun's orbital planes. Next print xdates:

In [181]: xdates
Out[181]:
DatetimeIndex(['2001-01-09', '2001-01-22', '2001-02-06', '2001-02-19',
               '2001-03-05', '2001-03-18', '2001-04-01', '2001-04-14',
               '2001-04-28', '2001-05-11', '2001-05-25', '2001-06-07',
               '2001-06-21', '2001-07-05', '2001-07-19', '2001-08-01',
               '2001-08-15', '2001-08-28', '2001-09-11', '2001-09-24',
               '2001-10-08', '2001-10-21', '2001-11-04', '2001-11-17',
               '2001-12-02', '2001-12-15', '2001-12-29'],
              dtype='datetime64[ns]', freq=None)

and compare the dates with a third party source, ie. Moon's crossing dates for the year 2001 can be found here: http://astropixels.com/ephemeris/moon/moonnodes2001.html . Note that the sequence of nodes is asc, desc, asc, desc, etc. and any two consecutive nodes are about 180 degrees apart.

This method works best if you need to quickly find the dates of the Moon's crossings with the Ecliptic, because the longitude values are computed for midnight and the Moon can travel up to ca. 15 degrees of longitude per day, so for a better accuracy you'll have to go through each minute of the crossing dates to find the actual time and longitude with 1-degree accuracy OR use this dirty little hack:

# Calculate lon and lat diffs on the main data frame:
df['lon_diff'] = df['lon'].diff()
df['lat_diff'] = df['lat'].diff()

# Export the nodes' rows into a separate data frame and for each row there
# find the longitude of the crossing by the formula:
#  lon_x = lon + (lat/lat_diff) * lon_diff
nodes = df.ix[xdates]
nodes['lon_x'] = nodes.apply(lambda r: r['lon'] + np.abs(r['lat']/r['lat_diff'])*r['lon_diff'], axis=1)

The diffs tell us how many degrees of lon or lat does the Moon travel per day, so the ratio lat:lat_diff tells how far was the Moon from the Ecliptic at midnight or in other words what fraction of 1 day did the Moon need to actually reach 0 degree lat. This value multiplied by lon_diff on the day gives the number of degrees of longitude between the Moon's position at midnight and the position at the crossing. I hope this method is accurate enough for you :)