A few more Rosetta Code entries to finish the year

[sarossell]

I hear ya. I've had plenty of crow myself (the trick is mustard).

I'm thinking it's a "one-way dead end' question. The kind that pops up to disprove mathematical theorems by proposing the ultimate reduction: 3^3, 2^2, 1^1, 0^0...Ummm. Now what? I just got involved because it was a simple Rosetta task and I'm learning BASIC with the easy ones first. It seemed odd, but eh, what the heck.

[Henko]

Explaining what a number to the 3th power means is easy, and that goes for all natural numbers 1,2,3,....ad inf.
The problem arizes when we ask for the meaning of a number to the power of 2.43 or 0.73.
Let us try a qualitative reasoning for the case X^2.43, X being non-zero.

X^2.43 means that X is multiplyed by itself more than 2 times, but less than 3 times.
Let us say then that this is equivalent to X multiplied 2 times with itself, and additional multiplied with a quantity less than itself, but 1 at minimum.
This means that when approaching the power of 1, say X^1.01, X is multiplied zero times with itself and addionally multiplied with a quantity very near to 1.
Now we pass the border and ask for the meaning of X^0.98. As we pass the border of doing something zero times, we can postulate that multiplying becomes division, and the quantity range 1 through X inverts also to X through 1
X^0.98 now means that X is divided zero times by itself, and addionally is divided by a quantity very near to 1. As the power lowers towards 0, the quantity by which X is divided grows near to X, and in the limit the division is X/X giving one.

Mathematically, the proof can be given by logarithms and applying a limit.

Speaking of cosmology, did you happen to find this program in the program section named "Newton"

Code: Select all

' 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=30 ! 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,scal)
  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 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)=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=17
  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=30
  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(),sc)
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))/sc,y_to_pix(yp(p))/sc,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()
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

[Joel]

As far as I can remember the best way to show the validity of e.g. e^0 is the development into taylor series around x=0 which is based on the characteristics of e-function (continuity, being its own derivative...)

[Henko]

This is indeed proof for e^x, having a special characteristic, hence that proof need not be valid for a^x, "a" being an arbitrary number. Try Taylor again for y=a^x and find yourself in a vicious circle problem .

[Joel]

Mhh, In that case I would extend a^x to e^(ln(a^x)) which leads to:
(e^x)^ln(a) where e^x is valid for all x resulting in b as shown before and b^ln(a) which is valid with ln(a) <> 0

[sarossell]

I like toast.

[Henko]

Mhh, In that case I would extend a^x to e^(ln(a^x)) which leads to:
(e^x)^ln(a) where e^x is valid for all x resulting in b as shown before and b^ln(a) which is valid with ln(a) <> 0

Nice trick! Clever

[GeorgeMcGinn]

Thanks, I'll try that next time - Goldens or the Grey Poop-on

I hear ya. I've had plenty of crow myself (the trick is mustard).

I'm thinking it's a "one-way dead end' question. The kind that pops up to disprove mathematical theorems by proposing the ultimate reduction: 3^3, 2^2, 1^1, 0^0...Ummm. Now what? I just got involved because it was a simple Rosetta task and I'm learning BASIC with the easy ones first. It seemed odd, but eh, what the heck.

[GeorgeMcGinn]

I like the program. It shows that eventually all the planets correct themselves and rotates around the main body counterclockwise.

I use a different method to calculating the raising of X by a decimal place. Either way is fine, I was just taught to convert into fractions.

For decimals, we convert them to a fraction and reduce. It makes doing the the math easier (at least for me, but when the numbers are more reasonable, it becomes easier). A decimal is nothing more than a fraction. I learned to do it this way first, and in a computer program it still works.

However, regardless of the method used, there is no way around the fact is that you will need to use some logic in guessing, as Henko stated, to determine where the answer falls.

To help understand how this works with fractions, whenever you raise a number to 1/2 power, that is the same as asking for the square root. So 144^1/2 = 12.

However, we have here an improper fraction: the decimal of 2.43 = 2 43/100 = 243/100 cannot be reduced further.

I was going to show an easier way to do this involving squaring fractions. However, due to the size of the numbers, i decided to use a scientific calculator as in both formulas below, doing this long hand will be nasty.

The rule is A^B/C = (A^B)^1/C also = (A^1/C)^B. The latter is usually the best way to solve this as the intermediate numbers remain small enough to work with.

Your problem is really 2^243/100, converting the decimal, which is an improper fraction, to a fraction.

I did this both ways and got the same answer, but I will do this the easy way – ( A^1/C)^B, or (2^1/100)^243

2 to the power of 1/100 = 1.00695555...
1.00695555^243=5.388934308

I had to use a calculator as even using the easier way, the numbers became too large to do long hand.

But if you use the formula to to 81^3/4, you get:
(81^1/4)^3.

First, 81^1/4 really means what is 81 cubed. So what X in 4^X = 81. The Answer is 3. 4 cubed = 81 (4x4=16x4=81)
So 81^1/4 = 3. Now we need to finish the formula, which is now 3^3, which is 27, the final answer.

So if you write a program to solve this as fractions instead of decimals, it will make a killer Rosetta code.

Yes, this is a little bit more than basic math, it's more or less algebra. But in my many years in computer science and writing programs in research, including cosmology, astronomy and physics, fractions are preferable over decimal places, and yes, they are harder to code (which is why I learned Assembler very early on).

But try this formula rule I gave you and try to write a basic program in it, and you should be able to write a function that can handle all situations.

May next time I will give you the rule for 2^-1/3 (yes, you can raise a number by a negative fraction, but not always).

Again, Henko, thanks for the program. I love examples showing Newton's law of motion and gravity.

George.


For example, 10^2.43 is the same as the improper fraction 2 43/100.

Explaining what a number to the 3th power means is easy, and that goes for all natural numbers 1,2,3,....ad inf.
The problem arizes when we ask for the meaning of a number to the power of 2.43 or 0.73.
Let us try a qualitative reasoning for the case X^2.43, X being non-zero.

X^2.43 means that X is multiplyed by itself more than 2 times, but less than 3 times.
Let us say then that this is equivalent to X multiplied 2 times with itself, and additional multiplied with a quantity less than itself, but 1 at minimum.
This means that when approaching the power of 1, say X^1.01, X is multiplied zero times with itself and addionally multiplied with a quantity very near to 1.
Now we pass the border and ask for the meaning of X^0.98. As we pass the border of doing something zero times, we can postulate that multiplying becomes division, and the quantity range 1 through X inverts also to X through 1
X^0.98 now means that X is divided zero times by itself, and addionally is divided by a quantity very near to 1. As the power lowers towards 0, the quantity by which X is divided grows near to X, and in the limit the division is X/X giving one.

Mathematically, the proof can be given by logarithms and applying a limit.

Speaking of cosmology, did you happen to find this program in the program section named "Newton"

Code: Select all

' 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=30 ! 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,scal)
  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 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)=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=17
  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=30
  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(),sc)
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))/sc,y_to_pix(yp(p))/sc,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()
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

[Joel]

again a wonderful small but highly efficient program worked out by Henko and a cool first approach to the n-body problem.
interessting: even with two objects, unfortunately the orbit will not remain stable in that simulation and the planet would spiral into the star (or move away from it?? can't remember how it was...) after some revolutions...would be fun if you could build in something like runge kutta

great control features btw... after a while i wished there was a lock-the-central-star-button though...

great work, thanks!!