feat(physics): Implement SDF-based physics engine and BVH

Completed Task #49.
- Implemented CPU-side SDF library (sphere, box, torus, plane).
- Implemented Dynamic BVH construction (rebuilt every frame).
- Implemented PhysicsSystem with semi-implicit Euler integration and collision resolution.
- Added visual debugging for BVH nodes.
- Created test_3d_physics interactive test and test_physics unit tests.
- Updated project docs and triaged new tasks.

13 files changed, 923 insertions, 2 deletions

diff --git a/src/3d/bvh.cc b/src/3d/bvh.cc
new file mode 100644
index 0000000..0c6bf9a
--- /dev/null
+++ b/src/3d/bvh.cc
@@ -0,0 +1,154 @@
+// This file is part of the 64k demo project.
+// It implements BVH construction and traversal.
+
+#include "3d/bvh.h"
+#include <algorithm>
+
+namespace {
+
+struct ObjectInfo {
+  int index;
+  AABB aabb;
+  vec3 centroid;
+};
+
+AABB get_world_aabb(const Object3D& obj) {
+  BoundingVolume local = obj.get_local_bounds();
+  mat4 model = obj.get_model_matrix();
+
+  vec3 corners[8] = {
+      {local.min.x, local.min.y, local.min.z},
+      {local.max.x, local.min.y, local.min.z},
+      {local.min.x, local.max.y, local.min.z},
+      {local.max.x, local.max.y, local.min.z},
+      {local.min.x, local.min.y, local.max.z},
+      {local.max.x, local.min.y, local.max.z},
+      {local.min.x, local.max.y, local.max.z},
+      {local.max.x, local.max.y, local.max.z},
+  };
+
+  AABB world;
+  for (int i = 0; i < 8; ++i) {
+    vec4 p = model * vec4(corners[i].x, corners[i].y, corners[i].z, 1.0f);
+    world.expand(p.xyz());
+  }
+  return world;
+}
+
+int build_recursive(std::vector<BVHNode>& nodes,
+                    std::vector<ObjectInfo>& obj_info, int start, int end) {
+  int node_idx = (int)nodes.size();
+  nodes.emplace_back();
+
+  AABB bounds;
+  for (int i = start; i < end; ++i) {
+    bounds.expand(obj_info[i].aabb);
+  }
+
+  int count = end - start;
+  if (count == 1) {
+    // Leaf node
+    nodes[node_idx].min_x = bounds.min.x;
+    nodes[node_idx].min_y = bounds.min.y;
+    nodes[node_idx].min_z = bounds.min.z;
+    nodes[node_idx].left_idx = -1;
+
+    nodes[node_idx].max_x = bounds.max.x;
+    nodes[node_idx].max_y = bounds.max.y;
+    nodes[node_idx].max_z = bounds.max.z;
+    nodes[node_idx].right_idx = obj_info[start].index;
+  } else {
+    // Internal node
+    // Find axis with largest variance (or just largest extent of centroids)
+    AABB centroid_bounds;
+    for (int i = start; i < end; ++i) {
+      centroid_bounds.expand(obj_info[i].centroid);
+    }
+
+    vec3 extent = centroid_bounds.max - centroid_bounds.min;
+    int axis = 0;
+    if (extent.y > extent.x)
+      axis = 1;
+    if (extent.z > (axis == 0 ? extent.x : extent.y))
+      axis = 2;
+
+    float split = (centroid_bounds.min[axis] + centroid_bounds.max[axis]) * 0.5f;
+
+    // Partition
+    int mid = start;
+    for (int i = start; i < end; ++i) {
+      if (obj_info[i].centroid[axis] < split) {
+        std::swap(obj_info[i], obj_info[mid]);
+        mid++;
+      }
+    }
+
+    // Fallback if partition failed
+    if (mid == start || mid == end) {
+      mid = start + count / 2;
+    }
+
+    int left = build_recursive(nodes, obj_info, start, mid);
+    int right = build_recursive(nodes, obj_info, mid, end);
+
+    nodes[node_idx].min_x = bounds.min.x;
+    nodes[node_idx].min_y = bounds.min.y;
+    nodes[node_idx].min_z = bounds.min.z;
+    nodes[node_idx].left_idx = left;
+
+    nodes[node_idx].max_x = bounds.max.x;
+    nodes[node_idx].max_y = bounds.max.y;
+    nodes[node_idx].max_z = bounds.max.z;
+    nodes[node_idx].right_idx = right;
+  }
+
+  return node_idx;
+}
+
+} // namespace
+
+void BVHBuilder::build(BVH& out_bvh, const std::vector<Object3D>& objects) {
+  out_bvh.nodes.clear();
+  if (objects.empty())
+    return;
+
+  std::vector<ObjectInfo> obj_info;
+  for (int i = 0; i < (int)objects.size(); ++i) {
+    if (objects[i].type == ObjectType::SKYBOX)
+      continue;
+    AABB aabb = get_world_aabb(objects[i]);
+    obj_info.push_back({i, aabb, aabb.center()});
+  }
+
+  if (obj_info.empty())
+    return;
+
+  out_bvh.nodes.reserve(obj_info.size() * 2);
+  build_recursive(out_bvh.nodes, obj_info, 0, (int)obj_info.size());
+}
+
+void BVH::query(const AABB& box, std::vector<int>& out_indices) const {
+  if (nodes.empty())
+    return;
+
+  std::vector<int> stack;
+  stack.push_back(0);
+
+  while (!stack.empty()) {
+    int idx = stack.back();
+    stack.pop_back();
+
+    const BVHNode& node = nodes[idx];
+    AABB node_aabb({node.min_x, node.min_y, node.min_z},
+                   {node.max_x, node.max_y, node.max_z});
+
+    if (node_aabb.intersects(box)) {
+      if (node.left_idx < 0) {
+        out_indices.push_back(node.right_idx);
+      } else {
+        stack.push_back(node.left_idx);
+        stack.push_back(node.right_idx);
+      }
+    }
+  }
+}
diff --git a/src/3d/bvh.h b/src/3d/bvh.h
new file mode 100644
index 0000000..97e9a06
--- /dev/null
+++ b/src/3d/bvh.h
@@ -0,0 +1,69 @@
+// This file is part of the 64k demo project.
+// It defines the Bounding Volume Hierarchy (BVH) for scene acceleration.
+
+#pragma once
+
+#include "3d/object.h"
+#include "util/mini_math.h"
+#include <cstdint>
+#include <vector>
+
+struct BVHNode {
+  float min_x, min_y, min_z;
+  int32_t left_idx; // If < 0, this is a leaf node.
+
+  float max_x, max_y, max_z;
+  int32_t right_idx; // If leaf, this holds the object_index.
+};
+
+struct AABB {
+  vec3 min;
+  vec3 max;
+
+  AABB() : min(1e10f, 1e10f, 1e10f), max(-1e10f, -1e10f, -1e10f) {
+  }
+  AABB(vec3 min, vec3 max) : min(min), max(max) {
+  }
+
+  void expand(vec3 p) {
+    if (p.x < min.x)
+      min.x = p.x;
+    if (p.y < min.y)
+      min.y = p.y;
+    if (p.z < min.z)
+      min.z = p.z;
+    if (p.x > max.x)
+      max.x = p.x;
+    if (p.y > max.y)
+      max.y = p.y;
+    if (p.z > max.z)
+      max.z = p.z;
+  }
+
+  void expand(const AABB& other) {
+    expand(other.min);
+    expand(other.max);
+  }
+
+  bool intersects(const AABB& other) const {
+    return (min.x <= other.max.x && max.x >= other.min.x) &&
+           (min.y <= other.max.y && max.y >= other.min.y) &&
+           (min.z <= other.max.z && max.z >= other.min.z);
+  }
+
+  vec3 center() const {
+    return (min + max) * 0.5f;
+  }
+};
+
+class BVH {
+ public:
+  std::vector<BVHNode> nodes;
+
+  void query(const AABB& box, std::vector<int>& out_indices) const;
+};
+
+class BVHBuilder {
+ public:
+  static void build(BVH& out_bvh, const std::vector<Object3D>& objects);
+};
diff --git a/src/3d/object.h b/src/3d/object.h
index 2099a5c..6d21393 100644
--- a/src/3d/object.h
+++ b/src/3d/object.h
@@ -31,9 +31,16 @@ class Object3D {
   // Material parameters could go here (color, roughness, etc.)
   vec4 color;
 
+  // Physics fields
+  vec3 velocity;
+  float mass;
+  float restitution;
+  bool is_static;
+
   Object3D(ObjectType t = ObjectType::CUBE)
       : position(0, 0, 0), rotation(0, 0, 0, 1), scale(1, 1, 1), type(t),
-        color(1, 1, 1, 1) {
+        color(1, 1, 1, 1), velocity(0, 0, 0), mass(1.0f), restitution(0.5f),
+        is_static(false) {
   }
 
   mat4 get_model_matrix() const {
@@ -47,6 +54,8 @@ class Object3D {
   // Returns the local-space AABB of the primitive (before transform)
   // Used to generate the proxy geometry for rasterization.
   BoundingVolume get_local_bounds() const {
+    if (type == ObjectType::TORUS)
+      return {{-1.5f, -0.5f, -1.5f}, {1.5f, 0.5f, 1.5f}};
     // Simple defaults for unit primitives
     return {{-1, -1, -1}, {1, 1, 1}};
   }
diff --git a/src/3d/physics.cc b/src/3d/physics.cc
new file mode 100644
index 0000000..351dd06
--- /dev/null
+++ b/src/3d/physics.cc
@@ -0,0 +1,144 @@
+// This file is part of the 64k demo project.
+// It implements a lightweight SDF-based physics engine.
+
+#include "3d/physics.h"
+#include "3d/bvh.h"
+#include "3d/sdf_cpu.h"
+#include <algorithm>
+
+namespace {
+// Helper to get world AABB (copied from bvh.cc or shared)
+AABB get_world_aabb(const Object3D& obj) {
+  BoundingVolume local = obj.get_local_bounds();
+  mat4 model = obj.get_model_matrix();
+
+  vec3 corners[8] = {
+      {local.min.x, local.min.y, local.min.z},
+      {local.max.x, local.min.y, local.min.z},
+      {local.min.x, local.max.y, local.min.z},
+      {local.max.x, local.max.y, local.min.z},
+      {local.min.x, local.min.y, local.max.z},
+      {local.max.x, local.min.y, local.max.z},
+      {local.min.x, local.max.y, local.max.z},
+      {local.max.x, local.max.y, local.max.z},
+  };
+
+  AABB world;
+  for (int i = 0; i < 8; ++i) {
+    vec4 p = model * vec4(corners[i].x, corners[i].y, corners[i].z, 1.0f);
+    world.expand(p.xyz());
+  }
+  return world;
+}
+} // namespace
+
+float PhysicsSystem::sample_sdf(const Object3D& obj, vec3 world_p) {
+  mat4 inv_model = obj.get_model_matrix().inverse();
+  vec4 local_p4 = inv_model * vec4(world_p.x, world_p.y, world_p.z, 1.0f);
+  vec3 q = local_p4.xyz();
+
+  float d = 1000.0f;
+  if (obj.type == ObjectType::SPHERE) {
+    d = q.len() - 1.0f;
+  } else if (obj.type == ObjectType::BOX || obj.type == ObjectType::CUBE) {
+    d = sdf::sdBox(q, vec3(1.0f, 1.0f, 1.0f));
+  } else if (obj.type == ObjectType::TORUS) {
+    d = sdf::sdTorus(q, vec2(1.0f, 0.4f));
+  } else if (obj.type == ObjectType::PLANE) {
+    d = sdf::sdPlane(q, vec3(0.0f, 1.0f, 0.0f), 0.0f);
+  }
+
+  // Extract scale from model matrix (assuming orthogonal with uniform or
+  // non-uniform scale)
+  mat4 model = obj.get_model_matrix();
+  float sx = vec3(model.m[0], model.m[1], model.m[2]).len();
+  float sy = vec3(model.m[4], model.m[5], model.m[6]).len();
+  float sz = vec3(model.m[8], model.m[9], model.m[10]).len();
+  float s = std::min(sx, std::min(sy, sz));
+
+  return d * s;
+}
+
+void PhysicsSystem::resolve_collision(Object3D& a, Object3D& b) {
+  if (a.is_static && b.is_static)
+    return;
+
+  // Probe points for 'a' (center and corners)
+  BoundingVolume local = a.get_local_bounds();
+  mat4 model_a = a.get_model_matrix();
+  vec3 probes[9] = {
+      {0, 0, 0}, // Center
+      {local.min.x, local.min.y, local.min.z},
+      {local.max.x, local.min.y, local.min.z},
+      {local.min.x, local.max.y, local.min.z},
+      {local.max.x, local.max.y, local.min.z},
+      {local.min.x, local.min.y, local.max.z},
+      {local.max.x, local.min.y, local.max.z},
+      {local.min.x, local.max.y, local.max.z},
+      {local.max.x, local.max.y, local.max.z},
+  };
+
+  for (int i = 0; i < 9; ++i) {
+    vec3 world_probe =
+        (model_a * vec4(probes[i].x, probes[i].y, probes[i].z, 1.0f)).xyz();
+    float d = sample_sdf(b, world_probe);
+
+    if (d < 0.0f) {
+      // Collision detected!
+      float penetration = -d;
+
+      // Calculate normal via gradient of b's SDF
+      auto b_sdf = [this, &b](vec3 p) { return sample_sdf(b, p); };
+      vec3 normal = sdf::calc_normal(world_probe, b_sdf);
+
+      // Resolution
+      if (!a.is_static) {
+        // Positional correction
+        a.position += normal * penetration;
+
+        // Velocity response
+        float v_dot_n = vec3::dot(a.velocity, normal);
+        if (v_dot_n < 0) {
+          a.velocity -= normal * (1.0f + a.restitution) * v_dot_n;
+        }
+      }
+    }
+  }
+}
+
+void PhysicsSystem::update(Scene& scene, float dt) {
+  if (dt <= 0)
+    return;
+
+  // 1. Integration
+  for (auto& obj : scene.objects) {
+    if (obj.is_static)
+      continue;
+    obj.velocity += gravity * dt;
+    obj.position += obj.velocity * dt;
+  }
+
+  // 2. Broad Phase
+  BVH bvh;
+  BVHBuilder::build(bvh, scene.objects);
+
+  // 3. Narrow Phase & Resolution
+  // We do multiple iterations for better stability? Just 1 for now.
+  for (int iter = 0; iter < 2; ++iter) {
+    for (int i = 0; i < (int)scene.objects.size(); ++i) {
+      Object3D& a = scene.objects[i];
+      if (a.is_static)
+        continue;
+
+      AABB query_box = get_world_aabb(a);
+      std::vector<int> candidates;
+      bvh.query(query_box, candidates);
+
+      for (int cand_idx : candidates) {
+        if (cand_idx == i)
+          continue;
+        resolve_collision(a, scene.objects[cand_idx]);
+      }
+    }
+  }
+}
diff --git a/src/3d/physics.h b/src/3d/physics.h
new file mode 100644
index 0000000..a014acd
--- /dev/null
+++ b/src/3d/physics.h
@@ -0,0 +1,18 @@
+// This file is part of the 64k demo project.
+// It implements a lightweight SDF-based physics engine.
+
+#pragma once
+
+#include "3d/scene.h"
+#include "util/mini_math.h"
+
+class PhysicsSystem {
+ public:
+  vec3 gravity = {0.0f, -9.81f, 0.0f};
+
+  void update(Scene& scene, float dt);
+
+ private:
+  void resolve_collision(Object3D& a, Object3D& b);
+  float sample_sdf(const Object3D& obj, vec3 world_p);
+};
diff --git a/src/3d/renderer.cc b/src/3d/renderer.cc
index e4320f4..64f6fc3 100644
--- a/src/3d/renderer.cc
+++ b/src/3d/renderer.cc
@@ -199,6 +199,13 @@ void Renderer3D::set_sky_texture(WGPUTextureView sky_view) {
   sky_texture_view_ = sky_view;
 }
 
+void Renderer3D::add_debug_aabb(const vec3& min, const vec3& max,
+                                const vec3& color) {
+#if !defined(STRIP_ALL)
+  visual_debug_.add_aabb(min, max, color);
+#endif
+}
+
 void Renderer3D::create_pipeline() {
   WGPUBindGroupLayoutEntry entries[5] = {};
   entries[0].binding = 0;
diff --git a/src/3d/renderer.h b/src/3d/renderer.h
index 57ad671..8068bdc 100644
--- a/src/3d/renderer.h
+++ b/src/3d/renderer.h
@@ -62,6 +62,8 @@ class Renderer3D {
   void set_noise_texture(WGPUTextureView noise_view);
   void set_sky_texture(WGPUTextureView sky_view);
 
+  void add_debug_aabb(const vec3& min, const vec3& max, const vec3& color);
+
   // Resize handler (if needed for internal buffers)
   void resize(int width, int height);
 
diff --git a/src/3d/sdf_cpu.h b/src/3d/sdf_cpu.h
new file mode 100644
index 0000000..1e461f6
--- /dev/null
+++ b/src/3d/sdf_cpu.h
@@ -0,0 +1,55 @@
+// This file is part of the 64k demo project.
+// It implements CPU-side Signed Distance Field (SDF) primitives
+// and utilities for physics and collision detection.
+
+#pragma once
+
+#include "util/mini_math.h"
+#include <algorithm>
+#include <cmath>
+
+namespace sdf {
+
+inline float sdSphere(vec3 p, float r) {
+  return p.len() - r;
+}
+
+inline float sdBox(vec3 p, vec3 b) {
+  vec3 q;
+  q.x = std::abs(p.x) - b.x;
+  q.y = std::abs(p.y) - b.y;
+  q.z = std::abs(p.z) - b.z;
+
+  vec3 max_q_0;
+  max_q_0.x = std::max(q.x, 0.0f);
+  max_q_0.y = std::max(q.y, 0.0f);
+  max_q_0.z = std::max(q.z, 0.0f);
+
+  return max_q_0.len() + std::min(std::max(q.x, std::max(q.y, q.z)), 0.0f);
+}
+
+inline float sdTorus(vec3 p, vec2 t) {
+  vec2 p_xz(p.x, p.z);
+  vec2 q(p_xz.len() - t.x, p.y);
+  return q.len() - t.y;
+}
+
+inline float sdPlane(vec3 p, vec3 n, float h) {
+  return vec3::dot(p, n) + h;
+}
+
+/**
+ * Calculates the normal of an SDF at point p using numerical differentiation.
+ * sdf_func must be a callable: float sdf_func(vec3 p)
+ */
+template <typename F>
+inline vec3 calc_normal(vec3 p, F sdf_func, float e = 0.001f) {
+  const float d2e = 2.0f * e;
+  return vec3(
+             sdf_func(vec3(p.x + e, p.y, p.z)) - sdf_func(vec3(p.x - e, p.y, p.z)),
+             sdf_func(vec3(p.x, p.y + e, p.z)) - sdf_func(vec3(p.x, p.y - e, p.z)),
+             sdf_func(vec3(p.x, p.y, p.z + e)) - sdf_func(vec3(p.x, p.y, p.z - e)))
+      .normalize();
+}
+
+} // namespace sdf
diff --git a/src/3d/visual_debug.cc b/src/3d/visual_debug.cc
index f8d9bed..ab4cb6c 100644
--- a/src/3d/visual_debug.cc
+++ b/src/3d/visual_debug.cc
@@ -182,6 +182,25 @@ void VisualDebug::add_box(const mat4& transform, const vec3& local_extent,
   }
 }
 
+void VisualDebug::add_aabb(const vec3& min, const vec3& max, const vec3& color) {
+  vec3 p[] = {{min.x, min.y, min.z}, {max.x, min.y, min.z}, {max.x, max.y, min.z},
+              {min.x, max.y, min.z}, {min.x, min.y, max.z}, {max.x, min.y, max.z},
+              {max.x, max.y, max.z}, {min.x, max.y, max.z}};
+
+  DebugLine edges[] = {
+      {p[0], p[1], color}, {p[1], p[2], color}, {p[2], p[3], color},
+      {p[3], p[0], color}, // Front
+      {p[4], p[5], color}, {p[5], p[6], color}, {p[6], p[7], color},
+      {p[7], p[4], color}, // Back
+      {p[0], p[4], color}, {p[1], p[5], color}, {p[2], p[6], color},
+      {p[3], p[7], color} // Connections
+  };
+
+  for (const auto& l : edges) {
+    lines_.push_back(l);
+  }
+}
+
 void VisualDebug::update_buffers(const mat4& view_proj) {
   // Update Uniforms
   wgpuQueueWriteBuffer(wgpuDeviceGetQueue(device_), uniform_buffer_, 0,
diff --git a/src/3d/visual_debug.h b/src/3d/visual_debug.h
index 456cb10..6173fc4 100644
--- a/src/3d/visual_debug.h
+++ b/src/3d/visual_debug.h
@@ -25,6 +25,8 @@ class VisualDebug {
   void add_box(const mat4& transform, const vec3& local_extent,
                const vec3& color);
 
+  void add_aabb(const vec3& min, const vec3& max, const vec3& color);
+
   // Render all queued primitives and clear the queue
   void render(WGPURenderPassEncoder pass, const mat4& view_proj);
 
diff --git a/src/tests/test_3d_physics.cc b/src/tests/test_3d_physics.cc
new file mode 100644
index 0000000..6d7f476
--- /dev/null
+++ b/src/tests/test_3d_physics.cc
@@ -0,0 +1,291 @@
+// This file is part of the 64k demo project.
+// Standalone "mini-demo" for testing the 3D physics engine.
+
+#include "3d/camera.h"
+#include "3d/object.h"
+#include "3d/renderer.h"
+#include "3d/scene.h"
+#include "3d/bvh.h"
+#include "3d/physics.h"
+#include "gpu/effects/shaders.h"
+#include "gpu/texture_manager.h"
+#include "platform.h"
+#include "procedural/generator.h"
+#include <cmath>
+#include <cstdio>
+#include <cstring>
+#include <vector>
+
+// Global State
+static Renderer3D g_renderer;
+static TextureManager g_textures;
+static Scene g_scene;
+static Camera g_camera;
+static PhysicsSystem g_physics;
+static WGPUDevice g_device = nullptr;
+static WGPUQueue g_queue = nullptr;
+static WGPUSurface g_surface = nullptr;
+static WGPUAdapter g_adapter = nullptr;
+static WGPUTextureFormat g_format = WGPUTextureFormat_Undefined;
+
+// ... (init_wgpu implementation same as before)
+void init_wgpu(PlatformState* platform_state) {
+  WGPUInstance instance = wgpuCreateInstance(nullptr);
+  if (!instance) {
+    fprintf(stderr, "Failed to create WGPU instance.\n");
+    exit(1);
+  }
+
+  g_surface = platform_create_wgpu_surface(instance, platform_state);
+  if (!g_surface) {
+    fprintf(stderr, "Failed to create WGPU surface.\n");
+    exit(1);
+  }
+
+  WGPURequestAdapterOptions adapter_opts = {};
+  adapter_opts.compatibleSurface = g_surface;
+  adapter_opts.powerPreference = WGPUPowerPreference_HighPerformance;
+
+#if defined(DEMO_CROSS_COMPILE_WIN32)
+  auto on_adapter = [](WGPURequestAdapterStatus status, WGPUAdapter adapter,
+                       const char* message, void* userdata) {
+    if (status == WGPURequestAdapterStatus_Success) {
+      *(WGPUAdapter*)userdata = adapter;
+    }
+  };
+  wgpuInstanceRequestAdapter(instance, &adapter_opts, on_adapter, &g_adapter);
+#else
+  auto on_adapter = [](WGPURequestAdapterStatus status, WGPUAdapter adapter,
+                       WGPUStringView message, void* userdata, void* user2) {
+    (void)user2;
+    if (status == WGPURequestAdapterStatus_Success) {
+      *(WGPUAdapter*)userdata = adapter;
+    }
+  };
+  WGPURequestAdapterCallbackInfo adapter_cb = {};
+  adapter_cb.mode = WGPUCallbackMode_WaitAnyOnly;
+  adapter_cb.callback = on_adapter;
+  adapter_cb.userdata1 = &g_adapter;
+  wgpuInstanceRequestAdapter(instance, &adapter_opts, adapter_cb);
+#endif
+
+  while (!g_adapter) {
+    platform_wgpu_wait_any(instance);
+  }
+
+  WGPUDeviceDescriptor device_desc = {};
+
+#if defined(DEMO_CROSS_COMPILE_WIN32)
+  auto on_device = [](WGPURequestDeviceStatus status, WGPUDevice device,
+                      const char* message, void* userdata) {
+    if (status == WGPURequestDeviceStatus_Success) {
+      *(WGPUDevice*)userdata = device;
+    }
+  };
+  wgpuAdapterRequestDevice(g_adapter, &device_desc, on_device, &g_device);
+#else
+  auto on_device = [](WGPURequestDeviceStatus status, WGPUDevice device,
+                      WGPUStringView message, void* userdata, void* user2) {
+    (void)user2;
+    if (status == WGPURequestDeviceStatus_Success) {
+      *(WGPUDevice*)userdata = device;
+    }
+  };
+  WGPURequestDeviceCallbackInfo device_cb = {};
+  device_cb.mode = WGPUCallbackMode_WaitAnyOnly;
+  device_cb.callback = on_device;
+  device_cb.userdata1 = &g_device;
+  wgpuAdapterRequestDevice(g_adapter, &device_desc, device_cb);
+#endif
+
+  while (!g_device) {
+    platform_wgpu_wait_any(instance);
+  }
+
+  g_queue = wgpuDeviceGetQueue(g_device);
+
+  WGPUSurfaceCapabilities caps = {};
+  wgpuSurfaceGetCapabilities(g_surface, g_adapter, &caps);
+  g_format = caps.formats[0];
+
+  WGPUSurfaceConfiguration config = {};
+  config.device = g_device;
+  config.format = g_format;
+  config.usage = WGPUTextureUsage_RenderAttachment;
+  config.width = platform_state->width;
+  config.height = platform_state->height;
+  config.presentMode = WGPUPresentMode_Fifo;
+  config.alphaMode = WGPUCompositeAlphaMode_Opaque;
+  wgpuSurfaceConfigure(g_surface, &config);
+}
+
+void setup_scene() {
+  g_scene.clear();
+  srand(12345); // Fixed seed
+
+  // Large floor, use BOX type (SDF) at index 0
+  Object3D floor(ObjectType::BOX);
+  floor.position = vec3(0, -2.0f, 0);
+  floor.scale = vec3(25.0f, 0.2f, 25.0f);
+  floor.color = vec4(0.8f, 0.8f, 0.8f, 1.0f);
+  floor.is_static = true;
+  g_scene.add_object(floor);
+
+  // Large center Torus (SDF)
+  Object3D center(ObjectType::TORUS);
+  center.position = vec3(0, 1.0f, 0);
+  center.scale = vec3(2.5f, 2.5f, 2.5f);
+  center.color = vec4(1, 0.2, 0.2, 1);
+  center.is_static = false;
+  center.restitution = 0.8f;
+  g_scene.add_object(center);
+
+  // Moving Sphere (SDF)
+  Object3D sphere(ObjectType::SPHERE);
+  sphere.position = vec3(4.0f, 2.0f, 0);
+  sphere.scale = vec3(1.5f, 1.5f, 1.5f);
+  sphere.color = vec4(0.2, 1, 0.2, 1);
+  sphere.is_static = false;
+  sphere.velocity = vec3(-2.0f, 5.0f, 1.0f);
+  g_scene.add_object(sphere);
+
+  // Random objects
+  for (int i = 0; i < 30; ++i) {
+    ObjectType type = ObjectType::SPHERE;
+    int r = rand() % 3;
+    if (r == 1)
+      type = ObjectType::TORUS;
+    if (r == 2)
+      type = ObjectType::BOX;
+
+    Object3D obj(type);
+    float angle = (rand() % 360) * 0.01745f;
+    float dist = 3.0f + (rand() % 100) * 0.05f;
+    float height = 5.0f + (rand() % 100) * 0.04f; // Start higher
+    obj.position = vec3(std::cos(angle) * dist, height, std::sin(angle) * dist);
+
+    // Random non-uniform scale for debugging
+    float s = 0.6f + (rand() % 100) * 0.008f;
+    obj.scale = vec3(s, s * 1.2f, s * 0.8f);
+
+    obj.color = vec4((rand() % 100) / 100.0f, (rand() % 100) / 100.0f,
+                     (rand() % 100) / 100.0f, 1.0f);
+    obj.is_static = false;
+    obj.velocity = vec3((rand() % 100 - 50) * 0.01f, 0, (rand() % 100 - 50) * 0.01f);
+    g_scene.add_object(obj);
+  }
+}
+
+// Wrapper to generate periodic noise
+bool gen_periodic_noise(uint8_t* buffer, int w, int h, const float* params,
+                        int num_params) {
+  if (!procedural::gen_noise(buffer, w, h, params, num_params))
+    return false;
+  float p_params[] = {0.1f}; // 10% overlap
+  return procedural::make_periodic(buffer, w, h, p_params, 1);
+}
+
+int main(int argc, char** argv) {
+  printf("Running 3D Physics Test...\n");
+
+#if !defined(STRIP_ALL)
+  for (int i = 1; i < argc; ++i) {
+    if (strcmp(argv[i], "--debug") == 0) {
+      Renderer3D::SetDebugEnabled(true);
+    }
+  }
+#else
+  (void)argc;
+  (void)argv;
+#endif
+
+  PlatformState platform_state = platform_init(false, 1280, 720);
+
+  // The test's own WGPU init sequence
+  init_wgpu(&platform_state);
+
+  InitShaderComposer();
+
+  g_renderer.init(g_device, g_queue, g_format);
+  g_renderer.resize(platform_state.width, platform_state.height);
+
+  g_textures.init(g_device, g_queue);
+  ProceduralTextureDef noise_def;
+  noise_def.width = 256;
+  noise_def.height = 256;
+  noise_def.gen_func = gen_periodic_noise;
+  noise_def.params.push_back(1234.0f);
+  noise_def.params.push_back(16.0f);
+  g_textures.create_procedural_texture("noise", noise_def);
+
+  g_renderer.set_noise_texture(g_textures.get_texture_view("noise"));
+
+  ProceduralTextureDef sky_def;
+  sky_def.width = 512;
+  sky_def.height = 256;
+  sky_def.gen_func = procedural::gen_perlin;
+  sky_def.params = {42.0f, 4.0f, 1.0f, 0.5f, 6.0f};
+  g_textures.create_procedural_texture("sky", sky_def);
+
+  g_renderer.set_sky_texture(g_textures.get_texture_view("sky"));
+
+  setup_scene();
+
+  g_camera.position = vec3(0, 5, 10);
+  g_camera.target = vec3(0, 0, 0);
+
+  while (!platform_should_close(&platform_state)) {
+    platform_poll(&platform_state);
+    float time = (float)platform_state.time;
+
+    float cam_radius = 10.0f + std::sin(time * 0.3f) * 4.0f;
+    float cam_height = 5.0f + std::cos(time * 0.4f) * 3.0f;
+    g_camera.set_look_at(vec3(std::sin(time * 0.5f) * cam_radius, cam_height,
+                              std::cos(time * 0.5f) * cam_radius),
+                         vec3(0, 0, 0), vec3(0, 1, 0));
+    g_camera.aspect_ratio = platform_state.aspect_ratio;
+
+    static double last_time = 0;
+    float dt = (float)(platform_state.time - last_time);
+    if (dt > 0.1f) dt = 0.1f; // Cap dt for stability
+    last_time = platform_state.time;
+
+    g_physics.update(g_scene, dt);
+
+    BVH bvh;
+    BVHBuilder::build(bvh, g_scene.objects);
+    for (const auto& node : bvh.nodes) {
+      g_renderer.add_debug_aabb({node.min_x, node.min_y, node.min_z},
+                                {node.max_x, node.max_y, node.max_z},
+                                {0.0f, 1.0f, 0.0f});
+    }
+
+#if !defined(STRIP_ALL)
+    Renderer3D::SetDebugEnabled(true);
+#endif
+
+    WGPUSurfaceTexture surface_tex;
+    wgpuSurfaceGetCurrentTexture(g_surface, &surface_tex);
+    if (surface_tex.status ==
+        WGPUSurfaceGetCurrentTextureStatus_SuccessOptimal) {
+      const WGPUTextureViewDescriptor view_desc = {
+          .format = g_format,
+          .dimension = WGPUTextureViewDimension_2D,
+          .mipLevelCount = 1,
+          .arrayLayerCount = 1,
+      };
+
+      const WGPUTextureView view =
+          wgpuTextureCreateView(surface_tex.texture, &view_desc);
+      g_renderer.render(g_scene, g_camera, time, view);
+      wgpuTextureViewRelease(view);
+      wgpuSurfacePresent(g_surface);
+      wgpuTextureRelease(surface_tex.texture);
+    }
+  }
+
+  g_renderer.shutdown();
+  g_textures.shutdown();
+  platform_shutdown(&platform_state);
+  return 0;
+}
\ No newline at end of file
diff --git a/src/tests/test_3d_render.cc b/src/tests/test_3d_render.cc
index 2e9b663..002cb55 100644
--- a/src/tests/test_3d_render.cc
+++ b/src/tests/test_3d_render.cc
@@ -276,4 +276,4 @@ int main(int argc, char** argv) {
   g_textures.shutdown();
   platform_shutdown(&platform_state);
   return 0;
-}
+}
\ No newline at end of file
diff --git a/src/tests/test_physics.cc b/src/tests/test_physics.cc
new file mode 100644
index 0000000..a59502c
--- /dev/null
+++ b/src/tests/test_physics.cc
@@ -0,0 +1,151 @@
+// This file is part of the 64k demo project.
+// It tests the CPU-side SDF library and BVH for physics and collision.
+
+#include "3d/bvh.h"
+#include "3d/physics.h"
+#include "3d/sdf_cpu.h"
+#include <cassert>
+#include <cmath>
+#include <iostream>
+
+bool near(float a, float b, float e = 0.001f) {
+  return std::abs(a - b) < e;
+}
+
+void test_sdf_sphere() {
+  std::cout << "Testing sdSphere..." << std::endl;
+  float r = 1.0f;
+  assert(near(sdf::sdSphere({0, 0, 0}, r), -1.0f));
+  assert(near(sdf::sdSphere({1, 0, 0}, r), 0.0f));
+  assert(near(sdf::sdSphere({2, 0, 0}, r), 1.0f));
+}
+
+void test_sdf_box() {
+  std::cout << "Testing sdBox..." << std::endl;
+  vec3 b(1, 1, 1);
+  assert(near(sdf::sdBox({0, 0, 0}, b), -1.0f));
+  assert(near(sdf::sdBox({1, 1, 1}, b), 0.0f));
+  assert(near(sdf::sdBox({2, 0, 0}, b), 1.0f));
+}
+
+void test_sdf_torus() {
+  std::cout << "Testing sdTorus..." << std::endl;
+  vec2 t(1.0f, 0.2f);
+  // Point on the ring: length(p.xz) = 1.0, p.y = 0
+  assert(near(sdf::sdTorus({1, 0, 0}, t), -0.2f));
+  assert(near(sdf::sdTorus({1.2f, 0, 0}, t), 0.0f));
+}
+
+void test_sdf_plane() {
+  std::cout << "Testing sdPlane..." << std::endl;
+  vec3 n(0, 1, 0);
+  float h = 1.0f; // Plane is at y = -1 (dot(p,n) + 1 = 0 => y = -1)
+  assert(near(sdf::sdPlane({0, 0, 0}, n, h), 1.0f));
+  assert(near(sdf::sdPlane({0, -1, 0}, n, h), 0.0f));
+}
+
+void test_calc_normal() {
+  std::cout << "Testing calc_normal..." << std::endl;
+
+  // Sphere normal at (1,0,0) should be (1,0,0)
+  auto sphere_sdf = [](vec3 p) { return sdf::sdSphere(p, 1.0f); };
+  vec3 n = sdf::calc_normal({1, 0, 0}, sphere_sdf);
+  assert(near(n.x, 1.0f) && near(n.y, 0.0f) && near(n.z, 0.0f));
+
+  // Box normal at side
+  auto box_sdf = [](vec3 p) { return sdf::sdBox(p, {1, 1, 1}); };
+  n = sdf::calc_normal({1, 0, 0}, box_sdf);
+  assert(near(n.x, 1.0f) && near(n.y, 0.0f) && near(n.z, 0.0f));
+
+  // Plane normal should be n
+  vec3 plane_n(0, 1, 0);
+  auto plane_sdf = [plane_n](vec3 p) { return sdf::sdPlane(p, plane_n, 1.0f); };
+  n = sdf::calc_normal({0, 0, 0}, plane_sdf);
+  assert(near(n.x, plane_n.x) && near(n.y, plane_n.y) && near(n.z, plane_n.z));
+}
+
+void test_bvh() {
+  std::cout << "Testing BVH..." << std::endl;
+  std::vector<Object3D> objects;
+
+  // Object 0: Left side
+  Object3D obj0(ObjectType::BOX);
+  obj0.position = {-10, 0, 0};
+  objects.push_back(obj0);
+
+  // Object 1: Right side
+  Object3D obj1(ObjectType::BOX);
+  obj1.position = {10, 0, 0};
+  objects.push_back(obj1);
+
+  BVH bvh;
+  BVHBuilder::build(bvh, objects);
+
+  assert(bvh.nodes.size() == 3); // 1 root + 2 leaves
+
+  // Query left side
+  std::vector<int> results;
+  bvh.query({{-12, -2, -2}, {-8, 2, 2}}, results);
+  assert(results.size() == 1);
+  assert(results[0] == 0);
+
+  // Query right side
+  results.clear();
+  bvh.query({{8, -2, -2}, {12, 2, 2}}, results);
+  assert(results.size() == 1);
+  assert(results[0] == 1);
+
+  // Query center (should miss both)
+  results.clear();
+  bvh.query({{-2, -2, -2}, {2, 2, 2}}, results);
+  assert(results.size() == 0);
+
+  // Query both
+    results.clear();
+    bvh.query({{-12, -2, -2}, {12, 2, 2}}, results);
+    assert(results.size() == 2);
+  }
+  
+  void test_physics_falling() {
+    std::cout << "Testing Physics falling..." << std::endl;
+    Scene scene;
+  
+    // Plane at y = -1
+    Object3D plane(ObjectType::PLANE);
+    plane.position = {0, -1, 0};
+    plane.is_static = true;
+    scene.add_object(plane);
+  
+    // Sphere at y = 5
+    Object3D sphere(ObjectType::SPHERE);
+    sphere.position = {0, 5, 0};
+    sphere.velocity = {0, 0, 0};
+    sphere.restitution = 0.0f; // No bounce for simple test
+    scene.add_object(sphere);
+  
+    PhysicsSystem physics;
+    float dt = 0.016f;
+    for (int i = 0; i < 100; ++i) {
+      physics.update(scene, dt);
+    }
+  
+    // Sphere should be above or at plane (y >= 0 because sphere radius is 1,
+    // plane is at -1)
+    assert(scene.objects[1].position.y >= -0.01f);
+    // Also should have slowed down
+    assert(scene.objects[1].velocity.y > -1.0f);
+  }
+  
+  int main() {
+    test_sdf_sphere();
+    test_sdf_box();
+    test_sdf_torus();
+    test_sdf_plane();
+    test_calc_normal();
+    test_bvh();
+    test_physics_falling();
+  
+    std::cout << "--- ALL PHYSICS TESTS PASSED ---" << std::endl;
+    return 0;
+  }
+  
\ No newline at end of file