| home | contents | send comment | send link | add bookmark |

Planet position calculation using mean orbital elements

by Stephen R. Schmitt

Enter your location:
latitude ° ' North South
longitude ° ' East West

today:

Enter the date and time:
mm/dd/yyyy hr:min:sec
/ /      : :


Contents

  1. About
  2. The source code
  3. Discussion

1. About

This Java Script calculator computes the right ascension and declination of the sun and the planets. It transforms these coordinates into the horizon coordinates of altitude and azimuth at your location. To operate the calculator, enter the latitude and longitude of the observing site and the local date and time. The program assumes that your computer's time and zone is set correctly. Press the Compute button to obtain the solution. On invalid entries, a popup window will display an error message.

Return to Contents


2. The source code

The Java Script source code for this program can be viewed by using the View|Source command of your web browser.

You may use or modify this source code in any way you find useful, provided that you agree that the author has no warranty, obligations or liability. You must determine the suitability of this source code for your use.

Return to Contents


3. Discussion

Mean orbital elements

The orbits of the major planets can be modeled as ellipses with the Sun at one focus. The effect of gravitational interactions between the planets perturbs these orbits so that an ellipse is not a exact match with a true orbit. Six numbers, the mean orbital elements, specify an elliptical orbit. Mean orbital elements average the effects of gravitational forces between planets. Calculation of a planets position based on these mean elements can be inaccurate by a few minutes of arc.

The position of a planet (the word originally meant wandering star) varies with time. The daily motion changes the mean longitude by the average number of degrees the planet moves in one (mean solar) day. The other elements change slowly with time. They are modeled using power series expansions of centuries from some fundamental epoch. Here, we use the elements with their linear rates of change from the epoch J2000 (12:00 UT, Jan 1, 2000).

Planet positions are computed in the Equatorial coordinate system as right ascension (RA) and declination (DEC). These give the coordinates of the planet with respect to the fixed stars. The origin for RA is the vernal equinox. Because the orientation of the Earth's axis is changing slowly with time, celestial coordinates must always be referred to an epoch, or date. By using orbital elements referred to epoch J2000, the orbits of the planets are described in a coordinate system that is based on the position the vernal equinox will have at J2000. The effect of nutation (the Earth's axis is nodding) is ignored since positions are relative to the mean ecliptic of J2000. The aberration effect caused by the finite speed of light is also ignored.

Computation steps

The main steps in calculating the RA and DEC of a planet from the mean elements are:

  1. Find the day number or time since the date of the elements
  2. Find the average (or mean) orbital elements of the planet
  3. Find the true anomaly, the angle to the planet from perihelion
    (helio means Sun)
  4. Find the heliocentric radius, distance to planet from sun
  5. Find the heliocentric ecliptic rectangular coordinates of the planet
  6. Find the heliocentric ecliptic rectangular coordinates of the Earth
  7. Transform the heliocentric coordinates to geocentric coordinates
    (geo means Earth)
  8. Transform the geocentric ecliptic coordinates to equatorial coordinates
  9. Calculate the RA and DEC and the distance to the planet from Earth

Notation

Elements

i - inclination (degrees), angle between the plane of the ecliptic (the plane of Earth's orbit about the Sun) and the plane of the planets orbit
O - longitude of ascending node (degrees), the position in the orbit where the elliptical path of the planet passes through the plane of the ecliptic, from below the plane to above the plane
w - longitude of perihelion (degrees), the position in the orbit where the planet is closest to the Sun
a - mean distance (AU), the value of the semi-major axis of the orbit (AU - Astronomical Unit - average Sun to Earth distance)
e - eccentricity of the ellipse which describes the orbit (dimensionless)
L - mean longitude (degrees), the position of the planet in the orbit

Calculated quantities

M - mean anomaly (degrees)
V - true anomaly (degrees)
R - heliocentric radius (AU)
X,Y,Z - rectangular coordinates (AU)
RA - right ascension (hour angle)
DEC - declination (degrees)

Position of the planet in its orbit

The day number

The day number is used to compute values of all time varying quantities. It is a measure of the number of days, hours, minutes, and seconds since epoch J2000. This can be calculated from the date using:

h = hour + minute/60 + second/3600
D = 367*year - 7*(year + (month + 9) div 12) div 4 + 275*month div 9 + day - 730531.5 + h/24
where div denotes integer division (for example, 17 div 4 = 4).

The average (or mean) orbital elements of the planet

Each of the orbital elements varies with time. The most rapidly changing element is the mean longitude, which changes due to the planets orbital motion about the Sun. The other elements change slowly. Gravitational forces between the planets cause these slow changes. For example, the orbital elements of the planet Mars are computed from:

a = 1.52366231 -  0.00007221*cy
e = 0.09341233 +  0.00011902*cy
i =    1.85061 -       25.47*cy/3600
O =   49.57854 -     1020.19*cy/3600
w =  336.04084 +     1560.78*cy/3600
L =  355.45332 + 68905103.78*cy/3600
where cy is the number of centuries since J2000 and is computed using:
cy = D/36525.0

The mean anomaly of the planet

The mean anomaly gives the planet's angular position for a circular orbit with radius equal to the semi major axis. It is computed directly from the elements using:

M = L - w 
To give more precise computer calculations, the value of M should be in the range 0...360.

The true anomaly of the planet

Kepler's second law states that the radius vector of a planet sweeps out equal areas in equal times. The planet must speed up and slow down in its orbit. The true anomaly gives the planets actual angular position in its orbit. It is the angle (at the Sun) between perihelion of the orbit and the current location of the planet. To obtain its value, first compute the eccentric anomaly, E, from M, and the eccentricity, e. As a first approximation:

E = M + e·sin(M)·(1.0 + e·cos(M))
Then iterate using:
E' = E
          E' - e·sin(E') - M
E  = E' - ------------------
            1 - e·cos(E')
until the magnitude of E - E' is sufficiently close to zero. Finally, convert the eccentric anomaly to the true anomaly using:
V  = 2·arctan(sqrt((1 + e)/(1 - e))·tan(0.5·E))
And ensure that the result is in the range 0...360

The heliocentric radius of the planet

The heliocentric radius is the distance from the focus of the ellipse (i.e. the Sun) to the planet. It is given by a formula based on the geometry of an ellipse:

     a·(1 - e2)
R = ------------
    1 + e·cos(V)

The heliocentric coordinates of the planet

By using the true anomaly, the heliocentric radius, and some of the orbital elements, the formulas below compute the heliocentric coordinates in the plane of the ecliptic.

XH = R·[cos(O)·cos(V + w - O) - sin(O)·sin(V + w - O)·cos(i)]
YH = R·[sin(O)·cos(V + w - O) + cos(O)·sin(V + w - O)·cos(i)]
ZH = R·[sin(V + w - O)·sin(i)]

The heliocentric coordinates of Earth

These are computed using the same method as above and are denoted as: XE YE ZE.

Geocentric ecliptic coordinates of the planet

The origin of the coordinate system is changed from the Sun to the Earth by subtracting the Earth's heliocentric coordinates from those of the planet:

XG = XH - XE
YG = YH - YE
ZG = ZH - ZE  

The geocentric equatorial coordinates of the planet

To change the coordinate system from geocentric ecliptic to geocentric equatorial, rotate around the X-axis by an angle equal to the obliquity of the ecliptic, ecl. The X-axis points towards the vernal equinox, where the Sun crosses the celestial equator in the spring.

ecl = 23.439281   the value of the obliquity of the ecliptic for J2000
XEQ = XG
YEQ = YG·cos(ecl) - ZG·sin(ecl)
ZEQ = YG·sin(ecl) + ZG·cos(ecl)

The RA and DEC and the distance

The geocentric equatorial coordinates are transformed into right ascension (RA) and declination (DEC) using the formulas:

       / arctan(YEQ / XEQ)
RA  = <  arctan(YEQ / XEQ) + 180°  if (XEQ < 0)
       \ arctan(YEQ / XEQ) + 360°  if (XEQ > 0) and (YEQ < 0)  

DEC = arctan( ZEQ / sqrt(XEQ2 + YEQ2) ) 
RA can be converted from degrees into hour angle by dividing by 15.

The distance from the planet to Earth is computed using:

Distance = sqrt(XEQ2 + YEQ2 + ZEQ2)          

Equatorial coordinates

By extending the lines of latitude and longitude outward from the Earth and onto the inside of the celestial sphere we get the equatorial coordinate system. The coordinates of stars, planets, and other celestial objects corresponding to latitude and longitude are declination (DEC) and right ascension (RA).

The declination of an object is its angle in degrees, minutes, and seconds of arc above or below the celestial equator. The right ascension is the angle between an object and the location of the vernal equinox (First Point in Aries) measured eastward along the celestial equator in hours, minutes, and seconds of sidereal time. Since the location of the vernal equinox changes due to the precession of the Earth's axis of rotation, coordinates must be given with reference to a date or epoch.

Right ascension is given in time units. One hour corresponds to 1/24 of a circle, or 15° of arc. As the Earth rotates, the sky moves to the West by about 1 hour of right ascension during each hour of clock time or exactly one hour of sidereal time. The Earth makes one full revolution in about 23 hours and 56 minutes of clock time or 24 hours of sidereal time. Sidereal time corresponds to the right ascension of the zenith, the point in the sky directly overhead.

For example, the coordinates of the star Regulus (Leo α) for epoch J2000 are:

RA:   10h 08m 22.3s
DEC: +11° 58' 02"
When the local sidereal time is 10h 08m 22.3s, it would be on the local meridian.

Horizon coordinates: azimuth and altitude

This is a local coordinate system to use for locating objects in the night sky as seen from a point on the Earth's surface. Azimuth is the angle of a celestial object around the sky from north. It is measure along the horizon in from North 0° through East 90°, South 180°, West 270° and back to North. Altitude is the complement of the zenith angle, which is the angle from the local meridian to the hour circle of object being observed. An object directly overhead would have an altitude of 90°. An object with a calculated altitude of 0° may not appear exactly on the horizon due to the refraction of light through the atmosphere.

Coordinate transformation

The azimuth (AZ) and altitude (ALT) of an object in the sky can be calculated easily using the date, universal time (UT), and the latitude (LAT) and longitude (LON) of the observing site and the right ascension (RA) and declination (DEC) of the object. All coordinates are expressed in degrees in the range 0° to 360°, so that trigonometric functions can be used for coordinate conversion.

Local Mean Sidereal Time

The mean sidereal time (MST) is calculated from a polynomial function of UT since epoch J2000. This formula gives MST, the sidereal time at the Greenwich meridian (at longitude 0°) in degrees. To get local mean sidereal time (LMST), add longitude if East or subtract longitude if West.

MST = f(UT)
LMST = MST + LON 

Hour Angle

The hour angle is the angle between an observer's meridian projected onto the celestial sphere and the right ascension of a celestial body. It is used in coordinate conversion.

HA = LMST - RA 

HA and DEC to ALT and AZ

Using the RA, DEC and HA for the object, and the Latitude (LAT) of the observing site, the following formulas give the ALT and AZ of the object at the time and longitude that was used to calculate HA.

sin(ALT) = sin(DEC)·sin(LAT) + cos(DEC)·cos(LAT)·cos(HA)

           sin(DEC) - sin(ALT)·sin(LAT)
cos(A)   = ----------------------------
                cos(ALT)·cos(LAT)
If sin(HA) is negative, then AZ = A, otherwise AZ = 360 - A

References

  1. US Naval Observatory, Explanatory Supplement to the Astronomical Almanac, 1992

Return to Contents


AbCd Classics - free on-line books


Copyright © 2004, Stephen R. Schmitt