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 9299 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$)