Ask
CoCalc Logo Icon
StoreFeaturesDocsShareSupport News Sign UpSign In
Views: 27
Visibility: Unlisted (only visible to those who know the link)
Image: ubuntu2004
Embed | Download | Raw |
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 x\mathbf{x} and v\mathbf{v} are determined from the previous x\mathbf{x}, v\mathbf{v}, and the calculated a\mathbf{a}.

v0=given\mathbf{v_0} = \text{given}x0=given\mathbf{x_0} = \text{given}a0=Gmx^x02\mathbf{a_0} = - G m \frac{\mathbf{\hat x}}{|\mathbf{x_0}|^2}\vdotsvn+1=vn+anΔt\mathbf{v_{n+1}} = \mathbf{v_n} + \mathbf{a_n} \Delta txn+1=xn+vn+1Δt\mathbf{x_{n+1}} = \mathbf{x_n} + \mathbf{v_{n+1}} \Delta tan+1=Gmxn+1^xn+12\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 x\mathbf{x} and v\mathbf{v}, and Newton's Universal Gravitation inverse square law for a\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()