I am working on project where I want to derive the orbital equation for a satellite of Mass M, orbiting around a planet of Mass P. I want to do this with these unknowns first to then apply the orbital equation to planets in our solar system. How would I go about doing this?

  • $\begingroup$ This is a tall order, and to answer everything would need a chapter of a book or more. I wonder if you can narrow this down a bit? Also, if you can add anything that you've tried so far, or a link to what you've read would really help. That way people who would like to answer could better understand which parts to focus on. I'll leave a short answer soon, to get things started at least. $\endgroup$ – uhoh Mar 4 '19 at 0:50
  • $\begingroup$ Which "planets"/primaries and moons? On Earth-Moon, you'll get errors on the order of 1% by assuming the moon has negligible mass. On Pluto-Charon it's more like 10%. If you try to avoid this by considering the moons to be nonneglible in mass, it gets very complicated unless there is just one moon. The way to an economical and yet accurate calculation depends on which system you consider and the accuracy you want! $\endgroup$ – Oscar Lanzi Mar 4 '19 at 0:52
  • $\begingroup$ Well I've been looking at the two-body problem to see if it would apply here. In general, I am not concerned with "accuracy". The assignment I am doing is more based on exploration. Do you have any suggestions of how I could narrow it down? This is what came my mind. $\endgroup$ – Dylan Mar 4 '19 at 0:57
  • 1
    $\begingroup$ youtube.com/watch?v=Yc3l1MEvkO0 I have been watching this video, which demonstrates a two-body problem using differential equations. Could this apply here? $\endgroup$ – Dylan Mar 4 '19 at 0:58
  • $\begingroup$ @Dylan you can use two bodies if you consider only one moon as nonneglible mass. That works in most places -- but not really on Jupiter and its four similar mass Galilean moons. There you might have to consider them all negligible. $\endgroup$ – Oscar Lanzi Mar 4 '19 at 1:03

From Quick facts #2: The two-body problem:

enter image description here

Masses $m_1$ and $m_2$ orbit around their center of mass at the same angular rate $\omega$, so

$$v_1 = \omega r_1, \ \ v_2 = \omega r_2$$


$$\omega = \sqrt{\frac{G(m_1 + m_2)}{r^3}} $$


$$m_1 r_1 = m_2 r_2.$$

The tutorial then derives:

$$ r_1 = r \frac{m_2}{m_1 + m_2}, \ \ r_2 = r \frac{m_1}{m_1 + m_2},$$

where $r = r_1 + r_2$. It then states the relationship between the distance $r$ and the period $T$ as:

$$\frac{r^3}{T^2} = \frac{G(m_1+m_2)}{4 \pi^2}$$

So to get the ball rolling (to start your orbit) with a planet of mass $M$ and satellite of mass $P$ separated by a distance $R$, position them at

$$\mathbf{x_M} = +R \frac{P}{M+P} \ \mathbf{\hat{x}}$$

$$\mathbf{x_P} = -R \frac{M}{M+P} \ \mathbf{\hat{x}}$$

and set the velocities to

$$\mathbf{v_M} = +\omega R \frac{P}{M+P} \ \mathbf{\hat{y}}$$

$$\mathbf{v_P} = -\omega R \frac{M}{M+P} \ \mathbf{\hat{y}}$$

Let's try it in Python.

The acceleration of each object is related to the mass of the other:

$$\mathbf{a_1} = -m_2 G \frac{1}{r_{12}^2} \mathbf{\hat{r}}_{12} = -m_2 G \frac{\mathbf{r_{12}}}{r_{12}^3}$$

$$\mathbf{a_2} = -m_1 G \frac{1}{r_{21}^2} \mathbf{\hat{r}}_{21} = -m_2 G \frac{\mathbf{r_{21}}}{r_{21}^3}$$

enter image description here

def deriv(X, t):
    x,   v   = X.reshape(2, -1)
    x1,  x2  = x.reshape(2, -1)
    x12, x21 = x1-x2, x2-x1

    acc1     = -x12 * m2 * G * ((x12**2).sum())**-1.5
    acc2     = -x21 * m1 * G * ((x21**2).sum())**-1.5

    return np.hstack((v, acc1, acc2))

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint as ODEint

halfpi, pi, twopi, fourpi = [f*np.pi for f in (0.5, 1, 2, 4)]

m1, m2 = pi, 42
R      = 5.
G      = 3.

# http://radio.astro.gla.ac.uk/a1dynamics/twobody.pdf
omega = np.sqrt(G * (m1 + m2) / R**3)
T     = np.sqrt(4 * pi**2 * R**3 / (G * (m1 + m2)))

print ("double check, are they equal? ", T, twopi/omega)

x1 = R * m2 / (m1 + m2)
x2 = R * m1 / (m1 + m2)

v1 = omega * R * m2 / (m1 + m2)
v2 = omega * R * m1 / (m1 + m2)

X0    = np.array([x1, 0] + [-x2, 0] + [0, v1] + [0, -v2])  # note the negative signs!
times = np.linspace(0, T, 201)

answer, info = ODEint(deriv, X0, times, full_output=True)
print answer.shape

pos1, pos2, vel1, vel2 = answer.T.reshape(4, 2, -1)
if True:
    for x, y in (pos1, pos2):
        plt.plot(x, y)
        plt.plot(x[:1], y[:1], 'ok')
  • $\begingroup$ In the video I linked above, differential equations are used to obtain an equation for motion at the end. The resource you posted does not include anything similar. Why is this? $\endgroup$ – Dylan Mar 4 '19 at 2:05
  • $\begingroup$ @Dylan re-read what you wrote in your two sentence question, then re-read my comment to it. If there is something specific in the video you need explained, then ask a new question, specify the time in the video (or show a screen shot) and ask something much more specific! You can ask as many questions as you like, but please write your questions carefully and explain better what it is that you are asking. $\endgroup$ – uhoh Mar 4 '19 at 2:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.