TWO - BODY - PROBLEM

[Link], IIC, UDSC

It is a laudable human pursuit to try to perceive order out of the apparent
randomness of nature; science is , after all , an attempt to make sense of
the world around us. Moving against the background of the xed stars,
the regularity of the Moon and planets demanded a dynamical explanation.
The history of astronomy is the history of a growing awareness of our posi-
tion ( or lack of it) in the universe. Observing , exploring , and ultimately
understanding our solar system is the rst step towards understanding the
rest of the universe. The key discovery in this process was Newtons formu-
lation of the universal law of gravitation; this made sense of the orbits of
planets , satellites, and comets , and their future motion could be predicted:
The Newtonian universe was a deterministic system. The Voyager missions
increased our knowledge of the outer solar system by several orders of mag-
nitude, and yet they would not have been possible without the knowledge of
Newtons laws and their consequences. However, advances in mathematics
and computer technology have now revealed that , even though our system
is deterministic, it is not necessarily predictable. The study of nonlinear
dynamics has revealed a solar system even more intricately structured than
Newton could have imagined.
MATHEMATICAL FORMULATION OF TWO-BODY PROB-
LEM
The two-body problem is perhaps the simplest, integrable problem in solar
system dynamics. It concerns the interaction of two point masses moving
under a mutual gravitational attraction described by Newtons universal
law of gravitation. The wide variety of masses in the solar system permits
the orbits of most planets and satellites to be approximated by two-body
motion, consisting of a smaller body moving around a much larger central
body. The eect of other bodies can usually be treated as perturbations to
the two-body system.
For example, the path of jupiter ( mass m
J
= 1.9 10
27
kg) around the
Sun ( mass m
Su
= 2.0 10
30
kg 1000m
J
) is basically an ellipse with
principal perturbation coming rom the other planets notably Saturn ( mass
m
Sa
= 5.7 10
26
kg). Consider two particles of mass m
1
and m
2
with
position vector r
1
and r
2
referred to some origin O xed in the inertial
1
space ( Fig.1). Then the gravitational forces and consequent acceleration
experienced are given by:
F
1
= G
m
1
m
2
r
3
r = m
1
r
1
F
2
= G
m
1
m
2
r
3
r = m
2
r
2
(1)

m
2

m
1

O
r

r
r

F
2

F
1

Fig.1
where r = r
2
r
1
and G = 6.67260 10
11
Nm
2
kg
2
is the universal grav-
itational constant. Observe that
m
1
r
1
+ m
2
r
2
= 0 (2)
which on integration results in
m
1
r
1
+ m
2
r
2
= a
m
1
r
1
+ m
2
r
2
= at +b (3)
where a and b are constant vectors. If the position vector of the center of
mass is R
G
then
R
G
=
m
1
r
1
+ m
2
r
2
m
1
+ m
2
=
at +b
m
1
+ m
2

R
G
=
m
1
r
1
+ m
2
r
2
m
1
+ m
2
=
a
m
1
+ m
2
(4)
2
Therefore the center of mass is either stationary (a = 0) or it is moving
with a constant velocity with respect to origin O. In terms of r = r
2
r
1
,
we may write (1) as the equation of relative motion.
r
2
r
1
=
G(m
1
+ m
2
)
r
3
r =

r
3
r
r +
r
r
3
= 0. (5)
In order to solve the foregoing equation and nd the path of m
2
relative to
m
1
, we must rst derive several constants of the motion.
Taking the vector product of r with (5) we nd that r r = 0 which on
integration gives
r r = h (6)
where h is a constant vector perpendicular to both r and r. Thus the motion
of m
2
about m
1
lies in a plane perpendicular to the direction dened by h
(Fig.2). Eq.(6) is known as angular momentum integral.













Fig.2
Since r and r always lie in the same plane ( the orbit plane) it is natural
that we now restrict ourselves to considering motion in that plane. We now
transform to a polar coordinate system ( r, ) referred to an origin centered
on the mass m
1
and an arbitrary reference line corresponding to = 0. If
we let r and

denote the unit vectors along and perpendicular to the radius
vector, then the position, velocity , and acceleration can be written in polar
3
coordinates as:
r = rr, r = rr + r

r = ( r r

2
)r +
_
1
r
d
dt
(r
2

)
_

(7)
Substituting the expression for r into eq.(6) gives
h = r
2

z
where z is a unit vector perpendicular to the plane of orbit forming right-
handed triad with r and

. Therefore
h = r
2

(8)
Consider the motion of m
2
during a time interval t (Fig.3). At time t = 0
it has polar coordinates (r, ), while at time t+t its polar co-ordinates have
changed to (r +r, +). The area swept out by the radius vector in time
t is
A
1
2
r(r + r) sin
1
2
r
2
(9)
Hence, by dividing each side by t and taking the limit as t 0 we have
dA
dt
=
1
2
r
2
d
dt
=
1
2
h (10)
Since h is a constant this implies that equal areas are swept out in equal
times and hence eq.(10) is the mathematical form of Keplers second law
of planetary motion. Observe that this does not require an inverse
square law of force , but only that the force is directed along the
line joining the two masses. Orbital Position and velocity
A scalar equation for the relative motion can be obtained by substituting
the expression for r from eq.(7) into (5). Comparing the r components gives
r r

2
=

r
2
(11)
To solve this equation and nd r as a function of we need to make the
substitution
u =
1
r
and to eliminate the time by making use of the relation
h = r
2

4























m
1

A
r + r
r
t = t
m
2

t=0
Fig.3
to obtain the following second order linear dierential equation
d
2
u
d
2
+ u =

h
2
(12)
The solution of the foregoing equation may be written as:
u =

h
2
[1 + e cos( )] (13)
where e- (an amplitude) and - (a phase) are two constant of integration.
In terms of r we may write
r =
p
1 + e cos( )
, p =
h
2

(14)
which is a general equation of conic in polar coordinate with e as its eccen-
tricity and p as semilatus rectum. The four possible conics are
circle : e = 0, p = a
ellipse 0 < e < 1, p = a(1 e
2
)
parabola : e = 1, p = 2q
hyperbola : e > 1, p = a(e
2
1)
(15)
where the constant a is the semimajor axis of the conic. In the special case
of the parabola p is dened in terms of q -the distance to the central mass
at closest approach.
5
In the context of two body problem the path of a planet about the Sun is
elliptical and closed in inertial space (Fig.4) and hence Keplers rst law of
planetary motion motion is a consequence of the inverse square law of force.
Observe that m
1
lls one focus of the ellipse while the other focus is empty.

= f +
Empty
focus
apocentr
e
pericentre
Reference
direction
b

m
1

ae
f
a
m
2

x
y
Fig.1
r
Fig.4
Although a large number of cometary orbits have e 1, most permanent
members of the solar system have e << 1. the notable exception among
the planets are Pluto (e = 0.25) and mercury (e = 0.21), while Nereid
(e = 0.75) a moon of Neptune has the largest eccentricity of any known
natural satellite. Consequently we consider only elliptical motion only such
that p = a(1 e
2
) and the quantities a and e are related by
b
2
= a
2
(1 e
2
)
where b is the semi-minor axis of the ellipse. We also have now
r =
a(1 e
2
)
1 + e cos( )
(16)
In celestial mechanics, it is customary to use the term longitude when refer-
ring to any angle that is measured with respect to a reference line xed in
inertial space. The angle is called the true longitude. It is readily observed
from above that the minimum and maximum of the orbital radius are
r
p
= a(1 e) r
a
= a(1 + e)
6
which occurs when = and = + respectively. These points in the or-
bit are called pericentre (or periapse) and apocentre (or apoapse)respectively.
The angle is termed as longitude of pericentre. Although this a constant for
the two-body problem, it can vary with time when additional perturbations
are introduced. It is usually more convenient to refer the angular coordinate
to the pericentre rather than to the arbitrary reference line. This leads to
the introduction of the angle f = which is called true anomaly. Since
is a constant the path is closed and the angular position is described by
f or , which are 2 periodic variables. Therefore the position of a particle
is described by
r =
a(1 e
2
)
1 + e cos f
(17)
Using a Cartesian coordinate system centred on the central mass with the x-
axis pointing towards the pericentre, the components of the position vector
are
x = r cos f, y = r sin f (18)
in one orbital period T the area swept out by a radius vector is simply the
area A = ab enclosed by the ellipse. From eq.(10) this area has to equal
hT/2 and hence, given that h
2
= a(1 e
2
),
T
2
=
4
2

a
3
(19)
which corresponds to Keplers third law of planetary motion. Note that the
period of the orbit is independent of e and is a function of and a only.
Consider the case of two objects of mass m and m

, orbiting a central
object of mass m
c
. Let the orbiting objects have semi-major axes a and a

and orbital period T and T

. Then
m
c
+ m
m
c
+ m

=
_
a
a

_
3
_
T

T
_
2
(20)
In the case of planets orbiting the Sun we have m, m

<< m
c
and hence
(a/a

)
3
(T/T

)
2
.
If any solar system object (e.g. asteroid or a comet) has a small natural or
articial satellite, then observations of the distance and period of the satellite
can be used with Keplers third law to derive an estimate of the mass of the
7
object. Consider the Sun-object and object- satellite system. Let m
c
, m,
and m

denote the masses of the Sun, object, and satellite respectively with
similar denitions of the semi-major axes and orbital periods.
Then
m + m

m
c
+ m

m
m
c
=
_
a

a
_
3
_
T
T

_
2
(21)
This means that the mass of the object (i.e., an asteroid or a comet) can
be estimated from the orbital properties of its satellite. Observation of
Dactyls motion [Dactyl is a moon of the asteroid Ida] taken by Galileo
spacecraft on its way to Jupiter have resulted in an estimated density of
2.6 [Link]
3
( cf. Belton [Link]., 1995, Nature, 374,785-788).
Since the angle covers 2 radians in one orbital period, we can dene the
average angular velocity, or the mean motion, n as
n =
2
T
(22)
and we can write
= n
2
a
3
, h = na
2
_
1 e
2
=
_
a(1 e
2
).
(23)
Although the mean motion is constant in the two-body problem, the actual
angular velocity

f of the orbiting body is a function of the longitude.
Now consider once again the eq.(5) derived earlier
r +
r
r
3
= 0 (24)
and take the scalar product with r to get
r.r +
r
r
2
= 0 (25)
which can be integrated to give
1
2
v
2

r
= C (26)
where v
2
= r. r is the square of the velocity and C is a constant of the motion.
The foregoing equation shows that the orbital energy per unit mass is con-
served. Thus the two-body problem has four constants of the motion: the
8
energy integral C and three components of the angular momentum integral,
h.
Expression for v
2
Substituting r from eq.()into the expression
v
2
= r. r
we get
v
2
= r. r = r
2
+ r
2

f
2
(27)
Further dierentiating eq.(20), we have
r =
r

fe sin f
1 + e cos f
(28)
Using r
2

f = h = na
2

1 e
2
, we can write
r =
na

1 e
2
e sinf (29)
and
r

f =
na

1 e
2
(1 + e cos f) (30)
Now substitute for r and r

f in the expression for v
2
, we get
v
2
=
n
2
a
2
1 e
2
_
1 + 2e cos f + e
2
_
=
n
2
a
2
1 e
2
_
2a(1 e
2
)
r
(1 e
2
)
_
(31)
Hence
v
2
=
_
2
r

1
a
_
(32)
where we have substituted = n
2
a
3
. Comparing with earlier expression
obtained for v
2
, we nd that
C =

2a
(33)
9
Hence for the elliptical orbit the energy is a function of its semi-major axis
alone and is independent of eccentricity, e.
Also the velocity of the orbiting body at pericentre(f = 0) and apocentre(f =
) is given by:
v
p
= na

1 + e
1 e
and v
a
= na

1 e
1 + e
(34)
Further since x = r cos f and y = r sin f, we nd the following expression
for x and y ( using the expressions for r and r

f)as:
x =
na

1 e
2
sin f (35)
y = +
na

1 e
2
[e + cos f] (36)
The Mean and Eccentric Anomalies
: Kepler Equation
The following equation summarizes the outcome of the discussions so far.
r =
a(1 e
2
)
1 + e cos f
x = r cos f, y = r sin f
r
p
= a(1 e) r
a
= a(1 + e)
v
2
=
_
2
r

1
a
_
(37)
It is readily observed that given the true anomaly f, we can calculate
the orbital radius and velocity of a body provided we know the eccentricity
and semi-major axis of its orbit.
However, in practice we usually want to calculate the location of a body at
a given time explicitly. Although f and r are function of t, we have not
shown the nature of this dependence.
In order to do so, we now introduce the term mean anomaly, M. We
dene M in terms of the average angular velocity or mean motion n as:
M = n(t ) (38)
10
where is the time of pericentric passage. Although M has the dimensions
of an angle, and it increases linearly with time at a constant rate equal to
the mean motion , it has no simple geometrical interpretation.
It may be noted in the foregoing denition of mean anomaly M that when
t = (pericentre passage), M = f = 0 and when t = + T/2 (apocentre
passage), M = f = .
Consider now a circumscribed circle of radius a that is concentric with an
orbital ellipse of semi-major axis a and eccentricity e (Fig.5).








F O
a
E f
r
( ) , x y
Circumscribe
d
circle
ellips
e
Fig.5
In this gure, a line perpendicular to the major axis of the ellipse is extended
through the point on the orbit and intersect the circle. We next dene the
parameter E- the eccentric anomaly to be the angle between the major
axis of the ellipse and the radius from the centre to the intersection
point on the circumscribed circle. Hence E = 0 corresponds to f = 0
and E = corresponds to f = .
Now the equation of a centered ellipse in rectangular coordinates is
_
x
a
_
2
+
_
y
b
_
2
= 1 (39)
11
From Fig.5, we observe that
x = a cos E, y
2
= b
2
sin
2
E y = a
_
1 e
2
sin E
Therefore the projection of r in the horizontal and vertical directions are:
x = a [cos E e]
y = a
_
1 e
2
sin E (40)
Using the foregoing equations, we nd the following expressions for r and
cos f:
r = a(1 e cos E), cos f =
cos E e
1 e cos E
(41)
a relation between E and f could be obtained by observing that
1 cos f =
(1 + e)(1 e cos E)
1 e cos E
2 sin
2
f
2
=
(1 + e)
1 e cos E
2 sin
2
E
2
(42)
Similarly writing the expression for 1 + cos f, we nd that
2 cos
2
f
2
=
(1 e)
1 e cos E
2 cos
2
E
2
(43)
Therefore using the foregoing equations, we obtain
tan
f
2
=

1 e
1 e
tan
E
2
(44)
Thus knowing eccentric anomaly E, we can determine r and f uniquely (
since E and f will always lie in the same half of the ellipse). However,
to locate a body in its orbit at some time t, we need to derive a
relationship between M and E.
Now consider the following equations:
v
2
= r
2
+ r
2

f
2
r =
na

1 e
2
e sin f
r

f =
na

1 e
2
(1 + e cos f)
v
2
=
_
2
r

1
a
_
12
Using these equations, we nd that
r
2
= n
2
a
3
_
2
r

1
a
_

n
2
a
4
(1 e
2
)
r
2
r
2
=
n
2
a
2
r
2
_
a
2
e
2
(r a)
2
_

dr
dt
=
na
r
_
a
2
e
2
(r a)
2
. (45)
Eq.(45) by substituting r a = ae cos E can be reduced to the form:
dE
dt
=
n
1 e cos E
(46)
which on integration gives
n(t ) = E e sin E (47)
where is a constant of integration such that at t = , E = 0. Further
using the denition of mean anomaly M, we have the Kepler equation:
M = E e sinE (48)
So far we have dened the true longitude (), true anomaly (f), the mean
anomaly (M), the eccentric anomaly (E), and the longitude of pericentre
(). To complete this set we dene the mean longitude() by
= M + . (49)
Therefore is a linear function of time and , since it is derived from M,
it has no simple geometrical interpretation, except in the special case of a
circular orbit. It is important to note that all longitudes (, , ) are dened
with respect to a common, arbitrary direction.
Solution of Keplers equation
Keplers equation cannot be solved directly because it is transcendental in
E and therefore, apart from the trivial solutions E = j when M = j for
integer j, it is impossible to express E as a simple function of time.
(a) Series solution of KE using iterative method
We may write the Keplers equation as:
E
i+1
= M e sin E
i
i = 0, 1,
(50)
13
where E
i
is the i-th approximation of the value E. In the foregoing , we
assume an initial guess or approximation such that E
0
= M. Therefore, we
nd that
E
1
= M + e sin M
E
2
= M + e sin(M + sin M)
M + e sin M +
1
2
e
2
sin 2M
E
3
= M + e sin(M + e sin M +
1
2
e
2
sin 2M)
E
3
M +
_
e
1
8
e
3
_
sinM
+
1
2
e
2
sin 2M +
3
8
e
3
sin 3M (51)
for the rst three steps, where we introduced only one additional term in e
at each step. It is clear from this approach that the nal series for E M
will have the form
E M =

s=1
b
s
(e) sin sM
where the lowest order term in b
s
(e) is O(e
s
).
Danbys (1988) method of solution of KE
By writing the Keplers equation as:
f(E) = E e sin E M (52)
we can use the Newton-Raphson method to nd the root of the nonlinear
equation f(E) = 0. The iteration scheme is
E
i+1
= E
i

f(E
i
)
f

(E
i
)
, i = 0, 1, 2,
(53)
Danby(1988) points out that the convergence of the Newton-Raphson scheme
is quadratic but the quartic convergence is also possible with a modied
scheme. Using Danbys notation and a Taylor series expansion we can write
0 = f(E
i
+
i
) = f(E
i
) +
i
f

(E
i
)
+
1
2

2
i
f

(E
i
) +
1
6

3
i
f

(E
i
) + O(
4
i
)
14
Neglecting higher order terms in , we can write
0 = f
i
+
i
f

i
+
1
2

2
i
f

i
+
1
6

3
i
f

i
where f
i
= f(E
i
), f

i
= f

(E
i
), etc. Hence

i
=
f
i
f

i
+
1
2

i
f

i
+
1
6

2
i
f

i
This can be solved for
i
by dening

i,1
=
f
i
f

i,2
=
f
i
f

i
+
1
2

i,1
f

i,3
=
f
i
f

i
+
1
2

i,2
f

i
+
1
6

2
i,2
f

i
and then using the iteration scheme
E
i+1
= E
i
+
i,3
(54)
Although this method has more arithmetic operations per iteration than
the standard Newton-Raphson scheme given above, it is more ecient since
(a) it can be programmed to make use of quantities that have already been
calculated at each iteration and (b) it will converge faster.
An important consideration in either of these numerical schemes is a suitable
starting value, E
0
. Obviously for small e we have E M and so E
0
= M
seems appropriate. However this guess is only correct in the cases where
e = 0 or M is a multiple of . Danby(1988) points out that, by rst reducing
M to the range 0 M 2, the initial guess
E
0
= M + sign(sin M)ke, 0 k 1
(55)
has a better chance of being correct and improves the convergence; the
recommended value is k = 0.85.
Solution for position and velocity
The solution of Keplers equation to nd E for a given value of M allows
the calculation of the position and velocity at any point t for an object in an
15
elliptical orbit. If the object has a position vector r
0
= r(t
0
) and a velocity
v
0
= v(t
0
) at time t
0
then this process can be simplied by introduction
of two special functions and their time derivatives. Provided that initial
vectors r
0
and v
0
are not parallel, r(t) can be written as:
r(t) = f(t, t
0
)r
0
+ g(t, t
0
)v
0
(56)
where f(t, t
0
) and g(t, t
0
) are referred to as f and g functions.
Separating the x and y components, we nd that
x = f(t, t
)
x
0
+ g(t, t
0
) x
0
y = f(t, t
)
y
0
+ g(t, t
0
) y
0
where we have taken r
0
= (x
0
, y
0
) and v
0
= ( x
0
, y
0
). This gives
f(t, t
0
) =
x y
0
y x
0
x
0
y
0
y
0
x
0
g(t, t
0
) =
y x
0
x y
0
x
0
y
0
y
0
x
0
Using the following relations derived earlier:
x = a(cos E e)
y = a
_
1 e
2
sin E
r = a(1 e cos E)
dE
dt
=
n
1 e cos E
we obtain the following expressions for f,

f and g, g functions:
f(t, t
0
) =
a
r
0
[cos(E E
0
) 1] + 1
g(t, t
0
) =
1
n
[sin(E E
0
) (E E
0
)]
+(t t
0
)

f(t, t
0
) =
a
2
rr
0
nsin(E E
0
)
g(t, t
0
) =
a
r
0
[cos(E E
0
) 1] + 1
The use of f and g functions means that once E is known from the solution
of Keplers equation, we can readily nd r and v.
16
Orbit in Space
So far we have shown that the position and velocity vectors of the mass m
2
with respect to mass m
1
always lie in a plane perpendicular to the angular
momentum vector. The values of r = (x, y) and r = ( x, y) or alternatively
(r, , r,

) of the mass m
2
with respect to m
1
at any time dene a unique orbit
and a location on that orbit by means of the three constants a, e,and and
the variable f. Our subsequent analysis was concerned with understanding
the motion in the orbital plane. However the motion in the solar system
is not conned to a single plane and we now consider the three- dimensional
representation of an orbit in space ( Fig.6).
In Fig.6 an arbitrary point has a position vector r = x x + y y + z z. The
x-axis is taken to lie along the major axis of the ellipse in the direction of
pericentre, the y-axis is perpendicular to the x-axis and lies in the orbital
plane, while the z-axis is mutually perpendicular to the x and y axes such
that all three form a right-handed triad.
We now wish to refer this orbital plane to a standard plane dened by
(

X,

Y,

Z) as shown in Fig.6. It may be noted that when considering the

motion of planets around the Sun, it is customary to use Sun-centered or
heliocentric coordinate system where the reference plane is the plane of the
Earths orbit (the ecliptic) and the reference line is in the direction of vernal
equinox, along the line of intersection of the plane of the earths equator and
the ecliptic.
In general the orbital plane will be inclined to the reference plane at an angle
I called the inclination of the orbit. The line of intersection between the
orbital plane and the standard reference frame is called the line of nodes.
The point in both planes where the orbit crosses the reference plane moving
from below to above the plane is called the ascending node while the angle
between the reference line and the radius vector to the ascending node is
called the longitude of ascending node,. The angle between this same radius
vector and the pericentre of the orbit is called the argument of pericentre,
.
The inclination is always in the range 0 I 180

. If I < 90

the motion
is said to be prograde whereas if I 90

the motion is retrograde.

In the limit as I 0 the orbital plane coincides with the reference plane
and we have
= + .
However , the denition of above is also used in the inclined case, despite
the fact that the angles and lie in dierent planes. Fig.7 shows the
17
relationship between the orbital plane coordinate system and the reference
plane system. It is clear that coordinates in one system can be expressed
in terms of the other by means of a series of three rotations about various
axes.
To transform from (x, y, z) orbital plane to the general (X, Y, Z) reference
system we have to carry out (i) a rotation about the z-axis through an angle
so that the x-axis coincides wit the line of nodes, (ii) a rotation about the
x-axis through an angle I so that the two plane are coincident, and nally
(iii) a rotation about the z-axis through an angle ( see Fig.6 and Fig.7).

Fig.6
We represent these transformation by three 33 rotation matrices, denoted
by P
1
, P
2
and P
3
respectively. Here
P
1
=
_
_
_
cos sin 0
sin cos 0
0 0 1
_
_
_
P
2
=
_
_
_
1 0 0
0 cos I sin I
0 sinI cos I
_
_
_
and
P
3
=
_
_
_
cos sin 0
sin cos 0
0 0 1
_
_
_
18








Fig.7
Consequently
_
_
_
X
Y
Z
_
_
_ = P
3
P
2
P
1
_
_
_
x
y
z
_
_
_
_
_
_
x
y
z
_
_
_ = P
1
1
P
1
2
P
1
3
_
_
_
X
Y
Z
_
_
_ (57)
If we now restrict ourselves to coordinates that lie in the orbital plane, we
have
_
_
_
X
Y
Z
_
_
_ = P
3
P
2
P
1
_
_
_
r cos f
r sin f
0
_
_
_ (58)
= r
_
_
_
cos cos( + f) sin sin( + f) cos I
sin cos( + f) + cos sin( + f) cos I
sin( + f) sin I
_
_
_
19
Observe that the value of a and e are unchanged by considering the ellipse
in this new coordinate system, since rotational transformations preserve
lengths.
We next attempt to use these formulae to nd the position of the plan-
ets at a given time, say September,25, 1994 at 5.32PM British
Summer time.
Julian Date
To calculate the position of a solar system body at a particular time it is
necessary to make use of the calender system based on a xed day, the
Julian day, consisting of 86,400s. Similarly , the Julian year consists of
365.25 Julian days and the Julian Century has 36,525 Julian days.
The Julian date (JD)is the number of Julian days since 12
h
Universal Time
( noon Greenwich) on 1 January 4713 B.C. The Julian date for any cal-
ender date can be calculated by using the following algorithm by Mon-
tenbruck(1989)( Practical Ephemeris Calculations, Springer-Verlag, Heidel-
berg). If Y , M, D and UT denote the year , month, day, and universal
Time, then the rst step is to dene the auxillary quantities y and m using
y = Y 1 and m = M + 12 if M 2,
y = Y and m + M if M > 2 (59)
and the quantity B using
B = 2, up to and including 4 Oct. 1582
B = Int[y/400] , from and including 15 Oct 1582
Int[y/100] (60)
where the function Int[x] denotes the largest whole number that is smaller
than or equal to x. The reason for the peculiar denition of B are to account
for the lost days in October 1582, when the Gregorian calender replaced
the Julian calender in Europe, and to deal with the introduction of a leap
day in the Gregorian calender.
The Julian date is given by
JD = Int[365.25y] + Int[30.6001(m + 1)]
+B + 1720996.5 + D + UT/24.
As an example, consider the calculation of the Julian date corresponding to
10
h
24
m
on 4 February, 1946. In this case Y = 1946 , M = 2, D = 4 and
20
UT = 10.4. The auxillary quantities are y = 1945 , m = 14 , and B = 15,
giving
JD = 2431855.933
Therefore for September 25, 1994 at 5.32 PM British Summer Time (it is
increased by an hour), we nd that Y = 1994, M = 9 , D = 25 and
UT = 16.32. Based on this the auxillary quantities y = 1993 , m = 9 and
B = 15. Putiing these values in the formula for JD , we get
JD = 2449256.18
Now for the J2000 epoch i.e., noon of Ist January, 2000, we observe that
JD = 2451545.0. Therefore the interval t in centuries between the given
date and the J2000 epoch date is 2288.811/36525 or 0.062664229 or
0.06266423.
For the orbital elements i.e. a, e I, , and , of a planet we use the
following table (cf. Table:1a). the orbital elements of the planets change
over time owing to their mutual perturbations. Table:1a and Table:1b give
the orbital elements of planets and their variations at the epoch of J2000 with
respect to the mean ecliptic and equinox of J2000 ( cf. Standish et al.(1992)).
To calculate the approximate elements at other times the following formulae
are used:
a = a
0
+ atAU,
e = e
0
+ et,
I = I
0
+ (

I/3600)t degrees
=
0
+ ( /3600)t degrees
=
0
+ (

/3600)t degrees
=
0
+ (

/3600 + 360N
r
)t degrees
Using the following table, we observe that for Jupiter at the epoch of J2000,
a
0
= 5.20336301, e
0
= 0.04839266, I
0
(

) = 1.30530,
0
(

) = 14.75385,

0
(

) = 100.55615 and
0
= 34.40438. Therefore using the data of Table:1b,
we observe that for Jupiter on September 25, 1994
a = 5.2033601 60737 0.06266423 10
8
= 5.20332
e = 0.04839266 +
(12880 0.06266423 10
8
)
21
= 0.04840073
I = 1.30530 + (4.15/3600) 0.06266423
= 1.305372

= 14.75385 +
(839.93/3600) 0.06266423
= 14.7392295

= 100.55615 +
(1217.17/3600) 0.06266423
= 100.53496305

= 34.40438 +
(557078.35/3600 + 360 8) 0.06266423
= 155.765515 = 204.234485

Therefore M = = 189.059

. The numerical solution of Kepler equation

gives E = 189.059

. We therefore obtain
x = a(cos E e) = 5.390263AU
y = a
_
1 e
2
sin E = 0.81831AU
z = 0.0AU
Substituting the values of I, and (= ) in eq.(57), we get
P = P
3
P
2
P
1
=
_
_
_
0.966838 0.254401 0.022397
0.254373 0.967097 0.004165
0.02272 0.00167 0.99974
_
_
_
for the transformation matrix and hence the coordinates of Jupiter in the
J2000 reference frame are
X = 0.996838x 0.254401y + 0.022397z
= 5.16504AU
Y = 0.254373x + 0.967097y + 0.004165z
= 2.162523AU
Z = 0.02272x + 0.00167y + 0.99974z
= 0.1211AU
22

Fig.8a
23
The procedure described here can be applied to nd the positions of the
other planets. The results are shown in Fig.8

Fig.8b
24
Table:1a*. Planetary orbital elements at the epoch
J2000(JD2451545.0)with respect to the mean ecliptic and
equinox of J2000.
Planet a
0
(AU) e
0
I
0
(

)
0
(

)
0
(

)
0
(

)
Mercury 0.38709893 0.20563069 7.00487 77.45645 48.33167 252.25084
Venus 0.72333199 0.00677323 3.39471 131.53298 76.68069 181.97973
Earth 1.00000011 0.01671022 0.00005 102.94719 348.73936 100.46435
Mars 1.52366231 0.09341233 1.85061 336.04084 49.57854 355.45332
Jupiter 5.20336301 0.04839266 1.30530 14.75385 100.55615 34.40438
Saturn 9.53707032 0.05415060 2.48446 92.43194 113.71504 49.94432
Uranus 19.19126393 0.04716771 0.76986 170.96424 74.22988 313.23218
Neptune 30.06896348 0.00858587 1.76917 44.97135 131.72169 304.88003
Pluto 39.48168677 0.24880766 17.14175 224.06676 110.30347 238.92881
*The data for the earth are for the Earth-Moon barycentre
Table:1b*. Rate of change of planetary orbital elements at the
epoch J2000(JD2451545.0)with respect to the mean ecliptic and
equinox of J2000.
Planet a
0
10
8
(AU) e
0
10
8

I
0
(

)
0
(

0
(

0
(

) N
r
Mercury 66 2527 -23.51 573.57 -446.30 261628.29 415
Venus 92 -4938 -2.86 -108.80 -996.89 712136.06 162
Earth -5 -3804 -46.94 1198.28 -18228.25 1293740.63 99
Mars -7221 11902 -25.47 1560.78 -1020.19 217103.78 53
Jupiter 60737 -12880 -4.15 839.93 1217.17 557078.35 8
Saturn -301530 -36762 6.11 -1948.89 -1591.05 513052.95 3
Uranus 152025 -19150 -2.09 1312.56 1681.40 246547.79 1
Neptune -125196 2514 -3.64 -844.43 -151.25 786449.21 0
Pluto -76912 6465 11.07 -132.25 -37.33 522747.90 0
*The inclination, longitude of perihelion, longitude of ascending node,
and mean longitude are measured in arcseconds per century (1

= 3600
arcseconds). The data for the earth are for the Earth-Moon barycentre.
25