2 //import java.util.Enumeration;
5 * A class that represents the root of the data structure used
6 * to represent the N-bodies in the Barnes-Hut algorithm.
15 * A reference to the root node.
19 * The complete list of bodies that have been created.
23 * The complete list of bodies that have been created - in reverse.
25 public Body bodyTabRev;
28 * Construct the root of the data structure that represents the N-bodies.
30 public Tree(double DTIME)
32 rmin = new MathVector();
46 * Return an enumeration of the bodies.
47 * @return an enumeration of the bodies.
49 /*final Enumeration bodies()
51 return bodyTab.elements();
55 * Return an enumeration of the bodies - in reverse.
56 * @return an enumeration of the bodies - in reverse.
58 /*final Enumeration bodiesRev()
60 return bodyTabRev.elementsRev();
64 * Random number generator used by the orignal BH benchmark.
65 * @param seed the seed to the generator
66 * @return a random number
68 public double myRand(double seed)
70 double t = 16807.0*seed + 1;
72 double iseed = t - (2147483647.0 * Math.floor(t / 2147483647.0f));
77 * Generate a doubleing point random number. Used by
78 * the original BH benchmark.
80 * @param xl lower bound
81 * @param xh upper bound
83 * @return a doubleing point randon number
85 public double xRand(double xl, double xh, double r)
87 double res = xl + (xh-xl)*r/2147483647.0;
92 * Create the testdata used in the benchmark.
93 * @param nbody the number of bodies to create
95 public void createTestData(int nbody)
97 MathVector cmr = new MathVector();
98 MathVector cmv = new MathVector();
100 Body head = new Body();
103 double rsc = 3.0 * Math.PI() / 16.0;
104 double vsc = Math.sqrt(1.0 / rsc);
106 //Random rand = new Random((long)seed);
107 //int max_int = ~(int)0x1+1;
109 for (int i = 0; i < nbody; i++) {
117 //seed = Math.abs((double)rand.nextInt()/max_int);
118 double t1 = xRand(0.0, 0.999, seed);
119 t1 = Math.pow(t1, (-2.0/3.0)) - 1.0;
120 double r = 1.0 / Math.sqrt(t1);
123 for (int k = 0; k < cmr.NDIM; k++) {
125 //seed = Math.abs((double)rand.nextInt()/max_int);
126 r = xRand(0.0, 0.999, seed);
127 p.pos.value(k, coeff*r);
135 //seed = Math.abs((double)rand.nextInt()/max_int);
136 x = xRand(0.0, 1.0, seed);
138 //seed = Math.abs((double)rand.nextInt()/max_int);
139 y = xRand(0.0, 0.1, seed);
140 } while (y > (x*x * Math.pow((1.0f - x*x), 3.5)));
142 double v = Math.sqrt(2.0) * x / Math.pow(1 + r*r, 0.25);
144 double rad = vsc * v;
147 for (int k = 0; k < cmr.NDIM; k++) {
149 //seed = Math.abs((double)rand.nextInt()/max_int);
150 p.vel.value(k, xRand(-1.0, 1.0, seed));
152 rsq = p.vel.dotProduct();
154 double rsc1 = rad / Math.sqrt(rsq);
155 p.vel.multScalar(rsc1);
161 // toss the dummy node at the beginning and set a reference to the first element
162 bodyTab = head.getNext();
164 cmr.divScalar(nbody);
165 cmv.divScalar(nbody);
168 Body b = this.bodyTab;
170 b.pos.subtraction(cmr);
171 b.vel.subtraction(cmv);
176 // set the reference to the last element
182 * Advance the N-body system one time-step.
183 * @param nstep the current time step
185 public void stepSystem(int nstep)
193 Body b = this.bodyTabRev;
195 b.hackGravity(this.rsize, this.root);
199 vp(this.bodyTabRev, nstep);
204 * Initialize the tree structure for hack force calculation.
205 * @param nsteps the current time step
207 public void makeTree(int nstep)
209 Body q = this.bodyTabRev;
212 q.expandBox(this, nstep);
213 MathVector xqic = intcoord(q.pos);
214 if (this.root == null) {
217 if(root instanceof Body) {
218 Body rootb = (Body) root;
219 this.root = rootb.loadTree(q, xqic, 1073741824/*Node.IMAX*/ >> 1, this);
220 } else if(root instanceof Cell) {
221 Cell rootc = (Cell)root;
222 this.root = rootc.loadTree(q, xqic, 1073741824/*Node.IMAX*/ >> 1, this);
228 if(root instanceof Body) {
229 Body rootb = (Body)root;
231 } else if(root instanceof Cell) {
232 Cell rootc = (Cell)root;
238 * Compute integerized coordinates.
239 * @return the coordinates or null if rp is out of bounds
241 public MathVector intcoord(MathVector vp)
243 MathVector xp = new MathVector();
245 double xsc = (vp.value(0) - rmin.value(0)) / rsize;
246 if (0.0 <= xsc && xsc < 1.0) {
247 xp.value(0, Math.floor(1073741824/*Node.IMAX*/ * xsc));
252 xsc = (vp.value(1) - rmin.value(1)) / rsize;
253 if (0.0 <= xsc && xsc < 1.0) {
254 xp.value(1, Math.floor(1073741824/*Node.IMAX*/ * xsc));
259 xsc = (vp.value(2) - rmin.value(2)) / rsize;
260 if (0.0 <= xsc && xsc < 1.0) {
261 xp.value(2, Math.floor(1073741824/*Node.IMAX*/ * xsc));
268 public void vp(Body p, int nstep)
270 MathVector dacc = new MathVector();
271 MathVector dvel = new MathVector();
272 double dthf = 0.5 * this.DTIME;
276 MathVector acc1 = (MathVector)b.newAcc.clone();
278 dacc.subtraction(acc1, b.acc);
279 dvel.multScalar(dacc, dthf);
280 dvel.addition(b.vel);
281 b.vel = (MathVector)dvel.clone();
283 b.acc = (MathVector)acc1.clone();
284 dvel.multScalar(b.acc, dthf);
286 MathVector vel1 = (MathVector)b.vel.clone();
288 MathVector dpos = (MathVector)vel1.clone();
289 dpos.multScalar(this.DTIME);
290 dpos.addition(b.pos);
291 b.pos = (MathVector)dpos.clone();
293 b.vel = (MathVector)vel1.clone();