Foundation settlement

/** Vertical stress increase under the corner of a uniformly loaded rectangle a × b at depth z (Newmark 1935): Δσ = q·I, I from the closed form. */
export function influenceCorner(a, b, z) {
  if (z <= 0) return 0.25;
  const m = a / z, n = b / z, m2 = m * m, n2 = n * n, v = m2 + n2 + 1, v1 = m2 * n2, den = v - v1;
  const ang = Math.atan2(2 * m * n * Math.sqrt(v), den);
  return ((2 * m * n * Math.sqrt(v)) / (v + v1) * ((v + 1) / v) + (den > 0 ? ang : ang + Math.PI * 0) ) / (4 * Math.PI);
}

/** Under the centre of a B × L footing: four corner rectangles of B/2 × L/2. Under a corner: one. */
export const stressRect = (q, B, L, z, at = 'centre') => (at === 'corner' ? q * influenceCorner(B, L, z) : 4 * q * influenceCorner(B / 2, L / 2, z));

/** 2:1 spread: Δσ = q·B·L / ((B+z)(L+z)). */
export const stress2to1 = (q, B, L, z) => (q * B * L) / ((B + z) * (L + z));

/** Elastic settlement of the corner of a flexible a × b rectangle carrying q on a half-space of modulus E and Poisson ν (Timoshenko and Goodier). */
export function elasticCorner(q, a, b, E, nu) {
  const d = Math.hypot(a, b);
  return ((q * (1 - nu * nu)) / (Math.PI * E)) * (a * Math.log((b + d) / a) + b * Math.log((a + d) / b));
}

/**
 * One-dimensional settlement of a stack of layers under a footing.
 *   layers: [{ top (elevation), bottom, soil kind: 'clay'|'sand', gamma, gammaSat, cc, cs, e0, sigmaC (kPa, effective), es, nu, cv, drain }]
 *   load: { q (kPa, net), B, L, depth (footing base elevation) , at: 'centre'|'corner', method: 'boussinesq' | '2:1' }
 *   water: elevation of the water table (null = none)
 * Compressible layers (cc given) are cut into sub-layers of at most `step` metres; granular layers (es given) give immediate settlement from
 * Δσ/Es·Δz (the strain method). Returns { consolidation, immediate, total, rows: [{ z, sigma0, dsigma, s, kind, layer }] }.
 */
export function settlement({ layers, load, water = null, ground = null, step = 1 }) {
  const rows = []; let sc = 0, si = 0;
  const base = load.depth ?? ground ?? layers[0].top;
  const sigmaAt = (z) => {   // effective vertical stress at elevation z, from the top of the first layer
    let s = 0, zt = ground ?? layers[0].top;
    for (const l of layers) {
      const hi = Math.min(zt, l.top), lo = Math.max(z, l.bottom);
      if (hi > lo) {
        const dry = water === null ? hi - lo : Math.max(0, Math.min(hi, Math.max(lo, water)) - lo) * 0 + Math.max(0, hi - Math.max(lo, water)), wet = hi - lo - dry;
        s += dry * l.gamma + wet * ((l.gammaSat ?? l.gamma) - G_W);
      }
      if (l.bottom <= z) break;
    }
    return s;
  };
  for (const [li, l] of layers.entries()) {
    const top = Math.min(l.top, base), bot = l.bottom;
    if (top <= bot) continue;
    const n = Math.max(1, Math.ceil((top - bot) / step)), dz = (top - bot) / n;
    for (let i = 0; i < n; i++) {
      const zm = top - (i + 0.5) * dz, depth = base - zm;
      const dsig = load.method === '2:1' ? stress2to1(load.q, load.B, load.L, depth) : stressRect(load.q, load.B, load.L, depth, load.at);
      if (l.cc !== undefined && l.e0 !== undefined) {
        const s0 = Math.max(sigmaAt(zm), 0.1), sf = s0 + dsig, sc0 = l.sigmaC ?? s0, cs = l.cs ?? l.cc / 5;
        let s;
        if (sc0 <= s0 * 1.0001) s = (l.cc * dz) / (1 + l.e0) * Math.log10(sf / s0);   // normally consolidated
        else if (sf <= sc0) s = (cs * dz) / (1 + l.e0) * Math.log10(sf / s0);   // stays below the preconsolidation stress
        else s = ((cs * dz) / (1 + l.e0)) * Math.log10(sc0 / s0) + ((l.cc * dz) / (1 + l.e0)) * Math.log10(sf / sc0);
        rows.push({ layer: li, z: zm, depth, sigma0: s0, dsigma: dsig, s, kind: 'consolidation' }); sc += s;
      } else if (l.es) {
        const s = (dsig / l.es) * dz * (1 - (l.nu ?? 0) * 0);
        rows.push({ layer: li, z: zm, depth, sigma0: Math.max(sigmaAt(zm), 0), dsigma: dsig, s, kind: 'immediate' }); si += s;
      }
    }
  }
  return { consolidation: sc, immediate: si, total: sc + si, rows };
}

/** Degree of consolidation for a time factor (Terzaghi, series). */
export function degreeOfConsolidation(Tv) {
  if (Tv <= 0) return 0;
  let s = 0; for (let m = 0; m < 200; m++) { const M = (Math.PI * (2 * m + 1)) / 2; s += (2 / (M * M)) * Math.exp(-M * M * Tv); if (M * M * Tv > 60) break; }
  return 1 - s;
}

/** Time factor for a degree of consolidation: by the usual closed forms below 60 %, bisection above. */
export function timeFactor(U) {
  if (U <= 0) return 0; if (U >= 1) return Infinity;
  let lo = 0, hi = 10; for (let i = 0; i < 80; i++) { const mid = (lo + hi) / 2; if (degreeOfConsolidation(mid) < U) lo = mid; else hi = mid; }
  return (lo + hi) / 2;
}

/** Time for U of the final consolidation settlement: t = Tv·Hdr² / cv. Hdr is the longest drainage path (m), cv in m²/year → years. */
export const timeFor = (U, Hdr, cv) => (timeFactor(U) * Hdr * Hdr) / cv;

/** Settlement curve: [{ t, U, s }] for a final consolidation settlement `sFinal`. */
export function timeCurve(sFinal, Hdr, cv, tMax, n = 40) {
  const out = []; for (let i = 0; i <= n; i++) { const t = (tMax * i) / n, U = degreeOfConsolidation((cv * t) / (Hdr * Hdr)); out.push({ t, U, s: sFinal * U }); }
  return out;
}

/** Secondary compression between t1 and t2 (Cα as e per log cycle, ep the void ratio at the end of primary, H the layer thickness). */
export const secondary = (Calpha, H, ep, t1, t2) => ((Calpha * H) / (1 + ep)) * Math.log10(t2 / t1);