Skip to content

Instantly share code, notes, and snippets.

@thquinn

thquinn/sagemath_juggling.py

Last active June 2, 2024 07:24
Show Gist options
  • Save thquinn/39c140d572dbc8d67eed702f3a8af836 to your computer and use it in GitHub Desktop.
Save thquinn/39c140d572dbc8d67eed702f3a8af836 to your computer and use it in GitHub Desktop.
Download ZIP
A not-very-fast symbolic physics simulation using SymPy.
Raw
sagemath_juggling.py
import os
import sage.all as sage
from itertools import combinations
from PIL import Image, ImageDraw
from typing import cast, Dict
from enum import Enum
os.system("clear")
class SimCircle:
center: sage.vector
radius: sage.Expression
velocity: sage.vector
def __init__(self, center: sage.vector, radius: sage.Expression, velocity: sage.vector = sage.vector([0, 0])):
self.center = center
self.radius = radius
self.velocity = velocity
self.first = False
def to_tuple(self):
return sage.vector([self.center[0], self.center[1], self.velocity[0], self.velocity[1]])
def __str__(self):
return f'SimCircle center: {self.center}, radius: {self.radius}, velocity: {self.velocity}'
class SimSegment:
p1: sage.vector
p2: sage.vector
def __init__(self, p1: sage.vector, p2: sage.vector):
self.p1 = p1
self.p2 = p2
def length(self):
return sage.sqrt((self.p2[0] - self.p1[0]) ** 2 + (self.p2[1] - self.p1[1]) ** 2)
def contains_box(self, p):
xs = sorted([self.p1[0], self.p2[0]])
ys = sorted([self.p1[1], self.p2[1]])
return bool(p[0] >= xs[0]) and bool(p[0] <= xs[1]) and bool(p[1] >= ys[0]) and bool(p[1] <= ys[1])
def projection(self, p):
AB = self.p2 - self.p1
AC = p - self.p1
AD = AB * (AB.dot_product(AC)) / AB.dot_product(AB)
return self.p1 + AD
def __str__(self):
return f'SimSegment p1: {self.p1}, p2: {self.p2}'
class SimState(Enum):
INCOMPLETE = 1
COMPLETE_WITH_PERIOD = 2
INVALID_MULTIPLE_COLLISIONS = 3
INVALID_SPAWN_OVERLAP = 4
INVALID_TIMEOUT = 5
class Sim:
SPAWN_RATE: sage.Expression = sage.Integer(1)
TIMEOUT: sage.Expression = SPAWN_RATE * sage.Integer(50)
SPAWN_HEIGHT: sage.Expression = sage.Integer(3)
DESPAWN_HEIGHT: sage.Expression = sage.Integer(-10)
CIRCLE_RADIUS: sage.Expression = sage.Rational('1/3')
GRAVITY: sage.Expression = sage.Integer(-3)
SEGMENT_LENGTH: sage.Expression = sage.Rational('3/2')
SEGMENT_SEPARATION: sage.Expression = sage.Integer(3)
COEFFICIENT_OF_RESTITUTION: sage.Expression = sage.Rational('1/1')
FRAME_TIME = sage.Rational('1/10')
t: sage.Expression
circles: list[SimCircle]
segments: list[SimSegment]
next_spawn_t: sage.Expression
previous_configurations: Dict[list, sage.Expression]
state: SimState
def __init__(self, theta1:sage.Expression, theta2:sage.Expression):
self.t = sage.Integer(0)
self.circles = []
left_segment_center = sage.vector([0, 0])
right_segment_center = sage.vector([Sim.SEGMENT_SEPARATION, 0])
self.segments = [
SimSegment(left_segment_center - sage.vector([sage.cos(theta1), -sage.sin(theta1)]) * Sim.SEGMENT_LENGTH / 2, left_segment_center + sage.vector([sage.cos(theta1), -sage.sin(theta1)]) * Sim.SEGMENT_LENGTH / 2),
SimSegment(right_segment_center - sage.vector([sage.cos(theta2), sage.sin(theta2)]) * Sim.SEGMENT_LENGTH / 2, right_segment_center + sage.vector([sage.cos(theta2), sage.sin(theta2)]) * Sim.SEGMENT_LENGTH / 2)
]
self.next_spawn_t = sage.Integer(0)
self.previous_configurations = dict()
self.state = SimState.INCOMPLETE
def iterate(self):
collisions = []
for circle_pair in combinations(self.circles, 2):
next_collision = self.t + circle_collision_next_time(*circle_pair)
collisions.append([next_collision, *circle_pair])
for segment in self.segments:
for circle in self.circles:
next_collision = self.t + circle_segment_collision_next_time(circle, segment, Sim.GRAVITY)
collisions.append([next_collision, circle, segment])
min_time_collisions = []
min_time: sage.Expression = sage.infinity
for collision in collisions:
if collision[0] < min_time:
min_time_collisions = []
min_time = collision[0]
if collision[0] <= min_time:
min_time_collisions.append(collision)
min_time = min_time.full_simplify() if min_time != sage.oo else sage.oo
print(f'min collision time: {min_time.n()}, next spawn: {self.next_spawn_t.n()}')
if min_time >= self.next_spawn_t:
self.advance_t_to(self.next_spawn_t)
self.next_spawn_t += Sim.SPAWN_RATE
# TODO: check for circles in spawn area
self.circles.append(SimCircle(sage.vector([0, Sim.SPAWN_HEIGHT]), Sim.CIRCLE_RADIUS))
if (len(self.circles) == 1):
self.circles[0].first = True
return -1
# TODO: check for circles involved in multiple collisions
self.advance_t_to(min_time)
for collision in min_time_collisions:
c1 = cast(SimCircle, collision[1])
if isinstance(collision[2], SimSegment):
seg = cast(SimSegment, collision[2])
perform_circle_point_collision(c1, seg.projection(c1.center), Sim.COEFFICIENT_OF_RESTITUTION)
else:
c2 = cast(SimCircle, collision[2])
perform_circle_collision(c1, c2, Sim.COEFFICIENT_OF_RESTITUTION)
self.circles = [c for c in self.circles if c.center[1] > Sim.DESPAWN_HEIGHT]
print('circles in the sim:', len(self.circles))
configuration = [c.to_tuple() for c in self.circles]
configuration.sort()
configuration = str(configuration)
if configuration in self.previous_configurations:
prev_config_time = self.previous_configurations[configuration]
config_delta = self.t - prev_config_time
self.state = SimState.COMPLETE_WITH_PERIOD
return config_delta / Sim.SPAWN_RATE
self.previous_configurations[configuration] = self.t
if self.t > Sim.TIMEOUT:
self.state = SimState.INVALID_TIMEOUT
return -1
def advance_t_to(self, next_t: sage.Expression):
if (self.t == 0):
self.save_img();
print(f'advancing to t={next_t.n()}')
next_frame = sage.ceil(self.t / Sim.FRAME_TIME) * Sim.FRAME_TIME
while (next_frame < next_t):
self.advance_t_to_impl(next_frame)
next_frame += Sim.FRAME_TIME
self.save_img();
self.advance_t_to_impl(next_t)
def advance_t_to_impl(self, next_t: sage.Expression):
delta: sage.Expression = next_t - self.t
if delta == 0:
return
for circle in self.circles:
circle.center = sage.vector([
circle.center[0] + circle.velocity[0] * delta,
circle.center[1] + circle.velocity[1] * delta + Sim.GRAVITY * delta ** 2 / 2
])
circle.velocity = sage.vector([
circle.velocity[0],
circle.velocity[1] + delta * Sim.GRAVITY
])
self.t = next_t
def save_img(self):
pixels_per_unit = 100
bounds_min = (-3, -3)
bounds_max = (6, 5)
def coor_convert(p: sage.vector):
px = (p[0].n() - bounds_min[0]) * pixels_per_unit
py = (bounds_max[1] - p[1].n()) * pixels_per_unit
return (px, py)
img = Image.new("RGBA", ((bounds_max[0] - bounds_min[0]) * pixels_per_unit, (bounds_max[1] - bounds_min[1]) * pixels_per_unit))
draw = ImageDraw.Draw(img)
draw.rectangle([(0,0), img.size], fill = (0, 0, 0) )
for segment in self.segments:
draw.line([*coor_convert(segment.p1), *coor_convert(segment.p2)], width=5, fill=(255, 255, 255))
for circle in self.circles:
draw.ellipse(
[
*coor_convert(circle.center - sage.vector([circle.radius, -circle.radius])),
*coor_convert(circle.center + sage.vector([circle.radius, -circle.radius]))
],
fill=(255, 255, 255)
)
img.save(f'img{self.t.n()}.png', 'PNG')
def circle_collision_next_time(c1: SimCircle, c2: SimCircle):
sage.forget()
print('circle collision next time in')
t = sage.var('t')
sage.assume(t > 0)
sage.assume(t, 'real')
eq = sage.sqrt( (c1.velocity[0] * t + c1.center[0] - c2.velocity[0] * t - c2.center[0]) ** 2 + (c1.velocity[1] * t + c1.center[1] - c2.velocity[1] * t - c2.center[1]) ** 2) == c1.radius + c2.radius
collision_times = eq.solve(t, to_poly_solve=True)
print('circle collision next time out')
return min([c.rhs() for c in collision_times]) if len(collision_times) > 0 else sage.oo
def circle_segment_collision_next_time(c: SimCircle, s: SimSegment, g: sage.Expression):
sage.forget()
if c.first:
print('circle-segment collision next time in')
print(c, s, g)
t = sage.var('t')
# sage.assume(t > 0)
# sage.assume(t, 'real')
p = sage.vector([c.center[0] + c.velocity[0] * t, c.center[1] + c.velocity[1] * t + g * t ** 2 / 2])
# Find collisions with the line defined by the segment.
eq = sage.abs_symbolic((s.p2[0] - s.p1[0]) * (p[1] - s.p1[1]) - (p[0] - s.p1[0]) * (s.p2[1] - s.p1[1])) / s.length() == c.radius
if c.first:
print('line equation', eq)
collision_times = eq.solve(t, to_poly_solve=True)
#####
# ^ TODO -- BUG! ^
#####
# to_poly_solve is missing solutions. e.g. the following:
# t = var('t')
# f = 2/3*abs(-1/245760000*sqrt(2)*((4*(5*sqrt(-32*sqrt(2) + 288) - 72)^2 - sqrt(2)*(sqrt(2)*((5*sqrt(-32*sqrt(2) + 288) - 72)^2 + 900*sqrt(2) + 720*sqrt(-32*sqrt(2) + 288) - 12384) + 1800) + 3600*sqrt(2) + 2880*sqrt(-32*sqrt(2) + 288) - 49536)*(sqrt(2)*(sqrt(2)*((5*sqrt(-32*sqrt(2) + 288) - 72)^2 + 900*sqrt(2) + 720*sqrt(-32*sqrt(2) + 288) - 12384) + 1800) - 3600*sqrt(2))*t*sqrt(-32*sqrt(2) + 288)/(abs(-1/2400*(5*sqrt(-32*sqrt(2) + 288) - 72)^2 + 1/9600*sqrt(2)*(sqrt(2)*((5*sqrt(-32*sqrt(2) + 288) - 72)^2 + 900*sqrt(2) + 720*sqrt(-32*sqrt(2) + 288) - 12384) + 1800) - 3/8*sqrt(2) - 3/10*sqrt(-32*sqrt(2) + 288) + 129/25)^2 + abs(-1/9600*sqrt(2)*(sqrt(2)*((5*sqrt(-32*sqrt(2) + 288) - 72)^2 + 900*sqrt(2) + 720*sqrt(-32*sqrt(2) + 288) - 12384) + 1800) + 3/8*sqrt(2))^2) + 69120000*sqrt(2) - 552960000) + 1/245760000*sqrt(2)*(((4*(5*sqrt(-32*sqrt(2) + 288) - 72)^2 - sqrt(2)*(sqrt(2)*((5*sqrt(-32*sqrt(2) + 288) - 72)^2 + 900*sqrt(2) + 720*sqrt(-32*sqrt(2) + 288) - 12384) + 1800) + 3600*sqrt(2) + 2880*sqrt(-32*sqrt(2) + 288) - 49536)^2*sqrt(-32*sqrt(2) + 288)/(abs(-1/2400*(5*sqrt(-32*sqrt(2) + 288) - 72)^2 + 1/9600*sqrt(2)*(sqrt(2)*((5*sqrt(-32*sqrt(2) + 288) - 72)^2 + 900*sqrt(2) + 720*sqrt(-32*sqrt(2) + 288) - 12384) + 1800) - 3/8*sqrt(2) - 3/10*sqrt(-32*sqrt(2) + 288) + 129/25)^2 + abs(-1/9600*sqrt(2)*(sqrt(2)*((5*sqrt(-32*sqrt(2) + 288) - 72)^2 + 900*sqrt(2) + 720*sqrt(-32*sqrt(2) + 288) - 12384) + 1800) + 3/8*sqrt(2))^2) - 46080000*sqrt(-32*sqrt(2) + 288))*t - 276480000*t^2 - 76800*(5*sqrt(-32*sqrt(2) + 288) - 72)^2 + 69120000*sqrt(2) - 55296000*sqrt(-32*sqrt(2) + 288) + 951091200)) == (1/3)
# solve(f, t, to_poly_solve=True)
# test on https://sagecell.sagemath.org/
# with to_poly_solve, returns [] — without, returns a simplified equation:
# 1/16*abs(18*sqrt(2)*t^2 + 3*sqrt(2)*t*sqrt(-32*sqrt(2) + 288) - 36*sqrt(2) - 8) == (1/2)
# that Wolfram can find real solutions to. adding assumptions doesn't help
#####
# maybe even worse, if you simplify the equation and square it on both sides you get:
# t = var('t')
# f = 1/256*(18*sqrt(2)*t^2 + 3*sqrt(2)*t*sqrt(-32*sqrt(2) + 288) - 36*sqrt(2) - 8)**2 == (1/4)
# s = solve(f, t, to_poly_solve=True)
# [s.rhs().n() for s in s]
# which returns four complex solutions. but if you plug the same equation into Wolfram, you get back the same 4 real solutions as the non-squared version, the smallest positive one is the one we want.
# does Sage sometimes return WRONG solutions, and exclude solutions even when to_poly_solve is False?
#####
# Filter out collisions where the collision point is not on the segment.
collision_times = [t_val for t_val in collision_times if s.contains_box(s.projection(p.subs(t = t_val.rhs())))]
if c.first:
print('collision times', collision_times)
# Add collisions with the endpoints themselves. This may include times where the circle contains part/all of the segment.
# We don't really care about that, though: as long as we always process the earliest event, we shouldn't be letting
# circles pass through segments to begin with.
for endpoint in [s.p1, s.p2]:
# TODO: can we square both sides of this kind of equation?
eq = sage.sqrt((endpoint[0] - p[0]) ** 2 + (endpoint[1] - p[1]) ** 2) == c.radius
if c.first:
print('endpoint equation', eq)
endpoint_solutions = eq.solve(t, to_poly_solve=True)
collision_times.extend(endpoint_solutions)
if c.first:
print('collision times w/ endpoints', collision_times)
collision_times = [c for c in collision_times if c.rhs() in sage.RR and c.rhs() > sage.Integer(0)]
if c.first:
print('collision times after filter', collision_times)
print('circle-segment collision next time out')
return min([c.rhs() for c in collision_times]) if len(collision_times) > 0 else sage.oo
def perform_circle_collision(c1: SimCircle, c2: SimCircle, e: sage.Expression):
dv = c2.velocity - c1.velocity
normal = (c2.center - c1.center).normalized()
vn = dv[0] * normal[0] + dv[1] * normal[1]
j = (1 + e) * vn / (1 / c1.radius + 1 / c2.radius)
c1.velocity += j * normal / c1.radius
c2.velocity -= j * normal / c2.radius
def perform_circle_point_collision(c: SimCircle, p: sage.vector, e: sage.Expression):
normal = (c.center - p).normalized()
vn = c.velocity[0] * normal[0] + c.velocity[1] * normal[1]
c.velocity -= normal * vn * (1 + e)
# c1 = SimCircle(sage.vector([0, 0]), 1, sage.vector([0, 0]))
# c2 = SimCircle(sage.vector([5, 0]), 1, sage.vector([-1, 0]))
# print(circle_collision_next_time(c1, c2))
# s = SimSegment(sage.vector([-2, -5]), sage.vector([20, -6]))
# print(circle_segment_collision_next_time(c1, s, -1))
sim = Sim(sage.pi / 4, sage.pi / 4)
while sim.state == SimState.INCOMPLETE:
print(f'iterating: {sim.t}')
period = sim.iterate()
if (sim.state == SimState.COMPLETE_WITH_PERIOD):
print(f'we actually found a period: {period}')
print(sim.state)
Raw
sympy_juggling.py
# rewrite in SymEngine? https://github.com/symengine/symengine
# seems to not have an analogue to solve()
# or: https://www.sagemath.org/
# we might be able to specify a range for t in sp.solve to reduce solve times?
# algo:
# - start at t=0 with two line segments and an empty list of circles, and repeat the following
# - find collision times between all pairs of circles and all circle/segment pairs
# - filter out all collisions that don't take place at the min time
# - if the min time is less than or equal to the next spawn time:
# - if any one circle is involved in multiple collisions, the result is INVALID: MULTIPLE COLLISIONS
# - process all collisions
# - advance time to the min collision time
# - if the min time is greater than the next spawn time:
# - advance time to the next spawn time
# - remove all circles with y < some negative number
# - if we are at a spawn time, spawn a new circle
# - if the new circle overlaps an existing circle, the result is INVALID: SPAWN OVERLAP
# - check if we've seen this configuration of circles before
# - if we have, return (elapsed time since we last saw this configuration) / (spawn period)
# - otherwise, record the time and configuration
# - if the time is over some maximum, the result is INVALID: TIMEOUT
import sympy as sp
from itertools import combinations
from PIL import Image, ImageDraw
from typing import cast, Dict
from enum import Enum
import os
os.system("clear")
class SimCircle:
circle: sp.Circle
velocity: sp.Point
def __init__(self, circle: sp.Circle):
self.circle = circle
self.velocity = sp.Point2D(0, 0)
def to_tuple(self):
return [self.circle.center.x, self.circle.center.y, self.velocity.x, self.velocity.y]
def __str__(self) -> str:
return f'SimCircle {self.circle}, {self.velocity}'
class SimState(Enum):
INCOMPLETE = 1
COMPLETE_WITH_PERIOD = 2
INVALID_MULTIPLE_COLLISIONS = 3
INVALID_SPAWN_OVERLAP = 4
INVALID_TIMEOUT = 5
class Sim:
SPAWN_RATE: sp.Number = 1
TIMEOUT: sp.Number = SPAWN_RATE * 50
SPAWN_HEIGHT: sp.Number = 3
DESPAWN_HEIGHT: sp.Number = -10
CIRCLE_RADIUS: sp.Number = sp.Rational(1, 3)
GRAVITY: sp.Number = -3
SEGMENT_LENGTH: sp.Number = sp.Rational(3, 2)
SEGMENT_SEPARATION: sp.Number = 3
COEFFICIENT_OF_RESTITUTION: sp.Number = sp.Rational(9, 10)
FRAME_TIME = sp.Rational(1, 10)
t: sp.Number
circles: list[SimCircle]
segments: list[sp.Segment2D]
next_spawn_t: sp.Number
previous_configurations: Dict[list, sp.Number]
state: SimState
def __init__(self, theta1:sp.Number, theta2:sp.Number):
self.t = sp.S.Zero
self.circles = []
left_segment_center = sp.Point2D(0, 0)
right_segment_center = sp.Point2D(Sim.SEGMENT_SEPARATION, 0)
self.segments = [
sp.Segment2D(left_segment_center - sp.Point2D(sp.cos(theta1), -sp.sin(theta1)) * Sim.SEGMENT_LENGTH / 2, left_segment_center + sp.Point2D(sp.cos(theta1), -sp.sin(theta1)) * Sim.SEGMENT_LENGTH / 2),
sp.Segment2D(right_segment_center - sp.Point2D(sp.cos(theta2), sp.sin(theta2)) * Sim.SEGMENT_LENGTH / 2, right_segment_center + sp.Point2D(sp.cos(theta2), sp.sin(theta2)) * Sim.SEGMENT_LENGTH / 2)
]
self.next_spawn_t = sp.S(0)
self.previous_configurations = dict()
self.state = SimState.INCOMPLETE
def iterate(self):
collisions = []
for circle_pair in combinations(self.circles, 2):
next_collision = self.t + sp.Min(*circle_collision_times(*circle_pair))
collisions.append([next_collision, *circle_pair])
for segment in self.segments:
for circle in self.circles:
next_collision = self.t + sp.Min(*circle_line_segment_collision_times(circle, segment, Sim.GRAVITY))
collisions.append([next_collision, circle, segment])
min_time_collisions = []
min_time: sp.Number = sp.oo
for collision in collisions:
if collision[0] < min_time:
min_time_collisions = []
min_time = collision[0]
if collision[0] <= min_time:
min_time_collisions.append(collision)
print(f'min collision time: {min_time.evalf()}, next spawn: {self.next_spawn_t.evalf()}')
if min_time >= self.next_spawn_t:
self.advance_t_to(self.next_spawn_t)
self.next_spawn_t += Sim.SPAWN_RATE
# TODO: check for circles in spawn area
self.circles.append(SimCircle(sp.Circle(sp.Point2D(0, Sim.SPAWN_HEIGHT), Sim.CIRCLE_RADIUS)))
return -1
# TODO: check for circles involved in multiple collisions
self.advance_t_to(min_time)
for collision in min_time_collisions:
c1 = cast(SimCircle, collision[1])
if isinstance(collision[2], sp.Segment2D):
seg = cast(sp.Segment2D, collision[2])
perform_circle_point_collision(c1, seg.projection(c1.circle.center), Sim.COEFFICIENT_OF_RESTITUTION)
else:
c2 = cast(SimCircle, collision[2])
perform_circle_collision(c1, c2, Sim.COEFFICIENT_OF_RESTITUTION)
self.circles = [c for c in self.circles if c.circle.center.y > Sim.DESPAWN_HEIGHT]
configuration = [c.to_tuple() for c in self.circles]
configuration.sort()
configuration = str(configuration)
if configuration in self.previous_configurations:
prev_config_time = self.previous_configurations[configuration]
config_delta = prev_config_time - self.t
self.state = SimState.COMPLETE_WITH_PERIOD
return config_delta / Sim.SPAWN_RATE
self.previous_configurations[configuration] = self.t
if self.t > Sim.TIMEOUT:
self.state = SimState.INVALID_TIMEOUT
return -1
def advance_t_to(self, next_t: sp.Number):
if (self.t == 0):
self.save_img();
print(f'advancing to t={next_t.evalf()}')
next_frame = sp.ceiling(self.t / Sim.FRAME_TIME) * Sim.FRAME_TIME
while (next_frame < next_t):
self.advance_t_to_impl(next_frame)
next_frame += Sim.FRAME_TIME
self.save_img();
self.advance_t_to_impl(next_t)
def advance_t_to_impl(self, next_t: sp.Number):
next_t = sp.simplify(next_t)
delta: sp.Number = next_t - self.t
if delta.is_zero:
return
for circle in self.circles:
circle.circle = sp.Circle(
sp.Point2D(sp.simplify(circle.circle.center.x + circle.velocity.x * delta),
sp.simplify(circle.circle.center.y + circle.velocity.y * delta + Sim.GRAVITY * delta ** 2 / 2)),
circle.circle.radius
)
circle.velocity = sp.Point2D(
circle.velocity.x,
sp.simplify(circle.velocity.y + delta * Sim.GRAVITY),
)
self.t = next_t
def save_img(self):
pixels_per_unit = 100
bounds_min = (-3, -3)
bounds_max = (6, 5)
def coor_convert(p: sp.Point2D):
px = (p.x.evalf() - bounds_min[0]) * pixels_per_unit
py = (bounds_max[1] - p.y.evalf()) * pixels_per_unit
return (px, py)
img = Image.new("RGBA", ((bounds_max[0] - bounds_min[0]) * pixels_per_unit, (bounds_max[1] - bounds_min[1]) * pixels_per_unit))
draw = ImageDraw.Draw(img)
draw.rectangle([(0,0), img.size], fill = (0, 0, 0) )
for segment in self.segments:
draw.line([*coor_convert(segment.p1), *coor_convert(segment.p2)], width=5, fill=(255, 255, 255))
for circle in self.circles:
draw.ellipse([*coor_convert(circle.circle.center - sp.Point2D(circle.circle.radius, -circle.circle.radius)),
*coor_convert(circle.circle.center + sp.Point2D(circle.circle.radius, -circle.circle.radius))],
fill=(255, 255, 255))
img.save(f'img{self.t.evalf()}.png', 'PNG')
def circle_collision_times(c1: SimCircle, c2: SimCircle):
t = sp.symbols('t', real=True, positive=True)
eq = sp.Eq(sp.sqrt( (c1.velocity.x * t + c1.circle.center.x - c2.velocity.x * t - c2.circle.center.x) ** 2 + (c1.velocity.y * t + c1.circle.center.y - c2.velocity.y * t - c2.circle.center.y) ** 2), c1.circle.radius + c2.circle.radius)
collision_times = sp.solve(eq, t)
collision_times = [time for time in collision_times if time.is_real and time > 0]
return collision_times
def circle_line_segment_collision_times(c: SimCircle, s: sp.Segment2D, g: sp.Number):
print('=== circle_line_segment_collision_times ===')
print(c)
print(s)
t = sp.symbols('t', real=True, positive=True)
p = sp.Point2D(c.circle.center.x + c.velocity.x * t, c.circle.center.y + c.velocity.y * t + g * t ** 2 / 2)
# Find collisions with the line defined by the segment.
line = sp.Line2D(s)
eq = sp.Eq(line.distance(p), c.circle.radius)
print('colliding with line')
collision_times = list(sp.solve(eq, t))
# Filter out collisions where the collision point is not on the segment.
collision_times = [t_val for t_val in collision_times if s.contains(line.projection(p.subs(t, t_val)))]
# Add collisions with the endpoints themselves. This may include times where the circle contains part/all of the segment.
# We don't really care about that, though: as long as we always process the earliest event, we shouldn't be letting
# circles pass through segments to begin with.
for endpoint in [s.p1, s.p2]:
eq = sp.Eq(sp.simplify(endpoint.distance(p)), c.circle.radius)
print(f'colliding with endpoint: {eq}')
endpoint_solutions = sp.solve(eq, t, manual=True)
print(f'got back: {endpoint_solutions}')
collision_times.extend(endpoint_solutions)
collision_times = [time for time in collision_times if time.is_real and time > 0]
print('=== done ===')
return collision_times
def perform_circle_collision(c1: SimCircle, c2: SimCircle, e: sp.Number):
dv = c2.velocity - c1.velocity
normal = (c2.circle.center - c1.circle.center).unit
vn = dv.x * normal.x + dv.y * normal.y
j = (1 + e) * vn / (1 / c1.circle.radius + 1 / c2.circle.radius)
c1.velocity += j * normal / c1.circle.radius
c2.velocity -= j * normal / c2.circle.radius
def perform_circle_point_collision(c: SimCircle, p: sp.Point2D, e: sp.Number):
normal = (c.circle.center - p).unit
vn = c.velocity.x * normal.x + c.velocity.y * normal.y
c.velocity -= normal * vn * (1 + e)
sim = Sim(sp.pi / 4, sp.pi / 4)
while sim.state == SimState.INCOMPLETE:
print(f'iterating: {sim.t}')
period = sim.iterate()
if (period > 0):
print(f'we actually found a period: {period}')
print(sim.state)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment