Celestial body simulation

Post Reply
Henko
Posts: 777
Joined: Tue Apr 09, 2013 12:23 pm
Flag: Netherlands

Celestial body simulation

Post by Henko » Wed Jul 25, 2018 6:19 pm

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
IMG_1538.PNG (127.76 KiB) Viewed 1300 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()
load_bodies()         ' starting system 
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
  read .mass(i),.diam(i), .xp(i),.yp(i),.vx(i),.vy(i)
  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 "load" text "Load" at 20,.xy+130 size 120,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
if bp("load") then
  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$)

Post Reply