// Fellwing flight model: a 6-degree-of-freedom rigid body with real aerodynamic coefficients (lift curve with stall
// break, induced drag, stability and damping derivatives, control power from airspeed and propwash), an ISA atmosphere,
// prop or jet thrust that falls with altitude, flaps, retractable gear, and spring-damper landing gear with brakes and
// nose-wheel steering, ground effect (less induced drag and more lift within a wingspan of the ground) and pitch trim.
// Pure maths, no DOM and no three.js, so node tests run it. Fixed 120 Hz steps.
// Body axes: x right, y up, z back (forward is -z). World: y up, -z north, +x east. SI units.
// (c) 2026 CyberMax. All rights reserved.

export const KT = 0.514444; // m/s per knot
export const FT = 0.3048;
const G = 9.80665;
const D2R = Math.PI / 180;

// ---------------------------------------------------------------- tiny vector / quaternion helpers
const v3 = (x = 0, y = 0, z = 0) => ({ x, y, z });
const add = (a, b) => v3(a.x + b.x, a.y + b.y, a.z + b.z);
const sub = (a, b) => v3(a.x - b.x, a.y - b.y, a.z - b.z);
const mul = (a, s) => v3(a.x * s, a.y * s, a.z * s);
const dot = (a, b) => a.x * b.x + a.y * b.y + a.z * b.z;
const cross = (a, b) => v3(a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x);
const len = (a) => Math.hypot(a.x, a.y, a.z);
const norm = (a) => { const l = len(a) || 1; return v3(a.x / l, a.y / l, a.z / l); };
/** Rotate vector v by unit quaternion q {x,y,z,w}. */
export function rot(q, v) {
  const tx = 2 * (q.y * v.z - q.z * v.y), ty = 2 * (q.z * v.x - q.x * v.z), tz = 2 * (q.x * v.y - q.y * v.x);
  return v3(v.x + q.w * tx + (q.y * tz - q.z * ty), v.y + q.w * ty + (q.z * tx - q.x * tz), v.z + q.w * tz + (q.x * ty - q.y * tx));
}
const conj = (q) => ({ x: -q.x, y: -q.y, z: -q.z, w: q.w });
const qmul = (a, b) => ({
  w: a.w * b.w - a.x * b.x - a.y * b.y - a.z * b.z,
  x: a.w * b.x + a.x * b.w + a.y * b.z - a.z * b.y,
  y: a.w * b.y - a.x * b.z + a.y * b.w + a.z * b.x,
  z: a.w * b.z + a.x * b.y - a.y * b.x + a.z * b.w,
});
const qnorm = (q) => { const l = Math.hypot(q.x, q.y, q.z, q.w) || 1; return { x: q.x / l, y: q.y / l, z: q.z / l, w: q.w / l }; };
/** Quaternion from heading (deg true, 0 = north), pitch and bank (deg). */
export function fromEuler(hdg, pitch = 0, bank = 0) {
  const ax = (axis, a) => { const s = Math.sin(a / 2); return { x: axis[0] * s, y: axis[1] * s, z: axis[2] * s, w: Math.cos(a / 2) }; };
  const yaw = ax([0, 1, 0], -hdg * D2R), pit = ax([1, 0, 0], pitch * D2R), rol = ax([0, 0, 1], -bank * D2R);
  return qnorm(qmul(qmul(yaw, pit), rol));
}

// ---------------------------------------------------------------- atmosphere (ISA, troposphere + lower stratosphere)
export function atmosphere(alt) {
  const h = Math.max(-500, Math.min(20000, alt));
  const T = Math.max(216.65, 288.15 - 0.0065 * h);
  const p = h <= 11000 ? 101325 * (T / 288.15) ** 5.25588 : 22632 * Math.exp(-(h - 11000) / 6341.6);
  return { T, p, rho: p / (287.053 * T) };
}
const RHO0 = 1.225;

// ---------------------------------------------------------------- aircraft (fictional types, real-world class data)
// I = [pitch (about x), yaw (about y), roll (about z)] kg m². Gear: [x, y, z, isNose]. Hull: points that must never touch.
export const AIRCRAFT = {
  kestrel: {
    id: 'kestrel', name: 'Kestrel T-3', kind: 'Turbocharged 4-seat trainer', free: true, seats: 4,
    m: 1250, S: 16.2, b: 11.0, c: 1.49, I: [1825, 2667, 1285],
    engine: 'prop', P: 230e3, turbo: 0.5, Tstatic: 3900, eta: 0.8,
    CL0: 0.28, CLa: 5.0, aStall: 15.5, CD0: 0.031, k: 0.056, flapCL: 0.55, flapCD: 0.045, flaps: [0, 0.5, 1], gearCD: 0, retract: false,
    Cma: -0.85, Cmq: -14, Cmde: 1.15, Clb: -0.09, Clp: -0.5, Clr: 0.1, Clda: 0.17, Cnb: 0.07, Cnr: -0.12, Cnda: -0.012, Cndr: 0.075, CYb: -0.4,
    gear: [[0, -1.25, -1.75, 1], [-1.25, -1.3, 0.35, 0], [1.25, -1.3, 0.35, 0]], hull: [[-5.5, 0.35, 0.2], [5.5, 0.35, 0.2], [0, 0.25, 6.2], [0, -0.1, -2.4], [0, 1.5, 5.9]],
    k_gear: 70000, c_gear: 7000, vr: 55, vcruise: 120, vne: 165, eye: [-0.3, 0.55, -0.4], spawnAgl: 1.32, maxG: 5.5,
  },
  rook: {
    id: 'rook', name: 'Rook 9 STOL', kind: 'Turboprop bush plane, short take-off and landing', free: false, seats: 8,
    m: 2000, S: 25.5, b: 15.9, c: 1.62, I: [4200, 7600, 3500],
    engine: 'prop', P: 500e3, turbo: 0.55, Tstatic: 9000, eta: 0.82,
    CL0: 0.34, CLa: 5.1, aStall: 16.5, CD0: 0.038, k: 0.05, flapCL: 0.95, flapCD: 0.07, flaps: [0, 0.4, 0.75, 1], gearCD: 0, retract: false,
    Cma: -0.9, Cmq: -16, Cmde: 1.25, Clb: -0.1, Clp: -0.55, Clr: 0.12, Clda: 0.16, Cnb: 0.075, Cnr: -0.13, Cnda: -0.012, Cndr: 0.08, CYb: -0.45,
    gear: [[0, -1.55, -2.2, 1], [-1.6, -1.6, 0.4, 0], [1.6, -1.6, 0.4, 0]], hull: [[-7.9, 0.6, 0.2], [7.9, 0.6, 0.2], [0, 0.3, 7.2], [0, -0.3, -3.2], [0, 1.9, 7.0]],
    k_gear: 110000, c_gear: 11000, vr: 45, vcruise: 150, vne: 190, eye: [-0.35, 0.7, -0.9], spawnAgl: 1.62, maxG: 4.5,
  },
  strix: {
    id: 'strix', name: 'Strix A2', kind: 'Aerobatic monoplane, +/-10 g', free: false, seats: 1,
    m: 760, S: 9.6, b: 7.4, c: 1.35, I: [950, 1550, 620],
    engine: 'prop', P: 250e3, turbo: 0, Tstatic: 4800, eta: 0.8,
    CL0: 0.0, CLa: 5.2, aStall: 17, CD0: 0.028, k: 0.075, flapCL: 0, flapCD: 0, flaps: [0], gearCD: 0, retract: false,
    Cma: -0.55, Cmq: -10, Cmde: 1.35, Clb: -0.03, Clp: -0.32, Clr: 0.06, Clda: 0.3, Cnb: 0.06, Cnr: -0.1, Cnda: -0.006, Cndr: 0.09, CYb: -0.45,
    gear: [[0, -1.0, -1.4, 1], [-1.1, -1.05, 0.2, 0], [1.1, -1.05, 0.2, 0]], hull: [[-3.7, 0.0, 0.2], [3.7, 0.0, 0.2], [0, 0.25, 5.0], [0, -0.1, -2.3], [0, 1.2, 4.8]],
    k_gear: 50000, c_gear: 5000, vr: 60, vcruise: 165, vne: 220, eye: [0, 0.65, 0.2], spawnAgl: 1.07, maxG: 10,
  },
  vesper: {
    id: 'vesper', name: 'Vesper SJ', kind: 'Light business jet, twin turbofan', free: false, seats: 8,
    m: 7400, S: 30, b: 15.8, c: 2.05, I: [42000, 80000, 38000],
    engine: 'jet', T: 2 * 16500, CL0: 0.22, CLa: 5.3, aStall: 14, CD0: 0.022, k: 0.048, flapCL: 0.85, flapCD: 0.06, flaps: [0, 0.35, 0.7, 1], gearCD: 0.022, retract: true,
    Cma: -1.0, Cmq: -18, Cmde: 1.2, Clb: -0.08, Clp: -0.45, Clr: 0.12, Clda: 0.12, Cnb: 0.09, Cnr: -0.16, Cnda: -0.01, Cndr: 0.08, CYb: -0.55,
    gear: [[0, -1.55, -5.2, 1], [-1.9, -1.6, 0.9, 0], [1.9, -1.6, 0.9, 0]], hull: [[-7.9, 0.7, 1.6], [7.9, 0.7, 1.6], [0, 0.0, 8.5], [0, -0.1, -7.6], [0, 3.9, 8.2]],
    k_gear: 400000, c_gear: 40000, vr: 115, vcruise: 400, vne: 450, eye: [-0.45, 0.85, -5.6], spawnAgl: 1.62, maxG: 3.5,
  },
};

// Pitch trim: each type is rigged to fly hands-off at 85 % of its cruise speed (IAS), clean, like a trimmed aircraft.
for (const A of Object.values(AIRCRAFT)) {
  const v = A.vcruise * 0.85 * KT, CLt = (2 * A.m * G) / (RHO0 * A.S * v * v);
  A.Cm0 = -A.Cma * ((CLt - A.CL0) / A.CLa);
}

// ---------------------------------------------------------------- state
export function createPlane(type, { pos = v3(), hdg = 0, pitch = 0, speed = 0, onGround = false } = {}) {
  const A = AIRCRAFT[type] || AIRCRAFT.kestrel;
  const q = fromEuler(hdg, pitch, 0);
  return {
    A, pos: { ...pos }, q, vel: rot(q, v3(0, 0, -speed)), w: v3(), // w = body angular velocity (rad/s)
    ctl: { pitch: 0, roll: 0, yaw: 0, throttle: onGround ? 0 : 0.72, flaps: 0, brake: onGround ? 1 : 0, gear: true, trim: 0 },
    surf: { e: 0, a: 0, r: 0, f: 0, gear: 1, thr: onGround ? 0 : 0.72 }, // actual deflections (rate-limited)
    onGround, crashed: null, t: 0, g: 1, alpha: 0, beta: 0, stall: 0, ias: speed, tas: speed, buffet: 0, contact: 0, prop: 0,
    maxVz: 0, landed: null,
  };
}

const approach = (cur, tgt, rate, dt) => cur + Math.max(-rate * dt, Math.min(rate * dt, tgt - cur));
const smooth = (a, b, x) => { const t = Math.max(0, Math.min(1, (x - a) / (b - a))); return t * t * (3 - 2 * t); };

/**
 * One fixed step. env = { ground(x, z) -> { h, water }, wind: {x,y,z} m/s world, turb: 0..1, rand: () => 0..1 }.
 * Returns the plane (mutated).
 */
export function step(P, dt, env) {
  if (P.crashed) return P;
  const A = P.A, c = P.ctl, s = P.surf;
  P.t += dt;
  // control surfaces and engine spool (real surfaces and engines are not instant)
  s.e = approach(s.e, Math.max(-1, Math.min(1, c.pitch + c.trim)), 3.2, dt);
  s.a = approach(s.a, c.roll, 4.0, dt);
  s.r = approach(s.r, c.yaw, 3.0, dt);
  s.f = approach(s.f, c.flaps, 0.35, dt);
  s.gear = approach(s.gear, A.retract ? (c.gear || P.onGround ? 1 : 0) : 1, 0.2, dt);
  s.thr = approach(s.thr, c.throttle, A.engine === 'jet' ? 0.35 : 1.6, dt);

  const alt = P.pos.y;
  const { rho } = atmosphere(alt);
  const wind = env.wind || v3();
  const turb = env.turb || 0;
  const gust = turb > 0 ? v3((env.rand() - 0.5) * 6 * turb, (env.rand() - 0.5) * 4 * turb, (env.rand() - 0.5) * 6 * turb) : v3();
  const air = sub(P.vel, add(wind, gust)); // velocity relative to the air mass, world
  const vb = rot(conj(P.q), air); // body frame
  const V = len(air);
  const right = rot(P.q, v3(1, 0, 0));
  let F = v3(0, -A.m * G, 0); // world forces
  let Tq = v3(); // body torques (x = pitch up, y = yaw left, z = roll left)

  // thrust
  let thrust = 0;
  if (A.engine === 'jet') thrust = s.thr * A.T * (rho / RHO0) ** 0.75;
  else {
    const pAvail = A.P * Math.min(1, (rho / RHO0) / Math.max(0.05, 1 - A.turbo));
    thrust = s.thr * Math.min(A.Tstatic, (A.eta * pAvail) / Math.max(V, 8));
  }
  F = add(F, rot(P.q, v3(0, 0, -thrust)));
  P.prop = s.thr;

  if (V > 0.5) {
    const qd = 0.5 * rho * V * V;
    const alpha = Math.atan2(-vb.y, -vb.z);
    const beta = Math.asin(Math.max(-1, Math.min(1, vb.x / V)));
    const aS = (A.aStall - 3 * s.f) * D2R;
    // ground effect (McCormick): within about one span of the ground the wingtip vortices are squashed, induced drag
    // falls and the lift slope rises a little, which is why aircraft float in the flare
    const hag = env.agl !== undefined ? Math.max(0.3, env.agl) : 1e9;
    const hb = hag / A.b, ge = hb < 1.1 ? (16 * hb) ** 2 / (1 + (16 * hb) ** 2) : 1;
    P.ge = 1 - ge;
    const lin = A.CL0 + A.flapCL * s.f + A.CLa * alpha * (1 + 0.1 * (1 - ge));
    const flat = 1.05 * Math.sin(2 * alpha);
    const st = smooth(aS, aS + 0.12, Math.abs(alpha));
    const linC = Math.max(-1.4, Math.min(A.CL0 + A.flapCL * s.f + A.CLa * aS, lin));
    const CL = linC * (1 - st) + flat * st;
    const CD = A.CD0 + A.k * ge * CL * CL * (1 - st) + A.flapCD * s.f * s.f + A.gearCD * s.gear + 1.25 * Math.sin(alpha) ** 2 * st + 0.4 * Math.sin(beta) ** 2;
    const vdir = mul(air, 1 / V);
    const liftDir = norm(cross(right, vdir));
    F = add(F, mul(liftDir, qd * A.S * CL));
    F = add(F, mul(vdir, -qd * A.S * CD));
    F = add(F, mul(right, qd * A.S * A.CYb * beta));
    // moments (conventional p roll-right, q nose-up, r nose-right)
    const p = -P.w.z, q = P.w.x, r = -P.w.y;
    const ph = (p * A.b) / (2 * V), qh = (q * A.c) / (2 * V), rh = (r * A.b) / (2 * V);
    const qc = qd + (0.45 * thrust) / A.S; // propwash over the tail gives control at low speed
    const stallRoll = st * (env.rand ? (env.rand() - 0.5) * 0.12 : 0) + st * 0.02 * Math.sign(beta || 1);
    const L = qd * A.S * A.b * (A.Clb * beta + A.Clp * ph + A.Clr * rh + stallRoll) + qd * A.S * A.b * A.Clda * s.a * 0.3 * (1 - 0.5 * st);
    const M = qd * A.S * A.c * (A.Cm0 - 0.04 * s.f + A.Cma * alpha + A.Cmq * qh) + qc * A.S * A.c * A.Cmde * s.e * 0.35;
    const N = qd * A.S * A.b * (A.Cnb * beta + A.Cnr * rh + A.Cnda * s.a) + qc * A.S * A.b * A.Cndr * s.r * 0.6;
    Tq = add(Tq, v3(M, -N, -L));
    P.alpha = alpha; P.beta = beta; P.stall = st; P.buffet = smooth(aS - 0.05, aS + 0.05, Math.abs(alpha));
  } else { P.alpha = 0; P.beta = 0; P.stall = 0; P.buffet = 0; }

  // landing gear + hull contact
  let onG = false, contacts = 0;
  const wWorld = rot(P.q, P.w);
  for (const g of A.gear) {
    if (s.gear < 0.98) break;
    const rb = v3(g[0], g[1], g[2]);
    const r = rot(P.q, rb);
    const pt = add(P.pos, r);
    const gr = env.ground(pt.x, pt.z);
    const gh = Math.max(gr.h, 0);
    const pen = gh - pt.y;
    if (pen <= 0) continue;
    const vp = add(P.vel, cross(wWorld, r));
    if (gr.h < 0.2 && gr.water) { P.crashed = 'Ditched in the water'; return P; }
    if (vp.y < -6.5 || pen > 1.2) { P.crashed = `Hard landing (${(-vp.y * 196.85).toFixed(0)} ft/min)`; return P; }
    onG = true; contacts++;
    P.maxVz = Math.min(P.maxVz, vp.y);
    const Nf = Math.max(0, A.k_gear * pen - A.c_gear * vp.y);
    let fwd = rot(P.q, v3(0, 0, -1));
    if (g[3]) { // nose wheel steering with rudder, less at speed
      const ang = -c.yaw * 0.45 * (1 - smooth(10, 30, len(P.vel)));
      const ca = Math.cos(ang), sa = Math.sin(ang);
      const rt = rot(P.q, v3(1, 0, 0));
      fwd = add(mul(fwd, ca), mul(rt, -sa));
    }
    fwd = norm(v3(fwd.x, 0, fwd.z));
    const lat = v3(-fwd.z, 0, fwd.x);
    const vl = dot(vp, fwd), vs = dot(vp, lat);
    const mu = g[3] ? 0.025 : 0.025 + c.brake * 0.65;
    const Fl = mul(fwd, -Math.max(-1, Math.min(1, vl * 4)) * mu * Nf);
    const Fs = mul(lat, -Math.max(-1, Math.min(1, vs * 2.5)) * 0.85 * Nf);
    const Fc = add(add(v3(0, Nf, 0), Fl), Fs);
    F = add(F, Fc);
    Tq = add(Tq, rot(conj(P.q), cross(r, Fc)));
  }
  for (const h of A.hull) {
    const pt = add(P.pos, rot(P.q, v3(h[0], h[1], h[2])));
    const gr = env.ground(pt.x, pt.z);
    if (pt.y < Math.max(gr.h, 0) - 0.05) { P.crashed = gr.h < 0.2 && gr.water ? 'Ditched in the water' : (Math.abs(h[0]) > 2 ? 'Wingtip strike' : 'Crashed into the terrain'); return P; }
  }
  if (!onG && s.gear < 0.98) { // gear up: the belly is a contact point
    const pt = add(P.pos, rot(P.q, v3(0, A.gear[1][1] * 0.7, 0)));
    if (pt.y < Math.max(env.ground(pt.x, pt.z).h, 0)) { P.crashed = 'Belly landing: the gear was up'; return P; }
  }
  if (P.onGround && !onG && !P.landed) P.landed = null;
  if (!P.onGround && onG) P.touch = { at: P.t, vz: P.maxVz };
  if (!onG) P.maxVz = 0;
  P.onGround = onG; P.contact = contacts;

  // integrate (semi-implicit Euler)
  const acc = mul(F, 1 / A.m);
  const bodyAcc = rot(conj(P.q), sub(acc, v3(0, -G, 0)));
  P.g = bodyAcc.y / G; // load factor felt by the pilot
  if (Math.abs(P.g) > A.maxG * 1.5 && !onG) { P.crashed = `Structural failure at ${P.g.toFixed(1)} g`; return P; }
  P.vel = add(P.vel, mul(acc, dt));
  P.pos = add(P.pos, mul(P.vel, dt));
  const I = A.I, w = P.w;
  const Iw = v3(I[0] * w.x, I[1] * w.y, I[2] * w.z);
  const gyro = cross(w, Iw);
  P.w = v3(w.x + ((Tq.x - gyro.x) / I[0]) * dt, w.y + ((Tq.y - gyro.y) / I[1]) * dt, w.z + ((Tq.z - gyro.z) / I[2]) * dt);
  if (onG) P.w = mul(P.w, 0.995); // tyre scrub
  const dq = { x: P.w.x * dt * 0.5, y: P.w.y * dt * 0.5, z: P.w.z * dt * 0.5, w: 0 };
  const qd = qmul(P.q, dq);
  P.q = qnorm({ x: P.q.x + qd.x, y: P.q.y + qd.y, z: P.q.z + qd.z, w: P.q.w + qd.w });
  P.tas = V;
  P.ias = V * Math.sqrt(rho / RHO0);
  if (P.ias > A.vne * KT * 1.25) { P.crashed = 'Overspeed: the airframe broke up'; return P; }
  return P;
}

/** Instruments: heading (deg true), pitch, bank (deg), vertical speed (m/s), IAS/TAS (kt), altitude (m), AoA (deg). */
export function readout(P) {
  const f = rot(P.q, v3(0, 0, -1)), up = rot(P.q, v3(0, 1, 0)), rt = rot(P.q, v3(1, 0, 0));
  const hdg = (Math.atan2(f.x, -f.z) / D2R + 360) % 360;
  const pitch = Math.asin(Math.max(-1, Math.min(1, f.y))) / D2R;
  const bank = Math.atan2(-rt.y, up.y) / D2R;
  const gs = Math.hypot(P.vel.x, P.vel.z);
  return { hdg, pitch, bank, vs: P.vel.y, ias: P.ias / KT, tas: P.tas / KT, gs: gs / KT, alt: P.pos.y, aoa: P.alpha / D2R, g: P.g, stall: P.stall > 0.25, buffet: P.buffet };
}

export const vec = { v3, add, sub, mul, dot, cross, len, norm, rot, conj, qmul, qnorm };