Gravity Wells

globals [ zoom-factor ] ;; used to "zoom" by scaling the displayed positions of particles
breed [ stayers stayer ] ;; used as points to form the mesh that displays the gravity wells
breed [ particles particle ] ;; particles that interact and exert gravitational forces upon each other
turtles-own [
  x-pos y-pos z-pos ;; used as the basis for physics calculations - xcor, ycor, zcor only used for display
]
stayers-own [
  gravity ;; the magnitude of the total gravitational field at a given point
]
particles-own [
  x-velocity y-velocity ;; components of 2d velocity, used for physics calculations
  mass ;; also used for physics; the greater the mass, the greater the gravitational pull
]

to setup
  clear-all

  ;; create stayers
  ask patches with [pzcor = 0] [
    sprout-stayers 1
    [
      set gravity 0
      create-links-with turtles-on neighbors
      set x-pos xcor
      set y-pos ycor
      hide-turtle
    ]
  ]

  set-default-shape turtles "circle"
  set zoom-factor 1
  create-particles number-of-objects [
    setxyz (random-float starting-box - random-float starting-box) (random-float starting-box - random-float starting-box) 0
    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
    ]
  compute-gravity-wells
  set-particle-display-positions
  reset-ticks
end 

to go
  compute-acceleration-changes
  check-boundaries
  move-particles
  check-slide
  compute-gravity-wells
  set-particle-display-positions
  tick
end 

to tick-once
  compute-acceleration-changes
  check-boundaries
  move-particles
  check-slide
  compute-gravity-wells
  set-particle-display-positions
  tick
end 

;; Update positions of stayers to correspond with gravitational field

to compute-gravity-wells

  ask stayers [ set gravity 0
    set x-pos xcor * zoom-factor
    set y-pos ycor * zoom-factor
  ]
  ask particles [
    ask stayers with [distancexy [xcor] of myself [ycor] of myself < SoI] [
      let xdist [x-pos] of self - [x-pos] of myself
      let ydist [y-pos] of self - [y-pos] of myself
      let scale-const max list ((xdist ^ 2 + ydist ^ 2) ^ .5) ([size] of myself / 2)

      set gravity gravity + [mass] of myself * G / scale-const
    ]
  ]
  ask stayers [
    set z-pos gravity * -1
    set zcor z-pos / zoom-factor
    ; set color scale-color red zcor 0 -50
  ]
end 

to set-particle-display-positions
  ask particles [
    set z-pos [z-pos] of min-one-of stayers [distancexy [xcor] of myself [ycor] of myself]
    set zcor z-pos / zoom-factor

    if not hidden? [
      set xcor x-pos / zoom-factor
      set ycor y-pos / zoom-factor
      set size mass ^ (1 / 3) / zoom-factor
    ]
  ]
end 

;; Check center of mass of the overall system, adjust view accordingly

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

;; Calculate the center of mass of the entire system

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))
end 

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
    ]
  ]

  ifelse flag [
    ;print "zooming out"
    set-zoom 1.1
  ]
  [
    ;print "zooming in"
    set-zoom .9
  ]

  ask particles [
    ifelse abs (x-pos + x-velocity) > max-pxcor or abs (y-pos + y-velocity) > max-pycor or abs (z-pos) > max-pzcor
    [ set hidden? true ]
    [ set hidden? false ]
  ]
end 

;; Have particles compute their acceleration based on gravitational interactions with all other particles

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

;; Move particles according to their velocity

to move-particles
  ask particles [
    set x-pos x-pos + x-velocity
    set y-pos y-pos + y-velocity
    ;set z-pos z-pos + z-velocity
  ]
end 

;; Compute gravitational pull between a pair of particles

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 xdist + ydist > 0 [
    set total-acc [mass] of affecting / (xdist ^ 2 + ydist ^ 2)
    let angle-heading atan xdist ydist
    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 + total-acc * (cos xyangle))
    set y-velocity (y-velocity + total-acc * (sin xyangle))
    ;set z-velocity (z-velocity + total-acc * (sin zangle))
  ]
end 

;; Move view laterally

to slide-view [deltax deltay deltaz]
  ask particles [
    set x-pos x-pos - deltax
    set y-pos y-pos - deltay
    set z-pos z-pos - deltaz
  ]
end 

;; Change zoom-factor to zoom in or out

to set-zoom [new-zoom]
  set zoom-factor zoom-factor * new-zoom
end