Make the tube actually cylindrical

src/tube/coords.rs

@@ -3,7 +3,7 @@ use glam::{vec3, Mat3, Vec3};
 use crate::riemann::Metric;
 use crate::types::{Location, Ray};
 
-use super::{Rect, Tube};
+use super::{Tube, YCylinder};
 
 pub trait FlatCoordinateSystem<T> {
 	fn flat_to_global(&self, v: T) -> T;
@@ -80,8 +80,9 @@ impl FlatCoordinateSystem<Vec3> for InnerCS {
 
 impl FlatRegion for InnerCS {
 	fn distance_to_boundary(&self, ray: Ray) -> Option<f32> {
-		Rect {
+		YCylinder {
-			size: vec3(self.0.inner_radius, self.0.internal_halflength, self.0.inner_radius),
+			radius: self.0.inner_radius,
+			half_length: self.0.internal_halflength,
 		}
 		.trace_out_of(ray)
 	}
@@ -97,12 +98,9 @@ impl MetricCS for OuterCS {
 
 impl FlatCoordinateSystem<Vec3> for OuterCS {
 	fn flat_to_global(&self, pos: Vec3) -> Vec3 {
-		let inner = Rect {
+		let inner = YCylinder {
-			size: vec3(
+			radius: self.0.inner_radius + 1.0,
-				self.0.inner_radius + 1.0,
+			half_length: self.0.external_halflength,
-				self.0.external_halflength,
-				self.0.inner_radius + 1.0,
-			),
 		};
 		if inner.is_inside(pos) {
 			let Vec3 { x, y: v, z } = pos;
@@ -116,12 +114,9 @@ impl FlatCoordinateSystem<Vec3> for OuterCS {
 	}
 
 	fn global_to_flat(&self, pos: Vec3) -> Vec3 {
-		let inner = Rect {
+		let inner = YCylinder {
-			size: vec3(
+			radius: self.0.inner_radius + 1.0,
-				self.0.inner_radius + 1.0,
+			half_length: self.0.external_halflength,
-				self.0.external_halflength,
-				self.0.inner_radius + 1.0,
-			),
 		};
 		if inner.is_inside(pos) {
 			let Vec3 { x: u, y, z: w } = pos; // в основной СК
@@ -135,8 +130,9 @@ impl FlatCoordinateSystem<Vec3> for OuterCS {
 
 impl FlatRegion for OuterCS {
 	fn distance_to_boundary(&self, ray: Ray) -> Option<f32> {
-		Rect {
+		YCylinder {
-			size: vec3(self.0.outer_radius, self.0.external_halflength, self.0.outer_radius),
+			radius: self.0.outer_radius,
+			half_length: self.0.external_halflength,
 		}
 		.trace_into(ray)
 	}
@@ -346,15 +342,16 @@ mod tests {
 			external_halflength: 300.0,
 		});
 		// TODO replace 200.20016 with something sane
+		// NOTE it can’t test −30..30 as that area intersects the boundary
 		test_flat_region(
 			&mapper,
-			(vec3(-30.0, -300.0, -30.0), vec3(30.0, -1.0, 30.0)),
+			(vec3(-20.0, -300.0, -20.0), vec3(20.0, -1.0, 20.0)),
-			(vec3(-30.0, -300.0, -30.0), vec3(30.0, -200.20016, 30.0)),
+			(vec3(-20.0, -300.0, -20.0), vec3(20.0, -200.20016, 20.0)),
 		);
 		test_flat_region(
 			&mapper,
-			(vec3(-30.0, 1.0, -30.0), vec3(30.0, 300.0, 30.0)),
+			(vec3(-20.0, 1.0, -20.0), vec3(20.0, 300.0, 20.0)),
-			(vec3(-30.0, 200.20016, -30.0), vec3(30.0, 300.0, 30.0)),
+			(vec3(-20.0, 200.20016, -20.0), vec3(20.0, 300.0, 20.0)),
 		);
 		test_flat_region(
 			&mapper,
@@ -451,9 +448,9 @@ mod tests {
 			}
 		}
 		// accelerating
-		for x in linspace(-29., 29., 20) {
+		for x in linspace(-21., 21., 20) {
 			for y in linspace(1., 299., 20) {
-				for z in linspace(-29., 29., 20) {
+				for z in linspace(-21., 21., 20) {
 					let v = mapper
 						.global_to_flat(Location {
 							pos: vec3(x, y, z),

src/tube/metric.rs

@@ -42,10 +42,11 @@ impl Tube {
 
 impl Metric for Tube {
 	fn sqrt_at(&self, pos: Vec3) -> Decomp3 {
-		let sx = self.fx().value(pos.x);
+		let r = f32::hypot(pos.x, pos.z);
+		let sr = self.fx().value(r);
 		let sy = self.fy().du(pos.y);
-		let s = sx + sy - sx * sy;
+		let s = sr + sy - sr * sy;
-		assert!(sx.is_finite());
+		assert!(sr.is_finite());
 		assert!(sy.is_finite());
 		assert!(sy > 0.0);
 		Decomp3 {
@@ -55,17 +56,23 @@ impl Metric for Tube {
 	}
 
 	fn part_derivs_at(&self, pos: Vec3) -> Tens3 {
-		let sx = self.fx().value(pos.x);
+		let r = f32::hypot(pos.x, pos.z);
+		let sr = self.fx().value(r);
 		let sy = self.fy().du(pos.y);
-		let s = sx + sy - sx * sy;
+		let s = sr + sy - sr * sy;
-		let dsx_dx = self.fx().derivative(pos.x);
+		let dsr_dr = self.fx().derivative(r);
 		let dsy_dy = self.fy().d2u(pos.y);
-		let ds2_dx = 2.0 * s * (1.0 - sy) * dsx_dx;
+		let ds2_dr = 2.0 * s * (1.0 - sy) * dsr_dr;
-		let ds2_dy = 2.0 * s * (1.0 - sx) * dsy_dy;
+		let ds2_dy = 2.0 * s * (1.0 - sr) * dsy_dy;
+		let (ds2_dx, ds2_dz) = if ds2_dr.abs() < 1e-5 {
+			(0., 0.)
+		} else {
+			(ds2_dr * pos.x / r, ds2_dr * pos.z / r)
+		};
 		[
 			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.]),
+			Mat3::from_cols_array(&[0., 0., 0., 0., ds2_dz, 0., 0., 0., 0.]),
 		]
 	}
 }

src/tube/mod.rs

@@ -28,10 +28,12 @@ pub enum Subspace {
 impl Space {
 	fn which_subspace(&self, pt: Vec3) -> Subspace {
 		if pt.y.abs() > self.tube.external_halflength {
+			return Outer;
+		}
+		let r = f32::hypot(pt.x, pt.z);
+		if r > self.tube.outer_radius {
 			Outer
-		} else if pt.x.abs() > self.tube.outer_radius {
+		} else if r > self.tube.inner_radius {
-			Outer
-		} else if pt.x.abs() > self.tube.inner_radius {
 			Boundary
 		} else {
 			Inner