This paper has been modified from the original form I received it in. I have
removed the double spacing for easier reading.
From pausch@electra.saaf.se Sun Jul 07 04:26:19 1996
Newsgroups: sci.astro.amateur
Subject: Re: REQ:Algorithm(s) for plotting the position of the planets at a given time.
From: pausch@electra.saaf.se (Paul Schlyter)
Date: 7 Jul 1996 11:26:19 +0200
In article <91Yp87AS8D3xEwwV@chrism.demon.co.uk>,
Chris Marriott wrote:
> In article <31db43b9.22110924@news.alt.net>, Greg Miller
> writes
>> Can anyone give me information regarding algorithms which give any
>>indication of the position of the planets in our solar system
>
> Computing accurate planetary positions is non-trivial, requiring the
> computation of HUNDREDS of trig terms. Meeus' "Astronomical Algorithms"
> published by Willman-Bell is the "bible" for all astronomical
> calculations of this type.
Why suggest an "overkill" method like that?
If you aim for the more modest accuracy of 1-2 arc minutes for any
date within a century or two from the present, about three dozen
perturbation terms will be enough: 19 for the moon, 7 each for
Jupiter and Saturn, and 2 for Uranus. None are needed for the inner
planets or Neptune where a two-body computation will be enough. The
text below describes the process.
Since this question is re-asked from time to time, perhaps this
should become a "FAQ"?
----------------------------------------------------------------------
HOW TO COMPUTE PLANETARY POSITIONS
==================================
By Paul Schlyter, Stockholm, Sweden
pausch@saaf.se
Below is a description of how to compute the positions of the Sun and
Moon and the major planets, as well as for comets and minor planets
from a set of orbital elements.
The algorithms have been simplified as much as possible while still
keeping a fairly good accuracy. The accuracy of the computed
positions is a fraction of an arc minute for the Sun and the inner
planets, about one arc minute for the outer planets, and 1-2 arc
minutes for the Moon. If we limit our accuracy demands to this
level, one can simplify further by e.g. ignoring the difference
between mean, true and apparent positions.
The positions computed below are for the 'equinox of the day', which
is suitable for computing rise/set times, but not for plotting the
position on a star map drawn for a fixed epoch. In the latter case,
correction for precession must be applied, which is most simply
performed as a rotation along the ecliptic.
These algortihms were developed by myself back in 1979, based on a
preprint from T. van Flandern's and K. Pulkkinen's paper "Low
precision formulae for planetary positions". It's basically a
simplification of these algorithms, while keeping a reasonable
accuracy. They were first implemented on a HP-41C programmable
pocket calculator, in 1979. Nowadays considerable more accurate
algorithms are available of course, as wel as more powerful
computers. Nevertheless I've retained these algorithms as what I
believe is the simplest way to compute solar/lunar positions with an
accuracy of one or a few arc minutes, or better.
1. Introduction
The text below describes how to compute the positions in the sky of
the Sun, Moon and the major planets out to Neptune. THe algorithm
for Pluto is taken from a fourier fit to Pluto's position as computed
by numerical integration at JPL. Positions of other heavenly bodies
as well (i.e. comets and asteroids) can also be computed, if their
orbital elements are available.
These formlulae may seem complicated, but I believe this is the
simplest method to compute planetary positions with the fairly good
accuracy of about one arc minute (=1/60 degree). Any furhter
simplifications will yield lower accuracy, but of course that may be
ok, depending on the application.
2. A few words about accuracy
The accuracy requirements are modest: a final position with an error
of no more than 1-2 arc minutes is aimed at (one arc minute = 1/60
degree). This accuracy is in one respect quite optimal: it is the
highest accuracy one can strive for, while still being available to
do many simplifications. The simplifications that is done here are:
1: Nutation and aberration are both ignored.
2: Planetary aberration (i.e. light travel time) is ignored.
3: The difference between Ephemeris Time (ET), (nowadays called
Terrestial Dynamic Time, or TDT) and Universal Time (UT) is ignored.
4: Precession is computed in a simplified way, by doing a simple
addition or subtraction to the ecliptic longitude.
5: Higher-order terms in the planetary orbital elements are ignored.
This will give an additional error of up to 2 arc min in 1000 years
from now. For the Moon, the error will be larger: 7 arc min 1000
years from now. This error will grow as the square of the time.
6: Most planetary perturbations are ignored. Only the major perturbation
terms for the Moon, Jupiter, Saturn, and Uranus, are included. If
less accuracy is acceptable, these perturbations can be ignored as
well.
7: The largest Uranus-Neptune perturbation is accounted for in the
orbital elements of these planets. Therefore, the orbital elements
of Uranus and Neptune are less accurace, especially in the distant
past and future. The elements for these planets should therefore
only be used for at most a few centuries into the past and the future.
3. The time scale
The time scale in these formulae are counted in days. Hours, minutes,
seconds are expressed as fractions of a day. Day 0.0 occurs at
2000 Jan 0.0 UT (or 1999 Dec 31, 0:00 UT). This "day number" d is
computed as follows (y=year, m=month, D=date, UT=UT in hours+decimals):
d = 367*y - 7 * ( y + (m+9)/12 ) / 4 + 275*m/9 + D - 730530
Note that ALL divisions here should be INTEGER divisions. In Pascal,
use "div" instead of "/", in MS-Basic, use "\" instead of "/". In
Fortran and C, "/" can be used if y and m are declared as integers.
Finally, include the time of the day, by adding:
d = d + UT/24.0 (this is a floating-point division)
4. The orbital elements
The usual orbital elements are here denoted as:
N = longitude of the ascending node
i = inclination
w = argument of perihelion
a = semi-major axis
e = eccentricity
M = mean anomaly
Another quantity that will be used as well is:
ecl = obliquity of the ecliptic
First, compute the "d" of the moment of interest (see section 3). Then,
compute the obliquity of the ecliptic:
ecl = 23.4393 - 3.563E-7 * d
Then compute the orbital elements of the planet of interest. If you
want the position of the Sun or the Moon, you only need to compute the
orbital elements for the Sun or the Moon. If you want the position of
any other planet, you must compute the orbital elements for that planet
AND for the Sun. This is necessary to be able to compute the geocentric
position of the planet.
Please note that a, the semi-major axis, is given in Earth radii for the
Moon, but in "astronomical units" for the Sun and all the planets. One
astronomical unit is the mean distance between the Earth and the Sun.
When computing M (and, for the Moon, when computing N and w as well), one
will quite often get a result that is larger than 360 degrees, or is
negative (all angles are computed in degrees!!). If negative, add 360
degrees until positive. If larger than 360 degrees, subtract 360 degrees
until the value is less than 360 degrees. Note that, in most programming
languages, one must then multiply these angles with pi/180 to convert them
to radians, before taking the sine or cosing of them.
Orbital elements of the Sun:
N = 0.0
i = 0.0
w = 282.9404 + 4.70935E-5 * d
a = 1.000000
e = 0.016709 - 1.151E-9 * d
M = 356.0470 + 0.9856002585 * d
Orbital elements of the Moon:
N = 125.1228 - 0.0529538083 * d
i = 5.1454
w = 318.0634 + 0.1643573223 * d
a = 60.2666 (Earth radii)
e = 0.054900
M = 115.3654 + 13.0649929509 * d
Orbital elements of Mercury:
N = 48.3313 + 3.24587E-5 * d
i = 7.0047 + 5.00E-8 * d
w = 29.1241 + 1.01444E-5 * d
a = 0.387098
e = 0.205635 + 5.59E-10 * d
M = 168.6562 + 4.0923344368 * d
Orbital elements of Venus:
N = 76.6799 + 2.46590E-5 * d
i = 3.3946 + 2.75E-8 * d
w = 54.8910 + 1.38374E-5 * d
a = 0.723330
e = 0.006773 - 1.302E-9 * d
M = 48.0052 + 1.6021302244 * d
Orbital elements of Mars:
N = 49.5574 + 2.11081E-5 * d
i = 1.8497 - 1.78E-8 * d
w = 286.5016 + 2.92961E-5 * d
a = 1.523688
e = 0.093405 + 2.516E-9 * d
M = 18.6021 + 0.5240207766 * d
Orbital elements of Jupiter:
N = 100.4542 + 2.76854E-5 * d
i = 1.3030 - 1.557E-7 * d
w = 273.8777 + 1.64505E-5 * d
a = 5.20256
e = 0.048498 + 4.469E-9 * d
M = 19.8950 + 0.0830853001 * d
Orbital elements of Saturn:
N = 113.6634 + 2.38980E-5 * d
i = 2.4886 - 1.081E-7 * d
w = 339.3939 + 2.97661E-5 * d
a = 9.55475
e = 0.055546 - 9.499E-9 * d
M = 316.9670 + 0.0334442282 * d
Orbital elements of Uranus:
N = 74.0005 + 1.3978E-5 * d
i = 0.7733 + 1.9E-8 * d
w = 96.6612 + 3.0565E-5 * d
a = 19.18171 - 1.55E-8 * d
e = 0.047318 + 7.45E-9 * d
M = 142.5905 + 0.011725806 * d
Orbital elements of Neptune:
N = 131.7806 + 3.0173E-5 * d
i = 1.7700 - 2.55E-7 * d
w = 272.8461 - 6.027E-6 * d
a = 30.05826 + 3.313E-8 * d
e = 0.008606 + 2.15E-9 * d
M = 260.2471 + 0.005995147 * d
Please note than the orbital elements of Uranus and Neptune as given here
are somewhat less accurate. They include a long period perturbation
between Uranus and Neptune. The period of the perturbation is about 4200
years. Therefore, these elements should not be expected to give results
within the stated accuracy for more than a few centuries in the past and
into the future.
5. The position of the Sun
The position of the Sun is computed just like the position of any other
planet, but since the Sun always is moving in the ecliptic and since
the eccentricity of the orbit is quite small, a few simplifications
can be made. Therefore, a separate presentation for the Sun is given.
Of course, we're here really computing the position of the Earth in its
orbit around the Sun, but since we're viewing the sky from an Earth-
centered perspective, we'll pretend that the Sun is in orbit around the
Earth instead.
First, compute the eccentric anomaly, E, from the mean anomaly, M, and
from the eccentricity, e, (E and M in degrees):
E = M + e*(180/pi) * sin(M) * ( 1.0 + e * cos(M) )
or (if E and M are expressed in radians):
E = M + e * sin(M) * ( 1.0 + e * cos(M) )
Then compute the Sun's distance, r, and its true anomaly, v, from:
xv = r * cos(v) = cos(E) - e
yv = r * sin(v) = sqrt(1.0 - e*e) * sin(E)
v = atan2( yv, xv )
r = sqrt( xv*xv + yv*yv )
atan2() is a function that converts an x,y coordinate pair to the
correct angle in all four quadrants. It is available as a library
function in C and Fortran. In other languages, one has to write one's
own atan2() function. It's not that difficult:
atan2( y, x ) = atan(y/x) if x > 0
atan2( y, x ) = atan(y/x) +- 180 degrees if x < 0
atan2( y, x ) = sign(y) * 90 degrees if x = 0
Now, compute the Sun's true longitude:
lonsun = v + w
Convert lonsun,r to xs,ys:
xs = r * cos(lonsun)
ys = r * sin(lonsun)
(zs is of course 0.0). xs,ys is the Sun's position in a coordinate
system in the plane of the ecliptic. To convert this to equatorial,
rectangular, geocentric coordinates, compute:
xe = xs
ye = ys * cos(ecl)
ze = ys * sin(ecl)
Finally, compute the Sun's Right Ascension (RA) and Declination (Dec):
RA = atan2( ye, xe )
Dec = atan2( ze, sqrt(xe*xe+ye*ye) )
6. The position of the Moon and of the planets
First, compute the eccentric anomaly, E, from M, the mean anomaly, and
e, the eccentricity. As a first approximation, do (E and M in degrees):
E = M + e*(180/pi) * sin(M) * ( 1.0 + e * cos(M) )
or, if E and M are in radians:
E = M + e * sin(M) * ( 1.0 + e * cos(M) )
If e, the eccentricity, is less than about 0.05-0.06, this approximation
is sufficiently accurate. If the eccentricity is larger, set E0=E and then
use this iteration formula (E and M in degrees):
E1 = E0 - ( E0 - e*(180/pi) * sin(E0) - M ) / ( 1 - e * cos(E0) )
or, if E and M are in radians:
E1 = E0 - ( E0 - e * sin(E0) - M ) / ( 1 - e * cos(E0) )
For each new iteration, replace E0 with E1. Iterate until E0 and E1 are
sufficiently close together (about 0.001 degrees. For comet orbits with
very high eccentricites, a difference of less than 1E-4 or 1E-5 degrees
should be required).
Now compute the planet's distance and true anomaly:
xv = r * cos(v) = a * ( cos(E) - e )
yv = r * sin(v) = a * ( sqrt(1.0 - e*e) * sin(E) )
v = atan2( yv, xv )
r = sqrt( xv*xv + yv*yv )
7. The position in space
Compute the planet's position in 3-dimensional space:
xh = r * ( cos(N) * cos(v+w) - sin(N) * sin(v+w) * cos(i) )
yh = r * ( sin(N) * cos(v+w) + cos(N) * sin(v+w) * cos(i) )
zh = r * ( sin(v+w) * sin(i) )
For the Moon, this is its geocentric (Earth-centered) position in
the ecliptic coordinate system. For the planets, this is the
heliocentric (Sun-centered) position, also in the ecliptic coordinate
system. If one wishes, one can compute the ecliptic longitude and
latitude (this must be done if one wishes to correct for perturbations,
or if one wants to precess the position to a standard epoch):
lonecl = atan2( yh, xh )
latecl = atan2( zh, sqrt(xh*xh+yh*yh) )
As a check one could also compute sqrt(xh*xh+yh*yh+zh*zh) -- this
should of course equal r (except for small round-off errors).
8. Precession
If one wishes to compute the planet's position for some standard epoch,
such as 1950.0 or 2000.0 (for instance to be able to plot the position
on a star atlas), one must add the correction below to lonecl. If a
planet's and not the Moon's position is computed, one must also add
the same correction to lonsun, the Sun's longitude. The desired Epoch
is expressed as a year with possibly a fraction.
lon_corr = 3.82394E-5 * ( 365.2422 * ( Epoch - 2000.0 ) - d )
If one wishes the position for today's epoch (useful when computing
rising/setting times and the like), no corrections need to be done.
9. Perturbations of the Moon
If the position of the Moon is computed, and one wishes a higher
accuracy than about 2 degrees, the most important perturbations has
to be taken into account. If one wishes 2 arc minute accuracy, all
the following terms should be accounted for. If less accuracy is
needed, some of the smaller terms can be omitted.
First compute:
Ms, Mm Mean Anomaly of the Sun and the Moon
Nm Longitude of the Moon's node
ws, wm Argument of perihelion for the Sun and the Moon
Ls = Ms + ws Mean Longitude of the Sun (Ns=0)
Lm = Mm + wm + Nm Mean longitude of the Moon
D = Lm - Ls Mean elongation of the Moon
F = Lm - Nm Argument of latitude for the Moon
Add these terms to the Moon's longitude (degrees):
-1.274 * sin(Mm - 2*D) (the Evection)
+0.658 * sin(2*D) (the Variation)
-0.186 * sin(Ms) (the Yearly Equation)
-0.059 * sin(2*Mm - 2*D)
-0.057 * sin(Mm - 2*D + Ms)
+0.053 * sin(Mm + 2*D)
+0.046 * sin(2*D - Ms)
+0.041 * sin(Mm - Ms)
-0.035 * sin(D) (the Parallactic Equation)
-0.031 * sin(Mm + Ms)
-0.015 * sin(2*F - 2*D)
+0.011 * sin(Mm - 4*D)
Add these terms to the Moon's latitude (degrees):
-0.173 * sin(F - 2*D)
-0.055 * sin(Mm - F - 2*D)
-0.046 * sin(Mm + F - 2*D)
+0.033 * sin(F + 2*D)
+0.017 * sin(2*Mm + F)
Add these terms to the Moon's distance (Earth radii):
-0.58 * cos(Mm - 2*D)
-0.46 * cos(2*D)
All perturbation terms that are smaller than 0.01 degrees in longitude
or latitude and smaller than 0.1 Earth radii in distance are omitted.
A few of the largest perturbation terms have their own names! The
Evection, the largest perturbation, was discovered by Ptolemy a few
thousand years ago. The Variation and the Yearly Equation were both
discovered by Tycho Brahe in the 16'th century.
The computations can be simplified by omitting the smaller perturbation
terms. The error introduced by this seldom exceeds the sum of the
amplitudes of the 4-5 largest omitted terms. If one only computes the
three largest perturbation terms in longitude and the largest term in
latitude, the error in longitude will rarley exceed 0.25 degrees, and
in latitude 0.15 degrees.
10. Perturbations of Jupiter, Saturn and Uranus
The only planets having perturbations larger than 0.01 degrees are
Jupiter, Saturn and Uranus. First compute:
Mj Mean anomaly of Jupiter
Ms Mean anomaly of Saturn
Mu Mean anomaly of Uranus (needed for Uranus only)
Perturbations for Jupiter. Add these terms to the longitude:
-0.332 * sin(2*Mj - 5*Ms - 67.6 degrees)
-0.056 * sin(2*Mj - 2*Ms + 21 degrees)
+0.042 * sin(3*Mj - 5*Ms + 21 degrees)
-0.036 * sin(Mj - 2*Ms)
+0.022 * cos(Mj - Ms)
+0.023 * sin(2*Mj - 3*Ms + 52 degrees)
-0.016 * sin(Mj - 5*Ms - 69 degrees)
Perturbations for Saturn. Add these terms to the longitude:
+0.812 * sin(2*Mj - 5*Ms - 67.6 degrees)
-0.229 * cos(2*Mj - 4*Ms - 2 degrees)
+0.119 * sin(Mj - 2*Ms - 3 degrees)
+0.046 * sin(2*Mj - 6*Ms - 69 degrees)
+0.014 * sin(Mj - 3*Ms + 32 degrees)
For Saturn: ALSO add these terms to the latitude:
-0.020 * cos(2*Mj - 4*Ms - 2 degrees)
+0.018 * sin(2*Mj - 6*Ms - 49 degrees)
For Uranus: Add these terms to the longitude:
+0.040 * sin(Ms - 2*Mu + 6 degrees)
+0.035 * sin(Ms - 3*Mu + 33 degrees)
-0.015 * sin(Mj - Mu + 20 degrees)
The largest perturbation, "the great Jupiter-Saturn term", is the largest
perturbation both for Jupiter and Saturn. It has a period of 918 years, and
an amplitude of 0.332 degrees for Jupiter and 0.812 degrees for Saturn.
These is also a "great Saturn-Uranus term", period 560 years, amplitude
0.035 degrees for Uranus, less than 0.01 degrees for Saturn (and therefore
omitted). The other perturbations have periods between 14 and 100 years.
One should also mention the "great Uranus-Neptune term", which has a period
of 4220 years and an amplitude of about one degree. It is not included
here. Instead it is accounted for in the orbital elements of Uranus and
Neptune.
If the planet is Mercury, Venus or Mars, all perturbations can be neglected.
For Neptune the only significant perturbation is already included in the
orbital elements, as mentioned above, and therefore no further perturbation
terms need to be accounted for.
11. Geocentric (Earth-centered) coordinates
Now we have computed the heliocentric (Sun-centered) coordinate of the
planet, and we have accounted for the most important perturbations.
We want to compute the geocentric (Earth-centerd) position. We should
now convert the perturbed lonecl, latecl, r to (perturbed) xh, yh, zh:
xh = r * cos(lonecl) * cos(latecl)
yh = r * sin(lonecl) * cos(latecl)
zh = r * sin(latecl)
If we are computing the Moon's position, this is already the geocentric
position, and thus we simply set xg=xh, yg=yh, zg=zh. Otherwise we must
also compute the Sun's position, and convert lonsun, rs to xs, ys:
xs = rs * cos(lonsun)
ys = rs * sin(lonsun)
(Of course, any correction for precession should be added to lonecl AND
lonsun BEFORE converting to xh,yh,zh and xs,ys).
Now convert from heliocentric to geocentric position:
xg = xh + xs
yg = yh + ys
zg = zh
Now we have the planet's geocentric (Earth centered) position in
rectangular, ecliptic coordinates.
12. Equatorial coordinates
Now let's convert our rectangular, ecliptic coordinates to rectangular,
equatorial coordinates: simply rotate the y-z-plane by ecl, the angle of
the obliquity of the ecliptic:
xe = xg
ye = yg * cos(ecl) - zg * sin(ecl)
ze = yg * sin(ecl) + zg * cos(ecl)
Finally, compute the planet's Right Ascension (RA) and Declination (Dec):
RA = atan2( ye, xe )
Dec = atan2( ze, sqrt(xe*xe+ye*ye) )
Compute the geocentric distance:
rg = sqrt(xg*xg+yg*yg+zg*zg) = sqrt(xe*xe+ye*ye+ze*ze)
13. The position of Pluto
Calculate d, our day number, as usual, Then calculate these three angles:
J = 34.23 + 0.083091190 * d
S = 50.03 + 0.033459652 * d
P = 238.95 + 0.003968789 * d
Next compute the heliocentroc ecliptic longitude and latitude (degrees),
and distance (a.u.):
lonecl = 238.9508 + 0.00400703 * d
- 19.799 * sin(P) + 19.848 * cos(P)
+ 0.897 * sin(2*P) - 4.956 * cos(2*P)
+ 0.610 * sin(3*P) + 1.211 * cos(3*P)
- 0.341 * sin(4*P) - 0.190 * cos(4*P)
+ 0.128 * sin(5*P) - 0.034 * cos(5*P)
- 0.038 * sin(6*P) + 0.031 * cos(6*P)
+ 0.020 * sin(S-P) - 0.010 * cos(S-P)
- 0.004 * sin(S) - 0.005 * cos(S)
- 0.006 * sin(S+P) - 0.003 * cos(S+P)
+ 0.007 * sin(J-P) + 0.001 * cos(J-P)
latecl = -3.9082
- 5.453 * sin(P) - 14.975 * cos(P)
+ 3.527 * sin(2*P) + 1.673 * cos(2*P)
- 1.051 * sin(3*P) + 0.328 * cos(3*P)
+ 0.179 * sin(4*P) - 0.292 * cos(4*P)
+ 0.019 * sin(5*P) + 0.100 * cos(5*P)
- 0.031 * sin(6*P) - 0.026 * cos(6*P)
+ 0.005 * sin(S-P) + 0.011 * cos(S-P)
r = 40.72
+ 6.68 * sin(P) + 6.90 * cos(P)
- 1.18 * sin(2*P) - 0.03 * cos(2*P)
+ 0.15 * sin(3*P) - 0.14 * cos(3*P)
+ 0.05 * cos(4*P)
- 0.01 * sin(5*P) - 0.01 * cos(5*P)
Now we know the heliocetric distance and ecliptic longitude/latitude
for Pluto. To convert to geocentric coordinates, do as for the other
planets.
14. Positions of asteroids
For asteroids, the orbital elements are often given as: N,i,w,a,e,M,
where N,i,w are valid for a specific epoch, usually 1950.0. In our
simplified computational scheme, the only significant changes with the
epoch occurs in N. To convert N_Epoch to the N (today's epoch) that we
want to use, simply do:
N = N_Epoch + 0.013967 * ( 2000.0 - Epoch ) + 3.82394E-5 * d
where Epoch is expressed as a year with fractions, e.g. 1950.0 or 2000.0
Most often M, the mean anomaly, is given for another day than the day
we want to compute the asteroid's position for. If the daily motion, n,
is given, simply add n * (time difference in days) to M. If n is not
given, but the period P (in days) is given, then n = 360.0/P. If P is
not given, it can be computed from:
P = 365.2568984 * a**1.5 (days) = 1.00004024 * a**1.5 (years)
** is the power-of operator. a**1.5 is the same as sqrt(a*a*a).
Now when all orbital elements has been computed, proceed as with the
other planets (goto section 6).
15. Position of comets.
For comets with elliptical orbits, M is usually not given. Instead T,
the time of perihelion, is given. At this moment, M is zero. To
compute M for any other moment, first compute the "day number" d of T,
let's call this dT. Then compute the "day number" d of the moment for
which you want to compute a position, let's call this d. Then M, the
mean anomaly, is computed like:
M = 360.0 * (d-dT)/P (degrees)
where P is given in days.
Also, a, the semi-major axis, is not given. Instead q, the perihelion
distance, is given. a can simply be computed from q and e:
a = q / (1.0 - e)
Then correct N for the epoch, if any, and proceed like for an asteroid or
a planet (goto section 6).
16. Parabolic orbits
If the comet has a parabolic orbit, a different method has to be used.
Now the comet has an infinite orbital period, and M, the mean anomaly,
is always zero. The eccentric anomaly, E, cannot be computed. The
eccentricity, e, is always exactly 1 and is therefore of no use. To
compute a parabolic orbit, do like this:
Compute the "day numbers" for T, the moment of perihelion, call it dT.
Compute d for the moment we want to compute a position, call it d.
Then compute:
H = (d-dT) / (82.21168627 * q**1.5)
where q**1.5 is the same as sqrt(q*q*q). Also compute:
h = 1.5 * H
g = sqrt( 1.0 + h*h )
s = cbrt( g + h ) - cbrt( g - h )
cbrt() is the cube root function. The formulae has been devised so that
both g+h and g-h always are positive. Therefore one can safely compute
cbrt(x) as exp( log(x)/3.0 ) . In general, cbrt(-x) = -cbrt(x) and of
course cbrt(0) = 0.
Compute the true anomaly and the heliocentric distance:
v = 2.0 * atan(s)
r = q * ( 1.0 + s*s )
Now that we know the true anomaly and the heliocentric distance, we can go
on to compute the position in space (goto section 7).
17. Near-parabolic orbits.
If the comet's orbit is elliptic, but so nearly parabolic that the
computation for an elliptic orbit never converges or becomes inaccurate,
or if the comet's orbit is slightly hyperbolic (as some comets' orbits
are), one can usually pretend that the orbit is exactly parabolic. The
errors in the geocentric position due to this assumption rarely exceeds
one or a few arc minutes while the comet resides in the inner solar
system. One can usually pretend that the comet's orbit is exactly
parabolic if the eccentricity is somewhere between about 0.998 and 1.002.
18. Validity of orbital elements.
Due to perturbations from mainly the heavy planets like Jupiter and Saturn,
the orbital elements of any heavenly body are constantly changing. The
orbital elements for the Sun, Moon and the major planets, as given here,
are valid for a long time period. However, orbital elemets given for
a comet or an asteroid (or for Pluto) are valid only for a limited time.
How long they are valid is hard to say generally. It depends on many
factors, such as the accuracy you need, and the magnitude of the
perturbations the comet or asteroid is subjected to from, say, Jupiter.
A comet might travel in roughly the same orbit several orbital periods,
experiencing only slight perturbations, but suddenly it might pass very
close to Jupiter and get its orbit changed drastically. To compute this
in a reliable way is quite complicated and completely out of scope for
this description. As a rule of thumb, though, one can assume that an
asteriod, if one uses the orbital elements for a certain moment, one
or a few revolutions away from that moment will have an error in its
computed position of at least one or a few arc minutes, and possibly
more. The errors will accumulate with time.
--
----------------------------------------------------------------
Paul Schlyter, Swedish Amateur Astronomer's Society (SAAF)
Grev Turegatan 40, S-114 38 Stockholm, SWEDEN
e-mail: pausch@saaf.se psr@home.ausys.se