diff --git a/bin/data/config.json b/bin/data/config.json
index 23d8ef51..aa261711 100644
--- a/bin/data/config.json
+++ b/bin/data/config.json
@@ -117,7 +117,7 @@
 					"gui": true,
 					"vsync": false, // vsync on vulkan side rather than engine-side
 					"hdr": true,
-					"vxgi": false, // true
+					"vxgi": true,
 					"culling": true,
 					"bloom": true,
 					"rt": false,
diff --git a/bin/data/entities/burger.json b/bin/data/entities/burger.json
index 5fe475e9..f7b52d49 100644
--- a/bin/data/entities/burger.json
+++ b/bin/data/entities/burger.json
@@ -29,10 +29,11 @@
 
 		"physics": {
 			"gravity": [ 0, -9.81, 0 ],
-			"inertia": [ 0, 0, 0 ],
 
 			"mass": 10,
-			"type": "bounding box",
+			"type": "sphere",
+			"radius": 1,
+			"static": false,
 			"recenter": false
 		},
 		"graph": {
diff --git a/bin/data/entities/prop.json b/bin/data/entities/prop.json
index b3ed5dee..130962ac 100644
--- a/bin/data/entities/prop.json
+++ b/bin/data/entities/prop.json
@@ -6,7 +6,7 @@
 	"metadata": {
 		"holdable": true,
 		"physics": {
-			"mass": 1000,
+			"mass": 0,
 			"inertia": false,
 			"type": "bounding box"
 		//	"type": "mesh"
diff --git a/bin/data/scenes/sourceengine/mds_mcdonalds.json b/bin/data/scenes/sourceengine/mds_mcdonalds.json
index fbece77f..4b3d2cfb 100644
--- a/bin/data/scenes/sourceengine/mds_mcdonalds.json
+++ b/bin/data/scenes/sourceengine/mds_mcdonalds.json
@@ -2,8 +2,8 @@
 	"import": "./base_sourceengine.json",
 	"assets": [	
 	//	{ "filename": "./models/mds_mcdonalds.glb" }
-		{ "filename": "./models/mds_mcdonalds/graph.json" }
-	//	{ "filename": "/burger.json", "delay": 1 }
+		{ "filename": "./models/mds_mcdonalds/graph.json" },
+		{ "filename": "/burger.json", "delay": 4 }
 	],
 	"metadata": {
 		"graph": {
diff --git a/bin/data/shaders/common/rt.h b/bin/data/shaders/common/rt.h
index 50b55238..1b3fface 100644
--- a/bin/data/shaders/common/rt.h
+++ b/bin/data/shaders/common/rt.h
@@ -71,76 +71,6 @@ vec4 traceStep( Ray ray ) {
 		eyeDepth = surface.position.eye.z;
 	}
 
-#if 0
-	// "transparency"
-	if ( payload.hit && surface.material.albedo.a < 0.999 ) {
-		const vec4 TRANSPARENCY_COLOR = vec4(1.0 - surface.material.albedo.a);
-		
-		if ( surface.material.albedo.a < 0.001 ) outFrag = vec4(0);
-
-		RayTracePayload surfacePayload = payload;
-		Ray transparency;
-		transparency.direction = ray.direction;
-		transparency.origin = surface.position.world;
-		fogRay = transparency;
-
-		trace( transparency, ubo.settings.rt.alphaTestOffset );
-		if ( payload.hit ) {
-			populateSurface( payload );
-			directLighting();
-		} else if (  0 <= ubo.settings.lighting.indexSkybox && ubo.settings.lighting.indexSkybox < CUBEMAPS ) {
-			surface.fragment = texture( samplerCubemaps[ubo.settings.lighting.indexSkybox], ray.direction );
-			surface.fragment.a = 4096;
-			surface.position.eye.z /= 8;
-		}
-	#if FOG
-		fog( transparency, surface.fragment.rgb, surface.fragment.a );
-	#endif
-		outFrag += TRANSPARENCY_COLOR * surface.fragment;
-		eyeDepth = surface.position.eye.z;
-
-		payload = surfacePayload;
-		populateSurface( payload );
-	}
-#if FOG
-	{
-	//	surface.position.eye.z = eyeDepth;
-	//	fog( fogRay, outFrag.rgb, outFrag.a );
-	//	fog( ray, surface.fragment.rgb, surface.fragment.a );
-	}
-#endif
-
-	// reflection
-	if ( payload.hit ) {
-		const float REFLECTIVITY = 1.0 - surface.material.roughness;
-		const vec4 REFLECTED_ALBEDO = surface.material.albedo * REFLECTIVITY;
-
-		if ( REFLECTIVITY > 0.001 ) {
-			RayTracePayload surfacePayload = payload;
-
-			Ray reflection;
-			reflection.origin = surface.position.world;
-			reflection.direction = reflect( ray.direction, surface.normal.world );
-
-			trace( reflection );
-
-			if ( payload.hit ) {
-				populateSurface( payload );
-				directLighting();
-			} else if (  0 <= ubo.settings.lighting.indexSkybox && ubo.settings.lighting.indexSkybox < CUBEMAPS ) {
-				surface.fragment = texture( samplerCubemaps[ubo.settings.lighting.indexSkybox], reflection.direction );
-				surface.fragment.a = 4096;
-			}
-		#if FOG
-			fog( reflection, surface.fragment.rgb, surface.fragment.a );
-		#endif
-			outFrag += REFLECTED_ALBEDO * surface.fragment;
-
-			payload = surfacePayload;
-			populateSurface( payload );
-		}
-	}
-#endif
 	surface = previousSurface;
 	return outFrag;
 }
@@ -170,7 +100,6 @@ void indirectLightingRT() {
 	const vec3 P = surface.position.world;
 	const vec3 N = surface.normal.world;
 
-#if 1
 	const vec3 right = normalize(orthogonal(N));
 	const vec3 up = normalize(cross(right, N));
 
@@ -183,26 +112,7 @@ void indirectLightingRT() {
 		normalize(N + -0.50937f * right + -0.7006629f * up),
 		normalize(N + -0.823639f * right + 0.267617f * up),
 	};
-#else
-	const vec3 ortho = normalize(orthogonal(N));
-	const vec3 ortho2 = normalize(cross(ortho, N));
 
-	const vec3 corner = 0.5f * (ortho + ortho2);
-	const vec3 corner2 = 0.5f * (ortho - ortho2);
-
-	const uint CONES_COUNT = 9;
-	const vec3 CONES[] = {
-		N,
-		normalize(mix(N, ortho, 0.5)),
-		normalize(mix(N, -ortho, 0.5)),
-		normalize(mix(N, ortho2, 0.5)),
-		normalize(mix(N, -ortho2, 0.5)),
-		normalize(mix(N, corner, 0.5)),
-		normalize(mix(N, -corner, 0.5)),
-		normalize(mix(N, corner2, 0.5)),
-		normalize(mix(N, -corner2, 0.5)),
-	};
-#endif
 
 	const float DIFFUSE_CONE_APERTURE = 2.0 * 0.57735f;
 	const float DIFFUSE_INDIRECT_FACTOR = 0; // 1.0f / float(CONES_COUNT) * 0.125f;
diff --git a/bin/data/shaders/common/shadows.h b/bin/data/shaders/common/shadows.h
index 34c74bd7..b8d1f27c 100644
--- a/bin/data/shaders/common/shadows.h
+++ b/bin/data/shaders/common/shadows.h
@@ -142,10 +142,11 @@ float shadowFactor( const Light light, float def ) {
 	const float bias = light.depthBias;
 	const float eyeDepth = positionClip.z;
 	const int samples = int(SHADOW_SAMPLES);
+	const float invResolution = 1.0 / textureSize( samplerTextures[nonuniformEXT(light.indexMap)], 0 ).x;
 	if ( samples < 1 ) return eyeDepth < texture(samplerTextures[nonuniformEXT(light.indexMap)], uv).r - bias ? 0.0 : factor;
 	for ( int i = 0; i < samples; ++i ) {
 		const int index = int( float(samples) * random(floor(surface.position.world.xyz * 1000.0), i)) % samples;
-		const float lightDepth = texture(samplerTextures[nonuniformEXT(light.indexMap)], uv + poissonDisk[index] / 700.0 ).r;
+		const float lightDepth = texture(samplerTextures[nonuniformEXT(light.indexMap)], uv + poissonDisk[index] * invResolution ).r;
 		if ( eyeDepth < lightDepth - bias ) factor -= 1.0 / samples;
 	}
 	return factor;
diff --git a/bin/data/shaders/common/vxgi.h b/bin/data/shaders/common/vxgi.h
index 3829baa3..c2c2e01a 100644
--- a/bin/data/shaders/common/vxgi.h
+++ b/bin/data/shaders/common/vxgi.h
@@ -135,11 +135,11 @@ void indirectLightingVXGI() {
 	};
 #endif
 
-	const float DIFFUSE_CONE_APERTURE = 2.0 * 0.57735f;
-	const float DIFFUSE_INDIRECT_FACTOR = 0; // 1.0f / float(CONES_COUNT) * 0.125f;
+	const float DIFFUSE_CONE_APERTURE = 2.0f * 0.57735f;
+	const float DIFFUSE_INDIRECT_FACTOR = 1.0f / float(CONES_COUNT) * 0.125f;
 	
 	const float SPECULAR_CONE_APERTURE = clamp(tan(PI * 0.5f * surface.material.roughness), 0.0174533f, PI); // tan( R * PI * 0.5f * 0.1f );
-	const float SPECULAR_INDIRECT_FACTOR = (1.0 - surface.material.metallic); // * 0.25; // 1.0f;
+	const float SPECULAR_INDIRECT_FACTOR = (1.0f - surface.material.metallic) * (1.0f - surface.material.roughness); // * 0.25; // 1.0f;
 	
 	vec4 indirectDiffuse = vec4(0);
 	vec4 indirectSpecular = vec4(0);
diff --git a/bin/data/shaders/display/vxgi/comp.glsl b/bin/data/shaders/display/vxgi/comp.glsl
index e9078a0e..2d49e312 100644
--- a/bin/data/shaders/display/vxgi/comp.glsl
+++ b/bin/data/shaders/display/vxgi/comp.glsl
@@ -10,9 +10,9 @@ layout (local_size_x = 8, local_size_y = 8, local_size_z = 8) in;
 
 #define VXGI 1
 #define MAX_CUBEMAPS CUBEMAPS
-#define GAMMA_CORRECT 1
-#define PBR 1
-#define LAMBERT 0
+#define GAMMA_CORRECT 0
+#define PBR 0
+#define LAMBERT 1
 
 layout (constant_id = 0) const uint TEXTURES = 512;
 layout (constant_id = 1) const uint CUBEMAPS = 128;
@@ -165,7 +165,6 @@ void main() {
 		#if 0
 			pbr();
 		#else
-		/*
 			surface.material.roughness *= 4.0;
 			const vec3 F0 = mix(vec3(0.04), surface.material.albedo.rgb, surface.material.metallic); 
 			const vec3 Lo = normalize( surface.position.world );
@@ -185,7 +184,7 @@ void main() {
 			//	if ( light.power * La * Ls <= LIGHT_POWER_CUTOFF ) continue;
 
 				const float cosLi = max(0.0, dot(surface.normal.world, Li));
-				const vec3 Lr = light.color.rgb * light.power * La * Ls;
+				const vec3 Lr = light.color.rgb * light.power * La * Ls / PI;
 			#if LAMBERT
 				const vec3 diffuse = surface.material.albedo.rgb;
 				const vec3 specular = vec3(0);
@@ -203,7 +202,6 @@ void main() {
 				surface.light.rgb += (diffuse + specular) * Lr * cosLi;
 				surface.light.a += light.power * La * Ls;
 			}
-		*/
 		#endif
 		}
 		surface.fragment.rgb += surface.light.rgb;
diff --git a/engine/inc/uf/utils/math/matrix.h b/engine/inc/uf/utils/math/matrix.h
index 14d68870..dc55a0c0 100644
--- a/engine/inc/uf/utils/math/matrix.h
+++ b/engine/inc/uf/utils/math/matrix.h
@@ -78,6 +78,9 @@ namespace uf {
 		template<typename T=NUM> /*FORCE_INLINE*/ pod::Vector3t<T> multiply( const pod::Matrix4t<T>& mat, const pod::Vector3t<T>& vector, T w = 1, bool = false );
 		template<typename T=NUM> /*FORCE_INLINE*/ pod::Vector4t<T> multiply( const pod::Matrix4t<T>& mat, const pod::Vector4t<T>& vector, bool = false );
 		template<typename T=NUM, size_t M, size_t N> /*FORCE_INLINE*/ pod::Matrix<T,M,N> /*UF_API*/ outerProduct( const pod::Vector<T,M>& a, const pod::Vector<T,N>& b );
+		
+		template<typename T=NUM, size_t R, size_t C> /*FORCE_INLINE*/ pod::Vector<T, R> diagonal(const pod::Matrix<T, R, C>& mat );
+		template<typename T=NUM, size_t N> /*FORCE_INLINE*/ pod::Matrix<T, N, N> diagonal(const pod::Vector<T, N>& vec );
 
 		template<typename T=pod::Matrix4> /*FORCE_INLINE*/ T /*UF_API*/ multiplyAll( const T& matrix, typename T::type_t scalar );
 		template<typename T=pod::Matrix4> /*FORCE_INLINE*/ T /*UF_API*/ add( const T& lhs, const T& rhs );
diff --git a/engine/inc/uf/utils/math/matrix/pod.inl b/engine/inc/uf/utils/math/matrix/pod.inl
index 75f8b7c2..0d6a1708 100644
--- a/engine/inc/uf/utils/math/matrix/pod.inl
+++ b/engine/inc/uf/utils/math/matrix/pod.inl
@@ -169,6 +169,21 @@ template<typename T, size_t M, size_t N> pod::Matrix<T,M,N> /*UF_API*/ uf::matri
     return m;
 }
 
+template<typename T, size_t R, size_t C> pod::Vector<T, R> /*UF_API*/ uf::matrix::diagonal(const pod::Matrix<T, R, C>& mat ) {
+	pod::Vector<T, R> vector;
+	FOR_EACH(R, {
+		vector[i] = m(i, i);
+	});
+	return vector;
+}
+template<typename T, size_t N> pod::Matrix<T, N, N> /*UF_API*/ uf::matrix::diagonal(const pod::Vector<T, N>& vector ) {
+	pod::Matrix<T, N, N> matrix;
+	FOR_EACH(N, {
+		matrix(i, i) = vector[i];
+	});
+	return matrix;
+}
+
 template<typename T> T /*UF_API*/ uf::matrix::multiplyAll( const T& m, typename T::type_t scalar ) {
 	T matrix;
 
diff --git a/engine/inc/uf/utils/math/quaternion.h b/engine/inc/uf/utils/math/quaternion.h
index 546bb76f..18fd38f0 100644
--- a/engine/inc/uf/utils/math/quaternion.h
+++ b/engine/inc/uf/utils/math/quaternion.h
@@ -40,6 +40,7 @@ namespace uf {
 		template<typename T> /*FORCE_INLINE*/ T /*UF_API*/ normalize( const T& vector );
 
 		template<typename T> /*FORCE_INLINE*/ pod::Matrix4t<T> matrix( const pod::Quaternion<T>& quaternion );
+		template<typename T> /*FORCE_INLINE*/ pod::Matrix3t<T> matrix3( const pod::Quaternion<T>& quaternion );
 		template<typename T> /*FORCE_INLINE*/ pod::Quaternion<T> axisAngle( const pod::Vector3t<T>& axis, T angle );
 		template<typename T> /*FORCE_INLINE*/ pod::Quaternion<T> unitVectors( const pod::Vector3t<T>& u, const pod::Vector3t<T>& v );
 		template<typename T> /*FORCE_INLINE*/ pod::Quaternion<T> lookAt( const pod::Vector3t<T>& source, const pod::Vector3t<T>& destination );
diff --git a/engine/inc/uf/utils/math/quaternion/pod.inl b/engine/inc/uf/utils/math/quaternion/pod.inl
index 60c4c9f2..312570ef 100644
--- a/engine/inc/uf/utils/math/quaternion/pod.inl
+++ b/engine/inc/uf/utils/math/quaternion/pod.inl
@@ -126,6 +126,27 @@ template<typename T> pod::Matrix4t<T> uf::quaternion::matrix( const pod::Quatern
 		0,             0,             0,         1
 	});
 }
+template<typename T> pod::Matrix3t<T> uf::quaternion::matrix3( const pod::Quaternion<T>& q ) {
+	auto normal = uf::quaternion::normalize(q);
+
+	const T xx = 2 * normal.x * normal.x;
+	const T xy = 2 * normal.x * normal.y;
+	const T xz = 2 * normal.x * normal.z;
+	const T xw = 2 * normal.x * normal.w;
+
+	const T yy = 2 * normal.y * normal.y;
+	const T yz = 2 * normal.y * normal.z;
+	const T yw = 2 * normal.y * normal.w;
+
+	const T zz = 2 * normal.z * normal.z;
+	const T zw = 2 * normal.z * normal.w;
+
+	return pod::Matrix3t<T>({
+		1 - yy - zz,   xy + zw,       xz - yw,
+		xy - zw,       1 - xx - zz,   yz + xw,
+		xz + yw,       yz - xw,       1 - xx - yy,
+	});
+}
 template<typename T> pod::Quaternion<T> uf::quaternion::axisAngle( const pod::Vector3t<T>& axis, T angle ) { 
 	pod::Quaternion<T> q;
 
diff --git a/engine/src/utils/math/physics/aabb.inl b/engine/src/utils/math/physics/aabb.inl
index 7dbd6ada..691b2811 100644
--- a/engine/src/utils/math/physics/aabb.inl
+++ b/engine/src/utils/math/physics/aabb.inl
@@ -139,11 +139,12 @@ namespace {
 		auto f1 = v2 - v1;
 		auto f2 = v0 - v2;
 
-		// sAT: test the 9 edge cross axes
+		// SAT: test the 9 edge cross axes
 		auto axisTest = [&]( const pod::Vector3f& axis ) {
-			if ( uf::vector::magnitude(axis) < eps ) return true;
+			float norm = uf::vector::norm( axis );
+			if ( norm < eps ) return true;
 			
-			auto a = uf::vector::normalize( axis );
+			auto a = axis / norm;
 			
 			float p0 = uf::vector::dot( v0, a );
 			float p1 = uf::vector::dot( v1, a );
diff --git a/engine/src/utils/math/physics/capsule.inl b/engine/src/utils/math/physics/capsule.inl
index 8d95b9ce..e31f34e4 100644
--- a/engine/src/utils/math/physics/capsule.inl
+++ b/engine/src/utils/math/physics/capsule.inl
@@ -9,7 +9,7 @@ namespace {
 		float r = a.collider.capsule.radius + b.collider.capsule.radius;
 		
 		auto delta = pB - pA;
-		float dist2 = uf::vector::dot(delta, delta);
+		float dist2 = uf::vector::magnitude( delta );
 		if ( dist2 > r * r ) return false;
 		float dist = std::sqrt( dist2 );
 		
@@ -63,18 +63,6 @@ namespace {
 		float penetration = r - dist;
 		manifold.points.emplace_back(pod::Contact{ contact, normal, penetration });
 
-		// build up to 4 contact points sampled around foot circle
-	#if 0
-		auto tangent1 = uf::vector::normalize( uf::vector::cross( normal, pod::Vector3f{1,0,0} ) );
-		if ( uf::vector::magnitude( tangent1 ) < EPS(1e-6f) ) tangent1 = uf::vector::normalize( uf::vector::cross( normal, pod::Vector3f{0,1,0} ) );
-		auto tangent2 = uf::vector::cross( normal, tangent1 );
-
-		// four directions around circle
-		manifold.points.emplace_back(pod::Contact{ foot + tangent1 * r - normal * d, normal, penetration });
-		manifold.points.emplace_back(pod::Contact{ foot - tangent1 * r - normal * d, normal, penetration });
-		manifold.points.emplace_back(pod::Contact{ foot + tangent2 * r - normal * d, normal, penetration });
-		manifold.points.emplace_back(pod::Contact{ foot - tangent2 * r - normal * d, normal, penetration });
-	#endif
 		return true;
 	}
 	bool capsuleSphere( const pod::PhysicsBody& a, const pod::PhysicsBody& b, pod::Manifold& manifold, float eps ) {
diff --git a/engine/src/utils/math/physics/gjk.inl b/engine/src/utils/math/physics/gjk.inl
index 86b60bd3..8e6011d9 100644
--- a/engine/src/utils/math/physics/gjk.inl
+++ b/engine/src/utils/math/physics/gjk.inl
@@ -124,9 +124,10 @@ namespace {
 	}
 
 	bool isDegenerate( const pod::Simplex& s, const pod::SupportPoint& newPt, float eps = EPS(1.0e-6f) ) {
+		const float eps2 = eps * eps;
 		// compare to existing
 		for ( auto& sp : s.pts ) {
-			if ( uf::vector::magnitude( sp.p - newPt.p ) < eps ) return true; // duplicate point
+			if ( uf::vector::magnitude( sp.p - newPt.p ) < eps2 ) return true; // duplicate point
 		}
 		// if simplex already has 3 points, check area
 		if ( s.pts.size() >= 2 ) {
@@ -135,7 +136,7 @@ namespace {
 			auto& C = newPt.p;
 			
 			float area = uf::vector::magnitude( uf::vector::cross( B - A, C - A ) );
-			if ( area < eps ) return true; // collinear with previous
+			if ( area < eps2 ) return true; // collinear with previous
 		}
 		if ( s.pts.size() >= 3 ) {
 			auto& A = s.pts[0].p;
@@ -150,8 +151,9 @@ namespace {
 	}
 
 	bool gjk( const pod::PhysicsBody& a, const pod::PhysicsBody& b, pod::Simplex& simplex, int maxIterations = 20, float eps = EPS(1e-6f) ) {
+		const float eps2 = eps * eps;
 		auto dir = ::getPosition( b ) - ::getPosition( a );
-		if ( uf::vector::magnitude(dir) < eps ) dir = {1,0,0}; // fallback direction
+		if ( uf::vector::magnitude( dir ) < eps2 ) dir = {1,0,0}; // fallback direction
 
 		// initial condition
 		simplex.pts.clear();
diff --git a/engine/src/utils/math/physics/helpers.inl b/engine/src/utils/math/physics/helpers.inl
index b4fef2bf..eb0e8bcf 100644
--- a/engine/src/utils/math/physics/helpers.inl
+++ b/engine/src/utils/math/physics/helpers.inl
@@ -58,10 +58,12 @@ namespace {
 namespace {
 	// create ID from pointers
 	uint64_t makePairKey( const pod::PhysicsBody& a, const pod::PhysicsBody& b ) {
-		auto idA = reinterpret_cast<uint64_t>(&a);
-		auto idB = reinterpret_cast<uint64_t>(&b);
-		if ( idA > idB ) std::swap(idA, idB); // ensure consistent order
-		return (idA << 32) ^ idB;
+		uint64_t lhs = reinterpret_cast<uint64_t>(&a);
+		uint64_t rhs = reinterpret_cast<uint64_t>(&b);
+		if (lhs > rhs) std::swap(lhs, rhs);
+
+		lhs ^= rhs + 0x9e3779b97f4a7c15 + (lhs << 6) + (lhs >> 2);
+		return lhs;
 	}
 
 	// marks a body as asleep
@@ -118,6 +120,15 @@ namespace {
 		return ::shouldCollide( a.collider, b.collider );
 	}
 
+	pod::Matrix3f computeWorldInverseInertia( const pod::PhysicsBody& b ) {
+		if ( b.isStatic || b.inverseMass == 0.0f ) return pod::Matrix3f{};
+
+		pod::Matrix3f invI_local = uf::matrix::diagonal( b.inverseInertiaTensor );
+		pod::Matrix3f R = uf::quaternion::matrix3(b.transform->orientation);
+
+		return R * invI_local * uf::matrix::transpose(R);
+	}
+
 	// normalizes the delta between two bodies / contacts by the distance (as it was already computed) if non-zero
 	// a lot of collider v colliders use this semantic
 	pod::Vector3f normalizeDelta( const pod::Vector3f& delta, float dist, float eps = EPS(1.0e-6), const pod::Vector3f& fallback = pod::Vector3f{0,1,0} ) {
@@ -183,16 +194,17 @@ namespace {
 	}
 
 	pod::Vector3f closestPointSegmentAabb( const pod::Vector3f& p1, const pod::Vector3f& p2, const pod::AABB& box, float eps = EPS(1e-6f) ) {
+		const float eps2 = eps * eps;
 		// AABB center and half extents
 		auto c = ( box.min + box.max ) * 0.5f;
 		auto e = ( box.max - box.min ) * 0.5f;
 
 		// direction of line segment
 		auto d = p2 - p1;
-		float len2 = uf::vector::dot( d, d );
+		float len2 = uf::vector::magnitude( d );
 		float t = 0.0f;
 
-		if ( len2 > eps ) {
+		if ( len2 > eps2 ) {
 			// parametric closest t from box center
 			t = uf::vector::dot( c - p1, d ) / len2;
 			t = std::clamp( t, 0.0f, 1.0f );
diff --git a/engine/src/utils/math/physics/impl.cpp b/engine/src/utils/math/physics/impl.cpp
index 1b6b4239..0e8e1982 100644
--- a/engine/src/utils/math/physics/impl.cpp
+++ b/engine/src/utils/math/physics/impl.cpp
@@ -3,6 +3,7 @@
 #include <uf/engine/scene/scene.h>
 #include <uf/utils/mesh/mesh.h>
 #include <uf/utils/memory/stack.h>
+#include <mutex>
 
 namespace {
 	bool warmupSolver = true; // cache manifold data to warm up the solver
@@ -32,6 +33,7 @@ namespace {
 	float baumgarteCorrectionPercent = 0.4f;
 	float baumgarteCorrectionSlop = 0.01f;
 	
+	std::mutex cacheMutex;
 	uf::stl::unordered_map<size_t, pod::Manifold> manifoldsCache;
 	uint32_t manifoldCacheLifetime = 6; // to-do: find a good value for this
 
@@ -48,7 +50,8 @@ namespace {
 	float groundedThreshold = 0.7f; // threshold before marking a body as grounded
 }
 
-#define EPS(x) x // 1.0e-6f
+#define EPS(x) 1.0e-6f
+#define EPS2 (EPS(1.0e-6) * EPS(1.0e-6))
 #define ASSERT_COLLIDER_TYPES( A, B ) UF_ASSERT( a.collider.type == pod::ShapeType::A && b.collider.type == pod::ShapeType::B );
 #define UF_PHYSICS_TEST 0
 
@@ -217,7 +220,10 @@ void uf::physics::impl::step( pod::World& world, float dt ) {
 		// do position correction
 		::solvePositions( manifolds, dt );
 		// cache manifold positions
-		if ( ::warmupSolver ) ::storeManifolds( manifolds, ::manifoldsCache );
+		if ( ::warmupSolver ) {
+            std::lock_guard<std::mutex> lock(::cacheMutex);
+            ::storeManifolds( manifolds, ::manifoldsCache );
+        }
 	}
 
 	for ( auto* b : bodies ) {
@@ -278,11 +284,29 @@ void uf::physics::impl::updateInertia( pod::PhysicsBody& body ) {
 
 	switch ( body.collider.type ) {
 		case pod::ShapeType::AABB: {
+			/*
+			// old
 			pod::Vector3f extents = (body.collider.aabb.max - body.collider.aabb.min);
 			extents *= extents; // square it;
 
 			body.inertiaTensor = extents * (body.mass / 12.0f);
 			body.inverseInertiaTensor = 1.0f / body.inertiaTensor;
+			*/
+			pod::Vector3f dims = (body.collider.aabb.max - body.collider.aabb.min);
+			pod::Vector3f dimsSq = dims * dims;
+
+			float massFactor = body.mass / 12.0f;
+			body.inertiaTensor = {
+				massFactor * (dimsSq.y + dimsSq.z),
+				massFactor * (dimsSq.x + dimsSq.z),
+				massFactor * (dimsSq.x + dimsSq.y),
+			};
+
+			body.inertiaTensor.x = std::max(body.inertiaTensor.x, EPS(1e-6f));
+			body.inertiaTensor.y = std::max(body.inertiaTensor.y, EPS(1e-6f));
+			body.inertiaTensor.z = std::max(body.inertiaTensor.z, EPS(1e-6f));
+
+			body.inverseInertiaTensor = 1.0f / body.inertiaTensor;
 		} break;
 		case pod::ShapeType::SPHERE: {
 			float I = 0.4f * body.mass * body.collider.sphere.radius * body.collider.sphere.radius;
@@ -309,13 +333,19 @@ void uf::physics::impl::updateInertia( pod::PhysicsBody& body ) {
 			float totalVolume = 0.0f;
 
 			// compute total volume
-			for ( auto& box : bvh.bounds ) {
+			for ( size_t i = 0; i < bvh.nodes.size(); ++i ) {
+				if ( bvh.nodes[i].getCount() == 0 ) continue;
+				const auto& box = bvh.bounds[i];
+
 				auto extents = box.max - box.min;
 				totalVolume += extents.x * extents.y * extents.z;
 			}
 
 			// accumulate inertia
-			for ( auto& box : bvh.bounds ) {
+			for ( size_t i = 0; i < bvh.nodes.size(); ++i ) {
+				if ( bvh.nodes[i].getCount() == 0 ) continue;
+				const auto& box = bvh.bounds[i];
+
 				auto extents = box.max - box.min;
 				float mass = body.mass * extents.x * extents.y * extents.z / totalVolume;
 
@@ -332,9 +362,9 @@ void uf::physics::impl::updateInertia( pod::PhysicsBody& body ) {
 				// parallel axis theorem
 				pod::Vector3f center = (box.min + box.max) * 0.5f;
 				pod::Vector3f d = center; // relative to mesh COM (assume COM at origin for now)
-				float d2 = uf::vector::dot(d, d);
+				float dist2 = uf::vector::magnitude( d );
 
-				pod::Matrix3f pat = uf::matrix::identityi<pod::Matrix3f>() * (mass * d2);
+				pod::Matrix3f pat = uf::matrix::identityi<pod::Matrix3f>() * (mass * dist2);
 				pat -= uf::matrix::outerProduct(d, d) * mass;
 
 				inertia += Ibox + pat;
diff --git a/engine/src/utils/math/physics/integration.inl b/engine/src/utils/math/physics/integration.inl
index 7758a2b3..c042a8b8 100644
--- a/engine/src/utils/math/physics/integration.inl
+++ b/engine/src/utils/math/physics/integration.inl
@@ -25,12 +25,10 @@ namespace {
 		if ( !a.isStatic ) {
 			a.velocity -= impulse * a.inverseMass;
 			a.angularVelocity -= (uf::vector::cross(rA, impulse)) * a.inverseInertiaTensor;
-		//	if ( uf::vector::magnitude( impulse ) > 1.0e-4 ) UF_MSG_DEBUG("aV delta={}", uf::vector::toString(impulse * a.inverseMass));
 		}
 		if ( !b.isStatic ) {
 			b.velocity += impulse * b.inverseMass;
 			b.angularVelocity += (uf::vector::cross(rB, impulse)) * b.inverseInertiaTensor;
-		//	if ( uf::vector::magnitude( impulse ) > 1.0e-4 ) UF_MSG_DEBUG("bV delta={}", uf::vector::toString(impulse * b.inverseMass));
 		}
 	}
 
@@ -38,8 +36,8 @@ namespace {
 		if ( body.isStatic ) return;
 
 		float rollingFriction = 0.02f; // to-do: derive from material
-		float angularSpeed = uf::vector::magnitude( body.angularVelocity );
-		if ( angularSpeed < EPS(1.0e-6f) ) return;
+		float angularSpeed2 = uf::vector::magnitude( body.angularVelocity );
+		if ( angularSpeed2 < EPS2 ) return;
 
 		body.angularVelocity += body.angularVelocity * body.mass * -rollingFriction * dt;
 		// body.angularVelocity *= -rollingFriction * dt;
@@ -274,14 +272,22 @@ namespace {
 		body.velocity += acceleration * dt;
 
 		// angular integration
-		body.angularVelocity += body.torqueAccumulator * body.inverseInertiaTensor * dt;
+		//body.angularVelocity += body.torqueAccumulator * body.inverseInertiaTensor * dt;
+		{
+			pod::Matrix3f R = uf::quaternion::matrix3(body.transform->orientation);
+			pod::Vector3f localTorque = uf::matrix::multiply( uf::matrix::transpose(R), body.torqueAccumulator );
+			pod::Vector3f localAngAccel = localTorque * body.inverseInertiaTensor; // element-wise
+			body.angularVelocity += uf::matrix::multiply( R, localAngAccel ) * dt;
+		}
 
 		// update position
 		body.transform/*.reference*/->position += body.velocity * dt;
 
 		// update orientation
-		if ( uf::vector::magnitude( body.angularVelocity ) > EPS(1.0e-8f) ) {
-			pod::Quaternion<> dq = uf::quaternion::axisAngle(uf::vector::normalize(body.angularVelocity), uf::vector::magnitude(body.angularVelocity)*dt);
+		float angularSpeed2 = uf::vector::magnitude( body.angularVelocity );
+		if ( angularSpeed2 > EPS2 ) {
+			float angularSpeed = std::sqrt( angularSpeed2 );
+			pod::Quaternion<> dq = uf::quaternion::axisAngle( body.angularVelocity / angularSpeed, angularSpeed * dt);
 			uf::transform::rotate( *body.transform/*.reference*/, dq );
 		}
 
diff --git a/engine/src/utils/math/physics/ray.inl b/engine/src/utils/math/physics/ray.inl
index 6c1d9e3b..5ea89686 100644
--- a/engine/src/utils/math/physics/ray.inl
+++ b/engine/src/utils/math/physics/ray.inl
@@ -146,7 +146,7 @@ namespace {
 
 		float b = dd * uf::vector::dot(m,n) - nd*md;
 		float discr = b*b - a*c;
-		if ( discr < 0.0f ) return false;
+		if ( fabs(a) < EPS(1.0e-6f) || discr < 0.0f ) return false;
 
 		float t = (-b - std::sqrt(discr)) / a;  // nearer hit
 
diff --git a/engine/src/utils/math/physics/solvers.inl b/engine/src/utils/math/physics/solvers.inl
index 9625e513..bd588956 100644
--- a/engine/src/utils/math/physics/solvers.inl
+++ b/engine/src/utils/math/physics/solvers.inl
@@ -38,9 +38,9 @@ namespace {
 
 		// tangent direction
 		pod::Vector3f tangent = rv - contact.normal * uf::vector::dot(rv, contact.normal);
-		float tangentMag = uf::vector::magnitude(tangent);
-		if (tangentMag > EPS(1e-6f)) {
-			tangent /= tangentMag;
+		float tangentMag2 = uf::vector::magnitude(tangent);
+		if ( tangentMag2 > EPS2 ) {
+			tangent /= std::sqrt( tangentMag2 );
 
 			// effective mass along tangent
 			float invMassT = ::computeEffectiveMass(a, b, rA, rB, tangent);
@@ -104,8 +104,11 @@ namespace {
 				pod::Vector3f raXnj = uf::vector::cross(rA_j, n_j);
 				pod::Vector3f rbXnj = uf::vector::cross(rB_j, n_j);
 
-				pod::Vector3f Ia_raXnj = a.inverseInertiaTensor * raXnj;
-				pod::Vector3f Ib_rbXnj = b.inverseInertiaTensor * rbXnj;
+				pod::Matrix3f invIa = computeWorldInverseInertia( a );
+				pod::Matrix3f invIb = computeWorldInverseInertia( b );
+
+				pod::Vector3f Ia_raXnj = uf::matrix::multiply( invIa, raXnj );
+				pod::Vector3f Ib_rbXnj = uf::matrix::multiply( invIb, rbXnj );
 
 				pod::Vector3f crossA = uf::vector::cross(Ia_raXnj, rA_i);
 				pod::Vector3f crossB = uf::vector::cross(Ib_rbXnj, rB_i);
@@ -141,6 +144,9 @@ namespace {
 			float penetrationBias = std::max(contact.penetration - ::baumgarteCorrectionSlop, 0.0f) * (::baumgarteCorrectionPercent / dt);
 			penetrationBias = std::min(penetrationBias, 2.0f / dt); // clamp
 
+			float maxPenetrationRecovery = 2.0f; // Limit to 2 units per second
+			if ( penetrationBias > maxPenetrationRecovery ) penetrationBias = maxPenetrationRecovery;
+
 			float cDot = vRel + penetrationBias;
 
 			rhs[i]	= (cDot < 0.0f) ? -cDot : 0.0f; // rHS is magnitude of correction needed
diff --git a/engine/src/utils/math/physics/sphere.inl b/engine/src/utils/math/physics/sphere.inl
index 4e28c57b..463cf1e2 100644
--- a/engine/src/utils/math/physics/sphere.inl
+++ b/engine/src/utils/math/physics/sphere.inl
@@ -5,7 +5,7 @@ namespace {
 		auto delta = ::getPosition( b ) - ::getPosition( a );
 		float r = a.collider.sphere.radius + b.collider.sphere.radius;
 
-		float dist2 = uf::vector::dot(delta, delta);
+		float dist2 = uf::vector::magnitude( delta );
 		if ( dist2 > r * r ) return false;
 		float dist = std::sqrt( dist2 );
 		
@@ -33,20 +33,13 @@ namespace {
 		};
 
 		auto delta = center - closest;
-		float dist2 = uf::vector::dot(delta, delta);
+		float dist2 = uf::vector::magnitude( delta );
 		if ( dist2 > r * r ) return false;
 		float dist = std::sqrt( dist2 );
 
 		float penetration = r - dist;
 		auto normal = ::normalizeDelta( delta, dist, eps );
 		auto contact = center - normal * r;
-
-		#if 0
-			pod::Vector3f dir = center - ::aabbCenter( bounds );
-			if (std::fabs(dir.x) > std::fabs(dir.y) && std::fabs(dir.x) > std::fabs(dir.z)) normal = { (dir.x > 0 ? 1 : -1), 0, 0 };
-			else if (std::fabs(dir.y) > std::fabs(dir.z)) normal = { 0, (dir.y > 0 ? 1 : -1), 0 };
-			else normal = { 0, 0, (dir.z > 0 ? 1 : -1) };
-		#endif
 		
 		manifold.points.emplace_back(pod::Contact{ contact, normal, penetration });
 		return true;
diff --git a/engine/src/utils/math/physics/triangle.inl b/engine/src/utils/math/physics/triangle.inl
index d6434651..9fddad07 100644
--- a/engine/src/utils/math/physics/triangle.inl
+++ b/engine/src/utils/math/physics/triangle.inl
@@ -16,6 +16,7 @@ namespace {
 	}
 
 	bool triangleTriangleIntersect( const pod::Triangle& a, const pod::Triangle& b, float eps = EPS(1e-6f) ) {
+		const float eps2 = eps * eps;
 		auto boxA = ::computeTriangleAABB( a );
 		auto boxB = ::computeTriangleAABB( b );
 
@@ -24,11 +25,11 @@ namespace {
 		// check vertices of a inside b or vice versa
 		FOR_EACH(3, {
 			auto q = ::closestPointOnTriangle( a.points[i], b );
-			if ( uf::vector::magnitude( q - a.points[i] ) < eps ) return true;
+			if ( uf::vector::magnitude( q - a.points[i] ) < eps2 ) return true;
 		});
 		FOR_EACH(3, {
 			auto q = ::closestPointOnTriangle( b.points[i], a );
-			if ( uf::vector::magnitude( q - b.points[i] ) < eps ) return true;
+			if ( uf::vector::magnitude( q - b.points[i] ) < eps2 ) return true;
 		});
 		return false;
 	}
@@ -239,6 +240,8 @@ namespace {
 	// i feel like it'd just be better to treat an AABB as a 12-triangle mesh and do a triangleTriangle collision instead of 3 + 1 + 3 * 3 axis tests
 	// but what do i know
 	bool triangleAabbSAT( const pod::TriangleWithNormal& tri, const pod::PhysicsBody& body, pod::Manifold& manifold, float eps = 1e-6f ) {
+		const float eps2 = eps * eps;
+
 		const auto& aabb = body.bounds;
 
 		// box center and half extents
@@ -259,7 +262,7 @@ namespace {
 		pod::Vector3f bestAxis;
 
 		auto testAxis = [&](const pod::Vector3f& axis) -> bool {
-			if ( uf::vector::magnitude(axis) < eps ) return true; // skip degenerate
+			if ( uf::vector::magnitude( axis ) < eps2 ) return true; // skip degenerate
 
 			pod::Vector3f n = uf::vector::normalize(axis);
 
@@ -297,17 +300,26 @@ namespace {
 		}
 
 		pod::Vector3f triNormal = uf::vector::normalize(uf::vector::cross(e0, e1));
+		float planeDist = uf::vector::dot(triNormal, v0);
+		if ( uf::vector::dot(bestAxis, triNormal) < 0.0f ) bestAxis = -bestAxis;
+		pod::Vector3f contact = boxCenter - bestAxis * (boxHalf.x * fabs(bestAxis.x) + boxHalf.y * fabs(bestAxis.y) + boxHalf.z * fabs(bestAxis.z));
+		//pod::Vector3f contact = boxCenter - triNormal * planeDist;
+	
+	/*
 		float d = uf::vector::dot(triNormal, v0);
 		float dist = uf::vector::dot(triNormal, -boxCenter) - d;
 		pod::Vector3f contact = boxCenter - triNormal * dist;
 
 		if ( uf::vector::dot(bestAxis, triNormal) < 0.0f ) bestAxis = -bestAxis;
+	*/
 
 		manifold.points.emplace_back( pod::Contact{ contact, bestAxis, minOverlap } );
 		return true;
 	}
 
 	bool triangleTriangle( const pod::TriangleWithNormal& a, const pod::TriangleWithNormal& b, pod::Manifold& manifold, float eps ) {
+		const float eps2 = eps * eps;
+
 		size_t axes = 0;
 		pod::Vector3f axesBuffer[12];
 		axesBuffer[axes++] = ::triangleNormal(a);
@@ -318,7 +330,7 @@ namespace {
 			for (int j = 0; j < 3; j++) {
 				auto eb = b.points[(j+1)%3] - b.points[j];
 				auto axis = uf::vector::cross(ea, eb);
-				if ( uf::vector::magnitude(axis) > eps ) axesBuffer[axes++] = uf::vector::normalize(axis);
+				if ( uf::vector::magnitude( axis ) > eps2 ) axesBuffer[axes++] = uf::vector::normalize(axis);
 			}
 		}
 
@@ -428,36 +440,8 @@ namespace {
 		return hit;
 	}
 	bool triangleAabb( const pod::TriangleWithNormal& tri, const pod::PhysicsBody& body, pod::Manifold& manifold, float eps ) {
-		const auto& aabb = body;
-
-	#if 1
-		return ::triangleAabbTri( tri, body, manifold, eps );
-		//return ::triangleAabbSAT( tri, body, manifold, eps );
-	#else
-		auto closest = ::closestPointOnTriangle( ::getPosition( aabb ), tri );
-		auto closestAabb = ::closestPointOnAABB( closest, aabb.bounds );
-
-		if ( !uf::vector::isValid( closest ) ) return false;
-		
-		// to-do: derive proper delta
-		auto delta = closestAabb - closest;
-		float dist2 = uf::vector::dot( delta, delta );
-		float tolerance = 1.0e-3;
-		if ( dist2 >= tolerance ) return false;
-		float dist = std::sqrt( dist2 );
-
-		// to-do: properly derive the contact information
-		auto contact = closest; // ( closest + closestAabb ) * 0.5f;
-		auto normal = ( dist > eps ) ? ( delta / dist ) : ::triangleNormal( tri );
-		float penetration = tolerance - dist;
-
-	#if REORIENT_NORMALS_ON_CONTACT
-		if ( uf::vector::dot( normal, delta ) < 0.0f ) normal = -normal;
-	#endif
-
-		manifold.points.emplace_back(pod::Contact{ contact, normal, penetration });
-		return true;
-	#endif
+		return ::triangleAabbSAT( tri, body, manifold, eps );
+		//return ::triangleAabbTri( tri, body, manifold, eps );
 	}
 	bool triangleSphere( const pod::TriangleWithNormal& tri, const pod::PhysicsBody& body, pod::Manifold& manifold, float eps ) {
 		const auto& sphere = body;
@@ -470,7 +454,7 @@ namespace {
 
 		// to-do: derive proper delta
 		auto delta = center - closest;
-		float dist2 = uf::vector::dot( delta, delta );
+		float dist2 = uf::vector::magnitude( delta );
 		if ( dist2 > r * r ) return false;
 		float dist = std::sqrt(dist2);