from dataclasses import dataclass from math import pi, tan, atan2, sqrt import numpy as np import matplotlib.pyplot as plt @dataclass class Line: x: float y: float angle: float def __repr__(self): return f"Line(x={self.x}, y={self.y}, angle={self.angle*180/pi})" def bounce_wall(line): if line.angle > pi: raise ValueError intercept_x = line.x - line.y / tan(line.angle) angle = -line.angle + 2*pi return Line(x=intercept_x, y=0, angle=angle) def bounce_circle(line): x0 = line.x t = tan(line.angle) t2 = t*t a = 1 + t2 b = 4*t - 2*t2*x0 c = t2*x0*x0 - 4*x0*t + 3 rad = b*b - 4*a*c if rad < 0: raise ValueError x1 = (-b + sqrt(rad)) / (2*a) y1 = (x1 - x0)*t x2 = (-b - sqrt(rad)) / (2*a) y2 = (x2 - x0)*t x, y = (x1, y1) if y1 > y2 else (x2, y2) t_tan = (atan2(y+2, x) - pi/2) % (2*pi) angle = (line.angle + 2*(t_tan - line.angle)) % (2*pi) return Line(x=x, y=y, angle=angle) def do_bounces(angle_deg, x=None): if x is not None: angle_deg = np.arctan2(2, x + 2) * 180 / pi lines = [Line(-2, -2, angle_deg*pi/180)] try: while True: lines.append(bounce_wall(lines[-1])) lines.append(bounce_circle(lines[-1])) except ValueError: pass return len(lines)-1, lines def draw_bounces(lines): xs = [line.x for line in lines] ys = [line.y for line in lines] plt.plot([-4, 4], [0, 0]) plt.plot(xs, ys) plt.gca().add_patch(plt.Circle((0, -2), radius=1)) plt.xlim((-1, 1)) plt.ylim((-1.8, 0.2)) plt.show() def count_range(start, end, count=100): xs = np.linspace(start, end, count) angles = np.arctan2(2, xs + 2) * 180 / pi bounces = [] for angle in angles: count, _ = do_bounces(angle) bounces.append(count) plt.plot(xs, bounces) plt.show() return xs, np.array(bounces) def iterative_zoom(start, end, n_iter): maxes = [] for _ in range(n_iter): xs, bounces = count_range(start, end, 1000) max_x = xs[np.argmax(bounces)] maxes.append(max_x) range_ = end - start start = max_x - 0.05*range_ end = max_x + 0.05*range_ print(maxes) plt.plot(maxes) plt.show() def main(): # count_range(-2, 0, 10000) # count_range(-0.824, -0.822, 10000) # count_range(-0.8225, -0.8223, 10000) # count_range(58, 62) # count_range(59, 60) # count_range(59.51, 59.515) # draw_bounces(do_bounces(None, x=-0.8224863249433897)[1]) iterative_zoom(-2, 0, 20) if __name__ == '__main__': main()