Skip to content

Two bodies meeting, live

Two bodies' velocities along the line between their centres after they meet: momentum kept, and the speed they part at e times the speed they met at (Newton's law of restitution), v_a′ = (p + m_b e (v_b − v_a)) / (m_a + m_b) and v_b′ = (p + m_a e (v_a − v_b)) / (m_a + m_b). An elastic hit, e = 1, keeps the kinetic energy; a plastic one, e = 0, leaves both moving together. Below, a 1 kg body at 2 m/s meets a 3 kg one at −1 m/s at a restitution swept from none to one, an eightieth a step, twenty steps a second, and the library's impact (lusine.integrator.impact) gives their velocities after on your GPU as you read. The share of their closing energy ½ μ u² (μ the reduced mass, u the closing speed) that the hit takes is 1 − e², which the program works out from the energies before and after. Beside it a pair whose masses and velocities change each step is held to the laws the impact must keep.

This page runs the program on your GPU through WebGPU, so it needs a browser with WebGPU (a recent Chrome, Edge, Safari or Firefox). Without it the program does not run and nothing is drawn below.

What it shows

  • The strip above the controls is green while every check holds, and turns red for good once one fails.
  • The swept pair: the chart draws the share of the closing energy lost against the restitution, a point a step over one whole sweep, and 1 − e² beside it: the two lie one over the other, from everything lost in a plastic hit to nothing in an elastic one. Below it, the restitution and the two velocities after.
  • The checks: equal masses of 1.25 kg at 3 and −0.5 m/s meeting elastically, the first leaving at −0.5 m/s; a 1 kg body at 2 m/s against one a million times heavier at rest, a wall, coming back at −2e m/s; energy lost, 1 − e², within 4 × 10⁻⁶ at every restitution; momentum kept and Newton's restitution, for the changing pair at every step, each within 16 roundings of its terms; elastic and plastic, equal masses swapping their velocities exactly, a plastic hit leaving both at exactly one speed, and an elastic one keeping the kinetic energy within 64 roundings; a wall, the light body coming back at −e times its speed within 20 µm/s and the wall unmoved.

Every number here is computed by the program on the GPU; the page only reads the values back and draws them. The pair and the sweep are the program's settings, which the host may write; this page offers no control to change them yet.

The program

The sweep, the impacts and the checks, the panels the page draws, and the one procedure, main, run each step: the sweep moves on, each pair meets, and the checks are taken.

// The impact's live page: two bodies meeting at a restitution swept from none to one, the library's impact giving
// their velocities after, the share of the closing energy lost beside 1 - e^2, and pairs that change each step held to
// momentum, Newton's restitution, an elastic swap, a plastic hit and a wall, with the page's checks computed here on
// the GPU and shown by the panels the program declares.

namespace impact;

lusine.version version = { { number = 1 } };

// The reader's settings: the swept pair, a 1 kg body at 2 m/s meeting a 3 kg one at -1 m/s; and how far each step
// moves the restitution, from none.
host struct {
  float<kg> light = 1.0 kg; float<m/s> light_speed = 2.0 m/s; float<kg> heavy = 3.0 kg;
  float<m/s> heavy_speed = -1.0 m/s;
  float stride = 0.0125;
} settings;

// The sweep: the step it is at, from 0 to 80 and round again, the restitution there, and an angle that changes the
// other pairs each step.
struct { integer step = -1; float restitution; float turn; } sweep;

// The swept pair after; the share of its closing energy (1/2 mu u^2, mu the reduced mass, u the closing speed) lost,
// and 1 - e^2; a pair that changes each step, its masses and velocities, after at the swept restitution, after at
// e = 0 and after at e = 1; and a light body against one a million times heavier, at rest.
struct {
  float<m/s> light; float<m/s> heavy; float lost; float closed_form;
  float<kg> mass_a; float<kg> mass_b; float<m/s> va; float<m/s> vb; float<m/s> a; float<m/s> b;
  float<m/s> plastic_a; float<m/s> plastic_b; float<m/s> elastic_a; float<m/s> elastic_b;
  float<m/s> wall_a; float<m/s> wall_b; float<m/s> swap_a; float<m/s> swap_b;
} hit;

// The checks: the share lost within 0.0004% of 1 - e^2 at every restitution; momentum kept and the bodies parting at
// e times the speed they met at, each within 16 roundings of its terms, for every pair; equal masses swapping their
// velocities exactly at e = 1, a plastic hit leaving both at one speed exactly, an elastic one keeping the kinetic
// energy within 64 roundings; and a body against a wall coming back at -e times its speed within 0.001%.
struct {
  bool energy = true; bool momentum = true; bool restitution = true; bool limits = true; bool wall = true;
} checks;

// A strip of lamps the page draws, lit green while every check holds and red once one fails.
struct lamp { vector<float, 4> colour; };
lattice<lamp, 2> lamps(spacing = 1 m, boundary = wall);

vector<float, 4> lamp_colour(bool holds)
{
  return holds ? vector<float, 4>(0.15, 0.6, 0.3, 1.0) : vector<float, 4>(0.75, 0.15, 0.15, 1.0);
}

float<J> kinetic(float<kg> mass, float<m/s> velocity) { return 0.5 * mass * velocity * velocity; }

// The share of the closing energy lost: the kinetic energy before less after, over 1/2 mu u^2.
float lost(float<kg> ma, float<m/s> va, float<kg> mb, float<m/s> vb, float<m/s> a, float<m/s> b)
{
  let mu = ma * mb / (ma + mb);
  let u = va - vb;
  return (kinetic(mass = ma, velocity = va) + kinetic(mass = mb, velocity = vb) - kinetic(mass = ma, velocity = a)
          - kinetic(mass = mb, velocity = b)) / (0.5 * mu * u * u);
}

// The scale of a pair's terms, against which its roundings are counted.
float<kg*m/s> scale(float<kg> ma, float<m/s> va, float<kg> mb, float<m/s> vb)
{
  return abs(ma * va) + abs(mb * vb) + (ma + mb) * abs(va - vb);
}

output table<lusine.widget> panels = {
  lusine.series { { panel = "the swept pair", chart = "share of the closing energy lost", quantity = hit.lost,
                    steps = 81, against = sweep.restitution, label = "the kernel" } },
  lusine.series { { panel = "the swept pair", chart = "share of the closing energy lost", quantity = hit.closed_form,
                    steps = 81, against = sweep.restitution, label = "1 - e^2" } },
  lusine.value { { panel = "the swept pair", quantity = sweep.restitution, label = "restitution" } },
  lusine.value { { panel = "the swept pair", quantity = hit.light, label = "1 kg after" } },
  lusine.value { { panel = "the swept pair", quantity = hit.heavy, label = "3 kg after" } },
  lusine.value { { panel = "the checks", quantity = hit.swap_a, label = "equal masses, 3 m/s becomes" } },
  lusine.value { { panel = "the checks", quantity = hit.wall_a, label = "at a wall, 2 m/s becomes" } },
  lusine.value { { panel = "the checks", quantity = checks.energy, label = "energy lost, 1 - e^2" } },
  lusine.value { { panel = "the checks", quantity = checks.momentum, label = "momentum kept" } },
  lusine.value { { panel = "the checks", quantity = checks.restitution, label = "Newton's restitution" } },
  lusine.value { { panel = "the checks", quantity = checks.limits, label = "elastic and plastic" } },
  lusine.value { { panel = "the checks", quantity = checks.wall, label = "a wall" } },
};

void main()
{
  sweep.step = sweep.step >= 80 ? 0 : sweep.step + 1;
  sweep.restitution = float(sweep.step) * settings.stride;
  sweep.turn = sweep.turn + 0.61;
  (hit.light, hit.heavy) = lusine.integrator.impact(mass_a = settings.light, velocity_a = settings.light_speed,
                                                    mass_b = settings.heavy, velocity_b = settings.heavy_speed,
                                                    sweep.restitution);
  hit.lost = lost(ma = settings.light, va = settings.light_speed, mb = settings.heavy, vb = settings.heavy_speed,
                  a = hit.light, b = hit.heavy);
  hit.closed_form = 1.0 - sweep.restitution * sweep.restitution;
  hit.mass_a = 0.1 kg + 2.45 kg * (1.0 + cos(sweep.turn));
  hit.mass_b = 0.1 kg + 2.45 kg * (1.0 + sin(1.7 * sweep.turn));
  hit.va = 3.0 m/s * cos(2.3 * sweep.turn);
  hit.vb = 3.0 m/s * sin(0.9 * sweep.turn);
  (hit.a, hit.b) = lusine.integrator.impact(hit.mass_a, velocity_a = hit.va, hit.mass_b, velocity_b = hit.vb,
                                            sweep.restitution);
  (hit.plastic_a, hit.plastic_b) = lusine.integrator.impact(hit.mass_a, velocity_a = hit.va, hit.mass_b,
                                                            velocity_b = hit.vb, restitution = 0.0);
  (hit.elastic_a, hit.elastic_b) = lusine.integrator.impact(hit.mass_a, velocity_a = hit.va, hit.mass_b,
                                                            velocity_b = hit.vb, restitution = 1.0);
  (hit.swap_a, hit.swap_b) = lusine.integrator.impact(mass_a = 1.25 kg, velocity_a = 3.0 m/s, mass_b = 1.25 kg,
                                                      velocity_b = -0.5 m/s, restitution = 1.0);
  (hit.wall_a, hit.wall_b) = lusine.integrator.impact(mass_a = 1.0 kg, velocity_a = 2.0 m/s, mass_b = 1000000.0 kg,
                                                      velocity_b = 0.0 m/s, sweep.restitution);
  checks.energy = checks.energy && abs(hit.lost - hit.closed_form) <= 0.000004;
  checks.momentum = checks.momentum
                    && abs(hit.mass_a * hit.a + hit.mass_b * hit.b - (hit.mass_a * hit.va + hit.mass_b * hit.vb))
                       <= 0.000002 * scale(ma = hit.mass_a, va = hit.va, mb = hit.mass_b, vb = hit.vb);
  checks.restitution = checks.restitution
                       && abs((hit.b - hit.a) - sweep.restitution * (hit.va - hit.vb))
                          <= 0.000002 * scale(ma = hit.mass_a, va = hit.va, mb = hit.mass_b, vb = hit.vb)
                             / (hit.mass_a < hit.mass_b ? hit.mass_a : hit.mass_b);
  checks.limits = checks.limits && hit.swap_a == -0.5 m/s && hit.swap_b == 3.0 m/s && hit.plastic_a == hit.plastic_b
                  && abs(kinetic(mass = hit.mass_a, velocity = hit.elastic_a)
                         + kinetic(mass = hit.mass_b, velocity = hit.elastic_b)
                         - kinetic(mass = hit.mass_a, velocity = hit.va)
                         - kinetic(mass = hit.mass_b, velocity = hit.vb))
                     <= 0.0000077 * scale(ma = hit.mass_a, va = hit.va, mb = hit.mass_b, vb = hit.vb) * 3.0 m/s;
  checks.wall = checks.wall && abs(hit.wall_a + sweep.restitution * 2.0 m/s) <= 0.00002 m/s
                && abs(hit.wall_b) <= 0.00002 m/s;
  lamps.colour = lamp_colour(holds = checks.energy && checks.momentum && checks.restitution && checks.limits
                                     && checks.wall);
}

The law's other proof (momentum and restitution over 500 impacts at random) is on its proof page.