Planetary positions for 300 years - Formulae and method by A. Ahad

Astronomical Observations & Research

			Formulae and method for calculating planetary positions from mean orbital elements. [A simplified, lower precision algorithm by A. Ahad] Article compiled: c. 1984-85, posted: 26 September, 2004 IMPORTANT: This page requires you to have Math fonts on your computer Copyright © 2004 Abdul Ahad. All rights reserved. This article provides a sequence of simply presented analytical steps, which when followed through, will enable one to calculate the apparent position of a planet at a given instant as seen from Earth, suitable for basic observational purposes. The formulae and method presented here can be evaluated either using a pocket calculator (with scientific functions like sines, cosines, tangents, etc) or you can program the steps into a home computer quite easily. The latter method is favoured, as it enables you to run a continuous ephemeris over a period of time. Obtaining the precise geocentric position of a planet or other solar system object at a given instant using analytical methods is not an exact science. This is because the so-called "many body problem" of classical celestial mechanics - of which our solar system is a prime example - does not have a 'closed' mathematical solution. The orbital elements of each planet are typically subjected to long periodic variations arising from the gravitational perturbations by each of the other planets. Hence, all ephemerides are calculated to a given order of approximation. The following is << my own >> algorithm which I painstakingly put together as long ago as between 1984 and 1985, so it is of historical importance to me! It is based on an order of analytical approximation which I have personal experience of, having used it quite extensively in the pre-home computer era through the mid- to late 1980s. Having compared the results from these equations and orbital elements against those published in official ephemeris such as the Astronomical Almanac over an albeit limited period of time, I would "project" errors in geocentric positions for the principal planets to not exceed 0.25 degrees in the 300 year time frame 1800-2100 A.D. The algorithm which follows may be used to compute geocentric positions of the Sun and principal planets for a given date and time for basic observational purposes. Step 1: Geocentric rectangular co-ordinates of the Sun from mean elements I shall refer to the instant for which the planetary position is to be calculated as the "ephemeris date". The first step is to translate the chosen ephemeris date into "Julian date". The Julian date is defined to be a continuous count of the number of days elapsed at Greenwich since noon on 1st January 4,713 B.C. and serves as a convenient basis for measuring the time interval between any two epochs. Here is an online source that lets you work out the Julian date for any given calendar ephemeris date. Next, from the mean elements page, evaluate the Sun's mean elements for the ephemeris date. The Sun's apparent yearly motion projected onto the geocentric celestial sphere is illustrated in Fig. 1 below:- Fig. 1 - Geocentric coordinates of the Sun [AA Institute] In Fig. 1 above consider a non-rotating 3D reference frame consisting of three orthogonal (x,y,z) axes passing through the centre of the Earth as origin. Further, imagine the positive direction of the x-axis to point towards the "First Point of Aries" (vernal equinox), the positive y-axis to be directed towards a point 90 degrees west on the celestial equator and the positive z-axis directed towards the north pole of the Earth. The xy-plane (shaded green) is the plane of the celestial equator. The instantaneous position of the Sun can then be specified by 3D equatorial rectangular co-ordinates (x,y,z) at any point on its yearly apparent travels along the ecliptic, as shown in the above diagram. We evaluate the Sun's geocentric equatorial rectangular coordinates (x,y,z) from its mean elements for the ephemeris date by working through the below sequence:- Ecliptic longitude (degrees): k = L+(360e/o) degrees sin g+0.0200sin 2g+0.0003sin 3g+... Geocentric distance (AUs): R = 1.000140 - ecos g - 0.000140cos 2g - 0.000 002cos 3g -... Equatorial rectangular coordinates (AUs): x = Rcos k y = Rsin k cos e z = Rsin k sin e The Sun's geocentric R.A. and Declination: a = arc tan (cos etan k) d = arc sin (sin esin k) Apparent semi-diameter = 961.18/R arc seconds Step 2: Heliocentric rectangular coordinates of the planet from mean elements Calculate mean orbital elements of the planet at the ephemeris date from the mean elements page. The planet's orbital path projected onto the heliocentric celestial sphere is shown in Fig. 2 below :- Fig. 2 - Heliocentric coordinates of the planet [AA Institute] In a similar fashion to the geocentric celestial reference frame described above, the planet's instantaneous heliocentric orbital position can be specified by a set of orthogonal (x0,y0,z0) co-ordinates as illustrated in Fig. 2 above. We obtain the planet's heliocentric equatorial rectangular coordinates (x0,y0,z0) on the ephemeris date by working through the below sequence. Calculate the Mean Anomaly (M) expressed in radian measure from:- M = (L- p)/57.2957796 Next, we have to solve Kepler's equation to get the Eccentric Anomaly (E) from the Mean Anomaly (M). Here is an iterative process for solving Kepler's equation which I have tested to the order of up to seven approximations, by taking the Mean Anomaly (M) as the first approximation (E1=M). The rate of convergence between successive approximations I have found to be inversely proportional to the orbital eccentricity (e), but even in the case of Mercury, which has the highest eccentricity of the principal planets, the 7-step iterations give results which are within limits to maintain a fair level of accuracy in the final geocentric position calculations:- E1=M+esin M E2=M+esin E1 E3=M+esin E2 E4=M+esin E3 E5=M+esin E4 E6=M+esin E5 E7=M+esin E6 Finally, we take the seventh order of approximation (E7) from above and work out Eccentric Anomaly (E) in degrees:- E = (M+esin E7)57.2957796 True Anomaly: m = 2arc tan (square root [(1+e)/(1-e)]tan (E/2)) degrees Ecliptic longitude: L = p + m degrees Ecliptic latitude: B = isin (L- W) degrees [where i = orbital inclination] Radius vector (AUs): r = a(1-ecos E) [where a = semi-major axis of orbit] Heliocentric equatorial rectangular coordinates (AUs): x0 = rcos Lcos B y0 = r[cos esin Lcos B - sine esin B] z0 = r[sin esin Lcos B + cos esin B] Step 3: Geocentric rectangular coordinates and R.A. and Declination of the planet The geocentric equatorial rectangular coordinates are derived by simply adding the planet's heliocentric coordinates to the Sun's geocentric coordinates, since both sets of coordinates are based on the same celestial frame of reference (the two orthogonal (x,y,z) and (x0,y0,z0) sets of axes merge into one and point in the same directions at infinity). Rectangular coordinates: n = x + x0 g = y + y0 f = z + z0 Right Ascension: a = arc tan (g/n) + K degrees ...... (1) In (1) above, the value of K is dependent upon the signs of the numerator and denominator in the fraction within the "arc tan" function: K = 0 if the numerator is positive and the denominator is positive K = 180 if numerator is positive, denominator is negative K = 180 if numerator is negative, denominator is negative K = 360 if numerator is negative, denominator is positive The value of K decided using the above criteria, when added as in (1) above, will fix the quadrant within the interval [0 to 360 degrees] in which the Right Ascension lies. Declination: d = arc sin (f /D) degrees Where D = geocentric distance = [square root (n^2 + g^2 + f^2)] (AUs) Additional formulae are available for calculating observational parameters such as apparent disc diameter, fraction illuminated (phase), elongation from the Sun, approximate visual magnitude, etc. which I have not detailed here. Illustration of Accuracy I plan to observe Mars at midnight on the 1st January, 2006. What are its geocentric coordinates predicted to be at that instant? Which constellation will Mars be in on that date? The ephemeris date is JD 2453736.5 and the coordinates derived from my given algorithm are found to be:- R.A. = 02 hrs 32.6 mins Dec. = +16º 37' Distance = 0.776 AU (116 million kilometres) App. Diam = 12.1" By way of comparison, Chris Marriott's "Sky Map Pro" gives: R.A. = 02 hrs 32.6 mins Dec. = +16º 37' Distance = 0.775 AU (116 million kilometres) App. Diam = 12.1" Hence, the positional errors from my simplified algorithm in this case are negligible, and just 0.001 AU out in the calculated geocentric distance. By referring to a star chart such as "Norton's", the predicted co-ordinates would place Mars in Aries (close to the sixth magnitude star 27 Arietis) on this night. Example of Practical Applications: Transit of Venus, June 8 2004 Using my then newly compiled algorithm along with a graphical technique of visualisation and linear interpolation << back in May 1985 >> I was able to show that the 2004 transit of Venus across the face of the sun would have the following basic circumstances:- Time of geocentric conjunction in Right Ascension between the center of the Sun's disk and Venus: June 8th at approx. 0900 E.T. (JD 245 3165.875) First contact (beginning of transit) would occur at approx. 0550 U.T. at P.A. 119 degrees east. Mid transit would occur at approx. 0830 hours U.T. Last contact (end of transit) would occur at approx. 1110 U.T. at P.A. 143 degrees west. Total duration of the transit would be about 5.3 hours At mid transit both Venus and the Sun would be in the zenith seen from latitude 22.7 degrees north, longitude 52.3 degrees east of Greenwich on the Earth's surface (in the United Arab Emirates). Summary Based on the above methodology and given orbital elements, it is easily possible to obtain the geocentric position of a planet, comet or asteroid at any given instant. However, the reverse problem facing 18th/19th century pioneers in celestial mechanics - deducing the orbital elements from 3 or more observed positions of a planet or comet - is one of far greater analytical complexity and is outside the scope of my article here. A vast amount of literature is available on the subject, and in particular the method of Gauss is favoured by many current day orbit computing professionals. R E F E R E N C E S I originally formulated the above methodology and calculations in their entirety back in 1984-85 << when I was just 15 >> by extracting and meticulously piecing together data and formulae from a number of widely scattered ephemeris sources. At that time, I also sought guidance from a book titled "Textbook on Spherical Astronomy" by W.M. Smart [first published in Leipzig, Germany, circa 1930s] which helped me in fully understanding the geocentric and heliocentric reference frames and their practical relationships with the formulae of planetary motion. I remember the book covering in detail the above mentioned iterative process for solving Kepler's equation. I recently found the below accompanying explanatory text which I wrote down in one of my old school exercise books from which I have taken the above formulae, methodology & orbital elements for this article:- Sources of orbital elements "All the orbital elements were obtained from expressions given in the 'Explanatory supplement to the Astronomical Ephemeris, American Ephemeris and Nautical Almamanc'. The geocentric orbital elements of the Sun are exactly the same as those which were used to evaluate 'Newcomb's Tables of the Sun' - except for the rate of change of Mean Anomaly, g, which I have modified to ensure it remains in line with [then] current values. The elements of Mercury, Venus and Mars are the same as those which were used in 'Newcomb's Theories of the Inner Planets', with Ross's corrections for Mars. For Jupiter and Saturn the values were obtained from 'Hill's Tables', by applying the variations of Le Verrier and Galliot and backdating to epoch 1900. In the case of Uranus and Neptune, the elements were taken from 'Newcomb's Tables of the Planets' and backdated to 1900." - A. Ahad [circa. 1985] NOTE: The original mathematicians behind much of this work perished many many years ago and my memory over the decades has faded considerably... but I will do my level best to provide answers in relation to this work if anyone has any questions. The Astronomical Almanac online Celestial Mechanics at Wolfram Research Celestial Mechanics: A Computational Guide for the Practitioner Methods of orbit determination for the microcomputer Orbit modeling online Copyright © 2004 Abdul Ahad. All rights reserved.

Hosted by www.Geocities.ws

R E F E R E N C E S