14. N-Particle


General Description

The gravitational attraction between any two masses is given by the formula:

  
 F = G * m1 * m2
  D*D

In the program the forces on each mass due to all the other masses are added up so that its new position after time 't' may be estimated. To minimise cumulative errors, all masses are moved at the same time - i.e. all the forces on all the masses and all their positions (calculated using the old positions), are calculated before they are all moved. (Thus, at least:
   Force of Mass A on Mass B = Force Mass B on Mass A
i.e. at least the momentum is conserved even for long sampling times.)
   A further refinement in the program is that rather than using the instantaneous acceleration (which will always be too little when masses are approaching and too large when they are separating), the mean of the two latest and the previous acceleration is used (except, of course, for the first calculation).

Detailed Description

The main program consists of three nested loops.
   The innermost loop takes the mass at a time (Ms, for example) and its data, and goes round the loop (N - 1) times to calculate the total forces exerted on Ms by the other (N - 1) masses. (This loop having, of course, the 'OLD' positions of all the masses.)
   The next loop out feeds masses one at a time to the inner loop, and when each emerges with its X & Y acceleration calculates its 'NEW' velocity (P & Q) components and coordinates (stored in ALSO(S, O) and ALSO(S, 1)).
   When all the masses have been processed, the outermost loop uses a mini-loop (360-370) to update all the information used by the other loops, and returns to the start.
   Lines 10-120 Descriptions, instructions and set up arrays etc.
   130-200 Accepts DATA and stores it in DAT(10, 5) (Notice that the sampling time may be changed if greater accuracy is required).
   210 Outer Loop.
   220 S=0 This resets the count for the centre loop - i.e. starts it all again with the new data (updated in the lower half of the loop).
   230 Centre Loop (to 350).
The top half of this loop takes one mass at a time, (S=0 to n)Ms, for example, gives the variables M, X, Y, U and V its coordinates and velocity from SAT, then feeds those variables to the inner loop (it also resets the count with H=0).
   250 Inner Loop.
This inner loop finds the forces due to the H mass on the mass Ms for all the masses - i.e. each time round it calculates the force on Ms due to another mass and adds it to the previous total.
   290 B & C hold the running total of the X & Y acceleration.
   320 ACC holds accelerations. The old values are transformed to E & F, before ACC is updated (with B & C). B & C are then recalculated using B & E or C & F (to get the average over the time interval 't').
   K is 0 the first time round and thereafter 0.5.
   330 P & Q are the new components of velocity calculated using

v = u + ft
They are then stored directly in DAT. The new positions calculated using
s = ut + ½ft ^ 2
are stored in the array ALSO until all particles have been treated.
   360-380 Mini-loop updates all masses by transferring ALSO to DAT.

Educational Notes

In the classroom this program will prove invaluable for the opportunities it offers teachers to provide students with a lucid demonstration of the laws of gravity. If restricted to two masses with the initial momentum at zero, the masses can orbit in quite stable ellipses about their common centre of gravity. This, of course, is the acid test!
   The program also has the useful capacity to demonstrate - or, more precisely, underline - the conservation of momentum even when it goes wrong - i.e. when a collision occurs.