Support 3D!

Yes this is code duplication. But glam is not dimension-generic.

src/bin/flat/main.rs

@@ -4,24 +4,15 @@ use flo_canvas::*;
 use flo_draw::*;
 use glam::*;
 
-use crate::ifaces::{DebugTraceable, Traceable};
-use crate::types::FlatTraceResult;
+use refraction::ifaces::{DebugTraceable, Traceable};
 use refraction::mathx::MatExt;
-use riemann::{trace_iter, Metric};
-use tube::metric::Tube;
-use tube::Space;
-use tube::Subspace::{Boundary, Inner, Outer};
-use types::{Location, Object, Ray};
+use refraction::riemann::{trace_iter, Metric};
+use refraction::tube::metric::Tube;
+use refraction::tube::Space;
+use refraction::types::{Location, Object, Ray};
+use refraction::DT;
 
-mod fns;
-mod ifaces;
-mod riemann;
-mod tube;
-mod types;
-
-const DT: f32 = 0.1;
-
-fn draw_loop(gc: &mut Vec<Draw>, mut pts: impl Iterator<Item = Vec2>) {
+fn draw_loop(gc: &mut Vec<Draw>, mut pts: impl Iterator<Item = Vec3>) {
 	gc.new_path();
 	let Some(first) = pts.next() else {
 		return;
@@ -52,23 +43,31 @@ pub fn main() {
 					id: k as i32,
 					loc: put_object(
 						&tube,
-						vec2(0.0, y * tube.external_halflength),
-						Mat2::from_angle(y),
+						vec3(0.0, y * tube.external_halflength, 0.0),
+						Mat3::from_mat2(Mat2::from_angle(y)),
 					),
 					r: 20.0,
 				})
 				.collect();
 			let space = Space { tube, objs };
-			let cam1 = put_object(&space.tube, vec2(-500., 0.), Mat2::IDENTITY);
+			let cam1 = put_object(&space.tube, vec3(-500., 0., 0.), Mat3::IDENTITY);
 			let cam2 = put_object(
 				&space.tube,
-				vec2(-2.5 * tube.outer_radius, 1.25 * tube.external_halflength),
-				mat2(vec2(1., -1.), vec2(1., 1.)),
+				vec3(
+					-2.5 * tube.outer_radius,
+					1.25 * tube.external_halflength,
+					0.,
+				),
+				mat3(vec3(1., -1., 0.), vec3(1., 1., 0.), vec3(0., 0., 1.)),
 			);
 			let cam3 = put_object(
 				&space.tube,
-				vec2(0.25 * tube.inner_radius, 0.25 * tube.external_halflength),
-				mat2(vec2(0., -1.), vec2(1., 0.)),
+				vec3(
+					0.25 * tube.inner_radius,
+					0.25 * tube.external_halflength,
+					0.,
+				),
+				mat3(vec3(0., -1., 0.), vec3(1., 0., 0.), vec3(0., 0., 1.)),
 			);
 
 			gc.canvas_height(500.0);
@@ -108,6 +107,7 @@ pub fn main() {
 						.skip(1)
 						.map(|φ| {
 							let dir = Vec2::from_angle(φ) * obj.r;
+							let dir = vec3(dir.x, dir.y, 0.);
 							let dir = obj.loc.rot * dir;
 							pos + dir
 						}),
@@ -119,6 +119,7 @@ pub fn main() {
 						.skip(1)
 						.map(|φ| {
 							let dir = Vec2::from_angle(φ) * obj.r;
+							let dir = vec3(dir.x, dir.y, 0.);
 							let dir = obj.loc.rot * dir;
 							space.trace_step(Ray { pos, dir }).pos
 						}),
@@ -132,6 +133,7 @@ pub fn main() {
 							let n = obj.r.floor();
 							let d = obj.r / n;
 							let dir = Vec2::from_angle(φ);
+							let dir = vec3(dir.x, dir.y, 0.);
 							let dir = obj.loc.rot * dir * d;
 							space
 								.trace_iter(Ray { pos, dir })
@@ -145,7 +147,7 @@ pub fn main() {
 	});
 }
 
-fn rel_to_abs(space: &impl Metric, base: &Location, rel: Vec2, steps: usize) -> Vec2 {
+fn rel_to_abs(space: &impl Metric, base: &Location, rel: Vec3, steps: usize) -> Vec3 {
 	let c = 1.0 / (steps as f32);
 	trace_iter(space, base.pos, base.rot * rel, c * rel.length())
 		.nth(steps - 1)
@@ -153,7 +155,7 @@ fn rel_to_abs(space: &impl Metric, base: &Location, rel: Vec2, steps: usize) ->
 }
 
 /// Converts a position and a rotation to a [Location]. Only the X direction is preserved from `rot` to ensure the resulting Location describes an orthonormal coordinate system.
-fn put_object(space: &impl Metric, pos: Vec2, rot: Mat2) -> Location {
+fn put_object(space: &impl Metric, pos: Vec3, rot: Mat3) -> Location {
 	let metric_sqrt = space.sqrt_at(pos);
 	let metric_inv_sqrt = space.sqrt_at(pos).inverse();
 	let rot = metric_inv_sqrt * (metric_sqrt * rot).orthonormalize();
@@ -165,27 +167,59 @@ fn test_put_object() {
 	use approx::assert_abs_diff_eq;
 
 	let ε = 1e-5;
-	let m = riemann::samples::ScaledMetric {
-		scale: vec2(3., 4.),
+	let m = refraction::riemann::samples::ScaledMetric {
+		scale: vec3(3., 4., 5.),
 	};
 
-	let loc = put_object(&m, vec2(1., 2.), mat2(vec2(1., 0.), vec2(0., 1.)));
-	assert_eq!(loc.pos, vec2(1., 2.));
-	assert_abs_diff_eq!(loc.rot * vec2(1., 0.), vec2(1. / 3., 0.), epsilon = ε);
-	assert_abs_diff_eq!(loc.rot * vec2(0., 1.), vec2(0., 1. / 4.), epsilon = ε);
+	let loc = put_object(
+		&m,
+		vec3(1., 2., 0.),
+		mat3(vec3(1., 0., 0.), vec3(0., 1., 0.), vec3(0., 0., 1.)),
+	);
+	assert_eq!(loc.pos, vec3(1., 2., 0.));
+	assert_abs_diff_eq!(
+		loc.rot * vec3(1., 0., 0.),
+		vec3(1. / 3., 0., 0.),
+		epsilon = ε
+	);
+	assert_abs_diff_eq!(
+		loc.rot * vec3(0., 1., 0.),
+		vec3(0., 1. / 4., 0.),
+		epsilon = ε
+	);
 
-	let loc = put_object(&m, vec2(1., 2.), mat2(vec2(0., 1.), vec2(-1., 0.)));
-	assert_eq!(loc.pos, vec2(1., 2.));
-	assert_abs_diff_eq!(loc.rot * vec2(1., 0.), vec2(0., 1. / 4.), epsilon = ε);
-	assert_abs_diff_eq!(loc.rot * vec2(0., 1.), vec2(-1. / 3., 0.), epsilon = ε);
+	let loc = put_object(
+		&m,
+		vec3(1., 2., 0.),
+		mat3(vec3(0., 1., 0.), vec3(-1., 0., 0.), vec3(0., 0., 1.)),
+	);
+	assert_eq!(loc.pos, vec3(1., 2., 0.));
+	assert_abs_diff_eq!(
+		loc.rot * vec3(1., 0., 0.),
+		vec3(0., 1. / 4., 0.),
+		epsilon = ε
+	);
+	assert_abs_diff_eq!(
+		loc.rot * vec3(0., 1., 0.),
+		vec3(-1. / 3., 0., 0.),
+		epsilon = ε
+	);
 
 	let c = 0.5 * std::f32::consts::SQRT_2;
-	let loc = put_object(&m, vec2(1., 2.), mat2(vec2(c, c), vec2(-c, c)));
-	assert_eq!(loc.pos, vec2(1., 2.));
-	assert_abs_diff_eq!(loc.rot * vec2(1., 0.), vec2(1. / 5., 1. / 5.), epsilon = ε);
+	let loc = put_object(
+		&m,
+		vec3(1., 2., 0.),
+		mat3(vec3(c, c, 0.), vec3(-c, c, 0.), vec3(0., 0., 1.)),
+	);
+	assert_eq!(loc.pos, vec3(1., 2., 0.));
 	assert_abs_diff_eq!(
-		loc.rot * vec2(0., 1.),
-		vec2(-4. / 15., 3. / 20.),
+		loc.rot * vec3(1., 0., 0.),
+		vec3(1. / 5., 1. / 5., 0.),
+		epsilon = ε
+	);
+	assert_abs_diff_eq!(
+		loc.rot * vec3(0., 1., 0.),
+		vec3(-4. / 15., 3. / 20., 0.),
 		epsilon = ε
 	);
 }
@@ -197,8 +231,8 @@ fn draw_cross(gc: &mut Vec<Draw>, pos: Vec2, r: f32) {
 	gc.line_to(pos.x + r, pos.y - r);
 }
 
-fn draw_ray_2(gc: &mut Vec<Draw>, space: &Space, camera: Location, dir: Vec2) {
-	let pos = vec2(0., 0.);
+fn draw_ray_2(gc: &mut Vec<Draw>, space: &Space, camera: Location, dir: Vec3) {
+	let pos = vec3(0., 0., 0.);
 	let (hits, path) = space.trace_dbg(camera, Ray { pos, dir });
 	let hits2 = space.trace(camera, Ray { pos, dir });
 	for (a, b) in hits.into_iter().zip(hits2.into_iter()) {
@@ -223,7 +257,7 @@ fn draw_ray_2(gc: &mut Vec<Draw>, space: &Space, camera: Location, dir: Vec2) {
 
 fn draw_fan_2(gc: &mut Vec<Draw>, space: &Space, camera: Location, spread: f32) {
 	for y in itertools_num::linspace(-spread, spread, 101) {
-		draw_ray_2(gc, space, camera, vec2(1., y));
+		draw_ray_2(gc, space, camera, vec3(1., y, 0.));
 	}
 }
 
@@ -234,10 +268,14 @@ fn draw_track(gc: &mut Vec<Draw>, space: &Space, start: Vec2, dir: Vec2) {
 	// 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)),
+		pos: vec3(start.x, start.y, 0.),
+		rot: mat3(
+			vec3(dir.x, dir.y, 0.),
+			vec3(-dir.y, dir.x, 0.),
+			vec3(0., 0., 1.),
+		),
 	};
-	let v = vec2(1.0, 0.0);
+	let v = vec3(1., 0., 0.);
 	let mut draw = |loc: &Location| {
 		let p = loc.pos;
 		let ax = p + loc.rot.x_axis * SCALE;

src/bin/flat/riemann.rs

@@ -1,281 +0,0 @@
-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<T> std::ops::Mul<T> for Decomp2
-where
-	Mat2: std::ops::Mul<T>,
-{
-	type Output = <Mat2 as std::ops::Mul<T>>::Output;
-
-	fn mul(self, rhs: T) -> Self::Output {
-		Mat2::from(self) * rhs
-	}
-}
-
-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)
-	}
-}
-
-pub struct TraceIter<'a, M: Metric> {
-	space: &'a M,
-	p: Vec2,
-	v: Vec2,
-}
-
-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 += a;
-		self.p += self.v;
-		Some(self.p)
-	}
-}
-
-pub fn trace_iter<M: Metric>(space: &M, base: Vec2, dir: Vec2, dt: f32) -> TraceIter<M> {
-	TraceIter {
-		space,
-		p: base,
-		v: dt * space.normalize_vec_at(base, dir),
-	}
-}
-
-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 {
-		/// Specifies unit size in each cardinal direction. E.g. with scale=(2, 3), vector (1, 0) has length 2 while a unit vector with the same direction is (1/2, 0).
-		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 super::*;
-	use approx::assert_abs_diff_eq;
-
-	use glam::{mat2, vec2, Mat2};
-	use rand::{Rng, SeedableRng};
-
-	#[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.)
-		);
-	}
-
-	#[test]
-	fn test_trace_iter() {
-		let metric = samples::ScaledMetric {
-			scale: vec2(2., 4.),
-		};
-		assert_eq!(
-			trace_iter(&metric, vec2(3., 5.), vec2(1., 0.), 1.)
-				.nth(7)
-				.unwrap(),
-			vec2(7., 5.)
-		);
-		assert_eq!(
-			trace_iter(&metric, vec2(3., 5.), vec2(2., 0.), 1.)
-				.nth(7)
-				.unwrap(),
-			vec2(7., 5.)
-		);
-		assert_eq!(
-			trace_iter(&metric, vec2(3., 5.), vec2(1., 0.), 0.5)
-				.nth(7)
-				.unwrap(),
-			vec2(5., 5.)
-		);
-		assert_eq!(
-			trace_iter(&metric, vec2(3., 5.), vec2(0., 1.), 1.)
-				.nth(9)
-				.unwrap(),
-			vec2(3., 7.5)
-		);
-		assert_eq!(
-			trace_iter(&metric, vec2(3., 5.), vec2(0., 4.), 1.)
-				.nth(9)
-				.unwrap(),
-			vec2(3., 7.5)
-		);
-		assert_eq!(
-			trace_iter(&metric, vec2(3., 5.), vec2(0., 1.), 0.5)
-				.nth(9)
-				.unwrap(),
-			vec2(3., 6.25)
-		);
-		assert_abs_diff_eq!(
-			trace_iter(
-				&metric,
-				vec2(3., 5.),
-				vec2(0.5, 0.25),
-				std::f32::consts::SQRT_2
-			)
-			.nth(7)
-			.unwrap(),
-			vec2(7., 7.),
-			epsilon = 1e-5
-		);
-	}
-}

src/bin/flat/tube/coords.rs

@@ -1,479 +0,0 @@
-use glam::{vec2, Mat2, 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::{assert_abs_diff_eq, AbsDiffEq};
-	use glam::{mat2, vec2, Mat2, 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),
-	) {
-		#[allow(non_upper_case_globals)]
-		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/fns.rs

@@ -1,4 +1,4 @@
-use refraction::mathx::FloatExt2;
+use crate::mathx::FloatExt2;
 
 pub trait Limiter {
 	fn value(&self, x: f32) -> f32;

src/ifaces.rs

@@ -1,5 +1,5 @@
 use crate::types::{Hit, Location, Ray};
-use glam::Vec2;
+use glam::Vec3;
 
 pub trait Traceable {
 	/// Traces a ray from a given starting point. `ray` is relative to the camera.
@@ -19,8 +19,8 @@ pub trait OptimizedTraceable: Traceable {
 }
 
 pub struct RayPath {
-	pub points: Vec<Vec2>,
-	pub end_dir: Vec2,
+	pub points: Vec<Vec3>,
+	pub end_dir: Vec3,
 }
 
 pub trait DebugTraceable: Traceable {

src/lib.rs

@@ -1,3 +1,10 @@
+mod fns;
+pub mod ifaces;
 pub mod mathx;
 pub mod mesh_loader;
 pub mod mesh_tracer;
+pub mod riemann;
+pub mod tube;
+pub mod types;
+
+pub const DT: f32 = 0.1;

src/mathx.rs

@@ -1,4 +1,4 @@
-use glam::{FloatExt, Mat2, Mat3};
+use glam::{FloatExt, Mat2, Mat3, Vec2, Vec3};
 
 mod bounds {
 	pub trait Pair<T> {}
@@ -44,6 +44,104 @@ impl MatExt for Mat3 {
 	}
 }
 
+/// Represents a 2×2 matrix decomposed as O^T D O, where O is orthogonal and D is diagonal.
+///
+/// Not every matrix can be decomposed like this, only that of a symmetric bilinear function.
+#[derive(Copy, Clone, Debug, PartialEq)]
+pub struct Decomp2 {
+	/// The orthogonal part.
+	///
+	/// Using a non-orthogonal matrix will yield to incorrect results (but no UB).
+	pub ortho: Mat2,
+
+	/// The diagonal part.
+	pub diag: Vec2,
+}
+
+impl Decomp2 {
+	/// Computes the square of this, more efficiently than doing that with a matrix.
+	pub fn square(&self) -> Self {
+		Self {
+			ortho: self.ortho,
+			diag: self.diag * self.diag,
+		}
+	}
+
+	/// Computes the inverse of this, more efficiently than doing that with a matrix.
+	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
+	}
+}
+
+impl<T> std::ops::Mul<T> for Decomp2
+where
+	Mat2: std::ops::Mul<T>,
+{
+	type Output = <Mat2 as std::ops::Mul<T>>::Output;
+
+	fn mul(self, rhs: T) -> Self::Output {
+		Mat2::from(self) * rhs
+	}
+}
+
+/// Represents a 3×3 matrix decomposed as O^T D O, where O is orthogonal and D is diagonal.
+///
+/// Not every matrix can be decomposed like this, only that of a symmetric bilinear function.
+#[derive(Copy, Clone, Debug, PartialEq)]
+pub struct Decomp3 {
+	/// The orthogonal part.
+	///
+	/// Using a non-orthogonal matrix will yield to incorrect results (but no UB).
+	pub ortho: Mat3,
+
+	/// The diagonal part.
+	pub diag: Vec3,
+}
+
+impl Decomp3 {
+	/// Computes the square of this, more efficiently than doing that with a matrix.
+	pub fn square(&self) -> Self {
+		Self {
+			ortho: self.ortho,
+			diag: self.diag * self.diag,
+		}
+	}
+
+	/// Computes the inverse of this, more efficiently than doing that with a matrix.
+	pub fn inverse(&self) -> Self {
+		Self {
+			ortho: self.ortho,
+			diag: Vec3::splat(1.0) / self.diag,
+		}
+	}
+}
+
+impl From<Decomp3> for Mat3 {
+	fn from(value: Decomp3) -> Self {
+		value.ortho.transpose() * Mat3::from_diagonal(value.diag) * value.ortho
+	}
+}
+
+impl<T> std::ops::Mul<T> for Decomp3
+where
+	Mat3: std::ops::Mul<T>,
+{
+	type Output = <Mat3 as std::ops::Mul<T>>::Output;
+
+	fn mul(self, rhs: T) -> Self::Output {
+		Mat3::from(self) * rhs
+	}
+}
+
 #[cfg(test)]
 mod tests {
 	use super::*;

src/riemann.rs

@@ -0,0 +1,248 @@
+use crate::mathx::Decomp3;
+use glam::*;
+
+pub type Tens3 = [Mat3; 3];
+
+pub trait Metric {
+	fn sqrt_at(&self, pos: Vec3) -> Decomp3;
+
+	fn at(&self, pos: Vec3) -> Mat3 {
+		self.sqrt_at(pos).square().into()
+	}
+
+	fn inverse_at(&self, pos: Vec3) -> Mat3 {
+		self.sqrt_at(pos).square().inverse().into()
+	}
+
+	fn part_derivs_at(&self, pos: Vec3) -> Tens3 {
+		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: Vec3, v: Vec3) -> f32 {
+		v.dot(self.at(at) * v).sqrt()
+	}
+
+	fn normalize_vec_at(&self, at: Vec3, v: Vec3) -> Vec3 {
+		v / self.vec_length_at(at, v)
+	}
+}
+
+pub struct TraceIter<'a, M: Metric> {
+	space: &'a M,
+	p: Vec3,
+	v: Vec3,
+}
+
+impl<'a, M: Metric> Iterator for TraceIter<'a, M> {
+	type Item = Vec3;
+
+	fn next(&mut self) -> Option<Self::Item> {
+		let a: Vec3 = -contract2(krist(self.space, self.p), self.v);
+		self.v += a;
+		self.p += self.v;
+		Some(self.p)
+	}
+}
+
+pub fn trace_iter<M: Metric>(space: &M, base: Vec3, dir: Vec3, dt: f32) -> TraceIter<M> {
+	TraceIter {
+		space,
+		p: base,
+		v: dt * space.normalize_vec_at(base, dir),
+	}
+}
+
+pub fn krist(space: &impl Metric, pos: Vec3) -> Tens3 {
+	// Γ^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_tens3(|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(Vec3) -> Mat3, pos: Vec3, delta: Vec3) -> Mat3 {
+	(f(pos + delta) - f(pos - delta)) / (2.0 * delta.length())
+}
+
+fn part_deriv(f: impl Fn(Vec3) -> Mat3, pos: Vec3, eps: f32) -> Tens3 {
+	[
+		dir_deriv(&f, pos, vec3(eps, 0.0, 0.0)),
+		dir_deriv(&f, pos, vec3(0.0, eps, 0.0)),
+		dir_deriv(&f, pos, vec3(0.0, 0.0, eps)),
+	]
+}
+
+/// Сворачивает тензор t с вектором u
+pub fn contract(t: Tens3, u: Vec3) -> Mat3 {
+	mat3(t[0] * u, t[1] * u, t[2] * u).transpose()
+}
+
+/// Сворачивает тензор t с вектором v дважды, по второму и третьему индексам.
+pub fn contract2(t: Tens3, v: Vec3) -> Vec3 {
+	contract(t, v) * v
+}
+
+fn make_vec3(f: impl Fn(usize) -> f32) -> Vec3 {
+	Vec3::from_array(std::array::from_fn(|i| f(i)))
+}
+
+fn make_mat3(f: impl Fn(usize, usize) -> f32) -> Mat3 {
+	Mat3::from_cols_array_2d(&std::array::from_fn(|i| std::array::from_fn(|j| f(i, j))))
+}
+
+fn make_tens3(f: impl Fn(usize, usize, usize) -> f32) -> Tens3 {
+	std::array::from_fn(|i| make_mat3(|j, k| f(i, j, k)))
+}
+
+#[test]
+fn m3() {
+	let m = make_mat3(|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 t3() {
+	let t = make_tens3(|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::{Mat3, Vec3};
+
+	use super::{Decomp3, Metric};
+
+	pub struct ScaledMetric {
+		/// Specifies unit size in each cardinal direction. E.g. with scale=(2, 3), vector (1, 0) has length 2 while a unit vector with the same direction is (1/2, 0).
+		pub scale: Vec3,
+	}
+
+	impl Metric for ScaledMetric {
+		fn sqrt_at(&self, _pos: Vec3) -> Decomp3 {
+			Decomp3 {
+				diag: self.scale,
+				ortho: Mat3::IDENTITY,
+			}
+		}
+	}
+}
+
+#[cfg(test)]
+mod tests {
+	use super::*;
+	use approx::assert_abs_diff_eq;
+
+	use glam::{mat3, vec3, Mat3};
+	use rand::{Rng, SeedableRng};
+
+	#[test]
+	fn uniform_scaled_metric() {
+		let mut rng = rand_pcg::Pcg64Mcg::seed_from_u64(17);
+		let metric = samples::ScaledMetric {
+			scale: vec3(3., 4., 5.),
+		};
+		assert_eq!(
+			metric.sqrt_at(rng.gen()),
+			Decomp3 {
+				ortho: Mat3::IDENTITY,
+				diag: vec3(3., 4., 5.)
+			}
+		);
+		assert_eq!(
+			metric.at(rng.gen()),
+			Mat3::from_cols_array(&[9., 0., 0., 0., 16., 0., 0., 0., 25.])
+		);
+		assert_eq!(
+			metric.inverse_at(rng.gen()),
+			Mat3::from_cols_array(&[1. / 9., 0., 0., 0., 1. / 16., 0., 0., 0., 1. / 25.])
+		);
+		assert_eq!(
+			metric.part_derivs_at(rng.gen()),
+			[Mat3::ZERO, Mat3::ZERO, Mat3::ZERO]
+		);
+		assert_eq!(metric.vec_length_at(rng.gen(), vec3(1., 0., 0.)), 3.);
+		assert_eq!(metric.vec_length_at(rng.gen(), vec3(0., 1., 0.)), 4.);
+		assert_eq!(metric.vec_length_at(rng.gen(), vec3(0., 0., 1.)), 5.);
+		assert_eq!(metric.vec_length_at(rng.gen(), vec3(1., 1., 0.)), 5.);
+		assert_eq!(
+			metric.normalize_vec_at(rng.gen(), vec3(1., 0., 0.)),
+			vec3(1. / 3., 0., 0.)
+		);
+		assert_eq!(
+			metric.normalize_vec_at(rng.gen(), vec3(0., 1., 0.)),
+			vec3(0., 1. / 4., 0.)
+		);
+		assert_eq!(
+			metric.normalize_vec_at(rng.gen(), vec3(1., 1., 0.)),
+			vec3(1. / 5., 1. / 5., 0.)
+		);
+	}
+
+	#[test]
+	fn test_trace_iter() {
+		let metric = samples::ScaledMetric {
+			scale: vec3(2., 4., 3.),
+		};
+		assert_eq!(
+			trace_iter(&metric, vec3(3., 5., 0.), vec3(1., 0., 0.), 1.)
+				.nth(7)
+				.unwrap(),
+			vec3(7., 5., 0.)
+		);
+		assert_eq!(
+			trace_iter(&metric, vec3(3., 5., 0.), vec3(2., 0., 0.), 1.)
+				.nth(7)
+				.unwrap(),
+			vec3(7., 5., 0.)
+		);
+		assert_eq!(
+			trace_iter(&metric, vec3(3., 5., 0.), vec3(1., 0., 0.), 0.5)
+				.nth(7)
+				.unwrap(),
+			vec3(5., 5., 0.)
+		);
+		assert_eq!(
+			trace_iter(&metric, vec3(3., 5., 0.), vec3(0., 1., 0.), 1.)
+				.nth(9)
+				.unwrap(),
+			vec3(3., 7.5, 0.)
+		);
+		assert_eq!(
+			trace_iter(&metric, vec3(3., 5., 0.), vec3(0., 4., 0.), 1.)
+				.nth(9)
+				.unwrap(),
+			vec3(3., 7.5, 0.)
+		);
+		assert_eq!(
+			trace_iter(&metric, vec3(3., 5., 0.), vec3(0., 1., 0.), 0.5)
+				.nth(9)
+				.unwrap(),
+			vec3(3., 6.25, 0.)
+		);
+		assert_abs_diff_eq!(
+			trace_iter(
+				&metric,
+				vec3(3., 5., 0.),
+				vec3(0.5, 0.25, 0.),
+				std::f32::consts::SQRT_2
+			)
+			.nth(7)
+			.unwrap(),
+			vec3(7., 7., 0.),
+			epsilon = 1e-5
+		);
+	}
+}

src/tube/coords.rs

@@ -0,0 +1,513 @@
+use glam::{vec3, Mat3, Vec3};
+
+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<Vec3> + FlatCoordinateSystem<Ray> + FlatCoordinateSystem<Location>
+{
+	// Измеряет расстояние до выхода за пределы области вдоль луча ray. Луч задаётся в плоской СК.
+	fn distance_to_boundary(&self, _ray: Ray) -> Option<f32> {
+		None
+	}
+}
+
+trait MetricCS: FlatCoordinateSystem<Vec3> {
+	fn global_metric(&self) -> &impl Metric;
+	fn flat_to_global_tfm_at(&self, pos: Vec3) -> Mat3 {
+		self.global_metric()
+			.sqrt_at(self.flat_to_global(pos))
+			.inverse()
+			.into()
+	}
+	fn global_to_flat_tfm_at(&self, pos: Vec3) -> Mat3 {
+		self.global_metric().sqrt_at(pos).into()
+	}
+}
+
+impl<T: FlatCoordinateSystem<Vec3> + 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<Vec3> + 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<Vec3> for InnerCS {
+	fn flat_to_global(&self, pos: Vec3) -> Vec3 {
+		vec3(pos.x, self.0.y(pos.y), pos.z)
+	}
+
+	// Работает только при |pos.x| ≤ inner_radius или |pos.y| ≥ external_halflength.
+	fn global_to_flat(&self, pos: Vec3) -> Vec3 {
+		vec3(pos.x, self.0.v(pos.y), pos.z)
+	}
+}
+
+impl FlatRegion for InnerCS {
+	fn distance_to_boundary(&self, ray: Ray) -> Option<f32> {
+		Rect {
+			size: vec3(
+				self.0.inner_radius,
+				self.0.internal_halflength,
+				self.0.inner_radius,
+			),
+		}
+		.trace_out_of(ray)
+	}
+}
+
+pub struct OuterCS(pub Tube);
+
+impl MetricCS for OuterCS {
+	fn global_metric(&self) -> &impl Metric {
+		&self.0
+	}
+}
+
+impl FlatCoordinateSystem<Vec3> for OuterCS {
+	fn flat_to_global(&self, pos: Vec3) -> Vec3 {
+		let inner = Rect {
+			size: vec3(
+				self.0.inner_radius + 1.0,
+				self.0.external_halflength,
+				self.0.inner_radius + 1.0,
+			),
+		};
+		if inner.is_inside(pos) {
+			let Vec3 { x, y: v, z } = pos;
+			let y = self
+				.0
+				.y(v - v.signum() * (self.0.external_halflength - self.0.internal_halflength));
+			vec3(x, y, z)
+		} else {
+			pos
+		}
+	}
+
+	fn global_to_flat(&self, pos: Vec3) -> Vec3 {
+		let inner = Rect {
+			size: vec3(
+				self.0.inner_radius + 1.0,
+				self.0.external_halflength,
+				self.0.inner_radius + 1.0,
+			),
+		};
+		if inner.is_inside(pos) {
+			let Vec3 { x: u, y, z: w } = pos; // в основной СК
+			let v = self.0.v(y)
+				+ y.signum() * (self.0.external_halflength - self.0.internal_halflength);
+			vec3(u, v, w) // в плоском продолжении СК Outer на область Inner
+		} else {
+			pos
+		}
+	}
+}
+
+impl FlatRegion for OuterCS {
+	fn distance_to_boundary(&self, ray: Ray) -> Option<f32> {
+		Rect {
+			size: vec3(
+				self.0.outer_radius,
+				self.0.external_halflength,
+				self.0.outer_radius,
+			),
+		}
+		.trace_into(ray)
+	}
+}
+
+#[cfg(test)]
+mod tests {
+	use approx::{assert_abs_diff_eq, AbsDiffEq};
+	use glam::{mat3, vec3, Mat3, Vec3};
+	use itertools_num::linspace;
+
+	use crate::riemann::samples;
+
+	use super::*;
+
+	#[test]
+	fn uniform_scaled_metric() {
+		struct Scaled(samples::ScaledMetric, Vec3);
+		impl FlatCoordinateSystem<Vec3> for Scaled {
+			fn flat_to_global(&self, pos: Vec3) -> Vec3 {
+				(pos - self.1) / self.0.scale
+			}
+			fn global_to_flat(&self, pos: Vec3) -> Vec3 {
+				pos * self.0.scale + self.1
+			}
+		}
+		impl MetricCS for Scaled {
+			fn global_metric(&self) -> &impl Metric {
+				&self.0
+			}
+		}
+		let cs = Scaled(
+			samples::ScaledMetric {
+				scale: vec3(3., 4., 5.),
+			},
+			vec3(2., 3., 7.),
+		);
+		assert_eq!(cs.global_to_flat(vec3(7., 3., 1.)), vec3(23., 15., 12.));
+		assert_eq!(cs.flat_to_global(vec3(8., 7., 17.)), vec3(2., 1., 2.));
+		assert_eq!(
+			cs.global_to_flat(Ray {
+				pos: vec3(7., 3., 0.),
+				dir: vec3(3., 2., 0.)
+			}),
+			Ray {
+				pos: vec3(23., 15., 7.),
+				dir: vec3(9., 8., 0.)
+			}
+		);
+		assert_eq!(
+			cs.flat_to_global(Ray {
+				pos: vec3(23., 15., 7.),
+				dir: vec3(9., 8., 0.)
+			}),
+			Ray {
+				pos: vec3(7., 3., 0.),
+				dir: vec3(3., 2., 0.)
+			}
+		);
+		assert_eq!(
+			cs.global_to_flat(Location {
+				pos: vec3(2., 1., 0.),
+				rot: mat3(vec3(0., 1., 0.), vec3(-1., 0., 0.), vec3(0., 0., 1.))
+			}),
+			Location {
+				pos: vec3(8., 7., 7.),
+				rot: mat3(vec3(0., 4., 0.), vec3(-3., 0., 0.), vec3(0., 0., 5.))
+			}
+		);
+		assert_eq!(
+			cs.flat_to_global(Location {
+				pos: vec3(2., 1., 7.),
+				rot: mat3(vec3(0., 1., 0.), vec3(-1., 0., 0.), vec3(0., 0., 1.))
+			}),
+			Location {
+				pos: vec3(0., -0.5, 0.),
+				rot: mat3(
+					vec3(0., 0.25, 0.),
+					vec3(-1. / 3., 0., 0.),
+					vec3(0., 0., 0.2)
+				)
+			}
+		);
+	}
+
+	fn test_flat_region(
+		region: &impl FlatRegion,
+		range_global: (Vec3, Vec3),
+		range_flat: (Vec3, Vec3),
+	) {
+		#[allow(non_upper_case_globals)]
+		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: Vec3,
+			range_a: (Vec3, Vec3),
+			name_b: &str,
+			b: Vec3,
+			range_b: (Vec3, Vec3),
+		) {
+			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) {
+				for z in linspace(range_global.0.z, range_global.1.z, 20) {
+					let pos_global = vec3(x, y, z);
+					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: Mat3::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: Mat3::IDENTITY
+							}))
+							.rot,
+						Mat3::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) {
+				for z in linspace(range_flat.0.z, range_flat.1.z, 20) {
+					let pos_flat = vec3(x, y, z);
+					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: Mat3::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: Mat3::IDENTITY
+							}))
+							.rot,
+						Mat3::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,
+			(vec3(-30.0, -300.0, -30.0), vec3(30.0, 300.0, 30.0)),
+			(vec3(-30.0, -100.0, -30.0), vec3(30.0, 100.0, 30.0)),
+		);
+		test_flat_region(
+			&mapper,
+			(vec3(-60.0, -400.0, -60.0), vec3(60.0, -300.0, 60.0)),
+			(vec3(-60.0, -200.0, -60.0), vec3(60.0, -100.0, 60.0)),
+		);
+		test_flat_region(
+			&mapper,
+			(vec3(-60.0, 300.0, -60.0), vec3(60.0, 400.0, 60.0)),
+			(vec3(-60.0, 100.0, -60.0), vec3(60.0, 200.0, 60.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,
+			(vec3(-30.0, -300.0, -30.0), vec3(30.0, -1.0, 30.0)),
+			(vec3(-30.0, -300.0, -30.0), vec3(30.0, -200.20016, 30.0)),
+		);
+		test_flat_region(
+			&mapper,
+			(vec3(-30.0, 1.0, -30.0), vec3(30.0, 300.0, 30.0)),
+			(vec3(-30.0, 200.20016, -30.0), vec3(30.0, 300.0, 30.0)),
+		);
+		test_flat_region(
+			&mapper,
+			(vec3(-60.0, -400.0, -60.0), vec3(60.0, -300.0, 60.0)),
+			(vec3(-60.0, -400.0, -60.0), vec3(60.0, -300.0, 60.0)),
+		);
+		test_flat_region(
+			&mapper,
+			(vec3(-60.0, 300.0, -60.0), vec3(60.0, 400.0, 60.0)),
+			(vec3(-60.0, 300.0, -60.0), vec3(60.0, 400.0, 60.0)),
+		);
+		// straight
+		for x in linspace(-60., 60., 20) {
+			for y in linspace(-320., 320., 20) {
+				for z in linspace(-60., 60., 20) {
+					assert_eq!(
+						mapper
+							.global_to_flat(Location {
+								pos: vec3(x, y, z),
+								rot: Mat3::IDENTITY
+							})
+							.pos
+							.x,
+						x
+					);
+				}
+			}
+		}
+		// symmetrical
+		for x in linspace(0., 60., 20) {
+			for y in linspace(0., 320., 20) {
+				for z in linspace(0., 60., 20) {
+					let pp = mapper
+						.global_to_flat(Location {
+							pos: vec3(x, y, z),
+							rot: Mat3::IDENTITY,
+						})
+						.pos;
+					let np = mapper
+						.global_to_flat(Location {
+							pos: vec3(-x, y, z),
+							rot: Mat3::IDENTITY,
+						})
+						.pos;
+					let pn = mapper
+						.global_to_flat(Location {
+							pos: vec3(x, -y, z),
+							rot: Mat3::IDENTITY,
+						})
+						.pos;
+					let nn = mapper
+						.global_to_flat(Location {
+							pos: vec3(-x, -y, z),
+							rot: Mat3::IDENTITY,
+						})
+						.pos;
+					assert_eq!(np, vec3(-pp.x, pp.y, pp.z));
+					assert_eq!(pn, vec3(pp.x, -pp.y, pp.z));
+					assert_eq!(nn, vec3(-pp.x, -pp.y, pp.z));
+				}
+			}
+		}
+		// clean boundary
+		for x in linspace(50., 60., 20) {
+			for y in linspace(0., 320., 20) {
+				for z in linspace(50., 60., 20) {
+					assert_eq!(
+						mapper
+							.global_to_flat(Location {
+								pos: vec3(x, y, z),
+								rot: Mat3::IDENTITY
+							})
+							.pos
+							.y,
+						y
+					);
+				}
+			}
+		}
+		for x in linspace(0., 60., 20) {
+			for y in linspace(300., 320., 20) {
+				for z in linspace(0., 60., 20) {
+					assert_eq!(
+						mapper
+							.global_to_flat(Location {
+								pos: vec3(x, y, z),
+								rot: Mat3::IDENTITY
+							})
+							.pos
+							.y,
+						y
+					);
+				}
+			}
+		}
+		// accelerating
+		for x in linspace(-29., 29., 20) {
+			for y in linspace(1., 299., 20) {
+				for z in linspace(-29., 29., 20) {
+					let v = mapper
+						.global_to_flat(Location {
+							pos: vec3(x, y, z),
+							rot: Mat3::IDENTITY,
+						})
+						.pos
+						.y;
+					assert!(v > 200.0);
+					assert!(v > y);
+				}
+			}
+		}
+	}
+}

src/tube/metric.rs

@@ -1,7 +1,8 @@
-use glam::{f32, vec2, Mat2, Vec2};
+use glam::{f32, vec3, Mat3, Vec3};
 
 use crate::fns::{self, Limiter};
-use crate::riemann::{Decomp2, Metric, Tens2};
+use crate::mathx::Decomp3;
+use crate::riemann::{Metric, Tens3};
 
 #[derive(Copy, Clone, Debug)]
 pub struct Tube {
@@ -40,20 +41,20 @@ impl Tube {
 }
 
 impl Metric for Tube {
-	fn sqrt_at(&self, pos: Vec2) -> Decomp2 {
+	fn sqrt_at(&self, pos: Vec3) -> Decomp3 {
 		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),
+		Decomp3 {
+			ortho: Mat3::IDENTITY,
+			diag: vec3(1.0, s, 1.0),
 		}
 	}
 
-	fn part_derivs_at(&self, pos: Vec2) -> Tens2 {
+	fn part_derivs_at(&self, pos: Vec3) -> Tens3 {
 		let sx = self.fx().value(pos.x);
 		let sy = self.fy().du(pos.y);
 		let s = sx + sy - sx * sy;
@@ -62,8 +63,9 @@ impl Metric for Tube {
 		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]),
+			Mat3::from_cols_array(&[0., 0., 0., 0., ds2_dx, 0., 0., 0., 0.]),
+			Mat3::from_cols_array(&[0., 0., 0., 0., ds2_dy, 0., 0., 0., 0.]),
+			Mat3::from_cols_array(&[0., 0., 0., 0., 0., 0., 0., 0., 0.]),
 		]
 	}
 }
@@ -71,10 +73,11 @@ impl Metric for Tube {
 #[cfg(test)]
 mod test {
 	use approx::assert_abs_diff_eq;
-	use glam::{vec2, Vec2};
+	use glam::{vec3, Vec3};
 	use itertools_num::linspace;
 
-	use crate::riemann::{Decomp2, Metric};
+	use crate::mathx::Decomp3;
+	use crate::riemann::Metric;
 	use crate::tube::Space;
 	use crate::types::Ray;
 
@@ -84,7 +87,7 @@ mod test {
 	fn test_tube_metric_derivs() {
 		struct Approx(Tube);
 		impl Metric for Approx {
-			fn sqrt_at(&self, pos: Vec2) -> Decomp2 {
+			fn sqrt_at(&self, pos: Vec3) -> Decomp3 {
 				self.0.sqrt_at(pos)
 			}
 		}
@@ -98,18 +101,25 @@ mod test {
 		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 x in itertools_num::linspace(-mul * testee.outer_radius, mul * testee.outer_radius, 20)
 		{
 			for y in itertools_num::linspace(
 				-mul * testee.external_halflength,
 				mul * testee.external_halflength,
-				100,
+				20,
 			) {
-				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");
+				for z in itertools_num::linspace(
+					-mul * testee.outer_radius,
+					mul * testee.outer_radius,
+					20,
+				) {
+					let pos = vec3(x, y, z);
+					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");
+				}
 			}
 		}
 	}
@@ -130,9 +140,9 @@ mod test {
 		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 a = vec3(ax, -(space.tube.external_halflength + off), 0.);
+				let b = vec3(bx, space.tube.external_halflength + off, 0.);
+				let Δ = vec3(bx - ax, 2.0 * (space.tube.internal_halflength + off), 0.);
 				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);
@@ -168,9 +178,9 @@ mod test {
 				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 a = vec3(ax, -(space.tube.external_halflength + off), 0.);
+				let b = vec3(bx, space.tube.external_halflength + off, 0.);
+				let Δ = vec3(bx - ax, 2.0 * (space.tube.internal_halflength + off), 0.);
 				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);
@@ -196,9 +206,9 @@ mod test {
 		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 a = vec3(x, -(space.tube.external_halflength + off), 0.);
+			let b = vec3(x, space.tube.external_halflength + off, 0.);
+			let Δ = vec3(0.0, 2.0 * (space.tube.internal_halflength + off), 0.);
 			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);
@@ -223,9 +233,9 @@ mod test {
 		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 a = vec3(x, -(space.tube.external_halflength + off), 0.);
+			let b = vec3(x, space.tube.external_halflength + off, 0.);
+			let Δ = vec3(0.0, 2.0 * (space.tube.external_halflength + off), 0.);
 			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);

src/tube/mod.rs

@@ -1,4 +1,4 @@
-use glam::{bool, f32, vec2, Mat2, Vec2};
+use glam::{bool, f32, vec3, Mat3, Vec3};
 
 use crate::ifaces::{DebugTraceable, RayPath, Traceable};
 use coords::{FlatCoordinateSystem, InnerCS, OuterCS};
@@ -26,7 +26,7 @@ pub enum Subspace {
 }
 
 impl Space {
-	fn which_subspace(&self, pt: Vec2) -> Subspace {
+	fn which_subspace(&self, pt: Vec3) -> Subspace {
 		if pt.y.abs() > self.tube.external_halflength {
 			Outer
 		} else if pt.x.abs() > self.tube.outer_radius {
@@ -41,7 +41,7 @@ impl Space {
 	/// Выполняет один шаг трассировки. Работает в любой части пространства, но вне Boundary доступны более эффективные методы.
 	/// ray задаётся в основной СК.
 	pub fn trace_step(&self, ray: Ray) -> Ray {
-		let a: Vec2 = -riemann::contract2(riemann::krist(&self.tube, ray.pos), ray.dir);
+		let a = -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 }
@@ -49,9 +49,9 @@ impl Space {
 
 	/// Выполняет один шаг перемещения. Работает в любой части пространства.
 	/// off задаётся в локальной СК. Рекомендуется считать небольшими шагами.
-	pub fn move_step(&self, loc: Location, off: Vec2) -> Location {
+	pub fn move_step(&self, loc: Location, off: Vec3) -> Location {
 		let corr =
-			Mat2::IDENTITY - riemann::contract(riemann::krist(&self.tube, loc.pos), loc.rot * off);
+			Mat3::IDENTITY - riemann::contract(riemann::krist(&self.tube, loc.pos), loc.rot * off);
 		let p = loc.pos + corr * loc.rot * off;
 		Location {
 			pos: p,
@@ -73,7 +73,7 @@ impl Space {
 		self.trace_flat(OuterCS(self.tube), ray)
 	}
 
-	fn obj_hitter(&self, pos: Vec2) -> Option<fn(&Self, ray: Ray) -> FlatTraceResult> {
+	fn obj_hitter(&self, pos: Vec3) -> Option<fn(&Self, ray: Ray) -> FlatTraceResult> {
 		match self.which_subspace(pos) {
 			Inner => Some(Self::trace_inner),
 			Outer => Some(Self::trace_outer),
@@ -113,7 +113,7 @@ impl Space {
 		objs: &[Object],
 		ray: Ray,
 		limit: Option<f32>,
-		globalize: impl Fn(Vec2) -> Vec2,
+		globalize: impl Fn(Vec3) -> Vec3,
 	) -> Vec<Hit> {
 		let limit = limit.unwrap_or(f32::INFINITY);
 		objs.iter()
@@ -146,7 +146,7 @@ impl Space {
 			.collect()
 	}
 
-	pub fn line(&self, a: Vec2, b: Vec2, step: f32) -> Vec<Vec2> {
+	pub fn line(&self, a: Vec3, b: Vec3, step: f32) -> Vec<Vec3> {
 		match self.which_subspace(a) {
 			Outer => vec![b],
 			Inner => {
@@ -210,7 +210,7 @@ impl DebugTraceable for Space {
 		let mut hits = vec![];
 		let mut ray = self.camera_ray_to_abs(camera, ray);
 
-		let trace_to_flat = |points: &mut Vec<Vec2>, ray| {
+		let trace_to_flat = |points: &mut Vec<Vec3>, ray| {
 			for ray in self.trace_iter(ray).skip(1) {
 				points.push(ray.pos);
 				if let Some(hitter) = self.obj_hitter(ray.pos) {
@@ -241,7 +241,7 @@ impl DebugTraceable for Space {
 }
 
 struct Rect {
-	pub size: Vec2,
+	pub size: Vec3,
 }
 
 impl Rect {
@@ -253,7 +253,7 @@ impl Rect {
 		}
 	}
 
-	fn is_inside(&self, pt: Vec2) -> bool {
+	fn is_inside(&self, pt: Vec3) -> bool {
 		pt.abs().cmplt(self.size).all()
 	}
 
@@ -285,146 +285,146 @@ impl Rect {
 fn test_rect() {
 	assert_eq!(
 		Rect::flip_ray(Ray {
-			pos: vec2(2.0, 3.0),
-			dir: vec2(4.0, 5.0)
+			pos: vec3(2.0, 3.0, 2.0),
+			dir: vec3(4.0, 5.0, 4.0)
 		}),
 		Ray {
-			pos: vec2(2.0, 3.0),
-			dir: vec2(4.0, 5.0)
+			pos: vec3(2.0, 3.0, 2.0),
+			dir: vec3(4.0, 5.0, 4.0)
 		}
 	);
 	assert_eq!(
 		Rect::flip_ray(Ray {
-			pos: vec2(2.0, 3.0),
-			dir: vec2(-4.0, 5.0)
+			pos: vec3(2.0, 3.0, 2.0),
+			dir: vec3(-4.0, 5.0, -4.0)
 		}),
 		Ray {
-			pos: vec2(-2.0, 3.0),
-			dir: vec2(4.0, 5.0)
+			pos: vec3(-2.0, 3.0, -2.0),
+			dir: vec3(4.0, 5.0, 4.0)
 		}
 	);
 	assert_eq!(
 		Rect::flip_ray(Ray {
-			pos: vec2(2.0, 3.0),
-			dir: vec2(4.0, -5.0)
+			pos: vec3(2.0, 3.0, 2.0),
+			dir: vec3(4.0, -5.0, 4.0)
 		}),
 		Ray {
-			pos: vec2(2.0, -3.0),
-			dir: vec2(4.0, 5.0)
+			pos: vec3(2.0, -3.0, 2.0),
+			dir: vec3(4.0, 5.0, 4.0)
 		}
 	);
 	assert_eq!(
 		Rect::flip_ray(Ray {
-			pos: vec2(2.0, 3.0),
-			dir: vec2(4.0, 0.0)
+			pos: vec3(2.0, 3.0, 2.0),
+			dir: vec3(4.0, 0.0, 4.0)
 		}),
 		Ray {
-			pos: vec2(2.0, 3.0),
-			dir: vec2(4.0, 0.0)
+			pos: vec3(2.0, 3.0, 2.0),
+			dir: vec3(4.0, 0.0, 4.0)
 		}
 	);
 
 	let r = Rect {
-		size: vec2(2.0, 3.0),
+		size: vec3(2.0, 3.0, 2.0),
 	};
 
 	assert_eq!(
 		r.trace_into(Ray {
-			pos: vec2(3.0, 3.0),
-			dir: vec2(1.0, 1.0)
+			pos: vec3(3.0, 3.0, 3.0),
+			dir: vec3(1.0, 1.0, 1.0)
 		}),
 		None
 	);
 	assert_eq!(
 		r.trace_into(Ray {
-			pos: vec2(-3.0, 2.0),
-			dir: vec2(1.0, 0.0)
+			pos: vec3(-3.0, 2.0, -3.0),
+			dir: vec3(1.0, 0.0, 1.0)
 		}),
 		Some(1.0)
 	);
 	assert_eq!(
 		r.trace_into(Ray {
-			pos: vec2(-3.0, 2.0),
-			dir: vec2(-1.0, 0.0)
+			pos: vec3(-3.0, 2.0, -3.0),
+			dir: vec3(-1.0, 0.0, -1.0)
 		}),
 		None
 	);
 	assert_eq!(
 		r.trace_into(Ray {
-			pos: vec2(-3.0, 1.0),
-			dir: vec2(2.0, 2.0)
+			pos: vec3(-3.0, 1.0, -3.0),
+			dir: vec3(2.0, 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)
+			pos: vec3(-3.0, 2.1, -3.0),
+			dir: vec3(2.0, 2.0, 2.0)
 		}),
 		None
 	);
 
 	assert_eq!(
 		r.trace_into(Ray {
-			pos: vec2(2.0, 3.0),
-			dir: vec2(1.0, 1.0)
+			pos: vec3(2.0, 3.0, 2.0),
+			dir: vec3(1.0, 1.0, 1.0)
 		}),
 		None
 	);
 	assert_eq!(
 		r.trace_into(Ray {
-			pos: vec2(-2.0, 3.0),
-			dir: vec2(-1.0, 1.0)
+			pos: vec3(-2.0, 3.0, -2.0),
+			dir: vec3(-1.0, 1.0, -1.0)
 		}),
 		None
 	);
 	assert_eq!(
 		r.trace_into(Ray {
-			pos: vec2(2.0, 3.0),
-			dir: vec2(-1.0, -1.0)
+			pos: vec3(2.0, 3.0, 2.0),
+			dir: vec3(-1.0, -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)
+			pos: vec3(2.0, -3.0, 2.0),
+			dir: vec3(-1.0, 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)
+			pos: vec3(0.0, 0.0, 0.0),
+			dir: vec3(1.0, 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)
+			pos: vec3(0.0, 0.0, 0.0),
+			dir: vec3(0.0, 1.0, 0.0)
 		}),
 		Some(3.0)
 	);
 	assert_eq!(
 		r.trace_out_of(Ray {
-			pos: vec2(0.0, 1.0),
-			dir: vec2(0.0, -1.0)
+			pos: vec3(0.0, 1.0, 0.0),
+			dir: vec3(0.0, -1.0, 0.0)
 		}),
 		Some(4.0)
 	);
 	assert_eq!(
 		r.trace_out_of(Ray {
-			pos: vec2(1.0, 1.0),
-			dir: vec2(0.0, -1.0)
+			pos: vec3(1.0, 1.0, 1.0),
+			dir: vec3(0.0, -1.0, 0.0)
 		}),
 		Some(4.0)
 	);
 	assert_eq!(
 		r.trace_out_of(Ray {
-			pos: vec2(2.0, 3.0),
-			dir: vec2(1.0, 1.0)
+			pos: vec3(2.0, 3.0, 2.0),
+			dir: vec3(1.0, 1.0, 1.0)
 		}),
 		Some(0.0)
 	);

src/types.rs

@@ -1,9 +1,9 @@
-use glam::{f32, i32, Mat2, Vec2};
+use glam::{f32, i32, Mat3, Vec3};
 
 #[derive(Copy, Clone, Debug, PartialEq)]
 pub struct Ray {
-	pub pos: Vec2,
-	pub dir: Vec2,
+	pub pos: Vec3,
+	pub dir: Vec3,
 }
 
 impl Ray {
@@ -15,7 +15,7 @@ impl Ray {
 	}
 }
 
-impl std::ops::Mul<Ray> for Mat2 {
+impl std::ops::Mul<Ray> for Mat3 {
 	type Output = Ray;
 
 	fn mul(self, rhs: Ray) -> Self::Output {
@@ -29,9 +29,9 @@ impl std::ops::Mul<Ray> for Mat2 {
 #[derive(Copy, Clone, Debug, PartialEq)]
 pub struct Location {
 	/// Положение в основной СК
-	pub pos: Vec2,
+	pub pos: Vec3,
 	/// Преобразование вектора из локальной ортонормированной в основную СК
-	pub rot: Mat2,
+	pub rot: Mat3,
 }
 
 #[derive(Copy, Clone, Debug)]
@@ -45,7 +45,7 @@ pub struct Object {
 pub struct Hit {
 	pub distance: f32,
 	pub id: i32,
-	pub pos: Vec2, // положение в основной СК
+	pub pos: Vec3, // положение в основной СК
 	pub rel: Ray,  // в локальной ортонормированной СК объекта
 }