Simulates the three-body problem in Scala using doodle for animation and breeze for numerical integration.

Gravity.scala

import breeze.linalg._
import breeze.linalg.LU.primitive.LU_DM_Impl_Double
import doodle.core._
import doodle.java2d._
import doodle.syntax.all._
import cats.effect.unsafe.implicits.global
import doodle.image._
import doodle.image.syntax._
import doodle.image.syntax.all.ImageOps
import doodle.image.syntax.image.ImageOps
import doodle.image.syntax.core._
import doodle.interact.animation._
import doodle.interact.syntax._
import doodle.reactor.Reactor
import doodle.core.PathElement._
import doodle.core.OpenPath._
import doodle.random._


@main def gravity: Unit =
  //conversion between vects and points
  def pointToVect(p : Point): breeze.linalg.DenseVector[Double] = DenseVector(p.x,p.y)
  def vectToPoint(v : breeze.linalg.DenseVector[Double]): Point = Point(v(0),v(1))
  
  //random angle
  val randomAngle: Random[Angle] =
  Random.double.map(x => x.turns)
  
  //random colors
  def randomColor(s: Double, l: Double): Random[Color] =
  randomAngle.map(hue => Color.hsl(hue, s, l))
  
  //random circles
  def randomCircle(r: Double, color: Random[Color]): Random[Image] =
  color.map(fill => Image.circle(r).fillColor(fill))
 
  //the length of a vector
  def norm( x : DenseVector[Double]) : Double = math.sqrt(math.pow(x(0),2) + math.pow(x(1),2))
  
  //ode for newton's law of gravitation  
  //inp = (x y vx vy), x' = vx, y' = vy, vx' = -k x/(x^2 +y^2), vy'= -ky/(x^2 + y^2)  
  def dxdt(inp : breeze.linalg.DenseVector[Double]): breeze.linalg.DenseVector[Double] = DenseVector(inp(2),inp(3),-inp(0)/norm(inp(0 to 1)),-inp(1)/norm(inp(0 to 1)))
  
  //for 3 bodies
  val nbodies: Int = 3
  //inp = (x1 y1 x2 y2 x3 y3 vx1 vy1 vx2 vy2 vx3 vy3) 
  def force(inp: breeze.linalg.DenseVector[Double]): breeze.linalg.DenseVector[Double] = DenseVector(-inp(0)/norm(inp), -inp(1)/norm(inp))
  
  def d3dt(inp : breeze.linalg.DenseVector[Double]): breeze.linalg.DenseVector[Double] = 
    val out = DenseVector.zeros[Double](12)
    val p1 = inp(0 to 1)
    val p2 = inp(2 to 3)
    val p3 = inp(4 to 5)
    //deriv of position is velocity
    out(0 to 1) := inp(6 to 7)
    out(2 to 3) := inp(8 to 9)
    out(4 to 5) := inp(10 to 11)
    // gravitation force
    out(6 to 7) := force(p1-p2) + force(p1-p3)
    out(8 to 9) := force(p2-p1) + force(p2-p3)
    out(10 to 11) := force(p3-p1) + force(p3-p2)
    out
    
  val init=(DenseVector(100.0,50.0,-200.0,-50.0,-100.0,1.0,0.0,1.0,0.0,1.0,-2.0,1.0),Color.crimson,Color.green,Color.white)
  

  val stepsize: Double = 0.8
  
  //euler's method of numerical integration
  def evolve(state : breeze.linalg.DenseVector[Double]) = state + (DenseVector.fill(12){stepsize} *:* d3dt(state))
  
  //
  def changestate(lst : List[(breeze.linalg.DenseVector[Double],Color,Color,Color)]): List[(breeze.linalg.DenseVector[Double],Color,Color,Color)] =
   lst match {
    case Nil => Nil
    case (x, c1,c2,c3) :: tail =>  (evolve(x),randomColor(0.5,0.5).run,randomColor(0.5,0.5).run,randomColor(0.5,0.5).run) +: lst
  }
  
  def drawpth(l:List[(breeze.linalg.DenseVector[Double],Color,Color,Color)]): Image = l match {
    case Nil => Image.empty
    case (x, c1,c2,c3) :: tail => Image.circle(10.0).fillColor(c1).at(vectToPoint(x(0 to 1))).on(
        Image.circle(10.0).fillColor(c2).at(vectToPoint(x(2 to 3))).on(
            Image.circle(10.0).fillColor(c3).at(vectToPoint(x(4 to 5))).on(drawpth(tail))
        )
    )
  }

  val travellingCircle =
  Reactor.init(List(init))
    .withOnTick(ptlst => changestate(ptlst))
    .withRender(ptlst => drawpth(ptlst))
    .withStop(ptlst => ptlst.length >= 300)
  
  travellingCircle.run(Frame.default.withSize(1000,700).withCenterOnPicture.withBackground(Color.midnightBlue))