From 98e81f8ff5068fc757217399798bd5aeac7d6f7b Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 28 Jul 2026 00:33:43 +0000 Subject: [PATCH] Add realistic heavy-symmetric-top dynamics on the nutation screen MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Screen 1 models a top with the idealized Ω = τ/(Iω) relation, which assumes the axis is handed exactly the steady-precession rate and holds its tilt forever. Screen 2 now integrates the same apparatus from its full Lagrangian instead, so the tilt is a dynamical variable and a released top behaves like a real one: the axis dips, picks up precession, rises again, and its tip traces cusps, loops, or smooth waves depending only on how it was released. Physics (src/common/rigid-body/HeavySymmetricTopPhysics.ts): - Euler-angle equations of motion for a heavy symmetric top, integrated with RK4 at a 0.5 ms internal substep - Turning points from the effective potential, by bisecting the cubic in cos θ - Steady-precession roots and the critical spin below which no steady solution exists, plus the sleeping-top stability criterion - Optional viscous friction (drag on the center of mass, spin friction at the pivot) that damps the nutation and spins the top down until it falls - Inelastic mechanical stop where the axle rests against its mount Screen: - Perspective view of the wheel with the tip path and the two dashed circles bounding the nutation band, drawn via an exact linear projection of the Euler frame (src/common/view/TopProjection.ts) - θ(t) graph with the turning points as reference lines - Controls for spin, release tilt, release mode, and friction, with readouts for the band, both frequencies, and the critical spin - Slow-motion option, since nutation and precession differ by roughly 5× in rate - Live accessible summary and localized strings in en/fr/es Apparatus constants are a demonstration gyroscope wheel: ω_nut · Ω_slow = M g l / I₁ is fixed by the apparatus alone, so the wheel is deliberately large and short-armed to bring both timescales on screen at once. Tests assert the integrator against analytic results — energy and both angular momenta conserved to 6+ decimals over 10 s, exact steady precession at the fixed point, motion confined to and reaching the predicted turning points, and the fast-top nutation frequency I₃ω₃/I₁. Also moves PlayAreaPanel to src/common/view so both screens can use it. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01H6aawcC32DhbbRsgmeugbk --- CLAUDE.md | 50 +- src/RigidBodyPrecessionColors.ts | 24 + src/RigidBodyPrecessionConstants.ts | 70 +++ .../rigid-body/HeavySymmetricTopPhysics.ts | 510 ++++++++++++++++++ src/common/rigid-body/TopTipTrace.ts | 99 ++++ .../view/PlayAreaPanel.ts | 2 +- src/common/view/TopProjection.ts | 163 ++++++ src/i18n/strings_en.json | 30 +- src/i18n/strings_es.json | 30 +- src/i18n/strings_fr.json | 30 +- src/nutation-screen/model/NutationModel.ts | 259 ++++++++- .../view/NutationAngleGraphNode.ts | 253 +++++++++ .../view/NutationControlPanel.ts | 254 +++++++++ .../view/NutationKeyboardHelpContent.ts | 12 +- .../view/NutationScreenSummaryContent.ts | 26 +- .../view/NutationScreenView.ts | 95 +++- src/nutation-screen/view/TopSceneNode.ts | 215 ++++++++ .../view/SteadyPrecessionScreenView.ts | 2 +- tests/HeavySymmetricTopPhysics.test.ts | 389 +++++++++++++ tests/NutationModel.test.ts | 142 +++++ 20 files changed, 2612 insertions(+), 43 deletions(-) create mode 100644 src/common/rigid-body/HeavySymmetricTopPhysics.ts create mode 100644 src/common/rigid-body/TopTipTrace.ts rename src/{steady-precession-screen => common}/view/PlayAreaPanel.ts (94%) create mode 100644 src/common/view/TopProjection.ts create mode 100644 src/nutation-screen/view/NutationAngleGraphNode.ts create mode 100644 src/nutation-screen/view/NutationControlPanel.ts create mode 100644 src/nutation-screen/view/TopSceneNode.ts create mode 100644 tests/HeavySymmetricTopPhysics.test.ts create mode 100644 tests/NutationModel.test.ts diff --git a/CLAUDE.md b/CLAUDE.md index a0ecec8..758ecca 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -13,20 +13,56 @@ multi-screen sims, see [`doc/multi-screen.md`](doc/multi-screen.md). | File | Purpose | |---|---| | `src/RigidBodyPrecessionColors.ts` | All `ProfileColorProperty` instances | -| `src/SimConstants.ts` | Named numeric constants (layout px, physics SI units) | +| `src/RigidBodyPrecessionConstants.ts` | Named numeric constants (layout px, physics SI units) | | `src/RigidBodyPrecessionNamespace.ts` | Namespace for color property names | | `src/i18n/StringManager.ts` | Singleton localized string accessor | -| `src/precession-screen/RigidBodyPrecessionScreen.ts` | Screen wrapper | -| `src/precession-screen/model/RigidBodyPrecessionModel.ts` | Simulation state and logic | -| `src/precession-screen/view/RigidBodyPrecessionScreenView.ts` | Visual nodes, layout, `screenSummaryContent` + `pdomOrder` | -| `src/precession-screen/view/RigidBodyPrecessionScreenSummaryContent.ts` | Accessible screen summary (reference a11y pattern) | -| `src/precession-screen/view/RigidBodyPrecessionKeyboardHelpContent.ts` | Keyboard-help dialog content | +| `src/steady-precession-screen/` | Screen 1 — idealized Ω = τ/(Iω) gyroscope | +| `src/nutation-screen/` | Screen 2 — heavy symmetric top integrated from its full Lagrangian | +| `src/torque-free-screen/` | Screen 3 — Euler/Poinsot tumbling (placeholder) | +| `src/common/rigid-body/SteadyPrecessionPhysics.ts` | Screen 1's closed-form relations | +| `src/common/rigid-body/HeavySymmetricTopPhysics.ts` | Screen 2's RK4 integrator, invariants, turning points | +| `src/common/rigid-body/TopTipTrace.ts` | Ring buffer of (t, θ, φ) samples behind the tip path and θ(t) graph | +| `src/common/view/TopProjection.ts` | Euler angles → oblique 2-D projection (axis, wheel rim, tilt circles) | +| `src/common/view/PlayAreaPanel.ts` | Titled panel wrapper shared by every screen's play area | | `src/common/SimPanel.ts` | Pre-themed `Panel` wrapper (uses `RigidBodyPrecessionColors` automatically) | | `src/common/SimButtonOptions.ts` | Flat button-appearance option bundles + light-control-surface combo-box options | | `src/common/TimeModel.ts` | Composable play/pause + elapsed-time model for animated sims | | `scripts/generate-icons.ts` | PNG icons from `public/icons/icon.svg` | | `scripts/rename-sim.ts` | Automated fork/rename across all files and folders | +## Physics + +### Screen 1 — steady precession + +Closed form only: Ω = τ / (I ω), with the tilt held fixed. Deliberately idealized so the +gyroscopic relation stands alone. + +### Screen 2 — nutation (`HeavySymmetricTopPhysics.ts`) + +The full heavy-symmetric-top Lagrangian in Euler angles (θ, φ, ψ), integrated with RK4 at +a fixed 0.5 ms internal substep. θ is a dynamical variable, so releasing the top produces +real nutation. Key entry points: + +- `stepHeavyTop` — the integrator; substeps `dt`, wraps ψ, applies the tilt limits +- `createReleaseState` — the four release modes, which differ only in φ̇(0): `cusp` (0), + `loop` (−Ω_slow), `smooth` (½Ω_slow), `steady` (Ω_slow) +- `nutationTurningPoints` — bisects the turning-point cubic u̇² = f(cos θ) for the band + the axis is confined to; drawn as the two dashed circles and the graph's reference lines +- `steadyPrecessionRates` / `criticalSpinRate` — roots of I₁cos θ Ω² − I₃ω₃Ω + Mgl = 0, + and the spin below which no steady precession exists +- `totalEnergy`, `verticalAngularMomentum`, `spinAngularMomentum` — invariants, exact + without friction; the test suite asserts they hold to 6+ decimal places over 10 s + +Friction is a phenomenological viscous model (`tipDrag` on the center of mass, `spinDrag` +on the spin), off by default. `maxTilt` is an inelastic mechanical stop where the axle +rests against its mount — the nutation screen sets it to 90°. + +Apparatus constants (`NUTATION_*` in `RigidBodyPrecessionConstants.ts`) are a +demonstration gyroscope wheel. Note that ω_nut · Ω_slow = M g l / I₁ is fixed by the +apparatus alone: a hand-sized top nutates at several hertz no matter how it is spun, so +the wheel is deliberately large and short-armed to bring both timescales on screen at +once (≈1.6 Hz nutation, ≈3 s per precession revolution), with a slow-motion option. + ## Common components ### SimPanel @@ -110,6 +146,8 @@ Fleet-standard Vitest layout (keep when forking): | `vitest.config.ts` | `happy-dom` environment; `setupFiles: ["./tests/setup.ts"]`; `execArgv: ["--expose-gc"]` | | `tests/setup.ts` | Canvas / AudioContext mocks + `init({ name: "…" })` before SceneryStack imports | | `tests/TimeModel.test.ts` | Sample model unit tests — replace with real physics tests | +| `tests/HeavySymmetricTopPhysics.test.ts` | Screen 2 integrator vs. analytic results (invariants, turning points, fast-top limits) | +| `tests/NutationModel.test.ts` | Screen 2 model wiring: release, re-release, trace bounds, friction | | `tests/memory-leak.test.ts` | WeakRef + `forceGC` dispose regression (fleet pattern) | | `tests/fuzz/fuzz.spec.ts` | Optional Playwright fuzz smoke via joist `?fuzz` | | `playwright.config.ts` | Chromium project + Vite webServer for fuzz | diff --git a/src/RigidBodyPrecessionColors.ts b/src/RigidBodyPrecessionColors.ts index 73f167f..58168c3 100644 --- a/src/RigidBodyPrecessionColors.ts +++ b/src/RigidBodyPrecessionColors.ts @@ -153,6 +153,30 @@ const RigidBodyPrecessionColors = { projector: new Color(153, 153, 153, 0.35), }), + /** Path traced by the axle tip on the nutation screen. */ + tipTraceColorProperty: new ProfileColorProperty(RigidBodyPrecessionNamespace, "tipTrace", { + default: "#ffb74d", + projector: "#e65100", + }), + + /** Turning-point circles bounding the nutation band. */ + nutationBandColorProperty: new ProfileColorProperty(RigidBodyPrecessionNamespace, "nutationBand", { + default: "#ce93d8", + projector: "#6a1b9a", + }), + + /** Face of the gyroscope wheel on the nutation screen (translucent so the trace shows through). */ + wheelFillColorProperty: new ProfileColorProperty(RigidBodyPrecessionNamespace, "wheelFill", { + default: new Color(79, 195, 247, 0.55), + projector: new Color(2, 119, 189, 0.4), + }), + + /** Readout color for a state that is unstable or out of range (spin below critical). */ + warningColorProperty: new ProfileColorProperty(RigidBodyPrecessionNamespace, "warning", { + default: "#ff8a65", + projector: "#bf360c", + }), + /** Top-view inset card fill (graph-background hue @ 70%). */ sceneInsetCardColorProperty: new ProfileColorProperty(RigidBodyPrecessionNamespace, "sceneInsetCard", { default: new Color(15, 26, 46, 0.7), diff --git a/src/RigidBodyPrecessionConstants.ts b/src/RigidBodyPrecessionConstants.ts index ffabe1c..f5f456f 100644 --- a/src/RigidBodyPrecessionConstants.ts +++ b/src/RigidBodyPrecessionConstants.ts @@ -20,6 +20,12 @@ export const STEADY_PRECESSION_PANEL_WIDTH = 300; /** Height of the precession-angle graph. */ export const PRECESSION_GRAPH_HEIGHT = 160; +/** Width of the control panel on the nutation screen. */ +export const NUTATION_PANEL_WIDTH = 300; + +/** Height of the nutation-angle graph. */ +export const NUTATION_GRAPH_HEIGHT = 130; + // ── Physics defaults (SI units) ─────────────────────────────────────────────── /** Gravitational acceleration (m/s²). */ @@ -52,12 +58,76 @@ export const DEFAULT_ARM_MASS_KG = 0.2; /** Default pivot-to-mass distance (m). */ export const DEFAULT_PIVOT_TO_MASS_DISTANCE_M = 0.35; +// ── Nutation screen: heavy symmetric top (SI units) ─────────────────────────── +// +// A pivoted gyroscope wheel, sized so both timescales are watchable. The product +// ω_nut · Ω_slow = M g l / I₁ is fixed by the apparatus alone, so a hand-sized top +// (which nutates at several hertz) cannot show both motions at once; this wheel is +// the classroom demonstration gyroscope, heavy and slow enough that nutation runs +// near 1.6 Hz while one precession revolution takes about 3 s. + +/** Mass of the gyroscope wheel (kg). */ +export const NUTATION_WHEEL_MASS_KG = 3.0; + +/** Rim radius of the gyroscope wheel (m). */ +export const NUTATION_WHEEL_RADIUS_M = 0.33; + +/** Distance from pivot to the wheel's center of mass along the axle (m). */ +export const NUTATION_COM_DISTANCE_M = 0.14; + +/** I₃ — spin-axis moment of inertia of a thin rim, M R² (kg·m²). */ +export const NUTATION_SPIN_INERTIA_KG_M2 = NUTATION_WHEEL_MASS_KG * NUTATION_WHEEL_RADIUS_M ** 2; + +/** I₁ — transverse moment about the pivot: ½ M R² for the rim plus the parallel-axis M l² (kg·m²). */ +export const NUTATION_TRANSVERSE_INERTIA_KG_M2 = + 0.5 * NUTATION_WHEEL_MASS_KG * NUTATION_WHEEL_RADIUS_M ** 2 + NUTATION_WHEEL_MASS_KG * NUTATION_COM_DISTANCE_M ** 2; + +/** Default spin about the symmetry axis (rad/s) — about 1.4× the critical spin at 45°. */ +export const DEFAULT_NUTATION_SPIN_RAD_S = 7; + +/** Default release tilt from the vertical (rad) — 45°. */ +export const DEFAULT_NUTATION_TILT_RAD = Math.PI / 4; + +/** Viscous drag on the center of mass when friction is enabled (N·m·s). */ +export const NUTATION_TIP_DRAG_N_M_S = 0.055; + +/** Spin friction at the pivot when friction is enabled (N·m·s). */ +export const NUTATION_SPIN_DRAG_N_M_S = 0.008; + +/** Tip-trace samples retained (8 s at 60 Hz), which also backs the θ(t) graph. */ +export const NUTATION_TRACE_CAPACITY = 480; + +/** + * Samples actually drawn as the tip path (about 4 s, a little over one precession + * revolution). Drawing the whole buffer overlays several revolutions and turns the + * cusps into a mesh. + */ +export const NUTATION_TRACE_DRAW_SAMPLES = 260; + +/** Sample interval for the tip trace and θ(t) graph (s). */ +export const NUTATION_SAMPLE_INTERVAL_S = 1 / 60; + +/** Visible time window of the θ(t) graph (s). */ +export const NUTATION_GRAPH_WINDOW_S = 8; + // ── Ranges ──────────────────────────────────────────────────────────────────── export const SPIN_RATE_RANGE = { min: 5 * 2 * Math.PI, max: 80 * 2 * Math.PI }; export const ARM_MASS_RANGE = { min: 0.05, max: 0.5 }; export const PIVOT_DISTANCE_RANGE = { min: 0.15, max: 0.5 }; +/** Spin range for the nutation screen (rad/s); spans the critical spin at 45°. */ +export const NUTATION_SPIN_RANGE = { min: 2, max: 20 }; + +/** + * Tilt at which the axle bottoms out on its mount (rad) — horizontal. A top that loses + * its spin falls only this far; past it the axle would swing through the support. + */ +export const NUTATION_MAX_TILT_RAD = Math.PI / 2; + +/** Release-tilt range for the nutation screen (rad) — 15° to 80°. */ +export const NUTATION_TILT_RANGE = { min: (15 * Math.PI) / 180, max: (80 * Math.PI) / 180 }; + RigidBodyPrecessionNamespace.register("RigidBodyPrecessionConstants", { SCREEN_VIEW_MARGIN, PANEL_CORNER_RADIUS, diff --git a/src/common/rigid-body/HeavySymmetricTopPhysics.ts b/src/common/rigid-body/HeavySymmetricTopPhysics.ts new file mode 100644 index 0000000..d86fb31 --- /dev/null +++ b/src/common/rigid-body/HeavySymmetricTopPhysics.ts @@ -0,0 +1,510 @@ +/** + * HeavySymmetricTopPhysics.ts + * + * Full Lagrangian dynamics of a heavy symmetric top pivoted at a fixed point — + * the realistic counterpart to Screen 1's idealized Ω = τ / (I ω) relation. + * + * ── Coordinates ─────────────────────────────────────────────────────────────── + * Euler angles in the z-x-z convention, with the pivot at the origin and gravity + * along −z: + * + * θ nutation angle — tilt of the symmetry axis from the upward vertical + * φ precession angle — azimuth of the symmetry axis about the vertical + * ψ spin angle — rotation of the body about its own symmetry axis + * + * The Lagrangian for a body with transverse moment I₁ (about the pivot) and + * spin moment I₃, mass M, and pivot-to-center-of-mass distance l is + * + * L = ½ I₁ (θ̇² + φ̇² sin²θ) + ½ I₃ (ψ̇ + φ̇ cos θ)² − M g l cos θ + * + * with ω₃ = ψ̇ + φ̇ cos θ the body spin about the symmetry axis. The + * Euler-Lagrange equations, with generalized damping forces Q, give + * + * I₁ θ̈ = I₁ φ̇² sin θ cos θ − I₃ ω₃ φ̇ sin θ + M g l sin θ + Q_θ + * ṗ_φ = Q_φ, p_φ = I₁ φ̇ sin²θ + I₃ ω₃ cos θ + * ṗ_ψ = Q_ψ, p_ψ = I₃ ω₃ + * + * ── Why this is "more realistic" ────────────────────────────────────────────── + * Screen 1 assumes the axis is handed exactly the steady-precession rate. A real + * top is released from rest, so θ is a dynamical variable: the axis dips, picks + * up precession, rises again, and the tip traces cusps, loops, or smooth waves + * depending on how it was released. Below a critical spin no steady precession + * exists at all and the top flops over. Optional viscous friction spins the top + * down, damps the nutation away (a real top's wobble dies out within seconds), + * and eventually drops the axis. + * + * ── Damping model ───────────────────────────────────────────────────────────── + * Phenomenological viscous drag on the center of mass (air resistance plus pivot + * friction), lumped into one coefficient c = b l². The center of mass moves with + * velocity components l θ̇ and l φ̇ sin θ, so a drag force −b v contributes + * + * Q_θ = −c θ̇ Q_φ = −c φ̇ sin²θ + * + * plus a separate spin friction Q_ψ = −c_s ω₃ at the pivot contact. With both + * coefficients zero the integrator conserves E, p_φ, and p_ψ. + */ + +export type HeavyTopParameters = { + /** I₁ — transverse moment of inertia about the pivot (kg·m²). */ + readonly transverseInertia: number; + /** I₃ — moment of inertia about the symmetry axis (kg·m²). */ + readonly spinInertia: number; + /** Total mass (kg). */ + readonly mass: number; + /** Gravitational acceleration (m/s²). */ + readonly gravity: number; + /** l — distance from pivot to center of mass along the symmetry axis (m). */ + readonly comDistance: number; + /** Viscous drag coefficient c on the center of mass (N·m·s); 0 disables it. */ + readonly tipDrag: number; + /** Spin friction coefficient c_s at the pivot (N·m·s); 0 disables it. */ + readonly spinDrag: number; + /** + * Largest tilt the axle can reach before it comes to rest against its mount (rad). + * Contact is inelastic: θ stops there while φ and ψ keep evolving. Defaults to the + * numerical limit just short of θ = π, i.e. no mechanical stop. + */ + readonly maxTilt?: number; +}; + +export type HeavyTopState = { + /** Nutation angle θ from the upward vertical (rad). */ + readonly theta: number; + /** θ̇ (rad/s). */ + readonly thetaDot: number; + /** Precession angle φ, unwrapped so its slope is the mean precession rate (rad). */ + readonly phi: number; + /** φ̇ (rad/s). */ + readonly phiDot: number; + /** Spin angle ψ about the symmetry axis, wrapped to [0, 2π) (rad). */ + readonly psi: number; + /** ω₃ = ψ̇ + φ̇ cos θ — body spin about the symmetry axis (rad/s). */ + readonly spin: number; +}; + +/** How the top is launched. Determines φ̇(0) and hence the shape of the tip trace. */ +export type ReleaseMode = + /** Released from rest: φ̇ = 0, so the tip traces cusps at the top of each dip. */ + | "cusp" + /** Pushed backwards against the precession: the tip traces retrograde loops. */ + | "loop" + /** Given a gentle forward nudge: smooth undulations with no cusps. */ + | "smooth" + /** Handed the exact steady-precession rate: the axis holds θ with no nutation. */ + | "steady"; + +export type NutationBand = { + /** Smallest θ reached (highest the axis rises), rad. */ + readonly thetaMin: number; + /** Largest θ reached (lowest the axis dips), rad. */ + readonly thetaMax: number; +}; + +export type SteadyPrecessionRoots = { + /** Whether ω₃ is large enough for steady precession at this θ to exist. */ + readonly exists: boolean; + /** Slow (physical) precession root Ω₋ (rad/s). */ + readonly slow: number; + /** Fast precession root Ω₊ (rad/s). */ + readonly fast: number; +}; + +/** sin θ is clamped to this magnitude in denominators to keep φ̈ finite near the poles. */ +const MIN_SIN_THETA = 1e-3; + +/** θ is confined to this open interval; a real top never reaches the coordinate poles. */ +const MIN_THETA = 1e-3; +const MAX_THETA = Math.PI - 1e-3; + +/** Largest internal RK4 step (s). Nutation is the fastest mode we must resolve. */ +const MAX_INTERNAL_STEP_S = 5e-4; + +/** Safety cap so a huge dt can never lock up the frame. */ +const MAX_SUBSTEPS = 400; + +const TURNING_POINT_ITERATIONS = 60; + +const TWO_PI = 2 * Math.PI; + +/** M g l — the gravitational torque coefficient (N·m). */ +export function gravityTorqueCoefficient(parameters: HeavyTopParameters): number { + return parameters.mass * parameters.gravity * parameters.comDistance; +} + +/** p_ψ = I₃ ω₃ — angular momentum about the symmetry axis (kg·m²/s). */ +export function spinAngularMomentum(parameters: HeavyTopParameters, state: HeavyTopState): number { + return parameters.spinInertia * state.spin; +} + +/** p_φ = I₁ φ̇ sin²θ + I₃ ω₃ cos θ — angular momentum about the vertical (kg·m²/s). */ +export function verticalAngularMomentum(parameters: HeavyTopParameters, state: HeavyTopState): number { + const sinTheta = Math.sin(state.theta); + return ( + parameters.transverseInertia * state.phiDot * sinTheta * sinTheta + + parameters.spinInertia * state.spin * Math.cos(state.theta) + ); +} + +/** Total mechanical energy, kinetic plus gravitational potential, about the pivot (J). */ +export function totalEnergy(parameters: HeavyTopParameters, state: HeavyTopState): number { + const sinTheta = Math.sin(state.theta); + const rotational = + 0.5 * + parameters.transverseInertia * + (state.thetaDot * state.thetaDot + state.phiDot * state.phiDot * sinTheta * sinTheta) + + 0.5 * parameters.spinInertia * state.spin * state.spin; + return rotational + gravityTorqueCoefficient(parameters) * Math.cos(state.theta); +} + +/** + * Effective one-dimensional potential governing θ, obtained by eliminating φ̇ and ψ̇ + * through the conserved momenta: + * + * V_eff(θ) = (p_φ − p_ψ cos θ)² / (2 I₁ sin²θ) + M g l cos θ + */ +export function effectivePotential( + parameters: HeavyTopParameters, + verticalMomentum: number, + spinMomentum: number, + theta: number, +): number { + const sinTheta = Math.max(Math.abs(Math.sin(theta)), MIN_SIN_THETA); + const numerator = verticalMomentum - spinMomentum * Math.cos(theta); + return ( + (numerator * numerator) / (2 * parameters.transverseInertia * sinTheta * sinTheta) + + gravityTorqueCoefficient(parameters) * Math.cos(theta) + ); +} + +/** + * The two steady-precession rates at a fixed tilt, from the θ̈ = 0 condition + * + * I₁ cos θ Ω² − I₃ ω₃ Ω + M g l = 0 + * + * Real roots require I₃²ω₃² ≥ 4 I₁ M g l cos θ — below that critical spin the top + * cannot precess at constant tilt and must nutate. For θ > 90° (cos θ < 0) the + * roots always exist and have opposite signs. + */ +export function steadyPrecessionRates( + parameters: HeavyTopParameters, + spin: number, + theta: number, +): SteadyPrecessionRoots { + const mgl = gravityTorqueCoefficient(parameters); + const pSpin = parameters.spinInertia * spin; + const cosTheta = Math.cos(theta); + const a = parameters.transverseInertia * cosTheta; + + // Degenerate at θ = 90°: the quadratic collapses to the gyroscopic relation Ω = M g l / p_ψ. + if (Math.abs(a) < 1e-9) { + const rate = pSpin === 0 ? 0 : mgl / pSpin; + return { exists: pSpin !== 0, slow: rate, fast: rate }; + } + + const discriminant = pSpin * pSpin - 4 * a * mgl; + if (discriminant < 0) { + return { exists: false, slow: 0, fast: 0 }; + } + + const root = Math.sqrt(discriminant); + const first = (pSpin - root) / (2 * a); + const second = (pSpin + root) / (2 * a); + const slow = Math.abs(first) <= Math.abs(second) ? first : second; + const fast = Math.abs(first) <= Math.abs(second) ? second : first; + return { exists: true, slow, fast }; +} + +/** + * Minimum spin for steady precession at this tilt: ω₃,min = 2 √(I₁ M g l cos θ) / I₃. + * Returns 0 for θ ≥ 90°, where steady precession is always possible. + */ +export function criticalSpinRate(parameters: HeavyTopParameters, theta: number): number { + const cosTheta = Math.cos(theta); + if (cosTheta <= 0) { + return 0; + } + return ( + (2 * Math.sqrt(parameters.transverseInertia * gravityTorqueCoefficient(parameters) * cosTheta)) / + parameters.spinInertia + ); +} + +/** + * The slow precession rate used when launching the top. Falls back to the + * gyroscopic approximation Ω ≈ M g l / (I₃ ω₃) below the critical spin, where the + * exact root does not exist. + */ +export function slowPrecessionRate(parameters: HeavyTopParameters, spin: number, theta: number): number { + const roots = steadyPrecessionRates(parameters, spin, theta); + if (roots.exists) { + return roots.slow; + } + const pSpin = parameters.spinInertia * spin; + return pSpin === 0 ? 0 : gravityTorqueCoefficient(parameters) / pSpin; +} + +/** + * Nutation angular frequency of a fast top, ω_nut ≈ I₃ ω₃ / I₁ (rad/s). + * Exact in the limit where the spin dominates the gravitational torque. + */ +export function nutationFrequency(parameters: HeavyTopParameters, spin: number): number { + return Math.abs((parameters.spinInertia * spin) / parameters.transverseInertia); +} + +/** + * Whether a top spinning upright would "sleep" (θ = 0 is stable), which requires + * I₃²ω₃² > 4 I₁ M g l. + */ +export function isSleepingTopStable(parameters: HeavyTopParameters, spin: number): boolean { + const pSpin = parameters.spinInertia * spin; + return pSpin * pSpin > 4 * parameters.transverseInertia * gravityTorqueCoefficient(parameters); +} + +/** + * The turning-point equation in u = cos θ. With a = p_ψ/I₁, b = p_φ/I₁, + * α = 2E′/I₁ (E′ the energy less the constant spin term) and β = 2Mgl/I₁, + * + * u̇² = (1 − u²)(α − β u) − (b − a u)² ≡ f(u) + * + * f is a cubic that is non-negative exactly on the interval of u the motion + * visits, and f(±1) = −(b ∓ a)² ≤ 0. + */ +function turningPointFunction(parameters: HeavyTopParameters, state: HeavyTopState): (u: number) => number { + const inertia = parameters.transverseInertia; + const pSpin = spinAngularMomentum(parameters, state); + const pVertical = verticalAngularMomentum(parameters, state); + const reducedEnergy = totalEnergy(parameters, state) - (pSpin * pSpin) / (2 * parameters.spinInertia); + + const a = pSpin / inertia; + const b = pVertical / inertia; + const alpha = (2 * reducedEnergy) / inertia; + const beta = (2 * gravityTorqueCoefficient(parameters)) / inertia; + + return (u: number): number => { + const linear = b - a * u; + return (1 - u * u) * (alpha - beta * u) - linear * linear; + }; +} + +/** Bisect f between a point where it is non-negative and one where it is non-positive. */ +function bisectRoot(f: (u: number) => number, insideU: number, outsideU: number): number { + let inside = insideU; + let outside = outsideU; + for (let i = 0; i < TURNING_POINT_ITERATIONS; i++) { + const mid = 0.5 * (inside + outside); + if (f(mid) >= 0) { + inside = mid; + } else { + outside = mid; + } + } + return 0.5 * (inside + outside); +} + +/** + * A point strictly inside the allowed band to bisect outward from. Released tops + * start exactly at a turning point, where f(u) is zero up to rounding and can come + * out slightly negative, so probe both sides for the interior. + */ +function findInteriorPoint(f: (u: number) => number, u: number): number | null { + if (f(u) >= 0) { + return u; + } + for (const probe of [1e-9, 1e-7, 1e-5, 1e-3, 1e-2]) { + if (u - probe >= -1 && f(u - probe) >= 0) { + return u - probe; + } + if (u + probe <= 1 && f(u + probe) >= 0) { + return u + probe; + } + } + return null; +} + +/** + * The band [θ_min, θ_max] the symmetry axis oscillates between, from the current + * energy and momenta. With friction enabled the invariants drift, so this is the + * band the top would settle into if friction were switched off right now. + */ +export function nutationTurningPoints(parameters: HeavyTopParameters, state: HeavyTopState): NutationBand { + const f = turningPointFunction(parameters, state); + const u = Math.cos(state.theta); + + const interior = findInteriorPoint(f, u); + if (interior === null) { + return { thetaMin: state.theta, thetaMax: state.theta }; + } + + const upper = f(1) >= 0 ? 1 : bisectRoot(f, interior, 1); + const lower = f(-1) >= 0 ? -1 : bisectRoot(f, interior, -1); + + // The mechanical stop can cut the band short of its analytic turning point. + return { + thetaMin: Math.acos(Math.min(1, Math.max(-1, upper))), + thetaMax: Math.min(maximumTilt(parameters), Math.acos(Math.min(1, Math.max(-1, lower)))), + }; +} + +type Derivatives = { + readonly theta: number; + readonly thetaDot: number; + readonly phi: number; + readonly phiDot: number; + readonly psi: number; + readonly spin: number; +}; + +/** + * Time derivatives of the full state from the Euler-Lagrange equations above. + * Exported for testing; {@link stepHeavyTop} is the integration entry point. + */ +export function heavyTopDerivatives(parameters: HeavyTopParameters, state: HeavyTopState): Derivatives { + const { transverseInertia: inertia1, spinInertia: inertia3 } = parameters; + const mgl = gravityTorqueCoefficient(parameters); + + const sinTheta = Math.sin(state.theta); + const cosTheta = Math.cos(state.theta); + const safeSinTheta = Math.sign(sinTheta || 1) * Math.max(Math.abs(sinTheta), MIN_SIN_THETA); + const sinSquared = safeSinTheta * safeSinTheta; + + // ṗ_ψ = −c_s ω₃ → spin decays exponentially with friction, is constant without it. + const spinAcceleration = (-parameters.spinDrag * state.spin) / inertia3; + + // I₁ θ̈ = I₁ φ̇² sin θ cos θ − I₃ ω₃ φ̇ sin θ + M g l sin θ − c θ̇ + const thetaAcceleration = + (inertia1 * state.phiDot * state.phiDot * sinTheta * cosTheta - + inertia3 * state.spin * state.phiDot * sinTheta + + mgl * sinTheta - + parameters.tipDrag * state.thetaDot) / + inertia1; + + // Differentiating p_φ = I₁ φ̇ sin²θ + I₃ ω₃ cos θ and setting ṗ_φ = −c φ̇ sin²θ. + const phiAcceleration = + (-parameters.tipDrag * state.phiDot * sinSquared - + 2 * inertia1 * state.phiDot * state.thetaDot * sinTheta * cosTheta - + inertia3 * spinAcceleration * cosTheta + + inertia3 * state.spin * state.thetaDot * sinTheta) / + (inertia1 * sinSquared); + + return { + theta: state.thetaDot, + thetaDot: thetaAcceleration, + phi: state.phiDot, + phiDot: phiAcceleration, + psi: state.spin - state.phiDot * cosTheta, + spin: spinAcceleration, + }; +} + +function addScaled(state: HeavyTopState, derivatives: Derivatives, scale: number): HeavyTopState { + return { + theta: state.theta + derivatives.theta * scale, + thetaDot: state.thetaDot + derivatives.thetaDot * scale, + phi: state.phi + derivatives.phi * scale, + phiDot: state.phiDot + derivatives.phiDot * scale, + psi: state.psi + derivatives.psi * scale, + spin: state.spin + derivatives.spin * scale, + }; +} + +function rungeKutta4(parameters: HeavyTopParameters, state: HeavyTopState, dt: number): HeavyTopState { + const k1 = heavyTopDerivatives(parameters, state); + const k2 = heavyTopDerivatives(parameters, addScaled(state, k1, dt / 2)); + const k3 = heavyTopDerivatives(parameters, addScaled(state, k2, dt / 2)); + const k4 = heavyTopDerivatives(parameters, addScaled(state, k3, dt)); + + const weight = dt / 6; + return { + theta: state.theta + weight * (k1.theta + 2 * k2.theta + 2 * k3.theta + k4.theta), + thetaDot: state.thetaDot + weight * (k1.thetaDot + 2 * k2.thetaDot + 2 * k3.thetaDot + k4.thetaDot), + phi: state.phi + weight * (k1.phi + 2 * k2.phi + 2 * k3.phi + k4.phi), + phiDot: state.phiDot + weight * (k1.phiDot + 2 * k2.phiDot + 2 * k3.phiDot + k4.phiDot), + psi: state.psi + weight * (k1.psi + 2 * k2.psi + 2 * k3.psi + k4.psi), + spin: state.spin + weight * (k1.spin + 2 * k2.spin + 2 * k3.spin + k4.spin), + }; +} + +/** The tilt at which the axle rests against its mount, defaulting to no mechanical stop. */ +export function maximumTilt(parameters: HeavyTopParameters): number { + return Math.min(MAX_THETA, parameters.maxTilt ?? MAX_THETA); +} + +/** + * Enforce the tilt limits: the mechanical stop at the top of the range, and the + * coordinate singularity at θ = 0 where the Euler angles break down. + */ +function applyTiltLimits(state: HeavyTopState, maxTheta: number): HeavyTopState { + if (state.theta >= MIN_THETA && state.theta <= maxTheta) { + return state; + } + if (state.theta > maxTheta) { + // Inelastic contact: the axle stops falling but is free to rise off the stop again. + return { ...state, theta: maxTheta, thetaDot: Math.min(0, state.thetaDot) }; + } + // Reflect θ̇ so the axis turns around at the pole rather than sticking to it. + return { ...state, theta: MIN_THETA, thetaDot: -state.thetaDot }; +} + +/** + * Advance the top by dt using RK4 with internal substeps small enough to resolve + * the nutation. ψ is wrapped to [0, 2π) each step; φ is left unwrapped so its + * slope is the mean precession rate. + */ +export function stepHeavyTop(parameters: HeavyTopParameters, state: HeavyTopState, dt: number): HeavyTopState { + if (dt <= 0) { + return state; + } + + const substeps = Math.min(MAX_SUBSTEPS, Math.max(1, Math.ceil(dt / MAX_INTERNAL_STEP_S))); + const h = dt / substeps; + + const maxTheta = maximumTilt(parameters); + let current = state; + for (let i = 0; i < substeps; i++) { + current = applyTiltLimits(rungeKutta4(parameters, current, h), maxTheta); + } + + const psi = current.psi % TWO_PI; + return { ...current, psi: psi < 0 ? psi + TWO_PI : psi }; +} + +/** + * Build the launch state for a given release mode. Every mode starts at rest in θ + * (θ̇ = 0) and differs only in the initial precession rate φ̇, which is what + * distinguishes cusps from loops from smooth undulations. + */ +export function createReleaseState( + parameters: HeavyTopParameters, + spin: number, + initialTilt: number, + mode: ReleaseMode, +): HeavyTopState { + const reference = slowPrecessionRate(parameters, spin, initialTilt); + + let phiDot: number; + switch (mode) { + case "cusp": + phiDot = 0; + break; + case "loop": + phiDot = -reference; + break; + case "smooth": + // Any 0 < φ̇(0) < Ω_slow keeps the top falling at first, and φ̇ only grows as + // θ increases — so the tip undulates forwards without ever stalling into a cusp. + phiDot = 0.5 * reference; + break; + case "steady": + phiDot = reference; + break; + } + + return { + theta: Math.min(maximumTilt(parameters), Math.max(MIN_THETA, initialTilt)), + thetaDot: 0, + phi: 0, + phiDot, + psi: 0, + spin, + }; +} diff --git a/src/common/rigid-body/TopTipTrace.ts b/src/common/rigid-body/TopTipTrace.ts new file mode 100644 index 0000000..2575818 --- /dev/null +++ b/src/common/rigid-body/TopTipTrace.ts @@ -0,0 +1,99 @@ +/** + * TopTipTrace.ts + * + * Fixed-length circular buffer of (t, θ, φ) samples for the nutation screen. + * Feeds both the tip path drawn on the unit sphere and the θ(t) graph. + */ + +export type TipSample = { + readonly time: number; + /** Nutation angle from the upward vertical (rad). */ + readonly theta: number; + /** Unwrapped precession angle (rad). */ + readonly phi: number; +}; + +export class TopTipTrace { + private readonly capacity: number; + private readonly times: Float64Array; + private readonly thetas: Float64Array; + private readonly phis: Float64Array; + private count = 0; + private startIndex = 0; + + public constructor(capacity: number) { + this.capacity = capacity; + this.times = new Float64Array(capacity); + this.thetas = new Float64Array(capacity); + this.phis = new Float64Array(capacity); + } + + public push(time: number, theta: number, phi: number): void { + const index = (this.startIndex + this.count) % this.capacity; + this.times[index] = time; + this.thetas[index] = theta; + this.phis[index] = phi; + if (this.count < this.capacity) { + this.count++; + } else { + this.startIndex = (this.startIndex + 1) % this.capacity; + } + } + + public clear(): void { + this.count = 0; + this.startIndex = 0; + } + + public getSampleCount(): number { + return this.count; + } + + private indexAt(offset: number): number { + return (this.startIndex + offset) % this.capacity; + } + + public toSamples(): TipSample[] { + const samples: TipSample[] = []; + for (let i = 0; i < this.count; i++) { + const index = this.indexAt(i); + samples.push({ + time: this.times[index] ?? 0, + theta: this.thetas[index] ?? 0, + phi: this.phis[index] ?? 0, + }); + } + return samples; + } + + /** θ(t) in the {x, y} form the Bamboo LinePlot consumes, with θ in degrees. */ + public toThetaDegreePoints(): Array<{ x: number; y: number }> { + const points: Array<{ x: number; y: number }> = []; + for (let i = 0; i < this.count; i++) { + const index = this.indexAt(i); + points.push({ + x: this.times[index] ?? 0, + y: ((this.thetas[index] ?? 0) * 180) / Math.PI, + }); + } + return points; + } + + /** + * Mean precession rate over the buffered window (rad/s), i.e. the slope of the + * unwrapped φ. This is the rate a stopwatch would measure — it averages the + * nutation ripple away, unlike the instantaneous φ̇. + */ + public estimateMeanPrecessionRate(minSamples = 8): number { + if (this.count < minSamples) { + return 0; + } + const first = this.indexAt(0); + const last = this.indexAt(this.count - 1); + const dt = (this.times[last] ?? 0) - (this.times[first] ?? 0); + if (Math.abs(dt) < 1e-9) { + return 0; + } + return ((this.phis[last] ?? 0) - (this.phis[first] ?? 0)) / dt; + } +} diff --git a/src/steady-precession-screen/view/PlayAreaPanel.ts b/src/common/view/PlayAreaPanel.ts similarity index 94% rename from src/steady-precession-screen/view/PlayAreaPanel.ts rename to src/common/view/PlayAreaPanel.ts index c815f34..507c32e 100644 --- a/src/steady-precession-screen/view/PlayAreaPanel.ts +++ b/src/common/view/PlayAreaPanel.ts @@ -8,8 +8,8 @@ import type { TReadOnlyProperty } from "scenerystack/axon"; import type { Node } from "scenerystack/scenery"; import { Text, VBox } from "scenerystack/scenery"; import { PhetFont } from "scenerystack/scenery-phet"; -import { SimPanel } from "../../common/SimPanel.js"; import RigidBodyPrecessionColors from "../../RigidBodyPrecessionColors.js"; +import { SimPanel } from "../SimPanel.js"; const TITLE_FONT = new PhetFont({ size: 13, weight: "bold" }); diff --git a/src/common/view/TopProjection.ts b/src/common/view/TopProjection.ts new file mode 100644 index 0000000..c50953c --- /dev/null +++ b/src/common/view/TopProjection.ts @@ -0,0 +1,163 @@ +/** + * TopProjection.ts + * + * Projection helpers for drawing a heavy symmetric top from its Euler angles. + * + * World coordinates are right-handed with +z up and the pivot at the origin. The + * view is a linear oblique projection, matching the schematic used on Screen 1: + * + * view.x = pivot.x + s·v_x + * view.y = pivot.y − s·v_z + s·k·v_y + * + * where s is the pixels-per-meter scale and k the depth foreshortening. Because + * the map is linear, the projection of the wheel's rim circle is exactly an + * ellipse spanned by the projections of the two in-plane basis vectors — no + * approximation is needed to draw the wheel at any tilt. + */ + +import { Vector2, Vector3 } from "scenerystack/dot"; + +export type TopProjection = { + /** Pivot location in view pixels. */ + readonly pivot: Vector2; + /** Scale from meters to view pixels. */ + readonly pxPerM: number; + /** Foreshortening applied to the depth (y) axis, 0–1. */ + readonly perspective: number; +}; + +/** + * Symmetry axis n̂ for a top at nutation angle θ and precession angle φ. + * θ is measured from the upward vertical, φ is the azimuth about it. + */ +export function symmetryAxis(theta: number, phi: number): Vector3 { + const sinTheta = Math.sin(theta); + return new Vector3(sinTheta * Math.cos(phi), sinTheta * Math.sin(phi), Math.cos(theta)); +} + +/** + * Orthonormal basis of the wheel plane (both ⊥ n̂). e₁ is the horizontal line of + * nodes, so the body's spin angle ψ is measured from it; e₂ completes the frame. + */ +export function wheelFrame(theta: number, phi: number): { e1: Vector3; e2: Vector3 } { + const sinPhi = Math.sin(phi); + const cosPhi = Math.cos(phi); + const sinTheta = Math.sin(theta); + const cosTheta = Math.cos(theta); + return { + e1: new Vector3(-sinPhi, cosPhi, 0), + e2: new Vector3(-cosTheta * cosPhi, -cosTheta * sinPhi, sinTheta), + }; +} + +/** Project a world point (in meters, relative to the pivot) into view pixels. */ +export function projectPoint(projection: TopProjection, point: Vector3): Vector2 { + const s = projection.pxPerM; + return new Vector2( + projection.pivot.x + s * point.x, + projection.pivot.y - s * point.z + s * projection.perspective * point.y, + ); +} + +/** Project a world direction — the same linear map without the pivot translation. */ +export function projectDirection(projection: TopProjection, direction: Vector3): Vector2 { + const s = projection.pxPerM; + return new Vector2(s * direction.x, -s * direction.z + s * projection.perspective * direction.y); +} + +/** View position of a point at `distance` meters along the symmetry axis. */ +export function projectAxisPoint(projection: TopProjection, theta: number, phi: number, distance: number): Vector2 { + return projectPoint(projection, symmetryAxis(theta, phi).timesScalar(distance)); +} + +/** + * Polygon approximating the projected wheel rim: a circle of the given radius, + * centered `comDistance` along the axis and lying in the plane ⊥ to the axis. + */ +export function wheelRimPoints( + projection: TopProjection, + theta: number, + phi: number, + comDistance: number, + radius: number, + pointCount = 48, +): Vector2[] { + const center = projectAxisPoint(projection, theta, phi, comDistance); + const { e1, e2 } = wheelFrame(theta, phi); + const a = projectDirection(projection, e1).timesScalar(radius); + const b = projectDirection(projection, e2).timesScalar(radius); + + const points: Vector2[] = []; + for (let i = 0; i < pointCount; i++) { + const angle = (i / pointCount) * 2 * Math.PI; + points.push( + new Vector2( + center.x + a.x * Math.cos(angle) + b.x * Math.sin(angle), + center.y + a.y * Math.cos(angle) + b.y * Math.sin(angle), + ), + ); + } + return points; +} + +/** View position of the rim point at body spin angle ψ, used to draw the spin marker. */ +export function projectRimPoint( + projection: TopProjection, + theta: number, + phi: number, + comDistance: number, + radius: number, + psi: number, +): Vector2 { + const center = projectAxisPoint(projection, theta, phi, comDistance); + const { e1, e2 } = wheelFrame(theta, phi); + const a = projectDirection(projection, e1).timesScalar(radius * Math.cos(psi)); + const b = projectDirection(projection, e2).timesScalar(radius * Math.sin(psi)); + return new Vector2(center.x + a.x + b.x, center.y + a.y + b.y); +} + +/** + * The two rim points on the projected silhouette — the ones furthest from the + * projected axis line. Drawing the axle-to-rim triangle through these gives the + * wheel a solid hub without a full 3-D renderer. + */ +export function wheelSilhouette( + projection: TopProjection, + theta: number, + phi: number, + comDistance: number, + radius: number, +): [Vector2, Vector2] { + const axis = projectDirection(projection, symmetryAxis(theta, phi)); + const { e1, e2 } = wheelFrame(theta, phi); + const a = projectDirection(projection, e1).timesScalar(radius); + const b = projectDirection(projection, e2).timesScalar(radius); + + // Maximize |(a cos α + b sin α) × â| over α. + const crossA = a.x * axis.y - a.y * axis.x; + const crossB = b.x * axis.y - b.y * axis.x; + const alpha = Math.atan2(crossB, crossA); + + const center = projectAxisPoint(projection, theta, phi, comDistance); + const offsetX = a.x * Math.cos(alpha) + b.x * Math.sin(alpha); + const offsetY = a.y * Math.cos(alpha) + b.y * Math.sin(alpha); + return [new Vector2(center.x + offsetX, center.y + offsetY), new Vector2(center.x - offsetX, center.y - offsetY)]; +} + +/** + * Projected outline of the horizontal circle traced by a point at `distance` + * along the axis when the tilt is held at θ — the boundary of the nutation band. + */ +export function projectTiltCircle( + projection: TopProjection, + theta: number, + distance: number, +): { center: Vector2; radiusX: number; radiusY: number } { + const s = projection.pxPerM; + const horizontal = distance * Math.sin(theta) * s; + return { + center: new Vector2(projection.pivot.x, projection.pivot.y - distance * Math.cos(theta) * s), + radiusX: horizontal, + radiusY: horizontal * projection.perspective, + }; +} diff --git a/src/i18n/strings_en.json b/src/i18n/strings_en.json index f8d12e2..b5603fa 100644 --- a/src/i18n/strings_en.json +++ b/src/i18n/strings_en.json @@ -19,7 +19,20 @@ "formulaHint": "Ω = τ / (I ω) — faster spin → slower precession" }, "nutation": { - "placeholder": "Nutation screen — release from a tilted angle to see wobble on top of steady precession. Coming soon." + "sceneTitle": "Heavy top — path traced by the axle tip", + "graphTitle": "Tilt θ vs. time", + "spinRate": "Spin ω₃ (Hz)", + "initialTilt": "Release tilt θ₀ (°)", + "releaseMode": "Release the top", + "releaseCusp": "From rest — cusps", + "releaseLoop": "Pushed backward — loops", + "releaseSmooth": "Pushed forward — waves", + "releaseSteady": "At the steady rate — no wobble", + "friction": "Friction (top spins down)", + "releaseAgain": "Release again", + "readoutsTitle": "Measurements", + "belowCritical": "Spin is below the critical value — no steady precession exists at this tilt, so the top must nutate.", + "insight": "Real tops are released, not handed the answer: the axis dips, picks up precession, and rises again." }, "torqueFree": { "placeholder": "Torque-free tumbling — asymmetric rigid body via Euler's equations. Coming soon." @@ -41,11 +54,18 @@ }, "nutation": { "screenSummary": { - "playArea": "The play area will show nutation — wobble on top of steady precession when the top is released off-equilibrium.", - "controlArea": "Controls for initial tilt and nutation kick will appear here.", - "interactionHint": "This screen is under construction." + "playArea": "The play area shows a heavy top pivoted at its base, the path traced by its axle tip, the two circles bounding the nutation band, and a graph of the tilt angle versus time.", + "controlArea": "The control area lets you set the spin rate and release tilt, choose how the top is released, turn friction on or off, and release the top again. Readouts give the nutation band, nutation and precession frequencies, and the critical spin.", + "interactionHint": "Release the top from rest to see cusps, then try pushing it backward or forward and compare the traced paths." }, - "currentDetails": "Nutation screen placeholder." + "currentDetailsPattern": "The axle is tilted {{tilt}} degrees from vertical, nutating between {{thetaMin}} and {{thetaMax}} degrees. Mean precession is {{precessionHz}} hertz and nutation is {{nutationHz}} hertz.", + "controls": { + "spinRate": "Spin rate about the axle", + "initialTilt": "Release tilt from vertical", + "releaseMode": "How the top is released", + "friction": "Friction", + "releaseAgain": "Release the top again" + } }, "torqueFree": { "screenSummary": { diff --git a/src/i18n/strings_es.json b/src/i18n/strings_es.json index 9050ee9..1111abb 100644 --- a/src/i18n/strings_es.json +++ b/src/i18n/strings_es.json @@ -19,7 +19,20 @@ "formulaHint": "Ω = τ / (I ω) — más giro → precesión más lenta" }, "nutation": { - "placeholder": "Pantalla de nutación — próximamente." + "sceneTitle": "Trompo pesado — trayectoria de la punta del eje", + "graphTitle": "Inclinación θ vs. tiempo", + "spinRate": "Giro propio ω₃ (Hz)", + "initialTilt": "Inclinación de suelta θ₀ (°)", + "releaseMode": "Soltar el trompo", + "releaseCusp": "En reposo — cúspides", + "releaseLoop": "Empujado hacia atrás — bucles", + "releaseSmooth": "Empujado hacia adelante — ondas", + "releaseSteady": "A la tasa estable — sin bamboleo", + "friction": "Fricción (el trompo se frena)", + "releaseAgain": "Soltar de nuevo", + "readoutsTitle": "Mediciones", + "belowCritical": "El giro está por debajo del valor crítico — no existe precesión estable en esta inclinación, así que el trompo debe nutar.", + "insight": "Un trompo real se suelta, no se lanza con la tasa exacta: el eje cae, gana precesión y vuelve a subir." }, "torqueFree": { "placeholder": "Volteo sin par — próximamente." @@ -41,11 +54,18 @@ }, "nutation": { "screenSummary": { - "playArea": "Próximamente: nutación.", - "controlArea": "Controles próximamente.", - "interactionHint": "En construcción." + "playArea": "El área de juego muestra un trompo pesado sobre su pivote, la trayectoria trazada por la punta del eje, los dos círculos que delimitan la banda de nutación y una gráfica de la inclinación frente al tiempo.", + "controlArea": "El área de control permite fijar el giro y la inclinación de suelta, elegir cómo se suelta el trompo, activar la fricción y volver a soltarlo. Las lecturas dan la banda de nutación, las frecuencias y el giro crítico.", + "interactionHint": "Suelte el trompo en reposo para ver las cúspides y luego pruebe a empujarlo hacia atrás o hacia adelante." }, - "currentDetails": "Marcador de posición de nutación." + "currentDetailsPattern": "El eje está inclinado {{tilt}} grados, nutando entre {{thetaMin}} y {{thetaMax}} grados. La precesión media es {{precessionHz}} hertz y la nutación {{nutationHz}} hertz.", + "controls": { + "spinRate": "Velocidad de giro propio", + "initialTilt": "Inclinación de suelta", + "releaseMode": "Modo de suelta", + "friction": "Fricción", + "releaseAgain": "Soltar el trompo de nuevo" + } }, "torqueFree": { "screenSummary": { diff --git a/src/i18n/strings_fr.json b/src/i18n/strings_fr.json index 71aee68..98cda9f 100644 --- a/src/i18n/strings_fr.json +++ b/src/i18n/strings_fr.json @@ -19,7 +19,20 @@ "formulaHint": "Ω = τ / (I ω) — plus de spin → précession plus lente" }, "nutation": { - "placeholder": "Écran nutation — bientôt disponible." + "sceneTitle": "Toupie pesante — trajectoire du bout de l'axe", + "graphTitle": "Inclinaison θ vs. temps", + "spinRate": "Rotation propre ω₃ (Hz)", + "initialTilt": "Inclinaison au lâcher θ₀ (°)", + "releaseMode": "Lâcher la toupie", + "releaseCusp": "Au repos — points de rebroussement", + "releaseLoop": "Poussée en arrière — boucles", + "releaseSmooth": "Poussée en avant — ondulations", + "releaseSteady": "À la vitesse régulière — sans oscillation", + "friction": "Frottement (la toupie ralentit)", + "releaseAgain": "Relâcher", + "readoutsTitle": "Mesures", + "belowCritical": "Rotation propre sous la valeur critique — aucune précession régulière n'existe à cette inclinaison, la toupie doit donc nuter.", + "insight": "Une vraie toupie est lâchée, pas lancée à la bonne vitesse : l'axe descend, gagne de la précession, puis remonte." }, "torqueFree": { "placeholder": "Rotation libre — bientôt disponible." @@ -41,11 +54,18 @@ }, "nutation": { "screenSummary": { - "playArea": "Bientôt : nutation.", - "controlArea": "Contrôles à venir.", - "interactionHint": "En construction." + "playArea": "La zone de jeu montre une toupie pesante sur son pivot, la trajectoire tracée par le bout de l'axe, les deux cercles délimitant la bande de nutation et un graphique de l'inclinaison en fonction du temps.", + "controlArea": "La zone de contrôle permet de régler la rotation propre et l'inclinaison au lâcher, de choisir le mode de lâcher, d'activer le frottement et de relâcher la toupie. Les affichages donnent la bande de nutation, les fréquences et la rotation critique.", + "interactionHint": "Lâchez la toupie au repos pour voir les points de rebroussement, puis essayez de la pousser en arrière ou en avant." }, - "currentDetails": "Écran nutation provisoire." + "currentDetailsPattern": "L'axe est incliné de {{tilt}} degrés, nutant entre {{thetaMin}} et {{thetaMax}} degrés. Précession moyenne {{precessionHz}} hertz, nutation {{nutationHz}} hertz.", + "controls": { + "spinRate": "Vitesse de rotation propre", + "initialTilt": "Inclinaison au lâcher", + "releaseMode": "Mode de lâcher", + "friction": "Frottement", + "releaseAgain": "Relâcher la toupie" + } }, "torqueFree": { "screenSummary": { diff --git a/src/nutation-screen/model/NutationModel.ts b/src/nutation-screen/model/NutationModel.ts index 81c4247..b82f255 100644 --- a/src/nutation-screen/model/NutationModel.ts +++ b/src/nutation-screen/model/NutationModel.ts @@ -1,15 +1,266 @@ /** - * NutationModel.ts — placeholder until Screen 2 is implemented. + * NutationModel.ts + * + * Model for Screen 2 — a heavy symmetric top integrated from its full Lagrangian + * equations of motion, instead of Screen 1's idealized steady-precession formula. + * + * The tilt θ is a genuine dynamical variable here, so releasing the top produces + * nutation: the axis dips, gains precession, rises, and repeats. Changing any + * launch parameter (spin, tilt, release mode) re-releases the top, since those are + * initial conditions rather than continuously variable properties of the motion. */ +import { BooleanProperty, DerivedProperty, EnumerationProperty, NumberProperty, Property } from "scenerystack/axon"; import type { TModel } from "scenerystack/joist"; +import { TimeSpeed } from "scenerystack/scenery-phet"; +import { + createReleaseState, + criticalSpinRate, + type HeavyTopParameters, + type HeavyTopState, + isSleepingTopStable, + type NutationBand, + nutationFrequency, + nutationTurningPoints, + type ReleaseMode, + slowPrecessionRate, + steadyPrecessionRates, + stepHeavyTop, + totalEnergy, + verticalAngularMomentum, +} from "../../common/rigid-body/HeavySymmetricTopPhysics.js"; +import { TopTipTrace } from "../../common/rigid-body/TopTipTrace.js"; +import { TimeModel } from "../../common/TimeModel.js"; +import { + DEFAULT_NUTATION_SPIN_RAD_S, + DEFAULT_NUTATION_TILT_RAD, + GRAVITY_MPS2, + NUTATION_COM_DISTANCE_M, + NUTATION_MAX_TILT_RAD, + NUTATION_SAMPLE_INTERVAL_S, + NUTATION_SPIN_DRAG_N_M_S, + NUTATION_SPIN_INERTIA_KG_M2, + NUTATION_TIP_DRAG_N_M_S, + NUTATION_TRACE_CAPACITY, + NUTATION_TRANSVERSE_INERTIA_KG_M2, + NUTATION_WHEEL_MASS_KG, +} from "../../RigidBodyPrecessionConstants.js"; + +/** Slow motion factor, so the ~1.6 Hz nutation can be followed by eye. */ +const SLOW_MOTION_FACTOR = 0.25; export class NutationModel implements TModel { + public readonly timer = new TimeModel(true); + + /** Normal or slow motion. Nutation and precession differ by roughly 5× in rate. */ + public readonly timeSpeedProperty = new EnumerationProperty(TimeSpeed.NORMAL); + + // ── Launch parameters (initial conditions) ────────────────────────────────── + + /** Body spin ω₃ about the symmetry axis at release (rad/s). */ + public readonly spinRateProperty = new NumberProperty(DEFAULT_NUTATION_SPIN_RAD_S); + /** Tilt from the vertical at release (rad). */ + public readonly initialTiltProperty = new NumberProperty(DEFAULT_NUTATION_TILT_RAD, { units: "radians" }); + /** How the top is launched — sets φ̇(0), which selects cusps, loops, or waves. */ + public readonly releaseModeProperty = new Property("cusp"); + /** Viscous drag and spin friction, which damp the nutation and spin the top down. */ + public readonly frictionEnabledProperty = new BooleanProperty(false); + + // ── Dynamical state ───────────────────────────────────────────────────────── + + public readonly thetaProperty = new NumberProperty(DEFAULT_NUTATION_TILT_RAD, { units: "radians" }); + public readonly thetaDotProperty = new NumberProperty(0); + public readonly phiProperty = new NumberProperty(0, { units: "radians" }); + public readonly phiDotProperty = new NumberProperty(0); + public readonly psiProperty = new NumberProperty(0, { units: "radians" }); + public readonly spinProperty = new NumberProperty(DEFAULT_NUTATION_SPIN_RAD_S); + + // ── Derived readouts ──────────────────────────────────────────────────────── + + /** [θ_min, θ_max] band the axis is confined to, from the current invariants. */ + public readonly nutationBandProperty; + /** Mean precession rate measured from the tip trace (rad/s). */ + public readonly meanPrecessionRateProperty; + /** Nutation frequency I₃ω₃/I₁ (rad/s). */ + public readonly nutationFrequencyProperty; + /** Minimum spin for steady precession at the current tilt (rad/s). */ + public readonly criticalSpinProperty; + /** Whether the current spin admits steady precession at the current tilt. */ + public readonly aboveCriticalSpinProperty; + /** Steady precession rate the top would need to hold its release tilt (rad/s). */ + public readonly steadyRateProperty; + /** Total mechanical energy about the pivot (J) — constant while friction is off. */ + public readonly energyProperty; + /** p_φ, angular momentum about the vertical (kg·m²/s) — also constant without friction. */ + public readonly verticalMomentumProperty; + + private readonly trace = new TopTipTrace(NUTATION_TRACE_CAPACITY); + private sampleAccumulator = 0; + + public constructor() { + this.nutationBandProperty = new DerivedProperty( + [this.thetaProperty, this.thetaDotProperty, this.phiDotProperty, this.spinProperty], + (): NutationBand => nutationTurningPoints(this.getParameters(), this.getState()), + ); + + this.meanPrecessionRateProperty = new DerivedProperty([this.phiProperty], () => + this.trace.estimateMeanPrecessionRate(), + ); + + this.nutationFrequencyProperty = new DerivedProperty([this.spinProperty], (spin) => + nutationFrequency(this.getParameters(), spin), + ); + + this.criticalSpinProperty = new DerivedProperty([this.thetaProperty], (theta) => + criticalSpinRate(this.getParameters(), theta), + ); + + this.aboveCriticalSpinProperty = new DerivedProperty( + [this.spinProperty, this.thetaProperty], + (spin, theta) => steadyPrecessionRates(this.getParameters(), spin, theta).exists, + ); + + this.steadyRateProperty = new DerivedProperty([this.spinProperty, this.initialTiltProperty], (spin, tilt) => + slowPrecessionRate(this.getParameters(), spin, tilt), + ); + + this.energyProperty = new DerivedProperty( + [this.thetaProperty, this.thetaDotProperty, this.phiDotProperty, this.spinProperty], + () => totalEnergy(this.getParameters(), this.getState()), + ); + + this.verticalMomentumProperty = new DerivedProperty( + [this.thetaProperty, this.phiDotProperty, this.spinProperty], + () => verticalAngularMomentum(this.getParameters(), this.getState()), + ); + + // Every launch parameter is an initial condition, so changing one re-releases the top. + this.spinRateProperty.lazyLink(() => this.release()); + this.initialTiltProperty.lazyLink(() => this.release()); + this.releaseModeProperty.lazyLink(() => this.release()); + + this.release(); + } + + public getParameters(): HeavyTopParameters { + const friction = this.frictionEnabledProperty.value; + return { + transverseInertia: NUTATION_TRANSVERSE_INERTIA_KG_M2, + spinInertia: NUTATION_SPIN_INERTIA_KG_M2, + mass: NUTATION_WHEEL_MASS_KG, + gravity: GRAVITY_MPS2, + comDistance: NUTATION_COM_DISTANCE_M, + tipDrag: friction ? NUTATION_TIP_DRAG_N_M_S : 0, + spinDrag: friction ? NUTATION_SPIN_DRAG_N_M_S : 0, + maxTilt: NUTATION_MAX_TILT_RAD, + }; + } + + public getState(): HeavyTopState { + return { + theta: this.thetaProperty.value, + thetaDot: this.thetaDotProperty.value, + phi: this.phiProperty.value, + phiDot: this.phiDotProperty.value, + psi: this.psiProperty.value, + spin: this.spinProperty.value, + }; + } + + /** Whether a top spun this fast would stay upright if stood vertically. */ + public isSleepingStable(): boolean { + return isSleepingTopStable(this.getParameters(), this.spinProperty.value); + } + + public getTraceSamples() { + return this.trace.toSamples(); + } + + public getThetaGraphPoints(): Array<{ x: number; y: number }> { + return this.trace.toThetaDegreePoints(); + } + + /** Re-launch the top from the current launch parameters and clear the history. */ + public release(): void { + const state = createReleaseState( + this.getParameters(), + this.spinRateProperty.value, + this.initialTiltProperty.value, + this.releaseModeProperty.value, + ); + + this.thetaProperty.value = state.theta; + this.thetaDotProperty.value = state.thetaDot; + this.phiProperty.value = state.phi; + this.phiDotProperty.value = state.phiDot; + this.psiProperty.value = state.psi; + this.spinProperty.value = state.spin; + + this.trace.clear(); + this.sampleAccumulator = 0; + this.timer.timeProperty.value = 0; + } + + public step(dt: number): void { + if (!this.timer.isPlayingProperty.value) { + return; + } + this.stepOnce(this.timeSpeedProperty.value === TimeSpeed.SLOW ? dt * SLOW_MOTION_FACTOR : dt); + } + + /** Advance the physics by dt regardless of the play/pause state (used by step-forward). */ + public stepOnce(dt: number): void { + this.timer.timeProperty.value += dt; + this.advance(dt); + } + + private advance(dt: number): void { + const next = stepHeavyTop(this.getParameters(), this.getState(), dt); + + this.thetaProperty.value = next.theta; + this.thetaDotProperty.value = next.thetaDot; + this.phiDotProperty.value = next.phiDot; + this.psiProperty.value = next.psi; + this.spinProperty.value = next.spin; + // Set φ last: the mean-precession readout derives from it and reads the trace. + this.sampleAccumulator += dt; + if (this.sampleAccumulator >= NUTATION_SAMPLE_INTERVAL_S) { + this.trace.push(this.timer.timeProperty.value, next.theta, next.phi); + this.sampleAccumulator = 0; + } + this.phiProperty.value = next.phi; + } + public reset(): void { - // Placeholder — Screen 2 will reset nutation state here. + this.timer.reset(); + this.timeSpeedProperty.reset(); + this.spinRateProperty.reset(); + this.initialTiltProperty.reset(); + this.releaseModeProperty.reset(); + this.frictionEnabledProperty.reset(); + this.release(); } - public step(_dt: number): void { - // Placeholder — Screen 2 integrator will advance state here. + public dispose(): void { + this.timer.dispose(); + this.timeSpeedProperty.dispose(); + this.spinRateProperty.dispose(); + this.initialTiltProperty.dispose(); + this.releaseModeProperty.dispose(); + this.frictionEnabledProperty.dispose(); + this.thetaProperty.dispose(); + this.thetaDotProperty.dispose(); + this.phiProperty.dispose(); + this.phiDotProperty.dispose(); + this.psiProperty.dispose(); + this.spinProperty.dispose(); + this.nutationBandProperty.dispose(); + this.meanPrecessionRateProperty.dispose(); + this.nutationFrequencyProperty.dispose(); + this.criticalSpinProperty.dispose(); + this.aboveCriticalSpinProperty.dispose(); + this.steadyRateProperty.dispose(); + this.energyProperty.dispose(); + this.verticalMomentumProperty.dispose(); } } diff --git a/src/nutation-screen/view/NutationAngleGraphNode.ts b/src/nutation-screen/view/NutationAngleGraphNode.ts new file mode 100644 index 0000000..d97b86f --- /dev/null +++ b/src/nutation-screen/view/NutationAngleGraphNode.ts @@ -0,0 +1,253 @@ +/** + * NutationAngleGraphNode.ts + * + * Rolling θ(t) plot. The trace oscillates between two dashed reference lines at + * θ_min and θ_max — the turning points predicted by the effective potential — so + * the graph and the scene's nutation band tell the same story two ways. + */ + +import { Multilink } from "scenerystack/axon"; +import { + ChartRectangle, + ChartTransform, + GridLineSet, + LinearEquationPlot, + LinePlot, + TickLabelSet, + TickMarkSet, +} from "scenerystack/bamboo"; +import { Bounds2, Range, toFixed, Vector2 } from "scenerystack/dot"; +import { Shape } from "scenerystack/kite"; +import { Orientation } from "scenerystack/phet-core"; +import { HBox, Node, Rectangle, Text } from "scenerystack/scenery"; +import { PhetFont } from "scenerystack/scenery-phet"; +import { calculateTickSpacing } from "../../common/view/BambooChartUtils.js"; +import RigidBodyPrecessionColors from "../../RigidBodyPrecessionColors.js"; +import { NUTATION_GRAPH_HEIGHT, NUTATION_GRAPH_WINDOW_S } from "../../RigidBodyPrecessionConstants.js"; +import type { NutationModel } from "../model/NutationModel.js"; + +const CHART_HEIGHT = NUTATION_GRAPH_HEIGHT; +const Y_AXIS_GUTTER = 56; +const X_AXIS_GUTTER = 34; + +const AXIS_FONT = new PhetFont({ size: 11 }); +const LEGEND_FONT = new PhetFont({ size: 10 }); +const TICK_LABEL_FONT = new PhetFont({ size: 10 }); +const TICK_EXTENT = 6; + +const RADIANS_TO_DEGREES = 180 / Math.PI; + +function legendSwatch(color: typeof RigidBodyPrecessionColors.tipTraceColorProperty, dashed = false): Node { + return new Rectangle(0, 0, 18, 3, { + fill: color, + stroke: color, + lineWidth: 1, + lineDash: dashed ? [4, 3] : [], + centerY: 0, + }); +} + +export class NutationAngleGraphNode extends Node { + private readonly chartTransform: ChartTransform; + private readonly dataPlot: LinePlot; + private readonly thetaMinPlot: LinearEquationPlot; + private readonly thetaMaxPlot: LinearEquationPlot; + private readonly verticalGrid: GridLineSet; + private readonly horizontalGrid: GridLineSet; + private readonly xTickMarks: TickMarkSet; + private readonly yTickMarks: TickMarkSet; + private readonly xTickLabels: TickLabelSet; + private readonly yTickLabels: TickLabelSet; + private readonly model: NutationModel; + + public constructor(model: NutationModel, width = 360) { + super(); + this.model = model; + + this.localBounds = new Bounds2(0, 0, width + Y_AXIS_GUTTER, CHART_HEIGHT + X_AXIS_GUTTER); + + this.chartTransform = new ChartTransform({ + viewWidth: width, + viewHeight: CHART_HEIGHT, + modelXRange: new Range(0, NUTATION_GRAPH_WINDOW_S), + modelYRange: new Range(0, 90), + }); + + const chartRectangle = new ChartRectangle(this.chartTransform, { + fill: RigidBodyPrecessionColors.graphBackgroundColorProperty, + stroke: RigidBodyPrecessionColors.panelBorderColorProperty, + lineWidth: 1, + cornerRadius: 4, + }); + + const initialXSpacing = calculateTickSpacing(NUTATION_GRAPH_WINDOW_S); + const initialYSpacing = 15; + + this.verticalGrid = new GridLineSet(this.chartTransform, Orientation.VERTICAL, initialXSpacing, { + stroke: RigidBodyPrecessionColors.graphGridColorProperty, + lineWidth: 0.5, + }); + this.horizontalGrid = new GridLineSet(this.chartTransform, Orientation.HORIZONTAL, initialYSpacing, { + stroke: RigidBodyPrecessionColors.graphGridColorProperty, + lineWidth: 0.5, + }); + + // Horizontal reference lines (slope 0) at the two turning points. + this.thetaMinPlot = new LinearEquationPlot(this.chartTransform, 0, 0, { + stroke: RigidBodyPrecessionColors.nutationBandColorProperty, + lineWidth: 1.5, + lineDash: [6, 4], + opacity: 0.85, + }); + this.thetaMaxPlot = new LinearEquationPlot(this.chartTransform, 0, 0, { + stroke: RigidBodyPrecessionColors.nutationBandColorProperty, + lineWidth: 1.5, + lineDash: [6, 4], + opacity: 0.85, + }); + + this.dataPlot = new LinePlot(this.chartTransform, [], { + stroke: RigidBodyPrecessionColors.tipTraceColorProperty, + lineWidth: 2.5, + }); + + const clippedChartContent = new Node({ + clipArea: Shape.bounds(new Bounds2(0, 0, width, CHART_HEIGHT)), + children: [ + chartRectangle, + this.verticalGrid, + this.horizontalGrid, + this.thetaMinPlot, + this.thetaMaxPlot, + this.dataPlot, + ], + }); + + this.xTickMarks = new TickMarkSet(this.chartTransform, Orientation.HORIZONTAL, initialXSpacing, { + edge: "min", + extent: TICK_EXTENT, + stroke: RigidBodyPrecessionColors.panelBorderColorProperty, + }); + this.yTickMarks = new TickMarkSet(this.chartTransform, Orientation.VERTICAL, initialYSpacing, { + edge: "min", + extent: TICK_EXTENT, + stroke: RigidBodyPrecessionColors.panelBorderColorProperty, + }); + this.xTickLabels = new TickLabelSet(this.chartTransform, Orientation.HORIZONTAL, initialXSpacing, { + edge: "min", + extent: TICK_EXTENT, + createLabel: (value: number) => + new Text(toFixed(value, 1), { + font: TICK_LABEL_FONT, + fill: RigidBodyPrecessionColors.textColorProperty, + maxWidth: 36, + }), + }); + this.yTickLabels = new TickLabelSet(this.chartTransform, Orientation.VERTICAL, initialYSpacing, { + edge: "min", + extent: TICK_EXTENT, + createLabel: (value: number) => + new Text(toFixed(value, 0), { + font: TICK_LABEL_FONT, + fill: RigidBodyPrecessionColors.textColorProperty, + maxWidth: 36, + }), + }); + + const chartNode = new Node({ + translation: new Vector2(Y_AXIS_GUTTER, 0), + children: [clippedChartContent, this.xTickMarks, this.yTickMarks, this.xTickLabels, this.yTickLabels], + }); + this.addChild(chartNode); + + this.addChild( + new Text("θ (°)", { + font: AXIS_FONT, + fill: RigidBodyPrecessionColors.textColorProperty, + rotation: -Math.PI / 2, + centerY: CHART_HEIGHT / 2, + // Clear of the y tick labels, which sit just left of the chart edge. + right: Y_AXIS_GUTTER - 24, + }), + ); + this.addChild( + new Text("t (s)", { + font: AXIS_FONT, + fill: RigidBodyPrecessionColors.textColorProperty, + centerX: Y_AXIS_GUTTER + width / 2, + top: CHART_HEIGHT + 18, + }), + ); + + this.addChild( + new HBox({ + spacing: 12, + align: "center", + children: [ + new HBox({ + spacing: 4, + align: "center", + children: [ + legendSwatch(RigidBodyPrecessionColors.tipTraceColorProperty), + new Text("θ(t)", { font: LEGEND_FONT, fill: RigidBodyPrecessionColors.textColorProperty }), + ], + }), + new HBox({ + spacing: 4, + align: "center", + children: [ + legendSwatch(RigidBodyPrecessionColors.nutationBandColorProperty, true), + new Text("turning points", { + font: LEGEND_FONT, + fill: RigidBodyPrecessionColors.textColorProperty, + }), + ], + }), + ], + // On the axis-title row, clear of the x tick labels directly under the chart. + right: Y_AXIS_GUTTER + width, + top: CHART_HEIGHT + 16, + }), + ); + + Multilink.multilink([model.thetaProperty, model.nutationBandProperty], () => this.updatePlot()); + this.updatePlot(); + } + + private updatePlot(): void { + const points = this.model.getThetaGraphPoints(); + const band = this.model.nutationBandProperty.value; + const thetaMinDegrees = band.thetaMin * RADIANS_TO_DEGREES; + const thetaMaxDegrees = band.thetaMax * RADIANS_TO_DEGREES; + + const latestTime = points[points.length - 1]?.x ?? 0; + const minTime = Math.max(0, latestTime - NUTATION_GRAPH_WINDOW_S); + const maxTime = Math.max(NUTATION_GRAPH_WINDOW_S, latestTime); + const xRange = new Range(minTime, maxTime); + this.chartTransform.setModelXRange(xRange); + + const visibleAngles = points.filter((point) => point.x >= minTime).map((point) => point.y); + const lowest = Math.min(...visibleAngles, thetaMinDegrees); + const highest = Math.max(...visibleAngles, thetaMaxDegrees); + const padding = Math.max(3, 0.12 * (highest - lowest)); + const yRange = new Range(Math.max(0, lowest - padding), Math.min(180, highest + padding)); + this.chartTransform.setModelYRange(yRange); + + const xSpacing = calculateTickSpacing(xRange.getLength()); + const ySpacing = calculateTickSpacing(yRange.getLength()); + this.verticalGrid.setSpacing(xSpacing); + this.horizontalGrid.setSpacing(ySpacing); + this.xTickMarks.setSpacing(xSpacing); + this.yTickMarks.setSpacing(ySpacing); + this.xTickLabels.setSpacing(xSpacing); + this.yTickLabels.setSpacing(ySpacing); + + this.dataPlot.setDataSet(points.map((point) => new Vector2(point.x, point.y))); + + this.thetaMinPlot.b = thetaMinDegrees; + this.thetaMaxPlot.b = thetaMaxDegrees; + const bandVisible = thetaMaxDegrees - thetaMinDegrees > 0.05; + this.thetaMinPlot.visible = bandVisible; + this.thetaMaxPlot.visible = bandVisible; + } +} diff --git a/src/nutation-screen/view/NutationControlPanel.ts b/src/nutation-screen/view/NutationControlPanel.ts new file mode 100644 index 0000000..e208bb4 --- /dev/null +++ b/src/nutation-screen/view/NutationControlPanel.ts @@ -0,0 +1,254 @@ +/** + * NutationControlPanel.ts + * + * Launch controls for the heavy top — spin, release tilt, release mode, friction — + * plus readouts of the nutation band, the two frequencies, and the critical spin. + * Every slider here is an initial condition, so moving one re-releases the top. + */ + +import { DerivedProperty, NumberProperty, type TReadOnlyProperty } from "scenerystack/axon"; +import { clamp, Dimension2, Range, toFixed } from "scenerystack/dot"; +import { type Color, HBox, Line, type Node, Text, VBox } from "scenerystack/scenery"; +import { NumberControl, PhetFont } from "scenerystack/scenery-phet"; +import { Checkbox, ComboBox, RectangularPushButton } from "scenerystack/sun"; +import type { ReleaseMode } from "../../common/rigid-body/HeavySymmetricTopPhysics.js"; +import { + FLAT_RECTANGULAR_BUTTON_OPTIONS, + LIGHT_SURFACE_TEXT_FILL, + SIM_COMBO_BOX_OPTIONS, +} from "../../common/SimButtonOptions.js"; +import { SimPanel } from "../../common/SimPanel.js"; +import { StringManager } from "../../i18n/StringManager.js"; +import RigidBodyPrecessionColors from "../../RigidBodyPrecessionColors.js"; +import { NUTATION_PANEL_WIDTH, NUTATION_SPIN_RANGE, NUTATION_TILT_RANGE } from "../../RigidBodyPrecessionConstants.js"; +import type { NutationModel } from "../model/NutationModel.js"; + +const TITLE_FONT = new PhetFont({ size: 13, weight: "bold" }); +const SECTION_FONT = new PhetFont({ size: 12, weight: "bold" }); +const READOUT_FONT = new PhetFont({ size: 13 }); +const READOUT_LABEL_FONT = new PhetFont({ size: 12 }); +const SLIDER_WIDTH = NUTATION_PANEL_WIDTH - 56; + +const HZ_PER_RAD_S = 1 / (2 * Math.PI); +const DEGREES_PER_RADIAN = 180 / Math.PI; + +/** + * A display-unit proxy for a model property (rad/s → Hz, rad → degrees). Edits flow + * both ways; the guard keeps the model's own resets from echoing back as user edits. + */ +function createUnitProxy(property: NumberProperty, scale: number, range: Range, units: "Hz" | "°"): NumberProperty { + const proxy = new NumberProperty(clamp(property.value * scale, range.min, range.max), { units }); + let suppress = false; + property.link((value) => { + if (!suppress) { + proxy.value = clamp(value * scale, range.min, range.max); + } + }); + proxy.lazyLink((value) => { + suppress = true; + property.value = value / scale; + suppress = false; + }); + return proxy; +} + +function createNumberControl( + title: TReadOnlyProperty, + property: NumberProperty, + range: Range, + delta: number, + decimalPlaces: number, + accessibleName: TReadOnlyProperty, +): NumberControl { + const titleNode = new Text(title, { font: TITLE_FONT, fill: RigidBodyPrecessionColors.textColorProperty }); + return new NumberControl(titleNode, property, range, { + delta, + layoutFunction: NumberControl.createLayoutFunction1({ align: "center", ySpacing: 2 }), + numberDisplayOptions: { + decimalPlaces, + textOptions: { font: READOUT_FONT }, + }, + sliderOptions: { + trackSize: new Dimension2(SLIDER_WIDTH, 4), + thumbSize: new Dimension2(14, 22), + }, + arrowButtonOptions: FLAT_RECTANGULAR_BUTTON_OPTIONS, + accessibleName, + }); +} + +function readoutRow( + label: string, + valueProperty: TReadOnlyProperty, + colorProperty: TReadOnlyProperty, +): Node { + return new HBox({ + spacing: 8, + children: [ + new Text(label, { font: READOUT_LABEL_FONT, fill: RigidBodyPrecessionColors.textColorProperty }), + new Text(valueProperty, { font: READOUT_FONT, fill: colorProperty }), + ], + }); +} + +export class NutationControlPanel extends SimPanel { + public constructor(model: NutationModel, listParent: Node) { + const strings = StringManager.getInstance().getNutationStrings(); + const a11y = StringManager.getInstance().getNutationA11yStrings(); + + const spinHzRange = new Range(NUTATION_SPIN_RANGE.min * HZ_PER_RAD_S, NUTATION_SPIN_RANGE.max * HZ_PER_RAD_S); + const spinControl = createNumberControl( + strings.spinRateStringProperty, + createUnitProxy(model.spinRateProperty, HZ_PER_RAD_S, spinHzRange, "Hz"), + spinHzRange, + 0.05, + 2, + a11y.controls.spinRateStringProperty, + ); + + const tiltDegreeRange = new Range( + Math.round(NUTATION_TILT_RANGE.min * DEGREES_PER_RADIAN), + Math.round(NUTATION_TILT_RANGE.max * DEGREES_PER_RADIAN), + ); + const tiltControl = createNumberControl( + strings.initialTiltStringProperty, + createUnitProxy(model.initialTiltProperty, DEGREES_PER_RADIAN, tiltDegreeRange, "°"), + tiltDegreeRange, + 1, + 0, + a11y.controls.initialTiltStringProperty, + ); + + const releaseLabel = new Text(strings.releaseModeStringProperty, { + font: TITLE_FONT, + fill: RigidBodyPrecessionColors.textColorProperty, + }); + + const releaseItems: Array<{ + value: ReleaseMode; + createNode: () => Node; + accessibleName: TReadOnlyProperty; + }> = [ + { + value: "cusp", + createNode: () => + new Text(strings.releaseCuspStringProperty, { font: READOUT_FONT, fill: LIGHT_SURFACE_TEXT_FILL }), + accessibleName: strings.releaseCuspStringProperty, + }, + { + value: "loop", + createNode: () => + new Text(strings.releaseLoopStringProperty, { font: READOUT_FONT, fill: LIGHT_SURFACE_TEXT_FILL }), + accessibleName: strings.releaseLoopStringProperty, + }, + { + value: "smooth", + createNode: () => + new Text(strings.releaseSmoothStringProperty, { font: READOUT_FONT, fill: LIGHT_SURFACE_TEXT_FILL }), + accessibleName: strings.releaseSmoothStringProperty, + }, + { + value: "steady", + createNode: () => + new Text(strings.releaseSteadyStringProperty, { font: READOUT_FONT, fill: LIGHT_SURFACE_TEXT_FILL }), + accessibleName: strings.releaseSteadyStringProperty, + }, + ]; + + const releaseComboBox = new ComboBox(model.releaseModeProperty, releaseItems, listParent, { + ...SIM_COMBO_BOX_OPTIONS, + accessibleName: a11y.controls.releaseModeStringProperty, + listPosition: "below", + }); + + const frictionCheckbox = new Checkbox( + model.frictionEnabledProperty, + new Text(strings.frictionStringProperty, { + font: READOUT_FONT, + fill: RigidBodyPrecessionColors.textColorProperty, + maxWidth: NUTATION_PANEL_WIDTH - 60, + }), + { + accessibleName: a11y.controls.frictionStringProperty, + boxWidth: 18, + }, + ); + + const releaseButton = new RectangularPushButton({ + ...FLAT_RECTANGULAR_BUTTON_OPTIONS, + content: new Text(strings.releaseAgainStringProperty, { font: READOUT_FONT, fill: LIGHT_SURFACE_TEXT_FILL }), + baseColor: RigidBodyPrecessionColors.controlSurfaceColorProperty, + accessibleName: a11y.controls.releaseAgainStringProperty, + listener: () => model.release(), + }); + + const bandValueProperty = new DerivedProperty([model.nutationBandProperty], (band) => { + const min = toFixed((band.thetaMin * 180) / Math.PI, 1); + const max = toFixed((band.thetaMax * 180) / Math.PI, 1); + return `${min}° – ${max}°`; + }); + const nutationValueProperty = new DerivedProperty( + [model.nutationFrequencyProperty], + (rate) => `${toFixed(rate * HZ_PER_RAD_S, 2)} Hz`, + ); + const precessionValueProperty = new DerivedProperty( + [model.meanPrecessionRateProperty], + (rate) => `${toFixed(rate * HZ_PER_RAD_S, 3)} Hz`, + ); + const criticalValueProperty = new DerivedProperty( + [model.criticalSpinProperty], + (rate) => `${toFixed(rate * HZ_PER_RAD_S, 2)} Hz`, + ); + const criticalColorProperty = new DerivedProperty([model.aboveCriticalSpinProperty], (above) => + above + ? RigidBodyPrecessionColors.precessionColorProperty.value + : RigidBodyPrecessionColors.warningColorProperty.value, + ); + + const separator = new Line(0, 0, NUTATION_PANEL_WIDTH - 40, 0, { + stroke: RigidBodyPrecessionColors.panelBorderColorProperty, + lineWidth: 1, + }); + + const readoutHeader = new Text(strings.readoutsTitleStringProperty, { + font: SECTION_FONT, + fill: RigidBodyPrecessionColors.accentColorProperty, + }); + + const criticalWarning = new Text(strings.belowCriticalStringProperty, { + font: new PhetFont({ size: 11 }), + fill: RigidBodyPrecessionColors.warningColorProperty, + maxWidth: NUTATION_PANEL_WIDTH - 40, + visibleProperty: new DerivedProperty([model.aboveCriticalSpinProperty], (above) => !above), + }); + + const insight = new Text(strings.insightStringProperty, { + font: new PhetFont({ size: 11 }), + fill: RigidBodyPrecessionColors.textColorProperty, + maxWidth: NUTATION_PANEL_WIDTH - 40, + opacity: 0.85, + }); + + const content = new VBox({ + spacing: 10, + align: "left", + children: [ + spinControl, + tiltControl, + new VBox({ spacing: 4, align: "left", children: [releaseLabel, releaseComboBox] }), + frictionCheckbox, + releaseButton, + separator, + readoutHeader, + readoutRow("θ range", bandValueProperty, RigidBodyPrecessionColors.nutationBandColorProperty), + readoutRow("f_nut", nutationValueProperty, RigidBodyPrecessionColors.tipTraceColorProperty), + readoutRow("Ω_mean", precessionValueProperty, RigidBodyPrecessionColors.precessionColorProperty), + readoutRow("ω₃ min", criticalValueProperty, criticalColorProperty), + criticalWarning, + insight, + ], + }); + + super(content, { minWidth: NUTATION_PANEL_WIDTH }); + } +} diff --git a/src/nutation-screen/view/NutationKeyboardHelpContent.ts b/src/nutation-screen/view/NutationKeyboardHelpContent.ts index ca51f9d..40dd2c9 100644 --- a/src/nutation-screen/view/NutationKeyboardHelpContent.ts +++ b/src/nutation-screen/view/NutationKeyboardHelpContent.ts @@ -2,10 +2,18 @@ * NutationKeyboardHelpContent.ts */ -import { BasicActionsKeyboardHelpSection, TwoColumnKeyboardHelpContent } from "scenerystack/scenery-phet"; +import { + BasicActionsKeyboardHelpSection, + ComboBoxKeyboardHelpSection, + SliderControlsKeyboardHelpSection, + TwoColumnKeyboardHelpContent, +} from "scenerystack/scenery-phet"; export class NutationKeyboardHelpContent extends TwoColumnKeyboardHelpContent { public constructor() { - super([new BasicActionsKeyboardHelpSection()], []); + super( + [new SliderControlsKeyboardHelpSection()], + [new ComboBoxKeyboardHelpSection(), new BasicActionsKeyboardHelpSection()], + ); } } diff --git a/src/nutation-screen/view/NutationScreenSummaryContent.ts b/src/nutation-screen/view/NutationScreenSummaryContent.ts index 70dfc47..7fb3d35 100644 --- a/src/nutation-screen/view/NutationScreenSummaryContent.ts +++ b/src/nutation-screen/view/NutationScreenSummaryContent.ts @@ -1,17 +1,39 @@ /** * NutationScreenSummaryContent.ts + * + * Live accessible summary of the top's state: where the axis is, how wide the + * nutation band is, and how fast it is precessing on average. */ +import { DerivedProperty, PatternStringProperty } from "scenerystack/axon"; +import { toFixed } from "scenerystack/dot"; import { ScreenSummaryContent } from "scenerystack/sim"; import { StringManager } from "../../i18n/StringManager.js"; +import type { NutationModel } from "../model/NutationModel.js"; + +const DEGREES_PER_RADIAN = 180 / Math.PI; +const HZ_PER_RAD_S = 1 / (2 * Math.PI); export class NutationScreenSummaryContent extends ScreenSummaryContent { - public constructor() { + public constructor(model: NutationModel) { const a11y = StringManager.getInstance().getNutationA11yStrings(); + + const currentDetails = new PatternStringProperty(a11y.currentDetailsPatternStringProperty, { + tilt: new DerivedProperty([model.thetaProperty], (theta) => toFixed(theta * DEGREES_PER_RADIAN, 0)), + thetaMin: new DerivedProperty([model.nutationBandProperty], (band) => + toFixed(band.thetaMin * DEGREES_PER_RADIAN, 0), + ), + thetaMax: new DerivedProperty([model.nutationBandProperty], (band) => + toFixed(band.thetaMax * DEGREES_PER_RADIAN, 0), + ), + precessionHz: new DerivedProperty([model.meanPrecessionRateProperty], (rate) => toFixed(rate * HZ_PER_RAD_S, 2)), + nutationHz: new DerivedProperty([model.nutationFrequencyProperty], (rate) => toFixed(rate * HZ_PER_RAD_S, 2)), + }); + super({ playAreaContent: a11y.screenSummary.playAreaStringProperty, controlAreaContent: a11y.screenSummary.controlAreaStringProperty, - currentDetailsContent: a11y.currentDetailsStringProperty, + currentDetailsContent: currentDetails, interactionHintContent: a11y.screenSummary.interactionHintStringProperty, }); } diff --git a/src/nutation-screen/view/NutationScreenView.ts b/src/nutation-screen/view/NutationScreenView.ts index 4b96b8e..b0378f0 100644 --- a/src/nutation-screen/view/NutationScreenView.ts +++ b/src/nutation-screen/view/NutationScreenView.ts @@ -1,22 +1,35 @@ /** - * NutationScreenView.ts — placeholder until Screen 2 is implemented. + * NutationScreenView.ts + * + * Layout for Screen 2: the top and its tip trace on the left, θ(t) below it, and + * the launch controls on the right. */ -import { Rectangle, Text } from "scenerystack/scenery"; -import { ResetAllButton } from "scenerystack/scenery-phet"; +import { Node, Rectangle, VBox } from "scenerystack/scenery"; +import { ResetAllButton, TimeControlNode, TimeSpeed } from "scenerystack/scenery-phet"; import type { ScreenViewOptions } from "scenerystack/sim"; import { ScreenView } from "scenerystack/sim"; -import { FLAT_RESET_ALL_BUTTON_OPTIONS } from "../../common/SimButtonOptions.js"; +import { + FLAT_PLAY_PAUSE_STEP_BUTTON_OPTIONS, + FLAT_RESET_ALL_BUTTON_OPTIONS, + TIME_CONTROL_SPEED_RADIO_OPTIONS, +} from "../../common/SimButtonOptions.js"; +import { PlayAreaPanel } from "../../common/view/PlayAreaPanel.js"; import { StringManager } from "../../i18n/StringManager.js"; import RigidBodyPrecessionColors from "../../RigidBodyPrecessionColors.js"; import { SCREEN_VIEW_MARGIN } from "../../RigidBodyPrecessionConstants.js"; import type { NutationModel } from "../model/NutationModel.js"; +import { NutationAngleGraphNode } from "./NutationAngleGraphNode.js"; +import { NutationControlPanel } from "./NutationControlPanel.js"; import { NutationScreenSummaryContent } from "./NutationScreenSummaryContent.js"; +import { TopSceneNode } from "./TopSceneNode.js"; + +const BOTTOM_CHROME_HEIGHT = 56; export class NutationScreenView extends ScreenView { public constructor(model: NutationModel, options?: ScreenViewOptions) { super({ - screenSummaryContent: new NutationScreenSummaryContent(), + screenSummaryContent: new NutationScreenSummaryContent(model), ...options, }); @@ -26,13 +39,56 @@ export class NutationScreenView extends ScreenView { }), ); - const placeholder = new Text(StringManager.getInstance().getNutationStrings().placeholderStringProperty, { - font: "24px sans-serif", - fill: RigidBodyPrecessionColors.textColorProperty, - center: this.layoutBounds.center, - maxWidth: this.layoutBounds.width - 80, + const strings = StringManager.getInstance().getNutationStrings(); + + // The combo-box list must sit above every other node in the view. + const comboBoxListParent = new Node(); + + const controlPanel = new NutationControlPanel(model, comboBoxListParent); + controlPanel.right = this.layoutBounds.maxX - SCREEN_VIEW_MARGIN; + controlPanel.top = SCREEN_VIEW_MARGIN; + + const playAreaRight = controlPanel.left - SCREEN_VIEW_MARGIN; + const playAreaWidth = playAreaRight - SCREEN_VIEW_MARGIN; + const playAreaBottom = this.layoutBounds.maxY - SCREEN_VIEW_MARGIN - BOTTOM_CHROME_HEIGHT; + + const scenePanel = new PlayAreaPanel(strings.sceneTitleStringProperty, new TopSceneNode(model)); + const graphPanel = new PlayAreaPanel( + strings.graphTitleStringProperty, + new NutationAngleGraphNode(model, Math.max(280, playAreaWidth - 60)), + ); + + const playColumn = new VBox({ + spacing: 8, + align: "left", + children: [scenePanel, graphPanel], + }); + + // Scale the column down if the two panels together overflow the play area. + const availableHeight = playAreaBottom - SCREEN_VIEW_MARGIN; + const widthScale = playColumn.width > playAreaWidth ? playAreaWidth / playColumn.width : 1; + const heightScale = playColumn.height > availableHeight ? availableHeight / playColumn.height : 1; + const scale = Math.min(widthScale, heightScale); + if (scale < 1) { + playColumn.setScaleMagnitude(scale); + } + playColumn.left = SCREEN_VIEW_MARGIN; + playColumn.top = SCREEN_VIEW_MARGIN; + + const timeControl = new TimeControlNode(model.timer.isPlayingProperty, { + ...TIME_CONTROL_SPEED_RADIO_OPTIONS, + timeSpeedProperty: model.timeSpeedProperty, + timeSpeeds: [TimeSpeed.NORMAL, TimeSpeed.SLOW], + playPauseStepButtonOptions: { + ...FLAT_PLAY_PAUSE_STEP_BUTTON_OPTIONS, + stepForwardButtonOptions: { + ...FLAT_PLAY_PAUSE_STEP_BUTTON_OPTIONS.stepForwardButtonOptions, + listener: () => model.stepOnce(1 / 60), + }, + }, + left: SCREEN_VIEW_MARGIN, + bottom: this.layoutBounds.maxY - SCREEN_VIEW_MARGIN, }); - this.addChild(placeholder); const resetAllButton = new ResetAllButton({ ...FLAT_RESET_ALL_BUTTON_OPTIONS, @@ -43,10 +99,25 @@ export class NutationScreenView extends ScreenView { right: this.layoutBounds.maxX - SCREEN_VIEW_MARGIN, bottom: this.layoutBounds.maxY - SCREEN_VIEW_MARGIN, }); + + this.addChild(playColumn); + this.addChild(controlPanel); + this.addChild(timeControl); this.addChild(resetAllButton); + this.addChild(comboBoxListParent); + + this.addChild( + new Node({ + pdomOrder: [controlPanel, timeControl, resetAllButton], + }), + ); } public reset(): void { - // No view-local state on the placeholder screen. + // View state is model-driven; nothing extra to reset. + } + + public override step(_dt: number): void { + // Animation updates via Multilink on model properties. } } diff --git a/src/nutation-screen/view/TopSceneNode.ts b/src/nutation-screen/view/TopSceneNode.ts new file mode 100644 index 0000000..2f8d1fa --- /dev/null +++ b/src/nutation-screen/view/TopSceneNode.ts @@ -0,0 +1,215 @@ +/** + * TopSceneNode.ts + * + * Perspective view of the heavy symmetric top: the pivoted wheel, the path its + * axle tip traces on the sphere, and the [θ_min, θ_max] band that path is + * confined to. The trace is the payoff of the full dynamics — cusps, loops, and + * smooth waves are all the same equations released three different ways. + */ + +import { Multilink } from "scenerystack/axon"; +import { Bounds2, Vector2 } from "scenerystack/dot"; +import { Shape } from "scenerystack/kite"; +import { Circle, Line, Node, Path, Rectangle, Text } from "scenerystack/scenery"; +import { PhetFont } from "scenerystack/scenery-phet"; +import { + projectAxisPoint, + projectRimPoint, + projectTiltCircle, + type TopProjection, + wheelRimPoints, + wheelSilhouette, +} from "../../common/view/TopProjection.js"; +import RigidBodyPrecessionColors from "../../RigidBodyPrecessionColors.js"; +import { + NUTATION_COM_DISTANCE_M, + NUTATION_TRACE_DRAW_SAMPLES, + NUTATION_WHEEL_RADIUS_M, +} from "../../RigidBodyPrecessionConstants.js"; +import type { NutationModel } from "../model/NutationModel.js"; + +export const TOP_SCENE_WIDTH = 420; +export const TOP_SCENE_HEIGHT = 345; + +/** + * Drawn length of the axle beyond the wheel (m). Visual only — the dynamics depend on + * the inertias and the center-of-mass distance, not on where the axle ends. It is long + * enough that the traced tip circle clears the wheel rim at every tilt in range. + */ +const AXLE_LENGTH_M = 0.85; + +const PROJECTION: TopProjection = { + pivot: new Vector2(210, 235), + pxPerM: 165, + perspective: 0.45, +}; + +const LABEL_FONT = new PhetFont({ size: 10 }); + +function bandShape(theta: number): Shape { + const circle = projectTiltCircle(PROJECTION, theta, AXLE_LENGTH_M); + return Shape.ellipse(circle.center.x, circle.center.y, Math.max(1, circle.radiusX), Math.max(1, circle.radiusY), 0); +} + +export class TopSceneNode extends Node { + public constructor(model: NutationModel) { + super(); + this.localBounds = new Bounds2(0, 0, TOP_SCENE_WIDTH, TOP_SCENE_HEIGHT); + + const pivot = PROJECTION.pivot; + + const ground = new Path(Shape.ellipse(pivot.x, pivot.y + 62, 150, 22, 0), { + fill: RigidBodyPrecessionColors.sceneGroundColorProperty, + stroke: RigidBodyPrecessionColors.panelBorderColorProperty, + lineWidth: 1, + }); + this.addChild(ground); + + // Vertical reference: θ is measured from this line. + const verticalAxis = new Line(pivot.x, pivot.y - 200, pivot.x, pivot.y, { + stroke: RigidBodyPrecessionColors.textColorProperty, + lineWidth: 1, + lineDash: [4, 5], + opacity: 0.45, + }); + this.addChild(verticalAxis); + + // Turning-point circles: the axis tip can never leave the band between them. + const bandMin = new Path(new Shape(), { + stroke: RigidBodyPrecessionColors.nutationBandColorProperty, + lineWidth: 1.5, + lineDash: [6, 5], + opacity: 0.9, + }); + const bandMax = new Path(new Shape(), { + stroke: RigidBodyPrecessionColors.nutationBandColorProperty, + lineWidth: 1.5, + lineDash: [6, 5], + opacity: 0.9, + }); + this.addChild(bandMin); + this.addChild(bandMax); + + const tipTrace = new Path(new Shape(), { + stroke: RigidBodyPrecessionColors.tipTraceColorProperty, + lineWidth: 2, + lineJoin: "round", + }); + this.addChild(tipTrace); + + const support = new Rectangle(pivot.x - 5, pivot.y, 10, 56, { + fill: RigidBodyPrecessionColors.gyroscopeColorProperty, + cornerRadius: 2, + }); + const supportBase = new Rectangle(pivot.x - 34, pivot.y + 52, 68, 10, { + fill: RigidBodyPrecessionColors.panelBorderColorProperty, + cornerRadius: 3, + }); + this.addChild(support); + this.addChild(supportBase); + + const axle = new Line(0, 0, 0, 0, { + stroke: RigidBodyPrecessionColors.gyroscopeColorProperty, + lineWidth: 4, + lineCap: "round", + }); + this.addChild(axle); + + const hub = new Path(new Shape(), { + fill: RigidBodyPrecessionColors.gyroscopeColorProperty, + opacity: 0.55, + }); + this.addChild(hub); + + const wheel = new Path(new Shape(), { + fill: RigidBodyPrecessionColors.wheelFillColorProperty, + stroke: RigidBodyPrecessionColors.angularMomentumColorProperty, + lineWidth: 3, + }); + this.addChild(wheel); + + const spinSpoke = new Line(0, 0, 0, 0, { + stroke: RigidBodyPrecessionColors.weightColorProperty, + lineWidth: 3, + lineCap: "round", + }); + this.addChild(spinSpoke); + + const pivotDot = new Circle(6, { + fill: RigidBodyPrecessionColors.accentColorProperty, + stroke: RigidBodyPrecessionColors.textColorProperty, + lineWidth: 1, + center: pivot, + }); + this.addChild(pivotDot); + + const tipDot = new Circle(5, { fill: RigidBodyPrecessionColors.tipTraceColorProperty }); + this.addChild(tipDot); + + const bandLabel = new Text("θ band", { + font: LABEL_FONT, + fill: RigidBodyPrecessionColors.nutationBandColorProperty, + opacity: 0.9, + }); + this.addChild(bandLabel); + + const update = (): void => { + const theta = model.thetaProperty.value; + const phi = model.phiProperty.value; + const psi = model.psiProperty.value; + const band = model.nutationBandProperty.value; + + const tip = projectAxisPoint(PROJECTION, theta, phi, AXLE_LENGTH_M); + axle.setPoint1(pivot.x, pivot.y); + axle.setPoint2(tip.x, tip.y); + tipDot.center = tip; + + const rim = wheelRimPoints(PROJECTION, theta, phi, NUTATION_COM_DISTANCE_M, NUTATION_WHEEL_RADIUS_M); + wheel.shape = Shape.polygon(rim); + + const [rimA, rimB] = wheelSilhouette(PROJECTION, theta, phi, NUTATION_COM_DISTANCE_M, NUTATION_WHEEL_RADIUS_M); + hub.shape = Shape.polygon([pivot, rimA, rimB]); + + const spokeTip = projectRimPoint(PROJECTION, theta, phi, NUTATION_COM_DISTANCE_M, NUTATION_WHEEL_RADIUS_M, psi); + const wheelCenter = projectAxisPoint(PROJECTION, theta, phi, NUTATION_COM_DISTANCE_M); + spinSpoke.setPoint1(wheelCenter.x, wheelCenter.y); + spinSpoke.setPoint2(spokeTip.x, spokeTip.y); + + bandMin.shape = bandShape(band.thetaMin); + bandMax.shape = bandShape(band.thetaMax); + const bandVisible = band.thetaMax - band.thetaMin > 1e-3; + bandMin.visible = bandVisible; + bandMax.visible = bandVisible; + bandLabel.visible = bandVisible; + if (bandVisible) { + const outer = projectTiltCircle(PROJECTION, band.thetaMax, AXLE_LENGTH_M); + bandLabel.right = outer.center.x - outer.radiusX - 4; + bandLabel.centerY = outer.center.y; + } + + const samples = model.getTraceSamples(); + const firstDrawn = Math.max(0, samples.length - NUTATION_TRACE_DRAW_SAMPLES); + if (samples.length - firstDrawn > 1) { + const traceShape = new Shape(); + for (let i = firstDrawn; i < samples.length; i++) { + const sample = samples[i]; + if (!sample) { + continue; + } + const point = projectAxisPoint(PROJECTION, sample.theta, sample.phi, AXLE_LENGTH_M); + if (i === firstDrawn) { + traceShape.moveToPoint(point); + } else { + traceShape.lineToPoint(point); + } + } + tipTrace.shape = traceShape; + } else { + tipTrace.shape = new Shape(); + } + }; + + Multilink.multilink([model.thetaProperty, model.phiProperty, model.psiProperty], update); + update(); + } +} diff --git a/src/steady-precession-screen/view/SteadyPrecessionScreenView.ts b/src/steady-precession-screen/view/SteadyPrecessionScreenView.ts index 0e2dbb0..5275f91 100644 --- a/src/steady-precession-screen/view/SteadyPrecessionScreenView.ts +++ b/src/steady-precession-screen/view/SteadyPrecessionScreenView.ts @@ -11,12 +11,12 @@ import { FLAT_RESET_ALL_BUTTON_OPTIONS, TIME_CONTROL_SPEED_RADIO_OPTIONS, } from "../../common/SimButtonOptions.js"; +import { PlayAreaPanel } from "../../common/view/PlayAreaPanel.js"; import { StringManager } from "../../i18n/StringManager.js"; import RigidBodyPrecessionColors from "../../RigidBodyPrecessionColors.js"; import { SCREEN_VIEW_MARGIN } from "../../RigidBodyPrecessionConstants.js"; import type { SteadyPrecessionModel } from "../model/SteadyPrecessionModel.js"; import { GyroscopeSceneNode } from "./GyroscopeSceneNode.js"; -import { PlayAreaPanel } from "./PlayAreaPanel.js"; import { PrecessionAngleGraphNode } from "./PrecessionAngleGraphNode.js"; import { SteadyPrecessionControlPanel } from "./SteadyPrecessionControlPanel.js"; import { SteadyPrecessionScreenSummaryContent } from "./SteadyPrecessionScreenSummaryContent.js"; diff --git a/tests/HeavySymmetricTopPhysics.test.ts b/tests/HeavySymmetricTopPhysics.test.ts new file mode 100644 index 0000000..74ae151 --- /dev/null +++ b/tests/HeavySymmetricTopPhysics.test.ts @@ -0,0 +1,389 @@ +/** + * HeavySymmetricTopPhysics.test.ts + * + * The integrator is only trustworthy if it reproduces the analytic results for the + * heavy symmetric top, so these tests check it against conservation laws, the + * steady-precession fixed point, the turning points of the effective potential, + * and the fast-top limits. + */ + +import { describe, expect, it } from "vitest"; +import { + createReleaseState, + criticalSpinRate, + gravityTorqueCoefficient, + type HeavyTopParameters, + type HeavyTopState, + isSleepingTopStable, + nutationFrequency, + nutationTurningPoints, + type ReleaseMode, + slowPrecessionRate, + spinAngularMomentum, + steadyPrecessionRates, + stepHeavyTop, + totalEnergy, + verticalAngularMomentum, +} from "../src/common/rigid-body/HeavySymmetricTopPhysics.js"; +import { + DEFAULT_NUTATION_SPIN_RAD_S, + DEFAULT_NUTATION_TILT_RAD, + GRAVITY_MPS2, + NUTATION_COM_DISTANCE_M, + NUTATION_SPIN_DRAG_N_M_S, + NUTATION_SPIN_INERTIA_KG_M2, + NUTATION_TIP_DRAG_N_M_S, + NUTATION_TRANSVERSE_INERTIA_KG_M2, + NUTATION_WHEEL_MASS_KG, +} from "../src/RigidBodyPrecessionConstants.js"; + +function baseParameters(overrides: Partial = {}): HeavyTopParameters { + return { + transverseInertia: NUTATION_TRANSVERSE_INERTIA_KG_M2, + spinInertia: NUTATION_SPIN_INERTIA_KG_M2, + mass: NUTATION_WHEEL_MASS_KG, + gravity: GRAVITY_MPS2, + comDistance: NUTATION_COM_DISTANCE_M, + tipDrag: 0, + spinDrag: 0, + ...overrides, + }; +} + +/** Integrate for `duration` seconds at 60 Hz, collecting the state at each frame. */ +function simulate( + parameters: HeavyTopParameters, + initial: HeavyTopState, + duration: number, + dt = 1 / 60, +): HeavyTopState[] { + const states: HeavyTopState[] = [initial]; + let state = initial; + const frames = Math.round(duration / dt); + for (let i = 0; i < frames; i++) { + state = stepHeavyTop(parameters, state, dt); + states.push(state); + } + return states; +} + +function release(mode: ReleaseMode, parameters = baseParameters(), spin = DEFAULT_NUTATION_SPIN_RAD_S): HeavyTopState { + return createReleaseState(parameters, spin, DEFAULT_NUTATION_TILT_RAD, mode); +} + +describe("HeavySymmetricTopPhysics", () => { + describe("conservation laws", () => { + it("conserves energy and both angular momenta without friction", () => { + const parameters = baseParameters(); + const initial = release("cusp"); + const states = simulate(parameters, initial, 10); + + const energy0 = totalEnergy(parameters, initial); + const vertical0 = verticalAngularMomentum(parameters, initial); + const spin0 = spinAngularMomentum(parameters, initial); + + for (const state of states) { + expect(totalEnergy(parameters, state)).toBeCloseTo(energy0, 6); + expect(verticalAngularMomentum(parameters, state)).toBeCloseTo(vertical0, 6); + expect(spinAngularMomentum(parameters, state)).toBeCloseTo(spin0, 8); + } + }); + + it("bleeds energy and spin away when friction is enabled", () => { + const parameters = baseParameters({ + tipDrag: NUTATION_TIP_DRAG_N_M_S, + spinDrag: NUTATION_SPIN_DRAG_N_M_S, + }); + const initial = release("cusp", parameters); + const states = simulate(parameters, initial, 8); + const final = states[states.length - 1]; + if (!final) { + throw new Error("no final state"); + } + + expect(totalEnergy(parameters, final)).toBeLessThan(totalEnergy(parameters, initial)); + + // ṗ_ψ = −c_s ω₃ gives exponential spin-down with time constant I₃/c_s. + const timeConstant = NUTATION_SPIN_INERTIA_KG_M2 / NUTATION_SPIN_DRAG_N_M_S; + expect(final.spin).toBeCloseTo(initial.spin * Math.exp(-8 / timeConstant), 4); + }); + }); + + describe("steady precession", () => { + it("holds the tilt exactly when released at the steady rate", () => { + const parameters = baseParameters(); + const initial = release("steady"); + const states = simulate(parameters, initial, 10); + + for (const state of states) { + expect(state.theta).toBeCloseTo(DEFAULT_NUTATION_TILT_RAD, 6); + expect(state.phiDot).toBeCloseTo(initial.phiDot, 6); + } + }); + + it("precesses through the steady rate over the simulated interval", () => { + const parameters = baseParameters(); + const initial = release("steady"); + const states = simulate(parameters, initial, 6); + const final = states[states.length - 1]; + if (!final) { + throw new Error("no final state"); + } + + const expected = slowPrecessionRate(parameters, DEFAULT_NUTATION_SPIN_RAD_S, DEFAULT_NUTATION_TILT_RAD); + expect((final.phi - initial.phi) / 6).toBeCloseTo(expected, 5); + }); + + it("has no steady solution below the critical spin", () => { + const parameters = baseParameters(); + const critical = criticalSpinRate(parameters, DEFAULT_NUTATION_TILT_RAD); + expect(critical).toBeGreaterThan(0); + + expect(steadyPrecessionRates(parameters, critical * 0.95, DEFAULT_NUTATION_TILT_RAD).exists).toBe(false); + expect(steadyPrecessionRates(parameters, critical * 1.05, DEFAULT_NUTATION_TILT_RAD).exists).toBe(true); + }); + + it("collapses to the gyroscopic rate Ω = Mgl/(I₃ω₃) for a fast top", () => { + const parameters = baseParameters(); + const spin = 400; + const roots = steadyPrecessionRates(parameters, spin, DEFAULT_NUTATION_TILT_RAD); + const gyroscopic = gravityTorqueCoefficient(parameters) / (NUTATION_SPIN_INERTIA_KG_M2 * spin); + expect(roots.slow / gyroscopic).toBeCloseTo(1, 2); + expect(roots.fast).toBeGreaterThan(roots.slow * 100); + }); + + it("requires a faster spin to precess steadily nearer the vertical", () => { + const parameters = baseParameters(); + expect(criticalSpinRate(parameters, Math.PI / 6)).toBeGreaterThan(criticalSpinRate(parameters, Math.PI / 3)); + expect(criticalSpinRate(parameters, Math.PI / 2)).toBeCloseTo(0, 6); + }); + + it("marks a vertical top as sleeping only above the stability threshold", () => { + const parameters = baseParameters(); + const threshold = criticalSpinRate(parameters, 0); + expect(isSleepingTopStable(parameters, threshold * 1.01)).toBe(true); + expect(isSleepingTopStable(parameters, threshold * 0.99)).toBe(false); + }); + }); + + describe("nutation", () => { + it("dips below the release tilt and returns to it when released from rest", () => { + const parameters = baseParameters(); + const initial = release("cusp"); + const states = simulate(parameters, initial, 6); + const thetas = states.map((state) => state.theta); + + expect(Math.min(...thetas)).toBeCloseTo(DEFAULT_NUTATION_TILT_RAD, 3); + expect(Math.max(...thetas)).toBeGreaterThan(DEFAULT_NUTATION_TILT_RAD + 0.05); + }); + + it("never precesses backwards when released from rest, and pauses at each cusp", () => { + const parameters = baseParameters(); + const states = simulate(parameters, release("cusp"), 6); + const phiDots = states.map((state) => state.phiDot); + + expect(Math.min(...phiDots)).toBeGreaterThan(-1e-6); + expect(Math.min(...phiDots)).toBeLessThan(1e-3); + expect(Math.max(...phiDots)).toBeGreaterThan(0); + }); + + it("precesses backwards for part of each cycle when pushed backwards", () => { + const parameters = baseParameters(); + const states = simulate(parameters, release("loop"), 6); + const phiDots = states.map((state) => state.phiDot); + + expect(Math.min(...phiDots)).toBeLessThan(0); + expect(Math.max(...phiDots)).toBeGreaterThan(0); + }); + + it("keeps precessing forwards, without cusps, when pushed forwards", () => { + const parameters = baseParameters(); + const states = simulate(parameters, release("smooth"), 6); + const phiDots = states.map((state) => state.phiDot); + + expect(Math.min(...phiDots)).toBeGreaterThan(0); + }); + + it("stays inside the turning points predicted by the effective potential", () => { + const parameters = baseParameters(); + const initial = release("cusp"); + const band = nutationTurningPoints(parameters, initial); + const states = simulate(parameters, initial, 12); + + expect(band.thetaMin).toBeCloseTo(DEFAULT_NUTATION_TILT_RAD, 4); + expect(band.thetaMax).toBeGreaterThan(band.thetaMin); + + for (const state of states) { + expect(state.theta).toBeGreaterThanOrEqual(band.thetaMin - 1e-4); + expect(state.theta).toBeLessThanOrEqual(band.thetaMax + 1e-4); + } + + // The motion must actually reach both turning points, not merely stay inside them. + const thetas = states.map((state) => state.theta); + expect(Math.max(...thetas)).toBeCloseTo(band.thetaMax, 3); + }); + + it("reports the same band from any point along the trajectory", () => { + const parameters = baseParameters(); + const initial = release("loop"); + const band = nutationTurningPoints(parameters, initial); + const states = simulate(parameters, initial, 5); + + for (const state of states) { + const current = nutationTurningPoints(parameters, state); + expect(current.thetaMin).toBeCloseTo(band.thetaMin, 4); + expect(current.thetaMax).toBeCloseTo(band.thetaMax, 4); + } + }); + + it("nutates at I₃ω₃/I₁ in the fast-top limit", () => { + const parameters = baseParameters(); + const spin = 120; + const initial = release("cusp", parameters, spin); + const dt = 1 / 2000; + const states = simulate(parameters, initial, 2, dt); + + // Time between successive cusps (θ̇ crossing zero from below) is the nutation period. + const cuspTimes: number[] = []; + for (let i = 1; i < states.length; i++) { + const previous = states[i - 1]; + const current = states[i]; + if (previous && current && previous.thetaDot > 0 && current.thetaDot <= 0) { + cuspTimes.push(i * dt); + } + } + expect(cuspTimes.length).toBeGreaterThan(2); + + const first = cuspTimes[0] ?? 0; + const last = cuspTimes[cuspTimes.length - 1] ?? 0; + const measuredPeriod = (last - first) / (cuspTimes.length - 1); + const predictedPeriod = (2 * Math.PI) / nutationFrequency(parameters, spin); + expect(measuredPeriod / predictedPeriod).toBeCloseTo(1, 1); + }); + + it("nutates further when the spin is lower", () => { + const parameters = baseParameters(); + const fast = nutationTurningPoints(parameters, release("cusp", parameters, 12)); + const slow = nutationTurningPoints(parameters, release("cusp", parameters, 6)); + expect(slow.thetaMax - slow.thetaMin).toBeGreaterThan(fast.thetaMax - fast.thetaMin); + }); + }); + + describe("friction", () => { + it("damps the nutation band down toward steady precession", () => { + const parameters = baseParameters({ + tipDrag: NUTATION_TIP_DRAG_N_M_S, + spinDrag: 0, + }); + const initial = release("cusp", parameters); + const initialBand = nutationTurningPoints(parameters, initial); + const states = simulate(parameters, initial, 20); + const final = states[states.length - 1]; + if (!final) { + throw new Error("no final state"); + } + + const finalBand = nutationTurningPoints(parameters, final); + expect(finalBand.thetaMax - finalBand.thetaMin).toBeLessThan( + 0.25 * (initialBand.thetaMax - initialBand.thetaMin), + ); + }); + + it("drops the axis once the spin decays past the critical value", () => { + const parameters = baseParameters({ + tipDrag: NUTATION_TIP_DRAG_N_M_S, + spinDrag: NUTATION_SPIN_DRAG_N_M_S, + }); + const initial = release("steady", parameters); + const states = simulate(parameters, initial, 60); + const final = states[states.length - 1]; + if (!final) { + throw new Error("no final state"); + } + + expect(final.spin).toBeLessThan(criticalSpinRate(parameters, DEFAULT_NUTATION_TILT_RAD)); + expect(final.theta).toBeGreaterThan(DEFAULT_NUTATION_TILT_RAD); + }); + }); + + describe("mechanical stop", () => { + it("rests the axle at maxTilt instead of swinging through it", () => { + const parameters = baseParameters({ + tipDrag: NUTATION_TIP_DRAG_N_M_S, + spinDrag: NUTATION_SPIN_DRAG_N_M_S, + maxTilt: Math.PI / 2, + }); + const states = simulate(parameters, release("steady", parameters), 90); + + for (const state of states) { + expect(state.theta).toBeLessThanOrEqual(Math.PI / 2 + 1e-9); + } + const final = states[states.length - 1]; + if (!final) { + throw new Error("no final state"); + } + expect(final.theta).toBeCloseTo(Math.PI / 2, 6); + }); + + it("reports a band that stops at the mount", () => { + const parameters = baseParameters({ maxTilt: Math.PI / 3 }); + const band = nutationTurningPoints(parameters, release("cusp", parameters, 4)); + expect(band.thetaMax).toBeCloseTo(Math.PI / 3, 9); + + // The same release without a stop nutates well past it. + const unbounded = nutationTurningPoints(baseParameters(), release("cusp", baseParameters(), 4)); + expect(unbounded.thetaMax).toBeGreaterThan(Math.PI / 3); + }); + + it("is a stop, not a trap — an axis already rising leaves it", () => { + const parameters = baseParameters({ maxTilt: Math.PI / 2 }); + // Gyroscopically supported at the stop, so gravity does not immediately slam it back. + const spin = 30; + const resting: HeavyTopState = { + theta: Math.PI / 2, + thetaDot: -0.5, + phi: 0, + phiDot: slowPrecessionRate(parameters, spin, Math.PI / 2), + psi: 0, + spin, + }; + // The contact zeroes a downward θ̇ but must leave an upward one alone, so the + // axis lifts off and nutates rather than being pinned at the stop forever. + const thetas = simulate(parameters, resting, 0.5, 1 / 2000).map((state) => state.theta); + expect(Math.min(...thetas)).toBeLessThan(Math.PI / 2 - 0.005); + }); + }); + + describe("stepHeavyTop", () => { + it("returns the same state for a non-positive timestep", () => { + const parameters = baseParameters(); + const initial = release("cusp"); + expect(stepHeavyTop(parameters, initial, 0)).toBe(initial); + expect(stepHeavyTop(parameters, initial, -0.1)).toBe(initial); + }); + + it("agrees with itself across timestep sizes", () => { + const parameters = baseParameters(); + const initial = release("cusp"); + const coarse = simulate(parameters, initial, 4, 1 / 30); + const fine = simulate(parameters, initial, 4, 1 / 240); + const coarseFinal = coarse[coarse.length - 1]; + const fineFinal = fine[fine.length - 1]; + if (!(coarseFinal && fineFinal)) { + throw new Error("no final state"); + } + + expect(coarseFinal.theta).toBeCloseTo(fineFinal.theta, 5); + expect(coarseFinal.phi).toBeCloseTo(fineFinal.phi, 5); + }); + + it("wraps the spin angle into [0, 2π)", () => { + const parameters = baseParameters(); + const states = simulate(parameters, release("cusp"), 5); + for (const state of states) { + expect(state.psi).toBeGreaterThanOrEqual(0); + expect(state.psi).toBeLessThan(2 * Math.PI); + } + }); + }); +}); diff --git a/tests/NutationModel.test.ts b/tests/NutationModel.test.ts new file mode 100644 index 0000000..dde6061 --- /dev/null +++ b/tests/NutationModel.test.ts @@ -0,0 +1,142 @@ +/** + * NutationModel.test.ts + */ + +import { describe, expect, it } from "vitest"; +import { criticalSpinRate, slowPrecessionRate } from "../src/common/rigid-body/HeavySymmetricTopPhysics.js"; +import { NutationModel } from "../src/nutation-screen/model/NutationModel.js"; +import { + DEFAULT_NUTATION_SPIN_RAD_S, + DEFAULT_NUTATION_TILT_RAD, + NUTATION_TRACE_CAPACITY, +} from "../src/RigidBodyPrecessionConstants.js"; + +function run(model: NutationModel, seconds: number): void { + const frames = Math.round(seconds * 60); + for (let i = 0; i < frames; i++) { + model.step(1 / 60); + } +} + +describe("NutationModel", () => { + it("starts released at the initial tilt with the clock running", () => { + const model = new NutationModel(); + expect(model.thetaProperty.value).toBeCloseTo(DEFAULT_NUTATION_TILT_RAD, 6); + expect(model.thetaDotProperty.value).toBe(0); + expect(model.spinProperty.value).toBe(DEFAULT_NUTATION_SPIN_RAD_S); + expect(model.timer.isPlayingProperty.value).toBe(true); + model.dispose(); + }); + + it("nutates within the reported band when released from rest", () => { + const model = new NutationModel(); + const band = model.nutationBandProperty.value; + let lowest = model.thetaProperty.value; + let highest = model.thetaProperty.value; + + for (let i = 0; i < 600; i++) { + model.step(1 / 60); + lowest = Math.min(lowest, model.thetaProperty.value); + highest = Math.max(highest, model.thetaProperty.value); + } + + expect(highest - lowest).toBeGreaterThan(0.05); + expect(lowest).toBeGreaterThanOrEqual(band.thetaMin - 1e-3); + expect(highest).toBeLessThanOrEqual(band.thetaMax + 1e-3); + model.dispose(); + }); + + it("measures a mean precession rate close to the steady rate for the cusp release", () => { + const model = new NutationModel(); + run(model, 6); + + // Averaged over whole nutation cycles the drift matches steady precession closely. + const steady = slowPrecessionRate(model.getParameters(), DEFAULT_NUTATION_SPIN_RAD_S, DEFAULT_NUTATION_TILT_RAD); + expect(model.meanPrecessionRateProperty.value).toBeGreaterThan(0); + expect(model.meanPrecessionRateProperty.value / steady).toBeGreaterThan(0.5); + expect(model.meanPrecessionRateProperty.value / steady).toBeLessThan(1.6); + model.dispose(); + }); + + it("holds the tilt with no nutation in steady release mode", () => { + const model = new NutationModel(); + model.releaseModeProperty.value = "steady"; + run(model, 8); + + const band = model.nutationBandProperty.value; + expect(model.thetaProperty.value).toBeCloseTo(DEFAULT_NUTATION_TILT_RAD, 4); + expect(band.thetaMax - band.thetaMin).toBeLessThan(1e-3); + model.dispose(); + }); + + it("re-releases the top when a launch parameter changes", () => { + const model = new NutationModel(); + run(model, 3); + expect(model.timer.timeProperty.value).toBeGreaterThan(0); + expect(model.getTraceSamples().length).toBeGreaterThan(0); + + model.spinRateProperty.value = 10; + + expect(model.timer.timeProperty.value).toBe(0); + expect(model.getTraceSamples().length).toBe(0); + expect(model.thetaProperty.value).toBeCloseTo(DEFAULT_NUTATION_TILT_RAD, 6); + expect(model.spinProperty.value).toBe(10); + model.dispose(); + }); + + it("keeps the trace bounded by its capacity", () => { + const model = new NutationModel(); + run(model, 30); + expect(model.getTraceSamples().length).toBe(NUTATION_TRACE_CAPACITY); + expect(model.getThetaGraphPoints().length).toBe(NUTATION_TRACE_CAPACITY); + model.dispose(); + }); + + it("does not advance while paused, but steps forward on demand", () => { + const model = new NutationModel(); + model.timer.isPlayingProperty.value = false; + const theta = model.thetaProperty.value; + + run(model, 2); + expect(model.thetaProperty.value).toBe(theta); + expect(model.timer.timeProperty.value).toBe(0); + + for (let i = 0; i < 60; i++) { + model.stepOnce(1 / 60); + } + expect(model.thetaProperty.value).not.toBe(theta); + expect(model.timer.timeProperty.value).toBeCloseTo(1, 6); + model.dispose(); + }); + + it("spins the top down and drops it when friction is enabled", () => { + const model = new NutationModel(); + model.releaseModeProperty.value = "steady"; + model.frictionEnabledProperty.value = true; + run(model, 60); + + // Steady precession at the release tilt is no longer possible, so the axis falls. + expect(model.spinProperty.value).toBeLessThan(criticalSpinRate(model.getParameters(), DEFAULT_NUTATION_TILT_RAD)); + expect(model.thetaProperty.value).toBeGreaterThan(DEFAULT_NUTATION_TILT_RAD); + model.dispose(); + }); + + it("reset restores the launch parameters and clears the history", () => { + const model = new NutationModel(); + model.releaseModeProperty.value = "loop"; + model.frictionEnabledProperty.value = true; + model.spinRateProperty.value = 12; + run(model, 4); + + model.reset(); + + expect(model.spinRateProperty.value).toBe(DEFAULT_NUTATION_SPIN_RAD_S); + expect(model.initialTiltProperty.value).toBeCloseTo(DEFAULT_NUTATION_TILT_RAD, 12); + expect(model.releaseModeProperty.value).toBe("cusp"); + expect(model.frictionEnabledProperty.value).toBe(false); + expect(model.timer.timeProperty.value).toBe(0); + expect(model.getTraceSamples().length).toBe(0); + expect(model.thetaProperty.value).toBeCloseTo(DEFAULT_NUTATION_TILT_RAD, 6); + model.dispose(); + }); +});