#**************************************
# Here, bodies do not start going in a random direction but rotate around the centre where a 'black hole' sits
# The black hole is the only body exerting gravity and devours stars that come to close. (no n-body problem) 
# As more stars get merged, the event horizon grows, gravitational effect of black hole increases
# Every time a star merges, a new one is created
# Biggest shortcoming: stars do not affect each other
#**************************************
#
#setup
gosub ScreenSetup
gosub ObjectDeclaration

#draw
while true
  gosub UpdateObjects
  gosub Show
end while

ScreenSetup:
  fastgraphics
  width=1000 
  height = 1000 
  graphsize (width, height)
  color black
  rect 0,0,width,height
  refresh
return


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

  #Set up n star bodies, each with 12 properties
  n = 1000
  dim Bodies[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 (x,y), distance to black hole (vector,x,y), velocity (x,y) and mass	
  for p = 0 to n-1
    Bodies[p,x_pos] = rand*width
    Bodies[p,y_pos] = rand*height
    Bodies[p, dx]   = px - Bodies[p,x_pos]
    Bodies[p, dy]   = py - Bodies[p,y_pos]
    Bodies[p, dist] = sqr(Bodies[p,dx]^2 + Bodies[p,dy]^2)
    if Bodies[p,dist] < 5 then Bodies[p,dist] = 5
    Bodies[p,speed] = 5 / sqr(Bodies[p,dist])
	Bodies[p,x_vel] = -Bodies[p, dy]/Bodies[p, dist] * Bodies[p,speed]
	Bodies[p,y_vel] =  Bodies[p, dx]/Bodies[p, dist] * Bodies[p,speed]
    Bodies[p,mass] = 1 + rand*4
  next p
return

UpdateObjects:
  for t = 0 to n-1
    #calculate the x and y components going from the black hole to the star
    Bodies[t,dx] = px - Bodies[t,x_pos]
    Bodies[t,dy] = py - Bodies[t,y_pos]
    #calculate distance of each star from the black hole (c² = x²  + y² so c = sqr(x² +y²)
    Bodies[t,dist] = sqr(Bodies[t,dx]^2 + Bodies[t,dy]^2)
    #if the distance gets too small the gravity will explode, so we cap the distance to 3
    if Bodies[t,dist] < 5 then Bodies[t,dist] = 5
    #gravity factor is mass over distance squared
    Bodies[t,gravity] = (Bodies[t,mass]*ev_hor)/((Bodies[t,dist]^2) + 25)
	#calculate the acceleration on x and y-axis based on the calculated gravity 
    Bodies[t,x_acc]= (Bodies[t,dx]/Bodies[t,dist]) * Bodies[t,gravity]
    Bodies[t,y_acc]= (Bodies[t,dy]/Bodies[t,dist]) * Bodies[t,gravity]
    #if the star has not merged with the black hole, update velocity and position
    if Bodies[t,crash]=0 then 
      Bodies[t,x_vel] = Bodies[t,x_vel]*0.999 + Bodies[t,x_acc] 
      Bodies[t,y_vel] = Bodies[t,y_vel]*0.999 + Bodies[t,y_acc]
      # add a tiny rotational bias to the velocity to creates differential rotation
      twist = 0.002
	  vx = Bodies[t,x_vel]
	  vy = Bodies[t,y_vel]
	  Bodies[t,x_vel] = vx - twist*vy
	  Bodies[t,y_vel] = vy + twist*vx
      Bodies[t,x_pos] = Bodies[t,x_pos] + Bodies[t,x_vel]
      Bodies[t,y_pos] = Bodies[t,y_pos] + Bodies[t,y_vel]
    end if
    #if the star enters the event horizon of the black hole and has not yet merged with it, 
    if Bodies[t,dist] < ev_hor and Bodies[t,crash]=0 then
		ev_hor = ev_hor + Bodies[t,mass] * 0.01
		if ev_hor > 10 then ev_hor = 10
	    Bodies[t,x_pos] = rand*width
	    Bodies[t,y_pos] = rand*height
	    Bodies[t, dx]   = px - Bodies[t,x_pos]
	    Bodies[t, dy]   = py - Bodies[t,y_pos]
	    Bodies[t, dist] = sqr(Bodies[t,dx]^2 + Bodies[t,dy]^2)
	    if Bodies[t,dist] < 5 then Bodies[t,dist] = 5
	    Bodies[t,speed] = 5 / sqr(Bodies[t,dist])
		Bodies[t,x_vel] = -Bodies[t, dy]/Bodies[t, dist] * Bodies[t,speed]
		Bodies[t,y_vel] =  Bodies[t, dx]/Bodies[t, dist] * Bodies[t,speed]
	    Bodies[t,mass] = 1 + rand*4
    endif
  next t  
return

Show:
  #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
    color rgb(Bodies[t,mass]*40, Bodies[t,mass]*30, Bodies[t,mass]*20)
    circle (Bodies[t,x_pos], Bodies[t,y_pos],Bodies[t,mass]/3 )
  next t    
  color white
  color rgb(ev_hor*12, 0, 0)
  circle px,py,ev_hor
  refresh
return