A simple educational particle simulator that captures laws of classical mechanics. In the current state, it will plot an exaggerated perihelion precision of Mercury.

SimpleParticleSimulator.py

# Author: Miladiouss (Milad Pourrahmani)
# Author's Website: Miladiouss.com
# Date: Sep 26, 2018
# Title: Simple Particle Simulator
# Description:
"""
I am hoping to assemble a series of brief educational computer programs that will capture the simplicity and elegance of laws of physics, both in nature and in apprehension.

This particular educational project (Simple Particle Simulator) is the first of this series. In less than 25 lines, plotting and comments aside, we can obtain the trajectory of a particle in ANY force field. The generated plot displays the simulated and exaggerated perihelion precision of Mercury. 

In future series, I will demonstrate how one can use this program to develop similar insights that are gained from most classical mechanics textbooks.

Unrelated, this project was developed on an Android device using the Pydroid 3 application. For some reason, what we can do with our portable devices still amazes me and makes me wonder what would Newton, Faraday, and Schrodinger do if they had access to the same technology.
"""

# Motivation:
"""
If the giants of physics somehow could stand on the shoulders of the computers, the formulation of physics would have been much different; it would have been computational rather than analytical. Since computers became accessible three centuries after major developments in physics, physics and applied analytical methods are considered synonymous as the result. For this reason, almost all educational resources in physics are dedicated to analytical methods. Aside from one course in numerical methods, computational methods are considered extracurricular and are not an integral part of all and every physics course. It is hard to deny the deep insights gained from an analytical approach, so is denying its limitations and arduousness. Perhaps numerical methods deserve as much if not more attention as analytical methods. All students regardless of their field of study, deserve to be exposed to computational thinking, especially in sciences.
"""

# Imports
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# In a Jupyter notebook, uncomment the following:
# matplotlib inline



# ========== Write a function to plot trajectory of a particle (IGNORE) ==========

def plot(r_vecs, scale = 1.5):
    """
    Creates a pretty rotating plot given a list of coordinates in the [..., (x, y, z, t), ...] format
    """

    # Convert coords to a NumPy array
    r_vecs = np.array(r_vecs)
    
    # Separate Coordinates
    xs, ys, zs, ts = r_vecs[:, 0], r_vecs[:, 1], r_vecs[:, 2], r_vecs[:, -1]

    # Calculate distance r from the origin
    rs = np.linalg.norm(r_vecs[:,:-1], axis=1)
    
    # For 2D plot, use: plt.plot(xs, ys)

    # Setup matplotlib fig and ax (for more flexibility)
    fig = plt.figure()
    fig.set_facecolor('black')
    ax = plt.axes(projection='3d')
    
    # Plot the trajectory
    ax.plot3D(xs, ys, zs, 'white')
    
    # Put a star at the origin
    ax.plot3D([0], [0], [0], "*")
    
    # Mark the r0 with a red dot
    # xs[:1] instead of xs[0] to create a list, i.e. [xs[0]], not just xs[0]
    ax.plot3D(xs[:1], ys[:1], zs[:1], "ro")
    
    # Mark the end with a blue dot
    ax.plot3D(xs[-1:], ys[-1:], zs[-1:], "bo")
    
    # Force all axis to have the same scale/range
    ax.axis('scaled')
    ax.set_xlim3d(-scale, scale)
    ax.set_ylim3d(-scale, scale)
    ax.set_zlim3d(-scale, scale)
    
    # Black rules!
    ax.set_facecolor('black')
    ax.w_xaxis.set_pane_color((0, 0, 0, 1.))
    ax.w_yaxis.set_pane_color((0, 0, 0, 1.))
    ax.w_zaxis.set_pane_color((0, 0, 0, 1.))
    
    # But axis should be white
    [t.set_color('white') for t in ax.xaxis.get_ticklabels()]
    [t.set_color('white') for t in ax.yaxis.get_ticklabels()]
    [t.set_color('white') for t in ax.zaxis.get_ticklabels()]
    
    # rotate the axes and update
    for angle in range(0, 90):
        ax.view_init(30, 4*angle)
        plt.draw()
        plt.pause(.05)



# ========== Create a particle class/object ==========

class Particle():

    def __init__(self, r0, v0, m0 = 1):
        """
        Position Notation: (r_x, r_y, r_z,   t) as NumPy array
        Velocity Notation: (v_x, v_y, v_z, v_t) as NumPy array
        """
        # r_vecs and v_vecs are lists of positions and velocities 
        # This is to keep track of the motion of the particle for plotting
        self.r_vecs = [r0]
        self.v_vecs = [v0]
        self.m      = [m0] # So you can do rocket problems :D


    def move(self, force_field, dt = 0.001):
        """
        Moves the particle one step
        """

        # 1. Update the position
        # Velocity is the rate of change of position (v = dr / dt)
        self.r_vecs.append( self.r_vecs[-1] + self.v_vecs[-1] * dt )
        # 2. Find acceleration at the new location
        a = force_field(self.r_vecs[-1]) / self.m[-1]
        # 3. Update the velocity at the new location 
        # Acceleration is the rate of change of velocity (a = dv / dt)
        self.v_vecs.append( self.v_vecs[-1] + a * dt )



# ========== Definition of a force field ==========

def force_field(r_vec, v_vec = None):
    """
    r_vec: (x, y, z, t)
    v_vec: Not implemented (for your exploration)
    Note: r_vec should be r_vecs[-1]
    Return: force at r_vec
    """

    # Cartesian Coords
    x, y, z, t = r_vec
    xyz_vec = np.array([x, y, z])

    # Spherical Coords
    r     = np.linalg.norm(xyz_vec)
    polar = np.arccos(z/r)
    azim  = np.arctan(y/x)

    # Newtons law:
    # f_vec = -1 * xyz_vec / r**(3)
    # It's r**3 since xyz_vec is not a unit vector and has length r

    # Use 3.05 instead of 3 to simulate and exaggerate the perihelion precession of Mercury
    # ||||||||||||||||||||||||||||||||||||
    f_vec = -1.5* xyz_vec / r**(3.1) # ||
    # ||||||||||||||||||||||||||||||||||||
    
    # Make f to have the form (fx, fy, fz, t), just like r and v
    # by appending "force on time", i.e. 0
    f_vec = np.append(f_vec, 0.)
        
    return f_vec

# ========== Run the simulation ==========

# Create a particle
# Notation: (r_x, r_y, r_z, t) for position and (v_x, v_y, v_z, v_t) for velocity, v_t should be 1.
p = Particle(
        r0 = np.array((1, 0, 0, 0)),
        v0 = np.array((0.5, 0.9998, 0.25, 1)), # 0.9998 for 1/r**5 cases, to demonstrate instability
        m0 = 1)

# Move the particle for certain number of steps
for step in range(5500):
    p.move(force_field, dt = 0.01)

# Plot
plot(p.r_vecs)