When Physics Broke My Brain, I Let Pygame Handle the Chaos.

The usual coding origin starts with a desire to a new carrer, financial freedom, or just a university course. But mine started with sheer passion to not only solve equations, but to see what’s actually happening behind the scenes. That’s when I wanted to learn Python, then discovered something called Pygame, which I can actually use to simulate the things which I wanted to see.

Pygame, which is a Python module used to make simple 2D video games, or maybe even at professional setups, became my toolbox to code whatever comes to my mind, either Spring-block systems, or Elastic Collisions of balls, or Planetary motion, or THE CHAOTIC DOUBLE PENDULUM. Suddenly physics wasn’t just homework for me; it became an Interactive Project.

The code: Where Physics meets Pygame.

Before we do all those simulations, calculations and stuff, we need a platform, a screen, to show them on. That’s where Pygame shows up. This is how literally each Pygame Project ever made starts off…

It does three key things: sets up the window, defines the simulation’s speed , and creates the game loop that keeps the simulation running.

The Boilerplate:

import pygame

pygame.init()

# Setting the screen dimensions
SCREEN_WIDTH = 800
SCREEN_HEIGHT = 600
screen = pygame.display.set_mode((SCREEN_WIDTH, SCREEN_HEIGHT))
pygame.display.set_caption("Physics Simulation")

# Setting the FPS
clock = pygame.time.Clock()
FPS = 60 

running = True
while running:  # The loop where the whole simulation lives in!

    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            running = False

    screen.fill((0, 0, 0)) # Actually putting something on to the screen(In this case, it's just a black screen))

    pygame.display.flip()
    clock.tick(FPS)

pygame.quit()

The essential part of this code, the heart of this code is the while running: loop. Everytime it runs, it updates the contents(we technicallly call ’em sprites in Pygame!) on the screen, like the velocity, acceleration, colour, size or whatever you can think of.

U can even make buffaloes dance. Yep! That’s possible!!

This is where and how we transform static math into continous fluid like simulations, giving it some life and soul, which breathes at 60 frames per second(w.r.t to the code above).

Project 1: Planetary Motion.

Once the basic canvas was set up, I could see what Newton saw when apple fell on his head, yeah.. the equations of Newton’s Laws Of Motion and could replicate it in a dynamic envioronment.

When simulating it, I quickly realised that I could change the fate of smaller, lighter bodies which revolve around larger, heavier bodies just by changing their velocities. I could make its orbit elliptical instead of a perfect circle, or make it spiral into the heavier body, or just kick it out of the frame.

Here’s the code I came up with…

import pygame
import random
from sys import exit
import math

earth_x = 600
earth_y = 200

moon_x = 884
moon_y = 200

earth_velocity = 0
moon_velocity = 0

G = 6.67 * 10**1 
mass_earth = 81.3
mass_moon = 1
r_earth = 3.67 *7
r_moon = 1 *7

# All the above variables are the initial conditions of the planets. 
# For this simulation, I took heavier body as The Sun, and lighter body as The Earth.
# For their masses, I took the ratio of their masses as force between two bodies is directly proportional to mass.

pygame.init()

earth_pos = pygame.Vector2(earth_x, earth_y)
moon_pos = pygame.Vector2(moon_x, moon_y)

direction = moon_pos - earth_pos

m_velocity = math.sqrt((G * mass_earth) / pygame.Vector2.magnitude(earth_pos - moon_pos))
e_velocity = math.sqrt((G * mass_moon) / pygame.Vector2.magnitude(earth_pos - moon_pos))

earth_velocity = pygame.Vector2(0, 0)
moon_velocity = pygame.Vector2(1, 2)

earth_surf = pygame.Surface((r_earth * 2, r_earth * 2), pygame.SRCALPHA)
pygame.draw.circle(earth_surf,"Orange",(r_earth, r_earth), r_earth)
earth_rect = earth_surf.get_rect(center = earth_pos)

moon_surf = pygame.Surface((r_moon * 2, r_moon * 2), pygame.SRCALPHA)
pygame.draw.circle(moon_surf,"green",(r_moon, r_moon),r_moon)
moon_rect = moon_surf.get_rect(center = moon_pos)

trail_moon = []
trail_earth = []

screen_center = pygame.Vector2(1200/2, 700/2)

screen = pygame.display.set_mode((1200, 700))

stars_world = [pygame.Vector2(random.randint(-2000, 2000), random.randint(-2000, 2000)) for _ in range(5000)]

clock = pygame.time.Clock()

while True:
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            pygame.quit()
            exit()

    screen.fill((30, 31, 32))

    earth_rect.center = earth_pos
    moon_rect.center = moon_pos

    screen_center = pygame.Vector2(1200//2, 700//2)
    camera_offset = screen_center - earth_pos

    direction = moon_pos - earth_pos
    direction_sq = direction.magnitude_squared()

    force_earth = (G * mass_earth * mass_moon) / direction_sq
    force_moon = -force_earth

    earth_acc = (force_earth/mass_earth) * direction.normalize()
    moon_acc = (force_moon/mass_moon) * direction.normalize()

    earth_velocity += earth_acc
    moon_velocity += moon_acc

    earth_pos += earth_velocity
    moon_pos += moon_velocity

    trail_moon.append(moon_pos.copy())


    for points in trail_moon:
        pygame.draw.circle(screen, "White", (int(points[0]), int(points[1])), 1)

    screen.blit(earth_surf, earth_rect)
    screen.blit(moon_surf, moon_rect)

    pygame.display.update()
    clock.tick(60)

But when we run this simulation, the planet actually kinda drifts away from the orbit and changes orbit every revolution. That’s because I’ve used Euler’s method of Numerical methods to calculate the position of the planet every frame. Even though the method is fairly simple, the error percentage for each iteration is what makes the planet drift away. To make the simulation a bit more stable, using methods like RK-4(Runge Kutta 4th Order) helps us achieve that stabilisation we aim for.

RK-4 methods use several intermediate steps to calculate the final position for the next frame, which is crucial to maintain stability in physics simulation. This helped me to learn that both coding and mathematical knowledge is required for physics simulations.

https://drive.google.com/file/d/1TQtmjFuYklns8gwSi5iFmkVxiNZE4ODq/view?usp=drive_link

The above video shows the Planetary Motion simulation.

Project 2: The chaotic DOUBLE PENDULUM.

When the planetary motion taught me presicion, this DOUBLE PENDULUM humbled me. The math behind this made me go NUTS! Solving the Lagrangian and later realising that the solution is Indeterministic was a bit painful, but finding the rate of change in angular accelaration in other initial terms like angular velocity, angular displacement, etc… made the whole thing easier, but it was still hard.

Here’s the code for the CHAOTIC DOUBLE PENDULUM:

import math
import pygame
import random
import matplotlib.pyplot as plt
from sys import exit

theta_1 = -(3.14)
theta_2 = 3.14
omega_1 = 0
omega_2 = 0

mass_1 = 10.00
mass_2 = 100.00

g =  9.80665

energy_values = []
time_values = []
time = 0
initial_energy = 0

# All the above variables are the initial conditions.

def heart(): # This is the main function where all the math goes into and spits out the final position each frame.
    global theta_1, theta_2, omega_1, omega_2, time

    an1 = -1 *  (mass_2 * l1 * (omega_1 ** 2) * math.sin(theta_1 - theta_2) * math.cos(theta_1 - theta_2)) + (mass_2 * g * math.sin(theta_2) * math.cos(theta_1 - theta_2)) - (mass_2 * l2 * (omega_2 ** 2) * math.sin(theta_1 - theta_2)) - ((mass_1 + mass_2) * g * math.sin(theta_1)) 
    ad1 = l1 * (mass_1 + (mass_2 * (math.sin(theta_1 - theta_2 )** 2)))

    an2 = (mass_1 + mass_2) * ((l1 * (omega_1 ** 2) * math.sin(theta_1 - theta_2)) - (g * math.sin(theta_2)) + (g * math.sin(theta_1) * math.cos(theta_1 - theta_2)) + (l2 * (omega_2 ** 2) * math.sin(theta_1 - theta_2) * math.cos(theta_1 - theta_2)))
    ad2 = l2 * (mass_1 + (mass_2 * (math.sin(theta_1 - theta_2 )** 2)))
    theta_1_acc = 1 * (an1/ad1)
    theta_2_acc = 1 * (an2/ad2)
    dt = 0.01  # Timestep

    omega_1 += theta_1_acc * dt
    omega_2 += theta_2_acc * dt

    theta_1 += omega_1
    theta_2 += omega_2

    ke1 = (0.5 * mass_1 *  (l1 ** 2) * (omega_1 ** 2))
    ke2 =  0.5 * mass_2 * ((l1 * omega_1) ** 2 + (l2 * omega_2) ** 2 + 2 * l1 * l2 * omega_1 * omega_2 * math.cos(theta_1 - theta_2))

    pe1 = (mass_1 * g * -1 * (l1 * math.cos(theta_1)))
    pe2 = (mass_2 * g *( -1 * (l1 * math.cos(theta_1)) - (l2 * math.cos(theta_2))))

    total_energy = ke1 + pe1 + ke2 + pe2
    energy_values.append(total_energy)
    time_values.append(time)
    time += dt

    initial_energy = energy_values[0]
    delta_e = total_energy - initial_energy

pygame.init()

state = [theta_1, theta_2, omega_1, omega_2]

pivot_x = 600
pivot_y = 350
ball_1_x = 600
ball_1_y = 450
ball_2_x = 600
ball_2_y = 550

pivot_pos = pygame.Vector2(pivot_x, pivot_y)
ball_1_pos = pygame.Vector2(ball_1_x, ball_1_y)
ball_2_pos = pygame.Vector2(ball_2_x, ball_2_y)

l1 = (pivot_pos - ball_1_pos).magnitude()
l2 = (ball_1_pos - ball_2_pos).magnitude()

pivot_surf = pygame.Surface((30, 30), pygame.SRCALPHA)
pivot = pygame.draw.circle(pivot_surf, "Blue", (15, 15), 15)
pivot_rect = pivot_surf.get_rect(center = pivot_pos)

ball_1_surf = pygame.Surface((20, 20), pygame.SRCALPHA)
ball_1 = pygame.draw.circle(ball_1_surf, "Orange", (10, 10), 10)
ball_1_rect = ball_1_surf.get_rect(center = ball_1_pos)

ball_2_surf = pygame.Surface((20, 20), pygame.SRCALPHA)
ball_2 = pygame.draw.circle(ball_2_surf, "Yellow", (10, 10), 10)
ball_2_rect = ball_2_surf.get_rect(center = ball_2_pos)

trail = []

clock = pygame.time.Clock()

screen = pygame.display.set_mode((1200, 700))

while True:
    for event in pygame.event.get():
        if event.type == pygame.QUIT: # This is where we plot the graph for Energy Conservation.
            plt.plot(time_values, energy_values)
            plt.xlabel('Time (s)')
            plt.ylabel('Total Mechanical Energy (Joules)')
            plt.title('Total Energy of Double Pendulum vs Time')
            plt.grid(True)
            plt.show()
            plt.axhline(y=energy_values[0], color='r', linestyle='--', label='Initial Energy')
            plt.legend()

            pygame.quit()
            exit()
    screen.fill("black")

    heart()

    # The below code converts polar coordinates of the bob of the pendulum to cartesian coordinates.

    ball_1_x = l1 * math.sin(theta_1)
    ball_1_y = -1 * (l1 * math.cos(theta_1))
    ball_2_x = ball_1_x + (l2 * math.sin(theta_2))
    ball_2_y = ball_1_y - (l2 * math.cos(theta_2))

    ball_1_rect.center = (pivot_x + ball_1_x, pivot_y - ball_1_y)
    ball_2_rect.center = (pivot_x + ball_2_x, pivot_y - ball_2_y)

    trail.append([pivot_x + ball_2_x, pivot_y - ball_2_y])

    for points in trail:
        pygame.draw.circle(screen, "yellow", (int(points[0]), int(points[1])), 1)

    pygame.draw.line(screen, "white", ball_1_rect.center, ball_2_rect.center, 2)
    pygame.draw.line(screen, "white", pivot_rect.center, ball_1_rect.center, 2)

    screen.blit(pivot_surf, pivot_rect)
    screen.blit(ball_1_surf, ball_1_rect)
    screen.blit(ball_2_surf, ball_2_rect)



    pygame.display.update()
    clock.tick(100)

This time, to even go crazier, I plotted the energy conservation graph to see if it actually conserves energy, but nah… it didn’t. The reason behind is, I again used Euler’s method for running the simulation which clearly messed up the energy conservation, which is accecptable cause.. It’s Euler’s method.

The Euler’s Method introduces cumulative error over the time, causing the energy conservation to fail. RK-4 will definetely do better in this case, but yeah, I chose to stick with Euler’s method cause, why not?!

Running this code was beautiful, it showed me chaos — not as a theory, but as a beautiful simulation which depicts the unpredictable reality of a simple physical system and infinite possibilites which arise from a subtle change like 0.0001 degrees which just follows a completely new path after just a few seconds.

https://drive.google.com/file/d/1H0SYVHf-or9SLgHWvmC6cAHSlk3ML3yk/view?usp=drive_link

The above video shows the Double Pendulum Experiment.

The Gift Of Scientific Coding…

This journey of myself into this coding, python, pygame and all the stuff that I’ve talked about above taught me an important lesson, simulations aren’t about precision, it is all about the problem-solving skill and application of knowledge.

I learnt to code because I genuinely wanted to answer a simple, yet computationally complex question. Debugging such a code is cruel and frustrating, but the dopamine hit you get when planets don’t act like buffaloes dancing and actually start behaving as planets, is just pure bliss.

Coding isn’t only about coding, You can create life, end the same life, maybe just let it grow to so that it ends your CPU. You can be the GOD of your universe where YOU determine the laws of physics…