diff --git a/src/bin/flat.rs b/src/bin/flat.rs index d601954..03cd373 100644 --- a/src/bin/flat.rs +++ b/src/bin/flat.rs @@ -2,6 +2,7 @@ use flo_draw::*; use flo_canvas::*; use glm::*; use num_traits::identities::Zero; +use riemann::{Decomp2, Metric, trace_iter}; pub fn main() { let space = Coil { @@ -75,126 +76,138 @@ impl Metric for Coil { } } -struct Decomp2 { - ortho: Mat2, - diag: Vec2, -} +mod riemann { + use glm::*; + use num_traits::identities::Zero; -type Tens2 = [Mat2; 2]; - -trait Metric { - fn halfmetric(&self, pos: Vec2) -> Decomp2; - - fn metric(&self, pos: Vec2) -> Mat2 { - let h = self.halfmetric(pos); - transpose(&h.ortho) * diagonal(h.diag * h.diag) * h.ortho + pub struct Decomp2 { + pub ortho: Mat2, + pub diag: Vec2, } - fn dmetric(&self, pos: Vec2) -> Tens2 { - part_deriv(|p| self.metric(p), pos, 1.0e-3) - } + type Tens2 = [Mat2; 2]; - fn length(&self, at: Vec2, v: Vec2) -> f32 { - sqrt(dot(v, self.metric(at) * v)) - } + pub trait Metric { + fn halfmetric(&self, pos: Vec2) -> Decomp2; - fn normalize(&self, at: Vec2, v: Vec2) -> Vec2 { - v / self.length(at, v) - } - fn globalize(&self, at: Vec2, v: Vec2) -> Vec2 { - let h = self.halfmetric(at); - transpose(&h.ortho) * diagonal( Vec2::from_s(1.0) / h.diag) * h.ortho * v - } -} + fn metric(&self, pos: Vec2) -> Mat2 { + let h = self.halfmetric(pos); + transpose(&h.ortho) * diagonal(h.diag * h.diag) * h.ortho + } -struct TraceIter<'a, M: Metric> { - space: &'a M, - p: Vec2, - v: Vec2, - dt: f32, -} + fn dmetric(&self, pos: Vec2) -> Tens2 { + part_deriv(|p| self.metric(p), pos, 1.0e-3) + } -impl<'a, M: Metric> Iterator for TraceIter<'a, M> { - type Item = Vec2; + fn length(&self, at: Vec2, v: Vec2) -> f32 { + sqrt(dot(v, self.metric(at) * v)) + } - fn next(&mut self) -> Option { - let a: Vec2 = -convolute(krist(self.space, self.p), self.v); - self.v = self.v + a * self.dt; - self.p = self.p + self.v * self.dt; - Some(self.p) - } -} - -fn trace(space: &impl Metric, base: Vec2, dir: Vec2, distance: f32, dt: f32) -> Vec { - trace_iter(space, base, dir, dt).take((distance / dt) as usize).collect() -} - -fn trace_iter(space: &M, base: Vec2, dir: Vec2, dt: f32) -> TraceIter { - TraceIter{ - space: space, - p: base, - v: space.normalize(base, dir), - dt: dt, - } -} - -#[test] -fn t_iter() { - let space = Coil { - coil_scale: 2.0, - coil_r: 300.0, - coil_w: 50.0, - coil_m: 10.0, - }; - let base = vec2(-500.0, 0.0); - let dir = vec2(1.0, 0.3); - let dt = 1.0; - let steps = 1000; - let a = trace(&space, base, dir, dt * (steps as f32), dt); - let b: Vec = trace_iter(&space, base, dir, dt).take(steps).collect(); - assert_eq!(a, b); -} - -fn krist(space: &impl Metric, pos: Vec2) -> Tens2 { - // Γ^i_k_l = .5 * g^i^m * (g_m_k,l + g_m_l,k - g_k_l,m) - let g = inverse(&space.metric(pos)); // с верхними индексами - let d = space.dmetric(pos); - let mut ret: Tens2 = [Mat2::zero(); 2]; - // ret[i][l][k] = sum((m) => .5f * g[m][i] * (d[k][l][m] + d[l][k][m] - d[m][k][l])) - for i in 0..2 { - for l in 0..2 { - for k in 0..2 { - let mut v = 0.0; - for m in 0..2 { - v += g[m][i] * (d[l][k][m] + d[k][m][l] - d[m][k][l]); - } - ret[i][l][k] = 0.5 * v; - } + fn normalize(&self, at: Vec2, v: Vec2) -> Vec2 { + v / self.length(at, v) + } + fn globalize(&self, at: Vec2, v: Vec2) -> Vec2 { + let h = self.halfmetric(at); + transpose(&h.ortho) * diagonal(Vec2::from_s(1.0) / h.diag) * h.ortho * v } } - ret -} -fn dir_deriv(f: impl Fn(Vec2) -> Mat2, pos: Vec2, delta: Vec2) -> Mat2 { - (f(pos + delta) - f(pos - delta)) / (2.0 * length(delta)) -} + pub struct TraceIter<'a, M: Metric> { + space: &'a M, + p: Vec2, + v: Vec2, + dt: f32, + } -fn part_deriv(f: impl Fn(Vec2) -> Mat2, pos: Vec2, eps: f32) -> Tens2 { - [ - dir_deriv(&f, pos, vec2(eps, 0.0)), - dir_deriv(&f, pos, vec2(0.0, eps)), - ] -} + impl<'a, M: Metric> Iterator for TraceIter<'a, M> { + type Item = Vec2; -fn convolute(G: Tens2, v: Vec2) -> Vec2 { - vec2( - dot(v, G[0] * v), - dot(v, G[1] * v) - ) -} + fn next(&mut self) -> Option { + let a: Vec2 = -convolute(krist(self.space, self.p), self.v); + self.v = self.v + a * self.dt; + self.p = self.p + self.v * self.dt; + Some(self.p) + } + } -fn diagonal(v: Vec2) -> Mat2 { - mat2(v.x, 0.0, 0.0, v.y) + fn trace(space: &impl Metric, base: Vec2, dir: Vec2, distance: f32, dt: f32) -> Vec { + trace_iter(space, base, dir, dt).take((distance / dt) as usize).collect() + } + + pub fn trace_iter(space: &M, base: Vec2, dir: Vec2, dt: f32) -> TraceIter { + TraceIter { + space: space, + p: base, + v: space.normalize(base, dir), + dt: dt, + } + } + + #[cfg(test)] + mod test { + use glm::*; + use crate::riemann::{trace, trace_iter}; + use crate::Coil; + + #[test] + fn t_iter() { + let space = Coil { + coil_scale: 2.0, + coil_r: 300.0, + coil_w: 50.0, + coil_m: 10.0, + }; + let base = vec2(-500.0, 0.0); + let dir = vec2(1.0, 0.3); + let dt = 1.0; + let steps = 1000; + let a = trace(&space, base, dir, dt * (steps as f32), dt); + let b: Vec = trace_iter(&space, base, dir, dt).take(steps).collect(); + assert_eq!(a, b); + } + } + + fn krist(space: &impl Metric, pos: Vec2) -> Tens2 { + // Γ^i_k_l = .5 * g^i^m * (g_m_k,l + g_m_l,k - g_k_l,m) + let g = inverse(&space.metric(pos)); // с верхними индексами + let d = space.dmetric(pos); + let mut ret: Tens2 = [Mat2::zero(); 2]; + // ret[i][l][k] = sum((m) => .5f * g[m][i] * (d[k][l][m] + d[l][k][m] - d[m][k][l])) + for i in 0..2 { + for l in 0..2 { + for k in 0..2 { + let mut v = 0.0; + for m in 0..2 { + v += g[m][i] * (d[l][k][m] + d[k][m][l] - d[m][k][l]); + } + ret[i][l][k] = 0.5 * v; + } + } + } + ret + } + + fn dir_deriv(f: impl Fn(Vec2) -> Mat2, pos: Vec2, delta: Vec2) -> Mat2 { + (f(pos + delta) - f(pos - delta)) / (2.0 * length(delta)) + } + + fn part_deriv(f: impl Fn(Vec2) -> Mat2, pos: Vec2, eps: f32) -> Tens2 { + [ + dir_deriv(&f, pos, vec2(eps, 0.0)), + dir_deriv(&f, pos, vec2(0.0, eps)), + ] + } + + fn convolute(G: Tens2, v: Vec2) -> Vec2 { + vec2( + dot(v, G[0] * v), + dot(v, G[1] * v) + ) + } + + fn diagonal(v: Vec2) -> Mat2 { + mat2(v.x, 0.0, 0.0, v.y) + } } fn sqr(x: f32) -> f32 {