Some concerns ain’t easily separated

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

@@ -5,6 +5,12 @@ use glam::*;
 use riemann::{Decomp2, Metric, trace_iter};
 use crate::boundary::Loop;
 
+#[cfg(test)]
+use approx::assert_abs_diff_eq;
+
+const RAYS_IN_FAN: usize = 101;
+const DT: f32 = 0.1;
+
 pub fn main() {
 	let space = Coil {
 		scale: 3.0,
@@ -13,15 +19,26 @@ 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");
 		canvas.draw(|gc| {
 			gc.canvas_height(1000.0);
-			space.render(gc);
+			gc.line_width(0.5);
+
+			// gc.canvas_height(100.0);
+			// gc.center_region(-100.0, -200.0, 100.0, 0.0);
+			// gc.line_width(0.1);
+
+			space.draw_background(gc);
+
+			// gc.stroke_color(Color::Rgba(0.5, 1.0, 0.0, 1.0));
+			// space.draw_fan(gc, vec2(0.0, -0.5 * space.external_halflength), vec2(1.0, 1.0), 1.0);
+			space.draw_rays(gc, &space);
 
 			let space = Rc::new(space);
 			let parts: Vec<Box<dyn SpacePart>> = vec![
@@ -29,40 +46,16 @@ pub fn main() {
 				Box::new(Wall { space: space.clone() }),
 				Box::new(Inside { space: space.clone() }),
 			];
+			space.draw_rays(gc, parts.as_slice());
+			// gc.stroke_color(Color::Rgba(1.0, 0.0, 1.0, 1.0));
+			// parts.draw_fan(gc, vec2(0.0, space.u(-0.5 * space.external_halflength)), vec2(1.0, 1.0), 1.0);
+		});
+	});
+}
 
-			gc.new_path();
+trait Interesting {
-			let dt = 1.0;
+	fn draw_background(&self, gc: &mut Vec<Draw>);
-			let mut s = boundary::Id(0);
+	fn draw_rays(&self, gc: &mut Vec<Draw>, tracer: &(impl Rayable + ?Sized));
-			let mut p = vec2(-500.0, 0.0);
-			let mut v = vec2(3.0, 1.0).normalize();
-			let part = &*parts[s.0 as usize];
-			gc.stroke_color(part.color());
-			gc.move_to(p.x, p.y);
-			for _ in 0..10000 {
-				let part = &*parts[s.0 as usize];
-				let a: Vec2 = -riemann::convolute(riemann::krist(part, p), v);
-				v = v + a * dt;
-				if let Some((id, base, dir)) = part.next(p, v, dt) {
-					gc.stroke();
-					gc.new_path();
-					let pt = part.globalize_loc(p);
-					gc.move_to(pt.x, pt.y);
-					s = id;
-					p = base;
-					v = dir;
-					let part = &*parts[s.0 as usize];
-					let pt = part.globalize_loc(p);
-					gc.stroke_color(part.color());
-					gc.line_to(pt.x, pt.y);
-				} else {
-					p = p + v * dt;
-					let pt = part.globalize_loc(p);
-					gc.line_to(pt.x, pt.y);
-				}
-			}
-			gc.stroke();
-		});
-	});
 }
 
 trait SpacePart: boundary::Boundary + Metric + SpaceVisual {}
@@ -100,37 +93,75 @@ impl SpaceVisual for Inside {
 	}
 
 	fn globalize_loc(&self, pos: Vec2) -> Vec2 {
-		vec2(pos.x, pos.y * self.space.scale)
+		vec2(pos.x, self.space.x(pos.y))
 	}
 }
 
-fn draw_ray(gc: &mut Vec<Draw>, space: &impl Metric, base: Vec2, dir: Vec2) {
+impl Rayable for [Box<dyn SpacePart>] {
-	let dir = space.globalize(base, dir);
+	fn draw_ray(&self, gc: &mut Vec<Draw>, base: Vec2, dir: Vec2) {
+		gc.new_path();
+		let dt = DT;
+		let mut s = boundary::Id(0);
+		let mut p = base;
+		let mut v = dir.normalize();
+		let part = &*self[s.0 as usize];
+		gc.stroke_color(part.color());
+		gc.move_to(p.x, p.y);
+		for _ in 0..10000 {
+			let part = &*self[s.0 as usize];
+			let a: Vec2 = -riemann::convolute(riemann::krist(part, p), v);
+			v = v + a * dt;
+			if let Some((id, base, dir)) = part.next(p, v, dt) {
+				gc.stroke();
+				gc.new_path();
+				let pt = part.globalize_loc(p);
+				gc.move_to(pt.x, pt.y);
+				s = id;
+				p = base;
+				v = dir;
+				let part = &*self[s.0 as usize];
+				let pt = part.globalize_loc(p);
+				gc.stroke_color(part.color());
+				gc.line_to(pt.x, pt.y);
+			} else {
+				p = p + v * dt;
+				let pt = part.globalize_loc(p);
+				gc.line_to(pt.x, pt.y);
+			}
+		}
+		gc.stroke();
+	}
+}
+
+trait Rayable {
+	fn draw_ray(&self, gc: &mut Vec<Draw>, base: Vec2, dir: Vec2);
+
+	fn draw_fan(&self, gc: &mut Vec<Draw>, base: Vec2, dir: Vec2, spread: f32) {
+		let dir = dir.normalize();
+		let v = vec2(-dir.y, dir.x);
+		for y in itertools_num::linspace(-spread, spread, RAYS_IN_FAN) {
+			self.draw_ray(gc, base, dir + y * v);
+		}
+	}
+}
+
+impl<T: Metric> Rayable for T {
+	fn draw_ray(&self, gc: &mut Vec<Draw>, base: Vec2, dir: Vec2) {
+		let dir = self.globalize(base, dir);
 		gc.new_path();
 		gc.move_to(base.x, base.y);
-	for pt in trace_iter(space, base, dir, 1.0).take(10000) {
+		for pt in trace_iter(self, base, dir, DT).take(10000) {
 			gc.line_to(pt.x, pt.y);
 			if pt.abs().cmpgt(Vec2::splat(1000.0)).any() {
 				break;
 			}
 		}
 		gc.stroke();
-}
-
-fn draw_fan(gc: &mut Vec<Draw>, space: &impl Metric, base: Vec2, dir: Vec2, spread: f32) {
-	let dir = dir.normalize();
-	let v = vec2(-dir.y, dir.x);
-	for y in itertools_num::linspace(-spread, spread, 101) {
-		draw_ray(gc, space, base, dir + y * v);
 	}
 }
 
-trait Renderable {
+impl Interesting for Coil {
-	fn render(&self, gc: &mut Vec<Draw>);
+	fn draw_background(&self, gc: &mut Vec<Draw>) {
-}
-
-impl Renderable for Coil {
-	fn render(&self, gc: &mut Vec<Draw>) {
 		gc.new_path();
 		gc.circle(0.0, 0.0, self.r + self.w + self.m);
 		gc.circle(0.0, 0.0, self.r + self.w);
@@ -139,27 +170,35 @@ impl Renderable for Coil {
 		gc.winding_rule(WindingRule::EvenOdd);
 		gc.fill_color(Color::Rgba(0.8, 0.8, 0.8, 1.0));
 		gc.fill();
-		gc.line_width(0.5);
+	}
+
+	fn draw_rays(&self, gc: &mut Vec<Draw>, tracer: &(impl Rayable + ?Sized)) {
 		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);
+		tracer.draw_fan(gc, 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, self.r), vec2(1.0, 0.0), 1.0);
+		tracer.draw_fan(gc, vec2(0.0, self.r), vec2(1.0, 0.0), 1.0);
 	}
 }
 
-impl Renderable for Rect {
+impl Interesting for Rect {
-	fn render(&self, gc: &mut Vec<Draw>) {
+	fn draw_background(&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();
-		gc.line_width(0.5);
+	}
+
+	fn draw_rays(&self, gc: &mut Vec<Draw>, tracer: &(impl Rayable + ?Sized)) {
 		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);
+		tracer.draw_fan(gc, vec2(-500.0, 0.0), vec2(1.0, 0.0), 1.0);
+		gc.stroke_color(Color::Rgba(0.5, 1.0, 0.0, 1.0));
+		tracer.draw_fan(gc, vec2(0.0, -0.5 * self.external_halflength), vec2(1.0, 1.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);
+		tracer.draw_fan(gc, vec2(0.0, -0.5 * self.external_halflength), vec2(1.0, 1.0), 1.0);
+		gc.stroke_color(Color::Rgba(0.2, 0.7, 0.0, 1.0));
+		tracer.draw_fan(gc, vec2(-0.5 * self.inner_radius, -1.2 * self.external_halflength), vec2(0.0, 1.0), 1.0);
 	}
 }
 
@@ -189,17 +228,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)),
 		}
 	}
 }
@@ -229,11 +328,20 @@ struct Inside {
 
 impl boundary::Boundary for Outside {
 	fn next(&self, base: Vec2, dir: Vec2, limit: f32) -> Option<(boundary::Id, Vec2, Vec2)> {
-		let size = self.space.r + self.space.m;
+		let or = self.space.outer_radius;
-		let bnd = Loop(vec![vec2(-size.x, -size.y), vec2(size.x, -size.y), vec2(size.x, size.y), vec2(-size.x, size.y)]);
+		let ir = self.space.inner_radius;
-		let (_, dist) = bnd.hit(base, dir)?;
+		let hl = self.space.external_halflength;
+		let bnd = Loop(vec![
+			vec2(-or, -hl), vec2(-ir, -hl), vec2(ir, -hl), vec2(or, -hl),
+			vec2(or, hl), vec2(ir, hl), vec2(-ir, hl), vec2(-or, hl),
+		]);
+		let (side, dist) = bnd.hit(base, dir)?;
 		if dist <= limit {
-			return Some((boundary::Id(1), base + dist * dir, dir));
+			let pt = base + dist * dir;
+			return match side {
+				1 | 5 => Some((boundary::Id(2), vec2(pt.x, self.space.u(pt.y)), vec2(dir.x, self.space.du(pt.y, dir.y)))),
+				_ => Some((boundary::Id(1), pt, dir)),
+			};
 		}
 		None
 	}
@@ -247,16 +355,17 @@ impl Metric for Outside {
 
 impl boundary::Boundary for Wall {
 	fn next(&self, base: Vec2, dir: Vec2, limit: f32) -> Option<(boundary::Id, Vec2, Vec2)> {
-		let osize = self.space.r + self.space.m;
+		let or = self.space.outer_radius;
-		let isize = self.space.r;
+		let ir = self.space.inner_radius;
-		let obnd = Loop(vec![vec2(-osize.x, -osize.y), vec2(-osize.x, osize.y), vec2(osize.x, osize.y), vec2(osize.x, -osize.y)]);
+		let hl = self.space.external_halflength;
-		let ibnd = Loop(vec![vec2(-isize.x, -isize.y), vec2(isize.x, -isize.y), vec2(isize.x, isize.y), vec2(-isize.x, isize.y)]);
+		let obnd = Loop(vec![vec2(-or, -hl), vec2(-or, hl), vec2(or, hl), vec2(or, -hl)]);
+		let ibnd = Loop(vec![vec2(-ir, -hl), vec2(ir, -hl), vec2(ir, hl), vec2(-ir, hl)]);
 		if let Some((_, dist)) = ibnd.hit(base, dir) {
 			if dist <= limit {
 				let p = base + dist * dir;
 				let v = dir;
-				let p = vec2(p.x, p.y / self.space.scale);
+				let v = vec2(v.x, self.space.du(p.y, v.y));
-				let v = vec2(v.x, v.y / self.space.scale);
+				let p = vec2(p.x, self.space.u(p.y));
 				return Some((boundary::Id(2), p, v));
 			}
 		}
@@ -277,16 +386,19 @@ impl Metric for Wall {
 
 impl boundary::Boundary for Inside {
 	fn next(&self, base: Vec2, dir: Vec2, limit: f32) -> Option<(boundary::Id, Vec2, Vec2)> {
-		let size = self.space.r;
+		let size = vec2(self.space.inner_radius, self.space.internal_halflength);
-		let size = vec2(size.x, size.y / self.space.scale);
 		let bnd = Loop(vec![vec2(-size.x, -size.y), vec2(-size.x, size.y), vec2(size.x, size.y), vec2(size.x, -size.y)]);
-		let (_, dist) = bnd.hit(base, dir)?;
+		let (side, dist) = bnd.hit(base, dir)?;
 		if dist <= limit {
 			let p = base + dist * dir;
 			let v = dir;
-			let p = vec2(p.x, p.y * self.space.scale);
+			let v = vec2(v.x, self.space.dx(p.y, v.y));
-			let v = vec2(v.x, v.y * self.space.scale);
+			let p = vec2(p.x, self.space.x(p.y));
-			return Some((boundary::Id(1), p, v));
+			return match side {
+				0 | 2 => Some((boundary::Id(1), p, v)),
+				1 | 3 => Some((boundary::Id(0), p, v)),
+				_ => panic!()
+			};
 		}
 		None
 	}
@@ -490,10 +602,15 @@ fn smoothstep(x: f32) -> f32 {
 	3.0 * x * x - 2.0 * x * x * x
 }
 
-/// 1.0 for val∈[range.x, range.y], 0.0 for val∉[range.x−pad, range.y+pad], smoothstep in-between.
+/// 1.0 for val∈[range.x, range.y], 0.0 for val∉[range.x−pad, range.y+pad], linear in-between.
-fn smoothbox(val: f32, range: Vec2, pad: f32) -> f32 {
+fn trapezoid(val: f32, range: Vec2, pad: f32) -> f32 {
 	let slope1 = 1.0 + (val - range.x) / pad;
 	let slope2 = 1.0 - (val - range.y) / pad;
 	let lin = slope1.min(slope2);
-	smoothstep(lin.clamp(0.0, 1.0))
+	lin.clamp(0.0, 1.0)
+}
+
+/// 1.0 for val∈[range.x, range.y], 0.0 for val∉[range.x−pad, range.y+pad], smoothstep in-between.
+fn smoothbox(val: f32, range: Vec2, pad: f32) -> f32 {
+	smoothstep(trapezoid(val, range, pad))
 }