 The magazine of the Melbourne PC User Group Dancing Stars — Part 1 Latest On The 3-Body Problem Ken Holmes |
 |
If you read the series on The Lagrange Points in PC Update, Jul/Aug/Sept 2000, (or were already aware), you will appreciate that these examples of equilibrium arrangements (stable or unstable) are of practical relevance in understanding the presence of the Trojan asteroids clustered on the orbit of Jupiter at the L4/L5 points, 60 degrees ahead of, or behind, the planet. Also, the L1/L2 points are of interest as sites for satellites, with the SOHO satellite already at the sun/earth L1, a million miles from earth, observing sunspots and the solar wind.
However, the 3-body problem is an irresistible challenge to mathematicians and the computer has finally provided the means to find other arrangements of more than two bodies in equilibrium under the influence of their mutual gravitational attractions. These articles are intended to introduce budding programmers to experimentation in this field; I'm sure many will find it fascinating.
Joseph Lagrange, in 1772, devised the situation, in Figure 1, of three stars of identical mass at the corners of an equilateral triangle. Provided they have exactly the right velocity, they will follow the circular orbit maintaining their relative positions - a bit boring - so we also show the effect of reducing the velocities by 20% to give three elliptical paths. The arrangement is, however, unstable; using single precision in QBasic (7 significant figures), accumulated errors allow the paths to diverge and Chaos ensues after about three orbits, with a pair going off in one direction as a binary and the other going in the opposite direction, maintaining the total nett momentum of zero. Since they don't achieve escape velocity they can return and sometimes swap partners. Using double precision (15 sig. figs) lets the original configuration last about seven orbits before errors accumulate. Raising the velocity of only star 1 by a mere 0.1% causes divergence in the second
orbit.
In fact, you can move progressively from an L4 arrangement to this one through a range of equilibrium patterns. Starting with the case of the earth and moon with an L4 body, increasing the moon's mass to that of the earth moves the centre of rotation to the midpoint of the earth/moon line. The L4 will move to a point on the bisector so that the `lead angle' of L4 ahead of the `moon' will be 90 degrees; however, at all stages the three bodies form an equilateral triangle. Ain't maths marvellous! If there seems to be a choice between a simple solution and a complex one, put your money on the simple one. Note that this is still a Lagrange L4 case and is still stable. If we now increase the mass of the L4 body to that of the earth, the centre of rotation must follow the centre of mass of the three bodies and we arrive at the
Figure 1 situation. Along the way we have moved from stable equilibrium to an unstable equilibrium.
|
|
|
What Next?
If you're still interested, look in New Scientist, Aug. 4, 2001. In the last decade some startling configurations have been devised. In 1993, Cristopher Moore at the Santa Fe Institute put three equal mass stars through a figure-of-eight circuit which is, amazingly, stable. It provides the basis for
Figure 2 and for Code Listing 1. In the last few years, several workers have contrived large numbers of paths with differing numbers of stars; these are all stars of equal mass, in two dimensions only, and are most unlikely to occur in reality but are nevertheless intriguing. About 20 are illustrated in New Scientist and it includes Web addresses where more info is available.
Our Program
The three-star figure-of-eight allows us to introduce the simplest possible version of coding to handle multiple star orbits. Having equal masses and only two dimensions will allow us to concentrate on the use of arrays to handle the stars' parameters, with the one set of code serving all stars. We won't attempt to relate our units to reality; the stars will perform on the VGA 640 x 480 screen with the origin at top left; all stars will have unit mass; our arbitrary 'gravitational constant' (set at 15) will, in fact, be the multiplication of a star mass by the Universal Gravitational Constant - in our own strange units. Our old friend, trial-and-error, will determine the appropriate magnitude of velocities to keep the stars on screen. See Listing 1.
'dstars1.bas (Listing 1 from Part 1, Dancing stars)
DEFINT I-K
DIM x(3), y(3), vx(3), vy(3) 'position & veloc. comps
DIM dsq(3, 3) 'distance squared i to j
DIM a(3, 3) 'gravitational acceleration of i to j
DIM ax(3, 3), ay(3, 3) 'x & y comps of acc. of i to j
DIM atx(3), aty(3) 'components of total accel. of i
SCREEN 12: g = 15: va = .108: dy = 94: ti = .1
DO
x(1) = 620: x(2) = 170: x(3) = x(2) 'star positions
y(1) = 200: y(2) = y(1) + dy: y(3) = y(1) - dy
r = (x(1) - x(2)) / (y(2) - y(3)) 'aspect ratio
vbx = va * r 'x (+-)velocity component of 2 & 3
CLS : LOCATE 2
PRINT " 1 = va+, 2 = va-, 3 = dy+, 4 = dy-"
PRINT " dy = "; dy; " ie. r = "; r
PRINT " va = "; va; " vbx = "; vbx
LOCATE 14, 77: PRINT "1": LOCATE 18, 22: PRINT "2"
LOCATE 8, 22: PRINT "3"
vx(1) = 0: vx(2) = vbx: vx(3) = -vbx
vy(1) = va: vy(2) = -va / 2: vy(3) = vy(2)
FOR i = 1 TO 3: CIRCLE (x(i), y(i)), 6, i
PAINT (x(i), y(i)), i, i: NEXT
FOR n = 1 TO 80000
FOR i = 1 TO 2: FOR j = 2 TO 3
IF i < j THEN 'i>j case is negative of this
xd = x(i) - x(j): yd = y(i) - y(j) 'x/y diffs
dsq(i, j) = xd * xd + yd * yd 'distance squared
d = SQR(dsq(i, j)) 'distance
a(i, j) = g / dsq(i, j) 'grav. accel. i to j
ax(i, j) = a(i, j) * (x(j) - x(i)) / d 'x comps..
ax(j, i) = -ax(i, j) '...of i to j and of j to i
ay(i, j) = a(i, j) * (y(j) - y(i)) / d 'y comps..
ay(j, i) = -ay(i, j) '...of i to j and of j to i
END IF
NEXT j: NEXT i
FOR i = 1 TO 3: atx(i) = 0: aty(i) = 0 'zero totals
FOR j = 1 TO 3 'add x & y acc. due to all others
IF i <> j THEN atx(i) = atx(i) + ax(i, j) 'cum. x
IF i <> j THEN aty(i) = aty(i) + ay(i, j) 'cum. y
NEXT j
vx(i) = vx(i) + atx(i) * ti 'change x vel. of i
x(i) = x(i) + vx(i) * ti 'change x pos. of i
vy(i) = vy(i) + aty(i) * ti 'change y vel. of i
y(i) = y(i) + vy(i) * ti 'change y pos. of i
PSET (x(i), y(i)), i 'mark position of i
NEXT i
NEXT n
DO: k$ = INKEY$: LOOP WHILE k$ = ""
SELECT CASE k$
CASE "1": va = va + .0001 'give a new plot with...
CASE "2": va = va - .0001 '....altered values
CASE "3": dy = dy + .5
CASE "4": dy = dy - .5
CASE ELSE 'any key gives a replot with same values
END SELECT
LOOP WHILE k$ <> " " 'SPACE key to exit program
Click to download DSTARS1.ZIP for source listing
and executable version of this program.
|
Getting Started - Basic Considerations
The path is symmetrical and we start the stars at the corners of an isosceles triangle, as shown, but we can only guess at the exact ratio of the lateral separation of star 1 from 2 & 3 to the vertical displacements of stars 2 & 3. We want the pattern to stay on screen, so the nett total momentum must be zero in both x and y directions. The upwards (y) velocity of 1 is called va; from symmetry the downwards velocities of 2 & 3 are equal and therefore must be va/2. The x velocity of 1 is zero, so the x velocities of 2 & 3 must be equal and opposite.
We don't want the pattern to rotate so the nett angular momentum of all stars must be zero. Consider angular momenta (velocity X offset) about the position of star 1; having zero offset, star 1's angular momentum is zero; now the velocity lines of both stars 2 & 3 must also pass through Star 1 since, if they were offset, their angular momenta would be additive giving nett rotation. As we know their y components (va/2), geometry will give their x components
(+-vbx).
In Listing 1, you may note that we provide the means to vary the y displacements of 2 & 3 (dy, set at 94) and va, set at 0.108, so that the program could be run repeatedly with these varied until we could get closest to the correct path. We have not achieved the exact single path since the true values of dy and va need many more decimal places, which would take a lot of trial-and-error. Mathematicians use more sophisticated techniques to get more accurate values. The arrangement is quite stable and alteration of the values and running extended orbits just gives thickened figure eights; you might wish to examine the range of stability. You may recall that we had L4 bodies wandering over large areas, just as the clusters of Trojan satellites do, without all disappearing.
Those Arrays
Using arrays can be initially daunting but they are the only way to go when handling large numbers of similar entities. We DIMension several arrays including some two dimensional ones. For example, dsq(i, j) represents the square of the distance between stars i & j - and, of course, equals dsq(j, i). Our coding uses these arrays and is enclosed in nested FOR...NEXT loops to deal with the gravitational acceleration of each star towards every other star. In further articles we will extend the code itself to stars of differing masses and to three dimensions, and we can also deal with any number of stars by simple increasing the array sizes and extending the range of the FOR...NEXT loops.
The Core Of the Code
Firstly, we consider every pair of stars but must avoid pairing any star with itself since dsq(i, i) is zero and will cause an error if we divide g by it. In this case, where all masses are equal, we can avoid calculating the acceleration of j to i since it is minus the acceleration of i to j. We have eliminated masses from our equations since they are all unity. Dividing g by dsq(i, j) gives the magnitude of the acceleration and we then resolve it into x and y components by multiplying by, for example, the x separation divided by the direct distance. The sign of the x or y separation determines the direction of the acceleration which is, of course, opposite for each star of a pair.
Next, in a separate nesting of FOR...NEXT loops, the acceleration components for each star are totalled. Then the most basic equations of motion are used to adjust each star's velocity and then its position. Note that we have a time interval, ti, of 0.1; we could make the code simpler by making it 1 to avoid many calculations. It is left in to remind first-timers that an acceleration should be multiplied by a time to get a velocity change and a velocity should be multiplied by a time to get a position change. You may eliminate it entirely from the motion equations, thus making it effectively equal to 1; then the 80000 iterations will give ten times as many orbits and rapidly fill in the thickened eight. A virtue of using ti = 0.1 is that the plotting of the orbit is slowed down to make it more interesting.
Next Month
I considered trying a five-star pattern from New Scientist but these involve more velocity and position variables which must be decided with great accuracy since the orbits are unstable;
I haven't as yet achieved a reasonable first orbit. We will try another tack to extend the versatility of the coding.
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 December 2001 issue of PC Update, the
magazine of Melbourne PC User Group, Australia
