[Leejohnston]

Code:

# orbital_simulator.py
# The Lumin Archive — Computer Science Section
# CC BY 4.0 | Written and Compiled by Lee Johnston

import matplotlib.pyplot as plt
import numpy as np

G = 6.674e-11  # gravitational constant

class Body:
    def __init__(self, name, mass, position, velocity, color):
        self.name = name
        self.mass = mass
        self.position = np.array(position, dtype=float)
        self.velocity = np.array(velocity, dtype=float)
        self.color = color
        self.trajectory = [self.position.copy()]

def gravitational_force(b1, b2):
    r_vec = b2.position - b1.position
    r = np.linalg.norm(r_vec)
    if r == 0:
        return np.zeros(2)
    return G * b1.mass * b2.mass * r_vec / r**3

def update_positions(bodies, dt):
    forces = {b: np.zeros(2) for b in bodies}
    for i, b1 in enumerate(bodies):
        for j, b2 in enumerate(bodies):
            if i != j:
                forces[b1] += gravitational_force(b1, b2)

    for b in bodies:
        a = forces[b] / b.mass
        b.velocity += a * dt
        b.position += b.velocity * dt
        b.trajectory.append(b.position.copy())

def simulate(bodies, steps=5000, dt=60):
    for _ in range(steps):
        update_positions(bodies, dt)

def plot_orbits(bodies):
    plt.style.use('dark_background')
    for b in bodies:
        traj = np.array(b.trajectory)
        plt.plot(traj[:,0], traj[:,1], color=b.color, label=b.name)
        plt.scatter(traj[-1,0], traj[-1,1], color=b.color, s=30)
    plt.xlabel("x position (m)")
    plt.ylabel("y position (m)")
    plt.legend()
    plt.axis('equal')
    plt.title("Orbital Simulator — The Lumin Archive")
    plt.show()

if __name__ == "__main__":
    sun = Body("Sun", 1.989e30, [0, 0], [0, 0], 'yellow')
    earth = Body("Earth", 5.972e24, [1.496e11, 0], [0, 29_780], 'cyan')
    mars = Body("Mars", 6.417e23, [2.279e11, 0], [0, 24_077], 'orange')

    bodies = [sun, earth, mars]
    simulate(bodies, steps=3000, dt=3600)
    plot_orbits(bodies)