Newton's two-body problem

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 ?

Framework – co-ordinates, quantification and invariants

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.

diagram showing co-ordinate system 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.

Dynamics

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

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

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

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

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,

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

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:

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

Digressions

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 sin2 +cos2 = 1 = cosh2 −sinh2, sinh' = cosh and cosh' = sinh.

For the special case of circular orbit,

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.

Collision Course (J = 0)

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:

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.

Energy (J=0)

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:

Detailed solutions (J=0)

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

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).

K = 0
(dR/dt)2 = 2.G.M/R

so

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

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

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 …

Scattering or orbiting with J > 0

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

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

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.

Energy (J>0)

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

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

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.

Detailed solution (J>0)

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

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

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

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,

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 hyperbolic trajectory

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

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

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.Sec2(θ/2).dθ/2/radian, so dθ/radian = 2.du.Cos2(θ)/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.k4/(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.k4/(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.Sec2(v).(dv/radian).√((1+e)/(1−e)) while k.k.(1+e) +(1−e).u.u = k.k.(1+e).(1 +Tan2(v)) = k.k.(1+e).Sec2(v), to get:

= integral(: 4.k.dv/radian/√(1−e.e)/(1+e)/Sec2(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(: Cos2(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 Cos2(v).dv/radian = (1 −Sin2(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)) −Cos2(v).v/radian, the last term in which is just what we started with negated, hence 2.Cos2(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.)

A little aside about units of angle

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.

Elliptical orbit minus its mean rotation (J>0, e<1)

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

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

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 ;^)

Valid CSSValid HTML 4.01 Written by Eddy.