Cargo.toml

@@ -16,6 +16,3 @@ 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,12 +5,6 @@ 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,
@@ -19,26 +13,15 @@ pub fn main() {
 		m: 10.0,
 	};
 	let space = Rect {
-		inner_radius: 30.0,
+		scale: 3.0,
-		outer_radius: 50.0,
+		r: vec2(30.0, 300.0),
-		internal_halflength: 100.0,
+		m: vec2(20.0, 50.0),
-		external_halflength: 300.0,
 	};
 	with_2d_graphics(move || {
 		let canvas = create_drawing_window("Refraction");
 		canvas.draw(|gc| {
 			gc.canvas_height(1000.0);
-			gc.line_width(0.5);
+			space.render(gc);
-
-			// 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![
@@ -46,16 +29,40 @@ 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);
-		});
-	});
-}
 
-trait Interesting {
+			gc.new_path();
-	fn draw_background(&self, gc: &mut Vec<Draw>);
+			let dt = 1.0;
-	fn draw_rays(&self, gc: &mut Vec<Draw>, tracer: &(impl Rayable + ?Sized));
+			let mut s = boundary::Id(0);
+			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 {}
@@ -93,75 +100,37 @@ impl SpaceVisual for Inside {
 	}
 
 	fn globalize_loc(&self, pos: Vec2) -> Vec2 {
-		vec2(pos.x, self.space.x(pos.y))
+		vec2(pos.x, pos.y * self.space.scale)
 	}
 }
 
-impl Rayable for [Box<dyn SpacePart>] {
+fn draw_ray(gc: &mut Vec<Draw>, space: &impl Metric, base: Vec2, dir: Vec2) {
-	fn draw_ray(&self, gc: &mut Vec<Draw>, base: Vec2, dir: Vec2) {
+	let dir = space.globalize(base, dir);
-		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(self, base, dir, DT).take(10000) {
+	for pt in trace_iter(space, base, dir, 1.0).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);
 	}
 }
 
-impl Interesting for Coil {
+trait Renderable {
-	fn draw_background(&self, gc: &mut Vec<Draw>) {
+	fn render(&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);
@@ -170,35 +139,27 @@ impl Interesting 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));
-		tracer.draw_fan(gc, vec2(-500.0, 0.0), vec2(1.0, 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));
-		tracer.draw_fan(gc, vec2(0.0, self.r), vec2(1.0, 0.0), 1.0);
+		draw_fan(gc, self, vec2(0.0, self.r), vec2(1.0, 0.0), 1.0);
 	}
 }
 
-impl Interesting for Rect {
+impl Renderable for Rect {
-	fn draw_background(&self, gc: &mut Vec<Draw>) {
+	fn render(&self, gc: &mut Vec<Draw>) {
 		gc.new_path();
-		gc.rect(-self.outer_radius, -self.external_halflength, self.outer_radius, self.external_halflength);
+		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.inner_radius, -self.external_halflength, self.inner_radius, self.external_halflength);
+		gc.rect(-self.r.x, -self.r.y, self.r.x, self.r.y);
 		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));
-		tracer.draw_fan(gc, vec2(-500.0, 0.0), vec2(1.0, 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.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));
-		tracer.draw_fan(gc, vec2(0.0, -0.5 * self.external_halflength), vec2(1.0, 1.0), 1.0);
+		draw_fan(gc, self, vec2(0.0, -0.5 * self.r.y), 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);
 	}
 }
 
@@ -228,77 +189,17 @@ impl Metric for Coil {
 }
 
 struct Rect {
-	outer_radius: f32,
+	scale: f32,
-	inner_radius: f32,
+	r: Vec2,
-	external_halflength: f32,
+	m: Vec2,
-	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 x = pos.y.abs();
+		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 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, sy.lerp(1.0, sx)),
+			diag: vec2(1.0, self.scale.powf(-s)),
 		}
 	}
 }
@@ -328,20 +229,11 @@ struct Inside {
 
 impl boundary::Boundary for Outside {
 	fn next(&self, base: Vec2, dir: Vec2, limit: f32) -> Option<(boundary::Id, Vec2, Vec2)> {
-		let or = self.space.outer_radius;
+		let size = self.space.r + self.space.m;
-		let ir = self.space.inner_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 hl = self.space.external_halflength;
+		let (_, dist) = bnd.hit(base, dir)?;
-		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 {
-			let pt = base + dist * dir;
+			return Some((boundary::Id(1), base + dist * dir, 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
 	}
@@ -355,17 +247,16 @@ impl Metric for Outside {
 
 impl boundary::Boundary for Wall {
 	fn next(&self, base: Vec2, dir: Vec2, limit: f32) -> Option<(boundary::Id, Vec2, Vec2)> {
-		let or = self.space.outer_radius;
+		let osize = self.space.r + self.space.m;
-		let ir = self.space.inner_radius;
+		let isize = self.space.r;
-		let hl = self.space.external_halflength;
+		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 obnd = Loop(vec![vec2(-or, -hl), vec2(-or, hl), vec2(or, hl), vec2(or, -hl)]);
+		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 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 v = vec2(v.x, self.space.du(p.y, v.y));
+				let p = vec2(p.x, p.y / self.space.scale);
-				let p = vec2(p.x, self.space.u(p.y));
+				let v = vec2(v.x, v.y / self.space.scale);
 				return Some((boundary::Id(2), p, v));
 			}
 		}
@@ -386,19 +277,16 @@ 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 = vec2(self.space.inner_radius, self.space.internal_halflength);
+		let size = self.space.r;
+		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 (side, dist) = bnd.hit(base, dir)?;
+		let (_, dist) = bnd.hit(base, dir)?;
 		if dist <= limit {
 			let p = base + dist * dir;
 			let v = dir;
-			let v = vec2(v.x, self.space.dx(p.y, v.y));
+			let p = vec2(p.x, p.y * self.space.scale);
-			let p = vec2(p.x, self.space.x(p.y));
+			let v = vec2(v.x, v.y * self.space.scale);
-			return match side {
+			return Some((boundary::Id(1), p, v));
-				0 | 2 => Some((boundary::Id(1), p, v)),
-				1 | 3 => Some((boundary::Id(0), p, v)),
-				_ => panic!()
-			};
 		}
 		None
 	}
@@ -602,15 +490,10 @@ 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], linear in-between.
+/// 1.0 for val∈[range.x, range.y], 0.0 for val∉[range.x−pad, range.y+pad], smoothstep in-between.
-fn trapezoid(val: f32, range: Vec2, pad: f32) -> f32 {
+fn smoothbox(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);
-	lin.clamp(0.0, 1.0)
+	smoothstep(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))
 }