1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125

use super::*;

impl<P> UnitHyperbola<P> {
    /// constructor
    #[inline]
    pub const fn new() -> UnitHyperbola<P> { UnitHyperbola(std::marker::PhantomData) }
}

impl ParametricCurve for UnitHyperbola<Point2> {
    type Point = Point2;
    type Vector = Vector2;
    #[inline]
    fn subs(&self, t: f64) -> Self::Point { Point2::new(f64::cosh(t), f64::sinh(t)) }
    #[inline]
    fn der(&self, t: f64) -> Self::Vector { Vector2::new(f64::sinh(t), f64::cosh(t)) }
    #[inline]
    fn der2(&self, t: f64) -> Self::Vector { Vector2::new(f64::cosh(t), f64::sinh(t)) }
}

impl ParametricCurve for UnitHyperbola<Point3> {
    type Point = Point3;
    type Vector = Vector3;
    #[inline]
    fn subs(&self, t: f64) -> Self::Point { Point3::new(f64::cosh(t), f64::sinh(t), 0.0) }
    #[inline]
    fn der(&self, t: f64) -> Self::Vector { Vector3::new(f64::sinh(t), f64::cosh(t), 0.0) }
    #[inline]
    fn der2(&self, t: f64) -> Self::Vector { Vector3::new(f64::cosh(t), f64::sinh(t), 0.0) }
}

impl<P> ParameterDivision1D for UnitHyperbola<P>
where
    UnitHyperbola<P>: ParametricCurve<Point = P>,
    P: EuclideanSpace<Scalar = f64> + MetricSpace<Metric = f64> + HashGen<f64>,
{
    type Point = P;
    fn parameter_division(&self, range: (f64, f64), tol: f64) -> (Vec<f64>, Vec<P>) {
        algo::curve::parameter_division(self, range, tol)
    }
}

impl SearchNearestParameter<D1> for UnitHyperbola<Point2> {
    type Point = Point2;
    fn search_nearest_parameter<H: Into<SPHint1D>>(
        &self,
        p: Point2,
        _: H,
        _: usize,
    ) -> Option<f64> {
        let a = -p.y;
        let b = (p.y * p.y - p.x * p.x) / 4.0 + 1.0;
        let c = -p.y;
        let d = p.y * p.y / 4.0;
        let y = solver::solve_quartic(a, b, c, d)
            .into_iter()
            .filter_map(|z| match z.im.so_small() {
                true => Some(z.re),
                false => None,
            })
            .min_by(|s, t| {
                p.distance2(self.subs(*s))
                    .partial_cmp(&p.distance2(self.subs(*t)))
                    .unwrap()
            })?;
        Some(f64::asinh(y))
    }
}

impl SearchNearestParameter<D1> for UnitHyperbola<Point3> {
    type Point = Point3;
    fn search_nearest_parameter<H: Into<SPHint1D>>(
        &self,
        p: Point3,
        _: H,
        _trials: usize,
    ) -> Option<f64> {
        UnitHyperbola::<Point2>::new().search_nearest_parameter(
            Point2::new(p.x, p.y),
            None,
            _trials,
        )
    }
}

impl SearchParameter<D1> for UnitHyperbola<Point2> {
    type Point = Point2;
    fn search_parameter<H: Into<SPHint1D>>(&self, p: Point2, _: H, _: usize) -> Option<f64> {
        let t = f64::asinh(p.y);
        match p.near(&self.subs(t)) {
            true => Some(t),
            false => None,
        }
    }
}

impl SearchParameter<D1> for UnitHyperbola<Point3> {
    type Point = Point3;
    fn search_parameter<H: Into<SPHint1D>>(&self, p: Point3, _: H, _: usize) -> Option<f64> {
        let t = f64::asinh(p.y);
        match p.near(&self.subs(t)) {
            true => Some(t),
            false => None,
        }
    }
}

#[test]
fn snp_test() {
    let curve = UnitHyperbola::<Point2>::new();
    let p = curve.subs(2.0);
    let q = p + Vector2::new(-p.x, p.y);
    let t = curve.search_nearest_parameter(q, None, 0).unwrap();
    assert_near!(t, 2.0);
}

#[test]
fn sp_test() {
    let curve = UnitHyperbola::<Point2>::new();
    let t = 100.0 * rand::random::<f64>() - 50.0;
    let p = curve.subs(t);
    assert_near!(curve.search_parameter(p, None, 0).unwrap(), t);

    let q = Point2::new(-1.0, 0.0);
    assert!(curve.search_parameter(q, None, 0).is_none());
}