Resources Numerical accuracy Calculation policy
Solves the magnetic field of a sketch's profiles (air, iron, magnets, coils) and reports the peak flux density and the force and torque on chosen bodies.
Everything is 64-bit: node coordinates, the potential, the band matrix and the per-element fields are Float64Array. The accuracy of the answer is set by the method: a planar model on linear triangles (one flux density per element), a mesh graded by a size field, sketch arcs entered as 48-segment polygons, and a Newton iteration stopped at a relative change of 1e-9. Forces come from the Maxwell stress on a box contour in air, with the field on each edge taken as the mean of its two elements. The CSV export writes 8 significant digits.
18 calculations in 3 source files. 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.
Tests tools/partcraft-test/test-sim08.mjs
The constant μ₀ used for air, linear materials and the Maxwell stress.
μ₀ = 4 × 10⁻⁷ × π
Written as 4e-7 × Math.PI, the exact value under the SI before 2019. Since then μ₀ is a measured quantity whose recommended value differs from this in about the tenth significant digit.
export const MU0 = 4e-7 * Math.PI;Turns a table of B and H points into ν(B²) and its derivative for the solver.
H(B) by a monotone cubic through the points; beyond the last point H = H_s + (B − B_s) ÷ μ₀; ν = H ÷ B; dν/dB² = (H′ ÷ B − H ÷ B²) ÷ (2B)
Below B = 1e-9 the first segment's slope is used for ν and 0 for its derivative, which avoids 0 ÷ 0. Where the curve is a straight line through the origin, H′ ÷ B and H ÷ B² are nearly equal and their difference is rounding noise; that value only steers the Newton iteration and does not enter the converged field.
/** A measured curve: points [[B (T), H (A/m)], …] with B and H rising; ν(B²) = H(B)/B through a monotone cubic (Fritsch–Carlson). */
export function table(points) {
const P = [[0, 0], ...points.filter((p) => p[0] > 0)].sort((a, b) => a[0] - b[0]);
for (let i = 1; i < P.length; i++) if (!(P[i][0] > P[i - 1][0] && P[i][1] > P[i - 1][1])) throw new Error('B–H table: B and H must both rise'); // i18n-ignore: internal check
const n = P.length, dx = [], sl = [];
for (let i = 0; i < n - 1; i++) { dx.push(P[i + 1][0] - P[i][0]); sl.push((P[i + 1][1] - P[i][1]) / dx[i]); }
const m = new Array(n);
m[0] = sl[0]; m[n - 1] = sl[n - 2];
for (let i = 1; i < n - 1; i++) m[i] = sl[i - 1] * sl[i] <= 0 ? 0 : 2 * sl[i - 1] * sl[i] / (sl[i - 1] + sl[i]);
for (let i = 0; i < n - 1; i++) { // keep it monotone
if (sl[i] === 0) { m[i] = m[i + 1] = 0; continue; }
const a = m[i] / sl[i], b = m[i + 1] / sl[i], s = a * a + b * b;
if (s > 9) { const t = 3 / Math.sqrt(s); m[i] = t * a * sl[i]; m[i + 1] = t * b * sl[i]; }
}
const Bs = P[n - 1][0], Hs = P[n - 1][1];
// H(B) and dH/dB
const HB = (B) => {
if (B >= Bs) return { h: Hs + (B - Bs) / MU0, d: 1 / MU0 };
let lo = 0, hi = n - 2; while (lo < hi) { const mid = (lo + hi + 1) >> 1; if (P[mid][0] <= B) lo = mid; else hi = mid - 1; }
const h = dx[lo], t = (B - P[lo][0]) / h, t2 = t * t, t3 = t2 * t;
const h00 = 2 * t3 - 3 * t2 + 1, h10 = t3 - 2 * t2 + t, h01 = -2 * t3 + 3 * t2, h11 = t3 - t2;
const val = h00 * P[lo][1] + h10 * h * m[lo] + h01 * P[lo + 1][1] + h11 * h * m[lo + 1];
const d = ((6 * t2 - 6 * t) * P[lo][1] + (3 * t2 - 4 * t + 1) * h * m[lo] + (-6 * t2 + 6 * t) * P[lo + 1][1] + (3 * t2 - 2 * t) * h * m[lo + 1]) / h;
return { h: val, d };
};
const mu0 = MU0, nu0 = P[1][1] / P[1][0]; // the first segment's slope for B → 0
return {
nu: (b2) => { const B = Math.sqrt(b2); return B < 1e-9 ? sl[0] : HB(B).h / B; },
dnu: (b2) => { const B = Math.sqrt(b2); if (B < 1e-9) return 0; const { h, d } = HB(B); return (d / B - h / (B * B)) / (2 * B); },
H: (B) => HB(B).h, mu0, nu0,
};
}Factors the tangent matrix and solves for the correction of the potential.
K = L Lᵀ in band storage; forward then back substitution
Float64Array band, factored in place with inner products subtracted in column order. A diagonal term that is not positive stops the solve with an error instead of taking a square root of it.
function bandSolve(K, n, bw, rhs) {
const w = bw + 1;
const L = K; // in place: lower band, L[i*w + (j − i + bw)]
for (let i = 0; i < n; i++) {
const j0 = Math.max(0, i - bw);
for (let j = j0; j <= i; j++) {
let s = L[i * w + (j - i + bw)];
for (let k = Math.max(j0, j - bw); k < j; k++) s -= L[i * w + (k - i + bw)] * L[j * w + (k - j + bw)];
if (i === j) { if (!(s > 0)) throw new Error('magnetostatic matrix is not positive definite (is some node left floating? fix A on the outer boundary)'); L[i * w + bw] = Math.sqrt(s); } // i18n-ignore: internal check
else L[i * w + (j - i + bw)] = s / L[j * w + bw];
}
}
const y = Float64Array.from(rhs);
for (let i = 0; i < n; i++) { let s = y[i]; for (let k = Math.max(0, i - bw); k < i; k++) s -= L[i * w + (k - i + bw)] * y[k]; y[i] = s / L[i * w + bw]; }
for (let i = n - 1; i >= 0; i--) { let s = y[i]; for (let k = i + 1; k <= Math.min(n - 1, i + bw); k++) s -= L[k * w + (i - k + bw)] * y[k]; y[i] = s / L[i * w + bw]; }
return y;
}The area, center and the three gradient vectors of every triangle.
d = (x₁ − x₀)(y₂ − y₀) − (x₂ − x₀)(y₁ − y₀); area = |d| ÷ 2; ∇Nᵢ = (bᵢ, cᵢ) ÷ d
Differences of node coordinates are formed first, then one 2 × 2 determinant. The division by d is not guarded; a triangle with no area would give infinite gradients.
/** Element geometry: area and the gradient vectors of the three shape functions. */
export function elementGeometry(mesh) {
const m = mesh.tris.length / 3, area = new Float64Array(m), gx = new Float64Array(3 * m), gy = new Float64Array(3 * m), cx = new Float64Array(m), cy = new Float64Array(m);
for (let e = 0; e < m; e++) {
const [a, b, c] = [mesh.tris[3 * e], mesh.tris[3 * e + 1], mesh.tris[3 * e + 2]], X = mesh.nodes;
const x = [X[2 * a], X[2 * b], X[2 * c]], y = [X[2 * a + 1], X[2 * b + 1], X[2 * c + 1]];
const d = (x[1] - x[0]) * (y[2] - y[0]) - (x[2] - x[0]) * (y[1] - y[0]);
area[e] = Math.abs(d) / 2; cx[e] = (x[0] + x[1] + x[2]) / 3; cy[e] = (y[0] + y[1] + y[2]) / 3;
const bI = [y[1] - y[2], y[2] - y[0], y[0] - y[1]], cI = [x[2] - x[1], x[0] - x[2], x[1] - x[0]];
for (let i = 0; i < 3; i++) { gx[3 * e + i] = bI[i] / d; gy[3 * e + i] = cI[i] / d; }
}
return { area, gx, gy, cx, cy };
}Assembles the out-of-balance load and, when asked, the tangent stiffness for the current potential.
R = load − Σ ν A (∇Nᵢ·∇A); K_ik = Σ A (ν ∇Nᵢ·∇N_k + 2ν′ (∇Nᵢ·∇A)(∇N_k·∇A))
Element contributions are added into Float64Array storage in element order. Fixed nodes get a zero residual, a unit diagonal and zeroed couplings.
/** The residual load − ∫ν∇A·∇N (zero at the fixed nodes) and, when asked, the tangent in band storage. */
const assemble = (Av, wantK) => {
const K = wantK ? new Float64Array(n * w) : null, R = Float64Array.from(load);
for (let e = 0; e < m; e++) {
const mat = mats[region[e]], nd = [mesh.tris[3 * e], mesh.tris[3 * e + 1], mesh.tris[3 * e + 2]], ar = G.area[e];
let ax = 0, ay = 0; for (let i = 0; i < 3; i++) { ax += G.gx[3 * e + i] * Av[nd[i]]; ay += G.gy[3 * e + i] * Av[nd[i]]; }
const b2 = ax * ax + ay * ay, nu = mat.nu(b2), dnu = wantK && nonlinear ? mat.dnu(b2) : 0;
for (let i = 0; i < 3; i++) {
const gix = G.gx[3 * e + i], giy = G.gy[3 * e + i];
R[nd[i]] -= nu * ar * (gix * ax + giy * ay);
if (!wantK) continue;
for (let k = 0; k < 3; k++) {
if (nd[k] > nd[i]) continue;
const gkx = G.gx[3 * e + k], gky = G.gy[3 * e + k];
K[nd[i] * w + (nd[k] - nd[i] + bw)] += ar * (nu * (gix * gkx + giy * gky) + 2 * dnu * (gix * ax + giy * ay) * (gkx * ax + gky * ay));
}
}
}
for (const [f] of fixed) {
R[f] = 0;
if (!wantK) continue;
for (let j = Math.max(0, f - bw); j <= Math.min(n - 1, f + bw); j++) {
if (j === f) continue;
const hi = Math.max(f, j), lo = Math.min(f, j); K[hi * w + (lo - hi + bw)] = 0;
}
K[f * w + bw] = 1;
}
return { R, K };
};Repeats solve and update until the potential stops changing.
A ← A + α ΔA, α halved until the residual norm falls; stop when |ΔA| ÷ |A| < tol or |R| < 1e-14 × |load|
A linear problem is solved in one step. Otherwise at most 40 iterations (an option), each with up to 30 halvings of the step; the default tolerance is 1e-9 on the relative change. If no halving lowers the residual the loop ends and the result is flagged as not converged.
let it = 0, change = Infinity, cur = assemble(A, true), rn = nrm(cur.R);
for (; it < (spec.maxIter || 40); it++) {
// the unknown is the correction dA, zero at the fixed nodes (their value is already in A)
const dA = bandSolve(cur.K, n, bw, cur.R);
let alpha = 1, next = null, nn = Infinity;
if (!nonlinear) { for (let i = 0; i < n; i++) A[i] += dA[i]; change = 0; it = 1; break; }
for (let k = 0; k < 30; k++, alpha /= 2) { // damped Newton: halve the step until the residual falls (and stays finite)
const At = Float64Array.from(A, (v, i) => v + alpha * dA[i]);
next = assemble(At, false); nn = nrm(next.R);
if (Number.isFinite(nn) && nn < rn) { A.set(At); break; }
next = null;
}
if (!next) break;
let dn = 0, an = 0; for (let i = 0; i < n; i++) { dn += dA[i] ** 2; an += A[i] * A[i]; }
change = Math.sqrt(dn / (an || 1));
cur = assemble(A, true); rn = nrm(cur.R);
if (change < (spec.tol || 1e-9) || rn < 1e-14 * (nrm(load) || 1)) { it++; break; }
}The flux density in each triangle and the stored energy per unit depth.
B = (∂A/∂y, −∂A/∂x); |B| = hypot; energy = Σ ½ ν(|B|²) |B|² × area
One constant B per triangle, so the reported peak is the largest element value and depends on the mesh. The energy sum is added in element order in a 64-bit number.
const Bx = new Float64Array(m), By = new Float64Array(m), Bmag = new Float64Array(m);
let energy = 0;
for (let e = 0; e < m; e++) {
let ax = 0, ay = 0; for (let i = 0; i < 3; i++) { const v = A[mesh.tris[3 * e + i]]; ax += G.gx[3 * e + i] * v; ay += G.gy[3 * e + i] * v; }
Bx[e] = ay; By[e] = -ax; Bmag[e] = Math.hypot(ax, ay);
energy += 0.5 * mats[region[e]].nu(Bmag[e] ** 2) * Bmag[e] ** 2 * G.area[e];
}The force per meter of depth and the torque about the origin on everything inside a set of elements.
t = ((B_x² − ½B²) n_x + B_x B_y n_y, B_x B_y n_x + (B_y² − ½B²) n_y) ÷ μ₀; F = Σ t × dl; τ = Σ (x t_y − y t_x) dl
B on each contour edge is the mean of the two triangles that share it. The force is a sum of signed edge terms in map order; where the net force is small beside the stresses on opposite sides, the terms cancel and the result depends on the mesh along the contour.
/**
* Force and torque per metre of depth on whatever lies inside the set of elements `inside(cx, cy) → bool` (it must be
* wholly surrounded by air): { fx, fy (N/m), torque (N, about the origin), edges }.
*/
export function maxwell(mesh, res, inside) {
const m = mesh.tris.length / 3, G = res.geom, X = mesh.nodes, side = new Uint8Array(m);
for (let e = 0; e < m; e++) side[e] = inside(G.cx[e], G.cy[e]) ? 1 : 0;
const edges = new Map();
for (let e = 0; e < m; e++) for (let k = 0; k < 3; k++) {
const a = mesh.tris[3 * e + k], b = mesh.tris[3 * e + (k + 1) % 3], key = a < b ? a * 4194304 + b : b * 4194304 + a;
const o = edges.get(key); if (o) o.push(e); else edges.set(key, [e]);
}
let fx = 0, fy = 0, tq = 0, cnt = 0;
for (const [key, es] of edges) {
if (es.length !== 2 || side[es[0]] === side[es[1]]) continue;
const a = Math.floor(key / 4194304), b = key % 4194304, ein = side[es[0]] ? es[0] : es[1];
const x1 = X[2 * a], y1 = X[2 * a + 1], x2 = X[2 * b], y2 = X[2 * b + 1], dl = Math.hypot(x2 - x1, y2 - y1);
let nx = (y2 - y1) / dl, ny = -(x2 - x1) / dl;
const mx = (x1 + x2) / 2, my = (y1 + y2) / 2;
if (nx * (mx - G.cx[ein]) + ny * (my - G.cy[ein]) < 0) { nx = -nx; ny = -ny; } // outward from the set
const bx = (res.Bx[es[0]] + res.Bx[es[1]]) / 2, by = (res.By[es[0]] + res.By[es[1]]) / 2, b2 = bx * bx + by * by;
const tx = ((bx * bx - 0.5 * b2) * nx + bx * by * ny) / MU0, ty = (bx * by * nx + (by * by - 0.5 * b2) * ny) / MU0;
fx += tx * dl; fy += ty * dl; tq += (mx * ty - my * tx) * dl; cnt++;
}
return { fx, fy, torque: tq, edges: cnt };
}Tests tools/partcraft-test/test-magnet-study.mjs
Converts each profile to meters and sets its material, magnet direction or coil current density.
x_m = x_mm × 1e-3; area_m² = area_mm² × 1e-3 × 1e-3; B_r = B (cos θ, sin θ); J = I ÷ area
Lengths are multiplied by 1e-3, which is not exactly representable, instead of divided by 1000; each coordinate can differ from the correctly rounded value in its last digit. Material inputs are clamped to fixed ranges. A coil with zero area gets no current.
const s = sp[i], id = i + 1, poly = cut[i].map(([x, y]) => [m(x), m(y)]), a = Math.abs(area(loops[i])) * MM * MM;
if (s.role === 'air') materials[id] = { type: 'linear', mur: 1 };
else if (s.role === 'iron') {
const c = parseCurve(s.curve);
if (String(s.curve || '').trim() && c.error) problems.push(t('study.mag.err.region', { n: id, why: c.error }));
if (c.points && c.points.length >= 2) materials[id] = { type: 'table', points: c.points.map(([B, H]) => [B, H]) };
else materials[id] = { type: 'linear', mur: clamp(s.mur, 1, LIMITS.mur, 2000) };
} else if (s.role === 'magnet') {
const Br = clamp(s.br, 0, LIMITS.br, 0), ang = clamp(s.angle, -3600, 3600, 0) * Math.PI / 180;
materials[id] = { type: 'linear', mur: clamp(s.murMag, 1, 100, 1.05) };
br[id] = [Br * Math.cos(ang), Br * Math.sin(ang)];
} else if (s.role === 'coil') {
materials[id] = { type: 'linear', mur: 1 };
const I = clamp(s.current, -LIMITS.current, LIMITS.current, 0);
if (a > 0) current[id] = I / a;The triangle size asked for in each body, in each air gap and in the surrounding air.
base = size × (1 + 2 pad) ÷ 14 ÷ q; body: narrowest side ÷ 5 ÷ q; gap: gap ÷ 4 ÷ q
These sizes decide how well the field and the force are resolved; q is the fineness setting (0.5, 1 or 2). A gap is refined only when it is under 0.3 of the model size.
const q = clamp(fine, 0.5, 2, 1), base = m(size) * (1 + 2 * clamp(pad, 0.25, 6, 1)) / 14 / q;
const zones = loops.map((p) => { const b = boxOf(p); return { box: b.map(m), h: Math.max(m(Math.max(1e-9, Math.min(b[2] - b[0], b[3] - b[1]))) / 5 / q, base / 40) }; });
for (const g of gaps(loops)) if (g.d < 0.3 * size) zones.push({ box: g.box.map(m), h: Math.max(m(g.d) / 4 / q, base / 120) });
const sizeSpec = { base, grow: 0.3, zones };The narrowest distance between two profiles that stand apart.
d = min over the corners of each profile of the distance to the other's edges
Corner-to-edge distances both ways, which gives the minimum between two polygons that do not cross. A distance not above 1e-6 (absolute, in sketch units) counts as touching. Sketch arcs enter as 48-segment polygons, so a gap next to an arc is measured to its chords.
/** Every pair of bodies that stand apart but close: [{ i, j, d (mm), box }] with `box` the region where both come within d of each other (the gap's own area). */
export function gaps(loops) {
const out = [], boxes = loops.map(boxOf);
for (let i = 0; i < loops.length; i++) for (let j = i + 1; j < loops.length; j++) {
let d = Infinity;
for (const [A, B] of [[loops[i], loops[j]], [loops[j], loops[i]]]) for (const [x, y] of A) for (let k = 0; k < B.length; k++) d = Math.min(d, segDist(x, y, B[k], B[(k + 1) % B.length]));
if (!(d > 1e-6)) continue;
const e = 1.5 * d, a = boxes[i], b = boxes[j];
const box = [Math.max(a[0], b[0]) - e, Math.max(a[1], b[1]) - e, Math.min(a[2], b[2]) + e, Math.min(a[3], b[3]) + e];
if (box[2] > box[0] && box[3] > box[1]) out.push({ i, j, d, box });
}
return out;
}The force in newtons and the torque about the body's own center, for the depth given.
F = F_per_m × depth; τ = (τ_origin − (c_x F_y − c_y F_x)) × depth
The torque is moved from the sketch origin to the center of the body's contour box by subtracting c × F, which cancels digits when the body is far from the origin. Depth is converted with × 1e-3.
// ---------------------------------------------------------------------------------------------------- the answer
export function summarize(res, meta) {
const d = meta.depth * MM; // metres of depth
const bodies = (res.forces || []).map((f) => {
const c = meta.contours.find((x) => x.id === f.id), L = meta.regions.find((r) => r.loop === c.loop);
const b = c.box, cx = (b[0] + b[2]) / 2, cy = (b[1] + b[3]) / 2; // the torque is asked about the body's own middle, not the sketch origin
const tau = f.torque - (cx * f.fy - cy * f.fx);
return { id: f.id, loop: c.loop + 1, role: L ? L.role : '', thin: !!c.thin, fx: f.fx * d, fy: f.fy * d, force: Math.hypot(f.fx, f.fy) * d, torque: tau * d, fxPerM: f.fx, fyPerM: f.fy, centre: [cx / MM, cy / MM] };
});
return { bmax: res.Bmax.value, bmaxAt: [res.Bmax.at[0] / MM, res.Bmax.at[1] / MM], iterations: res.iterations, tris: res.stats ? res.stats.tris : res.mesh.tris.length / 3, minAngle: res.stats ? res.stats.minAngle : null, bodies, depthMm: meta.depth, ms: res.ms, energy: res.energy };
}Writes the peak flux density, forces and torques to a CSV file.
number → toPrecision(8)
Every number in the exported file is rounded to 8 significant digits, fewer than the 15 to 17 the values hold.
export function resultCsv(sum) {
const q = (v) => (typeof v === 'number' ? String(+v.toPrecision(8)) : '"' + String(v).replace(/"/g, '""') + '"');
const rows = [['quantity', 'value', 'unit'], ['Bmax', q(sum.bmax), 'T'], ['Bmax at x', q(sum.bmaxAt[0]), 'mm'], ['Bmax at y', q(sum.bmaxAt[1]), 'mm'], ['depth', q(sum.depthMm), 'mm'], ['triangles', q(sum.tris), ''], ['iterations', q(sum.iterations), '']];
for (const b of sum.bodies) rows.push([`region ${b.loop} Fx`, q(b.fx), 'N'], [`region ${b.loop} Fy`, q(b.fy), 'N'], [`region ${b.loop} torque about its centre`, q(b.torque), 'N·m']);
return '' + rows.map((r) => r.join(',')).join('\r\n') + '\r\n';
}Tests tools/magnet-test/mesh.mjs
Which side of a directed line a point lies on, used to walk the triangulation.
(b_x − a_x)(c_y − a_y) − (b_y − a_y)(c_x − a_x)
A plain 64-bit determinant; its sign can be wrong for nearly collinear points and there is no exact or adaptive fallback. A walk that cycles because of this falls back to scanning every triangle.
const orient = (ax, ay, bx, by, cx, cy) => (bx - ax) * (cy - ay) - (by - ay) * (cx - ax);Whether a new point lies inside the circle through a triangle's corners, which decides the Delaunay cavity.
3 × 3 determinant of (a − d, b − d, c − d) with their squared lengths
Coordinates are taken relative to the tested point first, then the determinant is evaluated in plain 64-bit arithmetic with no exact fallback. The super triangle's corners are 1000 to 2000 box diagonals away, so tests that involve them work with much larger numbers than the model's.
function inCircle(ax, ay, bx, by, cx, cy, dx, dy) { // > 0 when d is inside the circle through counter-clockwise a, b, c
const adx = ax - dx, ady = ay - dy, bdx = bx - dx, bdy = by - dy, cdx = cx - dx, cdy = cy - dy;
return (adx * adx + ady * ady) * (bdx * cdy - cdx * bdy) + (bdx * bdx + bdy * bdy) * (cdx * ady - adx * cdy) + (cdx * cdx + cdy * cdy) * (adx * bdy - bdx * ady);
}A triangle's circumcenter, circumradius and shortest edge, used to decide whether to refine it.
u = (c_y |b|² − b_y |c|², b_x |c|² − c_x |b|²) ÷ (2 (b_x c_y − b_y c_x)), b and c relative to a
Corners are taken relative to the first corner before the squares are formed, which keeps digits for triangles far from the origin. A triangle whose denominator is under 1e-300 in magnitude is skipped.
const circ = (t) => {
const [a, b, c] = tri(t), ax = X[a], ay = Y[a], bx = X[b] - ax, by = Y[b] - ay, cx = X[c] - ax, cy = Y[c] - ay, d = 2 * (bx * cy - by * cx);
if (Math.abs(d) < 1e-300) return null;
const ux = (cy * (bx * bx + by * by) - by * (cx * cx + cy * cy)) / d, uy = (bx * (cx * cx + cy * cy) - cx * (bx * bx + by * by)) / d;
const sh = Math.min(Math.hypot(bx, by), Math.hypot(cx, cy), Math.hypot(bx - cx, by - cy));
return { x: ax + ux, y: ay + uy, R: Math.hypot(ux, uy), sh, gx: (X[a] + X[b] + X[c]) / 3, gy: (Y[a] + Y[b] + Y[c]) / 3 };
};Whether a point lies inside the circle whose diameter is a boundary segment, in which case the segment is split at its midpoint.
(p − mid)² < (|a − b|² ÷ 4) × (1 − 1e-12)
The factor 1 − 1e-12 keeps a point that is on the circle within rounding from counting as inside. The midpoint is (a + b) ÷ 2, with one rounding in the sum.
// ---- encroachment
const encroaches = (s, px, py) => {
const [a, b] = s; if (px === X[a] && py === Y[a]) return false; if (px === X[b] && py === Y[b]) return false;
const mx = (X[a] + X[b]) / 2, my = (Y[a] + Y[b]) / 2, r2 = ((X[a] - X[b]) ** 2 + (Y[a] - Y[b]) ** 2) / 4;
return (px - mx) ** 2 + (py - my) ** 2 < r2 * (1 - 1e-12);
};The mesh quality figures shown with the result.
angle = atan2(|u × v|, u·v) in degrees
atan2 of the cross and dot products keeps its accuracy for very small and very flat angles, where an acos of the normalized dot product would not.
const ang = Math.atan2(Math.abs(ux * vy - uy * vx), ux * vx + uy * vy) * 180 / Math.PI;
minA = Math.min(minA, ang); maxA = Math.max(maxA, ang);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.