Orbital Playground

breed [particles particle] ; turtles that represent planets, celestial bodies
breed [plotpoints plotpoint] ; turtles used to draw orbital paths
particles-own [ x-velocity y-velocity z-velocity ; used to store the 3-dimensional velocity of the particle
   x-pos y-pos z-pos ; used to store the 3d position (separate from actual xcor/ycor/zcor to allow the world to zoom and change perspective
   mass ; represents physical mass. Higher mass means greater gravitational pull
   density ; turtle size = mass / density (I wanted particles with greater mass to be larger, but didn't want something extremely massive to completely dominate the screen)
   identifier ; gives a simple way for the user to select specific planets
   copy? ; if true, the particle is a hidden copy used only for calculating
   ]
plotpoints-own [ number ; when an orbital path is plotted, the points are numbered for easy link creation
   x-pos y-pos z-pos ; same as particles - these are separate from actual x/y/zcor for zooming purposes
   ]
globals [ zoom-factor ; used to effectively zoom the world in and out; just scales positions and sizes accordingly
   max-plotpoint ; stores the next plotpoint number
   ]

to setup
  clear-all
  set-default-shape turtles "circle"
  set zoom-factor 1
  create-particles 1 [
    setxyz 0 0 0
    set x-velocity 0
    set y-velocity 0
    set z-velocity 0
    set mass 1000
    set density 1
    set size mass ^ (1 / 3) / density
    set x-pos xcor
    set y-pos ycor
    set z-pos zcor
    set color yellow
    set identifier "sun"
    set copy? false
    ]
  create-particles 1 [
    setxyz -20 0 0
    set x-velocity 0
    set y-velocity 1
    set z-velocity 0
    set mass 10
    set density 1
    set size mass ^ (1 / 3) / density
    set x-pos -40
    set y-pos ycor
    set z-pos zcor
    set identifier 1
    set copy? false
    set color red
    ]
  create-particles 1 [
    ;setxyz 45 0 0
    set x-velocity 0
    set y-velocity -.7
    set z-velocity 0
    set mass 40
    set density 1
    set size mass ^ (1 / 3) / density
    set x-pos 90
    set y-pos 0
    set z-pos 0
    set identifier 2
    set copy? false
    set color green
  ]
  create-particles 1 [
    ;setxyz 0 40 0
    set x-velocity .75
    set y-velocity 0
    set z-velocity 0
    set mass 10
    set density 1
    set size mass ^ (1 / 3) / density
    set x-pos 0
    set y-pos 80
    set z-pos 0
    set identifier 3
    set copy? false
    set color blue
  ]
  set-zoom 3
  set-visible-positions particles
  reset-ticks
end 

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

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

; checks center of mass of entire system, shifts positions accordingly to keep view centered

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

; calculates 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 position of turtles relative to boundaries of world, zooms in/out if needed

to check-boundaries
  let flag false
  ask particles [
    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 particles [
    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 

; top-level function to compute changes in acceleration due to gravitational force

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

; move all turtles based on their velocity

to move [turtleset]
  ask turtleset [
    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 

to set-visible-positions [turtleset]
  ask turtleset [
    if not hidden? [
      set xcor x-pos / zoom-factor
      set ycor y-pos / zoom-factor
      set zcor z-pos / zoom-factor
    ]
  ]
end 

; given two particles, calculate the gravitational pull of one upon the other

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 particles and plotpoints (effectively sliding "camera")

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
  ]
  set-visible-positions plotpoints
end 

; changes zoom factor, effectively zooming "camera" in or out

to set-zoom [new-zoom]
  set zoom-factor zoom-factor * new-zoom
  print zoom-factor
  ask particles [set size mass ^ (1 / 3) / density / zoom-factor]
  set-visible-positions plotpoints
end 

;; reports whichever turtle is selected in the "selected-body" dropdown menu

to-report selected-turtle
  let selected-set particles with [identifier = selected-body and not copy?]
  ;ifelse length selected-set = 1
  ;[
    report one-of selected-set
  ;]
  ;[ ; something is wrong
  ;  print "error selecting set"
  ;  report 0
  ;]
end 

;; Plots the projected path of the currently selected body
;; It does this by spawning hidden copies of every turtle, running physics on only the copies,
;; and hatching plotpoints every "tick"

to plot-orbital-path
  clear-orbital-path
  let selected selected-turtle
  if selected != 0 [
    ; spawn copies of all particles
    ask particles [
      hatch 1 [ set copy? true
        set hidden? true
      ]
    ]
    set max-plotpoint 0
    let temp path-plot-length
    while [temp > 0] [
      tick-copies
      set temp temp - 1
    ]

    ;delete copies (this is important)
    ask particles with [copy?] [die]
  ]
end 

;; Clears currently plotted path

to clear-orbital-path
  ask plotpoints [die]
end 

;; runs through physics only with copies

to tick-copies
  let copies particles with [copy?]
  compute-acceleration-changes copies
  move copies
  set-orbital-path-point
end 

;; puts down a new point on the projected orbital path

to set-orbital-path-point
  ask one-of particles with [copy? and identifier = selected-body]
  [
    hatch-plotpoints 1
    [
      set hidden? false
      set x-pos [x-pos] of myself
      set xcor x-pos / zoom-factor
      set y-pos [y-pos] of myself
      set ycor y-pos / zoom-factor
      set z-pos [z-pos] of myself
      set zcor z-pos / zoom-factor
      set size .1 / zoom-factor
      set number max-plotpoint
      if number > 0 [
        let prev-plotpoint one-of plotpoints with [number = max-plotpoint - 1]
        create-link-with prev-plotpoint
      ]
    ]
  ]
  set max-plotpoint max-plotpoint + 1
end 

;; A function to modify the velocity of the selected particle based on the corresponding slider
;; values in the interface.

to modify-velocity
  ask selected-turtle [
    set x-velocity x-velocity + x-velocity-mod
    set y-velocity y-velocity + y-velocity-mod
    set z-velocity z-velocity + z-velocity-mod
  ]
end