Make the tube actually cylindrical

rustfmt.toml

@@ -1 +1,2 @@
 hard_tabs = true
+max_width = 120

src/bin/flat/main.rs

@@ -52,20 +52,12 @@ pub fn main() {
 			let cam1 = put_object(&space.tube, vec3(-500., 0., 0.), Mat3::IDENTITY);
 			let cam2 = put_object(
 				&space.tube,
-				vec3(
+				vec3(-2.5 * tube.outer_radius, 1.25 * tube.external_halflength, 0.),
-					-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,
-				vec3(
+				vec3(0.25 * tube.inner_radius, 0.25 * tube.external_halflength, 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.)),
 			);
 
@@ -134,11 +126,7 @@ pub fn main() {
 							let dir = Vec2::from_angle(φ);
 							let dir = vec3(dir.x, dir.y, 0.);
 							let dir = obj.loc.rot * dir * d;
-							space
+							space.trace_iter(Ray { pos, dir }).nth(n as usize).unwrap().pos
-								.trace_iter(Ray { pos, dir })
-								.nth(n as usize)
-								.unwrap()
-								.pos
 						}),
 				);
 			}
@@ -169,10 +157,7 @@ fn draw_ray_2(gc: &mut Vec<Draw>, space: &Space, camera: Location, dir: Vec3) {
 	for pt in &path.points[1..] {
 		gc.line_to(pt.pos.x, pt.pos.y);
 	}
-	let end_pos = *path
+	let end_pos = *path.points.last().expect("the starting point is always in the path");
-		.points
-		.last()
-		.expect("the starting point is always in the path");
 	let dir_pos = end_pos.forward(1000. / DT).pos;
 	gc.line_to(dir_pos.x, dir_pos.y);
 	gc.stroke();
@@ -192,11 +177,7 @@ fn draw_track(gc: &mut Vec<Draw>, space: &Space, start: Vec2, dir: Vec2) {
 	// let v = space.tube.normalize(start, dir);
 	let mut loc = Location {
 		pos: vec3(start.x, start.y, 0.),
-		rot: mat3(
+		rot: mat3(vec3(dir.x, dir.y, 0.), vec3(-dir.y, dir.x, 0.), vec3(0., 0., 1.)),
-			vec3(dir.x, dir.y, 0.),
-			vec3(-dir.y, dir.x, 0.),
-			vec3(0., 0., 1.),
-		),
 	};
 	let v = vec3(1., 0., 0.);
 	let mut draw = |loc: &Location| {

src/bin/mesh/main.rs

@@ -81,9 +81,7 @@ fn render(mesh: &Mesh, camera: impl Fn(Vec2) -> (Vec3, Vec3)) -> Image {
 			} else {
 				bkg
 			};
-			let color = (color * 255.0)
+			let color = (color * 255.0).as_ivec3().clamp(IVec3::splat(0), IVec3::splat(255));
-				.as_ivec3()
-				.clamp(IVec3::splat(0), IVec3::splat(255));
 			img.put_pixel(x, y, Color(color.x as u8, color.y as u8, color.z as u8));
 		}
 	}
@@ -141,11 +139,7 @@ fn main() -> Result<(), Box<dyn Error>> {
 	loop {
 		for phi in 0..360 {
 			let proj = PROJS[PROJ_INDEX.load(Relaxed)];
-			let m_view = ypr_to_mat(vec3(
+			let m_view = ypr_to_mat(vec3((135.0 + phi as f32) * PI / 180.0, -30.0 * PI / 180.0, 0.0f32));
-				(135.0 + phi as f32) * PI / 180.0,
-				-30.0 * PI / 180.0,
-				0.0f32,
-			));
 			let m_camera = m_view.transpose();
 			let img = render(mesh.as_slice(), |off| {
 				let (base, ray) = proj(40., 20. * off);

src/bin/wireframe/camera.rs

@@ -64,7 +64,7 @@ impl Camera {
 			vec4(0., 0., 0., 1.),
 		);
 		let size = view_size.normalize() * std::f32::consts::SQRT_2;
-		let proj = make_proj_matrix(vec3(size.x, size.y, 2.), (1., 4096.)) * M;
+		let proj = make_proj_matrix(vec3(size.x, size.y, 2.), (1., (2f32).powi(16) + 1.)) * M;
 
 		let mvp = proj * view_mtx;
 		CameraUniform {

src/bin/wireframe/lines.rs

@@ -112,14 +112,12 @@ impl LineRenderer {
 							format: wgpu::VertexFormat::Float32x3,
 						},
 						wgpu::VertexAttribute {
-							offset: (mem::size_of::<Vertex>() + mem::offset_of!(Vertex, position))
+							offset: (mem::size_of::<Vertex>() + mem::offset_of!(Vertex, position)) as u64,
-								as u64,
 							shader_location: 2,
 							format: wgpu::VertexFormat::Float32x3,
 						},
 						wgpu::VertexAttribute {
-							offset: (mem::size_of::<Vertex>() + mem::offset_of!(Vertex, tangent))
+							offset: (mem::size_of::<Vertex>() + mem::offset_of!(Vertex, tangent)) as u64,
-								as u64,
 							shader_location: 3,
 							format: wgpu::VertexFormat::Float32x3,
 						},

src/bin/wireframe/main.rs

@@ -164,8 +164,7 @@ impl<'a> State<'a> {
 			.await
 			.unwrap();
 
-		let viewport =
+		let viewport = viewport::Viewport::new(&adapter, &device, surface, uvec2(size.width, size.height));
-			viewport::Viewport::new(&adapter, &device, surface, uvec2(size.width, size.height));
 
 		let kbd = keyctl::Keyboard::new();
 		let cam_loc = camctl::CameraLocation::new();
@@ -179,13 +178,7 @@ impl<'a> State<'a> {
 		};
 
 		let cam_obj = camera::Camera::new(&device);
-		let line_rend = lines::LineRenderer::new(
+		let line_rend = lines::LineRenderer::new(&device, cam_obj.bind_group_layout(), viewport.format(), depth, msaa);
-			&device,
-			cam_obj.bind_group_layout(),
-			viewport.format(),
-			depth,
-			msaa,
-		);
 
 		let scene = prepare_scene(&device);
 
@@ -226,10 +219,8 @@ impl<'a> State<'a> {
 		};
 		let size = self.viewport.size().as_vec2();
 
-		self.cam_loc
+		self.cam_loc.move_rel(100. * dt * self.kbd.control(&KEYS_MOVE));
-			.move_rel(100. * dt * self.kbd.control(&KEYS_MOVE));
+		self.cam_loc.rotate_rel_ypr(2. * dt * self.kbd.control(&KEYS_ROTATE));
-		self.cam_loc
-			.rotate_rel_ypr(2. * dt * self.kbd.control(&KEYS_ROTATE));
 		self.cam_obj.set(&self.queue, self.cam_loc.view_mtx(), size);
 	}
 
@@ -239,11 +230,8 @@ impl<'a> State<'a> {
 			.set_title(&format!("Space Refraction ({:.1} FPS)", self.fps.get()));
 		self.viewport
 			.render_single_pass(&self.device, &self.queue, |mut render_pass| {
-				self.line_rend.render(
+				self.line_rend
-					&mut render_pass,
+					.render(&mut render_pass, self.cam_obj.bind_group(), self.scene.iter());
-					self.cam_obj.bind_group(),
-					self.scene.iter(),
-				);
 			})
 	}
 }
@@ -262,10 +250,7 @@ pub async fn run() {
 	event_loop
 		.run(move |event, control_flow| {
 			match event {
-				Event::WindowEvent {
+				Event::WindowEvent { ref event, window_id } if window_id == state.window().id() => {
-					ref event,
-					window_id,
-				} if window_id == state.window().id() => {
 					match event {
 						WindowEvent::KeyboardInput {
 							device_id: _,

src/bin/wireframe/scene.rs

@@ -40,12 +40,7 @@ fn draw_rect(center: Vec3, u: Vec3, v: Vec3) -> Vec<Line> {
 	let b = center + u - v;
 	let c = center + u + v;
 	let d = center - u + v;
-	vec![
+	vec![draw_line(a, b), draw_line(b, c), draw_line(c, d), draw_line(d, a)]
-		draw_line(a, b),
-		draw_line(b, c),
-		draw_line(c, d),
-		draw_line(d, a),
-	]
 }
 
 fn draw_ellipse(center: Vec3, u: Vec3, v: Vec3) -> Line {
@@ -88,28 +83,21 @@ pub fn build() -> Vec<FancyLine> {
 	let cam1 = put_object(&space.tube, vec3(-500., 0., 0.), Mat3::IDENTITY);
 	let cam2 = put_object(
 		&space.tube,
-		vec3(
+		vec3(-2.5 * tube.outer_radius, 1.25 * tube.external_halflength, 0.),
-			-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,
-		vec3(
+		vec3(0.25 * tube.inner_radius, 0.25 * tube.external_halflength, 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.)),
 	);
 
 	let mut gc = vec![];
 	paint(&mut gc, vec3(0.6, 0.6, 0.6), tube.render());
-	paint(&mut gc, vec3(0.0, 0.6, 1.0), draw_fan_2(&space, cam3, 1.0));
+	paint(&mut gc, vec3(0.0, 0.6, 1.0), draw_fan_2(&space, cam3, vec3(0., 1., 0.)));
-	paint(&mut gc, vec3(0.2, 1.0, 0.0), draw_fan_2(&space, cam2, 1.0));
+	paint(&mut gc, vec3(0.2, 1.0, 0.0), draw_fan_2(&space, cam2, vec3(0., 1., 0.)));
-	paint(&mut gc, vec3(1.0, 0.2, 0.0), draw_fan_2(&space, cam1, 1.0));
+	paint(&mut gc, vec3(0.0, 1.0, 0.6), draw_fan_2(&space, cam2, vec3(0., 0., 1.)));
+	paint(&mut gc, vec3(1.0, 0.2, 0.0), draw_fan_2(&space, cam1, vec3(0., 1., 0.)));
 	gc
 }
 
@@ -124,18 +112,15 @@ fn draw_ray_2(gc: &mut Vec<Line>, space: &Space, camera: Location, dir: Vec3) {
 	}
 
 	let mut pts = path.points;
-	let end_pos = *pts
+	let end_pos = *pts.last().expect("the starting point is always in the path");
-		.last()
+	pts.extend(itertools::iterate(end_pos, |r| r.forward(100.0)).take(1000));
-		.expect("the starting point is always in the path");
-	let dir_pos = end_pos.forward(10000.0);
-	pts.push(dir_pos);
 	gc.push(Line::Strip(pts));
 }
 
-fn draw_fan_2(space: &Space, camera: Location, spread: f32) -> Vec<Line> {
+fn draw_fan_2(space: &Space, camera: Location, spread: Vec3) -> Vec<Line> {
 	let mut gc = vec![];
-	for y in itertools_num::linspace(-spread, spread, 101) {
+	for δ in itertools_num::linspace(-1., 1., 101) {
-		draw_ray_2(&mut gc, space, camera, vec3(1., y, 0.));
+		draw_ray_2(&mut gc, space, camera, vec3(1., 0., 0.) + δ * spread);
 	}
 	gc
 }

src/bin/wireframe/viewport.rs

@@ -8,12 +8,7 @@ pub struct Viewport<'a> {
 }
 
 impl<'a> Viewport<'a> {
-	pub fn new(
+	pub fn new(adapter: &wgpu::Adapter, device: &wgpu::Device, surface: wgpu::Surface<'a>, size: UVec2) -> Self {
-		adapter: &wgpu::Adapter,
-		device: &wgpu::Device,
-		surface: wgpu::Surface<'a>,
-		size: UVec2,
-	) -> Self {
 		let caps = surface.get_capabilities(adapter);
 		let format = wgpu::TextureFormat::Bgra8UnormSrgb;
 		let sample_count = adapter
@@ -73,9 +68,7 @@ impl<'a> Viewport<'a> {
 		f: impl FnOnce(wgpu::RenderPass),
 	) -> Result<(), wgpu::SurfaceError> {
 		let output = self.surface.get_current_texture()?;
-		let view = output
+		let view = output.texture.create_view(&wgpu::TextureViewDescriptor::default());
-			.texture
-			.create_view(&wgpu::TextureViewDescriptor::default());
 		let mut encoder = device.create_command_encoder(&wgpu::CommandEncoderDescriptor {
 			label: Some("Render CommandEncoder"),
 		});
@@ -110,12 +103,7 @@ struct Multisample {
 }
 
 impl Multisample {
-	fn new(
+	fn new(device: &wgpu::Device, format: wgpu::TextureFormat, size: UVec2, sample_count: u32) -> Multisample {
-		device: &wgpu::Device,
-		format: wgpu::TextureFormat,
-		size: UVec2,
-		sample_count: u32,
-	) -> Multisample {
 		let tex = device.create_texture(&wgpu::TextureDescriptor {
 			label: Some("Multisample texture"),
 			size: wgpu::Extent3d {

src/fns.rs

@@ -82,12 +82,7 @@ impl QuadraticAccelerator {
 	}
 
 	pub fn x(&self, u: f32) -> f32 {
-		extend_linear(
+		extend_linear(u, |u| (self.a() * u.abs() + self.b()) * u, self.internal, self.external)
-			u,
-			|u| (self.a() * u.abs() + self.b()) * u,
-			self.internal,
-			self.external,
-		)
 	}
 	pub fn u(&self, x: f32) -> f32 {
 		extend_linear(
@@ -98,12 +93,7 @@ impl QuadraticAccelerator {
 		)
 	}
 	pub fn dx(&self, u: f32) -> f32 {
-		extend_const(
+		extend_const(u, |u| 2.0 * self.a() * u.abs() + self.b(), self.internal, 1.0)
-			u,
-			|u| 2.0 * self.a() * u.abs() + self.b(),
-			self.internal,
-			1.0,
-		)
 	}
 	pub fn du(&self, x: f32) -> f32 {
 		extend_const(x, |x| 1.0 / self.root(x), self.external, 1.0)
@@ -148,28 +138,12 @@ mod test {
 
 	#[test]
 	fn test_smoothstep_limiter() {
-		test_limiter(
+		test_limiter(SmoothstepLimiter { min: 20.0, max: 30.0 }, 20.0, 30.0, 1.0 / 32.0);
-			SmoothstepLimiter {
-				min: 20.0,
-				max: 30.0,
-			},
-			20.0,
-			30.0,
-			1.0 / 32.0,
-		);
 	}
 
 	#[test]
 	fn test_smootherstep_limiter() {
-		test_limiter(
+		test_limiter(SmootherstepLimiter { min: 20.0, max: 30.0 }, 20.0, 30.0, 1.0 / 32.0);
-			SmootherstepLimiter {
-				min: 20.0,
-				max: 30.0,
-			},
-			20.0,
-			30.0,
-			1.0 / 32.0,
-		);
 	}
 
 	#[test]

src/mathx.rs

@@ -36,10 +36,8 @@ impl MatExt for Mat3 {
 	fn orthonormalize(&self) -> Self {
 		let fx = self.x_axis.normalize();
 		let fy = (self.y_axis - self.y_axis.project_onto_normalized(fx)).normalize();
-		let fz = (self.z_axis
+		let fz = (self.z_axis - self.z_axis.project_onto_normalized(fx) - self.z_axis.project_onto_normalized(fy))
-			- self.z_axis.project_onto_normalized(fx)
+			.normalize();
-			- self.z_axis.project_onto_normalized(fy))
-		.normalize();
 		Self::from_cols(fx, fy, fz)
 	}
 }

src/mesh_loader.rs

@@ -35,10 +35,7 @@ impl ObjMesh {
 	}
 
 	fn parse_fv(desc: &&str) -> ObjVertex {
-		let tokens: Vec<_> = desc
+		let tokens: Vec<_> = desc.split('/').map(|s| s.parse::<usize>().unwrap() - 1).collect();
-			.split('/')
-			.map(|s| s.parse::<usize>().unwrap() - 1)
-			.collect();
 		assert_eq!(tokens.len(), 3);
 		ObjVertex {
 			vertex: tokens[0],

src/mesh_tracer.rs

@@ -8,21 +8,13 @@ pub struct TraceResult {
 	pub normal: Vec3,
 }
 
-pub fn trace_to_mesh_all(
+pub fn trace_to_mesh_all(mesh: &Mesh, base: Vec3, ray: Vec3) -> impl Iterator<Item = TraceResult> + '_ {
-	mesh: &Mesh,
-	base: Vec3,
-	ray: Vec3,
-) -> impl Iterator<Item = TraceResult> + '_ {
 	mesh.iter().filter_map(move |f| {
 		let fs = (0..3).map(|k| edge_dist(f.vertices[k], f.vertices[(k + 1) % 3], base, ray));
 		if fs.into_iter().any(|f| f < 0.0) {
 			return None;
 		}
-		let m = mat3(
+		let m = mat3(f.vertices[1] - f.vertices[0], f.vertices[2] - f.vertices[0], -ray);
-			f.vertices[1] - f.vertices[0],
-			f.vertices[2] - f.vertices[0],
-			-ray,
-		);
 		let m = m.inverse();
 		let rel = m * (base - f.vertices[0]);
 		Some(TraceResult {

src/riemann.rs

@@ -15,7 +15,7 @@ pub trait Metric {
 	}
 
 	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
+		part_deriv(|p| self.at(p), pos, 1.0 / 64.0) // there just isn’t enough precision for a small ε.
 	}
 
 	fn vec_length_at(&self, at: Vec3, v: Vec3) -> f32 {
@@ -58,7 +58,7 @@ pub fn krist(space: &impl Metric, pos: Vec3) -> Tens3 {
 	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)
+		0.5 * (0..3)
 			.map(|m| g.col(m)[i] * (d[l].col(k)[m] + d[k].col(m)[l] - d[m].col(k)[l]))
 			.sum::<f32>()
 	})
@@ -161,10 +161,7 @@ mod tests {
 			metric.inverse_at(rng.gen()),
 			Mat3::from_cols_array(&[1. / 9., 0., 0., 0., 1. / 16., 0., 0., 0., 1. / 25.])
 		);
-		assert_eq!(
+		assert_eq!(metric.part_derivs_at(rng.gen()), [Mat3::ZERO, Mat3::ZERO, Mat3::ZERO]);
-			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.);
@@ -225,14 +222,9 @@ mod tests {
 			vec3(3., 6.25, 0.)
 		);
 		assert_abs_diff_eq!(
-			trace_iter(
+			trace_iter(&metric, vec3(3., 5., 0.), vec3(0.5, 0.25, 0.), std::f32::consts::SQRT_2)
-				&metric,
+				.nth(7)
-				vec3(3., 5., 0.),
+				.unwrap(),
-				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

@@ -3,16 +3,14 @@ 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;
 	fn global_to_flat(&self, v: T) -> T;
 }
 
-pub trait FlatRegion:
+pub trait FlatRegion: FlatCoordinateSystem<Vec3> + FlatCoordinateSystem<Ray> + FlatCoordinateSystem<Location> {
-	FlatCoordinateSystem<Vec3> + FlatCoordinateSystem<Ray> + FlatCoordinateSystem<Location>
-{
 	// Измеряет расстояние до выхода за пределы области вдоль луча ray. Луч задаётся в плоской СК.
 	fn distance_to_boundary(&self, _ray: Ray) -> Option<f32> {
 		None
@@ -22,10 +20,7 @@ pub trait FlatRegion:
 trait MetricCS: FlatCoordinateSystem<Vec3> {
 	fn global_metric(&self) -> &impl Metric;
 	fn flat_to_global_tfm_at(&self, pos: Vec3) -> Mat3 {
-		self.global_metric()
+		self.global_metric().sqrt_at(self.flat_to_global(pos)).inverse().into()
-			.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()
@@ -85,12 +80,9 @@ impl FlatCoordinateSystem<Vec3> for InnerCS {
 
 impl FlatRegion for InnerCS {
 	fn distance_to_boundary(&self, ray: Ray) -> Option<f32> {
-		Rect {
+		YCylinder {
-			size: vec3(
+			radius: self.0.inner_radius,
-				self.0.inner_radius,
+			half_length: self.0.internal_halflength,
-				self.0.internal_halflength,
-				self.0.inner_radius,
-			),
 		}
 		.trace_out_of(ray)
 	}
@@ -106,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;
@@ -125,17 +114,13 @@ 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; // в основной СК
-			let v = self.0.v(y)
+			let v = self.0.v(y) + y.signum() * (self.0.external_halflength - self.0.internal_halflength);
-				+ y.signum() * (self.0.external_halflength - self.0.internal_halflength);
 			vec3(u, v, w) // в плоском продолжении СК Outer на область Inner
 		} else {
 			pos
@@ -145,12 +130,9 @@ impl FlatCoordinateSystem<Vec3> for OuterCS {
 
 impl FlatRegion for OuterCS {
 	fn distance_to_boundary(&self, ray: Ray) -> Option<f32> {
-		Rect {
+		YCylinder {
-			size: vec3(
+			radius: self.0.outer_radius,
-				self.0.outer_radius,
+			half_length: self.0.external_halflength,
-				self.0.external_halflength,
-				self.0.outer_radius,
-			),
 		}
 		.trace_into(ray)
 	}
@@ -227,20 +209,12 @@ mod tests {
 			}),
 			Location {
 				pos: vec3(0., -0.5, 0.),
-				rot: mat3(
+				rot: mat3(vec3(0., 0.25, 0.), vec3(-1. / 3., 0., 0.), vec3(0., 0., 0.2))
-					vec3(0., 0.25, 0.),
-					vec3(-1. / 3., 0., 0.),
-					vec3(0., 0., 0.2)
-				)
 			}
 		);
 	}
 
-	fn test_flat_region(
+	fn test_flat_region(region: &impl FlatRegion, range_global: (Vec3, Vec3), range_flat: (Vec3, Vec3)) {
-		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 {
@@ -256,14 +230,7 @@ mod tests {
 				);
 			};
 		}
-		fn check_range(
+		fn check_range(name_a: &str, a: Vec3, range_a: (Vec3, Vec3), name_b: &str, b: Vec3, range_b: (Vec3, Vec3)) {
-			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, ε) {
@@ -284,14 +251,7 @@ mod tests {
 				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(
+					check_range("global", pos_global, range_global, "flat", pos_flat, range_flat);
-						"global",
-						pos_global,
-						range_global,
-						"flat",
-						pos_flat,
-						range_flat,
-					);
 					assert_eq_at!(
 						pos_global,
 						region
@@ -321,14 +281,7 @@ mod tests {
 				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(
+					check_range("flat", pos_flat, range_flat, "global", pos_global, range_global);
-						"flat",
-						pos_flat,
-						range_flat,
-						"global",
-						pos_global,
-						range_global,
-					);
 					assert_eq_at!(
 						pos_flat,
 						region
@@ -389,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,
@@ -494,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.]),
 		]
 	}
 }
@@ -101,24 +108,17 @@ 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, 20)
+		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, 20) {
-			for y in itertools_num::linspace(
+				for z in itertools_num::linspace(-mul * testee.outer_radius, mul * testee.outer_radius, 20) {
-				-mul * testee.external_halflength,
-				mul * testee.external_halflength,
-				20,
-			) {
-				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 =
+					let eq = (0..3).all(|coord| computed[coord].abs_diff_eq(reference[coord], epsilon));
-						(0..2).all(|coord| computed[coord].abs_diff_eq(reference[coord], epsilon));
+					assert!(
-					assert!(eq, "Bad derivative computation at {pos}:\n  explicit:  {computed:?}\n  numerical: {reference:?}\n");
+						eq,
+						"Bad derivative computation at {pos}:\n  explicit:  {computed:?}\n  numerical: {reference:?}\n"
+					);
 				}
 			}
 		}
@@ -168,16 +168,8 @@ mod test {
 		let ε = 1e-3;
 		let off = 10.0;
 		let steps = 4096;
-		for ax in linspace(
+		for ax in linspace(-space.tube.inner_radius + ε, space.tube.inner_radius - ε, 20) {
-			-space.tube.inner_radius + ε,
+			for bx in linspace(-space.tube.inner_radius + ε, space.tube.inner_radius - ε, 20) {
-			space.tube.inner_radius - ε,
-			20,
-		) {
-			for bx in linspace(
-				-space.tube.inner_radius + ε,
-				space.tube.inner_radius - ε,
-				20,
-			) {
 				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.);

src/tube/mod.rs

@@ -1,4 +1,4 @@
-use glam::{bool, f32, Mat3, Vec3};
+use glam::{bool, f32, Mat3, Vec3, Vec3Swizzles};
 
 use crate::ifaces::{DebugTraceable, RayPath, Traceable};
 use coords::{FlatCoordinateSystem, InnerCS, OuterCS};
@@ -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
@@ -50,8 +52,7 @@ impl Space {
 	/// Выполняет один шаг перемещения. Работает в любой части пространства.
 	/// off задаётся в локальной СК. Рекомендуется считать небольшими шагами.
 	pub fn move_step(&self, loc: Location, off: Vec3) -> Location {
-		let corr =
+		let corr = Mat3::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,
@@ -94,26 +95,16 @@ impl Space {
 	fn list_objects(&self, tfm: impl Fn(Location) -> Location) -> Vec<Object> {
 		self.objs
 			.iter()
-			.map(|&Object { id, loc, r }| Object {
+			.map(|&Object { id, loc, r }| Object { id, loc: tfm(loc), r })
-				id,
-				loc: tfm(loc),
-				r,
-			})
 			.collect()
 	}
 
-	fn hit_objects(
+	fn hit_objects(objs: &[Object], ray: Ray, limit: Option<f32>, globalize: impl Fn(Vec3) -> Vec3) -> Vec<Hit> {
-		objs: &[Object],
-		ray: Ray,
-		limit: Option<f32>,
-		globalize: impl Fn(Vec3) -> Vec3,
-	) -> Vec<Hit> {
 		let limit = limit.unwrap_or(f32::INFINITY);
 		objs.iter()
 			.filter_map(|obj| {
 				let rel = ray.pos - obj.loc.pos;
-				let diff = rel.dot(ray.dir).powi(2)
+				let diff = rel.dot(ray.dir).powi(2) - ray.dir.length_squared() * (rel.length_squared() - obj.r.powi(2));
-					- ray.dir.length_squared() * (rel.length_squared() - obj.r.powi(2));
 				if diff > 0.0 {
 					let t = (-rel.dot(ray.dir) - diff.sqrt()) / ray.dir.length_squared();
 					Some((obj, t))
@@ -277,154 +268,192 @@ impl Rect {
 	}
 }
 
-#[cfg(test)]
+fn solve_quadratic(a: f32, half_b: f32, c: f32) -> Option<(f32, f32)> {
-use glam::vec3;
+	let base = -half_b / a;
-
+	let d = base * base - c / a;
-#[test]
+	if d < 0. {
-fn test_rect() {
-	assert_eq!(
-		Rect::flip_ray(Ray {
-			pos: vec3(2.0, 3.0, 2.0),
-			dir: vec3(4.0, 5.0, 4.0)
-		}),
-		Ray {
-			pos: vec3(2.0, 3.0, 2.0),
-			dir: vec3(4.0, 5.0, 4.0)
-		}
-	);
-	assert_eq!(
-		Rect::flip_ray(Ray {
-			pos: vec3(2.0, 3.0, 2.0),
-			dir: vec3(-4.0, 5.0, -4.0)
-		}),
-		Ray {
-			pos: vec3(-2.0, 3.0, -2.0),
-			dir: vec3(4.0, 5.0, 4.0)
-		}
-	);
-	assert_eq!(
-		Rect::flip_ray(Ray {
-			pos: vec3(2.0, 3.0, 2.0),
-			dir: vec3(4.0, -5.0, 4.0)
-		}),
-		Ray {
-			pos: vec3(2.0, -3.0, 2.0),
-			dir: vec3(4.0, 5.0, 4.0)
-		}
-	);
-	assert_eq!(
-		Rect::flip_ray(Ray {
-			pos: vec3(2.0, 3.0, 2.0),
-			dir: vec3(4.0, 0.0, 4.0)
-		}),
-		Ray {
-			pos: vec3(2.0, 3.0, 2.0),
-			dir: vec3(4.0, 0.0, 4.0)
-		}
-	);
-
-	let r = Rect {
-		size: vec3(2.0, 3.0, 2.0),
-	};
-
-	assert_eq!(
-		r.trace_into(Ray {
-			pos: vec3(3.0, 3.0, 3.0),
-			dir: vec3(1.0, 1.0, 1.0)
-		}),
 		None
-	);
+	} else {
-	assert_eq!(
+		let δ = d.sqrt();
-		r.trace_into(Ray {
+		Some((base - δ, base + δ))
-			pos: vec3(-3.0, 2.0, -3.0),
+	}
-			dir: vec3(1.0, 0.0, 1.0)
+}
-		}),
+
-		Some(1.0)
+/// Цилиндр с центром в начале координат и осью вдоль оси Y.
-	);
+struct YCylinder {
-	assert_eq!(
+	pub half_length: f32,
-		r.trace_into(Ray {
+	pub radius: f32,
-			pos: vec3(-3.0, 2.0, -3.0),
+}
-			dir: vec3(-1.0, 0.0, -1.0)
+
-		}),
+impl YCylinder {
-		None
+	/// Отражает луч, чтобы все координаты направления были положительны (допустимо благодаря симметрии YCylinder).
-	);
+	fn flip_ray(ray: Ray) -> Ray {
-	assert_eq!(
+		Ray {
-		r.trace_into(Ray {
+			pos: ray.pos * ray.dir.signum(),
-			pos: vec3(-3.0, 1.0, -3.0),
+			dir: ray.dir.abs(),
-			dir: vec3(2.0, 2.0, 2.0)
+		}
-		}),
+	}
-		Some(0.5)
+
-	);
+	fn is_inside(&self, pt: Vec3) -> bool {
-	assert_eq!(
+		let r = f32::hypot(pt.x, pt.z);
-		r.trace_into(Ray {
+		pt.y.abs() < self.half_length && r < self.radius
-			pos: vec3(-3.0, 2.1, -3.0),
+	}
-			dir: vec3(2.0, 2.0, 2.0)
+
-		}),
+	fn trace_into(&self, ray: Ray) -> Option<f32> {
-		None
+		let ray = Self::flip_ray(ray);
-	);
+
-
+		// 1. ray.pos.y + t * ray.dir.y = −half_length
-	assert_eq!(
+		let t_cap_in = (-self.half_length - ray.pos.y) / ray.dir.y;
-		r.trace_into(Ray {
+		let t_cap_out = (self.half_length - ray.pos.y) / ray.dir.y;
-			pos: vec3(2.0, 3.0, 2.0),
+
-			dir: vec3(1.0, 1.0, 1.0)
+		// 2. (ray.pos.x + t * ray.dir.x)² + (ray.pos.z + t * ray.dir.z)² = radius²
-		}),
+		let pos = ray.pos.xz();
-		None
+		let dir = ray.dir.xz();
-	);
+		if dir.length_squared() < 1e-3 {
-	assert_eq!(
+			if pos.length_squared() >= self.radius.powi(2) {
-		r.trace_into(Ray {
+				return None;
-			pos: vec3(-2.0, 3.0, -2.0),
+			}
-			dir: vec3(-1.0, 1.0, -1.0)
+			return Some(t_cap_in).filter(|&t| t > 0.);
-		}),
+		}
-		None
+		let (t_side_in, t_side_out) = solve_quadratic(
-	);
+			dir.length_squared(),
-	assert_eq!(
+			pos.dot(dir),
-		r.trace_into(Ray {
+			pos.length_squared() - self.radius.powi(2),
-			pos: vec3(2.0, 3.0, 2.0),
+		)?;
-			dir: vec3(-1.0, -1.0, -1.0)
+		let t = f32::max(t_cap_in, t_side_in);
-		}),
+		if t < 0. {
-		Some(0.0)
+			return None;
-	);
+		}
-	assert_eq!(
+		if t >= t_cap_out || t >= t_side_out {
-		r.trace_into(Ray {
+			return None;
-			pos: vec3(2.0, -3.0, 2.0),
+		}
-			dir: vec3(-1.0, 1.0, -1.0)
+		Some(t)
-		}),
+	}
-		Some(0.0)
+
-	);
+	fn trace_out_of(&self, ray: Ray) -> Option<f32> {
-
+		let ray = Self::flip_ray(ray);
-	assert_eq!(
+		let t_cap_out = (self.half_length - ray.pos.y) / ray.dir.y;
-		r.trace_out_of(Ray {
+		let pos = ray.pos.xz();
-			pos: vec3(0.0, 0.0, 0.0),
+		let dir = ray.dir.xz();
-			dir: vec3(1.0, 1.0, 1.0)
+		if dir.length_squared() < 1e-3 {
-		}),
+			return Some(t_cap_out);
-		Some(2.0)
+		}
-	);
+		let (_t_side_in, t_side_out) = solve_quadratic(
-	assert_eq!(
+			dir.length_squared(),
-		r.trace_out_of(Ray {
+			pos.dot(dir),
-			pos: vec3(0.0, 0.0, 0.0),
+			pos.length_squared() - self.radius.powi(2),
-			dir: vec3(0.0, 1.0, 0.0)
+		)
-		}),
+		.expect("the ray starts inside and is not along the axis so *has* to cross the side");
-		Some(3.0)
+		Some(t_side_out)
-	);
+	}
-	assert_eq!(
+}
-		r.trace_out_of(Ray {
+
-			pos: vec3(0.0, 1.0, 0.0),
+#[cfg(test)]
-			dir: vec3(0.0, -1.0, 0.0)
+mod tests {
-		}),
+	use super::*;
-		Some(4.0)
+	use crate::types::ray;
-	);
+	use approx::assert_abs_diff_eq;
-	assert_eq!(
+	use glam::vec3;
-		r.trace_out_of(Ray {
+
-			pos: vec3(1.0, 1.0, 1.0),
+	#[test]
-			dir: vec3(0.0, -1.0, 0.0)
+	fn test_rect() {
-		}),
+		assert_eq!(
-		Some(4.0)
+			Rect::flip_ray(ray(vec3(2., 3., 2.), vec3(4., 5., 4.))),
-	);
+			ray(vec3(2., 3., 2.), vec3(4., 5., 4.)),
-	assert_eq!(
+		);
-		r.trace_out_of(Ray {
+		assert_eq!(
-			pos: vec3(2.0, 3.0, 2.0),
+			Rect::flip_ray(ray(vec3(2., 3., 2.), vec3(-4., 5., -4.))),
-			dir: vec3(1.0, 1.0, 1.0)
+			ray(vec3(-2., 3., -2.), vec3(4., 5., 4.)),
-		}),
+		);
-		Some(0.0)
+		assert_eq!(
-	);
+			Rect::flip_ray(ray(vec3(2., 3., 2.), vec3(4., -5., 4.))),
+			ray(vec3(2., -3., 2.), vec3(4., 5., 4.)),
+		);
+		assert_eq!(
+			Rect::flip_ray(ray(vec3(2., 3., 2.), vec3(4., 0., 4.))),
+			ray(vec3(2., 3., 2.), vec3(4., 0., 4.)),
+		);
+
+		let r = Rect { size: vec3(2., 3., 2.) };
+
+		assert_eq!(r.trace_into(ray(vec3(3., 3., 3.), vec3(1., 1., 1.))), None);
+		assert_eq!(r.trace_into(ray(vec3(-3., 2., -3.), vec3(1., 0., 1.))), Some(1.));
+		assert_eq!(r.trace_into(ray(vec3(-3., 2., -3.), vec3(-1., 0., -1.))), None);
+		assert_eq!(r.trace_into(ray(vec3(-3., 1., -3.), vec3(2., 2., 2.))), Some(0.5));
+		assert_eq!(r.trace_into(ray(vec3(-3., 2.1, -3.), vec3(2., 2., 2.))), None);
+
+		assert_eq!(r.trace_into(ray(vec3(2., 3., 2.), vec3(1., 1., 1.))), None);
+		assert_eq!(r.trace_into(ray(vec3(-2., 3., -2.), vec3(-1., 1., -1.))), None);
+		assert_eq!(r.trace_into(ray(vec3(2., 3., 2.), vec3(-1., -1., -1.))), Some(0.));
+		assert_eq!(r.trace_into(ray(vec3(2., -3., 2.), vec3(-1., 1., -1.))), Some(0.));
+
+		assert_eq!(r.trace_out_of(ray(vec3(0., 0., 0.), vec3(1., 1., 1.))), Some(2.));
+		assert_eq!(r.trace_out_of(ray(vec3(0., 0., 0.), vec3(0., 1., 0.))), Some(3.));
+		assert_eq!(r.trace_out_of(ray(vec3(0., 1., 0.), vec3(0., -1., 0.))), Some(4.));
+		assert_eq!(r.trace_out_of(ray(vec3(1., 1., 1.), vec3(0., -1., 0.))), Some(4.));
+		assert_eq!(r.trace_out_of(ray(vec3(2., 3., 2.), vec3(1., 1., 1.))), Some(0.));
+	}
+
+	#[test]
+	fn test_cylinder() {
+		assert_eq!(
+			YCylinder::flip_ray(ray(vec3(2., 3., 2.), vec3(4., 5., 4.))),
+			ray(vec3(2., 3., 2.), vec3(4., 5., 4.)),
+		);
+		assert_eq!(
+			YCylinder::flip_ray(ray(vec3(2., 3., 2.), vec3(-4., 5., -4.))),
+			ray(vec3(-2., 3., -2.), vec3(4., 5., 4.)),
+		);
+		assert_eq!(
+			YCylinder::flip_ray(ray(vec3(2., 3., 2.), vec3(4., -5., 4.))),
+			ray(vec3(2., -3., 2.), vec3(4., 5., 4.)),
+		);
+		assert_eq!(
+			YCylinder::flip_ray(ray(vec3(2., 3., 2.), vec3(4., 0., 4.))),
+			ray(vec3(2., 3., 2.), vec3(4., 0., 4.)),
+		);
+
+		let r = YCylinder {
+			half_length: 3.,
+			radius: 2.,
+		};
+
+		assert_eq!(r.trace_into(ray(vec3(3., 4., 3.), vec3(0., -1., 0.))), None);
+		assert_eq!(r.trace_into(ray(vec3(1., 4., 1.), vec3(0., -1., 0.))), Some(1.));
+		assert_eq!(r.trace_into(ray(vec3(3., 3., 3.), vec3(1., 1., 1.))), None);
+		assert_abs_diff_eq!(
+			r.trace_into(ray(vec3(-3., 2., -3.), vec3(1., 0., 1.))).unwrap(),
+			1.5857864
+		);
+		assert_eq!(r.trace_into(ray(vec3(-3., 2., -3.), vec3(-1., 0., -1.))), None);
+		assert_abs_diff_eq!(
+			r.trace_into(ray(vec3(-3., 1., -3.), vec3(2., 2., 2.))).unwrap(),
+			0.7928932
+		);
+		assert_eq!(r.trace_into(ray(vec3(-3., 2.1, -3.), vec3(2., 2., 2.))), None);
+
+		assert_eq!(r.trace_into(ray(vec3(2., 3., 2.), vec3(1., 1., 1.))), None);
+		assert_eq!(r.trace_into(ray(vec3(-2., 3., -2.), vec3(-1., 1., -1.))), None);
+		assert_eq!(
+			r.trace_into(ray(vec3(1.4142135, 3., 1.4142135), vec3(-1., -1., -1.))),
+			Some(0.)
+		);
+		assert_eq!(
+			r.trace_into(ray(vec3(1.4142135, -3., 1.4142135), vec3(-1., 1., -1.))),
+			Some(0.)
+		);
+
+		assert_abs_diff_eq!(
+			r.trace_out_of(ray(vec3(0., 0., 0.), vec3(1., 1., 1.))).unwrap(),
+			1.4142135
+		);
+		assert_eq!(r.trace_out_of(ray(vec3(0., 0., 0.), vec3(0., 1., 0.))), Some(3.));
+		assert_eq!(r.trace_out_of(ray(vec3(0., 1., 0.), vec3(0., -1., 0.))), Some(4.));
+		assert_eq!(r.trace_out_of(ray(vec3(1., 1., 1.), vec3(0., -1., 0.))), Some(4.));
+		assert_abs_diff_eq!(
+			r.trace_out_of(ray(vec3(1.4142135, 3., 1.4142135), vec3(1., 1., 1.)))
+				.unwrap(),
+			0.
+		);
+	}
 }

src/types.rs

@@ -6,6 +6,10 @@ pub struct Ray {
 	pub dir: Vec3,
 }
 
+pub fn ray(pos: Vec3, dir: Vec3) -> Ray {
+	Ray { pos, dir }
+}
+
 impl Ray {
 	pub fn forward(&self, dist: f32) -> Ray {
 		Ray {

src/utils.rs

@@ -43,16 +43,8 @@ mod tests {
 			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!(
+		assert_abs_diff_eq!(loc.rot * vec3(1., 0., 0.), vec3(1. / 3., 0., 0.), epsilon = ε);
-			loc.rot * vec3(1., 0., 0.),
+		assert_abs_diff_eq!(loc.rot * vec3(0., 1., 0.), vec3(0., 1. / 4., 0.), epsilon = ε);
-			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,
@@ -60,16 +52,8 @@ mod tests {
 			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!(
+		assert_abs_diff_eq!(loc.rot * vec3(1., 0., 0.), vec3(0., 1. / 4., 0.), epsilon = ε);
-			loc.rot * vec3(1., 0., 0.),
+		assert_abs_diff_eq!(loc.rot * vec3(0., 1., 0.), vec3(-1. / 3., 0., 0.), epsilon = ε);
-			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(
@@ -78,15 +62,7 @@ mod tests {
 			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!(
+		assert_abs_diff_eq!(loc.rot * vec3(1., 0., 0.), vec3(1. / 5., 1. / 5., 0.), epsilon = ε);
-			loc.rot * vec3(1., 0., 0.),
+		assert_abs_diff_eq!(loc.rot * vec3(0., 1., 0.), vec3(-4. / 15., 3. / 20., 0.), epsilon = ε);
-			vec3(1. / 5., 1. / 5., 0.),
-			epsilon = ε
-		);
-		assert_abs_diff_eq!(
-			loc.rot * vec3(0., 1., 0.),
-			vec3(-4. / 15., 3. / 20., 0.),
-			epsilon = ε
-		);
 	}
 }