Spaces:
Running
Running
| // 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 }; | |