Fixup the Cargo.lock

Apparently it wants more derivatives than are explicitly used

Cargo.toml

@@ -1,14 +1,25 @@
 [package]
-name = "hello"
+name = "refraction"
 version = "0.1.0"
 edition = "2021"
 
 # See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
 
 [profile.dev]
+panic = 'abort'
+
+[profile.dev.package."*"]
 opt-level = 3
 
 [dependencies]
 rand = "0.8.5"
-glm = "0.2.3"
+glam = { version = "0.27.0", features = ["approx", "fast-math", "rand"] }
 show-image = "0.14.0"
+flo_draw = "0.3.1"
+flo_canvas = "0.3.1"
+itertools-num = "0.1.3"
+
+[dev-dependencies]
+approx = "0.5.1"
+rand = "0.8.5"
+rand_pcg = "0.3.1"

src/bin/flat/float_fun.rs

@@ -0,0 +1,36 @@
+use glam::FloatExt;
+
+mod bounds {
+	pub trait Pair<T> {}
+
+	impl<T> Pair<T> for (T, T) {}
+}
+
+pub trait FloatExt2<T>: bounds::Pair<T> {
+	fn lerp(self, t: T) -> T;
+	fn inverse_lerp(self, y: T) -> T;
+}
+
+impl<F: FloatExt> FloatExt2<F> for (F, F) {
+	fn lerp(self, t: F) -> F { F::lerp(self.0, self.1, t) }
+	fn inverse_lerp(self, y: F) -> F { F::inverse_lerp(self.0, self.1, y) }
+}
+
+#[cfg(test)]
+mod test {
+	use super::FloatExt2;
+
+	#[test]
+	fn test_lerp() {
+		assert_eq!((3., 7.).lerp(-0.5), 1.);
+		assert_eq!((3., 7.).lerp(0.0), 3.);
+		assert_eq!((3., 7.).lerp(0.5), 5.);
+		assert_eq!((3., 7.).lerp(1.0), 7.);
+		assert_eq!((3., 7.).lerp(1.5), 9.);
+		assert_eq!((3., 7.).inverse_lerp(1.), -0.5);
+		assert_eq!((3., 7.).inverse_lerp(3.), 0.0);
+		assert_eq!((3., 7.).inverse_lerp(5.), 0.5);
+		assert_eq!((3., 7.).inverse_lerp(7.), 1.0);
+		assert_eq!((3., 7.).inverse_lerp(9.), 1.5);
+	}
+}

src/bin/flat/fns.rs

@@ -0,0 +1,150 @@
+use crate::float_fun::FloatExt2;
+
+pub trait Limiter {
+	fn value(&self, x: f32) -> f32;
+	fn derivative(&self, x: f32) -> f32;
+}
+
+pub struct SmoothstepLimiter {
+	pub min: f32,
+	pub max: f32,
+}
+
+impl Limiter for SmoothstepLimiter {
+	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
+	}
+	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 SmootherstepLimiter {
+	pub min: f32,
+	pub max: f32,
+}
+
+impl Limiter for SmootherstepLimiter {
+	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)
+	}
+	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,
+}
+
+/// Продолжает функцию f с [-lim, lim] линейно в предположении f(±lim) = ±val, f'(±lim) = 1.
+fn extend_linear(t: f32, f: impl FnOnce(f32) -> f32, lim: f32, val: f32) -> f32 {
+	if t.abs() <= lim { f(t) } else { t + t.signum() * (val - lim) }
+}
+
+/// Продолжает функцию f с [-lim, lim] константой в предположении f(±lim) = val, f'(±lim) = 0.
+fn extend_const(t: f32, f: impl FnOnce(f32) -> f32, lim: f32, val: f32) -> f32 {
+	if t.abs() <= lim { f(t) } else { val }
+}
+
+impl QuadraticAccelerator {
+	fn a(&self) -> f32 { -(self.external - self.internal) / self.internal.powi(2) }
+	fn b(&self) -> f32 { 2.0 * self.external / self.internal - 1.0 }
+	fn root(&self, x: f32) -> f32 { (self.b().powi(2) + 4.0 * self.a() * x.abs()).sqrt() }
+
+	pub fn x(&self, u: f32) -> f32 { extend_linear(u, |u| (self.a() * u.abs() + self.b()) * u, self.internal, self.external) }
+	pub fn u(&self, x: f32) -> f32 { extend_linear(x, |x| 0.5 * x.signum() * (-self.b() + self.root(x)) / self.a(), self.external, self.internal) }
+	pub fn dx(&self, u: f32) -> f32 { extend_const(u, |u| 2.0 * self.a() * u.abs() + self.b(), self.internal, 1.0) }
+	pub fn du(&self, x: f32) -> f32 { extend_const(x, |x| 1.0 / self.root(x), self.external, 1.0) }
+	pub fn d2u(&self, x: f32) -> f32 { extend_const(x, |x| -2.0 * x.signum() * self.a() * self.root(x).powi(-3), self.external, 0.0) }
+}
+
+#[cfg(test)]
+mod test {
+	use approx::{abs_diff_eq, AbsDiffEq, assert_abs_diff_eq};
+
+	use super::*;
+
+	fn test_limiter(testee: impl Limiter, min: f32, max: f32, δ: f32) {
+		let ε = 1.0e-4f32;
+		let margin = 1.0 / 16.0;
+		let mul = 1.0 + margin;
+		for x in itertools_num::linspace(0., 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(max, mul * 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 * max, mul * max, 100) {
+			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_smoothstep_limiter() {
+		test_limiter(SmoothstepLimiter { min: 20.0, max: 30.0 }, 20.0, 30.0, 1.0 / 32.0);
+	}
+
+	#[test]
+	fn test_smootherstep_limiter() {
+		test_limiter(SmootherstepLimiter { min: 20.0, max: 30.0 }, 20.0, 30.0, 1.0 / 32.0);
+	}
+
+	#[test]
+	fn test_quadratic_accelerator() {
+		let testee = super::QuadraticAccelerator { internal: 100.0, external: 150.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;
+		assert_abs_diff_eq!(testee.u(testee.external), testee.internal, epsilon = ε);
+		assert_abs_diff_eq!(testee.u(-testee.external), -testee.internal, epsilon = ε);
+		assert_abs_diff_eq!(testee.du(testee.external), 1., epsilon = ε);
+		assert_abs_diff_eq!(testee.du(-testee.external), 1., epsilon = ε);
+		for x in itertools_num::linspace(-mul * testee.external, mul * testee.external, 100) {
+			let ux = testee.u(x);
+			let xux = testee.x(ux);
+			assert!(abs_diff_eq!(x, xux, epsilon = ε), "At x={x}:\nu(x):    {ux}\nx(u(x)): {xux}\n");
+
+			let du_num = (testee.u(x + δ) - testee.u(x - δ)) / (2. * δ);
+			let du_expl = testee.du(x);
+			assert!(abs_diff_eq!(du_expl, du_num, epsilon = ε), "At x={x}, du/dx:\nnumerical: {du_num}\nexplicit:  {du_expl}\n");
+
+			let dudx = du_expl * testee.dx(ux);
+			assert!(abs_diff_eq!(dudx, 1.0, epsilon = ε), "At x={x}:\ndu/dx * dx/du: {dudx}\n");
+
+			let d2u_num = (testee.du(x + δ) - testee.du(x - δ)) / (2. * δ);
+			let d2u_expl = testee.d2u(x);
+			assert!(abs_diff_eq!(d2u_expl, d2u_num, epsilon = ε), "At x={x}, d^2u/dx^2:\nnumerical: {d2u_num}\nexplicit:  {d2u_expl}\n");
+		}
+		for u in itertools_num::linspace(-mul * testee.internal, mul * testee.internal, 100) {
+			let xu = testee.x(u);
+			let uxu = testee.u(xu);
+			assert!(abs_diff_eq!(u, uxu, epsilon = ε), "At u={u}:\nx(u):    {xu}\nu(x(u)): {uxu}\n");
+
+			let dx_num = (testee.x(u + δ) - testee.x(u - δ)) / (2. * δ);
+			let dx_expl = testee.dx(u);
+			assert!(abs_diff_eq!(dx_expl, dx_num, epsilon = ε), "At u={u}, dx/du:\nnumerical: {dx_num}\nexplicit:  {dx_expl}\n");
+
+			let dudx = testee.du(xu) * dx_expl;
+			assert!(abs_diff_eq!(dudx, 1.0, epsilon = ε), "At u={u}:\ndu/dx * dx/du: {dudx}\n");
+		}
+	}
+}

src/bin/flat/main.rs

@@ -0,0 +1,252 @@
+use std::f32::consts::{FRAC_PI_2, PI};
+
+use flo_canvas::*;
+use flo_draw::*;
+use glam::*;
+
+use riemann::{Metric, trace_iter};
+use tube::metric::Tube;
+use tube::Space;
+use tube::Subspace::{Boundary, Inner, Outer};
+use types::{Location, Object, Ray};
+
+mod riemann;
+mod fns;
+mod float_fun;
+mod tube;
+mod types;
+
+const DT: f32 = 0.1;
+
+fn draw_loop(gc: &mut Vec<Draw>, mut pts: impl Iterator<Item=Vec2>) {
+	gc.new_path();
+	let Some(first) = pts.next() else { return; };
+	gc.move_to(first.x, first.y);
+	for pt in pts {
+		gc.line_to(pt.x, pt.y);
+	}
+	gc.close_path();
+	gc.stroke();
+}
+
+pub fn main() {
+	with_2d_graphics(move || {
+		let canvas = create_drawing_window("Refraction");
+		canvas.draw(|gc| {
+			let tube = Tube {
+				inner_radius: 30.0,
+				outer_radius: 50.0,
+				internal_halflength: 100.0,
+				external_halflength: 300.0,
+			};
+
+			let objs: Vec<_> = [-1.25, -1.00, -0.85, -0.50, 0.00, 0.40, 0.70, 0.95, 1.05]
+				.iter()
+				.enumerate()
+				.map(|(k, &y)| Object {
+					id: k as i32,
+					loc: {
+						let pos = vec2(0.0, y * tube.external_halflength);
+						let adj: Mat2 = tube.sqrt_at(pos).inverse().into();
+						let rot = Mat2::from_angle(y);
+						Location {
+							pos,
+							rot: adj * rot,
+						}
+					},
+					r: 20.0,
+				})
+				.collect();
+			let space = Space { tube, objs };
+
+			gc.canvas_height(500.0);
+			gc.transform(Transform2D::rotate(FRAC_PI_2));
+			tube.render(gc);
+			gc.line_width(0.5);
+			// gc.stroke_color(Color::Rgba(1.0, 0.5, 0.0, 0.5));
+			// draw_fan(gc, &tube, vec2(-500.0, 0.0), vec2(1.0, 0.0), 1.0);
+			gc.stroke_color(Color::Rgba(1.0, 0.5, 0.0, 1.0));
+			draw_fan_2(gc, &space, 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));
+			draw_fan_2(gc, &space, vec2(-2.5 * tube.outer_radius, 1.25 * tube.external_halflength), vec2(1.0, -1.0), 1.0);
+			draw_track(gc, &space, vec2(-500.0, 0.0), vec2(1.0, 0.2));
+			draw_track(gc, &space, vec2(-500.0, 0.0), vec2(1.0, 0.5));
+			draw_track(gc, &space, vec2(-0.5 * tube.inner_radius, -1.25 * tube.external_halflength), vec2(0.1, 1.0));
+
+			let circle_segments = 47;
+			for obj in &space.objs {
+				let pos = obj.loc.pos;
+				gc.new_path();
+				gc.circle(pos.x, pos.y, 5.0);
+				gc.fill_color(Color::Rgba(0.0, 0.5, 1.0, 1.0));
+				gc.fill();
+
+				gc.stroke_color(Color::Rgba(0.0, 0.0, 0.0, 0.5));
+				draw_loop(gc, itertools_num::linspace(0.0, 2.0 * PI, circle_segments).skip(1).map(|φ| {
+					let dir = Vec2::from_angle(φ) * obj.r;
+					let dir = obj.loc.rot * dir;
+					pos + dir
+				}));
+				gc.stroke_color(Color::Rgba(0.0, 0.5, 1.0, 0.5));
+				draw_loop(gc, itertools_num::linspace(0.0, 2.0 * PI, circle_segments).skip(1).map(|φ| {
+					let dir = Vec2::from_angle(φ) * obj.r;
+					let dir = obj.loc.rot * dir;
+					space.trace_step(Ray { pos, dir }).pos
+				}));
+				gc.stroke_color(Color::Rgba(0.5, 0.0, 1.0, 1.0));
+				draw_loop(gc, itertools_num::linspace(0.0, 2.0 * PI, circle_segments).skip(1).map(|φ| {
+					let n = obj.r.floor();
+					let d = obj.r / n;
+					let dir = Vec2::from_angle(φ);
+					let dir = obj.loc.rot * dir * d;
+					space.trace_iter(Ray { pos, dir }).nth(n as usize).unwrap().pos
+				}));
+			}
+		});
+	});
+}
+
+fn draw_ray_2(gc: &mut Vec<Draw>, space: &Space, base: Vec2, dir: Vec2) {
+	let mut hits = Vec::<Draw>::new();
+	let dir = space.tube.globalize(base, dir);
+	gc.new_path();
+	gc.move_to(base.x, base.y);
+	let mut ray = Ray { pos: base, dir: space.tube.normalize_vec_at(base, dir) * DT };
+	for _ in 0..10000 {
+		ray = space.trace_step(ray);
+		gc.line_to(ray.pos.x, ray.pos.y);
+		if ray.pos.abs().cmpgt(Vec2::splat(1000.0)).any() {
+			break;
+		}
+		let sub = space.which_subspace(ray.pos);
+		if sub == Boundary {
+			continue;
+		}
+		gc.stroke();
+		gc.new_dash_pattern();
+		// gc.dash_length(6.0);
+		gc.new_path();
+		gc.move_to(ray.pos.x, ray.pos.y);
+		let ret = match sub {
+			Inner => space.trace_inner(ray),
+			Outer => space.trace_outer(ray),
+			Boundary => panic!(),
+		};
+		for hit in ret.objects {
+			let obj = space.objs[hit.id as usize];
+			hits.move_to(obj.loc.pos.x, obj.loc.pos.y);
+			for pt in trace_iter(&space.tube, obj.loc.pos, obj.loc.rot * hit.rel.pos, hit.rel.pos.length() / 100.0).take(100) {
+				hits.line_to(pt.x, pt.y);
+			}
+			hits.circle(hit.pos.x, hit.pos.y, 1.5);
+			let Ray { pos: rel, dir } = hit.rel;
+			let diff = rel.dot(dir).powi(2) - dir.length_squared() * (rel.length_squared() - obj.r.powi(2));
+			assert!(diff >= 0.0);
+			let t = (-rel.dot(dir) + diff.sqrt()) / dir.length_squared();
+			let rel2 = hit.rel.forward(t).pos;
+			let pos2 = trace_iter(&space.tube, obj.loc.pos, obj.loc.rot * rel2, rel2.length() / 100.0).nth(100).unwrap();
+			hits.move_to(pos2.x - 1.0, pos2.y - 1.0);
+			hits.line_to(pos2.x + 1.0, pos2.y + 1.0);
+			hits.move_to(pos2.x - 1.0, pos2.y + 1.0);
+			hits.line_to(pos2.x + 1.0, pos2.y - 1.0);
+		}
+		let a = ray.pos;
+		ray = match ret.end {
+			Some(r) => r,
+			None => {
+				ray = ray.forward(1000.0 / DT);
+				gc.line_to(ray.pos.x, ray.pos.y);
+				break;
+			}
+		};
+		for p in space.line(a, ray.pos, 10.0) {
+			gc.line_to(p.x, p.y);
+		}
+		gc.stroke();
+		gc.new_dash_pattern();
+		gc.new_path();
+		gc.move_to(ray.pos.x, ray.pos.y);
+	}
+	gc.stroke();
+	gc.new_path();
+	gc.new_dash_pattern();
+	gc.append(&mut hits);
+	gc.stroke();
+}
+
+fn draw_fan_2(gc: &mut Vec<Draw>, space: &Space, 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_2(gc, space, base, dir + y * v);
+	}
+}
+
+fn draw_ray(gc: &mut Vec<Draw>, space: &impl Metric, base: Vec2, dir: Vec2) {
+	let dir = space.globalize(base, dir);
+	gc.new_path();
+	gc.move_to(base.x, base.y);
+	for pt in trace_iter(space, 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_track(gc: &mut Vec<Draw>, space: &Space, start: Vec2, dir: Vec2) {
+	const SCALE: f32 = 5.0;
+	const STEP: f32 = 2.0 * SCALE;
+	// let mut loc = Location { pos: start, rot: Mat2::IDENTITY };
+	// let dir = space.tube.globalize(start, dir);
+	// let v = space.tube.normalize(start, dir);
+	let mut loc = Location { pos: start, rot: mat2(dir, vec2(-dir.y, dir.x)) };
+	let v = vec2(1.0, 0.0);
+	let mut draw = |loc: &Location| {
+		let p = loc.pos;
+		let ax = p + loc.rot.x_axis * SCALE;
+		let ay = p + loc.rot.y_axis * SCALE;
+		gc.new_path();
+		gc.stroke_color(Color::Rgba(0.7, 0.0, 0.0, 1.0));
+		gc.move_to(p.x, p.y);
+		gc.line_to(ax.x, ax.y);
+		gc.stroke();
+		gc.new_path();
+		gc.stroke_color(Color::Rgba(0.0, 0.7, 0.0, 1.0));
+		gc.move_to(p.x, p.y);
+		gc.line_to(ay.x, ay.y);
+		gc.stroke();
+	};
+	draw(&loc);
+	for _ in 0..1000 {
+		let N = (STEP / DT).floor() as i32;
+		for _ in 0..N {
+			loc = space.move_step(loc, v * DT);
+		}
+		draw(&loc);
+	}
+}
+
+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 {
+	fn render(&self, gc: &mut Vec<Draw>);
+}
+
+impl Renderable for Tube {
+	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.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();
+	}
+}

src/bin/flat/riemann.rs

@@ -0,0 +1,195 @@
+use glam::*;
+
+#[derive(Copy, Clone, Debug, PartialEq)]
+pub struct Decomp2 {
+	pub ortho: Mat2,
+	pub diag: Vec2,
+}
+
+impl Decomp2 {
+	fn square(&self) -> Self {
+		Self {
+			ortho: self.ortho,
+			diag: self.diag * self.diag,
+		}
+	}
+
+	pub fn inverse(&self) -> Self {
+		Self {
+			ortho: self.ortho,
+			diag: Vec2::splat(1.0) / self.diag,
+		}
+	}
+}
+
+impl From<Decomp2> for Mat2 {
+	fn from(value: Decomp2) -> Self {
+		value.ortho.transpose() * Mat2::from_diagonal(value.diag) * value.ortho
+	}
+}
+
+pub type Tens2 = [Mat2; 2];
+
+pub trait Metric {
+	fn sqrt_at(&self, pos: Vec2) -> Decomp2;
+
+	fn at(&self, pos: Vec2) -> Mat2 {
+		self.sqrt_at(pos).square().into()
+	}
+
+	fn inverse_at(&self, pos: Vec2) -> Mat2 {
+		self.sqrt_at(pos).square().inverse().into()
+	}
+
+	fn part_derivs_at(&self, pos: Vec2) -> Tens2 {
+		part_deriv(|p| self.at(p), pos, 1.0 / 1024.0) // division by such eps is exact which is good for overall precision
+	}
+
+	fn vec_length_at(&self, at: Vec2, v: Vec2) -> f32 {
+		v.dot(self.at(at) * v).sqrt()
+	}
+
+	fn normalize_vec_at(&self, at: Vec2, v: Vec2) -> Vec2 {
+		v / self.vec_length_at(at, v)
+	}
+
+	fn globalize(&self, at: Vec2, v: Vec2) -> Vec2 {
+		Mat2::from(self.sqrt_at(at).inverse()) * v
+	}
+}
+
+pub struct TraceIter<'a, M: Metric> {
+	space: &'a M,
+	p: Vec2,
+	v: Vec2,
+	dt: f32,
+}
+
+impl<'a, M: Metric> Iterator for TraceIter<'a, M> {
+	type Item = Vec2;
+
+	fn next(&mut self) -> Option<Self::Item> {
+		let a: Vec2 = -contract2(krist(self.space, self.p), self.v);
+		self.v = self.v + a * self.dt;
+		self.p = self.p + self.v * self.dt;
+		Some(self.p)
+	}
+}
+
+pub fn trace_iter<M: Metric>(space: &M, base: Vec2, dir: Vec2, dt: f32) -> TraceIter<M> {
+	TraceIter {
+		space,
+		p: base,
+		v: space.normalize_vec_at(base, dir),
+		dt,
+	}
+}
+
+pub fn krist(space: &impl Metric, pos: Vec2) -> Tens2 {
+	// Γ^i_k_l = .5 * g^i^m * (g_m_k,l + g_m_l,k - g_k_l,m)
+	let g = &space.inverse_at(pos); // с верхними индексами
+	let d = space.part_derivs_at(pos);
+	// ret[i][l][k] = sum((m) => .5f * g[m][i] * (d[k][l][m] + d[l][k][m] - d[m][k][l]))
+	make_tens2(|i, l, k| 0.5 * (0..2).map(|m| g.col(m)[i] * (d[l].col(k)[m] + d[k].col(m)[l] - d[m].col(k)[l])).sum::<f32>())
+}
+
+fn dir_deriv(f: impl Fn(Vec2) -> Mat2, pos: Vec2, delta: Vec2) -> Mat2 {
+	(f(pos + delta) - f(pos - delta)) / (2.0 * delta.length())
+}
+
+fn part_deriv(f: impl Fn(Vec2) -> Mat2, pos: Vec2, eps: f32) -> Tens2 {
+	[
+		dir_deriv(&f, pos, vec2(eps, 0.0)),
+		dir_deriv(&f, pos, vec2(0.0, eps)),
+	]
+}
+
+/// Сворачивает тензор t с вектором u
+pub fn contract(t: Tens2, u: Vec2) -> Mat2 {
+	mat2(t[0] * u, t[1] * u).transpose()
+}
+
+/// Сворачивает тензор t с вектором v дважды, по второму и третьему индексам.
+pub fn contract2(t: Tens2, v: Vec2) -> Vec2 {
+	contract(t, v) * v
+}
+
+fn make_vec2(f: impl Fn(usize) -> f32) -> Vec2 {
+	Vec2::from_array(std::array::from_fn(|i| f(i)))
+}
+
+fn make_mat2(f: impl Fn(usize, usize) -> f32) -> Mat2 {
+	Mat2::from_cols_array_2d(&std::array::from_fn(|i| std::array::from_fn(|j| f(i, j))))
+}
+
+fn make_tens2(f: impl Fn(usize, usize, usize) -> f32) -> Tens2 {
+	std::array::from_fn(|i| make_mat2(|j, k| f(i, j, k)))
+}
+
+#[test]
+fn m2() {
+	let m = make_mat2(|i, j| (i + 2 * j) as f32);
+	assert_eq!(m.col(0)[0], 0.0);
+	assert_eq!(m.col(1)[0], 1.0);
+	assert_eq!(m.col(0)[1], 2.0);
+	assert_eq!(m.col(1)[1], 3.0);
+}
+
+#[test]
+fn t2() {
+	let t = make_tens2(|i, j, k| (i + 2 * j + 4 * k) as f32);
+	assert_eq!(t[0].col(0)[0], 0.0);
+	assert_eq!(t[1].col(0)[0], 1.0);
+	assert_eq!(t[0].col(1)[0], 2.0);
+	assert_eq!(t[1].col(1)[0], 3.0);
+	assert_eq!(t[0].col(0)[1], 4.0);
+	assert_eq!(t[1].col(0)[1], 5.0);
+	assert_eq!(t[0].col(1)[1], 6.0);
+	assert_eq!(t[1].col(1)[1], 7.0);
+}
+
+pub mod samples {
+	use glam::{Mat2, Vec2};
+
+	use super::{Decomp2, Metric};
+
+	pub struct ScaledMetric {
+		pub scale: Vec2,
+	}
+
+	impl Metric for ScaledMetric {
+		fn sqrt_at(&self, _pos: Vec2) -> Decomp2 {
+			Decomp2 {
+				diag: self.scale,
+				ortho: Mat2::IDENTITY,
+			}
+		}
+	}
+}
+
+#[cfg(test)]
+mod tests {
+	use glam::{Mat2, mat2, vec2};
+	use rand::{Rng, SeedableRng};
+
+	use super::{Decomp2, Metric, samples};
+
+	#[test]
+	fn uniform_scaled_metric() {
+		let mut rng = rand_pcg::Pcg64Mcg::seed_from_u64(17);
+		let metric = samples::ScaledMetric { scale: vec2(3., 4.) };
+		assert_eq!(metric.sqrt_at(rng.gen()), Decomp2 { ortho: Mat2::IDENTITY, diag: vec2(3., 4.) });
+		assert_eq!(metric.at(rng.gen()), Mat2::from_cols_array(&[9., 0., 0., 16.]));
+		assert_eq!(metric.inverse_at(rng.gen()), Mat2::from_cols_array(&[1. / 9., 0., 0., 1. / 16.]));
+		assert_eq!(metric.part_derivs_at(rng.gen()), [Mat2::ZERO, Mat2::ZERO]);
+		assert_eq!(metric.vec_length_at(rng.gen(), vec2(1., 0.)), 3.);
+		assert_eq!(metric.vec_length_at(rng.gen(), vec2(0., 1.)), 4.);
+		assert_eq!(metric.vec_length_at(rng.gen(), vec2(1., 1.)), 5.);
+		assert_eq!(metric.normalize_vec_at(rng.gen(), vec2(1., 0.)), vec2(1. / 3., 0.));
+		assert_eq!(metric.normalize_vec_at(rng.gen(), vec2(0., 1.)), vec2(0., 1. / 4.));
+		assert_eq!(metric.normalize_vec_at(rng.gen(), vec2(1., 1.)), vec2(1. / 5., 1. / 5.));
+		assert_eq!(metric.globalize(rng.gen(), vec2(1., 0.)), vec2(1. / 3., 0.));
+		assert_eq!(metric.globalize(rng.gen(), vec2(0., 1.)), vec2(0., 1. / 4.));
+		assert_eq!(metric.globalize(rng.gen(), vec2(1., 1.)), vec2(1. / 3., 1. / 4.));
+	}
+}

src/bin/flat/tube/coords.rs

@@ -0,0 +1,248 @@
+use glam::{Mat2, Vec2, vec2};
+
+use crate::riemann::Metric;
+use crate::types::{Location, Ray};
+
+use super::{Rect, Tube};
+
+pub trait FlatCoordinateSystem<T> {
+	fn flat_to_global(&self, v: T) -> T;
+	fn global_to_flat(&self, v: T) -> T;
+}
+
+pub trait FlatRegion: FlatCoordinateSystem<Vec2> + FlatCoordinateSystem<Ray> + FlatCoordinateSystem<Location> {
+	// Измеряет расстояние до выхода за пределы области вдоль луча ray. Луч задаётся в плоской СК.
+	fn distance_to_boundary(&self, _ray: Ray) -> Option<f32> { None }
+}
+
+trait MetricCS: FlatCoordinateSystem<Vec2> {
+	fn global_metric(&self) -> &impl Metric;
+	fn flat_to_global_tfm_at(&self, pos: Vec2) -> Mat2 {
+		self.global_metric().sqrt_at(self.flat_to_global(pos)).inverse().into()
+	}
+	fn global_to_flat_tfm_at(&self, pos: Vec2) -> Mat2 {
+		self.global_metric().sqrt_at(pos).into()
+	}
+}
+
+impl<T: FlatCoordinateSystem<Vec2> + MetricCS> FlatCoordinateSystem<Ray> for T {
+	fn flat_to_global(&self, ray: Ray) -> Ray {
+		Ray {
+			pos: self.flat_to_global(ray.pos),
+			dir: self.flat_to_global_tfm_at(ray.pos) * ray.dir,
+		}
+	}
+
+	fn global_to_flat(&self, ray: Ray) -> Ray {
+		Ray {
+			pos: self.global_to_flat(ray.pos),
+			dir: self.global_to_flat_tfm_at(ray.pos) * ray.dir,
+		}
+	}
+}
+
+impl<T: FlatCoordinateSystem<Vec2> + MetricCS> FlatCoordinateSystem<Location> for T {
+	fn flat_to_global(&self, loc: Location) -> Location {
+		Location {
+			pos: self.flat_to_global(loc.pos),
+			rot: self.flat_to_global_tfm_at(loc.pos) * loc.rot,
+		}
+	}
+
+	fn global_to_flat(&self, loc: Location) -> Location {
+		Location {
+			pos: self.global_to_flat(loc.pos), // в плоской СК для Inner или её продолжении на Outer
+			rot: self.global_to_flat_tfm_at(loc.pos) * loc.rot,
+		}
+	}
+}
+
+pub struct InnerCS(pub Tube);
+
+impl MetricCS for InnerCS { fn global_metric(&self) -> &impl Metric { &self.0 } }
+
+impl FlatCoordinateSystem<Vec2> for InnerCS {
+	fn flat_to_global(&self, pos: Vec2) -> Vec2 {
+		vec2(pos.x, self.0.y(pos.y))
+	}
+
+	// Работает только при |pos.x| ≤ inner_radius или |pos.y| ≥ external_halflength.
+	fn global_to_flat(&self, pos: Vec2) -> Vec2 {
+		vec2(pos.x, self.0.v(pos.y))
+	}
+}
+
+impl FlatRegion for InnerCS {
+	fn distance_to_boundary(&self, ray: Ray) -> Option<f32> {
+		Rect { size: vec2(self.0.inner_radius, self.0.internal_halflength) }.trace_out_of(ray)
+	}
+}
+
+pub struct OuterCS(pub Tube);
+
+impl MetricCS for OuterCS { fn global_metric(&self) -> &impl Metric { &self.0 } }
+
+impl FlatCoordinateSystem<Vec2> for OuterCS {
+	fn flat_to_global(&self, pos: Vec2) -> Vec2 {
+		let inner = Rect { size: vec2(self.0.inner_radius + 1.0, self.0.external_halflength) };
+		if inner.is_inside(pos) {
+			let Vec2 { x, y: v } = pos;
+			let y = self.0.y(v - v.signum() * (self.0.external_halflength - self.0.internal_halflength));
+			vec2(x, y)
+		} else {
+			pos
+		}
+	}
+
+	fn global_to_flat(&self, pos: Vec2) -> Vec2 {
+		let inner = Rect { size: vec2(self.0.inner_radius + 1.0, self.0.external_halflength) };
+		if inner.is_inside(pos) {
+			let Vec2 { x: u, y } = pos; // в основной СК
+			let v = self.0.v(y) + y.signum() * (self.0.external_halflength - self.0.internal_halflength);
+			vec2(u, v) // в плоском продолжении СК Outer на область Inner
+		} else {
+			pos
+		}
+	}
+}
+
+impl FlatRegion for OuterCS {
+	fn distance_to_boundary(&self, ray: Ray) -> Option<f32> {
+		Rect { size: vec2(self.0.outer_radius, self.0.external_halflength) }.trace_into(ray)
+	}
+}
+
+#[cfg(test)]
+mod tests {
+	use approx::{AbsDiffEq, assert_abs_diff_eq};
+	use glam::{Mat2, mat2, vec2, Vec2};
+	use itertools_num::linspace;
+
+	use crate::riemann::samples;
+
+	use super::*;
+
+	#[test]
+	fn uniform_scaled_metric() {
+		struct Scaled(samples::ScaledMetric, Vec2);
+		impl FlatCoordinateSystem<Vec2> for Scaled {
+			fn flat_to_global(&self, pos: Vec2) -> Vec2 { (pos - self.1) / self.0.scale }
+			fn global_to_flat(&self, pos: Vec2) -> Vec2 { pos * self.0.scale + self.1 }
+		}
+		impl MetricCS for Scaled {
+			fn global_metric(&self) -> &impl Metric { &self.0 }
+		}
+		let cs = Scaled(samples::ScaledMetric { scale: vec2(3., 4.) }, vec2(2., 3.));
+		assert_eq!(cs.global_to_flat(vec2(7., 3.)), vec2(23., 15.));
+		assert_eq!(cs.flat_to_global(vec2(8., 7.)), vec2(2., 1.));
+		assert_eq!(cs.global_to_flat(Ray { pos: vec2(7., 3.), dir: vec2(3., 2.) }), Ray { pos: vec2(23., 15.), dir: vec2(9., 8.) });
+		assert_eq!(cs.flat_to_global(Ray { pos: vec2(23., 15.), dir: vec2(9., 8.) }), Ray { pos: vec2(7., 3.), dir: vec2(3., 2.) });
+		assert_eq!(cs.global_to_flat(Location { pos: vec2(2., 1.), rot: mat2(vec2(0., 1.), vec2(-1., 0.)) }), Location { pos: vec2(8., 7.), rot: mat2(vec2(0., 4.), vec2(-3., 0.)) });
+		assert_eq!(cs.flat_to_global(Location { pos: vec2(2., 1.), rot: mat2(vec2(0., 1.), vec2(-1., 0.)) }), Location { pos: vec2(0., -0.5), rot: mat2(vec2(0., 0.25), vec2(-1. / 3., 0.)) });
+	}
+
+	fn test_flat_region(region: &impl FlatRegion, range_global: (Vec2, Vec2), range_flat: (Vec2, Vec2)) {
+		const ε: f32 = 1e-3;
+		macro_rules! assert_eq_at {
+			($at: expr, $left: expr, $right: expr) => {
+				let at = $at;
+				let left = $left;
+				let right = $right;
+				assert!(left.abs_diff_eq(right, ε), "Assertion failed at {at}:\n  left: {left} = {}\n right: {right} = {}", stringify!($left), stringify!($right));
+			};
+		}
+		fn check_range(name_a: &str, a: Vec2, range_a: (Vec2, Vec2), name_b: &str, b: Vec2, range_b: (Vec2, Vec2)) {
+			assert!(b.cmpge(range_b.0 - ε).all() && b.cmple(range_b.1 + ε).all(), "Assertion failed:\nAt {name_a}: {a}, from range: {range_a:?}\nGot {name_b}: {b}, which is out of range {range_b:?}");
+			// TODO sort out when to check these conditions:
+			if a.x.abs_diff_eq(&range_a.0.x, ε) { assert_abs_diff_eq!(b.x, range_b.0.x, epsilon=ε); }
+			if a.y.abs_diff_eq(&range_a.0.y, ε) { assert_abs_diff_eq!(b.y, range_b.0.y, epsilon=ε); }
+			if a.x.abs_diff_eq(&range_a.1.x, ε) { assert_abs_diff_eq!(b.x, range_b.1.x, epsilon=ε); }
+			if a.y.abs_diff_eq(&range_a.1.y, ε) { assert_abs_diff_eq!(b.y, range_b.1.y, epsilon=ε); }
+		}
+		for x in linspace(range_global.0.x, range_global.1.x, 20) {
+			for y in linspace(range_global.0.y, range_global.1.y, 20) {
+				let pos_global = vec2(x, y);
+				let pos_flat = region.global_to_flat(pos_global);
+				check_range("global", pos_global, range_global, "flat", pos_flat, range_flat);
+				assert_eq_at!(pos_global, region.global_to_flat(Location { pos: pos_global, rot: Mat2::IDENTITY }).pos, pos_flat);
+				assert_eq_at!(pos_global, region.flat_to_global(pos_flat), pos_global);
+				assert_eq_at!(pos_global, region.flat_to_global(region.global_to_flat(Location { pos: pos_global, rot: Mat2::IDENTITY })).rot, Mat2::IDENTITY);
+			}
+		}
+		for x in linspace(range_flat.0.x, range_flat.1.x, 20) {
+			for y in linspace(range_flat.0.y, range_flat.1.y, 20) {
+				let pos_flat = vec2(x, y);
+				let pos_global = region.flat_to_global(pos_flat);
+				check_range("flat", pos_flat, range_flat, "global", pos_global, range_global);
+				assert_eq_at!(pos_flat, region.flat_to_global(Location { pos: pos_flat, rot: Mat2::IDENTITY }).pos, pos_global);
+				assert_eq_at!(pos_flat, region.global_to_flat(pos_global), pos_flat);
+				assert_eq_at!(pos_flat, region.global_to_flat(region.flat_to_global(Location { pos: pos_global, rot: Mat2::IDENTITY })).rot, Mat2::IDENTITY);
+			}
+		}
+	}
+
+	#[test]
+	fn test_mapper_inner() {
+		let mapper = InnerCS(Tube {
+			inner_radius: 30.0,
+			outer_radius: 50.0,
+			internal_halflength: 100.0,
+			external_halflength: 300.0,
+		});
+		test_flat_region(&mapper, (vec2(-30.0, -300.0), vec2(30.0, 300.0)), (vec2(-30.0, -100.0), vec2(30.0, 100.0)));
+		test_flat_region(&mapper, (vec2(-60.0, -400.0), vec2(60.0, -300.0)), (vec2(-60.0, -200.0), vec2(60.0, -100.0)));
+		test_flat_region(&mapper, (vec2(-60.0, 300.0), vec2(60.0, 400.0)), (vec2(-60.0, 100.0), vec2(60.0, 200.0)));
+	}
+
+	#[test]
+	fn test_mapper_outer() {
+		let mapper = OuterCS(Tube {
+			inner_radius: 30.0,
+			outer_radius: 50.0,
+			internal_halflength: 100.0,
+			external_halflength: 300.0,
+		});
+		// TODO replace 200.20016 with something sane
+		test_flat_region(&mapper, (vec2(-30.0, -300.0), vec2(30.0, -1.0)), (vec2(-30.0, -300.0), vec2(30.0, -200.20016)));
+		test_flat_region(&mapper, (vec2(-30.0, 1.0), vec2(30.0, 300.0)), (vec2(-30.0, 200.20016), vec2(30.0, 300.0)));
+		test_flat_region(&mapper, (vec2(-60.0, -400.0), vec2(60.0, -300.0)), (vec2(-60.0, -400.0), vec2(60.0, -300.0)));
+		test_flat_region(&mapper, (vec2(-60.0, 300.0), vec2(60.0, 400.0)), (vec2(-60.0, 300.0), vec2(60.0, 400.0)));
+		// straight
+		for x in linspace(-60., 60., 20) {
+			for y in linspace(-320., 320., 20) {
+				assert_eq!(mapper.global_to_flat(Location { pos: vec2(x, y), rot: Mat2::IDENTITY }).pos.x, x);
+			}
+		}
+		// symmetrical
+		for x in linspace(0., 60., 20) {
+			for y in linspace(0., 320., 20) {
+				let pp = mapper.global_to_flat(Location { pos: vec2(x, y), rot: Mat2::IDENTITY }).pos;
+				let np = mapper.global_to_flat(Location { pos: vec2(-x, y), rot: Mat2::IDENTITY }).pos;
+				let pn = mapper.global_to_flat(Location { pos: vec2(x, -y), rot: Mat2::IDENTITY }).pos;
+				let nn = mapper.global_to_flat(Location { pos: vec2(-x, -y), rot: Mat2::IDENTITY }).pos;
+				assert_eq!(np, vec2(-pp.x, pp.y));
+				assert_eq!(pn, vec2(pp.x, -pp.y));
+				assert_eq!(nn, vec2(-pp.x, -pp.y));
+			}
+		}
+		// clean boundary
+		for x in linspace(50., 60., 20) {
+			for y in linspace(0., 320., 20) {
+				assert_eq!(mapper.global_to_flat(Location { pos: vec2(x, y), rot: Mat2::IDENTITY }).pos.y, y);
+			}
+		}
+		for x in linspace(0., 60., 20) {
+			for y in linspace(300., 320., 20) {
+				assert_eq!(mapper.global_to_flat(Location { pos: vec2(x, y), rot: Mat2::IDENTITY }).pos.y, y);
+			}
+		}
+		// accelerating
+		for x in linspace(-29., 29., 20) {
+			for y in linspace(1., 299., 20) {
+				let v = mapper.global_to_flat(Location { pos: vec2(x, y), rot: Mat2::IDENTITY }).pos.y;
+				assert!(v > 200.0);
+				assert!(v > y);
+			}
+		}
+	}
+}

src/bin/flat/tube/metric.rs

@@ -0,0 +1,204 @@
+use glam::{f32, Mat2, Vec2, vec2};
+
+use crate::fns::{self, Limiter};
+use crate::riemann::{Decomp2, Metric, Tens2};
+
+#[derive(Copy, Clone, Debug)]
+pub struct Tube {
+	pub outer_radius: f32,
+	pub inner_radius: f32,
+	pub external_halflength: f32,
+	pub internal_halflength: f32,
+}
+
+impl Tube {
+	fn fx(&self) -> impl Limiter { 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) }
+	pub fn v(&self, y: f32) -> f32 { self.fy().u(y) }
+	pub fn dy(&self, v: f32) -> f32 { self.fy().dx(v) }
+	pub fn dv(&self, y: f32) -> f32 { self.fy().du(y) }
+}
+
+impl Metric for Tube {
+	fn sqrt_at(&self, pos: Vec2) -> Decomp2 {
+		let sx = self.fx().value(pos.x);
+		let sy = self.fy().du(pos.y);
+		let s = sx + sy - sx * sy;
+		assert!(sx.is_finite());
+		assert!(sy.is_finite());
+		assert!(sy > 0.0);
+		Decomp2 {
+			ortho: Mat2::IDENTITY,
+			diag: vec2(1.0, s),
+		}
+	}
+
+	fn part_derivs_at(&self, pos: Vec2) -> Tens2 {
+		let sx = self.fx().value(pos.x);
+		let sy = self.fy().du(pos.y);
+		let s = sx + sy - sx * sy;
+		let dsx_dx = self.fx().derivative(pos.x);
+		let dsy_dy = self.fy().d2u(pos.y);
+		let ds2_dx = 2.0 * s * (1.0 - sy) * dsx_dx;
+		let ds2_dy = 2.0 * s * (1.0 - sx) * dsy_dy;
+		[
+			Mat2::from_cols_array(&[0.0, 0.0, 0.0, ds2_dx]),
+			Mat2::from_cols_array(&[0.0, 0.0, 0.0, ds2_dy]),
+		]
+	}
+}
+
+#[cfg(test)]
+mod test {
+	use approx::assert_abs_diff_eq;
+	use glam::{Vec2, vec2};
+	use itertools_num::linspace;
+
+	use crate::riemann::{Decomp2, Metric};
+	use crate::tube::Space;
+	use crate::types::Ray;
+
+	use super::Tube;
+
+	#[test]
+	fn test_tube_metric_derivs() {
+		struct Approx(Tube);
+		impl Metric for Approx {
+			fn sqrt_at(&self, pos: Vec2) -> Decomp2 { self.0.sqrt_at(pos) }
+		}
+		let testee = Tube {
+			inner_radius: 30.0,
+			outer_radius: 50.0,
+			internal_halflength: 100.0,
+			external_halflength: 300.0,
+		};
+		let approx = Approx(testee);
+		let epsilon = 1.0e-3;
+		let margin = 1.0 / 16.0;
+		let mul = 1.0 + margin;
+		for x in itertools_num::linspace(-mul * testee.outer_radius, mul * testee.outer_radius, 100) {
+			for y in itertools_num::linspace(-mul * testee.external_halflength, mul * testee.external_halflength, 100) {
+				let pos = vec2(x, y);
+				let computed = testee.part_derivs_at(pos);
+				let reference = approx.part_derivs_at(pos);
+				let eq = (0..2).all(|coord| computed[coord].abs_diff_eq(reference[coord], epsilon));
+				assert!(eq, "Bad derivative computation at {pos}:\n  explicit:  {computed:?}\n  numerical: {reference:?}\n");
+			}
+		}
+	}
+
+	#[test]
+	fn test_accelerator() {
+		let space = Space {
+			tube: Tube {
+				inner_radius: 30.0,
+				outer_radius: 50.0,
+				internal_halflength: 100.0,
+				external_halflength: 300.0,
+			},
+			objs: vec![],
+		};
+		let ε = 1e-3;
+		let off = 10.0;
+		let steps = 1024;
+		for ax in [-30.0 + ε, -25.0, -3.0, 17.0, 30.0 - ε] {
+			for bx in [0.0, ε, 1.0, 7.0, 30.0 - ε] {
+				let a = vec2(ax, -(space.tube.external_halflength + off));
+				let b = vec2(bx, space.tube.external_halflength + off);
+				let Δ = vec2(bx - ax, 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.0e-2);
+				assert_abs_diff_eq!(traced.pos.y, b.y, epsilon=1.0e1);
+				assert_abs_diff_eq!(traced.dir.x, dir.x, epsilon=1.0e-3);
+				assert_abs_diff_eq!(traced.dir.y, dir.y, epsilon=1.0e-2);
+			}
+		}
+	}
+
+	#[test]
+	#[ignore]
+	fn test_accelerator_slow() {
+		let space = Space {
+			tube: Tube {
+				inner_radius: 30.0,
+				outer_radius: 50.0,
+				internal_halflength: 100.0,
+				external_halflength: 300.0,
+			},
+			objs: vec![],
+		};
+		let ε = 1e-3;
+		let off = 10.0;
+		let steps = 4096;
+		for ax in linspace(-space.tube.inner_radius + ε, space.tube.inner_radius - ε, 20) {
+			for bx in linspace(-space.tube.inner_radius + ε, space.tube.inner_radius - ε, 20) {
+				let a = vec2(ax, -(space.tube.external_halflength + off));
+				let b = vec2(bx, space.tube.external_halflength + off);
+				let Δ = vec2(bx - ax, 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.0e-2);
+				assert_abs_diff_eq!(traced.pos.y, b.y, epsilon=1.0e0);
+				assert_abs_diff_eq!(traced.dir.x, dir.x, epsilon=1.0e-3);
+				assert_abs_diff_eq!(traced.dir.y, dir.y, epsilon=1.0e-3);
+			}
+		}
+	}
+
+	#[test]
+	fn test_accelerator_inner_edge() {
+		let space = Space {
+			tube: Tube {
+				inner_radius: 30.0,
+				outer_radius: 50.0,
+				internal_halflength: 100.0,
+				external_halflength: 300.0,
+			},
+			objs: vec![],
+		};
+		let ε = 1e-3;
+		let off = 10.0;
+		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.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);
+		}
+	}
+
+	#[test]
+	fn test_accelerator_outer_edge() {
+		let space = Space {
+			tube: Tube {
+				inner_radius: 30.0,
+				outer_radius: 50.0,
+				internal_halflength: 100.0,
+				external_halflength: 300.0,
+			},
+			objs: vec![],
+		};
+		let ε = 1e-3;
+		let off = 10.0;
+		let steps = 4096;
+		for x in [space.tube.outer_radius + ε, space.tube.outer_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.external_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=2.0e0);
+			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);
+		}
+	}
+}

src/bin/flat/tube/mod.rs

@@ -0,0 +1,187 @@
+use glam::{bool, f32, Mat2, Vec2, vec2};
+
+use coords::{FlatCoordinateSystem, InnerCS, OuterCS};
+use metric::Tube;
+use Subspace::{Boundary, Inner, Outer};
+
+use crate::riemann;
+use crate::tube::coords::FlatRegion;
+use crate::types::{FlatTraceResult, Hit, Location, Object, Ray};
+
+pub mod metric;
+mod coords;
+
+pub struct Space {
+	pub tube: Tube,
+	pub objs: Vec<Object>,
+}
+
+#[derive(PartialEq, Eq, Debug)]
+pub enum Subspace {
+	Outer,
+	Boundary,
+	Inner,
+}
+
+impl Space {
+	pub fn which_subspace(&self, pt: Vec2) -> Subspace {
+		if pt.y.abs() > self.tube.external_halflength {
+			Outer
+		} else if pt.x.abs() > self.tube.outer_radius {
+			Outer
+		} else if pt.x.abs() > self.tube.inner_radius {
+			Boundary
+		} else {
+			Inner
+		}
+	}
+
+	/// Выполняет один шаг трассировки. Работает в любой части пространства, но вне Boundary доступны более эффективные методы.
+	/// ray задаётся в основной СК.
+	pub fn trace_step(&self, ray: Ray) -> Ray {
+		let a: Vec2 = -riemann::contract2(riemann::krist(&self.tube, ray.pos), ray.dir);
+		let v = ray.dir + a;
+		let p = ray.pos + v;
+		Ray { pos: p, dir: v }
+	}
+
+	/// Выполняет один шаг перемещения. Работает в любой части пространства.
+	/// off задаётся в локальной СК. Рекомендуется считать небольшими шагами.
+	pub fn move_step(&self, loc: Location, off: Vec2) -> Location {
+		let corr = Mat2::IDENTITY - riemann::contract(riemann::krist(&self.tube, loc.pos), loc.rot * off);
+		let p = loc.pos + corr * loc.rot * off;
+		Location { pos: p, rot: corr * loc.rot }
+	}
+
+	pub fn trace_iter(&self, ray: Ray) -> impl Iterator<Item=Ray> + '_ {
+		std::iter::successors(Some(ray), |&ray| Some(self.trace_step(ray)))
+	}
+
+	pub fn trace_inner(&self, ray: Ray) -> FlatTraceResult {
+		assert_eq!(self.which_subspace(ray.pos), Inner);
+		self.trace_flat(InnerCS(self.tube), ray)
+	}
+
+	pub fn trace_outer(&self, ray: Ray) -> FlatTraceResult {
+		assert_eq!(self.which_subspace(ray.pos), Outer);
+		self.trace_flat(OuterCS(self.tube), ray)
+	}
+
+	fn trace_flat(&self, cs: impl FlatRegion, ray: Ray) -> FlatTraceResult {
+		let ray = cs.global_to_flat(ray);
+		let dist = cs.distance_to_boundary(ray);
+		let objs = self.list_objects(|loc| cs.global_to_flat(loc));
+		FlatTraceResult {
+			end: dist.map(|dist| cs.flat_to_global(ray.forward(dist))),
+			objects: Self::hit_objects(objs.as_slice(), ray, dist, |pos| cs.flat_to_global(pos)),
+		}
+	}
+
+	fn trace_boundary(&self, ray: Ray) -> Ray {
+		assert_eq!(self.which_subspace(ray.pos), Boundary);
+		self.trace_iter(ray)
+			.find(|&ray| self.which_subspace(ray.pos) != Boundary)
+			.expect("Can't get outta the wall!")
+	}
+
+	fn list_objects(&self, tfm: impl Fn(Location) -> Location) -> Vec<Object> {
+		self.objs.iter().map(|&Object { id, loc, r }| Object { id, loc: tfm(loc), r }).collect()
+	}
+
+	fn hit_objects(objs: &[Object], ray: Ray, limit: Option<f32>, globalize: impl Fn(Vec2) -> Vec2) -> Vec<Hit> {
+		let limit = limit.unwrap_or(f32::INFINITY);
+		objs.iter()
+			.filter_map(|obj| {
+				let rel = ray.pos - obj.loc.pos;
+				let diff = rel.dot(ray.dir).powi(2) - ray.dir.length_squared() * (rel.length_squared() - obj.r.powi(2));
+				if diff > 0.0 {
+					let t = (-rel.dot(ray.dir) - diff.sqrt()) / ray.dir.length_squared();
+					Some((obj, t))
+				} else {
+					None
+				}
+			})
+			.filter(|&(_, t)| t >= 0.0 && t < limit)
+			.map(|(obj, t)| {
+				let pos = ray.forward(t).pos;
+				let rel = obj.loc.rot.inverse() * Ray { pos: pos - obj.loc.pos, dir: ray.dir };
+				Hit { id: obj.id, distance: t, pos: globalize(pos), rel }
+			})
+			.collect()
+	}
+
+	pub fn line(&self, a: Vec2, b: Vec2, step: f32) -> Vec<Vec2> {
+		match self.which_subspace(a) {
+			Outer => vec![b],
+			Inner => {
+				let cs = InnerCS(self.tube);
+				let n = ((b - a).length() / step) as usize + 1;
+				let a = cs.global_to_flat(a);
+				let b = cs.global_to_flat(b);
+				(1..=n).map(|k| cs.flat_to_global(a.lerp(b, k as f32 / n as f32))).collect()
+			}
+			Boundary => panic!("Can't draw a line here!"),
+		}
+	}
+}
+
+struct Rect {
+	pub size: Vec2,
+}
+
+impl Rect {
+	/// Отражает луч, чтобы все координаты направления были положительны (допустимо благодаря симметрии Rect).
+	fn flip_ray(ray: Ray) -> Ray {
+		Ray { pos: ray.pos * ray.dir.signum(), dir: ray.dir.abs() }
+	}
+
+	fn is_inside(&self, pt: Vec2) -> bool {
+		pt.abs().cmplt(self.size).all()
+	}
+
+	fn trace_into(&self, ray: Ray) -> Option<f32> {
+		let ray = Self::flip_ray(ray);
+		// ray.pos.x + t * ray.dir.x = −size.x
+		let ts = (-self.size - ray.pos) / ray.dir;
+		let t = ts.max_element();
+		let pt = ray.pos + t * ray.dir;
+		if t < 0.0 { return None; }
+		if pt.cmpgt(self.size).any() { return None; }
+		Some(t)
+	}
+
+	fn trace_out_of(&self, ray: Ray) -> Option<f32> {
+		let ray = Self::flip_ray(ray);
+		// ray.pos.x + t * ray.dir.x = +size.x
+		let ts = (self.size - ray.pos) / ray.dir;
+		let t = ts.min_element();
+		Some(t)
+	}
+}
+
+#[test]
+fn test_rect() {
+	assert_eq!(Rect::flip_ray(Ray { pos: vec2(2.0, 3.0), dir: vec2(4.0, 5.0) }), Ray { pos: vec2(2.0, 3.0), dir: vec2(4.0, 5.0) });
+	assert_eq!(Rect::flip_ray(Ray { pos: vec2(2.0, 3.0), dir: vec2(-4.0, 5.0) }), Ray { pos: vec2(-2.0, 3.0), dir: vec2(4.0, 5.0) });
+	assert_eq!(Rect::flip_ray(Ray { pos: vec2(2.0, 3.0), dir: vec2(4.0, -5.0) }), Ray { pos: vec2(2.0, -3.0), dir: vec2(4.0, 5.0) });
+	assert_eq!(Rect::flip_ray(Ray { pos: vec2(2.0, 3.0), dir: vec2(4.0, 0.0) }), Ray { pos: vec2(2.0, 3.0), dir: vec2(4.0, 0.0) });
+
+	let r = Rect { size: vec2(2.0, 3.0) };
+
+	assert_eq!(r.trace_into(Ray { pos: vec2(3.0, 3.0), dir: vec2(1.0, 1.0) }), None);
+	assert_eq!(r.trace_into(Ray { pos: vec2(-3.0, 2.0), dir: vec2(1.0, 0.0) }), Some(1.0));
+	assert_eq!(r.trace_into(Ray { pos: vec2(-3.0, 2.0), dir: vec2(-1.0, 0.0) }), None);
+	assert_eq!(r.trace_into(Ray { pos: vec2(-3.0, 1.0), dir: vec2(2.0, 2.0) }), Some(0.5));
+	assert_eq!(r.trace_into(Ray { pos: vec2(-3.0, 2.1), dir: vec2(2.0, 2.0) }), None);
+
+	assert_eq!(r.trace_into(Ray { pos: vec2(2.0, 3.0), dir: vec2(1.0, 1.0) }), None);
+	assert_eq!(r.trace_into(Ray { pos: vec2(-2.0, 3.0), dir: vec2(-1.0, 1.0) }), None);
+	assert_eq!(r.trace_into(Ray { pos: vec2(2.0, 3.0), dir: vec2(-1.0, -1.0) }), Some(0.0));
+	assert_eq!(r.trace_into(Ray { pos: vec2(2.0, -3.0), dir: vec2(-1.0, 1.0) }), Some(0.0));
+
+	assert_eq!(r.trace_out_of(Ray { pos: vec2(0.0, 0.0), dir: vec2(1.0, 1.0) }), Some(2.0));
+	assert_eq!(r.trace_out_of(Ray { pos: vec2(0.0, 0.0), dir: vec2(0.0, 1.0) }), Some(3.0));
+	assert_eq!(r.trace_out_of(Ray { pos: vec2(0.0, 1.0), dir: vec2(0.0, -1.0) }), Some(4.0));
+	assert_eq!(r.trace_out_of(Ray { pos: vec2(1.0, 1.0), dir: vec2(0.0, -1.0) }), Some(4.0));
+	assert_eq!(r.trace_out_of(Ray { pos: vec2(2.0, 3.0), dir: vec2(1.0, 1.0) }), Some(0.0));
+}

src/bin/flat/types.rs

@@ -0,0 +1,50 @@
+use glam::{f32, i32, Mat2, Vec2};
+
+#[derive(Copy, Clone, Debug, PartialEq)]
+pub struct Ray {
+	pub pos: Vec2,
+	pub dir: Vec2,
+}
+
+impl Ray {
+	pub fn forward(&self, dist: f32) -> Ray {
+		Ray { pos: self.pos + self.dir * dist, dir: self.dir }
+	}
+}
+
+impl std::ops::Mul<Ray> for Mat2 {
+	type Output = Ray;
+
+	fn mul(self, rhs: Ray) -> Self::Output {
+		Ray { pos: self * rhs.pos, dir: self * rhs.dir }
+	}
+}
+
+#[derive(Copy, Clone, Debug, PartialEq)]
+pub struct Location {
+	/// Положение в основной СК
+	pub pos: Vec2,
+	/// Преобразование вектора из локальной ортонормированной в основную СК
+	pub rot: Mat2,
+}
+
+#[derive(Copy, Clone, Debug)]
+pub struct Object {
+	pub id: i32,
+	pub loc: Location,
+	pub r: f32,
+}
+
+#[derive(Copy, Clone, Debug)]
+pub struct Hit {
+	pub distance: f32,
+	pub id: i32,
+	pub pos: Vec2, // положение в основной СК
+	pub rel: Ray, // в локальной ортонормированной СК объекта
+}
+
+#[derive(Clone, Debug)]
+pub struct FlatTraceResult {
+	pub end: Option<Ray>,
+	pub objects: Vec<Hit>,
+}

src/bin/mesh/main.rs

@@ -0,0 +1,151 @@
+mod mesh_loader;
+
+use std::fs::File;
+use std::{env};
+use std::error::Error;
+use std::f32::consts::PI;
+use std::io::{BufReader};
+use glam::*;
+use show_image::{ImageInfo, ImageView, WindowOptions};
+use crate::mesh_loader::{Face, load_mesh};
+
+const W: i32 = 320;
+const H: i32 = 240;
+
+#[derive(Copy, Clone)]
+struct Color(u8, u8, u8);
+
+struct Image {
+	w: i32,
+	h: i32,
+	data: Vec<u8>,
+}
+
+impl Image {
+	fn data(&self) -> &[u8] {
+		self.data.as_slice()
+	}
+
+	fn put_pixel(&mut self, x: i32, y: i32, color: Color) {
+		if x < 0 || x >= self.w || y < 0 || y > self.h {
+			return;
+		}
+		let index = 3 * (x + self.w * y) as usize;
+		self.data[index] = color.0;
+		self.data[index + 1] = color.1;
+		self.data[index + 2] = color.2;
+	}
+}
+
+fn ypr_to_mat(ypr: Vec3) -> Mat3 {
+	let Vec3 { x: yaw, y: pitch, z: roll } = ypr;
+	let m_roll = mat3(
+		vec3(roll.cos(), roll.sin(), 0.0),
+		vec3(-roll.sin(), roll.cos(), 0.0),
+		vec3(0.0, 0.0, 1.0));
+	let m_yaw = mat3(
+		vec3(yaw.cos(), 0.0, yaw.sin()),
+		vec3(0.0, 1.0, 0.0),
+		vec3(-yaw.sin(), 0.0, yaw.cos()));
+	let m_pitch = mat3(
+		vec3(1.0, 0.0, 0.0),
+		vec3(0.0, pitch.cos(), -pitch.sin()),
+		vec3(0.0, pitch.sin(), pitch.cos()));
+	m_roll * m_pitch * m_yaw
+}
+
+type Mesh = [Face];
+
+struct TraceResult {
+	distance: f32,
+	normal: Vec3,
+}
+
+fn trace_to_mesh(mesh: &Mesh, base: Vec3, ray: Vec3) -> Option<TraceResult> {
+	let mut ret: Option<TraceResult> = None;
+	let mut dist = f32::INFINITY;
+	for f in mesh {
+		let fs = (0..3).map(|k| edge_dist(f.vertices[k], f.vertices[(k + 1) % 3], base, ray));
+		if fs.into_iter().all(|f| f >= 0.0) {
+			let m = mat3(f.vertices[1] - f.vertices[0], f.vertices[2] - f.vertices[0], -ray);
+			let m = m.inverse();
+			let rel = m * (base - f.vertices[0]);
+			if rel.z > dist {
+				continue;
+			}
+			dist = rel.z;
+			ret = Some(TraceResult {
+				distance: rel.z,
+				normal: f.normal,
+			});
+		}
+	}
+	ret
+}
+
+struct Location {
+	pos: Vec3,
+	rot: Vec4,
+}
+
+fn render(mesh: &Mesh, camera: impl Fn(Vec2) -> (Vec3, Vec3)) -> Image {
+	let bkg = vec3(0.0, 0.0, 0.0);
+	let mut img = Image {
+		w: W,
+		h: H,
+		data: vec![0; (3 * W * H) as usize],
+	};
+	let img_size = vec2(W as f32, H as f32);
+	for y in 0..H {
+		for x in 0..W {
+			let img_coords = vec2(x as f32, y as f32);
+			let off = (img_coords - img_size * 0.5) / img_size.y;
+			let (base, ray) = camera(off);
+			let color = if let Some(r) = trace_to_mesh(mesh, base, ray.normalize()) {
+				// to_vec3(0.45) * dot(r.normal, normalize(vec3(-1.0, 1.0, -1.0))) + 0.50
+				r.normal * 0.45 + 0.50
+			} else {
+				bkg
+			};
+			let color = (color * 255.0).as_ivec3().clamp(IVec3::splat(0), IVec3::splat(255));
+			img.put_pixel(x, y, Color(color.x as u8, color.y as u8, color.z as u8));
+		}
+	}
+	img
+}
+
+#[show_image::main]
+fn main() -> Result<(), Box<dyn Error>> {
+	let args: Vec<String> = env::args().collect();
+	let mesh = {
+		let f = File::open(&args[1])?;
+		let mut f = BufReader::new(f);
+		load_mesh(&mut f)?
+	};
+	let window = show_image::create_window("Raytracing", WindowOptions::default())?;
+	loop {
+		for phi in 0..360 {
+			let m_view = ypr_to_mat(vec3((135.0 + phi as f32) * PI / 180.0, -30.0 * PI / 180.0, 0.0f32));
+			let m_camera = m_view.transpose();
+			let img = render(mesh.as_slice(), |off| {
+				// perspective projection
+				let base = vec3(0.0, 0.0, -40.0);
+				let ray = vec3(off.x, off.y, 2.0);
+
+				// orthographic projection
+				// let base = vec3(off.x, off.y, -10.0);
+				// let ray = vec3(0.0, 0.0, 1.0);
+
+				(m_camera * base, m_camera * ray)
+			});
+
+			let image = ImageView::new(ImageInfo::rgb8(W as u32, H as u32), img.data());
+			window.set_image("image", image)?;
+		}
+	}
+}
+
+fn edge_dist(a: Vec3, b: Vec3, base: Vec3, dir: Vec3) -> f32 {
+	// Note: given that the input is not arbitrary but comes from a cartesian product of certain (a, b) pairs and certain (base, dir) pairs, this can be optimized from Cnm to an+bm+cnm with c<C.
+	mat3(b - a, base - a, -dir).determinant()
+}

src/bin/mesh/mesh_loader.rs

@@ -1,5 +1,5 @@
 use std::io;
-use glm::{vec2, vec3, Vec2, Vec3};
+use glam::{vec2, vec3, Vec2, Vec3};
 
 #[derive(Copy, Clone, Debug)]
 struct ObjVertex {

src/main.rs

@@ -1,127 +0,0 @@
-mod mesh_loader;
-
-use std::fs::File;
-use std::{env};
-use std::error::Error;
-use std::f32::consts::PI;
-use std::io::{BufRead, BufReader};
-use std::io::Write;
-use rand::Rng;
-use glm::*;
-use show_image::{ImageInfo, ImageView, WindowOptions};
-use crate::mesh_loader::{Face, load_mesh};
-
-const W: i32 = 320;
-const H: i32 = 240;
-const SCALE: f32 = 30.0;
-
-#[derive(Copy, Clone)]
-struct Color(u8, u8, u8);
-
-struct Image {
-	w: i32,
-	h: i32,
-	data: Vec<u8>,
-}
-
-impl Image {
-	fn data(&self) -> &[u8] {
-		self.data.as_slice()
-	}
-
-	fn put_pixel(&mut self, x: i32, y: i32, color: Color) {
-		if x < 0 || x >= self.w || y < 0 || y > self.h {
-			return;
-		}
-		let index = 3 * (x + self.w * y) as usize;
-		self.data[index] = color.0;
-		self.data[index + 1] = color.1;
-		self.data[index + 2] = color.2;
-	}
-}
-
-fn ypr_to_mat(ypr: Vec3) -> Mat3 {
-	let Vec3 { x: yaw, y: pitch, z: roll } = ypr;
-	let m_roll = mat3(
-		roll.cos(), roll.sin(), 0.0,
-		-roll.sin(), roll.cos(), 0.0,
-		0.0, 0.0, 1.0);
-	let m_yaw = mat3(
-		yaw.cos(), 0.0, yaw.sin(),
-		0.0, 1.0, 0.0,
-		-yaw.sin(), 0.0, yaw.cos());
-	let m_pitch = mat3(
-		1.0, 0.0, 0.0,
-		0.0, pitch.cos(), -pitch.sin(),
-		0.0, pitch.sin(), pitch.cos());
-	m_roll * m_pitch * m_yaw
-}
-
-fn render(mesh: &[Face], camera: impl Fn(Vec2) -> (Vec3, Vec3)) -> Image {
-	let mut img = Image {
-		w: W,
-		h: H,
-		data: vec![0; (3 * W * H) as usize],
-	};
-	let img_size = vec2(W as f32, H as f32);
-	for y in 0..H {
-		for x in 0..W {
-			let img_coords = vec2(x as f32, y as f32);
-			let off = (img_coords - img_size * 0.5) / img_size.y;
-			let (base, ray) = camera(off);
-			let mut dist = f32::INFINITY;
-			for f in mesh {
-				let color = clamp(to_ivec3(f.normal * 120.0 + 128.0), ivec3(0, 0, 0), ivec3(255, 255, 255));
-				let fs = (0..3).map(|k| edge_dist(f.vertices[k], f.vertices[(k + 1) % 3], base, ray));
-				if fs.into_iter().all(|f| f >= 0.0) {
-					let m = Mat3 { c0: f.vertices[1] - f.vertices[0], c1: f.vertices[2] - f.vertices[0], c2: -ray };
-					if let Some(m) = m.inverse() {
-						let rel = m * (base - f.vertices[0]);
-						if rel.z > dist {
-							continue;
-						}
-						dist = rel.z;
-					} else {
-						continue;
-					}
-					img.put_pixel(x, y, Color(color.x as u8, color.y as u8, color.z as u8));
-				}
-			}
-		}
-	}
-	img
-}
-
-#[show_image::main]
-fn main() -> Result<(), Box<dyn Error>> {
-	let args: Vec<String> = env::args().collect();
-	let mesh = {
-		let f = File::open(&args[1])?;
-		let mut f = BufReader::new(f);
-		load_mesh(&mut f)?
-	};
-	let m_view = ypr_to_mat(vec3(135.0 * PI / 180.0, -30.0 * PI / 180.0, 0.0f32));
-	let m_camera = transpose(&m_view);
-	let img = render(mesh.as_slice(), |off| {
-		// perspective projection
-		let base = vec3(0.0, 0.0, -40.0);
-		let ray = vec3(off.x, off.y, 2.0);
-
-		// orthographic projection
-		// let base = vec3(off.x, off.y, -10.0);
-		// let ray = vec3(0.0, 0.0, 1.0);
-
-		(m_camera * base, m_camera * ray)
-	});
-
-	let window = show_image::create_window("Raytracing", WindowOptions::default())?;
-	let image = ImageView::new(ImageInfo::rgb8(W as u32, H as u32), img.data());
-	window.set_image("image", image)?;
-	window.wait_until_destroyed()?;
-	Ok(())
-}
-
-fn edge_dist(a: Vec3, b: Vec3, base: Vec3, dir: Vec3) -> f32 {
-	// Note: given that the input is not arbitrary but comes from a cartesian product of certain (a, b) pairs and certain (base, dir) pairs, this can be optimized from Cnm to an+bm+cnm with c<C.
-	Mat3 { c0: b - a, c1: base - a, c2: -dir }.determinant()
-}