The Magnus force¶

The lift on a spinning body, the acceleration it takes per unit mass from the air it moves through: a = (rho C_L A / 2 m) |u| (w^ x u), u = v - w_air its velocity through the air, w^ its angular velocity's unit axis, rho the air's density, C_L its lift coefficient (a sphere's from its spin ratio, lusine.laws.lift_coefficient) and A its frontal area.

Perpendicular to u, so it does no work: backspin lifts a body, topspin presses it down. Zero with no spin, or no motion through the air. Gravity, drag and any other force are the caller's to add. Source: Nathan, The effect of spin on the flight of a baseball, Am. J. Phys. 76 (2008) 119, eq. 2; Sawicki, Hubbard and Stronge, How to hit home runs, Am. J. Phys. 71 (2003) 1152. Checked against a reference in numpy, which has its own tests.

Source

Nathan, Am. J. Phys. 76 (2008) 119, eq. 2.

Proof

invariant, limit.

magnus.trait

namespace lusine.laws;

vector<float<m/s2>, 3> magnus(vector<float<m/s>, 3> velocity, vector<float<m/s>, 3> wind,
                              vector<float<rad/s>, 3> angular_velocity, float<kg/m3> density, float lift_coefficient,
                              float<m2> area, float<kg> mass)
{
  let u = velocity - wind;
  let spin = sqrt(dot(angular_velocity, angular_velocity));
  if (spin <= 0) return vector<float<m/s2>, 3>(0.0 m/s2, 0.0 m/s2, 0.0 m/s2);
  let axis = angular_velocity / spin;
  let speed = sqrt(dot(u, u));
  return cross(axis, u) * (speed * ((density * lift_coefficient * area) / (2.0 * mass)));
}

The lift on a spinning body, (rho C_L A / 2 m) |u| (w^ x u), perpendicular to its motion through the air: backspin lifts a lob, topspin presses it down, and it does no work.

In [1]:

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)

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.

In [2]:

print(trait_run.source('laws/magnus.trait'))

print(trait_run.source('laws/magnus.trait'))

// lusine.laws.magnus: the lift on a spinning body, the acceleration it takes per unit mass from the air it moves
// through: a = (rho C_L A / 2 m) |u| (w^ x u), u = v - w_air its velocity through the air, w^ its angular velocity's
// unit axis, rho the air's density, C_L its lift coefficient (a sphere's from its spin ratio,
// lusine.laws.lift_coefficient) and A its frontal area. Perpendicular to u, so it does no work: backspin lifts a body,
// topspin presses it down. Zero with no spin, or no motion through the air. Gravity, drag and any other force are the
// caller's to add.
// Source: Nathan, The effect of spin on the flight of a baseball, Am. J. Phys. 76 (2008) 119, eq. 2; Sawicki,
// Hubbard and Stronge, How to hit home runs, Am. J. Phys. 71 (2003) 1152.
// Checked against a reference in numpy, which has its own tests.

namespace lusine.laws;

vector<float<m/s2>, 3> magnus(vector<float<m/s>, 3> velocity, vector<float<m/s>, 3> wind,
                              vector<float<rad/s>, 3> angular_velocity, float<kg/m3> density, float lift_coefficient,
                              float<m2> area, float<kg> mass)
{
  let u = velocity - wind;
  let spin = sqrt(dot(angular_velocity, angular_velocity));
  if (spin <= 0) return vector<float<m/s2>, 3>(0.0 m/s2, 0.0 m/s2, 0.0 m/s2);
  let axis = angular_velocity / spin;
  let speed = sqrt(dot(u, u));
  return cross(axis, u) * (speed * ((density * lift_coefficient * area) / (2.0 * mass)));
}

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.

In [3]:

print(trait_run.law_program('laws/magnus.trait', 'magnus'))

print(trait_run.law_program('laws/magnus.trait', 'magnus'))

namespace probe;

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

struct item {
  host vector<float<m/s>, 3> velocity;
  host vector<float<m/s>, 3> wind;
  host vector<float<rad/s>, 3> angular_velocity;
  host float<kg/m3> density;
  host float lift_coefficient;
  host float<m2> area;
  host float<kg> mass;
  vector<float<m/s2>, 3> result_0;
};
lattice<item, 1> items(spacing = 1 m, boundary = wall);

void main() { items.result_0 = lusine.laws.magnus(items.velocity, items.wind, items.angular_velocity, items.density, items.lift_coefficient, items.area, items.mass); }

In [4]:

from lusine.kernels.laws.magnus import magnus as reference

from lusine.kernels.laws.magnus import magnus as reference

Reference: velocities and spins at random¶

The kernel's acceleration within binary32's roundings of the float64 reference's.

In [5]:

v = F32(rng.normal(0.0, 8.0, (500, 3))).astype(float)
w = F32(rng.normal(0.0, 3.0, (500, 3))).astype(float)
s = F32(rng.normal(0.0, 40.0, (500, 3))).astype(float)
area = float(F32(np.pi * 0.037 ** 2))
got = trait_run.law('laws/magnus.trait', 'magnus', velocity=v, wind=w, angular_velocity=s, density=1.204,
                    lift_coefficient=0.2, area=area, mass=0.7)
want = reference(v, w, s, float(F32(1.204)), float(F32(0.2)), area, float(F32(0.7)))
err = np.max(np.linalg.norm(got - want, axis=1) / np.maximum(np.linalg.norm(want, axis=1), 1e-6))
assert err <= 64 * U, err
print(f'{len(v)} cases within {err / U:.1f} unit roundoffs of the reference')

v = F32(rng.normal(0.0, 8.0, (500, 3))).astype(float) w = F32(rng.normal(0.0, 3.0, (500, 3))).astype(float) s = F32(rng.normal(0.0, 40.0, (500, 3))).astype(float) area = float(F32(np.pi * 0.037 ** 2)) got = trait_run.law('laws/magnus.trait', 'magnus', velocity=v, wind=w, angular_velocity=s, density=1.204, lift_coefficient=0.2, area=area, mass=0.7) want = reference(v, w, s, float(F32(1.204)), float(F32(0.2)), area, float(F32(0.7))) err = np.max(np.linalg.norm(got - want, axis=1) / np.maximum(np.linalg.norm(want, axis=1), 1e-6)) assert err <= 64 * U, err print(f'{len(v)} cases within {err / U:.1f} unit roundoffs of the reference')

500 cases within 8.1 unit roundoffs of the reference

In [6]:

FLIGHT = '''namespace flight;

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

struct body { vector<float<m>, 3> position; vector<float<m/s>, 3> velocity; vector<float<m/s2>, 3> acceleration; };
lattice<body, 1> bodies(spacing = 1 m, boundary = wall);
host struct {
  vector<float<m/s2>, 3> gravity; vector<float<m/s>, 3> wind; vector<float<rad/s>, 3> spin;
  float<kg/m3> density; float drag_coefficient; float<m2> area; float<m> radius; float<kg> mass; float<s> dt;
} settings;

// the body's acceleration: gravity, the wake's drag and the lift of its spin
vector<float<m/s2>, 3> pulled(vector<float<m/s>, 3> velocity, vector<float<m/s2>, 3> gravity, vector<float<m/s>, 3> wind,
                              vector<float<rad/s>, 3> spin, float<kg/m3> density, float drag_coefficient,
                              float<m2> area, float<m> radius, float<kg> mass)
{
  let lift = lusine.laws.lift_coefficient(spin_ratio = lusine.laws.spin_ratio(velocity = velocity, wind = wind,
                                                                               angular_velocity = spin, radius = radius));
  return (gravity + lusine.laws.quadratic_drag(velocity = velocity, wind = wind, density = density,
                                               drag_coefficient = drag_coefficient, area = area, mass = mass))
      + lusine.laws.magnus(velocity = velocity, wind = wind, angular_velocity = spin, density = density,
                           lift_coefficient = lift, area = area, mass = mass);
}

void start()
{
  bodies.acceleration = pulled(bodies.velocity, settings.gravity, settings.wind, settings.spin, settings.density,
                               settings.drag_coefficient, settings.area, settings.radius, settings.mass);
}

void main()
{
  (bodies.position, bodies.velocity) = lusine.integrator.kick_drift(bodies.position, bodies.velocity,
                                                                    bodies.acceleration, settings.dt);
  bodies.acceleration = pulled(bodies.velocity, settings.gravity, settings.wind, settings.spin, settings.density,
                               settings.drag_coefficient, settings.area, settings.radius, settings.mass);
  bodies.velocity = lusine.integrator.kick(bodies.velocity, bodies.acceleration, settings.dt);
}
'''
LIBRARY = ['laws/quadratic_drag.trait', 'laws/magnus.trait', 'laws/spin_ratio.trait', 'laws/lift_coefficient.trait',
           'integrator/kick.trait']
G = [0.0, -9.80665, 0.0]
BOULE = dict(density=1.204, drag_coefficient=0.47, area=float(np.pi * 0.037 ** 2), radius=0.037, mass=0.7)

def fly(velocity, steps, dt=0.001, spin=(0.0, 0.0, 0.0), wind=(0.0, 0.0, 0.0), position=(0.0, 0.0, 0.0), **over):
    run = trait_run.Program(LIBRARY, FLIGHT, {'flight.bodies': [1]})
    settings = dict(BOULE, gravity=G, wind=list(wind), spin=list(spin), dt=dt)
    settings.update(over)
    for name, value in settings.items():
        run.store('flight.settings', name, value)
    run.store('flight.bodies', 'position', [list(position)])
    run.store('flight.bodies', 'velocity', [list(velocity)])
    run.step('flight.start')
    p, v = run.steps('flight.main', steps, [('flight.bodies', 'position'), ('flight.bodies', 'velocity')])
    return p.reshape(steps, 3), v.reshape(steps, 3)

FLIGHT = '''namespace flight; lusine.version version = { { number = 1 } }; struct body { vector<float<m>, 3> position; vector<float<m/s>, 3> velocity; vector<float<m/s2>, 3> acceleration; }; lattice<body, 1> bodies(spacing = 1 m, boundary = wall); host struct { vector<float<m/s2>, 3> gravity; vector<float<m/s>, 3> wind; vector<float<rad/s>, 3> spin; float<kg/m3> density; float drag_coefficient; float<m2> area; float<m> radius; float<kg> mass; float<s> dt; } settings; // the body's acceleration: gravity, the wake's drag and the lift of its spin vector<float<m/s2>, 3> pulled(vector<float<m/s>, 3> velocity, vector<float<m/s2>, 3> gravity, vector<float<m/s>, 3> wind, vector<float<rad/s>, 3> spin, float<kg/m3> density, float drag_coefficient, float<m2> area, float<m> radius, float<kg> mass) { let lift = lusine.laws.lift_coefficient(spin_ratio = lusine.laws.spin_ratio(velocity = velocity, wind = wind, angular_velocity = spin, radius = radius)); return (gravity + lusine.laws.quadratic_drag(velocity = velocity, wind = wind, density = density, drag_coefficient = drag_coefficient, area = area, mass = mass)) + lusine.laws.magnus(velocity = velocity, wind = wind, angular_velocity = spin, density = density, lift_coefficient = lift, area = area, mass = mass); } void start() { bodies.acceleration = pulled(bodies.velocity, settings.gravity, settings.wind, settings.spin, settings.density, settings.drag_coefficient, settings.area, settings.radius, settings.mass); } void main() { (bodies.position, bodies.velocity) = lusine.integrator.kick_drift(bodies.position, bodies.velocity, bodies.acceleration, settings.dt); bodies.acceleration = pulled(bodies.velocity, settings.gravity, settings.wind, settings.spin, settings.density, settings.drag_coefficient, settings.area, settings.radius, settings.mass); bodies.velocity = lusine.integrator.kick(bodies.velocity, bodies.acceleration, settings.dt); } ''' LIBRARY = ['laws/quadratic_drag.trait', 'laws/magnus.trait', 'laws/spin_ratio.trait', 'laws/lift_coefficient.trait', 'integrator/kick.trait'] G = [0.0, -9.80665, 0.0] BOULE = dict(density=1.204, drag_coefficient=0.47, area=float(np.pi * 0.037 ** 2), radius=0.037, mass=0.7) def fly(velocity, steps, dt=0.001, spin=(0.0, 0.0, 0.0), wind=(0.0, 0.0, 0.0), position=(0.0, 0.0, 0.0), **over): run = trait_run.Program(LIBRARY, FLIGHT, {'flight.bodies': [1]}) settings = dict(BOULE, gravity=G, wind=list(wind), spin=list(spin), dt=dt) settings.update(over) for name, value in settings.items(): run.store('flight.settings', name, value) run.store('flight.bodies', 'position', [list(position)]) run.store('flight.bodies', 'velocity', [list(velocity)]) run.step('flight.start') p, v = run.steps('flight.main', steps, [('flight.bodies', 'position'), ('flight.bodies', 'velocity')]) return p.reshape(steps, 3), v.reshape(steps, 3)

Curve: a lob, with and without backspin¶

A boule lobbed at 7 m/s and 50 degrees, gravity, drag and the Magnus force stepped by the library's velocity Verlet: backspin of 40 rad/s about the axis across the throw (the ball's top moving back) holds it up a little longer and carries it farther; topspin brings it down short.

In [7]:

a = np.radians(50.0)
v0 = [7.0 * np.cos(a), 7.0 * np.sin(a), 0.0]
paths = {}
for name, spin in (('backspin', [0.0, 0.0, 40.0]), ('none', [0.0, 0.0, 0.0]), ('topspin', [0.0, 0.0, -40.0])):
    p, _ = fly(v0, 1600, dt=0.001, spin=spin, position=[0.0, 1.0, 0.0])
    land = np.argmax(p[:, 1] < 0.0)
    paths[name] = (p[:land], p[land, 0])
assert paths['backspin'][1] > paths['none'][1] > paths['topspin'][1]
fig, ax = plt.subplots(figsize=(7, 3))
for name, (p, x) in paths.items():
    ax.plot(p[:, 0], p[:, 1], label=f'{name}: lands at {x:.2f} m')
ax.set_xlabel('x (m)'); ax.set_ylabel('height (m)'); ax.set_aspect('equal'); ax.legend()
plt.show()

a = np.radians(50.0) v0 = [7.0 * np.cos(a), 7.0 * np.sin(a), 0.0] paths = {} for name, spin in (('backspin', [0.0, 0.0, 40.0]), ('none', [0.0, 0.0, 0.0]), ('topspin', [0.0, 0.0, -40.0])): p, _ = fly(v0, 1600, dt=0.001, spin=spin, position=[0.0, 1.0, 0.0]) land = np.argmax(p[:, 1] < 0.0) paths[name] = (p[:land], p[land, 0]) assert paths['backspin'][1] > paths['none'][1] > paths['topspin'][1] fig, ax = plt.subplots(figsize=(7, 3)) for name, (p, x) in paths.items(): ax.plot(p[:, 0], p[:, 1], label=f'{name}: lands at {x:.2f} m') ax.set_xlabel('x (m)'); ax.set_ylabel('height (m)'); ax.set_aspect('equal'); ax.legend() plt.show()

Invariant: the lift does no work¶

With the drag off, a spinning lob's energy, kinetic and potential, stays where it started (to the integrator's error, which the step's square sets), though its path bends.

In [8]:

p, v = fly(v0, 1200, dt=0.001, spin=[0.0, 0.0, 40.0], drag_coefficient=0.0, position=[0.0, 1.0, 0.0])
e = 0.5 * np.sum(v * v, axis=1) + 9.80665 * p[:, 1]
e0 = 0.5 * 49.0 + 9.80665
q, _ = fly(v0, 1200, dt=0.001, drag_coefficient=0.0, position=[0.0, 1.0, 0.0])
drift = np.max(np.abs(e - e0)) / e0
assert drift < 1e-4 and np.max(np.abs(p - q)) > 0.01, (drift, np.max(np.abs(p - q)))
print(f'energy within {drift:.1e} of its start; the path bent by up to {np.max(np.abs(p - q)) * 100:.1f} cm')

p, v = fly(v0, 1200, dt=0.001, spin=[0.0, 0.0, 40.0], drag_coefficient=0.0, position=[0.0, 1.0, 0.0]) e = 0.5 * np.sum(v * v, axis=1) + 9.80665 * p[:, 1] e0 = 0.5 * 49.0 + 9.80665 q, _ = fly(v0, 1200, dt=0.001, drag_coefficient=0.0, position=[0.0, 1.0, 0.0]) drift = np.max(np.abs(e - e0)) / e0 assert drift < 1e-4 and np.max(np.abs(p - q)) > 0.01, (drift, np.max(np.abs(p - q))) print(f'energy within {drift:.1e} of its start; the path bent by up to {np.max(np.abs(p - q)) * 100:.1f} cm')

energy within 5.0e-06 of its start; the path bent by up to 1.6 cm

Limit: no spin, no lift¶

With no spin the force is exactly zero, whatever the motion.

In [9]:

got = trait_run.law('laws/magnus.trait', 'magnus', velocity=v, wind=np.zeros((len(v), 3)),
                    angular_velocity=np.zeros((len(v), 3)), density=1.204, lift_coefficient=0.2, area=area, mass=0.7)
assert np.all(got == 0.0)
print('no spin, no lift')

got = trait_run.law('laws/magnus.trait', 'magnus', velocity=v, wind=np.zeros((len(v), 3)), angular_velocity=np.zeros((len(v), 3)), density=1.204, lift_coefficient=0.2, area=area, mass=0.7) assert np.all(got == 0.0) print('no spin, no lift')

no spin, no lift