Footing design checks

/** The concrete's f'c in MPa from the footing's text: "4000 psi", "28 MPa", "C30/37" (f'ck of the cylinder). */
export function fcOf(text) {
  const s = String(text || '');
  let m = s.match(/(\d+(?:\.\d+)?)\s*psi/i); if (m) return Number(m[1]) * 0.00689476;
  m = s.match(/C(\d+)\s*\/\s*\d+/i); if (m) return Number(m[1]);
  m = s.match(/(\d+(?:\.\d+)?)\s*MPa/i); if (m) return Number(m[1]);
  return 28;
}

/** A support reaction of the Analysis (N, N·mm, global axes, the force the support puts on the structure) as the load on the footing under it. */
export function toFootingAxes(R, rotationDeg = 0) {
  const a = rotationDeg * Math.PI / 180, c = Math.cos(a), s = Math.sin(a);
  const F = [-R[0], -R[1], -R[2]], M = [-R[3], -R[4], -R[5]];   // what the structure puts on the footing
  const loc = (u, v) => [c * u + s * v, -s * u + c * v];
  const [Vx, Vy] = loc(F[0], F[1]), [Mx, My] = loc(M[0], M[1]);
  return { P: -F[2] / 1000, Vx: Vx / 1000, Vy: Vy / 1000, Mx: Mx / 1e6, My: My / 1e6 };
}

const dX = H - c - bx.db / 2, dY = H - c - bx.db - by.db / 2, dAvg = (dX + dY) / 2;

const W = soil.gammaConc * A * h + soil.gammaSoil * A * soil.fill / 1000;
const N = sumP + W;
const moment = (key) => sv.reduce((s, l) => s + (l.service?.[key] || 0), 0);
const Vx = moment('Vx'), Vy = moment('Vy');
const ex = N > 0 ? (sv.reduce((s, l) => s + (l.service?.P || 0) * l.x / 1000, 0) + moment('My') + h * Vx) / N * 1000 : 0;
const ey = N > 0 ? (sv.reduce((s, l) => s + (l.service?.P || 0) * l.y / 1000, 0) - moment('Mx') - h * Vy) / N * 1000 : 0;
const field = pressureField(N, ex, ey, L, B);

/** q(x, y) = a + b·x + c·y (kPa, x and y in metres from the middle) wherever that is positive. */
export function pressureField(N, ex, ey, L, B) {
  const Lm = L / 1000, Bm = B / 1000, exm = ex / 1000, eym = ey / 1000;
  if (!(N > 0)) return { ok: false, reason: 'no-compression', N };
  if (Math.abs(exm) >= Lm / 2 || Math.abs(eym) >= Bm / 2) return { ok: false, reason: 'outside', N, ex: exm, ey: eym };
  const A = Lm * Bm;
  const full = { a: N / A, b: 12 * N * exm / (A * Lm * Lm), c: 12 * N * eym / (A * Bm * Bm) };
  const make = (p, partial) => {
    const corners = [[-Lm / 2, -Bm / 2], [Lm / 2, -Bm / 2], [-Lm / 2, Bm / 2], [Lm / 2, Bm / 2]].map(([x, y]) => p.a + p.b * x + p.c * y);
    const qmax = Math.max(...corners);
    return { ok: true, ...p, qmax, qmin: partial ? 0 : Math.min(...corners), partial, N, ex: exm, ey: eym, L: Lm, B: Bm, q: (x, y) => Math.max(0, p.a + p.b * x + p.c * y) };
  };
  if (6 * Math.abs(exm) / Lm + 6 * Math.abs(eym) / Bm <= 1 + 1e-12) return make(full, false);
  // part of the base lifts off: Newton on the plane so that the positive part carries N at (ex, ey)
  const n = 80, dx = Lm / n, dy = Bm / n, dA = dx * dy;
  const xs = Array.from({ length: n }, (_, i) => -Lm / 2 + (i + 0.5) * dx), ys = Array.from({ length: n }, (_, j) => -Bm / 2 + (j + 0.5) * dy);
  const F = (p) => {
    const f = [-N, -N * exm, -N * eym], J = Array.from({ length: 3 }, () => [0, 0, 0]);
    for (const y of ys) for (const x of xs) {
      const q = p.a + p.b * x + p.c * y;
      if (q <= 0) continue;
      const phi = [1, x, y];
      for (let i = 0; i < 3; i++) { f[i] += q * phi[i] * dA; for (let k = 0; k < 3; k++) J[i][k] += phi[i] * phi[k] * dA; }
    }
    return { f, J };
  };
  const norm = (f) => Math.hypot(f[0], f[1] / (Lm / 2), f[2] / (Bm / 2)) / N;
  let p = { ...full }, cur = F(p);
  for (let it = 0; it < 200 && norm(cur.f) > 1e-7; it++) {
    const J = cur.J.map((r, i) => r.map((v, k) => v + (i === k ? 1e-9 * (r[0] + 1e-12) : 0)));
    const d = solve3(J, cur.f.map((v) => -v));
    if (!d) break;
    let step = 1, next = null;
    for (let h = 0; h < 30; h++) {
      const q = { a: p.a + step * d[0], b: p.b + step * d[1], c: p.c + step * d[2] };
      const nf = F(q);
      if (norm(nf.f) < norm(cur.f)) { next = { q, nf }; break; }
      step /= 2;
    }
    if (!next) break;
    p = next.q; cur = next.nf;
  }
  return { ...make(p, true), residual: norm(cur.f) };
}

function solve3(A, b) {
  const M = A.map((r, i) => [...r, b[i]]);
  for (let c = 0; c < 3; c++) {
    let piv = c;
    for (let r = c + 1; r < 3; r++) if (Math.abs(M[r][c]) > Math.abs(M[piv][c])) piv = r;
    if (Math.abs(M[piv][c]) < 1e-14) return null;
    [M[c], M[piv]] = [M[piv], M[c]];
    for (let r = c + 1; r < 3; r++) { const k = M[r][c] / M[c][c]; for (let q = c; q < 4; q++) M[r][q] -= k * M[c][q]; }
  }
  const x = [0, 0, 0];
  for (let r = 2; r >= 0; r--) { let s = M[r][3]; for (let q = r + 1; q < 3; q++) s -= M[r][q] * x[q]; x[r] = s / M[r][r]; }
  return x;
}

/** ∫ q over the part of the base with x > xc (or y > yc), with the first moment about the cut: { V (kN), M (kN·m) }. */
function strip(field, axis, cut, n = 200) {
  const Lm = field.L, Bm = field.B;
  const len = axis === 'x' ? Lm : Bm, wid = axis === 'x' ? Bm : Lm;
  const t0 = Math.max(cut, -len / 2), t1 = len / 2;
  if (t1 <= t0) return { V: 0, M: 0 };
  const ds = (t1 - t0) / n, m = 80, dw = wid / m;
  let V = 0, M = 0;
  for (let i = 0; i < n; i++) {
    const t = t0 + (i + 0.5) * ds;
    let q = 0;
    for (let j = 0; j < m; j++) { const w = -wid / 2 + (j + 0.5) * dw; q += axis === 'x' ? field.q(t, w) : field.q(w, t); }
    const f = q * dw * ds;
    V += f; M += f * (t - cut);
  }
  return { V, M };
}

const mk = (id, clause, demand, capacity, unit, vars = {}, extra = {}) => {
  const ratio = capacity > 0 ? demand / capacity : (demand > 0 ? Infinity : 0);
  return { id, clause, demand, capacity, unit, ratio, status: ratio <= 1 + 1e-9 ? 'ok' : 'fail', vars, ...extra };
};

const H0 = Math.hypot(Vx, Vy);
if (H0 > 1e-9) {
  if (soil.mu > 0 && N > 0) checks.push(mk('sliding', 'ACI 318-19 13.2.6.1', soil.fsSlide * H0, soil.mu * N, 'kN', { fs: r2(soil.mu * N / H0), need: soil.fsSlide, H: r1(H0), N: r1(N) }));
  else checks.push(note('sliding', 'ACI 318-19 13.2.6.1', { H: r1(H0) }));
}
if (N > 0) {
  const eMax = Math.max(Math.abs(ex) / (L / 2), Math.abs(ey) / (B / 2));   // 1 = at the edge
  const fsOt = eMax > 0 ? 1 / eMax : Infinity;
  if (eMax > 1e-6) checks.push(mk('overturning', 'ACI 318-19 13.2.6.1', soil.fsOt, fsOt, '', { fs: r2(fsOt), need: soil.fsOt }));

const cutX = (xc) => {   // loads on the +x side of the cut
  const s = strip(fu, 'x', xc);
  let V = s.V, M = s.M;
  for (const k of cols) if (k.x > xc) { V -= k.P; M -= k.P * (k.x - xc) + k.My + h * k.Vx; }
  return { V, M };
};
const cutY = (yc) => {
  const s = strip(fu, 'y', yc);
  let V = s.V, M = s.M;
  for (const k of cols) if (k.y > yc) { V -= k.P; M -= k.P * (k.y - yc) - k.Mx - h * k.Vy; }
  return { V, M };
};
const insideX = (xc) => cols.some((k) => xc > k.x - k.cx / 2 + 1e-9 && xc < k.x + k.cx / 2 - 1e-9);
const insideY = (yc) => cols.some((k) => yc > k.y - k.cy / 2 + 1e-9 && yc < k.y + k.cy / 2 - 1e-9);
const scan = (axis) => {
  const len = (axis === 'x' ? L : B) / 1000, cut = axis === 'x' ? cutX : cutY, inside = axis === 'x' ? insideX : insideY;
  const faces = cols.flatMap((k) => (axis === 'x' ? [k.x - k.cx / 2, k.x + k.cx / 2] : [k.y - k.cy / 2, k.y + k.cy / 2]));
  const at = [...faces, ...Array.from({ length: 41 }, (_, i) => -len / 2 + (i + 0.5) * len / 41)].filter((v) => v > -len / 2 + 1e-6 && v < len / 2 - 1e-6 && !inside(v));
  let best = { M: 0, at: 0 };
  for (const v of at) { const r = cut(v); if (Math.abs(r.M) > Math.abs(best.M)) best = { M: r.M, at: v }; }
  return best;
};

const MuPerM = Math.abs(sc.M) / (b / 1000) * 1e6;   // N·mm per metre of width
const AsProv = lay.bar.area * 1000 / lay.s;
const AsReq = MuPerM > 0 ? requiredAs({ b: 1000, d: lay.d, fc, fy, Mu: MuPerM, phi: PHI_FLEX }) : 0;
const rhoMin = fy < 420 ? 0.0020 : Math.max(0.0014, 0.0018 * 420 / fy);
const AsMin = rhoMin * 1000 * H;

const cutsV = cols.flatMap((k) => (axis === 'x' ? [k.x - k.cx / 2 - lay.d / 1000, k.x + k.cx / 2 + lay.d / 1000] : [k.y - k.cy / 2 - lay.d / 1000, k.y + k.cy / 2 + lay.d / 1000])).filter((v) => Math.abs(v) < (axis === 'x' ? L : B) / 2000 - 1e-6 && !(axis === 'x' ? insideX(v) : insideY(v)));
let Vmax = 0;
for (const v of cutsV) { const r = axis === 'x' ? cutX(v) : cutY(v); Vmax = Math.max(Vmax, Math.abs(r.V)); }
if (cutsV.length) {
  const bs = beamShear({ bw: b, d: lay.d, h: H, fc, fyt: 420, Av: 0, s: 0, Vu: Vmax * 1000, As: AsProv * b / 1000, lambda: lam });
  checks.push(mk('shear-' + axis, 'ACI 318-19 13.2.6.2, 22.5.5.1', Vmax, PHI_SHEAR * bs.Vc / 1000, 'kN', { dir: axis, Vu: r1(Vmax), phiVc: r1(PHI_SHEAR * bs.Vc / 1000), d: Math.round(lay.d) }));

/** Straight development length of a bar in tension, ACI 318-19 Table 25.4.2.3 with ψt = ψe = 1 (bottom bars, uncoated), mm. `easy` = spacing and cover at least db. */
export function developmentLength(bar, fy, fc, lambda = 1, easy = true) {
  const small = bar.db <= 19.1;   // #6 and smaller
  const k = easy ? (small ? 2.1 : 1.7) : (small ? 1.4 : 1.1);
  return Math.max(300, fy / (k * lambda * Math.sqrt(fc)) * bar.db);
}

/** Hook development, ACI 318-19 25.4.3.1: 0.24·ψe·fy/(λ√f'c)·db, at least 8 db and 150 mm. */
export function hookDevelopment(bar, fy, fc, lambda = 1) { return Math.max(8 * bar.db, 150, 0.24 * fy / (lambda * Math.sqrt(fc)) * bar.db); }

const face = Math.max(...cols.map((k) => (axis === 'x' ? Math.abs(k.x) + k.cx / 2 : Math.abs(k.y) + k.cy / 2))) * 1000;
const avail = (axis === 'x' ? L : B) / 2 - face - c;
const ld = developmentLength(lay.bar, fy, fc, lam, lay.s - lay.bar.db >= lay.bar.db && c >= lay.bar.db);

/**
 * The critical section for punching at d/2 from a column face (ACI 318-19 22.6.4.1), as line segments of the slab depth d, with the sides that
 * fall outside the footing dropped and the others cut off at its edge. The footing runs from −L/2 to L/2 and −B/2 to B/2 (m). Returns the
 * segments, their total length bo (m), the centroid, and the section moduli J about each axis through the centroid (m⁴, ACI R22.6.5.3: the sides
 * along the bending direction add d·l³/12 + l·d³/12, the others only the transfer term d·l·arm²) with the farthest distance c to a segment end.
 */
export function punchSection(k, d, L, B) {
  const bx = k.cx + d, by = k.cy + d, x0 = k.x - bx / 2, x1 = k.x + bx / 2, y0 = k.y - by / 2, y1 = k.y + by / 2;
  const e = 1e-9, segs = [];
  const ya = Math.max(y0, -B / 2), yb = Math.min(y1, B / 2), xa = Math.max(x0, -L / 2), xb = Math.min(x1, L / 2);
  if (yb > ya + e) for (const x of [x0, x1]) if (x > -L / 2 + e && x < L / 2 - e) segs.push({ dir: 'y', len: yb - ya, mx: x, my: (ya + yb) / 2 });
  if (xb > xa + e) for (const y of [y0, y1]) if (y > -B / 2 + e && y < B / 2 - e) segs.push({ dir: 'x', len: xb - xa, mx: (xa + xb) / 2, my: y });
  const bo = segs.reduce((a, g) => a + g.len, 0);
  if (!(bo > e)) return { segs, bo: 0, x0, x1, y0, y1, bx, by };
  const xc = segs.reduce((a, g) => a + g.len * g.mx, 0) / bo, yc = segs.reduce((a, g) => a + g.len * g.my, 0) / bo;
  let Jx = 0, Jy = 0, cx = 0, cy = 0;
  for (const g of segs) {
    Jx += (g.dir === 'x' ? d * g.len ** 3 / 12 + g.len * d ** 3 / 12 : 0) + d * g.len * (g.mx - xc) ** 2;
    Jy += (g.dir === 'y' ? d * g.len ** 3 / 12 + g.len * d ** 3 / 12 : 0) + d * g.len * (g.my - yc) ** 2;
    cx = Math.max(cx, g.dir === 'x' ? Math.abs(g.mx - xc) + g.len / 2 : Math.abs(g.mx - xc));
    cy = Math.max(cy, g.dir === 'y' ? Math.abs(g.my - yc) + g.len / 2 : Math.abs(g.my - yc));
  }
  return { segs, bo, xc, yc, Jx, Jy, cx, cy, x0, x1, y0, y1, bx, by };
}

/**
 * Punching of one column (k: m and kN, kN·m) on the strength pressure field fu. The stress is Vu/(bo·d) plus γv·M·c/J in each direction, M being the
 * column's moment plus Vu times the distance from the column to the section's centroid (edge and corner columns), the two directions and the two
 * terms added without signs: slightly conservative. vc is the least of ACI 318-19 Table 22.6.5.2 (a)–(c) with the size effect λs.
 */
function punching(k, n, { fu, L, B, dAvg, fcr, lam }) {
  const d = dAvg / 1000, sec = punchSection(k, d, L / 1000, B / 1000);
  if (!(sec.bo > 0)) return [note('punching', 'ACI 318-19 22.6', { n, why: 'edge' })];
  const where = sec.segs.length >= 4 ? 'interior' : sec.segs.length === 3 ? 'edge' : 'corner';
  let inside = 0;
  const xa = Math.max(sec.x0, -L / 2000), xb = Math.min(sec.x1, L / 2000), ya = Math.max(sec.y0, -B / 2000), yb = Math.min(sec.y1, B / 2000);
  const m = 24, dx = (xb - xa) / m, dy = (yb - ya) / m;
  for (let a = 0; a < m; a++) for (let q = 0; q < m; q++) inside += fu.q(xa + (a + 0.5) * dx, ya + (q + 0.5) * dy) * dx * dy;
  const Vu = Math.max(0, k.P - inside);
  const gvx = 1 - 1 / (1 + (2 / 3) * Math.sqrt(sec.bx / sec.by)), gvy = 1 - 1 / (1 + (2 / 3) * Math.sqrt(sec.by / sec.bx));
  const eX = Math.abs(k.x - sec.xc), eY = Math.abs(k.y - sec.yc);
  const vu = Vu / (sec.bo * d) / 1000 + gvx * (Math.abs(k.My) + Vu * eX) * sec.cx / sec.Jx / 1000 + gvy * (Math.abs(k.Mx) + Vu * eY) * sec.cy / sec.Jy / 1000;   // MPa
  const beta = Math.max(k.cx, k.cy) / Math.min(k.cx, k.cy), lamS = Math.min(1, Math.sqrt(2 / (1 + 0.004 * dAvg)));
  const alphaS = where === 'interior' ? 40 : where === 'edge' ? 30 : 20;
  const vc = Math.min(0.33, 0.17 * (1 + 2 / beta), 0.083 * (2 + alphaS * d / sec.bo)) * lamS * lam * fcr;
  return [mk('punching', 'ACI 318-19 13.2.6.3, 22.6.5.2', vu, PHI_SHEAR * vc, 'MPa', { n, where, Vu: r1(Vu), vu: r2(vu), phivc: r2(PHI_SHEAR * vc), bo: Math.round(sec.bo * 1000), d: Math.round(dAvg) })];
}

const contactShare = (field) => {
  const n = 40; let k = 0;
  for (let i = 0; i < n; i++) for (let j = 0; j < n; j++) if (field.a + field.b * (-field.L / 2 + (i + 0.5) * field.L / n) + field.c * (-field.B / 2 + (j + 0.5) * field.B / n) > 0) k++;
  return k / (n * n);
};