Overhead conductor sag and tension

export function conductorLoad(c, lc) {
  const D = (c.diameter_mm || 0) / 1000;
  const t = Math.max(0, (lc?.ice_mm || 0) / 1000);
  const rhoI = lc?.iceDensity > 0 ? lc.iceDensity : 900;
  const iceArea = Math.PI * (sq(D / 2 + t) - sq(D / 2));
  const iceN = iceArea * rhoI * G;
  const wv = (c.mass_kg_m || 0) * G + iceN;
  const dWind = D + (lc?.windIce === false ? 0 : 2 * t);
  const wh = Math.max(0, lc?.wind_Pa || 0) * dWind;
  return { wv, wh, w: Math.hypot(wv, wh), iceN, dWind };
}

export function rulingSpan(spans) {
  let s3 = 0, s1 = 0;
  for (const s of spans) { const L = Math.abs(s.L); s3 += L * L * L; s1 += L; }
  return s1 > 0 ? Math.sqrt(s3 / s1) : 0;
}

const EA = (c) => (c.e_GPa || 0) * 1e9 * (c.area_mm2 || 0) * 1e-6;
const ALPHA = (c) => (c.alpha_e6 || 0) * 1e-6;

export function changeOfState(c, L, from, to) {
  const ea = EA(c), al = ALPHA(c);
  if (!(ea > 0) || !(from.H > 0) || !(L > 0)) return NaN;
  const K = ea * (sq(from.w * L) / (24 * sq(from.H)) + al * (to.temp - from.temp)) - from.H;
  const M = ea * sq(to.w * L) / 24;
  // f(H) = H³ + K·H² − M has exactly one positive root when M > 0
  const f = (H) => H * H * H + K * H * H - M;
  let H = Math.max(from.H, 1);
  for (let i = 0; i < 100; i++) {
    const fx = f(H), d = 3 * H * H + 2 * K * H;
    if (!(d > 0)) { H *= 2; continue; }
    const n = H - fx / d;
    if (!(n > 0)) { H /= 2; continue; }
    if (Math.abs(n - H) < 1e-9 * Math.max(1, H)) return n;
    H = n;
  }
  // bisection fallback
  let lo = 1e-6, hi = Math.max(1, from.H);
  while (f(hi) < 0) hi *= 2;
  for (let i = 0; i < 200; i++) { const m = (lo + hi) / 2; if (f(m) > 0) hi = m; else lo = m; }
  return (lo + hi) / 2;
}

export function catenary(span, H, w) {
  const L = Math.abs(span.L), dz = span.dz || 0;
  const a = H / w;
  // the lowest point's x (from the left end): dz = a·(cosh((L − x0)/a) − cosh(x0/a))
  const x0 = L / 2 - a * Math.asinh(dz / (2 * a * Math.sinh(L / (2 * a))));
  const y0 = -a * (Math.cosh(x0 / a) - 1);   // the lowest point's height below the left end
  const yAt = (x) => y0 + a * (Math.cosh((x - x0) / a) - 1);
  const xs = x0 + a * Math.asinh(dz / L);   // where the curve runs parallel to the chord
  const xm = Math.min(L, Math.max(0, xs));
  const sag = (dz * xm / L) - yAt(xm);
  const length = a * (Math.sinh((L - x0) / a) + Math.sinh(x0 / a));
  const Tleft = H * Math.cosh(x0 / a), Tright = H * Math.cosh((L - x0) / a);
  return { L, dz, a, x0, y0, xLow: x0, sag, length, Tleft, Tright, yAt };
}

/** Horizontal tension that gives a largest support tension T in the ruling span (level), by bisection on H. */
function hForT(T, w, Lr) {
  if (!(T > 0)) return NaN;
  let lo = 1e-6, hi = T;
  for (let i = 0; i < 200; i++) { const m = (lo + hi) / 2; if (m * Math.cosh(Lr / (2 * (m / w))) > T) hi = m; else lo = m; }
  return lo;
}

// a T limit is about the worst support of the real spans: scale down until every T limit holds (a few passes)
if (!fixed) {
  for (let pass = 0; pass < 20; pass++) {
    let worst = 1;
    for (const x of cands) {
      if (x.lim.kind !== 'T') continue;
      const N = x.lim.N > 0 ? x.lim.N : x.lim.pct / 100 * (c.rts_N || 0);
      const r = res[x.lc.id].Tmax / N;
      if (r > worst + 1e-9) { worst = r; gov = x; }
    }
    if (worst <= 1 + 1e-9) break;
    res = build(res[first.id].H / worst);
  }
}
// which limits hold
const checks = cands.map((x) => {
  const r = res[x.lc.id];
  const N = x.lim.fixed > 0 ? x.lim.fixed : x.lim.N > 0 ? x.lim.N : x.lim.pct / 100 * (c.rts_N || 0);
  const v = x.lim.kind === 'T' ? r.Tmax : r.H;
  return { caseId: x.lc.id, kind: x.lim.kind || 'H', limit: N, value: v, ok: x.lim.fixed > 0 || v <= N * (1 + 1e-6) };

export function stringingTable(c, spans, ref, temps) {
  const Lr = rulingSpan(spans);
  const w = conductorLoad(c, {}).w;
  return temps.map((temp) => {
    const H = changeOfState(c, Lr, ref, { w, temp });
    const sp = spans.map((s) => catenary(s, H, w));
    return { temp, H, Tmax: Math.max(...sp.map((x) => Math.max(x.Tleft, x.Tright))), length: sp.reduce((t, x) => t + x.length, 0),
      spans: sp.map((x) => ({ L: x.L, dz: x.dz, sag: x.sag, length: x.length, Tleft: x.Tleft, Tright: x.Tright })) };
  });
}

/** Each span on its own from the section's reference: tension with ice vs bare, at the ice case's temperature. */
export function unbalancedIce(c, spans, ref, iceCase, bareCase) {
  const wi = conductorLoad(c, iceCase).w, wb = conductorLoad(c, bareCase || { temp: iceCase.temp }).w;
  const per = spans.map((s) => {
    const Hi = changeOfState(c, Math.abs(s.L), ref, { w: wi, temp: iceCase.temp });
    const Hb = changeOfState(c, Math.abs(s.L), ref, { w: wb, temp: (bareCase || iceCase).temp });
    return { Hi, Hb };
  });
  const out = [];
  for (let i = 1; i < spans.length; i++) {
    const a = per[i - 1], b = per[i];
    out.push({ support: i, dT: Math.max(Math.abs(a.Hi - b.Hb), Math.abs(b.Hi - a.Hb)) });
  }
  return out;
}

function stranded({ key, name, category, wires, d, steelWires = 0, ds = 0, mat = 'al' }) {
  const M = { al: { rho: 2703, res: 0.028264, uts: 170e6, e: 59, a: 23 }, alloy: { rho: 2690, res: 0.0328, uts: 300e6, e: 59, a: 23 }, cu: { rho: 8890, res: 0.017241, uts: 400e6, e: 105, a: 17 } }[mat];
  const aA = wires * PI * d * d / 4, aS = steelWires * PI * ds * ds / 4;   // mm²
  const layers = (n) => (n <= 1 ? 1 : n <= 7 ? 3 : n <= 19 ? 5 : n <= 37 ? 7 : 9);
  const D = steelWires ? ds * layers(steelWires) + 2 * d * ((layers(wires + steelWires) - layers(steelWires)) / 2) : d * layers(wires);
  const mass = 1.02 * (aA * 1e-6 * M.rho + aS * 1e-6 * 7780);
  const rts = 0.9 * (aA * 1e-6 * M.uts + aS * 1e-6 * 1300e6);
  const e = steelWires ? (M.e * aA + 200 * aS) / (aA + aS) : M.e;
  const alpha = steelWires ? (M.a * M.e * aA + 11.5 * 200 * aS) / (M.e * aA + 200 * aS) : M.a;
  const r = 1.02 * M.res / aA * 1000;
  return { key, kind: 'conductor', source: 'starter', name, category, note: t('pl.st.note'),
    material: steelWires ? 'aluminium-steel' : mat === 'alloy' ? 'aluminium-alloy' : mat === 'cu' ? 'copper' : 'aluminium',
    area_mm2: +(aA + aS).toFixed(1), diameter_mm: +D.toFixed(2), mass_kg_m: +mass.toFixed(4), rts_N: Math.round(rts), e_GPa: +e.toFixed(1),
    alpha_e6: +alpha.toFixed(1), r20_ohm_km: +r.toFixed(4) };
}