Member axes and section properties

export const G_ACC = 9.80665;                       // m/s², standard gravity

/** Local axes of a member: ex along it, ez "up" (or +X for a vertical member), ey = ez × ex, then turned by roll. */
export function axesOf(pi, pj, rollDeg = 0) {
  const d = [pj[0] - pi[0], pj[1] - pi[1], pj[2] - pi[2]];
  const L = Math.hypot(d[0], d[1], d[2]);
  const ex = [d[0] / L, d[1] / L, d[2] / L];
  // "vertical" exactly as DetailCraft's geom.js frame() decides it (|x·Z| > 0.999999), so the analysis axes are the drawn ones
  let ez = Math.abs(ex[2]) > 0.999999 ? [1, 0, 0] : unit([-ex[2] * ex[0], -ex[2] * ex[1], 1 - ex[2] * ex[2]]);
  let ey = cross(ez, ex);
  if (rollDeg) {
    const r = rollDeg * Math.PI / 180, c = Math.cos(r), s = Math.sin(r);
    const y2 = [c * ey[0] + s * ez[0], c * ey[1] + s * ez[1], c * ey[2] + s * ez[2]];
    const z2 = [-s * ey[0] + c * ez[0], -s * ey[1] + c * ez[1], -s * ey[2] + c * ez[2]];
    ey = y2; ez = z2;
  }
  return { L, ex, ey, ez };
}

// a zero-length member (both ends on one spot) cannot be analysed
const at = new Map(nodes.map((n) => [n.id, n]));
for (let k = members.length - 1; k >= 0; k--) {
  const a = at.get(members[k].i), b = at.get(members[k].j);
  if (Math.hypot(b.x - a.x, b.y - a.y, b.z - a.z) < 1e-6) { warn('zero-length', { id: members[k].id }); memIds.delete(members[k].id); members.splice(k, 1); }
}

/**
 * Area and second moments of polygons (Green's theorem, exact for straight edges): loops in (u, v); the first is the
 * outline (any winding), the rest are holes (taken off whatever their winding). Iyy = ∫ (v − cv)² dA (about the u axis
 * through the centroid), Izz = ∫ (u − cu)² dA, Iyz = ∫ (u − cu)(v − cv) dA; I1 ≥ I2 principal, angle (deg) from +u to the
 * principal axis whose moment of inertia is I1 — i.e. the axis the section bends about most stiffly.
 */
export function polySection(loops) {
  let A = 0, Su = 0, Sv = 0, Iuu = 0, Ivv = 0, Iuv = 0;    // Iuu = ∫v², Ivv = ∫u², Iuv = ∫uv about the origin
  loops.forEach((loop, k) => {
    let a = 0, su = 0, sv = 0, iuu = 0, ivv = 0, iuv = 0;
    for (let i = 0, n = loop.length; i < n; i++) {
      const [u0, v0] = loop[i], [u1, v1] = loop[(i + 1) % n];
      const c = u0 * v1 - u1 * v0;
      a += c;
      su += (u0 + u1) * c; sv += (v0 + v1) * c;
      ivv += (u0 * u0 + u0 * u1 + u1 * u1) * c;
      iuu += (v0 * v0 + v0 * v1 + v1 * v1) * c;
      iuv += (u0 * v1 + 2 * u0 * v0 + 2 * u1 * v1 + u1 * v0) * c;
    }
    const sign = (k === 0 ? 1 : -1) * Math.sign(a || 1);
    A += sign * a / 2; Su += sign * su / 6; Sv += sign * sv / 6;
    Ivv += sign * ivv / 12; Iuu += sign * iuu / 12; Iuv += sign * iuv / 24;
  });
  const cu = Su / A, cv = Sv / A;
  const Iyy = Iuu - A * cv * cv, Izz = Ivv - A * cu * cu, Iyz = Iuv - A * cu * cv;
  const m = (Iyy + Izz) / 2, r = Math.hypot((Iyy - Izz) / 2, Iyz);
  // the axis about which I = I1: bending about an axis at angle φ from +u gives I(φ) = Iyy cos²φ + Izz sin²φ − 2 Iyz sinφ cosφ
  const angle = 0.5 * Math.atan2(-2 * Iyz, Iyy - Izz) * 180 / Math.PI;
  return { A, cu, cv, Iyy, Izz, Iyz, I1: m + r, I2: m - r, angle: Math.abs(Iyz) < 1e-9 * (Iyy + Izz) ? 0 : angle };
}

/**
 * Saint-Venant torsion constant J of a catalogue-style section (mm⁴):
 *   rolled I-shapes with root fillets: El Darwish & Johnston's formula (below), as steel tables give J (within about 1 %);
 *   other open sections (I without fillets, C, L): J = Σ b t³ / 3 over the plates, flanges b × tf, web (d − tf) × tw or
 *     legs (b − t/2) × t (S. P. Timoshenko & J. N. Goodier, Theory of Elasticity, 3rd ed., §109, thin-walled open sections)
 *   closed rectangular tube: J = 4 Am² t / p, Am and p of the wall's mid-line (Bredt; Timoshenko & Goodier §115)
 *   round tube: J = π (D⁴ − (D − 2t)⁴) / 32 (exact)
 *   flat bar / plate b × t (b ≥ t): J = b t³ [1/3 − 0.21 (t/b)(1 − t⁴ / (12 b⁴))] (Roark's Formulas for Stress and
 *     Strain, torsion of a solid rectangle, K = a b³ [16/3 − 3.36 (b/a)(1 − b⁴/12a⁴)] with sides 2a × 2b)
 */
export function torsionJ(p) {
  switch (p.kind) {
    case 'I':
      // with root fillets: I. A. El Darwish and B. G. Johnston, "Torsion of structural shapes", Journal of the Structural
      // Division (ASCE) 91 (ST1), 1965 — the formula steel tables use: J = 2/3 B T³ + 1/3 (D − 2T) t³ + 2 α D₁⁴ − 0.420 T⁴,
      // α = −0.042 + 0.2204 t/T + 0.1355 r/T − 0.0865 r t/T² − 0.0725 t²/T², D₁ = ((T + r)² + (r + t/4) t) / (2r + T)
      if (p.r > 0) {
        const B = p.bf, T = p.tf, D = p.d, t = p.tw, r = p.r;
        const a = -0.042 + 0.2204 * t / T + 0.1355 * r / T - 0.0865 * r * t / (T * T) - 0.0725 * t * t / (T * T);
        const D1 = ((T + r) ** 2 + (r + t / 4) * t) / (2 * r + T);
        return 2 / 3 * B * T ** 3 + (D - 2 * T) * t ** 3 / 3 + 2 * a * D1 ** 4 - 0.420 * T ** 4;
      }
      return (2 * p.bf * p.tf ** 3 + (p.d - p.tf) * p.tw ** 3) / 3;
    case 'C': return (2 * p.bf * p.tf ** 3 + (p.d - p.tf) * p.tw ** 3) / 3;
    case 'L': { const h = p.h || p.b; return ((p.b - p.t / 2) + (h - p.t / 2)) * p.t ** 3 / 3; }
    case 'RHS': { const t = p.t, Am = (p.h - t) * (p.b - t), per = 2 * ((p.h - t) + (p.b - t)); return 4 * Am * Am * t / per; }
    case 'CHS': return Math.PI * (p.D ** 4 - (p.D - 2 * p.t) ** 4) / 32;
    case 'PL': { const b = Math.max(p.b, p.t), t = Math.min(p.b, p.t); return b * t ** 3 * (1 / 3 - 0.21 * (t / b) * (1 - t ** 4 / (12 * b ** 4))); }
    case 'T': return (p.bf * p.tf ** 3 + (p.d - p.tf / 2) * p.tw ** 3) / 3;
    case 'POLY': {
      // a section known only by its outline (an IFC arbitrary profile): an estimate, said in the report. Compact shapes:
      // J ≈ A⁴ / (40 Ip) (Saint-Venant's approximation for solid sections; a circle comes out 1.3 % low); thin-walled ones
      // (A / perimeter² under 0.02): J ≈ A t² / 3 with t = A / (perimeter / 2), the open-section formula on the mid-line
      const s = polySection(p.loops);
      const per = p.loops.reduce((sum, l) => sum + l.reduce((q, pt, i) => q + Math.hypot(l[(i + 1) % l.length][0] - pt[0], l[(i + 1) % l.length][1] - pt[1]), 0), 0);
      if (s.A / (per * per) < 0.02) { const t = s.A / (per / 2); return s.A * t * t / 3; }
      return s.A ** 4 / (40 * (s.Iyy + s.Izz));
    }
    default: return NaN;
  }
}

/**
 * Section properties for the analysis from a section drawn in (u, v) — v the depth (DetailCraft's profiles; a member's
 * local y = u, local z = v). Symmetric sections: Iy = ∫ v² dA (the strong axis of a beam), Iz = ∫ u² dA, roll 0. A section
 * whose u, v are not principal (an angle) is analysed on its principal axes: Iy = I1, Iz = I2 and roll = the turn of the
 * principal axes, so loads in local directions follow those axes (said in the report).
 */
export function sectionOf(p, name = p.name) {
  const s = polySection(p.loops);
  const J = torsionJ(p);
  const A = Number.isFinite(p.area) && p.area > 0 ? p.area : s.A;     // the exact area (arcs exact) where the profile has one
  if (!s.angle) return { name, A, Iy: s.Iyy, Iz: s.Izz, J, roll: 0 };
  // the I1 axis is at `angle` from +u; local y must lie along it so that Iy (bending about y) = I1
  return { name, A, Iy: s.I1, Iz: s.I2, J, roll: s.angle };
}