ブラックホールに近づいた星は、なぜ回り続けられないのか

 // 円盤の内側と外側の半径(M)・天体の数の上限・軌跡の点の数・時間の速さ 1 倍のとき 1 秒に進む時間(M)const RIN = 6, ROUT = 16, MAXB = 12, TRAIL = 1600, GUIDE = 1400, RATE = 40;const KM = 1.4766, SEC = 4.9255e-6; // 太陽の質量のときの GM/c²(km)と GM/c³(秒)const MASS: Record<string, { m: number; label: string }> = {  stellar: { m: 10, label: "太陽の 10 倍" },  sgr: { m: 4.0e6, label: "いて座 A*" },  m87: { m: 6.5e9, label: "M87*" },};// 画質ごとの光の積分: 角の刻み h(ラジアン)と、1 画素で進める回数の上限 nconst QUAL = [{ h: 0.1, n: 120 }, { h: 0.07, n: 200 }, { h: 0.05, n: 320 }]; type Vec = [number, number, number];const cross = (a: Vec, b: Vec): Vec => [a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]];const dot = (a: Vec, b: Vec) => a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; interface Orb { r: number; pr: number; phi: number; E: number; L: number; tau: number; t: number }interface Info { kind: "circle" | "bound" | "plunge" | "escape"; rp: number; ra: number; prec: number; wall: boolean }interface Body extends Orb {  id: number; tx: Float32Array; ty: Float32Array; tn: number; ti: number; lx: number; ly: number;  rmin: number; rmax: number; lastPeri: number; meas: number; peri: number[];  info: Info | null; guide: Float32Array; gn: number; guideAt: number;}interface Box { x: number; y: number; w: number; h: number }interface Drag { mode: "cam" | "throw"; lx: number; ly: number; moved: boolean; px: number; py: number; vr: number; vt: number; v: number; snap: boolean; path: Float32Array; pn: number; info: Info | null }interface S {  gl: WebGL2RenderingContext | null; prog: Prog | null; quad: ReturnType<typeof fullscreen> | null; locB: WebGLUniformLocation | null;  bodies: Body[]; nextId: number; sel: Body | null; time: number; camPhi: number;  step: (dt: number, f: (h: number) => void) => number;  bArr: Float32Array; drag: Drag | null; mapR: number; map: Box | null; graph: Box | null;  mt: { cx: number; cy: number; k: number; cp: number; sp: number } | null;} // ---- 数字の表し方const num = (x: number, d: number) => {  if (!Number.isFinite(x)) return "-";  const s = Math.abs(x).toLocaleString("en-US", { minimumFractionDigits: d, maximumFractionDigits: d });  return (x < 0 && Number(s.replace(/,/g, "")) !== 0 ? "−" : "") + s;};const dg = (x: number) => (Math.abs(x) < 10 ? 1 : 0);function realTime(sec: number) {  const a = Math.abs(sec);  if (!Number.isFinite(a)) return "-";  if (a < 1e-3) return `${num(sec * 1e6, dg(sec * 1e6))} マイクロ秒`;  if (a < 1) return `${num(sec * 1e3, dg(sec * 1e3))} ミリ秒`;  if (a < 120) return `${num(sec, dg(sec))} 秒`;  if (a < 7200) return `${num(sec / 60, dg(sec / 60))} 分`;  if (a < 172800) return `${num(sec / 3600, dg(sec / 3600))} 時間`;  if (a < 6.3e7) return `${num(sec / 86400, dg(sec / 86400))} 日`;  return `${num(sec / 3.156e7, 1)} 年`;}function realLen(km: number) {  if (!Number.isFinite(km)) return "-";  if (km < 1e4) return `${num(km, dg(km))} km`;  if (km < 1e8) return `${num(km / 1e4, dg(km / 1e4))} 万 km`;  if (km < 1e12) return `${num(km / 1e8, dg(km / 1e8))} 億 km`;  return `${num(km / 1e12, 1)} 兆 km`;}// 文字の大きさ: 画面の幅に比例させ、大きな画面でも小さな画面でも読める範囲に収めるconst fs = (api: Api, n: number) => clamp(api.u(n), n * 0.95, n * 1.3);

// 天体の動きは、エネルギー E と角運動量 L(どちらも質量 1 あたり)が保存されることから決まる// E² の形の有効ポテンシャル: (dr/dτ)² = E² − V²(r) なので、E² が V² より小さい所には行けないconst V2 = (r: number, L: number) => (1 - 2 / r) * (1 + (L * L) / (r * r));// 動径の向きの加速度 d²r/dτ²。1 項目がニュートンの重力、2 項目が遠心力の項、3 項目が相対論で加わる引力(近いほど急に強くなる)const accR = (r: number, L2: number) => { const r2 = r * r; return -1 / r2 + L2 / (r2 * r) - (3 * L2) / (r2 * r2); };// 遠くの時計の進み dt/dτ = E/(1 − 2/r)。地平線に近づくほど大きくなる(遠くから見ると、天体の時計はゆっくり進む)const lapseInv = (E: number, r: number) => E / Math.max(1e-4, 1 - 2 / r);/** 4 次のルンゲ・クッタ法で、固有時を h だけ進める(r・dr/dτ・φ・t を同時に) */function rk4(o: Orb, h: number) {  const L = o.L, L2 = L * L, r = o.r, p = o.pr;  const a1 = accR(r, L2), f1 = L / (r * r);  const r2 = r + 0.5 * h * p, p2 = p + 0.5 * h * a1, a2 = accR(r2, L2), f2 = L / (r2 * r2);  const r3 = r + 0.5 * h * p2, p3 = p + 0.5 * h * a2, a3 = accR(r3, L2), f3 = L / (r3 * r3);  const r4 = r + h * p3, p4 = p + h * a3, a4 = accR(r4, L2), f4 = L / (r4 * r4);  o.t += (h / 6) * (lapseInv(o.E, r) + 2 * lapseInv(o.E, r2) + 2 * lapseInv(o.E, r3) + lapseInv(o.E, r4));  o.r = r + (h / 6) * (p + 2 * p2 + 2 * p3 + p4);  o.pr = p + (h / 6) * (a1 + 2 * a2 + 2 * a3 + a4);  o.phi += (h / 6) * (f1 + 2 * f2 + 2 * f3 + f4);  o.tau += h;}// 固有時の刻み: 円軌道の 1 周(2π r^1.5)のおよそ 500 分の 1。中心に近いほど細かくなるconst hOf = (r: number) => clamp(0.012 * r * Math.sqrt(r), 0.002, 0.8);/** 真上から見た図の軌跡に点を足す(前の点から少し離れたときだけ) */ function trail(b: Body) {  const x = b.r * Math.cos(b.phi), y = b.r * Math.sin(b.phi);  if (b.tn && Math.hypot(x - b.lx, y - b.ly) < 0.12 + 0.012 * b.r) return;  b.tx[b.ti] = x; b.ty[b.ti] = y; b.ti = (b.ti + 1) % TRAIL; b.tn = Math.min(TRAIL, b.tn + 1); b.lx = x; b.ly = y;}/** 遠くの時計で T まで進める。0 = 回っている・1 = 地平線へ落ちた・2 = 遠くへ飛び去った */function advance(b: Body, T: number) {  for (let n = 0; n < 4000 && T - b.t > 1e-9; n++) {    let h = hOf(b.r);    const rate = lapseInv(b.E, b.r);    if (h * rate > T - b.t) h = (T - b.t) / rate;    const prOld = b.pr, phiOld = b.phi;    rk4(b, h);    if (!(b.r > 2.03) || !Number.isFinite(b.r)) return 1;    if (b.r > 400) return 2;    b.rmin = Math.min(b.rmin, b.r); b.rmax = Math.max(b.rmax, b.r);    if (prOld < 0 && b.pr >= 0) {      // 近点を通った: 前の近点からの角 − 360° が近点の移動      const ph = phiOld + (b.phi - phiOld) * (-prOld / (b.pr - prOld));      if (b.rmax - b.rmin > 0.02 * b.r) {        if (Number.isFinite(b.lastPeri)) b.meas = ((Math.abs(ph - b.lastPeri) - TAU) * 180) / Math.PI;        b.lastPeri = ph; b.peri.push(ph); if (b.peri.length > 12) b.peri.shift();      }      b.rmin = b.r; b.rmax = b.r;    }    trail(b);  }  return 0;}

/** 軌道の形を、E と L だけから式で求める。E² = V²(r) を u = 1/r で書くと 2L²u³ − L²u² + 2u − (1 − E²) = 0 の 3 次式になり、 *  その解が遠点・近点・山の内側の点。L² < 12(L < √12 M)だと V² に山も谷も無く、天体は地平線へ落ちる */function classify(E: number, L: number, r: number, pr: number): Info {  const L2 = L * L, E2 = E * E;  const f = (u: number) => E2 - (1 - 2 * u) * (1 + L2 * u * u);  const root = (a: number, b: number) => { let fa = f(a); for (let k = 0; k < 60; k++) { const m = 0.5 * (a + b), fm = f(m); if (fm > 0 === fa > 0) { a = m; fa = fm; } else b = m; } return 0.5 * (a + b); };  const out: Info = { kind: "plunge", rp: NaN, ra: NaN, prec: NaN, wall: false };  const peak = L2 > 12, um = peak ? (1 + Math.sqrt(1 - 12 / L2)) / 6 : 0.5;  out.wall = peak && f(um) < 0; // ポテンシャルの山が、天体のエネルギーより高い  const u = 1 / r;  if (u >= um || !out.wall) { if (pr >= 0 && E >= 1 && !(out.wall && u >= um)) out.kind = "escape"; return out; }  if (E >= 1) { out.kind = "escape"; out.rp = 1 / root(Math.min(u, um * 0.999), um); return out; }  const ul = (1 - Math.sqrt(1 - 12 / L2)) / 6; // 谷の底(安定な円軌道の u)。山の頂上は um(不安定な円軌道)  if (f(ul) <= 1e-12) { out.kind = "circle"; out.rp = out.ra = 1 / ul; out.prec = (TAU / Math.sqrt(1 - 6 * ul) - TAU) * 180 / Math.PI; return out; }  const ua = root(0, ul), up = root(ul, um);  out.rp = 1 / up; out.ra = 1 / ua;  out.kind = (out.ra - out.rp) / (out.ra + out.rp) < 0.004 ? "circle" : "bound";  // 近点から遠点まで回る角: dφ/du = 1/√(2(u − ua)(up − u)(u3 − u))。u = ua + (up − ua)(1 − cosθ)/2 と置くと端の発散が消え、  // θ で 0〜π を積分するだけになる(シンプソンの公式・48 区間)。往復で 2 倍し、360° を引いたものが 1 周の近点の移動  const u3 = 0.5 - ua - up, N = 48;  let sum = 0;  for (let k = 0; k <= N; k++) { const th = (Math.PI * k) / N, uu = ua + ((up - ua) * (1 - Math.cos(th))) / 2; sum += (k === 0 || k === N ? 1 : k % 2 ? 4 : 2) / Math.sqrt(2 * (u3 - uu)); }  out.prec = (((2 * sum * Math.PI) / (3 * N) - TAU) * 180) / Math.PI;  return out;}/** これからの道筋(固有時で進める)。maxAng は回る角の上限 */function predict(E: number, L: number, r: number, pr: number, phi: number, out: Float32Array, rLim: number, maxAng: number) {  const o: Orb = { r, pr, phi, E, L, tau: 0, t: 0 };  let n = 0, lx = r * Math.cos(phi), ly = r * Math.sin(phi);  out[0] = lx; out[1] = ly; n = 1;  const cap = out.length >> 1;  for (let k = 0; k < 5000 && n < cap; k++) {    rk4(o, hOf(o.r) * 1.25);    const end = !(o.r > 2.02) || o.r > rLim || Math.abs(o.phi - phi) > maxAng || o.t > 6000;    const rr = Math.max(o.r, 2), x = rr * Math.cos(o.phi), y = rr * Math.sin(o.phi);    if (end || Math.hypot(x - lx, y - ly) > 0.2 + 0.012 * rr) { out[2 * n] = x; out[2 * n + 1] = y; n++; lx = x; ly = y; }    if (end) break;  }  return n;}

// ---- カメラ(止まった観測者)。i = 傾き(0 で真上・90° で真横)、camPhi = 赤道面の中のカメラの向き// 座標: 赤道面が x-z 面、円盤の軸が y。C = カメラの位置、F = 見る向き(ブラックホールの中心へ)、U = 画面の上、R = 画面の右function cam(api: Api, s: S) {  const i = (clamp(api.num("incl"), 0.5, 89.5) * Math.PI) / 180, D = api.num("dist"), f = s.camPhi;  const si = Math.sin(i), ci = Math.cos(i), cf = Math.cos(f), sf = Math.sin(f);  const C: Vec = [D * si * cf, D * ci, -D * si * sf], F: Vec = [-si * cf, -ci, si * sf], U: Vec = [-ci * cf, si, ci * sf];  return { C, F, U, R: cross(F, U), D };}/** 画角: 画面の横幅の半分に当たる tan。横長なら縦の半分を 0.4(約 22°)、縦長なら横の半分を 0.5 前後にして、どの形の画面でも真ん中に収める */function tanOf(api: Api) {  const a = api.w / Math.max(1, api.h);  return a >= 1 ? 0.4 * a : 0.4 + 0.1 * clamp((1 - a) / 0.4, 0, 1);}// 天体の見た目の半径(M)。本物の星よりずっと大きく描く(遠くから見ても見えるように距離で少し大きくする)const bodyR = (api: Api) => 0.45 + 0.008 * api.num("dist"); const FS = `#version 300 esprecision highp float;in vec2 vUv;out vec4 outColor;${DITHER}uniform vec2 uRes;uniform vec3 uCam;uniform vec3 uF;uniform vec3 uR;uniform vec3 uU;uniform float uTan, uTime, uDisk, uGrid, uH, uMax, uNB, uBR, uBRmax, uPix, uPixCss;uniform vec4 uB[12];const float RIN = 6.0;const float ROUT = 16.0;const float PI = 3.14159265; // 光子の軌道の式 u'' = 3u² − u(M = 1)void rk4(inout float u, inout float w, float h) {  float k1u = w, k1w = 3.0 * u * u - u;  float u2 = u + 0.5 * h * k1u, w2 = w + 0.5 * h * k1w;  float k2u = w2, k2w = 3.0 * u2 * u2 - u2;  float u3 = u + 0.5 * h * k2u, w3 = w + 0.5 * h * k2w;  float k3u = w3, k3w = 3.0 * u3 * u3 - u3;  float u4 = u + h * k3u, w4 = w + h * k3w;  float k4u = w4, k4w = 3.0 * u4 * u4 - u4;  u += h / 6.0 * (k1u + 2.0 * k2u + 2.0 * k3u + k4u);  w += h / 6.0 * (k1w + 2.0 * k2w + 2.0 * k3w + k4w);}// 刻みの両端の値と傾きから、途中(割合 f)の値を 3 次式で補間するfloat herm(float u0, float w0, float u1, float w1, float h, float f) {  float f2 = f * f, f3 = f2 * f;  return (2.0 * f3 - 3.0 * f2 + 1.0) * u0 + (f3 - 2.0 * f2 + f) * h * w0 + (3.0 * f2 - 2.0 * f3) * u1 + (f3 - f2) * h * w1;}// 線分 P0→P1 と球(中心 c・半径 R)の最初の交点を、線分の割合で返す(交わらなければ −1)float segSphere(vec3 P0, vec3 P1, vec3 c, float R) {  vec3 dv = P1 - P0; float lg = length(dv);  if (lg < 1e-6) return -1.0;  dv /= lg;  vec3 oc = P0 - c;  float bq = dot(oc, dv), cq = dot(oc, oc) - R * R, disc = bq * bq - cq;  if (disc < 0.0) return -1.0;  float sq = -bq - sqrt(disc);  if (sq < 0.0) sq = cq < 0.0 ? 0.0 : -1.0;  if (sq < 0.0 || sq > lg) return -1.0;  return sq / lg;}// 星空と円盤の関数(下の「星空と円盤」で定義する)float sky(vec3 d, float pw);float disk(vec3 Q, float r, float lam, float D, out float alpha);void main() {  // 1. この画素を通る視線の向き d  vec2 p = (gl_FragCoord.xy - 0.5 * uRes) / (0.5 * uRes.x);  vec3 d = normalize(uF + uTan * (p.x * uR + p.y * uU));  // 2. 光の通る平面: カメラの向き e1 と、それに直角で視線を含む向き e2。光の位置は r(e1 cosφ + e2 sinφ)  float D = length(uCam);  vec3 e1 = uCam / D;  float dr = dot(d, e1);  vec3 tv = d - dr * e1;  float dt = length(tv);  vec3 e2 = dt > 1e-6 ? tv / dt : normalize(cross(e1, vec3(0.0, 0.0, 1.0)));  dt = max(dt, 1e-6);  // 3. 始めの値: u = 1/D と、u' = −u√(1 − 2u)·cotψ(cotψ = dr/dt)。b は衝突パラメータ(光の角運動量 ÷ エネルギー)  float u = 1.0 / D;  float w = -u * sqrt(1.0 - 2.0 * u) * dr / dt;  float b = inversesqrt(w * w + u * u * (1.0 - 2.0 * u));  // λ: 光の角運動量の、円盤の軸の向きの成分(ドップラー効果に使う)。実際の光は視線と逆向きに進むので符号を返す  float lam = -b * cross(e1, e2).y;  // 赤道面からの高さは r(a1 cosφ + a2 sinφ)。この符号が変わった所で、光は赤道面を横切る  float a1 = e1.y, a2 = e2.y;  float phi = 0.0, s0 = a1, acc = 0.0, T = 1.0;  int fate = 0; // 0 = 決まらない(光子球のそばを回り続けた)・1 = 遠くへ届いた・2 = 地平線に落ちた・3 = 不透明な物に当たった  vec3 dirOut = vec3(0.0);  // 4. 角 φ を刻み h ずつ進める。遠い所(u が小さい)では刻みを 2.6 倍に広げて回数を減らす  for (int i = 0; i < 640; i++) {    if (float(i) >= uMax) break;    float h = uH * mix(2.6, 1.0, smoothstep(0.03, 0.16, u));    float u1 = u, w1 = w;    rk4(u1, w1, h);    // u が 0 を下回った = 無限遠に届いた。u = 0 になる角を直線で求め、その向きの星空を見る    if (u1 <= 0.0) {      float pf = phi + h * u / (u - u1);      dirOut = e1 * cos(pf) + e2 * sin(pf);      fate = 1;      break;    }    // この刻みの始めと終わりの位置 P0・P1(3 次元)。天体はこの線分と球が交わるかで調べる    float phi1 = phi + h, c0 = cos(phi), n0 = sin(phi), c1 = cos(phi1), n1 = sin(phi1);    float s1 = a1 * c1 + a2 * n1, r0 = 1.0 / u, r1 = 1.0 / u1;    vec3 P0 = r0 * (e1 * c0 + e2 * n0), P1 = r1 * (e1 * c1 + e2 * n1);    float fB = 2.0; int kB = -1;    if (uNB > 0.5 && min(r0, r1) < uBRmax) {      float z0 = r0 * s0, z1 = r1 * s1;      if (abs(z0) < uBR || abs(z1) < uBR || z0 * z1 < 0.0) {        for (int k = 0; k < 12; k++) {          if (float(k) >= uNB) break;          vec4 B = uB[k];          float f = segSphere(P0, P1, vec3(B.x, 0.0, -B.y), B.z);          if (f >= 0.0 && f < fB) { fB = f; kB = k; }        }      }    }    // 赤道面を横切ったか。横切った位置の半径 rc は、刻みの両端の u と u' からエルミート補間で求める    float fD = s0 * s1 < 0.0 ? s0 / (s0 - s1) : 2.0;    if (fD < fB) {      float rc = 1.0 / herm(u, w, u1, w1, h, fD), pc = phi + h * fD;      vec3 Q = rc * (e1 * cos(pc) + e2 * sin(pc));      if (uDisk > 0.5 && rc > RIN && rc < ROUT) { float al; float v = disk(Q, rc, lam, D, al); acc += T * al * v; T *= 1.0 - al; }      // 目印: 赤道面の 6M(ISCO)の円と、3M(光子球)の破線の円      if (uGrid > 0.5) {        float wm = max(length(Q - uCam) * uPix * 1.1, 0.02);        float ang = atan(-Q.z, Q.x);        float dash = step(0.4, fract(ang * 36.0 / (2.0 * PI)));        float ring = max(1.0 - smoothstep(0.5 * wm, wm, abs(rc - 6.0)), (1.0 - smoothstep(0.5 * wm, wm, abs(rc - 3.0))) * dash) * 0.95;        acc += T * ring; T *= 1.0 - ring;      }      if (T < 0.03) { fate = 3; break; }    }    // 天体に当たった: 中心のほうから照らされた球として濃淡を付け、縁取りをする(選んでいる天体は縞模様)    if (kB >= 0) {      vec4 B = uB[kB];      vec3 c = vec3(B.x, 0.0, -B.y), Hp = mix(P0, P1, fB), N = normalize(Hp - c), rd = normalize(P1 - P0);      float lit = max(dot(N, normalize(-c)), 0.0);      float rim = 1.0 - abs(dot(N, rd));      float vb = 0.08 + 0.92 * pow(lit, 0.7);      if (B.w > 0.5) vb = max(vb * 0.7, step(0.5, fract((N.y * 0.5 + 0.5) * 5.0)) * 0.95);      vb = max(vb, smoothstep(0.8, 0.96, rim));      acc += T * vb; T = 0.0; fate = 3;      break;    }    if (u1 > 0.5) { fate = 2; break; }    u = u1; w = w1; phi = phi1; s0 = s1;  }  // 5. 遠くへ届いた光は、まだ透けている分(T)だけ星空を足す。pw は隣の画素との向きの差(方眼の線の太さに使う)  float pw = length(fwidth(dirOut));  float v = acc;  if (fate == 1 && T > 0.02) v += T * sky(dirOut, pw);  // 6. 0〜1 の濃さを、白と黒の網点にする  outColor = vec4(mono(clamp(v, 0.0, 1.0), gl_FragCoord.xy), 1.0);}

float h12(vec2 p) { vec3 q = fract(vec3(p.xyx) * 0.1031); q += dot(q, q.yzx + 33.33); return fract((q.x + q.y) * q.z); }float h13(vec3 p) { p = fract(p * 0.1031); p += dot(p, p.zyx + 31.32); return fract((p.x + p.y) * p.z); }vec3 h33(vec3 p) { p = fract(p * vec3(0.1031, 0.1030, 0.0973)); p += dot(p, p.yxz + 33.33); return fract((p.xxy + p.yxx) * p.zyx); }// なめらかな雲のような模様(格子点ごとに向きをばらばらにした勾配のノイズ)。0〜1 のあたりの値vec2 g2(vec2 p) { float a = h12(p) * 6.2831853; return vec2(cos(a), sin(a)); }float gn(vec2 p) {  vec2 i = floor(p); vec2 f = fract(p); vec2 w = f * f * f * (f * (f * 6.0 - 15.0) + 10.0);  float a = dot(g2(i), f), b = dot(g2(i + vec2(1.0, 0.0)), f - vec2(1.0, 0.0));  float c = dot(g2(i + vec2(0.0, 1.0)), f - vec2(0.0, 1.0)), d = dot(g2(i + vec2(1.0, 1.0)), f - vec2(1.0, 1.0));  return 0.5 + 0.75 * mix(mix(a, b, w.x), mix(c, d, w.x), w.y);}// 星: 向き d を N 倍した 3 次元の格子のます目ごとに、確率 prob で星を 1 つ置く(近くの 8 ますを調べる)float stars(vec3 d, float N, float prob, float big) {  vec3 q = d * N, b0 = floor(q - 0.5);  float v = 0.0;  for (int k = 0; k < 8; k++) {    vec3 c = b0 + vec3(float(k & 1), float((k >> 1) & 1), float((k >> 2) & 1));    float h = h13(c);    if (h > prob) continue;    vec3 sp = normalize(c + 0.2 + 0.6 * h33(c));    float ang = length(d - sp), m = h / prob, br = pow(1.0 - m, 6.0);    float rad = max(uPixCss * (0.7 + big * br), uPix * 0.7);    float core = 1.0 - smoothstep(0.4 * rad, rad, ang);    float halo = (1.0 - smoothstep(rad, 4.5 * rad, ang)) * 0.2 * big * br;    v = max(v, max(core * (0.5 + 0.5 * (1.0 - m)), halo));  }  return v;}// 天球: 3 段の大きさの星と、15° ごとの方眼(赤道は太く)float sky(vec3 d, float pw) {  float v = max(max(stars(d, 7.0, 0.32, 2.4), stars(d, 19.0, 0.5, 0.5)), stars(d, 34.0, 0.35, 0.0) * 0.6);  if (uGrid > 0.5) {    float lat = asin(clamp(d.y, -1.0, 1.0)), lon = atan(d.z, d.x), st = PI / 12.0;    // 経線と緯線を 15° ごとに。線の太さは 1 画素ぶんの角の大きさ pw に合わせる(レンズで引き伸ばされた所は太くなりすぎないように上限)    float w = clamp(pw, uPix * 0.8, uPix * 2.5);    float dl = abs(fract(lat / st + 0.5) - 0.5) * st;    float dn = (abs(fract(lon / st + 0.5) - 0.5) * st * cos(lat)) / max(0.001, smoothstep(1.35, 1.1, abs(lat)));    float line = 1.0 - smoothstep(0.35 * w, w, min(dl, dn));    float eq = 1.0 - smoothstep(0.6 * w, 1.3 * w, abs(lat));    v = max(v, max(0.5 * line, 0.92 * eq));  }  return v;}// 円盤の明るさ(0〜1 の濃さ)と不透明さ alpha。放射 F に g⁴ を掛ける(g: 重力とドップラー効果による光の振動数の変化)float disk(vec3 Q, float r, float lam, float D, out float alpha) {  float x = RIN / r;  float F = (1.0 - sqrt(x)) * x * x * x / 0.05665; // r = 8.2M で最大 1 になるように割る  float Om = pow(r, -1.5);  float g = sqrt(1.0 - 3.0 / r) / (sqrt(1.0 - 2.0 / D) * (1.0 - Om * lam));  float g2 = g * g, I = F * g2 * g2;  // 模様: 半径 8% ごとの輪に分け、輪ごとにその半径の角速度 Ω で回す(内側ほど速い)  float kb = floor(log(r / RIN) / log(1.08));  float r0 = RIN * pow(1.08, kb), r1 = r0 * 1.08, rb = sqrt(r0 * r1), fb = (r - r0) / (r1 - r0);  float ph = atan(-Q.z, Q.x) - pow(rb, -1.5) * uTime;  vec2 q = vec2(cos(ph), sin(ph)) * rb * 1.2 + vec2(kb * 7.31, kb * 3.17);  float cl = clamp(0.5 * gn(q) + 0.32 * gn(q * 2.3 + 5.2) + 0.18 * gn(q * 5.3 + 1.7), 0.0, 1.0);  float groove = 0.5 + 0.5 * cos(fb * 4.0 * PI);  float edge = smoothstep(0.0, 0.07, fb) * smoothstep(1.0, 0.93, fb);  float outer = smoothstep(ROUT, ROUT - 4.0, r) * smoothstep(RIN, RIN + 0.2, r);  alpha = clamp((0.45 + 0.55 * smoothstep(0.25, 0.7, cl)) * mix(0.5, 1.0, edge) * outer, 0.0, 1.0);  float tex = (0.45 + 0.8 * cl) * (0.82 + 0.18 * groove);  return 1.0 - exp(-1.25 * I * tex); // 明るさを 0〜1 の濃さに(明るい所が飽和しすぎないように)}

/** 画面の配置: 真ん中はブラックホール。右下に真上から見た図、左上に有効ポテンシャルのグラフ(どの縦横比でも角に置く) */function layout(api: Api) {  const m = clamp(api.u(12), 8, 18);  const size = Math.round(clamp(Math.min(api.w, api.h) * 0.36, 110, 320));  const s2 = Math.min(size, Math.round(api.w * 0.46));  const map: Box = { x: api.w - s2 - m, y: api.h - s2 - m, w: s2, h: s2 };  const gw = Math.round(clamp(Math.min(api.w * 0.3, api.w - s2 - 3 * m), 150, 380));  const graph: Box = { x: m, y: m, w: gw, h: Math.round(clamp(gw * 0.56, 96, 210)) };  return { map, graph };}function project(api: Api, c: ReturnType<typeof cam>, P: Vec): [number, number] | null {  const d: Vec = [P[0] - c.C[0], P[1] - c.C[1], P[2] - c.C[2]], z = dot(d, c.F);  if (z <= 0) return null;  const t = tanOf(api);  return [api.w / 2 + ((dot(d, c.R) / z / t) * api.w) / 2, api.h / 2 - ((dot(d, c.U) / z / t) * api.w) / 2];} /** 文字の札が、右下の図や左上のグラフに重ならないときだけ描く(幅は全角 1 字 = 文字の大きさ、半角はその 0.55 倍で見積もる) */function labelIfFree(api: Api, s: S, t: string, x: number, y: number, size: number, align: "left" | "center" | "right") {  let w = 0;  for (const ch of t) w += ch.charCodeAt(0) < 0x2000 ? size * 0.55 : size;  const x0 = align === "center" ? x - w / 2 : align === "right" ? x - w : x, y0 = y - size, pad = size * 0.4;  for (const B of [s.map, s.graph]) if (B && x0 - pad < B.x + B.w && x0 + w + pad > B.x && y0 - pad < B.y + B.h && y + pad > B.y) return;  if (x0 < 0 || x0 + w > api.w) return;  api.label(t, x, y, { size, align });}function drawMain(api: Api, s: S) {  const c = cam(api, s), f1 = fs(api, 10.5), cx = api.w / 2, cy = api.h / 2;  if (api.bool("grid")) {    // 影の縁(b = √27 M の光が見える向き)と、光がまっすぐ進むとしたときの地平線の見かけの大きさ    const D = c.D, t = tanOf(api), as = Math.asin(Math.min(1, (Math.sqrt(27) * Math.sqrt(1 - 2 / D)) / D)), ah = Math.asin(2 / D);    const rs = ((Math.tan(as) / t) * api.w) / 2, rh = ((Math.tan(ah) / t) * api.w) / 2;    api.circle(cx, cy, rs, { w: api.u(1), dash: [api.u(5), api.u(4)] });    api.circle(cx, cy, rh, { w: api.u(1), dash: [api.u(1.5), api.u(3)] });    api.line(cx, cy - rs, cx, cy - rs - api.u(14), { w: api.u(0.8) });    labelIfFree(api, s, "影の縁  b = √27 M", cx, cy - rs - api.u(14) - f1 * 0.4, f1, "center");    labelIfFree(api, s, "光が曲がらない場合の地平線", cx, cy + rh + f1 * 1.5, f1, "center");  }  if (api.bool("disk") && api.num("incl") >= 50) {    // 円盤の左右の端(画面の右 = R の向き)。左が近づく側(ドップラー効果で明るい)    const rr = ROUT + 1.2, Rp: Vec = [c.R[0], 0, c.R[2]], l = Math.hypot(Rp[0], Rp[2]) || 1;    const L = project(api, c, [(-Rp[0] / l) * rr, 0, (-Rp[2] / l) * rr]), Rr = project(api, c, [(Rp[0] / l) * rr, 0, (Rp[2] / l) * rr]);    if (L) labelIfFree(api, s, "近づく側", L[0], L[1] + f1 * 0.35, f1, "right");    if (Rr) labelIfFree(api, s, "遠ざかる側", Rr[0], Rr[1] + f1 * 0.35, f1, "left");  }} function drawMap(api: Api, s: S, B: Box) {  const ctx = api.ctx, f1 = fs(api, 9.5), pad = api.u(10);  const k = (B.w / 2 - pad) / s.mapR, cx = B.x + B.w / 2, cy = B.y + B.h / 2 + f1 * 0.5;  const cp = Math.cos(s.camPhi), sp = Math.sin(s.camPhi);  s.mt = { cx, cy, k, cp, sp };  const P = (x: number, y: number): [number, number] => [cx + k * (-x * sp + y * cp), cy + k * (x * cp + y * sp)];  api.rect(B.x, B.y, B.w, B.h, { fill: "paper", w: api.u(1) });  ctx.save(); ctx.beginPath(); ctx.rect(B.x, B.y, B.w, B.h); ctx.clip();  if (api.bool("disk")) api.circle(cx, cy, ROUT * k, { w: api.u(0.6), dash: [api.u(1), api.u(3)] });  api.circle(cx, cy, 6 * k, { w: api.u(1), dash: [api.u(4), api.u(3)] });  api.circle(cx, cy, 3 * k, { w: api.u(0.9), dash: [api.u(1.2), api.u(2.2)] });  const hz: [number, number][] = [];  for (let i = 0; i < 28; i++) hz.push([cx + 2 * k * Math.cos((i / 28) * TAU), cy + 2 * k * Math.sin((i / 28) * TAU)]);  api.hatch(hz, { gap: api.u(2.4), w: api.u(0.7) }); api.circle(cx, cy, 2 * k, { w: api.u(1.2) });  const sel = s.sel;  // 選んだ天体のこれから(点線)と、同じ近点・遠点のニュートンの楕円  if (sel && sel.gn > 1) {    const pts: [number, number][] = [];    for (let i = 0; i < sel.gn; i++) pts.push(P(sel.guide[2 * i], sel.guide[2 * i + 1]));    api.poly(pts, { w: api.u(0.9), dash: [api.u(1), api.u(3)] });  }  const inf = sel?.info;  if (sel && inf && inf.kind === "bound" && Number.isFinite(sel.lastPeri)) {    const a = (inf.rp + inf.ra) / 2, e = (inf.ra - inf.rp) / (inf.ra + inf.rp), w0 = sel.lastPeri, pts: [number, number][] = [];    if (e > 0.08) {      for (let i = 0; i <= 120; i++) { const th = (i / 120) * TAU, rr = (a * (1 - e * e)) / (1 + e * Math.cos(th)); pts.push(P(rr * Math.cos(th + w0), rr * Math.sin(th + w0))); }      api.poly(pts, { w: api.u(0.7), dash: [api.u(5), api.u(2), api.u(1), api.u(2)] });    }  }  // 軌跡(古いほど細く)  for (const b of s.bodies) {    const n = b.tn;    if (n < 2) continue;    const start = (b.ti - n + TRAIL) % TRAIL, seg = 5, wm = api.u(b === sel ? 1.7 : 1.1);    ctx.strokeStyle = api.ink; ctx.setLineDash([]); ctx.lineJoin = "round";    for (let g = 0; g < seg; g++) {      const k0 = Math.floor((g * (n - 1)) / seg), k1 = Math.floor(((g + 1) * (n - 1)) / seg);      ctx.beginPath();      for (let kk = k0; kk <= k1; kk++) { const i = (start + kk) % TRAIL, [X, Y] = P(b.tx[i], b.ty[i]); if (kk === k0) ctx.moveTo(X, Y); else ctx.lineTo(X, Y); }      if (g === seg - 1) { const [X, Y] = P(b.r * Math.cos(b.phi), b.r * Math.sin(b.phi)); ctx.lineTo(X, Y); }      ctx.lineWidth = wm * (0.2 + (0.8 * (g + 1)) / seg); ctx.stroke();    }  }  // 近点の印  if (sel) for (let i = 0; i < sel.peri.length; i++) {    const ph = sel.peri[i], rp = sel.info?.rp ?? sel.rmin, r0 = Number.isFinite(rp) ? rp : sel.r;    const [X0, Y0] = P((r0 - 1.4) * Math.cos(ph), (r0 - 1.4) * Math.sin(ph)), [X1, Y1] = P((r0 + 1.4) * Math.cos(ph), (r0 + 1.4) * Math.sin(ph));    api.line(X0, Y0, X1, Y1, { w: api.u(i === sel.peri.length - 1 ? 1.6 : 0.8) });  }  // 天体  for (const b of s.bodies) {    const [X, Y] = P(b.r * Math.cos(b.phi), b.r * Math.sin(b.phi));    api.circle(X, Y, api.u(3.2), { fill: "ink", stroke: false });    if (b === sel) { api.circle(X, Y, api.u(7), { w: api.u(1), dash: [api.u(2), api.u(2.5)] }); api.label(`天体 ${b.id}`, X + api.u(9), Y - api.u(6), { size: f1 }); }    else api.text(String(b.id), X + api.u(5), Y - api.u(4), { size: f1 * 0.9, font: "mono" });  }  // 投げる前の矢印と予測  const d = s.drag;  if (d && d.mode === "throw") {    const [X0, Y0] = P(d.px, d.py);    if (d.moved && d.pn > 1) {      const pts: [number, number][] = [];      for (let i = 0; i < d.pn; i++) pts.push(P(d.path[2 * i], d.path[2 * i + 1]));      api.poly(pts, { w: api.u(1.2), dash: [api.u(2), api.u(3)] });    }    api.circle(X0, Y0, api.u(4), { fill: "ink", stroke: false });    if (d.moved) {      api.arrow(X0, Y0, api.pointer.x, api.pointer.y, { w: api.u(1.4), head: api.u(7) });      const kind = !d.info ? "" : d.snap ? "円軌道" : d.info.kind === "plunge" ? "落ちる" : d.info.kind === "escape" ? "飛び去る" : d.info.kind === "circle" ? "円に近い" : "回り続ける";      const tx = clamp(api.pointer.x + api.u(8), B.x + api.u(4), B.x + B.w - api.u(70));      const ty = clamp(api.pointer.y - api.u(8), B.y + f1 * 1.4, B.y + B.h - api.u(6));      api.label(`${fmt(d.v, 2)} c  ${kind}`, tx, ty, { size: f1, font: "mono" });    }  } else if (api.pointer.inside && !api.pointer.down && api.pointer.x >= B.x && api.pointer.x <= B.x + B.w && api.pointer.y >= B.y && api.pointer.y <= B.y + B.h) {    api.circle(api.pointer.x, api.pointer.y, api.u(9), { w: api.u(1), dash: [api.u(2), api.u(2.5)] });  }  ctx.restore();  // 視点の印(下が見ている側)・凡例・題  const ex = cx, ey = B.y + B.h - api.u(5);  api.poly([[ex, ey - api.u(8)], [ex - api.u(5), ey], [ex + api.u(5), ey]], { close: true, fill: "ink", stroke: false });  api.text("視点", ex + api.u(8), ey - api.u(1), { size: f1 * 0.9 });  const lx = B.x + api.u(7);  let ly = B.y + f1 * 1.3;  api.label(`真上から  半径 ${fmt(s.mapR, 0)} M`, lx, ly, { size: f1, pad: 1 });  const leg: [string, number[]][] = [["2M 地平線", []], ["3M 光子球", [api.u(1.2), api.u(2.2)]], ["6M ISCO", [api.u(4), api.u(3)]]];  for (const [t, dash] of leg) { ly += f1 * 1.25; api.line(lx, ly - f1 * 0.32, lx + api.u(14), ly - f1 * 0.32, { w: api.u(1), dash }); api.label(t, lx + api.u(18), ly, { size: f1 * 0.92, pad: 1 }); }  api.label("この図をドラッグすると天体を投げられます", B.x + B.w, B.y - api.u(7), { size: f1, align: "right" });} function drawGraph(api: Api, s: S, G: Box) {  const ctx = api.ctx, f1 = fs(api, 9.5), b = s.sel;  api.rect(G.x, G.y, G.w, G.h, { fill: "paper", w: api.u(1) });  api.text(b ? `有効ポテンシャル V(r)  天体 ${b.id}  L = ${fmt(Math.abs(b.L), 2)} M` : "有効ポテンシャル V(r)", G.x + api.u(7), G.y + f1 * 1.35, { size: f1 });  if (!b) { api.text("天体を選ぶと出ます", G.x + G.w / 2, G.y + G.h / 2 + f1 * 0.4, { size: f1, align: "center" }); return; }  const L2 = b.L * b.L, En = (b.E * b.E - 1) / 2, inf = b.info;  const V = (r: number) => -1 / r + L2 / (2 * r * r) - L2 / (r * r * r), VN = (r: number) => -1 / r + L2 / (2 * r * r);  const r1 = clamp(Math.max(24, (inf && Number.isFinite(inf.ra) ? inf.ra : b.r) * 1.3), 24, 90);  let lo = En;  for (let k = 0; k <= 60; k++) { const r = 4 + ((r1 - 4) * k) / 60; lo = Math.min(lo, V(r), VN(r)); }  const span = Math.max(0.02, Math.max(En, 0) - lo), hi = Math.max(En, 0) + 0.55 * span;  lo -= 0.3 * span;  const x0 = G.x + api.u(10), x1 = G.x + G.w - api.u(10), y0 = G.y + f1 * 2, y1 = G.y + G.h - f1 * 1.9;  const X = (r: number) => x0 + ((r - 2) / (r1 - 2)) * (x1 - x0), Y = (v: number) => y0 + ((hi - v) / (hi - lo)) * (y1 - y0);  ctx.save(); ctx.beginPath(); ctx.rect(x0, y0, x1 - x0, y1 - y0); ctx.clip();  const curve = (f: (r: number) => number, o: { w: number; dash?: number[] }) => {    const pts: [number, number][] = [];    for (let k = 0; k <= 160; k++) { const r = 2 + ((r1 - 2) * k) / 160, v = f(r); pts.push([X(r), clamp(Y(v), y0 - 400, y1 + 400)]); }    api.poly(pts, o);  };  api.line(x0, Y(0), x1, Y(0), { w: api.u(0.5), dash: [api.u(1), api.u(3)] });  api.line(X(6), y0, X(6), y1, { w: api.u(0.5), dash: [api.u(1), api.u(3)] });  curve(VN, { w: api.u(1), dash: [api.u(4), api.u(3)] });  curve(V, { w: api.u(1.8) });  api.line(x0, Y(En), x1, Y(En), { w: api.u(0.8) });  if (inf && inf.kind === "bound") api.line(X(inf.rp), Y(En), X(Math.min(inf.ra, r1)), Y(En), { w: api.u(2.6) });  api.circle(X(Math.min(b.r, r1)), Y(En), api.u(3.6), { fill: "ink", stroke: false });  if (L2 > 12) {    const um = (1 + Math.sqrt(1 - 12 / L2)) / 6, ul = (1 - Math.sqrt(1 - 12 / L2)) / 6;    const yp = Y(V(1 / um)), yv = Y(V(1 / ul));    if (yp > y0 + f1) api.text("山", X(1 / um), yp - api.u(5), { size: f1, align: "center" });    if (1 / ul < r1 && yv < y1 - f1 * 0.2) api.text("谷", X(1 / ul), yv + f1 * 1.1, { size: f1, align: "center" });  }  ctx.restore();  for (const r of [2, 6, 12, 24, 48]) if (r <= r1) {    api.line(X(r), y1, X(r), y1 + api.u(4), { w: api.u(0.8) });    api.text(`${r}M`, X(r), y1 + f1 * 1.15, { size: f1 * 0.9, font: "mono", align: "center" });  }  api.line(x0, y1, x1, y1, { w: api.u(0.8) });  const lgx = x1 - api.u(126), lgy = y0 + f1 * 0.9;  api.line(lgx, lgy - f1 * 0.3, lgx + api.u(14), lgy - f1 * 0.3, { w: api.u(1.8) }); api.text("相対論", lgx + api.u(18), lgy, { size: f1 * 0.9 });  api.line(lgx + api.u(62), lgy - f1 * 0.3, lgx + api.u(76), lgy - f1 * 0.3, { w: api.u(1), dash: [api.u(4), api.u(3)] }); api.text("ニュートン", lgx + api.u(80), lgy, { size: f1 * 0.9 });  api.label("天体のエネルギー", x0 + api.u(18), clamp(Y(En) - api.u(5), y0 + f1, y1 - api.u(2)), { size: f1 * 0.9 });} function readout(api: Api, s: S) {  const M = MASS[api.str("mass")] ?? MASS.sgr, b = s.sel, lines: string[] = [];  lines.push(`遠くの時計  t = ${num(s.time, 0)} M(${M.label}なら ${realTime(s.time * M.m * SEC)})`);  if (b) {    const lapse = 1 - 2 / b.r, vphi = (b.L * Math.sqrt(lapse)) / (b.E * b.r), vr = b.pr / b.E, v = Math.min(0.9999, Math.hypot(vr, vphi));    lines.push(`天体 ${b.id}  r = ${fmt(b.r, 1)} M(${realLen(b.r * M.m * KM)})  速さ ${fmt(v, 2)} c  時計の進み ${fmt(lapse / b.E, 2)} 倍`);    const o = b.info;    if (!o) lines.push("軌道  -");    else if (o.kind === "circle") lines.push(`軌道  円(1 周 ${num(TAU * Math.pow(b.r, 1.5), 0)} M = ${realTime(TAU * Math.pow(b.r, 1.5) * M.m * SEC)})`);    else if (o.kind === "bound") lines.push(`軌道  近点 ${fmt(o.rp, 1)} M・遠点 ${fmt(o.ra, 1)} M  近点の移動 1 周で ${fmt(o.prec, 0)}°${Number.isFinite(b.meas) ? `(計算では ${fmt(b.meas, 0)}°)` : ""}`);    else if (o.kind === "escape") lines.push("軌道  重力を振り切って遠くへ飛び去る");    else lines.push(Math.abs(b.L) < Math.sqrt(12) ? `軌道  地平線へ落ちる(L が √12 M より小さく、ポテンシャルに谷が無い)` : "軌道  地平線へ落ちる(ポテンシャルの山を越えるか、山の内側にいる)");  } else lines.push("天体  なし(右下の図をドラッグすると投げられます)");  const D = api.num("dist");  lines.push(`見ている所  傾き ${fmt(api.num("incl"), 0)}°  距離 ${fmt(D, 0)} M(${realLen(D * M.m * KM)})`);  if (api.bool("grid")) {    const as = Math.asin(Math.min(1, (Math.sqrt(27) * Math.sqrt(1 - 2 / D)) / D)), ah = Math.asin(2 / D);    lines.push(`影の見かけの半径  ${fmt((as * 180) / Math.PI, 1)}°(光がまっすぐ進むときの地平線の ${fmt(as / ah, 2)} 倍)`);  }  api.readout(lines);}

/** 地図の上の点(画面の px)→ 赤道面の座標(M) */function fromMap(s: S, X: number, Y: number): [number, number] {  const t = s.mt!; const a = (X - t.cx) / t.k, b = (Y - t.cy) / t.k;  return [-a * t.sp + b * t.cp, a * t.cp + b * t.sp];}const inBox = (B: Box | null, x: number, y: number) => !!B && x >= B.x && x <= B.x + B.w && y >= B.y && y <= B.y + B.h; function throwVel(api: Api, s: S, d: Drag) {  const B = s.map!, [qx, qy] = fromMap(s, api.pointer.x, api.pointer.y);  let dx = qx - d.px, dy = qy - d.py;  const lenPx = Math.hypot(api.pointer.x - api.pointer.x0, api.pointer.y - api.pointer.y0);  const v = Math.min(0.97, lenPx / (0.42 * B.w)), l = Math.hypot(dx, dy) || 1;  dx /= l; dy /= l;  const r = Math.hypot(d.px, d.py), ux = d.px / r, uy = d.py / r;  let vr = v * (dx * ux + dy * uy), vt = v * (-dx * uy + dy * ux);  d.snap = false;  if (r > 3.3) {    const vc = 1 / Math.sqrt(r - 2);    if (Math.abs(vr) < 0.08 * vc && Math.abs(Math.abs(vt) - vc) < 0.06 * vc) { vr = 0; vt = Math.sign(vt) * vc; d.snap = true; }  }  d.vr = vr; d.vt = vt; d.v = Math.hypot(vr, vt);} /** 指の操作: 右下の図の上なら天体を投げる(たたくと選ぶ)、ほかの所なら見る向きを変える */function onPointer(type: "down" | "move" | "up" | "cancel" | "wheel", api: Api, s: S) {  const P = api.pointer;  if (type === "down") {    if (inBox(s.map, P.x, P.y) && s.mt) {      const [px, py] = fromMap(s, P.x, P.y);      s.drag = { mode: "throw", lx: P.x, ly: P.y, moved: false, px, py, vr: 0, vt: 0, v: 0, snap: false, path: new Float32Array(2400), pn: 0, info: null };    } else s.drag = { mode: "cam", lx: P.x, ly: P.y, moved: false, px: 0, py: 0, vr: 0, vt: 0, v: 0, snap: false, path: new Float32Array(2), pn: 0, info: null };    return;  }  const d = s.drag;  if (!d) return;  if (type === "move") {    if (Math.hypot(P.x - P.x0, P.y - P.y0) > 5) d.moved = true;    if (d.mode === "cam") {      const dx = P.x - d.lx, dy = P.y - d.ly;      d.lx = P.x; d.ly = P.y;      s.camPhi -= dx * 0.008;      api.set("incl", Math.round(clamp(api.num("incl") - dy * 0.35, 0, 90)));    } else if (d.moved) {      throwVel(api, s, d);      const r = Math.hypot(d.px, d.py), phi = Math.atan2(d.py, d.px), tmp = local(s, Math.max(r, 2.3), phi, d.vr, d.vt);      s.nextId--;      d.info = classify(tmp.E, tmp.L, tmp.r, tmp.pr);      d.pn = predict(tmp.E, tmp.L, tmp.r, tmp.pr, phi, d.path, 3 * s.mapR, TAU * 3);    }    return;  }  s.drag = null;  if (type === "cancel" || d.mode === "cam") return;  const r = Math.hypot(d.px, d.py), phi = Math.atan2(d.py, d.px);  if (!d.moved) {    // たたいた: 近くの天体を選ぶ    let best: Body | null = null, bd = api.u(16);    const t = s.mt;    if (t) for (const b of s.bodies) {      const x = b.r * Math.cos(b.phi), y = b.r * Math.sin(b.phi);      const X = t.cx + t.k * (-x * t.sp + y * t.cp), Y = t.cy + t.k * (x * t.cp + y * t.sp), dd = Math.hypot(X - P.x, Y - P.y);      if (dd < bd) { bd = dd; best = b; }    }    if (best) { s.sel = best; best.guideAt = -1e9; }    return;  }  if (r < 2.3) { api.note("地平線に近すぎて投げられません。2.3 M より外から投げると動きます"); return; }  throwVel(api, s, d);  const b = local(s, r, phi, d.vr, d.vt);  addBody(s, b); s.sel = b; b.guideAt = -1e9;  api.emit("throw", { r: Math.round(r * 10) / 10, v: Math.round(d.v * 100) / 100 });}

const BC = Math.sqrt(27); // 影の縁になる衝突パラメータinterface Ray { pts: Float32Array; n: number; fate: "out" | "in" | "loop"; turn: number; rmin: number }/** 衝突パラメータ b の光をたどる。点は (x, y)(ブラックホールが原点、光は左から右へ) */function trace(b: number, out: Float32Array, rView: number): Ray {  const sb = Math.sign(b) || 1, ab = Math.max(Math.abs(b), 1e-6);  // 遠く(r = 1000 M)から始める。そこでの傾き du/dψ は、エネルギーの式 (du/dψ)² = 1/b² − u²(1 − 2u) から厳密に決める  let u = Math.min(0.001, 0.999 / ab), psi = Math.asin(clamp(u * ab, 0, 1)), w = Math.sqrt(Math.max(0, 1 / (ab * ab) - u * u * (1 - 2 * u)));  const f = (x: number) => 3 * x * x - x;  let n = 0, rmin = 1 / u;  const push = () => { if (2 * n + 1 < out.length) { const r = 1 / u; out[2 * n] = -r * Math.cos(psi); out[2 * n + 1] = sb * r * Math.sin(psi); n++; } };  push();  let fate: Ray["fate"] = "loop";  for (let k = 0; k < 30000; k++) {    // 中心に近いほど刻みを細かく(光子球のそばを回る光は、少しの差で行き先が変わる)    const h = u > 0.25 ? 0.002 : u > 0.08 ? 0.008 : 0.03;    const k1u = w, k1w = f(u), k2u = w + 0.5 * h * k1w, k2w = f(u + 0.5 * h * k1u);    const k3u = w + 0.5 * h * k2w, k3w = f(u + 0.5 * h * k2u), k4u = w + h * k3w, k4w = f(u + h * k3u);    u += (h / 6) * (k1u + 2 * k2u + 2 * k3u + k4u); w += (h / 6) * (k1w + 2 * k2w + 2 * k3w + k4w); psi += h;    if (u >= 0.5) { u = 0.5; push(); fate = "in"; break; }    if (u <= 1 / (3 * rView) && w < 0) { push(); fate = "out"; psi += Math.asin(clamp(u * ab, 0, 1)); break; } // 残りのまっすぐな部分の角を足す    rmin = Math.min(rmin, 1 / u);    if (k % 3 === 0 && 1 / u < 3 * rView) push();    if (psi > 7 * Math.PI) break;  }  return { pts: out, n, fate, turn: psi, rmin };}