Accuracy of planet positions using osculating elements

[Root ]

Contents

Overview

[Top]

"I was wondering why and how a sickle of just this thickness (0.00429) came into being. While this thought was driving me around.... I stumbled entirely by chance on the secant of the angle 5 degrees and 18 minutes, which is the measure of the greatest optical equation. When I realised that this secant equals 1.00429, I felt as if I had been awakened from a sleep..."
Kepler

This page is concerned with the accuracy of the Equatorial coordinates (RA and DEC) calculated using a method based on elliptical orbits using osculating elements from the Astronomical Almanac for 1997.

I have compared the positions calculated using a simple QBASIC program with those produced by the Planeph 4.2 program written by Francou and Chapront. For Mars, positions produced using the Osculating elements were found to be accurate to better than 4 seconds of time in RA and 20 arcseconds in DEC within 1 year of the date of the elements. The largest errors increase to 2 minutes of RA and 14 minutes of DEC for 10 years either side of the date of the elements. The error does not increase with time in a simple way - there are 'peaks' in the graph of error against time which recur at certain intervals.

Osculating elements and elliptical orbits

If you use an elliptical orbit to work out the position of a planet, you are assuming that the gravitational attraction between the planets themselves is very small in comparison to the attraction between the Sun and the planet. This is a good approximation! There are small, long term effects on the positions of the planets as a result of their gravitational interactions (known as perturbations of the orbit). The osculating elements of a planet's orbit take account of the perturbations for a given date, and the positions calculated from the osculating elements will be accurate around that date.

The elements used here are valid for the Julian date 2450680.5 (0h on 20th August 1997) and are taken from page E3 of the Astronomical Almanac for 1997. You can get the numerical values of the elements from the QBASIC listing below. Definitions of each of the elements and diagrams explaining their meaning are available from:

http://users.commkey.net/Braeunig/space/

or from Duffett-Smith section 53, or any celestial mechanics textbook.

Precession and the equinox and mean ecliptic of the elements

If you want to use the method here to find rough positions of the planets good to a few minutes either way, then this section is not essential reading. If you need to use the osculating elements for a more demanding application, then the information may be useful.

The osculating elements are referred to a given equinox and mean ecliptic. I have chosen to use the elements referred to the equinox and mean ecliptic of J2000.0 here (page E3 of the Astronomical Almanac). The values of RA and DEC calculated from these elements will also be referred to the J2000.0 epoch. This choice has several advantages and disadvantages;

You can choose to use the osculating elements referred to the equinox and mean ecliptic of the date of the elements. These elements are tabulated at intervals of 40 days for the year on pages E4 and E5 of the Astronomical Almanac.

Remember that the positions found using these elements will be referred to the date of the elements, not the date you want the positions for! If the date you want positions for is 4 years earlier than the date of the elements, then precession will amount to an error of 3 arc minutes in RA, and a 1 arc minute in DEC. See Duffett-Smith [section 34], or Meeus [Chapter 20] for low precision formulas for precession.

The 'mean ecliptic' does not take account of nutation, the 'nodding' of the Earth as the axis of rotation wobbles (precession). The maximum size of the nutation correction is about 18 arcseconds, and you do not want to correct for nutation if you are plotting positions on star charts!

See Meeus [Chapter 21] for a nutation theory based on J2000.0 with full accuracy. Duffett-Smith [section 35] uses an older theory to an accuracy of about half an arcsecond.

Procedure

[Top]

I wrote a simple QBASIC program based on the treatment of elliptical orbits found in Astronomy with your Calculator by Peter Duffett-Smith. The program includes an iterative solution to Kepler's equation, but does not include any correction for perturbations the orbit by other planets (other than those implied by the use of the osculating elements), and no correction for light travel time is made. As I use the osculating elements referred to the equinox and mean ecliptic of J2000.0, I can use a fixed value for the obliquity of the ecliptic when calculating the RA and DEC.

I modified the simple QBASIC program to produce positions for Mars from 4000 days before the date of the elements to 4000 days after the date of the elements, in steps of 40 days. This data was written to a file, and loaded into a spreadsheet.

I then used the Planeph 4.2 program to generate a series of values of RA and DEC over the same range of dates. The data can also be saved to a file, and loaded into a spreadsheet. The output has column headings added after every 50 lines, and these must be removed before the errors can be calculated!

The Planeph 4.2 program runs under DOS, and provides positions accurate to better than 0.01 arcsecond between the years 1900 and 2100. The program may be obtained from;

ftp://cdsarc.u-strasbg.fr/pub/cats/VI/87/

The positions found from Planeph 4.2 were geocentric astrometric positions (thus including the light time correction) referred to the fixed Equatorial frame for J2000.0, using the UT time scale. Planeph 4.2 includes nutation corrections for 'true' frames, and so I assume it does not correct for nutation in 'fixed' frames.

The spreadsheet was then used to calulate the errors in seconds of time for RA and in arcseconds for DEC. I always use the formula;

Error = Osculating element estimate - planeph 4.2 position
to define 'error'.

Accuracy

[Top]

The 'root mean square' errors for various ranges of position dates before and after the dates of the elements look quite impressive;

RMS errors for Mars, averaged over different ranges of time before
and after the date of the element.

                              RA Sec   DEC arcsec
+-1 year (2 year period)         2          8
+-3 years (6 year period)        5         24
+-10 years (20 year period)     26        145
The maximum errors within each time range either side of the date of the elements tell a different story;
Mars - Largest absolute error

Within;                        RA sec     DEC arcsec
+-1 year (2 year period)          4            17
+-3 years (6 year period)        15            80
+-10 years (20 year period)     130           832
With 'peak to mean ratios' like this, you might guess that a time series graph of error against position date would show a complex behaviour, not a simple monotonic rise with time from the date of the elements.

Error time series graphs for Mars

graph of RA error against year from date of
elements

The main features of the RA error time series for Mars are

graph of DEC error against year from date of
elements

The main features of the DEC error time series for Mars are;

The QBASIC program

[Top]
'*********************************************************
'   This program will calculate the positions of the
'   major planets using the current 'osculating elements'
'   from the Astronomical Almanac.
'
'   A simple elliptical orbit
'   is assumed for both the planet and the Earth - no
'   corrections are made from within the program as the
'   osculating elements will already take account of
'   perturbations.
'
'   The program is a simple 'straight through' version
'   of a manual calculation. Major results are printed
'   at each stage to allow comparison with manual
'   calculations. The method is based on Duffett-Smith's
'   book 'Practical Astronomy with your Calculator'
'   3rd edition Cambridge University Press, 1988
'   section 54.
'
'   QBASIC program by Keith Burnett (kburnett@geocity.com)
'
'
'   Work in double precision and define some constants
'
DEFDBL A-Z
pr$ = "\         \#####.####"
pi = 4 * ATN(1)
tpi = 2 * pi
twopi = tpi
degs = 180 / pi
rads = pi / 180
'
'   list of elements el()
'   List of the osculating elements of the 9 major
'   planets in the format used in the Astronomical
'   Ephemeris. Item el(64) in list is the Julian date
'   of the elements. Item el(65) is the epoch of the
'   mean ecliptic and equinox the elements are referred to.
'
'   If you want positions referred to
'   the mean equator and equinox of the date of the
'   osculating elements, then use the elements listed
'   on pages E4 and E5 of the AA. If you want the positions
'   referred to the mean equator and equinox of J2000
'   then use the elements found on page E3 of the AA.
'
'
DIM el(9 * 7 + 2)
'
'   below are the osculating elements for JD = 2450680.5
'   referred to mean ecliptic and equinox of J2000
'
'Mercury
el(1) = 7.00507# * rads
el(2) = 48.3339# * rads
el(3) = 77.45399999999999# * rads
el(4) = .3870978#
el(5) = 4.092353# * rads
el(6) = .2056324#
el(7) = 314.42369# * rads
'Venus
el(8) = 3.39472# * rads
el(9) = 76.6889# * rads
el(10) = 131.761# * rads
el(11) = .7233238#
el(12) = 1.602158# * rads
el(13) = .0067933#
el(14) = 236.94045# * rads
'Earth
el(15) = .00041# * rads
el(16) = 349.2# * rads
el(17) = 102.8517# * rads
el(18) = 1.00002#
el(19) = .9855796# * rads
el(20) = .0166967#
el(21) = 328.40353# * rads
'Mars
el(22) = 1.84992# * rads
el(23) = 49.5664# * rads
el(24) = 336.0882# * rads
el(25) = 1.5236365#
el(26) = .5240613# * rads
el(27) = .0934231#
el(28) = 262.42784# * rads
'Jupiter
el(29) = 1.30463# * rads
el(30) = 100.4713# * rads
el(31) = 15.6978# * rads
el(32) = 5.202597#
el(33) = 8.309618000000001D-02 * rads
el(34) = .0484646#
el(35) = 322.55983# * rads
'Saturn
el(36) = 2.48524# * rads
el(37) = 113.6358# * rads
el(38) = 88.863# * rads
el(39) = 9.571899999999999#
el(40) = .03328656# * rads
el(41) = .0531651#
el(42) = 20.95759# * rads
'Uranus
el(43) = .77343# * rads
el(44) = 74.0954# * rads
el(45) = 175.6807# * rads
el(46) = 19.30181#
el(47) = .01162295# * rads
el(48) = .0428959#
el(49) = 303.18967# * rads
'Neptune
el(50) = 1.7681# * rads
el(51) = 131.7925# * rads
el(52) = 7.206# * rads
el(53) = 30.26664#
el(54) = .005919282# * rads
el(55) = .0102981#
el(56) = 299.8641# * rads
'Pluto
el(57) = 17.12137# * rads
el(58) = 110.3833# * rads
el(59) = 224.8025# * rads
el(60) = 39.5804#
el(61) = .003958072# * rads
el(62) = .2501272#
el(63) = 235.7656# * rads
'Dates
el(64) = 2450680.5# 'date of elements
el(65) = 2451545#   'date of mean ecliptic and equinox of elements
'
'   Get the days to J2000
'   h is UT in decimal hours
'   FNday only works between 1901 to 2099 - see Meeus chapter 7
'
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
'
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 
'   use 8 for most situations, 12 for comets!
'   See Duffett-Smith section 47 or Meeus
'   Chapters 29 and 32 
'
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
'
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)
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 = 7 * (p - 1)
ip = el(q + 1)
op = el(q + 2)
pp = el(q + 3)
ap = el(q + 4)
np = el(q + 5)
ep = el(q + 6)
lp = el(q + 7)
ie = el(15)
oe = el(16)
pe = el(17)
ae = el(18)
ne = el(19)
ee = el(20)
le = el(21)
eldate = el(64) - 2451545#
'
'   now find position of Earth in orbit
'
me = FNrange(ne * (d - eldate) + le - pe)
PRINT USING pr$; " Earth M : "; me * degs
ve = FNkep(me, ee, 12)
PRINT USING pr$; " Earth V : "; ve * degs
le2 = FNrange(ve + pe)
PRINT USING pr$; " Earth L : "; le2 * degs
re = ae * (1 - ee * ee) / (1 + ee * COS(ve))
PRINT USING pr$; " Earth R : "; re
'
'   and position of planet in its orbit
'
mp = FNrange(np * (d - eldate) + lp - pp)
PRINT USING pr$; "Planet M : "; mp * degs
vp = FNkep(mp, ep, 12)
PRINT USING pr$; "Planet V : "; vp * degs
lp2 = FNrange(vp + pp)
PRINT USING pr$; "Planet L : "; lp2 * degs
rp = ap * (1 - ep * ep) / (1 + ep * COS(vp))
PRINT USING pr$; "Planet R : "; rp
'
'   project planet orbit onto ecliptic to get heliocentric
'   longitude, latitude and radius vector
'
phi = FNasin(SIN(lp2 - op) * SIN(ip))
PRINT USING pr$; "     phi : "; phi * degs
lp3 = FNatn2(SIN(lp2 - op) * COS(ip), COS(lp2 - op)) + op
PRINT USING pr$; " helio L : "; lp3 * degs
rp2 = rp * COS(phi)
PRINT USING pr$; " Helio R : "; rp2
'
'   calculate the geocentric ecliptic longitude
'
IF p > 2 THEN
    '
    'outer planet
    '
    lambda = FNrange(ATN(re * SIN(lp3 - le2) / (rp2 - re * COS(lp3 - le2)))
    + lp3)
ELSE
    '
    'inner planet
    '
    A = ATN(rp2 * SIN(le2 - lp3) / (re - rp2 * COS(le2 - lp3)))
    lambda = FNrange(pi + le2 + A)
END IF
PRINT USING pr$; "  lambda : "; lambda * degs
'
'   calculate the geocentric ecliptic latitiude
'
beta = ATN(rp2 * TAN(phi) * SIN(lambda - lp3) / (re * SIN(lp3 - le2)))
PRINT USING pr$; "    beta : "; beta * degs
'
'   find mean obliquity of ecliptic - just 23.439292 if J2000.0 elements
'   used. el(65) contains the equinox and mean ecliptic which elements
'   are referred to.
'
t = (el(65) - 2451545#) / 36525
e = 23.439292# + (-46.815 * t - .0006 * t ^ 2 + .00181 * t ^ 3) / 3600
e = rads * e
'
'   find RA and DEC
'
alpha = FNatn2(SIN(lambda) * COS(e) - TAN(beta) * SIN(e), COS(lambda))
alpha = FNrange(alpha)
PRINT USING pr$; "   alpha : "; alpha * degs / 15
delta = FNasin(SIN(beta) * COS(e) + COS(beta) * SIN(e) * SIN(lambda))
PRINT USING pr$; "   delta : "; delta * degs
'
'   print the results
'
END
'*********************************************************

Sample output of QBASIC program

Below is the output produced by the simple program for Mars (planet 4), on 1997 June 15 at 14:47 UT.
  year  : 1997
  month : 6
  day   : 15
hour UT : 14
 minute : 47
 planet : 4
    days :  -929.8840
 Earth M :   161.1107
 Earth V :   161.7181
 Earth L :   264.5698
 Earth R :     1.0158
Planet M :   252.0744
Planet V :   242.2900
Planet L :   218.3782
Planet R :     1.5789
     phi :     0.3589
 helio L :   218.3839
 Helio R :     1.5789
  lambda :   178.4491
    beta :     0.4962
   alpha :    11.9183
   delta :     1.0721
This compares with the following shortened output from the Planeph 4.2 program for the same date;
PLANETARY EPHEMERIDES                                PLANEPH 4.1
----------------------------------------------------------------

MARS             Astrometric Geocentric Positions
                 Fixed Equatorial Frame J2000.00 JD2451545.00000
                 Universal Time (UT)


Date         Hour       Right Ascension    Declination    

1997 Mar 27  14h47m00s    11.65381466 h   + 5.9785870 
1997 May 06  14h47m00s    11.26372419 h   + 6.6698851 
1997 Jun 15  14h47m00s    11.91811319 h   + 1.0738741 
1997 Jul 25  14h47m00s    13.12219611 h   - 7.4870472  
1997 Sep 03  14h47m00s    14.67884943 h   -16.5011926 

----------------------------------------------------------------
Here the errors are half a second of time on RA and 6 arc seconds on DEC, which is not typical (see Accuracy section above).

The spreadsheets

[Top]

I have set up an Excel 95 'workbook' with a page for each planet's orbit. The calculations are set out step by step, and the answers are available in decimal degrees for each major result. The workbook is available as a .ZIP file from:

http://www.btinternet.com/~kburnett/kepler/oscxls.zip (9 Kb)

If you do not have Microsoft Excel version 5 or 95 or later, then I have provided a set of 8 spreadsheets in Lotus .WKS format. Most spreadsheets can load files of this kind. These spreadsheets are available in a .ZIP file from:

http://www.btinternet.com/~kburnett/kepler/oscwks.zip (16 Kb)

Please note that the ATN2 function found on most spreadsheets works in a different way to that defined in C and added as a user function to other programming languages. The ATN2 function on Excel (and in Lotus 1-2-3, which everyone else seems to have copied) returns an angle in radians in the range -pi to +pi. Below is a copy of the 'help file' entry for atn2 in Excel 95;

Returns the arctangent of the specified x- and y- coordinates.
The arctangent is the angle from the x-axis to a line containing
the origin (0, 0)) and a point with coordinates (x_num, y_num).
The angle is given in radians between -p and p, excluding -p.

Syntax

ATAN2(x_num, y_num)

X_num    is the x-coordinate of the point.
Y_num    is the y-coordinate of the point.

Remarks

A positive result represents a counterclockwise angle from the
x-axis; a negative result represents a clockwise angle.ATAN2(a,b)
equals ATAN(b/a), except that a can equal 0 in ATAN2.
        
If both x_num and y_num are 0, ATAN2 returns the #DIV/0! error
value.
        
To express the arctangent in degrees, multiply the result by
180/PI( ).

Examples

ATAN2(1, 1) equals 0.785398 (p/4 radians)
ATAN2(-1, -1) equals -2.35619 (-3p/4 radians)
ATAN2(-1, -1)*180/PI() equals -135 (degrees)
I'm still scribbling on the backs of envelopes trying to work out how to use this ATN2 function to give the right answers for these spreadsheets - the values may be in the wrong quadrant at present!

References

[Top]

Below is a list of the books which I have used in compiling this page. These books should be found in most University libraries in English speaking countries.

Duffett-Smith, Peter
Practical Astronomy with your calculator
Cambridge University Press
3rd edition 1988
ISBN 0-521-35699-7

Meeus, Jean
Astronomical Algorithms
Willmann-Bell
1st English edition, 1991
ISBN 0-943396-35-2

[Root ]


Last Modified 8th July 1997
Keith Burnett
kburnett@btinternet.com