Orbital Motion.ipynb

Kernel: Python 3 (system-wide)

Orbital Motion Model

Goals:

• Use simple linear approximation to plot the orbits of bodies subject to gravitational force.

• Object oriented programming methods

  • including matplotlib

• Although this worksheet uses 3D vectors, it all works in 2D as well.

• Use physically reasonable physical parameters for solar system.

The model is an iterative function, where each $\mathbf{x}$ and $\mathbf{v}$ are determined from the previous $\mathbf{x}$, $\mathbf{v}$, and the calculated $\mathbf{a}$.

$$\mathbf{v_0} = \text{given}$$$$\mathbf{x_0} = \text{given}$$$$\mathbf{a_0} = - G m \frac{\mathbf{\hat x}}{|\mathbf{x_0}|^2}$$$$\vdots$$$$\mathbf{v_{n+1}} = \mathbf{v_n} + \mathbf{a_n} \Delta t$$$$\mathbf{x_{n+1}} = \mathbf{x_n} + \mathbf{v_{n+1}} \Delta t$$$$\mathbf{a_{n+1}} = - G m \frac{\mathbf{\hat{x_{n+1}}}}{|\mathbf{x_{n+1}}|^2}$$

The method is linear approximation, using kinematics equations for $\mathbf{x}$ and $\mathbf{v}$, and Newton's Universal Gravitation inverse square law for $\mathbf{a}$.

The vectors are represented in Python using NumPy arrays. NumPy is preferred for ease of manipulation as well as plotting via matplotlib.

import numpy as np
import matplotlib.pyplot as plt

from mpl_toolkits import mplot3d

# constant of Universal Gravitation
#   all units SI
G = 6.67430e-11

class Planet(object):
    """ A planet is a body subject to gravitational force, i.e. with mass, position, and velocity. """
    # Mass is a real number
    # Position, represented with an x, is a three-dimensional np array, i.e. a vector
    # Velocity, represented with a v, is a three-dimensional np array, i.e. a vector

    def __init__(self, mass, x, v):
        """ Sets the mass, position and velocity values """
        self.mass = mass
        self.x = x
        self.v = v

    def __str__(self):
        """ Returns a string representation of self """
        return "<" + str(self.mass) + " , " + str(self.x) + " , " + str(self.v) + ">"

    def distance(self, other):
        """ Returns the distance between two bodies """
        return np.linalg.norm(self.x - other.x)

    def unit(self, other):
        """ Returns a unit vector (np.array) pointing from self to other """
        return (other.x - self.x)/self.distance(other)

    def accel(self, other):
        """ Returns the acceleration of self due to other, gravitational constant = 1 """
        return G * other.mass /(self.distance(other))**2 * self.unit(other)

# some planet data - mass (kg), position (m), velocity (m/s)
sun = {"mass": 1.989e30, "position": np.array([0, 0, 0]), "velocity": np.array([0, 0, 0])}
mercury = {"mass": 3.285e23, "position":np.array([0, 57.9e9, 0]), "velocity": np.array([47850, 0, 0])}
venus = {"mass": 4.867e24, "position": np.array([0, 108.2e9, 0]), "velocity": np.array([35000, 0, 0])}
earth = {"mass": 5.972e24, "position": np.array([149.60e9, 0, 0]), "velocity": np.array([0, 29850, 0])}
mars = {"mass": 0.642e24, "position": np.array([227.9e9, 0, 0]), "velocity": np.array([0, 24120, 0])}
jupiter = {"mass": 1898e24, "position": np.array([0, 778.6e9, 0, 0]), "velocity": np.array([-13070, 0, 0])}
saturn = {"mass": 568e24, "position": np.array([-1433.5e9, 0, 0]), "velocity": np.array([-9690, 0, 0])}
comet = {"mass": 2.2e14, "position": np.array([-85.5e9, 0, 0]), "velocity": np.array([0, -54000, 10])}

# create objects, set initial conditions
p0 = Planet(sun['mass'], sun['position'], sun['velocity'])
p1 = Planet(earth['mass'], earth['position'], earth['velocity'])
p2 = Planet(venus['mass'], venus['position'], venus['velocity'])
p3 = Planet(comet['mass'], comet['position'], comet['velocity'])

# time increment
dt = 10000

# keep track of where we've been
p1_path = np.array([p1.x])
p2_path = np.array([p2.x])
p3_path = np.array([p3.x])

# orbit simulation loop
for i in range(2000):
    # p1 - Venus
    a1 = p1.accel(p0)
    p1.v = p1.v + a1 * dt
    p1.x = p1.x + p1.v * dt

    # p2 - Earth
    a2 = p2.accel(p0)
    p2.v = p2.v + a2 * dt
    p2.x = p2.x + p2.v * dt

    # p3 - Comet
    a3 = p3.accel(p0)
    p3.v = p3.v + a3 * dt
    p3.x = p3.x + p3.v * dt

    if i % 10 == 0: # plot every n_th point
        p1_path = np.append(p1_path, [p1.x], axis=0)
        p2_path = np.append(p2_path, [p2.x], axis=0)
        p3_path = np.append(p3_path, [p3.x], axis=0)

# comment out this line or use alternate style for light background or printing
# plt.style.use('dark_background')

# 3D plot
fig = plt.figure(figsize=(8, 8))
ax = plt.axes(projection='3d')

ax.plot3D([0], [0], [0], 'ro')
ax.scatter3D(p1_path[:,0], p1_path[:,1], p1_path[:,2], s=1, color='blue')
ax.scatter3D(p2_path[:,0], p2_path[:,1], p2_path[:,2], s=1, color='green')
ax.scatter3D(p3_path[:,0], p3_path[:,1], p3_path[:,2], s=1, color='red')

plt.show()

# create objects, set initial conditions for an imaginary system with inbteractions
p0 = Planet(sun['mass'], sun['position'], sun['velocity'])
p1 = Planet(venus['mass'], venus['position'], venus['velocity'])
p2 = Planet(earth['mass'], earth['position'], earth['velocity'])

p3 = Planet(earth['mass']/10, earth['position'] * 1.05, earth['velocity'] * 0.9)

# time increment
dt = 10000

# keep track of where we've been
p1_path = np.array([p1.x])
p2_path = np.array([p2.x])
p3_path = np.array([p3.x])

# orbit simulation loop
for i in range(4000):
    # p1 - Venus
    a1 = p1.accel(p0)
    p1.v = p1.v + a1 * dt
    p1.x = p1.x + p1.v * dt

    # p2 - Earth
    a2 = p2.accel(p0) + p2.accel(p3)
    p2.v = p2.v + a2 * dt
    p2.x = p2.x + p2.v * dt

    # p3 - co-planet
    a3 = p3.accel(p0) + p3.accel(p2)
    p3.v = p3.v + a3 * dt
    p3.x = p3.x + p3.v * dt

    if i % 10 == 0: # plot every n_th point
        p1_path = np.append(p1_path, [p1.x], axis=0)
        p2_path = np.append(p2_path, [p2.x], axis=0)
        p3_path = np.append(p3_path, [p3.x], axis=0)

# comment out this line or use alternate style for light background or printing
# plt.style.use('dark_background')

# 3D plot
fig = plt.figure(figsize=(8, 8))
ax = plt.axes(projection="3d")

ax.plot3D([0], [0], [0], 'ro')
ax.scatter3D(p1_path[:,0], p1_path[:,1], p1_path[:,2], s=1, color='blue')
ax.scatter3D(p2_path[:,0], p2_path[:,1], p2_path[:,2], s=1, color='green')
ax.scatter3D(p3_path[:,0], p3_path[:,1], p3_path[:,2], s=1, color='red')

plt.show()