Skip to content

The Coriolis term, live

A horizontal velocity, u east and v north, is turned by the Coriolis parameter f = 2 Ω sin φ, Ω the planet's rotation rate and φ the latitude: du/dt = f v and dv/dt = −f u, to the right of the motion in the north, to the left in the south, nothing at the equator. Left to that term alone, a parcel of air turns in an inertial circle at its own speed, round once in 2π / f. Below, a parcel at 45° north is set moving north at 10 m/s and stepped by the library's velocity Verlet (lusine.integrator.kick_drift, the term lusine.laws.coriolis, lusine.integrator.kick), thirty seconds a step and two hundred steps a second, on your GPU as you read: a seventeen-hour circle in ten seconds. Beside it the same program computes the closed form of the circle, u = U sin(f t), v = U cos(f t), and checks the one against the other 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 parcel: the first chart draws its velocity north against its velocity east, a point a step: the circle its velocity goes round (drawn as an ellipse, the chart being wider than it is tall). The second draws its velocity east at each of the last 2100 steps, stepped and closed form one over the other. Once round the circle the parcel starts again. Below them, its time, the inertial period, and its velocity east and north.
  • The checks: the chart draws the stepped velocity's distance from the closed form's at each step as a share of what agrees allows, 1% of the speed, on a log scale, so the check holds while it stays below 1; the integrator's error grows over each circle, to about half of that. does no work holds while the term has stayed at right angles to the velocity, u du + v dv zero to the rounding of its two products, at every step: it turns the parcel and never speeds it up. hemispheres holds while an east wind is turned south in the north, north in the south, and not at all at the equator, the first two shown above it.

Every number here is computed by the program on the GPU; the page only reads the values back and draws them. The rotation rate (the Earth's sidereal one), the latitude, the speed 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 parcel, the closed form and the checks, the panels the page draws, and the one procedure, main, run each step: one velocity Verlet step, the term taken at the step's start and at its middle, and the checks after it.

// The Coriolis term's live page: a parcel of air set moving north and left to the Coriolis term alone, stepped by the
// library's velocity Verlet, turning in an inertial circle beside the closed form of it, with the page's checks
// computed here on the GPU and shown by the panels the program declares.

namespace coriolis;

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

// The reader's settings: the planet's rotation rate (the Earth's sidereal one), the parcel's latitude (45 degrees
// north), its speed at the start, and the step.
host struct {
  float<rad/s> omega = 0.000072921159 rad/s; float<rad> latitude = 0.78539816 rad; float<m/s> speed = 10.0 m/s;
  float<s> dt = 30.0 s;
} settings;

// The parcel, its velocity east (u) and north (v) and where it has gone, stepped from moving north at the speed; the
// Coriolis term's acceleration of it; its time the sum of the steps it has taken.
struct {
  float<m> x; float<m> y; float<m/s> u; float<m/s> v; float<m/s2> du; float<m/s2> dv; float<s> time;
} parcel;

// The closed form at the parcel's time, an inertial circle, u = U sin(f t), v = U cos(f t) with f = 2 Omega sin(phi),
// and the period it closes in, 2 pi / f; the term at the parcel's velocity, and for a wind east at the speed in the
// north, in the south and at the equator.
struct {
  float<m/s> u; float<m/s> v; float<s> period; float<m/s2> du; float<m/s2> dv;
  float<m/s2> north; float<m/s2> south; float<m/s2> equator_u; float<m/s2> equator_v;
} closed;

// The checks: the stepped velocity within 1% of the speed of the closed form at every step (its miss also as a share
// of that, which the page charts); the term doing no work, u du + v dv zero to the rounding of its two products at
// every step; and the hemispheres, an east wind turned south in the north, north in the south, and not at all at the
// equator.
struct { float error_share; bool agrees = true; bool no_work = true; bool hemispheres = 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<s> inertial_period(float<rad/s> omega, float<rad> latitude)
{
  return 6.2831853 / (2.0 * omega * sin(latitude));
}

(float<m/s>, float<m/s>) circling(float<m/s> speed, float<rad/s> omega, float<rad> latitude, float<s> time)
{
  let turned = 2.0 * omega * sin(latitude) * time;
  return (speed * sin(turned), speed * cos(turned));
}

// How far the stepped velocity is from the closed form's, as a share of the 1% of the speed the check allows.
float miss(float<m/s> u, float<m/s> v, float<m/s> closed_u, float<m/s> closed_v, float<m/s> speed)
{
  let eu = u - closed_u;
  let ev = v - closed_v;
  return sqrt(eu * eu + ev * ev) / (0.01 * speed);
}

// Whether the work the term does on the velocity, u du + v dv, is within four roundings of f (u^2 + v^2) of zero.
bool workless(float<m/s> u, float<m/s> v, float<m/s2> du, float<m/s2> dv, float<rad/s> omega, float<rad> latitude)
{
  return abs(u * du + v * dv) <= 0.00000048 * abs(2.0 * omega * sin(latitude)) * (u * u + v * v);
}

output table<lusine.widget> panels = {
  lusine.series { { panel = "the parcel", chart = "its velocity", quantity = parcel.v, steps = 2100,
                    against = parcel.u, label = "north against east" } },
  lusine.series { { panel = "the parcel", chart = "eastward", quantity = parcel.u, steps = 2100, label = "stepped" } },
  lusine.series { { panel = "the parcel", chart = "eastward", quantity = closed.u, steps = 2100,
                    label = "closed form" } },
  lusine.value { { panel = "the parcel", quantity = parcel.time, label = "time" } },
  lusine.value { { panel = "the parcel", quantity = closed.period, label = "inertial period" } },
  lusine.value { { panel = "the parcel", quantity = parcel.u, label = "east" } },
  lusine.value { { panel = "the parcel", quantity = parcel.v, label = "north" } },
  lusine.series { { panel = "the checks", chart = "share of the tolerance", quantity = checks.error_share, steps = 2100,
                    log = true, label = "velocity error" } },
  lusine.value { { panel = "the checks", quantity = closed.north, label = "east wind, north" } },
  lusine.value { { panel = "the checks", quantity = closed.south, label = "east wind, south" } },
  lusine.value { { panel = "the checks", quantity = checks.agrees, label = "agrees" } },
  lusine.value { { panel = "the checks", quantity = checks.no_work, label = "does no work" } },
  lusine.value { { panel = "the checks", quantity = checks.hemispheres, label = "hemispheres" } },
};

void main()
{
  // Once round the circle the parcel starts again, moving north at the speed.
  closed.period = inertial_period(settings.omega, settings.latitude);
  parcel.u = parcel.time >= closed.period ? 0.0 m/s : parcel.u;
  parcel.v = parcel.time == 0.0 s || parcel.time >= closed.period ? settings.speed : parcel.v;
  parcel.x = parcel.time >= closed.period ? 0.0 m : parcel.x;
  parcel.y = parcel.time >= closed.period ? 0.0 m : parcel.y;
  parcel.time = parcel.time >= closed.period ? 0.0 s : parcel.time;
  (parcel.du, parcel.dv) = lusine.laws.coriolis(parcel.u, parcel.v, settings.latitude, settings.omega);
  (parcel.x, parcel.u) = lusine.integrator.kick_drift(position = parcel.x, velocity = parcel.u,
                                                      acceleration = parcel.du, dt = settings.dt);
  (parcel.y, parcel.v) = lusine.integrator.kick_drift(position = parcel.y, velocity = parcel.v,
                                                      acceleration = parcel.dv, dt = settings.dt);
  (parcel.du, parcel.dv) = lusine.laws.coriolis(parcel.u, parcel.v, settings.latitude, settings.omega);
  parcel.u = lusine.integrator.kick(velocity = parcel.u, acceleration = parcel.du, dt = settings.dt);
  parcel.v = lusine.integrator.kick(velocity = parcel.v, acceleration = parcel.dv, dt = settings.dt);
  parcel.time = parcel.time + settings.dt;
  (closed.u, closed.v) = circling(settings.speed, settings.omega, settings.latitude, time = parcel.time);
  checks.error_share = miss(parcel.u, parcel.v, closed_u = closed.u, closed_v = closed.v, settings.speed);
  checks.agrees = checks.agrees && checks.error_share <= 1.0;
  (closed.du, closed.dv) = lusine.laws.coriolis(parcel.u, parcel.v, settings.latitude, settings.omega);
  checks.no_work = checks.no_work && workless(parcel.u, parcel.v, du = closed.du, dv = closed.dv, settings.omega,
                                              settings.latitude);
  (closed.equator_u, closed.north) = lusine.laws.coriolis(u = settings.speed, v = 0.0 m/s, settings.latitude,
                                                           settings.omega);
  (closed.equator_u, closed.south) = lusine.laws.coriolis(u = settings.speed, v = 0.0 m/s,
                                                           latitude = -settings.latitude, settings.omega);
  (closed.equator_u, closed.equator_v) = lusine.laws.coriolis(u = settings.speed, v = 0.0 m/s, latitude = 0.0 rad,
                                                               settings.omega);
  checks.hemispheres = checks.hemispheres && closed.north < 0.0 m/s2 && closed.south > 0.0 m/s2
                       && closed.equator_u == 0.0 m/s2 && closed.equator_v == 0.0 m/s2;
  lamps.colour = lamp_colour(holds = checks.agrees && checks.no_work && checks.hemispheres);
}

The law's other checks (the deflection against 2 Ω sin(φ) u over every latitude) are on its proof page.