Celestial bodies simulation (Newton's gravitational law)

[Henko]

' This program simulates the movement of a system of celestial bodies.
' The movements are governed by Newtons law: F = g x m1 x m2 / r2
' All bodies attract each other following that law.
' The system may have a dominant object like the sun in our system.
' If you want to manipulate the proces other than with the buttons,
' you need to know the following:
' maxb=the dimension of the arrays
' nb=the number of objects, read by the read-data construct
' in the initprog subroutine. This is the default system
'
' Use of the buttons at the bottom of the screen:
' Center: put the centre of gravity of the system in the screen centre
' Center always: center after each iteration
' Grid: toggle grid on the screen on and off
' Traces: toggle tracing of objects on and off
' Random solar: generate random system with a central star
' Random: genrate random system without a central star
' Black hole: put a black hole in the system (not yet satisfactory)
' Meteor: let a meteor traverse the system
' Zoom buttons: zoom in and out
'
option base 1 ! option angle degrees ! randomize
maxb=15 ! nb=5 ! gravity=.05

dim dist(maxb,maxb), newt(maxb,maxb), gonio(maxb,maxb)
dim f_x(maxb,maxb), f_y(maxb,maxb)
dim mass(maxb), diam(maxb), col(maxb,3), fx(maxb), fy(maxb)
dim vx(maxb), vy(maxb), xp(maxb), yp(maxb), xold(maxb), yold(maxb)

gosub initprog

loop1:
if cent=0 then goto loop3
loop2:
mtot=0 ! zx=0 ! zy=0 ' calc centre of gravity of system
for i=1 to nb
mtot=mtot+mass(i)
zx=zx+mass(i)*xp(i) ! zy=zy+mass(i)*yp(i)
next i
zx=zx/mtot ! zy=zy/mtot
for i=1 to nb
xp(i)=xp(i)-zx ! yp(i)=yp(i)-zy ' translate planets
next i
fill rect 0,0 to maxx,maxx ! grid_on(maxx,gr)

loop3:
graphics lock ' paint the system
for i=1 to nb
if abs(xp(i))<105*scal and yp(i)<105*scal then
bpaint(i,xold,yold,xp(i),yp(i),diam(i),col,traces,scal)
end if
next i
graphics unlock

for i=1 to nb-1 ' calculate dx and dy
for j=i+1 to nb
dist(i,j)=yp(j)-yp(i)
dist(j,i)=xp(j)-xp(i)
next j
next i
min=100 ! p=0 ! q=0
for i=1 to nb-1 ' calculate (squared) distances and forces
for j=i+1 to nb
newt(i,j)=dist(i,j)*dist(i,j) + dist(j,i)*dist(j,i)
newt(j,i)=gravity*mass(i)*mass(j)/newt(i,j) ' Newton formula
newt(i,j)=sqrt(newt(i,j))
if newt(i,j)<(diam(i)+diam(j))/6 then
p=i ! q=j
end if
if newt(i,j)<min then min=newt(i,j)
next j
next i
if p then
nb=merge(p,q,nb,mass,diam,col,xp,yp,xold,yold,vx,vy)
p=0 ! q=0 ! fill rect 0,0 to maxx,maxx
grid_on(maxx,gr) ! goto loop3
end if
for i=1 to nb-1 ' calculate angles (sin and cos)
for j=i+1 to nb
gonio(i,j)=dist(i,j)/newt(i,j)
gonio(j,i)=dist(j,i)/newt(i,j)
next j
next i
for i=1 to nb-1 ' calculate x-force components
for j=i+1 to nb
f_x(j,i)=gonio(j,i)*newt(j,i) ! f_x(i,j)=-f_x(j,i)
next j
next i
for i=1 to nb-1 ' calculate y-force components
for j=i+1 to nb
f_y(j,i)=gonio(i,j)*newt(j,i) ! f_y(i,j)=-f_y(j,i)
next j
next i
for j=1 to nb ' calculate total forces per planet
fx(j)=0 ! fy(j)=0
for i=1 to nb
fx(j)=fx(j)+f_x(i,j) ! fy(j)=fy(j)+f_y(i,j)
next i
next j
if min>5 then dt=1 else dt=0.1+.03*min*min
if dt=1 then
for i=1 to nb ' calculate accel, velocity and position
acc=fx(i)/mass(i) ! vx(i)=vx(i)+acc ! xp(i)=xp(i)+vx(i)-acc/2
acc=fy(i)/mass(i) ! vy(i)=vy(i)+acc ! yp(i)=yp(i)+vy(i)-acc/2
next i
else
for i=1 to nb ' calculate accel, velocity and position
acc=fx(i)/mass(i) ! vx(i)=vx(i)+acc*dt
xp(i)=xp(i)+vx(i)*dt-acc*dt*dt/2
acc=fy(i)/mass(i) ! vy(i)=vy(i)+acc*dt
yp(i)=yp(i)+vy(i)*dt-acc*dt*dt/2
next i
end if

if button_pressed("grid") then
gr=1-gr ! grid_on(maxx,gr)
end if
if button_pressed("center") then goto loop2
if button_pressed("centeralw") then cent=1-cent
if button_pressed("trace") then traces=1-traces
if button_pressed("zoomplus") then
fill rect 0,0 to maxx,maxx ! scal=scal/2 ! grid_on(maxx,gr)
end if
if button_pressed("zoommin") then
fill rect 0,0 to maxx,maxx ! scal=scal*2 ! grid_on(maxx,gr)
end if
if button_pressed("solar") then
mass(1)=4900 ! diam(1)=10 ! gravity=0.01 ! nb=12
xp(1)=0 ! yp(1)=0 ! vx(1)=0 ! vy(1)=0 ! xold(1)=0 ! yold(1)=0
for i=2 to nb
diam(i)=2+rnd(7) ! mass(i)=diam(i)*diam(i)
col(i,1)=.5+rnd(.5) ! col(i,2)=.5+rnd(.5) ! col(i,3)=.5+rnd(.5)
rr=10*i-8-rnd(5) ! ang=60*i+rnd(5)
xp(i)=rr*cos(ang) ! yp(i)=rr*sin(ang)
xold(i)=xp(i) ! yold(i)=yp(i)
vel=7/sqrt(rr) ! dd=1-2*rnd(2)
vx(i)=vel*sin(ang)*dd ! vy(i)=-vel*cos(ang)*dd
next i
fill rect 0,0 to maxx,maxx ! grid_on(maxx,gr)
end if
if button_pressed("rand") then
gravity=0.01 ! nb=8
for i=1 to nb
diam(i)=2+rnd(7) ! mass(i)=diam(i)*diam(i)
col(i,1)=.5+rnd(.5) ! col(i,2)=.5+rnd(.5) ! col(i,3)=.5+rnd(.5)
xp(i)=90-rnd(180) ! yp(i)=90-rnd(180)
xold(i)=xp(i) ! yold(i)=yp(i)
vx(i)=.2-rnd(.3) ! vy(i)=.2-rnd(.3)
vx(i)=0 ! vy(i)=0
next i
fill rect 0,0 to maxx,maxx ! grid_on(maxx,gr)
end if
if button_pressed("meteor") then
if meteo=0 then
meteo=1 ! nb=nb+1
end if
mass(nb)=10 ! diam(nb)=3 ! col(nb,1)=1 ! col(nb,2)=1 ! col(nb,3)=1
xp(nb)=-110 ! yp(nb)=-rnd(100) ! xold(nb)=xp(nb)! yold(nb)=yp(nb)
vx(nb)=2 ! vy(nb)=1
end if
if button_pressed("hole") then
if meteo=0 then
meteo=1 ! nb=nb+1
end if
mass(nb)=10000 ! diam(nb)=2 ! col(nb,1)=0 ! col(nb,2)=0 ! col(nb,3)=0
xp(nb)=50-rnd(100) ! yp(nb)=50-rnd(100) ! vx(nb)=0 ! vy(nb)=0
end if

goto loop1
end

def bpaint(i,xold(),yold(),xn,yn,dia,col(,),tr,sc)
xpix=x_to_pix(xold(i)/sc) ! ypix=y_to_pix(yold(i)/sc)
if ypix<768-dia then
fill rect xpix,ypix size dia+tr
end if
xpix=x_to_pix(xn/sc) ! ypix=y_to_pix(yn/sc)
if ypix<768-dia then
fill color col(i,1),col(i,2),col(i,3)
fill circle xpix,ypix size dia
fillback()
end if
xold(i)=xn ! yold(i)=yn
end def

def merge(p,q,nb,mass(),diam(),col(,),xp(),yp(),xold(),yold(),vx(),vy())
mtot=mass(p)+mass(q)
r=diam(p) ! if diam(q)<r then r=diam(q)
diam(p)=sqrt(diam(p)^2+diam(q)^2)
if diam(p)>10 then diam(p)=10
for j=1 to 3
col(p,j)=(mass(p)*col(p,j)+mass(q)*col(q,j))/mtot
next j
vx(p)=(mass(p)*vx(p)+mass(q)*vx(q))/mtot
vy(p)=(mass(p)*vy(p)+mass(q)*vy(q))/mtot
if q<nb then
for k=q to nb-1
mass(k)=mass(k+1) ! diam(k)=diam(k+1)
for j=1 to 3 ! col(k,j)=col(k+1,j) ! next j
xp(k)=xp(k+1) ! yp(k)=yp(k+1)
xold(k)=xold(k+1) ! yold(k)=yold(k+1)
vx(k)=vx(k+1) ! vy(k)=vy(k+1)
next k
end if
nb=nb-1
blast(x_to_pix(xp(p)),y_to_pix(yp(p)),10*r)
merge=nb
end def

def x_to_pix(x) = 3.84*x+384
def y_to_pix(y) = 384-3.84*y

def fillback()
fill color .2,.2,.2
end def

def grid_on(mx,on)
draw size 1
if on then draw color .3,.3,.3 else draw color .2,.2,.2
for aa=96 to 672 step 96
draw line 0,aa to mx,aa ! draw line aa,0 to aa,mx
next aa
end def

def blast(x,y,rad)
for i=1 to 500
fill color 1,1,1
fill circle x,y size i*rad/500
next i
for i=1 to 200
fill color 1,1,1-i/250
fill circle x,y size rad
next i
for i=1 to 200
fill color 1,1-i/250,.2
fill circle x,y size rad
next i
for i=1 to 200
fill color 1-i/250,.2,.2
fill circle x,y size rad+1
next i
end def

initprog:
graphics ! graphics clear .2,.2,.2
maxx=screen_width() ! maxy=screen_height() ! set orientation top
fill color .8,.8,.8 ! fill rect 0,769 to maxx,maxy ! fillback()
xc=maxx/2 ! yc=xc ! traces=1 ! gr=0 ! cent=0 ! scal=1 ! meteo=0
randomize
for i=1 to nb
read mass(i),diam(i)
for j=1 to 3 ! read col(i,j) ! next j
read xp(i),yp(i),vx(i),vy(i)
xold(i)=xp(i) ! yold(i)=yp(i)
next i
data 100,8,1,1,0.8,0,0,0,0, 10,4,0.6,0.7,1,50,0,0,.33
data 7,3,1,0,0,-30,0,0,.4, 0.5,2,1,1,.6,45,0,0,.1
data 12,5,0,0,1,0,60,.4,0
button "center" title "Center" at 20,maxx+10 size 120,50
button "centeralw" title "Center always" at 160,maxx+10 size 120,50
button "grid" title "Grid" at 300,maxx+10 size 120,50
button "trace" title "traces" at 440,maxx+10 size 120,50
button "solar" title "random Solar" at 580,maxx+10 size 120,50
button "rand" title "random" at 580,maxx+80 size 120,50
button "zoomplus" title "zoom +" at 20,maxx+80 size 120,50
button "zoommin" title "zoom -" at 160,maxx+80 size 120,50
button "meteor" title "Meteor" at 300,maxx+80 size 120,50
button "hole" title "Black hole" at 440,maxx+80 size 120,50
return

[Mr. Kibernetik]

Very interesting!

I suggest you to force device rotation to vertical, because otherwise buttons can be missing from view...

[Henko]

Very interesting!

I suggest you to force device rotation to vertical, because otherwise buttons can be missing from view...

Indeed. I added "set orientation top" to the initprog subroutine

Very nice Basic now.

[Mr. Kibernetik]

Very nice Basic now.

Does it run normally on your device?

[Mr. Kibernetik]

It is more correct to put SET ORIENTATION before getting screen dimensions, because if device was initially horizontal then screen dimensions will not suit after change of orientation.

[Henko]

Very nice Basic now.

Does it run normally on your device?

I didn't use all new features yet. But so far so good (iPad3 with retina screen)

[Mr. Kibernetik]

I didn't use all new features yet. But so far so good (iPad3 with retina screen)

Very nice to know that you like it!