diff --git a/specs/t27f/particle_engine.t27 b/specs/t27f/particle_engine.t27 new file mode 100644 index 0000000000..72835332ed --- /dev/null +++ b/specs/t27f/particle_engine.t27 @@ -0,0 +1,1000 @@ +// SPDX-License-Identifier: Apache-2.0 +// specs/t27f/particle_engine.t27 -- the particle engine of T27F Live: real forces move sixteen +// particles, the agent's events steer them through a controller, and the picture is derived from +// the physical state. Nothing on the screen is drawn by hand; every pixel is a function of +// position, velocity, mass and charge. +// +// WHERE IT SITS: specs/t27f/particle_field.t27 is the conversation as a three-node graph (user, +// agent, tool). This spec is the matter that the graph is made of: sixteen bodies, one mesh-shader +// meshlet (16 quads = 64 vertices), simulated and then flattened into a fixed display list, the +// "display list -> paint" step of docs/nona-02-organism/t27f/genui-research.md section 5. The host +// draws the list as given; it holds no layout, no colour and no state. +// +// THE PHYSICS, each force with its law: +// gravity F = m g, g = 9.80665 m/s^2 (standard gravity, exact by definition, 3rd CGPM 1901) +// Coulomb F = KC qi qj r / (r^2 + eps^2)^(3/2), a Plummer-softened inverse square +// (H. C. Plummer, MNRAS 71, 460, 1911), charges are trits -1, 0, +1 +// body contact F = K d^(3/2) - c v_n while two spheres overlap by d: the Hertz law for +// elastic spheres (H. Hertz, J. reine angew. Math. 92, 156, 1882) with a dashpot, +// as in the discrete element method of P. A. Cundall and O. D. L. Strack +// (Geotechnique 29, 47, 1979). Bodies do not pass through each other, opposite +// charges bind into touching pairs, and a pile comes to rest. Two models were +// tried and measured first: a Lennard-Jones core soft enough not to swallow the +// controller let a falling body sink into the r^-12 wall, where omega H > 2 and +// Verlet diverges; a linear spring lost 10 to 20 percent of a bounce to the step +// where its force switches on. Hertz switches on smoothly and loses 0.5 percent. +// Stokes drag F = -b v, the air the particles float in +// floor, walls reflect a body with restitution e (a lacquered floor); an impact slower than +// V_REST is absorbed, and a body resting on the floor gets the normal force that +// cancels whatever pushes it down +// THE CONTROL: the agent never sets a position. An event changes a command (a ring centre, a ring +// radius, a spin, gravity on or off), and each body follows its slot on the ring through a +// critically damped PD controller, a = KP (target - x) - KD v with KD = 2 sqrt(KP), the fastest +// approach that never overshoots (K. Ogata, "Modern Control Engineering", 5th ed., 2010, ch. 5). +// IDLE ring at the agent, radius 0.25 m +// USER_MSG the ring tightens to 0.16 m: something came in +// THINK the ring spins at 1.2 rad/s +// TOOL_CALL the ring moves out to the tool and tightens +// TOOL_RESULT the ring comes back +// AGENT_MSG the ring moves down toward the user +// ERROR an outward radial kick of 1.2 m/s, then IDLE +// DONE control off, gravity on: the matter falls to the black floor and comes to rest +// +// THE NUMERICS: eight substeps of H = 1/480 s per displayed frame. Each substep is a Strang splitting: +// half a step of exact exponential damping, then velocity Verlet for every position-dependent force, +// then the other half of the damping. Verlet is symplectic and exact for a constant force; the +// exponential is exact for drag; the split is second order (G. Strang, SIAM J. Numer. Anal. 5, 506, +// 1968). Forces act pair by pair and are applied with opposite signs, so momentum is conserved to +// rounding. No clock and no random number: the same log gives the same frames bit for bit. +// +// THE PICTURE, derived from the state: +// glow a gaussian disc at the body; radius r = R_SCALE m^(1/3) (a sphere of fixed +// density); brightness 1 + v^2 / V_REF^2, so fast matter is hot; hue by charge +// (+1 gold 45, 0 violet 275, -1 cyan 190, the palette of particle_field.t27) +// reflection the mirror image in the black floor (y -> -y). Its strength is the Fresnel +// reflectance of lacquer, n = 1.5, in Schlick's form R = R0 + (1 - R0)(1 - cos t)^5, +// R0 = ((n - 1)/(n + 1))^2 = 0.04 (C. Schlick, Comput. Graph. Forum 13, 233, 1994), +// at the angle a camera 0.9 m high and 2.5 m away sees it: low bodies mirror +// strongly, high ones faintly, as on a real glossy floor. A glossy floor blurs with +// distance, so the image widens by ROUGH per metre of height. +// streak the path of the last TAU seconds, x - v TAU to x: motion blur from the velocity +// background true black. Light adds, and the exposure 1 - exp(-L) never clips, so a pixel far +// from every body is 0 of 255 and a crowded one saturates smoothly. +// Reduced motion (the OS setting) keeps the formations and the fall, and drops the spin, the kick +// and the streaks. +// +// WHAT IT DOES NOT CLAIM: the units are SI at a tabletop scale, but the charges, masses and +// strengths are design choices, labelled as such below; this is not a model of any real material. +// The reference raster here draws glows and reflections, not streaks. Writing frames to a file and +// the browser host are not in this spec; both are t27c work (see the PR). +// phi^2 + 1/phi^2 = 3 | TRINITY + +module particle_engine { + + pub const KIND : str = "widget-logic"; + pub const ID : str = "particle-engine"; + pub const VERSION : u8 = 1; + + pub const PI : f64 = 3.141592653589793; + pub const N : usize = 16; // one meshlet of quads + pub const FRAME_DT : f64 = 0.016666666666666666; + pub const SUB : i64 = 8; + pub const H : f64 = 0.0020833333333333333; // 1/480 s + + // Physical constant. + pub const G_STD : f64 = 9.80665; + + // Design constants (scene scale, chosen, not measured). + pub const W : f64 = 0.5; // half width of the box, m + pub const TOP : f64 = 1.6; // ceiling, m + pub const R_SCALE : f64 = 0.022; // radius of a 1 kg body, m + pub const KC : f64 = 0.005; // Coulomb strength, N m^2 + pub const EPS_C2 : f64 = 0.0025; // Plummer softening squared, m^2 + pub const K_CONTACT : f64 = 100000.0; // Hertz stiffness, N/m^1.5: 7.6 mm overlap at 1 m/s + pub const C_CONTACT : f64 = 20.0; // contact damping, N s/m + pub const V_REST : f64 = 0.14; // sqrt(2 g 1 mm): slower impacts do not bounce + pub const CONTACT : f64 = 0.0001; // a body this close to the floor rests on it + pub const DRAG : f64 = 0.6; // Stokes coefficient, kg/s + pub const RESTITUTION : f64 = 0.7; + pub const KP : f64 = 16.0; // controller stiffness, 1/s^2 (omega = 4 rad/s) + pub const KD : f64 = 8.0; // 2 sqrt(KP): critical + pub const SPIN : f64 = 1.2; // rad/s while thinking + pub const KICK : f64 = 1.2; // m/s, the error burst + pub const AGENT_X : f64 = 0.0; + pub const AGENT_Y : f64 = 0.9; + pub const TOOL_X : f64 = 0.25; + pub const TOOL_Y : f64 = 1.15; + pub const USER_Y : f64 = 0.55; + pub const RING : f64 = 0.25; + + // Optics (design: camera placement; physics: Fresnel). + pub const N_LACQUER : f64 = 1.5; + pub const CAM_H : f64 = 0.9; + pub const CAM_D : f64 = 2.5; + pub const ROUGH : f64 = 0.08; + pub const V_REF : f64 = 1.0; + pub const TAU : f64 = 0.05; // streak length in seconds of travel + pub const HORIZON : f64 = 0.06; // the faint glint where the floor meets the dark + + // Raster: 108 x 192, one tenth of a 1080 x 1920 reel frame. + pub const COLS : i64 = 108; + pub const ROWS : i64 = 192; + pub const PIX : f64 = 0.009259259259259259; // 2 W / COLS, m per pixel + pub const Y_TOP : f64 = 1.33; + + // Event kinds, the same numbers as particle_field.t27. + pub const USER_MSG : i64 = 0; + pub const THINK : i64 = 1; + pub const TOOL_CALL : i64 = 2; + pub const TOOL_RESULT : i64 = 3; + pub const AGENT_MSG : i64 = 4; + pub const ERROR : i64 = 5; + pub const DONE : i64 = 6; + + // Display list: 2 + 3 N ops of 8 fields: op, x, y, r, hue, alpha, x2, y2. + pub const OP_CLEAR : i64 = 1; + pub const OP_FLOOR : i64 = 2; + pub const OP_REFLECT : i64 = 3; + pub const OP_STREAK : i64 = 4; + pub const OP_GLOW : i64 = 5; + pub const OPS : i64 = 50; + pub const OP_BUDGET : i64 = 64; + + // Scenarios. + pub const SC_SHOW : i64 = 0; // the demo conversation + pub const SC_ERRSHOW : i64 = 1; // the same with an error in place of the tool result + pub const SC_FALL : i64 = 2; // one body, gravity only, no floor + pub const SC_BOUNCE : i64 = 3; // one body dropped on the floor, no drag + pub const SC_DRAG : i64 = 4; // one body coasting in air + pub const SC_TERMINAL : i64 = 5; // one body falling through air, no floor + pub const SC_NBODY : i64 = 6; // sixteen bodies, Coulomb + elastic contact, free space + pub const SC_DIMER : i64 = 7; // a +1 and a -1 body settling into a bond + pub const SC_PD : i64 = 8; // one body steered to a target 0.3 m away + pub const SC_SPIN : i64 = 9; // sixteen bodies thinking from frame 0 + pub const SC_COULOMB : i64 = 10; // SC_NBODY with the contact off: a smooth force, for energy + pub const SC_COLLIDE : i64 = 11; // two neutral 1 kg bodies meeting head-on at 1 m/s, no air + + // Metrics. + pub const M_X : i64 = 0; + pub const M_Y : i64 = 1; + pub const M_VX : i64 = 2; + pub const M_VY : i64 = 3; + pub const M_KE : i64 = 4; + pub const M_PX : i64 = 5; + pub const M_PY : i64 = 6; + pub const M_ENERGY : i64 = 7; + pub const M_APEX : i64 = 8; // highest y after the first floor contact + pub const M_DIST : i64 = 9; // distance between bodies 0 and 1 + pub const M_LZ : i64 = 10; // angular momentum about the agent + pub const M_OP : i64 = 11; // field `field` of op `idx` + pub const M_HASH : i64 = 12; // hash of the whole display list + pub const M_PIXEL : i64 = 13; // 0..255 at pixel idx = row * COLS + col + pub const M_MAXSPEED : i64 = 14; + pub const M_MAXX : i64 = 15; // largest x of body 0 over the run + pub const M_MAXY : i64 = 16; // largest y of any body + pub const M_MAXREST : i64 = 17; // largest (y - r) over the bodies + pub const M_COUNT : i64 = 18; // number of ops in the display list + + // ---- bodies ---- + + fn mass_of(i: i64) -> f64 { + return 1.0 + 0.5 * (@mod(i, 3) as f64); + } + + fn charge_of(i: i64) -> f64 { + if (i == 15) { + return 0.0; + } + const k = @mod(i, 3); + if (k == 0) { return 1.0; } + if (k == 1) { return 0.0; } + return 0.0 - 1.0; + } + + fn net_charge() -> f64 { + var s : f64 = 0.0; + var i : i64 = 0; + while (i < 16) { + s = s + charge_of(i); + i = i + 1; + } + return s; + } + + fn radius_of(m: f64) -> f64 { + return R_SCALE * @exp(@log(m) / 3.0); + } + + fn hue_of(q: f64) -> f64 { + if (q > 0.5) { return 45.0; } + if (q < 0.0 - 0.5) { return 190.0; } + return 275.0; + } + + // ---- the command an event sets ---- + + fn cmd_x(kind: i64) -> f64 { + if (kind == TOOL_CALL) { return TOOL_X; } + return AGENT_X; + } + + fn cmd_y(kind: i64) -> f64 { + if (kind == TOOL_CALL) { return TOOL_Y; } + if (kind == AGENT_MSG) { return USER_Y; } + return AGENT_Y; + } + + fn cmd_r(kind: i64) -> f64 { + if (kind == USER_MSG || kind == TOOL_CALL) { return 0.16; } + if (kind == AGENT_MSG) { return 0.18; } + return RING; + } + + fn cmd_spin(kind: i64) -> f64 { + if (kind == THINK) { return SPIN; } + return 0.0; + } + + fn cmd_control(kind: i64) -> f64 { + if (kind == DONE) { return 0.0; } + return 1.0; + } + + fn cmd_gravity(kind: i64) -> f64 { + if (kind == DONE) { return 1.0; } + return 0.0; + } + + // ---- the demo log ---- + + pub const LOG_LEN : i64 = 6; + + fn log_kind(i: i64, with_error: bool) -> i64 { + if (i == 0) { return USER_MSG; } + if (i == 1) { return THINK; } + if (i == 2) { return TOOL_CALL; } + if (i == 3) { + if (with_error) { return ERROR; } + return TOOL_RESULT; + } + if (i == 4) { return AGENT_MSG; } + return DONE; + } + + fn log_frame(i: i64) -> i64 { + if (i == 0) { return 120; } + if (i == 1) { return 150; } + if (i == 2) { return 210; } + if (i == 3) { return 260; } + if (i == 4) { return 300; } + return 330; + } + + // ---- optics ---- + + fn fresnel_r0(n: f64) -> f64 { + const a = (n - 1.0) / (n + 1.0); + return a * a; + } + + // Schlick reflectance at an angle of incidence whose cosine is c. + fn schlick(c: f64, n: f64) -> f64 { + const r0 = fresnel_r0(n); + const k = 1.0 - c; + return r0 + (1.0 - r0) * k * k * k * k * k; + } + + // Cosine of the incidence angle at which the camera sees, in the floor, a body at height y. + fn mirror_cos(y: f64) -> f64 { + const d = CAM_D * CAM_H / (CAM_H + y); + return CAM_H / @sqrt(CAM_H * CAM_H + d * d); + } + + fn mirror_alpha(y: f64, glow: f64) -> f64 { + return glow * schlick(mirror_cos(y), N_LACQUER); + } + + fn glow_of(vx: f64, vy: f64) -> f64 { + return 1.0 + (vx * vx + vy * vy) / (V_REF * V_REF); + } + + fn gauss(d2: f64, r: f64) -> f64 { + return @exp(0.0 - d2 / (2.0 * r * r)); + } + + // Exposure: light adds, the film never clips. + fn expose(l: f64) -> f64 { + return 255.0 * (1.0 - @exp(0.0 - l)); + } + + // The sum of a unit gaussian of width sigma pixels over the pixel grid, divided by 2 pi sigma^2: + // a glow keeps its light when it is rasterised. + fn kernel_sum(sigma: f64) -> f64 { + var s : f64 = 0.0; + var i : i64 = 0 - 40; + while (i <= 40) { + var j : i64 = 0 - 40; + while (j <= 40) { + const di = i as f64; + const dj = j as f64; + s = s + gauss(di * di + dj * dj, sigma); + j = j + 1; + } + i = i + 1; + } + return s / (2.0 * PI * sigma * sigma); + } + + // ---- contact and Coulomb, as force magnitude over r and as potential ---- + + // Hertz push between two elastic spheres whose centres are sqrt(r2) apart and touch at d0: + // F = K d^(3/2), U = (2/5) K d^(5/2) for an overlap d (H. Hertz, J. reine angew. Math. 92, 156, 1882). + fn contact_f_over_r(r2: f64, d0: f64) -> f64 { + if (r2 >= d0 * d0) { + return 0.0; + } + const r = @sqrt(r2); + const d = d0 - r; + return K_CONTACT * d * @sqrt(d) / r; + } + + fn contact_u(r2: f64, d0: f64) -> f64 { + if (r2 >= d0 * d0) { + return 0.0; + } + const d = d0 - @sqrt(r2); + return 0.4 * K_CONTACT * d * d * @sqrt(d); + } + + fn coulomb_f_over_r(qq: f64, r2: f64) -> f64 { + const s = r2 + EPS_C2; + return KC * qq / (s * @sqrt(s)); + } + + fn coulomb_u(qq: f64, r2: f64) -> f64 { + return KC * qq / @sqrt(r2 + EPS_C2); + } + + // The radial force on a +1 body from a -1 body of the same size, r apart (positive = apart). + fn dimer_force(r: f64) -> f64 { + const d0 = 2.0 * radius_of(1.0); + return contact_f_over_r(r * r, d0) * r + coulomb_f_over_r(0.0 - 1.0, r * r) * r; + } + + // Where the Coulomb pull and the contact push balance: Newton's method with a backward + // difference, so the slope is always taken inside the contact. + fn dimer_newton() -> f64 { + var r : f64 = 0.99 * 2.0 * radius_of(1.0); + var i : i64 = 0; + while (i < 60) { + const g = dimer_force(r); + const rh = r * 0.999999; + const gh = dimer_force(rh); + if (gh != g) { + r = r - g * (rh - r) / (gh - g); + } + i = i + 1; + } + return r; + } + + fn fnv(h: i64, v: f64) -> i64 { + const q = @as(i64, @intFromFloat(@round(v * 1000000.0))); + return @mod(h * 16777619 + q, 2147483647); + } + + // ---- the engine ---- + + // Run scenario `sc` for `frames` frames and report metric `metric` (with `idx` and `field` where + // the metric needs them) of the final state. + fn world(sc: i64, frames: i64, reduced: bool, metric: i64, idx: i64, field: i64) -> f64 { + var px : [16]f64 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]; + var py : [16]f64 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]; + var vx : [16]f64 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]; + var vy : [16]f64 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]; + var ax : [16]f64 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]; + var ay : [16]f64 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]; + var ms : [16]f64 = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]; + var qs : [16]f64 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]; + var rs : [16]f64 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]; + + // The switches a scenario sets. + var n : usize = 16; + var drag : f64 = DRAG; + var grav : f64 = 0.0; + var ctrl : f64 = 1.0; + var floor_on : bool = true; + var walls_on : bool = true; + var pair_on : bool = true; + var cdamp : bool = true; // contact dashpot; off where energy must be conserved + var contact_on : bool = true; + var log_on : bool = false; + var with_error : bool = false; + var cx : f64 = AGENT_X; + var cy : f64 = AGENT_Y; + var ring : f64 = RING; + var spin : f64 = 0.0; + var phase : f64 = 0.0; + var tx0 : f64 = 0.0; // a fixed target for SC_PD + var ty0 : f64 = 0.0; + var fixed_target : bool = false; + + // Default bodies: masses, charges and radii by index, placed on the ring at rest. + var i : usize = 0; + var ii : i64 = 0; + while (i < 16) { + ms[i] = mass_of(ii); + qs[i] = charge_of(ii); + rs[i] = radius_of(ms[i]); + const ang = 2.0 * PI * (ii as f64) / 16.0; + px[i] = AGENT_X + RING * @cos(ang); + py[i] = AGENT_Y + RING * @sin(ang); + i = i + 1; + ii = ii + 1; + } + + if (sc == SC_SHOW || sc == SC_ERRSHOW) { + log_on = true; + with_error = sc == SC_ERRSHOW; + } + if (sc == SC_FALL || sc == SC_BOUNCE || sc == SC_DRAG || sc == SC_TERMINAL || sc == SC_PD) { + n = 1; + ms[0] = 1.0; + rs[0] = radius_of(1.0); + px[0] = 0.0; + py[0] = 1.0; + pair_on = false; + walls_on = false; + ctrl = 0.0; + } + if (sc == SC_FALL) { + drag = 0.0; + grav = 1.0; + floor_on = false; + } + if (sc == SC_BOUNCE) { + drag = 0.0; + grav = 1.0; + py[0] = 0.5 + rs[0]; + } + if (sc == SC_DRAG) { + vx[0] = 1.0; + floor_on = false; + } + if (sc == SC_TERMINAL) { + grav = 1.0; + floor_on = false; + } + if (sc == SC_PD) { + drag = 0.0; + ctrl = 1.0; + fixed_target = true; + tx0 = 0.3; + ty0 = 1.0; + } + if (sc == SC_NBODY || sc == SC_COULOMB) { + drag = 0.0; + cdamp = false; + ctrl = 0.0; + floor_on = false; + walls_on = false; + } + if (sc == SC_DIMER) { + n = 2; + ctrl = 0.0; + floor_on = false; + walls_on = false; + ms[0] = 1.0; + ms[1] = 1.0; + qs[0] = 1.0; + qs[1] = 0.0 - 1.0; + rs[0] = radius_of(1.0); + rs[1] = radius_of(1.0); + px[0] = 0.0 - 0.05; + py[0] = 1.0; + px[1] = 0.05; + py[1] = 1.0; + } + if (sc == SC_COULOMB) { + contact_on = false; + } + if (sc == SC_COLLIDE) { + n = 2; + drag = 0.0; + cdamp = false; + ctrl = 0.0; + floor_on = false; + walls_on = false; + ms[0] = 1.0; + ms[1] = 1.0; + qs[0] = 0.0; + qs[1] = 0.0; + rs[0] = radius_of(1.0); + rs[1] = radius_of(1.0); + px[0] = 0.0 - 0.1; + py[0] = 1.0; + px[1] = 0.1; + py[1] = 1.0; + vx[0] = 0.5; + vx[1] = 0.0 - 0.5; + } + if (sc == SC_SPIN) { + spin = SPIN; + } + if (reduced) { + spin = 0.0; + } + + var apex : f64 = 0.0; + var touched : bool = false; + var maxspeed : f64 = 0.0; + var maxx : f64 = px[0]; + var maxy : f64 = 0.0; + + var f : i64 = 0; + while (f < frames) { + // Events change the command, never a position. + if (log_on) { + var e : i64 = 0; + while (e < LOG_LEN) { + if (log_frame(e) == f) { + const kind = log_kind(e, with_error); + if (kind == ERROR) { + if (reduced == false) { + var k : usize = 0; + while (k < n) { + const dx = px[k] - cx; + const dy = py[k] - cy; + const d = @sqrt(dx * dx + dy * dy) + 0.000001; + vx[k] = vx[k] + KICK * dx / d; + vy[k] = vy[k] + KICK * dy / d; + k = k + 1; + } + } + cx = AGENT_X; + cy = AGENT_Y; + ring = RING; + spin = 0.0; + } else { + cx = cmd_x(kind); + cy = cmd_y(kind); + ring = cmd_r(kind); + spin = cmd_spin(kind); + ctrl = cmd_control(kind); + grav = cmd_gravity(kind); + if (reduced) { + spin = 0.0; + } + } + } + e = e + 1; + } + } + + var s : i64 = 0; + while (s <= SUB) { + // Forces at the current positions. + var a : usize = 0; + var ai : i64 = 0; + while (a < n) { + var tx : f64 = cx + ring * @cos(phase + 2.0 * PI * (ai as f64) / 16.0); + var ty : f64 = cy + ring * @sin(phase + 2.0 * PI * (ai as f64) / 16.0); + if (fixed_target) { + tx = tx0; + ty = ty0; + } + ax[a] = ctrl * KP * (tx - px[a]); + ay[a] = ctrl * KP * (ty - py[a]) - grav * G_STD; + a = a + 1; + ai = ai + 1; + } + if (pair_on) { + var p : usize = 0; + while (p < n) { + var r : usize = p + 1; + while (r < n) { + const dx = px[p] - px[r]; + const dy = py[p] - py[r]; + const r2 = dx * dx + dy * dy; + const d0 = rs[p] + rs[r]; + var fr : f64 = coulomb_f_over_r(qs[p] * qs[r], r2); + if (contact_on) { + fr = fr + contact_f_over_r(r2, d0); + } + if (contact_on && cdamp && r2 < d0 * d0) { + const vn = (vx[p] - vx[r]) * dx + (vy[p] - vy[r]) * dy; + fr = fr - C_CONTACT * vn / r2; + } + ax[p] = ax[p] + fr * dx / ms[p]; + ay[p] = ay[p] + fr * dy / ms[p]; + ax[r] = ax[r] - fr * dx / ms[r]; + ay[r] = ay[r] - fr * dy / ms[r]; + r = r + 1; + } + p = p + 1; + } + } + // A body resting on the floor: the normal force cancels what pushes it down. + if (floor_on) { + var q : usize = 0; + while (q < n) { + if (py[q] <= rs[q] + CONTACT && ay[q] < 0.0) { + ay[q] = 0.0; + } + q = q + 1; + } + } + // Finish the previous substep: half kick, half damping. + if (s > 0) { + var b : usize = 0; + while (b < n) { + const dh = @exp(0.0 - 0.5 * H * (drag / ms[b] + ctrl * KD)); + vx[b] = (vx[b] + 0.5 * H * ax[b]) * dh; + vy[b] = (vy[b] + 0.5 * H * ay[b]) * dh; + b = b + 1; + } + } + // Start the next one: half damping, half kick, drift, contact. + if (s < SUB) { + var c : usize = 0; + while (c < n) { + const dh = @exp(0.0 - 0.5 * H * (drag / ms[c] + ctrl * KD)); + vx[c] = vx[c] * dh + 0.5 * H * ax[c]; + vy[c] = vy[c] * dh + 0.5 * H * ay[c]; + px[c] = px[c] + H * vx[c]; + py[c] = py[c] + H * vy[c]; + if (floor_on && py[c] < rs[c]) { + py[c] = 2.0 * rs[c] - py[c]; + vy[c] = 0.0 - RESTITUTION * vy[c]; + if (vy[c] < V_REST) { + vy[c] = 0.0; + } + touched = true; + } + if (walls_on) { + if (py[c] > TOP - rs[c]) { + py[c] = 2.0 * (TOP - rs[c]) - py[c]; + vy[c] = 0.0 - RESTITUTION * vy[c]; + } + if (px[c] > W - rs[c]) { + px[c] = 2.0 * (W - rs[c]) - px[c]; + vx[c] = 0.0 - RESTITUTION * vx[c]; + } + if (px[c] < rs[c] - W) { + px[c] = 2.0 * (rs[c] - W) - px[c]; + vx[c] = 0.0 - RESTITUTION * vx[c]; + } + } + c = c + 1; + } + } + s = s + 1; + } + phase = phase + spin * FRAME_DT; + + var t : usize = 0; + while (t < n) { + const sp = @sqrt(vx[t] * vx[t] + vy[t] * vy[t]); + if (sp > maxspeed) { maxspeed = sp; } + if (py[t] > maxy) { maxy = py[t]; } + t = t + 1; + } + if (touched && py[0] > apex) { apex = py[0]; } + if (px[0] > maxx) { maxx = px[0]; } + f = f + 1; + } + + // ---- metrics of the final state ---- + const i0 = @as(usize, @intCast(idx)); + if (metric == M_X) { return px[i0]; } + if (metric == M_Y) { return py[i0]; } + if (metric == M_VX) { return vx[i0]; } + if (metric == M_VY) { return vy[i0]; } + if (metric == M_APEX) { return apex; } + if (metric == M_MAXSPEED) { return maxspeed; } + if (metric == M_MAXX) { return maxx; } + if (metric == M_MAXY) { return maxy; } + if (metric == M_DIST) { + const dx = px[0] - px[1]; + const dy = py[0] - py[1]; + return @sqrt(dx * dx + dy * dy); + } + var ke : f64 = 0.0; + var pxs : f64 = 0.0; + var pys : f64 = 0.0; + var lz : f64 = 0.0; + var u : f64 = 0.0; + var rest : f64 = 0.0 - 1.0; + var k : usize = 0; + while (k < n) { + ke = ke + 0.5 * ms[k] * (vx[k] * vx[k] + vy[k] * vy[k]); + pxs = pxs + ms[k] * vx[k]; + pys = pys + ms[k] * vy[k]; + lz = lz + ms[k] * ((px[k] - AGENT_X) * vy[k] - (py[k] - AGENT_Y) * vx[k]); + u = u + grav * ms[k] * G_STD * py[k]; + if (py[k] - rs[k] > rest) { rest = py[k] - rs[k]; } + var j : usize = k + 1; + while (j < n) { + const dx = px[k] - px[j]; + const dy = py[k] - py[j]; + const r2 = dx * dx + dy * dy; + if (contact_on) { + u = u + contact_u(r2, rs[k] + rs[j]); + } + u = u + coulomb_u(qs[k] * qs[j], r2); + j = j + 1; + } + k = k + 1; + } + if (metric == M_KE) { return ke; } + if (metric == M_PX) { return pxs; } + if (metric == M_PY) { return pys; } + if (metric == M_LZ) { return lz; } + if (metric == M_ENERGY) { return ke + u; } + if (metric == M_MAXREST) { return rest; } + + // ---- the picture ---- + if (metric == M_PIXEL) { + const row = @divTrunc(idx, COLS); + const col = idx - row * COLS; + const wx = 0.0 - W + ((col as f64) + 0.5) * PIX; + const wy = Y_TOP - ((row as f64) + 0.5) * PIX; + var l : f64 = 0.0; + var g : usize = 0; + while (g < n) { + const glow = glow_of(vx[g], vy[g]); + const dx = wx - px[g]; + const dy = wy - py[g]; + l = l + glow * gauss(dx * dx + dy * dy, rs[g]); + if (wy < 0.0) { + const ry = wy + py[g]; + const rr = rs[g] + ROUGH * py[g]; + l = l + mirror_alpha(py[g], glow) * gauss(dx * dx + ry * ry, rr); + } + g = g + 1; + } + return expose(l); + } + + var h : i64 = 2166136261; + var count : i64 = 0; + var o : i64 = 0; + while (o < OPS) { + var v0 : f64 = 0.0; + var v1 : f64 = 0.0; + var v2 : f64 = 0.0; + var v3 : f64 = 0.0; + var v4 : f64 = 0.0 - 1.0; + var v5 : f64 = 0.0; + var v6 : f64 = 0.0; + var v7 : f64 = 0.0; + if (o == 0) { + v0 = OP_CLEAR as f64; + } + if (o == 1) { + v0 = OP_FLOOR as f64; + v1 = 0.0 - W; + v5 = HORIZON; + v6 = W; + } + if (o >= 2) { + const which = @divTrunc(o - 2, 3); + const part = o - 2 - 3 * which; + const b = @as(usize, @intCast(which)); + const glow = glow_of(vx[b], vy[b]); + v4 = hue_of(qs[b]); + if (part == 0) { + v0 = OP_REFLECT as f64; + v1 = px[b]; + v2 = 0.0 - py[b]; + v3 = rs[b] + ROUGH * py[b]; + v5 = mirror_alpha(py[b], glow); + } + if (part == 1) { + v0 = OP_STREAK as f64; + v1 = px[b] - vx[b] * TAU; + v2 = py[b] - vy[b] * TAU; + v3 = 0.5 * rs[b]; + if (reduced == false) { + v5 = 0.5 * @min(1.0, @sqrt(vx[b] * vx[b] + vy[b] * vy[b]) / V_REF); + } + v6 = px[b]; + v7 = py[b]; + } + if (part == 2) { + v0 = OP_GLOW as f64; + v1 = px[b]; + v2 = py[b]; + v3 = rs[b]; + v5 = glow; + } + } + if (metric == M_OP && o == idx) { + if (field == 0) { return v0; } + if (field == 1) { return v1; } + if (field == 2) { return v2; } + if (field == 3) { return v3; } + if (field == 4) { return v4; } + if (field == 5) { return v5; } + if (field == 6) { return v6; } + return v7; + } + h = fnv(h, v0); + h = fnv(h, v1); + h = fnv(h, v2); + h = fnv(h, v3); + h = fnv(h, v4); + h = fnv(h, v5); + h = fnv(h, v6); + h = fnv(h, v7); + count = count + 1; + o = o + 1; + } + if (metric == M_HASH) { + return h as f64; + } + return count as f64; + } + + // ---- tests ---- + + test "the bodies: sixteen trits that sum to zero, radius by the cube root of mass" { + assert net_charge() == 0.0; + assert charge_of(0) == 1.0; + assert charge_of(1) == 0.0; + assert charge_of(2) == 0.0 - 1.0; + assert hue_of(charge_of(0)) == 45.0; + assert hue_of(charge_of(1)) == 275.0; + assert hue_of(charge_of(2)) == 190.0; + assert @abs(radius_of(8.0) - 2.0 * R_SCALE) < 1.0e-15; + assert VERSION == 1; + assert N == 16; + } + + test "an event sets a command, never a position; DONE hands the matter to gravity" { + assert cmd_control(DONE) == 0.0; + assert cmd_gravity(DONE) == 1.0; + assert cmd_gravity(THINK) == 0.0; + assert cmd_spin(THINK) == SPIN; + assert cmd_x(TOOL_CALL) == TOOL_X; + assert cmd_y(AGENT_MSG) == USER_Y; + assert cmd_r(USER_MSG) < cmd_r(THINK); + } + + test "Verlet is exact in free fall: y = y0 - g t^2 / 2" { + const y = world(SC_FALL, 18, false, M_Y, 0, 0); + const t = 18.0 * FRAME_DT; + assert @abs(y - (1.0 - 0.5 * G_STD * t * t)) < 1.0e-12; + assert @abs(world(SC_FALL, 18, false, M_VY, 0, 0) + G_STD * t) < 1.0e-12; + } + + test "the exponential split is exact for drag: v = v0 exp(-b t / m)" { + const v = world(SC_DRAG, 120, false, M_VX, 0, 0); + assert @abs(v - @exp(0.0 - DRAG * 2.0)) < 1.0e-12; + } + + test "falling through air reaches the terminal speed m g / b" { + // 40 s is 24 time constants m / b. What is left is the splitting itself: its fixed point is + // (g / b)(1 - (b H)^2 / 24 + ...), 6.5e-8 of the speed for a 1 kg body. + const v = world(SC_TERMINAL, 2400, false, M_VY, 0, 0); + assert @abs(v / (0.0 - G_STD / DRAG) - 1.0) < 1.0e-6; + } + + test "a bounce on the lacquer floor returns e^2 of the height" { + const apex = world(SC_BOUNCE, 90, false, M_APEX, 0, 0); + const r = radius_of(1.0); + const ratio = (apex - r) / 0.5; + assert @abs(ratio - RESTITUTION * RESTITUTION) < 0.01; + } + + test "pair forces conserve momentum to rounding" { + assert @abs(world(SC_NBODY, 600, false, M_PX, 0, 0)) < 1.0e-12; + assert @abs(world(SC_NBODY, 600, false, M_PY, 0, 0)) < 1.0e-12; + } + + test "Coulomb plus elastic contact conserves energy under Verlet" { + const e0 = world(SC_NBODY, 0, false, M_ENERGY, 0, 0); + const e1 = world(SC_NBODY, 600, false, M_ENERGY, 0, 0); + assert @abs((e1 - e0) / e0) < 0.001; + assert e0 != 0.0; + } + + test "without the contact, the smooth Coulomb force keeps energy ten times better" { + const e0 = world(SC_COULOMB, 0, false, M_ENERGY, 0, 0); + const e1 = world(SC_COULOMB, 600, false, M_ENERGY, 0, 0); + assert @abs((e1 - e0) / e0) < 0.0001; + assert world(SC_COULOMB, 600, false, M_MAXSPEED, 0, 0) > 0.05; + } + + test "a head-on elastic collision gives back its energy: Hertz contact, 0.5 percent" { + // Two 1 kg bodies, 1 m/s closing, no air, no dashpot. The momentum is exact; the energy + // error is the step where the contact switches on, smaller for Hertz (force ~ d^1.5) than it + // was for a linear spring (measured 10 to 20 percent at 4 substeps). + assert @abs(world(SC_COLLIDE, 60, false, M_KE, 0, 0) / 0.25 - 1.0) < 0.005; + assert world(SC_COLLIDE, 60, false, M_VX, 0, 0) < 0.0 - 0.49; + assert @abs(world(SC_COLLIDE, 60, false, M_PX, 0, 0)) < 1.0e-12; + assert world(SC_COLLIDE, 60, false, M_DIST, 0, 0) > 2.0 * radius_of(1.0); + } + + test "a +1 and a -1 body bond where the Coulomb pull and the contact push balance" { + const sim = world(SC_DIMER, 900, false, M_DIST, 0, 0); + const root = dimer_newton(); + assert root < 2.0 * radius_of(1.0); + assert @abs(sim - root) < 1.0e-6; + } + + test "the controller is critically damped: it arrives and never overshoots" { + const x = world(SC_PD, 100, false, M_X, 0, 0); + assert @abs(x - 0.3) < 0.003; + assert world(SC_PD, 300, false, M_MAXX, 0, 0) < 0.3 + 0.000001; + } + + test "thinking spins the matter counterclockwise; reduced motion does not" { + assert world(SC_SPIN, 120, false, M_LZ, 0, 0) > 0.01; + assert @abs(world(SC_SPIN, 120, true, M_LZ, 0, 0)) < 0.001; + } + + test "an error is a radial burst: faster than anything the conversation does" { + // The kick is KICK = 1.2 m/s; the controller damps it by exp(-KD FRAME_DT) = 0.875 before the + // frame ends, so the frame sees more than 1.0 m/s. The calm log never passes 0.7 m/s. + assert world(SC_ERRSHOW, 262, false, M_MAXSPEED, 0, 0) > 1.0; + assert world(SC_SHOW, 262, false, M_MAXSPEED, 0, 0) < 0.7; + } + + test "after DONE the matter lies on the floor and is still" { + // 15 s after DONE: every body within 2 mm of the floor (none stacked), the kinetic energy + // below 0.1 mJ. + assert world(SC_SHOW, 330 + 900, false, M_MAXREST, 0, 0) < 0.002; + assert world(SC_SHOW, 330 + 900, false, M_KE, 0, 0) < 0.0001; + } + + test "Fresnel: lacquer reflects 4 percent head-on and everything at grazing" { + assert @abs(fresnel_r0(N_LACQUER) - 0.04) < 1.0e-15; + assert schlick(1.0, N_LACQUER) == fresnel_r0(N_LACQUER); + assert schlick(0.0, N_LACQUER) == 1.0; + assert @abs(schlick(0.5, N_LACQUER) - (0.04 + 0.96 / 32.0)) < 1.0e-15; + assert mirror_alpha(0.05, 1.0) > mirror_alpha(1.0, 1.0); + assert mirror_alpha(1.0, 1.0) > 0.04; + } + + test "the display list: 50 ops in budget, the reflection mirrors the glow in the floor" { + assert world(SC_SHOW, 200, false, M_COUNT, 0, 0) == (OPS as f64); + assert OPS <= OP_BUDGET; + assert OPS == 2 + 3 * 16; + const y = world(SC_SHOW, 200, false, M_OP, 4, 2); + assert world(SC_SHOW, 200, false, M_OP, 4, 0) == (OP_GLOW as f64); + assert world(SC_SHOW, 200, false, M_OP, 2, 0) == (OP_REFLECT as f64); + assert world(SC_SHOW, 200, false, M_OP, 2, 2) == 0.0 - y; + assert world(SC_SHOW, 200, false, M_OP, 2, 3) > world(SC_SHOW, 200, false, M_OP, 4, 3); + } + + test "reduced motion draws no streaks" { + assert world(SC_ERRSHOW, 262, true, M_OP, 3, 5) == 0.0; + assert world(SC_ERRSHOW, 262, true, M_OP, 3 + 3 * 7, 5) == 0.0; + assert world(SC_ERRSHOW, 262, false, M_OP, 3, 5) > 0.0; + } + + test "light is conserved when a glow is rasterised" { + assert @abs(kernel_sum(1.5) - 1.0) < 1.0e-6; + assert @abs(kernel_sum(3.0) - 1.0) < 1.0e-6; + } + + test "the background is true black and a body is bright" { + assert world(SC_SHOW, 60, false, M_PIXEL, 0, 0) < 1.0; + assert world(SC_SHOW, 60, false, M_PIXEL, 191 * 108 + 107, 0) < 1.0; + const col = @as(i64, @intFromFloat(@floor((AGENT_X + RING + W) / PIX))); + const row = @as(i64, @intFromFloat(@floor((Y_TOP - AGENT_Y) / PIX))); + assert world(SC_SHOW, 0, false, M_PIXEL, row * COLS + col, 0) > 100.0; + } + + test "replaying the same log gives the same picture; a different log does not" { + const h1 = world(SC_SHOW, 400, false, M_HASH, 0, 0); + const h2 = world(SC_SHOW, 400, false, M_HASH, 0, 0); + assert h1 == h2; + assert h1 != world(SC_ERRSHOW, 400, false, M_HASH, 0, 0); + } + + invariant "the controller gain is critical" { + assert KD * KD == 4.0 * KP; + } +}