The magazine of the Melbourne PC User Group Dancing Stars - Part 2 A Four-Star Configuration Ken Holmes

If you put three or more stars on the screen with random velocities and positions, they will almost certainly rapidly disappear. Two stars are no problem; unless they have escape velocity relative to each other, they will become a binary pair in perfectly elliptical orbits and, provided their Nett linear momentum is zero, they will stay on screen. In Figure 3, we show how they could be in concentric circular orbits and we can calculate the velocities required for this. This provides a starting point to then get a range of elliptic orbits by applying a factor to the velocities; here the factors are less than 1, but they could be greater than 1.

We will have to deal with stars of different masses and they will have a centre of gravity which is easily calculated. Here, we have a mass ratio of 2:1 and the CG will be a third of the way from the heavier star to the lighter. If we put them in circular orbits, the lighter star will orbit at twice the radius and at twice the velocity of the heavier. The mutual attractive force is proportional to the product of the masses and inversely proportional to the square of the separation distance. This force is divided by the mass of either to get its gravitational acceleration towards the other, ie towards the centre. This, of course, provides the centripetal acceleration and, to get a circle, a well known formula (a=v^2/r) allows us to find the required circumferential velocity. Since the velocities vary inversely with the masses, their linear momenta (mass x velocity) are equal and opposite, ie. Nett is zero and the orbits stay put on the screen. If you were to add a common velocity to both stars, in plane or out of plane, they would move off, tracing out cycloidal or helical paths. It is thought that more than half of stars in the galaxy are in a binary relationship; quite possibly, most of the rest have at some time had the experience before losing their partner to an interloper. For simplicity, we will adopt a time interval of 1 and also make the Universal Gravitational Constant equal to 1 (such omnipotence!).

Two Twos Are Four

To extend our code to more stars, of differing masses, and in three dimensions, we consider two pairs of binaries orbiting each other. In Figure 4, we see the pattern after they have completed one major circuit. In Figure 5, we have run the program four times as long until chaos developed, breaking up the cosy arrangement; Binary 1 & 2 broke up first and then disturbed binary 3 & 4 which stayed together but shifted their path.

Here we use a perspective presentation on the 2D screen, which is not very satisfactory but we are building up to the next article, where we will use a stereo method, the mirror, which is much more effective in letting us see the pattern more clearly in 3D. The main action occurs near the plane y = 150, extending into the screen and out in front, so we choose a viewing point at eyex = 320, eyey = -300, eyez = -600 to get a perspective view looking down at the plane so that it fills the screen. To plot each star position a simple geometric calculation determines the intersection of the line of sight with the screen surface. See Figure 6 for the side view to explain how screen y, Sy, is calculated. The screen x is similarly found.

In Listing 2, we put the four stars in line across the screen. Stars 1 and 2 form the first binary on the left, and stars 3 and 4 form the second binary on the right. With the separation between the pairs only about 5 times the separation within the binaries, as here, interactions will soon distort the binaries' orbits and, after a few sequences, chaos will take over. To maintain a stable pattern, the separation ratio should be more like 10; this program allows investigation of this subject. Once you have the program working, it is a simple matter to vary star masses or positions.

Our system of axes will have the origin at screen top-left with +x along the top, +y down the left side and +z into the screen. We will use circular orbits for the two binaries but a 0.8 factor for their large orbit, mainly to fit the picture to the screen. Firstly, we calculate the velocities to put stars 1 and 2 in circular orbits and allocate them parallel to the z axis. Then, we do the same for 3 and 4 and allocate the velocities parallel to the y axis. Next, we do the same for the two pairs and allocate the reduced velocities for each pair parallel to the z axis; this means that the pairs will orbit in a horizontal plane. Stars 1 and 2 will orbit each other in the same plane as the big circuit and so trace out epicycloidal paths; stars 3 and 4 will orbit in a vertical plane, which moves around the big circuit, so the paths will vary from cycloidal to helical each quarter revolution.

The Core Code

Extending to three dimensions is just a case of adding z equations identical to the x and y ones, and, to handle the four stars, the FOR...NEXT loops must cover 1 to 4. Now that star masses can differ we can't just calculate where i > j and take the i < j case accelerations as equal but negative,- as we did in Part 1. So using i <> j covers both, still excluding i = j. Excluding the time interval (ie. effectively = 1) simplifies the velocity and position change equations and speeds up the performance. The code will allow you to view an infinite number of situations by simply allotting different values to the masses and x coordinates of any stars. Also, if you change eyey to 150 and extend the n loop to 40000, you will see Figure from within the plane of the big orbit and see that chaos eventually moves the binary of 3 & 4 out of that plane.

You agree that a 2D presentation is not inspiring. In Part 3, we will look at the same situation in true stereo 3D; it is fascinating to see the orbits developing in 3D as the program runs, and extended chaotic runs can still be discerned, even though the screen is a complete mess in 2D.

  ‘Listing 2    dstars2.bas



DEFINT I-K

DIM m(4), x(4), y(4), z(4) 'masses & positions, 4 stars

DIM vx(4), vy(4), vz(4)  'and their velocity components

DIM dsq(4, 4), a(4, 4)   'distance sqrd & accel. i to j

DIM ax(4, 4), ay(4, 4), az(4, 4) 'x, y, z comps of acc.

DIM atx(4), aty(4), atz(4) 'x, y, z comps of total acc.

SCREEN 12           'using VGA screen, as is, no WINDOW

eyex = 320: eyey = -300: eyez = -600      'eye position

m(1) = 1: m(2) = 2: m(3) = 4: m(4) = 2  'easily changed

M12 = m(1) + m(2): M34 = m(3) + m(4)'masses of binaries

f12 = 1: f34 = 1: f14 = .8      ' "circularity" factors

x(1) = 20: x(2) = 105: x(3) = 540: x(4) = 620   'x posn

y(1) = 150: y(2) = 150: y(3) = 150: y(4) = 150  'y posn

z(1) = 0: z(2) = 0: z(3) = 0: z(4) = 0          'z posn

FOR i = 1 TO 4: CIRCLE (x(i), y(i)), 6, i

         PAINT (x(i), y(i)), i, i: COLOR i:

         LOCATE 12, x(i) / 8: PRINT CHR$(48 + i)

NEXT          'marking starting positions of four stars

      'Now get "internal" velocities of binary of 1 & 2

c12 = (m(1) * x(1) + m(2) * x(2)) / M12    'cg of 1 & 2

V1 = f12 * SQR(m(2) * (c12 - x(1)) / (x(2) - x(1)) ^ 2)

V2 = V1 * m(1) / m(2)         'velocity inverse to mass

vx(1) = 0: vx(2) = 0     'velocities parallel to y axis

vy(1) = V1: vy(2) = -V2: vz(1) = 0: vz(2) = 0

      'Now get "internal" velocities of binary of 3 & 4

c34 = (m(3) * x(3) + m(4) * x(4)) / M34    'cg of 3 & 4

V3 = f34 * SQR(m(4) * (c34 - x(3)) / (x(4) - x(3)) ^ 2)

V4 = V3 * m(3) / m(4)         'velocity inverse to mass

vx(3) = 0: vx(4) = 0     'velocities parallel to z axis

vy(3) = 0: vy(4) = 0: vz(3) = V3: vz(4) = -V4

       'Now get "external" velocities for both binaries

c = (M12 * c12 + M34 * c34) / (M12 + M34)    'cg of all

V12 = f14 * SQR(M34 * (c - c12) / (c34 - c12) ^ 2)

V34 = V12 * M12 / M34         'velocity inverse to mass

 'Now add "external" velocities to all, parallel z axis

vz(1) = vz(1) + V12: vz(2) = vz(2) + V12

vz(3) = vz(3) - V34: vz(4) = vz(4) - V34

FOR n = 1 TO 14000         'number of calculation loops

 FOR i = 1 TO 4: FOR j = 1 TO 4

  IF i <> j THEN        'for each star with every other

   xd = x(i) - x(j): yd = y(i) - y(j): zd = z(i) - z(j)

   dsq(i, j) = xd * xd + yd * yd + zd * zd  'dist. sqrd

   d = SQR(dsq(i, j))            'distance between pair

   a(i, j) = m(j) / dsq(i, j)    'gravity accel. i to j

   ax(i, j) = -a(i, j) * xd / d      'x comp. of accel.

   ay(i, j) = -a(i, j) * yd / d      'y comp. of accel.

   az(i, j) = -a(i, j) * zd / d      'z comp. of accel.

  END IF

 NEXT j: NEXT i

 FOR i = 1 TO 4           'Now for each individual star

  atx(i) = 0: aty(i) = 0: atz(i) = 0 'zero total accel.

  FOR j = 1 TO 4           'add x, y & z accel. due ...

   IF i <> j THEN             '... to all "other" stars

    atx(i) = atx(i) + ax(i, j)     'cumulative x accel.

    aty(i) = aty(i) + ay(i, j)     'cumulative y accel.

    atz(i) = atz(i) + az(i, j)     'cumulative z accel.

   END IF

  NEXT j

  vx(i) = vx(i) + atx(i)        'change x velocity of i

  x(i) = x(i) + vx(i)           'change x position of i

  vy(i) = vy(i) + aty(i)        'change y velocity of i

  y(i) = y(i) + vy(i)           'change y position of i

  vz(i) = vz(i) + atz(i)        'change z velocity of i

  z(i) = z(i) + vz(i)           'change z position of i

  zf = -eyez / (-eyez + z(i)) 'z factor for perspective

  Sx = eyex + (x(i) - eyex) * zf     'screen x position

  Sy = eyey + (y(i) - eyey) * zf     'screen y position

  PSET (Sx, Sy), i             'plot star in own colour

 NEXT i

NEXT n                                 'repeat the loop

DO: LOOP WHILE INKEY$ = ""                'There it is!

Note: This series consists of Parts 1 to 6 in PC Update issues from December 2001 to June 2002.  Click to download DSTARS.ZIP file (225 KB) which contains all text, source and executable copies of programs used in this series.

Reprinted from the February 2002 issue of PC Update, the magazine of Melbourne PC User Group, Australia

[About Melbourne PC User Group]