Author Archives: growescience

Halley’s comet: numerical calculation of orbit

Reference: Carroll, Bradley W. & Ostlie, Dale A. (2007), An Introduction to Modern Astrophysics, 2nd Edition; Pearson Education – Chapter 2, Problem 2.15.

We can use Kepler’s laws to generate numerical solutions for orbits on a computer. Carroll & Ostlie provide a program called Orbit (in Fortran or C++) that produces textual output that can be input to a spreadsheet or graphing program. However, this is a bit primitive, so I decided to implement a similar routine in Maple, since Maple has built-in plotting features.

The idea is to simulate a complete orbit of a planet by calculating its distance {r} from the focus of the ellipse (essentially, the distance from the Sun) as a function of the angle {\theta} from perihelion. The equation of an ellipse is:

\displaystyle r=\frac{a\left(1-e^{2}\right)}{1+e\cos\theta} \ \ \ \ \ (1)

Given the planet’s period {P} we can calculate its semimajor axis {a} from Kepler’s third law:

\displaystyle P^{2}=\frac{4\pi^{2}}{GM}a^{3} \ \ \ \ \ (2)

so the equation 1 becomes

\displaystyle r=\left(\frac{GMP^{2}}{4\pi^{2}}\right)^{1/3}\frac{\left(1-e^{2}\right)}{1+e\cos\theta} \ \ \ \ \ (3)

To use this formula in a numerical solution, we need to know how much {\theta} changes for a given time increment {dt}. We can use Kepler’s second law in the form

\displaystyle \frac{dA}{dt}=\frac{L}{2\mu} \ \ \ \ \ (4)

which gives the rate at which area is swept out in terms of the constant total angular momentum {L} and reduced mass {\mu}. The angular momentum is

\displaystyle L=\mu\sqrt{GMa\left(1-e^{2}\right)} \ \ \ \ \ (5)

The area increment {dA} is given in terms of {d\theta} by

\displaystyle dA=\frac{r\; d\theta}{2\pi r}\pi r^{2}=\frac{1}{2}r^{2}d\theta \ \ \ \ \ (6)

so we get

\displaystyle d\theta=\frac{L}{\mu r^{2}}dt=\frac{\sqrt{GMa\left(1-e^{2}\right)}}{r^{2}} \ \ \ \ \ (7)

The numerical solution divides the period {P} up into {n} equal time intervals {dt}, and starts with {r} at perihelion and {\theta=0}. We then use the value of {r} to calculate {d\theta}, add this to the current value of {\theta} and calculate the next value of {r} from 3. We then use the new value of {r} to get the next {d\theta} and so on until we’ve covered a complete period.

The Maple code for doing this is:

G := 0.6673e-10;
AU := 0.14959787066e12;
M__sun := 0.19891e31;
yr := 31558145.0;
rad2deg := 180/Pi; 
secsYear := 365.25*(3600*24);
orbit := proc (M__strsun, a__AU, e, n) 
local M__star, a, P, dt, t, theta, LoM, r, i, dtheta; 
M__star := M__strsun*M__sun; 
a := a__AU*AU; 
P := sqrt(4*Pi^2*a^3/(G*M__star)); 
dt := P/(n-1); 
t := Array(0 .. n, (i) -> i*dt); 
r := Array(0 .. n, datatype = float); 
theta := Array(0 .. n, datatype = float); 
theta[0] := 0.; 
LoM := sqrt(G*M__star*a*(1-e^2)); 
for i from 0 to n-1 do 
  r[i] := a*(1-e^2)/(1+e*cos(theta[i])); 
  dtheta := LoM*dt/r[i]^2; 
  theta[i+1] := evalf(theta[i]+dtheta); 
  if 0 < i and 1.0 < r[i]/AU and r[i-1]/AU < 1.0 then 
    print("Passes 1 AU at t = ", evalf(t[i]/secsYear)) 
  end if 
end do; 
r[n] := a*(1-e^2)/(1+e*cos(theta[n])); 
theta := theta*rad2deg; 
r := r/AU; 
print(polarplot(r, theta, angularunit = degrees)); 
print(plot(t/secsYear, r, labels = ["t (years)", "r (AU)"])) 
end proc;

The code pretty much just implements the algorithm given above. The Maple procedure orbit on line 8 takes as its arguments the mass of the central star in solar masses, the semimajor axis of the planet in AU, the eccentricity {e} and the number of time steps {n}. It then converts these quantities into SI units using the conversion factors given at the start, creates arrays for the time {t}, the radius {r} and the angle {\theta}, calculates {L/\mu} (as LoM), and then enters a for loop to calculate {r} and {d\theta} for each time increment. The ‘if’ statement finds the time at which the planet crosses from {r<1\mbox{ AU}} to {r>1\mbox{ AU}} and prints this out.

After the loop, we calculate the final value of {r}, convert {\theta} and {r} to degrees and AU, respectively, and print out a couple of plots. The first plot is a polar plot of {r} as a function of {\theta}, so it shows the elliptical orbit. The second plot graphs {r} as a function of {t} (the latter in years) for one complete period.

For Halley’s comet, {P=76\mbox{ years}} and {e=0.9673} from which we find {a=17.943\mbox{ AU}}. If we input these values (along with M\_\_strsun = 1 since the central star is the Sun), and choose {n=10000} time increments, we get the time at which the comet first crosses {r=1\mbox{ AU}} after perihelion as {t=0.1064\mbox{ years}} or around 39 days. The plots are as follows:

Carroll & Ostlie 02.15a

Carroll & Ostlie 02.15b

Velocity in an elliptical orbit

Reference: Carroll, Bradley W. & Ostlie, Dale A. (2007), An Introduction to Modern Astrophysics, 2nd Edition; Pearson Education – Chapter 2, Problem 2.3.

We did a sample calculation of the Sun-Jupiter system using Kepler’s laws by approximating the orbits of the Sun and Jupiter about the centre of mass by circles. We can, however, derive equations for the radial and tangential velocity components for the correct case of elliptical orbits. We can start with the polar equation of an ellipse:

\displaystyle  r=\frac{a\left(1-e^{2}\right)}{1+e\cos\theta} \ \ \ \ \ (1)

The velocity of an object in polar coordinates is

\displaystyle   \mathbf{v} \displaystyle  = \displaystyle  v_{r}\hat{\mathbf{r}}+v_{\theta}\hat{\boldsymbol{\theta}}\ \ \ \ \ (2)
\displaystyle  \displaystyle  = \displaystyle  \dot{r}\hat{\mathbf{r}}+r\dot{\theta}\hat{\boldsymbol{\theta}} \ \ \ \ \ (3)

Differentiating 1 we get

\displaystyle   \dot{r} \displaystyle  = \displaystyle  \frac{dr}{d\theta}\dot{\theta}\ \ \ \ \ (4)
\displaystyle  \displaystyle  = \displaystyle  \frac{ae\left(1-e^{2}\right)\sin\theta}{\left(1+e\cos\theta\right)^{2}}\dot{\theta} \ \ \ \ \ (5)

In deriving Kepler’s laws, we got an expression for the total angular momentum of the system:

\displaystyle  L=\mu\sqrt{GMa\left(1-e^{2}\right)} \ \ \ \ \ (6)

In vector form, the angular momentum can be expressed in terms of the velocity of the centre of mass as

\displaystyle   \mathbf{L} \displaystyle  = \displaystyle  \mu\mathbf{r}\times\mathbf{v}\ \ \ \ \ (7)
\displaystyle  \displaystyle  = \displaystyle  \mu\mathbf{r}\times\left(\dot{r}\hat{\mathbf{r}}+r\dot{\theta}\hat{\boldsymbol{\theta}}\right)\ \ \ \ \ (8)
\displaystyle  \displaystyle  = \displaystyle  \mu r^{2}\dot{\theta}\hat{\mathbf{z}} \ \ \ \ \ (9)


\displaystyle  \dot{\theta}=\frac{\sqrt{GMa\left(1-e^{2}\right)}}{r^{2}} \ \ \ \ \ (10)

Substituting 1 into this, we get

\displaystyle  \dot{\theta}=\frac{\sqrt{GMa\left(1-e^{2}\right)}\left(1+e\cos\theta\right)^{2}}{a^{2}\left(1-e^{2}\right)^{2}} \ \ \ \ \ (11)

From Kepler”s third law relating the period {P} of the orbit to the semimajor axis {a}:

\displaystyle  P^{2}=\frac{4\pi^{2}}{GM}a^{3} \ \ \ \ \ (12)

we can eliminate {GM} to get

\displaystyle   \dot{\theta} \displaystyle  = \displaystyle  \frac{\sqrt{4\pi^{2}a^{4}\left(1-e^{2}\right)}\left(1+e\cos\theta\right)^{2}}{a^{2}\left(1-e^{2}\right)^{2}P}\ \ \ \ \ (13)
\displaystyle  \displaystyle  = \displaystyle  \frac{2\pi\left(1+e\cos\theta\right)^{2}}{P\left(1-e^{2}\right)^{3/2}} \ \ \ \ \ (14)

From this we can get the velocity components:

\displaystyle   v_{\theta} \displaystyle  = \displaystyle  r\dot{\theta}\ \ \ \ \ (15)
\displaystyle  \displaystyle  = \displaystyle  \frac{a\left(1-e^{2}\right)}{1+e\cos\theta}\frac{2\pi\left(1+e\cos\theta\right)^{2}}{P\left(1-e^{2}\right)^{3/2}}\ \ \ \ \ (16)
\displaystyle  \displaystyle  = \displaystyle  \frac{2\pi a\left(1+e\cos\theta\right)}{P\sqrt{1-e^{2}}}\ \ \ \ \ (17)
\displaystyle  v_{r} \displaystyle  = \displaystyle  \frac{dr}{d\theta}\dot{\theta}\ \ \ \ \ (18)
\displaystyle  \displaystyle  = \displaystyle  \frac{ae\left(1-e^{2}\right)\sin\theta}{\left(1+e\cos\theta\right)^{2}}\frac{2\pi\left(1+e\cos\theta\right)^{2}}{P\left(1-e^{2}\right)^{3/2}}\ \ \ \ \ (19)
\displaystyle  \displaystyle  = \displaystyle  \frac{2\pi ae\sin\theta}{P\sqrt{1-e^{2}}} \ \ \ \ \ (20)

The square magnitude of the velocity is then

\displaystyle   v^{2} \displaystyle  = \displaystyle  v_{\theta}^{2}+v_{r}^{2}\ \ \ \ \ (21)
\displaystyle  \displaystyle  = \displaystyle  \frac{4\pi^{2}a^{2}}{p^{2}\left(1-e^{2}\right)}\left(1+2e\cos\theta+e^{2}\left(\sin^{2}\theta+\cos^{2}\theta\right)\right)\ \ \ \ \ (22)
\displaystyle  \displaystyle  = \displaystyle  \frac{4\pi^{2}a^{2}\left(1+e^{2}+2e\cos\theta\right)}{P^{2}\left(1-e^{2}\right)} \ \ \ \ \ (23)

From 12 we get

\displaystyle  v^{2}=\frac{GM\left(1+e^{2}+2e\cos\theta\right)}{a\left(1-e^{2}\right)} \ \ \ \ \ (24)

From 1 we have

\displaystyle   \frac{2}{r}-\frac{1}{a} \displaystyle  = \displaystyle  \frac{2\left(1+e\cos\theta\right)}{a\left(1-e^{2}\right)}-\frac{1}{a}\ \ \ \ \ (25)
\displaystyle  \displaystyle  = \displaystyle  \frac{2\left(1+e\cos\theta\right)-\left(1-e^{2}\right)}{a\left(1-e^{2}\right)}\ \ \ \ \ (26)
\displaystyle  \displaystyle  = \displaystyle  \frac{1+e^{2}+2e\cos\theta}{a\left(1-e^{2}\right)} \ \ \ \ \ (27)

Therefore, comparing with 24 we get

\displaystyle  v^{2}=GM\left(\frac{2}{r}-\frac{1}{a}\right) \ \ \ \ \ (28)

Kepler’s laws

Reference: Carroll, Bradley W. & Ostlie, Dale A. (2007), An Introduction to Modern Astrophysics, 2nd Edition; Pearson Education – Chapter 2, Problem 2.6.

The total angular momentum of a two-body system interacting under gravity is

\displaystyle  \mathbf{L}=\mu\mathbf{r}\times\mathbf{v} \ \ \ \ \ (1)

where {\mu} is the reduced mass

\displaystyle  \mu=\frac{m_{1}m_{2}}{m_{1}+m_{2}} \ \ \ \ \ (2)

and {\mathbf{r}=\mathbf{r}_{2}-\mathbf{r}_{1}} is the position of {m_{2}} relative to {m_{1}}, and {\mathbf{v}=d\mathbf{r}/dt}. An important conservation law follows from taking the derivative with respect to time:

\displaystyle   \frac{d\mathbf{L}}{dt} \displaystyle  = \displaystyle  \mu\frac{d\mathbf{r}}{dt}\times\mathbf{v}+\mu\mathbf{r}\times\frac{d\mathbf{v}}{dt}\ \ \ \ \ (3)
\displaystyle  \displaystyle  = \displaystyle  \mu\mathbf{v}\times\mathbf{v}+\mathbf{r}\times\mathbf{F}\ \ \ \ \ (4)
\displaystyle  \displaystyle  = \displaystyle  0 \ \ \ \ \ (5)

where {\mathbf{r}\times\mathbf{F}=0} because gravity acts along the direction separating the two masses, so {\mathbf{r}\parallel\mathbf{F}}. Thus for any two objects in orbit about their centre of mass, angular momentum is conserved.

Using this fact and Newton’s gravitational force law, Carroll & Ostlie derive (in section 2.3) a formula for {r} as a function of the angle {\theta} measured from the perihelion point on the major axis:

\displaystyle  \boxed{r=\frac{L^{2}/\mu^{2}}{GM\left(1+e\cos\theta\right)}} \ \ \ \ \ (6)

where {M=m_{1}+m_{2}}.

Comparing this with the polar equation for an ellipse:

\displaystyle  r=\frac{a\left(1-e^{2}\right)}{1+e\cos\theta} \ \ \ \ \ (7)

we see that if the angular momentum is

\displaystyle  L=\mu\sqrt{GMa\left(1-e^{2}\right)} \ \ \ \ \ (8)

then 6 is the equation of an ellipse. This is Kepler’s first law, derived from Newtonian mechanics.

To translate this to motions of the two masses separately, we can use the relations between the position vectors of the masses and {\mathbf{r}}.

\displaystyle   \mathbf{r}_{1} \displaystyle  = \displaystyle  -\frac{\mu}{m_{1}}\mathbf{r}\ \ \ \ \ (9)
\displaystyle  \mathbf{r}_{2} \displaystyle  = \displaystyle  \frac{\mu}{m_{2}}\mathbf{r} \ \ \ \ \ (10)

Substituting these into 6 we get

\displaystyle   r_{1} \displaystyle  = \displaystyle  \frac{L^{2}}{GMm_{1}\mu\left(1+e\cos\theta\right)}\ \ \ \ \ (11)
\displaystyle  r_{2} \displaystyle  = \displaystyle  \frac{L^{2}}{GMm_{2}\left(1+e\cos\theta\right)} \ \ \ \ \ (12)

These equations still have the form of ellipses, so we see that both masses are in elliptical orbits about the centre of mass, which itself is fixed at one focus of the ellipse. Since the eccentricity {e} is the same for both masses and {L} is constant (it’s the total angular momentum of the two masses together), we can get the semimajor axis for each mass from

\displaystyle   \frac{L^{2}}{GMm_{i}\mu} \displaystyle  = \displaystyle  a_{i}\left(1-e^{2}\right)\ \ \ \ \ (13)
\displaystyle  a_{i} \displaystyle  = \displaystyle  \frac{L^{2}}{GMm_{i}\mu\left(1-e^{2}\right)} \ \ \ \ \ (14)

That is, the semimajor axis is inversely proportional to the mass of the object.

The constancy of angular momentum can also be used to derive Kepler’s second law, which states that the rate at which the radius vector of each mass sweeps out area is a constant. Carroll & Ostlie give the details, with the result

\displaystyle  \boxed{\frac{dA}{dt}=\frac{L}{2\mu}} \ \ \ \ \ (15)

where {A} is the area of the sector of the ellipse between the perihelion and the vector {\mathbf{r}} at time {t}. This result pertains to the reduced mass, but given that

\displaystyle  \frac{L}{\mu}=\left|\mathbf{r}\times\mathbf{v}\right| \ \ \ \ \ (16)

we can follow through the same derivation in Carroll & Ostlie to get the equation for each mass:

\displaystyle  \frac{dA_{i}}{dt}=\frac{L_{i}}{2m_{i}} \ \ \ \ \ (17)


\displaystyle  \mathbf{L}_{i}=m_{i}\mathbf{r}_{i}\times\mathbf{v}_{i} \ \ \ \ \ (18)

By substituting 9 and 10 into 3 we can see that the angular momentum of each mass separately is also conserved, so the equal areas in equal times law applies to each mass separately as well.

At perihelion and aphelion, {\mathbf{r}_{i}} and {\mathbf{v}_{i}} are perpendicular so at these points (for both masses):

\displaystyle  L_{i}=m_{i}r_{i(a,p)}v_{i(a,p)} \ \ \ \ \ (19)

and for the centre of mass

\displaystyle  L=\mu r_{a,p}v_{a,p} \ \ \ \ \ (20)

For any ellipse

\displaystyle   r_{p} \displaystyle  = \displaystyle  a\left(1-e\right)\ \ \ \ \ (21)
\displaystyle  r_{a} \displaystyle  = \displaystyle  a\left(1+e\right) \ \ \ \ \ (22)

so for the centre of mass, using 8

\displaystyle   v_{a} \displaystyle  = \displaystyle  \frac{L}{\mu r_{a}}\ \ \ \ \ (23)
\displaystyle  \displaystyle  = \displaystyle  \frac{\mu\sqrt{GMa\left(1-e^{2}\right)}}{\mu a\left(1+e\right)}\ \ \ \ \ (24)
\displaystyle  \displaystyle  = \displaystyle  \sqrt{\frac{GM\left(1-e\right)}{a\left(1+e\right)}}\ \ \ \ \ (25)
\displaystyle  v_{p} \displaystyle  = \displaystyle  \sqrt{\frac{GM\left(1+e\right)}{a\left(1-e\right)}} \ \ \ \ \ (26)

Finally, we can get Kepler’s third law by integrating 15 over one orbital period {P} to get

\displaystyle  A=\frac{L}{2\mu}P \ \ \ \ \ (27)

Since the area of an ellipse is {A=\pi ab=\pi a^{2}\sqrt{1-e^{2}}} we have

\displaystyle   \pi a^{2}\sqrt{1-e^{2}} \displaystyle  = \displaystyle  \frac{1}{2}\sqrt{GMa\left(1-e^{2}\right)}P\ \ \ \ \ (28)
\displaystyle  P^{2} \displaystyle  = \displaystyle  \frac{4\pi^{2}}{GM}a^{3} \ \ \ \ \ (29)

The square of the period is proportional to the cube of the semimajor axis.

Example If we look at the Sun and Jupiter (ignoring all the other planets), we can plug in a few numbers. The eccentricity of Jupiter’s orbit is {e=0.048}, its period is {P=11.86\mbox{ yr}}, its mass is {0.0009546M_{S}} and its semimajor axis is {a=5.2\mbox{ AU}}. In these units, {4\pi^{2}/GM=1}, so from 8, the total angular momentum of the system is

\displaystyle  L=2\pi\mu\sqrt{a\left(1-e^{2}\right)} \ \ \ \ \ (30)

The reduced mass is

\displaystyle  \mu=\frac{0.0009546M_{S}^{2}}{1.0009546M_{S}}=0.000954M_{S} \ \ \ \ \ (31)


\displaystyle  L=0.01367M_{S}\mbox{ AU}^{2}\mbox{yr}^{-1} \ \ \ \ \ (32)

To find the contribution of the Sun to {L}, we need its semimajor axis, which we can get from 14

\displaystyle   a_{S} \displaystyle  = \displaystyle  \frac{\left(0.01365M_{S}\right)^{2}}{4\pi^{2}\left(0.000945\right)M_{S}^{2}\left(1-e^{2}\right)}\ \ \ \ \ (33)
\displaystyle  \displaystyle  = \displaystyle  0.005\mbox{ AU} \ \ \ \ \ (34)

In order to calculate {L_{S}} from 19, we need the velocity {v} at aphelion or perihelion, which at present we don’t know. However, if we approximate the sun’s orbit by a circle, then {v} is constant for the orbit and we have

\displaystyle  v_{S}=\frac{2\pi a_{S}}{P}=0.00265\mbox{ AU yr}^{-1} \ \ \ \ \ (35)

Then the Sun’s angular momentum is

\displaystyle  L_{S}=M_{S}a_{S}v_{S}=1.3\times10^{-5}M_{S}\mbox{ AU}^{2}\mbox{yr}^{-1} \ \ \ \ \ (36)

Doing the same calculation for Jupiter, we find

\displaystyle   v_{J} \displaystyle  = \displaystyle  \frac{2\pi a}{P}=2.755\mbox{ AU yr}^{-1}\ \ \ \ \ (37)
\displaystyle  L_{J} \displaystyle  = \displaystyle  0.0009546M_{S}\left(5.2\right)\left(2.755\right)\ \ \ \ \ (38)
\displaystyle  \displaystyle  = \displaystyle  0.01367M_{S}\mbox{ AU}^{2}\mbox{ yr}^{-1} \ \ \ \ \ (39)

Essentially all of the orbital angular momentum is due to Jupiter.

The rotational angular momentum of a solid sphere is

\displaystyle   L \displaystyle  = \displaystyle  I\omega\ \ \ \ \ (40)
\displaystyle  \displaystyle  = \displaystyle  \frac{2}{5}mr^{2}\omega \ \ \ \ \ (41)

where {I} is the moment of inertia, {r} is the radius of the sphere and {\omega} is the angular velocity. Although neither the Sun nor Jupiter is a solid sphere, we can get approximate values by assuming they are. The rotation period {D_{S}} of the Sun is around 26 days or 0.0712 years. Jupiter’s ‘day’ is only about 10 hours, or 0.00114 years. The angular velocities are

\displaystyle   \omega_{S} \displaystyle  = \displaystyle  \frac{2\pi}{D_{S}}=87.88\mbox{ yr}^{-1}\ \ \ \ \ (42)
\displaystyle  \omega_{J} \displaystyle  = \displaystyle  \frac{2\pi}{D_{J}}=5511\mbox{ yr}^{-1} \ \ \ \ \ (43)

and the radii are

\displaystyle   r_{S} \displaystyle  = \displaystyle  0.00465\mbox{ AU}\ \ \ \ \ (44)
\displaystyle  r_{J} \displaystyle  = \displaystyle  0.00047\mbox{ AU} \ \ \ \ \ (45)


\displaystyle   L_{S} \displaystyle  = \displaystyle  \frac{2}{5}M_{S}r_{S}^{2}\omega_{S}=7.6\times10^{-4}M_{S}\mbox{ AU}^{2}\mbox{yr}^{-1}\ \ \ \ \ (46)
\displaystyle  L_{J} \displaystyle  = \displaystyle  \frac{2}{5}\left(0.0009546M_{S}\right)r_{J}^{2}\omega_{J}=4.65\times10^{-7}M_{S}\mbox{ AU}^{2}\mbox{ yr}^{-1} \ \ \ \ \ (47)

Virtually all the angular momentum in the Sun-Jupiter system is due to the orbital motion of Jupiter.

Orbits in the centre of mass frame: energy and angular momentum

Reference: Carroll, Bradley W. & Ostlie, Dale A. (2007), An Introduction to Modern Astrophysics, 2nd Edition; Pearson Education – Chapter 2, Problems 2.4-2.5.

For two masses {m_{1}} and {m_{2}} interacting via gravity, it’s easiest to do calculations in the centre of mass frame. The position of the centre of mass is defined as

\displaystyle  \mathbf{R}\equiv\frac{m_{1}\mathbf{r}_{1}'+m_{2}\mathbf{r}_{2}'}{m_{1}+m_{2}} \ \ \ \ \ (1)

where {\mathbf{r}_{i}'} is the position of {m_{i}} in a coordinate system where the origin could be anywhere.

In the centre of mass frame, we take the centre of mass to be at the origin, so that {\mathbf{R}=0}, giving

\displaystyle  \frac{m_{1}\mathbf{r}_{1}+m_{2}\mathbf{r}_{2}}{m_{1}+m_{2}}=0 \ \ \ \ \ (2)

where now {\mathbf{r}_{i}} (without the prime) indicates the position relative to the centre of mass. If we introduce the relative position

\displaystyle  \mathbf{r}\equiv\mathbf{r}_{2}-\mathbf{r}_{1} \ \ \ \ \ (3)


\displaystyle   \frac{m_{1}\mathbf{r}_{1}+m_{2}\left(\mathbf{r}+\mathbf{r}_{1}\right)}{m_{1}+m_{2}} \displaystyle  = \displaystyle  0\ \ \ \ \ (4)
\displaystyle  \mathbf{r}_{1} \displaystyle  = \displaystyle  -\frac{m_{2}}{m_{1}+m_{2}}\mathbf{r}\ \ \ \ \ (5)
\displaystyle  \frac{m_{1}\left(\mathbf{r}_{2}-\mathbf{r}\right)+m_{2}\mathbf{r}_{2}}{m_{1}+m_{2}} \displaystyle  = \displaystyle  0\ \ \ \ \ (6)
\displaystyle  \mathbf{r}_{2} \displaystyle  = \displaystyle  \frac{m_{1}}{m_{1}+m_{2}}\mathbf{r} \ \ \ \ \ (7)

In terms of the reduced mass

\displaystyle  \mu\equiv\frac{m_{1}m_{2}}{m_{1}+m_{2}} \ \ \ \ \ (8)

we get

\displaystyle   \mathbf{r}_{1} \displaystyle  = \displaystyle  -\frac{\mu}{m_{1}}\mathbf{r}\ \ \ \ \ (9)
\displaystyle  \mathbf{r}_{2} \displaystyle  = \displaystyle  \frac{\mu}{m_{2}}\mathbf{r} \ \ \ \ \ (10)

The velocities of the two masses are just the derivatives of the positions, so that

\displaystyle  \mathbf{v}_{i}=\frac{d\mathbf{r}_{i}}{dt} \ \ \ \ \ (11)

and the rate of change of the relative position is

\displaystyle  \mathbf{v}=\frac{d\mathbf{r}}{dt} \ \ \ \ \ (12)

The total energy of the two mass system is the sum of the two kinetic energy terms and the gravitational potential energy term, so

\displaystyle  E=\frac{1}{2}m_{1}\left|\mathbf{v}_{1}\right|^{2}+\frac{1}{2}m_{2}\left|\mathbf{v}_{2}\right|^{2}-G\frac{m_{1}m_{2}}{r} \ \ \ \ \ (13)


\displaystyle  r=\left|\mathbf{r}_{2}-\mathbf{r}_{1}\right| \ \ \ \ \ (14)

The gravitational potential energy can be written in terms of the reduced mass and the total mass {M\equiv m_{1}+m_{2}} as

\displaystyle  -G\frac{m_{1}m_{2}}{r}=-G\left(m_{1}+m_{2}\right)\frac{\mu}{r}=-G\frac{M\mu}{r} \ \ \ \ \ (15)

The kinetic energy can be rewritten by using 9 and 10

\displaystyle   \mathbf{v}_{1} \displaystyle  = \displaystyle  \frac{d\mathbf{r}_{1}}{dt}\ \ \ \ \ (16)
\displaystyle  \displaystyle  = \displaystyle  -\frac{\mu}{m_{1}}\mathbf{v}\ \ \ \ \ (17)
\displaystyle  \mathbf{v}_{2} \displaystyle  = \displaystyle  \frac{\mu}{m_{2}}\mathbf{v}\ \ \ \ \ (18)
\displaystyle  \frac{1}{2}m_{1}\left|\mathbf{v}_{1}\right|^{2}+\frac{1}{2}m_{2}\left|\mathbf{v}_{2}\right|^{2} \displaystyle  = \displaystyle  \frac{1}{2}\frac{\mu^{2}}{m_{1}}v^{2}+\frac{1}{2}\frac{\mu^{2}}{m_{2}}v^{2}\ \ \ \ \ (19)
\displaystyle  \displaystyle  = \displaystyle  \frac{1}{2}\mu^{2}\frac{m_{1}+m_{2}}{m_{1}m_{2}}v^{2}\ \ \ \ \ (20)
\displaystyle  \displaystyle  = \displaystyle  \frac{1}{2}\mu v^{2} \ \ \ \ \ (21)

Thus the total energy is

\displaystyle  E=\frac{1}{2}\mu v^{2}-G\frac{M\mu}{r} \ \ \ \ \ (22)

which is the energy of a mass {\mu} moving around another mass {M}, the latter of which is fixed at the origin.

The orbital angular momentum (that is, the angular momentum due to the masses orbiting about each other, not including any angular momentum due to each mass rotating on its own axis) is

\displaystyle   \mathbf{L} \displaystyle  = \displaystyle  m_{1}\mathbf{r}_{1}\times\mathbf{v}_{1}+m_{2}\mathbf{r}_{2}\times\mathbf{v}_{2}\ \ \ \ \ (23)
\displaystyle  \displaystyle  = \displaystyle  -m_{1}\frac{\mu}{m_{1}}\mathbf{r}\times\left(-\frac{\mu}{m_{1}}\mathbf{v}\right)+m_{2}\frac{\mu}{m_{2}}\mathbf{r}\times\left(\frac{\mu}{m_{2}}\mathbf{v}\right)\ \ \ \ \ (24)
\displaystyle  \displaystyle  = \displaystyle  \mu^{2}\mathbf{r}\times\mathbf{v}\left(\frac{1}{m_{1}}+\frac{1}{m_{2}}\right)\ \ \ \ \ (25)
\displaystyle  \displaystyle  = \displaystyle  \frac{m_{1}+m_{2}}{m_{1}m_{2}}\mu^{2}\mathbf{r}\times\mathbf{v}\ \ \ \ \ (26)
\displaystyle  \displaystyle  = \displaystyle  \mu\mathbf{r}\times\mathbf{v} \ \ \ \ \ (27)

Thus the angular momentum is due to the reduced mass alone, which is consistent with this mass orbiting about a fixed mass at the origin.

Elliptical orbits

Reference: Carroll, Bradley W. & Ostlie, Dale A. (2007), An Introduction to Modern Astrophysics, 2nd Edition; Pearson Education – Chapter 2, Problems 2.1-2.2.

In a two-body solar system (that is, where we have only the Sun and a single planet or other object orbiting it), it turns out that the planet’s orbit is a conic section (ellipse, parabola or hyperbola) with the centre of mass of the system at one of the foci (or the single focus, in the case of a parabola) of the orbit. For any bound orbit (that is, one in which the planet never escapes to infinity), the orbital curve is an ellipse, so it’s useful to review a few properties of ellipses.

To define an ellipse, first specify two points which serve as the foci (see Fig. 2.4 in Carroll & Ostlie). The ellipse is then the set of points such that the sum of the distances of a given point from the two foci is a constant, namely {2a}, where {a} is called the semimajor axis (and {2a} is the major axis). If the distance from a point {P} on the ellipse to one focus is {r'} and to the other focus is {r}, then

\displaystyle  r'+r=2a \ \ \ \ \ (1)

If we place the ellipse in a rectangular coordinate system with the line between the foci on the {x} axis and the origin at the centre of the ellipse, then the distance along the {y} axis from the origin to the ellipse is known as {b}, the semiminor axis. The eccentricity {e} of the ellipse is defined as the distance between the foci divided by the major axis, so that the distance between the foci is {2ae}. If the two foci coincide, the distance between them is zero so {e=0}, {r'=r=a} and we get a circle. If the foci coincide with the ends of the major axis, then {ae=a} and {e=1}, giving just a line segment along the {x} axis extending from {x=-a} to {x=+a}. In the general case, if we draw a right angled triangle with vertices at one of the foci, the centre, and the point {\left(0,b\right)}, then the hypotenuse of this triangle is {a} so by Pythagoras we have

\displaystyle   a^{2}e^{2}+b^{2} \displaystyle  = \displaystyle  a^{2}\ \ \ \ \ (2)
\displaystyle  b^{2} \displaystyle  = \displaystyle  a^{2}\left(1-e^{2}\right)\ \ \ \ \ (3)
\displaystyle  e^{2} \displaystyle  = \displaystyle  1-\frac{b^{2}}{a^{2}} \ \ \ \ \ (4)

If we use a polar coordinate system with the origin at the right-hand focus, then the angle {\theta} between the {x} axis and the line {r} from that focus to the ellipse is the polar angle. To get the equation of an ellipse in polar coordinates, we can consider the triangle with its vertices at the polar origin, the point on the ellipse with coordinates {\left(r,\theta\right)} and the left-hand focus. The triangle’s angle at the origin is then {\pi-\theta} so by the law of cosines we have

\displaystyle   \left(r'\right)^{2} \displaystyle  = \displaystyle  4a^{2}e^{2}+r^{2}-4rae\cos\left(\pi-\theta\right)\ \ \ \ \ (5)
\displaystyle  \displaystyle  = \displaystyle  4a^{2}e^{2}+r^{2}+4rae\cos\theta \ \ \ \ \ (6)

From 1 we have

\displaystyle  \left(r'\right)^{2}=4a^{2}+r^{2}-4ar \ \ \ \ \ (7)


\displaystyle   4a^{2}-4ar \displaystyle  = \displaystyle  4a^{2}e^{2}+4rae\cos\theta\ \ \ \ \ (8)
\displaystyle  r \displaystyle  = \displaystyle  \frac{a\left(1-e^{2}\right)}{1+e\cos\theta}\ \ \ \ \ (9)
\displaystyle  \displaystyle  = \displaystyle  \frac{b^{2}}{a\left(1+e\cos\theta\right)} \ \ \ \ \ (10)

We can derive a rectangular coordinate version of the equation for an ellipse by two methods. We can start with 9 and shift the origin to the centre of the ellipse, so that a point on the ellipse has rectangular coordinates

\displaystyle   x \displaystyle  = \displaystyle  ae+r\cos\theta\ \ \ \ \ (11)
\displaystyle  y \displaystyle  = \displaystyle  r\sin\theta \ \ \ \ \ (12)

Now (since we know the answer we’re looking to prove) we can form the following quantity, using 4 and 10:

\displaystyle   \frac{x^{2}}{a^{2}}+\frac{y^{2}}{b^{2}} \displaystyle  = \displaystyle  \frac{1}{a^{2}}\left(a\sqrt{1-\frac{b^{2}}{a^{2}}}+r\cos\theta\right)^{2}+\frac{r^{2}\sin^{2}\theta}{b^{2}}\ \ \ \ \ (13)
\displaystyle  \displaystyle  = \displaystyle  \frac{1}{a^{2}}\left(a\sqrt{1-\frac{b^{2}}{a^{2}}}+\frac{b^{2}\cos\theta}{a\left(1+\sqrt{1-\frac{b^{2}}{a^{2}}}\cos\theta\right)}\right)^{2}+\frac{b^{2}\sin^{2}\theta}{a^{2}\left(1+\sqrt{1-\frac{b^{2}}{a^{2}}}\cos\theta\right)^{2}}\ \ \ \ \ (14)
\displaystyle  \displaystyle  = \displaystyle  \frac{1}{a^{4}\left(1+\sqrt{1-\frac{b^{2}}{a^{2}}}\cos\theta\right)^{2}}\left[\left(a^{2}\sqrt{1-\frac{b^{2}}{a^{2}}}\left(1+\sqrt{1-\frac{b^{2}}{a^{2}}}\cos\theta\right)+b^{2}\cos\theta\right)^{2}+a^{2}b^{2}\sin^{2}\theta\right]\ \ \ \ \ (15)
\displaystyle  \displaystyle  = \displaystyle  \frac{a^{4}}{a^{4}\left(1+\sqrt{1-\frac{b^{2}}{a^{2}}}\cos\theta\right)^{2}}\left[\left(1-\frac{b^{2}}{a^{2}}\right)\cos^{2}\theta+2\sqrt{1-\frac{b^{2}}{a^{2}}}\cos\theta+1\right]\ \ \ \ \ (16)
\displaystyle  \displaystyle  = \displaystyle  \frac{\left(1+\sqrt{1-\frac{b^{2}}{a^{2}}}\cos\theta\right)^{2}}{\left(1+\sqrt{1-\frac{b^{2}}{a^{2}}}\cos\theta\right)^{2}}\ \ \ \ \ (17)
\displaystyle  \displaystyle  = \displaystyle  1 \ \ \ \ \ (18)

The other method is to use the definition 1 and Pythagoras again. We have

\displaystyle   \left(r'\right)^{2} \displaystyle  = \displaystyle  \left(x+ae\right)^{2}+y^{2}\ \ \ \ \ (19)
\displaystyle  r^{2} \displaystyle  = \displaystyle  \left(x-ae\right)^{2}+y^{2} \ \ \ \ \ (20)


\displaystyle  \sqrt{\left(x+ae\right)^{2}+y^{2}}+\sqrt{\left(x-ae\right)^{2}+y^{2}}=2a \ \ \ \ \ (21)

Squaring both sides, we get

\displaystyle   \left(x+ae\right)^{2}+\left(x-ae\right)^{2}+2y^{2}+2\sqrt{\left(x+ae\right)^{2}+y^{2}}\sqrt{\left(x-ae\right)^{2}+y^{2}} \displaystyle  = \displaystyle  4a^{2}\ \ \ \ \ (22)
\displaystyle  x^{2}+y^{2}+a^{2}e^{2}+\sqrt{\left(x+ae\right)^{2}+y^{2}}\sqrt{\left(x-ae\right)^{2}+y^{2}} \displaystyle  = \displaystyle  2a^{2}\ \ \ \ \ (23)
\displaystyle  \sqrt{\left(x^{2}-a^{2}e^{2}\right)^{2}+y^{4}+y^{2}\left(\left(x+ae\right)^{2}+\left(x-ae\right)^{2}\right)} \displaystyle  = \displaystyle  a^{2}+b^{2}-x^{2}-y^{2}\ \ \ \ \ (24)
\displaystyle  \sqrt{y^{4}+2x^{2}y^{2}+x^{4}+2a^{2}e^{2}\left(y^{2}-x^{2}\right)+a^{4}e^{4}} \displaystyle  = \displaystyle  a^{2}+b^{2}-x^{2}-y^{2}\ \ \ \ \ (25)
\displaystyle  y^{4}+2x^{2}y^{2}+x^{4}+2a^{2}e^{2}\left(y^{2}-x^{2}\right)+a^{4}e^{4} \displaystyle  = \displaystyle  \left(a^{2}+b^{2}-x^{2}-y^{2}\right)^{2}\ \ \ \ \ (26)
\displaystyle  2x^{2}y^{2}+2\left(a^{2}-b^{2}\right)\left(y^{2}-x^{2}\right)+\left(a^{2}-b^{2}\right)^{2} \displaystyle  = \displaystyle  \left(a^{2}+b^{2}\right)^{2}-2\left(a^{2}+b^{2}\right)\left(x^{2}+y^{2}\right)+2x^{2}y^{2}\ \ \ \ \ (27)
\displaystyle  4b^{2}x^{2}+4a^{2}y^{2} \displaystyle  = \displaystyle  4a^{2}b^{2}\ \ \ \ \ (28)
\displaystyle  \frac{x^{2}}{a^{2}}+\frac{y^{2}}{b^{2}} \displaystyle  = \displaystyle  1 \ \ \ \ \ (29)

To find the area of an ellipse, we can again use either polar or rectangular coordinates. To find the area in polar coordinates, consider the area of a thin wedge of radius {r} and angular extent {d\theta}. If {r} were constant over the range {\theta\in\left[0,2\pi\right]}, we’d have the area of a circle, or {\pi r^{2}}. The wedge has a fraction {d\theta/2\pi} of this area, so the area of the wedge is {\frac{1}{2}r^{2}d\theta}, and the area of the ellipse is

\displaystyle  A=\int_{0}^{2\pi}\frac{1}{2}r^{2}d\theta \ \ \ \ \ (30)

Using 10 we get

\displaystyle  A=\frac{b^{4}}{2a^{2}}\int_{0}^{2\pi}\frac{d\theta}{\left(1+e\cos\theta\right)^{2}} \ \ \ \ \ (31)

I’m not sure how you would do this integral by hand, but for the indefinite integral Maple gives a complicated result involving arctans and tangents. Putting in the limits, however, and restricting {e} to {0\le e<1} gives a simple answer (using 3):

\displaystyle  A=\frac{\pi b^{4}}{a^{2}\left(1-e^{2}\right)^{3/2}}=\pi\frac{b^{2}}{\sqrt{1-e^{2}}}=\pi ab \ \ \ \ \ (32)

Using rectangular coordinates is somewhat easier. We can work out the area for the first quadrant and multiply by 4:

\displaystyle   A \displaystyle  = \displaystyle  4b\int_{0}^{a}\sqrt{1-\frac{x^{2}}{a^{2}}}dx\ \ \ \ \ (33)
\displaystyle  \displaystyle  = \displaystyle  2ab\left[\arctan\left(\frac{x}{\sqrt{a^{2}-x^{2}}}\right)+x\sqrt{a^{2}-x^{2}}\right]_{0}^{a}\ \ \ \ \ (34)
\displaystyle  \displaystyle  = \displaystyle  \pi ab \ \ \ \ \ (35)

Angular distances on the celestial sphere

Reference: Carroll, Bradley W. & Ostlie, Dale A. (2007), An Introduction to Modern Astrophysics, 2nd Edition; Pearson Education – Chapter 1, Problems 1.8 – 1.11.

Sometimes we need to know the angular distance between two points on the celestial sphere. By angular distance, I mean the angle subtended by lines drawn from the centre of the Earth to each of the two points. For example, the angular distance between the north pole and any point on the equator is {90^{\circ}}. The angular diameters of both the moon and the sun are around {30'} or half a degree. One practical application of such a calculation is to determine if two stars that are close to each other in the sky will both be visible within the field of view of a telescope using a particular eyepiece, since an eyepiece’s field of view is typically given as an angle.

Since the coordinates of objects are typically given in right ascension (RA) and declination (Dec), we’d like a formula that gives the angular distance between two sets of such coordinates. The derivation of the formula uses spherical trigonometry and is given in full in Carroll & Ostlie’s Section 1.3 (in particular Fig. 1.17 and surrounding text), so we won’t go through the details here. The exact formulas relating the angular distance {\Delta\theta} between a point {P} with coordinates {\left(RA,Dec\right)=\left(\alpha,\delta\right)} and an adjacent point {Q} with coordinates {\left(\alpha+\Delta\alpha,\delta+\Delta\delta\right)} are

\displaystyle   \sin\left(\Delta\alpha\right)\cos\left(\delta+\Delta\delta\right) \displaystyle  = \displaystyle  \sin\left(\Delta\theta\right)\sin\phi\ \ \ \ \ (1)
\displaystyle  \cos\left[90^{\circ}-\left(\delta+\Delta\delta\right)\right] \displaystyle  = \displaystyle  \cos\left(90^{\circ}-\delta\right)\cos\left(\Delta\theta\right)+\sin\left(90^{\circ}-\delta\right)\sin\left(\Delta\theta\right)\cos\phi \ \ \ \ \ (2)

where {\phi} is the angle, measured counterclockwise as we look at the sky, between a line due north from {P} and the line {PQ} (by ‘line’ I mean a segment of a great circle, since these lines are drawn on a sphere), measured . In principle, {\phi} is of course determined by {P} and {Q} so we should be able to eliminate it from these equations, but it turns out that {P} and {Q} are quite close together, there is an easier way to eliminate {\phi}, as we’ll see now.

If we assume that {\Delta\delta} and {\Delta\alpha} are both small enough, we can use the approximations {\sin x\approx x} and {\cos x\approx1} together with {\cos\left(90^{\circ}-x\right)=\sin x}, {\sin\left(90^{\circ}-x\right)=\cos x} and the trig formulas for the sum and difference of angles to simplify the above formulas. From 1 we get, by keeping only first order terms in {\Delta\delta} and {\Delta\alpha}:

\displaystyle   \Delta\alpha\left(\cos\delta-\sin\Delta\delta\sin\delta\right) \displaystyle  = \displaystyle  \Delta\theta\sin\phi\ \ \ \ \ (3)
\displaystyle  \Delta\alpha\cos\delta \displaystyle  = \displaystyle  \Delta\theta\sin\phi \ \ \ \ \ (4)

From 2 we get

\displaystyle   \sin\left(\delta+\Delta\delta\right) \displaystyle  = \displaystyle  \sin\delta+\Delta\theta\cos\delta\cos\phi\ \ \ \ \ (5)
\displaystyle  \sin\delta+\Delta\delta\cos\delta \displaystyle  = \displaystyle  \sin\delta+\Delta\theta\cos\delta\cos\phi\ \ \ \ \ (6)
\displaystyle  \Delta\delta\cos\delta \displaystyle  = \displaystyle  \Delta\theta\cos\delta\cos\phi\ \ \ \ \ (7)
\displaystyle  \Delta\delta \displaystyle  = \displaystyle  \Delta\theta\cos\phi \ \ \ \ \ (8)

We can now square 4 and 8 and add them to get

\displaystyle  \boxed{\left(\Delta\theta\right)^{2}=\left(\Delta\alpha\cos\delta\right)^{2}+\left(\Delta\delta\right)^{2}} \ \ \ \ \ (9)

Note that the approximation {\sin x\approx x} assumes that {x} is in radians, not degrees. However, because each term in this equation contains the square of an angle, the conversion factor from degrees to radians cancels out, so the formula is valid whether we specify angles in degrees or radians (as long as you remember to take {\delta} in the correct units when taking the cosine).

Example 1 The closest star system to the Sun is the {\alpha} Centauri triple star system. Of the three stars in this system, Proxima Centauri ({\alpha} Centauri C) is the closest to us. The coordinates of Proxima Centauri in the year 2000 (known as the epoch J2000.0) were {\left(\alpha,\delta\right)_{C}=\left(14^{h}29^{m}42.95^{s},-62^{\circ}40'46.1^{\prime\prime}\right)} (unfortunately for those of us in the far northern hemisphere, the {\alpha} Centauri system’s extreme southern declination means we can never see it). The brightest star in the system ({\alpha} Centauri A) has J2000.0 coordinates of {\left(\alpha,\delta\right)_{A}=\left(14^{h}39^{m}36.50^{s},-60^{\circ}50'02.3^{\prime\prime}\right)}. Since these two points are quite close to each other, we can use 9 to see how far apart they appear in the sky. We can first convert the coordinates into degrees in decimal form. The conversion for {\alpha} is {1^{h}=15^{\circ}}, so {1^{m}} is 1/60 of {15^{\circ}} or {15'=0.25^{\circ}}. Likewise, {1^{s}=15^{\prime\prime}=15/3600^{\circ}}. Therefore {14^{h}29^{m}42.95^{s}=14\times15+29\times0.25+42.95\times15/3600=217.42896^{\circ}}.

\displaystyle   \left(\alpha,\delta\right)_{C} \displaystyle  = \displaystyle  \left(217.42896^{\circ},-62.67947^{\circ}\right)\ \ \ \ \ (10)
\displaystyle  \left(\alpha,\delta\right)_{A} \displaystyle  = \displaystyle  \left(219.90208^{\circ},-60.83397^{\circ}\right) \ \ \ \ \ (11)

We therefore have

\displaystyle   \Delta\alpha \displaystyle  = \displaystyle  2.47312^{\circ}\ \ \ \ \ (12)
\displaystyle  \Delta\delta \displaystyle  = \displaystyle  1.8455^{\circ}\ \ \ \ \ (13)
\displaystyle  \cos\delta \displaystyle  = \displaystyle  \cos\left(-62.67947^{\circ}\right)\ \ \ \ \ (14)
\displaystyle  \displaystyle  = \displaystyle  0.45897\ \ \ \ \ (15)
\displaystyle  \Delta\theta \displaystyle  = \displaystyle  \sqrt{\left(\Delta\alpha\cos\delta\right)^{2}+\left(\Delta\delta\right)^{2}}\ \ \ \ \ (16)
\displaystyle  \displaystyle  = \displaystyle  2.1666^{\circ} \ \ \ \ \ (17)

If the distance to the {\alpha} Centauri system is {r=4.0\times10^{16}\mbox{ m}} (and we assume that all members of the system are the same distance from us, or at least that the radial distances between the three stars is negligible as a fraction of their distance from us), then Proxima Centauri and {\alpha} Centauri A are separated by

\displaystyle  d_{C\rightarrow A}=r\Delta\theta \ \ \ \ \ (18)

This time, we do have to use radians, so we get

\displaystyle  d_{C\rightarrow A}=\left(4.0\times10^{16}\right)\frac{\pi}{180}\left(2.1666\right)=1.512\times10^{15}\mbox{ m} \ \ \ \ \ (19)

This is about 3.78% of the system’s distance from us.

Example 2 Proper motion and precession. The position of a star as seen from Earth can change due to two effects: the star’s actual motion through space, known as proper motion, and the precession of the Earth’s axis. The latter effect is due to the direction of the Earth’s axis slowly changing so that the north celestial pole moves through a circle in the sky with a period of about 26,000 years. Both these effects cause a star’s coordinates to change. In the previous example, we gave the coordinates of the {\alpha} Centauri system for epoch J2000.0. Astronomical coordinates are typically updated for the effects of precession every 50 years, but of course it is sometimes necessary to know the precise coordinates of a star for a time between updates. Calculating these positions exactly is a complicated business, but there are approximate formulas that are usually good enough. The changes in coordinates due to precession relative to J2000.0 are given by

\displaystyle   \Delta\alpha \displaystyle  = \displaystyle  M+N\sin\alpha\tan\delta\ \ \ \ \ (20)
\displaystyle  \Delta\delta \displaystyle  = \displaystyle  N\cos\alpha \ \ \ \ \ (21)


\displaystyle   M \displaystyle  = \displaystyle  1^{\circ}.2812323T+0^{\circ}.0003879T^{2}+0^{\circ}.0000101T^{3}\ \ \ \ \ (22)
\displaystyle  N \displaystyle  = \displaystyle  0^{\circ}.5567530T-0^{\circ}.0001185T^{2}-0^{\circ}.0000116T^{3}\ \ \ \ \ (23)
\displaystyle  T \displaystyle  = \displaystyle  \left(t-2000.0\right)/100 \ \ \ \ \ (24)

and {t} is the current date, specified as the current year plus a fraction for the current date.

For Proxima Centauri in the year 2010, using these formulas gives

\displaystyle   \Delta\alpha \displaystyle  = \displaystyle  0^{\circ}.1936\ \ \ \ \ (25)
\displaystyle  \Delta\delta \displaystyle  = \displaystyle  -0^{\circ}.0442 \ \ \ \ \ (26)

Substituting into 9 we get

\displaystyle   \Delta\theta_{precession} \displaystyle  = \displaystyle  0^{\circ}.0992\ \ \ \ \ (27)
\displaystyle  \displaystyle  = \displaystyle  357^{\prime\prime}\approx6' \ \ \ \ \ (28)

The proper motion of Proxima Centauri is measured to be {3.84^{\prime\prime}} per year with a position angle of {\phi=282^{\circ}}. Over 10 years, proper motion results in a shift of {\Delta\theta_{pm}=38.4^{\prime\prime}}, so the precession effect is larger. We can get the shifts in RA and Dec from 4 and 8:

\displaystyle   \Delta\alpha_{pm} \displaystyle  = \displaystyle  \Delta\theta_{pm}\frac{\sin\phi}{\cos\delta}\ \ \ \ \ (29)
\displaystyle  \displaystyle  = \displaystyle  -81.84^{\prime\prime}=-5.46^{s}\ \ \ \ \ (30)
\displaystyle  \Delta\delta_{pm} \displaystyle  = \displaystyle  \Delta\theta_{pm}\cos\phi\ \ \ \ \ (31)
\displaystyle  \displaystyle  = \displaystyle  7.98^{\prime\prime} \ \ \ \ \ (32)

Julian dates

Reference: Carroll, Bradley W. & Ostlie, Dale A. (2007), An Introduction to Modern Astrophysics, 2nd Edition; Pearson Education – Chapter 1, Problem 1.7.

The various units commonly used to measure the passage of time, such as days, weeks, months and years, are not conducive to simple calculations since such quantities as determining the number of days between two dates are not easy to calculate. In astronomy, most time calculations refer to the time difference between two events (for example, the orbital period of a planet or the period of variability of a star), so a method of specifying time that makes such differences easy to calculate is preferred.

The most common such system is the Julian date, in which each day, starting from noon Universal Time (UT or Greenwich Mean Time (GMT), which is the time zone used in the UK during the winter months) on January 1, 4713 BC, is numbered sequentially. The historic reason for this choice of starting date relies on choosing some rather obscure periods (one of which is a 15 year cycle the Romans used for calculating land tax), and isn’t important in astronomy, but if you’re interested, I’ll refer you to the Wikipedia article.

The advantage of using Julian dates is, of course, that if the JD values for two events are known, the time between these events is just the difference in JDs. Of course, if we know the calendar date of an event, it’s a bit of a pain to calculate the corresponding JD, but online converters exist, such as this one provided by the US Naval Observatory.

In order to calculate a JD for a date near to the present, it’s useful to know the JD for a reference date in our current era. The JD of noon UT January 1, 2000 is

\displaystyle  JD2000=2451545.0 \ \ \ \ \ (1)

The JD of a time of day other than noon is calculated by adding on the fraction of a day since noon that has elapsed, so that midnight on January 2, 2000 (that is, 0:00 hours on January 2) is 2451545.5.

For a particular date such as July 14, 2006, 16:15 hours UT, we can work out the JD starting with JD2000 by counting the number of days since JD2000 (remembering to take leap years into account). The period from Jan 1, 2000 to Jan 1, 2006 includes 2 leap years (2000 and 2004) so there are {6\times365+2=2192} days between these dates. Between noon Jan 1, 2006 and noon July 14, 2006, there are {31+28+31+30+31+30+13=194} full days, so we’re now up to {2192+194=2386} complete days since Jan 1, 2000. The time from 12:00 to 16:15 is {4.25/24=0.177083} days, so the JD of July 14, 2006, 16:15 hours UT is

\displaystyle  JD=JD2000+2386+0.177083=2453931.177083 \ \ \ \ \ (2)

There are several variants of the Julian date that are sometimes used. The Modified Julian Date (MJD) is defined as

\displaystyle  MJD\equiv JD-2400000.5 \ \ \ \ \ (3)

which is equivalent to setting the starting date of {MJD=0.0} at 0h November 17, 1858 (thus the MJD starts at midnight UT rather than noon). Historically, the MJD was introduced in 1957 to allow the computers of the day (with very limited memory) to track the orbit of the first artificial satellite, Sputnik. Reducing the size of the numbers allowed the computer to handle the calculations. Thus the MJD of July 14, 2006, 16:15 hours UT is

\displaystyle  MJD=2453931.177083-2400000.5=53930.677083 \ \ \ \ \ (4)

Other variants include the Reduced JD, with starting date of 12:00 on November 16, 1858, Truncated JD (0:00, May 24, 1968), Dublin JD (12:00, December 31, 1899) and Chronological JD (0:00, January 1, 4713 BC, but adjusted for local time zone). All these variants, however, still number days sequentially from their starting date so they are all equally useful for calculating time differences between events.

Astronomical coordinates: declination and right ascension

Reference: Carroll, Bradley W. & Ostlie, Dale A. (2007), An Introduction to Modern Astrophysics, 2nd Edition; Pearson Education – Chapter 1, Problems 1.4 – 1.6.

Although objects visible in the sky are all at different distances from Earth, for the purposes of observing them from Earth-bound telescopes we can treat the objects in the sky as if they were all on the surface of a sphere whose centre is at the centre of the Earth. This sphere is known as the celestial sphere. The coordinates of an object in the sky can then be specified using only two coordinates, analogous to the latitude and longitude used to identify places on Earth. We can get celestial latitudes by simply projecting Earth’s latitudes onto the celestial sphere. The celestial latitude is called declension (Dec for short) and is measured in degrees ranging from {-90^{\circ}} at the south celestial pole (the point in the sky that is a projection of Earth’s south pole onto the celestial sphere) to {+90^{\circ}} at the north celestial pole.

When we try to specify celestial longitude we can’t just project Earth’s longitude onto the celestial sphere since Earth rotates relative to the stars so the longitude of a point on the Earth’s surface that is directly under some particular star changes continuously over the course of a day. A set of longitude lines is therefore specified that is fixed relative to the stars’ positions on the celestial sphere. This ‘celestial longitude’ is called right ascension and, rather than being measured in degrees as terrestrial longitude is, it is measured in hours, minutes and seconds, with values ranging from 0 hours around (in an easterly direction) to 24 hours.

Although right ascension (or RA as it is more commonly abbreviated) is measured in units of time, we need to be careful in relating RA time to terrestrial time. The 24 hours of RA correspond to one complete rotation of Earth relative to the stars, NOT the Sun! To see the difference, picture the Earth in its orbit about the Sun. Over the course of one complete orbit about the Sun, Earth goes through {360^{\circ}}. As the year is 365.25 days, Earth travels just under {1^{\circ}} of its orbit each day. That means that Earth must rotate about {361^{\circ}} between successive times at which the Sun is on the meridian (that is, at its highest point in the sky for a given day). However, Earth needs to rotate only {360^{\circ}} between successive times at which some particular distant star is on the meridian. The former time interval (which is what we think of as an ordinary day) is called a solar day, while the latter time interval is a sidereal day. A solar day is about 4 minutes longer than a sidereal day, as that’s how long it takes Earth to rotate through {1^{\circ}}. The units of RA measure sidereal time, not solar time, so the time it takes for one complete rotation of the celestial sphere (that is, the time taken to run through one complete cycle of RA from 0 hours up to the next 0 hours) is about 4 minutes less than 24 hours measured by your wristwatch. Another way of putting it is that any given star rises about 4 minutes earlier each day. It is this mismatch between sidereal and solar time that causes the constellations we see in the night sky to shift slowly over the course of a year.

The zero hour of RA is defined as the RA coordinate at the point where the Sun crosses the celestial equator heading north, which is the spring or vernal equinox in the northern hemisphere, and the autumnal equinox in the southern hemisphere and occurs around March 21. This point (0 hours RA and {0^{\circ}} Dec) is called the first point of Aries although confusingly it is actually located in Pisces. [The shift is due to Earth’s precession, or wobble in its axis, which we’ll get to eventually. A long time ago, this point actually was in Aries. Incidentally, if you believe in astrology, you should know that the dates of the various star signs are based on the positions of the constellations several thousand years ago. For example, the star sign of Aries runs from March 21 to April 19, during most of which time the Sun is actually in Pisces. The Sun doesn’t actually enter Aries until around April 17. Yet another reason, if one were needed, not to take horoscopes seriously.]

Due to the tilt of the Earth’s axis of around {23.5^{\circ}}, the path of the Sun across the sky (known as the ecliptic) varies between extremes of declination of {+23.5} (around June 21, at which point the Sun’s coordinates are {RA=6^{h};\; Dec=+23.5^{\circ}}) through the equinox around September 21 ({RA=12^{h};\; Dec=0^{\circ}}) to the winter (summer) solstice in the northern (southern) hemisphere around December 21 ({RA=18^{h};\; Dec=-23.5^{\circ}}).

The altitude (angle relative to the horizon) of an object depends on the latitude from which it is observed. At the north pole, the north celestial pole is directly overhead, so an object with {Dec=90^{\circ}} is directly overhead. In general, an object with Dec {\delta} has altitude {\delta} when seen from the north pole. At the north pole, the altitude of any given object is a constant; it just rotates around the sky in a circle around the zenith (the point directly overhead).

As we move south, the celestial north pole moves downwards, so that at latitude {L} the altitude of the celestial north pole is also {L}. Since objects rotate around the celestial north pole, an object with Dec {\delta} has an altitude that depends on where in the sky it is. When it’s on the meridian, its altitude is highest and at that point, the altitude is {90-L+\delta}. This means that an object with {L=\delta} is directly overhead on the meridian. The highest point in the sky reached by the Sun (on the summer solstice) at a given latitude is thus

\displaystyle A=90-L+23.5=113.5-L \ \ \ \ \ (1)

For latitudes less than 23.5, this gives a value greater than 90 which just indicates that the Sun appears past the zenith at an angle with the zenith of {A-90}. For my latitude of +56.5 north, the highest the Sun gets is {57^{\circ}}.

At the other extreme (the winter solstice), the Sun’s Dec is {-23.5} so the altitude is

\displaystyle A=90-L-23.5=66.5-L \ \ \ \ \ (2)

For my latitude, the winter sun thus never gets more than {10^{\circ}} above the horizon on December 21.

Since the north celestial pole has an altitude {L} at latitude {L}, any objects with declinations in the range {90-L<\delta<90} will never set as seen from that latitude and are known as circumpolar. Similarly, stars with declinations {-90<\delta<L-90} will never rise and are permanently invisible from that latitude. At the north pole, all stars with {\delta>0} are circumpolar, and only the northern half of the sky is ever visible. At the equator, no stars are circumpolar, but all stars are visible at some point.

Midnight sun at the summer solstice therefore occurs when the Sun, at declination {\delta=23.5}, is circumpolar, which occurs for latitudes {66.5<L<90}. At the vernal equinox the Sun’s declination is {\delta=0}, so it will never set only at the north and south poles, where it will just skim the horizon all the way round.

Sidereal and synodic periods of planets

Reference: Carroll, Bradley W. & Ostlie, Dale A. (2007), An Introduction to Modern Astrophysics, 2nd Edition; Pearson Education – Chapter 1, Problems 1.1 – 1.3.

As I recently got a shiny new telescope (for looking at the night sky, not the neighbours!), I thought it would be a good idea to revisit some astrophysics since it’s been something like 40 years since I studied it properly. The book I’ve chosen for this is the one referenced above, and it will take a while for me to work through it, since it has 30 chapters and comes close to 1500 pages. However, since I live in Scotland where clear nights are rare and the sky never gets truly dark between May and July (due to being so far north; the town of Monifieth where I live is almost as far north as Juneau, Alaska), having a form of astronomy that doesn’t depend on actually looking at the sky is very useful.

So let’s begin with some basics about planetary orbits. For most of history, people thought that the Earth was the centre of the universe. Looking at the sky, particularly at night, reinforced this view since all objects visible in the sky appear to rotate around the Earth and since the Earth felt solid an immovable (except during earthquakes), it was just common sense that the Earth was fixed and everything else revolved around it. This form of common sense is the same kind that got Aristotle in trouble by stating that the natural state of objects was at rest (since everything on Earth seemed to come to rest, given long enough) and that heavier objects fell faster than lighter ones (well, a cannonball falls faster than a feather anyway).

In 1543, Polish astronomer Nicolaus Copernicus published On the Revolution of the Celestial Sphere, in which he proposed that the Sun, not the Earth, was the true centre of the universe. Copernicus was still a bit of a mystic, however, so he assumed that the orbits of everything around the Sun were all geometrically ‘perfect’ circles. Since the planets’ orbits are all elliptical to some extent, this severely limited Copernicus’s theory’s ability to predict planetary positions, but it was still a giant leap forward in astronomical thought.

One consequence of the Copernican theory is that it is relatively easy to calculate the synodic period of a planet if we know its sidereal period. First, we need to define these two terms.

The sidereal period of a planet is simply the time it takes to make one complete orbit around the Sun. Thus Earth’s sidereal period is 1 year.

To understand the concept of a synodic period, we need to think about some observational astronomy. First, consider a superior planet (that is, a planet whose orbit is further from the Sun than Earth’s, such as Mars). When the positions of Earth and Mars are just right, Mars is directly opposite the Sun in the sky (that is, the angle between a line from Earth to Sun and a line from Earth to Mars is {\pi} or {180^{\circ}}). When this happens, Mars is said to be in opposition. After opposition, both planets continue in their orbits until some future time at which they line up again, giving another opposition. The time between two successive oppositions is the synodic period.

For inferior planets (planets whose orbits are closer to the Sun than Earth’s, of which there are only two: Venus and Mercury), opposition can never be achieved since it’s impossible for the Earth to get between the Sun and the planet. However, the Earth, the inferior planet and the Sun can still line up in two different ways. When the planet is directly between the Earth and the Sun, it is in inferior conjunction, at which point it’s at its closest position to Earth. When the planet is directly on the other side of the Sun, it is in superior conjunction, at which point it is at its furthest position from Earth. (Superior planets can be in superior conjunction as well, of course, but never in inferior conjunction. As a result, superior conjunction for a superior planet is usually called just ‘conjunction’.)

When Venus or Mercury is in inferior conjuction with the Earth, the Earth is in opposition as seen from Venus or Mercury.

We can work out the relation between the sidereal and synodic periods using a bit of geometry. Consider a superior planet such as Mars, and let’s suppose that at {t=0}, Mars is in opposition. How long will it be before the next opposition? Since Mars takes longer to orbit the Sun, its angular speed {\omega_{M}} is smaller than that of Earth {\omega_{E}}. After a synodic period {S} has elapsed, the position angles of Earth and Mars must be equal again, but since Earth is moving faster, it will have made one more complete orbit around the Sun than Mars has. The angle swept out by Earth in time {S} is {\omega_{E}S}, so

\displaystyle  \omega_{E}S=\omega_{M}S+2\pi \ \ \ \ \ (1)

In terms of the sidereal period of Earth {\omega_{E}=2\pi/P_{E}} (with a similar relation for Mars), so

\displaystyle   \frac{2\pi}{P_{E}}S \displaystyle  = \displaystyle  \frac{2\pi}{P_{M}}S+2\pi\ \ \ \ \ (2)
\displaystyle  \frac{1}{S} \displaystyle  = \displaystyle  \frac{1}{P_{E}}-\frac{1}{P_{M}} \ \ \ \ \ (3)

If all times are in years, then for any superior planet with period {P}

\displaystyle  \frac{1}{S}=1-\frac{1}{P} \ \ \ \ \ (4)

For an inferior planet, the roles of Earth and the planet are reversed, so

\displaystyle  \frac{1}{S}=\frac{1}{P}-1 \ \ \ \ \ (5)

All of this, of course, assumes perfectly circular orbits with planets moving at constant speeds, so the formulas aren’t exact. However, we can get an idea of how good they are by plugging in some actual numbers.

Planet {P} (years) {S} (calc) {S} (actual)
Mercury 0.241 0.3175 0.3176
Venus 0.616 1.604 1.599
Mars 1.9 2.111 2.136
Jupiter 11.9 1.092 1.092
Saturn 29.5 1.035 1.035
Uranus 84.0 1.012 1.013
Neptune 164.8 1.0061 1.0075
Pluto 248.5 1.0040 1.0048

The agreement is suprisingly good for such a simple formula. As you might expect, the further out the planet is, the closer {S} gets to 1 Earth year, since the outer planets don’t move that much in their own orbits over the course of a year. Thus the shortest synodic period for a superior planet is for the planet that is furthest away, namely Pluto (or, if you insist on degrading Pluto to the status of a dwarf planet, Neptune). For amateur astronomers, this means that each superior planet (Jupiter and beyond) is best placed for observing (which happens when it’s at opposition) roughly once per year. Mars is at closest approach to Earth about every two years.

Historically, the synodic period of a planet would be the available datum for the planet (since it’s easy to measure the time from one opposition or conjunction to the next), and the formula above would be inverted to give the planet’s actual period as

\displaystyle  P=\begin{cases} \left(1-\frac{1}{S}\right)^{-1} & \mbox{superior planet}\\ \left(1+\frac{1}{S}\right)^{-1} & \mbox{inferior planet} \end{cases} \ \ \ \ \ (6)

Given the observational data available to someone, such as Copernicus, in the pre-telescope age, we can work out the relative ordering of the planets from the Sun. Since Mercury and Venus are never observed in opposition, they must have orbits closer to the Sun than Earth, and since Mercury’s greatest elongation (angular separation from the Sun) is less than that of Venus, we can deduce that Mercury must be closer to the Sun than Venus.

For the superior planets, we can order them according to decreasing synodic period, which gives the correct ordering of Mars through Pluto (or up to Saturn, since Uranus, Neptune and Pluto were discovered only after the invention of the telescope).

Relativistic electromagnetic potentials

Reference: Griffiths, David J. (2007), Introduction to Electrodynamics, 3rd Edition; Pearson Education – Chapter 12, Problem 12.56.

Maxwell’s equations can be written in terms of the electromagnetic field tensor, its dual and the current four-vector as

\displaystyle   \partial_{j}F^{ij} \displaystyle  = \displaystyle  \mu_{0}J^{i}\ \ \ \ \ (1)
\displaystyle  \partial_{j}G^{ij} \displaystyle  = \displaystyle  0 \ \ \ \ \ (2)


\displaystyle   F^{ij} \displaystyle  = \displaystyle  \left[\begin{array}{cccc} 0 & E_{x} & E_{y} & E_{z}\\ -E_{x} & 0 & B_{z} & -B_{y}\\ -E_{y} & -B_{z} & 0 & B_{x}\\ -E_{z} & B_{y} & -B_{x} & 0 \end{array}\right]\ \ \ \ \ (3)
\displaystyle  G^{ij} \displaystyle  = \displaystyle  \left[\begin{array}{cccc} 0 & B_{x} & B_{y} & B_{z}\\ -B_{x} & 0 & -E_{z} & E_{y}\\ -B_{y} & E_{z} & 0 & -E_{x}\\ -B_{z} & -E_{y} & E_{x} & 0 \end{array}\right]\ \ \ \ \ (4)
\displaystyle  J^{i} \displaystyle  = \displaystyle  \left[\rho,J_{x},J_{y},J_{z}\right]\ \ \ \ \ (5)
\displaystyle  \displaystyle  = \displaystyle  \frac{\rho_{0}}{\sqrt{1-\beta^{2}}}\left[1,u_{x},u_{y},u_{z}\right] \ \ \ \ \ (6)

It turns out that an even more compact form for Maxwell’s equations can be written using the 4-vector potential

\displaystyle   A^{i} \displaystyle  = \displaystyle  \left(\frac{V}{c},A_{x},A_{y},A_{z}\right)\ \ \ \ \ (7)
\displaystyle  \displaystyle  = \displaystyle  \left(V,A_{x},A_{y},A_{z}\right) \ \ \ \ \ (8)

where the last line uses relativistic units where {c=1}.

Griffiths shows in section 12.3.5 that the field tensor can be written in terms of the potentials as

\displaystyle  F^{ij}=\partial^{i}A^{j}-\partial^{j}A^{i} \ \ \ \ \ (9)

Note that we’re using the contravariant gradient operator here, in order to get the signs right on the time components. Because of the form of {F^{ij}}, the gauge invariance shows up naturally, since if we replace {A^{i}} by

\displaystyle  A^{i}\rightarrow A^{i}+\partial^{i}\lambda \ \ \ \ \ (10)

where {\lambda} is any scalar function, {F^{ij}} is unchanged, as the order in which the partial derivatives of {\lambda} are taken doesn’t matter, so {\lambda} drops out of the equation for {F^{ij}}. The Lorentz gauge condition is

\displaystyle  \nabla\cdot\mathbf{A}=-\frac{\partial V}{\partial t}=-\partial_{0}V \ \ \ \ \ (11)

[Notice we’re back to using the covariant gradient operator.] This can be condensed to read

\displaystyle  \partial_{i}A^{i}=0 \ \ \ \ \ (12)

[Incidentally, Griffiths’s equation 12.135 is wrong; it should read {\partial A^{\mu}/\partial x^{\mu}=0}.]

Combining 1 with 9 gives

\displaystyle  \partial_{j}\partial^{i}A^{j}-\partial_{j}\partial^{j}A^{i}=\mu_{0}J^{i} \ \ \ \ \ (13)

If we use the Lorentz gauge, the first term is zero, so we get

\displaystyle   \partial_{j}\partial^{j}A^{i} \displaystyle  = \displaystyle  -\mu_{0}J^{i}\ \ \ \ \ (14)
\displaystyle  \Box^{2}A^{i} \displaystyle  = \displaystyle  -\mu_{0}J^{i} \ \ \ \ \ (15)

where the symbol {\Box^{2}} is the d’Alembertian operator, defined as

\displaystyle  \Box^{2}\equiv\partial_{j}\partial^{j} \ \ \ \ \ (16)

We can verify that the other Maxwell equation 2 is also satisfied by the potential formulation by using the earlier result

\displaystyle  \partial_{j}G^{ij}=\partial_{a}F_{bc}+\partial_{b}F_{ca}+\partial_{c}F_{ab} \ \ \ \ \ (17)

where {i}, {a}, {b} and {c} are all different.

Lowering the indexes on 9 we get

\displaystyle  F_{ij}=\partial_{i}A_{j}-\partial_{j}A_{i} \ \ \ \ \ (18)

Substituting 18 into 17 we get

\displaystyle  \partial_{j}G^{ij}=\partial_{a}\partial_{b}A_{c}-\partial_{a}\partial_{c}A_{b}+\partial_{b}\partial_{c}A_{a}-\partial_{b}\partial_{a}A_{c}+\partial_{c}\partial_{a}A_{b}-\partial_{c}\partial_{b}A_{a}=0 \ \ \ \ \ (19)


Get every new post delivered to your Inbox.

Join 514 other followers