A rolling sphere's spin¶
The angular velocity of a sphere rolling without slipping on a surface, w = (n x v) / r, n the surface's unit normal, v the sphere's velocity along it and r its radius: the point touching the surface, at -r n from the centre, then moves at v + w x (-r n) = 0.
Its component about n (a spin in place) is the caller's. Source: Taylor, Classical Mechanics (2005), 10.1 (rolling without slipping); Goldstein, Classical Mechanics, 3rd ed., 1.3. Checked against a reference in numpy, which has its own tests.
- Source
- Taylor, Classical Mechanics (2005), 10.1.
- Proof
- reference, invariant.
rolling_spin.trait
namespace lusine.kinematics;
vector<float<rad/s>, 3> rolling_spin(vector<float, 3> normal, vector<float<m/s>, 3> velocity, float<m> radius)
{
return cross(normal, velocity) / radius * 1.0 rad;
}
The angular velocity of a sphere rolling without slipping, (n x v) / r: the point touching the surface stands still. A boule rolling on the terrain turns so, and its grooves with it.
import os, sys
root = os.getcwd()
while not os.path.isfile(os.path.join(root, 'docsite.json')):
root = os.path.dirname(root)
sys.path.insert(0, os.path.join(root, 'tools'))
import numpy as np
import matplotlib.pyplot as plt
from lusine import trait_run
F32 = np.float32
U = float(np.finfo(F32).eps) / 2 # binary32's unit roundoff, 2^-24
rng = np.random.default_rng(2026)
The kernel¶
The law as the library states it.
print(trait_run.source('kinematics/rolling_spin.trait'))
// lusine.kinematics.rolling_spin: the angular velocity of a sphere rolling without slipping on a surface, w = (n x v)
// / r, n the surface's unit normal, v the sphere's velocity along it and r its radius: the point touching the surface,
// at -r n from the centre, then moves at v + w x (-r n) = 0. Its component about n (a spin in place) is the caller's.
// Source: Taylor, Classical Mechanics (2005), 10.1 (rolling without slipping); Goldstein, Classical Mechanics, 3rd
// ed., 1.3.
// Checked against a reference in numpy, which has its own tests.
namespace lusine.kinematics;
vector<float<rad/s>, 3> rolling_spin(vector<float, 3> normal, vector<float<m/s>, 3> velocity, float<m> radius)
{
return cross(normal, velocity) / radius * 1.0 rad;
}
The program¶
Each section below runs the kernel over a lattice of items, one a case: the host sets every parameter, and main() binds the kernel over the items.
print(trait_run.law_program('kinematics/rolling_spin.trait', 'rolling_spin'))
namespace probe;
lusine.version version = { { number = 1 } };
struct item {
host vector<float, 3> normal;
host vector<float<m/s>, 3> velocity;
host float<m> radius;
vector<float<rad/s>, 3> result_0;
};
lattice<item, 1> items(spacing = 1 m, boundary = wall);
void main() { items.result_0 = lusine.kinematics.rolling_spin(items.normal, items.velocity, items.radius); }
from lusine.kernels.kinematics.rolling_spin import rolling_spin as reference
Reference: rolls at random on a floor¶
Velocities along a floor with y up, a boule's radius: the kernel's spin within binary32's roundings of the reference's.
v = F32(rng.normal(0.0, 3.0, (400, 3)) * np.array([1.0, 0.0, 1.0])).astype(float)
up = np.array([0.0, 1.0, 0.0])
got = trait_run.law('kinematics/rolling_spin.trait', 'rolling_spin', normal=up, velocity=v, radius=0.037)
want = reference(up, v, float(F32(0.037)))
err = np.max(np.linalg.norm(got - want, axis=1) / np.linalg.norm(want, axis=1))
assert err <= 8 * U, err
print(f'{len(v)} rolls within {err / U:.1f} unit roundoffs')
400 rolls within 0.9 unit roundoffs
Invariant: no slip¶
The point touching the floor, at -r n from the centre, moves at v + w x (-r n): zero, to the roundings.
contact = v + np.cross(got, -0.037 * up)
assert np.max(np.linalg.norm(contact, axis=1) / np.linalg.norm(v, axis=1)) <= 1e-6
print('the touching point stands still')
the touching point stands still