diff --git a/.gitignore b/.gitignore
index cb46182..ac8a990 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,4 +1,6 @@
 Manifest.toml
 **/*.vtk
 **/*.stl
-**/*.csv
\ No newline at end of file
+**/*.csv
+**/*.txt
+**/*.gz
\ No newline at end of file
diff --git a/Project.toml b/Project.toml
index 5d1ea0c..9e4dd9a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -3,6 +3,7 @@ authors = ["Freddie Barter <fjbarter@outlook.com>"]
 uuid = "e89368bc-bc31-4f60-a9df-df669e2c49e6"
 
 [deps]
+AcceleratedKernels = "6a4ca0a5-0e36-4168-a932-d9be78d558f1"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
diff --git a/README.md b/README.md
index 1489dd6..2a707c3 100644
--- a/README.md
+++ b/README.md
@@ -10,13 +10,14 @@ using Pkg
 Pkg.develop(url="https://github.com/fjbarter/PRICK.jl")
 ```
 
-Then you can simply load in your sphere and triangle-wall data, define a start point and direction as SVectors, as below:
+Then you can load sphere and triangle-wall data, define one start point, and trace a batch of directions:
 ```julia
 using PRICK
 
 # Load mesh and tag each surface as a mirror or a sink
 vessel_tm = TriangleMesh("RAMCylinder.stl"; units=u"m")
 surfaces = [mirror(vessel_tm)]  # or sink(vessel_tm)
+surface_bvh = build_surface_bvh(surfaces)  # build once and reuse
 
 # Load spheres (example using Packing3D VTK reader)
 data = read_vtk_file("particles_0.vtk")
@@ -27,24 +28,19 @@ radii = Float64.(r)
 # Build the spheres together and construct BVH
 spheres = build_sphere_bvh(X, radii)
 
-# Trace a ray
-p0 = SVector(0.0079, 0.0, 0.022)
-d = SVector(1.0, 0.0, 0.0)
-res = trace_ray_geometric(p0, d, surfaces, spheres)
-```
-
-## Visualisation (optional)
-
-`visualise_trace` requires GLMakie to be loaded by the user:
-
-```julia
-using GLMakie
-fig = visualise_trace(
-    res;
-    X=X, radii=radii,
-    vessel_tm=vessel_tm
+# Trace N rays from a single start point p0 with a (3, N) direction matrix
+p0 = (0.0079, 0.0, 0.022)
+N = 10_000
+D = randn(3, N)  # directions are normalized internally
+
+res = trace_rays(
+    p0,
+    D,
+    surface_bvh,
+    spheres;
+    max_bounces=20_000,
+    max_length=1e6,
 )
-display(fig)
 ```
 
 ## Notes
@@ -52,3 +48,4 @@ display(fig)
 - `mirror` surfaces reflect specularly.
 - `sink` surfaces terminate rays.
 - The union mesh bounding box is used as an unmeshed sink and is expanded by 10%.
+- Batched tracing returns per-ray `escaped`, `total_length`, `nsteps`, and `termination`.
diff --git a/prototyping/cylinder/Project.toml b/prototyping/cylinder/Project.toml
index 021381b..d186a62 100644
--- a/prototyping/cylinder/Project.toml
+++ b/prototyping/cylinder/Project.toml
@@ -1,3 +1,4 @@
 [deps]
-GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+PProf = "e4faabce-9ead-11e9-39d9-4379958e3056"
 PRICK = "e89368bc-bc31-4f60-a9df-df669e2c49e6"
diff --git a/prototyping/cylinder/allocation_flamegraph.jl b/prototyping/cylinder/allocation_flamegraph.jl
new file mode 100644
index 0000000..bbc8867
--- /dev/null
+++ b/prototyping/cylinder/allocation_flamegraph.jl
@@ -0,0 +1,71 @@
+
+using Profile
+using PRICK
+
+using PProf
+
+function random_unit_directions(n::Integer)
+    n > 0 || throw(ArgumentError("n must be > 0"))
+    D = Matrix{Float64}(undef, 3, n)
+    @inbounds for i in 1:n
+        phi = 2pi * rand()
+        u = 2rand() - 1
+        s = sqrt(1 - u * u)
+        D[1, i] = s * cos(phi)
+        D[2, i] = s * sin(phi)
+        D[3, i] = u
+    end
+    return D
+end
+
+function main(;
+    n_rays::Int=1_000,
+    max_bounces::Int=20_000,
+    max_length::Float64=1e6,
+    sample_rate::Float64=0.01,
+    flat_report_path::String="alloc_profile_flat.txt",
+    tree_report_path::String="alloc_profile_tree.txt",
+)
+    println("Building geometry...")
+    vessel_trianglemesh = TriangleMesh("RAMCylinder.stl"; units=u"m")
+    surfaces = [sink(vessel_trianglemesh)]
+    surface_bvh = build_surface_bvh(surfaces)
+
+    println("Reading particles...")
+    data = read_vtk_file("particles_0.vtk")
+    x, y, z, r = retrieve_coordinates(data)
+    X = permutedims(Float64.(hcat(x, y, z)))
+    radii = Float64.(r)
+    spheres = build_sphere_bvh(X, radii)
+
+    p0 = (PRICK.mean(x), PRICK.mean(y), PRICK.mean(z))
+    D = random_unit_directions(n_rays)
+
+    println("Warmup run...")
+    trace_rays(p0, D, surface_bvh, spheres; max_bounces=max_bounces, max_length=max_length)
+
+    println("Profiling allocations (sample_rate=$sample_rate)...")
+    Profile.Allocs.clear()
+    Profile.Allocs.@profile sample_rate=sample_rate begin
+        trace_rays(p0, D, surface_bvh, spheres; max_bounces=max_bounces, max_length=max_length)
+    end
+
+    prof = Profile.Allocs.fetch()
+    println("Collected allocation samples: ", length(prof.allocs))
+
+    println("Writing text allocation reports...")
+    open(flat_report_path, "w") do io
+        Profile.Allocs.print(io, prof; format=:flat, C=false, maxdepth=12, mincount=1)
+    end
+    open(tree_report_path, "w") do io
+        Profile.Allocs.print(io, prof; format=:tree, C=false, maxdepth=20, mincount=1)
+    end
+    println("Wrote: ", abspath(flat_report_path))
+    println("Wrote: ", abspath(tree_report_path))
+
+    println("Opening pprof web UI...")
+    PProf.Allocs.pprof(prof; from_c=false)
+    println("If browser does not open automatically, use the URL printed by PProf.")
+end
+
+main()
diff --git a/prototyping/cylinder/cylinder.jl b/prototyping/cylinder/cylinder.jl
index ecc2725..95ace73 100644
--- a/prototyping/cylinder/cylinder.jl
+++ b/prototyping/cylinder/cylinder.jl
@@ -1,92 +1,63 @@
-using PRICK
-
 # ============================================================
-# example() demonstrating:
-# - provide X, radii
-# - provide vessel
-# - trace a ray and return results
+# Example batched ray tracing in a sphere bed
 # ============================================================
-function example()
-    # Vessel from stl file
-    vessel_trianglemesh = TriangleMesh("RAMCylinder.stl"; units=u"m")
-
-    # Construct surface (each is either a sink or mirror)
-    surfaces = [
-        sink(vessel_trianglemesh),
-    ]
-
-    # Spheres from VTK file
-    println("Reading file...")
-    file = "particles_0.vtk"
-    data = read_vtk_file(file)
-    x, y, z, r = retrieve_coordinates(data)
-
-    X = permutedims(Float64.(hcat(x, y, z)))   # 3×N
-    radii = Float64.(r)
-    spheres = build_sphere_bvh(X, radii)
-
-    # Build p0 from mean sphere position
-    p0x = PRICK.mean(x)
-    p0y = PRICK.mean(y)
-    p0z = PRICK.mean(z)
-    p0 = SVector{3,Float64}(p0x, p0y, p0z)
-
-    # p0 = SVector{3,Float64}(0.0079, 0.0, 0.022)
-
-    N = 1
-    path_lengths = zeros(Float64, N)
-    n_steps = zeros(Float64, N)
-    println("Tracing $N rays...")
-    for i in 1:N
-        # Fixed start point
-        
-        φ = 2π * rand()          # azimuth
-        u = 2rand() - 1          # u = cos(θ) uniformly in [-1,1]
-        s = sqrt(1 - u*u)
-        d = SVector(s*cos(φ), s*sin(φ), u)
-
-        res = trace_ray_geometric(p0, d, surfaces, spheres;
-                                  record_path=false,
-                                  max_steps=20000,
-                                  max_length=1e6)
-        
-        path_lengths[i] = res.escaped ? res.total_length : NaN
-        n_steps[i] = res.escaped ? res.nsteps : NaN
+using PRICK
+using BenchmarkTools
+
+function random_unit_directions(n::Integer)
+    n > 0 || throw(ArgumentError("n must be > 0"))
+    D = Matrix{Float64}(undef, 3, n)
+    @inbounds for i in 1:n
+        phi = 2pi * rand()
+        u = 2rand() - 1
+        s = sqrt(1 - u * u)
+        D[1, i] = s * cos(phi)
+        D[2, i] = s * sin(phi)
+        D[3, i] = u
     end
-
-    path_lengths_dense = filter(!isnan, path_lengths)
-
-    mean_length = PRICK.mean(path_lengths_dense)
-    println("Mean path length over $N rays (escaped: $(length(path_lengths_dense))): $(mean_length) m")
-
-    n_steps_dense = filter(!isnan, n_steps)
-    mean_number_of_steps = PRICK.mean(n_steps_dense)
-    println("Mean number of steps over escaped rays: $(mean_number_of_steps)")
-
-    # Single ray trace example
-
-    # Start point and direction
-    p0 = SVector{3,Float64}(0.0079, 0.0, 0.022)
-    φ = 2π * rand()          # azimuth
-        u = 2rand() - 1          # u = cos(θ) uniformly in [-1,1]
-        s = sqrt(1 - u*u)
-        d = SVector(s*cos(φ), s*sin(φ), u)
-
-    res = trace_ray_geometric(p0, d, surfaces, spheres;
-                              record_path=true,
-                              max_steps=20000,
-                              max_length=1e6)
-
-    println("escaped       = ", res.escaped)
-    println("total_length  = ", res.total_length)
-    println("nsteps        = ", res.nsteps)
-    println("path vertices = ", length(res.path))
-
-    # If you want to plot later:
-    P = path_to_matrix(res.path)  # 3xK
-    return (res=res, X=X, radii=r, vessel_tm=vessel_trianglemesh, path_matrix=P)
+    return D
 end
 
-out = example()
-# fig = visualise_trace(out.res; X=out.X, radii=out.radii, vessel_tm=out.vessel_tm)
-# display(fig)
+vessel_trianglemesh = TriangleMesh("RAMCylinder.stl"; units=u"m")
+surfaces = [sink(vessel_trianglemesh)]
+surface_bvh = build_surface_bvh(surfaces)
+
+println("Reading file...")
+data = read_vtk_file("particles_0.vtk")
+x, y, z, r = retrieve_coordinates(data)
+
+X = permutedims(Float64.(hcat(x, y, z)))
+radii = Float64.(r)
+spheres = build_sphere_bvh(X, radii)
+
+p0 = (PRICK.mean(x), PRICK.mean(y), PRICK.mean(z))
+n_rays = 1000
+D = random_unit_directions(n_rays)
+
+println("Tracing $n_rays rays...")
+t_trace0 = time()
+
+@benchmark trace_rays(
+    $p0,
+    $D,
+    $surface_bvh,
+    $spheres;
+    max_bounces=20_000,
+    max_length=1e6,
+)
+# t_trace = time() - t_trace0
+# println("Traced $n_rays rays in $t_trace seconds")
+
+# escaped_lengths = [res.total_length[i] for i in eachindex(res.total_length) if res.escaped[i]]
+# escaped_steps = [res.nsteps[i] for i in eachindex(res.nsteps) if res.escaped[i]]
+# escaped_term_dists = [res.p0_termination_distance[i] for i in eachindex(res.p0_termination_distance) if res.escaped[i]]
+
+# mean_length = isempty(escaped_lengths) ? NaN : PRICK.mean(escaped_lengths)
+# mean_steps = isempty(escaped_steps) ? NaN : PRICK.mean(escaped_steps)
+# mean_tortuosity = sum(escaped_lengths)/sum(escaped_term_dists)
+
+# println("Escaped rays: $(length(escaped_lengths)) / $n_rays")
+# println("Mean path length over escaped rays: $(mean_length) m")
+# println("Mean number of reflections over escaped rays: $(mean_steps)")
+# println("Mean tortuosity: $mean_tortuosity")
+
diff --git a/prototyping/polyh_bed/polyh_bed.jl b/prototyping/polyh_bed/polyh_bed.jl
index c4e2610..91291ce 100644
--- a/prototyping/polyh_bed/polyh_bed.jl
+++ b/prototyping/polyh_bed/polyh_bed.jl
@@ -1,24 +1,39 @@
 using PRICK
 using CSV
 using DataFrames
-using GLMakie
+using BenchmarkTools
+
+function random_unit_directions(n::Integer)
+    n > 0 || throw(ArgumentError("n must be > 0"))
+    D = Matrix{Float64}(undef, 3, n)
+    @inbounds for i in 1:n
+        phi = 2pi * rand()
+        u = 2rand() - 1
+        s = sqrt(1 - u * u)
+        D[1, i] = s * cos(phi)
+        D[2, i] = s * sin(phi)
+        D[3, i] = u
+    end
+    return D
+end
 
-function example2()
-    # Load geometry
-    container_tm = TriangleMesh("Extents.stl", units=u"m")
+function example(;
+    n_rays::Int=1_000,
+    max_bounces::Int=20_000,
+    max_length::Float64=1e6,
+)
+    println("Building container surface...")
+    container_tm = TriangleMesh("Extents.stl"; units=u"m")
     surfaces = [sink(container_tm)]
+    surface_bvh = build_surface_bvh(surfaces)
 
-    # Load particle mesh
+    println("Building particle template...")
     particle_tm = ParticleTriangleMesh("Particle.stl"; units=u"m")
 
-    # Load particle data
-    # Expected columns: 
-    #   x, y, z, r (radius)
-    #   rot_x, rot_y, rot_z (Euler angles)
+    println("Reading particle state...")
     df = CSV.read("ParticlePositions.csv", DataFrame)
     n_particles = nrow(df)
 
-    # Convert df Matrix & Vector
     X = Matrix{Float64}(undef, 3, n_particles)
     X[1, :] = df.x
     X[2, :] = df.y
@@ -30,58 +45,36 @@ function example2()
     orients[2, :] = df.rot_y
     orients[3, :] = df.rot_z
 
-    # Build polyhedral BVH
-    println("Building Polyhedral BVH...")
+    println("Building polyhedral BVH...")
     polyh_bvh = build_polyh_bvh(particle_tm, X, r, orients)
 
-    # Use a random restart hill-climb to find a void starting point near center
     p0 = find_void_rrhc(X, r)
+    D = random_unit_directions(n_rays)
 
-    # Random direction
-    φ = 2π * rand()
-    u = 2rand() - 1
-    s = sqrt(1 - u * u)
-    d = SVector(s * cos(φ), s * sin(φ), u)
-
-    # Run Trace
-    res = trace_ray_geometric(
+    println("Tracing $n_rays rays...")
+    res = trace_rays(
         p0,
-        d,
-        surfaces,
-        polyh_bvh
+        D,
+        surface_bvh,
+        polyh_bvh;
+        max_bounces=max_bounces,
+        max_length=max_length,
     )
 
-    # Results
-    println("--------------------------------")
-    println("Trace Result:")
-    println("  Escaped:      $(res.escaped)")
-    println("  Total Length: $(res.total_length) m")
-    println("  Steps:        $(res.nsteps)")
-    println("--------------------------------")
+    escaped_lengths = [res.total_length[i] for i in eachindex(res.total_length) if res.escaped[i]]
+    escaped_steps = [res.nsteps[i] for i in eachindex(res.nsteps) if res.escaped[i]]
+    escaped_term_dists = [res.p0_termination_distance[i] for i in eachindex(res.p0_termination_distance) if res.escaped[i]]
 
-    # Optional: Return data for visualization or analysis
-    P = isempty(res.path) ? Float64[] : PRICK.path_to_matrix(res.path)
+    mean_length = isempty(escaped_lengths) ? NaN : PRICK.mean(escaped_lengths)
+    mean_steps = isempty(escaped_steps) ? NaN : PRICK.mean(escaped_steps)
+    mean_tortuosity = isempty(escaped_term_dists) ? NaN : sum(escaped_lengths) / sum(escaped_term_dists)
 
-    return (
-        res=res,
-        path_matrix=P,
-        container=container_tm,
-        particles=polyh_bvh
-    )
-end
+    println("Escaped rays: $(length(escaped_lengths)) / $n_rays")
+    println("Mean path length over escaped rays: $(mean_length) m")
+    println("Mean number of reflections over escaped rays: $(mean_steps)")
+    println("Mean tortuosity: $mean_tortuosity")
 
-out = example2()
-
-fig = visualise_trace(
-    out.res,
-    out.particles;
-    vessel_tm=out.container,
-    show_vessel=true,
-    show_particles=false,
-)
-
-GLMakie.activate!()
-display(fig)
+    return res
+end
 
-println("Press Enter to exit...")
-readline()
\ No newline at end of file
+@benchmark example()
diff --git a/src/PRICK.jl b/src/PRICK.jl
index 3c42d61..cec536c 100644
--- a/src/PRICK.jl
+++ b/src/PRICK.jl
@@ -3,7 +3,9 @@
 module PRICK
 
 using StaticArrays
-import ImplicitBVH: BVH, BBox, BSphere, traverse_rays
+using ImplicitBVH
+import ImplicitBVH: BVH, BBox, BSphere
+import AcceleratedKernels as AK
 
 using Meshes
 using FileIO               # FileIO.load
@@ -15,6 +17,8 @@ using Unitful: ustrip, uconvert, m, @u_str
 include("utils.jl")
 include("mesh.jl")
 include("bvh.jl")
+include("ray_batch.jl")
+include("raytrace_active_lvt.jl")
 include("intersections.jl")
 include("tracer.jl")
 include("visualise.jl")
@@ -25,14 +29,11 @@ export @u_str
 # Packing3D
 export read_vtk_file, retrieve_coordinates
 
-# StaticArrays
-export SVector, SMatrix
-
 # PRICK
 export TriangleMesh, TriangleSurface, ParticleTriangleMesh, Mirror, Sink
 export mirror, sink
-export build_sphere_bvh, build_polyh_bvh
-export trace_ray_geometric, RayTraceResult
+export build_sphere_bvh, build_polyh_bvh, build_surface_bvh
+export trace_rays, RayTraceBatchResult, RayTermination
 export visualise_trace
 export path_to_matrix
 export find_void_rrhc
diff --git a/src/bvh.jl b/src/bvh.jl
index 70968c0..29ca917 100644
--- a/src/bvh.jl
+++ b/src/bvh.jl
@@ -3,22 +3,22 @@
 # ============================================================
 struct MeshBVH
     bvh::BVH
-    verts::Vector{SVector{3,Float64}}
+    verts::Vector{NTuple{3,Float64}}
     tris::Vector{NTuple{3,Int32}}
 end
 
 struct SurfaceBVH
     bvh::BVH
-    verts::Vector{SVector{3,Float64}}
+    verts::Vector{NTuple{3,Float64}}
     tris::Vector{NTuple{3,Int32}}
     kinds::Vector{SurfaceKind}
-    bbox_mins::SVector{3,Float64}
-    bbox_maxs::SVector{3,Float64}
+    bbox_mins::NTuple{3,Float64}
+    bbox_maxs::NTuple{3,Float64}
 end
 
 struct SphereBVH
     bvh::BVH
-    centres::Vector{SVector{3,Float64}}
+    centres::Vector{NTuple{3,Float64}}
     radii::Vector{Float64}
 end
 
@@ -27,25 +27,33 @@ struct PolyhedralBVH
     centres::Vector{SVector{3,Float64}}
     radii::Vector{Float64}
     orientations::Vector{SMatrix{3,3,Float64,9}}
+    scales::Vector{Float64}
+    local_verts::Vector{NTuple{3,Float64}}
+    local_tris::Vector{NTuple{3,Int32}}
+    local_normals::Vector{NTuple{3,Float64}}
     polyh_mesh::ParticleTriangleMesh
 end
 
 function build_mesh_bvh(tm::TriangleMesh)::MeshBVH
     V = tm.vertices
     C = tm.connectivity
-    verts = [SVector{3,Float64}(V[1, i], V[2, i], V[3, i]) for i in 1:size(V, 2)]
+    verts = [(V[1, i], V[2, i], V[3, i]) for i in 1:size(V, 2)]
     tris = [(Int32(C[1, j]), Int32(C[2, j]), Int32(C[3, j])) for j in 1:size(C, 2)]
 
-    tri_centres = Vector{SVector{3,Float64}}(undef, length(tris))
+    tri_centres = Vector{NTuple{3,Float64}}(undef, length(tris))
     tri_radii = Vector{Float64}(undef, length(tris))
     for (j, (i1, i2, i3)) in enumerate(tris)
         v0, v1, v2 = verts[i1], verts[i2], verts[i3]
-        c = (v0 + v1 + v2) / 3.0
+        c = (
+            (v0[1] + v1[1] + v2[1]) / 3.0,
+            (v0[2] + v1[2] + v2[2]) / 3.0,
+            (v0[3] + v1[3] + v2[3]) / 3.0,
+        )
         rtri = max(
-            sqrt(dot3(c - v0, c - v0)),
+            sqrt(dot3(sub3(c, v0), sub3(c, v0))),
             max(
-                sqrt(dot3(c - v1, c - v1)),
-                sqrt(dot3(c - v2, c - v2)),
+                sqrt(dot3(sub3(c, v1), sub3(c, v1))),
+                sqrt(dot3(sub3(c, v2), sub3(c, v2))),
             ),
         )
         tri_centres[j] = c
@@ -62,14 +70,23 @@ function build_surface_bvh(surfaces::AbstractVector{<:TriangleSurface})::Surface
 
     meshes = [s.mesh for s in surfaces]
     union_mesh = union_trianglemesh(meshes)
-    bbox_mins, bbox_maxs = compute_bounds(union_mesh)
-    center = 0.5 * (bbox_mins + bbox_maxs)
-    half = 0.5 * (bbox_maxs - bbox_mins)
-    half = 1.1 * half
-    bbox_mins = center - half
-    bbox_maxs = center + half
-
-    verts = Vector{SVector{3,Float64}}()
+    bbox_mins_raw, bbox_maxs_raw = compute_bounds(union_mesh)
+    mins = (Float64(bbox_mins_raw[1]), Float64(bbox_mins_raw[2]), Float64(bbox_mins_raw[3]))
+    maxs = (Float64(bbox_maxs_raw[1]), Float64(bbox_maxs_raw[2]), Float64(bbox_maxs_raw[3]))
+    center = (
+        0.5 * (mins[1] + maxs[1]),
+        0.5 * (mins[2] + maxs[2]),
+        0.5 * (mins[3] + maxs[3]),
+    )
+    half = (
+        0.55 * (maxs[1] - mins[1]),
+        0.55 * (maxs[2] - mins[2]),
+        0.55 * (maxs[3] - mins[3]),
+    )
+    bbox_mins = (center[1] - half[1], center[2] - half[2], center[3] - half[3])
+    bbox_maxs = (center[1] + half[1], center[2] + half[2], center[3] + half[3])
+
+    verts = Vector{NTuple{3,Float64}}()
     tris = Vector{NTuple{3,Int32}}()
     kinds = Vector{SurfaceKind}()
     nverts_total = 0
@@ -82,7 +99,7 @@ function build_surface_bvh(surfaces::AbstractVector{<:TriangleSurface})::Surface
         Nt = size(C, 2)
 
         @inbounds for i in 1:Nv
-            push!(verts, SVector{3,Float64}(V[1, i], V[2, i], V[3, i]))
+            push!(verts, (Float64(V[1, i]), Float64(V[2, i]), Float64(V[3, i])))
         end
 
         off = Int32(nverts_total)
@@ -94,16 +111,20 @@ function build_surface_bvh(surfaces::AbstractVector{<:TriangleSurface})::Surface
         nverts_total += Nv
     end
 
-    tri_centres = Vector{SVector{3,Float64}}(undef, length(tris))
+    tri_centres = Vector{NTuple{3,Float64}}(undef, length(tris))
     tri_radii = Vector{Float64}(undef, length(tris))
     @inbounds for (j, (i1, i2, i3)) in enumerate(tris)
         v0, v1, v2 = verts[i1], verts[i2], verts[i3]
-        c = (v0 + v1 + v2) / 3.0
+        c = (
+            (v0[1] + v1[1] + v2[1]) / 3.0,
+            (v0[2] + v1[2] + v2[2]) / 3.0,
+            (v0[3] + v1[3] + v2[3]) / 3.0,
+        )
         rtri = max(
-            sqrt(dot3(c - v0, c - v0)),
+            sqrt(dot3(sub3(c, v0), sub3(c, v0))),
             max(
-                sqrt(dot3(c - v1, c - v1)),
-                sqrt(dot3(c - v2, c - v2)),
+                sqrt(dot3(sub3(c, v1), sub3(c, v1))),
+                sqrt(dot3(sub3(c, v2), sub3(c, v2))),
             ),
         )
         tri_centres[j] = c
@@ -122,12 +143,12 @@ function build_sphere_bvh(X::AbstractMatrix{<:Real}, r::AbstractVector{<:Real}):
     n = size(X, 2)
     @assert length(r) == n "r must have length N"
 
-    centres = Vector{SVector{3,Float64}}(undef, n)
+    centres = Vector{NTuple{3,Float64}}(undef, n)
     radii = Vector{Float64}(undef, n)
     leaves = Vector{BSphere{Float64}}(undef, n)
 
     @inbounds for i in 1:n
-        c = SVector{3,Float64}(float(X[1, i]), float(X[2, i]), float(X[3, i]))
+        c = (float(X[1, i]), float(X[2, i]), float(X[3, i]))
         ri = float(r[i])
         centres[i] = c
         radii[i] = ri
@@ -155,11 +176,35 @@ function build_polyh_bvh(
     centres = Vector{SVector{3,Float64}}(undef, n)
     radii = Vector{Float64}(undef, n)
     orientations = Vector{SMatrix{3,3,Float64,9}}(undef, n)
+    scales = Vector{Float64}(undef, n)
     leaves = Vector{BBox{Float64}}(undef, n)
 
     # Pre-allocate bounds for loop efficiency
     V_local = polyh_mesh.vertices
+    C_local = polyh_mesh.connectivity
     Nv = size(V_local, 2)
+    Nt = size(C_local, 2)
+
+    local_verts = Vector{NTuple{3,Float64}}(undef, Nv)
+    @inbounds for i in 1:Nv
+        local_verts[i] = (Float64(V_local[1, i]), Float64(V_local[2, i]), Float64(V_local[3, i]))
+    end
+
+    local_tris = Vector{NTuple{3,Int32}}(undef, Nt)
+    local_normals = Vector{NTuple{3,Float64}}(undef, Nt)
+    @inbounds for j in 1:Nt
+        i1 = Int32(C_local[1, j])
+        i2 = Int32(C_local[2, j])
+        i3 = Int32(C_local[3, j])
+        local_tris[j] = (i1, i2, i3)
+
+        v0 = local_verts[i1]
+        v1 = local_verts[i2]
+        v2 = local_verts[i3]
+        nloc = cross3(sub3(v1, v0), sub3(v2, v0))
+        invn = inv(sqrt(dot3(nloc, nloc)))
+        local_normals[j] = mul3(nloc, invn)
+    end
 
     @inbounds for i in 1:n
         c = SVector{3,Float64}(float(X[1, i]), float(X[2, i]), float(X[3, i]))
@@ -172,6 +217,7 @@ function build_polyh_bvh(
         r_mat = get_rotation_matrix(view(orients, :, i), units=u"°")
         orientations[i] = r_mat
         sf = get_scale_factor(polyh_mesh, Float64(ri))
+        scales[i] = sf
 
         # Initialize bounds with the first transformed vertex
         v1 = SVector(V_local[1, 1], V_local[2, 1], V_local[3, 1])
@@ -192,5 +238,5 @@ function build_polyh_bvh(
     end
 
     bvh = BVH(leaves, BBox{Float64})
-    return PolyhedralBVH(bvh, centres, radii, orientations, polyh_mesh)
+    return PolyhedralBVH(bvh, centres, radii, orientations, scales, local_verts, local_tris, local_normals, polyh_mesh)
 end
diff --git a/src/intersections.jl b/src/intersections.jl
index 9f25d51..4d14562 100644
--- a/src/intersections.jl
+++ b/src/intersections.jl
@@ -1,23 +1,23 @@
-# ============================================================
 # Ray/geometry intersections
 # ============================================================
-function ray_triangle_intersect(
-    p::SVector{3,Float64},
-    d::SVector{3,Float64},
-    v0::SVector{3,Float64},
-    v1::SVector{3,Float64},
-    v2::SVector{3,Float64};
-    eps=1e-12
+
+@inline function ray_triangle_intersect(
+    p::NTuple{3,Float64},
+    d::NTuple{3,Float64},
+    v0::NTuple{3,Float64},
+    v1::NTuple{3,Float64},
+    v2::NTuple{3,Float64};
+    eps::Float64=1e-12,
 )
-    e1 = v1 - v0
-    e2 = v2 - v0
+    e1 = sub3(v1, v0)
+    e2 = sub3(v2, v0)
     h = cross3(d, e2)
     a = dot3(e1, h)
     if abs(a) < eps
         return nothing
     end
     f = 1.0 / a
-    s = p - v0
+    s = sub3(p, v0)
     u = f * dot3(s, h)
     if (u < 0.0) || (u > 1.0)
         return nothing
@@ -31,50 +31,24 @@ function ray_triangle_intersect(
     return t > eps ? t : nothing
 end
 
-function nearest_surface_hit(p::SVector{3,Float64},
-    d::SVector{3,Float64},
-    surface::SurfaceBVH;
-    eps=1e-12)
-    P = to_mat3x1(p)
-    D = to_mat3x1(d)
-    trav = traverse_rays(surface.bvh, P, D)
-
-    tmin = Inf
-    imin = 0
-
-    @inbounds for (leaf_idx, ray_idx) in trav.contacts
-        ray_idx == 1 || continue
-        i1, i2, i3 = surface.tris[leaf_idx]
-        v0, v1, v2 = surface.verts[i1], surface.verts[i2], surface.verts[i3]
-        t = ray_triangle_intersect(p, d, v0, v1, v2; eps=eps)
-        if t !== nothing && t < tmin
-            tmin = t
-            imin = leaf_idx
-        end
-    end
-
-    if imin == 0
-        return nothing
-    end
-
-    i1, i2, i3 = surface.tris[imin]
-    v0, v1, v2 = surface.verts[i1], surface.verts[i2], surface.verts[i3]
-    n = cross3(v1 - v0, v2 - v0)
-    invn = 1.0 / sqrt(dot3(n, n))
-    n = n * invn
-    kind = surface.kinds[imin]
-    return (tmin, imin, n, kind)
+@inline function triangle_unit_normal(
+    v0::NTuple{3,Float64},
+    v1::NTuple{3,Float64},
+    v2::NTuple{3,Float64},
+)
+    n = cross3(sub3(v1, v0), sub3(v2, v0))
+    invn = inv(sqrt(dot3(n, n)))
+    return mul3(n, invn)
 end
 
 @inline function ray_sphere_exit_distance(
-    p::SVector{3,Float64},
-    d::SVector{3,Float64},
-    c::SVector{3,Float64},
+    p::NTuple{3,Float64},
+    d::NTuple{3,Float64},
+    c::NTuple{3,Float64},
     r::Float64;
-    eps=1e-12
+    eps::Float64=1e-12,
 )
-    # Assumes d is unit. Returns the forward distance to exit the sphere.
-    m = p - c
+    m = sub3(p, c)
     b = dot3(m, d)
     c0 = dot3(m, m) - r * r
     disc = b * b - c0
@@ -87,11 +61,11 @@ end
 end
 
 function ray_aabb_intersect(
-    p::SVector{3,Float64},
-    d::SVector{3,Float64},
-    mins::SVector{3,Float64},
-    maxs::SVector{3,Float64};
-    eps=1e-12
+    p::NTuple{3,Float64},
+    d::NTuple{3,Float64},
+    mins::NTuple{3,Float64},
+    maxs::NTuple{3,Float64};
+    eps::Float64=1e-12,
 )
     tmin = -Inf
     tmax = Inf
@@ -121,14 +95,13 @@ function ray_aabb_intersect(
 end
 
 @inline function ray_sphere_intersect(
-    p::SVector{3,Float64},
-    d::SVector{3,Float64},
-    c::SVector{3,Float64},
+    p::NTuple{3,Float64},
+    d::NTuple{3,Float64},
+    c::NTuple{3,Float64},
     r::Float64;
-    eps=1e-12
+    eps::Float64=1e-12,
 )
-    # Solve ||p + t d - c||^2 = r^2
-    m = p - c
+    m = sub3(p, c)
     b = dot3(m, d)
     c0 = dot3(m, m) - r * r
     disc = b * b - c0
@@ -144,126 +117,180 @@ end
     return t2 > eps ? t2 : nothing
 end
 
-function nearest_sphere_hit(
-    p::SVector{3,Float64},
-    d::SVector{3,Float64},
-    sbvh::SphereBVH;
-    eps=1e-12
-)
-    P = to_mat3x1(p)
-    D = to_mat3x1(d)
-    trav = traverse_rays(sbvh.bvh, P, D)
+@inline function rotate3(R::SMatrix{3,3,Float64,9}, v::NTuple{3,Float64})
+    return (
+        R[1, 1] * v[1] + R[1, 2] * v[2] + R[1, 3] * v[3],
+        R[2, 1] * v[1] + R[2, 2] * v[2] + R[2, 3] * v[3],
+        R[3, 1] * v[1] + R[3, 2] * v[2] + R[3, 3] * v[3],
+    )
+end
 
-    tmin = Inf
-    imin = 0
+@inline function rotate3_transpose(R::SMatrix{3,3,Float64,9}, v::NTuple{3,Float64})
+    return (
+        R[1, 1] * v[1] + R[2, 1] * v[2] + R[3, 1] * v[3],
+        R[1, 2] * v[1] + R[2, 2] * v[2] + R[3, 2] * v[3],
+        R[1, 3] * v[1] + R[2, 3] * v[2] + R[3, 3] * v[3],
+    )
+end
 
-    @inbounds for (leaf_idx, ray_idx) in trav.contacts
-        ray_idx == 1 || continue
-        c = sbvh.centres[leaf_idx]
-        r = sbvh.radii[leaf_idx]
-        t = ray_sphere_intersect(p, d, c, r; eps=eps)
-        if t !== nothing && t < tmin
-            tmin = t
-            imin = leaf_idx
+@inline function ray_polyh_intersect(
+    p::NTuple{3,Float64},
+    d::NTuple{3,Float64},
+    pbvh::PolyhedralBVH,
+    idx::Int;
+    eps::Float64=1e-12,
+)
+    c = pbvh.centres[idx]
+    c_tuple = (c[1], c[2], c[3])
+    R = pbvh.orientations[idx]
+    s = pbvh.scales[idx]
+
+    invs = 1.0 / s
+    p_local = mul3(rotate3_transpose(R, sub3(p, c_tuple)), invs)
+    d_local = mul3(rotate3_transpose(R, d), invs)
+
+    t_best = Inf
+    n_local_best = (0.0, 0.0, 0.0)
+    tris = pbvh.local_tris
+    verts = pbvh.local_verts
+    normals = pbvh.local_normals
+
+    @inbounds for j in eachindex(tris)
+        i1, i2, i3 = tris[j]
+        v0 = verts[i1]
+        v1 = verts[i2]
+        v2 = verts[i3]
+        t = ray_triangle_intersect(p_local, d_local, v0, v1, v2; eps=eps)
+        t === nothing && continue
+        if t < t_best
+            t_best = t
+            n_local_best = normals[j]
         end
     end
 
-    return imin == 0 ? nothing : (tmin, imin)
-end
+    isfinite(t_best) || return nothing
 
+    n_world = normalize3(rotate3(R, n_local_best))
+    return (t_best, n_world)
+end
 
-function ray_polyh_intersect(
-    p::SVector{3,Float64},
-    d::SVector{3,Float64},
-    center::SVector{3,Float64},
-    orient::SMatrix{3,3,Float64,9},
-    scale::Float64,
-    template_mesh::ParticleTriangleMesh;
-    eps=1e-12
+function update_nearest_sphere_hits!(
+    rb::RayBatchBuffer{Float64},
+    sbvh::SphereBVH,
+    traversal::ActiveTraversalCache;
+    eps::Float64=1e-12,
 )
+    P = rb.positions
+    D = rb.directions
+    active = rb.active
+    sphere_t = rb.sphere_t
+    sphere_idx = rb.sphere_idx
+
+    @inbounds for icontact in 1:traversal.num_contacts
+        leaf_idx, ray_idx_raw = traversal.contacts[icontact]
+        ray_idx = Int(ray_idx_raw)
+        active[ray_idx] || continue
+
+        p = (P[1, ray_idx], P[2, ray_idx], P[3, ray_idx])
+        d = (D[1, ray_idx], D[2, ray_idx], D[3, ray_idx])
+        c = sbvh.centres[Int(leaf_idx)]
+        r = sbvh.radii[Int(leaf_idx)]
+        t = ray_sphere_intersect(p, d, c, r; eps=eps)
+        t === nothing && continue
 
-    # Transform global ray to local object space
-    # p_local = (1/s) * R' * (p - c)
-    # d_local = R' * d
-    p_rel = p - center
-    p_local = (orient' * p_rel) * (1.0 / scale)
-    d_local = orient' * d
-
-    # Intersect with template mesh in local space
-    V = template_mesh.vertices
-    C = template_mesh.connectivity
-
-    tmin_local = Inf
-    hit_found = false
-    n_local_min = SVector{3,Float64}(0, 0, 0)
-
-    # Brute force check against all template triangles
-    for i in 1:size(C, 2)
-        idx1, idx2, idx3 = C[1, i], C[2, i], C[3, i]
-        v0 = SVector(V[1, idx1], V[2, idx1], V[3, idx1])
-        v1 = SVector(V[1, idx2], V[2, idx2], V[3, idx2])
-        v2 = SVector(V[1, idx3], V[2, idx3], V[3, idx3])
-
-        t_loc = ray_triangle_intersect(p_local, d_local, v0, v1, v2; eps=eps)
-        if t_loc !== nothing && t_loc < tmin_local
-            tmin_local = t_loc
-            e1 = v1 - v0
-            e2 = v2 - v0
-            n_local_min = cross3(e1, e2)
-            hit_found = true
+        tcur = sphere_t[ray_idx]
+        idxcur = sphere_idx[ray_idx]
+        if (t < tcur) || ((t == tcur) && (idxcur == 0 || Int(leaf_idx) < idxcur))
+            sphere_t[ray_idx] = t
+            sphere_idx[ray_idx] = Int(leaf_idx)
         end
     end
+    return nothing
+end
 
-    if hit_found
-        t_global = tmin_local * scale
-
-        orient_n = orient * n_local_min
-        invlen = 1.0 / sqrt(dot3(orient_n, orient_n))
-        n_global = orient_n * invlen
-
-        return t_global, n_global
-    else
-        return nothing
+function update_nearest_polyh_hits!(
+    rb::RayBatchBuffer{Float64},
+    pbvh::PolyhedralBVH,
+    traversal::ActiveTraversalCache;
+    eps::Float64=1e-12,
+)
+    P = rb.positions
+    D = rb.directions
+    active = rb.active
+    polyh_t = rb.polyh_t
+    polyh_idx = rb.polyh_idx
+    polyh_nx = rb.polyh_nx
+    polyh_ny = rb.polyh_ny
+    polyh_nz = rb.polyh_nz
+
+    @inbounds for icontact in 1:traversal.num_contacts
+        leaf_idx_raw, ray_idx_raw = traversal.contacts[icontact]
+        ray_idx = Int(ray_idx_raw)
+        active[ray_idx] || continue
+        leaf_idx = Int(leaf_idx_raw)
+
+        p = (P[1, ray_idx], P[2, ray_idx], P[3, ray_idx])
+        d = (D[1, ray_idx], D[2, ray_idx], D[3, ray_idx])
+        hit = ray_polyh_intersect(p, d, pbvh, leaf_idx; eps=eps)
+        hit === nothing && continue
+        t, n = hit
+
+        tcur = polyh_t[ray_idx]
+        idxcur = polyh_idx[ray_idx]
+        if (t < tcur) || ((t == tcur) && (idxcur == 0 || leaf_idx < idxcur))
+            polyh_t[ray_idx] = t
+            polyh_idx[ray_idx] = leaf_idx
+            polyh_nx[ray_idx] = n[1]
+            polyh_ny[ray_idx] = n[2]
+            polyh_nz[ray_idx] = n[3]
+        end
     end
+    return nothing
 end
 
-function nearest_polyh_hit(
-    p::SVector{3,Float64},
-    d::SVector{3,Float64},
-    pbvh::PolyhedralBVH;
-    eps=1e-12
+function update_nearest_surface_hits!(
+    rb::RayBatchBuffer{Float64},
+    surface::SurfaceBVH,
+    traversal::ActiveTraversalCache;
+    eps::Float64=1e-12,
 )
-    P = to_mat3x1(p)
-    D = to_mat3x1(d)
+    P = rb.positions
+    D = rb.directions
+    active = rb.active
+    wall_t = rb.wall_t
+    wall_idx = rb.wall_idx
+    wall_nx = rb.wall_nx
+    wall_ny = rb.wall_ny
+    wall_nz = rb.wall_nz
+    wall_sink = rb.wall_sink
+
+    @inbounds for icontact in 1:traversal.num_contacts
+        leaf_idx_raw, ray_idx_raw = traversal.contacts[icontact]
+        ray_idx = Int(ray_idx_raw)
+        active[ray_idx] || continue
+        leaf_idx = Int(leaf_idx_raw)
 
-    trav = traverse_rays(pbvh.bvh, P, D)
-
-    tmin = Inf
-    imin = 0
-    nmin = SVector{3,Float64}(0, 0, 0)
-
-    template = pbvh.polyh_mesh
-
-    @inbounds for (leaf_idx, ray_idx) in trav.contacts
-        ray_idx == 1 || continue
-
-        c = pbvh.centres[leaf_idx]
-        r = pbvh.radii[leaf_idx]
-        orient = pbvh.orientations[leaf_idx]
-
-        sf = get_scale_factor(template, Float64(r))
-        res = ray_polyh_intersect(p, d, c, orient, Float64(sf), template; eps=eps / sf)
-
-        if res !== nothing
-            t, n = res
-            if t < tmin && t > eps
-                tmin = t
-                imin = leaf_idx
-                nmin = n
-            end
+        i1, i2, i3 = surface.tris[leaf_idx]
+        v0 = surface.verts[i1]
+        v1 = surface.verts[i2]
+        v2 = surface.verts[i3]
+        p = (P[1, ray_idx], P[2, ray_idx], P[3, ray_idx])
+        d = (D[1, ray_idx], D[2, ray_idx], D[3, ray_idx])
+        t = ray_triangle_intersect(p, d, v0, v1, v2; eps=eps)
+        t === nothing && continue
+
+        tcur = wall_t[ray_idx]
+        idxcur = wall_idx[ray_idx]
+        if (t < tcur) || ((t == tcur) && (idxcur == 0 || leaf_idx < idxcur))
+            n = triangle_unit_normal(v0, v1, v2)
+            wall_t[ray_idx] = t
+            wall_idx[ray_idx] = leaf_idx
+            wall_nx[ray_idx] = n[1]
+            wall_ny[ray_idx] = n[2]
+            wall_nz[ray_idx] = n[3]
+            wall_sink[ray_idx] = surface.kinds[leaf_idx] isa Sink
         end
     end
 
-    return imin == 0 ? nothing : (tmin, imin, nmin)
+    return nothing
 end
-
diff --git a/src/ray_batch.jl b/src/ray_batch.jl
new file mode 100644
index 0000000..2058925
--- /dev/null
+++ b/src/ray_batch.jl
@@ -0,0 +1,174 @@
+# Batched ray state and tracing results
+# ============================================================
+
+@enum RayTermination::UInt8 begin
+    term_sink = 1
+    term_bbox_sink = 2
+    term_max_bounces = 3
+    term_max_length = 4
+    term_no_hit = 5
+end
+
+struct RayTraceBatchResult
+    escaped::BitVector
+    total_length::Vector{Float64}
+    nsteps::Vector{Int}
+    termination::Vector{RayTermination}
+    p0_termination_distance::Vector{Float64}
+end
+
+mutable struct ActiveTraversalCache{I <: Integer}
+    contacts::Vector{ImplicitBVH.IndexPair{I}}
+    thread_ncontacts_raw::Vector{I}
+    num_contacts::Int
+end
+
+mutable struct RayBatchBuffer{T <: AbstractFloat}
+    positions::Matrix{T}
+    directions::Matrix{T}
+    active::BitVector
+    escaped::BitVector
+    total_length::Vector{T}
+    nsteps::Vector{Int}
+    termination::Vector{RayTermination}
+    p0_termination_distance::Vector{T}
+
+    sphere_t::Vector{T}
+    sphere_idx::Vector{Int}
+    polyh_t::Vector{T}
+    polyh_idx::Vector{Int}
+    polyh_nx::Vector{T}
+    polyh_ny::Vector{T}
+    polyh_nz::Vector{T}
+
+    wall_t::Vector{T}
+    wall_idx::Vector{Int}
+    wall_nx::Vector{T}
+    wall_ny::Vector{T}
+    wall_nz::Vector{T}
+    wall_sink::BitVector
+end
+
+function RayBatchBuffer(::Type{T}, nrays::Integer) where {T <: AbstractFloat}
+    nrays > 0 || throw(ArgumentError("nrays must be > 0"))
+
+    positions = Matrix{T}(undef, 3, nrays)
+    directions = Matrix{T}(undef, 3, nrays)
+    active = trues(nrays)
+    escaped = falses(nrays)
+    total_length = zeros(T, nrays)
+    nsteps = zeros(Int, nrays)
+    termination = fill(term_no_hit, nrays)
+    p0_termination_distance = zeros(T, nrays)
+
+    sphere_t = fill(T(Inf), nrays)
+    sphere_idx = zeros(Int, nrays)
+    polyh_t = fill(T(Inf), nrays)
+    polyh_idx = zeros(Int, nrays)
+    polyh_nx = zeros(T, nrays)
+    polyh_ny = zeros(T, nrays)
+    polyh_nz = zeros(T, nrays)
+
+    wall_t = fill(T(Inf), nrays)
+    wall_idx = zeros(Int, nrays)
+    wall_nx = zeros(T, nrays)
+    wall_ny = zeros(T, nrays)
+    wall_nz = zeros(T, nrays)
+    wall_sink = falses(nrays)
+
+    RayBatchBuffer(
+        positions, directions, active, escaped, total_length, nsteps, termination, p0_termination_distance,
+        sphere_t, sphere_idx, polyh_t, polyh_idx, polyh_nx, polyh_ny, polyh_nz,
+        wall_t, wall_idx, wall_nx, wall_ny, wall_nz, wall_sink,
+    )
+end
+
+function reset_ray_batch!(rb::RayBatchBuffer{T}, p0::NTuple{3,<:Real}, D0::AbstractMatrix{<:Real}) where {T}
+    size(D0, 1) == 3 || throw(ArgumentError("directions must be a (3, N) matrix"))
+    size(D0, 2) == size(rb.positions, 2) || throw(ArgumentError("direction count does not match buffer size"))
+
+    n = size(D0, 2)
+    px = T(p0[1]); py = T(p0[2]); pz = T(p0[3])
+
+    @inbounds for i in 1:n
+        rb.positions[1, i] = px
+        rb.positions[2, i] = py
+        rb.positions[3, i] = pz
+
+        dx = T(D0[1, i])
+        dy = T(D0[2, i])
+        dz = T(D0[3, i])
+        invn = inv(sqrt(dx * dx + dy * dy + dz * dz))
+        if !isfinite(invn)
+            throw(ArgumentError("direction column $i has near-zero norm"))
+        end
+        rb.directions[1, i] = dx * invn
+        rb.directions[2, i] = dy * invn
+        rb.directions[3, i] = dz * invn
+    end
+
+    fill!(rb.active, true)
+    fill!(rb.escaped, false)
+    fill!(rb.total_length, zero(T))
+    fill!(rb.nsteps, 0)
+    fill!(rb.termination, term_no_hit)
+    fill!(rb.p0_termination_distance, zero(T))
+
+    fill!(rb.sphere_t, T(Inf))
+    fill!(rb.sphere_idx, 0)
+    fill!(rb.polyh_t, T(Inf))
+    fill!(rb.polyh_idx, 0)
+    fill!(rb.polyh_nx, zero(T))
+    fill!(rb.polyh_ny, zero(T))
+    fill!(rb.polyh_nz, zero(T))
+    fill!(rb.wall_t, T(Inf))
+    fill!(rb.wall_idx, 0)
+    fill!(rb.wall_nx, zero(T))
+    fill!(rb.wall_ny, zero(T))
+    fill!(rb.wall_nz, zero(T))
+    fill!(rb.wall_sink, false)
+
+    return rb
+end
+
+function reset_bounce_hits!(rb::RayBatchBuffer{T}) where {T}
+    fill!(rb.sphere_t, T(Inf))
+    fill!(rb.sphere_idx, 0)
+    fill!(rb.polyh_t, T(Inf))
+    fill!(rb.polyh_idx, 0)
+    fill!(rb.polyh_nx, zero(T))
+    fill!(rb.polyh_ny, zero(T))
+    fill!(rb.polyh_nz, zero(T))
+    fill!(rb.wall_t, T(Inf))
+    fill!(rb.wall_idx, 0)
+    fill!(rb.wall_nx, zero(T))
+    fill!(rb.wall_ny, zero(T))
+    fill!(rb.wall_nz, zero(T))
+    fill!(rb.wall_sink, false)
+    return rb
+end
+
+@inline function finish_ray!(
+    rb::RayBatchBuffer,
+    i::Int,
+    reason::RayTermination,
+    escaped::Bool,
+)
+    rb.active[i] = false
+    rb.termination[i] = reason
+    rb.escaped[i] = escaped
+    return nothing
+end
+
+@inline function set_termination_distance!(
+    rb::RayBatchBuffer{T},
+    i::Int,
+    p0::NTuple{3,Float64},
+    p_term::NTuple{3,Float64},
+) where {T}
+    dx = p_term[1] - p0[1]
+    dy = p_term[2] - p0[2]
+    dz = p_term[3] - p0[3]
+    rb.p0_termination_distance[i] = T(sqrt(dx * dx + dy * dy + dz * dz))
+    return nothing
+end
diff --git a/src/raytrace_active_lvt.jl b/src/raytrace_active_lvt.jl
new file mode 100644
index 0000000..62043e1
--- /dev/null
+++ b/src/raytrace_active_lvt.jl
@@ -0,0 +1,126 @@
+# Active-ray LVT traversal wrapper (CPU only, in-place cache)
+# ============================================================
+
+function ensure_active_traversal_cache(
+    cache::Any,
+    bvh::BVH,
+    options,
+)
+    I = ImplicitBVH.get_index_type(bvh)
+    npadding = div(64, sizeof(I), RoundUp)
+    thread_ncontacts_length = npadding * options.num_threads
+
+    if cache isa ActiveTraversalCache{I}
+        if length(cache.thread_ncontacts_raw) < thread_ncontacts_length
+            resize!(cache.thread_ncontacts_raw, thread_ncontacts_length)
+        end
+        return cache, npadding
+    end
+
+    contacts = similar(bvh.nodes, ImplicitBVH.IndexPair{I}, 0)
+    thread_ncontacts_raw = similar(bvh.nodes, I, thread_ncontacts_length)
+    return ActiveTraversalCache{I}(contacts, thread_ncontacts_raw, 0), npadding
+end
+
+function traverse_rays_active_lvt!(
+    bvh::BVH,
+    points::AbstractMatrix,
+    directions::AbstractMatrix,
+    active::AbstractVector{Bool};
+    start_level::Int=ImplicitBVH.default_start_level(bvh, ImplicitBVH.LVTTraversal()),
+    narrow=(bv, p, d) -> true,
+    cache=nothing,
+    options=ImplicitBVH.BVHOptions(),
+)
+    size(points, 1) == 3 || throw(ArgumentError("points must have shape (3, N)"))
+    size(directions, 1) == 3 || throw(ArgumentError("directions must have shape (3, N)"))
+    size(points, 2) == size(directions, 2) || throw(ArgumentError("point and direction counts must match"))
+    length(active) == size(points, 2) || throw(ArgumentError("active mask length must match ray count"))
+    (bvh.built_level <= start_level <= bvh.tree.levels <= 32) || throw(ArgumentError("invalid start level"))
+
+    AK.get_backend(bvh.nodes) == AK.get_backend(points) == AK.get_backend(directions) ||
+        throw(ArgumentError("bvh, points and directions must share backend"))
+    AK.get_backend(bvh.nodes) isa typeof(AK.CPU_BACKEND) ||
+        throw(ArgumentError("traverse_rays_active_lvt! currently supports CPU backend only"))
+
+    cache_typed, npadding = ensure_active_traversal_cache(cache, bvh, options)
+
+    num_rays = size(points, 2)
+    if num_rays == 0
+        cache_typed.num_contacts = 0
+        return cache_typed
+    end
+
+    I = ImplicitBVH.get_index_type(bvh)
+    thread_ncontacts = view(reshape(cache_typed.thread_ncontacts_raw, npadding, :), 1, :)
+    thread_ncontacts .= I(0)
+
+    traverse_rays_lvt_active!(
+        bvh, points, directions, active, start_level,
+        thread_ncontacts, nothing, narrow, options,
+    )
+
+    AK.accumulate!(+, thread_ncontacts, init=I(0), max_tasks=1, block_size=options.block_size)
+
+    total_contacts = Int(thread_ncontacts[end])
+    cache_typed.num_contacts = total_contacts
+    if length(cache_typed.contacts) < total_contacts
+        resize!(cache_typed.contacts, total_contacts)
+    end
+
+    if total_contacts > 0
+        traverse_rays_lvt_active!(
+            bvh, points, directions, active, start_level,
+            thread_ncontacts, cache_typed.contacts, narrow, options,
+        )
+    end
+
+    return cache_typed
+end
+
+function traverse_rays_lvt_active!(
+    bvh::BVH,
+    points::AbstractMatrix,
+    directions::AbstractMatrix,
+    active::AbstractVector{Bool},
+    start_level::Int,
+    thread_ncontacts::AbstractVector,
+    contacts::Union{Nothing, AbstractVector},
+    narrow,
+    options,
+)
+    num_checks = size(points, 2)
+    AK.itask_partition(num_checks, options.num_threads, options.min_traversals_per_thread) do itask, irange
+        I = ImplicitBVH.get_index_type(bvh)
+        T = eltype(eltype(bvh.leaves))
+        stack = ImplicitBVH.SimpleMVector{32, I}(undef)
+        iwrite = if isnothing(contacts)
+            Ref(I(itask))
+        else
+            Ref(itask == 0x1 ? I(0x1) : I(thread_ncontacts[itask - 0x1] + 0x1))
+        end
+
+        for i in irange
+            @inbounds active[i] || continue
+            point = (
+                T(points[1, i]),
+                T(points[2, i]),
+                T(points[3, i]),
+            )
+            direction = (
+                T(directions[1, i]),
+                T(directions[2, i]),
+                T(directions[3, i]),
+            )
+            @inline ImplicitBVH.traverse_ray_lvt!(
+                point, direction, I(i),
+                bvh, stack,
+                I(start_level), iwrite,
+                thread_ncontacts, contacts,
+                narrow,
+            )
+        end
+    end
+
+    return nothing
+end
diff --git a/src/tracer.jl b/src/tracer.jl
index 57fe53f..3642c36 100644
--- a/src/tracer.jl
+++ b/src/tracer.jl
@@ -1,271 +1,474 @@
 # ============================================================
 # Tracer
 # ============================================================
-@inline function reflect_specular(d::SVector{3,Float64}, n::SVector{3,Float64})
-    return d - 2.0 * dot3(d, n) * n
-end
+
+@inline reflect_specular(d::NTuple{3,Float64}, n::NTuple{3,Float64}) = sub3(d, mul3(n, 2.0 * dot3(d, n)))
 
 const _bbox_warned = Ref(false)
+const _polyh_start_inside_warned = Ref(false)
 
+# Legacy visualisation container kept for backwards compatibility with visualise.jl.
 struct RayTraceResult
     escaped::Bool
     total_length::Float64
     nsteps::Int
-    path::Vector{SVector{3,Float64}}   # polyline vertices (start, hits..., exit)
+    path::Vector{SVector{3,Float64}}
 end
 
-function start_outside_spheres_aabb(p::SVector{3,Float64},
-                                   d::SVector{3,Float64},
-                                   spheres::SphereBVH;
-                                   tol::Float64 = 1e-12,
-                                   eps_hit::Float64 = 1e-12)
-
-    px, py, pz = p
-    tmax = 0.0
-    n_inside = 0
-
-    @inbounds for i in eachindex(spheres.radii)
-        c = spheres.centres[i]
-        r = spheres.radii[i]
-
-        dx = px - c[1];  adx = abs(dx)
-        adx > r && continue
-        dy = py - c[2];  ady = abs(dy)
-        ady > r && continue
-        dz = pz - c[3];  adz = abs(dz)
-        adz > r && continue
-
-        # Only now do the actual sphere check
-        r_eff = r - tol
-        r_eff <= 0 && continue
-        if (dx*dx + dy*dy + dz*dz) < (r_eff*r_eff)
-            n_inside += 1
-            t = ray_sphere_exit_distance(p, d, c, r; eps=eps_hit)
-            if t !== nothing && t > tmax
-                tmax = t
+function start_outside_spheres_aabb!(
+    rb::RayBatchBuffer{Float64},
+    spheres::SphereBVH;
+    tol::Float64=1e-12,
+    eps_hit::Float64=1e-12,
+    eps_shift::Float64=1e-10,
+    count_start_solid::Bool=false,
+)
+    P = rb.positions
+    D = rb.directions
+    T = rb.total_length
+
+    @inbounds for i in axes(P, 2)
+        p = (P[1, i], P[2, i], P[3, i])
+        d = (D[1, i], D[2, i], D[3, i])
+        tmax = 0.0
+        n_inside = 0
+
+        for j in eachindex(spheres.radii)
+            c = spheres.centres[j]
+            r = spheres.radii[j]
+
+            dx = p[1] - c[1]; adx = abs(dx)
+            adx > r && continue
+            dy = p[2] - c[2]; ady = abs(dy)
+            ady > r && continue
+            dz = p[3] - c[3]; adz = abs(dz)
+            adz > r && continue
+
+            r_eff = r - tol
+            r_eff <= 0 && continue
+            if (dx * dx + dy * dy + dz * dz) < (r_eff * r_eff)
+                n_inside += 1
+                t = ray_sphere_exit_distance(p, d, c, r; eps=eps_hit)
+                if t !== nothing && t > tmax
+                    tmax = t
+                end
             end
         end
-    end
-
-    return n_inside, tmax
-end
-
-function trace_ray_geometric(p0::SVector{3,Float64},
-                             d0::SVector{3,Float64},
-                             surfaces::AbstractVector{<:TriangleSurface},
-                             spheres::SphereBVH;
-                             eps::Float64=1e-10,
-                             tri_eps::Float64=1e-12,
-                             sph_eps::Float64=1e-12,
-                             max_steps::Int=1_000_000,
-                             max_length::Float64=Inf,
-                             record_path::Bool=true,
-                             allow_start_inside::Bool=true,
-                             count_start_solid::Bool=false)
-
-    surface_bvh = build_surface_bvh(surfaces)
-    ray_escaped = false
-
-    invn = 1.0 / sqrt(dot3(d0, d0))
-    d = d0 * invn
-    p = p0
 
-    path = record_path ? SVector{3,Float64}[p0] : SVector{3,Float64}[]
-    total = 0.0
-
-    # --- NEW: if starting inside spheres, move to outside all of them (no reflection) ---
-    if allow_start_inside
-        n_inside, tmax = start_outside_spheres_aabb(p, d, spheres; tol=1e-12, eps_hit=sph_eps)
         if n_inside > 0
-            p = p + (tmax + eps) * d
+            p_next = madd3(p, tmax + eps_shift, d)
+            P[1, i] = p_next[1]
+            P[2, i] = p_next[2]
+            P[3, i] = p_next[3]
             if count_start_solid
-                total += tmax
+                T[i] += tmax
             end
-            record_path && push!(path, p)
         end
     end
 
-    for step in 1:max_steps
-        sh = nearest_sphere_hit(p, d, spheres; eps=sph_eps)
-        mh = nearest_surface_hit(p, d, surface_bvh; eps=tri_eps)
-        tb = ray_aabb_intersect(p, d, surface_bvh.bbox_mins, surface_bvh.bbox_maxs; eps=eps)
-
-        ts = sh === nothing ? Inf : sh[1]
-        tm = mh === nothing ? Inf : mh[1]
-        tbox = tb === nothing ? Inf : tb
-
-        if isinf(ts) && isinf(tm) && isinf(tbox)
-            # If the vessel is truly closed and you're inside it, this shouldn't happen.
-            return RayTraceResult(false, total, step, path)
-        end
-
-        t = min(ts, tm, tbox)
-        p_hit = p + t * d
-        total += t
+    return nothing
+end
 
-        if record_path
-            push!(path, p_hit)
-        end
+function trace_rays(
+    p0::NTuple{3,<:Real},
+    directions::AbstractMatrix{<:Real},
+    surfaces::AbstractVector{<:TriangleSurface},
+    spheres::SphereBVH;
+    eps::Float64=1e-10,
+    tri_eps::Float64=1e-12,
+    sph_eps::Float64=1e-12,
+    max_bounces::Int=1_000_000,
+    max_length::Float64=Inf,
+    allow_start_inside::Bool=true,
+    count_start_solid::Bool=false,
+    options=ImplicitBVH.BVHOptions(),
+)
+    surface_bvh = build_surface_bvh(surfaces)
+    return trace_rays(
+        p0, directions, surface_bvh, spheres;
+        eps=eps,
+        tri_eps=tri_eps,
+        sph_eps=sph_eps,
+        max_bounces=max_bounces,
+        max_length=max_length,
+        allow_start_inside=allow_start_inside,
+        count_start_solid=count_start_solid,
+        options=options,
+    )
+end
 
-        if total > max_length
-            return RayTraceResult(false, total, step, path)
+function trace_rays(
+    p0::NTuple{3,<:Real},
+    directions::AbstractMatrix{<:Real},
+    surface_bvh::SurfaceBVH,
+    spheres::SphereBVH;
+    eps::Float64=1e-10,
+    tri_eps::Float64=1e-12,
+    sph_eps::Float64=1e-12,
+    max_bounces::Int=1_000_000,
+    max_length::Float64=Inf,
+    allow_start_inside::Bool=true,
+    count_start_solid::Bool=false,
+    options=ImplicitBVH.BVHOptions(),
+)
+    max_bounces >= 0 || throw(ArgumentError("max_bounces must be >= 0"))
+    size(directions, 1) == 3 || throw(ArgumentError("directions must be a (3, N) matrix"))
+    n_rays = size(directions, 2)
+    n_rays > 0 || throw(ArgumentError("directions must contain at least one ray"))
+
+    p0f = (Float64(p0[1]), Float64(p0[2]), Float64(p0[3]))
+    rb = RayBatchBuffer(Float64, n_rays)
+    reset_ray_batch!(rb, p0f, directions)
+
+    if max_bounces == 0
+        @inbounds for i in eachindex(rb.active)
+            set_termination_distance!(rb, i, p0f, p0f)
+            finish_ray!(rb, i, PRICK.term_max_bounces, false)
         end
+        return RayTraceBatchResult(
+            copy(rb.escaped),
+            copy(rb.total_length),
+            copy(rb.nsteps),
+            copy(rb.termination),
+            copy(rb.p0_termination_distance),
+        )
+    end
 
-        if tbox <= tm && tbox <= ts
-            # Hit bounding box sink
-            ray_escaped = true
-            if !_bbox_warned[]
-                @warn "Ray hit bounding box sink; mesh may not be watertight."
-                _bbox_warned[] = true
-            end
-            return RayTraceResult(true, total, step, path)
-        elseif tm < ts
-            # Hit triangle surface
-            _, _, nhat, kind = mh
-            if dot3(d, nhat) > 0.0
-                nhat = -nhat
-            end
-
-            if kind isa Sink
-                return RayTraceResult(true, total, step, path)
+    allow_start_inside && start_outside_spheres_aabb!(
+        rb, spheres; tol=1e-12, eps_hit=sph_eps, eps_shift=eps, count_start_solid=count_start_solid
+    )
+
+    sphere_traversal_cache = nothing
+    wall_traversal_cache = nothing
+
+    active_count = n_rays
+    while active_count > 0
+        reset_bounce_hits!(rb)
+
+        sphere_traversal = traverse_rays_active_lvt!(
+            spheres.bvh, rb.positions, rb.directions, rb.active;
+            cache=sphere_traversal_cache,
+            options=options,
+        )
+        sphere_traversal_cache = sphere_traversal
+        update_nearest_sphere_hits!(rb, spheres, sphere_traversal; eps=sph_eps)
+
+        wall_traversal = traverse_rays_active_lvt!(
+            surface_bvh.bvh, rb.positions, rb.directions, rb.active;
+            cache=wall_traversal_cache,
+            options=options,
+        )
+        wall_traversal_cache = wall_traversal
+        update_nearest_surface_hits!(rb, surface_bvh, wall_traversal; eps=tri_eps)
+
+        @inbounds for i in 1:n_rays
+            rb.active[i] || continue
+            p = (rb.positions[1, i], rb.positions[2, i], rb.positions[3, i])
+            d = (rb.directions[1, i], rb.directions[2, i], rb.directions[3, i])
+
+            ts = rb.sphere_t[i]
+            tm = rb.wall_t[i]
+            tbox = ray_aabb_intersect(p, d, surface_bvh.bbox_mins, surface_bvh.bbox_maxs; eps=eps)
+            tboxf = tbox === nothing ? Inf : tbox
+
+            if isinf(ts) && isinf(tm) && isinf(tboxf)
+                set_termination_distance!(rb, i, p0f, p)
+                finish_ray!(rb, i, term_no_hit, false)
+                active_count -= 1
+                continue
             end
 
-            d = reflect_specular(d, nhat)
-            d *= 1.0 / sqrt(dot3(d, d))
-            p = p_hit + eps * nhat
-        else
-            # Hit sphere => reflect and continue (robust)
-            _, idx = sh
-            c = spheres.centres[idx]
-            r = spheres.radii[idx]
+            t = min(ts, tm, tboxf)
+            p_hit = madd3(p, t, d)
+            rb.total_length[i] += t
 
-            # Robust unit normal from centre->hit
-            nhat = p_hit - c
-            invlen = 1.0 / sqrt(dot3(nhat, nhat))
-            nhat = nhat * invlen
+            if rb.total_length[i] > max_length
+                set_termination_distance!(rb, i, p0f, p_hit)
+                finish_ray!(rb, i, PRICK.term_max_length, false)
+                active_count -= 1
+                continue
+            end
 
-            # Specular reflection about unit normal
-            d = reflect_specular(d, nhat)
+            if tboxf <= tm && tboxf <= ts
+                if !_bbox_warned[]
+                    @warn "Ray hit bounding box sink; mesh may not be watertight."
+                    _bbox_warned[] = true
+                end
+                set_termination_distance!(rb, i, p0f, p_hit)
+                finish_ray!(rb, i, term_bbox_sink, true)
+                active_count -= 1
+                continue
+            elseif tm < ts
+                nhat = (rb.wall_nx[i], rb.wall_ny[i], rb.wall_nz[i])
+                if dot3(d, nhat) > 0.0
+                    nhat = mul3(nhat, -1.0)
+                end
 
-            # Renormalise direction to avoid drift
-            d *= 1.0 / sqrt(dot3(d, d))
+                if rb.wall_sink[i]
+                    set_termination_distance!(rb, i, p0f, p_hit)
+                    finish_ray!(rb, i, term_sink, true)
+                    active_count -= 1
+                    continue
+                end
 
-            # Reposition point to be just outside the sphere
-            eps_rel = max(eps, 1e-9 * r)
-            p = c + (r + eps_rel) * nhat
+                d_next = normalize3(reflect_specular(d, nhat))
+                p_next = madd3(p_hit, eps, nhat)
+                rb.directions[1, i] = d_next[1]
+                rb.directions[2, i] = d_next[2]
+                rb.directions[3, i] = d_next[3]
+                rb.positions[1, i] = p_next[1]
+                rb.positions[2, i] = p_next[2]
+                rb.positions[3, i] = p_next[3]
+
+                rb.nsteps[i] += 1
+                if rb.nsteps[i] >= max_bounces
+                    set_termination_distance!(rb, i, p0f, p_hit)
+                    finish_ray!(rb, i, PRICK.term_max_bounces, false)
+                    active_count -= 1
+                end
+            else
+                idx = rb.sphere_idx[i]
+                idx > 0 || throw(ArgumentError("internal error: missing sphere index for finite sphere hit"))
+
+                c = spheres.centres[idx]
+                r = spheres.radii[idx]
+                nhat = normalize3(sub3(p_hit, c))
+                d_next = normalize3(reflect_specular(d, nhat))
+                eps_rel = max(eps, 1e-9 * r)
+                p_next = add3(c, mul3(nhat, r + eps_rel))
+                rb.directions[1, i] = d_next[1]
+                rb.directions[2, i] = d_next[2]
+                rb.directions[3, i] = d_next[3]
+                rb.positions[1, i] = p_next[1]
+                rb.positions[2, i] = p_next[2]
+                rb.positions[3, i] = p_next[3]
+
+                rb.nsteps[i] += 1
+                if rb.nsteps[i] >= max_bounces
+                    set_termination_distance!(rb, i, p0f, p_hit)
+                    finish_ray!(rb, i, PRICK.term_max_bounces, false)
+                    active_count -= 1
+                end
+            end
         end
-
     end
 
-    return RayTraceResult(false, total, max_steps, path)
+    return RayTraceBatchResult(
+        copy(rb.escaped),
+        copy(rb.total_length),
+        copy(rb.nsteps),
+        copy(rb.termination),
+        copy(rb.p0_termination_distance),
+    )
 end
 
-
-function trace_ray_geometric(p0::SVector{3,Float64},
-    d0::SVector{3,Float64},
+function trace_rays(
+    p0::NTuple{3,<:Real},
+    directions::AbstractMatrix{<:Real},
     surfaces::AbstractVector{<:TriangleSurface},
-    particles::PolyhedralBVH;
+    polyh::PolyhedralBVH;
     eps::Float64=1e-10,
     tri_eps::Float64=1e-12,
-    part_eps::Float64=1e-12,
-    max_steps::Int=1_000_000,
+    polyh_eps::Float64=1e-12,
+    max_bounces::Int=1_000_000,
     max_length::Float64=Inf,
-    record_path::Bool=true,
     allow_start_inside::Bool=true,
-    count_start_solid::Bool=false)
-
+    count_start_solid::Bool=false,
+    options=ImplicitBVH.BVHOptions(),
+)
     surface_bvh = build_surface_bvh(surfaces)
+    return trace_rays(
+        p0, directions, surface_bvh, polyh;
+        eps=eps,
+        tri_eps=tri_eps,
+        polyh_eps=polyh_eps,
+        max_bounces=max_bounces,
+        max_length=max_length,
+        allow_start_inside=allow_start_inside,
+        count_start_solid=count_start_solid,
+        options=options,
+    )
+end
 
-    invn = 1.0 / sqrt(dot3(d0, d0))
-    d = d0 * invn
-    p = p0
-
-    path = record_path ? SVector{3,Float64}[p0] : SVector{3,Float64}[]
-    total = 0.0
-
-    # Note: start_outside_spheres logic is tricky for polyhedra. 
-    # For now we can assume ray is inside.
-    _ = allow_start_inside
-    _ = count_start_solid
-
-    for step in 1:max_steps
-        ph = nearest_polyh_hit(p, d, particles; eps=part_eps) # (t, idx, n)
-        mh = nearest_surface_hit(p, d, surface_bvh; eps=tri_eps)
-        tb = ray_aabb_intersect(p, d, surface_bvh.bbox_mins, surface_bvh.bbox_maxs; eps=eps)
-
-        ts = ph === nothing ? Inf : ph[1]
-        tm = mh === nothing ? Inf : mh[1]
-        tbox = tb === nothing ? Inf : tb
-
-        if isinf(ts) && isinf(tm) && isinf(tbox)
-            return RayTraceResult(false, total, step, path)
-        end
-
-        t = min(ts, tm, tbox)
-        p_hit = p + t * d
-        total += t
-
-        if record_path
-            push!(path, p_hit)
-        end
-
-        if total > max_length
-            return RayTraceResult(false, total, step, path)
+function trace_rays(
+    p0::NTuple{3,<:Real},
+    directions::AbstractMatrix{<:Real},
+    surface_bvh::SurfaceBVH,
+    polyh::PolyhedralBVH;
+    eps::Float64=1e-10,
+    tri_eps::Float64=1e-12,
+    polyh_eps::Float64=1e-12,
+    max_bounces::Int=1_000_000,
+    max_length::Float64=Inf,
+    allow_start_inside::Bool=true,
+    count_start_solid::Bool=false,
+    options=ImplicitBVH.BVHOptions(),
+)
+    max_bounces >= 0 || throw(ArgumentError("max_bounces must be >= 0"))
+    size(directions, 1) == 3 || throw(ArgumentError("directions must be a (3, N) matrix"))
+    n_rays = size(directions, 2)
+    n_rays > 0 || throw(ArgumentError("directions must contain at least one ray"))
+
+    p0f = (Float64(p0[1]), Float64(p0[2]), Float64(p0[3]))
+    rb = RayBatchBuffer(Float64, n_rays)
+    reset_ray_batch!(rb, p0f, directions)
+
+    if max_bounces == 0
+        @inbounds for i in eachindex(rb.active)
+            set_termination_distance!(rb, i, p0f, p0f)
+            finish_ray!(rb, i, PRICK.term_max_bounces, false)
         end
+        return RayTraceBatchResult(
+            copy(rb.escaped),
+            copy(rb.total_length),
+            copy(rb.nsteps),
+            copy(rb.termination),
+            copy(rb.p0_termination_distance),
+        )
+    end
 
-        if tbox <= tm && tbox <= ts
-            # Hit bounding box sink
-            return RayTraceResult(true, total, step, path)
-        elseif tm < ts
-            # Hit triangle surface (wall/mirror)
-            _, _, nhat, kind = mh
-            if dot3(d, nhat) > 0.0
-                nhat = -nhat
+    if allow_start_inside && !_polyh_start_inside_warned[]
+        @warn "allow_start_inside for polyhedral tracing is currently not implemented; start rays in void space."
+        _polyh_start_inside_warned[] = true
+    end
+    count_start_solid && @warn "count_start_solid has no effect for polyhedral tracing at present."
+
+    polyh_traversal_cache = nothing
+    wall_traversal_cache = nothing
+
+    active_count = n_rays
+    while active_count > 0
+        reset_bounce_hits!(rb)
+
+        polyh_traversal = traverse_rays_active_lvt!(
+            polyh.bvh, rb.positions, rb.directions, rb.active;
+            cache=polyh_traversal_cache,
+            options=options,
+        )
+        polyh_traversal_cache = polyh_traversal
+        update_nearest_polyh_hits!(rb, polyh, polyh_traversal; eps=polyh_eps)
+
+        wall_traversal = traverse_rays_active_lvt!(
+            surface_bvh.bvh, rb.positions, rb.directions, rb.active;
+            cache=wall_traversal_cache,
+            options=options,
+        )
+        wall_traversal_cache = wall_traversal
+        update_nearest_surface_hits!(rb, surface_bvh, wall_traversal; eps=tri_eps)
+
+        @inbounds for i in 1:n_rays
+            rb.active[i] || continue
+            p = (rb.positions[1, i], rb.positions[2, i], rb.positions[3, i])
+            d = (rb.directions[1, i], rb.directions[2, i], rb.directions[3, i])
+
+            tp = rb.polyh_t[i]
+            tm = rb.wall_t[i]
+            tbox = ray_aabb_intersect(p, d, surface_bvh.bbox_mins, surface_bvh.bbox_maxs; eps=eps)
+            tboxf = tbox === nothing ? Inf : tbox
+
+            if isinf(tp) && isinf(tm) && isinf(tboxf)
+                set_termination_distance!(rb, i, p0f, p)
+                finish_ray!(rb, i, term_no_hit, false)
+                active_count -= 1
+                continue
             end
 
-            if kind isa Sink
-                return RayTraceResult(true, total, step, path)
+            t = min(tp, tm, tboxf)
+            p_hit = madd3(p, t, d)
+            rb.total_length[i] += t
+
+            if rb.total_length[i] > max_length
+                set_termination_distance!(rb, i, p0f, p_hit)
+                finish_ray!(rb, i, PRICK.term_max_length, false)
+                active_count -= 1
+                continue
             end
 
-            d = reflect_specular(d, nhat)
-            # Normalize manually
-            d *= 1.0 / sqrt(dot3(d, d))
-            p = p_hit + eps * nhat
-        else
-            # Hit Polyhedron => reflect
-            _, idx, nhat = ph
+            if tboxf <= tm && tboxf <= tp
+                if !_bbox_warned[]
+                    @warn "Ray hit bounding box sink; mesh may not be watertight."
+                    _bbox_warned[] = true
+                end
+                set_termination_distance!(rb, i, p0f, p_hit)
+                finish_ray!(rb, i, term_bbox_sink, true)
+                active_count -= 1
+                continue
+            elseif tm < tp
+                nhat = (rb.wall_nx[i], rb.wall_ny[i], rb.wall_nz[i])
+                if dot3(d, nhat) > 0.0
+                    nhat = mul3(nhat, -1.0)
+                end
 
-            # Ensure normal faces against ray
-            if dot3(d, nhat) > 0.0
-                nhat = -nhat
-            end
+                if rb.wall_sink[i]
+                    set_termination_distance!(rb, i, p0f, p_hit)
+                    finish_ray!(rb, i, term_sink, true)
+                    active_count -= 1
+                    continue
+                end
 
-            # Specular reflection
-            d = reflect_specular(d, nhat)
-            # Normalize manually
-            d *= 1.0 / sqrt(dot3(d, d))
+                d_next = normalize3(reflect_specular(d, nhat))
+                p_next = madd3(p_hit, eps, nhat)
+                rb.directions[1, i] = d_next[1]
+                rb.directions[2, i] = d_next[2]
+                rb.directions[3, i] = d_next[3]
+                rb.positions[1, i] = p_next[1]
+                rb.positions[2, i] = p_next[2]
+                rb.positions[3, i] = p_next[3]
+
+                rb.nsteps[i] += 1
+                if rb.nsteps[i] >= max_bounces
+                    set_termination_distance!(rb, i, p0f, p_hit)
+                    finish_ray!(rb, i, PRICK.term_max_bounces, false)
+                    active_count -= 1
+                end
+            else
+                idx = rb.polyh_idx[i]
+                idx > 0 || throw(ArgumentError("internal error: missing polyhedral index for finite polyhedral hit"))
 
-            # Reposition point just outside
-            p = p_hit + eps * nhat
+                nhat = (rb.polyh_nx[i], rb.polyh_ny[i], rb.polyh_nz[i])
+                if dot3(d, nhat) > 0.0
+                    nhat = mul3(nhat, -1.0)
+                end
+                d_next = normalize3(reflect_specular(d, nhat))
+                eps_rel = max(eps, 1e-9 * polyh.radii[idx])
+                p_next = madd3(p_hit, eps_rel, nhat)
+                rb.directions[1, i] = d_next[1]
+                rb.directions[2, i] = d_next[2]
+                rb.directions[3, i] = d_next[3]
+                rb.positions[1, i] = p_next[1]
+                rb.positions[2, i] = p_next[2]
+                rb.positions[3, i] = p_next[3]
+
+                rb.nsteps[i] += 1
+                if rb.nsteps[i] >= max_bounces
+                    set_termination_distance!(rb, i, p0f, p_hit)
+                    finish_ray!(rb, i, PRICK.term_max_bounces, false)
+                    active_count -= 1
+                end
+            end
         end
     end
 
-    return RayTraceResult(false, total, max_steps, path)
+    return RayTraceBatchResult(
+        copy(rb.escaped),
+        copy(rb.total_length),
+        copy(rb.nsteps),
+        copy(rb.termination),
+        copy(rb.p0_termination_distance),
+    )
 end
 
 function find_void_rrhc(
     X::Matrix{Float64}, r::Vector{Float64};
-    restarts=50, steps=2_000, rel_tol=1e-2,
-    lambda=0.05
+    restarts::Int=50, steps::Int=2_000, rel_tol::Float64=1e-2,
+    lambda::Float64=0.05,
 )
-    # This is more robust than extrema for defining the "dense core"
     cx, cy, cz = mean(X[1, :]), mean(X[2, :]), mean(X[3, :])
     sx, sy, sz = stddev(X[1, :]), stddev(X[2, :]), stddev(X[3, :])
 
-    bed_center = SVector(cx, cy, cz)
-    # Restrict search to core 1-sigma region
-    bounds_min = SVector(cx - 0.5 * sx, cy - 0.5 * sy, cz - 0.5 * sz)
-    bounds_max = SVector(cx + 0.5 * sx, cy + 0.5 * sy, cz + 0.5 * sz)
+    bed_center = (cx, cy, cz)
+    bounds_min = (cx - 0.5 * sx, cy - 0.5 * sy, cz - 0.5 * sz)
+    bounds_max = (cx + 0.5 * sx, cy + 0.5 * sy, cz + 0.5 * sz)
     N = length(r)
     mean_r = mean(r)
     success_tol = rel_tol * mean_r
@@ -273,24 +476,23 @@ function find_void_rrhc(
     println("Searching for centralized void. Target gap: $(success_tol) m")
 
     for attempt in 1:restarts
-        # Start exactly at center or very close to it
-        if attempt == 1
-            p = bed_center
+        p = if attempt == 1
+            bed_center
         else
-            px = bounds_min[1] + rand() * (bounds_max[1] - bounds_min[1])
-            py = bounds_min[2] + rand() * (bounds_max[2] - bounds_min[2])
-            pz = bounds_min[3] + rand() * (bounds_max[3] - bounds_min[3])
-            p = SVector(px, py, pz)
+            (
+                bounds_min[1] + rand() * (bounds_max[1] - bounds_min[1]),
+                bounds_min[2] + rand() * (bounds_max[2] - bounds_min[2]),
+                bounds_min[3] + rand() * (bounds_max[3] - bounds_min[3]),
+            )
         end
 
         step_size = mean_r * 0.1
 
-        for s in 1:steps
-            # 1. Find nearest particle surface (Repulsive term)
+        for _ in 1:steps
             min_surf_dist = Inf
             nearest_idx = 0
 
-            for j in 1:N
+            @inbounds for j in 1:N
                 dx = p[1] - X[1, j]
                 dy = p[2] - X[2, j]
                 dz = p[3] - X[3, j]
@@ -302,41 +504,28 @@ function find_void_rrhc(
                 end
             end
 
-            # Check if outside of particles by stated tolerance
             if min_surf_dist > success_tol
                 println("Found valid start point via gradient ascent: $p (dist=$min_surf_dist m)")
                 return p
             end
 
-            # 2. Gradient of Repulsion (Away from surface)
-            c = SVector(X[1, nearest_idx], X[2, nearest_idx], X[3, nearest_idx])
-            vec = p - c
+            c = (X[1, nearest_idx], X[2, nearest_idx], X[3, nearest_idx])
+            vec = sub3(p, c)
             norm_vec = sqrt(dot3(vec, vec))
-
             dir_repulse = if norm_vec < 1e-12
-                d_rand = SVector(rand() - 0.5, rand() - 0.5, rand() - 0.5)
-                d_rand / sqrt(dot3(d_rand, d_rand))
+                d_rand = (rand() - 0.5, rand() - 0.5, rand() - 0.5)
+                normalize3(d_rand)
             else
-                vec * (1.0 / norm_vec)
+                mul3(vec, 1.0 / norm_vec)
             end
 
-            # 3. Gradient of Attraction (Towards Center)
-            vec_center = bed_center - p
+            vec_center = sub3(bed_center, p)
             dist_center = sqrt(dot3(vec_center, vec_center))
-            if dist_center > 1e-9
-                dir_attract = vec_center * (1.0 / dist_center)
-            else
-                dir_attract = SVector(0.0, 0.0, 0.0)
-            end
-
-            # Combine forces: push out of particle vs pull to center.
-            combined_dir = dir_repulse + lambda * dir_attract
+            dir_attract = dist_center > 1e-9 ? mul3(vec_center, 1.0 / dist_center) : (0.0, 0.0, 0.0)
 
-            # Re-normalise direction
-            len_comb = sqrt(dot3(combined_dir, combined_dir))
-            final_dir = combined_dir * (1.0 / len_comb)
+            combined_dir = add3(dir_repulse, mul3(dir_attract, lambda))
+            final_dir = normalize3(combined_dir)
 
-            # Adaptive step logic
             factor = 1.0
             if min_surf_dist < -0.5 * mean_r
                 factor = 2.0
@@ -344,22 +533,16 @@ function find_void_rrhc(
                 factor = 0.5
             end
 
-            move = final_dir * (step_size * factor)
-            p = p + move
-
-            # Clamp not strictly needed if attraction works, but good for safety
-            p = SVector(
+            move = mul3(final_dir, step_size * factor)
+            p = add3(p, move)
+            p = (
                 clamp(p[1], bounds_min[1], bounds_max[1]),
                 clamp(p[2], bounds_min[2], bounds_max[2]),
-                clamp(p[3], bounds_min[3], bounds_max[3])
+                clamp(p[3], bounds_min[3], bounds_max[3]),
             )
-
             step_size *= 0.995
         end
     end
 
     error("Gradient optimization failed to find void space > $(success_tol) m.")
 end
-
-
-
diff --git a/src/utils.jl b/src/utils.jl
index fb8d014..b3c6076 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -1,25 +1,44 @@
-# ------------------------------------------------------------
 # Small vector helpers
 # ------------------------------------------------------------
 @inline dot3(a, b) = @inbounds a[1] * b[1] + a[2] * b[2] + a[3] * b[3]
 
+@inline function cross3(a::NTuple{3,Float64}, b::NTuple{3,Float64})
+    @inbounds return (
+        a[2] * b[3] - a[3] * b[2],
+        a[3] * b[1] - a[1] * b[3],
+        a[1] * b[2] - a[2] * b[1],
+    )
+end
+
 @inline function cross3(a::SVector{3,Float64}, b::SVector{3,Float64})
     @inbounds return SVector{3,Float64}(
-        a[2]*b[3] - a[3]*b[2],
-        a[3]*b[1] - a[1]*b[3],
-        a[1]*b[2] - a[2]*b[1],
+        a[2] * b[3] - a[3] * b[2],
+        a[3] * b[1] - a[1] * b[3],
+        a[1] * b[2] - a[2] * b[1],
     )
 end
 
-@inline function to_mat3x1(v::SVector{3,Float64})
+@inline add3(a::NTuple{3,Float64}, b::NTuple{3,Float64}) = (a[1] + b[1], a[2] + b[2], a[3] + b[3])
+@inline sub3(a::NTuple{3,Float64}, b::NTuple{3,Float64}) = (a[1] - b[1], a[2] - b[2], a[3] - b[3])
+@inline mul3(a::NTuple{3,Float64}, s::Float64) = (a[1] * s, a[2] * s, a[3] * s)
+@inline madd3(a::NTuple{3,Float64}, s::Float64, b::NTuple{3,Float64}) = (a[1] + s * b[1], a[2] + s * b[2], a[3] + s * b[3])
+
+@inline function normalize3(v::NTuple{3,Float64})
+    invn = inv(sqrt(dot3(v, v)))
+    return mul3(v, invn)
+end
+
+@inline function to_mat3x1(v::NTuple{3,Float64})
     M = Matrix{Float64}(undef, 3, 1)
     @inbounds begin
-        M[1,1] = v[1]; M[2,1] = v[2]; M[3,1] = v[3]
+        M[1, 1] = v[1]
+        M[2, 1] = v[2]
+        M[3, 1] = v[3]
     end
     return M
 end
 
-function mean(values::AbstractVector{T}) where T <: Real
+function mean(values::AbstractVector{T}) where {T <: Real}
     if isempty(values)
         return T <: AbstractFloat ? T(NaN) : throw(ArgumentError("mean requires a non-empty vector"))
     end
@@ -30,7 +49,7 @@ function mean(values::AbstractVector{T}) where T <: Real
     return s / length(values)
 end
 
-function stddev(values::AbstractVector{T}) where T<:Real
+function stddev(values::AbstractVector{T}) where {T <: Real}
     if length(values) < 2
         return T <: AbstractFloat ? T(NaN) : throw(ArgumentError("stddev requires at least two elements"))
     end