Created
January 13, 2019 19:45
-
-
Save timakro/fd03263da485af7fd5a1a95b108ea429 to your computer and use it in GitHub Desktop.
Caculate Pi using a physics simulation of elastic collisions https://www.youtube.com/watch?v=HEfHFsfGXjs
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#!/usr/bin/python3 | |
import math | |
DIGITS = 4 | |
class Body: | |
def __init__(self, pos, vel, mass): | |
self.pos = float(pos) | |
self.vel = float(vel) | |
self.mass = float(mass) | |
def advance(self, time): | |
self.pos += self.vel * time | |
def get_velocity_after_collision_with(self, body): | |
if self.mass == math.inf: | |
return self.vel | |
if body.mass == math.inf: | |
return -self.vel | |
return ((self.mass * self.vel + body.mass * (2 * body.vel - self.vel)) | |
/ (self.mass + body.mass)) | |
def __str__(self): | |
return "pos={:.3f};vel={:.3f};mass={:.1f}".format( | |
self.pos, self.vel, self.mass) | |
a = Body(1, 0, 1) | |
b = Body(2, -1, 100**DIGITS) | |
wall = Body(0, 0, math.inf) | |
def time_until_collision(left_body, right_body): | |
vel_diff = left_body.vel - right_body.vel | |
if vel_diff > 0: # bodies moving towards each other | |
return (right_body.pos - left_body.pos) / vel_diff | |
return math.inf | |
def update_velocity_after_collision(body1, body2): | |
vel1 = body1.get_velocity_after_collision_with(body2) | |
vel2 = body2.get_velocity_after_collision_with(body1) | |
body1.vel = vel1 | |
body2.vel = vel2 | |
total_time = 0 | |
total_collisions = 0 | |
while True: | |
time_wall_a = time_until_collision(wall, a) | |
time_a_b = time_until_collision(a, b) | |
if time_wall_a == math.inf and time_a_b == math.inf: | |
break | |
if time_a_b < time_wall_a: | |
col_type = "a-b" # collision between body a and b happens next | |
time = time_a_b | |
else: | |
col_type = "w-a" # collision between wall and body a happens next | |
time = time_wall_a | |
total_time += time | |
total_collisions += 1 | |
a.advance(time) | |
b.advance(time) | |
wall.advance(time) | |
if col_type == "a-b": | |
update_velocity_after_collision(a, b) | |
elif col_type == "w-a": | |
update_velocity_after_collision(wall, a) | |
if DIGITS <= 4: | |
print("{} num={} time={:.3f} a({}) b({}) w({})".format( | |
col_type, total_collisions, total_time, a, b, wall)) | |
print("Pi is {}".format(total_collisions/10**DIGITS)) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment