電波は、アンテナからちぎれた電場の輪になって飛んでいく

export const C0 = 299792458; // 光の速さ m/sexport const EPS0 = 8.8541878188e-12; // 真空の誘電率 F/mexport const MU0 = 1.25663706127e-6; // 真空の透磁率 N/A²const F_DIP = 2.4e9; // 式の場面の周波数 Hz(Wi-Fi の 2.4 GHz 帯)const P_DIP = 0.1; // 式の場面のアンテナが出す電力 Wexport const I_ANT = 0.01; // 平面の場面のアンテナに流す電流の振幅 Aconst CPL = [6, 8, 10]; // 1 波長に入れる格子のます(画質 0・1・2)export const S_C = 0.5; // クーラン数 cΔt/Δx。2 次元では 1/√2 以下なら計算が発散しないexport const NPML = 10; // 端の吸収層(PML)の厚さ(ます)const MAXSTEP = [4, 6, 8]; // 1 コマに進める刻みの上限(画質ごと)const MAXANT = 12; // アンテナの数の上限/** 壁の材質。比誘電率 ε = a、導電率 σ = c·f^d(f は GHz・σ は S/m)。t = 描く厚さ(m)・tone = 描く濃さ */export interface MatDef { label: string; a: number; c: number; d: number; t: number; tone: number; pec?: boolean }export const MATS: Record<string, MatDef> = {  metal: { label: "金属", a: 1, c: 1e7, d: 0, t: 0.002, tone: 0.3, pec: true },  concrete: { label: "コンクリート", a: 5.31, c: 0.0326, d: 0.8095, t: 0.15, tone: 0.6 },  board: { label: "石膏ボード", a: 2.94, c: 0.0116, d: 0.7076, t: 0.0125, tone: 0.95 },  wood: { label: "木", a: 1.99, c: 0.0047, d: 1.0718, t: 0.03, tone: 0.8 },};const MKEY = ["", "metal", "concrete", "board", "wood"]; // 格子に入れる番号(0 = 空気)

/** 短いダイポール(双極子モーメント p = p0 cos ωt が縦に振動する)。K = p0/4πε0、KB = μ0 ω k p0/4π */export interface Dip { k: number; w: number; lam: number; K: number; KB: number }export function dipole(f: number, P: number): Dip {  // 出す電力の平均 P = μ0 p0² ω⁴ / (12π c) から p0 を決める  const w = TAU * f, k = w / C0, p0 = Math.sqrt((12 * Math.PI * C0 * P) / (MU0 * w ** 4));  return { k, w, lam: C0 / f, K: p0 / (4 * Math.PI * EPS0), KB: (MU0 * w * k * p0) / (4 * Math.PI) };}/** 距離 r・軸からの角 θ の場所の複素振幅。out = [Er, Eθ, Bφ] の実部と虚部。時刻 t の値は 実部·cos ωt + 虚部·sin ωt *  Er = K·2cosθ·(1/r³ − ik/r²)·e^{ikr} *  Eθ = K·sinθ·(1/r³ − ik/r² − k²/r)·e^{ikr} *  Bφ = −KB·sinθ·(1/r)·(1 + i/kr)·e^{ikr} */export function dipField(d: Dip, r: number, ct: number, st: number, out: Float64Array) {  const k = d.k, kr = k * r, c = Math.cos(kr), s = Math.sin(kr);  const a = 1 / (r * r * r), b = k / (r * r), f = (k * k) / r;  const nr = a * c + b * s, ni = a * s - b * c; // (1/r³ − ik/r²)·e^{ikr}  out[0] = d.K * 2 * ct * nr; out[1] = d.K * 2 * ct * ni;  out[2] = d.K * st * (nr - f * c); out[3] = d.K * st * (ni - f * s);  const g = (d.KB * st) / r;  out[4] = -g * (c - s / kr); out[5] = -g * (s + c / kr);}/** 電気力線の関数。この値が同じ点を結ぶと電気力線になる(遠くでは −1〜1 のあいだを動く) *  ψ = sin²θ·[sin(kr − ωt) + cos(kr − ωt)/kr] */export const psi = (kr: number, s2: number, wt: number) => s2 * (Math.sin(kr - wt) + Math.cos(kr - wt) / kr);

const DL = 0.34, PMAX = 10; // 線を引く ψ の間隔と、片側の本数の上限/** 縦の面(半径 R の円)の 4 分の 1(s ≥ 0・y ≥ 0)で ψ を格子の点ごとに求め、同じ値を結ぶ線の切れ端を作る(マーチング・スクエア法)。 *  場は上下と左右に対称なので、残りの 4 分の 3 は描くときに写す */function contours(D: DipState, wt: number) {  const { ni, nj, ps, seg } = D, hs = D.R / ni, hy = D.Y / nj, k = D.d.k, rmin = 0.06 * D.d.lam, W = ni + 1;  for (let j = 0; j <= nj; j++) for (let i = 0; i <= ni; i++) {    const s = i * hs, y = j * hy, r = Math.hypot(s, y);    ps[j * W + i] = r < rmin || r > D.R ? 1e9 : psi(k * r, (s * s) / (r * r), wt);  }  let n = 0;  const max = seg.length / 4;  for (let j = 0; j < nj; j++) for (let i = 0; i < ni; i++) {    const a = ps[j * W + i], b = ps[j * W + i + 1], c = ps[(j + 1) * W + i + 1], d = ps[(j + 1) * W + i];    if (a > 1e8 || b > 1e8 || c > 1e8 || d > 1e8) continue;    const lo = Math.max(-PMAX, Math.ceil(Math.min(a, b, c, d) / DL - 0.5)), hi = Math.min(PMAX - 1, Math.floor(Math.max(a, b, c, d) / DL - 0.5));    for (let q = lo; q <= hi && n < max - 2; q++) {      const L = (q + 0.5) * DL;      const id = (a > L ? 1 : 0) | (b > L ? 2 : 0) | (c > L ? 4 : 0) | (d > L ? 8 : 0);      if (id === 0 || id === 15) continue;      // 4 つの辺(0 = 下・1 = 右・2 = 上・3 = 左)の上で ψ = L になる点      const ex = (e: number) => (e === 0 ? i + (L - a) / (b - a) : e === 1 ? i + 1 : e === 2 ? i + (L - d) / (c - d) : i) * hs;      const ey = (e: number) => (e === 0 ? j : e === 1 ? j + (L - b) / (c - b) : e === 2 ? j + 1 : j + (L - a) / (d - a)) * hy;      const put = (e1: number, e2: number) => { const o = 4 * n++; seg[o] = ex(e1); seg[o + 1] = ey(e1); seg[o + 2] = ex(e2); seg[o + 3] = ey(e2); };      const mid = (a + b + c + d) / 4 > L;      switch (id) {        case 1: case 14: put(3, 0); break;        case 2: case 13: put(0, 1); break;        case 3: case 12: put(3, 1); break;        case 4: case 11: put(1, 2); break;        case 6: case 9: put(0, 2); break;        case 7: case 8: put(3, 2); break;        case 5: if (mid) { put(0, 1); put(2, 3); } else { put(3, 0); put(1, 2); } break;        case 10: if (mid) { put(3, 0); put(1, 2); } else { put(0, 1); put(2, 3); } break;      }    }  }  D.nseg = n;}/** 縦の面の点 (s, y) を画面の位置にする式(射影変換)。面の上の点は s·U + y·(上向き) なので、 *  カメラの行列を掛けた結果は s と y の 1 次式になる。h = [x の 3 係数, y の 3 係数, w の 3 係数] */function planeMap(v: View3, U: THREE.Vector3, h: Float64Array) {  const m = TMP_M.multiplyMatrices(v.camera.projectionMatrix, v.camera.matrixWorldInverse).elements;  for (let r = 0; r < 3; r++) {    const row = r === 2 ? 3 : r; // 同次座標の x・y・w の行    h[r * 3] = m[12 + row]; h[r * 3 + 1] = m[row] * U.x + m[8 + row] * U.z; h[r * 3 + 2] = m[4 + row];  }}const TMP_M = new THREE.Matrix4();/** 縦の面の点 → 画面(CSS px)。カメラの後ろなら false */function mapPt(api: Api, h: Float64Array, s: number, y: number, out: number[]): boolean {  const W = h[6] + h[7] * s + h[8] * y;  if (W <= 1e-6) return false;  out[0] = ((h[0] + h[1] * s + h[2] * y) / W + 1) * 0.5 * api.w;  out[1] = (1 - (h[3] + h[4] * s + h[5] * y) / W) * 0.5 * api.h;  return true;}

/** 水平な面の格子。縦向きの電場 E(ますの点)と、水平な磁場 Hx・Hz(点と点のあいだ)を、半刻みずらして交互に進める(Yee の格子)。 *  磁場には Z0 = √(μ0/ε0) を掛けて電場と同じ単位(V/m)にしてある。こうすると真空での係数が、どちらもクーラン数 S になる。 *  端の吸収層では、E を x の向きの分と z の向きの分に分けて、それぞれに吸収(導電率)を入れる(ベランジェの PML) */export class Grid {  readonly nx: number; readonly nz: number; readonly dx: number; readonly dt: number; readonly w: number; readonly k: number; readonly spp: number;  readonly x0: number; readonly z0: number;  readonly ey: Float32Array; readonly eyx: Float32Array; readonly eyz: Float32Array; readonly hx: Float32Array; readonly hz: Float32Array; readonly avg: Float32Array;  readonly cax: Float32Array; readonly cbx: Float32Array; readonly caz: Float32Array; readonly cbz: Float32Array; readonly mid: Uint8Array;  readonly dax: Float32Array; readonly dbx: Float32Array; readonly daz: Float32Array; readonly dbz: Float32Array;  private readonly aex: Float32Array; private readonly aez: Float32Array; private readonly al: number;  n = 0; src: { c: number; ph: number }[] = [];  constructor(readonly f: number, Lx: number, Lz: number, readonly cpl: number) {    const lam = C0 / f;    this.dx = lam / cpl; this.dt = (S_C * this.dx) / C0; this.w = TAU * f; this.k = TAU / lam; this.spp = cpl / S_C;    const ix = Math.max(8, Math.round(Lx / this.dx)), iz = Math.max(8, Math.round(Lz / this.dx));    const nx = (this.nx = ix + 1 + 2 * NPML), nz = (this.nz = iz + 1 + 2 * NPML), N = nx * nz;    this.x0 = -(ix / 2 + NPML) * this.dx; this.z0 = -(iz / 2 + NPML) * this.dx;    const F = () => new Float32Array(N);    this.ey = F(); this.eyx = F(); this.eyz = F(); this.hx = F(); this.hz = F(); this.avg = F();    this.cax = F(); this.cbx = F(); this.caz = F(); this.cbz = F(); this.mid = new Uint8Array(N);    this.dax = new Float32Array(nx); this.dbx = new Float32Array(nx); this.daz = new Float32Array(nz); this.dbz = new Float32Array(nz);    this.aex = new Float32Array(nx); this.aez = new Float32Array(nz);    // 吸収層の導電率: 内側の縁の 0 から外へ深さの 3 乗で増やす。正面から入った波が層の奥ではね返って戻ってきたとき、    // 振幅が 10⁻⁶ になる強さにする。    // 値は σΔt/2ε0 に直してある(磁場の側は σ* = σ·μ0/ε0 で、同じ値になる)    const amax = (S_C * 4 * Math.log(1e6)) / (4 * NPML);    const prof = (d: number) => amax * Math.min(1, Math.max(0, d)) ** 3;    const fill = (n: number, ae: Float32Array, da: Float32Array, db: Float32Array) => {      const hi = n - 1 - NPML;      for (let i = 0; i < n; i++) {        ae[i] = prof(Math.max(NPML - i, i - hi) / NPML);        const ah = prof(Math.max(NPML - i - 0.5, i + 0.5 - hi) / NPML);        da[i] = (1 - ah) / (1 + ah); db[i] = S_C / (1 + ah);      }    };    fill(nx, this.aex, this.dax, this.dbx); fill(nz, this.aez, this.daz, this.dbz);    for (let j = 0; j < nz; j++) for (let i = 0; i < nx; i++) this.vacuum(j * nx + i, i, j);    this.al = 1 / (3 * this.spp); // 強さの平均は 3 周期ほどの時定数でならす  }  private vacuum(c: number, i: number, j: number) {    const ax = this.aex[i], az = this.aez[j];    this.cax[c] = (1 - ax) / (1 + ax); this.cbx[c] = S_C / (1 + ax);    this.caz[c] = (1 - az) / (1 + az); this.cbz[c] = S_C / (1 + az);    this.mid[c] = 0;  }  /** 1 刻み(Δt)進める */  step() {    const { nx, nz, ey, eyx, eyz, hx, hz, cax, cbx, caz, cbz, dax, dbx, daz, dbz, avg, al } = this;    // 磁場: ∂Hx/∂t = (1/μ0)·∂E/∂z、∂Hz/∂t = −(1/μ0)·∂E/∂x    for (let j = 0; j < nz - 1; j++) { const a = daz[j], b = dbz[j]; for (let c = j * nx, e = c + nx; c < e; c++) hx[c] = a * hx[c] + b * (ey[c + nx] - ey[c]); }    for (let j = 0; j < nz; j++) for (let i = 0, c = j * nx; i < nx - 1; i++, c++) hz[c] = dax[i] * hz[c] - dbx[i] * (ey[c + 1] - ey[c]);    // 電場: ε·∂E/∂t + σE = ∂Hx/∂z − ∂Hz/∂x(いちばん外の点は 0 のまま)    for (let j = 1; j < nz - 1; j++) for (let c = j * nx + 1, e = j * nx + nx - 1; c < e; c++) {      const a = cax[c] * eyx[c] - cbx[c] * (hz[c] - hz[c - 1]);      const b = caz[c] * eyz[c] + cbz[c] * (hx[c] - hx[c - nx]);      eyx[c] = a; eyz[c] = b; const v = a + b; ey[c] = v; avg[c] += al * (v * v - avg[c]);    }    // アンテナの電流 I(t) = I0·sin(ωt + φ)。電流の密度 J = I/Δ² を電場の式に入れる(ΔE = −Δt·J/ε)。最初の 3 周期でゆっくり強める    const t = (this.n + 0.5) * this.dt, T3 = (3 * TAU) / this.w, ramp = t < T3 ? 0.5 - 0.5 * Math.cos((Math.PI * t) / T3) : 1;    const kJ = ((this.dt / (EPS0 * this.dx * this.dx)) * I_ANT * ramp) / S_C;    for (const s of this.src) { const c = s.c; eyx[c] -= kJ * cbx[c] * Math.sin(this.w * t + s.ph); ey[c] = eyx[c] + eyz[c]; }    this.n++;  }  /** 場所 (x, z) にいちばん近い格子の点(吸収層の外に収める) */  cell(x: number, z: number) {    const i = clamp(Math.round((x - this.x0) / this.dx), NPML, this.nx - 1 - NPML), j = clamp(Math.round((z - this.z0) / this.dx), NPML, this.nz - 1 - NPML);    return j * this.nx + i;  }  /** アンテナを置く。n 本目の電流の位相は −n·Δφ(右へ並べると、右ほど遅れる) */  setSources(ants: { x: number; z: number }[], dph: number) { this.src = ants.map((a, n) => ({ c: this.cell(a.x, a.z), ph: -n * dph })); }  /** 配列の値を (x, z) で読む(まわりの 4 点から直線で補う) */  sample(arr: Float32Array, x: number, z: number) {    const fi = clamp((x - this.x0) / this.dx, 0, this.nx - 1.001), fj = clamp((z - this.z0) / this.dx, 0, this.nz - 1.001);    const i = Math.floor(fi), j = Math.floor(fj), u = fi - i, w = fj - j, c = j * this.nx + i;    return (arr[c] * (1 - u) + arr[c + 1] * u) * (1 - w) + (arr[c + this.nx] * (1 - u) + arr[c + this.nx + 1] * u) * w;  }  /** 格子に入れる ε と σ を、格子の上の波が本物の材質の中と同じ波長と減り方で進むように直す。   *  コンクリートの中では波長が 2.3 分の 1 に縮み、1 波長が 3〜4 ます(画質が低いと 3 ます未満)しかなくなる。そのままでは格子の上の波が遅れすぎたり、   *  まったく進めなくなったりする。格子の上の波は「k の代わりに (2/Δ)·sin(kΔ/2)、ω の代わりに (2/Δt)·sin(ωΔt/2)」で   *  ふつうの式に従うので、本物の複素数の波数 k = (ω/c)·√(ε + iσ/ωε0) から逆に ε と σ を決める */  corrected(er: number, sg: number): [number, number] {    const ei = sg / (this.w * EPS0), mod = Math.hypot(er, ei), nr = Math.sqrt((mod + er) / 2), ni = Math.sqrt((mod - er) / 2);    const a = (Math.PI / this.cpl) * nr, b = (Math.PI / this.cpl) * ni; // kΔ/2(複素数)    const sr = Math.sin(a) * Math.cosh(b), si = Math.cos(a) * Math.sinh(b), q = Math.sin((Math.PI * S_C) / this.cpl);    const xr = (S_C * S_C * (sr * sr - si * si)) / (q * q), xi = (S_C * S_C * 2 * sr * si) / (q * q);    return [xr, (xi * EPS0 * (2 / this.dt) * q) / Math.cos((Math.PI * S_C) / this.cpl)];  }  /** 壁を塗る。線分 (x1, z1)-(x2, z2) から厚さの半分より近い点を、材質 key(空なら空気)にする。box なら長方形を塗りつぶす */  paint(x1: number, z1: number, x2: number, z2: number, key: string, th?: number, box = false) {    const m = key ? MATS[key] : null, mi = key ? MKEY.indexOf(key) : 0, dx = this.dx;    const t = th ?? m?.t ?? 0.12;    // 格子より薄い板は線から半ますの内側の点(ふつうは 1 列)に入れ、厚さの比 t/Δ で (ε − 1) と σ を薄める(薄い板の反射と吸収が同じになる)。    // 斜めに描いても、点が角でつながるので電波は抜けない    const thin = !!m && t < 1.5 * dx, half = box ? 0 : thin ? 0.5 * dx : t / 2;    let ca = 1, cb = S_C;    if (m?.pec) { ca = 0; cb = 0; }    else if (m) {      let er = m.a, sg = m.c * (this.f / 1e9) ** m.d;      if (thin) { er = 1 + (er - 1) * (t / dx); sg *= t / dx; }      else [er, sg] = this.corrected(er, sg);      const a = (sg * this.dt) / (2 * EPS0 * er); ca = (1 - a) / (1 + a); cb = S_C / (er * (1 + a));    }    const lim = !m || m.pec ? 1 : NPML; // 金属と空気だけは吸収層の中まで塗れる(どこまでも続く板)    const i0 = Math.max(lim, Math.floor((Math.min(x1, x2) - half - this.x0) / dx)), i1 = Math.min(this.nx - 1 - lim, Math.ceil((Math.max(x1, x2) + half - this.x0) / dx));    const j0 = Math.max(lim, Math.floor((Math.min(z1, z2) - half - this.z0) / dx)), j1 = Math.min(this.nz - 1 - lim, Math.ceil((Math.max(z1, z2) + half - this.z0) / dx));    const ux = x2 - x1, uz = z2 - z1, L2 = ux * ux + uz * uz, h2 = half * half * 1.0001 + 1e-12;    let changed = false;    for (let j = j0; j <= j1; j++) for (let i = i0; i <= i1; i++) {      const px = this.x0 + i * dx, pz = this.z0 + j * dx;      if (box) { if (px < Math.min(x1, x2) || px > Math.max(x1, x2) || pz < Math.min(z1, z2) || pz > Math.max(z1, z2)) continue; }      else {        const s = L2 > 0 ? clamp(((px - x1) * ux + (pz - z1) * uz) / L2, 0, 1) : 0;        const ex = px - x1 - s * ux, ez = pz - z1 - s * uz;        if (ex * ex + ez * ez > h2) continue;      }      const c = j * this.nx + i;      if (mi === 0) { if (this.mid[c] !== 0) { this.vacuum(c, i, j); changed = true; } continue; }      if (this.mid[c] !== mi) changed = true;      this.mid[c] = mi; this.cax[c] = this.caz[c] = ca; this.cbx[c] = this.cbz[c] = cb;      if (m?.pec) { this.ey[c] = this.eyx[c] = this.eyz[c] = 0; }    }    return changed;  }}

/** ベッセル関数 J0 と Y0(Numerical Recipes の有理式による近似・誤差は 10⁻⁷ ほど) */export function besJ0(x: number) {  const ax = Math.abs(x);  if (ax < 8) {    const y = x * x;    return (57568490574.0 + y * (-13362590354.0 + y * (651619640.7 + y * (-11214424.18 + y * (77392.33017 + y * -184.9052456))))) /      (57568490411.0 + y * (1029532985.0 + y * (9494680.718 + y * (59272.64853 + y * (267.8532712 + y)))));  }  const z = 8 / ax, y = z * z, xx = ax - 0.785398164;  const p1 = 1 + y * (-0.1098628627e-2 + y * (0.2734510407e-4 + y * (-0.2073370639e-5 + y * 0.2093887211e-6)));  const p2 = -0.1562499995e-1 + y * (0.1430488765e-3 + y * (-0.6911147651e-5 + y * (0.7621095161e-6 - y * 0.934935152e-7)));  return Math.sqrt(0.636619772 / ax) * (Math.cos(xx) * p1 - z * Math.sin(xx) * p2);}export function besY0(x: number) {  if (x < 8) {    const y = x * x;    return (-2957821389.0 + y * (7062834065.0 + y * (-512359803.6 + y * (10879881.29 + y * (-86327.92757 + y * 228.4622733))))) /      (40076544269.0 + y * (745249964.8 + y * (7189466.438 + y * (47447.2647 + y * (226.1030244 + y))))) + 0.636619772 * besJ0(x) * Math.log(x);  }  const z = 8 / x, y = z * z, xx = x - 0.785398164;  const p1 = 1 + y * (-0.1098628627e-2 + y * (0.2734510407e-4 + y * (-0.2073370639e-5 + y * 0.2093887211e-6)));  const p2 = -0.1562499995e-1 + y * (0.1430488765e-3 + y * (-0.6911147651e-5 + y * (0.7621095161e-6 - y * 0.934935152e-7)));  return Math.sqrt(0.636619772 / x) * (Math.sin(xx) * p1 + z * Math.cos(xx) * p2);}/** 壁が無いときの、場所 (x, z) の電場の実効値(V/m)。縦に長い電流 I が作る電場 E = −(ωμ0 I/4)·H0⁽²⁾(kρ) を、 *  アンテナごとの位相をつけて足し合わせる(H0⁽²⁾ = J0 − iY0)。one なら 1 本目だけ */export function freeRms(g: Grid, ants: { x: number; z: number }[], dph: number, x: number, z: number, one: boolean) {  let re = 0, im = 0;  const n = one ? Math.min(1, ants.length) : ants.length;  for (let a = 0; a < n; a++) {    const kr = Math.max(0.3, g.k * Math.hypot(x - ants[a].x, z - ants[a].z)), J = besJ0(kr), Y = besY0(kr), ph = -a * dph;    re += Math.cos(ph) * J + Math.sin(ph) * Y; im += Math.sin(ph) * J - Math.cos(ph) * Y;  }  return ((g.w * MU0 * I_ANT) / 4) * Math.hypot(re, im) / Math.SQRT2;}/** 並べたアンテナの電波の向き: 中心 (cx, cz) から半径 R の弧の上で、強さの平均がいちばん大きい角度(前 = −z から +x の側へ測る) */export function beamAngle(g: Grid, cx: number, cz: number, R: number) {  let best = -1, bi = 0;  const N = 161, v = new Float64Array(N);  for (let n = 0; n < N; n++) {    const th = ((-80 + (160 * n) / (N - 1)) * Math.PI) / 180;    v[n] = g.sample(g.avg, cx + R * Math.sin(th), cz - R * Math.cos(th));    if (v[n] > best) { best = v[n]; bi = n; }  }  // となりの 2 点と放物線でつないで、頂点の角度を細かく求める  let d = 0;  if (bi > 0 && bi < N - 1) { const a = v[bi - 1], b = v[bi], c = v[bi + 1], q = a - 2 * b + c; if (q < 0) d = clamp((0.5 * (a - c)) / q, -0.5, 0.5); }  return -80 + (160 * (bi + d)) / (N - 1);}

const HE = 0.24, HP = 0.5; // 面の高さ(波長に対する割合): 電場の山の高さ・強さの平均の高さ/** 電場 → 面の高さ。弱い所はそのまま比例し、強い所は tanh で縮める */const hOfE = (e: number, ed: number, lam: number) => HE * lam * Math.tanh(e / ed);/** 強さの平均(電場の 2 乗の平均) → 面の高さ。基準 er² の −25 dB を 0、+5 dB を 1 にして、そのあいだを dB に比例させる */const hOfP = (p: number, er2: number, lam: number) => HP * lam * clamp((10 * Math.log10(p / er2 + 1e-12) + 25) / 30, 0, 1.2); function teardown(s: S) {  s.root.traverse((o) => { const m = o as THREE.Mesh; if (m.geometry) m.geometry.dispose(); });  s.v.scene.remove(s.root);  s.root = new THREE.Group(); s.v.scene.add(s.root);  s.dip = null; s.fd = null;}/** 光: 波の斜面が読めるよう、太陽は奥の左から低めに当て、空からの光は弱くする */function setShadow(v: View3, half: number, dir: [number, number, number], sky: number) {  const c = v.sun.shadow.camera;  c.left = -half; c.right = half; c.top = half; c.bottom = -half; c.near = 0.05; c.far = half * 8; c.updateProjectionMatrix();  v.sun.target.position.set(0, 0, 0); v.sun.target.updateMatrixWorld();  v.sun.position.copy(new THREE.Vector3(...dir).normalize().multiplyScalar(half * 3));  v.sky.intensity = sky;} /** 1 本のアンテナの場面: 縦の面(電気力線)と、その奥の水平な半円(縦向きの電場の高さ)・放射の強さの形 */function buildDip(api: Api, s: S, donut: boolean) {  const d = dipole(F_DIP, P_DIP), lam = d.lam, q = api.quality;  const R = 3 * lam, Y = R, ni = [64, 100, 140][q], nj = ni;  const cut = new THREE.Group(); s.root.add(cut);  // アンテナ: 長さ 0.15 波長の 2 本の棒(真ん中のすき間から電流を送る)  const al = 0.15 * lam, rod = new THREE.CylinderGeometry(0.0022, 0.0022, al / 2 - 0.002, 12);  for (const sg of [1, -1]) { const m = new THREE.Mesh(rod, s.mats.ant); m.position.y = sg * (al / 4 + 0.001); m.castShadow = true; cut.add(m); }  let surf: THREE.Mesh | null = null;  const nr = [50, 80, 120][q], na = [24, 32, 48][q], r0 = 0.012;  if (!donut) {    const geo = new THREE.BufferGeometry(), nv = (nr + 1) * (na + 1);    geo.setAttribute("position", new THREE.BufferAttribute(new Float32Array(nv * 3), 3));    geo.setAttribute("normal", new THREE.BufferAttribute(new Float32Array(nv * 3), 3));    const idx: number[] = [];    for (let a = 0; a < na; a++) for (let r = 0; r < nr; r++) { const p = a * (nr + 1) + r; idx.push(p, p + 1, p + nr + 1, p + 1, p + nr + 2, p + nr + 1); }    geo.setIndex(idx);    surf = new THREE.Mesh(geo, s.mats.surf); surf.receiveShadow = true; surf.frustumCulled = false; cut.add(surf);  } else {    // 放射の強さの形: 中心からの距離を sin²θ に比例させた回転体の、奥の半分。切り口にはふたをする    const R0 = 2.4 * lam, nt = 64, nf = 24, pos: number[] = [], ix: number[] = [];    for (let f = 0; f <= nf; f++) for (let t = 0; t <= nt; t++) {      const th = (Math.PI * t) / nt, ph = (Math.PI * f) / nf, r = R0 * Math.sin(th) ** 2;      pos.push(r * Math.sin(th) * Math.cos(ph), r * Math.cos(th), -r * Math.sin(th) * Math.sin(ph));    }    for (let f = 0; f < nf; f++) for (let t = 0; t < nt; t++) { const p = f * (nt + 1) + t; ix.push(p, p + nt + 1, p + 1, p + 1, p + nt + 1, p + nt + 2); }    const g = new THREE.BufferGeometry(); g.setAttribute("position", new THREE.Float32BufferAttribute(pos, 3)); g.setIndex(ix); g.computeVertexNormals();    const shell = new THREE.Mesh(g, s.mats.donut); shell.castShadow = true; shell.receiveShadow = true; cut.add(shell);    const cp: number[] = [0, 0, 0], ci: number[] = [];    for (const sg of [1, -1]) {      const b = cp.length / 3;      for (let t = 0; t <= nt; t++) { const th = (Math.PI * t) / nt, r = R0 * Math.sin(th) ** 2; cp.push(sg * r * Math.sin(th), r * Math.cos(th), 0); }      for (let t = 0; t < nt; t++) ci.push(0, b + t, b + t + 1);    }    const cg = new THREE.BufferGeometry(); cg.setAttribute("position", new THREE.Float32BufferAttribute(cp, 3)); cg.setIndex(ci); cg.computeVertexNormals();    cut.add(new THREE.Mesh(cg, s.mats.cap));  }  s.dip = { d, R, Y, ni, nj, ps: new Float32Array((ni + 1) * (nj + 1)), seg: new Float32Array(4 * 60000), nseg: 0, cut, surf, nr, na, r0, hprof: new Float32Array(nr + 1), donut };  setShadow(s.v, R * 1.2, [0.4, 0.8, 0.6], 0.35);  s.v.orbit.minDist = 0.25; s.v.orbit.maxDist = 4;}/** 水平な半円の高さ(真横の面の縦向きの電場は、距離 ρ だけで決まる) */function updateDipSurface(api: Api, s: S, wt: number) {  const D = s.dip!, surf = D.surf; if (!surf) return;  const { nr, na, r0, R, d } = D, H = D.hprof, o = s.tmp, show = api.str("show"), cw = Math.cos(wt), sw = Math.sin(wt);  const ed = (d.K * d.k * d.k) / (1.5 * d.lam), lam = d.lam;  for (let r = 0; r <= nr; r++) {    const rho = r0 + ((R - r0) * r) / nr;    dipField(d, rho, 0, 1, o);    // 真横では θ の向きが真下なので、上向きの電場は −Eθ    H[r] = 2.4 * (show === "p" ? hOfP((o[2] * o[2] + o[3] * o[3]) / 2, (ed * ed) / 2, lam) : hOfE(-(o[2] * cw + o[3] * sw), ed, lam));  }  const pos = surf.geometry.attributes.position as THREE.BufferAttribute, nrm = surf.geometry.attributes.normal as THREE.BufferAttribute;  const P = pos.array as Float32Array, N = nrm.array as Float32Array, dr = (R - r0) / nr;  for (let a = 0; a <= na; a++) {    const ph = (Math.PI * a) / na, c = Math.cos(ph), sn = Math.sin(ph);    for (let r = 0; r <= nr; r++) {      const v = 3 * (a * (nr + 1) + r), rho = r0 + dr * r, sl = (H[Math.min(nr, r + 1)] - H[Math.max(0, r - 1)]) / (dr * (r > 0 && r < nr ? 2 : 1)), L = Math.hypot(1, sl);      P[v] = rho * c; P[v + 1] = H[r]; P[v + 2] = -rho * sn;      N[v] = (-sl * c) / L; N[v + 1] = 1 / L; N[v + 2] = (sl * sn) / L;    }  }  pos.needsUpdate = true; nrm.needsUpdate = true;}/** 画面に、電気力線・磁場の向き・観測点を描く */function drawDip(api: Api, s: S, wt: number, U: THREE.Vector3) {  const D = s.dip!, H = s.H, ctx = api.ctx, p = s.pt, q = [0, 0];  planeMap(s.v, U, H);  contours(D, wt);  // 4 つの象限に写して 1 本の道にまとめ、地の色の太い線の上に線の色の細い線を重ねる(網点の上でも読める)  ctx.beginPath();  for (let n = 0; n < D.nseg; n++) {    const o = 4 * n;    for (const sx of [1, -1]) for (const sy of [1, -1]) {      if (!mapPt(api, H, sx * D.seg[o], sy * D.seg[o + 1], p) || !mapPt(api, H, sx * D.seg[o + 2], sy * D.seg[o + 3], q)) continue;      ctx.moveTo(p[0], p[1]); ctx.lineTo(q[0], q[1]);    }  }  ctx.lineCap = "round"; ctx.setLineDash([]);  ctx.strokeStyle = api.paper; ctx.lineWidth = api.u(3.4); ctx.stroke();  ctx.strokeStyle = api.ink; ctx.lineWidth = api.u(1.3); ctx.stroke();  ctx.lineCap = "butt";  const o = s.tmp, cw = Math.cos(wt), sw = Math.sin(wt), lam = D.d.lam;  // 縦の面の縁(点線の円)  ctx.beginPath();  for (let n = 0; n <= 96; n++) { const a = (TAU * n) / 96; if (mapPt(api, H, D.R * Math.cos(a), D.R * Math.sin(a), p)) { if (n) ctx.lineTo(p[0], p[1]); else ctx.moveTo(p[0], p[1]); } }  ctx.strokeStyle = api.ink; ctx.lineWidth = api.u(1); ctx.setLineDash([api.u(2), api.u(5)]); ctx.stroke(); ctx.setLineDash([]);  // 磁場: 縦の面に垂直なので、手前へ向く所を ⊙、奥へ向く所を ⊗ で描く(アンテナの真横の線の上)。丸の大きさは、遠くの振幅に対する強さ  if (api.bool("bmark") && !D.donut) {    const gs = lam / 4;    for (let sx = -D.R + gs / 2; sx < D.R; sx += gs) {      const r = Math.abs(sx);      if (r < 0.3 * lam) continue;      dipField(D.d, r, 0, 1, o);      const B = o[4] * cw + o[5] * sw, rad = api.u(5.5) * Math.min(1, Math.abs(B) / (D.d.KB / r)) ** 0.7;      if (rad < api.u(1.4) || !mapPt(api, H, sx, 0, p)) continue;      bSym(api, p[0], p[1], rad, (B > 0) === (sx > 0));    }  }  // 観測点: 電場の矢印(長さは、その場所の振幅に対する今の強さ)と、磁場の記号  const [ps, py] = s.probe, r = Math.hypot(ps, py), st = Math.abs(ps) / r, ct = py / r, sg = Math.sign(ps) || 1;  dipField(D.d, r, ct, st, o);  const esR = sg * (o[0] * st + o[2] * ct), esI = sg * (o[1] * st + o[3] * ct), eyR = o[0] * ct - o[2] * st, eyI = o[1] * ct - o[3] * st;  const A = esR * esR + esI * esI + eyR * eyR + eyI * eyI, cR = esR * esR - esI * esI + eyR * eyR - eyI * eyI, cI = 2 * (esR * esI + eyR * eyI);  const eAmp = Math.sqrt((A + Math.hypot(cR, cI)) / 2), bAmp = Math.hypot(o[4], o[5]);  const es = esR * cw + esI * sw, ey = eyR * cw + eyI * sw, B = o[4] * cw + o[5] * sw, L = 0.55 * lam;  if (mapPt(api, H, ps, py, p)) {    api.circle(p[0], p[1], api.u(15), { w: api.u(1.3), dash: [api.u(2), api.u(3)], stroke: "ink" });    api.label("観測点", p[0] - api.u(20), p[1] + api.u(4), { size: api.u(11), align: "right" });    if (mapPt(api, H, ps + (L * es) / eAmp, py + (L * ey) / eAmp, q) && Math.hypot(q[0] - p[0], q[1] - p[1]) > api.u(4)) {      api.arrow(p[0], p[1], q[0], q[1], { w: api.u(4.5), stroke: "paper", head: api.u(13) });      api.arrow(p[0], p[1], q[0], q[1], { w: api.u(2.2), head: api.u(11) });      api.label("E", q[0] + api.u(6), q[1] - api.u(4), { size: api.u(13), font: "math" });    }    const rb = (api.u(10) * Math.abs(B)) / bAmp;    if (rb > api.u(2)) { bSym(api, p[0] + api.u(26), p[1] + api.u(22), rb, (B > 0) === (ps > 0)); api.label("B", p[0] + api.u(26) + rb + api.u(4), p[1] + api.u(26), { size: api.u(13), font: "math" }); }  }  // 目盛り(1 波長)と、アンテナの名前  const yb = -1.08 * D.R, sb = -lam / 2;  if (mapPt(api, H, sb, yb, p) && mapPt(api, H, sb + lam, yb, q)) {    api.line(p[0], p[1], q[0], q[1], { w: api.u(1.6) });    api.line(p[0], p[1] - api.u(5), p[0], p[1] + api.u(5), { w: api.u(1.6) }); api.line(q[0], q[1] - api.u(5), q[0], q[1] + api.u(5), { w: api.u(1.6) });    api.label(`1 波長 = ${fmt(lam * 100, 1)} cm`, q[0] + api.u(8), q[1] + api.u(4), { size: api.u(11.5) });  }  if (mapPt(api, H, 0, -0.1 * lam, p)) api.label("アンテナ", p[0], p[1] + api.u(22), { size: api.u(11.5), align: "center" });  return { r, st, ct, eAmp, bAmp, o, cR, cI };}/** 面に垂直な向きの記号。into = 奥へ(⊗)・そうでなければ手前へ(⊙) */function bSym(api: Api, x: number, y: number, r: number, into: boolean) {  api.circle(x, y, r + api.u(1.5), { stroke: false, fill: "paper" });  api.circle(x, y, r, { w: api.u(1.1) });  if (into) { const k = r * 0.62; api.line(x - k, y - k, x + k, y + k, { w: api.u(1.1) }); api.line(x - k, y + k, x + k, y - k, { w: api.u(1.1) }); }  else api.dot(x, y, Math.max(api.u(1), r * 0.28));} /** 平面の場面: 格子・面・壁・アンテナ・受信機を作る */function buildFd(api: Api, s: S, world: World) {  const g = new Grid(world.f, world.Lx, world.Lz, CPL[api.quality]);  for (const w of world.walls) g.paint(w.x1, w.z1, w.x2, w.z2, w.m, w.t, w.box);  g.setSources(world.ants, (api.num("phase") * Math.PI) / 180);  const ni = g.nx - 2 * NPML, nj = g.nz - 2 * NPML;  const surf = heightGrid(ni - 1, nj - 1, (ni - 1) * g.dx, (nj - 1) * g.dx, s.mats.surf);  surf.mesh.frustumCulled = false; s.root.add(surf.mesh);  const ah = world.wh * 1.3;  const antGeo = new THREE.CylinderGeometry(0.004, 0.004, ah, 10); antGeo.translate(0, ah / 2 - 0.005, 0);  const rxGeo = new THREE.BoxGeometry(0.04, 0.075, 0.01); rxGeo.translate(0, 0.0375, 0);  const rx = new THREE.Mesh(rxGeo, s.mats.rx); rx.castShadow = true; s.root.add(rx);  s.fd = { g, world, surf, hbuf: new Float32Array(ni * nj), ni, nj, walls: [], ants: [], rx, wallsDirty: true, meshDirty: true, lastShow: "", antGeo, antMat: s.mats.ant, beam: NaN, beamAt: -1 };  syncAnts(s);  const big = Math.max(world.Lx, world.Lz);  setShadow(s.v, big * 0.7, [-0.45, 0.55, -0.6], 0.2);  s.v.orbit.minDist = 0.4; s.v.orbit.maxDist = big * 4;}function syncAnts(s: S) {  const F = s.fd!;  while (F.ants.length < F.world.ants.length) { const m = new THREE.Mesh(F.antGeo, F.antMat); m.castShadow = true; s.root.add(m); F.ants.push(m); }  while (F.ants.length > F.world.ants.length) s.root.remove(F.ants.pop()!);  F.world.ants.forEach((a, i) => F.ants[i].position.set(a.x, 0, a.z));  F.rx.position.set(F.world.rx[0], 0, F.world.rx[1]);}/** 格子の材質から壁の立体を作り直す。同じ材質が横に並んだ列を 1 つの箱にし、上下の行で同じ列なら箱をつなぐ */function buildWalls(s: S) {  const F = s.fd!, g = F.g, nx = g.nx, dx = g.dx;  for (const m of F.walls) { s.root.remove(m); m.geometry.dispose(); }  F.walls = [];  for (let mi = 1; mi < MKEY.length; mi++) {    const boxes: number[] = []; // i0, i1, j0, j1    let open = new Map<number, number>(); // i0·4096 + i1 → boxes の中の位置    for (let j = NPML; j < g.nz - NPML; j++) {      const next = new Map<number, number>();      for (let i = NPML; i < nx - NPML; i++) {        if (g.mid[j * nx + i] !== mi) continue;        let e = i; while (e + 1 < nx - NPML && g.mid[j * nx + e + 1] === mi) e++;        const key = i * 4096 + e, b = open.get(key);        if (b !== undefined) { boxes[b + 3] = j; next.set(key, b); } else { next.set(key, boxes.length); boxes.push(i, e, j, j); }        i = e;      }      open = next;    }    if (!boxes.length) continue;    const nb = boxes.length / 4, P = new Float32Array(nb * 60), Nn = new Float32Array(nb * 60), I = new Uint32Array(nb * 30);    const top = F.world.wh * (mi === 3 ? 0.8 : 1);    for (let b = 0; b < nb; b++) {      const x0 = g.x0 + (boxes[4 * b] - 0.5) * dx, x1 = g.x0 + (boxes[4 * b + 1] + 0.5) * dx, z0 = g.z0 + (boxes[4 * b + 2] - 0.5) * dx, z1 = g.z0 + (boxes[4 * b + 3] + 0.5) * dx, y0 = -0.02;      // 上の面と 4 つの側面(下の面は見えないので作らない)      const faces: [number[], number[]][] = [        [[x0, top, z0, x0, top, z1, x1, top, z1, x1, top, z0], [0, 1, 0]],        [[x0, y0, z1, x1, y0, z1, x1, top, z1, x0, top, z1], [0, 0, 1]],        [[x1, y0, z0, x0, y0, z0, x0, top, z0, x1, top, z0], [0, 0, -1]],        [[x1, y0, z1, x1, y0, z0, x1, top, z0, x1, top, z1], [1, 0, 0]],        [[x0, y0, z0, x0, y0, z1, x0, top, z1, x0, top, z0], [-1, 0, 0]],      ];      faces.forEach(([v, n], f) => {        const o = b * 20 + f * 4;        for (let k = 0; k < 4; k++) { P.set(v.slice(3 * k, 3 * k + 3), 3 * (o + k)); Nn.set(n, 3 * (o + k)); }        I.set([o, o + 1, o + 2, o, o + 2, o + 3], b * 30 + f * 6);      });    }    const geo = new THREE.BufferGeometry();    geo.setAttribute("position", new THREE.BufferAttribute(P, 3)); geo.setAttribute("normal", new THREE.BufferAttribute(Nn, 3)); geo.setIndex(new THREE.BufferAttribute(I, 1));    const mesh = new THREE.Mesh(geo, s.mats[MKEY[mi]]); mesh.castShadow = true; mesh.receiveShadow = true; s.root.add(mesh); F.walls.push(mesh);  }}/** 高さの目安にする電場: 壁が無いとき、1 本のアンテナから 0.6 m の所の振幅 */const edOf = (g: Grid) => ((g.w * MU0 * I_ANT) / 4) * Math.sqrt(2 / (Math.PI * g.k * 0.6));/** 面の高さと向き(光の当たり方)を、格子の電場か強さの平均から作る。向きは、となりの点との高さの差から求める */function updateFdSurface(api: Api, s: S) {  const F = s.fd!, g = F.g, { ni, nj } = F, Hb = F.hbuf, show = api.str("show"), lam = C0 / g.f;  const ed = edOf(g), er2 = (ed * ed) / 2;  for (let j = 0, v = 0; j < nj; j++) for (let i = 0, c = (j + NPML) * g.nx + NPML; i < ni; i++, c++, v++) Hb[v] = g.mid[c] ? 0 : show === "p" ? hOfP(g.avg[c], er2, lam) : hOfE(g.ey[c], ed, lam);  const pos = F.surf.geo.attributes.position as THREE.BufferAttribute, nrm = F.surf.geo.attributes.normal as THREE.BufferAttribute;  const P = pos.array as Float32Array, N = nrm.array as Float32Array, k = 1 / (2 * g.dx);  for (let j = 0, v = 0; j < nj; j++) for (let i = 0; i < ni; i++, v++) {    const hx = (Hb[i < ni - 1 ? v + 1 : v] - Hb[i > 0 ? v - 1 : v]) * k, hz = (Hb[j < nj - 1 ? v + ni : v] - Hb[j > 0 ? v - ni : v]) * k, L = 1 / Math.sqrt(1 + hx * hx + hz * hz);    P[3 * v + 1] = Hb[v]; N[3 * v] = -hx * L; N[3 * v + 1] = L; N[3 * v + 2] = -hz * L;  }  pos.needsUpdate = true; nrm.needsUpdate = true;}function drawFd(api: Api, s: S) {  const F = s.fd!, W = F.world, v = s.v, lab = { size: api.u(11.5), align: "center" as const };  const pr = (x: number, y: number, z: number) => v.project([x, y, z]);  for (const [x, z, t] of W.labels) { const p = pr(x, W.wh * 1.05, z); if (p.on) api.label(t, p.x, p.y - api.u(8), lab); }  // アンテナの取っ手(四角)と受信機の取っ手(点線の輪)  const hs = api.u(W.ants.length > 4 ? 5 : 8);  for (const a of W.ants) { const p = pr(a.x, 0, a.z); if (p.on) api.rect(p.x - hs, p.y - hs, 2 * hs, 2 * hs, { w: api.u(1.4) }); }  const pRx = pr(W.rx[0], 0.04, W.rx[1]);  if (pRx.on) { api.circle(pRx.x, pRx.y, api.u(16), { w: api.u(1.3), dash: [api.u(2), api.u(3)] }); api.label("受信機", pRx.x, pRx.y + api.u(32), lab); }  // 目盛り 1 m(面の手前の縁)  const z = W.Lz / 2 - 0.06, x0 = -W.Lx / 2 + 0.1, a = pr(x0, 0, z), b = pr(x0 + 1, 0, z);  if (a.on && b.on) { api.line(a.x, a.y, b.x, b.y, { w: api.u(1.8) }); api.label("1 m", (a.x + b.x) / 2, (a.y + b.y) / 2 + api.u(18), lab); }  // 並べたアンテナの電波の向き: 測った向き(実線)と式の向き(破線)  if (W.beam) {    const B = W.beam, c = pr(B.cx, 0.01, B.cz), lam = C0 / W.f, sn = clamp(((api.num("phase") / 360) * lam) / B.d, -1, 1), th = Math.asin(sn);    const end = (t: number) => pr(B.cx + B.R * Math.sin(t), 0.01, B.cz - B.R * Math.cos(t));    const e1 = end(th);    if (c.on && e1.on) api.line(c.x, c.y, e1.x, e1.y, { w: api.u(1.6), dash: [api.u(6), api.u(5)] });    if (Number.isFinite(F.beam)) { const e2 = end((F.beam * Math.PI) / 180); if (c.on && e2.on) { api.line(c.x, c.y, e2.x, e2.y, { w: api.u(2.4), stroke: "paper" }); api.line(c.x, c.y, e2.x, e2.y, { w: api.u(1.4) }); api.label("測った向き", e2.x, e2.y - api.u(10), lab); } }  }}

const hitDist = (api: Api, x: number, y: number) => Math.hypot(api.pointer.x - x, api.pointer.y - y);/** 指の下の、水平な面(高さ 0)の点。格子の中(吸収層の外)なら返す */function floorAt(api: Api, s: S): [number, number] | null {  const W = s.world!, p = s.v.onPlane(api.pointer.x, api.pointer.y, new THREE.Plane(new THREE.Vector3(0, 1, 0), 0));  if (!p || Math.abs(p.x) > W.Lx / 2 || Math.abs(p.z) > W.Lz / 2) return null;  return [p.x, p.z];}/** 観測点を、指の下の縦の面の点へ動かす(アンテナに近づきすぎない) */function moveProbe(api: Api, s: S) {  const D = s.dip!, rot = D.cut.rotation.y, n = new THREE.Vector3(Math.sin(rot), 0, Math.cos(rot));  const p = s.v.onPlane(api.pointer.x, api.pointer.y, new THREE.Plane(n, 0));  if (!p) return;  let ps = clamp(p.x * Math.cos(rot) - p.z * Math.sin(rot), -D.R, D.R), py = clamp(p.y, -D.Y, D.Y);  const r = Math.hypot(ps, py), rmin = 0.08 * D.d.lam;  if (r < rmin) { if (r > 1e-9) { ps *= rmin / r; py *= rmin / r; } else ps = rmin; }  s.probe = [ps, py];}/** 指でなぞった線分に壁を塗る(「消す」なら空気に戻す) */function paintTo(s: S, x: number, z: number, key: string) {  const F = s.fd!, d = s.drag!, m = key === "erase" ? "" : key, t = key === "erase" ? 0.14 : undefined;  const w: Wall = { x1: d.lx, z1: d.lz, x2: x, z2: z, m, t };  if (F.g.paint(w.x1, w.z1, w.x2, w.z2, w.m, w.t)) F.wallsDirty = true;  s.world!.walls.push(w); d.lx = x; d.lz = z;}