Some concerns ain’t easily separated

This code works but boy it is ugly...

Cargo.toml

@@ -16,3 +16,6 @@ glam = "0.27.0"
 flo_draw = "0.3.1"
 flo_canvas = "0.3.1"
 itertools-num = "0.1.3"
+
+[dev-dependencies]
+approx = "0.5.1"

src/bin/flat.rs

@@ -1,7 +1,15 @@
+use std::rc::Rc;
 use flo_draw::*;
 use flo_canvas::*;
 use glam::*;
 use riemann::{Decomp2, Metric, trace_iter};
+use crate::boundary::Loop;
+
+#[cfg(test)]
+use approx::assert_abs_diff_eq;
+
+const RAYS_IN_FAN: usize = 101;
+const DT: f32 = 0.1;
 
 pub fn main() {
 	let space = Coil {
@@ -11,46 +19,149 @@ pub fn main() {
 		m: 10.0,
 	};
 	let space = Rect {
-		scale: 3.0,
+		inner_radius: 30.0,
-		r: vec2(30.0, 300.0),
+		outer_radius: 50.0,
-		m: vec2(20.0, 50.0),
+		internal_halflength: 100.0,
+		external_halflength: 300.0,
 	};
 	with_2d_graphics(move || {
 		let canvas = create_drawing_window("Refraction");
 		canvas.draw(|gc| {
 			gc.canvas_height(1000.0);
-			space.render(gc);
+			gc.line_width(0.5);
+
+			// gc.canvas_height(100.0);
+			// gc.center_region(-100.0, -200.0, 100.0, 0.0);
+			// gc.line_width(0.1);
+
+			space.draw_background(gc);
+
+			// gc.stroke_color(Color::Rgba(0.5, 1.0, 0.0, 1.0));
+			// space.draw_fan(gc, vec2(0.0, -0.5 * space.external_halflength), vec2(1.0, 1.0), 1.0);
+			space.draw_rays(gc, &space);
+
+			let space = Rc::new(space);
+			let parts: Vec<Box<dyn SpacePart>> = vec![
+				Box::new(Outside { space: space.clone() }),
+				Box::new(Wall { space: space.clone() }),
+				Box::new(Inside { space: space.clone() }),
+			];
+			space.draw_rays(gc, parts.as_slice());
+			// gc.stroke_color(Color::Rgba(1.0, 0.0, 1.0, 1.0));
+			// parts.draw_fan(gc, vec2(0.0, space.u(-0.5 * space.external_halflength)), vec2(1.0, 1.0), 1.0);
 		});
 	});
 }
 
-fn draw_ray(gc: &mut Vec<Draw>, space: &impl Metric, base: Vec2, dir: Vec2) {
+trait Interesting {
-	let dir = space.globalize(base, dir);
+	fn draw_background(&self, gc: &mut Vec<Draw>);
+	fn draw_rays(&self, gc: &mut Vec<Draw>, tracer: &(impl Rayable + ?Sized));
+}
+
+trait SpacePart: boundary::Boundary + Metric + SpaceVisual {}
+
+impl<T: boundary::Boundary + Metric + SpaceVisual> SpacePart for T {}
+
+trait SpaceVisual {
+	fn color(&self) -> Color;
+	fn globalize_loc(&self, pos: Vec2) -> Vec2;
+}
+
+impl SpaceVisual for Outside {
+	fn color(&self) -> Color {
+		Color::Rgba(0.7, 0.7, 0.7, 1.0)
+	}
+
+	fn globalize_loc(&self, pos: Vec2) -> Vec2 {
+		pos
+	}
+}
+
+impl SpaceVisual for Wall {
+	fn color(&self) -> Color {
+		Color::Rgba(1.0, 0.0, 0.0, 1.0)
+	}
+
+	fn globalize_loc(&self, pos: Vec2) -> Vec2 {
+		pos
+	}
+}
+
+impl SpaceVisual for Inside {
+	fn color(&self) -> Color {
+		Color::Rgba(0.0, 0.7, 0.0, 1.0)
+	}
+
+	fn globalize_loc(&self, pos: Vec2) -> Vec2 {
+		vec2(pos.x, self.space.x(pos.y))
+	}
+}
+
+impl Rayable for [Box<dyn SpacePart>] {
+	fn draw_ray(&self, gc: &mut Vec<Draw>, base: Vec2, dir: Vec2) {
 		gc.new_path();
-	gc.move_to(base.x, base.y);
+		let dt = DT;
-	for pt in trace_iter(space, base, dir, 1.0).take(10000) {
+		let mut s = boundary::Id(0);
+		let mut p = base;
+		let mut v = dir.normalize();
+		let part = &*self[s.0 as usize];
+		gc.stroke_color(part.color());
+		gc.move_to(p.x, p.y);
+		for _ in 0..10000 {
+			let part = &*self[s.0 as usize];
+			let a: Vec2 = -riemann::convolute(riemann::krist(part, p), v);
+			v = v + a * dt;
+			if let Some((id, base, dir)) = part.next(p, v, dt) {
+				gc.stroke();
+				gc.new_path();
+				let pt = part.globalize_loc(p);
+				gc.move_to(pt.x, pt.y);
+				s = id;
+				p = base;
+				v = dir;
+				let part = &*self[s.0 as usize];
+				let pt = part.globalize_loc(p);
+				gc.stroke_color(part.color());
+				gc.line_to(pt.x, pt.y);
+			} else {
+				p = p + v * dt;
+				let pt = part.globalize_loc(p);
 				gc.line_to(pt.x, pt.y);
-		if pt.abs().cmpgt(Vec2::splat(1000.0)).any() {
-			break
 			}
 		}
 		gc.stroke();
 	}
+}
 
-fn draw_fan(gc: &mut Vec<Draw>, space: &impl Metric, base: Vec2, dir: Vec2, spread: f32) {
+trait Rayable {
+	fn draw_ray(&self, gc: &mut Vec<Draw>, base: Vec2, dir: Vec2);
+
+	fn draw_fan(&self, gc: &mut Vec<Draw>, base: Vec2, dir: Vec2, spread: f32) {
 		let dir = dir.normalize();
 		let v = vec2(-dir.y, dir.x);
-	for y in itertools_num::linspace(-spread, spread, 101) {
+		for y in itertools_num::linspace(-spread, spread, RAYS_IN_FAN) {
-		draw_ray(gc, space, base, dir + y * v);
+			self.draw_ray(gc, base, dir + y * v);
+		}
 	}
 }
 
-trait Renderable {
+impl<T: Metric> Rayable for T {
-	fn render(&self, gc: &mut Vec<Draw>);
+	fn draw_ray(&self, gc: &mut Vec<Draw>, base: Vec2, dir: Vec2) {
+		let dir = self.globalize(base, dir);
+		gc.new_path();
+		gc.move_to(base.x, base.y);
+		for pt in trace_iter(self, base, dir, DT).take(10000) {
+			gc.line_to(pt.x, pt.y);
+			if pt.abs().cmpgt(Vec2::splat(1000.0)).any() {
+				break;
+			}
+		}
+		gc.stroke();
+	}
 }
 
-impl Renderable for Coil {
+impl Interesting for Coil {
-	fn render(&self, gc: &mut Vec<Draw>) {
+	fn draw_background(&self, gc: &mut Vec<Draw>) {
 		gc.new_path();
 		gc.circle(0.0, 0.0, self.r + self.w + self.m);
 		gc.circle(0.0, 0.0, self.r + self.w);
@@ -59,27 +170,35 @@ impl Renderable for Coil {
 		gc.winding_rule(WindingRule::EvenOdd);
 		gc.fill_color(Color::Rgba(0.8, 0.8, 0.8, 1.0));
 		gc.fill();
-		gc.line_width(0.5);
+	}
+
+	fn draw_rays(&self, gc: &mut Vec<Draw>, tracer: &(impl Rayable + ?Sized)) {
 		gc.stroke_color(Color::Rgba(1.0, 0.5, 0.0, 1.0));
-		draw_fan(gc, self, vec2(-500.0, 0.0), vec2(1.0, 0.0), 1.0);
+		tracer.draw_fan(gc, vec2(-500.0, 0.0), vec2(1.0, 0.0), 1.0);
 		gc.stroke_color(Color::Rgba(0.0, 0.5, 1.0, 1.0));
-		draw_fan(gc, self, vec2(0.0, self.r), vec2(1.0, 0.0), 1.0);
+		tracer.draw_fan(gc, vec2(0.0, self.r), vec2(1.0, 0.0), 1.0);
 	}
 }
 
-impl Renderable for Rect {
+impl Interesting for Rect {
-	fn render(&self, gc: &mut Vec<Draw>) {
+	fn draw_background(&self, gc: &mut Vec<Draw>) {
 		gc.new_path();
-		gc.rect(-self.r.x - self.m.x, -self.r.y - self.m.y, self.r.x + self.m.x, self.r.y + self.m.y);
+		gc.rect(-self.outer_radius, -self.external_halflength, self.outer_radius, self.external_halflength);
-		gc.rect(-self.r.x, -self.r.y, self.r.x, self.r.y);
+		gc.rect(-self.inner_radius, -self.external_halflength, self.inner_radius, self.external_halflength);
 		gc.winding_rule(WindingRule::EvenOdd);
 		gc.fill_color(Color::Rgba(0.8, 0.8, 0.8, 1.0));
 		gc.fill();
-		gc.line_width(0.5);
+	}
+
+	fn draw_rays(&self, gc: &mut Vec<Draw>, tracer: &(impl Rayable + ?Sized)) {
 		gc.stroke_color(Color::Rgba(1.0, 0.5, 0.0, 1.0));
-		draw_fan(gc, self, vec2(-500.0, 0.0), vec2(1.0, 0.0), 1.0);
+		tracer.draw_fan(gc, vec2(-500.0, 0.0), vec2(1.0, 0.0), 1.0);
+		gc.stroke_color(Color::Rgba(0.5, 1.0, 0.0, 1.0));
+		tracer.draw_fan(gc, vec2(0.0, -0.5 * self.external_halflength), vec2(1.0, 1.0), 1.0);
 		gc.stroke_color(Color::Rgba(0.0, 0.5, 1.0, 1.0));
-		draw_fan(gc, self, vec2(0.0, -0.5 * self.r.y), vec2(1.0, 1.0), 1.0);
+		tracer.draw_fan(gc, vec2(0.0, -0.5 * self.external_halflength), vec2(1.0, 1.0), 1.0);
+		gc.stroke_color(Color::Rgba(0.2, 0.7, 0.0, 1.0));
+		tracer.draw_fan(gc, vec2(-0.5 * self.inner_radius, -1.2 * self.external_halflength), vec2(0.0, 1.0), 1.0);
 	}
 }
 
@@ -109,21 +228,228 @@ impl Metric for Coil {
 }
 
 struct Rect {
-	scale: f32,
+	outer_radius: f32,
-	r: Vec2,
+	inner_radius: f32,
-	m: Vec2,
+	external_halflength: f32,
+	internal_halflength: f32,
+}
+
+impl Rect {
+	fn γ(&self) -> f32 { self.external_halflength / self.internal_halflength }
+	fn ri(&self) -> f32 { self.internal_halflength }
+	fn re(&self) -> f32 { self.external_halflength }
+	fn a(&self) -> f32 { (1.0 - self.γ()) / self.ri() }
+	fn b(&self) -> f32 { 2.0 * self.γ() - 1.0 }
+
+	fn root(&self, x: f32) -> f32 { ((2.0 * self.γ() - 1.0).powi(2) + 4.0 * (1.0 - self.γ()) * x / self.ri()).sqrt() }
+	fn d(&self, u: f32) -> f32 { 2.0 * self.a() * u + self.b() }
+	pub fn x(&self, u: f32) -> f32 { (self.a() * u.abs() + self.b()) * u }
+	pub fn u(&self, x: f32) -> f32 { 0.5 * self.ri() * (1.0 - 2.0 * self.γ() + self.root(x.abs())) / (1.0 - self.γ()) * x.signum() }
+	pub fn dx(&self, u: f32, du: f32) -> f32 { du * self.d(u.abs()) }
+	pub fn du(&self, x: f32, dx: f32) -> f32 { dx / self.root(x.abs()) }
+}
+
+#[test]
+fn test_rect() {
+	let r = Rect {
+		outer_radius: 50.0,
+		inner_radius: 20.0,
+		external_halflength: 100.0,
+		internal_halflength: 10.0,
+	};
+	assert_abs_diff_eq!(r.x(r.internal_halflength), r.external_halflength, epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.x(-r.internal_halflength), -r.external_halflength, epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.u(r.external_halflength), r.internal_halflength, epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.u(-r.external_halflength), -r.internal_halflength, epsilon = 1.0e-5);
+
+	assert_abs_diff_eq!(r.dx(r.internal_halflength, 3.0), 3.0, epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.dx(-r.internal_halflength, 3.0), 3.0, epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.du(r.external_halflength, 3.0), 3.0, epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.du(-r.external_halflength, 3.0), 3.0, epsilon = 1.0e-5);
+
+	assert_abs_diff_eq!(r.u(r.x(1.0)), 1.0, epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.u(r.x(5.0)), 5.0, epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.u(r.x(-5.0)), -5.0, epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.u(r.x(10.0)), 10.0, epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.x(r.u(10.0)), 10.0, epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.x(r.u(50.0)), 50.0, epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.x(r.u(-50.0)), -50.0, epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.x(r.u(100.0)), 100.0, epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.d(r.u(10.0)), r.root(10.0), epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.d(r.u(50.0)), r.root(50.0), epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.d(r.u(100.0)), r.root(100.0), epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.du(10.0, r.dx(r.u(10.0), 3.0)), 3.0, epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.du(50.0, r.dx(r.u(50.0), 3.0)), 3.0, epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.du(-50.0, r.dx(r.u(-50.0), 3.0)), 3.0, epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.du(100.0, r.dx(r.u(100.0), 3.0)), 3.0, epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.dx(1.0, r.du(r.x(1.0), 3.0)), 3.0, epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.dx(5.0, r.du(r.x(5.0), 3.0)), 3.0, epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.dx(-5.0, r.du(r.x(-5.0), 3.0)), 3.0, epsilon = 1.0e-5);
+	assert_abs_diff_eq!(r.dx(10.0, r.du(r.x(10.0), 3.0)), 3.0, epsilon = 1.0e-5);
 }
 
 impl Metric for Rect {
 	fn halfmetric(&self, pos: Vec2) -> Decomp2 {
-		let s = smoothbox(pos.x, vec2(-self.r.x, self.r.x), self.m.x) * smoothbox(pos.y, vec2(-self.r.y, self.r.y), self.m.y);
+		let x = pos.y.abs();
+		let sx = ((pos.x.abs() - self.inner_radius) / (self.outer_radius - self.inner_radius)).clamp(0.0, 1.0);
+		let sy = if x <= self.external_halflength { 1.0 / self.root(x) } else { 1.0 };
+		assert!(sx.is_finite());
+		assert!(sy.is_finite());
+		assert!(sy > 0.0);
 		Decomp2 {
 			ortho: Mat2::IDENTITY,
-			diag: vec2(1.0, self.scale.powf(-s)),
+			diag: vec2(1.0, sy.lerp(1.0, sx)),
 		}
 	}
 }
 
+struct Flat;
+
+impl Metric for Flat {
+	fn halfmetric(&self, pos: Vec2) -> Decomp2 {
+		Decomp2 {
+			ortho: Mat2::IDENTITY,
+			diag: Vec2::splat(1.0),
+		}
+	}
+}
+
+struct Outside {
+	space: Rc<Rect>,
+}
+
+struct Wall {
+	space: Rc<Rect>,
+}
+
+struct Inside {
+	space: Rc<Rect>,
+}
+
+impl boundary::Boundary for Outside {
+	fn next(&self, base: Vec2, dir: Vec2, limit: f32) -> Option<(boundary::Id, Vec2, Vec2)> {
+		let or = self.space.outer_radius;
+		let ir = self.space.inner_radius;
+		let hl = self.space.external_halflength;
+		let bnd = Loop(vec![
+			vec2(-or, -hl), vec2(-ir, -hl), vec2(ir, -hl), vec2(or, -hl),
+			vec2(or, hl), vec2(ir, hl), vec2(-ir, hl), vec2(-or, hl),
+		]);
+		let (side, dist) = bnd.hit(base, dir)?;
+		if dist <= limit {
+			let pt = base + dist * dir;
+			return match side {
+				1 | 5 => Some((boundary::Id(2), vec2(pt.x, self.space.u(pt.y)), vec2(dir.x, self.space.du(pt.y, dir.y)))),
+				_ => Some((boundary::Id(1), pt, dir)),
+			};
+		}
+		None
+	}
+}
+
+impl Metric for Outside {
+	fn halfmetric(&self, pos: Vec2) -> Decomp2 {
+		Flat {}.halfmetric(pos)
+	}
+}
+
+impl boundary::Boundary for Wall {
+	fn next(&self, base: Vec2, dir: Vec2, limit: f32) -> Option<(boundary::Id, Vec2, Vec2)> {
+		let or = self.space.outer_radius;
+		let ir = self.space.inner_radius;
+		let hl = self.space.external_halflength;
+		let obnd = Loop(vec![vec2(-or, -hl), vec2(-or, hl), vec2(or, hl), vec2(or, -hl)]);
+		let ibnd = Loop(vec![vec2(-ir, -hl), vec2(ir, -hl), vec2(ir, hl), vec2(-ir, hl)]);
+		if let Some((_, dist)) = ibnd.hit(base, dir) {
+			if dist <= limit {
+				let p = base + dist * dir;
+				let v = dir;
+				let v = vec2(v.x, self.space.du(p.y, v.y));
+				let p = vec2(p.x, self.space.u(p.y));
+				return Some((boundary::Id(2), p, v));
+			}
+		}
+		if let Some((_, dist)) = obnd.hit(base, dir) {
+			if dist <= limit {
+				return Some((boundary::Id(0), base + dist * dir, dir));
+			}
+		}
+		None
+	}
+}
+
+impl Metric for Wall {
+	fn halfmetric(&self, pos: Vec2) -> Decomp2 {
+		self.space.halfmetric(pos)
+	}
+}
+
+impl boundary::Boundary for Inside {
+	fn next(&self, base: Vec2, dir: Vec2, limit: f32) -> Option<(boundary::Id, Vec2, Vec2)> {
+		let size = vec2(self.space.inner_radius, self.space.internal_halflength);
+		let bnd = Loop(vec![vec2(-size.x, -size.y), vec2(-size.x, size.y), vec2(size.x, size.y), vec2(size.x, -size.y)]);
+		let (side, dist) = bnd.hit(base, dir)?;
+		if dist <= limit {
+			let p = base + dist * dir;
+			let v = dir;
+			let v = vec2(v.x, self.space.dx(p.y, v.y));
+			let p = vec2(p.x, self.space.x(p.y));
+			return match side {
+				0 | 2 => Some((boundary::Id(1), p, v)),
+				1 | 3 => Some((boundary::Id(0), p, v)),
+				_ => panic!()
+			};
+		}
+		None
+	}
+}
+
+impl Metric for Inside {
+	fn halfmetric(&self, pos: Vec2) -> Decomp2 {
+		Flat {}.halfmetric(pos)
+	}
+}
+
+mod boundary {
+	use glam::*;
+
+	pub struct Id(pub u8);
+
+	pub trait Boundary {
+		fn next(&self, base: Vec2, dir: Vec2, limit: f32) -> Option<(Id, Vec2, Vec2)>;
+	}
+
+	pub struct Loop(pub Vec<Vec2>);
+
+	impl Loop {
+		pub fn hit(&self, base: Vec2, dir: Vec2) -> Option<(usize, f32)> {
+			self.0.iter().enumerate().filter_map(|(k, &a)| {
+				let b = self.0[(k + 1) % self.0.len()];
+				let u = mat2(a - base, dir).determinant();
+				let v = mat2(b - base, dir).determinant();
+				if u < 0.0 && v > 0.0 {
+					let dist = mat2(a - base, b - a).determinant() / mat2(dir, b - a).determinant();
+					if dist >= 0.0 { Some((k, dist)) } else { None }
+				} else {
+					None
+				}
+			}).min_by(|(k1, dist1), (k2, dist2)| dist1.total_cmp(dist2))
+		}
+	}
+
+	#[test]
+	fn test_loop() {
+		let tri = Loop(vec![vec2(-1.0, -1.0), vec2(1.0, -1.0), vec2(0.0, 1.0)]);
+		assert_eq!(tri.hit(vec2(0.0, -2.0), vec2(0.0, 1.0)), Some((0, 1.0)));
+		assert_eq!(tri.hit(vec2(0.0, -2.0), vec2(0.0, 0.5)), Some((0, 2.0)));
+		assert_eq!(tri.hit(vec2(0.0, -2.0), vec2(1.0, 1.0)), None);
+		assert_eq!(tri.hit(vec2(0.0, 0.0), vec2(0.0, 1.0)), None);
+		assert_eq!(tri.hit(vec2(0.0, 0.0), vec2(0.0, -1.0)), None);
+		assert_eq!(tri.hit(vec2(-1.5, 0.5), vec2(2.0, -1.0)), Some((2, 0.5)));
+		assert_eq!(tri.hit(vec2(-1.5, 0.5), vec2(-2.0, 1.0)), None);
+	}
+}
+
 mod riemann {
 	use glam::*;
 
@@ -211,9 +537,9 @@ mod riemann {
 		}
 	}
 
-	fn krist(space: &impl Metric, pos: Vec2) -> Tens2 {
+	pub fn krist(space: &(impl Metric + ?Sized), 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.invmetric(pos); // с верхними индексами
+		let g = space.invmetric(pos); // с верхними индексами
 		let d = space.dmetric(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>())
@@ -230,10 +556,10 @@ mod riemann {
 		]
 	}
 
-	fn convolute(t: Tens2, v: Vec2) -> Vec2 {
+	pub fn convolute(t: Tens2, v: Vec2) -> Vec2 {
 		vec2(
 			v.dot(t[0] * v),
-			v.dot(t[1] * v)
+			v.dot(t[1] * v),
 		)
 	}
 
@@ -273,12 +599,18 @@ mod riemann {
 }
 
 fn smoothstep(x: f32) -> f32 {
-	return 3.0 * x * x - 2.0 * x * x * x;
+	3.0 * x * x - 2.0 * x * x * x
 }
 
-fn smoothbox(val: f32, range: Vec2, pad: f32) -> f32 {
+/// 1.0 for val∈[range.x, range.y], 0.0 for val∉[range.x−pad, range.y+pad], linear in-between.
+fn trapezoid(val: f32, range: Vec2, pad: f32) -> f32 {
 	let slope1 = 1.0 + (val - range.x) / pad;
 	let slope2 = 1.0 - (val - range.y) / pad;
 	let lin = slope1.min(slope2);
-	smoothstep(lin.clamp(0.0, 1.0))
+	lin.clamp(0.0, 1.0)
+}
+
+/// 1.0 for val∈[range.x, range.y], 0.0 for val∉[range.x−pad, range.y+pad], smoothstep in-between.
+fn smoothbox(val: f32, range: Vec2, pad: f32) -> f32 {
+	smoothstep(trapezoid(val, range, pad))
 }