n body simulation

turtles-own
[
  mass
  ;; x and y components of velocity
  vx
  vy
  ;; x and y components of force due to gravity
  gx
  gy
]

globals[
  ;; possible states for system: unary, binary, other
  state
  ;; gravitational constant
  g
  framerate
]

;; sets up a system with one central mass

to setup-unary
  reset
  crt 1
  [
    setxy 0 0
    set mass 10 ^ center-mass-exp
    set size log mass 10
    pen-down
  ]
end 

;; sets up a binary system

to setup-binary
  reset
  set-default-shape turtles "circle"

  crt 1
  [
    setxy 0 binary-r
    set mass 10 ^ binary-mass-exp
    pen-down
  ]

  crt 1
  [
    setxy 0 (- binary-r)
    set mass 10 ^ binary-mass-exp
    pen-down
  ]

  ;; setup orbits
  ask turtles
  [
    set size log mass 10
    let center max-one-of other turtles [mass]
    let m1+m2 mass + [mass] of center
    let r distance center
    let v orbital-velocity center 1
    let theta (180 - towards center) mod 360  ; <----- dont know why this works
    set vx (- v) * cos theta / 2
    set vy (- v) * sin theta / 2
  ]
end 

;; calculate forces, then move turtles

to go
  calculate-forces
  move
  tick
end 

;; adds random turtle

to add-random
  crt 1
  [
    ;; use exponent for greater density near center
    setxy (random-sign * 2 ^ (2 + random-float 6.5)) (random-sign * 2 ^ (2 + random-float 6.5))
    set mass 10 ^ (random 9)
    set size log mass 10

    orbit-strongest
    pen-down
  ]
end 

;; call from turtle context, attempts to orbit the turtle with the greatest pull

to orbit-strongest
  let center orbital-center
  let m1+m2 mass + [mass] of center
  let r distance center
  let v orbital-velocity center eccentricity
  let theta (180 - towards center) mod 360  ; <----- dont know why this works
  set vx (- v) * cos theta
  set vy (- v) * sin theta
  ;; if add-target-velocity is set, add centers velocity components
  if add-target-velocity
  [
    set vx vx + [vx] of center
    set vy vy + [vy] of center
  ]
end 

;; adds b's mass and momentum to a and kills b

to consume [a b]
  ;; momentum components for a
  let mvx1 [mass] of a * [vx] of a
  let mvy1 [mass] of a * [vy] of a
  ;; momentum components for b
  let mvx2 [mass] of b * [vx] of b
  let mvy2 [mass] of b * [vy] of b
  ;; modify caller turtle
  ask a
  [
    set vx (mvx1 + mvx2) / (mass + [mass] of b)
    set vy (mvy1 + mvy2) / (mass + [mass] of b)
    set mass mass + [mass] of b
    set size log mass 10
  ]
  ask b [ die ]
end 

;; iterate over every pair of turtles to sum the force due to gravity on each turtle

to calculate-forces
  ask turtles
  [
    let i 0
    let temp-gx 0
    let temp-gy 0
    let current self
    ask other turtles
    [
      ;; if close enough to a more massive turtle, be consumed
      ;; otherwise perform usual calculations
      ifelse distance current < [size] of current / 3 and
        [mass] of current >= mass
        [ consume current self ]
        [
          let m1m2 mass * [mass] of current
          let force g * m1m2 / (distance current ^ 2)
          let theta angle current self
          set temp-gx temp-gx + force * cos theta
          set temp-gy temp-gy + force * sin theta
        ]
    ]
    set gx temp-gx
    set gy temp-gy
  ]
end 

;; update turtle positions

to move
  ask turtles
  [
    let ax gx / mass
    let ay gy / mass
    set vx vx + ax
    set vy vy + ay
    let new-x xcor + (vx / 30) ; divide velocities by framerate
    let new-y ycor + (vy / 30)
    set gx 0
    set gy 0
    ifelse abs new-x > max-pxcor or abs new-y > max-pycor
      [ die ]
      [ setxy new-x new-y ]
  ]
end 

;; returns the angle between two turtles

to-report angle [turtle1 turtle2]
  let x1 [xcor] of turtle1
  let x2 [xcor] of turtle2
  let y1 [ycor] of turtle1
  let y2 [ycor] of turtle2
  let h atan (x2 - x1) (y2 - y1)
  report (90 - h) mod 360
end 

;; find turtle with greatest pull

to-report orbital-center
  let this self
  let center 0
  let maximum 0
  ask other turtles
  [
    let value mass / (distance this ^ 2)
    if value > maximum
      [
        set maximum value
        set center who
      ]
  ]
  report turtle center
end 

;; calculates rough velocity to orbit center with eccentricity ecc

to-report orbital-velocity [center ecc]
  let m1+m2 mass + [mass] of center
  let r distance center
  report ecc * sqrt (g * m1+m2 / r) * sqrt 30 ; needed to account for division by framerate
end 

;; returns 1 or -1 with equal probability

to-report random-sign
  report (-1) ^ random 2
end 

to reset
  clear-all
  reset-ticks
  set g 10 ^ g-exp
end