We now use the approximate masses of the sun and earth system.
m_sun <- 2.0e30
m_earth <- 6.0e24
df <- data.frame(m = c(m_sun, m_earth), x = c(100, 150), y = c(120, 50))Whoops, the earth disappears!
The sizes on the figure correspond to the masses, not the dimensions; we are treating these as point particles with no internal structure. This is just to say that the sun is massively more… massive than the earth, so we want to see how that works out for the physics of their motion.
Lets take a moment to review our previous results. We had a state vector of 12 quantities, \(S = [{\bf r}_0, {\bf r}_1, {\bf \dot{r}}_0, {\bf \dot{r}}_1]\), and a system of two differential equations describing the motion of the bodies:
\[ \begin{align} \ddot{\bf r}_0 &= -G m_1 \frac{{\bf r}_0 - {\bf r}_1}{r^3} \\ \ddot{\bf r}_1 &= -G m_0 \frac{{\bf r}_1 - {\bf r}_0}{r^3} \\ \end{align} \]
Note the signs chosen to show the similar form of the equations.
Now taking \({\bf r}\) as the position of the earth, \(m_1\) around the sun, \(m_0\), we have equations for its position, velocity, and acceleration:
\[ {\bf r} = {\bf r}_1 - {\bf r}_0 \\ \dot{\bf r} = \dot{\bf r}_1 - \dot{\bf r}_0 \\ \ddot{\bf r} = \ddot{\bf r}_1 - \ddot{\bf r}_0 \]
Substituting the accelerations derived from Newton’s law of gravitation, we develop the equations of motion.
\[ \begin{align} \ddot{\bf r} &= \ddot{\bf r}_1 - \ddot{\bf r}_0 \\ \ddot{\bf r} &= -G m_0 \frac{\bf r}{r^3} - G m_1 \frac{\bf r}{r^3} \\ \ddot{\bf r} &= -G(m_0 + m_1) \frac{\bf r}{r^3} \end{align} \]
we now have a single equation, with only six quantities in the state vector, \(S = [{\bf r}, {\bf \dot{r}}]\).
Writing the standard gravitational parameter, \(\mu = G(m_0 + m_1)\),
\[ \ddot{\bf r} = -\mu \frac{\bf r}{r^3} = -\mu \frac{\hat{\bf r}}{r^2} \]
So the acceleration is along the path between the two and decreases as the square of the distance. The equations are proportional to some negative constant, signalling an attractive force relating to a characteristic mass.
We want to recast this vector equation into some scalar form and see if we can tease out some of the properties of the motion.
We cross both sides of the equation by the specific angular momentum, \({\bf h}\). Note that this differs from \({\bf l}\) by a factor of the mass.
\[ \begin{align} {\bf l} &= {\bf r} \times {\bf p} \\ {\bf h} &= {\bf r} \times \dot{\bf r} = \frac{\bf l}{m} \end{align} \]
\[ \begin{align} \ddot{\bf r} \times {\bf h} &= -\left(\frac{\mu}{r^3}\right) {\bf r} \times {\bf h} \\ \frac{\mathrm{d}}{\mathrm{d} t}(\dot{\bf r} \times {\bf h}) - \dot{\bf r} \times \dot{\bf h} &= -\left(\frac{\mu}{r^3}\right) {\bf r} \times {\bf h} \\ \frac{\mathrm{d}}{\mathrm{d} t}(\dot{\bf r} \times {\bf h}) &= \mu \frac{\mathrm{d}}{\mathrm{d} t} \frac{\bf r}{r} \\ \frac{\mathrm{d}}{\mathrm{d} t}\left(\dot{\bf r} \times {\bf h} - \mu \frac{\bf r}{r} \right) &= {\bf 0} \end{align} \]
Where we have cancelled the term \(\dot{\bf r} \times \dot{\bf h}\) since the angular momentum is constant.
Integrating, we get
\[ \dot{\bf r} \times {\bf h} - \mu \frac{\bf r}{r} = \mu {\bf e} \]
where the eccentricity vector, \({\bf e}\), is a constant vector along the apse line, the line between the farthest and nearest points of the orbit. Rearranging:
\[ \frac{\bf r}{r} + {\bf e} = \frac{\dot{\bf r} \times {\bf h}}{\mu} \]
Now taking the dot product with \({\bf r}\):
\[ \begin{align} \frac{\bf r}{r} + {\bf e} &= \frac{\dot{\bf r} \times {\bf h}}{\mu} \\ \frac{\bf r}{r} \cdot {\bf r} &= \left( \frac{\dot{\bf r} \times {\bf h}}{\mu} - {\bf e} \right) \cdot {\bf r} \\ r &= \frac{{\bf r} \cdot \left( \dot{\bf r} \times {\bf h}\right)}{\mu} - {\bf r} \cdot {\bf e} \end{align} \]
Introducing the true anomaly, \(\nu\),
\[ \nu = \arccos \frac{{\bf e}\cdot{\bf r}}{er} \]
we find the orbit equation.
\[ r = \frac{h^2}{\mu}\frac{1}{1 + e cos(\nu)} \]
So here we have an equation for the radius of the motion with respect to one of the foci, and the angle from the major axis.
Let’s try to ballpark something more familiar. Using the relevant masses, eccentricities, and orbital periods, we can estimate the other parameters for Earth and Pluto.
E_earth <- -1.5501749e36 # kg m / s^2
l_earth <- 1.1e39 # kg m^2 / s
m_earth <- 5.972e24 # kg
G <- 6.67e-11 # N m^2 / kg^2
m_sun <- 1.99e30 # kg
earth_df <- data.frame(E = E_earth, l = l_earth, m = m_earth, G = G, M = m_sun)E_pluto <- -1.2405e32 # kg m / s^2
l_pluto <- 1.2140e37 # kg m^2 / s
m_pluto <- 1.303e22 # kg
G <- 6.67e-11 # N m^2 / kg^2
m_sun <- 1.99e30 # kg
pluto_df <- data.frame(E = E_pluto, l = l_pluto, m = m_pluto, G = G, M = m_sun)