Skip to content

Restitution and friction at a contact, live

What a contact does to a body's velocity, given the unit normal n out of the surface: the normal part, v_n = v · n, reversed and scaled by the restitution e (Newton's law of restitution); the tangential part, the slide, reduced by the friction coefficient μ times the normal impulse (1 + e) |v_n|, never past stopping (Coulomb's law as an impulse); and a body already leaving the surface untouched. Below, three balls are dropped from 2 m onto a floor, with restitutions of 0.5, 0.7 and 0.9, twenty steps a bounce and thirty steps a second: the library's contact (lusine.integrator.contact) turns each ball at the floor on your GPU as you read, and the flight between bounces is worked out exactly, the ball meeting the floor at √(2 g h) and climbing to v² / (2 g), so each bounce reaches e² of the last. Beside them, each step, a contact on a surface of another tilt, at another restitution and friction, is held to the two laws.

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 bounces: the chart draws each ball's height against the bounces it has made, twenty points a flight: an arch a bounce, each e² as high as the last, the stiffest ball's arches shrinking slowest. After eight bounces the balls are dropped again. Below it, the bounces, and the middle ball's height at the top of this flight beside e^(2n) of the drop.
  • A tilted surface: this step's contact, its normal speed in and out and its slide in and out.
  • The checks: each bounce e² of the last, every ball's height within 0.0008% of e^(2n) of the drop at every bounce; Newton's restitution, on every tilted surface the normal part after −e times before within 16 roundings; friction never reverses, the slide after neither faster than before nor turned against it; a body leaving untouched, a body moving away from the surface coming back exactly as it went in; the limits, on a floor, every value exact: an elastic bounce keeping its normal speed whole, a plastic one stopping it, no friction leaving the slide alone, a slide of 3 m/s losing μ (1 + e) 2 m/s = 0.6 m/s, and one of 0.5 m/s, less than its share, stopped, never reversed. Each is held at every step.

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

The program

The balls, the tilted contact and the checks, the panels the page draws, and the one procedure, main, run each step: one bounce of each ball, one contact on a tilted surface, and the checks.

// The contact's live page: three balls dropped on a floor, one bounce a step, the library's contact at each bounce and
// the free flight between worked out exactly, each bounce reaching e^2 of the last; and contacts on surfaces of every
// tilt held to Newton's restitution and Coulomb's friction, with the page's checks computed here on the GPU and shown
// by the panels the program declares.

namespace contact;

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

// The reader's settings: gravity; the height the balls are dropped from; their three restitutions; and how many
// bounces before they are dropped again.
host struct {
  float<m/s2> gravity = 9.80665 m/s2; float<m> drop = 2.0 m; float low = 0.5; float middle = 0.7; float high = 0.9;
  integer bounces = 8;
} settings;

// The balls: the bounce they are at, from 0 to 8 and round again, and the point of the flight after it, 0 to 19; the
// bounces and flights so far, a number for the chart; each one's height at the top of its flight; the heights e^(2n)
// times the drop gives; and each one's height at this point of its flight, h 4 f (1 - f) for f the share of the flight
// flown.
struct {
  integer bounce = 8; integer point = 19; float bounces; float<m> low; float<m> middle; float<m> high;
  float<m> low_closed; float<m> middle_closed; float<m> high_closed; float<m> low_now; float<m> middle_now;
  float<m> high_now;
} balls;

// A contact on a tilted surface, changing each step: its angle, its surface's normal, its velocity into the surface,
// its restitution and friction, and its velocity after; and the normal and tangential parts before and after.
struct {
  float turn; vector<float, 3> normal; vector<float<m/s>, 3> velocity; float restitution; float friction;
  vector<float<m/s>, 3> after; float<m/s> normal_before; float<m/s> normal_after; float<m/s> slide_before;
  float<m/s> slide_after;
} tilted;

// The checks: every bounce at e^(2n) of the drop within 0.0008%; on every tilted surface the normal part after -e times
// before within 16 roundings, and the slide after neither faster nor reversed; a body leaving the surface untouched,
// exactly; and the limits on a floor exact: an elastic bounce keeping the normal speed, a plastic one stopping it, no
// friction leaving the slide, a slide of 3 m/s losing mu (1 + e) 2 m/s, and one of 0.5 m/s stopped, never reversed.
struct {
  bool bounces = true; bool restitution = true; bool friction = true; bool leaving = true; bool limits = 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);
}

vector<float, 3> up() { return vector<float, 3>(0.0, 0.0, 1.0); }

// One bounce: the ball falls from its height to the floor at sqrt(2 g h), the contact turns it, and it climbs to
// v^2 / (2 g).
float<m> bounced(float<m> height, float restitution, float<m/s2> gravity)
{
  let falling = vector<float<m/s>, 3>(0.0 m/s, 0.0 m/s, -sqrt(2.0 * gravity * height));
  let rising = lusine.integrator.contact(velocity = falling, normal = up(), restitution, friction = 0.0);
  return rising.z * rising.z / (2.0 * gravity);
}

float<m> closed(float<m> drop, float restitution, float count) { return drop * pow(restitution, 2.0 * count); }

float<m> flying(float<m> top, integer point)
{
  let f = float(point) / 20.0;
  return top * 4.0 * f * (1.0 - f);
}

bool near(float<m> got, float<m> want) { return abs(got - want) <= 0.000008 * want; }

vector<float<m/s>, 3> floor_hit(float<m/s> slide, float restitution, float friction)
{
  return lusine.integrator.contact(velocity = vector<float<m/s>, 3>(slide, 0.0 m/s, -2.0 m/s), normal = up(),
                                   restitution, friction);
}

// The limits on a floor, every value exact in binary32 but the one slowed slide's, within a rounding.
bool limits()
{
  let elastic = floor_hit(slide = 1.5 m/s, restitution = 1.0, friction = 0.0);
  let plastic = floor_hit(slide = 1.5 m/s, restitution = 0.0, friction = 0.0);
  let smooth = floor_hit(slide = 1.5 m/s, restitution = 0.4, friction = 0.0);
  let slowed = floor_hit(slide = 3.0 m/s, restitution = 0.5, friction = 0.2);
  let stopped = floor_hit(slide = 0.5 m/s, restitution = 0.5, friction = 0.5);
  return elastic == vector<float<m/s>, 3>(1.5 m/s, 0.0 m/s, 2.0 m/s)
      && plastic == vector<float<m/s>, 3>(1.5 m/s, 0.0 m/s, 0.0 m/s)
      && smooth.x == 1.5 m/s && smooth.y == 0.0 m/s
      && abs(slowed.x - 2.4 m/s) <= 0.0000005 * 3.0 m/s && slowed.z == 1.0 m/s
      && stopped.x == 0.0 m/s && stopped.z == 1.0 m/s;
}

// A body leaving the surface comes back exactly as it went in.
bool leaves(vector<float, 3> normal, vector<float<m/s>, 3> velocity)
{
  let away = velocity - (2.0 * dot(velocity, normal)) * normal;
  return lusine.integrator.contact(velocity = away, normal, restitution = 0.5, friction = 0.5) == away;
}

output table<lusine.widget> panels = {
  lusine.series { { panel = "the bounces", chart = "height", quantity = balls.low_now, steps = 180,
                    against = balls.bounces, label = "e = 0.5" } },
  lusine.series { { panel = "the bounces", chart = "height", quantity = balls.middle_now, steps = 180,
                    against = balls.bounces, label = "e = 0.7" } },
  lusine.series { { panel = "the bounces", chart = "height", quantity = balls.high_now, steps = 180,
                    against = balls.bounces, label = "e = 0.9" } },
  lusine.value { { panel = "the bounces", quantity = balls.bounces, label = "bounces" } },
  lusine.value { { panel = "the bounces", quantity = balls.middle, label = "e = 0.7's height" } },
  lusine.value { { panel = "the bounces", quantity = balls.middle_closed, label = "e^(2n) of the drop" } },
  lusine.value { { panel = "a tilted surface", quantity = tilted.normal_before, label = "normal speed in" } },
  lusine.value { { panel = "a tilted surface", quantity = tilted.normal_after, label = "normal speed out" } },
  lusine.value { { panel = "a tilted surface", quantity = tilted.slide_before, label = "slide in" } },
  lusine.value { { panel = "a tilted surface", quantity = tilted.slide_after, label = "slide out" } },
  lusine.value { { panel = "the checks", quantity = checks.bounces, label = "each bounce e^2 of the last" } },
  lusine.value { { panel = "the checks", quantity = checks.restitution, label = "Newton's restitution" } },
  lusine.value { { panel = "the checks", quantity = checks.friction, label = "friction never reverses" } },
  lusine.value { { panel = "the checks", quantity = checks.leaving, label = "a body leaving untouched" } },
  lusine.value { { panel = "the checks", quantity = checks.limits, label = "the limits" } },
};

void main()
{
  // A step is a twentieth of a flight; at each flight's start the balls bounce, and every nine bounces they are
  // dropped again.
  balls.bounce = balls.point >= 19 ? (balls.bounce >= settings.bounces ? 0 : balls.bounce + 1) : balls.bounce;
  balls.low = balls.point < 19 ? balls.low : balls.bounce == 0 ? settings.drop
            : bounced(height = balls.low, restitution = settings.low, settings.gravity);
  balls.middle = balls.point < 19 ? balls.middle : balls.bounce == 0 ? settings.drop
               : bounced(height = balls.middle, restitution = settings.middle, settings.gravity);
  balls.high = balls.point < 19 ? balls.high : balls.bounce == 0 ? settings.drop
             : bounced(height = balls.high, restitution = settings.high, settings.gravity);
  balls.point = balls.point >= 19 ? 0 : balls.point + 1;
  balls.bounces = float(balls.bounce) + float(balls.point) / 20.0;
  balls.low_closed = closed(settings.drop, restitution = settings.low, count = float(balls.bounce));
  balls.middle_closed = closed(settings.drop, restitution = settings.middle, count = float(balls.bounce));
  balls.high_closed = closed(settings.drop, restitution = settings.high, count = float(balls.bounce));
  balls.low_now = flying(top = balls.low, balls.point);
  balls.middle_now = flying(top = balls.middle, balls.point);
  balls.high_now = flying(top = balls.high, balls.point);
  checks.bounces = checks.bounces && near(got = balls.low, want = balls.low_closed)
                   && near(got = balls.middle, want = balls.middle_closed)
                   && near(got = balls.high, want = balls.high_closed);
  // A surface of another tilt each step, and a velocity into it at 2 m/s along its normal.
  tilted.turn = tilted.turn + 0.73;
  tilted.normal = normalize(vector<float, 3>(cos(tilted.turn), sin(1.3 * tilted.turn), 0.2 + cos(0.7 * tilted.turn)));
  tilted.velocity = vector<float<m/s>, 3>(3.0 m/s * sin(2.1 * tilted.turn), 2.0 m/s * cos(0.4 * tilted.turn), 1.0 m/s);
  tilted.velocity = tilted.velocity - (dot(tilted.velocity, tilted.normal) + 2.0 m/s) * tilted.normal;
  tilted.restitution = 0.5 + 0.5 * sin(0.9 * tilted.turn);
  tilted.friction = 0.5 + 0.5 * cos(1.7 * tilted.turn);
  tilted.after = lusine.integrator.contact(tilted.velocity, tilted.normal, tilted.restitution, tilted.friction);
  tilted.normal_before = dot(tilted.velocity, tilted.normal);
  tilted.normal_after = dot(tilted.after, tilted.normal);
  tilted.slide_before = length(tilted.velocity - tilted.normal_before * tilted.normal);
  tilted.slide_after = length(tilted.after - tilted.normal_after * tilted.normal);
  checks.restitution = checks.restitution
                       && abs(tilted.normal_after + tilted.restitution * tilted.normal_before)
                          <= 0.000002 * length(tilted.velocity);
  checks.friction = checks.friction
                    && tilted.slide_after <= tilted.slide_before + 0.000002 * length(tilted.velocity)
                    && dot(tilted.after - tilted.normal_after * tilted.normal,
                           tilted.velocity - tilted.normal_before * tilted.normal)
                       >= -0.000002 * length(tilted.velocity) * tilted.slide_before;
  checks.leaving = checks.leaving && leaves(tilted.normal, tilted.velocity);
  checks.limits = checks.limits && limits();
  lamps.colour = lamp_colour(holds = checks.bounces && checks.restitution && checks.friction && checks.leaving
                                     && checks.limits);
}

The law's other proof (two hundred contacts at random) is on its proof page.