用 JavaScript 实现粒子系统和万有引力

export default class Vector2 {
  static zero = new Vector2(0, 0)
  static unit(x, y) {
    const u = new Vector2(x, y)
    const dvd = u.length()
    return u.divide(dvd)
  }
  static isNaV(vec) {
    return isNaN(vec.x) || isNaN(vec.y)
  }
  constructor (x, y) {
    this.x = x
    this.y = y
  }
  length() {
    return Math.sqrt(this.x * this.x + this.y * this.y)
  }
  normalize() {
    var inv = 1 / this.length()
    return new Vector2(this.x * inv, this.y * inv)
  }
  add(v) {
    return new Vector2(this.x + v.x, this.y + v.y)
  }
  subtract(v) {
    return new Vector2(this.x - v.x, this.y - v.y)
  }
  multiply(f) {
    return new Vector2(this.x * f, this.y * f)
  }
  divide(f) {
    var invf = 1 / f
    return new Vector2(this.x * invf, this.y * invf)
  }
  rotate(angle, center) {
    return new Vector2(
      (this.x - center.x) * Math.cos(angle) - (this.y - center.y) * Math.sin(angle) + center.x,
      (this.x - center.x) * Math.sin(angle) + (this.y - center.y) * Math.cos(angle) + center.y
    )
  }
}

import Vector2 from './vector2'

export default class Particle {
  static distancePow2(a, b) {
    return Math.pow(a.position.x - b.position.x, 2) + Math.pow(a.position.y - b.position.y, 2)
  }
  static massToRadius(mass) {
    return Math.pow(mass, 1 / 12)
  }
  static RocheLimitPow2(a, b) {
    return Math.pow((a.radius + b.radius) * 1.0134, 2)
  }
  constructor(position, velocity, color, mass, stasis = false, obWidth = 960, obHeight = 960) {
    this.position = position
    this.velocity = velocity
    this.acceleration = Vector2.zero
    this.color = color
    this.mass = mass
    this.radius = Particle.massToRadius(mass)

    this.stasis = stasis
    this.visible = true

    this.dead = false
    this.isSelected = false

    this.obWidth = obWidth
    this.obHeight = obHeight
  }
  devourOther(pOther) {
    if (pOther.mass < this.mass && !this.dead && !pOther.dead) {
      // m1 * v1 + m2 * v2 = (m1 + m2) * v'
      const newVelovity = this.velocity.multiply(this.mass).add(pOther.velocity.multiply(pOther.mass)).divide(this.mass + pOther.mass)
      this.velocity = newVelovity
      this.mass += pOther.mass
      this.radius = Particle.massToRadius(this.mass)
      this.color = this.color.add(pOther.color)
      pOther.dead = true
    }
  }
  outOfScreen() {
    return this.position.x + this.radius < 0 || this.position - this.radius > this.obWidth ||
      this.y + this.radius < 0 || this.y - this.radius > this.obHeight
  }

  render(ctx) {
    ctx.fillStyle = this.color.toRgba(1)

    ctx.beginPath()
    ctx.arc(this.position.x, this.position.y, this.radius, 0, Math.PI * 2, true)
    ctx.closePath()
    ctx.fill()
  }

  complexRender(ctx) {

    for (let i = 1; i <= Math.ceil(this.radius * 2); i += 1) {
      ctx.beginPath()
      ctx.arc(this.position.x, this.position.y, i, 0, Math.PI * 2)
      ctx.closePath()
      ctx.strokeStyle = this.color.toRgba(Math.random())
      ctx.stroke()
    }
  }

  renderPath() {
  }
}

import Particle from './particle'
import Vector2 from './vector2'

/**
 * 标准正态分布
 */
function standardNormalDistribution() {
  const numberPool = []
  return function () {
    if (numberPool.length > 0) {
      return numberPool.pop()
    }
    else {
      const u = Math.random(), v = Math.random()
      const p = Math.sqrt(-2 * Math.log(u)) * Math.cos(2 * Math.PI * v)
      const q = Math.sqrt(-2 * Math.log(u)) * Math.sin(2 * Math.PI * v)
      numberPool.push(q)
      return p
    }
  }()
}

/**
 * 正态分布
 * @param {number} off 期望分布顶点离开数轴中心的偏移量
 * @param {number} con 相对标准正态分布的系数
 */
function NormalDistribution(off, con) {
  const standard = standardNormalDistribution()
  return standard * con + off
}

export default class ParticleSystem {
  static G = 6.67408

  constructor(ctx, w, h) {
    this.w = w || ctx.canvas.width / window.devicePixelRatio
    this.h = h || ctx.canvas.height / window.devicePixelRatio

    this.workers = []
    /** @type {Particle[]} */
    this.particles = []
    /** @type {CanvasRenderingContext2D} */
    this.context = ctx
    this.effectors = []

    this.pauseSignal = false
    this.disableDevour = false
    this.disableGravitation = false

    this.lastRenderTime = new Float64Array(512)
    this.renderInterval = 0
  }

  emit(particle) {
    this.particles.push(particle)
  }
  emitMess(count, centerMass) {
    const __width = this.w
    const __height = this.h

    const O = new Vector2(__width / 2, __height / 2)
    const R = Math.min(__height, __width) / 2 - 50

    const ct = new Particle(O, Vector2.zero, Color.red, centerMass, true)
    ct.visible = false
    this.emit(ct)

    for (let i = 0; i < count; i++) {
      const randomX = random(__width / 2 - R, __width / 2 + R)
      const randomY = random(__height / 2 - Math.pow(R * R - Math.pow(randomX - __width / 2, 2), 0.5), Math.pow(R * R - Math.pow(randomX - __width / 2, 2), 0.5) + __height / 2)
      const P = new Vector2(randomX, randomY)
      this.emit(new Particle(
        P,
        O.subtract(P).rotate(Math.PI / 2, Vector2.zero).normalize().multiply(NormalDistribution(400, 200)),
        Color.random(),
        random(10, 1000)
      ))
    }
  }

  remove(index) {
    if (!this.particles[index] || !this.particles[index].dead) {
      return
    }
    this.particles[index] = this.particles[this.particles.length - 1]
    this.particles.pop()
  }

  kinematics(dt) {
    for (const particle of this.particles) {
      if (particle.stasis) {
        particle.acceleration = particle.velocity = Vector2.zero
        continue
      }
      particle.position = particle.position.add(particle.velocity.multiply(dt))
      particle.velocity = particle.velocity.add(particle.acceleration.multiply(dt))
    }
  }

  devour() {
    this.particles.sort((a, b) => b.mass - a.mass).forEach((p, _tmp, particlesRef) => {

      const foodCandidateIndex = []

      particlesRef.filter((pOther, innerIndex) => {

        if (p.mass <= pOther.mass) return
        if (p === pOther) return

        const canEat = Particle.distancePow2(p, pOther) < Particle.RocheLimitPow2(p, pOther)

        if (canEat) {
          foodCandidateIndex.push(innerIndex)
        }

        return canEat
      }).forEach(foodP => p.devourOther(foodP))
      
      foodCandidateIndex.forEach(beEatenIndex => this.remove(beEatenIndex))
    })
    
    this.particles.forEach((p, index) => {
      if (p.dead) this.remove(index)
    })
  }

  universalGravitation() {
    this.particles.forEach((p) => {

      const totalGravitation = this.particles.reduce((pv, cv) => {
        if (p === cv) {
          return pv
        }
        else {
          const gravAcc = ParticleSystem.G * cv.mass / Particle.distancePow2(p, cv)
          const gravVec = Vector2.unit(cv.position.x - p.position.x, cv.position.y - p.position.y).multiply(gravAcc)
          return pv.add(gravVec)
        }
      }, Vector2.zero)

      p.acceleration = totalGravitation
    })
  }

  applyEffectors() {
    for (const effector of this.effectors) {
      const apply = effector.apply
      for (const p of this.particles) apply(p)
    }
  }

  simulate(dt) {
    if (this.pauseSignal) return

    this.applyEffectors()
    if (!this.disableDevour) this.devour()
    if (!this.disableGravitation) this.universalGravitation()
    this.kinematics(dt)
  }

  recordRenderTime() {
    const now = performance.now()
    const zeroIndex = this.lastRenderTime.findIndex(v => v === 0)
    const actualLength = zeroIndex === -1 ? this.lastRenderTime.length : zeroIndex + 1

    if (zeroIndex === -1) {
      this.lastRenderTime.set(this.lastRenderTime.subarray(1))
      this.lastRenderTime.set([now], this.lastRenderTime.length - 1)
    }
    else {
      this.lastRenderTime.set([now], zeroIndex)
    }

    if (actualLength < 100) {
      this.renderInterval = (now - this.lastRenderTime[0]) / actualLength
    }
    else {
      this.renderInterval = (now - this.lastRenderTime[actualLength - 100]) / 100
    }
  }

  render() {
    for (const p of this.particles) {
      if (p.outOfScreen() || !p.visible) continue
      p.render(this.context)
    }
  }

  complexRender() {
    for (const p of this.particles) {
      p.complexRender(this.context)
    }

    this.recordRenderTime()
  }
}

import Vector2 from './vector2'

export class ChamberBox {
  constructor(x1, y1, x2, y2, elasticCoefficient) {
    elasticCoefficient = elasticCoefficient || 1
    this.apply = particle => {
      if (particle.position.x - particle.radius < x1 || particle.position.x + particle.radius > x2) {
        particle.velocity.x = -1 * elasticCoefficient * particle.velocity.x
      }
      if (particle.position.y - particle.radius < y1 || particle.position.y + particle.radius > y2) {
        particle.velocity.y = -1 * elasticCoefficient * particle.velocity.y
      }
    }
  }
}

export class Field {
  /**
   * @param {number} dt
   * @param {(pos: Vector2) => Vector2} fieldAcc
   */
  constructor(dt, fieldAcc) {
    fieldAcc = fieldAcc || (() => new Vector2(0, 9.8))
    this.apply = particle => {
      particle.velocity = particle.velocity.add(fieldAcc(particle.position).multiply(dt))
    }
  }
}

export function random(lower, upper) {
  return lower + (upper - lower) * Math.random()
}

export class Color {
  static black = new Color(0, 0, 0)
  static white = new Color(1, 1, 1)
  static red = new Color(1, 0, 0)
  static green = new Color(0, 1, 0)
  static blue = new Color(0, 0, 1)
  static yellow = new Color(1, 1, 0)
  static cyan = new Color(0, 1, 1)
  static purple = new Color(1, 0, 1)
  static random() {
    return new Color(random(.15, .9), random(.1, .9), random(.2, .9))
  }
  constructor(r, g, b) {
    if (r > 1 || g > 1 || b > 1) {
      r = r / 255
      g = g / 255
      b = b / 255
    }
    this.r = r
    this.g = g
    this.b = b
  }
  add(c) {
    return new Color(Math.min((this.r + c.r) / 2, 1), Math.min((this.g + c.g) / 2, 1), Math.min((this.b + c.b) / 2, 1))
  }
  toRgba(alpha) {
    return `rgba(${Math.floor(this.r * 255)},${Math.floor(this.g * 255)},${Math.floor(this.b * 255)},${alpha})`
  }
  toReversedRgba(alpha) {
    return `rgba(${Math.floor((1 - this.r) * 255)},${Math.floor((1 - this.g) * 255)},${Math.floor((1 - this.b) * 255)},${alpha})`
  }
}