Experiment 2: The Gravitational Force and Orbits of Comets¶

Aims¶

The aim of this experiment is to study basic celestial mechanics by computing the orbit of a comet about the Sun. We will explore the variation of the kinetic and potential energies of the comet with time, the angular momentum of the comet, and the eccentricity and curvature of its orbit.

Universal Law of Gravitation¶

To understand the interaction between a comet and the Sun, we can begin by looking at the simple case of two arbitrary objects. Sir Isaac Newton asserted that all objects exert a gravitational force on all other objects. He found that the magnitude of this force depended on three parameters: direction, mass, and distance.

The gravitational force between two objects is attractive, thus we can write that the force exerted on object 2 by object 1 (${\bf F}_{1\rightarrow 2}$) is directed in the opposite direction of the position vector between object 1 and 2. For example, if we think about the gravitational force that the earth (object 1) is exerting on us (object 2), the force is pulling us toward the earth (object 1). We can write this relationship

\begin{equation} {\bf \hat{F}}_{1\rightarrow2} = {\bf \hat{r}}_{2\rightarrow1} = -{\bf\hat{r}}_{1\rightarrow2}, \end{equation}

where the $\hat{}$ symbol denotes a unit vector. Furthermore, we know from Newton's Third Law that the gravitational force is dependent on the mass of both objects, thus

\begin{equation} {\bf F}_{1\rightarrow2} \propto m_{1}m_{2}. \end{equation}

Our daily experience gives us the intuition that the force between two objects becomes weaker as the distance between them becomes larger. Newton confirmed that this is indeed true. More specifically, he found that the force is inversely proportional to the square of the distance:

\begin{equation} {\bf F}_{1\rightarrow2} \propto \frac{1}{r^{2}}, \end{equation}

where $r^{2} = {\bf r}_{1\rightarrow2}^{2} = {\bf r}_{2\rightarrow1}^{2}$. Combining the previous three equations, we obtain the proportionality relationship \begin{equation} {\bf F}{1\rightarrow2} \propto -\frac{m{1}m{2}}{r^{2}}{\bf \hat{r}}{1\rightarrow2}. \end{equation}

We can then define a proportionality constant G, to write the Universal Law of Gravitation:

\begin{equation} {\bf F}_{1\rightarrow2} = -G \frac{m_{1}m_{2}}{r^{2}}{\bf \hat{r}}_{1\rightarrow2}. \end{equation}

The universal gravitation constant $G$, has been determined by experiment to be $6.67\times 10^{-11}$ $\mathrm{m}^3/\mathrm{kg}\,\mathrm{s}^2$.

In this experiment, we will consider the specific case of a comet and the Sun. The mass of the Sun is much larger than the mass of a comet, so we assume the Sun is stationary and simplify the previous equation to:

\begin{equation*} \mathbf{F} = - \frac{GMm}{\left\|\mathbf{r}\right\|^3} \, \mathbf{r} \end{equation*}

where $\mathbf{r}$ is the position of the comet, $m$ is its mass, $M$ ($=1.99\times 10^{30}$ kg) is the mass of the Sun, and we have substituted $\frac{{\bf r}}{||{\bf r}||}$ for ${\bf \hat{r}}$.

The natural units of length and time for this problem are not metres and seconds. Rather, as a unit of distance we will use the astronomical unit (1 AU $=1.496\times 10^{11}$ $\mathrm{m}$), which equals the mean Earth-Sun distance. The unit of time will be the AU year (the period of a circular orbit of radius 1 AU). In these units, the product $GM=4\pi^2$ $\mathrm{AU}^3/\mathrm{yr}^2$. We take the mass of the comet, $m$, to be unity; in MKS units the typical mass of a comet is $10^{15\pm 3}$ kg.

Simulating the Equations of Motion¶

In this experiment, we will use a Copernican coordinate system and fix the Sun at the origin. We assume that the trajectory of the comet lies in the $xy$-plane, so that we may write the position and velocity vectors as \begin{eqnarray} \mathbf{r} = \left( \begin{array}{c} x \ y \ 0 \end{array} \right) \;\;\;\; \textrm{and} \;\;\;\; \mathbf{v} = \left( \begin{array}{c} v_x \ v_y \ 0 \end{array} \right) . \end{eqnarray} The equations of motion of these two vectors are \begin{eqnarray*} \frac{d\mathbf{r}}{dt} = \mathbf{v} \;\;\;\; \textrm{and} \;\;\;\; \frac{d\mathbf{v}}{dt} = \mathbf{a} = \frac{\mathbf{F}}{m} =

\frac{GM}{\left|\mathbf{r}\right|^3} \, \mathbf{r} \, . \end{eqnarray*}

The script below -- comet -- solves the equations of motion using a higher-order Runga-Kutta integration method (the linear approximation, or Euler, method is not very suitable for this problem as it requires extremely small timesteps in order to avoid "blowing up" due to numerical instability).

The program stores the position and velocity vectors of the comet in the arrays rplot and vplot, respectively. In particular, the position vector after j timesteps is given by rplot[j,:], and the corresponding velocity vector is given by vplot[j,:]. The time after j timesteps is stored in tplot[j].

Note: The first and second indices of a 2D Numpy array correspond to the row and column of that array, respectively. For example, rplot[2,1] corresponds to the element in the third row, second column of the rplot array. The colon : denots all elements in a row or column. For example, rplot[1,:] will take every element in the second row, corresponding to the y-value of the position vector.

Run comet, and experiment with a variety of different input parameters defined by the sliders above the figure. For example, start with the following inputs, which correspond to a comet starting at the position $\mathbf{r}=(1,0,0)$ AU with velocity $\mathbf{v}=(0,4,0)$ AU/yr:

Initial radial distance, $r_{0}$ (AU) = 1

Initial tangential velocity, $v_{0}$ (AU/yr) = 4

Number of steps, $n$ = 500

Time step, $\tau$ (yr) = 0.001

You will see the comet's trajectory drawn in the $xy$-plane. Note its elliptical shape, with the Sun at one focus of the ellipse.

"""
comet - Program to compute the orbit of a comet.
"""
%matplotlib inline
%pylab inline
from math import pi
import numpy as np
from numpy import linalg as la
from scipy import integrate
import matplotlib.pyplot as plt

GM = 4*pi**2  # Grav. const. * Mass of Sun (au^3/yr^2)
mass = 1.0  # Mass of comet

def gravity(state, t, GM):
	"""Return RHS of Kepler ODE"""
	r = state[:3]
	v = state[3:]
	a = -GM*r/la.norm(r)**3
	return np.concatenate((v, a))
	
def orbit(r0, v0, n, tau, t0=0.0):
	# Set initial position and velocity of the comet.
	r1 = np.array([r0, 0., 0.])
	v1 = np.array([0., v0, 0.])
	state0 = np.concatenate((r1, v1))
	t = np.linspace(0.0, n*tau, n+1) + t0
	
	# Use SciPy standard ODE integrator (lsoda from odepack).
	states = integrate.odeint(gravity, state0, t, (GM,))
	r = states[:,:3]
	v = states[:,3:]		
	return t, r, v

def plot(r0,v0,n,tau):
# Set physical parameters.
	tplot, rplot, vplot = orbit(r0, v0, n, tau)
	
	# Graph the trajectory of the comet.
	plt.figure(1)
	plt.clf()
	plt.xlim((-5,5))
	plt.ylim((-5,5))
  
	plt.plot(rplot[:,0], rplot[:,1])
	plt.xlabel('Distance (AU)')
	plt.ylabel('Distance (AU)')
	plt.grid(True)
	width = 100
	height = 100
	plt.figure(figsize=(width, height))
	plt.show()
    
from ipywidgets import interactive
from IPython.display import display
import ipywidgets as widgets  

r0_slider   = widgets.FloatSlider(value = 1., min=0., max=2.5, step=0.1, description='$r_0$ (AU)')
v0_slider   = widgets.FloatSlider(value = 4., min=0., max=4.5, step=0.1, description='$v_0 (AU/yr)$')
n_slider   = widgets.IntSlider(value = 500, min=0, max=10000, step=500, description='$n$')
tau_fill   = widgets.FloatText(value = 0.001,description='$\\tau (yr)$')

w = interactive(plot,r0=r0_slider,v0=v0_slider,n=n_slider,tau = tau_fill)
display(w)

<matplotlib.figure.Figure at 0x7f23db811438>

1. Try different values of the initial tangential velocity. Describe what happens to the orbit as this initial velocity is made larger or smaller. (2 points)

Note: for small values of the initial tangential velocity, you may need to decrease the size of the timestep. Otherwise, numerical instability causes the trajectory to deviate from an ellipse. For large values of the initial velocity you may need to use more timesteps to see the full shape of the trajectory.

Kinetic Energy and Potential Energy¶

The kinetic energy $E_K$ and gravitational potential energy $E_P$ of the comet are given by \begin{eqnarray} E_K = \frac{1}{2}m\left|\mathbf{v}\right|^2 \;\;\;\; \textrm{and} \;\;\;\; E_P = - \frac{GMm}{\left|\mathbf{r}\right|} . \end{eqnarray} Note: The magnitude of a two-component vector is defined: $||{\bf a} || = \sqrt{a_{1}^{2}+a_{2}^{2}}$. Also, the comet script defines the constant GM = 4*pi**2 and defines mass = 1 to be the mass of the comet.

2. Complete the kepler script below to compute and plot $E_K$, $E_P$, and their sum $E_K+E_P$ versus time (tplot) for a particular trajectory (i.e., for a particular rplot and vplot generated by the comet script above). Label the axes of your plot, and describe in words what your results illustrate in a markdown cell. (3 points)

"""
kepler - Program to compute the orbit of a comet.
Write your code below.
"""
from pylab import *

# define the magnitude of the position and velocity vector, followed by the kinetic and potential energy
rlength = 
speed = 
Ek = 
Ep = 

# plot time versus Ek, Ep, and Ek + Ep

Angular Momentum¶

The angular momentum of the comet about the Sun is \begin{eqnarray} \mathbf{L} = m\,(\mathbf{r}\times\mathbf{v}) . \end{eqnarray}

3. Using the product rule for derivatives and the properties of the cross product and the gravitational force, show that the angular momentum of the comet is constant with time (i.e., calculate $d\mathbf{L}/dt$ and show that it is zero). The derivation is started for you in the cell below: (1 point)

$\frac{d{\bf L}}{dt} = m\frac{d}{dt}{\bf r} \times {\bf v}$

4. (a) Write a line or two of code below to compute and display the angular momentum vector of the comet at each timestep. Note that since the motion is entirely in the $x$-$y$ plane, $\mathbf{L}$ is in the $z$ direction. (1 point)

Hint: use the cross() function in the Numpy package to compute the cross product. Use help() if you are unsure how to use this function.

(b) Confirm for several different trajectories that $\mathbf{L}$ is a constant, and give its value for each case in the table below. Note that this will require running comet for a given set of parameters, and subsequently evaluating your code for $\mathbf{L}$|, and then repeating this process for additional trajectories until you are convinced $\mathbf{L}$ is a constant. (1 point)

$r_{0}$	$v_{0}$	$L_{z}$

5. Optional: Referring to Figure 1 below, the area swept out by the position vector in the (small) time interval $\Delta t$ is approximated by $\Delta A=\frac{1}{2}\left\|\mathbf{r}(t)\times\mathbf{r}(t+\Delta t)\right\|$. (Why?) Using this, show that the rate of change of the area, $\Delta A/\Delta t$, is a constant. This amounts to Kepler's Second Law of Planetary Motion -- the line joining the Sun and a planet (or orbiting comet) sweeps out equal areas in equal times. In other words, prove that $\Delta A/\Delta t$ is constant. ($\frac{1}{2}$ point)

Figure 1: Change in position vector in time interval $\Delta t$.

Eccentricity of the Orbit¶

The eccentricity of an elliptical orbit is given by the formula \begin{eqnarray} e = \sqrt{1-\frac{b^2}{a^2}} \, , \end{eqnarray} where $a$ and $b$ are the semimajor and semiminor axes of the ellipse, respectively (see Figure 2).

Figure 2: Elliptical orbit about the Sun.

Consider the following combinations of initial position and velocity vectors:

(i) $\mathbf{r}(t=0)=(1,0,0)$ AU, $\mathbf{v}(t=0)=(0,3,0)$ AU/yr

(ii) $\mathbf{r}(t=0)=(1,0,0)$ AU, $\mathbf{v}(t=0)=(0,4,0)$ AU/yr

(iii) $\mathbf{r}(t=0)=(1,0,0)$ AU, $\mathbf{v}(t=0)=(0,5,0)$ AU/yr

6. Run comet with the above inputs. For each set of initial conditions:

(a) Use rplot to determine $a$ and $b$ for each case and record in the table below. (1 point)

Hint: the min() and max() functions will find the minimum and maximum of an input vector or array, respectively. The absolute value function abs() may also be of use.

(b) Compute the eccentricity $e$ and value of $h=\left\|\mathbf{r}\times\mathbf{v}\right\| =\left\|\mathbf{L}/m\right\|$ for each case and type your results in the table below. (1 point)

	a	b	e	h
(i)
(ii)
(iii)

Polar Coordinates¶

Another way of representing the (two-dimensional) orbit is in terms of polar coordinates $r=\left\|\mathbf{r}\right\| =\sqrt{x^2+y^2}$ and $\theta = \tan^{-1}(y/x)$.

7.(a) Write a line of code to compute the angle $\theta$ at each point along the trajectory. Use the Python function arctan2(y,x) for the inverse tangent. (1 point)

(b) You can now use the command polar(theta,rlength) to generate a polar plot of a particular orbit. Try this for one set of initial conditions. (1 point)

Using calculus and geometry, one can show theoretically that the orbit of the comet in polar coordinates is in fact given by the expression \begin{eqnarray} r = \frac{h^2/GM}{1-e\cos (\theta )} \, , \end{eqnarray} which describes an ellipse of eccentricity $e$ with one focus at the origin (i.e., at the Sun). This, of course, is Kepler's First Law of Planetary Motion.

8. Set the initial conditions in the comet figure to the initial conditions (i) listed above Question 6. With the conditions set, click Cell > Run All Above to evaluate all of the cells in the notebook to this point. Then, use the following commands to superimpose the simulated and theoretical orbits on a polar graph. Repeat this process for initial conditions (ii) and (iii) and describe your results in a markdown cell. (2 points)

polar(theta,rlength)
hold(True)
phi = linspace(-pi,pi,20,endpoint=False)
polar(phi,(h**2/GM)/(1-e*cos(phi)),'*')
hold(False)

Curvature¶

The curvature of a curve at a given point is a measure of how quickly the curve changes direction at that point. Formally, it is defined as the magnitude of the rate of change of the unit tangent vector to the curve with respect to arc length, i.e., \begin{eqnarray} \kappa = \left| \frac{d\mathbf{T}}{ds} \right| , \;\;\;\; \textrm{where} \;\;\;\; \mathbf{T} = \frac{d\mathbf{r}/dt}{\left|d\mathbf{r}/dt\right|} = \frac{\mathbf{v}}{\left|\mathbf{v}\right|} , \end{eqnarray} and $s$ is the arc length. For a small timestep $\Delta t$, an approximation to the curvature is given by \begin{eqnarray} \kappa (t) = \frac{\left|\mathbf{T}(t+\Delta t)-\mathbf{T}(t)\right|} {\left|\mathbf{r}(t+\Delta t)-\mathbf{r}(t)\right|} . \end{eqnarray}

9.(a) Finish the incomplete lines in the following code cell (compute dT and ds) to calculate an estimate for the curvature $\kappa$ at each point along the trajectory according to the approximate formula given above. (2 point)

Hint: Examining the previous two equations for $\kappa$, it is clear that $d{\bf T}(t) = ||{\bf T}(t+\Delta t) - {\bf T}(t)||$. Furthermore, ${\bf T}(t+\Delta t) = \frac{{\bf v}(t+\Delta t)}{||{\bf v}(t+\Delta t)||}$.

(b) Add the code to plot curvature versus time and plot for initial conditions (i) and (ii). You will end up with one fewer values of curvature than time, so need to write plot(tplot[:-1],kappa) to plot the curvature versus time. Explain what you see in a markdown cell. (2 point)

# calculate kappa
kappa = np.zeros(np.shape(tplot))[:-1] # create an array filled with zeros that is the same length as tplot
for j in range(len(kappa)): # this line loops through *every value* in kappa
    dT = 
    ds = 
    kappa[j] = norm(dT)/norm(ds)
    
# plot kappa vs. time

10. Estimate the maximum and minimum values of the curvature, $\kappa_\textrm{max}$ and $\kappa_\textrm{min}$, for both cases, and identify the points (give $x$ and $y$ coordinates) on the trajectory where they occur. (2 points)

Hint: Use Cell > Run All Above and examine the figure produced by comet, and the previous figure produced in Question 9 for the initial conditions (i) and (ii). Changing the number of time steps may help to identify patterns in the graph of kappa.

Final Notes¶

Review your notebook and ensure you have answered all of the questions and all of your code, text, and figures are embedded. Then, select File > Download as > HTML (.html) and submit this file on Canvas.