Resources Numerical accuracy Calculation policy

Transition curves and sampled alignment geometry

Evaluates the six kinds of transition curve (clothoid, cubic parabola, Bloss, Helmert, sine half-wave, sine), finds the nearest point on one, and cuts curves into straight pieces for linework.

How the precision is kept

64-bit throughout, but the accuracy is set by the method: a 24-term power series for the clothoid, composite 8-point Gauss-Legendre quadrature for the others, Newton steps with a relative tolerance for the cubic parabola, and bisection for the nearest point. Linework from curves is a chain of chords limited by turn angle.

Calculations under this policy

12 calculations in 1 source file. The code is read from the released app when this page loads. It is the calculation people run, not a copy kept for this page.

alignment.js CivilCraft 0.8.1 64-bit 12 calculations Reviewed

Tests bin/test-civilcraft-spirals.mjs, bin/test-civilcraft-rail.mjs

Clothoid point

64-bit

The offsets along and across the tangent, and the turn so far, at a distance along a clothoid.

θ = l² ÷ (2·R·Ls); x = l·Σ (−1)ⁿ θ²ⁿ ÷ ((2n)!·(4n+1)); y = l·Σ (−1)ⁿ θ²ⁿ⁺¹ ÷ ((2n+1)!·(4n+3))

Always 24 terms, with no convergence test. The terms alternate in sign, so for a large turn angle big terms cancel; the comment gives about 2 radians as the range where the sum is exact to double precision.

  • alignment.js
  • lines 39 to 49
  • CivilCraft 0.8.1
/** A clothoid from zero curvature: at distance l along a spiral of length Ls that ends at radius R, the point (x along the starting
 *  tangent, y toward the turn) and the turn so far th = l² / (2·R·Ls). The power series is exact to double precision for th up to ~2 rad. */
export function clothoid(l, R, Ls) {
  const th = (l * l) / (2 * R * Ls), t2 = th * th;
  let x = 0, y = 0, fx = 1, fy = th;
  for (let n = 0; n < 24; n++) {
    x += fx / (4 * n + 1); y += fy / (4 * n + 3);
    fx *= -t2 / ((2 * n + 1) * (2 * n + 2)); fy *= -t2 / ((2 * n + 2) * (2 * n + 3));
  }
  return { x: l * x, y: l * y, th };
}

Gauss-Legendre quadrature

64-bit

The integral the other transition kinds and the cubic parabola's length are computed with.

∫ f ≈ (h ÷ 2)·Σ panels Σ₈ wₖ·f(mid + (h ÷ 2)·xₖ)

Eight nodes and weights per panel, written to 16 digits. All terms are added into one running sum in order. The number of panels is fixed by the caller, not adapted to an error estimate.

  • alignment.js
  • lines 70 to 77
  • CivilCraft 0.8.1
const GX = [-0.9602898564975363, -0.7966664774136267, -0.5255324099163290, -0.1834346424956498, 0.1834346424956498, 0.5255324099163290, 0.7966664774136267, 0.9602898564975363];
const GW = [0.1012285362903763, 0.2223810344533745, 0.3137066458778873, 0.3626837833783620, 0.3626837833783620, 0.3137066458778873, 0.2223810344533745, 0.1012285362903763];
function gauss(f, a, b, panels = 8) {
  let s = 0;
  const h = (b - a) / panels;
  for (let p = 0; p < panels; p++) { const m = a + h * (p + 0.5); for (let k = 0; k < 8; k++) s += GW[k] * f(m + (h / 2) * GX[k]); }
  return (s * h) / 2;
}

Curvature shapes

64-bit

The curvature law of each transition kind and its integral, which gives the turn angle.

κ = f(u) ÷ R; θ = Ls·F(u) ÷ R; u = l ÷ Ls

For the sine half-wave and the sine transition, F subtracts two terms that cancel to leading order near the start of the curve (u/2 − sin(πu)/2π, and u²/2 + (cos(2πu) − 1)/4π²), so the small turn angle there carries fewer digits. The absolute error stays at the rounding level of the terms.

  • alignment.js
  • lines 62 to 67
  • CivilCraft 0.8.1
const SHAPES = {
  bloss: { f: (u) => u * u * (3 - 2 * u), F: (u) => u ** 3 - u ** 4 / 2 },
  helmert: { f: (u) => (u <= 0.5 ? 2 * u * u : 1 - 2 * (1 - u) ** 2), F: (u) => (u <= 0.5 ? (2 * u ** 3) / 3 : u - 0.5 + (2 / 3) * (1 - u) ** 3) },
  sinehalf: { f: (u) => (1 - Math.cos(Math.PI * u)) / 2, F: (u) => u / 2 - Math.sin(Math.PI * u) / (2 * Math.PI) },
  sine: { f: (u) => u - Math.sin(2 * Math.PI * u) / (2 * Math.PI), F: (u) => (u * u) / 2 + (Math.cos(2 * Math.PI * u) - 1) / (4 * Math.PI * Math.PI) },
};

Length along a cubic parabola

64-bit

The true length of a cubic parabola transition from its start to a distance along the tangent.

s(x) = ∫₀ˣ √(1 + (v² ÷ (2·R·Ls))²) dv

Four panels of 8-point Gauss-Legendre, so 32 evaluations whatever the length. The integrand uses hypot, which does not overflow or lose the small term.

  • alignment.js
  • lines 78 to 79
  • CivilCraft 0.8.1
/** The cubic parabola's length along the curve from x = 0 to x (y' = x² / 2RLs). */
const cubicArc = (x, R, Ls) => gauss((v) => Math.hypot(1, (v * v) / (2 * R * Ls)), 0, x, 4);

Cubic parabola parameter from a length

64-bit

The tangent length Ls of the cubic parabola whose curve has a given length.

Newton steps on the length residual, stopping at 1e-12 × max(1, length) or after 30 steps. The slope used is the integrand at the end, which is not the full derivative because Ls changes too; the stop test is on the residual itself, so this affects the number of steps and not the answer when it converges. There is no report when 30 steps are used up.

  • alignment.js
  • lines 82 to 87
  • CivilCraft 0.8.1
/** The cubic parabola's Ls (length along the tangent) whose curve is `len` long. */
export function cubicParamFor(R, len) {
  let L = len;
  for (let i = 0; i < 30; i++) { const d = cubicArc(L, R, L) - len; if (Math.abs(d) < 1e-12 * Math.max(1, len)) break; L -= d / Math.hypot(1, L / (2 * R)); }
  return L;
}

Point on a transition of any kind

64-bit

The offsets and turn at a distance along a transition, for all six kinds.

x = ∫₀ˡ cos θ(s) ds; y = ∫₀ˡ sin θ(s) ds

Cubic parabola: Newton on arc length to 1e-13 × max(1, l), at most 30 steps, then y = x³ ÷ (6·R·Ls) and θ = atan(x² ÷ (2·R·Ls)). Other kinds: Gauss-Legendre with max(4, ⌈16·l ÷ Ls⌉) panels, so x and y come from two separate integrals of 32 or more points each.

  • alignment.js
  • lines 88 to 102
  • CivilCraft 0.8.1
/** A transition of `kind` from zero curvature: at distance l along it (its true length), the point (x along the starting tangent,
 *  y toward the turn) and the turn so far th. Clothoid: the exact series; cubic: on y = x³ / 6RLs; the rest by integration. */
export function transitionAt(kind, l, R, Ls) {
  if (kind === 'cubic') {
    // the x whose arc length is l (Newton on s(x) − l, ds/dx = √(1 + y'²))
    let x = l;
    for (let i = 0; i < 30; i++) { const d = cubicArc(x, R, Ls) - l; if (Math.abs(d) < 1e-13 * Math.max(1, l)) break; x -= d / Math.hypot(1, (x * x) / (2 * R * Ls)); }
    return { x, y: (x ** 3) / (6 * R * Ls), th: Math.atan((x * x) / (2 * R * Ls)) };
  }
  if (!SHAPES[kind]) return clothoid(l, R, Ls);
  const F = SHAPES[kind].F, th = (s) => (Ls * F(Math.min(1, Math.max(0, s / Ls)))) / R;
  if (!(l > 0)) return { x: 0, y: 0, th: 0 };
  const n = Math.max(4, Math.ceil(16 * (l / Ls)));
  return { x: gauss((s) => Math.cos(th(s)), 0, l, n), y: gauss((s) => Math.sin(th(s)), 0, l, n), th: th(l) };
}

Transition constants

64-bit

The end offsets Xc and Yc, the shift p and the tangent offset k of a transition, and its true length.

p = Yc − R·(1 − cos θ); k = Xc − R·sin θ

1 − cos θ is formed by subtraction. For a long radius and a short transition θ is small, and the shift p then carries fewer digits than Xc and Yc; the error is about R times the 64-bit rounding of the cosine.

  • alignment.js
  • lines 50 to 56
  • CivilCraft 0.8.1
/** The spiral's constants for radius R and length Ls: θs, Xc, Yc (the end's offsets), the shift p and the offset k (and `len`,
 *  its true length: Ls except for the cubic parabola). */
export function spiralConstants(R, Ls, kind = 'clothoid') {
  if (!(Ls > 0) || !(R > 0)) return { th: 0, Xc: 0, Yc: 0, p: 0, k: 0, len: 0 };
  const len = transitionLength(kind, R, Ls), c = transitionAt(kind, len, R, Ls);
  return { th: c.th, Xc: c.x, Yc: c.y, p: c.y - R * (1 - Math.cos(c.th)), k: c.x - R * Math.sin(c.th), len };
}

Transition element and its long and short tangents

64-bit

Places a transition on the alignment and works out its own point of intersection.

long tangent = Xc − Yc ÷ tan θ; short tangent = Yc ÷ sin θ

The far end of an entry transition is evaluated from its start; an exit transition is evaluated backward from its end. The long tangent is Xc minus Yc ÷ tan θ, about Ls minus a third of Ls, so little is lost in the subtraction. Both divisions are by functions of θ, which is above zero for any transition that has a length.

  • alignment.js
  • lines 204 to 216
  • CivilCraft 0.8.1
/** A spiral element. 'in': from p0 (tangent dir0) curvature 0 → 1/r; 'out': ends at p1 (tangent dir1), curvature 1/r → 0. */
function spiralElement(io, { p0, dir0, p1, dir1, r, Ls, ccw, sta, piIndex, tr = 'clothoid' }) {
  const sg = ccw ? 1 : -1, K = spiralConstants(r, Ls, tr), e = { type: 'spiral', io, ccw, r, Ls, Lp: Ls, tr, len: K.len, sta, piIndex, theta: K.th, delta: K.th, Xc: K.Xc, Yc: K.Yc, p: K.p, k: K.k };
  if (io === 'in') { e.p0 = p0; e.dir0 = dir0; e.dirIn = dir0; e.dir1 = dir0 + sg * K.th; }
  else { e.p1 = p1; e.dir1 = dir1; e.dir0 = dir1 - sg * K.th; e.dirIn = e.dir0; }
  const at = (s) => spiralAt(e, s);
  if (io === 'in') e.p1 = at(e.len); else e.p0 = at(0);
  // the spiral's own PI: where its two end tangents meet (tanLong from the tangent point, tanShort to the curve end)
  const tl = K.Xc - K.Yc / Math.tan(K.th), ts = K.Yc / Math.sin(K.th);
  e.tanLong = tl; e.tanShort = ts;
  e.pi = io === 'in' ? { x: e.p0.x + Math.cos(e.dir0) * tl, y: e.p0.y + Math.sin(e.dir0) * tl } : { x: e.p1.x - Math.cos(e.dir1) * tl, y: e.p1.y - Math.sin(e.dir1) * tl };
  return e;
}

Point on a placed transition

64-bit

The project coordinates and direction at a distance along a transition element.

p = p₀ + R(dir₀)·(x, ±y)

One rotation of the local offsets by the start direction and one addition of the anchor point. The local offsets are small and the anchor is a full-size coordinate, so the sum is limited by the coordinate's own rounding.

  • alignment.js
  • lines 217 to 227
  • CivilCraft 0.8.1
/** The point and forward direction on a spiral element at distance s from its start. */
export function spiralAt(e, s) {
  const sg = e.ccw ? 1 : -1;
  if (e.io === 'in') {
    const c = transitionAt(e.tr, s, e.r, e.Lp ?? e.Ls), co = Math.cos(e.dir0), si = Math.sin(e.dir0);
    return { x: e.p0.x + c.x * co - sg * c.y * si, y: e.p0.y + c.x * si + sg * c.y * co, dir: e.dir0 + sg * c.th };
  }
  // 'out': walked backwards from its end it is an 'in' spiral turning the other way
  const l = Math.max(0, e.len - s), c = transitionAt(e.tr, l, e.r, e.Lp ?? e.Ls), db = e.dir1 + Math.PI, co = Math.cos(db), si = Math.sin(db);
  return { x: e.p1.x + c.x * co + sg * c.y * si, y: e.p1.y + c.x * si - sg * c.y * co, dir: e.dir1 - sg * c.th };
}

Nearest point on a transition

64-bit

The station on a transition closest to a given point, for station and offset.

33 evenly spaced samples pick a bracket, then 70 bisections on the sign of the derivative of the squared distance. Seventy halvings shrink the bracket below 64-bit resolution; the loop always runs all 70. If the transition has two close local minima within one sample spacing, only the sampled one is found.

  • alignment.js
  • lines 228 to 240
  • CivilCraft 0.8.1
/** Nearest point on a spiral element to (x, y): { s, x, y, dir }. */
function spiralNearest(e, x, y) {
  const N = 32, P = [];
  let k = 0, best = Infinity;
  for (let i = 0; i <= N; i++) { const q = spiralAt(e, (e.len * i) / N), d = (q.x - x) ** 2 + (q.y - y) ** 2; P.push(q); if (d < best) { best = d; k = i; } }
  const g = (s) => { const q = spiralAt(e, s); return (q.x - x) * Math.cos(q.dir) + (q.y - y) * Math.sin(q.dir); };   // d(distance²)/2ds
  let lo = (e.len * Math.max(0, k - 1)) / N, hi = (e.len * Math.min(N, k + 1)) / N;
  const glo = g(lo), ghi = g(hi);
  let s;
  if (glo >= 0) s = lo; else if (ghi <= 0) s = hi;
  else { for (let i = 0; i < 70; i++) { const m = (lo + hi) / 2; if (g(m) < 0) lo = m; else hi = m; } s = (lo + hi) / 2; }
  return { s, ...spiralAt(e, s) };
}

Alignment as a polyline

64-bit

The alignment as vertices for linework and export: arcs and transitions become chords.

pieces = max(1, ⌈turn ÷ maxTurn⌉, ⌈length ÷ maxSeg⌉); default maxTurn = 0.5°

Vertices lie on the curve; the chords between them do not. With 0.5° of turn per chord the middle of a chord is R·(1 − cos 0.25°), about 0.0000095 R, inside an arc of radius R. Stations of the vertices are start + length × k ÷ n, not a running sum.

  • alignment.js
  • lines 302 to 311
  • CivilCraft 0.8.1
polyline(maxSeg = Infinity, maxTurn = Math.PI / 360) {
  const out = [];
  for (const e of this.elements) {
    if (!out.length) out.push({ x: e.p0.x, y: e.p0.y, sta: e.sta });
    if (e.type === 'line') { out.push({ x: e.p1.x, y: e.p1.y, sta: e.sta + e.len }); continue; }
    const n = Math.max(1, Math.ceil(e.delta / maxTurn), e.type === 'spiral' ? 4 : 1, Number.isFinite(maxSeg) ? Math.ceil(e.len / maxSeg) : 1);
    for (let k = 1; k <= n; k++) { const s = e.sta + (e.len * k) / n; const p = this.pointAt(s); out.push({ x: p.x, y: p.y, sta: s }); }
  }
  return out;
}

Transitions as straight pieces

64-bit

Each transition cut into short lines so that crossings with other linework can be found; lines and arcs are kept exact.

pieces = max(4, ⌈θ ÷ 0.25°⌉)

Each piece keeps the transition's true length share as its length and stores the chord length separately, so stations along the pieces match the transition while the geometry is the chord.

  • alignment.js
  • lines 333 to 349
  • CivilCraft 0.8.1
/** The elements with each spiral cut into straight pieces no more than `maxTurn` of turn apiece (lines and arcs as they are). */
pieces(maxTurn = Math.PI / 720) {
  if (this._pieces && this._pieces.maxTurn === maxTurn) return this._pieces.list;
  const list = [];
  for (const e of this.elements) {
    if (e.type !== 'spiral') { list.push(e); continue; }
    const n = Math.max(4, Math.ceil(e.delta / maxTurn));
    let q = this.pointAt(e.sta);
    for (let k = 1; k <= n; k++) {
      const sta = e.sta + (e.len * k) / n, p = this.pointAt(sta), len = Math.hypot(p.x - q.x, p.y - q.y);
      list.push({ type: 'line', p0: { x: q.x, y: q.y }, p1: { x: p.x, y: p.y }, dir: Math.atan2(p.y - q.y, p.x - q.x), sta: sta - e.len / n, len: e.len / n, chordLen: len });
      q = p;
    }
  }
  this._pieces = { maxTurn, list };
  return list;
}

In detail

Method
Clothoid: the power series for the Fresnel integrals, 24 terms. Bloss, Helmert, sine half-wave and sine: x and y are the integrals of cos θ and sin θ by 8-point Gauss-Legendre in at least 4 panels. Cubic parabola: y = x³ ÷ (6·R·Ls), with x found from arc length by Newton's method.
Number types
64-bit numbers. The eight Gauss nodes and weights are written to 16 digits.
Tolerances
Newton for the cubic parabola stops at 1e-13 × max(1, length) (point) or 1e-12 × max(1, length) (parameter from length), at most 30 steps. Nearest point: 33 samples, then 70 bisections. Chords: at most 0.5° of turn each in polyline, 0.25° in pieces.
Limits
The clothoid series has no convergence test; the code comment gives about 2 radians of turn as its range. The file header gives about 1e-12 of Ls for the integrated kinds; that figure is the comment's and was not re-derived in this review.

Approach: 64-bit, stated limit

The arithmetic is 64-bit, but the method or the data sets the accuracy: a solver tolerance, a mesh, a step size, a datum shift, a table. The limit is stated where the result is shown.

Every policy that takes this approach