From ca7f70284bbe6441f797671d85b7a5027cb58d6d Mon Sep 17 00:00:00 2001
From: Angel Wang <angel@angelcube.dev>
Date: Mon, 27 Jul 2026 15:19:41 -0600
Subject: [PATCH] feat(apply): physical model, basic functions

---
 TODO.md             |  6 ++---
 src/apply/colour.ts | 62 +++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 65 insertions(+), 3 deletions(-)
 create mode 100644 src/apply/colour.ts

diff --git a/TODO.md b/TODO.md
index 06950f5..65fb069 100644
--- a/TODO.md
+++ b/TODO.md
@@ -41,12 +41,12 @@ it's a small code change and noticeably more robust across sun angles.
 
 ### 1. Port the physical model (reuse, don't reinvent)
 
-- [ ] Port constants: Rayleigh/Mie/ozone coefficients, scale heights,
+- [x] Port constants: Rayleigh/Mie/ozone coefficients, scale heights,
       ground/atmosphere radii, `groundAlbedo` — same as `gradient.ts` /
       the `common` shader tab.
 - [ ] Port `getScatteringValues` (density + extinction at a given altitude).
-- [ ] Port the phase functions (`rayleighPhase`, `miePhase`).
-- [ ] Port `rayIntersectSphere` / `intersectSphere`.
+- [x] Port the phase functions (`rayleighPhase`, `miePhase`).
+- [x] Port `rayIntersectSphere` / `intersectSphere`.
 - [ ] Decide: **do you need a full LUT texture, or a direct function call?**
       Since you're only evaluating a handful of directions per frame (not a
       whole screen), you likely don't need GPU textures at all — a plain
diff --git a/src/apply/colour.ts b/src/apply/colour.ts
new file mode 100644
index 0000000..63c6ad1
--- /dev/null
+++ b/src/apply/colour.ts
@@ -0,0 +1,62 @@
+import { Colordx, colordx } from "@colordx/core"
+
+// constants from Hillaire, section 4 table 1
+const RAYLEIGH_SCATTER = [5.802e-6, 13.558e-6, 33.1e-6] // rgb
+const MIE_SCATTER = 3.996e-6
+const MIE_ABSORB = 4.44e-6
+const OZONE_SCATTER = 0
+const OZONE_ABSORB = [0.65e-6, 1.881e-6, 0.085e-6] // rgb
+
+const GROUND_ALBEDO = 0.3
+
+const RAYLEIGH_SCALE_HEIGHT = 8e3 // metres
+const MIE_SCALE_HEIGHT = 1.2e3 // metres
+
+const GROUND_RADIUS = 6_360e3 // radius of the earth, metres
+const TOP_RADIUS = 6_460e3 // radius of the atmosphere, metres
+
+type Vec3 = [number, number, number];
+
+function dot(v1: Vec3, v2: Vec3) {
+  return v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2];
+}
+
+// ACES tonemapper (Knarkowicz)
+function aces(colour: Colordx): Colordx {
+	Object.keys(colour.toRgb()).forEach((_key, colour) => {
+		const n = colour * (2.51 * colour + 0.03)
+		const d = colour * (2.43 * colour + 0.59) + 0.14
+		return Math.max(0, Math.min(1, n / d))
+	})
+	return colordx(colour)
+}
+
+function rayleighPhase(angle: number): number {
+  return (3 * (1 + Math.cos(angle) ** 2)) / (16 * Math.PI);
+}
+
+function miePhase(angle: number): number {
+  const g = 0.8;
+  const scale = 3 / (8 * Math.PI);
+  const num = (1 - g ** 2) * (1 + Math.cos(angle) ** 2);
+  const denom =
+    (2 + g ** 2) * (1 + g ** 2 - 2 * g * Math.cos(angle)) ** (3 / 2);
+  return (scale * num) / denom;
+}
+
+// From "5.3.2 Intersecting Ray or Segment Against Sphere" (Real-Time Collision Detection)
+function intersectSphere(rayOrigin: Vec3, direction: Vec3, radius: number): number | null {
+  // Sphere center is at origin (= 0), so rayOrigin - sphereCenter = rayOrigin
+  const b = dot(rayOrigin, direction);
+  const c = dot(rayOrigin, rayOrigin) - (radius ** 2);
+  const discr = (b ** 2) - c;
+  if (discr < 0)
+    // Ray misses sphere
+    return null;
+  const t = -b - Math.sqrt(discr);
+  if (t < 0)
+    // Ray inside sphere; Use far discriminant
+    return -b + Math.sqrt(discr);
+  return t;
+}
+
-- 
2.51.2