#!/usr/bin/env python3 """ Two-body problem — animation with an interactive mass-ratio slider. Two point masses interact through Newtonian gravity only. The centre of mass is fixed at the origin (equal and opposite momenta at t = 0), so both bodies trace similar conics about it, scaled by m2/M and m1/M respectively. A note on the "both bodies start stationary" case ------------------------------------------------- If the two bodies are genuinely released from rest, the angular momentum of the system is exactly zero and the motion is *rectilinear*: they fall straight at each other and collide. There is no orbit. To show the intended picture — bodies approaching each other and then swinging around the common centre of mass — the pair needs a non-zero initial transverse velocity. The second slider therefore sets f = (initial relative speed) / (circular-orbit speed at that separation), with the bodies starting at apoapsis, so that f = 0 both bodies at rest → head-on radial infall (the literal "stationary" case; the run stops when the surfaces touch), 0 < f < 1 ellipse of eccentricity e = 1 - f² (default), f = 1 circular orbit. Units are non-dimensional: G = 1, m1 + m2 = 1, initial separation = 1. Because the total mass and the initial separation are held fixed, the orbital period is independent of the mass ratio — only the *share* of the motion taken by each body changes. Usage ----- python two-body-animation.py # interactive, with sliders python two-body-animation.py --ratio 10 --speed 0.7 python two-body-animation.py --save img/two-body-animation.gif """ from __future__ import annotations import argparse import os import matplotlib.animation as animation import matplotlib.pyplot as plt import numpy as np from matplotlib.lines import Line2D from matplotlib.patches import Circle from matplotlib.widgets import Button, Slider from scipy.integrate import solve_ivp # ── Physical setup (non-dimensional) ────────────────────────────────────────── G = 1.0 # gravitational constant M_TOT = 1.0 # total mass m1 + m2, held fixed as the ratio is varied D0 = 1.0 # initial separation |R2 - R1| R_UNIT = 0.045 # plotted radius of a body of mass M_TOT (bodies scale as m^(1/3)) N_ORBITS = 2.0 # how many periods to integrate N_SAMPLES = 640 # animation frames per run (uniform in time) HOLD = 40 # frames the final state is held before the loop restarts Q_DEFAULT = 3.0 # m1 / m2 F_DEFAULT = 0.55 # → eccentricity 0.70 # ── Display ─────────────────────────────────────────────────────────────────── BG, FG, DIM, EDGE = '#0d1117', '#c9d1d9', '#8b949e', '#30363d' C1, C2 = '#e6194b', '#4363d8' # body 1 (m1), body 2 (m2) C_VEL, C_FOR = '#3cb44b', '#7e22ce' # velocity, gravitational force C_COM = '#f0b429' # centre of mass TRAIL_LEN = 320 # ── Dynamics ────────────────────────────────────────────────────────────────── def simulate(q, f, n_orbits=N_ORBITS, n_samples=N_SAMPLES): """Integrate the two-body problem for mass ratio ``q`` = m1/m2. The two-body problem reduces to a one-body problem in the relative coordinate r = r12 = R2 - R1 (the convention of the slides), which obeys r'' = -G M r / |r|³ with M = m1 + m2. The individual orbits follow from R1 = -(m2/M) r, R2 = (m1/M) r, which places the centre of mass at the origin for all time. Note the crossed weights: body 1 is scaled by m2, so the heavier body traces the smaller conic. """ m1 = M_TOT * q / (1.0 + q) m2 = M_TOT / (1.0 + q) R1 = R_UNIT * (m1 / M_TOT) ** (1.0 / 3.0) R2 = R_UNIT * (m2 / M_TOT) ** (1.0 / 3.0) mu = G * M_TOT # gravitational parameter of the relative orbit v_c = np.sqrt(mu / D0) # speed of a circular orbit of radius D0 e = abs(1.0 - f * f) # eccentricity (the start is apoapsis) a = D0 / (2.0 - f * f) # semi-major axis, from the vis-viva equation T = 2.0 * np.pi * np.sqrt(a ** 3 / mu) def rhs(t, y): x, yy, vx, vy = y r3 = (x * x + yy * yy) ** 1.5 return [vx, vy, -mu * x / r3, -mu * yy / r3] def contact(t, y): # stop cleanly if the surfaces touch return np.hypot(y[0], y[1]) - (R1 + R2) contact.terminal = True contact.direction = -1 t_end = n_orbits * T sol = solve_ivp( rhs, (0.0, t_end), [-D0, 0.0, 0.0, -f * v_c], t_eval=np.linspace(0.0, t_end, n_samples), events=contact, method='DOP853', rtol=1e-10, atol=1e-12, ) r, v = sol.y[:2], sol.y[2:] # relative position and velocity w1, w2 = -m2 / M_TOT, m1 / M_TOT d = np.hypot(r[0], r[1]) return dict( t=sol.t, d=d, F=G * m1 * m2 / d ** 2, p1=w1 * r, p2=w2 * r, v1=w1 * v, v2=w2 * v, m1=m1, m2=m2, R1=R1, R2=R2, q=q, f=f, e=e, a=a, T=T, collided=bool(sol.t_events[0].size), ) def view_limits(sims, pad=0.12, min_aspect=0.30): """Common axis limits covering every trajectory in ``sims``. The box hugs the trajectories (the axes keep an equal aspect ratio, so the drawn frame simply changes shape); ``min_aspect`` stops the frame from degenerating for the rectilinear case, where the motion is one-dimensional. """ xs = np.concatenate([np.concatenate([s['p1'][0], s['p2'][0]]) for s in sims]) ys = np.concatenate([np.concatenate([s['p1'][1], s['p2'][1]]) for s in sims]) rb = max(max(s['R1'], s['R2']) for s in sims) cx, cy = 0.5 * (xs.max() + xs.min()), 0.5 * (ys.max() + ys.min()) hx = (1.0 + pad) * (0.5 * (xs.max() - xs.min()) + rb) hy = (1.0 + pad) * (0.5 * (ys.max() - ys.min()) + rb) hx, hy = max(hx, min_aspect * hy), max(hy, min_aspect * hx) return (cx - hx, cx + hx), (cy - hy, cy + hy) def legend_handles(): """Proxy artists describing what is drawn, for a figure-level legend.""" return [ Line2D([], [], ls='none', marker='o', color=C1, ms=8, label='$m_1$'), Line2D([], [], ls='none', marker='o', color=C2, ms=8, label='$m_2$'), Line2D([], [], ls='none', marker='x', color=C_COM, ms=8, mew=2, label='centre of mass'), Line2D([], [], color=C_VEL, lw=2, label='velocity ($\\propto v$)'), Line2D([], [], color=C_FOR, lw=2, label='gravity ($\\propto F^{1/3}$)'), Line2D([], [], color=DIM, lw=1, ls=(0, (4, 3)), label='separation $d$'), ] # ── One animated view ───────────────────────────────────────────────────────── class Panel: """A single axes showing one two-body system.""" def __init__(self, ax): self.ax = ax self.sim = None ax.set_facecolor(BG) ax.set_aspect('equal') ax.tick_params(colors=DIM, labelsize=8) for sp in ax.spines.values(): sp.set_edgecolor(EDGE) self.trail1, = ax.plot([], [], '-', color=C1, lw=1.2, alpha=0.55, zorder=2) self.trail2, = ax.plot([], [], '-', color=C2, lw=1.2, alpha=0.55, zorder=2) self.sep, = ax.plot([], [], ls=(0, (4, 3)), color=DIM, lw=1.0, zorder=3) self.body1 = Circle((0, 0), R_UNIT, facecolor=C1, edgecolor=FG, lw=0.6, zorder=6) self.body2 = Circle((0, 0), R_UNIT, facecolor=C2, edgecolor=FG, lw=0.6, zorder=6) ax.add_patch(self.body1) ax.add_patch(self.body2) # centre of mass — fixed at the origin for all time ax.plot([0], [0], marker='x', color=C_COM, ms=9, mew=2.0, ls='none', zorder=7) arrow = dict(angles='xy', scale_units='xy', scale=1.0, width=0.006, zorder=5, capstyle='round') z = np.zeros(2) self.qvel = ax.quiver(z, z, z, z, color=C_VEL, **arrow) self.qfor = ax.quiver(z, z, z, z, color=C_FOR, **arrow) self.txt_info = ax.text(0.025, 0.975, '', transform=ax.transAxes, color=FG, fontsize=9, va='top', ha='left', family='monospace', zorder=8) self.txt_d = ax.text(0, 0, '', color=DIM, fontsize=8.5, ha='center', va='bottom', family='monospace', zorder=8) @property def n(self): return len(self.sim['t']) def set_sim(self, sim, limits=None): """Attach a trajectory and rescale the view and the arrow lengths.""" self.sim = sim self.body1.set_radius(sim['R1']) self.body2.set_radius(sim['R2']) (x0, x1), (y0, y1) = limits if limits else view_limits([sim]) self.ax.set_xlim(x0, x1) self.ax.set_ylim(y0, y1) # Velocity arrows are drawn to a linear scale, so Kepler's second law is # visible directly. The force varies as 1/d², far too wide a range to # draw linearly, so its length is compressed by a cube root. span = 0.5 * max(x1 - x0, y1 - y0) v_max = max(np.hypot(*sim['v1']).max(), np.hypot(*sim['v2']).max()) self.k_vel = 0.40 * span / v_max self.len_for = 0.26 * span self.off_for = 0.045 * span # keeps the force arrows off the d-line self.dy_lbl = 0.025 * (y1 - y0) self.F_max = sim['F'].max() self.ax.set_title( f"$m_1/m_2 = {sim['q']:.1f}$ ($m_1 = {sim['m1']:.3f}$, " f"$m_2 = {sim['m2']:.3f}$)", color=FG, fontsize=10, pad=8, ) self.draw(0) def draw(self, i): s = self.sim p1, p2 = s['p1'][:, i], s['p2'][:, i] d = s['d'][i] self.body1.center = tuple(p1) self.body2.center = tuple(p2) j = max(0, i - TRAIL_LEN) self.trail1.set_data(s['p1'][0, j:i + 1], s['p1'][1, j:i + 1]) self.trail2.set_data(s['p2'][0, j:i + 1], s['p2'][1, j:i + 1]) self.sep.set_data([p1[0], p2[0]], [p1[1], p2[1]]) self.qvel.set_offsets(np.array([p1, p2])) self.qvel.set_UVC(self.k_vel * np.array([s['v1'][0, i], s['v2'][0, i]]), self.k_vel * np.array([s['v1'][1, i], s['v2'][1, i]])) # equal and opposite (Newton's third law), each pointing at the other body; # drawn just off the line joining them so they stay legible u = (p2 - p1) / d nrm = np.array([-u[1], u[0]]) * self.off_for ell = self.len_for * (s['F'][i] / self.F_max) ** (1.0 / 3.0) self.qfor.set_offsets(np.array([p1 + nrm, p2 + nrm])) self.qfor.set_UVC(ell * np.array([u[0], -u[0]]), ell * np.array([u[1], -u[1]])) mid = 0.5 * (p1 + p2) self.txt_d.set_position((mid[0], mid[1] + self.dy_lbl)) self.txt_d.set_text(f'd = {d:.3f}') tail = '' if s['f'] == 0.0: tail = '\nzero angular momentum: radial infall' if s['collided'] and i == self.n - 1: tail = '\ncontact — the bodies collide' self.txt_info.set_text( f"t = {s['t'][i]:7.3f} (T = {s['T']:.3f})\n" f"e = {s['e']:.3f}{tail}" ) # ── Interactive mode ────────────────────────────────────────────────────────── def run_interactive(q0, f0): fig = plt.figure(figsize=(10.0, 8.0)) fig.patch.set_facecolor(BG) fig.suptitle('Two-Body Problem — motion about the common centre of mass', color=FG, fontsize=13, fontweight='bold', y=0.975) fig.legend(handles=legend_handles(), loc='upper center', bbox_to_anchor=(0.5, 0.945), ncol=6, fontsize=8.5, frameon=False, labelcolor=FG, columnspacing=1.6, handletextpad=0.6) ax = fig.add_axes([0.07, 0.23, 0.86, 0.62]) panel = Panel(ax) panel.set_sim(simulate(q0, f0)) def slider_axes(y): a = fig.add_axes([0.17, y, 0.55, 0.028], facecolor=EDGE) for sp in a.spines.values(): sp.set_edgecolor(EDGE) return a # mass ratio is swept logarithmically: the interesting range spans decades s_q = Slider(slider_axes(0.135), '$m_1/m_2$', 0.0, 2.0, valinit=np.log10(q0), color=C1, initcolor=DIM) s_f = Slider(slider_axes(0.075), '$v_0 / v_{circ}$', 0.0, 1.0, valinit=f0, color=C_VEL, initcolor=DIM) for s in (s_q, s_f): s.label.set_color(FG) s.label.set_fontsize(11) s.valtext.set_color(FG) s.valtext.set_family('monospace') s.track.set_facecolor(EDGE) # the unfilled part of the bar s_q.valtext.set_text(f'{q0:.1f}') state = dict(i=0, paused=False) def on_change(_): q = 10.0 ** s_q.val s_q.valtext.set_text(f'{q:.1f}') panel.set_sim(simulate(q, s_f.val)) state['i'] = 0 fig.canvas.draw_idle() s_q.on_changed(on_change) s_f.on_changed(on_change) def button(x, label, cb): a = fig.add_axes([x, 0.075, 0.08, 0.045]) b = Button(a, label, color=EDGE, hovercolor='#484f58') b.label.set_color(FG) b.on_clicked(cb) return b def toggle(_): state['paused'] = not state['paused'] b_play.label.set_text('Play' if state['paused'] else 'Pause') b_play = button(0.815, 'Pause', toggle) b_reset = button(0.900, 'Restart', lambda _: state.update(i=0)) fig._widgets = (s_q, s_f, b_play, b_reset) # widgets die if unreferenced fig.text(0.17, 0.012, 'G = 1, $m_1 + m_2 = 1$, initial separation = 1 — the period is ' 'therefore the same for every mass ratio;\nonly the share of the motion ' 'taken by each body changes. $v_0/v_{circ} = 0$ releases both bodies ' 'from rest.', color=DIM, fontsize=8.5, va='bottom') def tick(_): if not state['paused']: state['i'] += 1 if state['i'] >= panel.n + HOLD: state['i'] = 0 panel.draw(min(state['i'], panel.n - 1)) return () ani = animation.FuncAnimation(fig, tick, interval=25, blit=False, cache_frame_data=False, save_count=1) fig._ani = ani # keep a reference alive plt.show() # ── Rendering to a file ─────────────────────────────────────────────────────── def run_save(path, f0, ratios=(1.0, 10.0), fps=25, n_orbits=1.0, n_samples=200, dpi=80): """Render a side-by-side comparison of two mass ratios to a GIF/MP4.""" sims = [simulate(q, f0, n_orbits=n_orbits, n_samples=n_samples) for q in ratios] limits = view_limits(sims) fig, axes = plt.subplots(1, 2, figsize=(11.0, 5.4)) fig.patch.set_facecolor(BG) fig.suptitle('Two-Body Problem — both bodies orbit the common centre of mass', color=FG, fontsize=12, fontweight='bold', y=0.98) panels = [] for ax, sim in zip(axes, sims): p = Panel(ax) p.set_sim(sim, limits=limits) panels.append(p) fig.legend(handles=legend_handles(), loc='upper center', bbox_to_anchor=(0.5, 0.945), ncol=6, fontsize=8.5, frameon=False, labelcolor=FG, columnspacing=1.6, handletextpad=0.6) fig.text(0.5, 0.018, f'G = 1, $m_1 + m_2 = 1$, eccentricity e = {sims[0]["e"]:.2f} — ' 'as $m_1/m_2$ grows the heavier body barely moves and the centre of ' 'mass sinks towards it', color=DIM, fontsize=8.5, ha='center') fig.tight_layout(rect=(0, 0.05, 1, 0.90)) n = max(p.n for p in panels) def frame(k): for p in panels: p.draw(min(k, p.n - 1)) return () ani = animation.FuncAnimation(fig, frame, frames=n + HOLD, blit=False) os.makedirs(os.path.dirname(os.path.abspath(path)), exist_ok=True) writer = 'ffmpeg' if path.endswith('.mp4') else 'pillow' print(f'Writing {path} ({n + HOLD} frames) …') ani.save(path, writer=writer, fps=fps, dpi=dpi, savefig_kwargs=dict(facecolor=BG)) print('done.') # ── Entry point ─────────────────────────────────────────────────────────────── def main(): ap = argparse.ArgumentParser(description=__doc__.split('\n')[1]) ap.add_argument('--ratio', type=float, default=Q_DEFAULT, help='initial mass ratio m1/m2 (default: %(default)s)') ap.add_argument('--speed', type=float, default=F_DEFAULT, help='initial speed as a fraction of the circular speed; ' '0 = both bodies at rest (default: %(default)s)') ap.add_argument('--save', metavar='PATH', default=None, help='render a two-panel comparison to PATH (.gif or .mp4) ' 'instead of opening the interactive window') ap.add_argument('--save-ratios', type=float, nargs=2, default=(1.0, 10.0), metavar=('Q1', 'Q2'), help='mass ratios for the two saved panels') args = ap.parse_args() if args.save: run_save(args.save, args.speed, ratios=tuple(args.save_ratios)) else: import matplotlib if matplotlib.get_backend().lower() == 'agg': raise SystemExit( 'The interactive sliders need a GUI backend, but matplotlib is ' 'using "agg".\nInstall e.g. python3-tk, or use --save to render ' 'a file instead.' ) run_interactive(args.ratio, args.speed) if __name__ == '__main__': main()