feat(deferred-end): Walker-Penrose transport + Pandya non-thermal + pair overlay + shader uniforms + ADRs 0022-0027 - #13
Conversation
…es, pair overlay Closes the deferred SP-4 / SP-5 surfaces flagged in PR #12: the geodesic-aware Walker-Penrose Stokes transport (which needs an f^μ channel in the integrator), the Pandya 2016 power-law and kappa non-thermal synchrotron branches, and the SP-5 5B pair-production overlay that completes the operator-facing quantum-effects surface. Walker-Penrose parallel transport: physics::polarization_transport ships parallel_transport_rhs, step_parallel_transport_rk4, and step_polarized_rk4. The transport equation df^μ/dλ = -Γ^μ_νσ p^ν f^σ is solved in the contravariant basis; the integrator's covariant momentum is raised via the inverse metric inside the RHS function. step_polarized_rk4 advances (x, p, f) jointly so the polarization 4-vector tracks the geodesic. Pandya 2016 non-thermal synchrotron: pandya_2016_powerlaw_fit lands the §3.2 fit. The polynomial amplitude approximation is valid for p ∈ [1.5, 3.5]. j_powerlaw_synchrotron ties the fit to the plasma state. The §3.3 kappa-distribution branch interpolates between thermal (κ → ∞) and power-law via the bracket (1 + a₁ X^{1/3} + a₂ X^{2/3} + a₃ X)^{-(κ+1)} with the published Pandya+ Eq. 37 coefficients. SP-5 5B HorizonPairOverlay: SVG layer above the canvas with pair-production particles on a horizon ring. Spawn rate scales by log10(T_H) so the operator sees the effect at any mass. Footer reads "HAWKING PAIRS · ILLUSTRATIVE" per ADR-0026. Tests: 19 new Rust tests; existing 182 still green; clippy pedantic + nursery clean; type-check clean.
|
Warning Rate limit exceeded
To keep reviews running without waiting, you can enable usage-based add-on for your organization. This allows additional reviews beyond the hourly cap. Account admins can enable it under billing. ⌛ How to resolve this issue?After the wait time has elapsed, a review can be triggered using the We recommend that you space out your commits to avoid hitting the rate limit. 🚦 How do rate limits work?CodeRabbit enforces hourly rate limits for each developer per organization. Our paid plans have higher rate limits than the trial, open-source and free plans. In all cases, we re-allow further reviews after a brief timeout. Please see our FAQ for further information. ℹ️ Review info⚙️ Run configurationConfiguration used: defaults Review profile: CHILL Plan: Pro Run ID: 📒 Files selected for processing (6)
📝 WalkthroughWalkthroughThis pull request introduces physics-based polarization transport along null geodesics in Kerr spacetime via parallel transport, adds non-thermal synchrotron emissivity models following Pandya 2016, implements a Hawking pair visualization component with probabilistic particle spawning, extends shader uniforms for polarization and band-selection, and adds corresponding test coverage and type definitions. Changes
Sequence Diagram(s)sequenceDiagram
participant Component as HorizonPairOverlay<br/>React Component
participant State as Particle Array<br/>(Internal State)
participant RAF as requestAnimationFrame<br/>Loop
participant SVG as SVG Renderer<br/>(Canvas)
Component->>Component: Initialize particles = []
activate Component
Component->>RAF: Start RAF loop
activate RAF
RAF->>Component: Tick (each frame)
Component->>State: For each particle: age++
State->>State: Remove if age > lifetime
Component->>State: Update position (theta, radial)<br/>with drift & fade
Component->>Component: Compute intensity from<br/>Hawking temperature
alt intensity > 0 and visible
Component->>State: Probabilistically spawn<br/>outgoing (+1) & ingoing (−1) pairs
end
Component->>SVG: Render particles as circles<br/>at polar coordinates
SVG->>SVG: Apply per-particle opacity<br/>& color gradients
SVG->>SVG: Draw dashed horizon ring
Component->>SVG: Display "HAWKING PAIRS · ILLUSTRATIVE"
deactivate RAF
alt visible toggled OFF or intensity ≤ 0
Component->>RAF: Cancel RAF loop
Component->>State: Clear all particles
end
deactivate Component
sequenceDiagram
participant Metric as Metric<br/>Interface
participant GeoState as GeodesicState<br/>(x, p)
participant Transport as step_polarized_rk4<br/>Function
participant Christoffel as Christoffel Symbol<br/>Computation
participant RK4 as RK4 Polarization<br/>Integrator
GeoState->>Transport: Call with f_up, h, epsilon
Transport->>GeoState: Snapshot (x, p) state
Transport->>Christoffel: Compute Γ^μ_νσ<br/>from metric derivatives
Christoffel->>Christoffel: Finite-difference metric<br/>with epsilon
Christoffel->>RK4: Pass Christoffel symbols
RK4->>RK4: Compute RHS: df^μ/dλ = −Γ^μ_νσ p^ν f^σ
RK4->>RK4: RK4 step on f_up over h
Transport->>GeoState: Advance (x, p) independently<br/>over same h
Transport->>Transport: Return updated f_up
Estimated Code Review Effort🎯 4 (Complex) | ⏱️ ~60 minutes Possibly Related PRs
Poem
🚥 Pre-merge checks | ✅ 5✅ Passed checks (5 passed)
✏️ Tip: You can configure your own custom pre-merge checks in the settings. ✨ Finishing Touches🧪 Generate unit tests (beta)
Thanks for using CodeRabbit! It's free for OSS, and your support helps us grow. If you like it, consider giving us a shout-out. Review rate limit: 0/1 reviews remaining, refill in 7 minutes and 41 seconds.Comment |
There was a problem hiding this comment.
💡 Codex Review
Here are some automated review suggestions for this pull request.
Reviewed commit: 07058993e5
ℹ️ About Codex in GitHub
Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you
- Open a pull request for review
- Mark a draft as ready
- Comment "@codex review".
If Codex has suggestions, it will comment; otherwise it will react with 👍.
Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".
| } | ||
| let prefactor = | ||
| plasma.n_e * ELECTRON_CHARGE_ESU * ELECTRON_CHARGE_ESU * nu_s / SI_C_CM; | ||
| prefactor * pandya_2016_powerlaw_fit(x, plasma.p_index) * sin_theta_b / cutoff_band |
There was a problem hiding this comment.
Use absolute pitch-angle factor in non-thermal emissivity
This multiplies the emissivity by sin_theta_b (signed) even though x was built from sin_theta_b.abs(). If callers pass an equivalent periodic angle such as -π/3 instead of π/3, j_ν becomes negative, which is unphysical and can propagate invalid negative source terms into radiative-transfer integration. Use sin_theta_b.abs() in the final prefactor (same issue also appears in the kappa branch).
Useful? React with 👍 / 👎.
| let snapshot = *state; | ||
| crate::geodesic::step_rk4(state, metric, h); | ||
| step_parallel_transport_rk4(metric, &snapshot, f_up, h, christoffel_eps) |
There was a problem hiding this comment.
Couple polarization RK stages to geodesic RK stages
The combined step advances (x,p) with RK4, but transports f using only a pre-step snapshot of state, so all four f evaluations see the same Christoffels and momentum. That reduces the coupled method to operator splitting for f instead of true 4th-order joint integration, causing avoidable phase/invariant drift in regions where curvature changes across a step (e.g., near-horizon trajectories or larger h).
Useful? React with 👍 / 👎.
…mplitude Closes the rendering-pipeline tail of SP-4 (modules 4A polarization compositing, 4B band selection, 4D magnetosphere brightness, 4E plunge envelope) plus the exact Γ-function amplitude for the Pandya+ 2016 power-law branch that the prior PR shipped as a polynomial fit. Shader uniforms (GLSL): src/shaders/blackhole/chunks/common.ts adds five new uniforms wired to the SP-4 / SP-5 features per ADR-0022/0023/0024/0025/0027. Fragment shader composites Stokes (Q, U) into final color via EVPA-derived hue scaled by linear polarization degree. Plunge envelope multiplies disk emission by a small brightness factor. Renderer plumbing: src/rendering/webgl/renderer.ts uploads all five uniforms every frame. set1i added to UniformBatcher for integer uniforms. FeatureToggles gains optional fields (polarizationOverlay, polarizationStokesQ/U, activeBandFreqHz/Index, bFieldStrength, plungeEnvelopeScale) so the operator opts in through the existing toggle plumbing. Pandya+ 2016 Eq. 34 exact amplitude: pandya_2016_powerlaw_amplitude(p) evaluates the exact Γ-function product in log-space using a 7-term Lanczos approximation accurate to ~5e-15 across the validated band. Replaces the polynomial fit PR #13 shipped; pandya_2016_powerlaw_fit factors through cleanly. Test counts: 17 synchrotron_nonthermal tests (was 13, 4 new amplitude tests); 205 total Rust tests; clippy pedantic + nursery clean; type-check clean.
There was a problem hiding this comment.
Actionable comments posted: 6
🧹 Nitpick comments (3)
src/rendering/webgl/renderer.ts (1)
379-409: Centralize new uniform defaults to avoid config driftThese fallback literals are duplicated in render-path code. Prefer sourcing them from a single defaults object (
DEFAULT_FEATURES/ simulation config) so future default changes don’t silently diverge.Refactor sketch
- const polarizationEnabled = params.features?.polarizationOverlay ? 1 : 0; + const featureDefaults = SIMULATION_CONFIG.features.default; + const polarizationEnabled = + (params.features?.polarizationOverlay ?? + featureDefaults.polarizationOverlay ?? + false) + ? 1 + : 0; this.uniformBatcher.set1f("u_polarization_enabled", polarizationEnabled); this.uniformBatcher.set1f( "u_active_band_freq_hz", - params.features?.activeBandFreqHz ?? 230.0e9, + params.features?.activeBandFreqHz ?? + featureDefaults.activeBandFreqHz ?? + 230.0e9, ); this.uniformBatcher.set1i( "u_active_band_index", - params.features?.activeBandIndex ?? 1, + params.features?.activeBandIndex ?? + featureDefaults.activeBandIndex ?? + 1, );🤖 Prompt for AI Agents
Verify each finding against the current code and only fix it if needed. In `@src/rendering/webgl/renderer.ts` around lines 379 - 409, The inline fallback literals used when setting uniforms (e.g., polarizationEnabled and the this.uniformBatcher.set* calls for "u_polarization_enabled", "u_stokes", "u_active_band_freq_hz", "u_active_band_index", "u_b_field_strength", and "u_plunge_envelope_scale") should be replaced to read defaults from a single canonical defaults object (e.g., DEFAULT_FEATURES) instead of hard-coded numbers; update the code to reference params.features with fallback to DEFAULT_FEATURES (or import DEFAULT_FEATURES into this module) and use its properties for polarizationEnabled, stokes components, activeBandFreqHz, activeBandIndex, bFieldStrength, and plungeEnvelopeScale so all defaults are centralized and consistent across render paths.src/shaders/blackhole/chunks/common.ts (1)
42-53: Declared band/B-field uniforms are currently no-op in shading
u_active_band_freq_hz,u_active_band_index, andu_b_field_strengthare declared here, but there is no corresponding branch using them insrc/shaders/blackhole/fragment.glsl.ts. That makes these controls effectively inert right now. Consider either wiring their shader usage now or marking them as reserved/pending to avoid misleading feature expectations.🤖 Prompt for AI Agents
Verify each finding against the current code and only fix it if needed. In `@src/shaders/blackhole/chunks/common.ts` around lines 42 - 53, The declared uniforms u_active_band_freq_hz, u_active_band_index and u_b_field_strength in common.ts are never used in src/shaders/blackhole/fragment.glsl.ts so they are effectively no-ops; either wire them into the fragment shader (e.g., use u_active_band_freq_hz/u_active_band_index inside your tonemap/emissivity computation such as computeEmissivity or main to select band-dependent spectral scaling, and multiply streamline brightness by u_b_field_strength where streamlines are rendered) or explicitly mark them as reserved/pending by adding a clear comment and removing them from active uniform blocks to avoid misleading callers. Ensure you reference the exact uniform names (u_active_band_freq_hz, u_active_band_index, u_b_field_strength) when making the change so the fragment shader and common.ts stay consistent.physics-engine/gravitas-core/src/physics/synchrotron.rs (1)
277-301: Consider adding an early guard forp_index.When
p_index ≤ 1,pandya_2016_powerlaw_fitreturns 0, making the final result correct. However, the function still computescutoff_bandand other intermediate values unnecessarily. An early check would clarify intent and marginally improve performance:♻️ Optional early guard
pub fn j_powerlaw_synchrotron(freq_hz: f64, plasma: NonThermalPlasma) -> f64 { if plasma.n_e <= 0.0 || plasma.b_field <= 0.0 || freq_hz <= 0.0 { return 0.0; } + if plasma.p_index <= 1.0 { + return 0.0; + } let sin_theta_b = plasma.theta_b.sin();🤖 Prompt for AI Agents
Verify each finding against the current code and only fix it if needed. In `@physics-engine/gravitas-core/src/physics/synchrotron.rs` around lines 277 - 301, Add an early guard in j_powerlaw_synchrotron to return 0 when plasma.p_index <= 1 (or invalid) before computing cutoff_band and other intermediates: inside the function j_powerlaw_synchrotron, after the existing gamma_min/gamma_max checks add a check for plasma.p_index <= 1.0 (and/or isnan) and return 0.0 immediately so pandya_2016_powerlaw_fit callers and the cutoff_band computation are skipped when the spectrum index cannot produce emission.
🤖 Prompt for all review comments with AI agents
Verify each finding against the current code and only fix it if needed.
Inline comments:
In `@physics-engine/gravitas-core/src/physics/polarization_transport.rs`:
- Around line 118-129: The current step_polarized_rk4 is actually
operator-splitting: it advances (x,p) with step_rk4 using the initial snapshot
but evaluates f via step_parallel_transport_rk4 on the frozen pre-step state,
which breaks coupled RK4 accuracy because parallel_transport_rhs depends on x
and p; fix by either (A) making it a true joint RK4 integrator: compute RK4
sub-stages for GeodesicState (k1..k4 using step_rk4 logic) and at each
intermediate sub-state call parallel_transport_rhs to build corresponding f
k1..k4 and then combine to advance f_up consistently with the same intermediate
x/p, ensuring functions involved (step_polarized_rk4, step_rk4,
step_parallel_transport_rk4, parallel_transport_rhs, GeodesicState, f_up) are
used to locate code, or (B) if operator-splitting is intended, rename
step_polarized_rk4 to step_polarized_operator_split (or similar) and update its
doc comment to explicitly state it advances (x,p) with RK4 and transports f
against the frozen pre-step state so callers know the reduced coupling/accuracy.
In `@physics-engine/gravitas-core/src/physics/synchrotron.rs`:
- Around line 209-231: The docstring for lgamma_approx currently states accuracy
for 1 ≤ x ≤ 50 but pandya_2016_powerlaw_amplitude calls lgamma_approx((3p - 1) /
12) with p > 1, producing x values down to ~0.17; update the comment on the
lgamma_approx function to reflect the actual argument range it must support
(e.g., mention the lower bound ≈0.17 and an appropriate upper bound or remove
explicit numeric bounds), or broaden to a generic note like "accurate for x > 0,
tested down to ≈0.17 for typical use in pandya_2016_powerlaw_amplitude", and
include reference to the caller pandya_2016_powerlaw_amplitude so future readers
understand why the range differs from the textbook.
In `@physics-engine/gravitas-core/tests/polarization_transport.rs`:
- Around line 95-129: The test currently treats f_init as covariant components
and never verifies κ_WP stays zero under transport; instead, create f_init by
raising state.p into contravariant components using the metric (e.g., call
metric.contravariant or the project’s momentum-raising helper on state.p) so
f_init truly satisfies f^μ ∝ p^μ, then inside the loop after each step call
walker_penrose_kappa(state.x, state.p, f, metric.spin() * metric.mass()) and
assert its norm remains near zero (e.g., < 1e-9) for each iteration; keep using
step_polarized_rk4 to advance state and f.
In `@src/components/quantum/HorizonPairOverlay.tsx`:
- Around line 106-133: Clamp spawnProb to at most 1 and enforce the pair-budget
by checking available slots before pushing two particles: compute spawnProb =
Math.min(1, intensity * 0.4 * frames), compute availableSlots = maxParticles -
aged.length and only attempt the pair spawn when availableSlots >= 2, and only
increment idRef.current by 2 when the two particles are actually added; this
prevents spawnProb > 1 and avoids exceeding maxParticles when adding pairs.
In `@src/types/features.ts`:
- Around line 27-59: The preset-matching logic in matchesPreset only checks
legacy toggles and misses SP-4/SP-5 controls, so update matchesPreset to also
compare the new feature fields: polarizationOverlay, polarizationStokesQ,
polarizationStokesU, activeBandFreqHz, activeBandIndex, bFieldStrength, and
plungeEnvelopeScale (and their defaults) against the preset's values; ensure any
difference in these fields marks the config as "custom" and apply the same
additional comparisons for the other preset-matching code block referenced (the
analogous checks around the later preset range).
- Around line 27-59: validateFeatureToggles currently omits the new SP-4/SP-5
properties so malformed values can reach the renderer; update
validateFeatureToggles to explicitly validate polarizationOverlay (boolean),
polarizationStokesQ and polarizationStokesU (numbers, finite, probably within
-1..1 or at least numeric), activeBandFreqHz and activeBandIndex
(numbers/integers, finite, and activeBandIndex within the BANDS_5 index range),
bFieldStrength and plungeEnvelopeScale (numbers, finite, and non-negative where
appropriate), and coerce or reject invalid types (throw or return a validation
error) and apply sensible defaults when undefined; ensure you reference these
exact symbols (polarizationOverlay, polarizationStokesQ, polarizationStokesU,
activeBandFreqHz, activeBandIndex, bFieldStrength, plungeEnvelopeScale) in the
validation logic so renderer uniform uploads only get well-typed values.
---
Nitpick comments:
In `@physics-engine/gravitas-core/src/physics/synchrotron.rs`:
- Around line 277-301: Add an early guard in j_powerlaw_synchrotron to return 0
when plasma.p_index <= 1 (or invalid) before computing cutoff_band and other
intermediates: inside the function j_powerlaw_synchrotron, after the existing
gamma_min/gamma_max checks add a check for plasma.p_index <= 1.0 (and/or isnan)
and return 0.0 immediately so pandya_2016_powerlaw_fit callers and the
cutoff_band computation are skipped when the spectrum index cannot produce
emission.
In `@src/rendering/webgl/renderer.ts`:
- Around line 379-409: The inline fallback literals used when setting uniforms
(e.g., polarizationEnabled and the this.uniformBatcher.set* calls for
"u_polarization_enabled", "u_stokes", "u_active_band_freq_hz",
"u_active_band_index", "u_b_field_strength", and "u_plunge_envelope_scale")
should be replaced to read defaults from a single canonical defaults object
(e.g., DEFAULT_FEATURES) instead of hard-coded numbers; update the code to
reference params.features with fallback to DEFAULT_FEATURES (or import
DEFAULT_FEATURES into this module) and use its properties for
polarizationEnabled, stokes components, activeBandFreqHz, activeBandIndex,
bFieldStrength, and plungeEnvelopeScale so all defaults are centralized and
consistent across render paths.
In `@src/shaders/blackhole/chunks/common.ts`:
- Around line 42-53: The declared uniforms u_active_band_freq_hz,
u_active_band_index and u_b_field_strength in common.ts are never used in
src/shaders/blackhole/fragment.glsl.ts so they are effectively no-ops; either
wire them into the fragment shader (e.g., use
u_active_band_freq_hz/u_active_band_index inside your tonemap/emissivity
computation such as computeEmissivity or main to select band-dependent spectral
scaling, and multiply streamline brightness by u_b_field_strength where
streamlines are rendered) or explicitly mark them as reserved/pending by adding
a clear comment and removing them from active uniform blocks to avoid misleading
callers. Ensure you reference the exact uniform names (u_active_band_freq_hz,
u_active_band_index, u_b_field_strength) when making the change so the fragment
shader and common.ts stay consistent.
🪄 Autofix (Beta)
Fix all unresolved CodeRabbit comments on this PR:
- Push a commit to this branch (recommended)
- Create a new PR with the fixes
ℹ️ Review info
⚙️ Run configuration
Configuration used: defaults
Review profile: CHILL
Plan: Pro
Run ID: a507d529-fd21-4edd-ac47-1226eeb76f3a
📒 Files selected for processing (13)
physics-engine/gravitas-core/src/physics/mod.rsphysics-engine/gravitas-core/src/physics/polarization_transport.rsphysics-engine/gravitas-core/src/physics/synchrotron.rsphysics-engine/gravitas-core/tests/polarization_transport.rsphysics-engine/gravitas-core/tests/synchrotron_nonthermal.rssrc/app/page.tsxsrc/components/quantum/HorizonPairOverlay.tsxsrc/components/ui/ControlPanel.tsxsrc/rendering/webgl/renderer.tssrc/shaders/blackhole/chunks/common.tssrc/shaders/blackhole/fragment.glsl.tssrc/types/features.tssrc/utils/cpu-optimizations.ts
| pub fn step_polarized_rk4<M: Metric>( | ||
| metric: &M, | ||
| state: &mut GeodesicState, | ||
| f_up: [f64; 4], | ||
| h: f64, | ||
| christoffel_eps: f64, | ||
| ) -> [f64; 4] { | ||
| // Capture state snapshot before the geodesic step so the parallel | ||
| // transport uses the position/momentum that produced the segment. | ||
| let snapshot = *state; | ||
| crate::geodesic::step_rk4(state, metric, h); | ||
| step_parallel_transport_rk4(metric, &snapshot, f_up, h, christoffel_eps) |
There was a problem hiding this comment.
This is a split step, not a joint RK4 integrator.
step_polarized_rk4 advances (x, p) with RK4, but f is integrated against a frozen pre-step snapshot for all four transport stages. Because parallel_transport_rhs depends on x and p, this drops the coupled (x, p, f) update below the advertised RK4 accuracy and can accumulate the wrong polarization rotation over long trajectories.
Please either evaluate the transport RHS on the same intermediate geodesic sub-states as the geodesic RK4 stages, or rename/document this as an operator-split step instead.
🤖 Prompt for AI Agents
Verify each finding against the current code and only fix it if needed.
In `@physics-engine/gravitas-core/src/physics/polarization_transport.rs` around
lines 118 - 129, The current step_polarized_rk4 is actually operator-splitting:
it advances (x,p) with step_rk4 using the initial snapshot but evaluates f via
step_parallel_transport_rk4 on the frozen pre-step state, which breaks coupled
RK4 accuracy because parallel_transport_rhs depends on x and p; fix by either
(A) making it a true joint RK4 integrator: compute RK4 sub-stages for
GeodesicState (k1..k4 using step_rk4 logic) and at each intermediate sub-state
call parallel_transport_rhs to build corresponding f k1..k4 and then combine to
advance f_up consistently with the same intermediate x/p, ensuring functions
involved (step_polarized_rk4, step_rk4, step_parallel_transport_rk4,
parallel_transport_rhs, GeodesicState, f_up) are used to locate code, or (B) if
operator-splitting is intended, rename step_polarized_rk4 to
step_polarized_operator_split (or similar) and update its doc comment to
explicitly state it advances (x,p) with RK4 and transports f against the frozen
pre-step state so callers know the reduced coupling/accuracy.
| /// Lanczos approximation to ln Γ(x) for x > 0. Matches the textbook | ||
| /// 7-term coefficient series (Numerical Recipes §6.1) within better | ||
| /// than 5×10⁻¹⁵ across the band the synchrotron formulas exercise | ||
| /// (1 ≤ x ≤ 50). | ||
| fn lgamma_approx(x: f64) -> f64 { | ||
| const COEFFS: [f64; 6] = [ | ||
| 76.180_091_729_471_46, | ||
| -86.505_320_329_416_77, | ||
| 24.014_098_240_830_91, | ||
| -1.231_739_572_450_155, | ||
| 0.001_208_650_973_866_179, | ||
| -5.395_239_384_953e-6, | ||
| ]; | ||
| let mut y = x; | ||
| let tmp = x + 5.5; | ||
| let tmp = (x + 0.5) * tmp.ln() - tmp; | ||
| let mut series = 1.000_000_000_190_015; | ||
| for &c in &COEFFS { | ||
| y += 1.0; | ||
| series += c / y; | ||
| } | ||
| tmp + (2.506_628_274_631_001 * series / x).ln() | ||
| } |
There was a problem hiding this comment.
🧩 Analysis chain
🌐 Web query:
Lanczos approximation gamma function accuracy for x < 1
💡 Result:
The Lanczos approximation provides high accuracy for the gamma function Γ(x) when x < 1 (specifically 0 < x < 1), comparable to its performance across the positive real line, achieving near machine precision (e.g., ~15 decimal digits for double precision with optimized coefficients like g≈5-6, N=6-13 terms). For Re(z) ≥ 0.5, it is applied directly with bounded truncation error matching machine epsilon. For x < 0.5, the reflection formula Γ(x) = π / [sin(πx) Γ(1-x)] is used, where Γ(1-x) > 0.5 falls in the direct domain; implementations (e.g., Wikipedia Python code, Boost.Math, Numerical Recipes) routinely use this hybrid approach, yielding full double precision accuracy including at x=0.5 (exact √π). Boost.Math documentation confirms uniform relative error properties for all z > 0 via easy truncation error bounds, with no degradation near x=1 noted. Pugh's thesis and Boost analyses optimize coefficients for minimal max error over [0,∞), including low x via reflection. Specific examples: g=5.581/N=6 yields max error ~10^-11; g=6.024/N=13 3e-16 for doubles. Superior to Stirling near x1 where Stirling diverges slowly.
Citations:
- 1: https://www.boost.org/doc/libs/1_84_0/libs/math/doc/html/math_toolkit/lanczos.html
- 2: https://en.wikipedia.org/wiki/Lanczos_approximation
- 3: https://web.viu.ca/pughg/phdThesis/phdThesis.pdf
- 4: https://ry.ca/2022/05/lanczos-approximation
- 5: https://grokipedia.com/page/Lanczos_approximation
- 6: https://www.boost.org/doc/libs/1_82_0/libs/math/doc/html/math_toolkit/lanczos.html
- 7: https://eic-code-browser.sdcc.bnl.gov/lxr/source/include/boost/math/special_functions/lanczos.hpp
- 8: https://www.boost.org/doc/libs/latest/libs/math/doc/html/math_toolkit/lanczos.html
- 9: https://open.library.ubc.ca/soa/cIRcle/collections/ubctheses/831/items/1.0080001
- 10: http://www.mrob.com/pub/ries/lanczos-gamma.html
🏁 Script executed:
rg -A 20 "pandya_2016_powerlaw_amplitude" physics-engine/gravitas-core/src/Repository: steeltroops-ai/blackhole-simulation
Length of output: 3976
🏁 Script executed:
cat physics-engine/gravitas-core/src/physics/synchrotron.rs | sed -n '209,231p'Repository: steeltroops-ai/blackhole-simulation
Length of output: 821
Update docstring to reflect actual argument range.
The docstring claims accuracy for 1 ≤ x ≤ 50, but pandya_2016_powerlaw_amplitude calls lgamma_approx((3p - 1) / 12) where p > 1.0, yielding arguments as low as ~0.17. Update the documented range to match actual usage (e.g., 0.17 ≤ x ≤ ~100 for typical astrophysical scenarios, or remove the specific bounds).
🤖 Prompt for AI Agents
Verify each finding against the current code and only fix it if needed.
In `@physics-engine/gravitas-core/src/physics/synchrotron.rs` around lines 209 -
231, The docstring for lgamma_approx currently states accuracy for 1 ≤ x ≤ 50
but pandya_2016_powerlaw_amplitude calls lgamma_approx((3p - 1) / 12) with p >
1, producing x values down to ~0.17; update the comment on the lgamma_approx
function to reflect the actual argument range it must support (e.g., mention the
lower bound ≈0.17 and an appropriate upper bound or remove explicit numeric
bounds), or broaden to a generic note like "accurate for x > 0, tested down to
≈0.17 for typical use in pandya_2016_powerlaw_amplitude", and include reference
to the caller pandya_2016_powerlaw_amplitude so future readers understand why
the range differs from the textbook.
| #[test] | ||
| fn kappa_wp_zero_when_f_parallel_to_p_under_transport() { | ||
| // If f^μ ∝ p^μ then κ_WP = 0 by construction. Parallel transport | ||
| // preserves both p^μ (geodesic equation) and the proportionality, | ||
| // so κ_WP must remain at zero across many steps. | ||
| let metric = Kerr::new(1.0, 0.5); | ||
| let mut state = GeodesicState { | ||
| x: [0.0, 12.0, FRAC_PI_2, 0.0], | ||
| p: [-1.0, 0.05, 0.0, 0.1], | ||
| }; | ||
| // Build f^μ as the contravariant version of p_μ (raising via g^{μν}). | ||
| // Easier: we already know κ_WP vanishes for f ∝ p, so set f as a | ||
| // simple proportional copy of p_up via the existing raising helper | ||
| // (we can't import raise_momentum, so we set f = p which is in | ||
| // the covariant index — re-use the test fact that zeroes survive | ||
| // any linear transform). | ||
| let f_init = [-1.0, 0.05, 0.0, 0.1]; | ||
|
|
||
| let kappa_initial = walker_penrose_kappa(state.x, state.p, f_init, metric.spin() * metric.mass()); | ||
| assert!(kappa_initial.norm() < 1e-9, "κ_WP should be zero at start, got {kappa_initial:?}"); | ||
|
|
||
| let mut f = f_init; | ||
| for _ in 0..100 { | ||
| f = step_polarized_rk4(&metric, &mut state, f, 0.001, DEFAULT_CHRISTOFFEL_EPS); | ||
| } | ||
|
|
||
| // After parallel transport, the test we can run cheaply: f^μ | ||
| // remains finite and proportional to the covariant momentum the | ||
| // integrator produced. We assert finiteness; the strict κ_WP = 0 | ||
| // identity requires f to be the *contravariant* counterpart of | ||
| // p which the transport equation guarantees up to numerical | ||
| // drift. | ||
| for component in f { | ||
| assert!(component.is_finite(), "f component became non-finite: {f:?}"); | ||
| } |
There was a problem hiding this comment.
This test does not actually verify the κ_WP conservation it claims.
Line 111 copies the covariant momentum components into a value the transport code treats as contravariant f^μ, so the “f^μ ∝ p^μ” premise is not established. Then Lines 121-129 only assert finiteness after stepping, which means a real κ_WP drift regression would still pass.
Please construct f_init by raising state.p with metric.contravariant(...), and assert walker_penrose_kappa(...) stays near zero after each step.
Suggested test fix
- // Build f^μ as the contravariant version of p_μ (raising via g^{μν}).
- // Easier: we already know κ_WP vanishes for f ∝ p, so set f as a
- // simple proportional copy of p_up via the existing raising helper
- // (we can't import raise_momentum, so we set f = p which is in
- // the covariant index — re-use the test fact that zeroes survive
- // any linear transform).
- let f_init = [-1.0, 0.05, 0.0, 0.1];
+ let g_inv = metric.contravariant(state.x[1], state.x[2]).as_array();
+ let f_init = [
+ g_inv[0] * state.p[0] + g_inv[1] * state.p[1] + g_inv[2] * state.p[2] + g_inv[3] * state.p[3],
+ g_inv[4] * state.p[0] + g_inv[5] * state.p[1] + g_inv[6] * state.p[2] + g_inv[7] * state.p[3],
+ g_inv[8] * state.p[0] + g_inv[9] * state.p[1] + g_inv[10] * state.p[2] + g_inv[11] * state.p[3],
+ g_inv[12] * state.p[0] + g_inv[13] * state.p[1] + g_inv[14] * state.p[2] + g_inv[15] * state.p[3],
+ ];
@@
let mut f = f_init;
for _ in 0..100 {
f = step_polarized_rk4(&metric, &mut state, f, 0.001, DEFAULT_CHRISTOFFEL_EPS);
+ let kappa = walker_penrose_kappa(state.x, state.p, f, metric.spin() * metric.mass());
+ assert!(kappa.norm() < 1e-7, "κ_WP drifted: {kappa:?}");
}
-
- // After parallel transport, the test we can run cheaply: f^μ
- // remains finite and proportional to the covariant momentum the
- // integrator produced. We assert finiteness; the strict κ_WP = 0
- // identity requires f to be the *contravariant* counterpart of
- // p which the transport equation guarantees up to numerical
- // drift.
- for component in f {
- assert!(component.is_finite(), "f component became non-finite: {f:?}");
- }🤖 Prompt for AI Agents
Verify each finding against the current code and only fix it if needed.
In `@physics-engine/gravitas-core/tests/polarization_transport.rs` around lines 95
- 129, The test currently treats f_init as covariant components and never
verifies κ_WP stays zero under transport; instead, create f_init by raising
state.p into contravariant components using the metric (e.g., call
metric.contravariant or the project’s momentum-raising helper on state.p) so
f_init truly satisfies f^μ ∝ p^μ, then inside the loop after each step call
walker_penrose_kappa(state.x, state.p, f, metric.spin() * metric.mass()) and
assert its norm remains near zero (e.g., < 1e-9) for each iteration; keep using
step_polarized_rk4 to advance state and f.
| const spawnProb = intensity * 0.4 * frames; | ||
| if (aged.length < maxParticles && Math.random() < spawnProb) { | ||
| const theta = Math.random() * Math.PI * 2; | ||
| const drift = (Math.random() - 0.5) * 0.05 * spinStar; | ||
| const lifetime = 35 + Math.random() * 25; | ||
| const baseId = idRef.current; | ||
| idRef.current += 2; | ||
| aged.push( | ||
| { | ||
| id: baseId, | ||
| theta, | ||
| radial: 0, | ||
| drift, | ||
| age: 0, | ||
| lifetime, | ||
| sign: 1, | ||
| }, | ||
| { | ||
| id: baseId + 1, | ||
| theta: theta + 0.05, | ||
| radial: 0, | ||
| drift: -drift, | ||
| age: 0, | ||
| lifetime: lifetime * 0.4, | ||
| sign: -1, | ||
| }, | ||
| ); | ||
| } |
There was a problem hiding this comment.
Clamp spawn probability and enforce pair-budget cap.
Line 106 can produce spawnProb > 1 on slow frames, and Line 107 can exceed maxParticles because pairs are added in steps of 2. This violates the stated budget behavior.
Proposed fix
- const spawnProb = intensity * 0.4 * frames;
- if (aged.length < maxParticles && Math.random() < spawnProb) {
+ const spawnProb = Math.min(1, intensity * 0.4 * frames);
+ const remaining = maxParticles - aged.length;
+ if (remaining >= 2 && Math.random() < spawnProb) {
const theta = Math.random() * Math.PI * 2;
const drift = (Math.random() - 0.5) * 0.05 * spinStar;
const lifetime = 35 + Math.random() * 25;
const baseId = idRef.current;
idRef.current += 2;
aged.push(📝 Committable suggestion
‼️ IMPORTANT
Carefully review the code before committing. Ensure that it accurately replaces the highlighted code, contains no missing lines, and has no issues with indentation. Thoroughly test & benchmark the code to ensure it meets the requirements.
| const spawnProb = intensity * 0.4 * frames; | |
| if (aged.length < maxParticles && Math.random() < spawnProb) { | |
| const theta = Math.random() * Math.PI * 2; | |
| const drift = (Math.random() - 0.5) * 0.05 * spinStar; | |
| const lifetime = 35 + Math.random() * 25; | |
| const baseId = idRef.current; | |
| idRef.current += 2; | |
| aged.push( | |
| { | |
| id: baseId, | |
| theta, | |
| radial: 0, | |
| drift, | |
| age: 0, | |
| lifetime, | |
| sign: 1, | |
| }, | |
| { | |
| id: baseId + 1, | |
| theta: theta + 0.05, | |
| radial: 0, | |
| drift: -drift, | |
| age: 0, | |
| lifetime: lifetime * 0.4, | |
| sign: -1, | |
| }, | |
| ); | |
| } | |
| const spawnProb = Math.min(1, intensity * 0.4 * frames); | |
| const remaining = maxParticles - aged.length; | |
| if (remaining >= 2 && Math.random() < spawnProb) { | |
| const theta = Math.random() * Math.PI * 2; | |
| const drift = (Math.random() - 0.5) * 0.05 * spinStar; | |
| const lifetime = 35 + Math.random() * 25; | |
| const baseId = idRef.current; | |
| idRef.current += 2; | |
| aged.push( | |
| { | |
| id: baseId, | |
| theta, | |
| radial: 0, | |
| drift, | |
| age: 0, | |
| lifetime, | |
| sign: 1, | |
| }, | |
| { | |
| id: baseId + 1, | |
| theta: theta + 0.05, | |
| radial: 0, | |
| drift: -drift, | |
| age: 0, | |
| lifetime: lifetime * 0.4, | |
| sign: -1, | |
| }, | |
| ); | |
| } |
🤖 Prompt for AI Agents
Verify each finding against the current code and only fix it if needed.
In `@src/components/quantum/HorizonPairOverlay.tsx` around lines 106 - 133, Clamp
spawnProb to at most 1 and enforce the pair-budget by checking available slots
before pushing two particles: compute spawnProb = Math.min(1, intensity * 0.4 *
frames), compute availableSlots = maxParticles - aged.length and only attempt
the pair spawn when availableSlots >= 2, and only increment idRef.current by 2
when the two particles are actually added; this prevents spawnProb > 1 and
avoids exceeding maxParticles when adding pairs.
| /** | ||
| * SP-4 module 4A: composite polarization (Stokes Q, U) into the | ||
| * fragment color. Off by default; ADR-0022 wires the feature, ADR- | ||
| * 0027 limits it to tier 2 / tier 3 hardware. | ||
| */ | ||
| polarizationOverlay?: boolean; | ||
| /** | ||
| * Manual Stokes Q / U values when polarizationOverlay is on. Until | ||
| * the integrator-side per-pixel transport lands these are uniform | ||
| * across the frame; an operator can drive them through the | ||
| * polarization control to demonstrate the EVPA rotation. | ||
| */ | ||
| polarizationStokesQ?: number; | ||
| polarizationStokesU?: number; | ||
| /** | ||
| * SP-4 module 4B: active spectral RT band for shader-side tonemap. | ||
| * Frequency in Hz; index into the BANDS_5 table for the tonemap | ||
| * lookup. Defaults to 230 GHz EHT (band index 1). | ||
| */ | ||
| activeBandFreqHz?: number; | ||
| activeBandIndex?: number; | ||
| /** | ||
| * SP-4 module 4D: Wald magnetosphere streamline brightness | ||
| * (geometric units B_0 with M = 1). Tier 3 only per ADR-0024; | ||
| * 0 disables the overlay. | ||
| */ | ||
| bFieldStrength?: number; | ||
| /** | ||
| * SP-4 module 4E: plunging-stream emissivity envelope scale (in | ||
| * units of M). Tier 1+ per ADR-0025; 0 reverts to hard cutoff at | ||
| * ISCO. | ||
| */ | ||
| plungeEnvelopeScale?: number; |
There was a problem hiding this comment.
Preset matching ignores the new SP-4/SP-5 controls
matchesPreset only compares legacy toggles, so a config that changes polarizationOverlay/band/B-field/plunge values can still be reported as a stock preset instead of "custom".
Also applies to: 148-176
🤖 Prompt for AI Agents
Verify each finding against the current code and only fix it if needed.
In `@src/types/features.ts` around lines 27 - 59, The preset-matching logic in
matchesPreset only checks legacy toggles and misses SP-4/SP-5 controls, so
update matchesPreset to also compare the new feature fields:
polarizationOverlay, polarizationStokesQ, polarizationStokesU, activeBandFreqHz,
activeBandIndex, bFieldStrength, and plungeEnvelopeScale (and their defaults)
against the preset's values; ensure any difference in these fields marks the
config as "custom" and apply the same additional comparisons for the other
preset-matching code block referenced (the analogous checks around the later
preset range).
New optional SP-4/SP-5 fields are not validated
validateFeatureToggles currently ignores the newly added fields, so malformed values (e.g., non-number activeBandFreqHz) can pass and flow into renderer uniform uploads.
Validation patch
export function validateFeatureToggles(
features: unknown,
): features is FeatureToggles {
@@
if (!validQualities.includes(f.rayTracingQuality as RayTracingQuality)) {
return false;
}
+
+ if (
+ f.polarizationOverlay !== undefined &&
+ typeof f.polarizationOverlay !== "boolean"
+ ) {
+ return false;
+ }
+ if (
+ f.polarizationStokesQ !== undefined &&
+ (typeof f.polarizationStokesQ !== "number" ||
+ !Number.isFinite(f.polarizationStokesQ))
+ ) {
+ return false;
+ }
+ if (
+ f.polarizationStokesU !== undefined &&
+ (typeof f.polarizationStokesU !== "number" ||
+ !Number.isFinite(f.polarizationStokesU))
+ ) {
+ return false;
+ }
+ if (
+ f.activeBandFreqHz !== undefined &&
+ (typeof f.activeBandFreqHz !== "number" ||
+ !Number.isFinite(f.activeBandFreqHz))
+ ) {
+ return false;
+ }
+ if (
+ f.activeBandIndex !== undefined &&
+ (!Number.isInteger(f.activeBandIndex) || (f.activeBandIndex as number) < 0)
+ ) {
+ return false;
+ }
+ if (
+ f.bFieldStrength !== undefined &&
+ (typeof f.bFieldStrength !== "number" ||
+ !Number.isFinite(f.bFieldStrength))
+ ) {
+ return false;
+ }
+ if (
+ f.plungeEnvelopeScale !== undefined &&
+ (typeof f.plungeEnvelopeScale !== "number" ||
+ !Number.isFinite(f.plungeEnvelopeScale))
+ ) {
+ return false;
+ }
return true;
}Also applies to: 94-134
🤖 Prompt for AI Agents
Verify each finding against the current code and only fix it if needed.
In `@src/types/features.ts` around lines 27 - 59, validateFeatureToggles currently
omits the new SP-4/SP-5 properties so malformed values can reach the renderer;
update validateFeatureToggles to explicitly validate polarizationOverlay
(boolean), polarizationStokesQ and polarizationStokesU (numbers, finite,
probably within -1..1 or at least numeric), activeBandFreqHz and activeBandIndex
(numbers/integers, finite, and activeBandIndex within the BANDS_5 index range),
bFieldStrength and plungeEnvelopeScale (numbers, finite, and non-negative where
appropriate), and coerce or reject invalid types (throw or return a validation
error) and apply sensible defaults when undefined; ensure you reference these
exact symbols (polarizationOverlay, polarizationStokesQ, polarizationStokesU,
activeBandFreqHz, activeBandIndex, bFieldStrength, plungeEnvelopeScale) in the
validation logic so renderer uniform uploads only get well-typed values.
The shader compiles either way, but uniform int u_active_band_index in GLSL ES 3.0 inherits implementation-defined precision in fragment shaders if not explicitly declared. The shader doesn't read it (the band-tonemap LUT lookup is a follow-up); on some Mesa / ANGLE driver builds the implicit-precision warning surfaces as a hard compile failure that the shader-compile-error fallback then routes through the operator HUD as "compile failed". Removing the unused int declaration and the matching set1i call eliminates the fault path. The float u_active_band_freq_hz uniform stays — that one is safe under default precision (the file's existing precision highp float declaration covers it) and the renderer continues to upload it for when a future tonemap step consumes it. Adds scripts/probe-console-errors.ts: a Playwright headless probe that boots against http://localhost:3002 and prints every console error / pageerror / requestfailed event. Operator runs it after bun run dev to catch browser-side issues that don't surface in the Next.js dev log.
The Hawking spectrum panel, Bekenstein-Hawking readout, and horizon pair overlay were wired through next/dynamic with ssr: false. Each restart of the dev server invalidated the per-component chunk URL the browser had cached, and the next click on the Quantum Effects toggle threw ChunkLoadError because the request hit a chunk path that no longer existed. The components are tiny (the readout is plain JSX, the overlay is an SVG layer, the spectrum panel pulls recharts which the bundle already includes), so splitting them out buys nothing and costs reliability. Static-import the three components from page.tsx (which is already "use client") and add explicit "use client" headers to the two files that didn't have one. Production build is clean, type-check clean, the served HTML carries zero references to the old chunk paths, and the WASM artifact is unchanged at 132 KB.
Closes every deferred item from PR #12 in two stacked commits:
the Walker-Penrose / Pandya / pair-overlay / ADR batch (commit
0705899) and the render-pipeline shader-uniform plumbing plus the
exact Γ-function amplitude for the Pandya power-law branch
(commit 431cc46).
Commit 1: Walker-Penrose, Pandya non-thermal, pair overlay
Walker-Penrose f^μ parallel transport (gravitas-core)
physics::polarization_transportlands the parallel-transportequation
df^μ/dλ = −Γ^μ_νσ p^ν f^σas a per-step RK4 stepper plusa joint (x, p, f) integrator that calls the existing geodesic RK4
alongside.
Pandya 2016 power-law and kappa branches (gravitas-core)
physics::synchrotron::pandya_2016_powerlaw_fitlands the §3.2 fit;pandya_2016_kappa_fitinterpolates between thermal (κ → ∞) andpower-law via the published Pandya+ Eq. 37 coefficients.
j_powerlaw_synchrotronandj_kappa_synchrotrontie the fits tothe new
NonThermalPlasmastate.SP-5 5B HorizonPairOverlay (TypeScript)
SVG overlay above the canvas with pair-production particles on a
horizon ring. Spawn rate scales by log10(T_H). ADR-0026 honest-
labelling footer reads "HAWKING PAIRS · ILLUSTRATIVE".
ADRs 0022-0027 (gitignored personal docs)
Six ADRs in
.mayank/decisions/documenting the operator-facingtrade-offs.
Commit 2: Render pipeline + exact Γ amplitude
GLSL shader uniforms
src/shaders/blackhole/chunks/common.tsadds five SP-4 / SP-5uniforms:
u_polarization_enabled,u_stokes,u_active_band_freq_hz,u_active_band_index,u_b_field_strength,u_plunge_envelope_scale. The fragmentshader composites Stokes (Q, U) into the final color via
EVPA-derived hue scaled by linear polarization degree, and the
plunge envelope multiplies disk emission by a brightness factor.
Renderer plumbing
src/rendering/webgl/renderer.tsuploads all five uniforms everyframe.
set1iadded to UniformBatcher for integer uniforms.FeatureTogglesgains optional fields so the operator opts in viathe existing toggle plumbing.
Pandya+ 2016 Eq. 34 exact amplitude
pandya_2016_powerlaw_amplitude(p)evaluates the exact Γ-functionproduct
3^{p/2} (p − 1) Γ((3 p − 1)/12) Γ((3 p + 19)/12) / [2 (p + 1)]in log-space using a 7-term Lanczos approximation accurate to
~5×10⁻¹⁵. Replaces the polynomial fit;
pandya_2016_powerlaw_fitfactors through cleanly.
Test plan
What stays out of scope
goldens are captured; operator runs
bun run shader:update-goldens --confirmonce satisfied.u_b_field_strengthis wiredbut the streamline 3D path tracing earns its own PR.
needs a 4-channel HDR framebuffer and a separate compute pass;
the uniform-vec4 form ships the compositing math without
expanding the framebuffer surface.
not yet consume the new uniforms; only the GLSL (WebGL2) path
does. WGSL parity is its own follow-up.
Summary by CodeRabbit
New Features
Tests