4. The orbital elements
The primary orbital elements are here denoted as:
N = longitude of the ascending node
i = inclination to the ecliptic (plane of the Earth's orbit)
w = argument of perihelion
a = semi-major axis, or mean distance from Sun
e = eccentricity (0=circle, 0-1=ellipse, 1=parabola)
M = mean anomaly (0 at perihelion; increases uniformly with time)
Related orbital elements are:
w1 = N + w = longitude of perihelion
L = M + w1 = mean longitude
q = a*(1-e) = perihelion distance
Q = a*(1+e) = aphelion distance
P = a ^ 1.5 = orbital period (years if a is in AU, astronomical units)
T = Epoch_of_M - (M(deg)/360_deg) / P = time of perihelion
v = true anomaly (angle between position and perihelion)
E = eccentric anomaly
One Astronomical Unit (AU) is the Earth's mean distance to the
Sun, or 149.6 million km. When closest to the Sun, a planet is in
perihelion, and when most distant from the Sun it's in
aphelion. For the Moon, an artificial satellite, or any other
body orbiting the Earth, one says perigee and apogee
instead, for the points in orbit least and most distant from Earth.
To describe the position in the orbit, we use three angles: Mean
Anomaly, True Anomaly, and Eccentric Anomaly. They are all zero when the
planet is in perihelion:
Mean Anomaly (M): This angle increases uniformly over time, by
360 degrees per orbital period. It's zero at perihelion. It's
easily computed from the orbital period and the time since last
perihelion.
True Anomaly (v): This is the actual angle between the planet
and the perihelion, as seen from the central body (in this case the
Sun). It increases non-uniformly with time, changing most rapidly at
perihelion.
Eccentric Anomaly (E): This is an auxiliary angle used in
Kepler's Equation, when computing the True Anomaly from the Mean
Anomaly and the orbital eccentricity.
Note that for a circular orbit (eccentricity=0), these three angles
are all equal to each other.
Another quantity we will need is ecl, the obliquity of the
ecliptic, i.e. the "tilt" of the Earth's axis of rotation
(currently ca 23.4 degrees and slowly decreasing). First, compute the
"d" of the moment of interest (
Orbital elements of the Sun:
N = 0.0
i = 0.0
w = 282.9404 + 4.70935E-5 * d
a = 1.000000 (AU)
e = 0.016709 - 1.151E-9 * d
M = 356.0470 + 0.9856002585 * d
Orbital elements of the Moon:
N = 125.1228 - 0.0529538083 * d
i = 5.1454
w = 318.0634 + 0.1643573223 * d
a = 60.2666 (Earth radii)
e = 0.054900
M = 115.3654 + 13.0649929509 * d
Orbital elements of Mercury:
N = 48.3313 + 3.24587E-5 * d
i = 7.0047 + 5.00E-8 * d
w = 29.1241 + 1.01444E-5 * d
a = 0.387098 (AU)
e = 0.205635 + 5.59E-10 * d
M = 168.6562 + 4.0923344368 * d
Orbital elements of Venus:
N = 76.6799 + 2.46590E-5 * d
i = 3.3946 + 2.75E-8 * d
w = 54.8910 + 1.38374E-5 * d
a = 0.723330 (AU)
e = 0.006773 - 1.302E-9 * d
M = 48.0052 + 1.6021302244 * d
Orbital elements of Mars:
N = 49.5574 + 2.11081E-5 * d
i = 1.8497 - 1.78E-8 * d
w = 286.5016 + 2.92961E-5 * d
a = 1.523688 (AU)
e = 0.093405 + 2.516E-9 * d
M = 18.6021 + 0.5240207766 * d
Orbital elements of Jupiter:
N = 100.4542 + 2.76854E-5 * d
i = 1.3030 - 1.557E-7 * d
w = 273.8777 + 1.64505E-5 * d
a = 5.20256 (AU)
e = 0.048498 + 4.469E-9 * d
M = 19.8950 + 0.0830853001 * d
Orbital elements of Saturn:
N = 113.6634 + 2.38980E-5 * d
i = 2.4886 - 1.081E-7 * d
w = 339.3939 + 2.97661E-5 * d
a = 9.55475 (AU)
e = 0.055546 - 9.499E-9 * d
M = 316.9670 + 0.0334442282 * d
Orbital elements of Uranus:
N = 74.0005 + 1.3978E-5 * d
i = 0.7733 + 1.9E-8 * d
w = 96.6612 + 3.0565E-5 * d
a = 19.18171 - 1.55E-8 * d (AU)
e = 0.047318 + 7.45E-9 * d
M = 142.5905 + 0.011725806 * d
Orbital elements of Neptune:
N = 131.7806 + 3.0173E-5 * d
i = 1.7700 - 2.55E-7 * d
w = 272.8461 - 6.027E-6 * d
a = 30.05826 + 3.313E-8 * d (AU)
e = 0.008606 + 2.15E-9 * d
M = 260.2471 + 0.005995147 * d
Please note than the orbital elements of Uranus and Neptune as given
here are somewhat less accurate. They include a long period
perturbation between Uranus and Neptune. The period of the
perturbation is about 4200 years. Therefore, these elements should
not be expected to give results within the stated accuracy for more
than a few centuries in the past and into the future.
6. The position of the Moon and of the planets
First, compute the eccentric anomaly, E, from M, the mean anomaly,
and e, the eccentricity. As a first approximation, do (E and M in
degrees):
E = M + e*(180/pi) * sin(M) * ( 1.0 + e * cos(M) )
or, if E and M are in radians:
E = M + e * sin(M) * ( 1.0 + e * cos(M) )
If e, the eccentricity, is less than about 0.05-0.06, this
approximation is sufficiently accurate. If the eccentricity is
larger, set E0=E and then use this iteration formula (E and M in
degrees):
E1 = E0 - ( E0 - e*(180/pi) * sin(E0) - M ) / ( 1 - e * cos(E0) )
or (E and M in radians):
E1 = E0 - ( E0 - e * sin(E0) - M ) / ( 1 - e * cos(E0) )
For each new iteration, replace E0 with E1. Iterate until E0 and E1
are sufficiently close together (about 0.001 degrees). For comet
orbits with eccentricites close to one, a difference of less than
1E-4 or 1E-5 degrees should be required.
If this iteration formula won't converge, the eccentricity is
probably too close to one. Then you should instead use the formulae
for 9. Perturbations of the Moon
If the position of the Moon is computed, and one wishes a better
accuracy than about 2 degrees, the most important perturbations has
to be taken into account. If one wishes 2 arc minute accuracy, all
the following terms should be accounted for. If less accuracy is
needed, some of the smaller terms can be omitted.
First compute:
Ms, Mm Mean Anomaly of the Sun and the Moon
Nm Longitude of the Moon's node
ws, wm Argument of perihelion for the Sun and the Moon
Ls = Ms + ws Mean Longitude of the Sun (Ns=0)
Lm = Mm + wm + Nm Mean longitude of the Moon
D = Lm - Ls Mean elongation of the Moon
F = Lm - Nm Argument of latitude for the Moon
Add these terms to the Moon's longitude (degrees):
-1.274 * sin(Mm - 2*D) (the Evection)
+0.658 * sin(2*D) (the Variation)
-0.186 * sin(Ms) (the Yearly Equation)
-0.059 * sin(2*Mm - 2*D)
-0.057 * sin(Mm - 2*D + Ms)
+0.053 * sin(Mm + 2*D)
+0.046 * sin(2*D - Ms)
+0.041 * sin(Mm - Ms)
-0.035 * sin(D) (the Parallactic Equation)
-0.031 * sin(Mm + Ms)
-0.015 * sin(2*F - 2*D)
+0.011 * sin(Mm - 4*D)
Add these terms to the Moon's latitude (degrees):
-0.173 * sin(F - 2*D)
-0.055 * sin(Mm - F - 2*D)
-0.046 * sin(Mm + F - 2*D)
+0.033 * sin(F + 2*D)
+0.017 * sin(2*Mm + F)
Add these terms to the Moon's distance (Earth radii):
-0.58 * cos(Mm - 2*D)
-0.46 * cos(2*D)
All perturbation terms that are smaller than 0.01 degrees in
longitude or latitude and smaller than 0.1 Earth radii in distance
have been omitted here. A few of the largest perturbation terms even
have their own names! The Evection (the largest perturbation) was
discovered already by Ptolemy a few thousand years ago (the Evection
was one of Ptolemy's epicycles). The Variation and the Yearly
Equation were both discovered by Tycho Brahe in the 16'th
century.
The computations can be simplified by omitting the smaller
perturbation terms. The error introduced by this seldom exceeds the
sum of the amplitudes of the 4-5 largest omitted terms. If one only
computes the three largest perturbation terms in longitude and the
largest term in latitude, the error in longitude will rarley exceed
0.25 degrees, and in latitude 0.15 degrees.
10. Perturbations of Jupiter, Saturn and Uranus
The only planets having perturbations larger than 0.01 degrees are
Jupiter, Saturn and Uranus. First compute:
Mj Mean anomaly of Jupiter
Ms Mean anomaly of Saturn
Mu Mean anomaly of Uranus (needed for Uranus only)
Perturbations for Jupiter. Add these terms to the longitude:
-0.332 * sin(2*Mj - 5*Ms - 67.6 degrees)
-0.056 * sin(2*Mj - 2*Ms + 21 degrees)
+0.042 * sin(3*Mj - 5*Ms + 21 degrees)
-0.036 * sin(Mj - 2*Ms)
+0.022 * cos(Mj - Ms)
+0.023 * sin(2*Mj - 3*Ms + 52 degrees)
-0.016 * sin(Mj - 5*Ms - 69 degrees)
Perturbations for Saturn. Add these terms to the longitude:
+0.812 * sin(2*Mj - 5*Ms - 67.6 degrees)
-0.229 * cos(2*Mj - 4*Ms - 2 degrees)
+0.119 * sin(Mj - 2*Ms - 3 degrees)
+0.046 * sin(2*Mj - 6*Ms - 69 degrees)
+0.014 * sin(Mj - 3*Ms + 32 degrees)
For Saturn: also add these terms to the latitude:
-0.020 * cos(2*Mj - 4*Ms - 2 degrees)
+0.018 * sin(2*Mj - 6*Ms - 49 degrees)
Perturbations for Uranus: Add these terms to the longitude:
+0.040 * sin(Ms - 2*Mu + 6 degrees)
+0.035 * sin(Ms - 3*Mu + 33 degrees)
-0.015 * sin(Mj - Mu + 20 degrees)
The "great Jupiter-Saturn term" is the largest perturbation for both
Jupiter and Saturn. Its period is 918 years, and its amplitude is
0.332 degrees for Jupiter and 0.812 degrees for Saturn. These is
also a "great Saturn-Uranus term", period 560 years, amplitude 0.035
degrees for Uranus, less than 0.01 degrees for Saturn (and therefore
omitted). The other perturbations have periods between 14 and 100
years. One should also mention the "great Uranus-Neptune term",
which has a period of 4220 years and an amplitude of about one
degree. It is not included here, instead it is included in the
orbital elements of Uranus and Neptune.
For Mercury, Venus and Mars we can ignore all perturbations. For
Neptune the only significant perturbation is already included in the
orbital elements, as mentioned above, and therefore no further
perturbation terms need to be accounted for.
13. The Moon's topocentric position
The Moon's position, as computed earlier, is geocentric, i.e. as seen
by an imaginary observer at the center of the Earth. Real observers
dwell on the surface of the Earth, though, and they will see a
different position - the topocentric position. This position can
differ by more than one degree from the geocentric position. To
compute the topocentric positions, we must add a correction to the
geocentric position.
Let's start by computing the Moon's parallax, i.e. the apparent
size of the (equatorial) radius of the Earth, as seen from the Moon:
mpar = asin( 1/r )
where r is the Moon's distance in Earth radii. It's simplest to apply
the correction in horizontal coordinates (azimuth and altitude):
within our accuracy aim of 1-2 arc minutes, no correction need to be
applied to the azimuth. One need only apply a correction to the
altitude above the horizon:
alt_topoc = alt_geoc - mpar * cos(alt_geoc)
Sometimes one need to correct for topocentric position directly in
equatorial coordinates though, e.g. if one wants to draw on a star
map how the Moon passes in front of the Pleiades, as seen from some
specific location. Then we need to know the Moon's geocentric
Right Ascension and Declination (RA, Decl), the Local Sidereal
Time (LST), and our latitude (lat).
Our astronomical latitude (lat) must first be converted to a
geocentric latitude (gclat), and distance from the center of the Earth
(rho) in Earth equatorial radii. If we only want an approximate
topocentric position, it's simplest to pretend that the Earth is
a perfect sphere, and simply set:
gclat = lat, rho = 1.0
However, if we do wish to account for the flattening of the Earth,
we instead compute:
gclat = lat - 0.1924_deg * sin(2*lat)
rho = 0.99883 + 0.00167 * cos(2*lat)
Next we compute the Moon's geocentric Hour Angle (HA):
HA = LST - RA
where LST is our Local Sidereal Time, computed as:
LST = GMST0 + UT + LON/15
where UT is the Universal Time in hours, LON is the observer's
geographical longitude (east longitude positive, west negative).
GMST0 is the Greenwich Mean Sidereal Time at 0h UT, and is most
easily computed by adding, or subtracting, 180 degrees to Ls, the
Sun's mean longitude, and then converting from degrees to hours by
duvuding by 15:
GMST0 = ( Ls + 180_deg ) / 15
We also need an auxiliary angle, g:
g = atan( tan(gclat) / cos(HA) )
Now we're ready to convert the geocentric Right Ascention and
Declination (RA, Decl) to their topocentric values (topRA, topDecl):
topRA = RA - mpar * rho * cos(gclat) * sin(HA) / cos(Decl)
topDecl = Decl - mpar * rho * sin(gclat) * sin(g - Decl) / sin(g)
(Note that if decl is exactly 90 deg, cos(Decl) becomes zero and we get a
division by zero when computing topRA, but that formula breaks down
very close to the celestial poles anyway. Also if gclat is precisely
zero, g becomes zero too, and we get a division by zero when computing
topDecl. In that case, replace the formula for topDecl with
topDecl = Decl - mpar * rho * sin(-Decl) * cos(HA)
which is valid for gclat equal to zero; it can also be used for gclat
extremely close to zero).
This correction to topocentric position can also be applied to the
Sun and the planets. But since they're much farther away, the
correction becomes much smaller. It's largest for Venus at inferior
conjunction, when Venus' parallax is somewhat larger than 32 arc
seconds. Within our aim of obtaining a final accuracy of 1-2 arc
minutes, it might barely be justified to correct to topocentric
position when Venus is close to inferior conjunction, and perhaps
also when Mars is at a favourable opposition. But in all other cases
this correction can safely be ignored within our accuracy aim. We
only need to worry about the Moon in this case.
If you want to compute topocentric coordinates for the planets anyway,
you do it the same way as for the Moon, with one exception: the
parallax of the planet (ppar) is computed from this formula:
ppar = (8.794/3600)_deg / r
where r is the distance of the planet from the Earth, in astronomical
units.
15. The elongation and physical ephemerides of the planets
When we finally have completed our computation of the heliocentric
and geocentric coordinates of the planets, it could also be
interesting to know what the planet will look like. How large will
it appear? What's its phase and magnitude (brightness)? These
computations are much simpler than the computations of the
positions.
Let's start by computing the apparent diameter of the planet:
d = d0 / R
R is the planet's geocentric distance in astronomical units, and
d is the planet's apparent diameter at a distance of 1 astronomical
unit. d0 is of course different for each planet. The values
below are given in seconds of arc. Some planets have different
equatorial and polar diameters:
Mercury 6.74"
Venus 16.92"
Earth 17.59" equ 17.53" pol
Mars 9.36" equ 9.28" pol
Jupiter 196.94" equ 185.08" pol
Saturn 165.6" equ 150.8" pol
Uranus 65.8" equ 62.1" pol
Neptune 62.2" equ 60.9" pol
The Sun's apparent diameter at 1 astronomical unit is 1919.26". The
Moon's apparent diameter is:
d = 1873.7" * 60 / r
where r is the Moon's distance in Earth radii.
Two other quantities we'd like to know are the phase angle and the
elongation.
The phase angle tells us the phase: if it's zero the planet appears
"full", if it's 90 degrees it appears "half", and if it's 180 degrees
it appears "new". Only the Moon and the inferior planets (i.e.
Mercury and Venus) can have phase angles exceeding about 50 degrees.
The elongation is the apparent angular distance of the planet from
the Sun. If the elongation is smaller than ca 20 degrees, the planet
is hard to observe, and if it's smaller than ca 10 degrees it's
usually not possible to observe the planet.
To compute phase angle and elongation we need to know the planet's
heliocentric distance, r, its geocentric distance, R, and the
distance to the Sun, s. Now we can compute the phase angle, FV,
and the elongation, elong:
elong = acos( ( s*s + R*R - r*r ) / (2*s*R) )
FV = acos( ( r*r + R*R - s*s ) / (2*r*R) )
When we know the phase angle, we can easily compute the phase:
phase = ( 1 + cos(FV) ) / 2 = hav(180_deg - FV)
hav(FV) is the "haversine" of the phase angle. The "haversine" (or
"half versine") is an old and now obsolete trigonometric function;
it's defined as:
hav(x) = ( 1 - cos(x) ) / 2 = sin^2 (x/2)
As usual we must use a different procedure for the Moon. Since the
Moon is so close to the Earth, the procedure above would introduce
too big errors. Instead we use the Moon's ecliptic longitude and
latitude, mlon and mlat, and the Sun's ecliptic longitude, mlon, to
compute first the elongation, then the phase angle, of the Moon:
elong = acos( cos(slon - mlon) * cos(mlat) )
FV = 180_deg - elong
Finally we'll compute the magnitude (or brightness) of the planets.
Here we need to use a formula that's different for each planet. FV
is the phase angle (in degrees), r is the heliocentric and R the
geocentric distance (both in AU):
Mercury: -0.36 + 5*log10(r*R) + 0.027 * FV + 2.2E-13 * FV**6
Venus: -4.34 + 5*log10(r*R) + 0.013 * FV + 4.2E-7 * FV**3
Mars: -1.51 + 5*log10(r*R) + 0.016 * FV
Jupiter: -9.25 + 5*log10(r*R) + 0.014 * FV
Saturn: -9.0 + 5*log10(r*R) + 0.044 * FV + ring_magn
Uranus: -7.15 + 5*log10(r*R) + 0.001 * FV
Neptune: -6.90 + 5*log10(r*R) + 0.001 * FV
Moon: +0.23 + 5*log10(r*R) + 0.026 * FV + 4.0E-9 * FV**4
** is the power operator, thus FV**6 is the phase angle (in degrees)
raised to the sixth power. If FV is 150 degrees then FV**6 becomes
ca 1.14E+13, which is a quite large number.
For the Moon, we also need the heliocentric distance, r, and
geocentric distance, R, in AU (astronomical units). Here r can be set
equal to the Sun's geocentric distance in AU. The Moon's geocentric
distance, R, previously computed i Earth radii, must be converted
to AU's - we do this by multiplying by sin(17.59"/2) = 1/23450. Or
we could modify the magnitude formula for the Moon so it uses
r in AU's and R in Earth radii:
Moon: -21.62 + 5*log10(r*R) + 0.026 * FV + 4.0E-9 * FV**4
Saturn needs special treatment due to its rings: when Saturn's
rings are "open" then Saturn will appear much brighter than when
we view Saturn's rings edgewise. We'll compute ring_mang like
this:
ring_magn = -2.6 * sin(abs(B)) + 1.2 * (sin(B))**2
Here B is the tilt of Saturn's rings which we also need to compute.
Then we start with Saturn's geocentric ecliptic longitude and
latitude (los, las) which we've already computed. We also need the
tilt of the rings to the ecliptic, ir, and the "ascending node" of
the plane of the rings, Nr:
ir = 28.06_deg
Nr = 169.51_deg + 3.82E-5_deg * d
Here d is our "day number" which we've used so many times before.
Now let's compute the tilt of the rings:
B = asin( sin(las) * cos(ir) - cos(las) * sin(ir) * sin(los-Nr) )
This concludes our computation of the magnitudes of the planets.