diff --git a/src/bin/flat/main.rs b/src/bin/flat/main.rs index 7981c61..a5baffc 100644 --- a/src/bin/flat/main.rs +++ b/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, mut pts: impl Iterator) { +fn draw_loop(gc: &mut Vec, mut pts: impl Iterator) { 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), - Mat2::from_angle(y), + vec3(0.0, y * tube.external_halflength, 0.0), + Mat3::from_mat2(Mat2::from_angle(y)), ), r: 20.0, }) .collect(); let space = Space { tube, objs }; - let cam1 = put_object(&space.tube, vec2(-500., 0.), Mat2::IDENTITY); + let cam1 = put_object(&space.tube, vec3(-500., 0., 0.), Mat3::IDENTITY); let cam2 = put_object( &space.tube, - vec2(-2.5 * tube.outer_radius, 1.25 * tube.external_halflength), - mat2(vec2(1., -1.), vec2(1., 1.)), + vec3( + -2.5 * tube.outer_radius, + 1.25 * tube.external_halflength, + 0., + ), + mat3(vec3(1., -1., 0.), vec3(1., 1., 0.), vec3(0., 0., 1.)), ); let cam3 = put_object( &space.tube, - vec2(0.25 * tube.inner_radius, 0.25 * tube.external_halflength), - mat2(vec2(0., -1.), vec2(1., 0.)), + vec3( + 0.25 * tube.inner_radius, + 0.25 * tube.external_halflength, + 0., + ), + mat3(vec3(0., -1., 0.), vec3(1., 0., 0.), vec3(0., 0., 1.)), ); gc.canvas_height(500.0); @@ -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.))); - assert_eq!(loc.pos, vec2(1., 2.)); - assert_abs_diff_eq!(loc.rot * vec2(1., 0.), vec2(1. / 3., 0.), epsilon = ε); - assert_abs_diff_eq!(loc.rot * vec2(0., 1.), vec2(0., 1. / 4.), epsilon = ε); + let loc = put_object( + &m, + vec3(1., 2., 0.), + mat3(vec3(1., 0., 0.), vec3(0., 1., 0.), vec3(0., 0., 1.)), + ); + assert_eq!(loc.pos, vec3(1., 2., 0.)); + assert_abs_diff_eq!( + loc.rot * vec3(1., 0., 0.), + vec3(1. / 3., 0., 0.), + epsilon = ε + ); + assert_abs_diff_eq!( + loc.rot * vec3(0., 1., 0.), + vec3(0., 1. / 4., 0.), + epsilon = ε + ); - let loc = put_object(&m, vec2(1., 2.), mat2(vec2(0., 1.), vec2(-1., 0.))); - assert_eq!(loc.pos, vec2(1., 2.)); - assert_abs_diff_eq!(loc.rot * vec2(1., 0.), vec2(0., 1. / 4.), epsilon = ε); - assert_abs_diff_eq!(loc.rot * vec2(0., 1.), vec2(-1. / 3., 0.), epsilon = ε); + let loc = put_object( + &m, + vec3(1., 2., 0.), + mat3(vec3(0., 1., 0.), vec3(-1., 0., 0.), vec3(0., 0., 1.)), + ); + assert_eq!(loc.pos, vec3(1., 2., 0.)); + assert_abs_diff_eq!( + loc.rot * vec3(1., 0., 0.), + vec3(0., 1. / 4., 0.), + epsilon = ε + ); + assert_abs_diff_eq!( + loc.rot * vec3(0., 1., 0.), + vec3(-1. / 3., 0., 0.), + epsilon = ε + ); let c = 0.5 * std::f32::consts::SQRT_2; - let loc = put_object(&m, vec2(1., 2.), mat2(vec2(c, c), vec2(-c, c))); - assert_eq!(loc.pos, vec2(1., 2.)); - assert_abs_diff_eq!(loc.rot * vec2(1., 0.), vec2(1. / 5., 1. / 5.), epsilon = ε); + let loc = put_object( + &m, + vec3(1., 2., 0.), + mat3(vec3(c, c, 0.), vec3(-c, c, 0.), vec3(0., 0., 1.)), + ); + assert_eq!(loc.pos, vec3(1., 2., 0.)); assert_abs_diff_eq!( - loc.rot * vec2(0., 1.), - vec2(-4. / 15., 3. / 20.), + loc.rot * vec3(1., 0., 0.), + vec3(1. / 5., 1. / 5., 0.), + epsilon = ε + ); + assert_abs_diff_eq!( + loc.rot * vec3(0., 1., 0.), + vec3(-4. / 15., 3. / 20., 0.), epsilon = ε ); } @@ -188,8 +231,8 @@ fn draw_cross(gc: &mut Vec, pos: Vec2, r: f32) { gc.line_to(pos.x + r, pos.y - r); } -fn draw_ray_2(gc: &mut Vec, space: &Space, camera: Location, dir: Vec2) { - let pos = vec2(0., 0.); +fn draw_ray_2(gc: &mut Vec, space: &Space, camera: Location, dir: Vec3) { + let pos = vec3(0., 0., 0.); let (hits, path) = space.trace_dbg(camera, Ray { pos, dir }); let hits2 = space.trace(camera, Ray { pos, dir }); for (a, b) in hits.into_iter().zip(hits2.into_iter()) { @@ -214,7 +257,7 @@ fn draw_ray_2(gc: &mut Vec, space: &Space, camera: Location, dir: Vec2) { fn draw_fan_2(gc: &mut Vec, 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, space: &Space, start: Vec2, dir: Vec2) { // let dir = space.tube.globalize(start, dir); // let v = space.tube.normalize(start, dir); let mut loc = Location { - pos: start, - rot: mat2(dir, vec2(-dir.y, dir.x)), + pos: vec3(start.x, start.y, 0.), + rot: mat3( + vec3(dir.x, dir.y, 0.), + vec3(-dir.y, dir.x, 0.), + vec3(0., 0., 1.), + ), }; - let v = vec2(1.0, 0.0); + let v = vec3(1., 0., 0.); let mut draw = |loc: &Location| { let p = loc.pos; let ax = p + loc.rot.x_axis * SCALE; diff --git a/src/ifaces.rs b/src/ifaces.rs index d28be75..d8d3e4c 100644 --- a/src/ifaces.rs +++ b/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, - pub end_dir: Vec2, + pub points: Vec, + pub end_dir: Vec3, } pub trait DebugTraceable: Traceable { diff --git a/src/riemann.rs b/src/riemann.rs index 961bfa4..d2b90a7 100644 --- a/src/riemann.rs +++ b/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, - v: Vec2, + p: Vec3, + v: Vec3, } impl<'a, M: Metric> Iterator for TraceIter<'a, M> { - type Item = Vec2; + type Item = Vec3; fn next(&mut self) -> Option { - 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(space: &M, base: Vec2, dir: Vec2, dt: f32) -> TraceIter { +pub fn trace_iter(space: &M, base: Vec3, dir: Vec3, dt: f32) -> TraceIter { TraceIter { space, p: base, @@ -52,54 +52,55 @@ pub fn trace_iter(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::() }) } -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, vec2(0.0, eps)), + dir_deriv(&f, pos, vec3(eps, 0.0, 0.0)), + dir_deriv(&f, pos, vec3(0.0, eps, 0.0)), + dir_deriv(&f, pos, vec3(0.0, 0.0, eps)), ] } /// Сворачивает тензор t с вектором u -pub fn contract(t: Tens2, u: Vec2) -> Mat2 { - mat2(t[0] * u, t[1] * u).transpose() +pub fn contract(t: Tens3, u: Vec3) -> Mat3 { + mat3(t[0] * u, t[1] * u, t[2] * u).transpose() } /// Сворачивает тензор t с вектором v дважды, по второму и третьему индексам. -pub fn contract2(t: 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 { - Vec2::from_array(std::array::from_fn(|i| f(i))) +fn make_vec3(f: impl Fn(usize) -> f32) -> Vec3 { + Vec3::from_array(std::array::from_fn(|i| f(i))) } -fn make_mat2(f: impl Fn(usize, usize) -> f32) -> Mat2 { - Mat2::from_cols_array_2d(&std::array::from_fn(|i| std::array::from_fn(|j| f(i, j)))) +fn make_mat3(f: impl Fn(usize, usize) -> f32) -> Mat3 { + Mat3::from_cols_array_2d(&std::array::from_fn(|i| std::array::from_fn(|j| f(i, j)))) } -fn make_tens2(f: impl Fn(usize, usize, usize) -> f32) -> Tens2 { - std::array::from_fn(|i| make_mat2(|j, k| f(i, j, k))) +fn make_tens3(f: impl Fn(usize, usize, usize) -> f32) -> Tens3 { + std::array::from_fn(|i| make_mat3(|j, k| f(i, j, k))) } #[test] -fn m2() { - let m = make_mat2(|i, j| (i + 2 * j) as f32); +fn m3() { + let m = make_mat3(|i, j| (i + 2 * j) as f32); assert_eq!(m.col(0)[0], 0.0); assert_eq!(m.col(1)[0], 1.0); assert_eq!(m.col(0)[1], 2.0); @@ -107,8 +108,8 @@ fn m2() { } #[test] -fn t2() { - let t = make_tens2(|i, j, k| (i + 2 * j + 4 * k) as f32); +fn t3() { + let t = make_tens3(|i, j, k| (i + 2 * j + 4 * k) as f32); assert_eq!(t[0].col(0)[0], 0.0); assert_eq!(t[1].col(0)[0], 1.0); assert_eq!(t[0].col(1)[0], 2.0); @@ -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 { - Decomp2 { + fn sqrt_at(&self, _pos: Vec3) -> Decomp3 { + 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 { - ortho: Mat2::IDENTITY, - diag: vec2(3., 4.) + Decomp3 { + ortho: Mat3::IDENTITY, + 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.]) - ); - 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.) + Mat3::from_cols_array(&[1. / 9., 0., 0., 0., 1. / 16., 0., 0., 0., 1. / 25.]) ); assert_eq!( - metric.normalize_vec_at(rng.gen(), vec2(0., 1.)), - vec2(0., 1. / 4.) + metric.part_derivs_at(rng.gen()), + [Mat3::ZERO, Mat3::ZERO, Mat3::ZERO] + ); + assert_eq!(metric.vec_length_at(rng.gen(), vec3(1., 0., 0.)), 3.); + assert_eq!(metric.vec_length_at(rng.gen(), vec3(0., 1., 0.)), 4.); + assert_eq!(metric.vec_length_at(rng.gen(), vec3(0., 0., 1.)), 5.); + assert_eq!(metric.vec_length_at(rng.gen(), vec3(1., 1., 0.)), 5.); + assert_eq!( + metric.normalize_vec_at(rng.gen(), vec3(1., 0., 0.)), + vec3(1. / 3., 0., 0.) ); assert_eq!( - metric.normalize_vec_at(rng.gen(), vec2(1., 1.)), - vec2(1. / 5., 1. / 5.) + metric.normalize_vec_at(rng.gen(), vec3(0., 1., 0.)), + vec3(0., 1. / 4., 0.) + ); + assert_eq!( + metric.normalize_vec_at(rng.gen(), vec3(1., 1., 0.)), + vec3(1. / 5., 1. / 5., 0.) ); } #[test] fn test_trace_iter() { let metric = samples::ScaledMetric { - scale: 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.), - vec2(0.5, 0.25), + vec3(3., 5., 0.), + vec3(0.5, 0.25, 0.), std::f32::consts::SQRT_2 ) .nth(7) .unwrap(), - vec2(7., 7.), + vec3(7., 7., 0.), epsilon = 1e-5 ); } diff --git a/src/tube/coords.rs b/src/tube/coords.rs index edcc0f1..a410fc9 100644 --- a/src/tube/coords.rs +++ b/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 { } pub trait FlatRegion: - FlatCoordinateSystem + FlatCoordinateSystem + FlatCoordinateSystem + FlatCoordinateSystem + FlatCoordinateSystem + FlatCoordinateSystem { // Измеряет расстояние до выхода за пределы области вдоль луча ray. Луч задаётся в плоской СК. fn distance_to_boundary(&self, _ray: Ray) -> Option { @@ -19,20 +19,20 @@ pub trait FlatRegion: } } -trait MetricCS: FlatCoordinateSystem { +trait MetricCS: FlatCoordinateSystem { 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 + MetricCS> FlatCoordinateSystem for T { +impl + MetricCS> FlatCoordinateSystem for T { fn flat_to_global(&self, ray: Ray) -> Ray { Ray { pos: self.flat_to_global(ray.pos), @@ -48,7 +48,7 @@ impl + MetricCS> FlatCoordinateSystem for T { } } -impl + MetricCS> FlatCoordinateSystem for T { +impl + MetricCS> FlatCoordinateSystem 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 for InnerCS { - fn flat_to_global(&self, pos: Vec2) -> Vec2 { - vec2(pos.x, self.0.y(pos.y)) +impl FlatCoordinateSystem for InnerCS { + fn flat_to_global(&self, pos: Vec3) -> Vec3 { + vec3(pos.x, self.0.y(pos.y), pos.z) } // Работает только при |pos.x| ≤ inner_radius или |pos.y| ≥ external_halflength. - fn global_to_flat(&self, pos: Vec2) -> Vec2 { - vec2(pos.x, self.0.v(pos.y)) + fn global_to_flat(&self, pos: Vec3) -> Vec3 { + vec3(pos.x, self.0.v(pos.y), pos.z) } } impl FlatRegion for InnerCS { fn distance_to_boundary(&self, ray: Ray) -> Option { 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 for OuterCS { - fn flat_to_global(&self, pos: Vec2) -> Vec2 { +impl FlatCoordinateSystem for OuterCS { + 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 for OuterCS { impl FlatRegion for OuterCS { fn distance_to_boundary(&self, ray: Ray) -> Option { 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); - impl FlatCoordinateSystem for Scaled { - fn flat_to_global(&self, pos: Vec2) -> Vec2 { + struct Scaled(samples::ScaledMetric, Vec3); + impl FlatCoordinateSystem for Scaled { + 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.flat_to_global(vec2(8., 7.)), vec2(2., 1.)); + assert_eq!(cs.global_to_flat(vec3(7., 3., 1.)), vec3(23., 15., 12.)); + assert_eq!(cs.flat_to_global(vec3(8., 7., 17.)), vec3(2., 1., 2.)); assert_eq!( cs.global_to_flat(Ray { - pos: vec2(7., 3.), - dir: vec2(3., 2.) + pos: vec3(7., 3., 0.), + dir: vec3(3., 2., 0.) }), Ray { - pos: vec2(23., 15.), - dir: vec2(9., 8.) + pos: vec3(23., 15., 7.), + dir: vec3(9., 8., 0.) } ); assert_eq!( cs.flat_to_global(Ray { - pos: vec2(23., 15.), - dir: vec2(9., 8.) + pos: vec3(23., 15., 7.), + dir: vec3(9., 8., 0.) }), Ray { - pos: vec2(7., 3.), - dir: vec2(3., 2.) + pos: vec3(7., 3., 0.), + dir: vec3(3., 2., 0.) } ); assert_eq!( cs.global_to_flat(Location { - pos: vec2(2., 1.), - rot: mat2(vec2(0., 1.), vec2(-1., 0.)) + pos: vec3(2., 1., 0.), + rot: mat3(vec3(0., 1., 0.), vec3(-1., 0., 0.), vec3(0., 0., 1.)) }), Location { - pos: vec2(8., 7.), - rot: mat2(vec2(0., 4.), vec2(-3., 0.)) + pos: vec3(8., 7., 7.), + rot: mat3(vec3(0., 4., 0.), vec3(-3., 0., 0.), vec3(0., 0., 5.)) } ); assert_eq!( cs.flat_to_global(Location { - pos: vec2(2., 1.), - rot: mat2(vec2(0., 1.), vec2(-1., 0.)) + pos: vec3(2., 1., 7.), + rot: mat3(vec3(0., 1., 0.), vec3(-1., 0., 0.), vec3(0., 0., 1.)) }), Location { - pos: vec2(0., -0.5), - rot: mat2(vec2(0., 0.25), vec2(-1. / 3., 0.)) + pos: vec3(0., -0.5, 0.), + rot: mat3( + vec3(0., 0.25, 0.), + vec3(-1. / 3., 0., 0.), + vec3(0., 0., 0.2) + ) } ); } fn test_flat_region( region: &impl FlatRegion, - range_global: (Vec2, Vec2), - range_flat: (Vec2, Vec2), + range_global: (Vec3, Vec3), + 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, - range_a: (Vec2, Vec2), + a: Vec3, + range_a: (Vec3, Vec3), name_b: &str, - b: Vec2, - range_b: (Vec2, Vec2), + 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: @@ -261,72 +281,76 @@ 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); - let pos_flat = region.global_to_flat(pos_global); - check_range( - "global", - pos_global, - range_global, - "flat", - pos_flat, - range_flat, - ); - assert_eq_at!( - pos_global, - region - .global_to_flat(Location { - pos: pos_global, - rot: Mat2::IDENTITY - }) - .pos, - pos_flat - ); - assert_eq_at!(pos_global, region.flat_to_global(pos_flat), pos_global); - assert_eq_at!( - pos_global, - region - .flat_to_global(region.global_to_flat(Location { - pos: pos_global, - rot: Mat2::IDENTITY - })) - .rot, - Mat2::IDENTITY - ); + for z in linspace(range_global.0.z, range_global.1.z, 20) { + let pos_global = vec3(x, y, z); + let pos_flat = region.global_to_flat(pos_global); + check_range( + "global", + pos_global, + range_global, + "flat", + pos_flat, + range_flat, + ); + assert_eq_at!( + pos_global, + region + .global_to_flat(Location { + pos: pos_global, + rot: Mat3::IDENTITY + }) + .pos, + pos_flat + ); + assert_eq_at!(pos_global, region.flat_to_global(pos_flat), pos_global); + assert_eq_at!( + pos_global, + region + .flat_to_global(region.global_to_flat(Location { + pos: pos_global, + rot: Mat3::IDENTITY + })) + .rot, + Mat3::IDENTITY + ); + } } } for x in linspace(range_flat.0.x, range_flat.1.x, 20) { for y in linspace(range_flat.0.y, range_flat.1.y, 20) { - let pos_flat = vec2(x, y); - let pos_global = region.flat_to_global(pos_flat); - check_range( - "flat", - pos_flat, - range_flat, - "global", - pos_global, - range_global, - ); - assert_eq_at!( - pos_flat, - region - .flat_to_global(Location { - pos: pos_flat, - rot: Mat2::IDENTITY - }) - .pos, - pos_global - ); - assert_eq_at!(pos_flat, region.global_to_flat(pos_global), pos_flat); - assert_eq_at!( - pos_flat, - region - .global_to_flat(region.flat_to_global(Location { - pos: pos_global, - rot: Mat2::IDENTITY - })) - .rot, - Mat2::IDENTITY - ); + for z in linspace(range_flat.0.z, range_flat.1.z, 20) { + let pos_flat = vec3(x, y, z); + let pos_global = region.flat_to_global(pos_flat); + check_range( + "flat", + pos_flat, + range_flat, + "global", + pos_global, + range_global, + ); + assert_eq_at!( + pos_flat, + region + .flat_to_global(Location { + pos: pos_flat, + rot: Mat3::IDENTITY + }) + .pos, + pos_global + ); + assert_eq_at!(pos_flat, region.global_to_flat(pos_global), pos_flat); + assert_eq_at!( + pos_flat, + region + .global_to_flat(region.flat_to_global(Location { + pos: pos_global, + rot: Mat3::IDENTITY + })) + .rot, + Mat3::IDENTITY + ); + } } } } @@ -341,18 +365,18 @@ mod tests { }); test_flat_region( &mapper, - (vec2(-30.0, -300.0), vec2(30.0, 300.0)), - (vec2(-30.0, -100.0), vec2(30.0, 100.0)), + (vec3(-30.0, -300.0, -30.0), vec3(30.0, 300.0, 30.0)), + (vec3(-30.0, -100.0, -30.0), vec3(30.0, 100.0, 30.0)), ); test_flat_region( &mapper, - (vec2(-60.0, -400.0), vec2(60.0, -300.0)), - (vec2(-60.0, -200.0), vec2(60.0, -100.0)), + (vec3(-60.0, -400.0, -60.0), vec3(60.0, -300.0, 60.0)), + (vec3(-60.0, -200.0, -60.0), vec3(60.0, -100.0, 60.0)), ); test_flat_region( &mapper, - (vec2(-60.0, 300.0), vec2(60.0, 400.0)), - (vec2(-60.0, 100.0), vec2(60.0, 200.0)), + (vec3(-60.0, 300.0, -60.0), vec3(60.0, 400.0, 60.0)), + (vec3(-60.0, 100.0, -60.0), vec3(60.0, 200.0, 60.0)), ); } @@ -367,112 +391,122 @@ mod tests { // TODO replace 200.20016 with something sane test_flat_region( &mapper, - (vec2(-30.0, -300.0), vec2(30.0, -1.0)), - (vec2(-30.0, -300.0), vec2(30.0, -200.20016)), + (vec3(-30.0, -300.0, -30.0), vec3(30.0, -1.0, 30.0)), + (vec3(-30.0, -300.0, -30.0), vec3(30.0, -200.20016, 30.0)), ); test_flat_region( &mapper, - (vec2(-30.0, 1.0), vec2(30.0, 300.0)), - (vec2(-30.0, 200.20016), vec2(30.0, 300.0)), + (vec3(-30.0, 1.0, -30.0), vec3(30.0, 300.0, 30.0)), + (vec3(-30.0, 200.20016, -30.0), vec3(30.0, 300.0, 30.0)), ); test_flat_region( &mapper, - (vec2(-60.0, -400.0), vec2(60.0, -300.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)), + (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)), - (vec2(-60.0, 300.0), vec2(60.0, 400.0)), + (vec3(-60.0, 300.0, -60.0), vec3(60.0, 400.0, 60.0)), + (vec3(-60.0, 300.0, -60.0), vec3(60.0, 400.0, 60.0)), ); // straight for x in linspace(-60., 60., 20) { for y in linspace(-320., 320., 20) { - assert_eq!( - mapper - .global_to_flat(Location { - pos: vec2(x, y), - rot: Mat2::IDENTITY - }) - .pos - .x, - x - ); + for z in linspace(-60., 60., 20) { + assert_eq!( + mapper + .global_to_flat(Location { + pos: vec3(x, y, z), + rot: Mat3::IDENTITY + }) + .pos + .x, + x + ); + } } } // symmetrical for x in linspace(0., 60., 20) { for y in linspace(0., 320., 20) { - let pp = mapper - .global_to_flat(Location { - pos: vec2(x, y), - rot: Mat2::IDENTITY, - }) - .pos; - let np = mapper - .global_to_flat(Location { - pos: vec2(-x, y), - rot: Mat2::IDENTITY, - }) - .pos; - let pn = mapper - .global_to_flat(Location { - pos: vec2(x, -y), - rot: Mat2::IDENTITY, - }) - .pos; - let nn = mapper - .global_to_flat(Location { - pos: vec2(-x, -y), - rot: Mat2::IDENTITY, - }) - .pos; - assert_eq!(np, vec2(-pp.x, pp.y)); - assert_eq!(pn, vec2(pp.x, -pp.y)); - assert_eq!(nn, vec2(-pp.x, -pp.y)); + for z in linspace(0., 60., 20) { + let pp = mapper + .global_to_flat(Location { + pos: vec3(x, y, z), + rot: Mat3::IDENTITY, + }) + .pos; + let np = mapper + .global_to_flat(Location { + pos: vec3(-x, y, z), + rot: Mat3::IDENTITY, + }) + .pos; + let pn = mapper + .global_to_flat(Location { + pos: vec3(x, -y, z), + rot: Mat3::IDENTITY, + }) + .pos; + let nn = mapper + .global_to_flat(Location { + pos: vec3(-x, -y, z), + rot: Mat3::IDENTITY, + }) + .pos; + assert_eq!(np, vec3(-pp.x, pp.y, pp.z)); + assert_eq!(pn, vec3(pp.x, -pp.y, pp.z)); + assert_eq!(nn, vec3(-pp.x, -pp.y, pp.z)); + } } } // clean boundary for x in linspace(50., 60., 20) { for y in linspace(0., 320., 20) { - assert_eq!( - mapper - .global_to_flat(Location { - pos: vec2(x, y), - rot: Mat2::IDENTITY - }) - .pos - .y, - y - ); + for z in linspace(50., 60., 20) { + assert_eq!( + mapper + .global_to_flat(Location { + pos: vec3(x, y, z), + rot: Mat3::IDENTITY + }) + .pos + .y, + y + ); + } } } for x in linspace(0., 60., 20) { for y in linspace(300., 320., 20) { - assert_eq!( - mapper - .global_to_flat(Location { - pos: vec2(x, y), - rot: Mat2::IDENTITY - }) - .pos - .y, - y - ); + for z in linspace(0., 60., 20) { + assert_eq!( + mapper + .global_to_flat(Location { + pos: vec3(x, y, z), + rot: Mat3::IDENTITY + }) + .pos + .y, + y + ); + } } } // accelerating for x in linspace(-29., 29., 20) { for y in linspace(1., 299., 20) { - let v = mapper - .global_to_flat(Location { - pos: vec2(x, y), - rot: Mat2::IDENTITY, - }) - .pos - .y; - assert!(v > 200.0); - assert!(v > y); + for z in linspace(-29., 29., 20) { + let v = mapper + .global_to_flat(Location { + pos: vec3(x, y, z), + rot: Mat3::IDENTITY, + }) + .pos + .y; + assert!(v > 200.0); + assert!(v > y); + } } } } diff --git a/src/tube/metric.rs b/src/tube/metric.rs index 7f923ab..2b463d0 100644 --- a/src/tube/metric.rs +++ b/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::riemann::{Metric, Tens2}; +use crate::mathx::Decomp3; +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 { - ortho: Mat2::IDENTITY, - diag: vec2(1.0, s), + Decomp3 { + ortho: Mat3::IDENTITY, + diag: vec3(1.0, s, 1.0), } } - fn part_derivs_at(&self, pos: Vec2) -> Tens2 { + fn part_derivs_at(&self, pos: Vec3) -> Tens3 { let sx = self.fx().value(pos.x); let sy = self.fy().du(pos.y); let s = sx + sy - sx * sy; @@ -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]), - Mat2::from_cols_array(&[0.0, 0.0, 0.0, ds2_dy]), + Mat3::from_cols_array(&[0., 0., 0., 0., ds2_dx, 0., 0., 0., 0.]), + Mat3::from_cols_array(&[0., 0., 0., 0., ds2_dy, 0., 0., 0., 0.]), + Mat3::from_cols_array(&[0., 0., 0., 0., 0., 0., 0., 0., 0.]), ] } } @@ -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,18 +101,25 @@ mod test { let epsilon = 1.0e-3; let margin = 1.0 / 16.0; let mul = 1.0 + margin; - for x in itertools_num::linspace(-mul * testee.outer_radius, mul * testee.outer_radius, 100) + for x in itertools_num::linspace(-mul * testee.outer_radius, mul * testee.outer_radius, 20) { for y in itertools_num::linspace( -mul * testee.external_halflength, mul * testee.external_halflength, - 100, + 20, ) { - let pos = vec2(x, y); - let computed = testee.part_derivs_at(pos); - let reference = approx.part_derivs_at(pos); - let eq = (0..2).all(|coord| computed[coord].abs_diff_eq(reference[coord], epsilon)); - assert!(eq, "Bad derivative computation at {pos}:\n explicit: {computed:?}\n numerical: {reference:?}\n"); + for z in itertools_num::linspace( + -mul * testee.outer_radius, + mul * testee.outer_radius, + 20, + ) { + let pos = vec3(x, y, z); + let computed = testee.part_derivs_at(pos); + let reference = approx.part_derivs_at(pos); + let eq = + (0..2).all(|coord| computed[coord].abs_diff_eq(reference[coord], epsilon)); + assert!(eq, "Bad derivative computation at {pos}:\n explicit: {computed:?}\n numerical: {reference:?}\n"); + } } } } @@ -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 b = vec2(bx, space.tube.external_halflength + off); - let Δ = vec2(bx - ax, 2.0 * (space.tube.internal_halflength + off)); + let a = vec3(ax, -(space.tube.external_halflength + off), 0.); + let b = vec3(bx, space.tube.external_halflength + off, 0.); + let Δ = vec3(bx - ax, 2.0 * (space.tube.internal_halflength + off), 0.); let dir = Δ / (steps as f32); let traced = space.trace_iter(Ray { pos: a, dir }).nth(steps).unwrap(); assert_abs_diff_eq!(traced.pos.x, b.x, epsilon = 1.0e-2); @@ -170,9 +178,9 @@ mod test { space.tube.inner_radius - ε, 20, ) { - let a = vec2(ax, -(space.tube.external_halflength + off)); - let b = vec2(bx, space.tube.external_halflength + off); - let Δ = vec2(bx - ax, 2.0 * (space.tube.internal_halflength + off)); + let a = vec3(ax, -(space.tube.external_halflength + off), 0.); + let b = vec3(bx, space.tube.external_halflength + off, 0.); + let Δ = vec3(bx - ax, 2.0 * (space.tube.internal_halflength + off), 0.); let dir = Δ / (steps as f32); let traced = space.trace_iter(Ray { pos: a, dir }).nth(steps).unwrap(); assert_abs_diff_eq!(traced.pos.x, b.x, epsilon = 1.0e-2); @@ -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 b = vec2(x, space.tube.external_halflength + off); - let Δ = vec2(0.0, 2.0 * (space.tube.internal_halflength + off)); + let a = vec3(x, -(space.tube.external_halflength + off), 0.); + let b = vec3(x, space.tube.external_halflength + off, 0.); + let Δ = vec3(0.0, 2.0 * (space.tube.internal_halflength + off), 0.); let dir = Δ / (steps as f32); let traced = space.trace_iter(Ray { pos: a, dir }).nth(steps).unwrap(); assert_abs_diff_eq!(traced.pos.x, b.x, epsilon = 1.0e-1); @@ -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 b = vec2(x, space.tube.external_halflength + off); - let Δ = vec2(0.0, 2.0 * (space.tube.external_halflength + off)); + let a = vec3(x, -(space.tube.external_halflength + off), 0.); + let b = vec3(x, space.tube.external_halflength + off, 0.); + let Δ = vec3(0.0, 2.0 * (space.tube.external_halflength + off), 0.); let dir = Δ / (steps as f32); let traced = space.trace_iter(Ray { pos: a, dir }).nth(steps).unwrap(); assert_abs_diff_eq!(traced.pos.x, b.x, epsilon = 2.0e0); diff --git a/src/tube/mod.rs b/src/tube/mod.rs index 846a7a0..e7e31af 100644 --- a/src/tube/mod.rs +++ b/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 FlatTraceResult> { + fn obj_hitter(&self, pos: Vec3) -> Option 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, - globalize: impl Fn(Vec2) -> Vec2, + globalize: impl Fn(Vec3) -> Vec3, ) -> Vec { 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 { + pub fn line(&self, a: Vec3, b: Vec3, step: f32) -> Vec { 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, ray| { + let trace_to_flat = |points: &mut Vec, ray| { for ray in self.trace_iter(ray).skip(1) { points.push(ray.pos); if let Some(hitter) = self.obj_hitter(ray.pos) { @@ -241,7 +241,7 @@ impl DebugTraceable for Space { } struct Rect { - pub size: Vec2, + pub size: Vec3, } impl Rect { @@ -253,7 +253,7 @@ impl Rect { } } - fn is_inside(&self, pt: Vec2) -> bool { + fn is_inside(&self, pt: Vec3) -> bool { pt.abs().cmplt(self.size).all() } @@ -285,146 +285,146 @@ impl Rect { fn test_rect() { assert_eq!( Rect::flip_ray(Ray { - pos: vec2(2.0, 3.0), - dir: vec2(4.0, 5.0) + pos: vec3(2.0, 3.0, 2.0), + dir: vec3(4.0, 5.0, 4.0) }), Ray { - pos: vec2(2.0, 3.0), - dir: vec2(4.0, 5.0) + pos: vec3(2.0, 3.0, 2.0), + dir: vec3(4.0, 5.0, 4.0) } ); assert_eq!( Rect::flip_ray(Ray { - pos: vec2(2.0, 3.0), - dir: vec2(-4.0, 5.0) + pos: vec3(2.0, 3.0, 2.0), + dir: vec3(-4.0, 5.0, -4.0) }), Ray { - pos: vec2(-2.0, 3.0), - dir: vec2(4.0, 5.0) + pos: vec3(-2.0, 3.0, -2.0), + dir: vec3(4.0, 5.0, 4.0) } ); assert_eq!( Rect::flip_ray(Ray { - pos: vec2(2.0, 3.0), - dir: vec2(4.0, -5.0) + pos: vec3(2.0, 3.0, 2.0), + dir: vec3(4.0, -5.0, 4.0) }), Ray { - pos: vec2(2.0, -3.0), - dir: vec2(4.0, 5.0) + pos: vec3(2.0, -3.0, 2.0), + dir: vec3(4.0, 5.0, 4.0) } ); assert_eq!( Rect::flip_ray(Ray { - pos: vec2(2.0, 3.0), - dir: vec2(4.0, 0.0) + pos: vec3(2.0, 3.0, 2.0), + dir: vec3(4.0, 0.0, 4.0) }), Ray { - pos: vec2(2.0, 3.0), - dir: vec2(4.0, 0.0) + pos: vec3(2.0, 3.0, 2.0), + dir: vec3(4.0, 0.0, 4.0) } ); let r = Rect { - size: vec2(2.0, 3.0), + size: vec3(2.0, 3.0, 2.0), }; assert_eq!( r.trace_into(Ray { - pos: vec2(3.0, 3.0), - dir: vec2(1.0, 1.0) + pos: vec3(3.0, 3.0, 3.0), + dir: vec3(1.0, 1.0, 1.0) }), None ); assert_eq!( r.trace_into(Ray { - pos: vec2(-3.0, 2.0), - dir: vec2(1.0, 0.0) + pos: vec3(-3.0, 2.0, -3.0), + dir: vec3(1.0, 0.0, 1.0) }), Some(1.0) ); assert_eq!( r.trace_into(Ray { - pos: vec2(-3.0, 2.0), - dir: vec2(-1.0, 0.0) + pos: vec3(-3.0, 2.0, -3.0), + dir: vec3(-1.0, 0.0, -1.0) }), None ); assert_eq!( r.trace_into(Ray { - pos: vec2(-3.0, 1.0), - dir: vec2(2.0, 2.0) + pos: vec3(-3.0, 1.0, -3.0), + dir: vec3(2.0, 2.0, 2.0) }), Some(0.5) ); assert_eq!( r.trace_into(Ray { - pos: vec2(-3.0, 2.1), - dir: vec2(2.0, 2.0) + pos: vec3(-3.0, 2.1, -3.0), + dir: vec3(2.0, 2.0, 2.0) }), None ); assert_eq!( r.trace_into(Ray { - pos: vec2(2.0, 3.0), - dir: vec2(1.0, 1.0) + pos: vec3(2.0, 3.0, 2.0), + dir: vec3(1.0, 1.0, 1.0) }), None ); assert_eq!( r.trace_into(Ray { - pos: vec2(-2.0, 3.0), - dir: vec2(-1.0, 1.0) + pos: vec3(-2.0, 3.0, -2.0), + dir: vec3(-1.0, 1.0, -1.0) }), None ); assert_eq!( r.trace_into(Ray { - pos: vec2(2.0, 3.0), - dir: vec2(-1.0, -1.0) + pos: vec3(2.0, 3.0, 2.0), + dir: vec3(-1.0, -1.0, -1.0) }), Some(0.0) ); assert_eq!( r.trace_into(Ray { - pos: vec2(2.0, -3.0), - dir: vec2(-1.0, 1.0) + pos: vec3(2.0, -3.0, 2.0), + dir: vec3(-1.0, 1.0, -1.0) }), Some(0.0) ); assert_eq!( r.trace_out_of(Ray { - pos: vec2(0.0, 0.0), - dir: vec2(1.0, 1.0) + pos: vec3(0.0, 0.0, 0.0), + dir: vec3(1.0, 1.0, 1.0) }), Some(2.0) ); assert_eq!( r.trace_out_of(Ray { - pos: vec2(0.0, 0.0), - dir: vec2(0.0, 1.0) + pos: vec3(0.0, 0.0, 0.0), + dir: vec3(0.0, 1.0, 0.0) }), Some(3.0) ); assert_eq!( r.trace_out_of(Ray { - pos: vec2(0.0, 1.0), - dir: vec2(0.0, -1.0) + pos: vec3(0.0, 1.0, 0.0), + dir: vec3(0.0, -1.0, 0.0) }), Some(4.0) ); assert_eq!( r.trace_out_of(Ray { - pos: vec2(1.0, 1.0), - dir: vec2(0.0, -1.0) + pos: vec3(1.0, 1.0, 1.0), + dir: vec3(0.0, -1.0, 0.0) }), Some(4.0) ); assert_eq!( r.trace_out_of(Ray { - pos: vec2(2.0, 3.0), - dir: vec2(1.0, 1.0) + pos: vec3(2.0, 3.0, 2.0), + dir: vec3(1.0, 1.0, 1.0) }), Some(0.0) ); diff --git a/src/types.rs b/src/types.rs index 72e88dc..77ec8c4 100644 --- a/src/types.rs +++ b/src/types.rs @@ -1,9 +1,9 @@ -use glam::{f32, i32, Mat2, Vec2}; +use glam::{f32, i32, Mat3, Vec3}; #[derive(Copy, Clone, Debug, PartialEq)] pub struct Ray { - pub pos: Vec2, - pub dir: Vec2, + pub pos: Vec3, + pub dir: Vec3, } impl Ray { @@ -15,7 +15,7 @@ impl Ray { } } -impl std::ops::Mul for Mat2 { +impl std::ops::Mul for Mat3 { type Output = Ray; fn mul(self, rhs: Ray) -> Self::Output { @@ -29,9 +29,9 @@ impl std::ops::Mul 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, // в локальной ортонормированной СК объекта }