Is there a way to get the longitude of the Moon's ascending/descending nodes (North/South nodes) in pyEphem?
Thanks
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 :)