axiolid_curve/elevation.rs
1//! Elevation laws and the composition of a planar curve with one.
2//!
3//! A road or rail centreline is authored as two independent laws: a planar
4//! layout carrying the transition spirals, and a vertical profile giving
5//! height as a function of distance along that layout. The 3D centreline is
6//! their composition. This module holds the vertical half and the pairing;
7//! the planar half is an ordinary [`Curve2`].
8//!
9//! # Height is a function of PLAN distance, not 3D arc length
10//!
11//! The vertical profile a surveyor writes is indexed by chainage measured on
12//! the horizontal layout, so `height_at` takes plan distance. The two differ
13//! whenever the grade is non-zero, because `ds3 = sqrt(1 + g^2) ds_plan`.
14//! Over a 120 m curve at a 2% entry grade the 3D length exceeds the plan
15//! length by 12.5 mm, which reads the profile 0.35 mm off if the distances
16//! are confused. That is small but it is a systematic misreading, not noise,
17//! and it compounds along a chain of segments.
18//!
19//! Naming the convention here means a consumer never has to guess which
20//! distance a law is written in.
21//!
22//! # Closed forms here, quadrature in an evaluator
23//!
24//! A polynomial and a circular arc in the `(d, z)` plane have closed-form
25//! heights and grades, so [`ElevationLaw::height_at`] and
26//! [`ElevationLaw::grade_at`] answer them here. A profile given by its
27//! curvature against its own arc length ([`ElevationLaw::Intrinsic`], the
28//! clothoid between two grades) does not: its height is the Fresnel-type
29//! integral the [`CurvatureLaw`] module refuses to compute, followed by an
30//! inversion of plan distance against arc length. This crate stores that
31//! law exactly and reports `None` for it; an evaluator that can state its
32//! quadrature tolerance reads it (`axiolid-evaluate`'s `elevation` module).
33
34use axiolid_core::Scalar;
35
36use crate::{CurvatureLaw, Curve2};
37
38/// Height as a function of distance along the plan.
39///
40/// Mirrors [`CurvatureLaw`] in shape so the two halves of
41/// an alignment read the same way, but stays a separate type: a curvature law
42/// is a property of a planar curve and an elevation law is not, and sharing one
43/// enum would make a meaningless pairing representable.
44///
45/// Dirty imported data stays representable, as everywhere else in this crate.
46/// Mismatched piece lists report `None` rather than guessing.
47#[non_exhaustive]
48#[derive(Debug, Clone, PartialEq)]
49#[cfg_attr(
50 feature = "serde",
51 derive(serde::Serialize, serde::Deserialize),
52 serde(deny_unknown_fields)
53)]
54pub enum ElevationLaw {
55 /// `z(d) = coefficients[0] + coefficients[1] * d + coefficients[2] * d^2 + ...`
56 ///
57 /// Degree 1 is a constant gradient, degree 2 the parabolic vertical curve
58 /// used to join two grades. Those are the two vertical segment kinds that
59 /// carry most alignment data, and both are exact here.
60 ///
61 /// An empty coefficient list is the zero polynomial: height zero.
62 Polynomial {
63 /// Coefficients in ascending powers of plan distance.
64 coefficients: Vec<Scalar>,
65 },
66 /// Pieces laid end to end along plan distance, each with its own law.
67 ///
68 /// `breaks` holds the INTERIOR seam positions measured from the start, so
69 /// `laws.len() == breaks.len() + 1` and piece `i` spans
70 /// `breaks[i - 1] .. breaks[i]`.
71 ///
72 /// Each piece's law is written in its OWN distance, restarting at zero at
73 /// its seam, so moving a piece never rewrites its coefficients. This
74 /// matches [`CurvatureLaw::Piecewise`]
75 /// deliberately: a vertical profile is authored as a run of segments and
76 /// the seams are observable data.
77 Piecewise {
78 /// Interior seam positions in plan distance, ascending.
79 breaks: Vec<Scalar>,
80 /// One law per piece; `laws.len() == breaks.len() + 1`.
81 laws: Vec<ElevationLaw>,
82 },
83 /// A circular arc in the `(d, z)` plane: the grade angle changes at a
84 /// constant rate along the profile's own arc length.
85 ///
86 /// With `t0 = atan(grade)` and the signed `radius` `R` (positive is a
87 /// sag, turning counter-clockwise in `(d, z)`; negative is a crest):
88 ///
89 /// - `sin t(d) = sin t0 + d / R`
90 /// - `z(d) = height + R (cos t0 - cos t(d))`
91 /// - `grade(d) = tan t(d)`
92 ///
93 /// This is the circle itself, not the parabola `d^2 / (2R)` that
94 /// approximates it near the vertex; the two differ by `O(d^4 / R^3)`.
95 /// The arc has heights only while `|sin t0 + d / R| < 1`: past that the
96 /// circle turns vertical, and `height_at` and `grade_at` report `None`.
97 ///
98 /// Evaluated without cancellation: `R (cos t0 - cos t(d))` is computed
99 /// as `d (sin t(d) + sin t0) / (cos t0 + cos t(d))`, which removes the
100 /// difference of nearly equal cosines that `d / R -> 0` would otherwise
101 /// multiply by a large `R`.
102 CircularArc {
103 /// Height at the piece start.
104 height: Scalar,
105 /// Grade `dz/dd` at the piece start.
106 grade: Scalar,
107 /// Signed radius: positive sag, negative crest. Finite and non-zero.
108 radius: Scalar,
109 },
110 /// A profile given by its curvature against its OWN arc length.
111 ///
112 /// The vertical counterpart of [`Curve2::Intrinsic`]: a planar curve in
113 /// `(d, z)` starting at `(0, height)` in the direction `atan(grade)`,
114 /// whose curvature is `curvature(s)` with `s` the arc length of the
115 /// profile (not plan distance), positive counter-clockwise in `(d, z)`
116 /// -- a sag. A linear law is the clothoid joining two grades; a constant
117 /// law is the circular arc that [`Self::CircularArc`] holds in closed
118 /// form.
119 ///
120 /// Height at a plan distance has no closed form (see the module
121 /// documentation), so `height_at` and `grade_at` report `None`, and an
122 /// evaluator integrates the curve and inverts `d(s)` to a stated
123 /// tolerance. Only `d >= 0` is defined, since the law is integrated
124 /// forward from its start, and only while `d(s)` keeps increasing.
125 Intrinsic {
126 /// Height at the piece start.
127 height: Scalar,
128 /// Grade `dz/dd` at the piece start.
129 grade: Scalar,
130 /// Curvature against the profile's own arc length from its start.
131 curvature: CurvatureLaw,
132 },
133}
134
135impl ElevationLaw {
136 /// A constant height.
137 #[must_use]
138 pub fn level(height: Scalar) -> Self {
139 Self::Polynomial {
140 coefficients: vec![height],
141 }
142 }
143
144 /// A constant gradient: `z(d) = height + grade * d`.
145 #[must_use]
146 pub fn constant_grade(height: Scalar, grade: Scalar) -> Self {
147 Self::Polynomial {
148 coefficients: vec![height, grade],
149 }
150 }
151
152 /// A parabolic vertical curve joining `entry_grade` to `exit_grade` over
153 /// `length`.
154 ///
155 /// The rate of change of grade is `(exit - entry) / length`, so
156 /// `z(d) = height + entry * d + (exit - entry) / (2 * length) * d^2`.
157 /// A non-positive length is storable; naming it is a validator's job.
158 #[must_use]
159 pub fn parabolic(
160 height: Scalar,
161 entry_grade: Scalar,
162 exit_grade: Scalar,
163 length: Scalar,
164 ) -> Self {
165 Self::Polynomial {
166 coefficients: vec![
167 height,
168 entry_grade,
169 (exit_grade - entry_grade) / (2.0 * length),
170 ],
171 }
172 }
173
174 /// A circular vertical arc from `height` at `grade`, of signed `radius`
175 /// (positive sag, negative crest). See [`Self::CircularArc`].
176 #[must_use]
177 pub const fn circular_arc(height: Scalar, grade: Scalar, radius: Scalar) -> Self {
178 Self::CircularArc {
179 height,
180 grade,
181 radius,
182 }
183 }
184
185 /// A profile from `height` at `grade` whose curvature runs along its own
186 /// arc length. See [`Self::Intrinsic`].
187 #[must_use]
188 pub const fn intrinsic(height: Scalar, grade: Scalar, curvature: CurvatureLaw) -> Self {
189 Self::Intrinsic {
190 height,
191 grade,
192 curvature,
193 }
194 }
195
196 /// Whether the piece lists agree and the seams ascend, and every stored
197 /// number is finite (an arc's radius also non-zero).
198 #[must_use]
199 pub fn is_well_formed(&self) -> bool {
200 match self {
201 Self::Polynomial { coefficients } => coefficients.iter().all(|c| c.is_finite()),
202 Self::CircularArc {
203 height,
204 grade,
205 radius,
206 } => height.is_finite() && grade.is_finite() && radius.is_finite() && *radius != 0.0,
207 Self::Intrinsic {
208 height,
209 grade,
210 curvature,
211 } => height.is_finite() && grade.is_finite() && curvature.is_well_formed(),
212 Self::Piecewise { breaks, laws } => {
213 laws.len() == breaks.len() + 1
214 && breaks.iter().all(|b| b.is_finite())
215 && breaks.windows(2).all(|w| w[0] <= w[1])
216 && laws.iter().all(Self::is_well_formed)
217 }
218 }
219 }
220
221 /// Height at `distance` along the plan.
222 ///
223 /// Returns `None` when the law is malformed or the distance is not finite,
224 /// rather than extrapolating off a piece that does not exist.
225 #[must_use]
226 pub fn height_at(&self, distance: Scalar) -> Option<Scalar> {
227 if !distance.is_finite() {
228 return None;
229 }
230 match self {
231 Self::Polynomial { coefficients } => {
232 Some(horner(coefficients, distance)).filter(|z| z.is_finite())
233 }
234 Self::Piecewise { breaks, laws } => {
235 let (law, local) = piece_at(breaks, laws, distance)?;
236 law.height_at(local)
237 }
238 Self::CircularArc {
239 height,
240 grade,
241 radius,
242 } => {
243 let arc = ArcState::at(*grade, *radius, distance)?;
244 // R (cos t0 - cos t) = R (sin^2 t - sin^2 t0) / (cos t0 + cos t)
245 // and R (sin t - sin t0) = d, so the large R cancels exactly.
246 Some(height + distance * (arc.sin + arc.sin0) / (arc.cos0 + arc.cos))
247 .filter(|z| z.is_finite())
248 }
249 // No closed form: an evaluator integrates it (module docs).
250 Self::Intrinsic { .. } => None,
251 }
252 }
253
254 /// Grade -- `dz/dd` -- at `distance` along the plan.
255 ///
256 /// This is the slope the 3D tangent needs, and it is why the polynomial is
257 /// differentiated exactly rather than differenced.
258 #[must_use]
259 pub fn grade_at(&self, distance: Scalar) -> Option<Scalar> {
260 if !distance.is_finite() {
261 return None;
262 }
263 match self {
264 Self::Polynomial { coefficients } => {
265 let derivative: Vec<Scalar> = coefficients
266 .iter()
267 .enumerate()
268 .skip(1)
269 .map(|(power, c)| c * power as Scalar)
270 .collect();
271 Some(horner(&derivative, distance)).filter(|g| g.is_finite())
272 }
273 Self::Piecewise { breaks, laws } => {
274 let (law, local) = piece_at(breaks, laws, distance)?;
275 law.grade_at(local)
276 }
277 Self::CircularArc { grade, radius, .. } => {
278 let arc = ArcState::at(*grade, *radius, distance)?;
279 Some(arc.sin / arc.cos).filter(|g| g.is_finite())
280 }
281 // No closed form: an evaluator integrates it (module docs).
282 Self::Intrinsic { .. } => None,
283 }
284 }
285
286 /// The innermost piece covering `distance`, and that distance rebased
287 /// to the piece's own start.
288 ///
289 /// Descends through nested [`Self::Piecewise`] laws, so the returned law
290 /// is never piecewise; any other law is its own piece at the unchanged
291 /// distance. Seams belong to the piece that starts there, as in
292 /// `height_at` ([`SeamSide::Outgoing`](crate::SeamSide::Outgoing)). `None` when the distance is not finite or a piece list
293 /// is malformed.
294 ///
295 /// This is what an evaluator needs to read a law it evaluates itself
296 /// (an [`Self::Intrinsic`] piece) inside a composed profile.
297 #[must_use]
298 pub fn piece_at(&self, distance: Scalar) -> Option<(&Self, Scalar)> {
299 if !distance.is_finite() {
300 return None;
301 }
302 match self {
303 Self::Piecewise { breaks, laws } => {
304 let (law, local) = piece_at(breaks, laws, distance)?;
305 law.piece_at(local)
306 }
307 _ => Some((self, distance)),
308 }
309 }
310}
311
312/// The trigonometry of a circular vertical arc at one plan distance.
313struct ArcState {
314 /// `sin t0` of the start direction.
315 sin0: Scalar,
316 /// `cos t0` of the start direction.
317 cos0: Scalar,
318 /// `sin t(d)` at the distance.
319 sin: Scalar,
320 /// `cos t(d)` at the distance, positive inside the domain.
321 cos: Scalar,
322}
323
324impl ArcState {
325 /// `None` outside the arc's domain `|sin t0 + d / R| < 1`, or for a
326 /// non-finite grade or a zero or non-finite radius.
327 fn at(grade: Scalar, radius: Scalar, distance: Scalar) -> Option<Self> {
328 if !(grade.is_finite() && radius.is_finite() && radius != 0.0) {
329 return None;
330 }
331 // sin and cos of atan(grade) without forming the angle: no rounding
332 // through atan and back, and exact at grade zero.
333 let norm = grade.hypot(1.0);
334 let sin0 = grade / norm;
335 let cos0 = 1.0 / norm;
336 let sin = sin0 + distance / radius;
337 if sin.is_nan() || sin.abs() >= 1.0 {
338 return None;
339 }
340 // (1 - s)(1 + s) keeps its relative accuracy as |s| -> 1, where
341 // 1 - s^2 would not.
342 let cos = ((1.0 - sin) * (1.0 + sin)).sqrt();
343 if cos <= 0.0 {
344 return None;
345 }
346 Some(Self {
347 sin0,
348 cos0,
349 sin,
350 cos,
351 })
352 }
353}
354
355/// Evaluate ascending-power coefficients at `x`.
356fn horner(coefficients: &[Scalar], x: Scalar) -> Scalar {
357 coefficients
358 .iter()
359 .rev()
360 .fold(0.0, |accumulated, c| accumulated * x + c)
361}
362
363/// The piece covering `distance`, and that distance rebased to the piece start.
364///
365/// Pieces are half-open so a seam belongs to the piece that starts there; the
366/// final piece is closed at its far end so the curve's last point evaluates.
367fn piece_at<'a>(
368 breaks: &[Scalar],
369 laws: &'a [ElevationLaw],
370 distance: Scalar,
371) -> Option<(&'a ElevationLaw, Scalar)> {
372 if laws.len() != breaks.len() + 1 {
373 return None;
374 }
375 let index = breaks.partition_point(|b| *b <= distance);
376 let start = if index == 0 { 0.0 } else { breaks[index - 1] };
377 laws.get(index).map(|law| (law, distance - start))
378}
379
380/// A planar curve carrying an independent elevation law.
381///
382/// This is the composition an alignment centreline actually is: the plan is
383/// exact -- including the transition spirals a [`Curve2::Intrinsic`] holds --
384/// and the vertical profile is exact, and neither is approximated to pair
385/// them. Evaluation is parameterised by PLAN DISTANCE, which is the parameter
386/// both halves are authored against.
387///
388/// Composition rather than re-encoding: a `Curve3::BSpline` fitted through the
389/// pair would lose both. Storing the two laws keeps each one's own exactness
390/// and lets a consumer recover either half unchanged.
391///
392/// Only a plan whose parameter is arc length (line, circle, intrinsic,
393/// [`Chain2`](crate::Chain2)) can carry a law. A B-spline's parameter is not
394/// a distance, so evaluators refuse an elevated B-spline plan rather than
395/// reinterpret the law (ADR 0060); a B-spline read by arc length is a
396/// chain piece.
397#[derive(Debug, Clone, PartialEq)]
398#[cfg_attr(
399 feature = "serde",
400 derive(serde::Serialize, serde::Deserialize),
401 serde(deny_unknown_fields)
402)]
403pub struct Elevated3 {
404 /// Horizontal layout. Boxed to keep [`Curve3`](crate::Curve3) small.
405 pub plan: Box<Curve2>,
406 /// Height as a function of distance along `plan`.
407 pub elevation: ElevationLaw,
408}
409
410impl Elevated3 {
411 /// Pair a planar curve with an elevation law.
412 #[must_use]
413 pub fn new(plan: Curve2, elevation: ElevationLaw) -> Self {
414 Self {
415 plan: Box::new(plan),
416 elevation,
417 }
418 }
419}