nx_softdyn.nx source

1// nx_softdyn.nx -- LIB: SOVEREIGN SECONDARY DYNAMICS ("jiggle physics"), the general solver every animated 2// thing inherits. A driven ANCHOR (a bone from the rig, a weapon mount, a hair root) pulls a mass point through 3// a spring-damper; the mass LAGS under acceleration, OVERSHOOTS when the anchor stops, oscillates, and damps to 4// rest. The SAME solver serves: soft-tissue inertia on a moving body, firearm RECOIL (impulse), cloth and hair 5// sway, vehicle suspension. 100% integer + deterministic -> the same motion replays bit-identically. 6// 7// State stride 6 per point: [px,py,pz, vx,vy,vz] all Q8 (model-units * 256). 8// Semi-implicit Euler: a = K*(anchor-p)/1024 - C*v/1024 ; v += a ; p += v 9// K = stiffness per-1024 (how hard tissue is bound to bone), C = damping per-1024 (how fast wobble dies). 10// Underdamped (C < 2*sqrt(K)) => visible jiggle. Critically/over-damped => rigid, no wobble. 11// A hard displacement CLAMP makes divergence impossible BY CONSTRUCTION (tissue can never fly off the body). 12// license_tier: ORIGINAL 13import "nx_syscalls.nx" 14import "nx_vecmath.nx" 15 16const SD_Q8: i64 = 256 17const SD_G: i64 = 1024 18const SD_STRIDE: i64 = 6 19// stability envelope (enforced in sd_step): C must stay under one full per-tick velocity reversal, 20// K under omega*dt < 2. Values outside diverge -- so they are clamped, never honoured. 21const SD_CAPC: i64 = 1000 22const SD_CAPK: i64 = 3600 23 24func sd_abs(v: i64) -> i64 { if v < 0 { return 0 - v } return v } 25func sd_isqrt(v: i64) -> i64 { return vm_isqrt(v) } 26 27// STABILITY ENVELOPE as one callable each, so the per-axis path cannot drift from the scalar one. 28// Extracted from sd_step's own inline clamps (nx_softdyn_gate T7, 2026-07-23) -- same arithmetic, 29// one definition. A config outside the envelope stays IMPOSSIBLE TO REQUEST rather than trusted. 30func sd_clampk(K: i64) -> i64 { var k: i64 = K; if k > SD_CAPK { k = SD_CAPK } if k < 0 { k = 0 } return k } 31func sd_clampc(C: i64) -> i64 { var c: i64 = C; if c > SD_CAPC { c = SD_CAPC } if c < 0 { c = 0 } return c } 32 33// allocate state for n points 34func sd_alloc(n: i64) -> *i64 { return sys_mmap(n*SD_STRIDE*8) as *i64 } 35 36// seat point i exactly at its anchor, at rest (no startup transient) 37func sd_seat(st: *i64, i: i64, ax: i64, ay: i64, az: i64) -> i64 { 38 let b: i64 = i*SD_STRIDE 39 st[b] = ax*SD_Q8; st[b+1] = ay*SD_Q8; st[b+2] = az*SD_Q8 40 st[b+3] = 0; st[b+4] = 0; st[b+5] = 0 41 return 0 42} 43 44// inject an impulse (velocity kick) -- RECOIL, a hit, a footfall. Units: model-units/tick * 256. 45func sd_impulse(st: *i64, i: i64, ix: i64, iy: i64, iz: i64) -> i64 { 46 let b: i64 = i*SD_STRIDE 47 st[b+3] = st[b+3] + ix; st[b+4] = st[b+4] + iy; st[b+5] = st[b+5] + iz 48 return 0 49} 50 51// ★TRUNCATION DEAD-ZONE SNAP (found 2026-08-03 by nx_softtissue_conserve_gate T1 + nx_cg_restprobe). 52// Both force terms are integer divisions by SD_G that TRUNCATE toward zero, so the damping term is 53// exactly 0 whenever |v|*C < SD_G and the spring term is exactly 0 whenever |e|*K < SD_G. Inside BOTH 54// dead zones the solver applies NO force at all and the point coasts ballistically forever: measured 55// as a perfect period-30 limit cycle at +-52 q8 (+-3.2 mm at 2.0 Hz) still running undecayed after 56// 200,000 ticks with a stationary anchor. A damped system that cannot reach its own rest state is a 57// perpetual-motion bug, not a rounding detail -- the tissue never stops moving on a standing character. 58// The remedy acts ONLY where the model already computes zero force, so it cannot alter any trajectory 59// that carries real force: it replaces "coast forever" with the one physically correct force-free 60// outcome, rest. Thresholds are DERIVED from K, C and SD_G (no magic numbers). Guarded on kk>0 and 61// cc>0 because with no spring there is no rest pose to snap to, and with no damping rest is genuinely 62// NOT an attractor -- an undamped oscillator must keep oscillating. 63func sd_deadsnap(st: *i64, b: i64, axis: i64, q: i64, kk: i64, cc: i64) -> i64 { 64 if kk <= 0 { return 0 } 65 if cc <= 0 { return 0 } 66 let e: i64 = sd_abs(st[b+axis] - q) 67 let v: i64 = sd_abs(st[b+3+axis]) 68 if e*kk < SD_G { if v*cc < SD_G { st[b+axis] = q; st[b+3+axis] = 0 } } 69 return 0 70} 71 72// one integration tick against anchor (ax,ay,az) in MODEL UNITS. maxd = displacement clamp (model units). 73func sd_step(st: *i64, i: i64, ax: i64, ay: i64, az: i64, K: i64, C: i64, maxd: i64) -> i64 { 74 return sd_step_ax(st, i, ax, ay, az, K, K, K, K, C, C, C, C, maxd) 75} 76// ANISOTROPIC + SIGN-ASYMMETRIC step -- the ONE integrator. sd_step above is exactly this with 77// equal axes and a symmetric vertical pair, so every existing caller stays bit-identical BY 78// CONSTRUCTION rather than by promise. Real tissue is neither isotropic nor sign-symmetric, and 79// a single scalar K/C cannot express either. Both facts are CITED ROWS, never tuned here: 80// knowledge/gamefeel_oracle.conf tissue_aniso_ap_permil / tissue_aniso_ml_permil 81// (mills2025-ejss, CC-BY): vertical (SI) carries the largest excursion; anterior-posterior 82// and medio-lateral run ~610 / ~600 per-mille of it. One K for all three axes asserts 1000. 83// knowledge/gamefeel_oracle.conf tissue_k_up_n_per_m 70..78 vs tissue_k_down_n_per_m 620..700 84// (cai2018-j-biomech-67-137, model RMSE <= 2.6pct vs real running): a piecewise mass-spring 85// -damper roughly 9x STIFFER BELOW static equilibrium than above it. The conf also states 86// the implementation: ONE SIGN TEST ON DISPLACEMENT SELECTS THE BRANCH, exact in integers. 87func sd_step_ax(st: *i64, i: i64, ax: i64, ay: i64, az: i64, 88 kx: i64, kyu: i64, kyd: i64, kz: i64, 89 cx: i64, cyu: i64, cyd: i64, cz: i64, maxd: i64) -> i64 { 90 return sd_step_axq(st, i, ax*SD_Q8, ay*SD_Q8, az*SD_Q8, 91 kx, kyu, kyd, kz, cx, cyu, cyd, cz, maxd) 92} 93// SUB-UNIT ANCHOR ENTRY POINT. Identical arithmetic, but the anchor arrives ALREADY IN Q8 so a 94// caller that derives an anchor from another point's CURRENT position -- a HAIR CHAIN, where each 95// segment hangs off the one above it -- does not lose that position to model-unit truncation at 96// every hop. A chain built on the model-unit entry point quantises its own coupling away: the 97// parent's sub-unit motion is floored before the child ever sees it, so the tip of a long strand 98// reads as still while the root is visibly moving. sd_step_ax above is exactly this with the 99// multiply applied, so every existing caller stays ARITHMETICALLY identical BY CONSTRUCTION: the 100// old body's first act was `let qx = ax*SD_Q8`, and that multiply is now the argument expression. 101// THE BINARY IS NOT BYTE-IDENTICAL AND SAYING SO WOULD BE A FALSE CLAIM -- extracting a function 102// adds a symbol and a delegating frame. MEASURED across this exact edit: nx_softbind_gate went 103// 35,485 -> 35,915 B (+430), sha a5d872b0 -> 0adef123. A const-rename refactor rebuilds byte- 104// identical; a function extraction cannot, so the proof owed here is BEHAVIOURAL, not byte, and 105// it is nx_softbind_gate staying GREEN across the change. 106func sd_step_axq(st: *i64, i: i64, qx: i64, qy: i64, qz: i64, 107 kx: i64, kyu: i64, kyd: i64, kz: i64, 108 cx: i64, cyu: i64, cyd: i64, cz: i64, maxd: i64) -> i64 { 109 let b: i64 = i*SD_STRIDE 110 // ★STABILITY ENVELOPE, ENFORCED (found by nx_softdyn_gate T7, 2026-07-23): with semi-implicit Euler the 111 // damping term is a PER-TICK velocity multiplier -- C >= SD_G inverts and amplifies velocity (divergence), 112 // and K > SD_CAPK breaks omega*dt < 2. Clamping here makes an unstable config IMPOSSIBLE TO REQUEST rather 113 // than trusting every caller to know the envelope. 114 // THE SIGN TEST: the vertical branch is selected by where the tissue currently sits relative 115 // to its anchor. Below it (dy < 0) the suspensory ligament and skin carry tension -- the stiff 116 // row. This is the whole of the cai2018 piecewise model and it is exact in integers. 117 var ky: i64 = kyu 118 var cy: i64 = cyu 119 if st[b+1] - qy < 0 { ky = kyd; cy = cyd } 120 let kkx: i64 = sd_clampk(kx); let kky: i64 = sd_clampk(ky); let kkz: i64 = sd_clampk(kz) 121 let ccx: i64 = sd_clampc(cx); let ccy: i64 = sd_clampc(cy); let ccz: i64 = sd_clampc(cz) 122 // spring toward anchor + viscous damping, per axis 123 let accx: i64 = ((qx - st[b])*kkx)/SD_G - (st[b+3]*ccx)/SD_G 124 let accy: i64 = ((qy - st[b+1])*kky)/SD_G - (st[b+4]*ccy)/SD_G 125 let accz: i64 = ((qz - st[b+2])*kkz)/SD_G - (st[b+5]*ccz)/SD_G 126 st[b+3] = st[b+3] + accx; st[b+4] = st[b+4] + accy; st[b+5] = st[b+5] + accz 127 st[b] = st[b] + st[b+3]; st[b+1] = st[b+1] + st[b+4]; st[b+2] = st[b+2] + st[b+5] 128 // HARD CLAMP: tissue can never separate further than maxd from bone -- divergence impossible by construction 129 var dx: i64 = st[b] - qx; var dy: i64 = st[b+1] - qy; var dz: i64 = st[b+2] - qz 130 let lim: i64 = maxd*SD_Q8 131 let d2: i64 = dx*dx + dy*dy + dz*dz 132 if d2 > lim*lim { 133 var d: i64 = sd_isqrt(d2) 134 if d < 1 { d = 1 } 135 st[b] = qx + dx*lim/d; st[b+1] = qy + dy*lim/d; st[b+2] = qz + dz*lim/d 136 // kill the outward velocity component so it rests on the limit instead of grinding 137 st[b+3] = st[b+3]/2; st[b+4] = st[b+4]/2; st[b+5] = st[b+5]/2 138 } 139 sd_deadsnap(st, b, 0, qx, kkx, ccx) 140 sd_deadsnap(st, b, 1, qy, kky, ccy) 141 sd_deadsnap(st, b, 2, qz, kkz, ccz) 142 return 0 143} 144 145// current displacement of point i from its anchor, in MODEL UNITS (the jiggle amplitude) 146func sd_disp(st: *i64, i: i64, ax: i64, ay: i64, az: i64) -> i64 { 147 let b: i64 = i*SD_STRIDE 148 let dx: i64 = st[b] - ax*SD_Q8; let dy: i64 = st[b+1] - ay*SD_Q8; let dz: i64 = st[b+2] - az*SD_Q8 149 return sd_isqrt(dx*dx + dy*dy + dz*dz)/SD_Q8 150} 151// signed displacement along one axis (for overshoot / sign-change detection), model units 152func sd_disp_y(st: *i64, i: i64, ay: i64) -> i64 { return (st[i*SD_STRIDE+1] - ay*SD_Q8)/SD_Q8 } 153// ★SUB-UNIT accessors: jiggle amplitude is often a FRACTION of a model unit -- measuring in truncated 154// integer units hides real oscillation (that bug made gate T2 read "no overshoot" when it was oscillating). 155func sd_disp_y_q8(st: *i64, i: i64, ay: i64) -> i64 { return st[i*SD_STRIDE+1] - ay*SD_Q8 } 156func sd_disp_q8(st: *i64, i: i64, ax: i64, ay: i64, az: i64) -> i64 { 157 let b: i64 = i*SD_STRIDE 158 let dx: i64 = st[b] - ax*SD_Q8; let dy: i64 = st[b+1] - ay*SD_Q8; let dz: i64 = st[b+2] - az*SD_Q8 159 return sd_isqrt(dx*dx + dy*dy + dz*dz) 160} 161// kinetic energy proxy (sum v^2) -- must DECAY for a damped system 162func sd_energy(st: *i64, i: i64) -> i64 { 163 let b: i64 = i*SD_STRIDE 164 return (st[b+3]*st[b+3] + st[b+4]*st[b+4] + st[b+5]*st[b+5])/SD_Q8 165}