632 lines
19 KiB
Rust
632 lines
19 KiB
Rust
use flo_draw::*;
|
||
use flo_canvas::*;
|
||
use glam::*;
|
||
use riemann::{Decomp2, Metric, trace_iter};
|
||
use shape::Shape;
|
||
use Subspace::{Boundary, Inner, Outer};
|
||
|
||
#[cfg(test)]
|
||
use approx::assert_abs_diff_eq;
|
||
|
||
const DT: f32 = 0.1;
|
||
|
||
pub fn main() {
|
||
with_2d_graphics(move || {
|
||
let canvas = create_drawing_window("Refraction");
|
||
canvas.draw(|gc| {
|
||
let tube = Rect {
|
||
inner_radius: 30.0,
|
||
outer_radius: 50.0,
|
||
internal_halflength: 100.0,
|
||
external_halflength: 300.0,
|
||
};
|
||
|
||
let space = Space { rect: tube };
|
||
|
||
gc.canvas_height(1000.0);
|
||
tube.render(gc);
|
||
gc.line_width(0.5);
|
||
gc.stroke_color(Color::Rgba(1.0, 0.5, 0.0, 0.5));
|
||
draw_fan(gc, &tube, vec2(-500.0, 0.0), vec2(1.0, 0.0), 1.0);
|
||
gc.stroke_color(Color::Rgba(1.0, 0.5, 0.0, 1.0));
|
||
draw_fan_2(gc, &space, vec2(-500.0, 0.0), vec2(1.0, 0.0), 1.0);
|
||
gc.stroke_color(Color::Rgba(0.0, 0.5, 1.0, 0.5));
|
||
draw_fan(gc, &tube, vec2(0.0, -0.5 * tube.internal_halflength), vec2(1.0, 1.0), 1.0);
|
||
gc.stroke_color(Color::Rgba(0.0, 0.5, 1.0, 1.0));
|
||
draw_fan_2(gc, &space, vec2(0.0, -0.5 * tube.internal_halflength), vec2(1.0, 1.0), 1.0);
|
||
draw_track(gc, &space, vec2(-500.0, 0.0), vec2(1.0, 0.2));
|
||
draw_track(gc, &space, vec2(-500.0, 0.0), vec2(1.0, 0.5));
|
||
});
|
||
});
|
||
}
|
||
|
||
struct Location {
|
||
pos: Vec2,
|
||
rot: Mat2,
|
||
}
|
||
|
||
struct Space {
|
||
rect: Rect,
|
||
}
|
||
|
||
#[derive(PartialEq, Eq)]
|
||
enum Subspace {
|
||
Outer,
|
||
Boundary,
|
||
Inner,
|
||
}
|
||
|
||
impl Space {
|
||
fn which_subspace(&self, pt: Vec2) -> Subspace {
|
||
if pt.y.abs() > self.rect.external_halflength {
|
||
Outer
|
||
} else if pt.x.abs() > self.rect.outer_radius {
|
||
Outer
|
||
} else if pt.x.abs() > self.rect.inner_radius {
|
||
Boundary
|
||
} else {
|
||
Inner
|
||
}
|
||
}
|
||
|
||
fn flat_to_global(&self, at: Vec2, v: Vec2) -> Vec2 {
|
||
Mat2::from(self.rect.sqrt_at(at).inverse()) * v
|
||
}
|
||
|
||
fn global_to_flat(&self, at: Vec2, v: Vec2) -> Vec2 {
|
||
Mat2::from(self.rect.sqrt_at(at)) * v
|
||
}
|
||
|
||
/// Выполняет один шаг трассировки. Работает в любой части пространства, но вне Boundary доступны более эффективные методы.
|
||
fn trace_step(&self, ray: Ray) -> Ray {
|
||
let a: Vec2 = -riemann::contract2(riemann::krist(&self.rect, ray.pos), ray.dir);
|
||
let v = ray.dir + a * DT;
|
||
let p = ray.pos + v * DT;
|
||
Ray { pos: p, dir: v }
|
||
}
|
||
|
||
/// Выполняет один шаг перемещения. Работает в любой части пространства, но вне Boundary доступны более эффективные методы.
|
||
fn move_step(&self, loc: Location, off: Vec2) -> Location {
|
||
let corr = Mat2::IDENTITY - riemann::contract(riemann::krist(&self.rect, loc.pos), loc.rot * off);
|
||
let p = loc.pos + corr * loc.rot * off;
|
||
Location { pos: p, rot: corr * loc.rot }
|
||
}
|
||
}
|
||
|
||
fn draw_ray_2(gc: &mut Vec<Draw>, space: &Space, base: Vec2, dir: Vec2) {
|
||
let dir = space.rect.globalize(base, dir);
|
||
gc.new_path();
|
||
gc.move_to(base.x, base.y);
|
||
let dt = DT;
|
||
let mut p = base;
|
||
let mut v = space.rect.normalize_vec_at(base, dir);
|
||
for _ in 0..10000 {
|
||
let a: Vec2 = -riemann::contract2(riemann::krist(&space.rect, p), v);
|
||
v = v + a * dt;
|
||
p = p + v * dt;
|
||
gc.line_to(p.x, p.y);
|
||
if p.abs().cmpgt(Vec2::splat(1000.0)).any() {
|
||
break;
|
||
}
|
||
let sub = space.which_subspace(p);
|
||
if sub == Boundary {
|
||
continue;
|
||
}
|
||
gc.stroke();
|
||
gc.new_dash_pattern();
|
||
gc.dash_length(6.0);
|
||
gc.new_path();
|
||
gc.move_to(p.x, p.y);
|
||
match sub {
|
||
Inner => {
|
||
let cell = RectInside { rect: space.rect };
|
||
let ray = cell.ray_to_local(Ray { pos: p, dir: v });
|
||
let Some(ray) = cell.to_boundary(ray) else {
|
||
panic!("Can't get outta here!");
|
||
};
|
||
Ray { pos: p, dir: v } = cell.ray_to_global(ray);
|
||
}
|
||
Outer => {
|
||
let cell = basic_shapes::Rect { size: vec2(space.rect.outer_radius, space.rect.external_halflength) };
|
||
let Some(dist) = cell.trace_into(Ray { pos: p, dir: v }) else {
|
||
p += v * 1000.0;
|
||
gc.line_to(p.x, p.y);
|
||
gc.stroke();
|
||
break;
|
||
};
|
||
p += v * dist;
|
||
}
|
||
Boundary => panic!(),
|
||
}
|
||
gc.line_to(p.x, p.y);
|
||
gc.stroke();
|
||
gc.new_dash_pattern();
|
||
gc.new_path();
|
||
gc.move_to(p.x, p.y);
|
||
}
|
||
gc.stroke();
|
||
gc.new_dash_pattern();
|
||
}
|
||
|
||
fn draw_fan_2(gc: &mut Vec<Draw>, space: &Space, 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) {
|
||
draw_ray_2(gc, space, base, dir + y * v);
|
||
}
|
||
}
|
||
|
||
fn draw_ray(gc: &mut Vec<Draw>, space: &impl Metric, base: Vec2, dir: Vec2) {
|
||
let dir = space.globalize(base, dir);
|
||
gc.new_path();
|
||
gc.move_to(base.x, base.y);
|
||
for pt in trace_iter(space, base, dir, DT).take(10000) {
|
||
gc.line_to(pt.x, pt.y);
|
||
if pt.abs().cmpgt(Vec2::splat(1000.0)).any() {
|
||
break;
|
||
}
|
||
}
|
||
gc.stroke();
|
||
}
|
||
|
||
fn draw_track(gc: &mut Vec<Draw>, space: &Space, start: Vec2, dir: Vec2) {
|
||
const SCALE: f32 = 5.0;
|
||
const STEP: f32 = 2.0 * SCALE;
|
||
// let mut loc = Location { pos: start, rot: Mat2::IDENTITY };
|
||
// let dir = space.rect.globalize(start, dir);
|
||
// let v = space.rect.normalize(start, dir);
|
||
let mut loc = Location { pos: start, rot: mat2(dir, vec2(-dir.y, dir.x)) };
|
||
let v = vec2(1.0, 0.0);
|
||
let mut draw = |loc: &Location| {
|
||
let p = loc.pos;
|
||
let ax = p + loc.rot.x_axis * SCALE;
|
||
let ay = p + loc.rot.y_axis * SCALE;
|
||
gc.new_path();
|
||
gc.stroke_color(Color::Rgba(0.7, 0.0, 0.0, 1.0));
|
||
gc.move_to(p.x, p.y);
|
||
gc.line_to(ax.x, ax.y);
|
||
gc.stroke();
|
||
gc.new_path();
|
||
gc.stroke_color(Color::Rgba(0.0, 0.7, 0.0, 1.0));
|
||
gc.move_to(p.x, p.y);
|
||
gc.line_to(ay.x, ay.y);
|
||
gc.stroke();
|
||
};
|
||
draw(&loc);
|
||
for _ in 0..1000 {
|
||
let N = (STEP / DT).floor() as i32;
|
||
for _ in 0..N {
|
||
loc = space.move_step(loc, v * DT);
|
||
}
|
||
draw(&loc);
|
||
}
|
||
}
|
||
|
||
fn draw_fan(gc: &mut Vec<Draw>, space: &impl Metric, 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) {
|
||
draw_ray(gc, space, base, dir + y * v);
|
||
}
|
||
}
|
||
|
||
trait Renderable {
|
||
fn render(&self, gc: &mut Vec<Draw>);
|
||
}
|
||
|
||
impl Renderable for Rect {
|
||
fn render(&self, gc: &mut Vec<Draw>) {
|
||
gc.new_path();
|
||
gc.rect(-self.outer_radius, -self.external_halflength, self.outer_radius, self.external_halflength);
|
||
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();
|
||
}
|
||
}
|
||
|
||
#[derive(Debug, PartialEq)]
|
||
struct Ray {
|
||
pos: Vec2,
|
||
dir: Vec2,
|
||
}
|
||
|
||
mod basic_shapes {
|
||
use glam::{Vec2, vec2};
|
||
use crate::Ray;
|
||
use crate::shape::Shape;
|
||
|
||
pub struct Rect {
|
||
pub size: Vec2,
|
||
}
|
||
|
||
impl Rect {
|
||
/// Отражает луч, чтобы все координаты направления были положительны (допустимо благодаря симметрии Rect).
|
||
fn flip_ray(ray: Ray) -> Ray {
|
||
Ray { pos: ray.pos * ray.dir.signum(), dir: ray.dir.abs() }
|
||
}
|
||
}
|
||
|
||
impl Shape for Rect {
|
||
fn is_inside(&self, pt: Vec2) -> bool {
|
||
pt.abs().cmplt(self.size).all()
|
||
}
|
||
|
||
fn trace_into(&self, ray: Ray) -> Option<f32> {
|
||
let ray = Self::flip_ray(ray);
|
||
// ray.pos.x + t * ray.dir.x = −size.x
|
||
let ts = (-self.size - ray.pos) / ray.dir;
|
||
let t = ts.max_element();
|
||
let pt = ray.pos + t * ray.dir;
|
||
if t < 0.0 { return None; }
|
||
if pt.cmpgt(self.size).any() { return None; }
|
||
Some(t)
|
||
}
|
||
|
||
fn trace_out_of(&self, ray: Ray) -> Option<f32> {
|
||
let ray = Self::flip_ray(ray);
|
||
// ray.pos.x + t * ray.dir.x = +size.x
|
||
let ts = (self.size - ray.pos) / ray.dir;
|
||
let t = ts.min_element();
|
||
Some(t)
|
||
}
|
||
|
||
fn visualise(&self) -> Vec<Vec2> {
|
||
vec![vec2(-self.size.x, -self.size.y), vec2(self.size.x, -self.size.y), vec2(self.size.x, self.size.y), vec2(-self.size.x, self.size.y)]
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn test_rect() {
|
||
assert_eq!(Rect::flip_ray(Ray { pos: vec2(2.0, 3.0), dir: vec2(4.0, 5.0) }), Ray { pos: vec2(2.0, 3.0), dir: vec2(4.0, 5.0) });
|
||
assert_eq!(Rect::flip_ray(Ray { pos: vec2(2.0, 3.0), dir: vec2(-4.0, 5.0) }), Ray { pos: vec2(-2.0, 3.0), dir: vec2(4.0, 5.0) });
|
||
assert_eq!(Rect::flip_ray(Ray { pos: vec2(2.0, 3.0), dir: vec2(4.0, -5.0) }), Ray { pos: vec2(2.0, -3.0), dir: vec2(4.0, 5.0) });
|
||
assert_eq!(Rect::flip_ray(Ray { pos: vec2(2.0, 3.0), dir: vec2(4.0, 0.0) }), Ray { pos: vec2(2.0, 3.0), dir: vec2(4.0, 0.0) });
|
||
|
||
let r = Rect { size: vec2(2.0, 3.0) };
|
||
|
||
assert_eq!(r.trace_into(Ray { pos: vec2(3.0, 3.0), dir: vec2(1.0, 1.0) }), None);
|
||
assert_eq!(r.trace_into(Ray { pos: vec2(-3.0, 2.0), dir: vec2(1.0, 0.0) }), Some(1.0));
|
||
assert_eq!(r.trace_into(Ray { pos: vec2(-3.0, 2.0), dir: vec2(-1.0, 0.0) }), None);
|
||
assert_eq!(r.trace_into(Ray { pos: vec2(-3.0, 1.0), dir: vec2(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) }), None);
|
||
|
||
assert_eq!(r.trace_into(Ray { pos: vec2(2.0, 3.0), dir: vec2(1.0, 1.0) }), None);
|
||
assert_eq!(r.trace_into(Ray { pos: vec2(-2.0, 3.0), dir: vec2(-1.0, 1.0) }), None);
|
||
assert_eq!(r.trace_into(Ray { pos: vec2(2.0, 3.0), dir: vec2(-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) }), Some(0.0));
|
||
|
||
assert_eq!(r.trace_out_of(Ray { pos: vec2(0.0, 0.0), dir: vec2(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) }), Some(3.0));
|
||
assert_eq!(r.trace_out_of(Ray { pos: vec2(0.0, 1.0), dir: vec2(0.0, -1.0) }), Some(4.0));
|
||
assert_eq!(r.trace_out_of(Ray { pos: vec2(1.0, 1.0), dir: vec2(0.0, -1.0) }), Some(4.0));
|
||
assert_eq!(r.trace_out_of(Ray { pos: vec2(2.0, 3.0), dir: vec2(1.0, 1.0) }), Some(0.0));
|
||
}
|
||
}
|
||
|
||
mod shape {
|
||
use glam::Vec2;
|
||
use crate::Ray;
|
||
|
||
pub trait Shape {
|
||
fn is_inside(&self, pt: Vec2) -> bool;
|
||
|
||
/// Ищет ближайшее пересечение луча с границей в направлении внутрь контура. Возвращает расстояние (в ray.dir).
|
||
fn trace_into(&self, ray: Ray) -> Option<f32>;
|
||
/// Ищет ближайшее пересечение луча с границей в направлении вовне контура. Возвращает расстояние (в ray.dir).
|
||
fn trace_out_of(&self, ray: Ray) -> Option<f32>;
|
||
|
||
/// Возвращает визуальное представление контура, для отладки.
|
||
fn visualise(&self) -> Vec<Vec2>;
|
||
}
|
||
}
|
||
|
||
trait FlatCell: std::fmt::Debug {
|
||
fn pos_to_global(&self, pos: Vec2) -> Vec2;
|
||
fn pos_to_local(&self, pos: Vec2) -> Vec2;
|
||
fn ray_to_global(&self, ray: Ray) -> Ray;
|
||
fn ray_to_local(&self, ray: Ray) -> Ray;
|
||
|
||
fn is_inside(&self, pos: Vec2) -> bool {
|
||
let bnd = self.local_bounds();
|
||
pos.cmpge(bnd.0).all() && pos.cmple(bnd.1).all()
|
||
}
|
||
fn local_bounds(&self) -> (Vec2, Vec2);
|
||
|
||
fn to_boundary(&self, ray: Ray) -> Option<Ray> {
|
||
assert!(self.is_inside(ray.pos));
|
||
let sgn = ray.dir.signum();
|
||
let p = ray.pos * sgn;
|
||
let v = ray.dir * sgn;
|
||
let mut bnd = self.local_bounds();
|
||
if sgn.x < 0.0 {
|
||
(bnd.0.x, bnd.1.x) = (-bnd.1.x, -bnd.0.x);
|
||
}
|
||
if sgn.y < 0.0 {
|
||
(bnd.0.y, bnd.1.y) = (-bnd.1.y, -bnd.0.y);
|
||
}
|
||
let t = if (bnd.1.x - p.x) * v.y <= (bnd.1.y - p.y) * v.x {
|
||
(bnd.1.x - p.x) / v.x
|
||
} else {
|
||
(bnd.1.y - p.y) / v.y
|
||
};
|
||
if t <= 1000.0 {
|
||
Some(Ray { pos: sgn * (p + v * t), dir: sgn * v })
|
||
} else {
|
||
None
|
||
}
|
||
}
|
||
}
|
||
|
||
#[derive(Debug)]
|
||
struct RectInside {
|
||
rect: Rect,
|
||
}
|
||
|
||
impl FlatCell for RectInside {
|
||
fn pos_to_global(&self, pos: Vec2) -> Vec2 {
|
||
vec2(pos.x, self.rect.x(pos.y))
|
||
}
|
||
|
||
fn pos_to_local(&self, pos: Vec2) -> Vec2 {
|
||
vec2(pos.x, self.rect.u(pos.y))
|
||
}
|
||
|
||
fn ray_to_global(&self, ray: Ray) -> Ray {
|
||
Ray {
|
||
pos: self.pos_to_global(ray.pos),
|
||
dir: vec2(ray.dir.x, self.rect.dx(ray.pos.y, ray.dir.y)),
|
||
}
|
||
}
|
||
|
||
fn ray_to_local(&self, ray: Ray) -> Ray {
|
||
Ray {
|
||
pos: self.pos_to_local(ray.pos),
|
||
dir: vec2(ray.dir.x, self.rect.du(ray.pos.y, ray.dir.y)),
|
||
}
|
||
}
|
||
|
||
fn local_bounds(&self) -> (Vec2, Vec2) {
|
||
(vec2(-self.rect.inner_radius, -self.rect.internal_halflength), vec2(self.rect.inner_radius, self.rect.internal_halflength))
|
||
}
|
||
}
|
||
|
||
#[derive(Copy, Clone, Debug)]
|
||
struct Rect {
|
||
outer_radius: f32,
|
||
inner_radius: f32,
|
||
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 sqrt_at(&self, pos: Vec2) -> Decomp2 {
|
||
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, sy.lerp(1.0, sx)),
|
||
}
|
||
}
|
||
}
|
||
|
||
mod riemann {
|
||
use glam::*;
|
||
|
||
pub struct Decomp2 {
|
||
pub ortho: Mat2,
|
||
pub diag: Vec2,
|
||
}
|
||
|
||
impl Decomp2 {
|
||
fn square(&self) -> Self {
|
||
Self {
|
||
ortho: self.ortho,
|
||
diag: self.diag * self.diag,
|
||
}
|
||
}
|
||
|
||
pub(crate) fn inverse(&self) -> Self {
|
||
Self {
|
||
ortho: self.ortho,
|
||
diag: Vec2::splat(1.0) / self.diag,
|
||
}
|
||
}
|
||
}
|
||
|
||
impl From<Decomp2> for Mat2 {
|
||
fn from(value: Decomp2) -> Self {
|
||
value.ortho.transpose() * Mat2::from_diagonal(value.diag) * value.ortho
|
||
}
|
||
}
|
||
|
||
type Tens2 = [Mat2; 2];
|
||
|
||
pub trait Metric {
|
||
fn sqrt_at(&self, pos: Vec2) -> Decomp2;
|
||
|
||
fn at(&self, pos: Vec2) -> Mat2 {
|
||
self.sqrt_at(pos).square().into()
|
||
}
|
||
|
||
fn inverse_at(&self, pos: Vec2) -> Mat2 {
|
||
self.sqrt_at(pos).square().inverse().into()
|
||
}
|
||
|
||
fn part_derivs_at(&self, pos: Vec2) -> Tens2 {
|
||
part_deriv(|p| self.at(p), pos, 1.0e-3)
|
||
}
|
||
|
||
fn vec_length_at(&self, at: Vec2, v: Vec2) -> f32 {
|
||
v.dot(self.at(at) * v).sqrt()
|
||
}
|
||
|
||
fn normalize_vec_at(&self, at: Vec2, v: Vec2) -> Vec2 {
|
||
v / self.vec_length_at(at, v)
|
||
}
|
||
|
||
fn globalize(&self, at: Vec2, v: Vec2) -> Vec2 {
|
||
Mat2::from(self.sqrt_at(at).inverse()) * v
|
||
}
|
||
}
|
||
|
||
pub struct TraceIter<'a, M: Metric> {
|
||
space: &'a M,
|
||
p: Vec2,
|
||
v: Vec2,
|
||
dt: f32,
|
||
}
|
||
|
||
impl<'a, M: Metric> Iterator for TraceIter<'a, M> {
|
||
type Item = Vec2;
|
||
|
||
fn next(&mut self) -> Option<Self::Item> {
|
||
let a: Vec2 = -contract2(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)
|
||
}
|
||
}
|
||
|
||
pub fn trace_iter<M: Metric>(space: &M, base: Vec2, dir: Vec2, dt: f32) -> TraceIter<M> {
|
||
TraceIter {
|
||
space,
|
||
p: base,
|
||
v: space.normalize_vec_at(base, dir),
|
||
dt,
|
||
}
|
||
}
|
||
|
||
pub 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 = &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| 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>())
|
||
}
|
||
|
||
fn dir_deriv(f: impl Fn(Vec2) -> Mat2, pos: Vec2, delta: Vec2) -> Mat2 {
|
||
(f(pos + delta) - f(pos - delta)) / (2.0 * delta.length())
|
||
}
|
||
|
||
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)),
|
||
]
|
||
}
|
||
|
||
/// Сворачивает тензор t с вектором u
|
||
pub fn contract(t: Tens2, u: Vec2) -> Mat2 {
|
||
mat2(t[0] * u, t[1] * u).transpose()
|
||
}
|
||
|
||
/// Сворачивает тензор t с вектором v дважды, по второму и третьему индексам.
|
||
pub fn contract2(t: Tens2, v: Vec2) -> Vec2 {
|
||
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_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_tens2(f: impl Fn(usize, usize, usize) -> f32) -> Tens2 {
|
||
std::array::from_fn(|i| make_mat2(|j, k| f(i, j, k)))
|
||
}
|
||
|
||
#[test]
|
||
fn m2() {
|
||
let m = make_mat2(|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);
|
||
assert_eq!(m.col(1)[1], 3.0);
|
||
}
|
||
|
||
#[test]
|
||
fn t2() {
|
||
let t = make_tens2(|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);
|
||
assert_eq!(t[1].col(1)[0], 3.0);
|
||
assert_eq!(t[0].col(0)[1], 4.0);
|
||
assert_eq!(t[1].col(0)[1], 5.0);
|
||
assert_eq!(t[0].col(1)[1], 6.0);
|
||
assert_eq!(t[1].col(1)[1], 7.0);
|
||
}
|
||
}
|
||
|
||
fn smoothstep(x: f32) -> f32 {
|
||
return 3.0 * x * x - 2.0 * x * x * x;
|
||
}
|
||
|
||
fn smoothbox(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))
|
||
}
|