A few more Rosetta Code entries to finish the year

[Henko]

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!!

There is : the "center always" button. Be it, that it does not center the central star, but the mass center of the entire system. In that case, systems without a "central star" can also be more or less "fixed" on the screen.

[Joel]

Oh, yeah, right...Didn't realize that that button would practically do the same!

One thing though: It is so much fun to see, how you deal with the collision of two planets, hehehe !!!

[GeorgeMcGinn]

Henko is correct.

Newton's Law, without getting into too much "over the head" science, states that two bodies in motion around each other is known as the "gravitational two-body problem" or from the Ancient Greek name the Barycenter (https://en.m.wikipedia.org/wiki/Barycenter).

The program got it right. When two bodies of unequal mass rotate around each other, several things are happening. The dominant body (with the greater mass) controls the orbit of the smaller object, and the smaller object does cause the larger to wobble. Examples are the orbit of Mercury and the technique used to detect exoplanets using the wobble effect. It is the wobble of the star that is being detected, not the planet's.

However, to see the effect on a planet's orbit, check the website. Mercury follows the same kind of orbit, Due to its closeness to the Sun, Mercury goes around in concentric circles.

With the exception of Mercury's orbit, once you calculate the center of gravity between the two bodies, you'll see that the orbit's are not only elliptical, but due to the wobble the bodies do not stay in a center like the one in the program. The program accurately shows how the two bodies drift because the true center of gravity is a point between the two where the pull is equal.

That's about as good as I can put it in layman terms without having to get into the triginomitry and calculus of motion. And yes, I may have gotten some things wrong by simplifying it, but the trade off so many can understand it far outweighs any small mistakes or omissions I may have made.

Besides, I hate eating crow!

George.

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!!

[Joel]

Hey George, yeah... the wobbling effect is very cool... and correct of course. actio= reactio.

But the other thing, namely the instability of a simulated and relatively small planet turning around a huge star is due to very interesting effects in mathematics - numerical mathematics to be precise.

Of course Henko is right by applying newton here. The latter, undoubtedly, is right too (at least, if you neglect relativistic effects, hehe ). Don't get me wrong!! The simulation is very cool!! I love it, there is something happening there, but the simulation only approaches the mathematical solution ( which is ok, as I said)

Just try this out: you will see, that the planet doesn't stay on his elliptical orbit. It gets more and more distant from his central star.

This is due to the linearity of the simulation which uses discrete mathematics to depict a continous world... For a better result you could try for example a 2nd order of Runge Kutta method...

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=2 ! 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
dt=5'IF MIN>5 THEN dt=1 ELSE dt=0.1+.03*MIN*MIN
IF dt=1 THEN
  FOR i=2 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=2 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

[Henko]

I think it is not too difficult to prove that a stationnary orbit of an object around a planet or star is inherently unstable. The radial force on the object with speed "v" is:
F = m*v^2/R - f*m*M/R^2
Equilibrium: F=0 gives the needed speed v' to remain in orbit: v' = sqrt(f*M/R)
If no disturbances take place, the object may remain in orbit for a long time.
Now suppose that your kid is playing with the buttons in the cockpit, and start the rocket engine for a short time. The speed will become a little larger. It's obvious that the radial force F will increase (the partial derivate wilt respect to v is positive). The object will start to accelerate away from the planet, hence R will increase.
Given that the speed is about v', dF/dR is also positive, hence the increase of F will become still stronger.
The same goes for the inverse situation. If the velocity gets smaller then v', the resulting F will be directed towards the planet and the object will eventually crash or burn.
This has nothing to do with the discretization applied in the program. It's systematic instability.
Recently i made a SB simulation of bringing an object in a stationary orbit around the earth, taking into account the specifics of the earth, the drag in the standard atmosphere, and a very critical amount of rocket fuel. So far i did not manage to achieve the target. The reason for that (amongst others) is the instability property around v'. I have to think about (practical) help tools for the poor astronaut (assuming that "Houston" suffers a total blackout).
The "fly-by-wire" solution is of course not yet invented in that little game. That would be to simple. The simulation also needs some graphic upgrading for posting it.