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