N体問題をブラウザで — 万有引力シミュレーションの実装

type CelestialBody = {
  x: number
  y: number
  vx: number
  vy: number
  mass: number
  radius: number
  color: string
  trail: { x: number; y: number }[]
}

const G = 6.674 // スケーリング済みの重力定数

function computeGravity(bodies: CelestialBody[]) {
  for (let i = 0; i < bodies.length; i++) {
    for (let j = i + 1; j < bodies.length; j++) {
      const a = bodies[i]
      const b = bodies[j]

      const dx = b.x - a.x
      const dy = b.y - a.y
      const distSq = dx * dx + dy * dy
      const dist = Math.sqrt(distSq)

      // ソフトニング: 距離が近すぎると力が発散する
      const softening = 10
      const force = (G * a.mass * b.mass) / (distSq + softening * softening)

      const fx = force * dx / dist
      const fy = force * dy / dist

      a.vx += fx / a.mass
      a.vy += fy / a.mass
      b.vx -= fx / b.mass
      b.vy -= fy / b.mass
    }
  }
}

function integrateVerlet(bodies: CelestialBody[], dt: number) {
  // 半ステップの速度更新（加速度は前ステップの値）
  for (const b of bodies) {
    b.vx += b.ax * dt * 0.5
    b.vy += b.ay * dt * 0.5
  }

  // 位置更新
  for (const b of bodies) {
    b.x += b.vx * dt
    b.y += b.vy * dt
  }

  // 新しい加速度を計算
  computeGravity(bodies)

  // 残り半ステップの速度更新
  for (const b of bodies) {
    b.vx += b.ax * dt * 0.5
    b.vy += b.ay * dt * 0.5
  }
}

const SUN: CelestialBody = {
  x: 400, y: 300,
  vx: 0, vy: 0,
  mass: 10000,
  radius: 20,
  color: '#FFD700',
  trail: [],
}

// 円軌道の初速 v = sqrt(GM/r)
function orbitalVelocity(centralMass: number, distance: number): number {
  return Math.sqrt(G * centralMass / distance)
}

const EARTH: CelestialBody = {
  x: 400 + 150, y: 300,
  vx: 0, vy: orbitalVelocity(10000, 150),
  mass: 10,
  radius: 6,
  color: '#4169E1',
  trail: [],
}

function drawTrail(ctx: CanvasRenderingContext2D, trail: { x: number; y: number }[]) {
  if (trail.length < 2) return

  ctx.beginPath()
  ctx.moveTo(trail[0].x, trail[0].y)

  for (let i = 1; i < trail.length; i++) {
    ctx.lineTo(trail[i].x, trail[i].y)
  }

  ctx.strokeStyle = 'rgba(255, 255, 255, 0.3)'
  ctx.lineWidth = 1
  ctx.stroke()
}

function updateTrail(body: CelestialBody, maxLength: number) {
  body.trail.push({ x: body.x, y: body.y })
  if (body.trail.length > maxLength) {
    body.trail.shift()
  }
}

function attractToPointer(
  particles: CelestialBody[],
  pointerX: number,
  pointerY: number,
  strength: number,
) {
  for (const p of particles) {
    const dx = pointerX - p.x
    const dy = pointerY - p.y
    const distSq = dx * dx + dy * dy
    const dist = Math.sqrt(distSq)

    if (dist < 1) continue

    const force = strength / (distSq + 100)
    p.vx += (force * dx) / dist
    p.vy += (force * dy) / dist
  }
}