gravity and orbits

[Antoni Gual]

Mango: You could have limited the minimum distance when calculating the force

Webberboy: Fine code!

I did a blend of all ideas. Webberboy multiple balls + mango interactions from all objects+ minimum distance limitation:

Code:

SCREEN 12
TYPE balltype
x AS SINGLE
y AS SINGLE
speedx AS SINGLE
speedy AS SINGLE
accelx AS SINGLE
accely AS SINGLE
mass AS SINGLE
rad AS INTEGER
END TYPE

CONST nballs = 4
'play with these constants
CONST k = -100
CONST dt = .2

DIM b(nballs) AS balltype
DIM indx(nballs, nballs) AS INTEGER

DATA 320,240,0,0,100
DATA 320,140,.5,.3,10
DATA 320,350,-1,-.0,20
DATA 350,200,.1,.3,10

'get ball's parameters
FOR i% = 1 TO nballs
READ b(i%).x, b(i%).y, b(i%).speedx, b(i%).speedy, b(i%).mass
b(i%).rad = SQR(b(i%).mass)
NEXT

'erase accelerations
DO
ERASE indx
FOR i% = 1 TO nballs
b(i%).accelx = 0
b(i%).accely = 0
'a ball can't accelerate itself
indx(i%, i%) = -1
NEXT

'calculate interactions between every pair of balls
FOR i% = 1 TO nballs
  FOR j% = 1 TO nballs
    'if this pair is not already calculated
    IF NOT indx(i%, j%) THEN
      'we're calculating it now!
      indx(i%, j%) = -1
    
      dx = b(i%).x - b(j%).x
      dy = b(i%).y - b(j%).y
      'square of distance
      dist2 = (dx * dx + dy * dy)
      'distance
      dist = SQR(dist2)
      'limit minimum distance to avoid force becoming infinite
      IF dist < 10 THEN dist = 10: dist2 = 100
      'calculate force between balls i and j
      force = k * (b(i%).mass + b(j%).mass) / dist2
      'add this force's acceleration to both balls
      b(i%).accelx = b(i%).accelx + (force / b(i%).mass * dx / dist)
      b(i%).accely = b(i%).accely + (force / b(i%).mass * dy / dist)
      b(j%).accelx = b(j%).accelx - (force / b(j%).mass * dx / dist)
      b(j%).accely = b(j%).accely - (force / b(j%).mass * dy / dist)
    END IF
  NEXT
NEXT

'update velocities and positions and draw balls
WAIT &H3DA, 8
FOR i% = 1 TO nballs
    b(i%).speedx = b(i%).speedx + b(i%).accelx * dt
    b(i%).speedy = b(i%).speedy + b(i%).accely * dt
    CIRCLE (b(i%).x, b(i%).y), b(i%).rad, 0
    b(i%).x = b(i%).x + b(i%).speedx * dt
    b(i%).y = b(i%).y + b(i%).speedy * dt
    CIRCLE (b(i%).x, b(i%).y), b(i%).rad, i%
NEXT
LOOP UNTIL LEN(INKEY$)

[Mango]

Mango: You could have limited the minimum distance when calculating the force

Yep...I think I actually made it so dist couldn't be zero (as that would cause a div by 0 error), but left it with the ability to get too close. I altered your code to give random initial conditions, and allow for an arbtrary # of particles. Now...if you included a zoom factor...

I think that this change makes your code more flexible since it allows you to add or remove balls much more easily:

Code:

CONST nballs = 16  'however many you want.
'play with these constants
CONST k = -100
CONST dt = .2

DIM b(nballs) AS balltype
DIM indx(nballs, nballs) AS INTEGER

'get ball's parameters
randomize timer
FOR i% = 1 TO nballs
  b(i%).x = RND * 640
  b(i%).y = RND * 480
  b(i%).speedx = 0  'keeps center of mass on-screen
  b(i%).speedy = 0  'keeps center of mass on-screen
  b(i%).mass = 5 + (RND * 95)
  b(i%).rad = SQR(b(i%).mass) / 4
NEXT i%

Then, I let the balls get a littel closer...

Code:

'limit minimum distance to avoid force becoming infinite
IF dist < 4 THEN dist = 4: dist2 = 16

Anyway...nice code. more polished than mine was by a lot. Mine was in '96 when I first got a Pentium (90)...it was about 100 times faster than my old PC, and I was thrilled to be able to write this type of program with 1000 elements and run it in real time!!! Totally unfeasible on my old machine(IBM PS2/30). Now with GHz machines, yours works great with circles instead of points....and a wait to keep it from running too fast ;-)....how things change.

[Antoni Gual]

In fact I used DATA for ball data to be able to find and then hard-code a set of initial conditions where balls remain onscreen for a long time. With random initial conditions probably most balls go out of screen in the first ten seconds...
To build a stable system that behaves as our Solar system is not trivial. We live in a little wonder!

[Mango]

With random initial conditions probably most balls go out of screen in the first ten seconds...!

That's why you need zoom-out!!

Code:

RANDOMIZE TIMER
SCREEN 12
TYPE balltype
x AS SINGLE
y AS SINGLE
speedx AS SINGLE
speedy AS SINGLE
accelx AS SINGLE
accely AS SINGLE
mass AS SINGLE
rad AS INTEGER
END TYPE
scale% = 1
LOCATE 1: PRINT "(L)eft (R)ight (U)p (D)own Zoom(+\-) = "; scale%
CONST nballs = 30
'play with these constants
CONST k = -100
CONST dt = .2

DIM b(nballs) AS balltype
DIM indx(nballs, nballs) AS INTEGER
centerx = 320
centery = 240
'get ball's parameters
FOR i% = 1 TO nballs
  b(i%).x = (RND * 640) - centerx
  b(i%).y = (RND * 480) - centery
  b(i%).speedx = 0  'keeps center of mass on-screen
  b(i%).speedy = 0  'keeps center of mass on-screen
  b(i%).mass = 5 + (RND * 95)
  b(i%).rad = b(i%).mass / 10
NEXT i%


'erase accelerations


DO
ERASE indx
FOR i% = 1 TO nballs
b(i%).accelx = 0
b(i%).accely = 0
'a ball can't accelerate itself
indx(i%, i%) = -1
NEXT

'calculate interactions between every pair of balls
FOR i% = 1 TO nballs
  FOR j% = 1 TO nballs
    'if this pair is not already calculated
    IF NOT indx(i%, j%) THEN
      'we're calculating it now!
      indx(i%, j%) = -1
      
      dx = b(i%).x - b(j%).x
      dy = b(i%).y - b(j%).y
      'square of distance
      dist2 = (dx * dx + dy * dy)
      'distance
      dist = SQR(dist2)
      'limit minimum distance to avoid force becoming infinite
      IF dist < 10 THEN dist = 10: dist2 = 100
      'calculate force between balls i and j
      force = k * (b(i%).mass + b(j%).mass) / dist2
      'add this force's acceleration to both balls
      b(i%).accelx = b(i%).accelx + (force / b(i%).mass * dx / dist)
      b(i%).accely = b(i%).accely + (force / b(i%).mass * dy / dist)
      b(j%).accelx = b(j%).accelx - (force / b(j%).mass * dx / dist)
      b(j%).accely = b(j%).accely - (force / b(j%).mass * dy / dist)
    END IF
  NEXT
NEXT

a$ = LCASE$(INKEY$)
IF a$ <> "" THEN
   SELECT CASE a$
   CASE "-"    'zoom out
      scale% = scale% * 2
      centerx = 320: centery = 240
   CASE "="    'zoom in
      IF scale% >= 2 THEN
          scale% = scale% \ 2
      ELSE      'recenter
          centerx = 320
          centery = 240
      END IF
   CASE "q"    'quit
      END
   CASE "u"    'move center up
      centery = centery - 50
   CASE "d"    'move center down
      centery = centery + 50
   CASE "l"    'left
      centerx = centerx - 50
   CASE "r"    'right
      centerx = centerx + 50
   END SELECT
   CLS
   LOCATE 1: PRINT "(L)eft (R)ight (U)p (D)own Zoom(+\-) = "; scale%
END IF
'update velocities and positions and draw balls
WAIT &H3DA, 8

FOR i% = 1 TO nballs
    b(i%).speedx = b(i%).speedx + b(i%).accelx * dt
    b(i%).speedy = b(i%).speedy + b(i%).accely * dt
    CIRCLE (centerx + b(i%).x \ scale%, centery + b(i%).y \ scale%), b(i%).rad \ scale%, 0
    b(i%).x = b(i%).x + b(i%).speedx * dt
    b(i%).y = b(i%).y + b(i%).speedy * dt
    CIRCLE (centerx + b(i%).x \ scale%, centery + b(i%).y \ scale%), b(i%).rad \ scale%, 1 + (i% MOD 16)
NEXT
LOOP

PS: Why does nballs *need* to be a const? I changed it to a user-defined variable so the num of "planets" could be specified by the user at runtime, and it choked. I don't understand.

[Antoni Gual]

nballs does not need to be a constant. Check your %....

Tonight at home i will check your Zoom. We're doing a planetarium online.... What about 3d?

[Mango]

nballs does not need to be a constant.

I dunno...when I make it a non-const, I get a subscript out of range error on this line:

indx(i%, i%) = -1

it may have to do with zero-indexing the arrays, but I can't figure it out...

Tonight at home i will check your Zoom. We're doing a planetarium online.... What about 3d?

3d should be pretty straight forward...one cavat, scaled z-position might be visually confused with (size~mass).

[Antoni Gual]

indx(i%, i%) = -1

I found out. By changing nballs to a variable, the array indx becomes DYNAMIC, so ERASE does actually erase it.
ERASE should be a REDIM in this case...

But you can get rid of all lines where indx appears. by doing:

Code:

FOR i% = 1 TO nballs
  FOR j% = i%+1 TO nballs

[oracle]

You guys should check out plntwars.bas from qbasic.com or QBNZ. It's a whole game based around gravity of planets with different sizes bending the shots you fire etc. It's an ancient game, and the algo for working out where the shot should go next is so ancient the program doesn't need a delay, but it shows what you can do

[Mango]

You guys should check out plntwars.bas from qbasic.com or QBNZ...

The thing is, we weren't really looking for the answer...we were simply enjoying the process of rediscovery...cheers

Antonio...I was looking at the code...it can be simplified significantly, I believe, but I'm into a bottle of wine at the moment, and my algebra focus appears lost. However, I think that the total number of equations can be substantially reduced, and the whole "index" thing can be eliminated by doing

Code:

for i = 0 to numcirc-1
   for j = i + 1 to numcirc
     do calcs between i and j here
   next j
next i

edit...lol...sorry antonio...I see that this is the exact change you proposed. I read your post on my main computer, then went to another (that was off-line) to understand the code. I posted the previous once I understood your code, without re-reading your earlier post...again...rediscovery is always a pleasure...lol

[Antoni Gual]

Oracle: I will check it!

Mango: I hope you enjoyed your wine!

I tried the 3d stuff, but with only circles there 's no 3d effect at all. I will try to add some shadowing...