### Notes from Under Sky

#### How to figure out a long-period comet's orbit based on just three observations

Ever wonder how astronomers can figure out a comet's orbit in three dimensions so quickly, even though you can't necessarily see how far the comet is at first glance? It turns out that the mathematics behind it are elementary, but complex. By that I mean that nothing beyond the first weeks of trigonometry is required, but the derivation is not all that easy to follow.

We will begin by assuming that our observations can be reckoned as if they were taken from the sun. As long as the comet is still far away, this will be a harmless enough assumption for us to obtain an approximate orbit.

We will also assume that the comet is approaching the sun on a parabolic orbit. Technically, Kepler's first law of planetary motion states that all planets and gravitationally bound comets travel in ellipses, but the long-period comets that venture into the inner solar system are so eccentric (e very close to 1) that they are for practical purposes parabolic (e equal to 1). In fact, some of them are perturbed sufficiently that they end up with a hyperbolic orbit (where e is greater than 1).

Kepler's second law of planetary motion states that a planet or comet sweeps our equal areas in equal times. Another way of saying the same thing is to say that the comet's angular velocity is inversely proportional to the square of its distance from the sun. For example, if it is twice as far, it must cover one-fourth the angle per unit of time, in order to sweep out an equal area. In algebraic terms, we can write this as

w (r) = K / r2

where r is the comet's distance from the sun, w is the comet's angular velocity, and K is some constant which we have to determine. Now, consider our three observations. They will reveal that the comet is not moving at a constant angular velocity, and this is the key to the enterprise. If it were travelling at a constant angular velocity, then it would have to be orbiting the sun in a circle, and that is counter to our original assumption. Instead, it will gradually be increasing its angular velocity, and the rate at which it does this will indicate how fast the comet is approaching, and determine the parameters of its orbit.

Suppose that the first observation O1 is taken at time t1 = 0 and at position angle ß1 = 0, and the second observation O2 is taken at time t2 and at position angle ß2. We can then obtain a derived observation D1.5 at time t1.5 at position angle ß1.5, where

t1.5 = t2 / 2
ß1.5 = ß2 / 2

In other words, we pretend that we saw the comet at a place exactly halfway between O1 and O2, at a time exactly halfway between as well. Similarly, if we say that the third observation O3 is taken at time t3 at position angle ß3, we can obtain a second derived observation D2.5 at time s2.5 at position angle b2.5, where

t2.5 = (t2 + t3) / 2
ß2.5 = (ß2 + ß3) / 2

We note that the average angular velocity between the first two observations is

w1.5 = ß2 / t2

and the the average angular velocity between the second and third observations is

w2.5 = (ß3 - ß2) / (t3 - t2)

and finally we also observe that the angular separation between the two derived observations is just

ßD = ß3 / 2

If we denote by r1.5 and r2.5 the distances of the comet from the sun at the moments of these two derived observations, then taking into account our first equation, we can write

 w2.5 / w1.5 = (K / r22.5) / (K / r21.5) = (r1.5 / r2.5)2

Consider now the triangle made by the sun S at one corner, and our two derived observations D1.5 and D2.5 at the others. If we denote by d the linear distance between the two derived observations, the law of cosines give us

d2 = r21.5 + r22.5 - 2r1.5r2.5 cos ßD

Dividing both sides of this equation by r22.5, we get

 (d / r2.5)2 = 1 + (r1.5 / r2.5)2 - 2 (r1.5 / r2.5) cos ßD = 1 + u - 2 sqrt (u) cos ßD

where u = w2.5 / w1.5. Now, the component of the motion along D1.5 D2.5 that is "sideways" is given by ßDr2.5. So, the angle which the segment D1.5D2.5 makes with the segment SD1.5 is

 v = sin–1 (ßDr2.5 / d) = sin–1 (ßD / sqrt (1 + u - 2 sqrt (u) cos ßD) )

From simple algebra, we know that the angle which SO2 makes with D1.5D2.5 is

vadj = v + ß2 / 2

Finally, the segment D1.5D2.5 is approximately parallel to the parabola at O2, and using the geometry of the parabola, we can determine that the perihelion point of the comet's orbit is located at position angle

 ßq = ß2 + pi - 2vadj = pi - 2v

Now, we must determine the actual perihelion distance q. We know that the perihelion is separated from the observation O1 by an angle of pi - 2v. We can find the actual distance r1 at O1 by solving the simultaneous equations

x = y tan 2v
y = (x2 / 4) - 1

This can be solved using the quadratic equation and some simple trigonometric identities to yield

y = (2 + 2 sec 2v) / tan2 2v

The distance proportion p is the ratio between the distance at O1 and the perihelion distance q, and is given by

 p = sqrt (x2 + y2) = y sec 2v

at the proper intersection solution point defined by the line and the parabola. Again, using Kepler's second law, we know that

 wq = K / q2 = K / (r1 / p)2 = p2 (K / r21) = p2 w1

where w1 is the derived angular velocity at O1

w1 = [t3 (ß2 / t2) - t2 (ß3 / t3) ] / (t3 - t2)

By equating observed centripetal acceleration with the predicted acceleration due to gravity, we can write

qw2q = GM / q2

from which we can easily derive

q = cbrt (GM / w2q)

where G is the universal gravitation constant, M is the mass of the sun, and cbrt is the cube root function. The orbit is now completely determined, with the perihelion at distance q, coplanar with the three observations, located at the angle ßq.

### Correcting for the Earth's Offset

Here is an approximate but simple way of accounting for the earth's offset in computing the orbit. Using the parabolic orbit determined above, find the predicted location of the comet in three-dimensional space at the three observation points, O1, O2, and O3.

Next, determine the earth's position relative to the sun at the time of the observations, and imagine an anti-earth E' such that the sun is precisely between the earth and E'. Compute the new angles ß'1 = 0, ß'2 = O1E'O2, and ß'3 = O2E'O3. Redo the orbit determination with the new angle values, and remember to reckon the orbital axis to the new ß'1. The new orbit should be approximately corrected for the earth's offset.