Update the metric

Cargo.toml

@@ -16,3 +16,6 @@ glam = "0.27.0"
 flo_draw = "0.3.1"
 flo_canvas = "0.3.1"
 itertools-num = "0.1.3"
+
+[dev-dependencies]
+approx = "0.5.1"

src/bin/flat.rs

@@ -3,6 +3,9 @@ use flo_canvas::*;
 use glam::*;
 use riemann::{Decomp2, Metric, trace_iter};
 
+#[cfg(test)]
+use approx::assert_abs_diff_eq;
+
 pub fn main() {
 	let space = Coil {
 		scale: 3.0,
@@ -11,9 +14,10 @@ pub fn main() {
 		m: 10.0,
 	};
 	let space = Rect {
-		scale: 3.0,
+		inner_radius: 30.0,
-		r: vec2(30.0, 300.0),
+		outer_radius: 50.0,
-		m: vec2(20.0, 50.0),
+		internal_halflength: 100.0,
+		external_halflength: 300.0,
 	};
 	with_2d_graphics(move || {
 		let canvas = create_drawing_window("Refraction");
@@ -70,8 +74,8 @@ impl Renderable for Coil {
 impl Renderable for Rect {
 	fn render(&self, gc: &mut Vec<Draw>) {
 		gc.new_path();
-		gc.rect(-self.r.x - self.m.x, -self.r.y - self.m.y, self.r.x + self.m.x, self.r.y + self.m.y);
+		gc.rect(-self.outer_radius, -self.external_halflength, self.outer_radius, self.external_halflength);
-		gc.rect(-self.r.x, -self.r.y, self.r.x, self.r.y);
+		gc.rect(-self.inner_radius, -self.external_halflength, self.inner_radius, self.external_halflength);
 		gc.winding_rule(WindingRule::EvenOdd);
 		gc.fill_color(Color::Rgba(0.8, 0.8, 0.8, 1.0));
 		gc.fill();
@@ -79,7 +83,7 @@ impl Renderable for Rect {
 		gc.stroke_color(Color::Rgba(1.0, 0.5, 0.0, 1.0));
 		draw_fan(gc, self, vec2(-500.0, 0.0), vec2(1.0, 0.0), 1.0);
 		gc.stroke_color(Color::Rgba(0.0, 0.5, 1.0, 1.0));
-		draw_fan(gc, self, vec2(0.0, -0.5 * self.r.y), vec2(1.0, 1.0), 1.0);
+		draw_fan(gc, self, vec2(0.0, -0.5 * self.internal_halflength), vec2(1.0, 1.0), 1.0);
 	}
 }
 
@@ -109,17 +113,77 @@ impl Metric for Coil {
 }
 
 struct Rect {
-	scale: f32,
+	outer_radius: f32,
-	r: Vec2,
+	inner_radius: f32,
-	m: Vec2,
+	external_halflength: f32,
+	internal_halflength: f32,
+}
+
+impl Rect {
+	fn γ(&self) -> f32 { self.external_halflength / self.internal_halflength }
+	fn ri(&self) -> f32 { self.internal_halflength }
+	fn re(&self) -> f32 { self.external_halflength }
+	fn a(&self) -> f32 { (1.0 - self.γ()) / self.ri() }
+	fn b(&self) -> f32 { 2.0 * self.γ() - 1.0 }
+
+	fn root(&self, x: f32) -> f32 { ((2.0 * self.γ() - 1.0).powi(2) + 4.0 * (1.0 - self.γ()) * x / self.ri()).sqrt() }
+	fn d(&self, u: f32) -> f32 { 2.0 * self.a() * u + self.b() }
+	pub fn x(&self, u: f32) -> f32 { (self.a() * u.abs() + self.b()) * u }
+	pub fn u(&self, x: f32) -> f32 { 0.5 * self.ri() * (1.0 - 2.0 * self.γ() + self.root(x.abs())) / (1.0 - self.γ()) * x.signum() }
+	pub fn dx(&self, u: f32, du: f32) -> f32 { du * self.d(u.abs()) }
+	pub fn du(&self, x: f32, dx: f32) -> f32 { dx / self.root(x.abs()) }
+}
+
+#[test]
+fn test_rect() {
+	let r = Rect {
+		outer_radius: 50.0,
+		inner_radius: 20.0,
+		external_halflength: 100.0,
+		internal_halflength: 10.0,
+	};
+	assert_abs_diff_eq!(r.x(r.internal_halflength), r.external_halflength, epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.x(-r.internal_halflength), -r.external_halflength, epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.u(r.external_halflength), r.internal_halflength, epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.u(-r.external_halflength), -r.internal_halflength, epsilon = 1.0e-5);
+
+	assert_abs_diff_eq!(r.dx(r.internal_halflength, 3.0), 3.0, epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.dx(-r.internal_halflength, 3.0), 3.0, epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.du(r.external_halflength, 3.0), 3.0, epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.du(-r.external_halflength, 3.0), 3.0, epsilon = 1.0e-5);
+
+	assert_abs_diff_eq!(r.u(r.x(1.0)), 1.0, epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.u(r.x(5.0)), 5.0, epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.u(r.x(-5.0)), -5.0, epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.u(r.x(10.0)), 10.0, epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.x(r.u(10.0)), 10.0, epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.x(r.u(50.0)), 50.0, epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.x(r.u(-50.0)), -50.0, epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.x(r.u(100.0)), 100.0, epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.d(r.u(10.0)), r.root(10.0), epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.d(r.u(50.0)), r.root(50.0), epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.d(r.u(100.0)), r.root(100.0), epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.du(10.0, r.dx(r.u(10.0), 3.0)), 3.0, epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.du(50.0, r.dx(r.u(50.0), 3.0)), 3.0, epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.du(-50.0, r.dx(r.u(-50.0), 3.0)), 3.0, epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.du(100.0, r.dx(r.u(100.0), 3.0)), 3.0, epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.dx(1.0, r.du(r.x(1.0), 3.0)), 3.0, epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.dx(5.0, r.du(r.x(5.0), 3.0)), 3.0, epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.dx(-5.0, r.du(r.x(-5.0), 3.0)), 3.0, epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.dx(10.0, r.du(r.x(10.0), 3.0)), 3.0, epsilon = 1.0e-5);
 }
 
 impl Metric for Rect {
 	fn halfmetric(&self, pos: Vec2) -> Decomp2 {
-		let s = smoothbox(pos.x, vec2(-self.r.x, self.r.x), self.m.x) * smoothbox(pos.y, vec2(-self.r.y, self.r.y), self.m.y);
+		let x = pos.y.abs();
+		let sx = ((pos.x.abs() - self.inner_radius) / (self.outer_radius - self.inner_radius)).clamp(0.0, 1.0);
+		let sy = if x <= self.external_halflength { 1.0 / self.root(x) } else { 1.0 };
+		assert!(sx.is_finite());
+		assert!(sy.is_finite());
+		assert!(sy > 0.0);
 		Decomp2 {
 			ortho: Mat2::IDENTITY,
-			diag: vec2(1.0, self.scale.powf(-s)),
+			diag: vec2(1.0, sy.lerp(1.0, sx)),
 		}
 	}
 }