Solar System Simulation

I built this as an A-Level Computer Science project — a 3D Newtonian gravity simulator, written from scratch in Python and PyGame, designed to be used as a teaching tool by GCSE and A-Level physics students. It renders the Sun and the eight planets as shaded 3D meshes, integrates their motion under gravity in real time, and lets a student drag sliders to change a planet's mass or the Sun's mass and watch the orbits respond. There's no 3D or physics library underneath any of it: the vector maths, the renderer, the depth sorting and the GUI are all mine.

The physics

Every frame, each planet's acceleration is recalculated from Newton's law of gravitation and integrated forward:

I deliberately do not model gravity between planets — only Sun-to-planet. In our solar system that cross-term is negligible, and skipping it turns an $O(n^2)$ force calculation into $O(n)$, which matters once meteors are in the picture.

Getting the scale right

The solar system does not fit on a screen at a consistent scale. If the Sun is drawn 50 pixels across, Earth — at 1/109th the Sun's diameter — should be $0.46$ pixels wide, and would vanish. Distances are worse: on a linear scale that fits Neptune on screen, the inner four planets collapse into a single point. I scaled both radius and orbital distance logarithmically instead: a planet's real size grows exponentially outward from Mercury to Jupiter, so a log scale turns that into a linear, readable progression on screen. It's not physically honest — Earth and Venus end up closer together on screen than they should be, and at times planets visually overlap when in reality they're millions of kilometres apart — but it's the only version of this simulation a student can actually look at.

A 3D engine with no 3D library

PyGame draws 2D primitives and nothing else, so every planet is a hand-built triangle mesh: five rings of manually-placed points — pole, tropic, equator, tropic, pole — joined into faces. I settled on 36 triangles per planet after testing denser meshes generated by an icosphere subdivision algorithm; they looked better close up but cost more to sort and shade for a difference I decided wasn't worth it at this viewing distance, and increasing the face count later is an $O(n)$ change if I ever need it.

Hiding the far side of each planet — back-face culling — is a single dot product: for every triangle, if its outward normal points away from the camera, don't draw it. I tested the technique on a plain cube before trusting it on a planet.

The first real bug came from where I'd put the Sun: at the origin, $(0,0,0)$ — which is also where the camera's projection maths measures distances from. The unit-vector function divided a zero vector by its own (zero) magnitude, and rather than raising an error, it silently produced a vector of $\mathrm{NaN}$s. Every triangle that touched that calculation stopped drawing. Nothing crashed; the Sun and its light source just quietly disappeared, along with anything whose shading depended on it.

def calcUnitVector(vector):
    magnitude = calcMagnitude(vector)
    if magnitude != 0:
        return vector[0]/magnitude, vector[1]/magnitude, vector[2]/magnitude
    else:
        return 0, 0, 0

Saturn's rings broke the same back-face culling I'd just proven worked. A ring isn't convex — from some angles you can see the underside of the far edge and the topside of the near edge at once — so treating it as ordinary faces on the planet mesh always drew it wrong from certain angles.

Depth, without a depth buffer

There's no z-buffer here — planets are simply sorted by distance from the camera every frame and drawn back-to-front, so a nearer planet is painted over a farther one. Python's built-in sort() (Timsort, $O(n\log n)$) handles the ordering; a lambda pulls the distance out of each planet as the sort key. It's the classic painter's algorithm, and it's the same trick used to fix Saturn's rings above.

Meteors and the inverse-square law

The Meteor class reuses the Planet physics almost unchanged, with one difference: a planet only feels the Sun, but a meteor sums the pull of every body in the system before it moves. Each meteor also carries a visible line pointing along its current acceleration vector — not scaled to the force's magnitude (the range is too extreme to draw sensibly at both a planet's surface and the edge of the system), just its direction, so a student can see the field, not only the motion it produces.

To make that pull visible, the simulation lets you increase a planet's mass live. Radius and density sliders combine as $m \propto r^{3}$, so a 10× radius increase is a 1000× mass increase — enough to visibly bend meteor paths toward whichever planet you've inflated.

With enough meteors, some eventually drift far enough from the system that their acceleration lines stretch across the entire screen, which looks like a rendering bug even though the physics is correct. I clean these up with a Manhattan-distance check — cheaper than the true Euclidean distance since it needs no square root, and I only need to know a meteor is far away, not exactly how far.

for meteor in meteors:
    meteor.accelerate(planetsToDraw)
    meteor.move()
    meteorPosition = meteor.position
    meteor.rotateAndScale(totalxrotation, totalyrotation)
    meteor.draw()
    if abs(meteorPosition[0]) + abs(meteorPosition[1]) + abs(meteorPosition[2]) > 10**13:
        meteors.remove(meteor)
        del meteor

Validating against NASA data

To check the physics was actually right and not just plausible-looking, I seeded the simulation with Mercury's real position, velocity and mass from NASA's planetary fact sheet, ran it for three orbits, and compared the measured orbital period and mean distance against the published values. My threshold for a pass was within 5%.

96.6%

orbital period accuracy

0.233 yr measured vs. 0.241 yr actual, averaged over 3 orbits of Mercury

97.9%

orbital distance accuracy

56.66M km measured mean distance vs. 57.9M km actual

±0.8%

spread from float rounding

measured by seeding many Mercury-like orbits at random starting angles

169 orbits

vs. 168.37 expected

after 168.37 virtual years — about 30 real minutes of runtime

The remaining error traces to three places: the initial conditions use averaged rather than instantaneous position and velocity, acceleration is only recalculated 100 times a second rather than continuously, and 64-bit floats round differently depending on a planet's starting angle. I isolated that last one by seeding dozens of Mercury-like orbits at random angles around the Sun and plotting them — if the physics were exact they'd land on exactly the same point on a $T^2$ vs. $R^3$ graph; the small scatter that appears instead is rounding error, not a modelling mistake.

Kepler's third law says $T^2 \propto R^3$ for every planet orbiting the same star, with the constant of proportionality $\dfrac{GM_\odot}{4\pi^2}$. Plotting the simulation's own output for all eight planets is the real end-to-end test: it doesn't just check one planet against one fact sheet, it checks whether the whole system's dynamics are self-consistent. On a linear scale the four inner planets collapse into an unreadable cluster in the corner, so I plotted it again on a log-log scale.

The graph above only proves Kepler's law — it doesn't let a student ask their own question of the data. So every run also writes each planet's live position, velocity and orbital data out to an Excel spreadsheet via Openpyxl, overwriting the previous run rather than leaving a folder of near-duplicate files behind. A student who wants to test something the built-in graph plotter doesn't cover — checking $F=ma$ by hand, say — can just open the numbers themselves.

Performance

The frame rate is capped at 100 FPS. With just the Sun and eight planets it holds that cap comfortably; the cost that actually matters is sorting and shading planets live as the user adds more.

100 FPS

baseline, capped

Sun + 8 planets, steady

~40 FPS

with 50 user-added planets

depth-sorted and shaded every frame

~20 FPS

with 500 meteors

meteors skip depth-sorting — too small for draw order to matter

7.5%

CPU load at 100 FPS

measured in Activity Monitor on the development machine

Designing it with the people who'd use it

Before writing any code I interviewed my A-level physics teacher, a Year 9 student about to start GCSE physics, and an A-level classmate applying to study planetary science. Astrophysics is the one part of school physics where students never get to run their own experiment — there's no lab kit for it — and that's the gap this was built to fill: not a demo to watch, but a system a student could perturb and question. The classmate's own frustration, not being able to visualise Lagrange points from a textbook diagram, was part of what pushed me toward letting the user change any variable live rather than only pressing play on a fixed scenario.