#**************************************
# Single black hole devouring stars that come to close.
# As more stars get merged, the event horizon grows, gravitational effect of black hole increases and stars stop orbiting
# Biggest shortcoming: stars do not affect each other
#**************************************
#screen setup
#*************
fastgraphics
width=1000
height = 1000
graphsize (width, height)
color black
rect 0,0,width,height
Outputvisible 0
Editvisible 0
refresh

#place black hole initialisation
#*******************************************
ev_hor = 5.1
px=width/2 : py=height/2

#Set up n star systems with 12 properties
#******************************************
n = 2000
dim sol[n,13]
x_pos = 0
y_pos = 1
x_vel = 2
y_vel = 3
dist = 4
dx = 5
dy = 6
x_acc = 7
y_acc = 8
gravity = 9
crash = 10
mass = 11
speed = 12
# initialize position and speed/direction   
for p = 0 to n-1
    sol[p,x_pos] = rand*width
    sol[p,y_pos] = rand*height
    sol[p, dx]   = px - sol[p,x_pos]
    sol[p, dy]   = py - sol[p,y_pos]
    sol[p, dist] = sqr(sol[p,dx]^2 + sol[p,dy]^2)
    if sol[p,dist] < 5 then sol[p,dist] = 5
    sol[p,speed] = 5 / sqr(sol[p,dist])
    sol[p,x_vel] = -sol[p, dy]/sol[p, dist] * sol[p,speed]
    sol[p,y_vel] =  sol[p, dx]/sol[p, dist] * sol[p,speed]
    sol[p,mass] = 1 + rand*4
next p


while true
  #paint the canvas black with a small transparency to see the trails
  color rgb(0,0,0,15)
  rect 0,0,width,height
  for t = 0 to n-1
    #calculate the x and y components going from the black hole to the star
    sol[t,dx] = px - sol[t,x_pos]
    sol[t,dy] = py - sol[t,y_pos]
    #calculate distance of each star from the black hole (c² = x²  + y² so c = sqr(x² +y²)
    sol[t,dist] = sqr(sol[t,dx]^2 + sol[t,dy]^2)
    #if the distance gets too small the gravity will explode, so we cap the distance to 3
    if sol[t,dist] < 5 then sol[t,dist] = 5
    #gravity factor is mass over distance squared
    sol[t,gravity] = (sol[t,mass]*ev_hor)/((sol[t,dist]^2) + 25)
    #calculate the acceleration on x and y-axis based on the calculated gravity
    sol[t,x_acc]= (sol[t,dx]/sol[t,dist]) * sol[t,gravity]
    sol[t,y_acc]= (sol[t,dy]/sol[t,dist]) * sol[t,gravity]
    #if the star has not merged with the black hole, update velocity and position
    if sol[t,crash]=0 then
      sol[t,x_vel] = sol[t,x_vel]*0.999 + sol[t,x_acc]
      sol[t,y_vel] = sol[t,y_vel]*0.999 + sol[t,y_acc]
      twist = 0.002
      vx = sol[t,x_vel]
      vy = sol[t,y_vel]
      sol[t,x_vel] = vx - twist*vy
      sol[t,y_vel] = vy + twist*vx
      sol[t,x_pos] = sol[t,x_pos] + sol[t,x_vel]
      sol[t,y_pos] = sol[t,y_pos] + sol[t,y_vel]
    end if
    #if the star enters the event horizon of the black hole and has not yet merged with it,
    if sol[t,dist] < ev_hor and sol[t,crash]=0 then
       ev_hor = ev_hor + sol[t,mass] * 0.01
       if ev_hor > 10 then ev_hor = 10
       sol[t,x_pos] = rand*width
       sol[t,y_pos] = rand*height
       sol[t, dx]   = px - sol[t,x_pos]
       sol[t, dy]   = py - sol[t,y_pos]
       sol[t, dist] = sqr(sol[t,dx]^2 + sol[t,dy]^2)
       if sol[t,dist] < 5 then sol[t,dist] = 5
       sol[t,speed] = 5 / sqr(sol[t,dist])
       sol[t,x_vel] = -sol[t, dy]/sol[t, dist] * sol[t,speed]
       sol[t,y_vel] =  sol[t, dx]/sol[t, dist] * sol[t,speed]
       sol[t,mass] = 1 + rand*4
    endif
    color rgb(sol[t,mass]*40, sol[t,mass]*30, sol[t,mass]*20)
    circle (sol[t,x_pos], sol[t,y_pos],sol[t,mass]/3 )
  next t   
  color white
  color rgb(ev_hor*12, 0, 0)
  circle px,py,ev_hor
  refresh
end while