From 3fb8bf42b76f44b4faae1cf906f00c1268c7c414 Mon Sep 17 00:00:00 2001
From: numzero <numzer0@yandex.ru>
Date: Fri, 28 Jun 2024 23:50:01 +0300
Subject: [PATCH] Use smmothstep to smoothen the metric

---
 src/bin/flat/fns.rs         | 46 +++++++++++++++++++++++++++++++++++++
 src/bin/flat/tube/metric.rs | 10 ++++----
 2 files changed, 51 insertions(+), 5 deletions(-)

diff --git a/src/bin/flat/fns.rs b/src/bin/flat/fns.rs
index 526e6d7..9c84bb3 100644
--- a/src/bin/flat/fns.rs
+++ b/src/bin/flat/fns.rs
@@ -10,6 +10,26 @@ impl LinearLimiter {
 	pub fn derivative(&self, x: f32) -> f32 { if x.abs() > self.min && x.abs() < self.max { x.signum() / (self.max - self.min) } else { 0.0 } }
 }
 
+pub struct SmoothstepLimiter {
+	pub min: f32,
+	pub max: f32,
+}
+
+impl SmoothstepLimiter {
+	pub fn value(&self, x: f32) -> f32 {
+		let y = (self.min, self.max).inverse_lerp(x.abs()).clamp(0.0, 1.0);
+		3.0 * y * y - 2.0 * y * y * y
+	}
+	pub fn derivative(&self, x: f32) -> f32 {
+		if x.abs() > self.min && x.abs() < self.max {
+			let t = (self.min, self.max).inverse_lerp(x.abs());
+			6.0 * x.signum() * t * (1.0 - t) / (self.max - self.min)
+		} else {
+			0.0
+		}
+	}
+}
+
 pub struct QuadraticAccelerator {
 	pub internal: f32,
 	pub external: f32,
@@ -67,6 +87,32 @@ mod test {
 		}
 	}
 
+	#[test]
+	fn test_smoothsteo_limiter() {
+		let testee = super::SmoothstepLimiter { min: 20.0, max: 30.0 };
+		let ε = 1.0e-4f32;
+		let δ = 1.0 / 8.0; // Mathematically, you want this to be small. Computationally, you don’t.
+		let margin = 1.0 / 16.0;
+		let mul = 1.0 + margin;
+		for x in itertools_num::linspace(0., testee.min, 10) {
+			assert_abs_diff_eq!(testee.value(x), 0., epsilon = ε);
+			assert_abs_diff_eq!(testee.value(-x), 0., epsilon = ε);
+		}
+		for x in itertools_num::linspace(testee.max, mul * testee.max, 10) {
+			assert_abs_diff_eq!(testee.value(x), 1., epsilon = ε);
+			assert_abs_diff_eq!(testee.value(-x), 1., epsilon = ε);
+		}
+		for x in itertools_num::linspace(-mul * testee.max, mul * testee.max, 100) {
+			// Currently, the derivative is discontinuous at ±min and ±max... let’s just skip these for now.
+			if x.abs().abs_diff_eq(&testee.min, δ) || x.abs().abs_diff_eq(&testee.max, δ) {
+				continue;
+			}
+			let df_num = (testee.value(x + δ) - testee.value(x - δ)) / (2. * δ);
+			let df_expl = testee.derivative(x);
+			assert!(abs_diff_eq!(df_expl, df_num, epsilon = ε), "At x={x}, df/dx:\nnumerical: {df_num}\nexplicit:  {df_expl}\n");
+		}
+	}
+
 	#[test]
 	fn test_quadratic_accelerator() {
 		let testee = super::QuadraticAccelerator { internal: 100.0, external: 150.0 };
diff --git a/src/bin/flat/tube/metric.rs b/src/bin/flat/tube/metric.rs
index 96c6fe2..78111fb 100644
--- a/src/bin/flat/tube/metric.rs
+++ b/src/bin/flat/tube/metric.rs
@@ -11,7 +11,7 @@ pub struct Tube {
 }
 
 impl Tube {
-	fn fx(&self) -> fns::LinearLimiter { fns::LinearLimiter { min: self.inner_radius, max: self.outer_radius } }
+	fn fx(&self) -> fns::SmoothstepLimiter { fns::SmoothstepLimiter { min: self.inner_radius, max: self.outer_radius } }
 	fn fy(&self) -> fns::QuadraticAccelerator { fns::QuadraticAccelerator { internal: self.internal_halflength, external: self.external_halflength } }
 
 	pub fn y(&self, v: f32) -> f32 { self.fy().x(v) }
@@ -146,7 +146,6 @@ mod test {
 	}
 
 	#[test]
-	#[ignore] // TODO make the metric smooth so that this passes
 	fn test_accelerator_inner_edge() {
 		let space = Space {
 			tube: Tube {
@@ -157,16 +156,17 @@ mod test {
 			},
 			objs: vec![],
 		};
-		let ε = 1e-3;
+		let ε = 1e-5;
 		let off = 10.0;
-		let steps = 4096;
+		let steps = 10000;
 		for x in [space.tube.inner_radius - ε, space.tube.inner_radius + ε] {
 			let a = vec2(x, -(space.tube.external_halflength + off));
 			let b = vec2(x, space.tube.external_halflength + off);
 			let Δ = vec2(0.0, 2.0 * (space.tube.internal_halflength + off));
 			let dir = Δ / (steps as f32);
 			let traced = space.trace_iter(Ray { pos: a, dir }).nth(steps).unwrap();
-			assert_abs_diff_eq!(traced.pos.x, b.x, epsilon=1.0e0);
+			println!("{traced:?}");
+			assert_abs_diff_eq!(traced.pos.x, b.x, epsilon=1.0e1);
 			assert_abs_diff_eq!(traced.pos.y, b.y, epsilon=1.0e0);
 			assert_abs_diff_eq!(traced.dir.x, dir.x, epsilon=1.0e-2);
 			assert_abs_diff_eq!(traced.dir.y, dir.y, epsilon=1.0e-2);