You can describe the orbit of a major planet as an ellipse, with the Sun at one focus. This simple picture ignores the effect of gravitational interactions between the planets, which are described as perturbing the orbit. The elliptical orbit is specified by six numbers called the orbital elements. These orbital elements can be osculating elements or mean elements. Using the mean elements and the method discussed here will result in worst case positions of Mars correct to within 6 arc minutes, of Jupiter to within 12 arc minutes, and of Neptune to within two arcminutes on the sky for 50 years either side of J2000. 'Most of the time' the errors will be much lower.
Osculating elements are accurate for some date, and have the effects of the perturbations of the orbit built into them for that date. Positions calculated from the ellipse defined by a set of osculating elements become less accurate at dates further from the date of the elements. The 'planet positions using elliptical orbit' page on this site discusses the accuracy of the osculating elements in some detail.
Mean elements are designed to reflect the average orbit of the
planet over a range of time, and do not include adjustments for
perturbations based on any single date. Slow changes in the parameters can be
modelled using power series expansions of the time (usually in Julian
centuaries) from some fundamental epoch. Here, I have used the mean elements
given in the Explanatory Supplement to the Astronomical Almanac,
1992, page 318, table 5.8.1. The elements are given with their rates of change,
i.e. the linear term in time from J2000. The values and the daily rates of
change can be found in the QBASIC program below, or in a table exported
from my spreadsheet as a graphic.
The method of calulation is identical to that using the osculating elements. Please refer to the 'planet positions using elliptical orbit' page on this site for details and numerical examples.
' This program finds the geocentric positions of the major
' and planets for the mean equinox and equator of
' J2000. The program uses a simple Kepler ellipse, but with
' mean elements which change slightly with the time since
' J2000. See 'Explanatory supplement to the Astronomical
' Almanac' 1992, page 316, table 5.8.1
'
' Work in double precision and define some constants
DEFDBL A-Z
pr$ = "\ \#####.##"
pr2$ = "\ \###.#######"
pi = 4 * ATN(1)
tpi = 2 * pi
degs = 180 / pi
rads = pi / 180
DIM pname$(9)
pname$(1) = "Mercury"
pname$(2) = "Venus"
pname$(3) = "Earth"
pname$(4) = "Mars"
pname$(5) = "Jupiter"
pname$(6) = "Saturn"
pname$(7) = "Uranus"
pname$(8) = "Neptune"
pname$(9) = "Pluto"
' list of elements el()
DIM el(9 * 6 + 2)
' define some functions
'
' days to J2000 - FNday may only work from 1900 to 2100 - see Meeus
DEF FNday (y, m, d, h) = 367 * y - 7 * (y + (m + 9) \ 12) \ 4 + 275 * m \ 9 + d - 730531.5 + h / 24
'
' define some arc cos and arc sin functions and a modified inverse
' tangent function
DEF FNacos (x)
s = SQR(1 - x * x)
FNacos = ATN(s / x)
END DEF
DEF FNasin (x)
c = SQR(1 - x * x)
FNasin = ATN(x / c)
END DEF
' the atn2 function below returns an angle in the range 0 to two pi
' depending on the signs of x and y.
DEF FNatn2 (y, x)
a = ATN(y / x)
IF x < 0 THEN a = a + pi
IF (y < 0) AND (x > 0) THEN a = a + tpi
FNatn2 = a
END DEF
' the function below returns the true integer part,
' even for negative numbers
DEF FNipart (x) = SGN(x) * INT(ABS(x))
' the function below returns an angle in the range
' 0 to two pi
DEF FNrange (x)
b = x / tpi
a = tpi * (b - FNipart(b))
IF a < 0 THEN a = tpi + a
FNrange = a
END DEF
' The next function finds the elements of the orbits of all the
' planets. The function just returns 'true' but alters the values
' in the global array el(). Inelegant but it allows me to compile
' using the FirstBas basic compiler.
DEF FNelement (d)
'
' Mercury
'
el(1) = (7.00487 - .000000178797# * d) * rads
el(2) = (48.33167 - .0000033942# * d) * rads
el(3) = (77.45645 + .00000436208# * d) * rads
el(4) = .38709893# + .0000000000180698# * d
el(5) = .20563069# + .000000000691855# * d
el(6) = FNrange(rads * (252.25084# + 4.092338796# * d))
'
' Venus
'
el(7) = (3.39471 - .0000000217507# * d) * rads
el(8) = (76.68069 - .0000075815# * d) * rads
el(9) = (131.53298# - .000000827439# * d) * rads
el(10) = .72333199# + .0000000000251882# * d
el(11) = .00677323# - .00000000135195# * d
el(12) = FNrange(rads * (181.97973# + 1.602130474# * d))
'
' Earth
'
el(13) = (.00005 - .000000356985# * d) * rads
el(14) = (-11.26064 - .00013863# * d) * rads
el(15) = (102.94719# + .00000911309# * d) * rads
el(16) = 1.00000011# - 1.36893D-12 * d
el(17) = .01671022# - .00000000104148# * d
el(18) = FNrange(rads * (100.46435# + .985609101# * d))
'
' Mars
'
el(19) = (1.85061 - .000000193703# * d) * rads
el(20) = (49.57854 - 7.758700000000001D-06 * d) * rads
el(21) = (336.04084# + .00001187# * d) * rads
el(22) = 1.52366231# - .000000001977# * d
el(23) = .09341233# - .00000000325859# * d
el(24) = FNrange(rads * (355.45332# + .524033035# * d))
'
' Jupiter
'
el(25) = (1.3053 - .0000000315613# * d) * rads
el(26) = (100.55615# + 9.256750000000001D-06 * d) * rads
el(27) = (14.75385# + .00000638779# * d) * rads
el(28) = 5.20336301# + .0000000166289# * d
el(29) = .04839266# - .00000000352635# * d
el(30) = FNrange(rads * (34.40438# + 8.308676199999999D-02 * d))
'
' Saturn
'
el(31) = (2.48446 + .0000000464674# * d) * rads
el(32) = (113.71504# - .0000121# * d) * rads
el(33) = (92.43194# - .0000148216# * d) * rads
el(34) = 9.53707032# - 8.255439999999999D-08 * d
el(35) = .0541506# - .0000000100649# * d
el(36) = FNrange(rads * (49.94432# + .033470629# * d))
'
' Uranus
'
el(37) = (.76986 - .0000000158947# * d) * rads
el(38) = (74.22987999999999# + .0000127873# * d) * rads
el(39) = (170.96424# + 9.982200000000001D-06 * d) * rads
el(40) = 19.19126393# + .0000000416222# * d
el(41) = .04716771# - .00000000524298# * d
el(42) = FNrange(rads * (313.23218# + .011731294# * d))
'
' Neptune
'
el(43) = (1.76917 - .0000000276827# * d) * rads
el(44) = (131.72169# - .0000011503# * d) * rads
el(45) = (44.97135 - .00000642201# * d) * rads
el(46) = 30.06896348# - .0000000342768# * d
el(47) = 8.585870000000001D-03 + .000000000688296# * d
el(48) = FNrange(rads * (304.88003# + 0.0059810572# * d))
'
' Pluto
'
el(49) = (17.14175 + 8.418889999999999D-08 * d) * rads
el(50) = (110.30347# - .0000002839# * d) * rads
el(51) = (224.06676# - .00000100578# * d) * rads
el(52) = 39.48168677# - .0000000210574# * d
el(53) = .24880766# + .00000000177002# * d
el(54) = FNrange(rads * (238.92881# + .013831834# * d))
'
' return a dummy value to finish function
FNelement = 1
END DEF
'
DEF FNkep (m, ecc, eps)
'
' returns the true anomaly given
' m - the mean anomaly in radians
' ecc - the eccentricity of the orbit
' eps - the convergence paramter (8 or 9 is usually fine
' 12 or 14 for very accurate work)
'
e = m
delta = .05#
DO WHILE ABS(delta) >= 10 ^ -eps
delta = e - ecc * SIN(e) - m
e = e - delta / (1 - ecc * COS(e))
LOOP
v = 2 * ATN(((1 + ecc) / (1 - ecc)) ^ .5 * TAN(.5 * e))
IF v < 0 THEN v = v + tpi
FNkep = v
END DEF
'
DEF FNdegmin (x)
' cosmetic function returns angular values as a made up decimal
' number - ddd.mm - the digits after the decimal point are the
' minutes.
a = FNipart(x)
b = x - a
e = FNipart(60 * b)
' deal with carry on minutes
IF e >= 60 THEN
e = 0
a = a + 1
END IF
FNdegmin = a + e / 100
END DEF
'
CLS
'
' get the date and planet number from the user
'
INPUT " year : ", y
INPUT " month : ", m
INPUT " day : ", day
INPUT "hour UT : ", h
INPUT " minute : ", mins
h = h + mins / 60
INPUT " planet : ", p
d = FNday(y, m, day, h)
scrap = FNelement(d)
PRINT USING pr$; " days : "; d
'
' get the osculating elements from the list
' using letters instead of the array element
' makes the program easier to read.
'
q = 6 * (p - 1)
ip = el(q + 1)
op = el(q + 2)
pp = el(q + 3)
ap = el(q + 4)
ep = el(q + 5)
lp = el(q + 6)
ie = el(13)
oe = el(14)
pe = el(15)
ae = el(16)
ee = el(17)
le = el(18)
'
' now find position of Earth in orbit
'
me = FNrange(le - pe)
ve = FNkep(me, ee, 12)
re = ae * (1 - ee * ee) / (1 + ee * COS(ve))
xe = re * COS(ve + pe)
ye = re * SIN(ve + pe)
ze = 0
'
' and position of planet in its orbit
'
mp = FNrange(lp - pp)
vp = FNkep(mp, ep, 12)
rp = ap * (1 - ep * ep) / (1 + ep * COS(vp))
'
' heliocentric rectangular coordinates of planet
'
xh = rp * (COS(op) * COS(vp + pp - op) - SIN(op) * SIN(vp + pp - op) * COS(ip))
yh = rp * (SIN(op) * COS(vp + pp - op) + COS(op) * SIN(vp + pp - op) * COS(ip))
zh = rp * (SIN(vp + pp - op) * SIN(ip))
'
' convert to geocentric rectangular coordinates
'
xg = xh - xe
yg = yh - ye
zg = zh
'
' rotate around x axis from ecliptic to equatorial coords
'
ecl = 23.429292# * rads# 'value for J2000.0 frame
xeq = xg
yeq = yg * COS(ecl) - zg * SIN(ecl)
zeq = yg * SIN(ecl) + zg * COS(ecl)
'
' find the RA and DEC from the rectangular equatorial coords
'
ra = FNatn2(yeq, xeq)
dec = ATN(zeq / SQR(xeq * xeq + yeq * yeq))
rvec = SQR(xeq * xeq + yeq * yeq + zeq * zeq)
PRINT USING pr$; " RA : "; FNdegmin(ra * degs / 15)
PRINT USING pr$; " DEC : "; FNdegmin(dec * degs)
PRINT USING pr2$; "Distance : "; rvec
END
I used Manfred Ding's excellent and convenient Ephemeris Tool 2.2 (not to be confused with the asteroid and comet package of similar name) to check the accuracy of the mean element positions. Ephtool can be used with an add-in DLL which calculates planetary positions to the full accuracy of the VSOP-87 series. Ephtool is programmed in Delphi, and uses a rather nice 'spreadsheet control' to present data. The spreadsheet knows a useful range of formulas (including such esoterica as '=countif()'), including basic statistical formulas, and has 'autofill' and a few other Excel like features.
I implemented the mean element position calculation as a series of 'user defined functions' in an Excel VBA 'module' (using a method very similar to the QBASIC code above), and I constructed an Excel spreadsheet which calculated the RA and DEC using the mean element method. Accurate positions for the same time period were copied from Ephtool 2.2 (be sure to select the J2000 option from the 'ephemeris setup' menu). It takes about a minute on an ageing Gateway 75 Mhz Pentium to calculate 7300 planet positions using this VBA module.
The error was defined as;
Ephtool position - mean element position
in each case. I calculated errors for roughly 18000 days either side of J2000 (about 50 years) and positions were calculated at 10 day intervals. Errors are stated in units of arcseconds both for RA and for DEC. I also calculated the standard deviation of the errors - this gives a 'root mean square' type average error. The results for Mars, Jupiter and Neptune are shown below;
Mars Jupiter Neptune
ra dec ra dec ra dec
max + 223 58 412 185 55 35
max - -358 -129 -683 -226 -58 -38
stdev 50 30 203 69 27 20
Errors in arcseconds both for RA and DEC
As you can see from these figures, Mars has an assymetrical error curve, the negative errors (ie when the mean element position is ahead of the VSOP position) are larger than the positive errors. The ratio of the maximum error to the 'root mean square' error is around 7:1 for RA, from which we might expect a 'corrugated' error graph.
Jupiter also has an assymetrical graph, but the 'peak to mean' ratio is lower at about 7:2. Neptune has more symmetrical figures and a less corrugated curve, until we remember that the whole 100 year period here is less than 2/3 of an orbit!
Screenshots of the error graphs are available;
I may be deluding myself, but the appearence of these graphs tends to suggest a large amplitude term which depends on the planet period, and a smaller amplitude term which depends on the Earth's period. In the graphs for Jupiter, the time axis has 'major ticks' of 4333 days and minor ticks of 365 days. The Neptune time axis has 10 and 1 year periods marked.
Another way of summarising these results is to calculate the fraction of the results out of the 3600 or so which lie within a certain range of the VSOP position. The tables below show the fraction of results within the specified range of the VSOP position;
Mars
ra dec dist on sky
± 5" 0.082 0.092 0.005
± 10" 0.165 0.176 0.028
± 20" 0.318 0.371 0.099
± 30" 0.473 0.630 0.246
± 40" 0.586 0.862 0.444
± 60" 0.774 0.974 0.711
±120" 0.940 0.999 0.931
±300" 0.999 1.000 0.998
±360" 1.000 1.000 1.000
dist on sky = sqrt(ra*ra + dec*dec)
(measure of total angular difference
between mean element and VSOP position)
Jupiter
ra dec dist on sky
± 20" 0.080 0.161 0.014
± 60" 0.223 0.639 0.174
±120" 0.442 0.918 0.410
±240" 0.762 1.000 0.746
±300" 0.855 1.000 0.839
±600" 0.996 1.000 0.994
±720" 1.000 1.000 1.000
Neptune
ra dec dist on sky
± 5" 0.114 0.071 0.006
± 10" 0.230 0.148 0.034
± 20" 0.445 0.405 0.178
± 30" 0.692 0.843 0.402
± 40" 0.840 1.000 0.711
± 60" 1.000 1.000 0.978
±120" 1.000 1.000 1.000
The 'distance on the sky' is simply a measure of the angular distance between the mean element position and the VSOP position. I believe that the mean element position is good enough to plot positions on finder charts for most telescopes used by amateurs! Rise and set times should not be too much in error.
A conservative summary of these results would be to say that the positions of Mars using the mean elements are correct to within 6 arc minutes, of Jupiter to within 12 arc minutes, and of Neptune to within two arcminutes on the sky for 50 years either side of J2000.
A more meaningful (or optimistic) way of looking at the results might be to say that the position for Mars will be within 1 arc minute 'most of the time' (i.e. 71.1% of results), and Jupiter within 4 arc minutes (74.6%), and Neptune within 40 arc seconds (71.1%) by a similar criteria. The Explanatory Supplement states a 'approximate maximum error' for Jupiter of 300" in RA and 100" in DEC, which is consistent with a criteria of 'right 85% of the time'.
If you compare these results with the corresponding results for Mars using the osculating orbits, you will see that the error does not worsen with time from J2000 to the same extent as with the osculating elements. On the other hand, using the mean elements gives worse errors than using the osculating elements for dates within few years either side of the date of the elements. The moral is clear: use mean elements for 'reasonable' accuracy over a spread of time - use osculating elements for higher accuracy over short times either side of the date of the elements. It may even be worth allowing for the light travel time (sometimes called 'planetary abberation') when calculating results from the osculating elements.
I expect the error using these mean elements to be more apparent in the case of Jupiter, Saturn and the outer planets as a result of perturbations of the mean orbit. Unlike the osculating elements, you can correct the positions found using the mean elements for perturbations - see the extensive page provided by Paul Schlyter for an example of perturbation corrections for Jupiter, Saturn and Neptune. Paul's page is at;
http://hotel04.ausys.se/pausch/comp/ppcomp.html
[Root ]
Error in Neptune elements found and corrected 19980829
Last Modified 1st
September 1998
Keith Burnett