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.
The main steps in calculating the RA and DEC of a planet from the mean elements are:
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)
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:
where div denotes integer division (for example, 17 div 4 = 4).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
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:
where cy is the number of centuries since J2000 and is computed using: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
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:
To give more precise computer calculations, the value of M should be in the range 0...360.M = L - w
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:
Then iterate using:E = M + e·sin(M)·(1.0 + e·cos(M))
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:
And ensure that the result is in the range 0...360V = 2·arctan(sqrt((1 + e)/(1 - e))·tan(0.5·E))
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)
Zeno 1.2 is an interpreter for the Zeno programming language. It is an easy to learn and is suitable for educational purposes.
const pi : real := 4 * arctan(1)
const degs : real := 180 / pi
const rads : real := pi / 180
const eps : real := 1.0e-12
type elem : record
a : real % semi-major axis [AU]
e : real % eccentricity of orbit
i : real % inclination of orbit [deg]
O : real % longitude of the ascending node [deg]
w : real % longitude of perihelion [deg]
L : real % mean longitude [deg]
end record
type coord : record
ra, dec, rvec : real
end record
program
var year, month, day, hour, mins, p : int
% get the date and universal time
year := get_int( "year" )
month := get_int( "month" )
day := get_int( "day" )
hour := get_int( "hour" )
mins := get_int( "minute" )
% compute day number for date/time
var d : real := day_number( year, month, day, hour, mins )
put "Date: ", month, "/", day, "/", year, " UT: ", hour, ":", mins
put "days since J2000: ", d
put ""
put "Object RA DEC Distance"
put "----------------------------------------"
% compute location of objects
for p := 1...9 do
var obj : coord
get_coord( obj, p, d )
put pname(p)...
put " ", ha2str(obj.ra*degs)...
put " ", dec2str(obj.dec*degs)...
put " ", obj.rvec:10:6
end for
put ""
end program
% day number to/from J2000 (Jan 1.5, 2000)
function day_number( y, m, d, hour, mins : int ) : real
var rv : real
var h : real := hour + mins / 60
rv := 367*y - 7*(y + (m + 9) div 12) div 4
+ 275*m div 9 + d - 730531.5 + h / 24
return rv
end function
% compute RA, DEC, and distance of planet-p for day number-d
% result returned in structure obj
procedure get_coord( var obj : coord, p : int, d : real )
var planet : elem
mean_elements( planet, p, d )
var ap : real := planet.a
var ep : real := planet.e
var ip : real := planet.i
var op : real := planet.O
var pp : real := planet.w
var lp : real := planet.L
var earth : elem
mean_elements( earth, 3, d )
var ae : real := earth.a
var ee : real := earth.e
var ie : real := earth.i
var oe : real := earth.O
var pe : real := earth.w
var le : real := earth.L
% position of Earth in its orbit
var me : real := mod2pi(le - pe)
var ve : real := true_anomaly(me, ee)
var re : real := ae*(1 - ee*ee) / (1 + ee*cos(ve))
% heliocentric rectangular coordinates of Earth
var xe : real := re*cos(ve + pe)
var ye : real := re*sin(ve + pe)
var ze : real := 0
% position of planet in its orbit
var mp : real := mod2pi(lp - pp)
var vp : real := true_anomaly(mp, planet.e)
var rp : real := ap*(1 - ep*ep) / (1 + ep*cos(vp))
% heliocentric rectangular coordinates of planet
var xh : real := rp*(cos(op)*cos(vp + pp - op) - sin(op)*sin(vp + pp - op)*cos(ip))
var yh : real := rp*(sin(op)*cos(vp + pp - op) + cos(op)*sin(vp + pp - op)*cos(ip))
var zh : real := rp*(sin(vp + pp - op)*sin(ip))
if p = 3 then % earth --> compute sun
xh := 0
yh := 0
zh := 0
end if
% convert to geocentric rectangular coordinates
var xg : real := xh - xe
var yg : real := yh - ye
var zg : real := zh - ze
% rotate around x axis from ecliptic to equatorial coords
var ecl : real := 23.439281*rads %value for J2000.0 frame
var xeq : real := xg
var yeq : real := yg*cos(ecl) - zg*sin(ecl)
var zeq : real := yg*sin(ecl) + zg*cos(ecl)
% find the RA and DEC from the rectangular equatorial coords
obj.ra := mod2pi( arctanxy(xeq, yeq) )
obj.dec := arctan(zeq / sqrt(xeq*xeq + yeq*yeq))
obj.rvec := sqrt(xeq*xeq + yeq*yeq + zeq*zeq)
end procedure
% Compute the elements of the orbit for planet-i at day number-d
% result is returned in structure p
procedure mean_elements( var p : elem, i : int, d : real )
var cy : real := d / 36525 % centuries since J2000
case i of
value 1: % Mercury
p.a := 0.38709893 + 0.00000066*cy
p.e := 0.20563069 + 0.00002527*cy
p.i := ( 7.00487 - 23.51*cy/3600)*rads
p.O := (48.33167 - 446.30*cy/3600)*rads
p.w := (77.45645 + 573.57*cy/3600)*rads
p.L := mod2pi((252.25084 + 538101628.29*cy/3600)*rads)
value 2: % Venus
p.a := 0.72333199 + 0.00000092*cy
p.e := 0.00677323 - 0.00004938*cy
p.i := ( 3.39471 - 2.86*cy/3600)*rads
p.O := ( 76.68069 - 996.89*cy/3600)*rads
p.w := (131.53298 - 108.80*cy/3600)*rads
p.L := mod2pi((181.97973 + 210664136.06*cy/3600)*rads)
value 3: % Earth/Sun
p.a := 1.00000011 - 0.00000005*cy
p.e := 0.01671022 - 0.00003804*cy
p.i := ( 0.00005 - 46.94*cy/3600)*rads
p.O := (-11.26064 - 18228.25*cy/3600)*rads
p.w := (102.94719 + 1198.28*cy/3600)*rads
p.L := mod2pi((100.46435 + 129597740.63*cy/3600)*rads)
value 4: % Mars
p.a := 1.52366231 - 0.00007221*cy
p.e := 0.09341233 + 0.00011902*cy
p.i := ( 1.85061 - 25.47*cy/3600)*rads
p.O := ( 49.57854 - 1020.19*cy/3600)*rads
p.w := (336.04084 + 1560.78*cy/3600)*rads
p.L := mod2pi((355.45332 + 68905103.78*cy/3600)*rads)
value 5: % Jupiter
p.a := 5.20336301 + 0.00060737*cy
p.e := 0.04839266 - 0.00012880*cy
p.i := ( 1.30530 - 4.15*cy/3600)*rads
p.O := (100.55615 + 1217.17*cy/3600)*rads
p.w := ( 14.75385 + 839.93*cy/3600)*rads
p.L := mod2pi((34.40438 + 10925078.35*cy/3600)*rads)
value 6: % Saturn
p.a := 9.53707032 - 0.00301530*cy
p.e := 0.05415060 - 0.00036762*cy
p.i := ( 2.48446 + 6.11*cy/3600)*rads
p.O := (113.71504 - 1591.05*cy/3600)*rads
p.w := ( 92.43194 - 1948.89*cy/3600)*rads
p.L := mod2pi((49.94432 + 4401052.95*cy/3600)*rads)
value 7: % Uranus
p.a := 19.19126393 + 0.00152025*cy
p.e := 0.04716771 - 0.00019150*cy
p.i := ( 0.76986 - 2.09*cy/3600)*rads
p.O := ( 74.22988 - 1681.40*cy/3600)*rads
p.w := (170.96424 + 1312.56*cy/3600)*rads
p.L := mod2pi((313.23218 + 1542547.79*cy/3600)*rads)
value 8: % Neptune
p.a := 30.06896348 - 0.00125196*cy
p.e := 0.00858587 + 0.00002510*cy
p.i := ( 1.76917 - 3.64*cy/3600)*rads
p.O := (131.72169 - 151.25*cy/3600)*rads
p.w := ( 44.97135 - 844.43*cy/3600)*rads
p.L := mod2pi((304.88003 + 786449.21*cy/3600)*rads)
value 9: % Pluto
p.a := 39.48168677 - 0.00076912*cy
p.e := 0.24880766 + 0.00006465*cy
p.i := ( 17.14175 + 11.07*cy/3600)*rads
p.O := (110.30347 - 37.33*cy/3600)*rads
p.w := (224.06676 - 132.25*cy/3600)*rads
p.L := mod2pi((238.92881 + 522747.90*cy/3600)*rads)
value:
assert false
end case
end procedure
% compute the true anomaly from mean anomaly using iteration
% M - mean anomaly in radians
% e - orbit eccentricity
function true_anomaly( M, e : real ) : real
var V, E, E1 : real
% initial approximation of eccentric anomaly
E := M + e*sin(M)*(1.0 + e*cos(M))
% iterate to improve accuracy
repeat
E1 := E
E := E1 - (E1 - e*sin(E1) - M)/(1 - e*cos(E1))
until abs( E - E1 ) < eps
% convert eccentric anomaly to true anomaly
V := 2*arctan(sqrt((1 + e)/(1 - e))*tan(0.5*E))
% modulo 2pi
if (V < 0) then
V := V + (2*pi)
end if
return V
end function
% converts hour angle in degrees into hour angle string
function ha2str( x : real ) : string
assert (0 <= x) and (x < 360) % assure valid range
var ra : real := x / 15 % degrees to hours
var h : int := floor( ra )
var m : real := 60 * ( ra - h )
return intstr( h, 3 ) & "h " & frealstr( m, 4, 1 ) & "m"
end function
% converts declination angle in degrees into string
function dec2str( x : real ) : string
assert (-90 <= x) and (x <= +90) % assure valid range
var dec : real := abs( x ) % work with absolute value
var d : int := floor( dec )
var m : real := 60 * ( dec - d )
return intstr( sgn(x)*d, 3 ) & "° " & frealstr( m, 4, 1 ) & "'"
end function
% get object name
function pname( i : int ) : string
case i of
value 1 : return "Mercury"
value 2 : return "Venus "
value 3 : return "Sun "
value 4 : return "Mars "
value 5 : return "Jupiter"
value 6 : return "Saturn "
value 7 : return "Uranus "
value 8 : return "Neptune"
value 9 : return "Pluto "
end case
return "not in range"
end function
% prompts for input of an int
function get_int( prompt : string ) : int
var n : int
put prompt...
get n
return n
end function
% return the integer part of a number
function abs_floor( x : real ) : int
var r : int
if x >= 0.0 then
r := floor( x )
else
r := ceil( x )
end if
return r
end function
% return an angle in the range 0 to 2pi
function mod2pi( x : real ) : real % radians
var b : real := x / (2*pi)
var a : real := (2*pi)*(b - abs_floor(b))
if a < 0 then
a := (2*pi) + a
end if
return a
end function
% returns absolute value of real number
function abs( x : real ) : real
if x < 0 then
return -x
end if
return x
end function
% returns sign of real number as int
function sgn( x : real ) : int
if x < 0 then
return -1
end if
return +1
end function
year? 2004 month? 5 day? 1 hour? 0 minute? 0 Date: 5/1/2004 UT: 0:0 days since J2000: 1581.5 Object RA DEC Distance ---------------------------------------- Mercury 1h 20.5m 6° 34.4' 0.633271 Venus 5h 20.1m 27° 43.9' 0.462291 Sun 2h 33.9m 15° 5.9' 1.007611 Mars 5h 42.0m 24° 36.1' 2.166172 Jupiter 10h 44.1m 9° 28.4' 4.879948 Saturn 6h 37.7m 22° 45.8' 9.527284 Uranus 22h 32.5m -9° 58.4' 20.458105 Neptune 21h 11.0m -16° 18.3' 30.133788 Pluto 17h 26.7m -14° 17.4' 30.032601