Accuracy of an Epicycle Model for Predicting Eclipses

One of the first geometric models for predicting the motion of the sun and moon was devised by Hipparchus over 2000 years ago (in 150 BC). In this model, the moon moves on a small circle called the epicycle, and this epicycle moves on a larger circle called the deferent. Both these motions are uniform. The main reason for introducing this model was to explain the non-uniform motion of the that was observed. The ancient Babylonian astronomers observed that the moon was sometimes faster in its orbit. Indeed, this has effect of making the length of a synodic month variable. In the previous article, I looked at how the Babylonians deduced this, and also how they calculated the length of a month.

When the moon moves forwards on its epicycle, it appears to accelerate, and in the opposite direction it appears to decelerate. The mean speed is represented by the motion of the center of the epicycle on the deferent. The point on the epicycle, where the moon is closer to the Earth, is called perigee. And the point, where the moon is farther away, is called the apogee. In Hipparchus’ model, the moon moves forwards when at perigee, and backwards at apogee. This makes the moon slower at apogee and faster at perigee.

In addition to the epicycle, astronomers were also aware that the path of the moon and the sun were not exactly the same in the sky but were inclined at a small angle. The two points where the moon intersects the sun’s path in the sky are called the nodes. At one node, when the moon crosses, it goes further northwards, this is called the ascending node. The other node where the moon moves southwards, is called the descending node. These nodes are not fixed but were observed to move in the opposite direction to the moon’s motion at a rate which Hipparchus and others judged to be uniform.

This mechanism allowed astronomers to explain eclipses, as eclipses can only occur when the moon crosses the path of the sun at one of the nodes. The addition of the epicycle to account for the moon’s variable speed allowed Hipparchus to predict the time for when an eclipse would occur, as well as to determine how much of the moon or sun is eclipsed. My goal in this article is to determine how accurate this epicycle model is in predicting eclipses, using the best modern values for the parameters.

How to Determine the Parameters for the Model

Hipparchus had found a way to geometrically determine the size of the epicycle for the sun and moon, using just three observations. For the sun, these observations are taken to be the equinoxes and solstices (or equivalently the length of the seasons). For the moon, these observations are taken at lunar eclipses, where the sun is opposite the moon, and therefore the position of the sun allows us to locate the position of the moon.

The details of the method used by Hipparchus can be found in a paper by Duke ¹. I have written the following program in python to implement the procedure:

# Hipparchus' method

from math import *

def angle(x):
    return x%360

# Input the values
T = float(input("Enter the period of orbit (true):"))
T_a = float(input("Enter the period of orbit (anomaly):"))

l1 = float(input("long. of Observation 1:"))
t1 = float(input("time of Observation 1:"))

l2 = float(input("long. of Observation 2:"))
t2 = float(input("time of Observation 2:"))

l3 = float(input("long. of Observation 3:"))
t3 = float(input("time of Observation 3:"))

w = 360 / T  # Mean Angular Speed
w_a = 360/T_a 

w_d = w - w_a #Difference in angular speed

# Convert the angles to radians
l1 = radians(l1)
l2 = radians(l2)
l3 = radians(l3)

w = radians(w)
w_d = radians(w_d)

# Mean motion of the apogee
p = w_d*(t2 - t1)
q = w_d*(t3 - t2)

# Mean Angles, from the center of the deferent
a = w * (t2 - t1) - p
b = w * (t3 - t2) - q

# Longitude, angles from the observer
g = l2 - l1 - p
d = l3 - l2 - q

# Solve the trigonometric problem
y = sin(d) * sin(a - g) - sin(g) * sin(b - d)

x = sin(a + b - g - d) - cos(g) * sin(b - d) - cos(d) * sin(a - g)

z = sin(g) * cos(b - d) + sin(d) * cos(a - g) - sin(g + d)

p = atan2(y, x)
s = atan2(y, z)

# Calculate the radius of the epicycle
r = sin(s) / sin(p)

# Convert the angles to degrees
p = degrees(p)
s = degrees(s)

phi = p + s

# Find the apogee at the second observation
apo = angle(degrees(l2) - p)

print("The radius of the epicycle is:", r)
print("The apogee is at", apo)
print("Daily motion of apogee (Apsidal precession):", degrees(w_d))

For the sun, the input values are the time from vernal equinox to summer solstice (or the length of spring), and the time from summer solstice to autumnal equinox (the length of summer). The values for the longitudes represent the position of the sun at each time. Because of how equinoxes and solstices are defined, these are simply $90\degree$ apart. The time values are given in days. For simplicity, I take the vernal equinox as day 0, and from there use the lengths of the season to determining the other time values. For the length of the seasons, I found the following values from a web search: spring $92.75$ days, summer $93.65$ days. I ran the python program and inputted these values. Note: for the third observation simply add $92.75+93.65=186.4$.

 
Enter the period of orbit (true):365.25
Enter the period of orbit (anomaly):365.25
long. of Observation 1:0
time of Observation 1:0
long. of Observation 2:90
time of Observation 2:92.75
long. of Observation 3:180
time of Observation 3:186.4
The radius of the epicycle is: 0.03337403397527309
The apogee is at 103.41173814032092
Daily motion of apogee (Apsidal precession): 0.0
 

For an epicycle, the apogee is the point on the epicycle where the sun reaches its furthest point from Earth, and the size of the epicycle is measured relative to a deferent of unit radius. Notice the apogee corresponds to the aphelion in our modern theory. Also notice that the radius of the epicycle is very nearly twice the value of the eccentricity of Earth’s orbit $0.03337/2 \approx 0.01669$. So, Hipparchus’ procedure from 2000 years ago, can give us a reasonably accurate estimate for the elliptical orbit of the Earth.

Sticking with the epicycle model, our next step is to determine the epoch of the sun. That is, we choose some arbitrary date and time, which will represent our starting position, then determine the exact position of the sun at that time. For my epoch date, I will choose January 1, 2001, at 0:00 UTC. To find the position of the sun at this time, I simply need to find the time between the nearest vernal equinox and the start date of our epoch. For this, I looked up the date of the vernal equinox online and found March 20, 2001, 13:19 UTC to be the closest equinox date. The difference between this equinox and our start date is $78.555$ days. To calculate the mean longitude of the sun (the position of the center of its epicycle), we simply multiply by 360 and divide by the length of the year: $78.555/365.25\times 360 = 77.425$. The result is then subtracted from the position of the vernal equinox $360-77.426 =282.57\degree$. This is the position of the sun at the starting point of our epoch. The angle between this and the apogee is called the mean anomaly. This is $282.57-103.412 = 179.159\degree$. I wrote a simple epoch calculator to find the position

# Epoch Calculator

from math import *

def mod(a,b): #Works for negative numbers
	return a%b

T = float(input("Period of orbit (true):"))
T_a = float(input("Period of orbit (anomaly):"))
radius = float(input("Radius of epicycle:"))
days = float(input("Time of observation since epoch (in days) :"))
apogee = float(input("Longitude of apogee:"))
longt = float(input("Longitude of observation:"))

true_ano = mod(longt - apogee,360) #True Anomaly - Angle between observation and apogee

#Determine the equation of anomaly (i.e difference between true and mean anomaly) 
x = 1 + radius*cos(radians(true_ano))
y = radius*sin(radians(true_ano))
eq = degrees(atan(y/x))

mean_ano = mod(true_ano + eq,360) #Mean anomaly - difference between mean position or apogee

longm = mod(mean_ano + apogee,360) #Mean longitude of observation

#Mean longitude of epoch
longm_ep = mod(longm - (days/T)*360,360)

#Mean anomaly of epoch
mean_ano_ep = mod(longm_ep - apogee,360)

# Equation of anomaly (at epoch)
x = 1 + radius*cos(radians(mean_ano_ep))
y = radius*sin(radians(mean_ano_ep))
eq = degrees(atan(y/x))

#True anomaly at epoch
ano_ep = mod(mean_ano_ep - eq,360)

#Longitude of epoch
long_epoch = mod(ano_ep + apogee,360)

#Motion of the apogee
w_d = 360/T - 360/T_a

apogee_ep = apogee - w_d*days 

print(f"Mean longitude at epoch:{longm_ep}")
print(f"Longitude at epoch:{long_epoch}")
print(f"Longitude of apogee at epoch:{apogee_ep}")

The following is the result of running the program for the sun:

 
Period of orbit (true):365.25
Period of orbit (anomaly):365.25
Radius of epicycle:0.0334
Time of observation since epoch (in days) :78.555
Longitude of apogee:103.41
Longitude of observation:0
Mean longitude at epoch:280.69876207754385
Longitude at epoch:280.60511633619694
Longitude of apogee at epoch:103.41
 

To continue on to the next step, I have written a generic epicycle calculator in python:

# Epicycle Calculator

from math import *

def angle(a):
	return a%360

T = float(input("Enter period of orbit (true):"))
T_a = float(input("Enter period of orbit (anomaly):"))
radius = float(input("radius:"))
epoch = float(input("Mean Longitude at epoch:"))
apogee_ep = float(input("Longitude of apogee at epoch:"))

while True: #Repeat
	days = float(input("\nTime of observation (days since epoch):"))

    #Motion of the apogee
    w_d = 360/T - 360/T_a
    apogee = apogee_ep + w_d*days
    
	mean_long = epoch + days/T*360
	mean_long = angle(mean_long)

	print(f"Mean Longitude:{mean_long}")

	epi_ano = mean_long - apogee
	print(f"Epicyclic anomaly:{epi_ano}")


	x = 1 + radius*cos(radians(epi_ano))
	y = radius*sin(radians(epi_ano))

	#Equation of anomaly (i.e difference between them)
	eq = degrees(atan(y/x))
	
	#True anomaly
	ano = angle(epi_ano - eq)
	print(f"True Anomaly:{ano}")

	#True longitude
	longt = angle(ano + apogee)
	print(f"Longitude:{longt}")

With the parameters for the sun now established, the next step is to determine the parameters for the moon. We proceed by finding three lunar eclipses and determining the time of mid-eclipse for each. At this time, the moon and the sun are $180\degree$ apart in longitude. So, we find the exact position of the moon by finding the position of the sun and reflecting this by $180\degree$. For data on lunar eclipses, I found three observations close to our epoch from the NASA Eclipse Catalogue.

The following is the list of eclipses:

• The first eclipse was on January 9, 2001, 20:22 UTC. The number of days since the start of our epoch is $8.85$ days.
• The second eclipse was on July 5, 2001, 14:56 UTC. The number of days since the start of our epoch is $185.62$ days
• The third eclipse was on December 30, 2001, 10:30 UTC. The number of days since the start of our epoch is $363.44$ days

Now we can calculate the position of the sun using our epicycle calculator. The result is as follows:

 
Enter period of orbit (true) :365.25
Enter period of orbit (anomaly) :365.25
radius:0.0334
Mean Longitude at epoch:280.7
Longitude of apogee:103.41

Time of observation (days since epoch):8.85
Mean Longitude:289.4227926078029
Epicyclic anomaly:186.01279260780288
True Anomaly:186.2201378195743
Longitude:289.6301378195743

Time of observation (days since epoch):185.62
Mean Longitude:103.65195071868584
Epicyclic anomaly:0.2419507186858425
True Anomaly:0.23413077261939347
Longitude:103.64413077261939

Time of observation (days since epoch):363.44
Mean Longitude:278.91601642710475
Epicyclic anomaly:175.50601642710475
True Anomaly:175.35090686500425
Longitude:278.76090686500424

Time of observation (days since epoch): 
 

From the values, we see that the longitude of the sun at each eclipse are given as $289.63\degree$, $103.64\degree$ and $278.76\degree$. To find the position of the moon, we simply need to flip these values $180\degree$; we can either add or subtract $180\degree$.

So, the longitude of the moon at each eclipse is: $109.63\degree$, $283.64$ and $98.76\degree$. We can now input these values into our Hipparchus’ method calculator from earlier to obtain the parameters for the moon. The results are as follows:

 
Enter the period of orbit (true):27.3216
Enter the period of orbit (anomaly):27.5546
long. of Observation 1:109.63
time of Observation 1:8.85
long. of Observation 2:283.64
time of Observation 2:185.62
long. of Observation 3:98.76
time of Observation 3:363.44
The radius of the epicycle is: 0.08573619368442606
The apogee is at 326.67908028855817
Daily motion of apogee (Apsidal precession): 0.11141872441248958
 

Now our final step is to determine the epoch of the moon, using our epoch calculator. For this I will use the second observation, since the longitude of apogee calculated from our program is based on this value. The result for our epoch is:

 
Period of orbit (true):27.3216
Period of orbit (anomaly):27.5546
Radius of epicycle:0.086
Time of observation since epoch (in days) :185.62
Longitude of apogee:326.68
Longitude of observation:283.64
Mean longitude at epoch:354.67795053849795
Longitude at epoch:352.5290740672692
Longitude of apogee at epoch:305.9984563745537
 

Predicting Lunar Eclipses

There are two conditions required for an eclipse to occur:

• The angle between the sun and the moon — referred to as the elongation — must be either $0\degree$ or $180\degree$. If the elongation is $0\degree$, then the eclipse is a solar eclipse. And if the elongation is $180\degree$, the eclipse is a lunar eclipse.
• The sun must be sufficiently close to one of the nodes. These are the two points where the moon’s orbit intersects the path of the sun. The angle between the node, required for an eclipse to occur, is referred to as an ecliptic limit.
  • For a lunar eclipse, the limit is $18\degree$ for penumbral, $12\degree$ for partial and $5\degree$ for total.
  • For a solar eclipse, the limits are $18\degree$ for partial, and $11\degree$ for total or annular.

These values can vary based on the distance to the moon.

Another aspect of eclipse prediction is the magnitude of the eclipse, which refers to the fraction of the moon or sun that is covered. The closer the sun is to the node, the greater the magnitude of the eclipse. The magnitude is also affected by the relative sizes of the eclipsed body and the eclipsing body. The sun and moon have approximately the same angular diameter ($0.5 \degree$), and the size of the earth’s shadow on the moon is approximately $2.5$ times the diameter of the moon ( $1.25\degree$).

We have the following formula for determining the magnitude of an eclipse. $R_S$ is the smaller of the radius (i.e. half the angular diameter), $R_L$ is the larger radius, and $d$ is the distance between the centers. $$M = \frac{R_L + R_S - d}{2R_S}$$

Before predicting eclipses, there is one final parameter we must determine: the position of the lunar nodes at epoch. For this we will need an eclipse of known magnitude, from which we will calculate the angle from the node. We consult the NASA Eclipse Catalogue. For this we will select the second eclipse, which occurred on July 5, 2001, because this is a partial eclipse. The magnitude of the eclipse is given as $0.4947$. From this we can calculate the distance between the center of the moon, and the center of the Earth’s shadow. $$0.4947 = =\frac{0.625+0.25 - d}{2\times0.25}$$

From which we get $d=0.627$

Now this distance represents the latitude of the moon. $\sin(0.627) = \sin(5)\times\sin(\varphi_L)$

We get the angle of the observation from the nodes as $\varphi_L = 7.2\degree$.

From observation, the shadow covers the northern half of the moon’s disk. This mean’s the moon’s latitude is southerly. And since the eclipse occurs at a descending node, this means that it is heading south. We can therefore say that the moon is ahead of the node. And thus, we subtract $7.2$ from the longitude of the moon to determine the position of the descending node. Now earlier, we had determined the longitude of the moon at this eclipse to be $283.64\degree$. And so the longitude of the descending node is $276.44\degree$. And if we prefer, we can take the ascending node, which is $276.44-180=96.44\degree$. We can reduce these values to epoch by using the value for the draconic month, which is $27.212$ days, and the time of our observation from epoch:

$$ 276.44 - \left(\frac{360}{27.321} - \frac{360}{27.212} \right)\times 185.62 = 284.65\degree$$.

So, at the time of our epoch, the longitude of the descending node is $284.65\degree$ and the longitude of the ascending node is $96.64\degree$.

Now that we have determined all the parameters for the sun and the moon, the next step is to predict the eclipses. The following is a simple algorithm for predicting lunar eclipses:

• Start with the known parameters for the sun and moon at some initial time, such as our epoch.
• Determine the difference between the longitude of the sun and moon; this is the elongation. If the elongation is $180\degree$, then we already have a full moon. Otherwise, we add the difference to the mean position of the moon. This will bring us close to a full moon, but not exactly, because the moon must also have moved on its epicycle. But in that case, we just add the difference again until the elongation is close enough to $180\degree$, then we have a full moon.
• After we have determined the time of full moon, we then find the difference between moon’s longitude and the nearest node; this is the argument of latitude. We can quickly tell if the eclipse is total or partial by the value of the argument of latitude. If it is less than $12\degree$ we have a partial eclipse. If it is less than $5\degree$, we have a total eclipse, and if it is greater than $18\degree$, no eclipse. We can also compute the magnitude of the eclipse.
• We can then progress the moon’s position forward $360\degree$ and repeat the procedure to predict the next eclipse. And we can even go backwards to determine past eclipses, before our epoch date.

I wrote the following code in python to predict lunar eclipses:

# Lunar Eclipse Prediction

from math import *
from datetime import datetime, date, timedelta

def angle(x):
        return x%360
def angle_diff(x,y):
    diff = abs(x - y) %360
    return min(diff, 360 - diff)
        
#Parameters for the sun
sun_r = 0.033 #Radius of epicycle
sun_T = 365.25 # Orbital period
sun_Ta = 365.25 # Period of anomaly
sun_w = 360/sun_T
sun_wa = 360/sun_Ta

# Epoch positions for the sun
sun_lm_ep = 280.6  #Mean Longitude
sun_apo_ep = 103.41 #  Apogee

#Parameters for the moon
moon_r = 0.086 #Radius of epicycle
moon_T = 27.321582 #Orbital period
moon_Ta = 27.554550 #Period of anomaly
moon_Td = 27.212221 # Period of latitude (draconic month)
moon_w = 360/moon_T
moon_wa = 360/moon_Ta
moon_wd = 360/moon_Td

#Epoch positions for the moon
moon_lm_ep = 354.67 #Mean longitude
moon_apo_ep = 306 #Apogee
dnode_ep = 284.65 #Longitude of descending node
moon_i = 5 #Orbital inclination of the moon

# Corrections
dnode_ep += 1.5532
moon_lm_ep -= 0.227
sun_lm_ep += 0.05

#moon_apo_ep  +=  8.5

#Copied from our epicycle calculator and modified
def epi_calc(T,T_a,radius,epoch,apogee_ep,days):
    w_d = 360/T - 360/T_a
    apogee = apogee_ep + w_d*days

    mean_long = epoch + days/T*360
    mean_long = angle(mean_long)
    
    epi_ano = mean_long - apogee

    x = 1 + radius*cos(radians(epi_ano))
    y = radius*sin(radians(epi_ano))

    eq = degrees(atan(y/x))
    
    ano = angle(epi_ano - eq)

    longt = angle(ano + apogee)
    return longt

# Calculate the date in a printable format
def calc_datetime(days):
    epoch = datetime.strptime("2001-01-01 00:00:00","%Y-%m-%d %H:%M:%S")
    epoch += timedelta(days)
    return epoch

#Find the exact time of the full moon
def find_moon(n):
    # n is the lunation number, i.e the number of full moons
    # since epoch.
    # Negative values of n give eclipses before epoch.

    # Point opposite the moon
    opp_lm_ep = angle(moon_lm_ep - 180)
    
    #Determine the angle reuired to reach opposite the moon
    elong_m = angle(sun_lm_ep - opp_lm_ep)

    # Difference in speed between sun and moon
    w_e = moon_w - sun_w

    T_syn = 360/w_e # Synodic month (approx 29.5)


    #Calculate the number of days to the nth mean full moon
    days_m = T_syn * n +  elong_m/w_e

    # Repeat until the angle between sun and moon is sufficiently close to 180
    for i in range(5):
        #Calculate the mean longitude for the sun and moon
        moon_lm = angle(moon_lm_ep + moon_w * days_m)
        sun_lm = angle(sun_lm_ep + sun_w *  days_m)

        # Calculate the true longitude for the sun and moon
        moon_lt = epi_calc(moon_T, moon_Ta, moon_r,  moon_lm_ep, moon_apo_ep, days_m)
        sun_lt = epi_calc(sun_T, sun_Ta, sun_r, sun_lm_ep, sun_apo_ep, days_m)

        #Point opposite the moon
        opp_lt = angle(moon_lt - 180)
        
        # Angle between sun and moon
        elong_t = angle(sun_lt - opp_lt)

        elong_t = -(360 - elong_t) if elong_t > 180 else elong_t

        # Approximate number of days left for the moon to reach opposite the sun
        days_x = elong_t/w_e

        days_m += days_x
    
    return days_m

# Tell whether a given full moon is an ecipse and predict its magnitude
def predict_eclipse(n):
    
    #Time of full moon from epoch
    days = find_moon(n)

    #Motion of the node
    w_d = (moon_w - moon_wd)

    # Calculate the longitude of descending node
    dnode = angle(dnode_ep + w_d * days)

    #Ascending node
    anode = angle(dnode - 180)

    #Calculate the longitude for the sun
    sun_lt = epi_calc(sun_T, sun_Ta, sun_r, sun_lm_ep, sun_apo_ep, days)

    #Calculate the argument of latitude
    arg_d = angle_diff(sun_lt , dnode)
    arg_a = angle_diff(sun_lt , anode)
    arg_lat = min(arg_d,arg_a)

    # Moon Latitude - also distance from center of shadow
    moon_lat = degrees(asin(sin(radians(arg_lat))*sin(radians(moon_i))))
    
    # Eclipse Magnitude
    Mag = 2*(0.875 - moon_lat)

    # Determine whether there is an eclipse, and type
    if(arg_lat > 18):
        return False # No eclipse

    ecl_type = "N" #Default penumbral

    if(Mag > 0):
       ecl_type = "P" #Partial

    if(Mag > 1):
        ecl_type = "T" #Total
    
    Lun = n + 12 #Lunation number, to match the standard epoch (J2000) used by NASA

    show_date = calc_datetime(days)

    print(f"{show_date}\t{Lun}\t{Mag}\t{ecl_type}")
    return True

print("Date and Time\t\t   Lunation\t Magnitude\t \tType")

for i in range(223):
    predict_eclipse(i)

The following is the result of running our program:

 
Date and Time              Lunation      Magnitude              Type
2001-01-09 20:23:33.640631      12      1.0809695425549712      T
2001-07-05 14:50:43.257775      18      0.4947041996342709      P
2001-12-30 10:34:10.150156      24      -0.2845878029645501     N
2002-05-26 11:17:00.870877      29      -0.7024032429161089     N
2002-06-24 21:14:54.720935      30      -0.9190359275820748     N
2002-11-20 00:48:22.741719      35      -0.3836926271353569     N
2003-05-16 03:03:29.463578      41      0.7702240704747019      P
2003-11-09 00:56:25.261840      47      0.8963103449460311      P
2004-05-04 20:12:29.122011      53      1.2328396685035106      T
2004-10-28 03:09:46.188322      59      1.2826309707381098      T
2005-04-24 09:06:58.265271      65      -0.21292774979164752    N
2005-10-17 12:00:43.609386      71      -0.09124302054516287    N
2006-03-14 22:16:01.319825      76      -0.25504311654414025    N
2006-09-07 19:11:13.121962      82      -0.21257941094613186    N
2007-03-03 22:32:51.940187      88      1.041565703102378       T
2007-08-28 10:51:43.503824      94      1.2604580651885184      T
2008-02-21 02:41:33.456188      100     1.1311343495494386      T
2008-08-16 20:33:08.502797      106     0.7916847882242458      P
2009-02-09 13:47:32.099062      112     -0.22875730503424507    N
2009-07-07 08:43:09.960390      117     -1.0520738917936905     N
2009-08-05 23:53:08.070463      118     -0.5975521083329518     N
2009-12-31 19:22:40.714607      123     -0.33516781163293885    N
2010-06-26 11:14:41.934383      129     0.3296586519817353      P
2010-12-21 07:57:01.640007      135     1.017702413512806       T
2011-06-15 19:34:34.753541      141     1.7194986286965974      T
2011-12-10 13:30:51.507130      147     1.1671098273897351      T
2012-06-04 10:25:49.077593      153     0.2253696233695679      P
2012-11-28 13:55:15.569668      159     -0.104286407661633      N
2013-04-25 19:32:14.266452      164     -0.40544005300912156    N
2013-05-25 03:54:15.742093      165     -1.2483695255682008     N
2013-10-18 23:39:14.308246      170     -0.512836398747341      N
2014-04-15 06:26:54.301449      176     1.0129092116063738      T
2014-10-08 10:39:50.772621      182     0.8790673519806677      P
2015-04-04 10:22:36.393058      188     1.1170692993842082      T
2015-09-28 02:59:50.460096      194     1.1605820624172454      T
2016-03-23 10:50:14.583030      200     -0.2036310027225292     N
2016-09-16 19:18:12.422130      206     -0.2974179149950271     N
2017-02-10 23:41:18.104347      211     -0.3843456973892252     N
2017-08-07 17:12:29.246199      217     0.02691962948749027     P
2018-01-31 12:35:10.425051      223     0.979916870079237       P
2018-07-27 19:24:40.166134      229     1.4246767158897369      T
 

The data is also presented in the following table:

We can compare this with the value from the NASA Catalogue

Predicting Solar Eclipses

We might think that we can modify the program above to predict solar eclipses, by modifying the function find_moon() to determine the time of new moon instead and adjusting the calculations for magnitude and eclipse type. This idea would be correct, however we run into the issue of parallax. Since the moon is close to the Earth, different observers do not see the moon in exactly the same position in the sky at the same time. For a solar eclipse to occur, the center of the moon must line up with the sun, as viewed from the observer.

When we calculate the position of the moon from theory, the position we determine corresponds to what would be seen by a hypothetical observer, viewing from the center of Earth. This is the geocentric position of the moon — or we may say the true position. But the position of the moon that we actually see depends on our location on the Earth — this is the apparent position. An observer who sees the moon directly overhead — at zenith, sees the moon at its true position, because a straight line can be drawn from the center of the Earth, through the observer and the moon. The lower the moon appears in the sky, the further it is from the true position. The angle between the true and apparent position is what we refer to as parallax. The standard value for parallax is usually given for apparent position at the horizon (i.e $90\degree$ from zenith), where the value is maximum.

There are a few ways to measure parallax from observations. One way, suggested by Hipparchus, is to observe a solar eclipse from two locations on the Earth. The difference in the observed magnitude indicates the different positions observed. This method however requires us to assume something about the parallax of the sun. The more modern method involves measuring the angle between the moon and zenith, when it reaches its highest point — at meridian, from two different locations on the Earth on the same meridian.

When we measure the parallax, the value comes out to be about $0.95\degree$. We can, calculate the distance to the moon from this:

$$d = \frac {R_E}{\sin(0.95\degree)} = 60.3\times 6371 \mathrm{km} $$ $$d = 384, 171 \mathrm{km}$$

The effect of parallax means that the suggested modification of our program would only be able to determine the time of a solar eclipse for an observer who sees the sun and moon directly overhead — these observers are said to be at the subsolar and sublunar points. Our program would therefore only predict eclipses which occur at noon for an observer in the tropics. This is however not a bad starting point.

Let us modify our program in the following manner to predict the time of new moon.

def find_moon(n):
    # n is the lunation number, i.e the number of new moons
    # since epoch.
    # Negative values of n give eclipses before epoch.

    #Determine the angle required for  the sun and moon to align
    elong_m = angle(sun_lm_ep - moon_lm_ep)

    # Difference in speed between sun and moon
    w_e = moon_w - sun_w

    T_syn = 360/w_e # Synodic month (approx 29.5)

    #Calculate the number of days to the nth mean full moon
    days_m = T_syn * n +  elong_m/w_e

    # Repeat until the angle between sun and moon is sufficiently close to 180
    for i in range(5):
        #Calculate the mean longitude for the sun and moon
        moon_lm = angle(moon_lm_ep + moon_w * days_m)
        sun_lm = angle(sun_lm_ep + sun_w *  days_m)

        # Calculate the true longitude for the sun and moon
        moon_lt = epi_calc(moon_T, moon_Ta, moon_r,  moon_lm_ep, moon_apo_ep, days_m)
        sun_lt = epi_calc(sun_T, sun_Ta, sun_r, sun_lm_ep, sun_apo_ep, days_m)


        # Angle between sun and moon
        elong_t = angle(sun_lt - moon_lt)

        elong_t = -(360 - elong_t) if elong_t > 180 else elong_t

        # Approximate number of days left for the moon to reach opposite the sun
        days_x = elong_t/w_e

        days_m += days_x
    
    return days_m

Get new posts via email:

Subscribe

Powered by

Reading Time: 23 minutes

Word Count: 4814