void Interact(body *b1,vector b2pos,double b2mass)
{
double maga; /* magnitude of acceleration */
double r; /* distance from b1 to b2 */
double r2; /* r*r */
vector b1b2; /* vector from b1 to b2 (b2.pos - b1.pos) */
vector da; /* the change that will be applied to the acceleration vector (delta a) */
/* b1 is the body that will be moved, b2 is the body causing the change in acceleration */
/* Given that F = (G m1 m2) / (R^2), and F = m1 a, we can calculate
|a| = G m2 / R^2, and then use that magnitude change to adjust the acceleration
along the direction of the interaction */
SUBV(b1b2,b2pos,b1->pos);
DOTVP(r2,b1b2,b1b2);
r2 += constants.epssq; /* Weaken closeness of the two objects */
r = sqrt(r2);
maga = G * b2mass / (r2*r); /* Account for b1b2 not being a unit vector
in maga */
MULVS(da,b1b2,maga); /* make da the right length */
ADDV(b1->acc,b1->acc,da); /* Accumulate the acceleration */
}
This function will correctly calculate the interaction of any two bodies.
The tree-walk optimization pre-computes the center of mass of a group of
bodies and calculates the interaction of a single body with that entire
group all in a single step. The criteria for whether or not to recurse
down the oct-tree or perform the interaction with the entire group looks like:
void TW_DownWalk(node *global root,body *particle,double dsq)
{
SUBV(dr,root->pos,particle->pos);
DOTVP(drsq,dr,dr);
if (constants.tolsq * drsq < dsq) {
/* need to open node, i.e. recurse down the tree */
} else {
/* perform the interaction */
}
}
The important piece to note is that the Interaction function uses
the mass of the second body in the calculation to determine the
amount of acceleration, but the opening criteria (whether or not to
recurse) does not take into account the accumulated mass at the
node, but only the tolerance. If the simulation is sufficiently well
behaved to keep large clumps from occurring, then having a single
tolerance is sufficient. However, since most of these simulations
move from randomly distributed particles to clustered particles, it
is important to take the mass into account. Our implementation does not
do so, and as a result has the same problem with unbounded error as many
other implementations. An area for further work would be to either use
the opening criteria in Warren and Salmon's paper, or work out another
criteria for opening the node. Further reduction in the amount of
work for the Barnes-Hut algorithm could be found by performing node-node
interactions when the two nodes are sufficiently far apart as to minimize
the error.
-----
maintained by <hodes@cs.berkeley.edu>