diff --git a/src/bin/flat/fns.rs b/src/bin/flat/fns.rs index 9c84bb3..a4ed940 100644 --- a/src/bin/flat/fns.rs +++ b/src/bin/flat/fns.rs @@ -30,6 +30,26 @@ impl SmoothstepLimiter { } } +pub struct SmootherstepLimiter { + pub min: f32, + pub max: f32, +} + +impl SmootherstepLimiter { + pub fn value(&self, x: f32) -> f32 { + let y = (self.min, self.max).inverse_lerp(x.abs()).clamp(0.0, 1.0); + 6.0 * y.powi(5) - 15.0 * y.powi(4) + 10.0 * y.powi(3) + } + 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()); + 30.0 * (t * (1.0 - t)).powi(2) * x.signum() / (self.max - self.min) + } else { + 0.0 + } + } +} + pub struct QuadraticAccelerator { pub internal: f32, pub external: f32, @@ -113,6 +133,32 @@ mod test { } } + #[test] + fn test_smootherstep_limiter() { + let testee = super::SmootherstepLimiter { min: 20.0, max: 30.0 }; + let ε = 1.0e-4f32; + let δ = 1.0 / 32.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 78111fb..689ace5 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::SmoothstepLimiter { fns::SmoothstepLimiter { min: self.inner_radius, max: self.outer_radius } } + fn fx(&self) -> fns::SmootherstepLimiter { fns::SmootherstepLimiter { 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) } @@ -156,7 +156,7 @@ mod test { }, objs: vec![], }; - let ε = 1e-5; + let ε = 1e-3; let off = 10.0; let steps = 10000; for x in [space.tube.inner_radius - ε, space.tube.inner_radius + ε] { @@ -166,7 +166,7 @@ mod test { let dir = Δ / (steps as f32); let traced = space.trace_iter(Ray { pos: a, dir }).nth(steps).unwrap(); println!("{traced:?}"); - assert_abs_diff_eq!(traced.pos.x, b.x, epsilon=1.0e1); + assert_abs_diff_eq!(traced.pos.x, b.x, epsilon=1.0e-1); 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);