Support 3D!

src/bin/flat/main.rs

@@ -12,7 +12,7 @@ use refraction::tube::Space;
 use refraction::types::{Location, Object, Ray};
 use refraction::DT;
 
-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;
@@ -43,23 +43,31 @@ pub fn main() {
 					id: k as i32,
 					loc: put_object(
 						&tube,
-						vec2(0.0, y * tube.external_halflength),
+						vec3(0.0, y * tube.external_halflength, 0.0),
-						Mat2::from_angle(y),
+						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),
+				vec3(
-				mat2(vec2(1., -1.), vec2(1., 1.)),
+					-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),
+				vec3(
-				mat2(vec2(0., -1.), vec2(1., 0.)),
+					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);
@@ -99,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
 						}),
@@ -110,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
 						}),
@@ -123,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 })
@@ -136,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)
@@ -144,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();
@@ -157,26 +168,58 @@ fn test_put_object() {
 
 	let ε = 1e-5;
 	let m = refraction::riemann::samples::ScaledMetric {
-		scale: vec2(3., 4.),
+		scale: vec3(3., 4., 5.),
 	};
 
-	let loc = put_object(&m, vec2(1., 2.), mat2(vec2(1., 0.), vec2(0., 1.)));
+	let loc = put_object(
-	assert_eq!(loc.pos, vec2(1., 2.));
+		&m,
-	assert_abs_diff_eq!(loc.rot * vec2(1., 0.), vec2(1. / 3., 0.), epsilon = ε);
+		vec3(1., 2., 0.),
-	assert_abs_diff_eq!(loc.rot * vec2(0., 1.), vec2(0., 1. / 4.), epsilon = ε);
+		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.)));
+	let loc = put_object(
-	assert_eq!(loc.pos, vec2(1., 2.));
+		&m,
-	assert_abs_diff_eq!(loc.rot * vec2(1., 0.), vec2(0., 1. / 4.), epsilon = ε);
+		vec3(1., 2., 0.),
-	assert_abs_diff_eq!(loc.rot * vec2(0., 1.), vec2(-1. / 3., 0.), epsilon = ε);
+		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)));
+	let loc = put_object(
-	assert_eq!(loc.pos, vec2(1., 2.));
+		&m,
-	assert_abs_diff_eq!(loc.rot * vec2(1., 0.), vec2(1. / 5., 1. / 5.), epsilon = ε);
+		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.),
+		loc.rot * vec3(1., 0., 0.),
-		vec2(-4. / 15., 3. / 20.),
+		vec3(1. / 5., 1. / 5., 0.),
+		epsilon = ε
+	);
+	assert_abs_diff_eq!(
+		loc.rot * vec3(0., 1., 0.),
+		vec3(-4. / 15., 3. / 20., 0.),
 		epsilon = ε
 	);
 }
@@ -188,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) {
+fn draw_ray_2(gc: &mut Vec<Draw>, space: &Space, camera: Location, dir: Vec3) {
-	let pos = vec2(0., 0.);
+	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()) {
@@ -214,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.));
 	}
 }
 
@@ -225,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,
+		pos: vec3(start.x, start.y, 0.),
-		rot: mat2(dir, vec2(-dir.y, dir.x)),
+		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/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 points: Vec<Vec3>,
-	pub end_dir: Vec2,
+	pub end_dir: Vec3,
 }
 
 pub trait DebugTraceable: Traceable {

src/riemann.rs

@@ -1,50 +1,50 @@
-use crate::mathx::Decomp2;
+use crate::mathx::Decomp3;
 use glam::*;
 
-pub type Tens2 = [Mat2; 2];
+pub type Tens3 = [Mat3; 3];
 
 pub trait Metric {
-	fn sqrt_at(&self, pos: Vec2) -> Decomp2;
+	fn sqrt_at(&self, pos: Vec3) -> Decomp3;
 
-	fn at(&self, pos: Vec2) -> Mat2 {
+	fn at(&self, pos: Vec3) -> Mat3 {
 		self.sqrt_at(pos).square().into()
 	}
 
-	fn inverse_at(&self, pos: Vec2) -> Mat2 {
+	fn inverse_at(&self, pos: Vec3) -> Mat3 {
 		self.sqrt_at(pos).square().inverse().into()
 	}
 
-	fn part_derivs_at(&self, pos: Vec2) -> Tens2 {
+	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: Vec2, v: Vec2) -> f32 {
+	fn vec_length_at(&self, at: Vec3, v: Vec3) -> f32 {
 		v.dot(self.at(at) * v).sqrt()
 	}
 
-	fn normalize_vec_at(&self, at: Vec2, v: Vec2) -> Vec2 {
+	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: Vec2,
+	p: Vec3,
-	v: Vec2,
+	v: Vec3,
 }
 
 impl<'a, M: Metric> Iterator for TraceIter<'a, M> {
-	type Item = Vec2;
+	type Item = Vec3;
 
 	fn next(&mut self) -> Option<Self::Item> {
-		let a: Vec2 = -contract2(krist(self.space, self.p), self.v);
+		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: Vec2, dir: Vec2, dt: f32) -> TraceIter<M> {
+pub fn trace_iter<M: Metric>(space: &M, base: Vec3, dir: Vec3, dt: f32) -> TraceIter<M> {
 	TraceIter {
 		space,
 		p: base,
@@ -52,54 +52,55 @@ pub fn trace_iter<M: Metric>(space: &M, base: Vec2, dir: Vec2, dt: f32) -> Trace
 	}
 }
 
-pub fn krist(space: &impl Metric, pos: Vec2) -> Tens2 {
+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_tens2(|i, l, k| {
+	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(Vec2) -> Mat2, pos: Vec2, delta: Vec2) -> Mat2 {
+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(Vec2) -> Mat2, pos: Vec2, eps: f32) -> Tens2 {
+fn part_deriv(f: impl Fn(Vec3) -> Mat3, pos: Vec3, eps: f32) -> Tens3 {
 	[
-		dir_deriv(&f, pos, vec2(eps, 0.0)),
+		dir_deriv(&f, pos, vec3(eps, 0.0, 0.0)),
-		dir_deriv(&f, pos, vec2(0.0, eps)),
+		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: Tens2, u: Vec2) -> Mat2 {
+pub fn contract(t: Tens3, u: Vec3) -> Mat3 {
-	mat2(t[0] * u, t[1] * u).transpose()
+	mat3(t[0] * u, t[1] * u, t[2] * u).transpose()
 }
 
 /// Сворачивает тензор t с вектором v дважды, по второму и третьему индексам.
-pub fn contract2(t: Tens2, v: Vec2) -> Vec2 {
+pub fn contract2(t: Tens3, v: Vec3) -> Vec3 {
 	contract(t, v) * v
 }
 
-fn make_vec2(f: impl Fn(usize) -> f32) -> Vec2 {
+fn make_vec3(f: impl Fn(usize) -> f32) -> Vec3 {
-	Vec2::from_array(std::array::from_fn(|i| f(i)))
+	Vec3::from_array(std::array::from_fn(|i| f(i)))
 }
 
-fn make_mat2(f: impl Fn(usize, usize) -> f32) -> Mat2 {
+fn make_mat3(f: impl Fn(usize, usize) -> f32) -> Mat3 {
-	Mat2::from_cols_array_2d(&std::array::from_fn(|i| std::array::from_fn(|j| f(i, j))))
+	Mat3::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 {
+fn make_tens3(f: impl Fn(usize, usize, usize) -> f32) -> Tens3 {
-	std::array::from_fn(|i| make_mat2(|j, k| f(i, j, k)))
+	std::array::from_fn(|i| make_mat3(|j, k| f(i, j, k)))
 }
 
 #[test]
-fn m2() {
+fn m3() {
-	let m = make_mat2(|i, j| (i + 2 * j) as f32);
+	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);
@@ -107,8 +108,8 @@ fn m2() {
 }
 
 #[test]
-fn t2() {
+fn t3() {
-	let t = make_tens2(|i, j, k| (i + 2 * j + 4 * k) as f32);
+	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);
@@ -120,20 +121,20 @@ fn t2() {
 }
 
 pub mod samples {
-	use glam::{Mat2, Vec2};
+	use glam::{Mat3, Vec3};
 
-	use super::{Decomp2, Metric};
+	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: Vec2,
+		pub scale: Vec3,
 	}
 
 	impl Metric for ScaledMetric {
-		fn sqrt_at(&self, _pos: Vec2) -> Decomp2 {
+		fn sqrt_at(&self, _pos: Vec3) -> Decomp3 {
-			Decomp2 {
+			Decomp3 {
 				diag: self.scale,
-				ortho: Mat2::IDENTITY,
+				ortho: Mat3::IDENTITY,
 			}
 		}
 	}
@@ -144,99 +145,103 @@ mod tests {
 	use super::*;
 	use approx::assert_abs_diff_eq;
 
-	use glam::{mat2, vec2, Mat2};
+	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: vec2(3., 4.),
+			scale: vec3(3., 4., 5.),
 		};
 		assert_eq!(
 			metric.sqrt_at(rng.gen()),
-			Decomp2 {
+			Decomp3 {
-				ortho: Mat2::IDENTITY,
+				ortho: Mat3::IDENTITY,
-				diag: vec2(3., 4.)
+				diag: vec3(3., 4., 5.)
 			}
 		);
 		assert_eq!(
 			metric.at(rng.gen()),
-			Mat2::from_cols_array(&[9., 0., 0., 16.])
+			Mat3::from_cols_array(&[9., 0., 0., 0., 16., 0., 0., 0., 25.])
 		);
 		assert_eq!(
 			metric.inverse_at(rng.gen()),
-			Mat2::from_cols_array(&[1. / 9., 0., 0., 1. / 16.])
+			Mat3::from_cols_array(&[1. / 9., 0., 0., 0., 1. / 16., 0., 0., 0., 1. / 25.])
-		);
-		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.)),
+			metric.part_derivs_at(rng.gen()),
-			vec2(0., 1. / 4.)
+			[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(), vec2(1., 1.)),
+			metric.normalize_vec_at(rng.gen(), vec3(0., 1., 0.)),
-			vec2(1. / 5., 1. / 5.)
+			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: vec2(2., 4.),
+			scale: vec3(2., 4., 3.),
 		};
 		assert_eq!(
-			trace_iter(&metric, vec2(3., 5.), vec2(1., 0.), 1.)
+			trace_iter(&metric, vec3(3., 5., 0.), vec3(1., 0., 0.), 1.)
 				.nth(7)
 				.unwrap(),
-			vec2(7., 5.)
+			vec3(7., 5., 0.)
 		);
 		assert_eq!(
-			trace_iter(&metric, vec2(3., 5.), vec2(2., 0.), 1.)
+			trace_iter(&metric, vec3(3., 5., 0.), vec3(2., 0., 0.), 1.)
 				.nth(7)
 				.unwrap(),
-			vec2(7., 5.)
+			vec3(7., 5., 0.)
 		);
 		assert_eq!(
-			trace_iter(&metric, vec2(3., 5.), vec2(1., 0.), 0.5)
+			trace_iter(&metric, vec3(3., 5., 0.), vec3(1., 0., 0.), 0.5)
 				.nth(7)
 				.unwrap(),
-			vec2(5., 5.)
+			vec3(5., 5., 0.)
 		);
 		assert_eq!(
-			trace_iter(&metric, vec2(3., 5.), vec2(0., 1.), 1.)
+			trace_iter(&metric, vec3(3., 5., 0.), vec3(0., 1., 0.), 1.)
 				.nth(9)
 				.unwrap(),
-			vec2(3., 7.5)
+			vec3(3., 7.5, 0.)
 		);
 		assert_eq!(
-			trace_iter(&metric, vec2(3., 5.), vec2(0., 4.), 1.)
+			trace_iter(&metric, vec3(3., 5., 0.), vec3(0., 4., 0.), 1.)
 				.nth(9)
 				.unwrap(),
-			vec2(3., 7.5)
+			vec3(3., 7.5, 0.)
 		);
 		assert_eq!(
-			trace_iter(&metric, vec2(3., 5.), vec2(0., 1.), 0.5)
+			trace_iter(&metric, vec3(3., 5., 0.), vec3(0., 1., 0.), 0.5)
 				.nth(9)
 				.unwrap(),
-			vec2(3., 6.25)
+			vec3(3., 6.25, 0.)
 		);
 		assert_abs_diff_eq!(
 			trace_iter(
 				&metric,
-				vec2(3., 5.),
+				vec3(3., 5., 0.),
-				vec2(0.5, 0.25),
+				vec3(0.5, 0.25, 0.),
 				std::f32::consts::SQRT_2
 			)
 			.nth(7)
 			.unwrap(),
-			vec2(7., 7.),
+			vec3(7., 7., 0.),
 			epsilon = 1e-5
 		);
 	}

src/tube/coords.rs

@@ -1,4 +1,4 @@
-use glam::{vec2, Mat2, Vec2};
+use glam::{vec3, Mat3, Vec3};
 
 use crate::riemann::Metric;
 use crate::types::{Location, Ray};
@@ -11,7 +11,7 @@ pub trait FlatCoordinateSystem<T> {
 }
 
 pub trait FlatRegion:
-	FlatCoordinateSystem<Vec2> + FlatCoordinateSystem<Ray> + FlatCoordinateSystem<Location>
+	FlatCoordinateSystem<Vec3> + FlatCoordinateSystem<Ray> + FlatCoordinateSystem<Location>
 {
 	// Измеряет расстояние до выхода за пределы области вдоль луча ray. Луч задаётся в плоской СК.
 	fn distance_to_boundary(&self, _ray: Ray) -> Option<f32> {
@@ -19,20 +19,20 @@ pub trait FlatRegion:
 	}
 }
 
-trait MetricCS: FlatCoordinateSystem<Vec2> {
+trait MetricCS: FlatCoordinateSystem<Vec3> {
 	fn global_metric(&self) -> &impl Metric;
-	fn flat_to_global_tfm_at(&self, pos: Vec2) -> Mat2 {
+	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: Vec2) -> Mat2 {
+	fn global_to_flat_tfm_at(&self, pos: Vec3) -> Mat3 {
 		self.global_metric().sqrt_at(pos).into()
 	}
 }
 
-impl<T: FlatCoordinateSystem<Vec2> + MetricCS> FlatCoordinateSystem<Ray> for T {
+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),
@@ -48,7 +48,7 @@ impl<T: FlatCoordinateSystem<Vec2> + MetricCS> FlatCoordinateSystem<Ray> for T {
 	}
 }
 
-impl<T: FlatCoordinateSystem<Vec2> + MetricCS> FlatCoordinateSystem<Location> for T {
+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),
@@ -72,21 +72,25 @@ impl MetricCS for InnerCS {
 	}
 }
 
-impl FlatCoordinateSystem<Vec2> for InnerCS {
+impl FlatCoordinateSystem<Vec3> for InnerCS {
-	fn flat_to_global(&self, pos: Vec2) -> Vec2 {
+	fn flat_to_global(&self, pos: Vec3) -> Vec3 {
-		vec2(pos.x, self.0.y(pos.y))
+		vec3(pos.x, self.0.y(pos.y), pos.z)
 	}
 
 	// Работает только при |pos.x| ≤ inner_radius или |pos.y| ≥ external_halflength.
-	fn global_to_flat(&self, pos: Vec2) -> Vec2 {
+	fn global_to_flat(&self, pos: Vec3) -> Vec3 {
-		vec2(pos.x, self.0.v(pos.y))
+		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: vec2(self.0.inner_radius, self.0.internal_halflength),
+			size: vec3(
+				self.0.inner_radius,
+				self.0.internal_halflength,
+				self.0.inner_radius,
+			),
 		}
 		.trace_out_of(ray)
 	}
@@ -100,31 +104,39 @@ impl MetricCS for OuterCS {
 	}
 }
 
-impl FlatCoordinateSystem<Vec2> for OuterCS {
+impl FlatCoordinateSystem<Vec3> for OuterCS {
-	fn flat_to_global(&self, pos: Vec2) -> Vec2 {
+	fn flat_to_global(&self, pos: Vec3) -> Vec3 {
 		let inner = Rect {
-			size: vec2(self.0.inner_radius + 1.0, self.0.external_halflength),
+			size: vec3(
+				self.0.inner_radius + 1.0,
+				self.0.external_halflength,
+				self.0.inner_radius + 1.0,
+			),
 		};
 		if inner.is_inside(pos) {
-			let Vec2 { x, y: v } = pos;
+			let Vec3 { x, y: v, z } = pos;
 			let y = self
 				.0
 				.y(v - v.signum() * (self.0.external_halflength - self.0.internal_halflength));
-			vec2(x, y)
+			vec3(x, y, z)
 		} else {
 			pos
 		}
 	}
 
-	fn global_to_flat(&self, pos: Vec2) -> Vec2 {
+	fn global_to_flat(&self, pos: Vec3) -> Vec3 {
 		let inner = Rect {
-			size: vec2(self.0.inner_radius + 1.0, self.0.external_halflength),
+			size: vec3(
+				self.0.inner_radius + 1.0,
+				self.0.external_halflength,
+				self.0.inner_radius + 1.0,
+			),
 		};
 		if inner.is_inside(pos) {
-			let Vec2 { x: u, y } = 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);
-			vec2(u, v) // в плоском продолжении СК Outer на область Inner
+			vec3(u, v, w) // в плоском продолжении СК Outer на область Inner
 		} else {
 			pos
 		}
@@ -134,7 +146,11 @@ impl FlatCoordinateSystem<Vec2> for OuterCS {
 impl FlatRegion for OuterCS {
 	fn distance_to_boundary(&self, ray: Ray) -> Option<f32> {
 		Rect {
-			size: vec2(self.0.outer_radius, self.0.external_halflength),
+			size: vec3(
+				self.0.outer_radius,
+				self.0.external_halflength,
+				self.0.outer_radius,
+			),
 		}
 		.trace_into(ray)
 	}
@@ -143,7 +159,7 @@ impl FlatRegion for OuterCS {
 #[cfg(test)]
 mod tests {
 	use approx::{assert_abs_diff_eq, AbsDiffEq};
-	use glam::{mat2, vec2, Mat2, Vec2};
+	use glam::{mat3, vec3, Mat3, Vec3};
 	use itertools_num::linspace;
 
 	use crate::riemann::samples;
@@ -152,12 +168,12 @@ mod tests {
 
 	#[test]
 	fn uniform_scaled_metric() {
-		struct Scaled(samples::ScaledMetric, Vec2);
+		struct Scaled(samples::ScaledMetric, Vec3);
-		impl FlatCoordinateSystem<Vec2> for Scaled {
+		impl FlatCoordinateSystem<Vec3> for Scaled {
-			fn flat_to_global(&self, pos: Vec2) -> Vec2 {
+			fn flat_to_global(&self, pos: Vec3) -> Vec3 {
 				(pos - self.1) / self.0.scale
 			}
-			fn global_to_flat(&self, pos: Vec2) -> Vec2 {
+			fn global_to_flat(&self, pos: Vec3) -> Vec3 {
 				pos * self.0.scale + self.1
 			}
 		}
@@ -168,58 +184,62 @@ mod tests {
 		}
 		let cs = Scaled(
 			samples::ScaledMetric {
-				scale: vec2(3., 4.),
+				scale: vec3(3., 4., 5.),
 			},
-			vec2(2., 3.),
+			vec3(2., 3., 7.),
 		);
-		assert_eq!(cs.global_to_flat(vec2(7., 3.)), vec2(23., 15.));
+		assert_eq!(cs.global_to_flat(vec3(7., 3., 1.)), vec3(23., 15., 12.));
-		assert_eq!(cs.flat_to_global(vec2(8., 7.)), vec2(2., 1.));
+		assert_eq!(cs.flat_to_global(vec3(8., 7., 17.)), vec3(2., 1., 2.));
 		assert_eq!(
 			cs.global_to_flat(Ray {
-				pos: vec2(7., 3.),
+				pos: vec3(7., 3., 0.),
-				dir: vec2(3., 2.)
+				dir: vec3(3., 2., 0.)
 			}),
 			Ray {
-				pos: vec2(23., 15.),
+				pos: vec3(23., 15., 7.),
-				dir: vec2(9., 8.)
+				dir: vec3(9., 8., 0.)
 			}
 		);
 		assert_eq!(
 			cs.flat_to_global(Ray {
-				pos: vec2(23., 15.),
+				pos: vec3(23., 15., 7.),
-				dir: vec2(9., 8.)
+				dir: vec3(9., 8., 0.)
 			}),
 			Ray {
-				pos: vec2(7., 3.),
+				pos: vec3(7., 3., 0.),
-				dir: vec2(3., 2.)
+				dir: vec3(3., 2., 0.)
 			}
 		);
 		assert_eq!(
 			cs.global_to_flat(Location {
-				pos: vec2(2., 1.),
+				pos: vec3(2., 1., 0.),
-				rot: mat2(vec2(0., 1.), vec2(-1., 0.))
+				rot: mat3(vec3(0., 1., 0.), vec3(-1., 0., 0.), vec3(0., 0., 1.))
 			}),
 			Location {
-				pos: vec2(8., 7.),
+				pos: vec3(8., 7., 7.),
-				rot: mat2(vec2(0., 4.), vec2(-3., 0.))
+				rot: mat3(vec3(0., 4., 0.), vec3(-3., 0., 0.), vec3(0., 0., 5.))
 			}
 		);
 		assert_eq!(
 			cs.flat_to_global(Location {
-				pos: vec2(2., 1.),
+				pos: vec3(2., 1., 7.),
-				rot: mat2(vec2(0., 1.), vec2(-1., 0.))
+				rot: mat3(vec3(0., 1., 0.), vec3(-1., 0., 0.), vec3(0., 0., 1.))
 			}),
 			Location {
-				pos: vec2(0., -0.5),
+				pos: vec3(0., -0.5, 0.),
-				rot: mat2(vec2(0., 0.25), vec2(-1. / 3., 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: (Vec2, Vec2),
+		range_global: (Vec3, Vec3),
-		range_flat: (Vec2, Vec2),
+		range_flat: (Vec3, Vec3),
 	) {
 		#[allow(non_upper_case_globals)]
 		const ε: f32 = 1e-3;
@@ -238,11 +258,11 @@ mod tests {
 		}
 		fn check_range(
 			name_a: &str,
-			a: Vec2,
+			a: Vec3,
-			range_a: (Vec2, Vec2),
+			range_a: (Vec3, Vec3),
 			name_b: &str,
-			b: Vec2,
+			b: Vec3,
-			range_b: (Vec2, Vec2),
+			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:
@@ -261,7 +281,8 @@ mod tests {
 		}
 		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);
+				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",
@@ -276,7 +297,7 @@ mod tests {
 						region
 							.global_to_flat(Location {
 								pos: pos_global,
-							rot: Mat2::IDENTITY
+								rot: Mat3::IDENTITY
 							})
 							.pos,
 						pos_flat
@@ -287,16 +308,18 @@ mod tests {
 						region
 							.flat_to_global(region.global_to_flat(Location {
 								pos: pos_global,
-							rot: Mat2::IDENTITY
+								rot: Mat3::IDENTITY
 							}))
 							.rot,
-					Mat2::IDENTITY
+						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) {
-				let pos_flat = vec2(x, y);
+				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",
@@ -311,7 +334,7 @@ mod tests {
 						region
 							.flat_to_global(Location {
 								pos: pos_flat,
-							rot: Mat2::IDENTITY
+								rot: Mat3::IDENTITY
 							})
 							.pos,
 						pos_global
@@ -322,14 +345,15 @@ mod tests {
 						region
 							.global_to_flat(region.flat_to_global(Location {
 								pos: pos_global,
-							rot: Mat2::IDENTITY
+								rot: Mat3::IDENTITY
 							}))
 							.rot,
-					Mat2::IDENTITY
+						Mat3::IDENTITY
 					);
 				}
 			}
 		}
+	}
 
 	#[test]
 	fn test_mapper_inner() {
@@ -341,18 +365,18 @@ mod tests {
 		});
 		test_flat_region(
 			&mapper,
-			(vec2(-30.0, -300.0), vec2(30.0, 300.0)),
+			(vec3(-30.0, -300.0, -30.0), vec3(30.0, 300.0, 30.0)),
-			(vec2(-30.0, -100.0), vec2(30.0, 100.0)),
+			(vec3(-30.0, -100.0, -30.0), vec3(30.0, 100.0, 30.0)),
 		);
 		test_flat_region(
 			&mapper,
-			(vec2(-60.0, -400.0), vec2(60.0, -300.0)),
+			(vec3(-60.0, -400.0, -60.0), vec3(60.0, -300.0, 60.0)),
-			(vec2(-60.0, -200.0), vec2(60.0, -100.0)),
+			(vec3(-60.0, -200.0, -60.0), vec3(60.0, -100.0, 60.0)),
 		);
 		test_flat_region(
 			&mapper,
-			(vec2(-60.0, 300.0), vec2(60.0, 400.0)),
+			(vec3(-60.0, 300.0, -60.0), vec3(60.0, 400.0, 60.0)),
-			(vec2(-60.0, 100.0), vec2(60.0, 200.0)),
+			(vec3(-60.0, 100.0, -60.0), vec3(60.0, 200.0, 60.0)),
 		);
 	}
 
@@ -367,32 +391,33 @@ mod tests {
 		// TODO replace 200.20016 with something sane
 		test_flat_region(
 			&mapper,
-			(vec2(-30.0, -300.0), vec2(30.0, -1.0)),
+			(vec3(-30.0, -300.0, -30.0), vec3(30.0, -1.0, 30.0)),
-			(vec2(-30.0, -300.0), vec2(30.0, -200.20016)),
+			(vec3(-30.0, -300.0, -30.0), vec3(30.0, -200.20016, 30.0)),
 		);
 		test_flat_region(
 			&mapper,
-			(vec2(-30.0, 1.0), vec2(30.0, 300.0)),
+			(vec3(-30.0, 1.0, -30.0), vec3(30.0, 300.0, 30.0)),
-			(vec2(-30.0, 200.20016), vec2(30.0, 300.0)),
+			(vec3(-30.0, 200.20016, -30.0), vec3(30.0, 300.0, 30.0)),
 		);
 		test_flat_region(
 			&mapper,
-			(vec2(-60.0, -400.0), vec2(60.0, -300.0)),
+			(vec3(-60.0, -400.0, -60.0), vec3(60.0, -300.0, 60.0)),
-			(vec2(-60.0, -400.0), vec2(60.0, -300.0)),
+			(vec3(-60.0, -400.0, -60.0), vec3(60.0, -300.0, 60.0)),
 		);
 		test_flat_region(
 			&mapper,
-			(vec2(-60.0, 300.0), vec2(60.0, 400.0)),
+			(vec3(-60.0, 300.0, -60.0), vec3(60.0, 400.0, 60.0)),
-			(vec2(-60.0, 300.0), vec2(60.0, 400.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: vec2(x, y),
+								pos: vec3(x, y, z),
-							rot: Mat2::IDENTITY
+								rot: Mat3::IDENTITY
 							})
 							.pos
 							.x,
@@ -400,46 +425,50 @@ mod tests {
 					);
 				}
 			}
+		}
 		// 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: vec2(x, y),
+							pos: vec3(x, y, z),
-						rot: Mat2::IDENTITY,
+							rot: Mat3::IDENTITY,
 						})
 						.pos;
 					let np = mapper
 						.global_to_flat(Location {
-						pos: vec2(-x, y),
+							pos: vec3(-x, y, z),
-						rot: Mat2::IDENTITY,
+							rot: Mat3::IDENTITY,
 						})
 						.pos;
 					let pn = mapper
 						.global_to_flat(Location {
-						pos: vec2(x, -y),
+							pos: vec3(x, -y, z),
-						rot: Mat2::IDENTITY,
+							rot: Mat3::IDENTITY,
 						})
 						.pos;
 					let nn = mapper
 						.global_to_flat(Location {
-						pos: vec2(-x, -y),
+							pos: vec3(-x, -y, z),
-						rot: Mat2::IDENTITY,
+							rot: Mat3::IDENTITY,
 						})
 						.pos;
-				assert_eq!(np, vec2(-pp.x, pp.y));
+					assert_eq!(np, vec3(-pp.x, pp.y, pp.z));
-				assert_eq!(pn, vec2(pp.x, -pp.y));
+					assert_eq!(pn, vec3(pp.x, -pp.y, pp.z));
-				assert_eq!(nn, vec2(-pp.x, -pp.y));
+					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: vec2(x, y),
+								pos: vec3(x, y, z),
-							rot: Mat2::IDENTITY
+								rot: Mat3::IDENTITY
 							})
 							.pos
 							.y,
@@ -447,13 +476,15 @@ mod tests {
 					);
 				}
 			}
+		}
 		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: vec2(x, y),
+								pos: vec3(x, y, z),
-							rot: Mat2::IDENTITY
+								rot: Mat3::IDENTITY
 							})
 							.pos
 							.y,
@@ -461,13 +492,15 @@ mod tests {
 					);
 				}
 			}
+		}
 		// 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: vec2(x, y),
+							pos: vec3(x, y, z),
-						rot: Mat2::IDENTITY,
+							rot: Mat3::IDENTITY,
 						})
 						.pos
 						.y;
@@ -476,4 +509,5 @@ mod tests {
 				}
 			}
 		}
+	}
 }

src/tube/metric.rs

@@ -1,8 +1,8 @@
-use glam::{f32, vec2, Mat2, Vec2};
+use glam::{f32, vec3, Mat3, Vec3};
 
 use crate::fns::{self, Limiter};
-use crate::mathx::Decomp2;
+use crate::mathx::Decomp3;
-use crate::riemann::{Metric, Tens2};
+use crate::riemann::{Metric, Tens3};
 
 #[derive(Copy, Clone, Debug)]
 pub struct Tube {
@@ -41,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 {
+		Decomp3 {
-			ortho: Mat2::IDENTITY,
+			ortho: Mat3::IDENTITY,
-			diag: vec2(1.0, s),
+			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;
@@ -63,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]),
+			Mat3::from_cols_array(&[0., 0., 0., 0., ds2_dx, 0., 0., 0., 0.]),
-			Mat2::from_cols_array(&[0.0, 0.0, 0.0, ds2_dy]),
+			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.]),
 		]
 	}
 }
@@ -72,10 +73,10 @@ 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::mathx::Decomp2;
+	use crate::mathx::Decomp3;
 	use crate::riemann::Metric;
 	use crate::tube::Space;
 	use crate::types::Ray;
@@ -86,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)
 			}
 		}
@@ -100,21 +101,28 @@ 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);
+				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));
+					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() {
@@ -132,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 a = vec3(ax, -(space.tube.external_halflength + off), 0.);
-				let b = vec2(bx, space.tube.external_halflength + off);
+				let b = vec3(bx, space.tube.external_halflength + off, 0.);
-				let Δ = vec2(bx - ax, 2.0 * (space.tube.internal_halflength + off));
+				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);
@@ -170,9 +178,9 @@ mod test {
 				space.tube.inner_radius - ε,
 				20,
 			) {
-				let a = vec2(ax, -(space.tube.external_halflength + off));
+				let a = vec3(ax, -(space.tube.external_halflength + off), 0.);
-				let b = vec2(bx, space.tube.external_halflength + off);
+				let b = vec3(bx, space.tube.external_halflength + off, 0.);
-				let Δ = vec2(bx - ax, 2.0 * (space.tube.internal_halflength + off));
+				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);
@@ -198,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 a = vec3(x, -(space.tube.external_halflength + off), 0.);
-			let b = vec2(x, space.tube.external_halflength + off);
+			let b = vec3(x, space.tube.external_halflength + off, 0.);
-			let Δ = vec2(0.0, 2.0 * (space.tube.internal_halflength + off));
+			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);
@@ -225,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 a = vec3(x, -(space.tube.external_halflength + off), 0.);
-			let b = vec2(x, space.tube.external_halflength + off);
+			let b = vec3(x, space.tube.external_halflength + off, 0.);
-			let Δ = vec2(0.0, 2.0 * (space.tube.external_halflength + off));
+			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),
+			pos: vec3(2.0, 3.0, 2.0),
-			dir: vec2(4.0, 5.0)
+			dir: vec3(4.0, 5.0, 4.0)
 		}),
 		Ray {
-			pos: vec2(2.0, 3.0),
+			pos: vec3(2.0, 3.0, 2.0),
-			dir: vec2(4.0, 5.0)
+			dir: vec3(4.0, 5.0, 4.0)
 		}
 	);
 	assert_eq!(
 		Rect::flip_ray(Ray {
-			pos: vec2(2.0, 3.0),
+			pos: vec3(2.0, 3.0, 2.0),
-			dir: vec2(-4.0, 5.0)
+			dir: vec3(-4.0, 5.0, -4.0)
 		}),
 		Ray {
-			pos: vec2(-2.0, 3.0),
+			pos: vec3(-2.0, 3.0, -2.0),
-			dir: vec2(4.0, 5.0)
+			dir: vec3(4.0, 5.0, 4.0)
 		}
 	);
 	assert_eq!(
 		Rect::flip_ray(Ray {
-			pos: vec2(2.0, 3.0),
+			pos: vec3(2.0, 3.0, 2.0),
-			dir: vec2(4.0, -5.0)
+			dir: vec3(4.0, -5.0, 4.0)
 		}),
 		Ray {
-			pos: vec2(2.0, -3.0),
+			pos: vec3(2.0, -3.0, 2.0),
-			dir: vec2(4.0, 5.0)
+			dir: vec3(4.0, 5.0, 4.0)
 		}
 	);
 	assert_eq!(
 		Rect::flip_ray(Ray {
-			pos: vec2(2.0, 3.0),
+			pos: vec3(2.0, 3.0, 2.0),
-			dir: vec2(4.0, 0.0)
+			dir: vec3(4.0, 0.0, 4.0)
 		}),
 		Ray {
-			pos: vec2(2.0, 3.0),
+			pos: vec3(2.0, 3.0, 2.0),
-			dir: vec2(4.0, 0.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),
+			pos: vec3(3.0, 3.0, 3.0),
-			dir: vec2(1.0, 1.0)
+			dir: vec3(1.0, 1.0, 1.0)
 		}),
 		None
 	);
 	assert_eq!(
 		r.trace_into(Ray {
-			pos: vec2(-3.0, 2.0),
+			pos: vec3(-3.0, 2.0, -3.0),
-			dir: vec2(1.0, 0.0)
+			dir: vec3(1.0, 0.0, 1.0)
 		}),
 		Some(1.0)
 	);
 	assert_eq!(
 		r.trace_into(Ray {
-			pos: vec2(-3.0, 2.0),
+			pos: vec3(-3.0, 2.0, -3.0),
-			dir: vec2(-1.0, 0.0)
+			dir: vec3(-1.0, 0.0, -1.0)
 		}),
 		None
 	);
 	assert_eq!(
 		r.trace_into(Ray {
-			pos: vec2(-3.0, 1.0),
+			pos: vec3(-3.0, 1.0, -3.0),
-			dir: vec2(2.0, 2.0)
+			dir: vec3(2.0, 2.0, 2.0)
 		}),
 		Some(0.5)
 	);
 	assert_eq!(
 		r.trace_into(Ray {
-			pos: vec2(-3.0, 2.1),
+			pos: vec3(-3.0, 2.1, -3.0),
-			dir: vec2(2.0, 2.0)
+			dir: vec3(2.0, 2.0, 2.0)
 		}),
 		None
 	);
 
 	assert_eq!(
 		r.trace_into(Ray {
-			pos: vec2(2.0, 3.0),
+			pos: vec3(2.0, 3.0, 2.0),
-			dir: vec2(1.0, 1.0)
+			dir: vec3(1.0, 1.0, 1.0)
 		}),
 		None
 	);
 	assert_eq!(
 		r.trace_into(Ray {
-			pos: vec2(-2.0, 3.0),
+			pos: vec3(-2.0, 3.0, -2.0),
-			dir: vec2(-1.0, 1.0)
+			dir: vec3(-1.0, 1.0, -1.0)
 		}),
 		None
 	);
 	assert_eq!(
 		r.trace_into(Ray {
-			pos: vec2(2.0, 3.0),
+			pos: vec3(2.0, 3.0, 2.0),
-			dir: vec2(-1.0, -1.0)
+			dir: vec3(-1.0, -1.0, -1.0)
 		}),
 		Some(0.0)
 	);
 	assert_eq!(
 		r.trace_into(Ray {
-			pos: vec2(2.0, -3.0),
+			pos: vec3(2.0, -3.0, 2.0),
-			dir: vec2(-1.0, 1.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),
+			pos: vec3(0.0, 0.0, 0.0),
-			dir: vec2(1.0, 1.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),
+			pos: vec3(0.0, 0.0, 0.0),
-			dir: vec2(0.0, 1.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),
+			pos: vec3(0.0, 1.0, 0.0),
-			dir: vec2(0.0, -1.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),
+			pos: vec3(1.0, 1.0, 1.0),
-			dir: vec2(0.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),
+			pos: vec3(2.0, 3.0, 2.0),
-			dir: vec2(1.0, 1.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 pos: Vec3,
-	pub dir: Vec2,
+	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,  // в локальной ортонормированной СК объекта
 }