Boring logs and the soil profile between borings

const layers = (Array.isArray(o.layers) ? o.layers : []).slice(0, 80).map(layerOf).filter(Boolean).sort((a, b) => a.from - b.from);
// layers must be end to end: a gap is closed by stretching the layer above; an overlap is cut at the layer below's top
for (let i = 1; i < layers.length; i++) { if (layers[i].from > layers[i - 1].to) layers[i - 1].to = layers[i].from; else if (layers[i].from < layers[i - 1].to) layers[i - 1].to = layers[i].from; }
const kept = layers.filter((l) => l.to > l.from);

export const boringRows = (b) => b.layers.map((l) => ({ ...l, top: b.z - l.from, bottom: b.z - l.to, thickness: l.to - l.from }));

// SPT corrections: Pa = 101.325 kPa. Defaults are the common US safe-hammer values; every factor is the user's.
export function sptCorrect(n, { er = 60, cb = 1, cs = 1, cr = 1, sigma = 100 } = {}) {
  const n60 = n * (er / 60) * cb * cs * cr, cn = Math.min(1.7, Math.sqrt(101.325 / Math.max(sigma, 1)));
  return { n60, cn, n160: n60 * cn };
}

/**
 * The layers between borings along a line. Each boring's layers are named by `unit` (a unit is one stratum that runs between borings); a unit
 * found in two neighbouring borings has its top and bottom joined by straight lines; a unit found in only one pinches out half way to the next boring.
 * Borings are projected onto the line (station = along, offset = across); those further than `width` from the line are left out.
 * line: [{ x, y }, …] (a polyline); returns { length, borings: [{ id, name, sta, off }], strata: [{ unit, soil, color, polygons: [[{ sta, z }]] }], ground: [{ sta, z }] }
 */
export function projectToLine(line, p) {
  let best = null, run = 0;
  for (let i = 0; i + 1 < line.length; i++) {
    const a = line[i], b = line[i + 1], dx = b.x - a.x, dy = b.y - a.y, L2 = dx * dx + dy * dy, L = Math.sqrt(L2);
    if (L2 === 0) continue;
    const t = clamp(((p.x - a.x) * dx + (p.y - a.y) * dy) / L2, 0, 1), qx = a.x + t * dx, qy = a.y + t * dy, off = Math.hypot(p.x - qx, p.y - qy);
    if (!best || off < best.off) best = { sta: run + t * L, off };
    run += L;
  }
  return best;
}

export function sectionProfile(borings, line, { width = Infinity } = {}) {
  const pts = [];
  for (const b of borings) { const q = projectToLine(line, b); if (q && q.off <= width) pts.push({ b, sta: q.sta, off: q.off }); }
  pts.sort((p, q) => p.sta - q.sta);
  let length = 0; for (let i = 0; i + 1 < line.length; i++) length += Math.hypot(line[i + 1].x - line[i].x, line[i + 1].y - line[i].y);
  if (pts.length < 1) return { length, borings: [], strata: [], ground: [], problems: [{ code: 'geo.profile.none' }] };
  // the order of the units, top to bottom, over all borings: first appearance fixes the position
  const order = [];
  for (const p of pts) for (const l of p.b.layers) if (!order.includes(l.unit)) order.push(l.unit);
  const top = (p, unit) => { const l = p.b.layers.find((q) => q.unit === unit); return l ? { top: p.b.z - l.from, bottom: p.b.z - l.to, soil: l.soil } : null; };
  const problems = [];
  const strata = order.map((unit) => ({ unit, soil: null, polygons: [] }));
  for (const s of strata) {
    // where this unit exists: a run of consecutive borings
    const at = pts.map((p) => top(p, s.unit));
    s.soil = at.find(Boolean)?.soil || null;
    // each pair of neighbours (and the free ends): a polygon between the two top lines and the two bottom lines
    const edge = (i) => ({ lo: i > 0 ? (pts[i - 1].sta + pts[i].sta) / 2 : pts[i].sta, hi: i + 1 < pts.length ? (pts[i].sta + pts[i + 1].sta) / 2 : pts[i].sta });
    for (let i = 0; i < pts.length; i++) {
      // from this boring to the next
      if (i + 1 < pts.length) {
        const a = at[i], c = at[i + 1], sa = pts[i].sta, sb = pts[i + 1].sta, mid = (sa + sb) / 2;
        if (a && c) s.polygons.push([{ sta: sa, z: a.top }, { sta: sb, z: c.top }, { sta: sb, z: c.bottom }, { sta: sa, z: a.bottom }]);
        else if (a) s.polygons.push([{ sta: sa, z: a.top }, { sta: mid, z: a.top - (a.top - a.bottom) * 0 + 0 }, { sta: sa, z: a.bottom }].map((q, k) => (k === 1 ? { sta: mid, z: (a.top + a.bottom) / 2 } : q)));
        else if (c) s.polygons.push([{ sta: sb, z: c.top }, { sta: mid, z: (c.top + c.bottom) / 2 }, { sta: sb, z: c.bottom }]);
      }
    }
    void edge;
  }
  const ground = pts.map((p) => ({ sta: p.sta, z: p.b.z }));
  const water = pts.filter((p) => p.b.water !== null).map((p) => ({ sta: p.sta, z: p.b.z - p.b.water }));
  return { length, borings: pts.map((p) => ({ id: p.b.id, name: p.b.name, sta: p.sta, off: p.off, depth: p.b.depth, z: p.b.z })), strata, ground, water, problems };
}

/** The elevation of the top of every unit at one station of a profile: [{ unit, soil, top, bottom }] top to bottom (the stack the slope code reads). */
export function stackAt(profile, sta) {
  const out = [];
  for (const s of profile.strata) {
    for (const poly of s.polygons) {
      const xs = poly.map((q) => q.sta), lo = Math.min(...xs), hi = Math.max(...xs);
      if (sta < lo - 1e-9 || sta > hi + 1e-9) continue;
      const cuts = [];
      for (let i = 0; i < poly.length; i++) {
        const a = poly[i], c = poly[(i + 1) % poly.length];
        if (Math.abs(a.sta - c.sta) < 1e-12) continue;
        const t = (sta - a.sta) / (c.sta - a.sta);
        if (t >= -1e-9 && t <= 1 + 1e-9) cuts.push(a.z + t * (c.z - a.z));
      }
      if (cuts.length >= 2) { out.push({ unit: s.unit, soil: s.soil, top: Math.max(...cuts), bottom: Math.min(...cuts) }); break; }
    }
  }
  return out.sort((p, q) => q.top - p.top);
}