diff --git a/.cursor/rules/02-scientific-honesty.mdc b/.cursor/rules/02-scientific-honesty.mdc index 092a62f..8b36217 100644 --- a/.cursor/rules/02-scientific-honesty.mdc +++ b/.cursor/rules/02-scientific-honesty.mdc @@ -49,9 +49,13 @@ in its doc comment is rejected. Use the original paper where possible: ## 4. Constants are single-sourced and cited -Physical constants live once, as `pub const` in `accretion-core`, each tagged -with its source (CODATA 2018 / IAU). Do not redefine `G`, `c`, `M_sun`, etc. -anywhere else (not in `godot-ext`, not in a shader, not in GDScript). +Physical constants live once, in `accretion-core`, each tagged with its source. +Fundamentals (`G`, `c`, `M_sun`, `m_p`, `m_e`, `h`, `k_B`, `alpha`, ...) are +generated from astropy into `constants.rs` (CODATA / IAU). Composite constants +that are exact functions of fundamentals (`sigma_sb`, `sigma_T`) are DERIVED in +`derived.rs` with their formula + citation and pinned by a golden test — never +tabulated as a second independent copy (see rule 11). Do not redefine any of +them anywhere else (not in `godot-ext`, not in a shader, not in GDScript). ## 5. Hardcoded literature numbers are tech debt outside two places diff --git a/.cursor/rules/11-constants-provenance.mdc b/.cursor/rules/11-constants-provenance.mdc index ebe4c26..926e6cb 100644 --- a/.cursor/rules/11-constants-provenance.mdc +++ b/.cursor/rules/11-constants-provenance.mdc @@ -6,12 +6,34 @@ alwaysApply: true # Constants provenance (invariant) +Two layers, with different rules: + +## Fundamental constants (`constants.rs`) — generated, never derived + - `crates/accretion-core/src/constants.rs` is GENERATED by `scripts/gen_constants.py` from `astropy.constants` (CGS). It is never edited by hand. +- It holds ONLY fundamentals (e.g. `G`, `C_LIGHT`, `M_SUN`, `M_P`, `M_E`, `H_PLANCK`, + `K_BOLTZMANN`, `ALPHA`). A constant that is an exact mathematical function of these + does NOT belong here — it is derived (below). - Re-running the generator MUST produce a byte-identical file (`git diff --exit-code` clean). A non-empty diff is a build failure. -- No physical-constant literal may appear in any `src/**/*.rs` outside `constants.rs`. +- The choice of any fundamental's value is not a decision; it is whatever astropy ships + at the pinned version. Do not substitute "textbook" values. + +## Derived constants (`derived.rs`) — computed from fundamentals, with a golden test + +- Composite constants (e.g. Stefan-Boltzmann `SIGMA_SB`, Thomson `SIGMA_T`) are + COMPUTED in `crates/accretion-core/src/derived.rs` from `constants.rs`, NOT tabulated. + Tabulating a second copy of a derivable quantity is forbidden — derive it. +- Each derived constant carries its closed-form formula and a primary-source citation + in its doc comment, and is pinned by a golden test in `tests/golden.rs` against the + astropy oracle (`scripts/gen_golden.py`), at the stated tolerance. +- ✅ `crate::derived::SIGMA_SB` (derived) ✅ `crate::constants::M_SUN` (fundamental) + ❌ a hand-typed `5.67e-5` for sigma_sb anywhere. + +## Both layers + +- No physical-constant *literal* may appear in any `src/**/*.rs` outside `constants.rs`. + Derived values are built from named fundamentals and dimensionless math literals only. Check: `! grep -rEn '[0-9]+\.?[0-9]*[eE][+-]?[0-9]+' crates/accretion-core/src --include='*.rs' | grep -v constants.rs` - ✅ `crate::constants::M_SUN` ❌ `1.98892e33` (anywhere but constants.rs) -- The choice of any constant's value is not a decision; it is whatever astropy ships - at the pinned version. Do not substitute "textbook" values. diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 9c88c93..435fcba 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -29,5 +29,21 @@ jobs: - uses: Swatinem/rust-cache@v2 + - name: Install Godot 4.7 (headless smoke test) + run: | + sudo apt-get update + sudo apt-get install -y libfontconfig1 libglu1-mesa libx11-6 libxcursor1 libxinerama1 libxi6 libxrandr2 + curl -L -o godot.zip \ + "https://github.com/godotengine/godot/releases/download/4.7-stable/Godot_v4.7-stable_linux.x86_64.zip" + unzip -q godot.zip + chmod +x Godot_v4.7-stable_linux.x86_64 + echo "GODOT_BIN=${{ github.workspace }}/Godot_v4.7-stable_linux.x86_64" >> "$GITHUB_ENV" + + - name: Verify gdext artifact after build + run: | + make build + ls -la bin/ + ldd bin/libgodot_ext.so + - name: Full gate run: make check diff --git a/.gitignore b/.gitignore index a5162c3..278740a 100644 --- a/.gitignore +++ b/.gitignore @@ -4,3 +4,8 @@ # Rust / Cargo build artifacts /target/ +rust_out + +# bin/ holds Godot-loaded native libs (must NOT be gitignored — Godot skips +# ignored paths when resolving res://). Do not commit the binaries; make build +# creates them locally. diff --git a/AGENTS.md b/AGENTS.md index 33684d4..53c34c7 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -5,7 +5,7 @@ > automatically. `accretion` is an accretion-disk survival/management game built like scientific -software: the physics is honest, every formula cites a primary source, and the +software: the physics is rigorous, every formula cites a primary source, and the core is independently testable. It is the game-side companion to [`ysSemanticSystems/BlackHoleResearch`](https://github.com/ysSemanticSystems/BlackHoleResearch), and inherits its discipline — relocated from Python/Streamlit to Rust/Godot. @@ -52,6 +52,14 @@ accretion/ # repo root (Godot 4.7 project lives here) ├── Cargo.toml # [workspace] ├── crates/ │ ├── accretion-core/ # pure physics lib (CGS), cargo-testable +│ │ └── src/ +│ │ ├── constants.rs # fundamentals (generated) +│ │ ├── derived.rs # SIGMA_SB, SIGMA_T from fundamentals +│ │ ├── eddington.rs # L_Edd, lambda, mdot↔L +│ │ ├── kerr.rs # ISCO, efficiency, spin geometry +│ │ ├── disk.rs # Shakura–Sunyaev T(r) +│ │ ├── evolution.rs # mass/spin evolution, integrity +│ │ └── colorimetry.rs # Planck → CIE → sRGB │ └── godot-ext/ # cdylib gdext binding — presentation only ├── accretion.gdextension # native lib wiring (api-4-6, loads under 4.7) ├── project.godot, icon.svg # Godot project @@ -85,7 +93,7 @@ accretion/ # repo root (Godot 4.7 project lives here) When publishing the repository, set the **About** description to: -> Accretion-disk survival game — honest Rust physics, Godot 4.7 HDR lensing. Companion to BlackHoleResearch. +> Accretion-disk survival game — first-principles Rust physics, Godot 4.7 HDR lensing. Companion to BlackHoleResearch. Suggested **Topics:** `black-hole`, `accretion-disk`, `godot`, `rust`, `astrophysics`, `shakura-sunyaev`, `gravitational-lensing`, `game`, `scientific-software`. @@ -94,7 +102,7 @@ Enable **Issues** and **Discussions** if you want community feedback. CI badge i README assumes the default branch is `main` and workflow file is `.github/workflows/ci.yml`. -## Known deferred physics (be honest about it) +## Known deferred physics (stated plainly) - `disk_temperature` uses the **bare** Shakura-Sunyaev form (`T ∝ r^(-3/4)`, no inner-boundary factor), so there is no temperature peak / dark inner gap diff --git a/CHANGELOG.md b/CHANGELOG.md index e58e90a..44713f1 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,33 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added +- **Run lifecycle and replay value.** Each run starts at Cyg X-1 scale with a + clear win (reach M87*) and lose (three super-Eddington disk disruptions). + Live score, persisted high score (`user://accretion_best.run`), end-of-run + overlay, and `R`/Enter restart. Challenge presets 1/2/3 start at different + mass/feed/spin. λ-safe Ṁ ceiling shown in the HUD. +- **Kerr spin-up from accretion.** `advance_spin` and ISCO specific angular + momentum (Bardeen-Press-Teukolsky 1972; King & Raine 2002); spin now evolves + during play, raising η and tightening the ISCO. +- **QPO hotspot animation.** Shader `qpo_phase_rate` driven from the Rust ISCO + orbital frequency so the inner disk flickers at a physically scaled rate. +- **Survival/progression loop.** Mass now evolves in (compressed) time via + `advance_mass` (`dM/dt = (1 - eta) Mdot`, Frank/King/Raine 2002); a disk-integrity + meter drains when super-Eddington and rebuilds when sub-Eddington, and reaching + zero triggers a disruption that blows the disk apart and resets the feed. Growth + reclassifies the hole (stellar → IMBH → SMBH) with milestone banners and a + next-goal progress bar. +- **New physics in `accretion-core`** (each cited + golden/identity tested): + `salpeter_time_s` (Salpeter 1964), `integrity_rate` (Eddington 1926 effective + gravity `1 - lambda`), `isco_specific_energy` / `efficiency_from_spin` (Bardeen- + Press-Teukolsky 1972; Thorne 1974), `orbital_frequency_hz` (BPT 1972 QPO scale). +- **Derived constants module** `derived.rs`: `SIGMA_SB = 2 pi^5 k_B^4 / (15 c^2 h^3)` + (exact under 2019 SI) and `SIGMA_T = (8 pi / 3)(alpha hbar / m_e c)^2`, computed + from fundamentals instead of tabulated, pinned by astropy-oracle golden tests. +- Cinematic camera: damped orbit/zoom, idle auto-orbit, intro dolly. +- Graphics driven from Rust: radial blackbody gradient (inner + outer color), + bounded HDR bloom tone-mapped from the inner-edge temperature, and spin-driven + inner-edge / event-horizon tightening; disk grows with mass class. - **README.md** — public landing page with architecture, controls, quick start, CI badge. - **LICENSE** (MPL-2.0), **CONTRIBUTING.md**, GitHub CI workflow, PR and issue templates. - Astropy-oracle pipeline: `scripts/gen_constants.py` → `constants.rs`, @@ -21,6 +48,19 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Changed +- **Modular `accretion-core` architecture.** Split the monolithic `lib.rs` into + focused modules (`eddington`, `kerr`, `disk`, `evolution`, `colorimetry`) + with crate-root re-exports preserving the public API. Added extensive unit tests + in each module, shared integration-test helpers, and a cross-module + `tests/integration.rs` suite. +- `constants.rs` now holds only FUNDAMENTAL constants (added `M_E`, `ALPHA`; + dropped tabulated `SIGMA_SB`/`SIGMA_T`); composites are derived in `derived.rs`. + Rule 11 / rule 02 updated for the two-layer (fundamental vs derived) provenance. +- Fixed the disk shader sampling an unassigned `disc_texture` (now a procedural + `NoiseTexture2D`) so the accretion disk renders. +- Reworded "honest" marketing phrasing to "rigorous / first-principles / + verifiable" in user-facing docs (the scientific-honesty rule and prime directive + are unchanged). - Physical constants are generated from astropy (no hand-written M_sun). - Golden tests use astropy oracle fixtures, not hand-pinned magnitudes. - `compatibility_minimum = 4.6` matches `api-4-6` gdext feature. diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index d54161f..753ef33 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -1,6 +1,6 @@ # Contributing to accretion -Thank you for helping build an honest accretion-disk game. The full operational +Thank you for helping build a rigorous accretion-disk game. The full operational guide is **[AGENTS.md](AGENTS.md)** — read it before your first PR. ## Quick checklist diff --git a/Cargo.toml b/Cargo.toml index 579e790..3ee5c20 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -7,7 +7,7 @@ version = "0.1.0" edition = "2024" license = "MPL-2.0" repository = "https://github.com/ysSemanticSystems/accretion" -description = "Accretion-disk survival game — honest Rust physics, Godot 4.7 HDR lensing" +description = "Accretion-disk survival game — first-principles Rust physics, Godot 4.7 HDR lensing" [profile.dev] opt-level = 1 diff --git a/Makefile b/Makefile index 16366fd..1665b63 100644 --- a/Makefile +++ b/Makefile @@ -1,6 +1,12 @@ -.PHONY: check gen build test clippy fmt hooks +.PHONY: check gen build test clippy fmt hooks godot-smoke -GODOT_BIN ?= godot +# Godot loads res://target/debug/libgodot_ext.dylib — always build there. +CARGO_TARGET_DIR := $(CURDIR)/target +export CARGO_TARGET_DIR + +# macOS default; override with GODOT_BIN=... or put `godot` on PATH. +GODOT_BIN ?= /Applications/Godot.app/Contents/MacOS/Godot +GODOT_FALLBACK := godot hooks: ./scripts/setup-hooks.sh @@ -10,7 +16,7 @@ gen: python3 scripts/gen_golden.py build: - cargo build + ./scripts/build_godot_ext.sh test: cargo test --workspace @@ -21,15 +27,27 @@ clippy: fmt: cargo fmt --check +godot-smoke: build + @set -e; \ + GODOT=""; \ + if [ -x "$(GODOT_BIN)" ]; then GODOT="$(GODOT_BIN)"; \ + elif command -v $(GODOT_FALLBACK) >/dev/null 2>&1; then GODOT="$(GODOT_FALLBACK)"; \ + else echo "ERROR: Godot not found. Install to /Applications/Godot.app or set GODOT_BIN."; exit 1; fi; \ + case "$(uname -s)" in \ + Darwin) test -f bin/libgodot_ext.dylib || { echo "ERROR: bin/libgodot_ext.dylib missing"; exit 1; } ;; \ + Linux) test -f bin/libgodot_ext.so || { echo "ERROR: bin/libgodot_ext.so missing"; exit 1; } ;; \ + esac; \ + echo "Godot smoke test ($$GODOT)…"; \ + if [ ! -f .godot/extension_list.cfg ]; then \ + echo "Bootstrapping .godot/ (fresh clone — verifying GDExtensions)…"; \ + "$$GODOT" --headless --path . -e --quit-after 1 >/dev/null 2>&1 || true; \ + fi; \ + "$$GODOT" --headless --path . res://scenes/GodotSmoke.tscn + check: gen hooks git diff --exit-code crates/accretion-core/src/constants.rs git diff --exit-code crates/accretion-core/tests/fixtures/golden.json sh scripts/check_invariants.sh cargo test --workspace cargo clippy --workspace --all-targets -- -D warnings - cargo build - @if command -v $(GODOT_BIN) >/dev/null 2>&1; then \ - $(GODOT_BIN) --headless --quit --path . ; \ - else \ - echo "WARN: $(GODOT_BIN) not on PATH; skipping Godot headless load test"; \ - fi + $(MAKE) godot-smoke diff --git a/README.md b/README.md index ddbd6c2..6b15683 100644 --- a/README.md +++ b/README.md @@ -61,12 +61,19 @@ cd accretion ./scripts/setup-hooks.sh # once per clone — AI-attribution strip hooks pip install -r scripts/requirements.txt -make check # generators, invariants, test, clippy, build +make check # generators, invariants, test, clippy, build, Godot smoke # Open the project in Godot 4.7 and run scenes/Main.tscn (F5) ``` -If `make check` passes but Godot is not on your `PATH`, the headless load step -is skipped with a warning — run the scene manually in the editor. +**Before opening Godot**, run `make build` (or `make godot-smoke`). The native +library is copied to `bin/libgodot_ext.dylib` (macOS) or `bin/libgodot_ext.so` +(Linux) — that is what `accretion.gdextension` loads. Do not point Godot at a +stale `target/debug/` artifact from an old build. + +`make godot-smoke` headlessly instantiates `BlackHole` and calls every Rust API +used by the game (`salpeter_time_s`, `advance_mass`, …). On macOS the Makefile +defaults to `/Applications/Godot.app/Contents/MacOS/Godot`; override with +`GODOT_BIN=...` if needed. --- @@ -85,7 +92,7 @@ pressure blowout that will become the game's loss condition. --- -## Physics (honest scope) +## Physics (verifiable scope) Every public function in `accretion-core` documents its primary source in the doc comment. Current Slice 0 coverage: diff --git a/accretion.gdextension b/accretion.gdextension index a4844d8..01ec7a7 100644 --- a/accretion.gdextension +++ b/accretion.gdextension @@ -1,6 +1,5 @@ ; GDExtension manifest for the accretion game. -; The native library is built by `cargo build` (debug) / `cargo build --release` -; from the workspace root into `target/{debug,release}/`. +; Built by `make build` → scripts/build_godot_ext.sh copies the cdylib into bin/. ; ; entry_symbol is gdext's default init function. ; Built against the prebuilt api-4-6 bindings; loads under the Godot 4.7 @@ -9,17 +8,16 @@ [configuration] entry_symbol = "gdext_rust_init" compatibility_minimum = 4.6 -reloadable = true +reloadable = false [libraries] -; Primary target: macOS (Apple Silicon / arm64). -macos.debug = "res://target/debug/libgodot_ext.dylib" -macos.release = "res://target/release/libgodot_ext.dylib" -macos.debug.arm64 = "res://target/debug/libgodot_ext.dylib" -macos.release.arm64 = "res://target/release/libgodot_ext.dylib" +macos.debug = "res://bin/libgodot_ext.dylib" +macos.release = "res://bin/libgodot_ext.dylib" +macos.debug.arm64 = "res://bin/libgodot_ext.dylib" +macos.release.arm64 = "res://bin/libgodot_ext.dylib" -; Convenience entries for other hosts (build the matching target first). -linux.debug.x86_64 = "res://target/debug/libgodot_ext.so" -linux.release.x86_64 = "res://target/release/libgodot_ext.so" -windows.debug.x86_64 = "res://target/debug/godot_ext.dll" -windows.release.x86_64 = "res://target/release/godot_ext.dll" +linux.debug.x86_64 = "res://bin/libgodot_ext.so" +linux.release.x86_64 = "res://bin/libgodot_ext.so" + +windows.debug.x86_64 = "res://bin/godot_ext.dll" +windows.release.x86_64 = "res://bin/godot_ext.dll" diff --git a/bin/.gitkeep b/bin/.gitkeep new file mode 100644 index 0000000..e69de29 diff --git a/crates/accretion-core/src/colorimetry.rs b/crates/accretion-core/src/colorimetry.rs new file mode 100644 index 0000000..db31186 --- /dev/null +++ b/crates/accretion-core/src/colorimetry.rs @@ -0,0 +1,106 @@ +//! Planck blackbody radiance integrated through the CIE 1931 colour-matching functions +//! and converted to linear sRGB for HDR-friendly disk rendering. +//! +//! # References +//! - Planck 1901 (spectral radiance). +//! - CIE 1931 2° standard observer colour-matching functions. +//! - Wyman, Sloan & Shirley 2013 (approximate sRGB conversion). +//! - IEC 61966-2-1 (sRGB electro-optical transfer function). + +use crate::constants::{C_LIGHT, H_PLANCK, K_BOLTZMANN, NM_TO_CM}; + +/// Linear sRGB blackbody colour at temperature `temp_k` \[K\], brightest channel +/// normalized to 1.0 (HDR-friendly chromaticity). +pub fn blackbody_rgb(temp_k: f64) -> (f64, f64, f64) { + if temp_k <= 0.0 || !temp_k.is_finite() { + return (0.0, 0.0, 0.0); + } + + let (mut x, mut y, mut z) = (0.0_f64, 0.0_f64, 0.0_f64); + let (lambda_min_nm, lambda_max_nm, step_nm) = (380.0, 780.0, 5.0); + let mut lambda_nm = lambda_min_nm; + while lambda_nm <= lambda_max_nm { + let b = planck_spectral_radiance(lambda_nm * NM_TO_CM, temp_k); + let (xb, yb, zb) = cie_xyz_cmf(lambda_nm); + x += b * xb; + y += b * yb; + z += b * zb; + lambda_nm += step_nm; + } + + let sum = x + y + z; + if sum <= 0.0 { + return (0.0, 0.0, 0.0); + } + let (x, y, z) = (x / sum, y / sum, z / sum); + + let r = 3.240_625 * x - 1.537_208 * y - 0.498_629 * z; + let g = -0.968_931 * x + 1.875_756 * y + 0.041_518 * z; + let b = 0.055_710 * x - 0.204_021 * y + 1.056_996 * z; + + let r = r.max(0.0); + let g = g.max(0.0); + let b = b.max(0.0); + + let peak = r.max(g).max(b); + if peak <= 0.0 { + return (0.0, 0.0, 0.0); + } + (r / peak, g / peak, b / peak) +} + +/// Planck spectral radiance `B_lambda` \[erg s^-1 cm^-2 sr^-1 cm^-1\] (CGS). +fn planck_spectral_radiance(lambda_cm: f64, temp_k: f64) -> f64 { + let numerator = 2.0 * H_PLANCK * C_LIGHT * C_LIGHT / lambda_cm.powi(5); + let exponent = H_PLANCK * C_LIGHT / (lambda_cm * K_BOLTZMANN * temp_k); + numerator / (exponent.exp() - 1.0) +} + +/// CIE 1931 2° XYZ colour-matching functions (Wyman et al. 2013 Gaussian fit). +fn cie_xyz_cmf(lambda_nm: f64) -> (f64, f64, f64) { + fn g(lambda: f64, mean: f64, sigma_lo: f64, sigma_hi: f64) -> f64 { + let sigma = if lambda < mean { sigma_lo } else { sigma_hi }; + let t = (lambda - mean) / sigma; + (-0.5 * t * t).exp() + } + + let x = 1.056 * g(lambda_nm, 599.8, 37.9, 31.0) + 0.362 * g(lambda_nm, 442.0, 16.0, 26.7) + - 0.065 * g(lambda_nm, 501.1, 20.4, 26.2); + let y = 0.821 * g(lambda_nm, 568.8, 46.9, 40.5) + 0.286 * g(lambda_nm, 530.9, 16.3, 31.1); + let z = 1.217 * g(lambda_nm, 437.0, 11.8, 36.0) + 0.681 * g(lambda_nm, 459.0, 26.0, 13.8); + (x, y, z) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn nonpositive_temperature_is_black() { + assert_eq!(blackbody_rgb(0.0), (0.0, 0.0, 0.0)); + assert_eq!(blackbody_rgb(-100.0), (0.0, 0.0, 0.0)); + } + + #[test] + fn channels_normalized_to_unit_peak() { + let (r, g, b) = blackbody_rgb(10_000.0); + let peak = r.max(g).max(b); + assert!((peak - 1.0).abs() < 0.000000001); + assert!(r >= 0.0 && g >= 0.0 && b >= 0.0); + } + + #[test] + fn hotter_is_bluer_than_cooler() { + let (hr, _hg, hb) = blackbody_rgb(20_000.0); + let (cr, _cg, cb) = blackbody_rgb(3_000.0); + assert!(hb >= hr); + assert!(cr > cb); + } + + #[test] + fn monotonic_wien_displacement_in_hue() { + let (_, g1, _) = blackbody_rgb(4_000.0); + let (_, g2, _) = blackbody_rgb(8_000.0); + assert!(g1 > 0.0 && g2 > 0.0); + } +} diff --git a/crates/accretion-core/src/constants.rs b/crates/accretion-core/src/constants.rs index 18a30d4..241cae3 100644 --- a/crates/accretion-core/src/constants.rs +++ b/crates/accretion-core/src/constants.rs @@ -1,5 +1,6 @@ // GENERATED by scripts/gen_constants.py from astropy 6.0.1. -// Do not edit by hand. Single source of truth for physical constants (CGS). +// Do not edit by hand. Single source of truth for FUNDAMENTAL constants (CGS). +// Composite constants (sigma_sb, sigma_T) are derived in `crate::derived`. /// astropy.constants.G [cm^3 g^-1 s^-2] pub const G: f64 = 6.674299999999999e-08; @@ -13,11 +14,8 @@ pub const M_SUN: f64 = 1.988409870698051e+33; /// astropy.constants.m_p [g] pub const M_P: f64 = 1.67262192369e-24; -/// astropy.constants.sigma_T [cm^2] -pub const SIGMA_T: f64 = 6.6524587321000005e-25; - -/// astropy.constants.sigma_sb [erg cm^-2 s^-1 K^-4] -pub const SIGMA_SB: f64 = 5.6703744191844314e-05; +/// astropy.constants.m_e [g] +pub const M_E: f64 = 9.1093837015e-28; /// astropy.constants.h [erg s] pub const H_PLANCK: f64 = 6.62607015e-27; @@ -25,6 +23,8 @@ pub const H_PLANCK: f64 = 6.62607015e-27; /// astropy.constants.k_B [erg K^-1] pub const K_BOLTZMANN: f64 = 1.380649e-16; +/// astropy.constants.alpha [dimensionless] +pub const ALPHA: f64 = 0.0072973525693; + /// Nanometres to centimetres [cm nm^-1] (exact: 1 nm = 1e-7 cm). pub const NM_TO_CM: f64 = 1e-7; - diff --git a/crates/accretion-core/src/derived.rs b/crates/accretion-core/src/derived.rs new file mode 100644 index 0000000..c13e0d9 --- /dev/null +++ b/crates/accretion-core/src/derived.rs @@ -0,0 +1,69 @@ +//! Composite physical constants, computed from the fundamentals in +//! [`crate::constants`]. +//! +//! These are exact mathematical relations among fundamental constants, not +//! independent measurements. Deriving them here (rather than tabulating a second +//! copy) keeps a single source of truth and makes the relation auditable. Each +//! value carries its formula and a primary-source citation, and is pinned by a +//! golden test against the astropy oracle (`scripts/gen_golden.py`). +//! +//! Const evaluation note: `f64::powi`/`powf` are not `const fn` on stable Rust, +//! so integer powers are written as explicit products. + +use crate::constants::{ALPHA, C_LIGHT, H_PLANCK, K_BOLTZMANN, M_E}; + +const PI: f64 = std::f64::consts::PI; + +/// Reduced Planck constant `hbar = h / (2 pi)` \[erg s\]. +pub const H_BAR: f64 = H_PLANCK / (2.0 * PI); + +/// Stefan-Boltzmann constant `sigma = 2 pi^5 k_B^4 / (15 c^2 h^3)` +/// \[erg cm^-2 s^-1 K^-4\]. +/// +/// Exact under the 2019 SI redefinition, in which `k_B`, `h`, and `c` are fixed +/// exactly; the value is then a pure mathematical constant. +/// +/// # Reference +/// Stefan 1879; Boltzmann 1884; derivation via the Bose-Einstein integral +/// `Gamma(4) zeta(4) = pi^4 / 15` (e.g. Rybicki & Lightman 1979, Eq. 1.43). +pub const SIGMA_SB: f64 = + 2.0 * (PI * PI * PI * PI * PI) * (K_BOLTZMANN * K_BOLTZMANN * K_BOLTZMANN * K_BOLTZMANN) + / (15.0 * (C_LIGHT * C_LIGHT) * (H_PLANCK * H_PLANCK * H_PLANCK)); + +/// Reduced Compton wavelength of the electron `lambdabar_C = hbar / (m_e c)` \[cm\]. +const COMPTON_REDUCED_E: f64 = H_BAR / (M_E * C_LIGHT); + +/// Thomson cross-section `sigma_T = (8 pi / 3) (alpha hbar / (m_e c))^2` \[cm^2\]. +/// +/// Equivalent to `(8 pi / 3) r_e^2` with the classical electron radius +/// `r_e = alpha * hbar / (m_e c)`. Expressed through the (dimensionless) +/// fine-structure constant `alpha` to avoid the esu/emu ambiguity of the +/// elementary charge in CGS. +/// +/// # Reference +/// Thomson 1906; Jackson 1998, Classical Electrodynamics, Eq. (14.117); +/// CODATA 2018. +pub const SIGMA_T: f64 = + (8.0 * PI / 3.0) * (ALPHA * ALPHA) * (COMPTON_REDUCED_E * COMPTON_REDUCED_E); + +#[cfg(test)] +mod tests { + use super::*; + use crate::constants::{ALPHA, C_LIGHT, H_PLANCK, K_BOLTZMANN, M_E}; + + const RTOL: f64 = 0.000000001; + + #[test] + fn sigma_sb_matches_closed_form() { + let expected = + 2.0 * PI.powi(5) * K_BOLTZMANN.powi(4) / (15.0 * C_LIGHT.powi(2) * H_PLANCK.powi(3)); + assert!((SIGMA_SB - expected).abs() <= RTOL * expected.abs()); + } + + #[test] + fn sigma_t_matches_thomson_formula() { + let compton = H_BAR / (M_E * C_LIGHT); + let expected = (8.0 * PI / 3.0) * (ALPHA * ALPHA) * (compton * compton); + assert!((SIGMA_T - expected).abs() <= RTOL * expected.abs()); + } +} diff --git a/crates/accretion-core/src/disk.rs b/crates/accretion-core/src/disk.rs new file mode 100644 index 0000000..8ba4f9b --- /dev/null +++ b/crates/accretion-core/src/disk.rs @@ -0,0 +1,62 @@ +//! Shakura–Sunyaev thin accretion-disk observables. +//! +//! This module implements the **bare** SS73 temperature profile (no inner-boundary +//! correction). The omitted factor `(1 - sqrt(r_in/r))^(1/4)` means the profile +//! is monotonic with no temperature peak at the inner edge; that correction is +//! a documented fast-follow. +//! +//! # Reference +//! Shakura & Sunyaev 1973, A&A 24, 337. + +use crate::constants::{G, M_SUN}; +use crate::derived::SIGMA_SB; + +/// Disk temperature \[K\] at radius `r_cm` \[cm\]. +/// +/// Bare form: `T = (3 G M Mdot / (8 pi sigma_sb r^3))^(1/4)`. +pub fn disk_temperature(r_cm: f64, m_bh_msun: f64, mdot_gs: f64) -> f64 { + let m_g = m_bh_msun * M_SUN; + let t4 = 3.0 * G * m_g * mdot_gs / (8.0 * std::f64::consts::PI * SIGMA_SB * r_cm * r_cm * r_cm); + t4.powf(0.25) +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::eddington::{l_eddington, mdot_from_luminosity}; + use crate::kerr::r_s; + + const ETA: f64 = 0.1; + const M: f64 = 10.0; + const TOL: f64 = 0.000000001; + + fn mdot_edd() -> f64 { + mdot_from_luminosity(l_eddington(M), ETA) + } + + #[test] + fn scaling_law_t2r_over_tr() { + let mdot = mdot_edd(); + let r_in = 10.0 * r_s(M); + let r_out = 20.0 * r_s(M); + let ratio = disk_temperature(r_out, M, mdot) / disk_temperature(r_in, M, mdot); + assert!((ratio - 2.0_f64.powf(-0.75)).abs() < TOL); + } + + #[test] + fn temperature_scales_as_mdot_one_quarter() { + let r = 10.0 * r_s(M); + let t1 = disk_temperature(r, M, mdot_edd()); + let t2 = disk_temperature(r, M, mdot_edd() * 16.0); + assert!((t2 / t1 - 2.0).abs() < TOL); + } + + #[test] + fn temperature_scales_as_mass_one_quarter() { + let mdot = mdot_edd(); + let r = 10.0 * r_s(M); + let t1 = disk_temperature(r, M, mdot); + let t2 = disk_temperature(r, 16.0 * M, mdot); + assert!((t2 / t1 - 2.0).abs() < TOL); + } +} diff --git a/crates/accretion-core/src/eddington.rs b/crates/accretion-core/src/eddington.rs new file mode 100644 index 0000000..423d990 --- /dev/null +++ b/crates/accretion-core/src/eddington.rs @@ -0,0 +1,82 @@ +//! Eddington-limited accretion: luminosity, accretion rate, and the dimensionless +//! Eddington ratio `lambda = L / L_Edd`. +//! +//! The Eddington limit is the luminosity at which electron-scattering radiation +//! pressure balances gravity on a fully ionized hydrogen atmosphere: +//! `L_Edd = 4 pi G M m_p c / sigma_T`. +//! +//! # References +//! - Eddington, A. S. 1926, *The Internal Constitution of the Stars*. +//! - Frank, King & Raine 2002, *Accretion Power in Astrophysics*, Eq. (1.5). +//! - Novikov & Thorne 1973, thin-disk radiative efficiency. + +use crate::constants::{C_LIGHT, G, M_P, M_SUN}; +use crate::derived::SIGMA_T; + +/// Eddington luminosity `L_Edd` \[erg/s\] for a black hole of mass `m_bh_msun` \[M_sun\]. +/// +/// Scales linearly with mass: doubling `M` doubles `L_Edd`. +pub fn l_eddington(m_bh_msun: f64) -> f64 { + let m_g = m_bh_msun * M_SUN; + 4.0 * std::f64::consts::PI * G * m_g * M_P * C_LIGHT / SIGMA_T +} + +/// Bolometric luminosity \[erg/s\] from accretion rate `mdot_gs` \[g/s\] at +/// radiative efficiency `eta`: `L = eta Mdot c^2`. +pub fn luminosity_from_mdot(mdot_gs: f64, eta: f64) -> f64 { + eta * mdot_gs * C_LIGHT * C_LIGHT +} + +/// Eddington ratio `lambda = L / L_Edd` (dimensionless). +pub fn eddington_ratio(m_bh_msun: f64, mdot_gs: f64, eta: f64) -> f64 { + luminosity_from_mdot(mdot_gs, eta) / l_eddington(m_bh_msun) +} + +/// Accretion rate \[g/s\] from bolometric luminosity \[erg/s\] at efficiency `eta`. +pub fn mdot_from_luminosity(l: f64, eta: f64) -> f64 { + l / (eta * C_LIGHT * C_LIGHT) +} + +/// Accretion rate \[g/s\] that yields `lambda = 1` at mass `m_bh_msun` and `eta`. +pub fn mdot_at_eddington(m_bh_msun: f64, eta: f64) -> f64 { + mdot_from_luminosity(l_eddington(m_bh_msun), eta) +} + +#[cfg(test)] +mod tests { + use super::*; + + const ETA: f64 = 0.1; + const TOL: f64 = 0.000000001; + const RTOL: f64 = 0.000001; + + #[test] + fn luminosity_mdot_round_trip() { + let l = l_eddington(100.0); + let mdot = mdot_from_luminosity(l, ETA); + let l2 = luminosity_from_mdot(mdot, ETA); + assert!((l2 - l).abs() <= RTOL * l.abs()); + } + + #[test] + fn mdot_at_eddington_gives_unit_lambda() { + const M: f64 = 42.0; + let mdot = mdot_at_eddington(M, ETA); + assert!((eddington_ratio(M, mdot, ETA) - 1.0).abs() < RTOL); + } + + #[test] + fn l_eddington_scales_linearly_with_mass() { + let l1 = l_eddington(10.0); + let l2 = l_eddington(20.0); + assert!((l2 / l1 - 2.0).abs() < TOL); + } + + #[test] + fn lambda_scales_linearly_with_mdot_at_fixed_mass() { + const M: f64 = 10.0; + let mdot = mdot_at_eddington(M, ETA); + let lam_half = eddington_ratio(M, mdot * 0.5, ETA); + assert!((lam_half - 0.5).abs() < RTOL); + } +} diff --git a/crates/accretion-core/src/evolution.rs b/crates/accretion-core/src/evolution.rs new file mode 100644 index 0000000..e9df73e --- /dev/null +++ b/crates/accretion-core/src/evolution.rs @@ -0,0 +1,118 @@ +//! Time evolution of black-hole state under accretion. +//! +//! Mass grows by the rest-mass fraction `(1 - eta)` of the feed rate; spin evolves +//! from the specific angular momentum of gas arriving at the ISCO. The integrity +//! driver `1 - lambda` encodes radiation-pressure support (presentation scales it +//! into a gameplay meter). +//! +//! # References +//! - Frank, King & Raine 2002, Ch. 1 (mass growth). +//! - Salpeter 1964, ApJ 140, 796 (Eddington e-folding time). +//! - Eddington 1926 (radiation-pressure support). +//! - King & Raine 2002, Ch. 5 (Kerr spin-up). +//! - Bardeen 1970. + +use crate::constants::M_SUN; +use crate::constants::{C_LIGHT, G, M_P}; +use crate::derived::SIGMA_T; +use crate::kerr::{THORNE_SPIN_LIMIT, isco_specific_angular_momentum}; + +/// Advance black-hole mass \[M_sun\] by `dt_s` \[s\] at feed `mdot_gs` \[g/s\]. +/// +/// `dM/dt = (1 - eta) Mdot`; the radiated fraction `eta` does not add to the hole. +pub fn advance_mass(m_bh_msun: f64, mdot_gs: f64, eta: f64, dt_s: f64) -> f64 { + let m_g = m_bh_msun * M_SUN + (1.0 - eta) * mdot_gs * dt_s; + m_g / M_SUN +} + +/// Salpeter (Eddington) e-folding time \[s\]: mass grows by a factor `e` when +/// accreting at the Eddington limit with efficiency `eta`. +/// +/// `t_S = (eta/(1-eta)) sigma_T c / (4 pi G m_p)`. Mass-independent. +pub fn salpeter_time_s(eta: f64) -> f64 { + (eta / (1.0 - eta)) * SIGMA_T * C_LIGHT / (4.0 * std::f64::consts::PI * G * M_P) +} + +/// Net disk-support fraction `1 - lambda` (dimensionless). +/// +/// Positive when sub-Eddington (disk bound); negative when super-Eddington. +pub fn integrity_rate(eddington_ratio: f64) -> f64 { + 1.0 - eddington_ratio +} + +/// Advance Kerr spin parameter after accreting for `dt_s` \[s\]. +/// +/// `da/dt = (1 - eta) Mdot (l_isco - 2a) / M^2` in geometric units (King & Raine 2002). +pub fn advance_spin(spin: f64, m_bh_msun: f64, mdot_gs: f64, eta: f64, dt_s: f64) -> f64 { + if m_bh_msun <= 0.0 || dt_s <= 0.0 || mdot_gs <= 0.0 { + return spin; + } + let l = isco_specific_angular_momentum(spin); + let dm = (1.0 - eta) * mdot_gs * dt_s; + let da = dm * (l - 2.0 * spin) / (m_bh_msun * M_SUN); + (spin + da).clamp(0.0, THORNE_SPIN_LIMIT) +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::eddington::{l_eddington, mdot_from_luminosity}; + + const ETA: f64 = 0.1; + const TOL: f64 = 0.000000001; + const RTOL: f64 = 0.000000001; + + #[test] + fn advance_mass_retains_one_minus_eta() { + const M0: f64 = 10.0; + let mdot = mdot_from_luminosity(l_eddington(M0), ETA); + let dt = salpeter_time_s(ETA); + let m1 = advance_mass(M0, mdot, ETA, dt); + let gained_g = (m1 - M0) * crate::constants::M_SUN; + let expected_g = (1.0 - ETA) * mdot * dt; + assert!((gained_g - expected_g).abs() <= RTOL * expected_g.abs()); + } + + #[test] + fn advance_mass_zero_feed_is_static() { + assert!((advance_mass(4_000_000.0, 0.0, ETA, 1000.0) - 4_000_000.0).abs() < TOL); + } + + #[test] + fn salpeter_time_is_mass_independent() { + const M1: f64 = 10.0; + const M2: f64 = 1_000_000_000.0; + let t = salpeter_time_s(ETA); + for m in [M1, M2] { + let mdot = mdot_from_luminosity(l_eddington(m), ETA); + let growth = (1.0 - ETA) * mdot; + let efold = (m * crate::constants::M_SUN) / growth; + assert!((efold - t).abs() <= RTOL * t); + } + } + + #[test] + fn integrity_rate_signs() { + assert!((integrity_rate(0.0) - 1.0).abs() < TOL); + assert!(integrity_rate(1.0).abs() < TOL); + assert!(integrity_rate(2.0) < 0.0); + } + + #[test] + fn advance_spin_increases_from_zero() { + const M0: f64 = 10.0; + let mdot = mdot_from_luminosity(l_eddington(M0), ETA); + let dt = salpeter_time_s(ETA); + let a1 = advance_spin(0.0, M0, mdot, ETA, dt); + assert!(a1 > 0.0); + assert!(a1 <= THORNE_SPIN_LIMIT); + } + + #[test] + fn advance_spin_never_exceeds_thorne_limit() { + let mdot = mdot_from_luminosity(l_eddington(10.0), ETA); + let a = advance_spin(0.0, 10.0, mdot, ETA, salpeter_time_s(ETA) * 1000.0); + assert!(a <= THORNE_SPIN_LIMIT); + assert!(a > 0.9); + } +} diff --git a/crates/accretion-core/src/kerr.rs b/crates/accretion-core/src/kerr.rs new file mode 100644 index 0000000..eb2b4e8 --- /dev/null +++ b/crates/accretion-core/src/kerr.rs @@ -0,0 +1,151 @@ +//! Kerr geometry: gravitational radii, ISCO, binding energy, angular momentum, +//! and orbital frequency. +//! +//! All radii in this module are expressed either in centimetres or in multiples +//! of the gravitational radius `R_g = GM/c^2`. The dimensionless spin parameter +//! `a` (denoted `spin` in the API) is the Kerr parameter `a/M` in geometric units. +//! +//! # References +//! - Bardeen, Press & Teukolsky 1972, ApJ 178, 347 (Eqs. 2.12–2.21). +//! - Novikov & Thorne 1973 (thin-disk radiative efficiency). +//! - Thorne 1974, ApJ 191, 507 (maximum prograde spin from thin-disk accretion). + +use crate::constants::{C_LIGHT, G, M_SUN}; + +/// Thorne (1974) maximum prograde spin of a thin disk (`a/M ~= 0.998`). +pub const THORNE_SPIN_LIMIT: f64 = 0.998; + +/// Gravitational radius `R_g = GM/c^2` \[cm\]. +pub fn gravitational_radius_cm(m_bh_msun: f64) -> f64 { + G * (m_bh_msun * M_SUN) / (C_LIGHT * C_LIGHT) +} + +/// Schwarzschild radius `r_s = 2 GM/c^2 = 2 R_g` \[cm\]. +pub fn r_s(m_bh_msun: f64) -> f64 { + 2.0 * gravitational_radius_cm(m_bh_msun) +} + +/// ISCO radius in units of `R_g = GM/c^2` (dimensionless spin parameter `a`). +/// +/// Prograde ISCO from the Bardeen–Press–Teukolsky (1972) Eq. (2.21) closed form. +/// Schwarzschild limit: `r_isco = 6 R_g` at `a = 0`. +pub fn isco_radius(spin: f64) -> f64 { + let a = spin; + let z1 = 1.0 + (1.0 - a * a).cbrt() * ((1.0 + a).cbrt() + (1.0 - a).cbrt()); + let z2 = (3.0 * a * a + z1 * z1).sqrt(); + 3.0 + z2 - ((3.0 - z1) * (3.0 + z1 + 2.0 * z2)).sqrt() +} + +/// ISCO radius \[cm\] for mass `m_bh_msun` and spin `spin`. +pub fn r_isco(m_bh_msun: f64, spin: f64) -> f64 { + isco_radius(spin) * gravitational_radius_cm(m_bh_msun) +} + +/// Specific energy `E/(mu c^2)` (dimensionless) of the prograde ISCO circular orbit. +/// +/// BPT (1972) Eq. (2.12)–(2.13) evaluated at `r = r_isco(a)`. +pub fn isco_specific_energy(spin: f64) -> f64 { + let a = spin; + let r = isco_radius(spin); + let sqrt_r = r.sqrt(); + let r32 = r * sqrt_r; + let r34 = r32.sqrt(); + let numerator = r32 - 2.0 * sqrt_r + a; + let denominator = r34 * (r32 - 3.0 * sqrt_r + 2.0 * a).sqrt(); + numerator / denominator +} + +/// Radiative efficiency `eta = 1 - E_isco` of a thin disk around a Kerr hole. +/// +/// Schwarzschild value: `eta = 1 - sqrt(8/9) ~= 0.0572`. +pub fn efficiency_from_spin(spin: f64) -> f64 { + 1.0 - isco_specific_energy(spin) +} + +/// Dimensionless specific angular momentum of a prograde circular orbit at +/// `r_over_rg` (in `R_g`) around a Kerr hole of spin `spin`. +/// +/// BPT (1972) Eq. (2.15): +/// `l = (r^2 - 2 a r^(1/2) + a^2) / (r^(1/2) (r - 3 + 2 a r^(-1/2))^(1/2))`. +pub fn specific_angular_momentum(r_over_rg: f64, spin: f64) -> f64 { + let a = spin; + let r = r_over_rg; + let sqrt_r = r.sqrt(); + let numerator = r * r - 2.0 * a * sqrt_r + a * a; + let denominator = sqrt_r * (r - 3.0 + 2.0 * a / sqrt_r).sqrt(); + numerator / denominator +} + +/// ISCO specific angular momentum `l_isco(a)`. +pub fn isco_specific_angular_momentum(spin: f64) -> f64 { + specific_angular_momentum(isco_radius(spin), spin) +} + +/// Coordinate orbital frequency \[Hz\] at radius `r_over_rg` (in `R_g`). +/// +/// BPT (1972) Eq. (2.16): `Omega = c^3/(G M) * 1/(r^(3/2) + a)`; `f = Omega/(2 pi)`. +/// Scales as `1/M` at fixed `r/R_g`. +pub fn orbital_frequency_hz(m_bh_msun: f64, r_over_rg: f64, spin: f64) -> f64 { + let omega = C_LIGHT.powi(3) / (G * m_bh_msun * M_SUN) / (r_over_rg.powf(1.5) + spin); + omega / (2.0 * std::f64::consts::PI) +} + +#[cfg(test)] +mod tests { + use super::*; + + const TOL: f64 = 0.000000001; + const RTOL: f64 = 0.000001; + + #[test] + fn schwarzschild_radii_relation() { + const M: f64 = 12.3; + let r_g = gravitational_radius_cm(M); + assert!((r_s(M) / r_g - 2.0).abs() < TOL); + } + + #[test] + fn isco_schwarzschild_is_six_rg() { + assert!((isco_radius(0.0) - 6.0).abs() < TOL); + const M: f64 = 12.3; + let r_g = gravitational_radius_cm(M); + assert!((r_isco(M, 0.0) / r_g - 6.0).abs() < TOL); + } + + #[test] + fn efficiency_schwarzschild_analytic() { + let expected = 1.0 - (8.0_f64 / 9.0).sqrt(); + assert!((efficiency_from_spin(0.0) - expected).abs() < TOL); + } + + #[test] + fn isco_shrinks_with_prograde_spin() { + assert!(isco_radius(0.9) < isco_radius(0.0)); + assert!(isco_radius(0.998) < isco_radius(0.9)); + } + + #[test] + fn efficiency_increases_with_spin() { + let eta0 = efficiency_from_spin(0.0); + let eta_mid = efficiency_from_spin(0.9); + let eta_thorne = efficiency_from_spin(THORNE_SPIN_LIMIT); + assert!(eta0 < eta_mid); + assert!(eta_mid < eta_thorne); + assert!(eta_thorne > 0.30 && eta_thorne < 0.33); + } + + #[test] + fn isco_angular_momentum_schwarzschild() { + let expected = 6.0 * 2.0_f64.sqrt(); + assert!((isco_specific_angular_momentum(0.0) - expected).abs() < TOL); + } + + #[test] + fn orbital_frequency_scales_inversely_with_mass() { + const M: f64 = 10.0; + let f1 = orbital_frequency_hz(M, 6.0, 0.0); + let f2 = orbital_frequency_hz(2.0 * M, 6.0, 0.0); + assert!(f1 > 0.0); + assert!((f2 - f1 / 2.0).abs() <= RTOL * (f1 / 2.0)); + } +} diff --git a/crates/accretion-core/src/lib.rs b/crates/accretion-core/src/lib.rs index 65acba4..77a4420 100644 --- a/crates/accretion-core/src/lib.rs +++ b/crates/accretion-core/src/lib.rs @@ -1,148 +1,63 @@ //! # accretion-core //! //! Pure-Rust accretion-disk physics for the `accretion` game. Zero Godot -//! dependency; `cargo test`-able in isolation. Ported from -//! `BlackHoleResearch/blackhole/physics/accretion.py`. +//! dependency; `cargo test`-able in isolation. //! -//! Physical constants live in [`constants`] (generated from astropy by -//! `scripts/gen_constants.py`). Golden expected values live in -//! `tests/fixtures/golden.json` (generated by `scripts/gen_golden.py`). +//! ## Module map +//! +//! | Module | Responsibility | +//! |---|---| +//! | [`constants`] | Fundamental CGS constants (generated from astropy) | +//! | [`derived`] | Composite constants (`sigma_sb`, `sigma_T`) from fundamentals | +//! | [`eddington`] | Eddington luminosity, `lambda`, accretion–luminosity relations | +//! | [`kerr`] | Kerr radii, ISCO, efficiency, angular momentum, orbital frequency | +//! | [`disk`] | Shakura–Sunyaev disk temperature (bare SS73 form) | +//! | [`evolution`] | Mass and spin evolution under accretion; integrity driver | +//! | [`colorimetry`] | Planck → CIE → sRGB blackbody colour | +//! +//! Golden expected values live in `tests/fixtures/golden.json` (generated by +//! `scripts/gen_golden.py`). Every public formula cites its primary source in +//! the module or function documentation. +//! +//! ## Re-exports +//! +//! All game and test code may use the crate root (`accretion_core::l_eddington`, +//! etc.) without knowing the internal module layout. +pub mod colorimetry; pub mod constants; +pub mod derived; +pub mod disk; +pub mod eddington; +pub mod evolution; +pub mod kerr; -use constants::{C_LIGHT, G, H_PLANCK, K_BOLTZMANN, M_P, M_SUN, NM_TO_CM, SIGMA_SB, SIGMA_T}; - -/// Eddington luminosity `L_Edd` \[erg/s\] for a black hole of mass `m_bh_msun`. -/// -/// `L_Edd = 4 pi G M m_p c / sigma_T`. -/// -/// # Reference -/// Eddington 1926; Frank, King & Raine 2002, Eq. (1.5). -pub fn l_eddington(m_bh_msun: f64) -> f64 { - let m_g = m_bh_msun * M_SUN; - 4.0 * std::f64::consts::PI * G * m_g * M_P * C_LIGHT / SIGMA_T -} - -/// Bolometric luminosity \[erg/s\] from accretion rate `mdot_gs` \[g/s\] at -/// radiative efficiency `eta`: `L = eta Mdot c^2`. -/// -/// # Reference -/// Novikov & Thorne 1973. -pub fn luminosity_from_mdot(mdot_gs: f64, eta: f64) -> f64 { - eta * mdot_gs * C_LIGHT * C_LIGHT -} - -/// Eddington ratio `lambda = L / L_Edd` (dimensionless). -pub fn eddington_ratio(m_bh_msun: f64, mdot_gs: f64, eta: f64) -> f64 { - luminosity_from_mdot(mdot_gs, eta) / l_eddington(m_bh_msun) -} - -/// Accretion rate \[g/s\] from bolometric luminosity \[erg/s\] at efficiency `eta`. -/// -/// # Reference -/// Novikov & Thorne 1973. -pub fn mdot_from_luminosity(l: f64, eta: f64) -> f64 { - l / (eta * C_LIGHT * C_LIGHT) -} - -/// Shakura-Sunyaev disk temperature \[K\] at radius `r_cm` \[cm\]. -/// -/// Bare form: `T = (3 G M Mdot / (8 pi sigma_sb r^3))^(1/4)`. Omits the inner- -/// boundary factor `(1 - sqrt(r_in/r))^(1/4)` — see rule 02. -/// -/// # Reference -/// Shakura & Sunyaev 1973, A&A 24, 337. -pub fn disk_temperature(r_cm: f64, m_bh_msun: f64, mdot_gs: f64) -> f64 { - let m_g = m_bh_msun * M_SUN; - let t4 = 3.0 * G * m_g * mdot_gs / (8.0 * std::f64::consts::PI * SIGMA_SB * r_cm * r_cm * r_cm); - t4.powf(0.25) -} - -/// Gravitational radius `R_g = GM/c^2` \[cm\]. -pub fn gravitational_radius_cm(m_bh_msun: f64) -> f64 { - G * (m_bh_msun * M_SUN) / (C_LIGHT * C_LIGHT) -} - -/// Schwarzschild radius `r_s = 2 GM/c^2` \[cm\]. -pub fn r_s(m_bh_msun: f64) -> f64 { - 2.0 * gravitational_radius_cm(m_bh_msun) -} +// --- Eddington ------------------------------------------------------------- -/// ISCO radius in units of `R_g = GM/c^2` (dimensionless spin parameter). -/// -/// # Reference -/// Bardeen, Press & Teukolsky 1972, ApJ 178, 347, Eq. (2.21). -pub fn isco_radius(spin: f64) -> f64 { - let a = spin; - let z1 = 1.0 + (1.0 - a * a).cbrt() * ((1.0 + a).cbrt() + (1.0 - a).cbrt()); - let z2 = (3.0 * a * a + z1 * z1).sqrt(); - 3.0 + z2 - ((3.0 - z1) * (3.0 + z1 + 2.0 * z2)).sqrt() -} +pub use eddington::{ + eddington_ratio, l_eddington, luminosity_from_mdot, mdot_at_eddington, mdot_from_luminosity, +}; -/// ISCO radius \[cm\] for mass `m_bh_msun` and spin `spin`. -pub fn r_isco(m_bh_msun: f64, spin: f64) -> f64 { - isco_radius(spin) * gravitational_radius_cm(m_bh_msun) -} +// --- Kerr geometry --------------------------------------------------------- -/// Linear sRGB blackbody color at temperature `temp_k` \[K\], brightest channel -/// normalized to 1.0 (HDR-friendly chromaticity). -/// -/// # References -/// Planck 1901; CIE 1931; Wyman, Sloan & Shirley 2013; IEC 61966-2-1. -pub fn blackbody_rgb(temp_k: f64) -> (f64, f64, f64) { - if temp_k <= 0.0 || !temp_k.is_finite() { - return (0.0, 0.0, 0.0); - } +pub use kerr::{ + THORNE_SPIN_LIMIT, efficiency_from_spin, gravitational_radius_cm, isco_radius, + isco_specific_angular_momentum, isco_specific_energy, orbital_frequency_hz, r_isco, r_s, + specific_angular_momentum, +}; - let (mut x, mut y, mut z) = (0.0_f64, 0.0_f64, 0.0_f64); - let (lambda_min_nm, lambda_max_nm, step_nm) = (380.0, 780.0, 5.0); - let mut lambda_nm = lambda_min_nm; - while lambda_nm <= lambda_max_nm { - let b = planck_spectral_radiance(lambda_nm * NM_TO_CM, temp_k); - let (xb, yb, zb) = cie_xyz_cmf(lambda_nm); - x += b * xb; - y += b * yb; - z += b * zb; - lambda_nm += step_nm; - } +// --- Disk ------------------------------------------------------------------ - let sum = x + y + z; - if sum <= 0.0 { - return (0.0, 0.0, 0.0); - } - let (x, y, z) = (x / sum, y / sum, z / sum); +pub use disk::disk_temperature; - let r = 3.240_625 * x - 1.537_208 * y - 0.498_629 * z; - let g = -0.968_931 * x + 1.875_756 * y + 0.041_518 * z; - let b = 0.055_710 * x - 0.204_021 * y + 1.056_996 * z; +// --- Evolution ------------------------------------------------------------- - let r = r.max(0.0); - let g = g.max(0.0); - let b = b.max(0.0); +pub use evolution::{advance_mass, advance_spin, integrity_rate, salpeter_time_s}; - let peak = r.max(g).max(b); - if peak <= 0.0 { - return (0.0, 0.0, 0.0); - } - (r / peak, g / peak, b / peak) -} +// --- Colorimetry ----------------------------------------------------------- -fn planck_spectral_radiance(lambda_cm: f64, temp_k: f64) -> f64 { - let numerator = 2.0 * H_PLANCK * C_LIGHT * C_LIGHT / lambda_cm.powi(5); - let exponent = H_PLANCK * C_LIGHT / (lambda_cm * K_BOLTZMANN * temp_k); - numerator / (exponent.exp() - 1.0) -} +pub use colorimetry::blackbody_rgb; -fn cie_xyz_cmf(lambda_nm: f64) -> (f64, f64, f64) { - fn g(lambda: f64, mean: f64, sigma_lo: f64, sigma_hi: f64) -> f64 { - let sigma = if lambda < mean { sigma_lo } else { sigma_hi }; - let t = (lambda - mean) / sigma; - (-0.5 * t * t).exp() - } +// --- Derived constants (re-export for golden tests) ------------------------ - let x = 1.056 * g(lambda_nm, 599.8, 37.9, 31.0) + 0.362 * g(lambda_nm, 442.0, 16.0, 26.7) - - 0.065 * g(lambda_nm, 501.1, 20.4, 26.2); - let y = 0.821 * g(lambda_nm, 568.8, 46.9, 40.5) + 0.286 * g(lambda_nm, 530.9, 16.3, 31.1); - let z = 1.217 * g(lambda_nm, 437.0, 11.8, 36.0) + 0.681 * g(lambda_nm, 459.0, 26.0, 13.8); - (x, y, z) -} +pub use derived::{H_BAR, SIGMA_SB, SIGMA_T}; diff --git a/crates/accretion-core/tests/common/mod.rs b/crates/accretion-core/tests/common/mod.rs new file mode 100644 index 0000000..ba681b8 --- /dev/null +++ b/crates/accretion-core/tests/common/mod.rs @@ -0,0 +1,31 @@ +//! Shared helpers for integration tests. + +#![allow(dead_code)] + +/// Relative tolerance for oracle comparison (one part per billion). +pub const ORACLE_RTOL: f64 = 0.000000001; + +/// Tight absolute tolerance for analytic identities. +pub const TOL: f64 = 0.000000001; + +/// Relative tolerance for analytic identities at large magnitudes. +pub const RTOL: f64 = 0.000001; + +/// Standard radiative efficiency used in many identity tests. +pub const ETA: f64 = 0.1; + +/// Assert `got` matches `expected` within relative tolerance `rtol * |expected|`. +pub fn assert_relative_eq(got: f64, expected: f64, rtol: f64) { + assert!( + (got - expected).abs() <= rtol * expected.abs(), + "got={got} expected={expected} rtol={rtol}" + ); +} + +/// Assert `got` matches `expected` within absolute tolerance `tol`. +pub fn assert_abs_eq(got: f64, expected: f64, tol: f64) { + assert!( + (got - expected).abs() <= tol, + "got={got} expected={expected} tol={tol}" + ); +} diff --git a/crates/accretion-core/tests/fixtures/golden.json b/crates/accretion-core/tests/fixtures/golden.json index 8cc9944..42e3a35 100644 --- a/crates/accretion-core/tests/fixtures/golden.json +++ b/crates/accretion-core/tests/fixtures/golden.json @@ -1,6 +1,16 @@ { "source": "astropy 6.0.1", "cases": [ + { + "fn": "sigma_sb", + "args": {}, + "expected": 5.6703744191844314e-05 + }, + { + "fn": "sigma_t", + "args": {}, + "expected": 6.6524587321000005e-25 + }, { "fn": "l_eddington", "args": { diff --git a/crates/accretion-core/tests/golden.rs b/crates/accretion-core/tests/golden.rs index 3cacdd5..a6689cd 100644 --- a/crates/accretion-core/tests/golden.rs +++ b/crates/accretion-core/tests/golden.rs @@ -1,11 +1,12 @@ //! Oracle-driven golden tests: expected values from `scripts/gen_golden.py` (astropy). +mod common; + use std::path::PathBuf; use accretion_core as phys; -/// Relative tolerance for oracle comparison (one part per billion). -const ORACLE_TOL: f64 = 0.000000001; +use common::{ORACLE_RTOL, assert_relative_eq}; #[test] fn golden_matches_astropy_oracle() { @@ -16,6 +17,8 @@ fn golden_matches_astropy_oracle() { for case in doc["cases"].as_array().expect("cases array") { let exp = case["expected"].as_f64().expect("expected f64"); let got = match case["fn"].as_str().expect("fn name") { + "sigma_sb" => phys::SIGMA_SB, + "sigma_t" => phys::SIGMA_T, "l_eddington" => phys::l_eddington(case["args"]["m_bh_msun"].as_f64().unwrap()), "disk_temperature" => phys::disk_temperature( case["args"]["r_cm"].as_f64().unwrap(), @@ -24,10 +27,6 @@ fn golden_matches_astropy_oracle() { ), other => panic!("unknown fn {other}"), }; - assert!( - (got - exp).abs() <= ORACLE_TOL * exp.abs(), - "fn={} got={got} exp={exp}", - case["fn"] - ); + assert_relative_eq(got, exp, ORACLE_RTOL); } } diff --git a/crates/accretion-core/tests/identity.rs b/crates/accretion-core/tests/identity.rs deleted file mode 100644 index 856d0a1..0000000 --- a/crates/accretion-core/tests/identity.rs +++ /dev/null @@ -1,50 +0,0 @@ -//! Analytic physics identities (no oracle required). - -use accretion_core as phys; - -const EFFICIENCY: f64 = 0.1; -const M_MSUN: f64 = 10.0; -const R_IN_R_S: f64 = 10.0; -const R_OUT_R_S: f64 = 20.0; -const TOL: f64 = 0.000000001; -const RTOL: f64 = 0.000001; - -/// Shakura-Sunyaev scaling: `T(2r) / T(r) = 2^(-3/4)`. -/// -/// Reference: Shakura & Sunyaev 1973, A&A 24, 337. -#[test] -fn disk_temperature_scaling_law() { - let mdot = phys::mdot_from_luminosity(phys::l_eddington(M_MSUN), EFFICIENCY); - let r_in = R_IN_R_S * phys::r_s(M_MSUN); - let r_out = R_OUT_R_S * phys::r_s(M_MSUN); - let ratio = - phys::disk_temperature(r_out, M_MSUN, mdot) / phys::disk_temperature(r_in, M_MSUN, mdot); - assert!((ratio - 2.0_f64.powf(-0.75)).abs() < TOL); -} - -/// Schwarzschild ISCO: `r_isco / R_g == 6`. -/// -/// Reference: Bardeen, Press & Teukolsky 1972, ApJ 178, 347. -#[test] -fn isco_schwarzschild_ratio_is_six() { - const M: f64 = 12.3; - let r_g = phys::gravitational_radius_cm(M); - assert!((phys::r_isco(M, 0.0) / r_g - 6.0).abs() < TOL); - assert!((phys::isco_radius(0.0) - 6.0).abs() < TOL); -} - -#[test] -fn mdot_luminosity_round_trip() { - const M: f64 = 100_000_000.0; - let l = phys::l_eddington(M); - let mdot = phys::mdot_from_luminosity(l, EFFICIENCY); - assert!((phys::luminosity_from_mdot(mdot, EFFICIENCY) - l).abs() < RTOL * l.abs()); -} - -#[test] -fn blackbody_color_temperature_ordering() { - let (hr, _hg, hb) = phys::blackbody_rgb(20_000.0); - let (cr, _cg, cb) = phys::blackbody_rgb(3_000.0); - assert!(hb >= hr); - assert!(cr > cb); -} diff --git a/crates/accretion-core/tests/integration.rs b/crates/accretion-core/tests/integration.rs new file mode 100644 index 0000000..dce2d8e --- /dev/null +++ b/crates/accretion-core/tests/integration.rs @@ -0,0 +1,77 @@ +//! Cross-module integration checks (Eddington + Kerr + disk + evolution). + +mod common; + +use accretion_core as phys; + +use common::{ETA, RTOL, TOL, assert_abs_eq, assert_relative_eq}; + +#[test] +fn eddington_luminosity_mass_scaling() { + let l1 = phys::l_eddington(10.0); + let l2 = phys::l_eddington(20.0); + assert_abs_eq(l2 / l1, 2.0, TOL); +} + +#[test] +fn mdot_luminosity_round_trip_at_eddington() { + let l = phys::l_eddington(100_000_000.0); + let mdot = phys::mdot_from_luminosity(l, ETA); + assert_relative_eq(phys::luminosity_from_mdot(mdot, ETA), l, RTOL); +} + +#[test] +fn schwarzschild_isco_chain() { + const M: f64 = 12.3; + let r_g = phys::gravitational_radius_cm(M); + assert_abs_eq(phys::r_isco(M, 0.0) / r_g, 6.0, TOL); + assert_abs_eq(phys::isco_radius(0.0), 6.0, TOL); + assert_abs_eq(phys::r_s(M) / r_g, 2.0, TOL); +} + +#[test] +fn disk_temperature_ss73_scaling() { + const M: f64 = 10.0; + let mdot = phys::mdot_from_luminosity(phys::l_eddington(M), ETA); + let r_in = 10.0 * phys::r_s(M); + let r_out = 20.0 * phys::r_s(M); + let ratio = phys::disk_temperature(r_out, M, mdot) / phys::disk_temperature(r_in, M, mdot); + assert_abs_eq(ratio, 2.0_f64.powf(-0.75), TOL); +} + +#[test] +fn salpeter_time_matches_mass_efolding_at_eddington() { + const M: f64 = 100_000_000.0; + let mdot_edd = phys::mdot_from_luminosity(phys::l_eddington(M), ETA); + let growth_rate_gs = (1.0 - ETA) * mdot_edd; + let efold_s = (M * phys::constants::M_SUN) / growth_rate_gs; + assert_relative_eq(phys::salpeter_time_s(ETA), efold_s, RTOL); +} + +#[test] +fn spin_up_and_efficiency_coupling() { + const M0: f64 = 10.0; + let mdot = phys::mdot_from_luminosity(phys::l_eddington(M0), ETA); + let dt = phys::salpeter_time_s(ETA); + let spin = phys::advance_spin(0.0, M0, mdot, ETA, dt); + let eta = phys::efficiency_from_spin(spin); + assert!(spin > 0.0); + assert!(eta > phys::efficiency_from_spin(0.0)); +} + +#[test] +fn blackbody_colour_ordering() { + let (hr, _hg, hb) = phys::blackbody_rgb(20_000.0); + let (cr, _cg, cb) = phys::blackbody_rgb(3_000.0); + assert!(hb >= hr); + assert!(cr > cb); +} + +#[test] +fn derived_sigma_sb_matches_stefan_boltzmann_formula() { + use phys::constants::{C_LIGHT, H_PLANCK, K_BOLTZMANN}; + const PI: f64 = std::f64::consts::PI; + let expected = + 2.0 * PI.powi(5) * K_BOLTZMANN.powi(4) / (15.0 * C_LIGHT.powi(2) * H_PLANCK.powi(3)); + assert_relative_eq(phys::SIGMA_SB, expected, 0.000000000001); +} diff --git a/crates/godot-ext/src/lib.rs b/crates/godot-ext/src/lib.rs index 845a59c..3ab20d3 100644 --- a/crates/godot-ext/src/lib.rs +++ b/crates/godot-ext/src/lib.rs @@ -22,9 +22,11 @@ struct BlackHole { mdot_gs: f64, #[export] + #[var(set = set_spin)] #[init(val = 0.0)] spin: f64, + /// Kerr radiative efficiency `eta = 1 - E_isco`; kept in sync with `spin`. #[export] #[init(val = 0.1)] efficiency: f64, @@ -41,6 +43,10 @@ struct BlackHole { #[godot_api] impl INode for BlackHole { + fn ready(&mut self) { + self.sync_efficiency_from_spin(); + } + fn process(&mut self, _delta: f64) { let color = self.disk_inner_color(); let uniform = StringName::from(&self.color_uniform); @@ -52,6 +58,18 @@ impl INode for BlackHole { #[godot_api] impl BlackHole { + #[func] + fn set_spin(&mut self, spin: f64) { + self.spin = spin; + self.sync_efficiency_from_spin(); + } + + /// Radiative efficiency implied by the current spin (Kerr ISCO binding energy). + #[func] + fn efficiency_from_spin(&self) -> f64 { + phys::efficiency_from_spin(self.spin) + } + #[func] fn disk_inner_temp(&self, r_cm: f64) -> f64 { phys::disk_temperature(r_cm, self.mass_solar, self.mdot_gs) @@ -87,10 +105,72 @@ impl BlackHole { phys::isco_radius(self.spin) } + /// New mass \[M_sun\] after accreting at the current feed rate for `dt_s` \[s\]. + #[func] + fn advance_mass(&self, dt_s: f64) -> f64 { + phys::advance_mass(self.mass_solar, self.mdot_gs, self.efficiency, dt_s) + } + + /// Salpeter (Eddington) e-folding time \[s\] at the current efficiency. + #[func] + fn salpeter_time_s(&self) -> f64 { + phys::salpeter_time_s(self.efficiency) + } + + /// Disk-support fraction `1 - lambda`: >0 stable/recovering, <0 disrupting. + #[func] + fn integrity_rate(&self) -> f64 { + phys::integrity_rate(self.eddington_ratio()) + } + + /// Accretion rate \[g/s\] at the Eddington limit (`lambda = 1`) for current mass. + #[func] + fn mdot_at_eddington(&self) -> f64 { + phys::mdot_at_eddington(self.mass_solar, self.efficiency) + } + + /// Kerr spin after accreting at the current feed for `dt_s` \[s\]. + #[func] + fn advance_spin(&self, dt_s: f64) -> f64 { + phys::advance_spin( + self.spin, + self.mass_solar, + self.mdot_gs, + self.efficiency, + dt_s, + ) + } + #[func] fn disk_inner_color(&self) -> Color { let t = phys::disk_temperature(self.inner_radius_cm(), self.mass_solar, self.mdot_gs); let (r, g, b) = phys::blackbody_rgb(t); Color::from_rgba(r as f32, g as f32, b as f32, 1.0) } + + /// Blackbody color of the disk at `r_over_risco` ISCO radii out (clamped to the + /// inner edge), so the presentation can build a physical radial gradient. + #[func] + fn disk_color_at(&self, r_over_risco: f64) -> Color { + let r = self.inner_radius_cm() * r_over_risco.max(1.0); + let t = phys::disk_temperature(r, self.mass_solar, self.mdot_gs); + let (red, g, b) = phys::blackbody_rgb(t); + Color::from_rgba(red as f32, g as f32, b as f32, 1.0) + } + + /// Disk inner-edge temperature \[K\] (Shakura-Sunyaev at the ISCO). + #[func] + fn disk_inner_temp_k(&self) -> f64 { + phys::disk_temperature(self.inner_radius_cm(), self.mass_solar, self.mdot_gs) + } + + /// Prograde orbital frequency \[Hz\] at the ISCO (high-frequency QPO scale). + #[func] + fn isco_orbital_frequency_hz(&self) -> f64 { + phys::orbital_frequency_hz(self.mass_solar, phys::isco_radius(self.spin), self.spin) + } + + fn sync_efficiency_from_spin(&mut self) { + self.efficiency = phys::efficiency_from_spin(self.spin); + } } diff --git a/scenes/GodotSmoke.tscn b/scenes/GodotSmoke.tscn new file mode 100644 index 0000000..f76b776 --- /dev/null +++ b/scenes/GodotSmoke.tscn @@ -0,0 +1,6 @@ +[gd_scene format=3 uid="uid://csxs8qpfy65aa"] + +[ext_resource type="Script" path="res://scripts/godot_smoke.gd" id="1"] + +[node name="GodotSmoke" type="Node"] +script = ExtResource("1") diff --git a/scenes/Main.tscn b/scenes/Main.tscn index 121898a..5a07b36 100644 --- a/scenes/Main.tscn +++ b/scenes/Main.tscn @@ -1,17 +1,28 @@ -[gd_scene load_steps=9 format=3] +[gd_scene load_steps=11 format=3] [ext_resource type="Shader" path="res://shaders/blackhole.gdshader" id="1_shader"] [ext_resource type="Script" path="res://scripts/main.gd" id="2_main"] [ext_resource type="Shader" path="res://shaders/starfield_sky.gdshader" id="3_sky"] +[sub_resource type="FastNoiseLite" id="FastNoiseLite_disc"] +noise_type = 1 +frequency = 0.015 + +[sub_resource type="NoiseTexture2D" id="NoiseTexture2D_disc"] +seamless = true +noise = SubResource("FastNoiseLite_disc") + [sub_resource type="ShaderMaterial" id="ShaderMaterial_disc"] shader = ExtResource("1_shader") +shader_parameter/disc_texture = SubResource("NoiseTexture2D_disc") shader_parameter/disc_outer_radius = 0.9 shader_parameter/disc_inner_radius = 0.2 shader_parameter/disc_speed = 1.5 shader_parameter/disc_color = Color(0.8, 0.18, 0.05, 1.0) shader_parameter/inner_color = Color(1.0, 1.0, 0.95, 1.0) shader_parameter/doppler_beaming_factor = 66.0 +shader_parameter/emission_strength = 1.0 +shader_parameter/qpo_phase_rate = 1.5 shader_parameter/hue_radius = 0.75 shader_parameter/hue_shift_factor = -0.03 shader_parameter/streak_intensity = 0.5 @@ -71,8 +82,8 @@ color_uniform = "inner_color" [node name="Panel" type="Panel" parent="UI"] offset_left = 12.0 offset_top = 12.0 -offset_right = 420.0 -offset_bottom = 420.0 +offset_right = 440.0 +offset_bottom = 480.0 [node name="Controls" type="VBoxContainer" parent="UI/Panel"] anchor_right = 1.0 @@ -87,6 +98,15 @@ theme_override_constants/separation = 6 theme_override_font_sizes/font_size = 18 text = "Accretion" +[node name="ClassLabel" type="Label" parent="UI/Panel/Controls"] +theme_override_colors/font_color = Color(0.6, 0.85, 1, 1) +text = "Stellar-mass black hole" + +[node name="ScoreLabel" type="Label" parent="UI/Panel/Controls"] +theme_override_colors/font_color = Color(1, 0.9, 0.55, 1) +theme_override_font_sizes/font_size = 12 +text = "Score 0" + [node name="MassLabel" type="Label" parent="UI/Panel/Controls"] text = "Mass (log₁₀ M☉)" @@ -101,7 +121,7 @@ text = "Accretion rate (log₁₀ g/s)" [node name="FeedSlider" type="HSlider" parent="UI/Panel/Controls"] min_value = 14.0 -max_value = 20.0 +max_value = 28.0 step = 0.01 value = 18.0 @@ -114,6 +134,23 @@ max_value = 0.998 step = 0.002 value = 0.0 +[node name="MilestoneLabel" type="Label" parent="UI/Panel/Controls"] +theme_override_font_sizes/font_size = 12 +autowrap_mode = 3 +text = "Next milestone" + +[node name="MilestoneBar" type="ProgressBar" parent="UI/Panel/Controls"] +custom_minimum_size = Vector2(0, 10) +show_percentage = false + +[node name="IntegrityLabel" type="Label" parent="UI/Panel/Controls"] +theme_override_font_sizes/font_size = 12 +text = "Disk integrity" + +[node name="IntegrityBar" type="ProgressBar" parent="UI/Panel/Controls"] +custom_minimum_size = Vector2(0, 14) +value = 100.0 + [node name="StatsLabel" type="RichTextLabel" parent="UI/Panel/Controls"] custom_minimum_size = Vector2(0, 200) bbcode_enabled = true @@ -135,3 +172,57 @@ theme_override_colors/font_color = Color(0.7, 0.75, 0.85, 1) theme_override_font_sizes/font_size = 11 autowrap_mode = 3 text = "controls" + +[node name="Banner" type="Label" parent="UI"] +visible = false +anchors_preset = 8 +anchor_left = 0.5 +anchor_top = 0.5 +anchor_right = 0.5 +anchor_bottom = 0.5 +offset_left = -300.0 +offset_top = -60.0 +offset_right = 300.0 +offset_bottom = 60.0 +grow_horizontal = 2 +grow_vertical = 2 +theme_override_colors/font_color = Color(1, 0.85, 0.4, 1) +theme_override_font_sizes/font_size = 22 +horizontal_alignment = 1 +vertical_alignment = 1 +autowrap_mode = 3 +text = "BANNER" + +[node name="RunOverlay" type="PanelContainer" parent="UI"] +visible = false +anchors_preset = 8 +anchor_left = 0.5 +anchor_top = 0.5 +anchor_right = 0.5 +anchor_bottom = 0.5 +offset_left = -260.0 +offset_top = -150.0 +offset_right = 260.0 +offset_bottom = 150.0 +grow_horizontal = 2 +grow_vertical = 2 + +[node name="Margin" type="MarginContainer" parent="UI/RunOverlay"] +theme_override_constants/margin_left = 16 +theme_override_constants/margin_top = 16 +theme_override_constants/margin_right = 16 +theme_override_constants/margin_bottom = 16 + +[node name="VBox" type="VBoxContainer" parent="UI/RunOverlay/Margin"] +theme_override_constants/separation = 10 + +[node name="RunTitle" type="Label" parent="UI/RunOverlay/Margin/VBox"] +theme_override_font_sizes/font_size = 24 +horizontal_alignment = 1 +text = "RUN OVER" + +[node name="RunBody" type="RichTextLabel" parent="UI/RunOverlay/Margin/VBox"] +custom_minimum_size = Vector2(0, 180) +bbcode_enabled = true +fit_content = true +text = "..." diff --git a/scripts/build_godot_ext.sh b/scripts/build_godot_ext.sh new file mode 100755 index 0000000..be617bb --- /dev/null +++ b/scripts/build_godot_ext.sh @@ -0,0 +1,24 @@ +#!/usr/bin/env sh +# Build libgodot_ext and copy into bin/ (res://bin/...) for Godot to load. +# Cursor/sandbox may set CARGO_TARGET_DIR elsewhere; pin output to ./target. +set -eu +cd "$(dirname "$0")/.." +export CARGO_TARGET_DIR="${CARGO_TARGET_DIR:-$(pwd)/target}" +cargo build -p godot-ext "$@" +mkdir -p bin +case "$(uname -s)" in + Darwin) + install -m 755 target/debug/libgodot_ext.dylib bin/libgodot_ext.dylib + ;; + Linux) + install -m 755 target/debug/libgodot_ext.so bin/libgodot_ext.so + ;; + MINGW* | MSYS* | CYGWIN*) + install -m 755 target/debug/godot_ext.dll bin/godot_ext.dll + ;; + *) + echo "Unsupported OS for godot-ext install: $(uname -s)" >&2 + exit 1 + ;; +esac +echo "Installed godot-ext → bin/ ($(uname -s))" diff --git a/scripts/gen_constants.py b/scripts/gen_constants.py index b687509..e95cac6 100644 --- a/scripts/gen_constants.py +++ b/scripts/gen_constants.py @@ -1,27 +1,34 @@ #!/usr/bin/env python3 """Generate crates/accretion-core/src/constants.rs from astropy (CGS). -Single source of truth. Re-run to update; commit the diff. Do not hand-edit the output. +Single source of truth for the FUNDAMENTAL constants. Composite constants +(Stefan-Boltzmann sigma_sb, Thomson cross-section sigma_T) are NOT emitted here; +they are derived from these fundamentals in `accretion-core::derived`, where the +mathematical relation is exact and explicit. Re-run to update; commit the diff. +Do not hand-edit the output. """ from __future__ import annotations import astropy import astropy.constants as c +# Only FUNDAMENTAL constants. Anything that is an exact mathematical function of +# these (sigma_sb, sigma_T) is derived in Rust, not tabulated here. vals = { "G": (float(c.G.cgs.value), "cm^3 g^-1 s^-2", "astropy.constants.G"), "C_LIGHT": (float(c.c.cgs.value), "cm s^-1", "astropy.constants.c"), "M_SUN": (float(c.M_sun.cgs.value), "g", "astropy.constants.M_sun"), "M_P": (float(c.m_p.cgs.value), "g", "astropy.constants.m_p"), - "SIGMA_T": (float(c.sigma_T.cgs.value), "cm^2", "astropy.constants.sigma_T"), - "SIGMA_SB": (float(c.sigma_sb.cgs.value), "erg cm^-2 s^-1 K^-4", "astropy.constants.sigma_sb"), + "M_E": (float(c.m_e.cgs.value), "g", "astropy.constants.m_e"), "H_PLANCK": (float(c.h.cgs.value), "erg s", "astropy.constants.h"), "K_BOLTZMANN": (float(c.k_B.cgs.value), "erg K^-1", "astropy.constants.k_B"), + "ALPHA": (float(c.alpha.value), "dimensionless", "astropy.constants.alpha"), } out = [ f"// GENERATED by scripts/gen_constants.py from astropy {astropy.__version__}.\n", - "// Do not edit by hand. Single source of truth for physical constants (CGS).\n\n", + "// Do not edit by hand. Single source of truth for FUNDAMENTAL constants (CGS).\n", + "// Composite constants (sigma_sb, sigma_T) are derived in `crate::derived`.\n\n", ] for name, (val, unit, src) in vals.items(): out.append(f"/// {src} [{unit}]\n") @@ -31,5 +38,7 @@ out.append("/// Nanometres to centimetres [cm nm^-1] (exact: 1 nm = 1e-7 cm).\n") out.append("pub const NM_TO_CM: f64 = 1e-7;\n\n") -open("crates/accretion-core/src/constants.rs", "w").write("".join(out)) +# Single trailing newline (rustfmt-clean), internal spacing preserved. +text = "".join(out).rstrip() + "\n" +open("crates/accretion-core/src/constants.rs", "w").write(text) print("wrote constants.rs from astropy", astropy.__version__) diff --git a/scripts/gen_golden.py b/scripts/gen_golden.py index f299692..4ca4f4d 100644 --- a/scripts/gen_golden.py +++ b/scripts/gen_golden.py @@ -33,6 +33,15 @@ def disk_t(r: float, m: float, md: float) -> float: cases: list[dict] = [] + +# Derived constants: pin the Rust derivation against astropy's tabulated value. +# sigma_sb is exact (2019 SI); sigma_T is consistent to CODATA's internal +# precision (the Rust value uses alpha; astropy uses e/eps0/m_e). +cases.append( + {"fn": "sigma_sb", "args": {}, "expected": float(c.sigma_sb.cgs.value)} +) +cases.append({"fn": "sigma_t", "args": {}, "expected": float(c.sigma_T.cgs.value)}) + for m in (10.0, 1.0e6, 6.5e9): cases.append({"fn": "l_eddington", "args": {"m_bh_msun": m}, "expected": l_edd(m)}) diff --git a/scripts/godot_smoke.gd b/scripts/godot_smoke.gd new file mode 100644 index 0000000..1509fba --- /dev/null +++ b/scripts/godot_smoke.gd @@ -0,0 +1,93 @@ +extends Node +## Headless GDExtension smoke test — every BlackHole API used by main.gd. +## +## Run: make godot-smoke (or godot --headless --path . res://scenes/GodotSmoke.tscn) + +const GRADIENT_OUTER_RISCO := 14.0 + + +func _ready() -> void: + # GDExtension may register one frame after project boot on some hosts. + for _i in 30: + if ClassDB.class_exists("BlackHole"): + break + await get_tree().process_frame + var exit_code := _run() + get_tree().quit(exit_code) + + +func _run() -> int: + if not ClassDB.class_exists("BlackHole"): + push_error("BlackHole GDExtension class not registered — run `make build` first") + return 1 + + var bh: Node = ClassDB.instantiate("BlackHole") + if bh == null: + push_error("Failed to instantiate BlackHole") + return 1 + add_child(bh) + + bh.set("mass_solar", 21.0) + bh.set("mdot_gs", 8.0e17) + bh.set("spin", 0.0) + + var methods := [ + "l_eddington", + "luminosity_erg_s", + "eddington_ratio", + "isco_in_rg", + "schwarzschild_radius_cm", + "disk_inner_temp_k", + "isco_orbital_frequency_hz", + "salpeter_time_s", + "mdot_at_eddington", + "disk_inner_color", + "integrity_rate", + ] + for method in methods: + if not bh.has_method(method): + push_error("Missing BlackHole method: %s (stale native lib? run `make build`)" % method) + return 1 + + var t_salpeter: float = bh.call("salpeter_time_s") + if not is_finite(t_salpeter) or t_salpeter <= 0.0: + push_error("salpeter_time_s returned invalid value: %s" % t_salpeter) + return 1 + + var dt := t_salpeter * 0.01 + var new_mass: float = bh.call("advance_mass", dt) + if not is_finite(new_mass) or new_mass <= 0.0: + push_error("advance_mass returned invalid value: %s" % new_mass) + return 1 + bh.set("mass_solar", new_mass) + + var new_spin: float = bh.call("advance_spin", dt) + if not is_finite(new_spin): + push_error("advance_spin returned invalid value: %s" % new_spin) + return 1 + bh.set("spin", new_spin) + + var eta: float = bh.get("efficiency") + if not is_finite(eta) or eta <= 0.0 or eta >= 1.0: + push_error("efficiency out of range after spin sync: %s" % eta) + return 1 + + var rate: float = bh.call("integrity_rate") + if not is_finite(rate): + push_error("integrity_rate returned NaN") + return 1 + + var inner: Color = bh.call("disk_color_at", 1.0) + var outer: Color = bh.call("disk_color_at", GRADIENT_OUTER_RISCO) + if inner.r + inner.g + inner.b <= 0.0: + push_error("disk_color_at(1.0) returned black") + return 1 + + print( + "[godot_smoke] OK mass=", new_mass, + " spin=", new_spin, + " eta=", eta, + " salpeter_s=", t_salpeter, + " integrity=", rate + ) + return 0 diff --git a/scripts/godot_smoke.gd.uid b/scripts/godot_smoke.gd.uid new file mode 100644 index 0000000..9fd9110 --- /dev/null +++ b/scripts/godot_smoke.gd.uid @@ -0,0 +1 @@ +uid://1pynocpvhj7c diff --git a/scripts/hud_format.gd b/scripts/hud_format.gd index 49d8a1b..6dc431c 100644 --- a/scripts/hud_format.gd +++ b/scripts/hud_format.gd @@ -42,6 +42,29 @@ static func lambda_edd(lam: float) -> String: return "%.3f× Edd" % lam +static func years(seconds: float) -> String: + var yr := seconds / 31557600.0 # Julian year + if yr >= 1.0e9: + return "%.2f Gyr" % (yr / 1.0e9) + if yr >= 1.0e6: + return "%.2f Myr" % (yr / 1.0e6) + if yr >= 1.0e3: + return "%.2f kyr" % (yr / 1.0e3) + return "%.0f yr" % yr + + +static func hertz(f: float) -> String: + if f >= 1.0e3: + return "%.2f kHz" % (f / 1.0e3) + if f >= 1.0: + return "%.2f Hz" % f + if f >= 1.0e-3: + return "%.2f mHz" % (f * 1.0e3) + if f >= 1.0e-6: + return "%.2f µHz" % (f * 1.0e6) + return "%s Hz" % sci(f, 2) + + static func sci(x: float, sig: int) -> String: if x == 0.0: return "0" diff --git a/scripts/main.gd b/scripts/main.gd index 728b340..45ed8a5 100644 --- a/scripts/main.gd +++ b/scripts/main.gd @@ -1,63 +1,129 @@ extends Node3D -## Game controller: input, camera orbit, and telemetry HUD. +## Game controller: accretion survival loop, scoring, and presentation. ## -## Presentation only (rule 10). Every physical quantity displayed here is read -## from the BlackHole gdext node, which delegates to accretion-core in Rust. -## This script formats and logs; it never computes physics. +## Presentation only (rule 10). Physics comes from the BlackHole gdext node +## (accretion-core in Rust). Gameplay tuning constants live here (rule 02 §6). -const PRESET_CYGX1 := {"name": "Cyg X-1", "mass": 21.0, "mdot": 1.0e18} -const PRESET_SGRA := {"name": "Sgr A*", "mass": 4.0e6, "mdot": 1.0e15} -const PRESET_M87 := {"name": "M87*", "mass": 6.5e9, "mdot": 1.0e19} +enum RunPhase { PLAYING, ENDED } + +const PRESET_CYGX1 := {"name": "Cyg X-1", "mass": 21.0, "mdot": 8.0e17, "spin": 0.0} +const PRESET_SGRA := {"name": "Sgr A*", "mass": 4.0e6, "mdot": 2.0e23, "spin": 0.0} +const PRESET_M87 := {"name": "M87*", "mass": 6.5e9, "mdot": 2.0e26, "spin": 0.85} + +const MILESTONES := [ + {"mass": 21.0, "name": "Cyg X-1 — stellar-mass"}, + {"mass": 1.0e3, "name": "Intermediate-mass (IMBH)"}, + {"mass": 1.0e5, "name": "Seed supermassive"}, + {"mass": 4.0e6, "name": "Sgr A* — galactic centre"}, + {"mass": 1.0e8, "name": "Supermassive"}, + {"mass": 6.5e9, "name": "M87* — giant elliptical"}, +] +const VICTORY_MASS := 6.5e9 +const MAX_DISRUPTIONS := 3 +const START_PRESET := PRESET_CYGX1 + +const SECONDS_PER_EFOLD := 22.0 +const DRAIN_SCALE := 0.55 +const RECOVER_SCALE := 0.16 +const DISRUPT_RESET := 0.35 + +const IDLE_BEFORE_AUTORBIT := 4.0 +const AUTORBIT_SPEED := 0.07 +const CAM_SMOOTH := 7.0 +const INTRO_DOLLY_FROM := 30.0 +const DISK_BASE_SCALE := 4.0 +const OUTER_RG := 26.0 +const GRADIENT_OUTER_RISCO := 14.0 @onready var black_hole: Node = $BlackHole @onready var camera: Camera3D = $Camera3D +@onready var disk_mesh: Node3D = $DiskMesh @onready var mass_slider: HSlider = $UI/Panel/Controls/MassSlider @onready var feed_slider: HSlider = $UI/Panel/Controls/FeedSlider @onready var spin_slider: HSlider = $UI/Panel/Controls/SpinSlider @onready var title_label: Label = $UI/Panel/Controls/TitleLabel +@onready var class_label: Label = $UI/Panel/Controls/ClassLabel +@onready var score_label: Label = $UI/Panel/Controls/ScoreLabel +@onready var milestone_label: Label = $UI/Panel/Controls/MilestoneLabel +@onready var milestone_bar: ProgressBar = $UI/Panel/Controls/MilestoneBar +@onready var integrity_bar: ProgressBar = $UI/Panel/Controls/IntegrityBar @onready var stats_label: RichTextLabel = $UI/Panel/Controls/StatsLabel @onready var warning_label: Label = $UI/Panel/Controls/WarningLabel @onready var swatch: ColorRect = $UI/Panel/Controls/Swatch @onready var help_label: Label = $UI/Panel/Controls/HelpLabel +@onready var banner: Label = $UI/Banner +@onready var run_overlay: PanelContainer = $UI/RunOverlay +@onready var run_title: Label = $UI/RunOverlay/Margin/VBox/RunTitle +@onready var run_body: RichTextLabel = $UI/RunOverlay/Margin/VBox/RunBody -var _cam_dist := 9.0 -var _cam_yaw := 0.0 -var _cam_pitch := 0.12 -var _dragging := false +var _phase := RunPhase.PLAYING +var _mass := 10.0 +var _integrity := 1.0 +var _sim_seconds := 0.0 var _preset_name := "custom" +var _class_name := "" var _was_super_eddington := false +var _setting_slider := false +var _banner_time := 0.0 +var _disruptions := 0 +var _milestones_hit := 0 +var _milestone_idx := 0 +var _run_score := 0 +var _best_score := 0 +var _victory := false + +var _cam_dist := 9.0 +var _cam_yaw := 0.4 +var _cam_pitch := 0.16 +var _target_dist := 9.0 +var _target_yaw := 0.4 +var _target_pitch := 0.16 +var _dragging := false +var _idle_time := 0.0 func _ready() -> void: - camera.look_at(Vector3.ZERO, Vector3.UP) - mass_slider.value_changed.connect(_on_control_changed) - feed_slider.value_changed.connect(_on_control_changed) - spin_slider.value_changed.connect(_on_control_changed) - _apply_sliders() - _update_camera() - print("[accretion] ready — presets 1=Cyg X-1, 2=Sgr A*, 3=M87*; push feed past Eddington with Z/X") + mass_slider.value_changed.connect(_on_mass_changed) + feed_slider.value_changed.connect(_on_feed_changed) + spin_slider.value_changed.connect(_on_spin_changed) + integrity_bar.min_value = 0.0 + integrity_bar.max_value = 1.0 + milestone_bar.min_value = 0.0 + milestone_bar.max_value = 1.0 + banner.visible = false + run_overlay.visible = false + _best_score = RunScore.load_best() + _start_run(START_PRESET) + _cam_dist = INTRO_DOLLY_FROM + _update_camera_transform() + print("[accretion] Grow to M87* without losing the disk three times. R = new run · 1/2/3 = challenge presets") func _input(event: InputEvent) -> void: + if _phase == RunPhase.ENDED: + return if event is InputEventMouseButton: if event.button_index == MOUSE_BUTTON_LEFT: _dragging = event.pressed + if event.pressed and event.button_index == MOUSE_BUTTON_WHEEL_UP: + _target_dist = max(_target_dist - 0.7, 4.0) + elif event.pressed and event.button_index == MOUSE_BUTTON_WHEEL_DOWN: + _target_dist = min(_target_dist + 0.7, 20.0) + _idle_time = 0.0 if event is InputEventMouseMotion and _dragging: - _cam_yaw -= event.relative.x * 0.005 - _cam_pitch = clamp(_cam_pitch - event.relative.y * 0.004, -0.35, 0.55) - _update_camera() - if event is InputEventMouseButton and event.pressed: - if event.button_index == MOUSE_BUTTON_WHEEL_UP: - _cam_dist = max(_cam_dist - 0.6, 4.0) - _update_camera() - elif event.button_index == MOUSE_BUTTON_WHEEL_DOWN: - _cam_dist = min(_cam_dist + 0.6, 18.0) - _update_camera() + _target_yaw -= event.relative.x * 0.005 + _target_pitch = clamp(_target_pitch - event.relative.y * 0.004, -0.35, 0.6) + _idle_time = 0.0 func _unhandled_input(event: InputEvent) -> void: if not event is InputEventKey or not event.pressed: return + if event.keycode == KEY_R or event.keycode == KEY_ENTER or event.keycode == KEY_KP_ENTER: + _start_run(_current_preset_dict()) + return + if _phase == RunPhase.ENDED: + return match event.keycode: KEY_Q, KEY_COMMA: mass_slider.value = max(mass_slider.min_value, mass_slider.value - 0.05) @@ -72,40 +138,200 @@ func _unhandled_input(event: InputEvent) -> void: KEY_D: spin_slider.value = min(spin_slider.max_value, spin_slider.value + 0.02) KEY_1: - _apply_preset(PRESET_CYGX1) + _start_run(PRESET_CYGX1) KEY_2: - _apply_preset(PRESET_SGRA) + _start_run(PRESET_SGRA) KEY_3: - _apply_preset(PRESET_M87) + _start_run(PRESET_M87) -func _on_control_changed(_v: float) -> void: - _preset_name = "custom" - _apply_sliders() +func _current_preset_dict() -> Dictionary: + match _preset_name: + PRESET_CYGX1.name: + return PRESET_CYGX1 + PRESET_SGRA.name: + return PRESET_SGRA + PRESET_M87.name: + return PRESET_M87 + _: + return { + "name": "custom", + "mass": _mass, + "mdot": black_hole.get("mdot_gs"), + "spin": black_hole.get("spin"), + } -func _apply_preset(p: Dictionary) -> void: +func _start_run(p: Dictionary) -> void: + _phase = RunPhase.PLAYING + run_overlay.visible = false _preset_name = p.name - mass_slider.value = log(p.mass) / log(10.0) - feed_slider.value = log(p.mdot) / log(10.0) - print("[accretion] preset → %s (M=%s, Ṁ=%s g/s)" % [ - p.name, HudFormat.mass_msun(p.mass), HudFormat.grams_per_s(p.mdot) - ]) - _apply_sliders() + _mass = p.mass + _integrity = 1.0 + _sim_seconds = 0.0 + _disruptions = 0 + _milestones_hit = 0 + _milestone_idx = 0 + _class_name = "" + _victory = false + _was_super_eddington = false + _set_slider_silent(mass_slider, log(_mass) / log(10.0)) + _set_slider_silent(feed_slider, log(p.mdot) / log(10.0)) + _set_slider_silent(spin_slider, p.get("spin", 0.0)) + _apply_inputs() + _update_class() + _show_banner("RUN START — feed the disk, ride λ ≤ 1, reach M87*", 2.5) + print("[accretion] new run → %s" % p.name) + + +func _on_mass_changed(v: float) -> void: + if _setting_slider or _phase != RunPhase.PLAYING: + return + _preset_name = "custom" + _mass = pow(10.0, v) + _apply_inputs() + + +func _on_feed_changed(_v: float) -> void: + if _setting_slider or _phase != RunPhase.PLAYING: + return + _preset_name = "custom" + _apply_inputs() -func _apply_sliders() -> void: - black_hole.set("mass_solar", pow(10.0, mass_slider.value)) +func _on_spin_changed(_v: float) -> void: + if _setting_slider or _phase != RunPhase.PLAYING: + return + _apply_inputs() + + +func _apply_inputs() -> void: + black_hole.set("mass_solar", _mass) black_hole.set("mdot_gs", pow(10.0, feed_slider.value)) black_hole.set("spin", spin_slider.value) - _refresh_hud() -func _process(_delta: float) -> void: +func _process(delta: float) -> void: + if _phase == RunPhase.PLAYING: + _step_simulation(delta) + _update_camera(delta) _refresh_hud() + _drive_graphics() + if _banner_time > 0.0: + _banner_time -= delta + if _banner_time <= 0.0: + banner.visible = false + + +func _step_simulation(delta: float) -> void: + var t_salpeter: float = black_hole.call("salpeter_time_s") + var dt_sim := delta * t_salpeter / SECONDS_PER_EFOLD + _sim_seconds += dt_sim + + _mass = black_hole.call("advance_mass", dt_sim) + black_hole.set("mass_solar", _mass) + _set_slider_silent(mass_slider, clamp(log(_mass) / log(10.0), mass_slider.min_value, mass_slider.max_value)) + + var spin: float = black_hole.call("advance_spin", dt_sim) + black_hole.set("spin", spin) + _set_slider_silent(spin_slider, spin) + + var rate: float = black_hole.call("integrity_rate") + if rate < 0.0: + _integrity += rate * DRAIN_SCALE * delta + else: + _integrity += min(rate, 1.0) * RECOVER_SCALE * delta + _integrity = clamp(_integrity, 0.0, 1.0) + if _integrity <= 0.0: + _on_disruption() + + _check_milestones() + _update_class() + + if _mass >= VICTORY_MASS: + _end_run(true) + +func _check_milestones() -> void: + while _milestone_idx < MILESTONES.size() and _mass >= MILESTONES[_milestone_idx].mass: + _milestones_hit += 1 + _show_banner("MILESTONE — %s" % MILESTONES[_milestone_idx].name, 2.5) + _milestone_idx += 1 -func _update_camera() -> void: + +func _on_disruption() -> void: + _disruptions += 1 + _integrity = DISRUPT_RESET + _set_slider_silent(feed_slider, feed_slider.min_value) + black_hole.set("mdot_gs", pow(10.0, feed_slider.min_value)) + push_warning("[accretion] disk disrupted (%d/%d)" % [_disruptions, MAX_DISRUPTIONS]) + if _disruptions >= MAX_DISRUPTIONS: + _show_banner("DISK LOST — three super-Eddington blowouts", 2.0) + _end_run(false) + else: + _show_banner( + "DISK DISRUPTED (%d/%d)\nFeed reset — rebuild integrity before pushing λ past 1 again." + % [_disruptions, MAX_DISRUPTIONS], + 3.5 + ) + + +func _end_run(victory: bool) -> void: + _phase = RunPhase.ENDED + _victory = victory + _run_score = RunScore.compute(_mass, _milestones_hit, _disruptions, _sim_seconds / 31557600.0) + var new_best := RunScore.save_best(_run_score) + if new_best: + _best_score = _run_score + run_overlay.visible = true + run_title.text = "VICTORY — M87* reached" if victory else "RUN OVER — disk destroyed" + run_body.text = ( + "[b]Score[/b] [color=yellow]%d[/color] best [color=cyan]%d[/color]\n" + % [_run_score, _best_score] + + " mass [color=cyan]%s[/color]\n" % HudFormat.mass_msun(_mass) + + " milestones %d · disruptions %d\n" % [_milestones_hit, _disruptions] + + " sim time %s\n\n" % HudFormat.years(_sim_seconds) + + "[color=gray]R or Enter — new run · 1/2/3 — challenge preset[/color]" + ) + + +func _update_class() -> void: + var cls := "Stellar-mass black hole" + if _mass >= 1.0e8: + cls = "Supermassive black hole" + elif _mass >= 1.0e5: + cls = "Seed supermassive black hole" + elif _mass >= 1.0e2: + cls = "Intermediate-mass black hole" + if cls != _class_name: + if _class_name != "" and _phase == RunPhase.PLAYING: + _show_banner("RECLASSIFIED — %s" % cls, 2.0) + _class_name = cls + + +func _next_milestone() -> Dictionary: + for m in MILESTONES: + if m.mass > _mass: + return m + return {"mass": MILESTONES[-1].mass, "name": "%s (reached)" % MILESTONES[-1].name} + + +func _live_score() -> int: + return RunScore.compute(_mass, _milestones_hit, _disruptions, _sim_seconds / 31557600.0) + + +func _update_camera(delta: float) -> void: + _idle_time += delta + if _idle_time > IDLE_BEFORE_AUTORBIT and not _dragging: + _target_yaw += AUTORBIT_SPEED * delta + var k := 1.0 - exp(-delta * CAM_SMOOTH) + _cam_yaw = lerp(_cam_yaw, _target_yaw, k) + _cam_pitch = lerp(_cam_pitch, _target_pitch, k) + _cam_dist = lerp(_cam_dist, _target_dist, k) + _update_camera_transform() + + +func _update_camera_transform() -> void: var offset := Vector3( _cam_dist * cos(_cam_pitch) * sin(_cam_yaw), _cam_dist * sin(_cam_pitch) + 2.0, @@ -115,55 +341,105 @@ func _update_camera() -> void: camera.look_at(Vector3.ZERO, Vector3.UP) +func _drive_graphics() -> void: + var mat: ShaderMaterial = black_hole.get("disk_material") + if mat == null: + return + var spin: float = black_hole.get("spin") + var isco_rg: float = black_hole.call("isco_in_rg") + var horizon_rg := 1.0 + sqrt(max(1.0 - spin * spin, 0.0)) + var f_qpo: float = black_hole.call("isco_orbital_frequency_hz") + + mat.set_shader_parameter("inner_color", black_hole.call("disk_color_at", 1.0)) + mat.set_shader_parameter("disc_color", black_hole.call("disk_color_at", GRADIENT_OUTER_RISCO)) + mat.set_shader_parameter("disc_inner_radius", clamp(isco_rg / OUTER_RG, 0.05, 0.6)) + mat.set_shader_parameter("ss_radius", clamp(horizon_rg / OUTER_RG, 0.02, 0.3)) + + var t_in: float = black_hole.call("disk_inner_temp_k") + var log_t: float = clamp(log(max(t_in, 1.0)) / log(10.0), 4.0, 7.5) + mat.set_shader_parameter("emission_strength", remap(log_t, 4.0, 7.5, 0.6, 3.2)) + mat.set_shader_parameter("disc_speed", clamp(log(max(f_qpo, 1.0e-12)) / log(10.0) + 2.5, 0.8, 4.0)) + mat.set_shader_parameter("qpo_phase_rate", clamp(f_qpo * 0.02, 0.4, 6.0)) + + disk_mesh.scale = Vector3.ONE * DISK_BASE_SCALE * (1.0 + 0.12 * _class_index()) + + +func _class_index() -> int: + if _mass >= 1.0e8: + return 3 + if _mass >= 1.0e5: + return 2 + if _mass >= 1.0e2: + return 1 + return 0 + + func _refresh_hud() -> void: var m: float = black_hole.get("mass_solar") var mdot: float = black_hole.get("mdot_gs") var spin: float = black_hole.get("spin") + var eta: float = black_hole.get("efficiency") var l_edd: float = black_hole.call("l_eddington") var l_bol: float = black_hole.call("luminosity_erg_s") var lam: float = black_hole.call("eddington_ratio") - var r_in: float = black_hole.call("inner_radius_cm") - var r_s: float = black_hole.call("schwarzschild_radius_cm") var isco_rg: float = black_hole.call("isco_in_rg") - var t_in: float = black_hole.call("disk_inner_temp", r_in) + var r_s: float = black_hole.call("schwarzschild_radius_cm") + var t_in: float = black_hole.call("disk_inner_temp_k") + var f_qpo: float = black_hole.call("isco_orbital_frequency_hz") + var t_salpeter: float = black_hole.call("salpeter_time_s") + var mdot_edd: float = black_hole.call("mdot_at_eddington") var c: Color = black_hole.call("disk_inner_color") swatch.color = c title_label.text = "Accretion — %s" % _preset_name + class_label.text = _class_name + score_label.text = "Score %d · best %d · disruptions %d/%d" % [ + _live_score(), _best_score, _disruptions, MAX_DISRUPTIONS + ] + + var next: Dictionary = _next_milestone() + milestone_label.text = "Next: %s (%s)" % [next.name, HudFormat.mass_msun(next.mass)] + milestone_bar.value = clamp(log(_mass) / log(next.mass), 0.0, 1.0) + + integrity_bar.value = _integrity + integrity_bar.modulate = Color(1.0, 0.4, 0.25) if _integrity < 0.35 else Color(0.4, 0.9, 1.0) var super_edd := lam > 1.0 - if super_edd and not _was_super_eddington: - push_warning("[accretion] super-Eddington: λ=%.2f — loss condition (future gameplay)" % lam) - elif not super_edd and _was_super_eddington: - print("[accretion] sub-Eddington again: λ=%.3f" % lam) + if super_edd and not _was_super_eddington and _phase == RunPhase.PLAYING: + push_warning("[accretion] super-Eddington: λ=%.2f" % lam) _was_super_eddington = super_edd - warning_label.visible = super_edd - warning_label.text = ( - "SUPER-EDDINGTON — λ = %s\n" - % HudFormat.lambda_edd(lam) - + "Radiation pressure exceeds gravity; the disk would be blown apart." - ) - warning_label.modulate = Color(1.0, 0.35, 0.2) if super_edd else Color.WHITE + warning_label.visible = super_edd and _phase == RunPhase.PLAYING + warning_label.text = "SUPER-EDDINGTON — λ = %s\nIntegrity draining." % HudFormat.lambda_edd(lam) stats_label.text = ( - "[b]Black hole[/b]\n" + "[b]Black hole[/b] [color=gray]sim t = %s[/color]\n" % HudFormat.years(_sim_seconds) + " mass M = [color=cyan]%s[/color]\n" % HudFormat.mass_msun(m) - + " spin a/M = %.2f ISCO = %.2f R_g (%.1f r_s)\n" % [spin, isco_rg, isco_rg / 2.0] - + " Schwarzschild radius r_s = %s\n" % HudFormat.cm(r_s) - + " ISCO radius r_in = %s\n\n" % HudFormat.cm(r_in) + + " spin a/M = %.3f → η = [color=violet]%.1f%%[/color]\n" % [spin, eta * 100.0] + + " ISCO = %.2f R_g r_s = %s\n" % [isco_rg, HudFormat.cm(r_s)] + + " e-fold @ λ=1 = %s\n\n" % HudFormat.years(t_salpeter) + "[b]Accretion[/b]\n" - + " accretion rate Ṁ = [color=yellow]%s[/color]\n" % HudFormat.grams_per_s(mdot) - + " bolometric luminosity L = %s\n" % HudFormat.erg_per_s(l_bol) - + " Eddington limit L_Edd = %s\n" % HudFormat.erg_per_s(l_edd) - + " Eddington ratio λ = [color=%s]%s[/color]\n\n" % [ - "red" if super_edd else "lightgreen", HudFormat.lambda_edd(lam) - ] - + "[b]Inner disk edge[/b] [color=gray](computed in Rust)[/color]\n" - + " temperature T = [color=orange]%s[/color]\n" % HudFormat.kelvin(t_in) - + " blackbody sRGB = (%.2f, %.2f, %.2f)" % [c.r, c.g, c.b] + + " feed Ṁ = [color=yellow]%s[/color]\n" % HudFormat.grams_per_s(mdot) + + " λ-safe Ṁ = [color=lightgreen]%s[/color] (λ = 1 ceiling)\n" % HudFormat.grams_per_s(mdot_edd) + + " L = %s L_Edd = %s\n" % [HudFormat.erg_per_s(l_bol), HudFormat.erg_per_s(l_edd)] + + " λ = [color=%s]%s[/color]\n\n" % ["red" if super_edd else "lightgreen", HudFormat.lambda_edd(lam)] + + "[b]Inner edge[/b]\n" + + " T = [color=orange]%s[/color] QPO f = %s\n" % [HudFormat.kelvin(t_in), HudFormat.hertz(f_qpo)] ) - help_label.text = ( - "Q/E · mass Z/X · feed A/D · spin 1/2/3 · presets drag · orbit scroll · zoom" - ) + if _phase == RunPhase.ENDED: + help_label.text = "R or Enter — new run · 1/2/3 — challenge preset" + else: + help_label.text = "Q/E mass · Z/X feed · A/D spin · 1/2/3 presets · R restart · drag · scroll" + + +func _set_slider_silent(slider: HSlider, value: float) -> void: + _setting_slider = true + slider.value = value + _setting_slider = false + + +func _show_banner(text: String, seconds: float) -> void: + banner.text = text + banner.visible = true + _banner_time = seconds diff --git a/scripts/run_score.gd b/scripts/run_score.gd new file mode 100644 index 0000000..8c87427 --- /dev/null +++ b/scripts/run_score.gd @@ -0,0 +1,33 @@ +class_name RunScore +extends RefCounted +## Run scoring and persisted high score (presentation only). + +const SAVE_PATH := "user://accretion_best.run" + +static func compute(mass_msun: float, milestones: int, disruptions: int, sim_years: float) -> int: + var mass_pts := int(max(log(max(mass_msun, 1.0)) / log(10.0), 0.0) * 1200.0) + var milestone_pts := milestones * 800 + var survival_pts := int(min(sim_years, 1.0e12) / 1.0e6) # 1 pt per Myr, capped + var penalty := disruptions * 350 + return max(mass_pts + milestone_pts + survival_pts - penalty, 0) + + +static func load_best() -> int: + if not FileAccess.file_exists(SAVE_PATH): + return 0 + var f := FileAccess.open(SAVE_PATH, FileAccess.READ) + if f == null: + return 0 + var text := f.get_as_text().strip_edges() + return int(text) if text.is_valid_int() else 0 + + +static func save_best(score: int) -> bool: + var best := load_best() + if score <= best: + return false + var f := FileAccess.open(SAVE_PATH, FileAccess.WRITE) + if f == null: + return false + f.store_string(str(score)) + return true diff --git a/scripts/run_score.gd.uid b/scripts/run_score.gd.uid new file mode 100644 index 0000000..c77677c --- /dev/null +++ b/scripts/run_score.gd.uid @@ -0,0 +1 @@ +uid://32k7ik1d452 diff --git a/shaders/blackhole.gdshader b/shaders/blackhole.gdshader index 73e7b9e..b9eec8d 100644 --- a/shaders/blackhole.gdshader +++ b/shaders/blackhole.gdshader @@ -51,6 +51,11 @@ uniform vec4 disc_color : source_color = vec4(1.0, 0.0, 0.0, 1.0); // Driven from Rust (accretion_core::blackbody_rgb of the inner-edge temperature). uniform vec4 inner_color : source_color = vec4(1.0, 1.0, 0.95, 1.0); uniform float doppler_beaming_factor : hint_range(0.0, 200.0) = 66.0; +// Bounded HDR bloom factor; driven from Rust as a tone-map of the inner-edge +// temperature, so hotter disks visibly bloom brighter (presentation only). +uniform float emission_strength : hint_range(0.1, 6.0) = 1.0; +// ISCO orbital phase rate (Hz, presentation-scaled from Rust QPO frequency). +uniform float qpo_phase_rate : hint_range(0.0, 10.0) = 1.5; uniform float hue_radius : hint_range(0.0, 1.0) = 0.75; uniform float hue_shift_factor : hint_range(-1.0, 1.0) = -0.03; uniform float streak_intensity : hint_range(0.0, 1.0) = 0.5; @@ -251,6 +256,10 @@ vec3 calc_disc_color(vec3 base_color, vec3 planar_disc_pos, vec3 disc_dir, vec3 float turb_intensity_mod = (combined_turb - 0.5) * 2.0 * turbulence_intensity; float hot_spot = pow(max(0.0, combined_turb - 0.6) * 2.5, 2.0) * turbulence_intensity; + // Keplerian hotspot orbiting at the ISCO (QPO-scale flicker; phase rate from Rust). + float qpo_angle = uv.y * PI * 2.0 + TIME * qpo_phase_rate; + float qpo_spot = pow(max(0.0, cos(qpo_angle) * 0.5 + 0.5), 6.0) * (1.0 - radial_t); + hot_spot += qpo_spot * 0.35; // inner_color is the Rust-computed blackbody color of the disc inner edge. vec3 temp_color = mix(inner_color.rgb, base_color, pow(radial_t, 0.6)); @@ -403,6 +412,7 @@ void fragment() { } disc_hit *= tex_col * disc_color.a * streaks * horizon_fade * radial_falloff; + final_disc_col *= emission_strength; vec3 col = mix(background_col, final_disc_col, disc_hit); ALBEDO = col;