The Magnus force, live¶
The lift on a spinning body, the acceleration it takes per unit mass from the air it moves through: a = (ρ C_L A / 2m) |u| (ω̂ × u), u its velocity through the air, ω̂ its spin's axis, ρ the air's density, C_L its lift coefficient (a sphere's from its spin ratio) and A its frontal area. It is at right angles to the motion, so it does no work: backspin holds a ball up, topspin presses it down. Below, four tennis balls are lobbed alike at 20 m/s and 50° and stepped by the library's velocity Verlet (lusine.integrator.kick_drift, lusine.integrator.kick) under gravity, the drag of the wake (lusine.laws.quadratic_drag) and the lift of the spin (lusine.laws.magnus, its coefficient lusine.laws.lift_coefficient of lusine.laws.spin_ratio), two milliseconds a step and four hundred steps a second, on your GPU as you read: one with 100 rad/s of backspin, one with none, one with topspin, and a fourth with backspin and no drag, whose energy the same program adds up at every step.
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 lobs: the chart draws the three lobbed balls' heights against their time, a point a step: they rise alike, and the backspun ball comes down last, the topspun one first. Every 3.2 seconds the balls are lobbed again. Below it, where each landed, the first step its height was below the ground.
- The checks: the chart draws the ball with no drag's energy drift, kinetic and potential per unit mass, at each step as a share of what the lift does no work allows, 0.02% of its start, on a log scale, so the check holds while it stays below 1; the drift is the integrator's. Below it, how far that ball's path has bent from the parabola it would fly with no spin: by most of a metre, though its energy has not changed. in order holds while every lob has landed backspin farthest and topspin shortest; no spin, no lift while a ball with no spin has felt no lift at all.
Every number here is computed by the program on the GPU; the page only reads the values back and draws them. Gravity, the air's density, the ball, the lob, the spins and the step are the program's settings, which the host may write; this page offers no control to change them yet.
The program¶
The balls, the throw and the checks, the panels the page draws, and the one procedure, main, run each step: one velocity Verlet step for every ball, its acceleration taken at the step's start and at its middle, and the checks after it.
// The Magnus force's live page: three tennis balls lobbed alike, one with backspin, one with none and one with
// topspin, and a fourth with backspin and no drag, each stepped by the library's velocity Verlet under gravity, the
// drag of its wake and the lift of its spin, with the page's checks computed here on the GPU and shown by the panels
// the program declares.
namespace magnus;
lusine.version version = { { number = 1 } };
// The reader's settings: gravity, the air's density, a tennis ball's frontal area, radius and mass, the lob's speed and
// its angle's cosine and sine (50 degrees), the height it leaves the hand at, the step, and how long each lob is
// followed.
host struct {
float<m/s2> gravity = 9.80665 m/s2; float<kg/m3> density = 1.204 kg/m3; float<m2> area = 0.0034211944 m2;
float<m> radius = 0.033 m; float<kg> mass = 0.057 kg; float<m/s> speed = 20.0 m/s; float cosine = 0.64278761;
float sine = 0.76604444; float<m> hand = 1.0 m; float<s> dt = 0.002 s; float<s> flight = 3.2 s;
} settings;
// The balls: each one's number, its spin about the axis across the throw (positive the ball's top moving back, a
// backspin), its drag coefficient, where it is and how it moves, and where it landed, once it has: backspin, no spin,
// topspin, and backspin with no drag.
struct ball {
integer number; float<rad/s> spin; float drag_coefficient; vector<float<m>, 3> position;
vector<float<m/s>, 3> velocity; vector<float<m/s2>, 3> acceleration; bool landed; float<m> landing;
};
table<ball> balls(capacity = 4) = {
{ number = 0, spin = 100.0 rad/s, drag_coefficient = 0.55 },
{ number = 1, spin = 0.0 rad/s, drag_coefficient = 0.55 },
{ number = 2, spin = -100.0 rad/s, drag_coefficient = 0.55 },
{ number = 3, spin = 100.0 rad/s, drag_coefficient = 0.0 },
};
// The throw: its time since the balls left the hand; each lobbed ball's height and where it landed; the ball with
// no drag's energy per unit mass, kinetic and potential, against what it started with, and the farthest its path has
// bent from the path it takes with no spin, a closed form, a parabola; and the lift on a ball with no spin.
struct {
float<s> time; float<m> back; float<m> plain; float<m> top; float<m> back_lands; float<m> plain_lands;
float<m> top_lands; float<J/kg> energy; float<J/kg> start; float drift_share; float<m> bent;
vector<float<m/s2>, 3> unspun;
} throw;
// The checks: the landings in order, backspin carrying a ball farthest and topspin bringing it down shortest; the
// ball with no drag's energy within 0.02% of its start at every step, though its path bends by most of a metre; and
// no spin, no lift at all.
struct { bool order = true; bool no_work = true; bool no_spin = 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<m/s>, 3> still() { return vector<float<m/s>, 3>(0.0 m/s, 0.0 m/s, 0.0 m/s); }
vector<float<rad/s>, 3> spun(float<rad/s> spin) { return vector<float<rad/s>, 3>(0.0 rad/s, 0.0 rad/s, spin); }
vector<float<m/s>, 3> lobbed(float<m/s> speed, float cosine, float sine)
{
return vector<float<m/s>, 3>(speed * cosine, speed * sine, 0.0 m/s);
}
vector<float<m>, 3> in_hand(float<m> hand) { return vector<float<m>, 3>(0.0 m, hand, 0.0 m); }
bool nothing(vector<float<m/s2>, 3> a) { return a == vector<float<m/s2>, 3>(0.0 m/s2, 0.0 m/s2, 0.0 m/s2); }
// A ball's acceleration: gravity, the wake's drag and the lift of its spin, its lift coefficient from its spin ratio.
vector<float<m/s2>, 3> pulled(vector<float<m/s>, 3> velocity, float<rad/s> spin, float drag_coefficient,
float<m/s2> gravity, float<kg/m3> density, float<m2> area, float<m> radius,
float<kg> mass)
{
let lift = lusine.laws.lift_coefficient(spin_ratio = lusine.laws.spin_ratio(velocity = velocity, wind = still(),
angular_velocity = spun(spin),
radius = radius));
let down = vector<float<m/s2>, 3>(0.0 m/s2, -gravity, 0.0 m/s2);
return (down + lusine.laws.quadratic_drag(velocity = velocity, wind = still(), density = density,
drag_coefficient = drag_coefficient, area = area, mass = mass))
+ lusine.laws.magnus(velocity = velocity, wind = still(), angular_velocity = spun(spin), density = density,
lift_coefficient = lift, area = area, mass = mass);
}
// Where a ball lands: where it is at the first step its height is below the ground.
bool landing_now(vector<float<m>, 3> position, bool landed) { return !landed && position.y < 0.0 m; }
// One ball's height, or where it landed, picked from the table by its number.
float<m> height_of(table<vector<float<m>, 3>> position, table<integer> number, integer which)
{
return sum(position.y) where number == which;
}
float<m> landing_of(table<float<m>> landing, table<integer> number, integer which)
{
return sum(landing) where number == which;
}
float<J/kg> energy_of(table<vector<float<m/s>, 3>> velocity, table<vector<float<m>, 3>> position, table<integer> number,
float<m/s2> gravity)
{
return sum(0.5 * dot(velocity, velocity) + gravity * position.y) where number == 3;
}
// How far the ball with no drag is from where it would be with no spin either: a parabola from the hand.
float<m> bent(table<vector<float<m>, 3>> position, table<integer> number, float<m/s> speed, float cosine, float sine,
float<m> hand, float<m/s2> gravity, float<s> time)
{
let x = speed * cosine * time;
let y = hand + speed * sine * time - 0.5 * gravity * time * time;
return sum(sqrt((position.x - x) * (position.x - x) + (position.y - y) * (position.y - y))) where number == 3;
}
output table<lusine.widget> panels = {
lusine.series { { panel = "the lobs", chart = "height", quantity = throw.back, steps = 1600, against = throw.time,
label = "backspin" } },
lusine.series { { panel = "the lobs", chart = "height", quantity = throw.plain, steps = 1600, against = throw.time,
label = "no spin" } },
lusine.series { { panel = "the lobs", chart = "height", quantity = throw.top, steps = 1600, against = throw.time,
label = "topspin" } },
lusine.value { { panel = "the lobs", quantity = throw.back_lands, label = "backspin lands at" } },
lusine.value { { panel = "the lobs", quantity = throw.plain_lands, label = "no spin lands at" } },
lusine.value { { panel = "the lobs", quantity = throw.top_lands, label = "topspin lands at" } },
lusine.series { { panel = "the checks", chart = "share of the tolerance", quantity = throw.drift_share, steps = 1600,
log = true, label = "energy drift with no drag" } },
lusine.value { { panel = "the checks", quantity = throw.bent, label = "its path bent by" } },
lusine.value { { panel = "the checks", quantity = checks.order, label = "in order" } },
lusine.value { { panel = "the checks", quantity = checks.no_work, label = "the lift does no work" } },
lusine.value { { panel = "the checks", quantity = checks.no_spin, label = "no spin, no lift" } },
};
void main()
{
// Every lob's flight the four balls are lobbed again.
throw.time = throw.time >= settings.flight ? 0.0 s : throw.time;
balls.position = throw.time == 0.0 s ? in_hand(settings.hand) : balls.position;
balls.velocity = throw.time == 0.0 s ? lobbed(settings.speed, settings.cosine, settings.sine) : balls.velocity;
balls.landed = throw.time == 0.0 s ? false : balls.landed;
throw.start = throw.time == 0.0 s ? energy_of(balls.velocity, balls.position, balls.number, settings.gravity)
: throw.start;
balls.acceleration = pulled(balls.velocity, balls.spin, balls.drag_coefficient, settings.gravity, settings.density,
settings.area, settings.radius, settings.mass);
(balls.position, balls.velocity) = lusine.integrator.kick_drift(balls.position, balls.velocity, balls.acceleration,
settings.dt);
balls.acceleration = pulled(balls.velocity, balls.spin, balls.drag_coefficient, settings.gravity, settings.density,
settings.area, settings.radius, settings.mass);
balls.velocity = lusine.integrator.kick(balls.velocity, balls.acceleration, settings.dt);
throw.time = throw.time + settings.dt;
balls.landing = landing_now(balls.position, balls.landed) ? balls.position.x : balls.landing;
balls.landed = balls.landed || balls.position.y < 0.0 m;
throw.back = height_of(balls.position, balls.number, which = 0);
throw.plain = height_of(balls.position, balls.number, which = 1);
throw.top = height_of(balls.position, balls.number, which = 2);
throw.back_lands = landing_of(balls.landing, balls.number, which = 0);
throw.plain_lands = landing_of(balls.landing, balls.number, which = 1);
throw.top_lands = landing_of(balls.landing, balls.number, which = 2);
throw.energy = energy_of(balls.velocity, balls.position, balls.number, settings.gravity);
throw.drift_share = abs(throw.energy - throw.start) / (0.0002 * throw.start);
throw.bent = bent(balls.position, balls.number, settings.speed, settings.cosine, settings.sine, settings.hand,
settings.gravity, throw.time);
checks.order = checks.order && (throw.time < settings.flight || throw.back_lands > throw.plain_lands
&& throw.plain_lands > throw.top_lands);
checks.no_work = checks.no_work && throw.drift_share <= 1.0;
throw.unspun = lusine.laws.magnus(velocity = lobbed(settings.speed, settings.cosine, settings.sine), wind = still(),
angular_velocity = spun(spin = 0.0 rad/s), settings.density,
lift_coefficient = 0.2, settings.area, settings.mass);
checks.no_spin = checks.no_spin && nothing(a = throw.unspun);
lamps.colour = lamp_colour(holds = checks.order && checks.no_work && checks.no_spin);
}
The law's other proof (the lift against a reference over velocities and spins at random) is on its proof page.