Blackhole

turtles-own
[ fx    ;; x-component of force vector
  fy    ;; y-component of force vector
  vx    ;; x-component of velocity vector
  vy    ;; y-component of velocity vector
  xc    ;; real x coordinate (in case the particle leaves the screen)
  yc    ;; real y coordinate (in case the particle leaves the screen)
  mass  ;; the particle mass
]

breed [ stars star ]
breed [planets planet]

globals

[  center-of-mass-yc   ;; y-coordinate of the center of mass
   center-of-mass-xc   ;; x-coordinate of the center of mass
   g  ;; Gravitational constant
]

;; Setup Proceedures

to setup 
  clear-all
  set g 20
  set-default-shape turtles "circle"
  crt number
  ifelse symmetrical-setup?
  [ zero-sum-initial-setup ]
  [ random-initial-setup ]
  if keep-centered?
  [ recenter ]
  reset-ticks
end 

to random-initial-setup
  ask turtles
  [  set vx (( random-float (( 2 * max-initial-speed) - 1 )) - max-initial-speed )
     set vy (( random-float (( 2 * max-initial-speed) - 1 )) - max-initial-speed )
     set mass ( random-float max-initial-mass ) + 1
     set size sqrt mass 
     set heading ( random-float 360 )
     jump ( random-float ( max-pxcor - 10))
     set xc xcor
     set yc ycor    
     ]
end 

to zero-sum-initial-setup
    ;;First we set up the initial velocities of the first half of the particles.  
    ask turtles with [ who < ( number / 2 ) ]
    [  set vx ( random-float ((( 2 * max-initial-speed ) - 1 )) - max-initial-speed )
       set vy ( random-float ((( 2 * max-initial-speed ) - 1 )) - max-initial-speed )
       setxy random-xcor random-ycor
       set xc xcor
       set yc ycor  
       set mass ( random-float max-initial-mass ) + 1 
       set size sqrt mass
    ] 
  
  ;;Now, as we're zero-summing, we set the velocities of the second half of the 
  ;;particles to be the opposites of their counterparts in the first half. 
  
  ask turtles with  [ who >= (number / 2) ]
  [ set vx (- ([vx] of turtle (who - ( number / 2 ))))
    set vy (- ([vy] of turtle (who - ( number / 2 ))))
    set xc (- ([xc] of turtle (who - ( number / 2 ))))
    set yc (- ([yc] of turtle (who - ( number / 2 ))))
    setxy xc yc 
    set mass [mass] of turtle (who - (number / 2 ))
    set size sqrt mass 
  ]
  set center-of-mass-xc 0
  set center-of-mass-yc 0 
end  

;  This is the creation of the Sun in this model.  This would eventually become our black hole.

to create-particle 
  if mouse-down?
  [ let mx mouse-xcor
    let my mouse-ycor
    if (not any? turtles-on patch mx my)
    [
      crt 1
      [ set xc mx ;initial-psition-x 
        set xc yc ;initial-position-y
        set vx initial-velocity-x
        set vy initial-velocity-y
        set mass initial-mass
        set size sqrt mass
        set color particle-color
      ]
      display
    ]
  ]
    while [mouse-down?]
    [ ]
    if keep-centered?
    [
      recenter
      display
    ]
end 

to setup-two-planet
    set number 0 
    setup
    crt 1
    [ set color yellow
      set mass 200
      set size sqrt mass
      set xc -50
      set yc 0
      setxy xc yc
      set vx 0
      set vy -5
    ] 
    ;;crt 1
   ;; [ set color blue
    ;;  set mass 5 
     ;; set size sqrt mass
     ;; set xc 50
     ;; set yc 0
     ;; setxy xc yc 
     ;; set vx 0
     ;; set vy 9
    ;; ]
    crt 1
    [ set color Yellow
      set mass 200
      set size sqrt mass
      set xc 50
      set yc 0
      setxy xc yc
      set vx 0
      set vy 5
    ]
end 

;; Runtime Proceedures

to go 
  ask turtles
  [ set fx 0 
    set fy 0
  ]
  ;; must do all of these steps separately to get correct results
  ;; since all trutles interact with one another
  ask turtles [ check-for-collisions ] 
  ask turtles [ update-force ] 
  ask turtles [ update-velocity ] 
  ask turtles [ update-position ] 
  set-color
  if keep-centered?
  [ recenter ] 
  fade-patches 
  tick
end 

to check-for-collisions
  if any? other turtles-here
  [
    ask other turtles-here
    [
      set vx vx + [vx] of myself
      set vy vy + [vy] of myself
      set mass mass + [ mass ] of myself
      set size sqrt mass
    ]
    die
  ]
end   

to update-force ;; Turtle Procedure
  ;; This recrusive over all the turtles, each turtle asks this of all others turtles
  ask other turtles [ sum-its-force-on-me myself ]
end 

to sum-its-force-on-me [it] ;; Turtles Procedure
    let xd xc - [xc] of it
    let yd yc - [yc] of it 
    let d sqrt ((xd * xd) + (yd * yd)) 
    set fx fx + (cos (atan (- yd) (- xd))) * ([mass] of it * mass) / (d * d)
    set fy fy + (sin (atan (- yd) (- xd))) * ([mass] of it * mass) / (d * d)
end 

to update-velocity ;; Turtle Proceedure 
    ;; Now update each particle's velocity, by taking last time-step's velocity
    ;; and adding the effect of the force to it.  
    set vx (vx + (fx * g / mass))
    set vy (vy + (fy * g / mass))
end 

to update-position ;; Turtle Proceedure
    ;; As our system is closed, we can safely recenter of mass to the origin.
    set xc (xc + vx)
    set yc (yc + vy)
    adjust-position
end 

to adjust-position ;; Turtle Procedure
    ;; If we're in the visible world (the world inside the view)
    ;; update our x and y coordinates.
    ;; if there is no patch at the xc yc that means it is outside the world
    ;; and the turtale should just be hidden until it returns to the viewable world.
  setxy xc yc  
;    ifelse patch-at (xc - xcor) (yc - ycor ) != nobody
;    [ setxy xc yc
;      show-turtle
;      if (fade-rate != 100 )
;      [ set pcolor color + 3 ]
;    ]
;    [ hide-turtle ]
end 

to set-color
  ask turtle 1 [
    ifelse unhide  [
      set pcolor yellow
      ht
    ]
    [ st ]
  ]
end 
  
  ;; Center of Mass

to recenter
    find-center-of-mass
    ask turtles
    [ set xc (xc - center-of-mass-xc)
      set yc (yc - center-of-mass-yc)
      adjust-position
    ]
end 

to find-center-of-mass
    if any? turtles 
    [ set center-of-mass-xc sum [mass * xc] of turtles / sum [mass] of turtles
      set center-of-mass-yc sum [mass * yc] of turtles / sum [mass] of turtles
    ]
end 

to fade-patches
    ask patches with [pcolor != black]
    [ ifelse (fade-rate = 100)
      [ set pcolor black ]
      [ if (fade-rate != 0)
         [ fade ]
      ]
    ]
end 

to fade ;; Patch Proceedure
      let new-color pcolor - 8 * fade-rate / 100
      ;; if the new-color is no longer the same shade then it's faded to black.
      ifelse (shade-of? pcolor new-color)
      [ set pcolor new-color ]
      
      [ set pcolor black ]
end