Investigating Collisions in Phase Space

Using Newton’s Law of Restitution and the Law of Conservation of Momentum one can calculate the final velocities of two bodies after a collision, given the masses of each body; their initial velocities and the coefficient of restitution between them. I wanted to display this visually, with controllable initial conditions, so I coded something up in python.

The program plots the displacements of the two bodies every second, with one body on the x axis, and the other on the y axis.

Collisions in Phase Space

Each point is the position at an integer number of seconds. The starting point is one that lies at (x2, x1) – the position of the second body is the x-coordinate, the position of the first body is the y-coordinate.
e – coefficient of restitution
u1 – initial velocity of body 1
u2 – initial velocity of body 2
m1 – mass of body 1
m2 – mass of body 2
x1 – initial displacement of body 1
x2 – initial displacement of body 2

The code:
It’s a tad messy… but it works. Requires matplotlib

#!/usr/bin/env python

import numpy as np
import matplotlib.path as mpath
import matplotlib.patches as mpatches
import matplotlib.pyplot as plt

def v1(u1, u2, m1, m2, e):
    return (m2*e*(2*u2 - u1))/(m1+m2)

def v2(u1, u2, m1, m2, e, v1):
    return float(m1*(u1 - v1) + m1*u2)/m2

def t_collide(u1, u2, x1, x2):
    return float(x2 - x1)/(u1 - u2)

def produce_phase_space_graph(ax, u1, u2, m1, m2, e, x1, x2):
    points = []
    collision_time = t_collide(u1, u2, x1, x2)

    max_time = int(2*collision_time)

    for t in range(0, int(collision_time)+1):
        current_x1 = x1 + t*u1
        current_x2 = x2 + t*u2
        points.append((current_x1, current_x2))

    v_1 = v1(u1, u2, m1, m2, e)
    v_2 = v2(u1, u2, m1, m2, e, v_1)

    print v_1
    print v_2

    collision_point = x1 + collision_time*u1

    points.append((collision_point, collision_point))

    for t in range(int(collision_time)+1, max_time+1):
        current_x1 = collision_point + (t-collision_time)*v_1
        current_x2 = collision_point + (t-collision_time)*v_2
        points.append((current_x1, current_x2))

    print points

    x = []
    y = []

    for point in points:
        x.append(point[0])
        y.append(point[1])

    ax.clear()
    ax.plot(y, x, 'ro-')

    ax.grid()
    ax.set_xlim(-12, 12)
    ax.set_ylim(-12, 12)
    plt.show(False)
    plt.draw()

#################
#INTERFACE
#################
import gtk

class Interface():
    params = { #(u1, u2, m1, m2, e, x1, x2)
        'u1': 2,
        'u2': -4,
        'm1': 1,
        'm2': 1,
        'e': 0.8,
        'x1': 0,
        'x2': 12,
        }

    fig = None
    ax = None

    def __init__(self):
        self.init_gui()
        self.update()

    def init_gui(self):
        self.dialog = gtk.Dialog("Collisions")
        for param,default in self.params.items():
            self.dialog.vbox.pack_start(self.create_parameter_control(param))
        self.dialog.show_all()

    def create_parameter_control(self, parameter):
        label = gtk.Label(parameter)
        if parameter != 'e':
            adjustment = gtk.Adjustment(
                value=self.params[parameter],
                lower=-100,
                upper=100,
                step_incr=0.5,
                page_incr=1,
                page_size=0
                )
        else:
            adjustment = gtk.Adjustment(
                value=self.params[parameter],
                lower = 0,
                upper = 1,
                step_incr = 0.05
                )
        spinbutton = gtk.SpinButton(adjustment=adjustment, digits=2)
        spinbutton.connect('value_changed', self.update_parameter, parameter)
        hbox = gtk.HBox()
        hbox.pack_start(label)
        hbox.pack_start(spinbutton)
        hbox.show_all()
        return hbox

    def update_parameter(self, spinbutton, parameter):
        print "in callback"
        value = spinbutton.get_value()
        self.params[parameter] = value
        self.update()
        return True

    def update(self):
        if not self.fig:
            self.fig = plt.figure()
        if not self.ax:
            self.ax = self.fig.add_subplot(111)
        produce_phase_space_graph(self.ax, **self.params)

#produce_phase_space_graph(2, -4, 1, 1, 0.8, 0, 12)

interface = Interface()

gtk.main()