766 lines
25 KiB
Rust
766 lines
25 KiB
Rust
use std::f32::consts::{FRAC_PI_2, PI};
|
||
use flo_draw::*;
|
||
use flo_canvas::*;
|
||
use glam::*;
|
||
|
||
mod riemann;
|
||
|
||
use riemann::{Tens2, Decomp2, Metric, trace_iter};
|
||
use shape::Shape;
|
||
use Subspace::{Boundary, Inner, Outer};
|
||
|
||
const DT: f32 = 0.1;
|
||
|
||
fn draw_loop(gc: &mut Vec<Draw>, mut pts: impl Iterator<Item=Vec2>) {
|
||
gc.new_path();
|
||
let Some(first) = pts.next() else { return; };
|
||
gc.move_to(first.x, first.y);
|
||
for pt in pts {
|
||
gc.line_to(pt.x, pt.y);
|
||
}
|
||
gc.close_path();
|
||
gc.stroke();
|
||
}
|
||
|
||
pub fn main() {
|
||
with_2d_graphics(move || {
|
||
let canvas = create_drawing_window("Refraction");
|
||
canvas.draw(|gc| {
|
||
let tube = Tube {
|
||
inner_radius: 30.0,
|
||
outer_radius: 50.0,
|
||
internal_halflength: 100.0,
|
||
external_halflength: 300.0,
|
||
};
|
||
|
||
let objs: Vec<_> = [-1.25, -1.00, -0.85, -0.50, 0.00, 0.40, 0.70, 0.95, 1.05]
|
||
.iter()
|
||
.enumerate()
|
||
.map(|(k, &y)| Object {
|
||
id: k as i32,
|
||
loc: {
|
||
let pos = vec2(0.0, y * tube.external_halflength);
|
||
let adj: Mat2 = tube.sqrt_at(pos).inverse().into();
|
||
let rot = Mat2::from_angle(y);
|
||
Location {
|
||
pos,
|
||
rot: adj * rot,
|
||
}
|
||
},
|
||
r: 20.0,
|
||
})
|
||
.collect();
|
||
let space = Space { tube, objs };
|
||
|
||
gc.canvas_height(500.0);
|
||
gc.transform(Transform2D::rotate(FRAC_PI_2));
|
||
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.5, 1.0, 0.0, 1.0));
|
||
draw_fan_2(gc, &space, vec2(-2.5 * tube.outer_radius, 1.25 * tube.external_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));
|
||
draw_track(gc, &space, vec2(-0.5 * tube.inner_radius, -1.25 * tube.external_halflength), vec2(0.1, 1.0));
|
||
|
||
let circle_segments = 47;
|
||
for obj in &space.objs {
|
||
let pos = obj.loc.pos;
|
||
gc.new_path();
|
||
gc.circle(pos.x, pos.y, 5.0);
|
||
gc.fill_color(Color::Rgba(0.0, 0.5, 1.0, 1.0));
|
||
gc.fill();
|
||
|
||
gc.stroke_color(Color::Rgba(0.0, 0.0, 0.0, 0.5));
|
||
draw_loop(gc, itertools_num::linspace(0.0, 2.0 * PI, circle_segments).skip(1).map(|φ| {
|
||
let dir = Vec2::from_angle(φ) * obj.r;
|
||
let dir = obj.loc.rot * dir;
|
||
pos + dir
|
||
}));
|
||
gc.stroke_color(Color::Rgba(0.0, 0.5, 1.0, 0.5));
|
||
draw_loop(gc, itertools_num::linspace(0.0, 2.0 * PI, circle_segments).skip(1).map(|φ| {
|
||
let dir = Vec2::from_angle(φ) * obj.r;
|
||
let dir = obj.loc.rot * dir;
|
||
space.trace_step(Ray { pos, dir }).pos
|
||
}));
|
||
gc.stroke_color(Color::Rgba(0.5, 0.0, 1.0, 1.0));
|
||
draw_loop(gc, itertools_num::linspace(0.0, 2.0 * PI, circle_segments).skip(1).map(|φ| {
|
||
let n = obj.r.floor();
|
||
let d = obj.r / n;
|
||
let dir = Vec2::from_angle(φ);
|
||
let dir = obj.loc.rot * dir * d;
|
||
space.trace_iter(Ray { pos, dir }).nth(n as usize).unwrap().pos
|
||
}));
|
||
}
|
||
});
|
||
});
|
||
}
|
||
|
||
#[derive(Copy, Clone, Debug)]
|
||
struct Location {
|
||
/// Положение в основной СК
|
||
pos: Vec2,
|
||
/// Преобразование вектора из локальной ортонормированной в основную СК
|
||
rot: Mat2,
|
||
}
|
||
|
||
#[derive(Copy, Clone, Debug)]
|
||
struct Object {
|
||
id: i32,
|
||
loc: Location,
|
||
r: f32,
|
||
}
|
||
|
||
struct Space {
|
||
tube: Tube,
|
||
objs: Vec<Object>,
|
||
}
|
||
|
||
#[derive(PartialEq, Eq, Debug)]
|
||
enum Subspace {
|
||
Outer,
|
||
Boundary,
|
||
Inner,
|
||
}
|
||
|
||
struct Hit {
|
||
distance: f32,
|
||
id: i32,
|
||
pos: Vec2, // положение в основной СК
|
||
rel: Ray, // в локальной ортонормированной СК объекта
|
||
}
|
||
|
||
struct FlatTraceResult {
|
||
end: Option<Ray>,
|
||
objects: Vec<Hit>,
|
||
}
|
||
|
||
impl Space {
|
||
fn which_subspace(&self, pt: Vec2) -> Subspace {
|
||
if pt.y.abs() > self.tube.external_halflength {
|
||
Outer
|
||
} else if pt.x.abs() > self.tube.outer_radius {
|
||
Outer
|
||
} else if pt.x.abs() > self.tube.inner_radius {
|
||
Boundary
|
||
} else {
|
||
Inner
|
||
}
|
||
}
|
||
|
||
fn flat_to_global(&self, at: Vec2) -> Mat2 {
|
||
Mat2::from(self.tube.sqrt_at(at).inverse())
|
||
}
|
||
|
||
fn global_to_flat(&self, at: Vec2) -> Mat2 {
|
||
Mat2::from(self.tube.sqrt_at(at))
|
||
}
|
||
|
||
/// Выполняет один шаг трассировки. Работает в любой части пространства, но вне Boundary доступны более эффективные методы.
|
||
/// ray задаётся в основной СК.
|
||
fn trace_step(&self, ray: Ray) -> Ray {
|
||
let a: Vec2 = -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 }
|
||
}
|
||
|
||
/// Выполняет один шаг перемещения. Работает в любой части пространства.
|
||
/// off задаётся в локальной СК. Рекомендуется считать небольшими шагами.
|
||
fn move_step(&self, loc: Location, off: Vec2) -> Location {
|
||
let corr = Mat2::IDENTITY - riemann::contract(riemann::krist(&self.tube, loc.pos), loc.rot * off);
|
||
let p = loc.pos + corr * loc.rot * off;
|
||
Location { pos: p, rot: corr * loc.rot }
|
||
}
|
||
|
||
fn trace_iter(&self, ray: Ray) -> impl Iterator<Item=Ray> + '_ {
|
||
std::iter::successors(Some(ray), |&ray| Some(self.trace_step(ray)))
|
||
}
|
||
|
||
fn trace_inner(&self, ray: Ray) -> FlatTraceResult {
|
||
assert_eq!(self.which_subspace(ray.pos), Inner);
|
||
let cell = TubeInside { tube: self.tube };
|
||
let ray = cell.ray_to_local(ray);
|
||
let objs = self.list_objects_inner();
|
||
let dist = cell.to_boundary(ray).expect("Can't get outta here!");
|
||
FlatTraceResult {
|
||
end: Some(cell.ray_to_global(ray.forward(dist))),
|
||
objects: Self::hit_objects(objs.as_slice(), ray, dist, |pos| cell.pos_to_global(pos)),
|
||
}
|
||
}
|
||
|
||
fn trace_outer(&self, ray: Ray) -> FlatTraceResult {
|
||
assert_eq!(self.which_subspace(ray.pos), Outer);
|
||
let cell = basic_shapes::Rect { size: vec2(self.tube.outer_radius, self.tube.external_halflength) };
|
||
let objs = self.list_objects_outer();
|
||
let lim = cell.trace_into(ray);
|
||
let dist = lim.unwrap_or(f32::INFINITY);
|
||
FlatTraceResult {
|
||
end: lim.map(|dist| ray.forward(dist)),
|
||
objects: Self::hit_objects(objs.as_slice(), ray, dist, |pos| pos),
|
||
}
|
||
}
|
||
|
||
fn trace_boundary(&self, ray: Ray) -> Ray {
|
||
assert_eq!(self.which_subspace(ray.pos), Boundary);
|
||
self.trace_iter(ray)
|
||
.find(|&ray| self.which_subspace(ray.pos) != Boundary)
|
||
.expect("Can't get outta the wall!")
|
||
}
|
||
|
||
fn list_objects(&self, tfm: impl Fn(Location) -> Location) -> Vec<Object> {
|
||
self.objs.iter().map(|&Object { id, loc, r }| Object { id, loc: tfm(loc), r }).collect()
|
||
}
|
||
|
||
fn list_objects_outer(&self) -> Vec<Object> {
|
||
self.list_objects(|loc|
|
||
match self.which_subspace(loc.pos) {
|
||
Outer => loc,
|
||
Inner => {
|
||
let Vec2 { x: u, y } = loc.pos; // в основной СК
|
||
let v = self.tube.v(y) + y.signum() * (self.tube.external_halflength - self.tube.internal_halflength);
|
||
Location {
|
||
pos: vec2(u, v), // в плоском продолжении СК Outer на область Inner
|
||
rot: self.global_to_flat(loc.pos) * loc.rot,
|
||
}
|
||
}
|
||
Boundary => panic!("Object at {} was destroyed by the space curvature", loc.pos),
|
||
})
|
||
}
|
||
|
||
fn list_objects_inner(&self) -> Vec<Object> {
|
||
self.list_objects(|Location { pos, rot }|
|
||
match self.which_subspace(pos) {
|
||
Inner | Outer => {
|
||
// NB: не работает для частей Outer с |y| < external_halflength. Но они и не нужны.
|
||
Location {
|
||
pos: vec2(pos.x, self.tube.v(pos.y)), // в плоской СК для Inner или её продолжении на Outer
|
||
rot: self.global_to_flat(pos) * rot,
|
||
}
|
||
}
|
||
Boundary => panic!("Object at {} was destroyed by the space curvature", pos),
|
||
})
|
||
}
|
||
|
||
fn hit_objects(objs: &[Object], ray: Ray, limit: f32, globalize: impl Fn(Vec2) -> Vec2) -> Vec<Hit> {
|
||
objs.iter()
|
||
.filter_map(|obj| {
|
||
let rel = ray.pos - obj.loc.pos;
|
||
let diff = rel.dot(ray.dir).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))
|
||
} else {
|
||
None
|
||
}
|
||
})
|
||
.filter(|&(_, t)| t >= 0.0 && t < limit)
|
||
.map(|(obj, t)| {
|
||
let pos = ray.forward(t).pos;
|
||
let rel = obj.loc.rot.inverse() * Ray { pos: pos - obj.loc.pos, dir: ray.dir };
|
||
Hit { id: obj.id, distance: t, pos: globalize(pos), rel }
|
||
})
|
||
.collect()
|
||
}
|
||
|
||
fn line(&self, a: Vec2, b: Vec2, step: f32) -> Vec<Vec2> {
|
||
match self.which_subspace(a) {
|
||
Outer => vec![b],
|
||
Inner => {
|
||
let cell = TubeInside { tube: self.tube };
|
||
let n = ((b - a).length() / step) as usize + 1;
|
||
let a = cell.pos_to_local(a);
|
||
let b = cell.pos_to_local(b);
|
||
(1..=n).map(|k| cell.pos_to_global(a.lerp(b, k as f32 / n as f32))).collect()
|
||
}
|
||
Boundary => panic!("Can't draw a line here!"),
|
||
}
|
||
}
|
||
}
|
||
|
||
fn draw_ray_2(gc: &mut Vec<Draw>, space: &Space, base: Vec2, dir: Vec2) {
|
||
let mut hits = Vec::<Draw>::new();
|
||
let dir = space.tube.globalize(base, dir);
|
||
gc.new_path();
|
||
gc.move_to(base.x, base.y);
|
||
let mut ray = Ray { pos: base, dir: space.tube.normalize_vec_at(base, dir) * DT };
|
||
for _ in 0..10000 {
|
||
ray = space.trace_step(ray);
|
||
gc.line_to(ray.pos.x, ray.pos.y);
|
||
if ray.pos.abs().cmpgt(Vec2::splat(1000.0)).any() {
|
||
break;
|
||
}
|
||
let sub = space.which_subspace(ray.pos);
|
||
if sub == Boundary {
|
||
continue;
|
||
}
|
||
gc.stroke();
|
||
gc.new_dash_pattern();
|
||
// gc.dash_length(6.0);
|
||
gc.new_path();
|
||
gc.move_to(ray.pos.x, ray.pos.y);
|
||
let ret = match sub {
|
||
Inner => space.trace_inner(ray),
|
||
Outer => space.trace_outer(ray),
|
||
Boundary => panic!(),
|
||
};
|
||
for hit in ret.objects {
|
||
let obj = space.objs[hit.id as usize];
|
||
hits.move_to(obj.loc.pos.x, obj.loc.pos.y);
|
||
for pt in trace_iter(&space.tube, obj.loc.pos, obj.loc.rot * hit.rel.pos, hit.rel.pos.length() / 100.0).take(100) {
|
||
hits.line_to(pt.x, pt.y);
|
||
}
|
||
hits.circle(hit.pos.x, hit.pos.y, 1.5);
|
||
let Ray { pos: rel, dir } = hit.rel;
|
||
let diff = rel.dot(dir).powi(2) - dir.length_squared() * (rel.length_squared() - obj.r.powi(2));
|
||
assert!(diff >= 0.0);
|
||
let t = (-rel.dot(dir) + diff.sqrt()) / dir.length_squared();
|
||
let rel2 = hit.rel.forward(t).pos;
|
||
let pos2 = trace_iter(&space.tube, obj.loc.pos, obj.loc.rot * rel2, rel2.length() / 100.0).nth(100).unwrap();
|
||
hits.move_to(pos2.x - 1.0, pos2.y - 1.0);
|
||
hits.line_to(pos2.x + 1.0, pos2.y + 1.0);
|
||
hits.move_to(pos2.x - 1.0, pos2.y + 1.0);
|
||
hits.line_to(pos2.x + 1.0, pos2.y - 1.0);
|
||
}
|
||
let a = ray.pos;
|
||
ray = match ret.end {
|
||
Some(r) => r,
|
||
None => {
|
||
ray = ray.forward(1000.0 / DT);
|
||
gc.line_to(ray.pos.x, ray.pos.y);
|
||
break;
|
||
}
|
||
};
|
||
for p in space.line(a, ray.pos, 10.0) {
|
||
gc.line_to(p.x, p.y);
|
||
}
|
||
gc.stroke();
|
||
gc.new_dash_pattern();
|
||
gc.new_path();
|
||
gc.move_to(ray.pos.x, ray.pos.y);
|
||
}
|
||
gc.stroke();
|
||
gc.new_path();
|
||
gc.new_dash_pattern();
|
||
gc.append(&mut hits);
|
||
gc.stroke();
|
||
}
|
||
|
||
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.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)) };
|
||
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 Tube {
|
||
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(Copy, Clone, Debug, PartialEq)]
|
||
struct Ray {
|
||
pos: Vec2,
|
||
dir: Vec2,
|
||
}
|
||
|
||
impl Ray {
|
||
fn forward(&self, dist: f32) -> Ray {
|
||
Ray { pos: self.pos + self.dir * dist, dir: self.dir }
|
||
}
|
||
}
|
||
|
||
impl std::ops::Mul<Ray> for Mat2 {
|
||
type Output = Ray;
|
||
|
||
fn mul(self, rhs: Ray) -> Self::Output {
|
||
Ray { pos: self * rhs.pos, dir: self * rhs.dir }
|
||
}
|
||
}
|
||
|
||
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<f32> {
|
||
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 <= 100000.0 {
|
||
Some(t)
|
||
} else {
|
||
None
|
||
}
|
||
}
|
||
}
|
||
|
||
#[derive(Debug)]
|
||
struct TubeInside {
|
||
tube: Tube,
|
||
}
|
||
|
||
impl FlatCell for TubeInside {
|
||
fn pos_to_global(&self, pos: Vec2) -> Vec2 {
|
||
vec2(pos.x, self.tube.y(pos.y))
|
||
}
|
||
|
||
fn pos_to_local(&self, pos: Vec2) -> Vec2 {
|
||
vec2(pos.x, self.tube.v(pos.y))
|
||
}
|
||
|
||
fn ray_to_global(&self, ray: Ray) -> Ray {
|
||
Ray {
|
||
pos: self.pos_to_global(ray.pos),
|
||
dir: vec2(ray.dir.x, self.tube.dy(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.tube.dv(ray.pos.y) * ray.dir.y),
|
||
}
|
||
}
|
||
|
||
fn local_bounds(&self) -> (Vec2, Vec2) {
|
||
(vec2(-self.tube.inner_radius, -self.tube.internal_halflength), vec2(self.tube.inner_radius, self.tube.internal_halflength))
|
||
}
|
||
}
|
||
|
||
#[derive(Copy, Clone, Debug)]
|
||
struct Tube {
|
||
outer_radius: f32,
|
||
inner_radius: f32,
|
||
external_halflength: f32,
|
||
internal_halflength: f32,
|
||
}
|
||
|
||
impl Tube {
|
||
fn fx(&self) -> fns::TubeX { fns::TubeX { min: self.inner_radius, max: self.outer_radius } }
|
||
fn fy(&self) -> fns::TubeY { fns::TubeY { internal: self.internal_halflength, external: self.external_halflength } }
|
||
|
||
pub fn y(&self, v: f32) -> f32 { self.fy().x(v) }
|
||
pub fn v(&self, y: f32) -> f32 { self.fy().u(y) }
|
||
pub fn dy(&self, v: f32) -> f32 { self.fy().dx(v) }
|
||
pub fn dv(&self, y: f32) -> f32 { self.fy().du(y) }
|
||
}
|
||
|
||
mod fns {
|
||
use crate::FloatExt2;
|
||
|
||
#[cfg(test)]
|
||
use approx::abs_diff_eq;
|
||
|
||
pub struct TubeX {
|
||
pub min: f32,
|
||
pub max: f32,
|
||
}
|
||
|
||
impl TubeX {
|
||
pub fn value(&self, x: f32) -> f32 { (self.min, self.max).inverse_lerp(x.abs()).clamp(0.0, 1.0) }
|
||
pub fn derivative(&self, x: f32) -> f32 { if x.abs() > self.min && x.abs() < self.max { x.signum() / (self.max - self.min) } else { 0.0 } }
|
||
}
|
||
|
||
pub struct TubeY {
|
||
pub internal: f32,
|
||
pub external: f32,
|
||
}
|
||
|
||
/// Продолжает функцию f с [-lim, lim] линейно в предположении f(±lim) = ±val, f'(±lim) = 1.
|
||
fn extend_linear(t: f32, f: impl FnOnce(f32) -> f32, lim: f32, val: f32) -> f32 {
|
||
if t.abs() <= lim { f(t) } else { t + t.signum() * (val - lim) }
|
||
}
|
||
|
||
/// Продолжает функцию f с [-lim, lim] константой в предположении f(±lim) = val, f'(±lim) = 0.
|
||
fn extend_const(t: f32, f: impl FnOnce(f32) -> f32, lim: f32, val: f32) -> f32 {
|
||
if t.abs() <= lim { f(t) } else { val }
|
||
}
|
||
|
||
impl TubeY {
|
||
fn a(&self) -> f32 { -(self.external - self.internal) / self.internal.powi(2) }
|
||
fn b(&self) -> f32 { 2.0 * self.external / self.internal - 1.0 }
|
||
fn root(&self, x: f32) -> f32 { (self.b().powi(2) + 4.0 * self.a() * x.abs()).sqrt() }
|
||
|
||
pub fn x(&self, u: f32) -> f32 { extend_linear(u, |u| (self.a() * u.abs() + self.b()) * u, self.internal, self.external) }
|
||
pub fn u(&self, x: f32) -> f32 { extend_linear(x, |x| 0.5 * x.signum() * (-self.b() + self.root(x)) / self.a(), self.external, self.internal) }
|
||
pub fn dx(&self, u: f32) -> f32 { extend_const(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) }
|
||
pub fn d2u(&self, x: f32) -> f32 { extend_const(x, |x| -2.0 * x.signum() * self.a() * self.root(x).powi(-3), self.external, 0.0) }
|
||
}
|
||
|
||
#[test]
|
||
fn test_tube_y() {
|
||
let testee = TubeY { internal: 100.0, external: 150.0 };
|
||
let ε = 1.0e-4f32;
|
||
let δ = 1.0 / 8.0; // Mathematically, you want this to be small. Computationally, you don’t.
|
||
let margin = 1.0 / 16.0;
|
||
let mul = 1.0 + margin;
|
||
for x in itertools_num::linspace(-mul * testee.external, mul * testee.external, 100) {
|
||
let u = testee.u(x);
|
||
assert!(abs_diff_eq!(x, testee.x(u), epsilon = ε), "At x={}:\nu(x): {}\nx(u(x)): {}\n", x, u, testee.x(u));
|
||
|
||
let du = (testee.u(x + δ) - testee.u(x - δ)) / (2. * δ);
|
||
assert!(abs_diff_eq!(du, testee.du(x), epsilon = ε), "At x={}, u':\nexpected: {}\nactual: {}\n", x, du, testee.du(x));
|
||
|
||
assert!(abs_diff_eq!(testee.du(x) * testee.dx(u), 1.0, epsilon = ε), "At x={}:\nu'x': {}\n", x, testee.du(x) * testee.dx(u));
|
||
|
||
let d2u = (testee.du(x + δ) - testee.du(x - δ)) / (2. * δ);
|
||
assert!(abs_diff_eq!(d2u, testee.d2u(x), epsilon = ε), "At x={}, u'':\nexpected: {}\nactual: {}\n", x, d2u, testee.d2u(x));
|
||
}
|
||
for u in itertools_num::linspace(-mul * testee.internal, mul * testee.internal, 100) {
|
||
let x = testee.x(u);
|
||
assert!(abs_diff_eq!(x, testee.x(u), epsilon = ε), "At u={}:\nx(u): {}\nu(x(u)): {}\n", u, x, testee.u(x));
|
||
|
||
let dx = (testee.x(u + δ) - testee.x(u - δ)) / (2. * δ);
|
||
assert!(abs_diff_eq!(dx, testee.dx(u), epsilon = ε), "At u={}, x':\nexpected: {}\nactual: {}\n", u, dx, testee.dx(u));
|
||
|
||
assert!(abs_diff_eq!(testee.du(x) * testee.dx(u), 1.0, epsilon = ε), "At u={}:\nu'x': {}\n", u, testee.du(x) * testee.dx(u));
|
||
}
|
||
}
|
||
}
|
||
|
||
trait FloatExt2 {
|
||
fn lerp(&self, t: f32) -> f32;
|
||
fn inverse_lerp(&self, y: f32) -> f32;
|
||
}
|
||
|
||
impl FloatExt2 for (f32, f32) {
|
||
fn lerp(&self, t: f32) -> f32 { f32::lerp(self.0, self.1, t) }
|
||
fn inverse_lerp(&self, y: f32) -> f32 { f32::inverse_lerp(self.0, self.1, y) }
|
||
}
|
||
|
||
impl Metric for Tube {
|
||
fn sqrt_at(&self, pos: Vec2) -> Decomp2 {
|
||
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),
|
||
}
|
||
}
|
||
|
||
fn part_derivs_at(&self, pos: Vec2) -> Tens2 {
|
||
let sx = self.fx().value(pos.x);
|
||
let sy = self.fy().du(pos.y);
|
||
let s = sx + sy - sx * sy;
|
||
let dsx_dx = self.fx().derivative(pos.x);
|
||
let dsy_dy = self.fy().d2u(pos.y);
|
||
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]),
|
||
]
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn test_tube_metric_derivs() {
|
||
struct Approx(Tube);
|
||
impl Metric for Approx {
|
||
fn sqrt_at(&self, pos: Vec2) -> Decomp2 { self.0.sqrt_at(pos) }
|
||
}
|
||
let testee = Tube {
|
||
inner_radius: 30.0,
|
||
outer_radius: 50.0,
|
||
internal_halflength: 100.0,
|
||
external_halflength: 300.0,
|
||
};
|
||
let approx = Approx(testee);
|
||
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 y in itertools_num::linspace(-mul * testee.external_halflength, mul * testee.external_halflength, 100) {
|
||
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 {}:\n explicit: {:?}\n numerical: {:?}\n", pos, computed, reference);
|
||
}
|
||
}
|
||
}
|