## Celestial body simulation

Henko
Posts: 797
Joined: Tue Apr 09, 2013 12:23 pm
Windows
Location: Groningen, Netherlands
Flag:

### Celestial body simulation

A remake of a program some years ago. This time with the use of sprites and some other simplification. The dynamics are based upon Newtons well known gravitation law. In this sim the units are focussed on the usage of the screen. Moreover, the simulation is 2-dimensional, in x and y.
A 3-D simulation would not be more complex, but requires more calculations, and eventually, the 3-D results have to be projected on the 2-D screen after all.
You may enter your own stellar system data in the red code sections. For each object that would be:
- mass
- diameter
- x and y of starting positio
- vx and vy of starting velocity vector
Look for the example for about proper magnitudes.
The calculation method is such that every object influences and is influenced by all other objects, as should be the case.
IMG_1538.PNG (127.76 KiB) Viewed 7649 times

Code: Select all

``````'r'
maxb=30               ' maximum # of celestial bodies
gravity=.09           ' gravitational constant (in this sim)
col_sense=.2          ' sensitivity for collision ( 0>= x <=1)
deffn\$="celest_test"  ' default data file
''
option base 1
dim mass(maxb),diam(maxb),vx(maxb),vy(maxb),xp(maxb),yp(maxb)
init_prog()
make_sprites()
do
if rest=0 then gravity_calc(mass,diam,vx,vy,xp,yp,dt)
show_sprites()
user_action()
until forever
end

def load_bodies()  '  fill / change .nb and the data statements
'r'
.nb=3
''
mm=0
for i=1 to .nb
mm=max(mm,.mass(i)) ! .max_m=mm
next i
'r'
data 100,8,-100,-100,0,0,   100,8,-100,200,.1,0
data 1,3,300,-100,-1,-.0058
''
end def

def init_prog()
graphics ! graphics clear .2,.2,.2
set orientation portrait ! set toolbar off
option sprite pos central
randomize
maxx=screen_width() ! maxy=screen_height() ! .xy=min(maxx,maxy)
fill color .7,.7,.7 ! fill rect 0,769 to .xy,maxy
.sc=.xy/2
.trail=1 ! .centron=0 ! .scale=1 ! .rest=1 ! .dt=0.05
set buttons custom
draw color 1,0,0 ! fill color 0,1,0 ! fill alpha .3
set buttons font size 30
switch "trail" state 1 at maxx-220,.xy+10
switch "collo" state 0 at maxx-220,.xy+60
switch "centr" state 0 at maxx-220,.xy+110
button "zoom+" text "+" at 20,.xy+10 size 60,40
button "zoom-" text "-" at 230,.xy+10 size 60,40
button "speed+" text "+" at 230,.xy+70 size 60,40
button "speed-" text "-" at 20,.xy+70 size 60,40
button "save" text "Save" at 20,.xy+190 size 120,40
button "stop" text "Quit" at maxx-180,maxy-60 size 150,40
button "cent" text "Centre" at 170,.xy+130 size 120,40
button "rand" text "Random" at 170,.xy+190 size 120,40
button "start" text "Start" at 355,.xy+70 size 120,40
button "clear" text "Clear" at 355,.xy+190 size 120,40
field "fname" text .deffn\$ at 310,.xy-130 size 200,60
field "fname" hide
draw color 0,0,0
draw text "<- Zoom ->" at 95,.xy+20
draw text "<- Speed ->" at 90,.xy+80
draw text "Trail on" at maxx-160,.xy+15
draw text "Collisions on" at maxx-160,.xy+65
draw text "Centre always" at maxx-160,.xy+115
fill alpha 1
sprite "nova" begin 56,56
fill color 1,1,.8 ! fill circle 28,28 size 25
sprite end
sprite "nova" hide
end def

def make_sprites()
for i=1 to .nb
d=.diam(i) ! d2=2.2*d
sprite i begin d2,d2
set_fill_color(.mass(i))
fill circle d2/2,d2/2 size d
sprite end
next i
end def

def show_sprites()
for i= 1 to .nb
xs=.sc+.scale*.xp(i) ! ys=.sc-.scale*.yp(i) ! d=.diam(i)
sprite i at xs,ys
if ys<.xy-2*d then
sprite i show
if .trail then draw pixel 2*(xs),2*(ys) color .4,.4,.4
else ! sprite i hide
end if
next i
end def

def user_action()
if bp("start") then
.rest=1-.rest
if .rest then ! button "start" text "Start"
else ! button "start" text "Pause"
end if
end if
if bp("zoom+") then .scale*=1.5
if bp("zoom-") then .scale/=1.5
if bp("speed+") then .dt*=1.5
if bp("speed-") then .dt/=1.5
if bp("cent") then centre()
if bp("rand") then
.nb=rnd(15) ! .gravity=.1+rnd(.1)
for i=1 to .nb
.mass(i)=4+rnd(10) ! .diam(i)=2+rnd(7)
.xp(i)=200-rnd(400) ! .yp(i)=200-rnd(400)
.vx(i)=0 ! .vy(i)=0
next i
make_sprites()
end if
field "fname" show ! field "fname" select
do until field_changed("fname")
deffn\$=field_text\$("fname")
if not file_exists(deffn\$) then
field "fname" hide ! return
end if
file deffn\$ reset
file deffn\$ input .nb, .gravity, .dt
for i=1 to .nb
file deffn\$ input .mass(i),.diam(i),.xp(i),.yp(i),.vx(i),.vy(i)
next i
field "fname" hide
make_sprites()
end if
if bp("save") then
if .nb=0 then return
field "fname" show ! field "fname" select
do until field_changed("fname")
deffn\$=field_text\$("fname")
if file_exists(deffn\$) then file deffn\$ delete
file deffn\$ print .nb, .gravity, .dt
for i=1 to .nb
file deffn\$ print .mass(i),.diam(i),.xp(i),.yp(i),.vx(i),.vy(i)
next i
field "fname" hide
end if
if switch_state("collo") then
for i=1 to .nb-1 ! i\$=i
for j=i+1 to .nb ! j\$=j
if sprites_collide(i\$,j\$) then
s=sqrt((.xp(i)-.xp(j))^2+(.yp(i)-.yp(j))^2)
if s<.col_sense*(.diam(i)+.diam(j)) then merge(i,j)
end if
next j
next i
end if
.centron=switch_state("centr")
.trail=switch_state("trail")
if bp("clear") then
fill color .2,.2,.2
fill rect 0,0 to .xy,769

end if
if bp("stop") then stop
end def

def centre()
mtot=0 ! zx=0 ! zy=0
for i=1 to .nb
mtot+=.mass(i) ! zx+=.mass(i)*.xp(i) ! zy+=.mass(i)*.yp(i)
next i
zx=zx/mtot ! zy=zy/mtot
for i=1 to .nb ! .xp(i)-=zx ! .yp(i)-=zy ! next i
end def

def merge (p,q)    ' collision !
super_nova(.xp(p),.yp(p))
mn=.mass(p)+.mass(q) ! .mass(p)=mn
.diam(p)=(.diam(p)^3+.diam(q)^3)^(1/3)
.xp(p)=(.xp(p)+.xp(q))/2 ! .yp(p)=(.yp(p)+.yp(q))/2
.vx(p)=(.mass(p)*.vx(p)+.mass(q)*.vx(q))/mn
.vy(p)=(.mass(p)*.vy(p)+.mass(q)*.vy(q))/mn
sprite p delete ! sprite q delete
d=.diam(p) ! d2=2.2*d
sprite p begin d2,d2
set_fill_color(.mass(p))
fill circle d2/2,d2/2 size d
sprite end
.nb-=1
for i=q to .nb
.mass(i)=.mass(i+1) ! .diam(i)=.diam(i+1)
.xp(i)=.xp(i+1) ! .yp(i)=.yp(i+1)
.vx(i)=.vx(i+1) ! .vy(i)=.vy(i+1)
sprite str\$(i+1) rename str\$(i)
next i
end def

def super_nova(x,y)
xs=.sc+.scale*x ! ys=.sc-.scale*y
sprite "nova" at xs,ys ! sprite "nova" show
for s=.5 to 2.5 step .1
sprite "nova" at xs,ys scale s ! pause .01
next s
for a=1 to .1 step -.1
sprite "nova" alpha a ! pause .1
next a
sprite "nova" hide
end def

def set_fill_color(m)
if m=.max_m then ! fill color .9,.9,.9
else
t=max(1,100*rnd(4)+rnd(40)-20) ! pal(t)
fill color pal.r, pal.g, pal.b
end if
end def

def gravity_calc(mass(),diam(),vx(),vy(),xp(),yp(),dt)
n=.maxb ! nb=.nb
dim dist(n,n),newt(n,n),gonio(n,n),f_x(n,n),f_y(n,n)
dim fx(n),fy(n)
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
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))
next j
next i
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 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)
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)+=f_x(i,j) ! fy(j)+=f_y(i,j) ! next i
next j
if dt=1 then
for i=1 to nb         ' calculate accel, velocity and position
acc=fx(i)/mass(i) ! vx(i)+=acc ! xp(i)+=vx(i)-acc/2
acc=fy(i)/mass(i) ! vy(i)+=acc ! 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)+=acc*dt
xp(i)+=vx(i)*dt-acc*dt*dt/2
acc=fy(i)/mass(i) ! vy(i)+=acc*dt
yp(i)+=vy(i)*dt-acc*dt*dt/2
next i
end if
if .centron then centre()
end def

def pal(t)
r=0 ! g=0 ! b=0
if t<120 or t>240 then r=palsub(abs(t-360*floor(t/240)))
if t<240 then g=palsub(abs(t-120))
if t>120 then b=palsub(abs(t-240))
end def
def palsub(e)
f=.5   ' 0<=f<=1 better balance between prim. and sec. colors
if e<60 then c=1 else ! x=(120-e)/60 ! c=x*(1+f-f*x) ! end if
return c
end def

'  atan for option degrees
'
def deg_atan(x,y)
if x=0 then return 90*(2-sign(y)) else ta=atan(y/x)
if x<0 then ! ta+=180 ! else ! if y<0 then ! ta+=360 ! end if ! end if
return ta
end def

def bp(a\$) = button_pressed(a\$)
``````