Gravitation

turtles-own [ x-velocity y-velocity z-velocity ; 3d velocity of the particle
  x-pos y-pos z-pos ; physical position of the particle
  mass ; represents physical mass. More mass means a greater gravitational pull
  ]
globals [
  zoom-factor ; used to relate x, y, z-pos variables to actual xcor, ycor, zcor
  ]

to setup
  clear-all
  set-default-shape turtles "circle"
  set zoom-factor 1
  create-turtles number-of-objects [
    setxyz (random-float starting-box - random-float starting-box) (random-float starting-box - random-float starting-box) (random-float starting-box - random-float starting-box)
    set x-velocity ((random-float max-starting-velocity) - (random-float max-starting-velocity))
    set y-velocity ((random-float max-starting-velocity) - (random-float max-starting-velocity))
    set z-velocity ((random-float max-starting-velocity) - (random-float max-starting-velocity))
    set mass random-float max-starting-mass
    set size mass ^ (1 / 3)

    set x-pos xcor
    set y-pos ycor
    set z-pos zcor
    ]
  reset-ticks
end 

;; Sets up a symmetric system of two particles orbiting each other

to setup-two-orbits
  clear-all
  set-default-shape turtles "circle"
  set zoom-factor 1
  create-turtles 1 [
    setxyz 5 0 0
    set x-velocity .03
    set y-velocity .1
    set z-velocity 0
    set mass 3
    set size mass ^ (1 / 3)
    set x-pos xcor
    set y-pos ycor
    set z-pos zcor
    ]
  create-turtles 1 [
    setxyz -5 0 0
    set x-velocity 0
    set y-velocity -.1
    set z-velocity 0
    set mass 3
    set size mass ^ (1 / 3)
    set x-pos xcor
    set y-pos ycor
    set z-pos zcor
    ]
  reset-ticks
end 

; Sets up a system with a "sun" and a smaller particle orbiting it

to setup-heliocentric-orbit
  clear-all
  set-default-shape turtles "circle"
  set zoom-factor 1
  create-turtles 1 [
    setxyz 0 0 0
    set x-velocity 0
    set y-velocity 0
    set z-velocity 0
    set mass 1000
    set size mass ^ (1 / 3)
    set x-pos xcor
    set y-pos ycor
    set z-pos zcor
    ]
  create-turtles 1 [
    setxyz -20 0 0
    set x-velocity 0
    set y-velocity 2
    set z-velocity 0
    set mass 1
    set size mass
    set x-pos xcor
    set y-pos ycor
    set z-pos zcor
  ]
  reset-ticks
end 

; Sets up a system of two pairs of particles - each particle orbits the other in its pair, and the pairs orbit each other

to setup-4
  clear-all
  set-default-shape turtles "circle"
  set zoom-factor 1
  create-turtles 1 [
    setxyz 23 5 0
    set x-velocity .07
    set y-velocity -.05
    set z-velocity 0
    set mass 2
    set size mass ^ (1 / 3)
    set x-pos xcor
    set y-pos ycor
    set z-pos zcor
  ]
  create-turtles 1 [
    setxyz 23 -5 0
    set x-velocity -.07
    set y-velocity -.05
    set z-velocity 0
    set mass 2
    set size mass ^ (1 / 3)
    set x-pos xcor
    set y-pos ycor
    set z-pos zcor
  ]
  create-turtles 1 [
    setxyz -23 5 0
    set x-velocity .07
    set y-velocity .05
    set z-velocity 0
    set mass 2
    set size mass ^ (1 / 3)
    set x-pos xcor
    set y-pos ycor
    set z-pos zcor
  ]
  create-turtles 1 [
    setxyz -23 -5 0
    set x-velocity -.07
    set y-velocity .05
    set z-velocity 0
    set mass 2
    set size mass ^ (1 / 3)
    set x-pos xcor
    set y-pos ycor
    set z-pos zcor
  ]
  reset-ticks
end 

to go
  compute-acceleration-changes
  check-boundaries
  move
  check-slide
  tick
end 

; Keeps the world centered on the center of mass of all particles

to check-slide
  let CoM center-of-mass turtles
  slide-view (item 0 CoM) (item 1 CoM) (item 2 CoM)
end 

; Calculates the center of mass of a turtleset

to-report center-of-mass [turtleset]
  ; R = 1/M sigma (mi * ri)
  let M 0
  let totalx 0
  let totaly 0
  let totalz 0
  ask turtleset [
    let mi [mass] of self
    set M M + mi
    set totalx totalx + [x-pos] of self * mi
    set totaly totaly + [y-pos] of self * mi
    set totalz totalz + [z-pos] of self * mi
  ]

  report (list (totalx / M) (totaly / M) (totalz / M))
end 

; Checks turtle positions in the world, zooms out if any turtle is too close to the edge

to check-boundaries
  let flag false
  ask turtles [
    if abs (x-pos + x-velocity) > max-pxcor * zoom-factor * .8 or
    abs (y-pos + y-velocity) > max-pycor * zoom-factor * .8 or
    abs (z-pos + z-velocity) > max-pzcor * zoom-factor * .8
    [
      set flag true
    ]
  ]

  if flag [
    print "zooming out"
    set-zoom 1.02
  ]

  set flag false
  ask turtles [
    if abs (x-pos + x-velocity) > max-pxcor * zoom-factor * .5 or
    abs (y-pos + y-velocity) > max-pycor * zoom-factor * .5 or
    abs (z-pos + z-velocity) > max-pzcor * zoom-factor * .5
    [
      set flag true
    ]
  ]
  if (not flag) [
    ;print "zooming in"
    ;set-zoom .98
  ]
end 

to tick-once
  compute-acceleration-changes
  check-boundaries
  move
  check-slide
  tick
end 

; Asks turtles to gravitationally affect each other

to compute-acceleration-changes
  ask turtles [
    let turtle1 self
    ask other turtles [
      gravitationally-affect turtle1 self
    ]
  ]
end 

; Update turtle positions based on velocity

to move
  ask turtles [
    set x-pos x-pos + x-velocity
    set y-pos y-pos + y-velocity
    set z-pos z-pos + z-velocity

    if not hidden? [
      set xcor x-pos / zoom-factor
      set ycor y-pos / zoom-factor
      set zcor z-pos / zoom-factor
    ]
  ]
end 

; Calculates gravitational pull of one turtle upon another

to gravitationally-affect [ affected affecting ]
  ; g = m1 * m2 / r^2
  ; r = sqrt(xdist^2 + ydist^2) -> r^2 = xdist^2 + ydist^2
  let xdist [x-pos] of affecting - [x-pos] of affected
  ;print xdist
  let ydist [y-pos] of affecting - [y-pos] of affected
  ;print ydist
  let zdist [z-pos] of affecting - [z-pos] of affected
  let total-acc 0
  let xyangle 90
  let zangle 90
  let xydist (xdist ^ 2 + ydist ^ 2) ^ .5
  if abs(xdist) + abs(ydist) + abs(zdist) > 0 [
    set total-acc [mass] of affecting / max (list (xdist ^ 2 + ydist ^ 2 + zdist ^ 2) ((([size] of affected + [size] of affecting) / 2) ^ 2))
    let angle-heading atan xdist ydist
    let z-angle-heading atan xydist zdist
    set xyangle (90 - angle-heading) mod 360
    set zangle (90 - z-angle-heading) mod 360
  ]

  if total-acc > 50 [ set total-acc 50 ]

  set total-acc total-acc / 20
  let xyacc total-acc * (cos zangle)

  ask affected [
    set x-velocity (x-velocity + xyacc * (cos xyangle))
    set y-velocity (y-velocity + xyacc * (sin xyangle))
    set z-velocity (z-velocity + total-acc * (sin zangle))
  ]
end 

; Shifts positions of all turtles, used to stay centered around center of mass

to slide-view [deltax deltay deltaz]
  ask turtles [
    set x-pos x-pos - deltax
    set y-pos y-pos - deltay
    set z-pos z-pos - deltaz
  ]
  ;print "slide"
end 

; Zooms scale of the world in or out

to set-zoom [new-zoom]
  set zoom-factor zoom-factor * new-zoom
  print zoom-factor
  ask turtles [set size mass ^ (1 / 3) / zoom-factor * particle-scale]
end