axiolid_curve/
quadric_section.rs

1//! Where a quadric cuts a cylinder or cone, read in the ruled surface's own
2//! parameters (ADR 0076).
3//!
4//! A cylinder, elliptical cylinder or cone is RULED: at a fixed angle `u`
5//! its points run along a straight line in `v`,
6//!
7//! ```text
8//! P(u, v) = O + (rx + s v) cos(u) X + (ry + s v) sin(u) Y + v Z,
9//! ```
10//!
11//! so substituting it into any quadric's equation gives, for each `u`, a
12//! quadratic in `v`:
13//!
14//! ```text
15//! a(u) v^2 + b(u) v + c(u) = 0,
16//! ```
17//!
18//! with `a`, `b`, `c` trigonometric polynomials of degree at most two. The
19//! section curve is the graph of one root of that quadratic over the angles
20//! where its discriminant is not negative. [`QuadraticGraph2`] holds that
21//! graph as a pcurve and [`RuledSection3`] the same curve in space. Both are
22//! exact: nothing is sampled or fitted, and a plane's cut ([`Sinusoid2`]'s
23//! case, `a = 0`) is the degenerate member of the family.
24//!
25//! [`Sinusoid2`]: crate::Sinusoid2
26
27use axiolid_core::{Frame3, Point3, Scalar, Vec3};
28
29/// `constant + cos cos(t) + sin sin(t) + cos2 cos(2t) + sin2 sin(2t)`.
30#[derive(Debug, Clone, Copy, PartialEq, Default)]
31pub struct Trig2 {
32    /// Constant term.
33    pub constant: Scalar,
34    /// Coefficient of `cos(t)`.
35    pub cos: Scalar,
36    /// Coefficient of `sin(t)`.
37    pub sin: Scalar,
38    /// Coefficient of `cos(2t)`.
39    pub cos2: Scalar,
40    /// Coefficient of `sin(2t)`.
41    pub sin2: Scalar,
42}
43
44impl Trig2 {
45    /// Value at `t`.
46    #[must_use]
47    pub fn value(&self, t: Scalar) -> Scalar {
48        let (s, c) = t.sin_cos();
49        let (s2, c2) = (2.0 * t).sin_cos();
50        self.constant + self.cos * c + self.sin * s + self.cos2 * c2 + self.sin2 * s2
51    }
52
53    /// First derivative at `t`.
54    #[must_use]
55    pub fn derivative(&self, t: Scalar) -> Scalar {
56        let (s, c) = t.sin_cos();
57        let (s2, c2) = (2.0 * t).sin_cos();
58        -self.cos * s + self.sin * c - 2.0 * self.cos2 * s2 + 2.0 * self.sin2 * c2
59    }
60
61    /// Second derivative at `t`.
62    #[must_use]
63    pub fn second(&self, t: Scalar) -> Scalar {
64        let (s, c) = t.sin_cos();
65        let (s2, c2) = (2.0 * t).sin_cos();
66        -self.cos * c - self.sin * s - 4.0 * self.cos2 * c2 - 4.0 * self.sin2 * s2
67    }
68
69    /// Whether every coefficient is finite.
70    #[must_use]
71    pub fn is_finite(&self) -> bool {
72        self.constant.is_finite()
73            && self.cos.is_finite()
74            && self.sin.is_finite()
75            && self.cos2.is_finite()
76            && self.sin2.is_finite()
77    }
78}
79
80/// Which root of the quadratic a graph follows.
81#[derive(Debug, Clone, Copy, PartialEq, Eq)]
82pub enum Branch {
83    /// `v = (-b + sqrt(b^2 - 4ac)) / 2a`.
84    Plus,
85    /// `v = (-b - sqrt(b^2 - 4ac)) / 2a`.
86    Minus,
87}
88
89impl Branch {
90    /// `+1` or `-1`.
91    #[must_use]
92    pub fn sign(self) -> Scalar {
93        match self {
94            Self::Plus => 1.0,
95            Self::Minus => -1.0,
96        }
97    }
98}
99
100/// The graph of one root of `a(t) v^2 + b(t) v + c(t) = 0`, parameterised by
101/// its first coordinate: the point at `t` is `(t, v(t))`.
102///
103/// Where `a` vanishes the quadratic is linear and the branch that stays
104/// finite is `v = -c / b`; the evaluation below is the rationalised form
105/// `v = 2c / (-b - sign sqrt(D))`, which is finite there and free of the
106/// cancellation `-b + sqrt(D)` suffers when `4ac` is small.
107///
108/// The graph is only defined where `D = b^2 - 4ac >= 0`. Which spans those
109/// are is decided when the graph is constructed (exactly, from the
110/// operands' own numbers) and carried by the edge interval that uses it;
111/// evaluating outside them is refused, not extrapolated.
112#[derive(Debug, Clone, Copy, PartialEq)]
113pub struct QuadraticGraph2 {
114    /// Coefficient of `v^2`.
115    pub a: Trig2,
116    /// Coefficient of `v`.
117    pub b: Trig2,
118    /// Constant coefficient.
119    pub c: Trig2,
120    /// Which root.
121    pub branch: Branch,
122}
123
124impl QuadraticGraph2 {
125    /// The discriminant `b^2 - 4ac` at `t`.
126    #[must_use]
127    pub fn discriminant(&self, t: Scalar) -> Scalar {
128        let b = self.b.value(t);
129        b * b - 4.0 * self.a.value(t) * self.c.value(t)
130    }
131
132    /// Height at `t`, or `None` where the root does not exist or diverges.
133    #[must_use]
134    pub fn height(&self, t: Scalar) -> Option<Scalar> {
135        let (a, b, c) = (self.a.value(t), self.b.value(t), self.c.value(t));
136        let d = b * b - 4.0 * a * c;
137        // A discriminant a few ulps below zero at a branch end is rounding
138        // of an exact zero; the caller's span says the point is on the
139        // curve, so it is read as zero rather than refused. At the end both
140        // terms are themselves near zero, so the bound is taken from the
141        // coefficients' sizes, not from the values there.
142        let size =
143            |t: &Trig2| t.constant.abs() + t.cos.abs() + t.sin.abs() + t.cos2.abs() + t.sin2.abs();
144        let scale = size(&self.b).powi(2) + 4.0 * size(&self.a) * size(&self.c);
145        if d < -1e-12 * scale {
146            return None;
147        }
148        let root = d.max(0.0).sqrt();
149        let denominator = -b - self.branch.sign() * root;
150        if denominator != 0.0 {
151            let v = 2.0 * c / denominator;
152            if v.is_finite() {
153                return Some(v);
154            }
155        }
156        // `-b - sign sqrt(D)` vanishes only where this branch's root is the
157        // other formula's: `(-b + sign sqrt(D)) / 2a`.
158        if a != 0.0 {
159            let v = (-b + self.branch.sign() * root) / (2.0 * a);
160            return v.is_finite().then_some(v);
161        }
162        None
163    }
164
165    /// `dv/dt` at `t`, by differentiating the quadratic implicitly; `None`
166    /// where the height is undefined or the tangent is vertical (a branch
167    /// end, where `2 a v + b = 0`).
168    #[must_use]
169    pub fn slope(&self, t: Scalar) -> Option<Scalar> {
170        let v = self.height(t)?;
171        let f_t = self.a.derivative(t) * v * v + self.b.derivative(t) * v + self.c.derivative(t);
172        let f_v = 2.0 * self.a.value(t) * v + self.b.value(t);
173        let slope = -f_t / f_v;
174        slope.is_finite().then_some(slope)
175    }
176
177    /// `d2v/dt2` at `t`.
178    #[must_use]
179    pub fn bend(&self, t: Scalar) -> Option<Scalar> {
180        let v = self.height(t)?;
181        let p = self.slope(t)?;
182        let f_v = 2.0 * self.a.value(t) * v + self.b.value(t);
183        let f_tt = self.a.second(t) * v * v + self.b.second(t) * v + self.c.second(t);
184        let f_tv = 2.0 * self.a.derivative(t) * v + self.b.derivative(t);
185        let f_vv = 2.0 * self.a.value(t);
186        let bend = -(f_tt + 2.0 * f_tv * p + f_vv * p * p) / f_v;
187        bend.is_finite().then_some(bend)
188    }
189
190    /// Whether every coefficient is finite.
191    #[must_use]
192    pub fn is_finite(&self) -> bool {
193        self.a.is_finite() && self.b.is_finite() && self.c.is_finite()
194    }
195}
196
197/// A ruled carrier surface, as the curve needs it: a cylinder
198/// (`rx = ry`, `slope = 0`), an elliptical cylinder (`slope = 0`) or a cone
199/// (`rx = ry`, `slope = tan(semi-angle)`), in the same parameterisation as
200/// the matching `axiolid_surface` family.
201#[derive(Debug, Clone, Copy, PartialEq)]
202pub struct RuledCarrier {
203    /// Local frame; `z` is the axis.
204    pub frame: Frame3,
205    /// Radius along local `x` at `v = 0`.
206    pub x_radius: Scalar,
207    /// Radius along local `y` at `v = 0`.
208    pub y_radius: Scalar,
209    /// How fast both radii grow with `v`.
210    pub slope: Scalar,
211}
212
213impl RuledCarrier {
214    /// The carrier's point at `(u, v)`.
215    #[must_use]
216    pub fn point(&self, u: Scalar, v: Scalar) -> Point3 {
217        let (s, c) = u.sin_cos();
218        self.frame.origin
219            + self.frame.x * ((self.x_radius + self.slope * v) * c)
220            + self.frame.y * ((self.y_radius + self.slope * v) * s)
221            + self.frame.z * v
222    }
223
224    /// `(dP/du, dP/dv)` at `(u, v)`.
225    #[must_use]
226    pub fn partials(&self, u: Scalar, v: Scalar) -> (Vec3, Vec3) {
227        let (s, c) = u.sin_cos();
228        (
229            self.frame.x * (-(self.x_radius + self.slope * v) * s)
230                + self.frame.y * ((self.y_radius + self.slope * v) * c),
231            self.frame.x * (self.slope * c) + self.frame.y * (self.slope * s) + self.frame.z,
232        )
233    }
234
235    /// Whether every number is finite.
236    #[must_use]
237    pub fn is_finite(&self) -> bool {
238        self.frame.origin.is_finite()
239            && self.frame.x.is_finite()
240            && self.frame.y.is_finite()
241            && self.frame.z.is_finite()
242            && self.x_radius.is_finite()
243            && self.y_radius.is_finite()
244            && self.slope.is_finite()
245    }
246}
247
248/// A quadric's section of a ruled carrier, in space: the point at `t` is the
249/// carrier's point at `(t, graph.height(t))`.
250#[derive(Debug, Clone, Copy, PartialEq)]
251pub struct RuledSection3 {
252    /// The ruled surface the curve lies on.
253    pub carrier: RuledCarrier,
254    /// The curve in the carrier's parameters.
255    pub graph: QuadraticGraph2,
256}
257
258impl RuledSection3 {
259    /// Point at `t`, or `None` outside the graph's spans.
260    #[must_use]
261    pub fn point(&self, t: Scalar) -> Option<Point3> {
262        Some(self.carrier.point(t, self.graph.height(t)?))
263    }
264
265    /// Tangent at `t`: `P_u + P_v dv/dt`.
266    #[must_use]
267    pub fn tangent(&self, t: Scalar) -> Option<Vec3> {
268        let v = self.graph.height(t)?;
269        let (pu, pv) = self.carrier.partials(t, v);
270        Some(pu + pv * self.graph.slope(t)?)
271    }
272
273    /// Second derivative at `t`: `P_uu + 2 P_uv v' + P_v v''` (the carrier
274    /// is linear in `v`, so `P_vv = 0`).
275    #[must_use]
276    pub fn bend(&self, t: Scalar) -> Option<Vec3> {
277        let v = self.graph.height(t)?;
278        let slope = self.graph.slope(t)?;
279        let bend = self.graph.bend(t)?;
280        let (s, c) = t.sin_cos();
281        let k = &self.carrier;
282        let p_uu = k.frame.x * (-(k.x_radius + k.slope * v) * c)
283            + k.frame.y * (-(k.y_radius + k.slope * v) * s);
284        let p_uv = k.frame.x * (-k.slope * s) + k.frame.y * (k.slope * c);
285        let (_, p_v) = k.partials(t, v);
286        Some(p_uu + p_uv * (2.0 * slope) + p_v * bend)
287    }
288}
axiolid_curve/quadric_section.rs

axiolid_curve/
quadric_section.rs
