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

axiolid_curve/
quadric_section.rs
