Consider two gravitating bodies far from any other gravitating bodies or
causes of other forces (but let's *not* assume either of them is
negligible compared to the other). How do they move ?

Chose a frame of reference (for some neighbourhood of their entwined space-time trajectories) in which the sum of their momenta is entirely time-wards; i.e. their spatial momenta are equal and opposite.

Let the total mass of the bodies be M; let their several masses be a.M and
A.M respectively; a+A = 1. When they are a distance R apart, sub-divide the
line between them in the proportions a.R at A.M's end and A.R at a.M's end;
the point of sub-division is known as their centre of gravity

. Because
the total spatial momentum is zero, yet equal to M times the velocity of the
centre of gravity, the centre of gravity isn't moving; so take it as origin
for non-spinning spatial co-ordinates in conjunction with the time-like
co-ordinate choice which made the total (spatial) momentum zero. Then the
centre of gravity's trajectory in space-time is the time axis of the
co-ordinate system. I shall presume that M, A and a are constant; R, and the
direction of the line between the masses, will vary with time, t.

Let w be the result of
dividing a.M's velocity by A.R (i.e. a.M's distance from the origin); because
the centre of gravity (our spatial origin) always lies on the line from a.M to
A.M, dividing it in the same proportions, dividing A.M's velocity by a.R
(A.M's distance from the origin) must give −w, the same but negated
(because we measured both radii as positive, despite their being in opposite
directions; while the velocities encoded opposite direction by
negation). This fractional velocity

, w, encodes both the fractional
rate of change of R, as its radial component, and the angular velocity of the
bodies about their common centre, as a circumferential component. Its
location on a diagram should be at unit

radius along the line from the
origin to a.M, extended past a.M if the unit

exceeds A.R; crucially,
it's on a.M's side of the origin; the corresponding vector to display on A.M's
side of the unit

circle would be −w.

Since a.M has velocity A.R.w, its momentum is a.A.M.R.w, while that of A.M will be −a.A.M.R.w. As a, A and M are constant, the rate of change of a.M's momentum is thus a.A.M.d(R.w)/dt; in accord with Newton's second law (for each action there is an equal and opposite reaction), that of A.M is the same but negated.

Newton tells us that the rate of change, a.A.M.d(R.w)/dt, of a.M's momentum is towards the origin with magnitude G.A.a.M.M/R/R, as long as the only force involved is their mutual gravitation (and neither R.w nor √(G.M/R) is near enough to light's speed to oblige us to worry about relativistic effects). The bodies' gravitational potential energies, each in the field of the other, are then −G.M.M/R times a.A.A and times A.a.a for a sum of −a.A.G.M.M/R.

Space-time has a metric, g; it's constant

and symmetric, so g(x,y)
= g(y,x), and I'll try not to be distracted by the fact that our co-ordinates
have been chosen to diagonalise it (with one component of one sign, three of
the opposite sign and all off-diagonals zero). Since w is the spatial
components of either body's velocity divided by the body's radius, it's
entirely spatial so g(w,w) only sees the spatial components of g, which I'll
presume to be positive-definite (making the time component negative).

Now, a.M's position vector divided by A.R is a unit vector, u. The velocity of a.M is A.R.w, yet is equally A.d(R.u)/dt = A.(u.dR/dt + R.du/dt), whence R.w = u.dR/dt + R.du/dt; so, in particular, w is in the plane spanned by u and du/dt. Now

- g(u,u) = 1 so 0 = d(g(u,u))/dt = 2.g(u,du/dt) whence g(u,du/dt) = 0 so
- g(R.u,w) = R.g(u,w) = g(u,R.w) = g(u,u).dR/dt = dR/dt, whence
- du/dt = w −u.dR/dt/R = w −u.g(u,w)

while a.M's acceleration divided by A, or A.M's acceleration divided by a, is - d(R.w)/dt = R.dw/dt + w.dR/dt = R.(dw/dt +w.g(u,w))

and Newton tells us (see above) that this last is −u.G.M/R/R; i.e.

- −u.G.M/R/R = d(R.w)/dt = R.(dw/dt +w.g(u,w)), so
- 0 = dw/dt +w.g(u,w) +u.G.M/R/R/R.

Introduce a unit vector, n, in the direction of du/dt, which is
orthogonal to u: this changes direction in the same manner

(subject to
the rotation which maps u to n) as u with dn/dt = −u.g(n,du/dt) and
du/dt = n.g(n,du/dt). We already know du/dt = w−u.g(u,w) and g(n,u) =
0, so we can infer g(n,du/dt) = g(n,w), whence

- du/dt = n.g(n,w) and
- dn/dt = −u.g(n,w).

Our force-law is central, i.e. g(n,w.g(u,w)+dw/dt) = −g(n,u).G.M/R/R/R = 0 as g(n,u) = 0, so

- g(n,w).g(u,w) +g(n,dw/dt) = 0 and
- g(n,dw/dt) = −g(n,w).g(u,w) = −g(u.g(n,w),w) = g(dn/dt,w), making
- d(g(n,w))/dt = −2.g(n,w).g(u,w) = −2.g(n,w).dR/dt/R, so
- d(R.R.g(n,w))/dt = 2.R.g(n,w).dR/dt +R.R.d(g(n,w))/dt = 0

making J = g(R.n,R.w) constant: multiplied by suitable masses, this gives the angular momenta of the bodies; a.M has angular momentum a.M.g(A.R.n,A.R.w) = a.A.A.M.J while A.M's is A.a.a.M.J and, as a+A = 1, the sum is S = a.A.M.J. As 1/(a.A.M) = (a.M+A.M)/(a.M)/(A.M) = 1/(a.M) +1/(A.M), J is the total angular momentum multiplied by the sum of inverses of the masses. [Thus S = J/(1/(a.M) +1/(A.M)), so 2.S/J = 2/(1/(a.M) +1/(A.M)), which is called the harmonic mean of the two masses.] J has dimensions area/time and, thanks to our choice of n's direction, J is non-negative.

In particular, g(n,w) = J/(R.R) can never change sign. Note that, for non-zero R, J = 0 iff g(n,w) = 0 iff w is parallel to u iff du/dt is zero iff u and n are constant. In any case,

- u.J/R/R = u.g(n,w) = −dn/dt, giving
- d(J.w.R)/dt = −u.G.M.J/R/R = d(G.M.n)/dt, whence
- J.w.R −G.M.n is constant

Let Q be the dimensionless vector w.R.J/G/M −n, which is constant; we then have

- 1 +g(n,Q) = g(n, n+Q) = g(n, J.w.R/G/M) = J.g(n.R, w.R)/R/G/M = (J.J/G/M)/R,

so introduce the length L = J.J/G/M, the semi-latus rectum

(which is all Latin to me) of the conic section

(the curve of
intersection between a cone and a plane) our trajectory will turn out to be
following. Our trajectory can now be written as:

- L/R = 1+g(n,Q), with
- Q = J.w.R/G/M −n = w.R.L/J −n constant; only w, R and n vary.

Discussion must, given that J is non-negative (by construction), now break into two cases: J = 0 and J > 0.

The following discussion involves the trigonometric functions sine and
cosine; following habits from elsewhere in
this web site, I'll use the Capitalised names Sin and Cos for the functions
which map angles (in units of radians, turns, degrees or whatever) to scalars;
and their lower-case equivalents, sin = (: Sin(x.radian) ←x :) and cos = (:
Cos(x.radian) ←x :), for the closely related pure-scalar functions which
satisfy cos' = −sin, sin' = cos and, in the complex plane, exp(i.z) =
cos(z) +i.sin(z). The J=0 case will also exploit sinh(z) = (exp(z)
−exp(−z))/2 and cosh(z) = (exp(z) +exp(−z))/2, making
cosh(i.z) = cos(z) and sinh(i.z) = i.sin(z), with sin^{2}
+cos^{2} = 1 = cosh^{2} −sinh^{2}, sinh' = cosh
and cosh' = sinh.

For the special case of circular orbit,

- R is constant, so 0 = dJ/dt = d(R.R.g(n,w))/dt implies j = g(n,w) is constant,
- 0 = dR/dt = R.g(u,w) and R > 0, so g(u,w) = 0,
- w = n.g(n,w) + u.g(u,w) = j.n,
- du/dt = w = j.n, dn/dt = −j.u, ddu/dt/dt = −j.j.u
- −u.G.M/R/R/R = dw/dt + w.g(u,w) = j.dn/dt + 0 = −j.j.u
- j.j = G.M/R/R/R

so, supposing the A.M body has radius C.R and density B, while the a.M body has radius c.R and density b, so that A.M = (4.π/3).B.C.C.C.R.R.R, a.M = (4.π/3).b.c.c.c.R.R.R, we find

- j.j/G
- = M/R/R/R
- = (4.π/3).B.C.C.C/A
- = (4.π/3).b.c.c.c/a

enabling us to infer the densities – b from j, c and a, B from
j, C and A. If we are observing the mutually-orbiting system from afar,
possibly without knowledge of our distance from the system (hence without
knowledge of the absolute value of R), we'll still be able to determine j
(because 2.π/j will be the period of the mutual orbit, given ddu/dt/dt =
−j.j.u), while the remaining quantities (A, a, C and c)
are *ratios* between lengths, which we can infer from the angles we
see, so insensitive to our ignorance of R. Thus we can determine
the *densities* of mutually orbiting bodies from afar, at least when
the trajectory is a circular orbit.

Furthermore, though one of the bodies may be too small to distinguish its
radius from zero, if the other is large enough that we can see what fraction
of their separation *its* radius is, then we can still
determine *its* density, even though we cannot determine that of its
partner. Thus Newton's contemporaries presumably had enough data to determine
Jupiter's average density, if not especially accurately; at 1.33 g/cc it's
only slightly denser than water, on average – less than a quarter of
Earth's average density. This should at least have given some clue to the
effect that Jupiter is a lot bigger than Earth. As moons of other planets
showed up, the densities of Mars (3.94 g/cc, Earth / 1.4), Saturn (0.7 g/cc,
Earth / 7.9), Uranus (1.3 g/cc, Earth / 4.2) and Neptune (1.76 g/cc, Earth /
3.1) must have come to light … oh, and Pluto (2.0 g/cc, Earth / 2.8),
back when it was counted as a planet, but only relatively recently did we
realize it was a binary system.

When J is zero, hence so is L and the above description of the trajectory becomes fatuous, 0 = 0, we still have Q = w.R.J/G/M −n = −n constant, so conclude that n and u are constant: the bodies are confined to a fixed line through the origin, with no angular velocity; their momenta are collinear and there is no mutual angular momentum. Now, with u and n constant, w is simply u.g(u,w) = u.dR/dt/R. Our dynamical equation reduces to:

- 0 = ddR/dt/dt +G.M/R/R,
- from which we may infer:
- d(dR/dt.dR/dt)/dt
- = 2.dR/dt.ddR/dt/dt
- = −2.G.M.dR/dt/R/R
- = d(2.G.M/R)/dt,

- so (dR/dt)
^{2}−2.G.M/R is constant.

Multiplying this last by a.A.M/2 yields the total energy; dividing it by 2.G.M yields a constant with dimensions 1/length; call this K; it might be zero, positive or negative. We obtain ±dR/dt = √(2.G.M.(K+1/R)) so integrating dt = ±dR/√(2.G.M.(K+1/R)) will complete our solution.

At this point we can examine energies: the kinetic energy of a.M is
(1/2).a.M.(d(A.R)/dt)^{2} = (1/2).a.A.A.M.(dR/dt)^{2} and that
of A.M is (1/2).A.a.a.M.(dR/dt)^{2} for a sum of
(1/2).a.A.M.(dR/dt)^{2} = (1/2).a.A.M.2.G.M.(K+1/R); as observed
above, the total gravitational potential energy is −a.A.G.M.M/R; adding
this to the kinetic energy we get a.A.M.((dR/dt)^{2} −2.G.M/R)/2
= G.(a.M).(A.M).K which is constant; it's the gravitational potential energy
of either body in the field of the other at separation 1/K. So the (constant)
total energy for the collision course is:

- E = a.A.G.M.M.K = a.A.M.((1/2).(dR/dt)
^{2}−G.M/R)

Now, to solve (dR/dt)^{2} = 2.G.M.(K+1/R), consider three
branches, depending on the sign of K:

- K < 0
with K+1/R > 0, hence R < −1/K, so introduce an angular parameter s with −K.R = (1 −cos(s))/2 = (sin(s/2))

^{2}so that- (dR/dt)
^{2}/ (−2.G.M.K) = −(1+1/K/R) = 1/(sin(s/2))^{2}−1 = (cos(s/2)/sin(s/2))^{2}, - so tan(s/2).dR/dt = ±√(−2.G.M.K), which is constant, while
- −K.dR/dt = sin(s).ds/dt/2 = sin(s/2).cos(s/2).ds/dt, so
- ±√(−2.G.M.K) = −tan(s/2).sin(s/2).cos(s/2).ds/dt/K =
−sin(s/2)
^{2}.ds/dt/K = R.ds/dt, giving - ±K.√(−2.G.M.K).dt = K.R.ds = −ds.sin(s/2)
^{2}= −ds.(1−cos(s))/2 = d(sin(s)−s)/2, so - s −sin(s) ±2.K.t.√(−2.G.M.K) is constant

then, choosing t = 0 at a moment when s = 0, hence R = 0, the given constant is zero so we have sin(s) −s = ±2.K.t.√(−2.G.M.K) with R = (cos(s) −1)/2/K describing our solution.

Now, R = 0 happens when cos(s) is 1, i.e. when s.radian is a whole number of turns. Since sin(s) −s changes by 2.π between such moments, at each of which sin(s) is 0, t changes by −π/K/√(−2.K.M.G), giving this as the time it takes the bodies to fly apart after one collision, then return to collide once more. Of course, in practice, the bodies shall interact in some non-gravitational manner when they collide, so you'll only see part of one cycle of this process; the analysis ignores this, effectively pretending that they either pass through one another, unaffected, or bounce off one another perfectly elastically.

Writing H for −1/K, the maximum value R attains, we get a total energy of a.A.G.M.M/H and sin(s) −s = ±2.(t/H).√(2.G.M/H) with a period of H.π.√(H/2/G/M).

- (dR/dt)
- K = 0
- (dR/dt)
^{2}= 2.G.M/Rso

- d(√(R.R.R))/dt = (3/2).√(R).dR/dt = ±3.√(G.M/2), i.e.
- ±3.t.√(G.M/2) −√(R.R.R) is constant

By taking t = 0 when R = 0, at the moment of collision, this becomes

- t = ±√(2.R.R.R/G/M)/3, so
- 9.G.M.t.t = 2.R.R.R

for which the bodies fall together, pass through or bounce off one another (the given solution doesn't tell you which, or care) at t=0 then fly apart in the time-reverse (possibly

centrally inverted

spatially) of their earlier collision trajectories. Both long before and long after the collision, the bodies are very far apart and essentially stationary; in the limit as long tends to infinitely, they are infinitely far apart and actually stationary. - K > 0
By analogy with the negative case, consider a parameter s for which K.R = sinh(s/2)

^{2}= (cosh(s) −1)/2 so that- K.dR/dt = sinh(s).ds/dt/2 = sinh(s/2).cosh(s/2).ds/dt, yet
- (dR/dt)
^{2}/ (2.G.M.K) = 1+1/sinh(s/2)^{2}= (cosh(s/2)/sinh(s/2))^{2}whence - ds/dt = K.dR/dt/sinh(s/2)/cosh(s/2) =
±K.√(2.G.M.K)/sinh(s/2)
^{2}, or - ±√(2.G.M.K).K.dt = ds.sinh(s/2)
^{2}= ds.(cosh(s) −1)/2 = d(sinh(s) −s)/2, so - s −sinh(s) ±2.K.t.√(2.G.M.K) is constant

again – choosing t = 0 when s = 0, hence R = 0 – we can make the constant zero, giving sinh(s) −s = ±2.K.t.√(2.G.M.K) while R = (cosh(s) −1)/2/K. For large values of s, cosh(s) and sinh(s) are roughly equal and s (despite being large) is negligible compared to them, yielding R = t.√(2.G.M.K). This solution describes bodies intially falling towards one another at relative speed √(2.G.M.K) and ultimately flying apart at the same relative speed.

Writing V = √(2.G.M.K) for this initial and final relative speed, yielding 2.K = V.V/(G.M), we find that the solution is sinh(s) −s = ±2.V.V.V.t/G/M.

The similarities between the K<0 and K>0 cases arise from the broken duality between [sin, cos] and [sinh, cosh]. Still, enough of the unrealistic special case with exactly no angular momentum …

If J = √(L.G.M) is non-zero (in which case our choice of n makes it positive), we can divide G.M.Q = J.w.R −G.M.n by J to obtain a constant velocity V = G.M.Q/J = w.R −n.G.M/J = w.R −n.J/L, giving Q = J.V/G/M = V.L/J, so expressing our trajectory as

- L/R = 1+g(L.n,V)/J, with
- w.R = V +n.J/L as the tangent to the trajectory.

If ever u is parallel to V, g(n,V) will be zero, making R = L; i.e. L is the radius at which the trajectory crosses the line, though the origin, in the direction of V. [V is perpendicular to the focal axis of the conic that R.u traces out; V is of order 36 km/s for Earth's orbit around the Sun; one thousandth of G.M/J]. In any case, we now have

- dR/dt = g(u,R.w) = g(u,V).

and R.w, our bodies' normalised

velocity, consists of a
constant, V, plus a rotating cirumferential

component of constant
magnitude, J/L. Any resemblance to the epicycles

of the ancient Greeks
can be construed as evidence they were less wrong that they're sometimes
accused of ;^)

In the following analysis, I'll make use of a constant speed, v, which is
V's magnitude, i.e. v is a positive scalar with v.v = g(V,V). Likewise, let e
be the magnitude of our dimensionless vector Q = V.L/J, i.e. e = v.L/J,
yielding v.v = e.e.G.M/L. This e is known as the eccentricity

of the
trajectory. Since n was so chosen as to make J non-negative, and we're dealing
with the case of non-zero J, e cannot be negative.

The kinetic energy of a.M is a.M.g(A.w.R,A.w.R)/2, that of A.M is A.M.g(a.w.R,a.w.R)/2, for a sum of a.A.M.g(w.R,w.R)/2, so

- w.R = V +n.G.M/J = V +n.J/L makes
- g(w.R,w.R) = g(V,V) + 2.g(n,V).J/L +J.J/L/L, and
- g(n,V) = J.(1/R −1/L), so
- g(w.R,w.R) = v.v +J.J/L/L +2.J.J.(1/R −1/L)/L = v.v +J.J.(2/R−1/L)/L, and J.J/L = G.M, so the kinetic energy is
- a.A.M.v.v/2 +a.A.G.M.M.(2/R −1/L)/2, or
- a.A.M.(v.v/2 −G.M/L/2 +G.M/R) and v.v/e/e = G.M/L, so this is
- (e.e −1 +2.L/R).a.A.G.M.M/L/2 = a.A.M.((e.e−1).G.M/L/2 +G.M/R)

in which only the last term varies; and its value, a.A.G.M.M/R, is just −1 times the total potential energy, so the sum of potential and kinetic energy is

- E = (e.e−1).G.a.M.A.M/L/2, which is constant.

Thus the total energy is positive iff e.e > 1, i.e. the eccentricity is bigger than 1. When R = L, the kinetic energy is a.A.(e.e+1).G.M.M/L/2, which is −(e.e+1)/2 times the potential energy.

The rotating frame of reference induced by n and u sees w.R = V +n.J/L as
the constant

circumferential velocity n.J/L plus a rotating

offset, V, which adds a sinusoidally varying component g(n,V) to the angular
velocity and contributes a complementary sinusoidally varying component g(u,V)
as the radial velocity; but note that the sinusoid's angular velocity is the
angular velocity to which g(V,n) is contributing, so the sinusoid's phase
doesn't vary linearly with time.

Introduce an angle θ specified by: g(n,V) = v.Cos(θ), g(u,V) = v.Sin(θ) – i.e. chose (positive) y-axis parallel to V and look at the plane of the trajectory from whichever (spatial) side makes the angular velocity positive (anti-clockwise), to infer the x-axis (a quarter turn clockwise from the y-axis), then use the usual polar co-ordinate measured anti-clockwise with θ zero on the positive x-axis.

Since V = n.v.Cos(θ) +u.v.Sin(θ) is constant, so are both its magnitude, v, and its direction, V/v = n.Cos(θ) +u.Sin(θ). Our trajectory is now expressed as

- L/R = 1 +g(n,V).L/J = 1 +(v.L/J).Cos(θ) = 1 + e.Cos(θ)

with dR/dt = g(u,V) = v.Sin(θ). We can now re-write the total kinetic energy as

- (e.e −1 +2.L/R).G.a.M.A.M/L/2 = (e.e +1 +2.e.Cos(θ)).G.a.M.A.M/L/2

If we write X = R.Cos(θ), Y = R.Sin(θ), we now have L = R + e.X whence

- X.X + Y.Y = R.R = (L−e.X)
^{2}= L.L −2.e.L.X +e.e.X.X, whence - Y.Y = L.L −2.e.L.X +(e.e −1).X.X = (e.e−1).(X
−e.L/(e.e−1))
^{2}+L.L.(1 −e.e/(e.e−1)), so - Y.Y.(1−e.e) +(X.(1−e.e) +e.L)
^{2}= L.L,

which describes an ellipse for e.e<1 and a hyperbola for e.e>1. In either of these cases, substituting X=0 yields Y.Y = L.L so Y = ±L (as already anticipated); while substituting Y=0 yields X = L.(±1 −e)/(1−e.e) = L/(e±1). For the case e.e = 1, we get Y.Y = L.L −2.L.X, or 2.L.X = L.L − Y.Y, which is a parabola.

Our dynamical equation now reduces to θ's time-dependence,

- dθ/dt/radian = g(n,w) = J/R/R = (J/L/L).(1+e.Cos(θ))
^{2}

so time is the integral of L.L/J/(1+e.Cos(θ))^{2} with
respect to θ/radian, plus some arbitrary offset – albeit this
integral is no picnic ! Again, we break into three cases, according as
the total energy is positive, zero or negative; which is to say, according as
e is greater than, equal to or less than 1. The X,Y analysis above revealed
that these cases give hyperbolic, parabolic and elliptical orbits,
respectively.

- e > 1, hyperbola
From L/R = 1 + e.Cos(θ) with L, R and e positive, we can infer Cos(θ) > −1/e. This gives us a wedge-shaped portion of the plane from which the (R,θ) representation of our bodies' trajectories are excluded; the apex of the wedge is our origin and its two bounding rays point outwards in the directions with Cos(θ) = −1/e; these are in the −ve X half-plane, one in its +ve Y quadrant, the other in its −ve Y quadrant. The portion between these, straddling the −ve X-axis, is the wedge-shaped region from which the orbit is excluded.

Our hyperbola is

centred on

the point X = e.L/(e.e−1), Y = 0, i.e. R = e.L/(e.e−1), θ = 0. If we translate the exclusion-wedge in the X-direction through e.L/(e.e−1), the hyperbola lies entirely within the translated wedge, though outside the un-translated wedge; it is confined to the >-shaped trench between the two wedges' boundaries. In its (R,θ) representation, our bodies' trajectories come from far out along the −ve Y branch of the translated wedge's edge; move steadily towards the apex of the > and away from this leading edge, speeding up as they move inwards; cross the line θ = −turn/4 at R=L; swing round through the angle in the trench, crossing θ=zero at R=L/(1+e) < L/2; then escape in like manner, though time-reversed and reflected in the X-axis, to ultimately approach the translated wedge's edge's +ve Y branch.The relative speed of the two particles while very far apart is just v.√(1−1/e/e) = v.√(e.e−1)/e = J.√(e.e−1)/L. The maximum relative speed is J.(e+1)/L, which happens when θ=zero, R=L/(1+e). The ratio of these two speeds, dividing the former by the latter, is √((e−1)/(e+1)); it tends to 1 for large e but tends to zero as e approaches 1. When far from each other, the bodies' velocities are almost exactly either towards or away from one another; the magnitude increases with e; near e = 1, when the internal angle of the wedges is tiny (and the outer vertex is far from the origin, to keep the Y-wards width still L), they are barely moving (relative to one another).

- e = 1, parabola
This will swing round the origin in very similar manner to the hyperbola, from far away in the −X direction, through the −ve X, −ve Y quadrant, crossing into the +ve X, −ve Y quadrant at R = L on the −ve Y-axis, swinging through this quadrant to pass through the X-axis at X = L/2 and a top speed of 2.J/L before followinng a time-reversed mirror image of its infalling trajectory. Its speed tends to zero at large radius. This can be thought of as the limiting case of the hyperbola as e tends to 1 from above, albeit the wedge has become the whole plane (its outer apex is infinitely far away and the only place excluded is the negative X-axis); θ's value is arbitrarily close above −turn/2 for large negative t and tends to +turn/2 from below as t tends to infinity. It can also be thought of as the limiting case of an ellipse (see below) with the distant end of the major axis

at infinity

.We can exploit 1+Cos(θ) = 2.Cos(θ/2)

^{2}to obtain J.dt/L/L = 4.dθ/Cos(θ/2)^{4}/radian. Substituting z = Tan(θ/2), we obtain 2.radian.dz/dθ = Sec(θ/2)^{2}= 1+z.z whence- J.dt/L/L/8 = dz.(1+z.z) = d(z +z.z.z/3)

whence z.(1 +z.z/3) = J.t/L/L/8 = 3.t.(G.M)

^{2}/(2.J)^{3}plus a constant which we can make zero by taking t = 0 when z = 0, i.e. when θ is zero. We start with z large and negative (θ just increasing from −turn/2) and end with z large and positive (θ tending towards turn/2). For large t, Tan(θ/2) = z is then approximately half the cube root of 3.J.t/L/L or, equally, approximately 1/2/J times the cube root of 3.t.(G.M)^{2}.- e < 1, ellipse
with (1−e.e).Y.Y +((1−e.e).X +e.L)

^{2}= L.L, R.u's orbit encloses an area π.L.L/(1−e.e)/√(1−e.e). The rate at which the R.u ray sweeps out area is g(R.n,R.w)/2 = J/2, which is constant, and this ray sweeps out the entire area of the ellipse in each orbit, so the period of the orbit is- period = 2.π.L.L/J/(1−e.e)/√(1−e.e)

Substituting J = √(L.G.M) into this gives a formula in which L/(1−e.e) appears with power 3/2; and L/(1−e.e) is just the semi-major axis of our ellipse, since L/(1−e) +L/(1+e) = 2.L/(1−e.e). Thus, writing H for the semi-major axis, we get period 2.π.H.√(H/G/M) regardless of eccentricity. Furthermore, E = (e.e−1).G.a.M.A.M/L/2 = −a.A.G.M.M/H/2, so (for a given pair of bodies) both the total energy and the period depend only on the semi-major axis of the orbit. The product of energy and period, −π.a.A.M.√(G.M.H) has the same dimensions as angular momentum, S; so introduce

- f = −E.period/S
- = 2.π.L.L/J/(1−e.e)/√(1−e.e) . a.A.(1−e.e).G.M.M/L/2 / (a.A.M.J)
- = π.L.G.M/J/J/√(1−e.e)
- = π/√(1−e.e), since L.G.M = J.J.

The orbit is circular when e is 0, making f = π so −E.period = π.S with L equal to the semi-major axis (which, in this case, is the radius).

Since R, u and n are expressed in terms of θ, relating θ
to time (i.e. integrating up the equation dt.J/L/L =
dθ/radian/(1+e.Cos(θ))^{2}) is all that remains for a
full, general solution of Newton's two-body problem.

- dt.J/L/L
- = dθ/radian/(1+e.Cos(θ))
^{2} Substitute u = k.Tan(θ/2), for some constant k to be chosen in due course: this gives Cos(θ) = (k.k −u.u)/(k.k +u.u), du = k.Sec

^{2}(θ/2).dθ/2/radian, so dθ/radian = 2.du.Cos^{2}(θ)/k and we get- = 2.du.(k.k −u.u)
^{2}/(k.k +u.u)^{2}/k/(1+e.(k.k −u.u)/(k.k +u.u))^{2} - = 2.du.(k.k −u.u)
^{2}/k/((k.k +u.u) +e.(k.k −u.u))^{2} - = (2.du/k).((k.k −u.u)/(k.k.(1+e)
+(1 −e).u.u))
^{2} - = (2.du/k).((k.k −u.u)/(1−e)/(k.k.(1+e)/(1−e) +u.u))
^{2} - = (2.du/k).(( (k.k.(1 +(1+e)/(1−e)) −(k.k.(1+e)/(1−e)
+u.u))/(k.k.(1+e)/(1−e) +u.u) )/(1−e))
^{2} - = (2.du/k/(1−e)
^{2}).(2.k.k/(k.k.(1+e) +(1−e).u.u) −1)^{2} - = (2.du/k/(1−e)
^{2}).(4.k^{4}/(k.k.(1+e) +(1−e).u.u)^{2}−4.k.k/(k.k.(1+e) +(1−e).u.u) +1) whence, picking t = 0 at some moment when θ is 0 or a half turn, hence u is zero,

- k.t.J.(1−e).(1−e)/L/L/2
- = integral(: 4.du.k
^{4}/(k.k.(1+e) +(1−e).u.u)^{2}−4.du.k.k/(k.k.(1+e) +(1−e).u.u) +du ←u :{reals between 0 and k.Tan(θ/2)}) The final du is easy; in the other two terms, substitute k.Tan(v).√(1+e) = u.√(1−e), i.e. v = arcTan(Tan(θ/2).√((1−e)/(1+e))), so that du = k.Sec

^{2}(v).(dv/radian).√((1+e)/(1−e)) while k.k.(1+e) +(1−e).u.u = k.k.(1+e).(1 +Tan^{2}(v)) = k.k.(1+e).Sec^{2}(v), to get:- = integral(:
4.k.dv/radian/√(1−e.e)/(1+e)/Sec
^{2}(v) −4.k.dv/radian/√(1−e.e) ←v :{angles between 0 and arcTan(Tan(θ/2).√((1−e)/(1+e)))}) +k.Tan(θ/2) - = 4.k.integral(: Cos
^{2}(v).dv/radian ←v :{angles between 0 and arcTan(Tan(θ/2).√((1−e)/(1+e)))})/√(1−e.e)/(1+e) −4.k.arcTan(Tan(θ/2).√((1−e)/(1+e)))/radian/√(1−e.e) +k.Tan(θ/2) in which Cos

^{2}(v).dv/radian = (1 −Sin^{2}(v)).dv/radian = dv/radian +d(Cos(v)).Sin(v) = dv/radian +d(Cos(v).Sin(v)) −Cos(v).d(Sin(v)) = dv/radian +d(Cos(v).Sin(v)) −Cos^{2}(v).v/radian, the last term in which is just what we started with negated, hence 2.Cos^{2}(v).dv/radian = dv/radian +d(Cos(v).Sin(v)), so:- = 2.k.integral(: dv/radian +d(Cos(v).Sin(v)) ←v :{angles between 0 and arcTan(Tan(θ/2).√((1−e)/(1+e)))})/√(1−e.e)/(1+e) −4.k.arcTan(Tan(θ/2).√((1−e)/(1+e)))/radian/√(1−e.e) +k.Tan(θ/2)
- = 2.k.arcTan(Tan(θ/2).√((1−e)/(1+e)))/radian/√(1−e.e)/(1+e) +2.k.Sin(2.arcTan(Tan(θ/2).√((1−e)/(1+e))))/√(1−e.e)/(1+e) −4.k.arcTan(Tan(θ/2).√((1−e)/(1+e)))/radian/√(1−e.e) +k.Tan(θ/2)
- = 2.k.(1/(1+e) −2).arcTan(Tan(θ/2).√((1−e)/(1+e)))/radian/√(1−e.e) +2.k.Sin(2.arcTan(Tan(θ/2).√((1−e)/(1+e))))/√(1−e.e)/(1+e) +k.Tan(θ/2)
- = −2.k.(2.e +1).arcTan(Tan(θ/2).√((1−e)/(1+e)))/radian/√(1−e.e)/(e+1) +2.k.Sin(2.arcTan(Tan(θ/2).√((1−e)/(1+e))))/√(1−e.e)/(1+e) +k.Tan(θ/2)

which should at least let us compute t from θ, albeit going the other way is apt to be rather tricky !

Note that, since the energy (for a given pair of bodies) and the period only depend on the semi-major axis of the orbit, in each case as a monotonic function of H, we can infer that orbits which pass through a given point at equal speeds must necessarily have equal periods – and vice versa. (I first met this detail as a school-boy in the form of a question specifying that an object, instantaneously at rest in a gravitational field, is blown up by a bomb which causes every part of the object to get a share of the released energy, in proportion to its mass; which amounts to specifying equal speeds: all the fragments shall return to their start-point at a common later moment – ignoring any which run into something else along the way.)

The period of anything is an amount of time per cycle; so, properly, the
period of our elliptical orbit (J>0, e<1) should be construed as an
inverse rate of orbit

of
2.π.L.L/J/(1−e.e)/√(1−e.e) *per turn*, hence
equally as L.L/J/(1−e.e)/√(1−e.e) per radian. This turns f
into E/S/(rate of orbit) and makes 1/f equal to 2.(1−e.e).radian unless
we give S units that include a 1/angle term to cancel the angle term arising
here in f. The proper rôle and place of units of angle, in particular
in connection with angular

momentum, is an area that causes me some
doubts: notably, Planck's constant is an energy times period, like f.S, so
should really involve the unit of angle – with the result that
multiplying it by turn yields the quantity usually known as Planck's constant,
while multiplying it by radian yields Dirac's constant, which is a factor of
2.π smaller.

Note, also, that L.L/J/(1−e.e)/√(1−e.e)/radian = J.J.J/G/G/M/M/(1−e.e)/√(1−e.e) =

- L.L/J/(1−e.e)/√(1−e.e)
- = J.J.J/G/G/M/M/(1−e.e)/√(1−e.e)
- = (J/√(1−e.e))
^{3}/(G.M)^{2} - = (L/(1−e.e))
^{3/2}/(G.M)^{1/2} - = (H.H.H/G.M)
^{1/2}

which encourages a review of the above with attention to J/√(1 −e.e) or L/( 1−e.e) as the terms to which to pay attention.

In all of the above, I have used a non-rotating frame of reference to
describe the system; however, since the elliptical orbit has an over-all
average angular velocity, it makes sense to subtract off

this rotation
and describe the orbit from the point of view of
a rotating frame of reference. The rate of
orbit

(or inverse period) alluded to above is
J.radian.(1−e.e).√(1−e.e)/L/L, so this is the angular
velocity to use for our rotating frame. The resulting frame will be useful
for considering the motions of smaller bodies in the
neighbourhood of such a mutually orbiting pair.

Our natural point of departure from the above is the pair of equations

- J.radian = R.R.dθ/dt is constant
- L = R.(1 + e.Cos(θ)) is also constant, as is e

and we are given that J>0 and 0≤e<1. We are changing to a co-ordinate system in which R is unchanged but θ is replaced by

- φ = θ − t.J.radian.(1−e.e).√(1−e.e)/L/L

so that, introducing D = L.L/(1−e.e)/√(1−e.e) for the result of dividing our ellipse's area by π,

- dφ/dt/radian
- = J/R/R −J/D
- L/R
- = 1 + e.Cos(θ)
- = 1 + e.Cos(φ + t.radian.J/D)
- = 1 + e.(Cos(φ).cos(t.J/D) −Sin(φ).sin(t.J/D))
For circular orbit, e = 0, we have R = L and φ constant

I think the right approach will be to assume small e, derive the dynamical equations in the rotating frame, apply some suitable approximations and solve by perturbation – expect another epicycle ;^)

Written by Eddy.