diff --git a/.github/dependabot.yml b/.github/dependabot.yml
new file mode 100644
index 0000000..0c2a409
--- /dev/null
+++ b/.github/dependabot.yml
@@ -0,0 +1,19 @@
+version: 2
+updates:
+  - package-ecosystem: github-actions
+    directory: /
+    schedule:
+      interval: weekly
+    open-pull-requests-limit: 5
+
+  - package-ecosystem: pip
+    directory: /python
+    schedule:
+      interval: weekly
+    open-pull-requests-limit: 5
+
+  - package-ecosystem: cargo
+    directory: /native/core
+    schedule:
+      interval: weekly
+    open-pull-requests-limit: 5
diff --git a/.github/workflows/arch-report.yml b/.github/workflows/arch-report.yml
deleted file mode 100644
index a2b5776..0000000
--- a/.github/workflows/arch-report.yml
+++ /dev/null
@@ -1,34 +0,0 @@
-name: Architecture Report
-
-on:
-  push:
-    branches: [main]
-  workflow_dispatch:
-  pull_request:
-
-jobs:
-  report:
-    name: arch-${{ matrix.os }}
-    runs-on: ${{ matrix.os }}
-    strategy:
-      fail-fast: false
-      matrix:
-        include:
-          - os: ubuntu-24.04
-          - os: ubuntu-24.04-arm
-          - os: windows-2025
-          - os: windows-11-arm
-          - os: macos-15-intel
-          - os: macos-15
-    defaults:
-      run:
-        working-directory: python
-    steps:
-      - uses: actions/checkout@v4
-      - uses: actions/setup-python@v5
-        with:
-          python-version: "3.12"
-      - name: Install package
-        run: python -m pip install -U pip && python -m pip install -e .[dev]
-      - name: Doctor
-        run: blackhole-accelerators doctor --json --fail-on-emulation
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
new file mode 100644
index 0000000..1ed5ab9
--- /dev/null
+++ b/.github/workflows/ci.yml
@@ -0,0 +1,210 @@
+name: CI
+
+on:
+  pull_request:
+  push:
+    branches: [main]
+  workflow_dispatch:
+
+permissions:
+  contents: read
+
+concurrency:
+  group: ci-${{ github.workflow }}-${{ github.ref }}
+  cancel-in-progress: true
+
+jobs:
+  hygiene:
+    name: hygiene
+    runs-on: ubuntu-24.04
+    timeout-minutes: 10
+    steps:
+      - uses: actions/checkout@v4
+      - uses: actions/setup-python@v5
+        with:
+          python-version: "3.12"
+      - name: Install YAML parser
+        run: python -m pip install -U pip pyyaml
+      - name: Validate workflow YAML and repository artifact policy
+        run: |
+          python - <<'PY'
+          from pathlib import Path
+          import sys
+          import yaml
+
+          for path in sorted(Path(".github/workflows").glob("*.yml")):
+              with path.open("r", encoding="utf-8") as f:
+                  yaml.safe_load(f)
+
+          import subprocess
+
+          blocked_suffixes = {".npz", ".h5", ".hdf5", ".png", ".tiff", ".whl"}
+          tracked = subprocess.check_output(["git", "ls-files"], text=True).splitlines()
+          blocked = [p for p in tracked if Path(p).suffix.lower() in blocked_suffixes]
+          if blocked:
+              print("Generated or binary artifacts must not be committed:")
+              print("\n".join(sorted(blocked)))
+              sys.exit(1)
+          PY
+      - name: Git whitespace check
+        run: git diff --check
+
+  python:
+    name: python-${{ matrix.os }}-${{ matrix.python-version }}
+    needs: hygiene
+    runs-on: ${{ matrix.os }}
+    timeout-minutes: 30
+    strategy:
+      fail-fast: false
+      matrix:
+        include:
+          - os: ubuntu-24.04
+            python-version: "3.11"
+          - os: ubuntu-24.04
+            python-version: "3.12"
+          - os: ubuntu-24.04
+            python-version: "3.13"
+          - os: windows-2025
+            python-version: "3.12"
+          - os: macos-15
+            python-version: "3.12"
+    defaults:
+      run:
+        working-directory: python
+    steps:
+      - uses: actions/checkout@v4
+      - uses: actions/setup-python@v5
+        with:
+          python-version: ${{ matrix.python-version }}
+          cache: pip
+          cache-dependency-path: python/pyproject.toml
+      - name: Install package
+        run: python -m pip install -U pip && python -m pip install -e ".[dev,hdf5]"
+      - name: Compile
+        run: python -m compileall blackhole_sim
+      - name: Test
+        run: python -m pytest -q
+      - name: Accelerator doctor
+        run: python -m blackhole_sim.accelerator_cli doctor --json --fail-on-emulation
+
+  science-contracts:
+    name: science-contracts
+    needs: hygiene
+    runs-on: ubuntu-24.04
+    timeout-minutes: 25
+    defaults:
+      run:
+        working-directory: python
+    steps:
+      - uses: actions/checkout@v4
+      - uses: actions/setup-python@v5
+        with:
+          python-version: "3.12"
+          cache: pip
+          cache-dependency-path: python/pyproject.toml
+      - name: Install package
+        run: python -m pip install -U pip && python -m pip install -e ".[dev,hdf5]"
+      - name: Focused physics contracts
+        run: >
+          python -m pytest -q
+          tests/test_kerr.py
+          tests/test_radiative_transfer.py
+          tests/test_synchrotron_polarized_transfer.py
+          tests/test_public_dumps.py
+          tests/test_external_validation.py
+          tests/test_native_kernel_assets.py
+      - name: Benchmark parity smoke
+        run: python -m blackhole_sim.benchmark_cli --json --nr 3 --ntheta 3 --nphi 4 --points 7 --iterations 1
+      - name: Accelerated renderer smoke
+        run: python -m blackhole_sim.accelerated_cli --width 32 --height 18 --max-steps 64 --output out/ci_stokes_smoke.npz
+
+  rust:
+    name: rust-${{ matrix.os }}
+    needs: hygiene
+    runs-on: ${{ matrix.os }}
+    timeout-minutes: 20
+    strategy:
+      fail-fast: false
+      matrix:
+        os: [ubuntu-24.04, windows-2025, macos-15]
+    steps:
+      - uses: actions/checkout@v4
+      - uses: dtolnay/rust-toolchain@stable
+      - name: Format
+        run: cargo fmt --check --manifest-path native/core/Cargo.toml
+      - name: Test
+        run: cargo test --manifest-path native/core/Cargo.toml
+
+  native-wheel:
+    name: native-wheel-${{ matrix.os }}
+    needs: [python, rust, science-contracts]
+    runs-on: ${{ matrix.os }}
+    timeout-minutes: 45
+    strategy:
+      fail-fast: false
+      matrix:
+        include:
+          - os: ubuntu-24.04
+          - os: ubuntu-24.04-arm
+          - os: windows-2025
+          - os: windows-11-arm
+          - os: macos-15-intel
+          - os: macos-15
+    steps:
+      - uses: actions/checkout@v4
+      - uses: actions/setup-python@v5
+        with:
+          python-version: "3.12"
+          cache: pip
+          cache-dependency-path: python/pyproject.toml
+      - uses: dtolnay/rust-toolchain@stable
+      - name: Install maturin
+        run: python -m pip install -U pip maturin
+      - name: Build native wheel
+        run: maturin build --manifest-path native/core/Cargo.toml --release --out dist --features extension-module
+      - name: Install simulator and native wheel
+        shell: python
+        run: |
+          from pathlib import Path
+          import subprocess
+          import sys
+
+          subprocess.check_call([sys.executable, "-m", "pip", "install", "-e", "./python[dev,hdf5]"])
+          wheels = sorted(Path("dist").glob("blackhole_native-*.whl"))
+          assert wheels, "no blackhole_native wheel produced"
+          subprocess.check_call([sys.executable, "-m", "pip", "install", "--force-reinstall", str(wheels[-1])])
+      - name: Native wheel doctor
+        run: |
+          python -c "import blackhole_native, json; info=blackhole_native.detect_arch(); print(json.dumps({'version': blackhole_native.core_version(), 'arch': info}, sort_keys=True)); assert info.get('arch') in {'arm64', 'x86_64', 'x86'}"
+          python -m blackhole_sim.accelerator_cli doctor --json --fail-on-emulation
+      - name: Native parity regression
+        run: python -m pytest -q python/tests/test_native_stokes_parity.py python/tests/test_native_sampler_parity.py python/tests/test_benchmark.py
+      - name: Native benchmark smoke
+        run: python -m blackhole_sim.benchmark_cli --json --nr 3 --ntheta 3 --nphi 4 --points 7 --iterations 1
+      - uses: actions/upload-artifact@v4
+        with:
+          name: native-wheel-${{ matrix.os }}
+          path: dist/*
+          if-no-files-found: error
+          retention-days: 14
+
+  required:
+    name: ci-required
+    if: always()
+    needs: [hygiene, python, science-contracts, rust, native-wheel]
+    runs-on: ubuntu-24.04
+    timeout-minutes: 5
+    steps:
+      - name: Check required jobs
+        run: |
+          python - <<'PY'
+          import json
+          import sys
+
+          results = json.loads(r'''${{ toJson(needs) }}''')
+          failed = {name: data["result"] for name, data in results.items() if data["result"] != "success"}
+          if failed:
+              print(failed)
+              sys.exit(1)
+          print("all required CI jobs passed")
+          PY
diff --git a/.github/workflows/native.yml b/.github/workflows/native.yml
deleted file mode 100644
index a84410f..0000000
--- a/.github/workflows/native.yml
+++ /dev/null
@@ -1,63 +0,0 @@
-name: Native Core
-
-on:
-  push:
-    branches: [main]
-  pull_request:
-
-jobs:
-  rust:
-    name: rust-${{ matrix.os }}
-    runs-on: ${{ matrix.os }}
-    strategy:
-      fail-fast: false
-      matrix:
-        os: [ubuntu-24.04, windows-2025, macos-15]
-    steps:
-      - uses: actions/checkout@v4
-      - uses: actions/setup-python@v5
-        with:
-          python-version: "3.12"
-      - uses: dtolnay/rust-toolchain@stable
-      - name: Format
-        run: cargo fmt --check --manifest-path native/core/Cargo.toml
-      - name: Test
-        run: cargo test --manifest-path native/core/Cargo.toml --no-default-features
-
-  maturin-smoke:
-    name: maturin-${{ matrix.os }}
-    runs-on: ${{ matrix.os }}
-    strategy:
-      fail-fast: false
-      matrix:
-        include:
-          - os: ubuntu-24.04
-          - os: ubuntu-24.04-arm
-          - os: windows-2025
-          - os: windows-11-arm
-          - os: macos-15-intel
-          - os: macos-15
-    steps:
-      - uses: actions/checkout@v4
-      - uses: actions/setup-python@v5
-        with:
-          python-version: "3.12"
-      - uses: dtolnay/rust-toolchain@stable
-      - name: Install maturin
-        run: python -m pip install -U pip maturin
-      - name: Build native wheel
-        run: maturin build --manifest-path native/core/Cargo.toml --release --out dist --features extension-module
-      - name: Install simulator and native wheel
-        run: |
-          python -m pip install -e "./python[dev]"
-          python -c "from pathlib import Path; import subprocess, sys; wheels=sorted(Path('dist').glob('blackhole_native-*.whl')); assert wheels, 'no blackhole_native wheel produced'; subprocess.check_call([sys.executable, '-m', 'pip', 'install', '--force-reinstall', str(wheels[-1])])"
-      - name: Native wheel doctor
-        run: |
-          python -c "import blackhole_native, json; info=blackhole_native.detect_arch(); print(json.dumps({'version': blackhole_native.core_version(), 'arch': info}, sort_keys=True)); assert info.get('arch') in {'arm64', 'x86_64', 'x86'}"
-          blackhole-accelerators doctor --json --fail-on-emulation
-      - name: Native parity regression
-        run: python -m pytest -q python/tests/test_native_stokes_parity.py python/tests/test_native_sampler_parity.py
-      - uses: actions/upload-artifact@v4
-        with:
-          name: native-wheel-${{ matrix.os }}
-          path: dist/*
diff --git a/.github/workflows/python.yml b/.github/workflows/python.yml
deleted file mode 100644
index ddaa05d..0000000
--- a/.github/workflows/python.yml
+++ /dev/null
@@ -1,32 +0,0 @@
-name: Python
-
-on:
-  push:
-    branches: [main]
-  pull_request:
-
-jobs:
-  tests:
-    name: python-${{ matrix.os }}-${{ matrix.python-version }}
-    runs-on: ${{ matrix.os }}
-    strategy:
-      fail-fast: false
-      matrix:
-        os: [ubuntu-24.04, windows-2025, macos-15]
-        python-version: ["3.11", "3.12"]
-    defaults:
-      run:
-        working-directory: python
-    steps:
-      - uses: actions/checkout@v4
-      - uses: actions/setup-python@v5
-        with:
-          python-version: ${{ matrix.python-version }}
-      - name: Install package
-        run: python -m pip install -U pip && python -m pip install -e .[dev]
-      - name: Compile
-        run: python -m compileall blackhole_sim
-      - name: Test
-        run: python -m pytest -q
-      - name: Accelerator doctor
-        run: blackhole-accelerators doctor --json --fail-on-emulation
diff --git a/.github/workflows/release-artifacts.yml b/.github/workflows/release-artifacts.yml
new file mode 100644
index 0000000..6774fd6
--- /dev/null
+++ b/.github/workflows/release-artifacts.yml
@@ -0,0 +1,109 @@
+name: Release Artifacts
+
+on:
+  workflow_dispatch:
+  push:
+    tags:
+      - "v*"
+
+permissions:
+  contents: read
+
+concurrency:
+  group: release-artifacts-${{ github.ref }}
+  cancel-in-progress: false
+
+jobs:
+  python-package:
+    name: python-package
+    runs-on: ubuntu-24.04
+    timeout-minutes: 20
+    steps:
+      - uses: actions/checkout@v4
+      - uses: actions/setup-python@v5
+        with:
+          python-version: "3.12"
+          cache: pip
+          cache-dependency-path: python/pyproject.toml
+      - name: Build Python distribution
+        working-directory: python
+        run: |
+          python -m pip install -U pip build
+          python -m build --sdist --wheel --outdir ../dist/python
+      - name: Smoke install wheel
+        run: |
+          python -m venv .venv-smoke
+          . .venv-smoke/bin/activate
+          python -m pip install -U pip
+          python -m pip install dist/python/blackhole_sim-*.whl
+          python -m blackhole_sim.accelerator_cli doctor --json --fail-on-emulation
+      - uses: actions/upload-artifact@v4
+        with:
+          name: python-package
+          path: dist/python/*
+          if-no-files-found: error
+          retention-days: 30
+
+  native-wheels:
+    name: native-wheel-${{ matrix.os }}
+    runs-on: ${{ matrix.os }}
+    timeout-minutes: 45
+    strategy:
+      fail-fast: false
+      matrix:
+        include:
+          - os: ubuntu-24.04
+          - os: ubuntu-24.04-arm
+          - os: windows-2025
+          - os: windows-11-arm
+          - os: macos-15-intel
+          - os: macos-15
+    steps:
+      - uses: actions/checkout@v4
+      - uses: actions/setup-python@v5
+        with:
+          python-version: "3.12"
+      - uses: dtolnay/rust-toolchain@stable
+      - name: Install maturin
+        run: python -m pip install -U pip maturin
+      - name: Build native wheel
+        run: maturin build --manifest-path native/core/Cargo.toml --release --out dist/native --features extension-module
+      - name: Smoke install native wheel
+        shell: python
+        run: |
+          from pathlib import Path
+          import subprocess
+          import sys
+
+          subprocess.check_call([sys.executable, "-m", "pip", "install", "-e", "./python[dev]"])
+          wheels = sorted(Path("dist/native").glob("blackhole_native-*.whl"))
+          assert wheels, "no native wheel produced"
+          subprocess.check_call([sys.executable, "-m", "pip", "install", "--force-reinstall", str(wheels[-1])])
+          subprocess.check_call([sys.executable, "-m", "blackhole_sim.accelerator_cli", "doctor", "--json", "--fail-on-emulation"])
+      - uses: actions/upload-artifact@v4
+        with:
+          name: native-wheel-${{ matrix.os }}
+          path: dist/native/*
+          if-no-files-found: error
+          retention-days: 30
+
+  artifact-summary:
+    name: artifact-summary
+    if: always()
+    needs: [python-package, native-wheels]
+    runs-on: ubuntu-24.04
+    timeout-minutes: 5
+    steps:
+      - name: Check artifact jobs
+        run: |
+          python - <<'PY'
+          import json
+          import sys
+
+          results = json.loads(r'''${{ toJson(needs) }}''')
+          failed = {name: data["result"] for name, data in results.items() if data["result"] != "success"}
+          if failed:
+              print(failed)
+              sys.exit(1)
+          print("release artifacts built and uploaded; no package or GitHub release was published")
+          PY
diff --git a/README.md b/README.md
index 809cf5d..5e7aaf5 100644
--- a/README.md
+++ b/README.md
@@ -11,6 +11,24 @@ web/          WebGPU renderer and guided browser showcase
 docs/         build state, native roadmap, data validation, performance notes
 ```
 
+## Physics Contract
+
+Python is the authoritative physics reference path. The repo models Kerr
+spacetime, local ZAMO camera rays, null geodesics, GRMHD snapshot sampling,
+invariant redshift hooks, and polarized Stokes transfer, but it is still a
+research reference implementation. WebGPU and native kernels are acceleration
+targets and must pass deterministic parity gates before any physics-equivalence
+claim is made.
+
+The unpolarized preview path defaults to `educational_proxy` mode. `validated`
+mode requires metric-aware local plasma calculations plus a non-proxy
+coefficient model or external baseline evidence.
+
+See `docs/PHYSICS_VALIDATION.md` for the current grounding, approximation
+boundaries, validation commands, and the WebGPU-versus-game-engine policy. See
+`docs/SCIENTIFIC_ROADMAP.md` for the release gates that turn the project into a
+validated GRRT platform.
+
 ## Install
 
 ```bash
@@ -72,7 +90,10 @@ for local adapter diagnostics.
 The CUDA, Metal, HIP, OpenCL, and WGSL kernel assets share the same staged
 contract: invalid nonperiodic brick samples are guarded before Stokes transfer
 updates, and no GPU physics parity claim is made until a small deterministic
-render passes against the CPU reference envelope on the target hardware.
+render passes against the CPU reference envelope on the target hardware. The
+Stokes kernel source now derives camera rays from a ZAMO launch helper rather
+than fixed photon momenta, but readback parity is still required before claiming
+GPU renderer parity.
 
 ## Showcase
 
diff --git a/ci/README.md b/ci/README.md
index 436f2ff..6a53681 100644
--- a/ci/README.md
+++ b/ci/README.md
@@ -1,7 +1,30 @@
-# CI Notes
+# CI/CD Notes
 
-CI is defined in `.github/workflows/`.
+CI/CD is defined in `.github/workflows/`.
 
-The current workflows are smoke gates for the native-foundation milestone. They upload build artifacts where useful but do not publish packages, create releases, sign apps, or notarize macOS bundles.
+## Pull Request And Main CI
 
-Native wheel jobs must install the wheel they just built and run `blackhole-accelerators doctor --json --fail-on-emulation` on the target runner before uploading artifacts.
+`ci.yml` is the required validation pipeline for pull requests and pushes to
+`main`. It runs:
+
+- repository hygiene and workflow YAML parsing,
+- Python compile/test/accelerator-doctor gates across Linux, Windows, and macOS,
+- focused science-contract tests plus benchmark and accelerated-render smoke,
+- Rust format and native-core tests,
+- native wheel builds on x64/ARM Linux, Windows, and macOS,
+- native wheel install, architecture doctor, parity tests, and benchmark smoke.
+
+The final `ci-required` job is the branch-protection target. It fails if any
+upstream validation job fails or is skipped.
+
+## Artifact Delivery
+
+`release-artifacts.yml` runs manually or on `v*` tags. It builds Python package
+artifacts and native wheels, smoke-installs them, and uploads artifacts for
+review. It does not publish to PyPI, create GitHub releases, sign binaries, or
+notarize macOS bundles.
+
+## Dependency Maintenance
+
+`.github/dependabot.yml` opens weekly update PRs for GitHub Actions, Python
+package metadata, and the Rust native core.
diff --git a/docs/BUILD_STATE.md b/docs/BUILD_STATE.md
index e80e690..741593f 100644
--- a/docs/BUILD_STATE.md
+++ b/docs/BUILD_STATE.md
@@ -1,6 +1,6 @@
 # Build State
 
-Date: 2026-06-29
+Date: 2026-07-01
 
 ## Current Milestone
 
@@ -36,6 +36,11 @@ Date: 2026-06-29
 - WebGPU accelerator discovery now reports OS-visible GPU adapters on Windows so direct browser GPU execution can be audited alongside shader availability.
 - The WebGPU browser renderer now reports the granted adapter and labels compute/fragment paths explicitly.
 - CUDA, Metal, HIP, OpenCL, and WebGPU Stokes kernel assets now guard invalid nonperiodic coefficient-brick samples before applying Stokes transfer updates.
+- Added a science-grade validation direction: explicit `educational_proxy` versus `validated` transfer mode, metric-aware local plasma magnetic-field and pitch-angle helpers, and a release-gate roadmap for GRRT validation.
+- Coefficient-brick precomputation and polarized transfer now use the shared metric-aware local plasma helper when snapshot spin is available.
+- WebGPU Stokes compute and static CUDA/Metal/HIP/OpenCL kernel assets now derive launch state and conserved photon momenta through a ZAMO camera helper instead of fixed `p_t`/`p_phi` constants.
+- Installed seven local Codex science skills under `C:\Users\mkang\.codex\skills` for Kerr validation, GRMHD data validation, GRRT transfer, numerical methods, GPU parity, scientific Python quality, and scientific visualization review.
+- Replaced separate Python/native/architecture smoke workflows with a unified CI pipeline plus a manual/tag artifact-delivery workflow and Dependabot update configuration. The artifact workflow uploads review artifacts only; it does not publish packages or releases.
 
 ## Release Boundary
 
@@ -142,6 +147,29 @@ v0.9.0 direct GPU contract local evidence:
 - `maturin build --manifest-path native/core/Cargo.toml --release --out native/core/target/wheels-gpu-contract`: passed and produced `blackhole_native-0.9.0-cp310-abi3-win_arm64.whl`.
 - In-app browser smoke against `http://127.0.0.1:8800/index.html?shader=stokes`: loaded the default WebGPU page, status reported `WebGPU direct compute: Stokes coefficient bricks on qualcomm / adreno-7xx`, and console warnings/errors were empty. Follow-up automation for the opt-in `diagnostics=1` readback path timed out before DOM readback, so GPU readback verification is not claimed here.
 
+v0.9.0 science-grade redirection local evidence:
+
+- `python -m pytest -q tests/test_synchrotron_polarized_transfer.py tests/test_radiative_transfer.py tests/test_native_kernel_assets.py tests/test_kerr.py`: passed.
+- Local Codex skill validation with `quick_validate.py`: passed for `blackhole-gr-kerr-validation`, `blackhole-grmhd-data-validation`, `blackhole-grrt-polarized-transfer`, `blackhole-numerical-methods`, `blackhole-hpc-gpu-parity`, `scientific-python-quality`, and `scientific-visualization-review`.
+- `python -m compileall blackhole_sim`: passed.
+- `python -m pytest -q`: passed with 3 optional skips.
+- `cargo fmt --check --manifest-path native/core/Cargo.toml`: passed.
+- `cargo test --manifest-path native/core/Cargo.toml`: passed; 6 Rust unit tests passed.
+- `python -m blackhole_sim.accelerator_cli doctor --json --fail-on-emulation`: passed; `process_arch=arm64`, `python_arch=arm64`, `native_core_loaded=true`, `native_core_arch=arm64`, `native_core_version=0.9.0`, `gpu_backend=webgpu`, and `emulation_detected=false`.
+- `python -m blackhole_sim.benchmark_cli --json --nr 3 --ntheta 3 --nphi 4 --points 7 --iterations 1`: passed; native sampler and Stokes RK2 parity both reported `allclose=true`, `max_abs_diff=0.0`, and `max_rel_diff=0.0`.
+- `python -m blackhole_sim.accelerated_cli --width 32 --height 18 --max-steps 64 --output out\science_redirection_smoke.npz`: passed; output is ignored and not committed.
+- Direct WebGPU browser readback and concrete Illinois v3 selected-dump/ipole validation were not run in this pass, so GPU image parity and external GRRT parity are still not claimed.
+
+v0.9.0 CI/CD pipeline local evidence:
+
+- `.github/workflows/ci.yml` added as the pull-request and `main` validation pipeline, including hygiene, Python, science-contract, Rust, native-wheel, and final required aggregation jobs.
+- `.github/workflows/release-artifacts.yml` added for manual or `v*` tag artifact builds; it uploads artifacts but intentionally does not publish PyPI packages or GitHub releases.
+- `.github/dependabot.yml` added for weekly GitHub Actions, Python, and Cargo update PRs.
+- Local workflow YAML parsing with Python/PyYAML: passed for `.github/workflows/ci.yml`, `.github/workflows/release-artifacts.yml`, and `.github/dependabot.yml`.
+- Tracked-artifact policy check using `git ls-files`: passed; no tracked `.npz`, `.h5`, `.hdf5`, `.png`, `.tiff`, or `.whl` artifacts.
+- `git diff --check`: passed; Windows line-ending warnings are non-blocking.
+- `actionlint` was not available on PATH, so GitHub Actions execution is the next validation layer after push.
+
 ## GitHub CI Evidence
 
 Validated on commit `46f20630f41556dce75cb6142a54ff816185c54a`:
diff --git a/docs/PHYSICS_VALIDATION.md b/docs/PHYSICS_VALIDATION.md
new file mode 100644
index 0000000..232a766
--- /dev/null
+++ b/docs/PHYSICS_VALIDATION.md
@@ -0,0 +1,117 @@
+# Physics Validation Contract
+
+## Purpose
+
+This repository should be treated as a research reference implementation, not as
+a production EHT/GRRT pipeline. The code is grounded in known physics where the
+model and tests explicitly say so, and it must not claim scientific accuracy for
+paths that are only visual, approximate, or parity-unproven.
+
+## Physics Anchors
+
+The authoritative physics path is the Python reference implementation.
+
+- Kerr spacetime is represented in Boyer-Lindquist coordinates with
+  geometrized units `G = c = M = 1`.
+- Photon rays are integrated as null geodesics with Hamiltonian equations,
+  `H = 0.5 g^{mu nu} p_mu p_nu = 0`.
+- Camera rays are launched from a local ZAMO orthonormal tetrad so screen-space
+  directions map to local physical directions before converting to coordinate
+  momenta.
+- The code includes Kerr horizon, static-limit, ISCO, Keplerian orbit, circular
+  four-velocity, and redshift helpers.
+- GRMHD snapshots are represented as explicit `r, theta, phi` grids with
+  density, electron temperature, pressure, magnetic-field four-vector, and
+  fluid four-velocity fields.
+- Radiative transfer uses invariant redshift hooks and Stokes `I, Q, U, V`
+  transport with exact matrix stepping when SciPy is available and RK2 fallback
+  otherwise.
+- Validated plasma-frame helpers compute magnetic-field magnitude with the Kerr
+  metric invariant `sqrt(b_mu b^mu)` when snapshot spin is available, and compute
+  photon/B pitch angle from the fluid-frame invariant
+  `(p_mu b^mu)/(E_fluid |B|)`.
+
+## Current Validation Gates
+
+The minimum physics checks are:
+
+```bash
+cd python
+python -m pytest tests/test_geodesics.py tests/test_kerr.py tests/test_grmhd.py tests/test_radiative_transfer.py tests/test_synchrotron_polarized_transfer.py -q
+```
+
+These cover:
+
+- Schwarzschild critical impact parameter and capture/escape behavior.
+- Kerr metric inverse identity and Schwarzschild-limit metric components.
+- Local ZAMO tetrad orthonormality.
+- Null camera-ray launch and short-run Hamiltonian conservation.
+- Kerr horizon and ISCO limit checks.
+- Circular orbit four-velocity normalization.
+- GRMHD fixture schema, interpolation, and four-velocity normalization.
+- Basic radiative-transfer and polarized-transfer invariants.
+
+The broader repository gate remains:
+
+```bash
+cd python
+python -m compileall blackhole_sim
+python -m pytest -q
+cd ..
+cargo fmt --check --manifest-path native/core/Cargo.toml
+cargo test --manifest-path native/core/Cargo.toml
+python -m blackhole_sim.accelerator_cli doctor --json --fail-on-emulation
+```
+
+## Approximation Boundaries
+
+The following are intentionally not full scientific validation:
+
+- The analytic torus generator is a deterministic fixture. It is not GRMHD
+  solver output.
+- The unpolarized thermal coefficient model is a compact testing fit. Absolute
+  Jy-scale interpretation requires pinned unit conversion and calibration. It is
+  the default `educational_proxy` mode and is rejected when `validated` mode is
+  requested without a non-proxy coefficient model.
+- Coefficient-brick precomputation now uses the same metric-aware local plasma
+  helper as polarized transfer, but image-level scientific validation still
+  requires a selected dump and external baseline.
+- The public-data manifest is currently a collection-level discovery anchor.
+  Release-grade validation requires one concrete selected dump, SHA-256,
+  expected field map, accepted numeric ranges, and an external-code comparison.
+- The `ipole` path provides comparison plumbing, but parity is not established
+  until the same dump, camera, frequency, mass, unit normalization, and image
+  convention are pinned and the comparison report passes.
+
+## WebGPU And Native Policy
+
+WebGPU, CUDA, Metal, HIP, OpenCL, and Rust native paths are acceleration targets,
+not alternate sources of physics truth. They must reproduce the Python reference
+within an accepted regression envelope before any physics parity claim.
+
+Current policy:
+
+- Keep Python as the reference model and validation oracle.
+- Use WebGPU compute for browser-interactive kernels where direct buffer access,
+  shader auditability, and deterministic readback are needed.
+- Use native CPU/GPU kernels only after small hot-loop parity tests pass.
+- Do not move geodesic integration or Stokes transport into a game engine's
+  built-in physics system.
+- A game engine may be used later as a presentation shell only if the physical
+  state is still produced by the validated reference/native/WebGPU kernels.
+
+The current WebGPU Stokes compute shader and static CUDA/Metal/HIP/OpenCL source
+assets derive camera rays through a ZAMO launch helper instead of fixed photon
+momenta. They still must not be described as full physics-equivalent renderers
+until geodesic stepping, coefficient sampling, Stokes transport, and readback
+image regression pass against the CPU reference on target hardware.
+
+## Next Accuracy Milestones
+
+1. Select one public GRMHD dump and record checksum, field map, unit convention,
+   camera convention, and accepted numeric ranges.
+2. Run the same selected dump through this renderer and a local `ipole` baseline.
+3. Add a signed comparison report for Stokes image shape, orientation, flux
+   scale, and tolerance envelope.
+4. Add GPU/native readback tests that compare sampled rays and small rendered
+   Stokes images against the Python reference before claiming parity.
diff --git a/docs/SCIENTIFIC_ROADMAP.md b/docs/SCIENTIFIC_ROADMAP.md
new file mode 100644
index 0000000..cb6f106
--- /dev/null
+++ b/docs/SCIENTIFIC_ROADMAP.md
@@ -0,0 +1,37 @@
+# Scientific Roadmap
+
+## Direction
+
+BlackHole Sim should move from an accelerator-backed visualization prototype to
+a validated GRRT research platform. Python remains the canonical physics
+implementation. Native Rust and WebGPU are acceleration layers that must prove
+parity before they can support correctness, performance, or visual claims.
+
+## Release Gates
+
+- `educational_proxy`: fast visual and fixture paths are allowed, but UI, docs,
+  and CLI output must identify them as approximations.
+- `validated`: requires metric-aware plasma-frame calculations, pinned data
+  provenance, non-proxy coefficient models or coefficient tables, and regression
+  evidence.
+- GPU/native parity: requires Python reference comparison for ray launch,
+  geodesic stepping, coefficient sampling, Stokes stepping, and direct readback
+  or compact comparison reports.
+- Public-data validation: requires one selected dump with SHA-256, field map,
+  accepted ranges, unit/camera/frequency conventions, and external baseline
+  comparison.
+- External baseline: `ipole` comparison reports must record shape, orientation,
+  flux scale, Stokes component metrics, and tolerance envelope.
+
+## Near-Term Work
+
+1. Finish selected Illinois v3 SANE `a=+0.5` dump acquisition metadata without
+   committing the dump.
+2. Add a small ipole baseline report for the selected dump and a matching local
+   Stokes cube report.
+3. Add WebGPU diagnostics automation that proves readback for a tiny Stokes
+   render on the local adapter.
+4. Port the ZAMO launch numeric parity test from static source checks to runtime
+   shader/native checks once vendor/WebGPU test harnesses are available.
+5. Keep game engines out of the physics core. Revisit them only as presentation
+   shells after validated kernels produce the physical state.
diff --git a/native/cuda/kerr_stokes_kernel.cu b/native/cuda/kerr_stokes_kernel.cu
index a4592db..9f2a3b7 100644
--- a/native/cuda/kerr_stokes_kernel.cu
+++ b/native/cuda/kerr_stokes_kernel.cu
@@ -2,12 +2,14 @@
 // One thread = one pixel. Coefficients are flattened as [nr][nt][np][11].
 #include <math.h>
 struct State6 { float t,r,th,ph,pr,pth; };
+struct RayLaunch { State6 y; float p_t; float p_phi; };
 struct Mat4 { float m[16]; };
 __device__ float mget(const Mat4& m,int i,int j){return m.m[i*4+j];}
 __device__ void mset(Mat4& m,int i,int j,float v){m.m[i*4+j]=v;}
 __device__ float clampf(float x,float lo,float hi){return fminf(fmaxf(x,lo),hi);} 
 __device__ float wrap_phi(float x){float y=fmodf(x,6.28318530718f);return y<0?y+6.28318530718f:y;}
 __device__ Mat4 metric_contravariant(float r,float theta,float a){float ct=cosf(theta),st=sinf(theta),s2=fmaxf(st*st,1e-8f);float sig=r*r+a*a*ct*ct,dlt=r*r-2*r+a*a;float A=(r*r+a*a)*(r*r+a*a)-a*a*dlt*s2;Mat4 g={0};mset(g,0,0,-A/(sig*dlt));mset(g,0,3,-2*a*r/(sig*dlt));mset(g,3,0,mget(g,0,3));mset(g,1,1,dlt/sig);mset(g,2,2,1/sig);mset(g,3,3,(dlt-a*a*s2)/(sig*dlt*s2));return g;}
+__device__ RayLaunch zamo_camera_ray_initial_state(int x,int y,int width,int height,float camera_r,float camera_theta,float fov_y,float spin_a){float ndcx=2*((x+.5f)/width)-1,ndcy=1-2*((y+.5f)/height);float aspect=((float)width)/fmaxf((float)height,1.0f),tan_y=tanf(.5f*fov_y);float nr=-1.0f,nth=-ndcy*tan_y,nph=ndcx*aspect*tan_y;float invn=rsqrtf(nr*nr+nth*nth+nph*nph);nr*=invn;nth*=invn;nph*=invn;float r=camera_r,th=clampf(camera_theta,1e-6f,3.1415926f-1e-6f),a=spin_a;float ct=cosf(th),st=sinf(th),s2=fmaxf(st*st,1e-8f),sig=r*r+a*a*ct*ct,dlt=r*r-2*r+a*a,A=(r*r+a*a)*(r*r+a*a)-a*a*dlt*s2;float gtt=-(1-2*r/sig),gtp=-2*a*r*s2/sig,grr=sig/dlt,gth=sig,gpp=A*s2/sig,lapse=sqrtf(sig*dlt/A),omega=2*a*r/A;float ptcon=1/lapse,prcon=nr*sqrtf(dlt/sig),pthcon=nth/sqrtf(sig),pphcon=omega/lapse+nph/sqrtf(gpp);float p_t=gtt*ptcon+gtp*pphcon,p_phi=gtp*ptcon+gpp*pphcon;return {{0,r,th,0,grr*prcon,gth*pthcon},p_t,p_phi};}
 __device__ void metric_derivative_numeric(float r,float th,float a,Mat4& dr,Mat4& dth){float er=1e-3f,et=1e-4f;Mat4 gr0=metric_contravariant(r-er,th,a),gr1=metric_contravariant(r+er,th,a),gt0=metric_contravariant(r,th-et,a),gt1=metric_contravariant(r,th+et,a);for(int i=0;i<4;i++)for(int j=0;j<4;j++){mset(dr,i,j,(mget(gr1,i,j)-mget(gr0,i,j))/(2*er));mset(dth,i,j,(mget(gt1,i,j)-mget(gt0,i,j))/(2*et));}}
 __device__ State6 kerr_rhs(State6 y,float p_t,float p_phi,float a){Mat4 g=metric_contravariant(y.r,y.th,a),dgdr,dgdth;metric_derivative_numeric(y.r,y.th,a,dgdr,dgdth);float p[4]={p_t,y.pr,y.pth,p_phi},xd[4]={0,0,0,0};for(int mu=0;mu<4;mu++)for(int nu=0;nu<4;nu++)xd[mu]+=mget(g,mu,nu)*p[nu];float pr=0,pth=0;for(int mu=0;mu<4;mu++)for(int nu=0;nu<4;nu++){pr+=-0.5f*p[mu]*mget(dgdr,mu,nu)*p[nu];pth+=-0.5f*p[mu]*mget(dgdth,mu,nu)*p[nu];}return {xd[0],xd[1],xd[2],xd[3],pr,pth};}
 __device__ State6 rk2_geodesic_step(State6 y,float h,float p_t,float p_phi,float a){State6 k1=kerr_rhs(y,p_t,p_phi,a);State6 mid={y.t+.5f*h*k1.t,y.r+.5f*h*k1.r,y.th+.5f*h*k1.th,y.ph+.5f*h*k1.ph,y.pr+.5f*h*k1.pr,y.pth+.5f*h*k1.pth};State6 k2=kerr_rhs(mid,p_t,p_phi,a);return {y.t+h*k2.t,y.r+h*k2.r,clampf(y.th+h*k2.th,1e-6f,3.1415926f-1e-6f),y.ph+h*k2.ph,y.pr+h*k2.pr,y.pth+h*k2.pth};}
@@ -17,4 +19,4 @@ __device__ void bracket_phi(const float* g,int n,float ph,int& i0,int& i1,float&
 __device__ int sample_brick_trilinear(const float* coeffs,const float* rg,const float* tg,const float* pg,int nr,int nt,int np,float r,float th,float ph,float outv[11]){int r0,r1,t0,t1,p0,p1;float wr,wt,wp;if(!bracket_linear(rg,nr,r,r0,r1,wr)||!bracket_linear(tg,nt,th,t0,t1,wt)){for(int c=0;c<11;c++)outv[c]=NAN;return 0;}bracket_phi(pg,np,ph,p0,p1,wp);for(int c=0;c<11;c++){float c000=coeffs[cidx(r0,t0,p0,c,nt,np)],c001=coeffs[cidx(r0,t0,p1,c,nt,np)],c010=coeffs[cidx(r0,t1,p0,c,nt,np)],c011=coeffs[cidx(r0,t1,p1,c,nt,np)],c100=coeffs[cidx(r1,t0,p0,c,nt,np)],c101=coeffs[cidx(r1,t0,p1,c,nt,np)],c110=coeffs[cidx(r1,t1,p0,c,nt,np)],c111=coeffs[cidx(r1,t1,p1,c,nt,np)];float c00=c000+wp*(c001-c000),c01=c010+wp*(c011-c010),c10=c100+wp*(c101-c100),c11=c110+wp*(c111-c110);outv[c]=(c00+wt*(c01-c00))+wr*((c10+wt*(c11-c10))-(c00+wt*(c01-c00)));}return 1;}
 __device__ void rhs_stokes(const float S[4],const float c[11],float out[4]){out[0]=c[0]-(c[4]*S[0]+c[5]*S[1]+c[6]*S[2]+c[7]*S[3]);out[1]=c[1]-(c[5]*S[0]+c[4]*S[1]+c[8]*S[2]-c[9]*S[3]);out[2]=c[2]-(c[6]*S[0]-c[8]*S[1]+c[4]*S[2]+c[10]*S[3]);out[3]=c[3]-(c[7]*S[0]+c[9]*S[1]-c[10]*S[2]+c[4]*S[3]);}
 __device__ void stokes_step_rk2(float S[4],const float c[11],float ds){float k1[4],mid[4],k2[4];rhs_stokes(S,c,k1);for(int i=0;i<4;i++)mid[i]=S[i]+0.5f*ds*k1[i];rhs_stokes(mid,c,k2);for(int i=0;i<4;i++)S[i]+=ds*k2[i];}
-extern "C" __global__ void kerr_stokes_render_kernel(float* out_stokes,const float* coeffs,const float* r_grid,const float* theta_grid,const float* phi_grid,int width,int height,int nr,int nt,int np,float spin_a,int max_steps,float step){int x=blockIdx.x*blockDim.x+threadIdx.x,y=blockIdx.y*blockDim.y+threadIdx.y;if(x>=width||y>=height)return;int pix=y*width+x;float ndcx=2*((x+.5f)/width)-1,ndcy=1-2*((y+.5f)/height);State6 s={0,55,1.134464f,0,-1,-0.22f*ndcy};float p_t=-1,p_phi=0.12f*ndcx,S[4]={0,0,0,0};float rplus=1+sqrtf(fmaxf(1-spin_a*spin_a,0));for(int n=0;n<max_steps;n++){State6 prev=s;s=rk2_geodesic_step(s,step,p_t,p_phi,spin_a);float c[11];int valid=sample_brick_trilinear(coeffs,r_grid,theta_grid,phi_grid,nr,nt,np,.5f*(prev.r+s.r),.5f*(prev.th+s.th),.5f*(prev.ph+s.ph),c);if(valid){float ds=sqrtf((s.r-prev.r)*(s.r-prev.r)+(s.th-prev.th)*(s.th-prev.th)+(s.ph-prev.ph)*(s.ph-prev.ph));stokes_step_rk2(S,c,ds);}if(s.r<=rplus*1.0002f||s.r>220||!isfinite(s.r))break;}out_stokes[4*pix]=S[0];out_stokes[4*pix+1]=S[1];out_stokes[4*pix+2]=S[2];out_stokes[4*pix+3]=S[3];}
+extern "C" __global__ void kerr_stokes_render_kernel(float* out_stokes,const float* coeffs,const float* r_grid,const float* theta_grid,const float* phi_grid,int width,int height,int nr,int nt,int np,float spin_a,int max_steps,float step){int x=blockIdx.x*blockDim.x+threadIdx.x,y=blockIdx.y*blockDim.y+threadIdx.y;if(x>=width||y>=height)return;int pix=y*width+x;RayLaunch launch=zamo_camera_ray_initial_state(x,y,width,height,55.0f,1.134464f,0.5934119f,spin_a);State6 s=launch.y;float p_t=launch.p_t,p_phi=launch.p_phi,S[4]={0,0,0,0};float rplus=1+sqrtf(fmaxf(1-spin_a*spin_a,0));for(int n=0;n<max_steps;n++){State6 prev=s;s=rk2_geodesic_step(s,step,p_t,p_phi,spin_a);float c[11];int valid=sample_brick_trilinear(coeffs,r_grid,theta_grid,phi_grid,nr,nt,np,.5f*(prev.r+s.r),.5f*(prev.th+s.th),.5f*(prev.ph+s.ph),c);if(valid){float ds=sqrtf((s.r-prev.r)*(s.r-prev.r)+(s.th-prev.th)*(s.th-prev.th)+(s.ph-prev.ph)*(s.ph-prev.ph));stokes_step_rk2(S,c,ds);}if(s.r<=rplus*1.0002f||s.r>220||!isfinite(s.r))break;}out_stokes[4*pix]=S[0];out_stokes[4*pix+1]=S[1];out_stokes[4*pix+2]=S[2];out_stokes[4*pix+3]=S[3];}
diff --git a/native/metal/kerr_stokes_kernel.metal b/native/metal/kerr_stokes_kernel.metal
index f6da3bf..5ad1d67 100644
--- a/native/metal/kerr_stokes_kernel.metal
+++ b/native/metal/kerr_stokes_kernel.metal
@@ -3,6 +3,7 @@ using namespace metal;
 
 struct RenderParams { uint width; uint height; uint nr; uint ntheta; uint nphi; float spin_a; uint max_steps; float step; };
 struct State6 { float t; float r; float th; float ph; float pr; float pth; };
+struct RayLaunch { State6 y; float p_t; float p_phi; };
 
 inline float wrap_phi(float x) { float y = fmod(x, 6.28318530718f); return y < 0.0f ? y + 6.28318530718f : y; }
 inline uint cidx(uint ir, uint it, uint ip, uint c, uint nt, uint np) { return (((ir * nt) + it) * np + ip) * 11u + c; }
@@ -17,6 +18,25 @@ float4x4 metric_contravariant(float r, float theta, float a) {
     return g;
 }
 
+RayLaunch zamo_camera_ray_initial_state(uint x, uint y, uint width, uint height, float camera_r, float camera_theta, float fov_y, float spin_a) {
+    float ndcx = 2.0f * ((float(x) + .5f) / float(width)) - 1.0f;
+    float ndcy = 1.0f - 2.0f * ((float(y) + .5f) / float(height));
+    float aspect = float(width) / max(float(height), 1.0f), tan_y = tan(.5f * fov_y);
+    float nr = -1.0f, nth = -ndcy * tan_y, nph = ndcx * aspect * tan_y;
+    float invn = rsqrt(nr * nr + nth * nth + nph * nph);
+    nr *= invn; nth *= invn; nph *= invn;
+    float r = camera_r, th = clamp(camera_theta, 1.0e-6f, 3.1415926f - 1.0e-6f), a = spin_a;
+    float ct = cos(th), st = sin(th), s2 = max(st * st, 1.0e-8f);
+    float sig = r * r + a * a * ct * ct, dlt = r * r - 2.0f * r + a * a;
+    float A = (r * r + a * a) * (r * r + a * a) - a * a * dlt * s2;
+    float gtt = -(1.0f - 2.0f * r / sig), gtp = -2.0f * a * r * s2 / sig;
+    float grr = sig / dlt, gth = sig, gpp = A * s2 / sig;
+    float lapse = sqrt(sig * dlt / A), omega = 2.0f * a * r / A;
+    float ptcon = 1.0f / lapse, prcon = nr * sqrt(dlt / sig), pthcon = nth / sqrt(sig), pphcon = omega / lapse + nph / sqrt(gpp);
+    float p_t = gtt * ptcon + gtp * pphcon, p_phi = gtp * ptcon + gpp * pphcon;
+    return RayLaunch{State6{0.0f, r, th, 0.0f, grr * prcon, gth * pthcon}, p_t, p_phi};
+}
+
 void metric_derivative_numeric(float r, float th, float a, thread float4x4& dr, thread float4x4& dth) {
     float er = 1.0e-3f, et = 1.0e-4f;
     float4x4 gr0 = metric_contravariant(r-er, th, a), gr1 = metric_contravariant(r+er, th, a);
@@ -64,8 +84,8 @@ void stokes_step_rk2(thread float S[4], thread const float c[11], float ds) { fl
 
 kernel void kerr_stokes_render_kernel(device float4* out_stokes [[buffer(0)]], device const float* coeffs [[buffer(1)]], device const float* r_grid [[buffer(2)]], device const float* theta_grid [[buffer(3)]], device const float* phi_grid [[buffer(4)]], constant RenderParams& params [[buffer(5)]], uint2 gid [[thread_position_in_grid]]) {
     if (gid.x >= params.width || gid.y >= params.height) return; uint pix = gid.y * params.width + gid.x;
-    float ndcx = 2.0f * ((float(gid.x) + .5f) / float(params.width)) - 1.0f; float ndcy = 1.0f - 2.0f * ((float(gid.y) + .5f) / float(params.height));
-    State6 s = State6{0.0f,55.0f,1.134464f,0.0f,-1.0f,-.22f*ndcy}; float p_t=-1.0f,p_phi=.12f*ndcx,S[4]={0,0,0,0}; float rplus=1.0f+sqrt(max(1.0f-params.spin_a*params.spin_a,0.0f));
+    RayLaunch launch = zamo_camera_ray_initial_state(gid.x, gid.y, params.width, params.height, 55.0f, 1.134464f, 0.5934119f, params.spin_a);
+    State6 s = launch.y; float p_t=launch.p_t,p_phi=launch.p_phi,S[4]={0,0,0,0}; float rplus=1.0f+sqrt(max(1.0f-params.spin_a*params.spin_a,0.0f));
     for(uint n=0;n<params.max_steps;++n){ State6 prev=s; s=rk2_geodesic_step(s,params.step,p_t,p_phi,params.spin_a); float c[11]; bool valid=sample_brick_trilinear(coeffs,r_grid,theta_grid,phi_grid,params.nr,params.ntheta,params.nphi,.5f*(prev.r+s.r),.5f*(prev.th+s.th),.5f*(prev.ph+s.ph),c); if(valid){ float ds=length(float3(s.r-prev.r,s.th-prev.th,s.ph-prev.ph)); stokes_step_rk2(S,c,ds); } if(s.r<=rplus*1.0002f||s.r>220.0f)break; }
     out_stokes[pix] = float4(S[0],S[1],S[2],S[3]);
 }
diff --git a/native/opencl/kerr_stokes_kernel.cl b/native/opencl/kerr_stokes_kernel.cl
index 1bacda3..40ad612 100644
--- a/native/opencl/kerr_stokes_kernel.cl
+++ b/native/opencl/kerr_stokes_kernel.cl
@@ -1,8 +1,10 @@
 // v0.7 OpenCL/SYCL-style hot loop for Intel/AMD/ARM targets.
 typedef struct { float t,r,th,ph,pr,pth; } State6;
+typedef struct { State6 y; float p_t; float p_phi; } RayLaunch;
 inline float wrap_phi(float x){float y=fmod(x,6.28318530718f);return y<0?y+6.28318530718f:y;}
 inline int cidx(int ir,int it,int ip,int c,int nt,int np){return (((ir*nt)+it)*np+ip)*11+c;}
 float16 metric_contravariant(float r,float theta,float a){float ct=cos(theta),st=sin(theta),s2=fmax(st*st,1e-8f);float sig=r*r+a*a*ct*ct,dlt=r*r-2*r+a*a;float A=(r*r+a*a)*(r*r+a*a)-a*a*dlt*s2;float16 g=(float16)(0);g.s0=-A/(sig*dlt);g.s3=-2*a*r/(sig*dlt);g.sc=g.s3;g.s5=dlt/sig;g.sa=1/sig;g.sf=(dlt-a*a*s2)/(sig*dlt*s2);return g;}
+RayLaunch zamo_camera_ray_initial_state(int x,int y,int width,int height,float camera_r,float camera_theta,float fov_y,float spin_a){float ndcx=2*((x+.5f)/width)-1,ndcy=1-2*((y+.5f)/height);float aspect=((float)width)/fmax((float)height,1.0f),tan_y=tan(.5f*fov_y);float nr=-1.0f,nth=-ndcy*tan_y,nph=ndcx*aspect*tan_y;float invn=rsqrt(nr*nr+nth*nth+nph*nph);nr*=invn;nth*=invn;nph*=invn;float r=camera_r,th=clamp(camera_theta,1e-6f,3.1415926f-1e-6f),a=spin_a;float ct=cos(th),st=sin(th),s2=fmax(st*st,1e-8f),sig=r*r+a*a*ct*ct,dlt=r*r-2*r+a*a,A=(r*r+a*a)*(r*r+a*a)-a*a*dlt*s2;float gtt=-(1-2*r/sig),gtp=-2*a*r*s2/sig,grr=sig/dlt,gth=sig,gpp=A*s2/sig,lapse=sqrt(sig*dlt/A),omega=2*a*r/A;float ptcon=1/lapse,prcon=nr*sqrt(dlt/sig),pthcon=nth/sqrt(sig),pphcon=omega/lapse+nph/sqrt(gpp);float p_t=gtt*ptcon+gtp*pphcon,p_phi=gtp*ptcon+gpp*pphcon;RayLaunch out;out.y=(State6)(0,r,th,0,grr*prcon,gth*pthcon);out.p_t=p_t;out.p_phi=p_phi;return out;}
 inline float mget(float16 m,int i,int j){return m.s[i*4+j];}
 void metric_derivative_numeric(float r,float th,float a,__private float16* dr,__private float16* dt){float er=1e-3f,et=1e-4f;float16 r0=metric_contravariant(r-er,th,a),r1=metric_contravariant(r+er,th,a),t0=metric_contravariant(r,th-et,a),t1=metric_contravariant(r,th+et,a);*dr=(r1-r0)/(2*er);*dt=(t1-t0)/(2*et);}
 State6 kerr_rhs(State6 y,float p_t,float p_phi,float a){float16 g=metric_contravariant(y.r,y.th,a),dr,dt;metric_derivative_numeric(y.r,y.th,a,&dr,&dt);float p[4]={p_t,y.pr,y.pth,p_phi},xd[4]={0,0,0,0};for(int mu=0;mu<4;mu++)for(int nu=0;nu<4;nu++)xd[mu]+=mget(g,mu,nu)*p[nu];float pr=0,pth=0;for(int mu=0;mu<4;mu++)for(int nu=0;nu<4;nu++){pr+=-.5f*p[mu]*mget(dr,mu,nu)*p[nu];pth+=-.5f*p[mu]*mget(dt,mu,nu)*p[nu];}return (State6)(xd[0],xd[1],xd[2],xd[3],pr,pth);}
@@ -12,4 +14,4 @@ void bracket_phi(__global const float* g,int n,float ph,__private int* i0,__priv
 int sample_brick_trilinear(__global const float* coeffs,__global const float* rg,__global const float* tg,__global const float* pg,int nr,int nt,int np,float r,float th,float ph,__private float outv[11]){int r0,r1,t0,t1,p0,p1;float wr,wt,wp;if(!bracket_linear(rg,nr,r,&r0,&r1,&wr)||!bracket_linear(tg,nt,th,&t0,&t1,&wt)){for(int c=0;c<11;c++)outv[c]=nan((uint)0);return 0;}bracket_phi(pg,np,ph,&p0,&p1,&wp);for(int c=0;c<11;c++){float c000=coeffs[cidx(r0,t0,p0,c,nt,np)],c001=coeffs[cidx(r0,t0,p1,c,nt,np)],c010=coeffs[cidx(r0,t1,p0,c,nt,np)],c011=coeffs[cidx(r0,t1,p1,c,nt,np)],c100=coeffs[cidx(r1,t0,p0,c,nt,np)],c101=coeffs[cidx(r1,t0,p1,c,nt,np)],c110=coeffs[cidx(r1,t1,p0,c,nt,np)],c111=coeffs[cidx(r1,t1,p1,c,nt,np)];float c00=mix(c000,c001,wp),c01=mix(c010,c011,wp),c10=mix(c100,c101,wp),c11=mix(c110,c111,wp);outv[c]=mix(mix(c00,c01,wt),mix(c10,c11,wt),wr);}return 1;}
 void rhs_stokes(__private const float S[4],__private const float c[11],__private float o[4]){o[0]=c[0]-(c[4]*S[0]+c[5]*S[1]+c[6]*S[2]+c[7]*S[3]);o[1]=c[1]-(c[5]*S[0]+c[4]*S[1]+c[8]*S[2]-c[9]*S[3]);o[2]=c[2]-(c[6]*S[0]-c[8]*S[1]+c[4]*S[2]+c[10]*S[3]);o[3]=c[3]-(c[7]*S[0]+c[9]*S[1]-c[10]*S[2]+c[4]*S[3]);}
 void stokes_step_rk2(__private float S[4],__private const float c[11],float ds){float k1[4],m[4],k2[4];rhs_stokes(S,c,k1);for(int i=0;i<4;i++)m[i]=S[i]+.5f*ds*k1[i];rhs_stokes(m,c,k2);for(int i=0;i<4;i++)S[i]+=ds*k2[i];}
-__kernel void kerr_stokes_render_kernel(__global float4* out_stokes,__global const float* coeffs,__global const float* r_grid,__global const float* theta_grid,__global const float* phi_grid,const int width,const int height,const int nr,const int ntheta,const int nphi,const float spin_a,const int max_steps,const float step){int x=get_global_id(0),y=get_global_id(1);if(x>=width||y>=height)return;int pix=y*width+x;float ndcx=2*((x+.5f)/width)-1,ndcy=1-2*((y+.5f)/height);State6 s=(State6)(0,55,1.134464f,0,-1,-.22f*ndcy);float p_t=-1,p_phi=.12f*ndcx,S[4]={0,0,0,0};float rplus=1+sqrt(fmax(1-spin_a*spin_a,0.0f));for(int n=0;n<max_steps;n++){State6 prev=s;s=rk2_geodesic_step(s,step,p_t,p_phi,spin_a);float c[11];int valid=sample_brick_trilinear(coeffs,r_grid,theta_grid,phi_grid,nr,ntheta,nphi,.5f*(prev.r+s.r),.5f*(prev.th+s.th),.5f*(prev.ph+s.ph),c);if(valid){float ds=sqrt((s.r-prev.r)*(s.r-prev.r)+(s.th-prev.th)*(s.th-prev.th)+(s.ph-prev.ph)*(s.ph-prev.ph));stokes_step_rk2(S,c,ds);}if(s.r<=rplus*1.0002f||s.r>220)break;}out_stokes[pix]=(float4)(S[0],S[1],S[2],S[3]);}
+__kernel void kerr_stokes_render_kernel(__global float4* out_stokes,__global const float* coeffs,__global const float* r_grid,__global const float* theta_grid,__global const float* phi_grid,const int width,const int height,const int nr,const int ntheta,const int nphi,const float spin_a,const int max_steps,const float step){int x=get_global_id(0),y=get_global_id(1);if(x>=width||y>=height)return;int pix=y*width+x;RayLaunch launch=zamo_camera_ray_initial_state(x,y,width,height,55.0f,1.134464f,0.5934119f,spin_a);State6 s=launch.y;float p_t=launch.p_t,p_phi=launch.p_phi,S[4]={0,0,0,0};float rplus=1+sqrt(fmax(1-spin_a*spin_a,0.0f));for(int n=0;n<max_steps;n++){State6 prev=s;s=rk2_geodesic_step(s,step,p_t,p_phi,spin_a);float c[11];int valid=sample_brick_trilinear(coeffs,r_grid,theta_grid,phi_grid,nr,ntheta,nphi,.5f*(prev.r+s.r),.5f*(prev.th+s.th),.5f*(prev.ph+s.ph),c);if(valid){float ds=sqrt((s.r-prev.r)*(s.r-prev.r)+(s.th-prev.th)*(s.th-prev.th)+(s.ph-prev.ph)*(s.ph-prev.ph));stokes_step_rk2(S,c,ds);}if(s.r<=rplus*1.0002f||s.r>220)break;}out_stokes[pix]=(float4)(S[0],S[1],S[2],S[3]);}
diff --git a/native/rocm/kerr_stokes_kernel.hip b/native/rocm/kerr_stokes_kernel.hip
index a4592db..9f2a3b7 100644
--- a/native/rocm/kerr_stokes_kernel.hip
+++ b/native/rocm/kerr_stokes_kernel.hip
@@ -2,12 +2,14 @@
 // One thread = one pixel. Coefficients are flattened as [nr][nt][np][11].
 #include <math.h>
 struct State6 { float t,r,th,ph,pr,pth; };
+struct RayLaunch { State6 y; float p_t; float p_phi; };
 struct Mat4 { float m[16]; };
 __device__ float mget(const Mat4& m,int i,int j){return m.m[i*4+j];}
 __device__ void mset(Mat4& m,int i,int j,float v){m.m[i*4+j]=v;}
 __device__ float clampf(float x,float lo,float hi){return fminf(fmaxf(x,lo),hi);} 
 __device__ float wrap_phi(float x){float y=fmodf(x,6.28318530718f);return y<0?y+6.28318530718f:y;}
 __device__ Mat4 metric_contravariant(float r,float theta,float a){float ct=cosf(theta),st=sinf(theta),s2=fmaxf(st*st,1e-8f);float sig=r*r+a*a*ct*ct,dlt=r*r-2*r+a*a;float A=(r*r+a*a)*(r*r+a*a)-a*a*dlt*s2;Mat4 g={0};mset(g,0,0,-A/(sig*dlt));mset(g,0,3,-2*a*r/(sig*dlt));mset(g,3,0,mget(g,0,3));mset(g,1,1,dlt/sig);mset(g,2,2,1/sig);mset(g,3,3,(dlt-a*a*s2)/(sig*dlt*s2));return g;}
+__device__ RayLaunch zamo_camera_ray_initial_state(int x,int y,int width,int height,float camera_r,float camera_theta,float fov_y,float spin_a){float ndcx=2*((x+.5f)/width)-1,ndcy=1-2*((y+.5f)/height);float aspect=((float)width)/fmaxf((float)height,1.0f),tan_y=tanf(.5f*fov_y);float nr=-1.0f,nth=-ndcy*tan_y,nph=ndcx*aspect*tan_y;float invn=rsqrtf(nr*nr+nth*nth+nph*nph);nr*=invn;nth*=invn;nph*=invn;float r=camera_r,th=clampf(camera_theta,1e-6f,3.1415926f-1e-6f),a=spin_a;float ct=cosf(th),st=sinf(th),s2=fmaxf(st*st,1e-8f),sig=r*r+a*a*ct*ct,dlt=r*r-2*r+a*a,A=(r*r+a*a)*(r*r+a*a)-a*a*dlt*s2;float gtt=-(1-2*r/sig),gtp=-2*a*r*s2/sig,grr=sig/dlt,gth=sig,gpp=A*s2/sig,lapse=sqrtf(sig*dlt/A),omega=2*a*r/A;float ptcon=1/lapse,prcon=nr*sqrtf(dlt/sig),pthcon=nth/sqrtf(sig),pphcon=omega/lapse+nph/sqrtf(gpp);float p_t=gtt*ptcon+gtp*pphcon,p_phi=gtp*ptcon+gpp*pphcon;return {{0,r,th,0,grr*prcon,gth*pthcon},p_t,p_phi};}
 __device__ void metric_derivative_numeric(float r,float th,float a,Mat4& dr,Mat4& dth){float er=1e-3f,et=1e-4f;Mat4 gr0=metric_contravariant(r-er,th,a),gr1=metric_contravariant(r+er,th,a),gt0=metric_contravariant(r,th-et,a),gt1=metric_contravariant(r,th+et,a);for(int i=0;i<4;i++)for(int j=0;j<4;j++){mset(dr,i,j,(mget(gr1,i,j)-mget(gr0,i,j))/(2*er));mset(dth,i,j,(mget(gt1,i,j)-mget(gt0,i,j))/(2*et));}}
 __device__ State6 kerr_rhs(State6 y,float p_t,float p_phi,float a){Mat4 g=metric_contravariant(y.r,y.th,a),dgdr,dgdth;metric_derivative_numeric(y.r,y.th,a,dgdr,dgdth);float p[4]={p_t,y.pr,y.pth,p_phi},xd[4]={0,0,0,0};for(int mu=0;mu<4;mu++)for(int nu=0;nu<4;nu++)xd[mu]+=mget(g,mu,nu)*p[nu];float pr=0,pth=0;for(int mu=0;mu<4;mu++)for(int nu=0;nu<4;nu++){pr+=-0.5f*p[mu]*mget(dgdr,mu,nu)*p[nu];pth+=-0.5f*p[mu]*mget(dgdth,mu,nu)*p[nu];}return {xd[0],xd[1],xd[2],xd[3],pr,pth};}
 __device__ State6 rk2_geodesic_step(State6 y,float h,float p_t,float p_phi,float a){State6 k1=kerr_rhs(y,p_t,p_phi,a);State6 mid={y.t+.5f*h*k1.t,y.r+.5f*h*k1.r,y.th+.5f*h*k1.th,y.ph+.5f*h*k1.ph,y.pr+.5f*h*k1.pr,y.pth+.5f*h*k1.pth};State6 k2=kerr_rhs(mid,p_t,p_phi,a);return {y.t+h*k2.t,y.r+h*k2.r,clampf(y.th+h*k2.th,1e-6f,3.1415926f-1e-6f),y.ph+h*k2.ph,y.pr+h*k2.pr,y.pth+h*k2.pth};}
@@ -17,4 +19,4 @@ __device__ void bracket_phi(const float* g,int n,float ph,int& i0,int& i1,float&
 __device__ int sample_brick_trilinear(const float* coeffs,const float* rg,const float* tg,const float* pg,int nr,int nt,int np,float r,float th,float ph,float outv[11]){int r0,r1,t0,t1,p0,p1;float wr,wt,wp;if(!bracket_linear(rg,nr,r,r0,r1,wr)||!bracket_linear(tg,nt,th,t0,t1,wt)){for(int c=0;c<11;c++)outv[c]=NAN;return 0;}bracket_phi(pg,np,ph,p0,p1,wp);for(int c=0;c<11;c++){float c000=coeffs[cidx(r0,t0,p0,c,nt,np)],c001=coeffs[cidx(r0,t0,p1,c,nt,np)],c010=coeffs[cidx(r0,t1,p0,c,nt,np)],c011=coeffs[cidx(r0,t1,p1,c,nt,np)],c100=coeffs[cidx(r1,t0,p0,c,nt,np)],c101=coeffs[cidx(r1,t0,p1,c,nt,np)],c110=coeffs[cidx(r1,t1,p0,c,nt,np)],c111=coeffs[cidx(r1,t1,p1,c,nt,np)];float c00=c000+wp*(c001-c000),c01=c010+wp*(c011-c010),c10=c100+wp*(c101-c100),c11=c110+wp*(c111-c110);outv[c]=(c00+wt*(c01-c00))+wr*((c10+wt*(c11-c10))-(c00+wt*(c01-c00)));}return 1;}
 __device__ void rhs_stokes(const float S[4],const float c[11],float out[4]){out[0]=c[0]-(c[4]*S[0]+c[5]*S[1]+c[6]*S[2]+c[7]*S[3]);out[1]=c[1]-(c[5]*S[0]+c[4]*S[1]+c[8]*S[2]-c[9]*S[3]);out[2]=c[2]-(c[6]*S[0]-c[8]*S[1]+c[4]*S[2]+c[10]*S[3]);out[3]=c[3]-(c[7]*S[0]+c[9]*S[1]-c[10]*S[2]+c[4]*S[3]);}
 __device__ void stokes_step_rk2(float S[4],const float c[11],float ds){float k1[4],mid[4],k2[4];rhs_stokes(S,c,k1);for(int i=0;i<4;i++)mid[i]=S[i]+0.5f*ds*k1[i];rhs_stokes(mid,c,k2);for(int i=0;i<4;i++)S[i]+=ds*k2[i];}
-extern "C" __global__ void kerr_stokes_render_kernel(float* out_stokes,const float* coeffs,const float* r_grid,const float* theta_grid,const float* phi_grid,int width,int height,int nr,int nt,int np,float spin_a,int max_steps,float step){int x=blockIdx.x*blockDim.x+threadIdx.x,y=blockIdx.y*blockDim.y+threadIdx.y;if(x>=width||y>=height)return;int pix=y*width+x;float ndcx=2*((x+.5f)/width)-1,ndcy=1-2*((y+.5f)/height);State6 s={0,55,1.134464f,0,-1,-0.22f*ndcy};float p_t=-1,p_phi=0.12f*ndcx,S[4]={0,0,0,0};float rplus=1+sqrtf(fmaxf(1-spin_a*spin_a,0));for(int n=0;n<max_steps;n++){State6 prev=s;s=rk2_geodesic_step(s,step,p_t,p_phi,spin_a);float c[11];int valid=sample_brick_trilinear(coeffs,r_grid,theta_grid,phi_grid,nr,nt,np,.5f*(prev.r+s.r),.5f*(prev.th+s.th),.5f*(prev.ph+s.ph),c);if(valid){float ds=sqrtf((s.r-prev.r)*(s.r-prev.r)+(s.th-prev.th)*(s.th-prev.th)+(s.ph-prev.ph)*(s.ph-prev.ph));stokes_step_rk2(S,c,ds);}if(s.r<=rplus*1.0002f||s.r>220||!isfinite(s.r))break;}out_stokes[4*pix]=S[0];out_stokes[4*pix+1]=S[1];out_stokes[4*pix+2]=S[2];out_stokes[4*pix+3]=S[3];}
+extern "C" __global__ void kerr_stokes_render_kernel(float* out_stokes,const float* coeffs,const float* r_grid,const float* theta_grid,const float* phi_grid,int width,int height,int nr,int nt,int np,float spin_a,int max_steps,float step){int x=blockIdx.x*blockDim.x+threadIdx.x,y=blockIdx.y*blockDim.y+threadIdx.y;if(x>=width||y>=height)return;int pix=y*width+x;RayLaunch launch=zamo_camera_ray_initial_state(x,y,width,height,55.0f,1.134464f,0.5934119f,spin_a);State6 s=launch.y;float p_t=launch.p_t,p_phi=launch.p_phi,S[4]={0,0,0,0};float rplus=1+sqrtf(fmaxf(1-spin_a*spin_a,0));for(int n=0;n<max_steps;n++){State6 prev=s;s=rk2_geodesic_step(s,step,p_t,p_phi,spin_a);float c[11];int valid=sample_brick_trilinear(coeffs,r_grid,theta_grid,phi_grid,nr,nt,np,.5f*(prev.r+s.r),.5f*(prev.th+s.th),.5f*(prev.ph+s.ph),c);if(valid){float ds=sqrtf((s.r-prev.r)*(s.r-prev.r)+(s.th-prev.th)*(s.th-prev.th)+(s.ph-prev.ph)*(s.ph-prev.ph));stokes_step_rk2(S,c,ds);}if(s.r<=rplus*1.0002f||s.r>220||!isfinite(s.r))break;}out_stokes[4*pix]=S[0];out_stokes[4*pix+1]=S[1];out_stokes[4*pix+2]=S[2];out_stokes[4*pix+3]=S[3];}
diff --git a/python/blackhole_sim/__init__.py b/python/blackhole_sim/__init__.py
index 22fdbc3..49c23a4 100644
--- a/python/blackhole_sim/__init__.py
+++ b/python/blackhole_sim/__init__.py
@@ -37,7 +37,14 @@
 from .calibration import PhysicalScaling, calibrate_flux_scale
 from .grmhd_adapters import load_harm_hdf5, load_bhac_hdf5, load_koral_hdf5
 from .polarized_transfer import integrate_polarized_kerr_grrt, stokes_step_exact
-from .synchrotron import HybridSynchrotronCoefficients, ThermalSynchrotronCoefficients, NonthermalPowerLawSynchrotronCoefficients
+from .synchrotron import (
+    HybridSynchrotronCoefficients,
+    NonthermalPowerLawSynchrotronCoefficients,
+    ThermalSynchrotronCoefficients,
+    local_plasma_from_sample,
+    magnetic_field_strength_code,
+    magnetic_pitch_cosine,
+)
 
 from .accelerated_renderer import (
     AcceleratedRenderConfig,
@@ -98,6 +105,9 @@
     "HybridSynchrotronCoefficients",
     "ThermalSynchrotronCoefficients",
     "NonthermalPowerLawSynchrotronCoefficients",
+    "local_plasma_from_sample",
+    "magnetic_field_strength_code",
+    "magnetic_pitch_cosine",
     "AcceleratedRenderConfig",
     "render_stokes_image_bricks",
     "render_progressive_stokes_bricks",
diff --git a/python/blackhole_sim/coefficient_bricks.py b/python/blackhole_sim/coefficient_bricks.py
index 13d038a..93f1535 100644
--- a/python/blackhole_sim/coefficient_bricks.py
+++ b/python/blackhole_sim/coefficient_bricks.py
@@ -16,7 +16,7 @@
 
 from .calibration import PhysicalScaling
 from .grmhd import GRMHDSnapshot, FluidSample
-from .synchrotron import HybridSynchrotronCoefficients, LocalPlasmaFrame
+from .synchrotron import HybridSynchrotronCoefficients, LocalPlasmaFrame, local_plasma_from_sample
 
 COEFF_NAMES: tuple[str, ...] = (
     "j_i", "j_q", "j_u", "j_v",
@@ -107,13 +107,8 @@ def _cell_sample(snapshot: GRMHDSnapshot, i: int, j: int, k: int) -> FluidSample
     )
 
 
-def _fast_frame(sample: FluidSample, scaling: PhysicalScaling) -> LocalPlasmaFrame:
-    b_vec = np.asarray(sample.b_con[1:], dtype=float)
-    b_code = float(np.linalg.norm(b_vec))
-    b_gauss = float(scaling.magnetic_field_gauss(b_code))
-    n_e = float(scaling.electron_number_density_cm3(sample.rho))
-    evpa = float(np.arctan2(b_vec[2], b_vec[1])) if b_code > 0.0 else 0.0
-    return LocalPlasmaFrame(n_e, max(sample.theta_e, 1e-8), max(b_gauss, 0.0), 0.0, evpa)
+def _fast_frame(sample: FluidSample, scaling: PhysicalScaling, spin_a: float) -> LocalPlasmaFrame:
+    return local_plasma_from_sample(sample, scaling, spin_a=spin_a)
 
 
 def precompute_coefficient_bricks(
@@ -141,7 +136,7 @@ def precompute_coefficient_bricks(
         for jj, j in enumerate(t_idx):
             for kk, k in enumerate(p_idx):
                 sample = _cell_sample(snapshot, int(i), int(j), int(k))
-                frame = _fast_frame(sample, scaling)
+                frame = _fast_frame(sample, scaling, float(snapshot.spin_a))
                 c = coeff_model.coefficients(frame, nu_hz)
                 coeffs[ii, jj, kk] = np.asarray([
                     c.j_i, c.j_q, c.j_u, c.j_v,
diff --git a/python/blackhole_sim/polarized_transfer.py b/python/blackhole_sim/polarized_transfer.py
index 07d2fc1..96bc9b5 100644
--- a/python/blackhole_sim/polarized_transfer.py
+++ b/python/blackhole_sim/polarized_transfer.py
@@ -117,7 +117,7 @@ def integrate_polarized_kerr_grrt(
         if g_shift < pcfg.min_redshift:
             continue
         nu_emit = pcfg.observing_frequency_hz / g_shift
-        frame = local_plasma_from_sample(sample, scaling, p_cov)
+        frame = local_plasma_from_sample(sample, scaling, p_cov, spin_a=float(snapshot.spin_a))
         coeff = model.coefficients(frame, nu_emit)
         # Invariant transfer: j/nu^2 and alpha*nu are invariant; for a compact
         # screen-frame Stokes integral we use g^2 on emission and g^-1 on K.
diff --git a/python/blackhole_sim/radiative_transfer.py b/python/blackhole_sim/radiative_transfer.py
index 688003e..9ea1c9f 100644
--- a/python/blackhole_sim/radiative_transfer.py
+++ b/python/blackhole_sim/radiative_transfer.py
@@ -11,17 +11,24 @@
 
 from dataclasses import dataclass
 import math
-from typing import Protocol
+from typing import Literal, Protocol
 
 import numpy as np
 
 from .grmhd import FluidSample, GRMHDSnapshot
 from .kerr import KerrTraceResult
 from .physics import C_SI
+from .synchrotron import magnetic_field_strength_code
 
 
 class CoefficientModel(Protocol):
-    def coefficients(self, sample: FluidSample, nu_emit_hz: float, p_cov: np.ndarray) -> tuple[float, float, np.ndarray]:
+    def coefficients(
+        self,
+        sample: FluidSample,
+        nu_emit_hz: float,
+        p_cov: np.ndarray,
+        spin_a: float | None = None,
+    ) -> tuple[float, float, np.ndarray]:
         """Return emission j_nu, absorption alpha_nu, and RGB source color."""
 
 
@@ -42,10 +49,16 @@ class ThermalSynchrotronFit:
     spectral_index: float = 1.2
     color_temperature_bias: float = 0.18
 
-    def coefficients(self, sample: FluidSample, nu_emit_hz: float, p_cov: np.ndarray) -> tuple[float, float, np.ndarray]:
+    def coefficients(
+        self,
+        sample: FluidSample,
+        nu_emit_hz: float,
+        p_cov: np.ndarray,
+        spin_a: float | None = None,
+    ) -> tuple[float, float, np.ndarray]:
         if not sample.valid or sample.rho <= 0.0 or sample.theta_e <= 0.0:
             return 0.0, 0.0, np.zeros(3)
-        b_mag = _magnetic_magnitude_proxy(sample.b_con)
+        b_mag = magnetic_field_strength_code(sample, spin_a=spin_a)
         theta_e = max(sample.theta_e, 1.0e-8)
         # Critical-frequency proxy in dimensionless field units. Real datasets
         # must provide the physical B and density scale for absolute flux.
@@ -69,6 +82,7 @@ class TransferConfig:
     max_optical_depth: float = 18.0
     min_redshift: float = 1.0e-6
     intensity_floor: float = 0.0
+    physics_mode: Literal["educational_proxy", "validated"] = "educational_proxy"
 
 
 @dataclass(frozen=True)
@@ -81,13 +95,6 @@ class TransferResult:
     redshift_max: float
 
 
-def _magnetic_magnitude_proxy(b_con: np.ndarray) -> float:
-    # Positive definite proxy for local field strength. A production adapter can
-    # replace this with b^2 = b_mu b^mu in the local tetrad frame.
-    b = np.asarray(b_con, dtype=float)
-    return float(np.linalg.norm(b[1:]) + 1.0e-30)
-
-
 def photon_energy_in_fluid_frame(p_cov: np.ndarray, u_con: np.ndarray) -> float:
     return -float(np.asarray(p_cov, dtype=float) @ np.asarray(u_con, dtype=float))
 
@@ -99,6 +106,21 @@ def invariant_redshift(p_cov: np.ndarray, u_emit: np.ndarray, observed_energy: f
     return observed_energy / e_emit
 
 
+def _coefficient_values(
+    model: CoefficientModel,
+    sample: FluidSample,
+    nu_emit_hz: float,
+    p_cov: np.ndarray,
+    spin_a: float,
+) -> tuple[float, float, np.ndarray]:
+    try:
+        return model.coefficients(sample, nu_emit_hz, p_cov, spin_a=spin_a)
+    except TypeError as exc:
+        if "spin_a" not in str(exc):
+            raise
+        return model.coefficients(sample, nu_emit_hz, p_cov)  # type: ignore[call-arg]
+
+
 def integrate_kerr_grrt(
     trace: KerrTraceResult,
     snapshot: GRMHDSnapshot,
@@ -114,6 +136,8 @@ def integrate_kerr_grrt(
     """
     model = coeffs or ThermalSynchrotronFit()
     tc = cfg or TransferConfig()
+    if tc.physics_mode == "validated" and isinstance(model, ThermalSynchrotronFit):
+        raise ValueError("validated transfer mode requires a non-proxy CoefficientModel")
     states = trace.states
     if len(states) < 2:
         return TransferResult(np.zeros(3), 0.0, 0, 0, math.inf, 0.0)
@@ -143,7 +167,7 @@ def integrate_kerr_grrt(
         if g_shift < tc.min_redshift:
             continue
         nu_emit = tc.observing_frequency_hz / g_shift
-        j, alpha, color = model.coefficients(sample, nu_emit, p_cov)
+        j, alpha, color = _coefficient_values(model, sample, nu_emit, p_cov, float(snapshot.spin_a))
         if j <= 0.0 and alpha <= 0.0:
             continue
         dlambda = float(np.linalg.norm(cur[1:4] - prev[1:4]))
diff --git a/python/blackhole_sim/synchrotron.py b/python/blackhole_sim/synchrotron.py
index c74779a..83c65f2 100644
--- a/python/blackhole_sim/synchrotron.py
+++ b/python/blackhole_sim/synchrotron.py
@@ -22,6 +22,7 @@
 from .calibration import C_CGS, E_CHARGE_ESU, K_BOLTZMANN_CGS, M_ELECTRON_CGS
 from .calibration import PhysicalScaling
 from .grmhd import FluidSample
+from .kerr import kerr_metric_covariant, zamo_tetrad
 
 
 @dataclass(frozen=True)
@@ -101,19 +102,72 @@ def scaled(self, factor: float) -> "PolarizedCoefficients":
         )
 
 
-def local_plasma_from_sample(sample: FluidSample, scaling: PhysicalScaling, p_cov: np.ndarray | None = None) -> LocalPlasmaFrame:
-    b_vec = np.asarray(sample.b_con[1:], dtype=float)
-    b_code = float(np.linalg.norm(b_vec))
+def magnetic_field_strength_code(sample: FluidSample, spin_a: float | None = None) -> float:
+    """Return local magnetic-field magnitude in code units.
+
+    When ``spin_a`` is available this uses the Kerr metric invariant
+    ``sqrt(b_mu b^mu)`` for the supplied contravariant magnetic four-vector.
+    The Euclidean spatial norm is retained only as a compatibility fallback for
+    callers that do not yet have a metric context.
+    """
+
+    b_con = np.asarray(sample.b_con, dtype=float)
+    if spin_a is not None and sample.valid:
+        try:
+            g_cov = kerr_metric_covariant(float(sample.r), float(sample.theta), float(spin_a))
+            b2 = float(b_con @ g_cov @ b_con)
+            if np.isfinite(b2) and b2 > 0.0:
+                return math.sqrt(b2)
+        except (FloatingPointError, ValueError, OverflowError):
+            pass
+    return float(np.linalg.norm(b_con[1:]))
+
+
+def magnetic_pitch_cosine(sample: FluidSample, p_cov: np.ndarray, spin_a: float | None = None) -> float:
+    """Return cos(angle) between photon direction and magnetic field.
+
+    The validated path uses the invariant fluid-frame relation
+    ``cos(alpha) = (p_mu b^mu) / (E_fluid |B|)`` where
+    ``E_fluid = -p_mu u^mu``. This avoids coordinate-basis Euclidean dot
+    products, which are not physically meaningful in Kerr coordinates.
+    """
+
+    b_code = magnetic_field_strength_code(sample, spin_a=spin_a)
+    if b_code <= 0.0:
+        return 0.0
+    p = np.asarray(p_cov, dtype=float)
+    u = np.asarray(sample.u_con, dtype=float)
+    e_fluid = -float(p @ u)
+    if e_fluid <= 0.0 or not np.isfinite(e_fluid):
+        return 0.0
+    numerator = float(p @ np.asarray(sample.b_con, dtype=float))
+    return float(np.clip(numerator / max(e_fluid * b_code, 1.0e-300), -1.0, 1.0))
+
+
+def _zamo_spatial_components(sample: FluidSample, spin_a: float | None) -> np.ndarray:
+    b_con = np.asarray(sample.b_con, dtype=float)
+    if spin_a is None or not sample.valid:
+        return b_con[1:].copy()
+    try:
+        g_cov = kerr_metric_covariant(float(sample.r), float(sample.theta), float(spin_a))
+        tetrad = zamo_tetrad(float(sample.r), float(sample.theta), float(spin_a))
+        return np.array([float(b_con @ g_cov @ tetrad[i]) for i in (1, 2, 3)], dtype=float)
+    except (FloatingPointError, ValueError, OverflowError):
+        return b_con[1:].copy()
+
+
+def local_plasma_from_sample(
+    sample: FluidSample,
+    scaling: PhysicalScaling,
+    p_cov: np.ndarray | None = None,
+    spin_a: float | None = None,
+) -> LocalPlasmaFrame:
+    b_hat = _zamo_spatial_components(sample, spin_a)
+    b_code = magnetic_field_strength_code(sample, spin_a=spin_a)
     b_gauss = float(scaling.magnetic_field_gauss(b_code))
     n_e = float(scaling.electron_number_density_cm3(sample.rho))
-    cos_los = 0.0
-    if p_cov is not None and b_code > 0.0:
-        # Coordinate-basis proxy for angle between B and photon spatial covector.
-        k = np.asarray(p_cov[1:], dtype=float)
-        nk = float(np.linalg.norm(k))
-        if nk > 0.0:
-            cos_los = float(np.clip(np.dot(b_vec, k) / (b_code * nk), -1.0, 1.0))
-    evpa = math.atan2(float(b_vec[2]), float(b_vec[1])) if b_code > 0.0 else 0.0
+    cos_los = magnetic_pitch_cosine(sample, p_cov, spin_a=spin_a) if p_cov is not None else 0.0
+    evpa = math.atan2(float(b_hat[2]), float(b_hat[1])) if b_code > 0.0 else 0.0
     return LocalPlasmaFrame(n_e, max(float(sample.theta_e), 1.0e-8), max(b_gauss, 0.0), cos_los, evpa)
 
 
diff --git a/python/tests/test_kerr.py b/python/tests/test_kerr.py
index 86c275a..86276e1 100644
--- a/python/tests/test_kerr.py
+++ b/python/tests/test_kerr.py
@@ -12,18 +12,36 @@
     kerr_metric_covariant,
     kerr_metric_contravariant,
     orbital_period_code,
+    static_limit_radius,
     trace_kerr_null_geodesic,
+    zamo_tetrad,
 )
 
 
 def test_kerr_horizon_and_isco_limits():
     assert np.isclose(horizon_radius(0.0), 2.0)
     assert np.isclose(horizon_radius(0.9), 1.0 + math.sqrt(1.0 - 0.9**2))
+    assert np.isclose(static_limit_radius(math.pi / 2.0, 0.9), 2.0)
+    assert np.isclose(static_limit_radius(0.0, 0.9), horizon_radius(0.9))
     assert np.isclose(isco_radius(0.0, prograde=True), 6.0)
     assert isco_radius(0.9, prograde=True) < 3.0
     assert isco_radius(0.9, prograde=False) > 8.0
 
 
+def test_kerr_metric_reduces_to_schwarzschild_at_zero_spin():
+    r, th, a = 11.0, 1.2, 0.0
+    f = 1.0 - 2.0 / r
+    s2 = math.sin(th) ** 2
+
+    gcov = kerr_metric_covariant(r, th, a)
+    expected_cov = np.diag([-f, 1.0 / f, r * r, r * r * s2])
+    assert np.allclose(gcov, expected_cov, atol=1e-12)
+
+    gcon = kerr_metric_contravariant(r, th, a)
+    expected_con = np.diag([-1.0 / f, f, 1.0 / (r * r), 1.0 / (r * r * s2)])
+    assert np.allclose(gcon, expected_con, atol=1e-12)
+
+
 def test_metric_inverse_identity():
     r, th, a = 8.0, 1.1, 0.72
     gcov = kerr_metric_covariant(r, th, a)
@@ -31,6 +49,14 @@ def test_metric_inverse_identity():
     assert np.allclose(gcov @ gcon, np.eye(4), atol=1e-11)
 
 
+def test_zamo_tetrad_is_orthonormal():
+    r, th, a = 12.0, 1.0, 0.7
+    tetrad = zamo_tetrad(r, th, a)
+    gcov = kerr_metric_covariant(r, th, a)
+    gram = tetrad @ gcov @ tetrad.T
+    assert np.allclose(gram, np.diag([-1.0, 1.0, 1.0, 1.0]), atol=1e-11)
+
+
 def test_camera_ray_is_null():
     cam = LocalCamera.from_degrees(r=40.0, inclination_degrees=62.0, fov_y_degrees=30.0)
     y0, pt, pphi = camera_ray_initial_state(cam, a=0.6, ndc_x=0.1, ndc_y=-0.2, aspect=16/9)
diff --git a/python/tests/test_native_kernel_assets.py b/python/tests/test_native_kernel_assets.py
index 4776a05..c359715 100644
--- a/python/tests/test_native_kernel_assets.py
+++ b/python/tests/test_native_kernel_assets.py
@@ -55,3 +55,22 @@ def test_webgpu_stokes_invalid_samples_use_finite_shader_values():
     assert "@compute @workgroup_size(8, 8, 1)\nfn main" in text
     assert "@compute @workgroup_size(8, 8, 1)\nfn kerr_stokes_render_kernel" not in text
     assert "fn kerr_stokes_render_kernel(gid: vec3<u32>)" in text
+
+
+def test_gpu_kernel_assets_use_zamo_camera_launch_not_hardcoded_momenta():
+    root = Path(__file__).resolve().parents[2]
+    expected = [
+        root / "native/cuda/kerr_stokes_kernel.cu",
+        root / "native/metal/kerr_stokes_kernel.metal",
+        root / "native/opencl/kerr_stokes_kernel.cl",
+        root / "native/rocm/kerr_stokes_kernel.hip",
+        root / "web/webgpu/src/stokes_brick_compute.wgsl",
+    ]
+    stale_tokens = ("p_t=-1", "p_t = -1", "p_phi=.12", "p_phi = 0.12", "-0.22")
+    for path in expected:
+        text = path.read_text()
+        assert "zamo_camera_ray_initial_state" in text, f"ZAMO launch missing from {path}"
+        assert "RayLaunch" in text, f"conserved-momentum launch payload missing from {path}"
+        assert "lapse" in text and "omega" in text, f"ZAMO lapse/frame-dragging terms missing from {path}"
+        for tok in stale_tokens:
+            assert tok not in text, f"stale hardcoded ray momentum token {tok!r} remains in {path}"
diff --git a/python/tests/test_radiative_transfer.py b/python/tests/test_radiative_transfer.py
index 3150c25..3004c08 100644
--- a/python/tests/test_radiative_transfer.py
+++ b/python/tests/test_radiative_transfer.py
@@ -1,4 +1,5 @@
 import numpy as np
+import pytest
 
 from blackhole_sim.grmhd import generate_analytic_grmhd_torus
 from blackhole_sim.kerr import LocalCamera, camera_ray_initial_state, trace_kerr_null_geodesic
@@ -27,3 +28,12 @@ def test_grrt_integrator_accumulates_nonnegative_intensity():
     assert result.valid_steps >= 0
     assert result.optical_depth >= 0.0
     assert np.all(result.intensity_rgb >= 0.0)
+
+
+def test_validated_mode_rejects_default_proxy_coefficients():
+    snap = generate_analytic_grmhd_torus(spin_a=0.4, nr=8, ntheta=6, nphi=5)
+    cam = LocalCamera.from_degrees(r=35, inclination_degrees=62, fov_y_degrees=28)
+    y0, pt, pph = camera_ray_initial_state(cam, snap.spin_a, 0.0, 0.0, 16 / 9)
+    trace = trace_kerr_null_geodesic(y0, pt, pph, snap.spin_a, step=0.12, max_steps=8, escape_radius=80)
+    with pytest.raises(ValueError, match="validated transfer mode"):
+        integrate_kerr_grrt(trace, snap, cfg=TransferConfig(physics_mode="validated"))
diff --git a/python/tests/test_synchrotron_polarized_transfer.py b/python/tests/test_synchrotron_polarized_transfer.py
index 3c67afb..8c33097 100644
--- a/python/tests/test_synchrotron_polarized_transfer.py
+++ b/python/tests/test_synchrotron_polarized_transfer.py
@@ -1,6 +1,16 @@
 import numpy as np
 
-from blackhole_sim.synchrotron import LocalPlasmaFrame, ThermalSynchrotronCoefficients, NonthermalPowerLawSynchrotronCoefficients, PolarizedCoefficients
+from blackhole_sim.calibration import PhysicalScaling
+from blackhole_sim.grmhd import FluidSample
+from blackhole_sim.synchrotron import (
+    LocalPlasmaFrame,
+    NonthermalPowerLawSynchrotronCoefficients,
+    PolarizedCoefficients,
+    ThermalSynchrotronCoefficients,
+    local_plasma_from_sample,
+    magnetic_field_strength_code,
+    magnetic_pitch_cosine,
+)
 from blackhole_sim.polarized_transfer import stokes_step_exact
 
 
@@ -22,3 +32,61 @@ def test_faraday_rotation_converts_q_to_u():
     assert np.isclose(out[0], 0.0, atol=1e-12)
     assert abs(out[2]) > 0.05
     assert np.isclose(np.linalg.norm(out[1:3]), 1.0, atol=1e-12)
+
+
+def test_metric_magnetic_strength_replaces_coordinate_norm():
+    r = 10.0
+    f = 1.0 - 2.0 / r
+    sample = FluidSample(
+        r=r,
+        theta=np.pi / 2.0,
+        phi=0.0,
+        rho=1.0,
+        theta_e=10.0,
+        pressure=1.0,
+        b_con=np.array([0.0, np.sqrt(f), 0.0, 0.0]),
+        u_con=np.array([1.0 / np.sqrt(f), 0.0, 0.0, 0.0]),
+        valid=True,
+    )
+    assert np.linalg.norm(sample.b_con[1:]) != 1.0
+    assert np.isclose(magnetic_field_strength_code(sample, spin_a=0.0), 1.0, atol=1e-12)
+
+
+def test_pitch_angle_uses_invariant_fluid_frame_projection():
+    r = 10.0
+    f = 1.0 - 2.0 / r
+    sample = FluidSample(
+        r=r,
+        theta=np.pi / 2.0,
+        phi=0.0,
+        rho=1.0,
+        theta_e=10.0,
+        pressure=1.0,
+        b_con=np.array([0.0, np.sqrt(f), 0.0, 0.0]),
+        u_con=np.array([1.0 / np.sqrt(f), 0.0, 0.0, 0.0]),
+        valid=True,
+    )
+    radial_photon = np.array([-np.sqrt(f), 1.0 / np.sqrt(f), 0.0, 0.0])
+    theta_photon = np.array([-np.sqrt(f), 0.0, r, 0.0])
+    assert np.isclose(magnetic_pitch_cosine(sample, radial_photon, spin_a=0.0), 1.0, atol=1e-12)
+    assert np.isclose(magnetic_pitch_cosine(sample, theta_photon, spin_a=0.0), 0.0, atol=1e-12)
+
+
+def test_local_plasma_frame_applies_metric_field_before_scaling():
+    r = 10.0
+    f = 1.0 - 2.0 / r
+    sample = FluidSample(
+        r=r,
+        theta=np.pi / 2.0,
+        phi=0.0,
+        rho=2.0,
+        theta_e=4.0,
+        pressure=1.0,
+        b_con=np.array([0.0, np.sqrt(f), 0.0, 0.0]),
+        u_con=np.array([1.0 / np.sqrt(f), 0.0, 0.0, 0.0]),
+        valid=True,
+    )
+    scaling = PhysicalScaling(1.0, 1.0, rho_cgs_per_code=2.0, b_gauss_per_code=7.0)
+    frame = local_plasma_from_sample(sample, scaling, spin_a=0.0)
+    assert np.isclose(frame.b_gauss, 7.0, atol=1e-12)
+    assert frame.n_e_cm3 > 0.0
diff --git a/web/webgpu/src/stokes_brick_compute.wgsl b/web/webgpu/src/stokes_brick_compute.wgsl
index d72c8de..0ef3953 100644
--- a/web/webgpu/src/stokes_brick_compute.wgsl
+++ b/web/webgpu/src/stokes_brick_compute.wgsl
@@ -20,6 +20,7 @@ struct RenderParams {
   pad0: u32,
 };
 struct State6 { t: f32, r: f32, th: f32, ph: f32, pr: f32, pth: f32 };
+struct RayLaunch { y: State6, p_t: f32, p_phi: f32 };
 struct BrickSample {
   valid: f32,
   coeffs: array<f32, 11>,
@@ -53,6 +54,40 @@ fn metric_contravariant(r: f32, theta: f32, a: f32) -> mat4x4<f32> {
   return g;
 }
 
+fn zamo_camera_ray_initial_state(px: u32, py: u32) -> RayLaunch {
+  let ndcx = 2.0 * ((f32(px) + 0.5) / f32(params.width)) - 1.0;
+  let ndcy = 1.0 - 2.0 * ((f32(py) + 0.5) / f32(params.height));
+  let aspect = f32(params.width) / max(f32(params.height), 1.0);
+  let tan_y = tan(0.5 * params.fov_y);
+  let n = normalize(vec3<f32>(-1.0, -ndcy * tan_y, ndcx * aspect * tan_y));
+
+  let r = params.camera_r;
+  let th = clamp_theta(params.camera_theta);
+  let a = params.spin_a;
+  let ct = cos(th);
+  let st = sin(th);
+  let s2 = max(st * st, 1.0e-8);
+  let sig = r * r + a * a * ct * ct;
+  let dlt = r * r - 2.0 * r + a * a;
+  let A = (r * r + a * a) * (r * r + a * a) - a * a * dlt * s2;
+  let gtt = -(1.0 - 2.0 * r / sig);
+  let gtphi = -2.0 * a * r * s2 / sig;
+  let grr = sig / dlt;
+  let gthth = sig;
+  let gphiphi = A * s2 / sig;
+  let lapse = sqrt(sig * dlt / A);
+  let omega = 2.0 * a * r / A;
+
+  let pcon_t = 1.0 / lapse;
+  let pcon_r = n.x * sqrt(dlt / sig);
+  let pcon_th = n.y / sqrt(sig);
+  let pcon_ph = omega / lapse + n.z / sqrt(gphiphi);
+  let p_t = gtt * pcon_t + gtphi * pcon_ph;
+  let p_phi = gtphi * pcon_t + gphiphi * pcon_ph;
+  let y = State6(0.0, r, th, 0.0, grr * pcon_r, gthth * pcon_th);
+  return RayLaunch(y, p_t, p_phi);
+}
+
 fn metric_derivative_r(r: f32, th: f32, a: f32) -> mat4x4<f32> {
   let e = 1.0e-3;
   return (metric_contravariant(r + e, th, a) - metric_contravariant(r - e, th, a)) * (1.0 / (2.0 * e));
@@ -152,18 +187,16 @@ fn stokes_step_rk2(S: vec4<f32>, c: array<f32, 11>, ds: f32) -> vec4<f32> {
 }
 
 fn ray_initial_state(px: u32, py: u32) -> State6 {
-  let ndcx = 2.0 * ((f32(px) + 0.5) / f32(params.width)) - 1.0;
-  let ndcy = 1.0 - 2.0 * ((f32(py) + 0.5) / f32(params.height));
-  return State6(0.0, params.camera_r, params.camera_theta, 0.0, -1.0, -0.22 * ndcy);
+  return zamo_camera_ray_initial_state(px, py).y;
 }
 
 fn kerr_stokes_render_kernel(gid: vec3<u32>) {
   if (gid.x >= params.width || gid.y >= params.height) { return; }
   let pix = gid.y * params.width + gid.x;
-  let ndcx = 2.0 * ((f32(gid.x) + 0.5) / f32(params.width)) - 1.0;
-  var y = ray_initial_state(gid.x, gid.y);
-  let p_t = -1.0;
-  let p_phi = 0.12 * ndcx;
+  let launch = zamo_camera_ray_initial_state(gid.x, gid.y);
+  var y = launch.y;
+  let p_t = launch.p_t;
+  let p_phi = launch.p_phi;
   let r_plus = 1.0 + sqrt(max(1.0 - params.spin_a * params.spin_a, 0.0));
   var S = vec4<f32>(0.0);
   let maxSteps = min(params.max_steps, 520u);