Last active
November 22, 2023 18:33
-
-
Save Jademaster/be6aa20492a0a9ee48178eca13b60870 to your computer and use it in GitHub Desktop.
Simulates the three-body problem in Scala using doodle for animation and breeze for numerical integration.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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)) | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment