1 #include <omp.h>
    2 /******************************************************************************/
    3 /*                                                                            */
    4 /*  Copyright (c) 1996 University of California at Santa Barbara              */
    5 /*                                                                            */
    6 /*  All rights reserved.                                                      */
    7 /*                                                                            */
    8 /*  Permission is given to use, copy, and modify this software for any        */
    9 /*  non-commercial purpose as long as this copyright notice is not removed    */
   10 /*  All other uses, including redistribution in whole or in part, are         */
   11 /*  forbidden without prior written permission.                               */
   12 /*                                                                            */
   13 /*  This software is provided with absolutely no warranty and no support.     */
   14 /*                                                                            */
   15 /*  For further information contact: pedro@cs.ucsb.edu or martin@cs.ucsb.edu  */
   16 /*                                                                            */
   17 /******************************************************************************/
   18 
   19 #include <stdio.h>
   20 #include <string.h>
   21 #include <math.h>
   22 #include <stdlib.h>
   23 #include <iostream>
   24 #include <fstream>
   25 #include <limits.h>
   26 
   27 #define CHECK
   28 using namespace std;
   29 
   30 // const int NULL = 0;
   31 // const int  FALSE = 0;
   32 // const int  TRUE  = 1;
   33 
   34 const double    PI = 3.14159265358979323846;
   35 // const double    TWO_PI = 6.28318530717958647693;
   36 // const double    FOUR_PI = 12.56637061435917295385;
   37 // const double    HALF_PI = 1.57079632679489661923;
   38 // const double    FRTHRD_PI = 4.18879020478639098462;
   39 
   40 const int  MAX_NAME = 80;
   41 
   42 extern "C" void reset_stat();
   43 extern "C" void print_stat();
   44 
   45 extern "C" void         get_time(int *r);
   46 extern "C" void         get_ticks_time(int *r);
   47 
   48 extern "C" double       prand();
   49 extern "C" double       xrand(double xl, double xh);
   50 extern "C" void         pranset(int seed);
   51 
   52 extern "C" int          scanbind(char **bvec, char *name);
   53 extern "C" int          matchname(char *bind, char *name);
   54 extern "C" char*        extrvalue(char *arg);
   55 #define _OPENMP_
   56 #include "water.h"
   57 //const int NDIM   = 3;
   58 inline void in_vector(ifstream& f, double *v){
   59   double x1, x2, x3;
   60   f >> x1 >> x2 >> x3;
   61   v[0] = x1; v[1] = x2; v[2] = x3;
   62 }
   63 
   64 inline void tvecClr(double *v) { int i; for(i=0; i < NDIM; i++) v[i] = 0.0; }
   65 
   66 inline void tvecAdd(double *v1, double *v2) { int i;
   67 for(i=0; i < NDIM; i++) v1[i] += v2[i]; }
   68 
   69 inline void tvecSub(double *v1, double *v2) { int i;
   70 for(i=0; i < NDIM; i++) v1[i] -= v2[i]; }
   71 
   72 inline void tvecLoad(double *v1, double *v2) { int i;
   73 for(i=0; i < NDIM; i++) v1[i] = v2[i]; }
   74 
   75 inline void tvecStore(double *v1, double *v2) { int i;
   76 for(i=0; i < NDIM; i++) v2[i] = v1[i]; }
   77 
   78 inline void tvecDiv(double *v1, double *v2) { int i;
   79 for(i=0; i < NDIM; i++) v1[i] /= v2[i]; }
   80 
   81 inline void tvecProd(double *v1, double *v2) { int i;
   82 for(i=0; i < NDIM; i++) v1[i] *= v2[i]; }
   83 
   84 inline void tvecScale(double *v, double s) { int i;
   85 for(i=0; i < NDIM; i++) v[i] *= s; }
   86 
   87 inline double tvecDotProd(double *v1,double *v2) { int i; double d = 0.0;
   88   for(i=0; i < NDIM; i++) d += (v1[i] * v2[i]); return d; }
   89 
   90 inline void tvecCrossProd(double *v, double *u, double *w){
   91   v[0] = u[1]*w[2] - u[2]*w[1];
   92   v[1] = u[2]*w[0] - u[0]*w[2];
   93   v[2] = u[0]*w[1] - u[1]*w[0];
   94 }
   95 
   96 inline void tvecMin(double *v1, double  *v2){ int i;
   97   for(i=0; i < NDIM; i++) if(v1[i] > v2[i]) v1[i] = v2[i]; }
   98 
   99 inline double tvecNorm1(double *v){ int i; double d;
  100 d = 0.0; for(i=0; i < NDIM; i++) d += v[i]; return d; }
  101 
  102 
  103 // ----------- FILE: body.h ----------
  104 // BH-tree capacity constants
  105 const int LEAFMAXBODIES = 10;
  106 const int NSUB = 8; 			// (1 << NDIM)
  107 const int IMAX	= 1073741824;		// (1 << MAXLEVEL)
  108 
  109 const int CELL = 01;
  110 const int LEAF = 02;
  111 const int BODY = 04;
  112 
  113 class node {	// inherited by all CELL, LEAF and BODY classes
  114 public:
  115    int type;
  116    double mass;      /* total mass of node */
  117    vector pos;       /* position of node */
  118    int level;
  119    node *parent;     /* ptr to parent of this node in tree */
  120    int child_num;    /* Index that this node should be put at in parent cell */
  121 	node()		{ type = 0; mass = 0.0; pos.vecClr(); level = 0; parent = NULL; child_num =0;}
  122 	~node()		{ }
  123 int 	Type()		{ return type; }
  124 double	Mass()		{ return mass; }
  125 vector  *Pos()  	{ return &pos; }
  126 int	Level()  	{ return level; }
  127 node 	*Parent()  	{ return parent; }
  128 int 	ChildNum()  	{ return child_num; }
  129 void    setChildNum(int n) { child_num = n; }
  130 void	setKind(int t)	{ type = t; }
  131 void	setMass(double m) { mass = m; }
  132 void	setLevel(int d)	{ level =d; }
  133 void	setParent(node *p) { parent = p; }
  134 void	PosvecStore(double *p) { pos.vecStore(p); }
  135 void	PosvecLoad(double *p) { pos.vecLoad(p); }
  136 void	PosvecAdd(double *p) { pos.vecAdd(p); }
  137 void	clear()		{ mass= 0.0; pos.vecClr(); level = 0;
  138 			  parent = NULL; child_num =0; }
  139 void	adjLevel(node *p) {
  140 			  if (p == NULL) level = (IMAX >> 1);
  141 			  else level = (p->Level() >> 1);
  142 			}
  143 int	subindex(int *x, int lev);
  144 void	computecofm();
  145 };
  146 
  147 // ----------------------------------------------------------------------------
  148 // BODY: data structure used to represent particles.
  149 // ----------------------------------------------------------------------------
  150 class body : public node {
  151    vector vel;          /* velocity of body                             */
  152    vector acc;          /* acceleration of body                         */
  153    double phi;          /* potential at body                            */
  154    vector tmp;          /* Skratch vector		                */
  155    double  subdiv(node *p, double *ai);
  156 public:
  157         body()		{ type = BODY; vel.vecClr(); acc.vecClr(); tmp.vecClr(); phi = 0.0; }
  158         ~body()         { }
  159 vector  *Vel()          { return &vel; }
  160 vector  *Acc()          { return &acc; }
  161 double  Phi()           { return phi; }
  162 vector  *Tmp()          { return &tmp; }
  163 void 	BodyClear()	{ node::clear();
  164 			  vel.vecClr(); acc.vecClr(); tmp.vecClr();
  165 			  phi = 0.0;
  166 			}
  167 void	VelvecAdd(double *p)	{ vel.vecAdd(p); }
  168 void	VelvecStore(double *p)	{ vel.vecStore(p); }
  169 void	VelvecLoad(double *p)	{ vel.vecLoad(p); }
  170 void	AccvecStore(double *p)	{ acc.vecStore(p); }
  171 
  172 double  computeInter(node *p, double epsSq, double *res);
  173 void	openCell(node *n, double tolsq, double dsq, double epssq);
  174 void	openLeaf(node *n, double epssq, int iter);
  175 void 	gravsub(node *p, double epssq);
  176 void 	updatePhi(double inc);
  177 void    fastgravsub(node *p, double *ai, double drsq, double epssq);
  178 void   	walksub(node *n, double tolsq, double dsq, double epssq);
  179 void	advance(double hts, double ts);
  180 void	swapAcc();
  181 void 	startVel(double f);
  182 #ifdef CHECK
  183 void dump();
  184 #endif
  185 };
  186 
  187 // ----------------------------------------------------------------------------
  188 // LEAF: structure used to represent leaf nodes of tree.
  189 // ----------------------------------------------------------------------------
  190 class cell;
  191 class leaf : public node {
  192    int num_bodies;
  193    body *bodyp[LEAFMAXBODIES];    // bodies of leaf
  194 public:
  195         leaf()          { type = LEAF; }
  196         ~leaf()         { }
  197 body    *Bodyp(int idx) { return bodyp[idx]; }
  198 void    setBodyp(int idx, body *p) { bodyp[idx] = p; }
  199 int     NumBody()       { return num_bodies; }
  200 void    setNumBody(int n) { num_bodies = n; }
  201 int     full()          { return (num_bodies == LEAFMAXBODIES); }
  202 void    LeafClear()         { int i;
  203                           node::clear();
  204                           num_bodies = 0;
  205                           for (i = 0; i < LEAFMAXBODIES; i++)
  206                             bodyp[i] = NULL;
  207                         }
  208 
  209 void    Subdivide(cell *ce, int idx);
  210 void    LeafAddBody(body *bp);
  211 };
  212 
  213 // ----------------------------------------------------------------------------
  214 // CELL: structure used to represent internal nodes of tree.
  215 // ----------------------------------------------------------------------------
  216 class cell : public node {
  217    node *subp[NSUB];    // descendents of cell
  218 public:
  219         cell()          { type = CELL; }
  220         ~cell()         { }
  221 node    *Subp(int idx)  { return subp[idx]; }
  222 void    setSubp(int idx, node *p)  { subp[idx] = p; }
  223 void    CellClear()         { int i;
  224                           node::clear();
  225                           for (i = 0; i < NSUB; i++)
  226                             subp[i] = NULL;
  227                         }
  228 void CellAddBody(body *p, int *coords);
  229 };
  230 
  231 // ----------------------------------------------------------------------------
  232 // CELLSET: set of cells.
  233 // ----------------------------------------------------------------------------
  234 class cellset {
  235   int	numcell;       	/* no. cells  used in celltab 	*/
  236   int	maxcell;       	/* max no. of cells  allocated 	*/
  237   cell	*celltab;	/* array of cells allocated 	*/
  238 public:
  239   cellset(int max)	{ int i;
  240 			  celltab = new cell[max]; maxcell=max; numcell=0;
  241 			  for(i=0; i < maxcell; i++)
  242 			    celltab[i].clear();
  243 			}
  244   ~cellset()		{ delete[] celltab; maxcell=0; numcell=0; }
  245 void	CellSetclear()	{ int i;
  246 			  for(i=0; i < maxcell; i++)
  247 			    celltab[i].clear();
  248 			}
  249 void 	CellSetReset() 	{ numcell = 0; }
  250  cell	*newcell();
  251  cell	*makecell(node *parent);
  252  cell	*getcell(int idx) { return &celltab[idx]; }
  253 };
  254 
  255 // ----------------------------------------------------------------------------
  256 // LEAFSET: set of leaves.
  257 // ----------------------------------------------------------------------------
  258 class leafset {
  259   int	numleaf;       	/* no. leaves used in leafctab 	*/
  260   int	maxleaf;       	/* max no. of leaves allocated 	*/
  261   leaf	*leaftab;	/* array of leaves allocated 	*/
  262 public:
  263   leafset(int max)	{ leaftab = new leaf[max]; maxleaf=max; numleaf=0; LeafSetClear(); }
  264   ~leafset()		{ delete[] leaftab; maxleaf=0; numleaf=0;}
  265 void	LeafSetClear()	{ int i; for(i=0; i < maxleaf; i++)
  266 			    leaftab[i].clear();
  267 			}
  268 void 	LeafSetReset() 	{ numleaf = 0; }
  269 leaf* newleaf();
  270 leaf* makeleaf(node *parent);
  271 leaf *getleaf(int idx)	{ return &leaftab[idx]; }
  272 // void  add_body(body *bp);
  273 };
  274 
  275 // ----------------------------------------------------------------------------
  276 // BODYSET: set of particle.
  277 // ----------------------------------------------------------------------------
  278 class bodyset {
  279   int   numbody;        /* no. bodies used in bodytab   */
  280   int   maxbody;        /* max no. of bodies allocated  */
  281   body  *bodytab;       /* array of bodies   allocated  */
  282 public:
  283   bodyset(int m)      { int i; bodytab = new body[m]; maxbody=m; numbody=0; for(i=0; i < maxbody; i++) bodytab[i].clear(); }
  284   ~bodyset()            { delete[] bodytab; maxbody=0; numbody=0;}
  285 void  BodySetClear()    { int i; for(i=0; i < maxbody; i++)
  286                             bodytab[i].clear();
  287                         }
  288 void  BodySetReset()         { numbody = 0; }
  289 body* getbody(int idx)  { return &bodytab[idx]; }
  290 body* newbody();
  291   int size() { return numbody; }
  292 };
  293 
  294 class parms;
  295 class NBodySystem {
  296   int nbody;
  297   cell* BH_tree;
  298   cellset* cells;
  299   leafset* leaves;
  300   bodyset* bodies;
  301 
  302   acc_double mtot;
  303   vector *etot;
  304   matrix *keten;
  305   matrix *peten;
  306   vector *cmphasePos;
  307   vector *cmphaseVel;
  308   vector *amvec;
  309 
  310   vector *mincorner;
  311   vector *maxcorner;
  312 
  313   double size;
  314   double tnow;
  315   double tout;
  316 
  317   int nstep;
  318   int n2bcalc;
  319   int nbccalc;
  320   int computestart;
  321   int computeend;
  322   int tracktime;
  323   int partitiontime;
  324   int treebuildtime;
  325   int forcecalctime;
  326   int updatetime;
  327   void pickshell(double *v, double rad);
  328 public:
  329   NBodySystem() {  nbody = 0; cells = NULL; leaves = NULL; bodies = NULL; }
  330   ~NBodySystem() { delete cells; delete leaves; delete bodies; nbody = 0; }
  331   void init();
  332   void intcoord(int *xp, double *rp);
  333   void updateInter(int t) { if(t == BODY) n2bcalc++; else nbccalc++; }
  334 bodyset *getbodies()    { return bodies; }
  335 leafset *getleaves()    { return leaves; }
  336 cellset *getcells()     { return cells; }
  337   void loadData(char *name);
  338   void loadTestData();
  339   void clearStats();
  340   void updateStats();
  341   void outputStats();
  342   void updateEnergy();
  343   void outputEnergy();
  344   void clearTiming();
  345   void outputTiming();
  346   void initOutput(parms *p);
  347   void output(parms *p);
  348   void startrun(parms *p);
  349   void stepsystem(parms *p);
  350   void setbound();
  351   void maketree();
  352   void SwapAccs();
  353   void Advance(double hts, double ts);
  354   void ComputeAccels(double tolerance, double eps);
  355   void StartVels(double halftimestep);
  356   void computeForces(double tol, double hts, double eps);
  357 #ifdef CHECK
  358 void dump();
  359 #endif
  360 };
  361 
  362 // ----------- FILE: simulation.h ----------
  363 class parms {
  364 char    headline[MAX_NAME];/* message describing calculation */
  365 char    infile[MAX_NAME];  /* file name for snapshot input */
  366 char    outfile[MAX_NAME]; /* file name for snapshot output */
  367 double  dtime;          /* timestep for leapfrog integrator */
  368 double  dtout;          /* time between data outputs */
  369 double  tstop;          /* time to stop calculation */
  370 double  fcells;         /* ratio of cells/leaves allocated */
  371 double  fleaves;        /* ratio of leaves/bodies allocated */
  372 double  tol;            /* accuracy parameter: 0.0 => exact */
  373 double  tolsq;          /* square of previous */
  374 double  eps;            /* potential softening parameter */
  375 double  epssq;          /* square of previous */
  376 double  dthf;           /* half time step */
  377 char    *defaults[12];  /* pairs of name = value */
  378 
  379 public:
  380         parms();
  381         ~parms(){ }
  382   int     getiparam(char *name);
  383   double  getdparam(char *name);
  384   void    getparam(char *name);
  385 
  386   void    setparms(char **fnames, double *args);
  387   void    loadInfile();
  388   void    loadParms();
  389 char    *getheadline()  { return headline; }
  390 char    *getinfile()    { return infile; }
  391 char    *getoutfile()   { return outfile; }
  392 double  getdtime()      { return dtime; }
  393 double  getdtout()      { return dtout; }
  394 double  gettstop()      { return tstop; }
  395 double  getfcells()     { return fcells; }
  396 double  getfleaves()    { return fleaves; }
  397 double  gettol()        { return tol; }
  398 double  gettolsq()      { return tolsq; }
  399 double  geteps()        { return eps; }
  400 double  getepssq()      { return epssq; }
  401 double  getdthf()       { return dthf; }
  402 };
  403 
  404 // --- Global variables ---------------------------------------------------
  405 NBodySystem Nbody;
  406 parms simparms;
  407 char    buf[128];
  408 
  409 // ### CODE BELOW  ###################################################
  410 // ### CODE BELOW  ###################################################
  411 // ### CODE BELOW  ###################################################
  412 
  413 
  414 #ifdef CHECK
  415 // ----------------------------------------------------------------------------
  416 void NBodySystem::dump(){
  417   int i;
  418   body *p;
  419   printf ("## mass\t (position)\t (velocity)\t (acceleration)\n");
  420   for (i = 0; i < nbody ; i++){
  421     p = bodies->getbody(i);
  422     p->dump();
  423   }
  424 }
  425 #endif
  426 
  427 #ifdef CHECK
  428 void body::dump(){
  429  printf("%8.4f ",mass);
  430  pos.vecPrint();
  431  vel.vecPrint();
  432  acc.vecPrint();
  433  printf("\n");
  434 }
  435 #endif
  436 
  437 
  438 
  439 // ---------------------------------------------------------------------------
  440 void NBodySystem::init() {
  441   etot = new vector;
  442   keten = new matrix;
  443   peten = new matrix;
  444   cmphasePos = new vector;
  445   cmphaseVel = new vector;
  446   amvec = new vector;
  447   mincorner = new vector;
  448   maxcorner = new vector;
  449 }
  450 
  451 // ---------------------------------------------------------------------------
  452 void NBodySystem::clearStats(){
  453   nstep = 0;
  454   n2bcalc = 0;
  455   nbccalc = 0;
  456 }
  457 
  458 // --- NBodySystem ---------------------------------------------------
  459 void NBodySystem::intcoord(int *xp, double *rp){
  460   int k;
  461   double xsc, v[NDIM];
  462 
  463   mincorner->vecStore(v);
  464   for(k = 0; k < NDIM; k++) {
  465     xsc = (rp[k] - v[k]) / size;
  466       xp[k] = (int)(floor(IMAX * xsc));
  467   }
  468 }
  469 
  470 // ---------------------------------------------------------------------------
  471 void NBodySystem::updateStats(){
  472   int nttot, nbavg, ncavg;
  473   nttot = n2bcalc + nbccalc;
  474   nbavg = (int) ((double) n2bcalc / (double) nbody);
  475   ncavg = (int) ((double) nbccalc / (double) nbody);
  476 }
  477 
  478 void NBodySystem::outputStats(){
  479   cout << "\n\nINTERACTIONS\tTOTAL\tITER\tAVG\n";
  480   cout << "body-body\t\t"<< n2bcalc << "\t" << n2bcalc/nstep<< "\t" <<
  481     (int)sqrt((double)(n2bcalc/nstep)) << "\n";
  482   cout << "body-cell\t\t"<< nbccalc << "\t" << nbccalc/nstep << "\t" <<
  483     (int)sqrt((double)(nbccalc/nstep)) << " \n";;
  484 }
  485 
  486 // ----------------------------------------------------------------------------
  487 void NBodySystem::updateEnergy(){
  488   int k;
  489   double m, velsq;
  490   body *p;
  491 
  492   mtot.writeval(0.0);
  493   etot->vecClr();
  494   keten->matClr();
  495   peten->matClr();
  496   cmphasePos->vecClr();
  497   cmphaseVel->vecClr();
  498   amvec->vecClr();
  499 
  500    #pragma omp parallel for private (velsq, m, p) schedule(static)
  501   for(k = 0; k < nbody; k++) {
  502     p = bodies->getbody(k);
  503     matrix m1;
  504     double v0[3], v1[3], v2[3], v3[3], v4[3];
  505     m = p->Mass();
  506     mtot.addvalRepl(m);
  507     p->VelvecStore(v1);
  508     velsq = tvecDotProd(v1,v1);
  509     v4[1] = 0.5 * m * velsq;
  510     v4[2] = 0.5 * m * p->Phi();
  511     v4[0] = v4[1] + v4[2];
  512     etot->vecAddRepl(v4);
  513     p->PosvecStore(v1);
  514     p->VelvecStore(v2);
  515     p->AccvecStore(v3);
  516     tvecLoad(v0,v2);
  517     tvecScale(v0,0.5*m);
  518     m1.matOutProd(v0,v2);
  519     keten->matAddRepl(&m1);
  520     tvecLoad(v0,v1);
  521     tvecScale(v0,m);
  522     m1.matOutProd(v0,v3);
  523     peten->matAddRepl(&m1);
  524     tvecLoad(v0,v1);
  525     tvecScale(v0,m);
  526     cmphasePos->vecAddRepl(v0);
  527     tvecLoad(v0,v2);
  528     tvecScale(v0,m);
  529     cmphaseVel->vecAddRepl(v0);
  530     tvecCrossProd(v0,v1,v2);
  531     tvecScale(v0,m);
  532     amvec->vecAddRepl(v0);
  533   }
  534 
  535    if (mtot.readval() != 0.0){
  536       cmphasePos->vecScale(1.0/mtot.readval());
  537       cmphaseVel->vecScale(1.0/mtot.readval());
  538    }
  539 }
  540 
  541 
  542 // ----------------------------------------------------------------------------
  543 void NBodySystem::outputEnergy(){
  544   cout << "\n\nENERGY TOTALS\n";
  545   cout << " MASS:        "<< mtot.readval() << "\n";
  546   cout << " ENERGY:     "; etot->vecPrint(); cout << "\n";
  547   cout << " KIN ENERGY: "; keten->matPrint(); cout << "\n";
  548   cout << " POT ENERGY: "; peten->matPrint(); cout << "\n";
  549   cout << " AGGREGATE CM POS: "; cmphasePos->vecPrint(); cout << "\n";
  550   cout << " AGGREGATE CM VEL: "; cmphaseVel->vecPrint(); cout << "\n";
  551   cout << " ANGULAR MOMENTUM: "; amvec->vecPrint(); cout << "\n\n";
  552 }
  553 
  554 // ----------------------------------------------------------------------------
  555 void NBodySystem::output(parms *sim) {
  556   if( (tout - 0.01 * sim->getdtime()) <= tnow)
  557     tout += sim->getdtout();
  558 }
  559 
  560 // ----------------------------------------------------------------------------
  561 void NBodySystem::clearTiming(){
  562    tracktime = 0;
  563    updatetime = 0;
  564    partitiontime = 0;
  565    treebuildtime = 0;
  566    forcecalctime = 0;
  567 }
  568 
  569 // ----------------------------------------------------------------------------
  570 void NBodySystem::initOutput(parms *p){
  571   cout << "\n\t\tOutput: " << p->getheadline() << "\n\n";
  572   cout << "nbody   dtime   eps   tol   dtout   tstop   fcells\n";
  573   cout << nbody << "     " << p->getdtime() << "   " << p->geteps()
  574       << "   " << p->gettol() << "   " << p->getdtout() << "    "
  575       << p->gettstop() << "   " << p->getfcells() << "\n\n";
  576 }
  577 
  578 // ----------------------------------------------------------------------------
  579 void NBodySystem::outputTiming(){
  580   int ticks;
  581   double dticks;
  582 
  583   get_ticks_time(&ticks);
  584   dticks = ticks;
  585   cout << "COMPUTESTART  = " <<
  586     (computestart / dticks) << "\n";
  587   cout << "COMPUTEEND    = " <<
  588     (computeend / dticks) << "\n";
  589   cout << "COMPUTETIME   = " <<
  590     ((computeend - computestart) / dticks) << "\n";
  591   cout << "TRACKTIME     = " <<
  592     (tracktime / dticks) << "\n";
  593   cout << "PARTITIONTIME = " <<
  594     (partitiontime / dticks) << "\t"
  595         << ((float)partitiontime)/tracktime << "\n";
  596   cout << "TREEBUILDTIME = " <<
  597     (treebuildtime / dticks) << "\t"
  598         << ((float)treebuildtime)/tracktime << "\n";
  599   cout << "FORCECALCTIME = " <<
  600     (forcecalctime / dticks) << "\t"
  601         << ((float)forcecalctime)/tracktime << "\n";
  602   cout << "UPDATEPOSTIME = " <<
  603     (updatetime / dticks) << "\t"
  604         << ((float)updatetime)/tracktime << "\n";
  605   cout << "RESTTIME      = " <<
  606     ((tracktime - partitiontime -
  607      updatetime - treebuildtime - forcecalctime)
  608     / dticks) << "\t"
  609         << ((float)(tracktime-partitiontime- updatetime - treebuildtime-
  610         forcecalctime))/tracktime<<"\n";
  611   cout.flush();
  612 }
  613 
  614 // ----------------------------------------------------------------------------
  615 void NBodySystem::loadData(char *name){
  616   int i, ndim;
  617   body *p;
  618   double m, tmp[3];
  619   ifstream instr(name);
  620 
  621   cerr << "Reading input file : " << name << '\n';
  622   cerr.flush();
  623   if (!instr)
  624      cerr << "loadData: cannot find file " << name << '\n';
  625   instr >> nbody;
  626   if (nbody < 1)
  627      cerr << "loadData: nbody = " << nbody << " is absurd\n";
  628   instr >> ndim;
  629   if(ndim != NDIM)
  630     cerr << "inputdata: NDIM = " << NDIM << "ndim = " << ndim << " is absurd\n";
  631   instr >> tnow;
  632 
  633   for(i = 0; i < nbody; i++){
  634     p = bodies->getbody(i);
  635     instr >> m;
  636     p->setMass(m);
  637   }
  638   for(i = 0; i < nbody; i++){
  639     p = bodies->getbody(i);
  640     in_vector(instr,tmp);
  641     p->PosvecLoad(tmp);
  642   }
  643   for(i = 0; i < nbody; i++){
  644     p = bodies->getbody(i);
  645     in_vector(instr,tmp);
  646     p->VelvecLoad(tmp);
  647   }
  648   instr.close();
  649 }
  650 
  651 // ----------------------------------------------------------------------------
  652 void NBodySystem::loadTestData(){
  653   double rsc, vsc, r, v, x, y, offset, tmp;
  654   double cmr[NDIM], cmv[NDIM], tmpv[NDIM], tmpv2[NDIM];
  655   body *p, *cp;
  656   int halfnbody, i, k;
  657 
  658   tnow = 0.0;
  659   rsc = 9 * PI / 16;
  660   vsc = sqrt(1.0 / rsc);
  661 
  662   tvecClr(cmr);
  663   tvecClr(cmv);
  664 
  665   srand(123);
  666 
  667 
  668 
  669   halfnbody = nbody / 2;
  670   if (nbody % 2 != 0) halfnbody++;
  671 
  672   for(k = 0; k < halfnbody; k++) {
  673     p = bodies->newbody();
  674     //p->setMass((double)(1.0/nbody));
  675     //p->setMass(0.1 + 0.9*((double)rand()/INT_MAX));
  676     //p->setMass(0.5);
  677 
  678     /* reject radii greater than 10 */
  679     do {
  680       tmp = xrand(0.0, 0.999);
  681       r = 1.0 / sqrt(pow(tmp, -2.0/3.0) - 1);
  682     } while (r > 9.0);
  683 
  684     pickshell(tmpv, rsc * r);
  685     p->PosvecAdd(tmpv);
  686     tvecAdd(cmr,tmpv);
  687     do{
  688       x = xrand(0.0, 1.0);
  689       y = xrand(0.0, 0.1);
  690     } while (y > x*x * pow(1 - x*x, 3.5));
  691     v = sqrt(2.0) * x / pow(1 + r*r, 0.25);
  692     pickshell(tmpv, vsc * v);
  693     p->VelvecAdd(tmpv);
  694     tvecAdd(cmv,tmpv);
  695   }
  696 
  697   offset = 4.0;
  698   for(k = halfnbody; k < nbody; k++) {
  699     p = bodies->newbody();
  700     //p->setMass((double)(1.0/nbody));
  701     //p->setMass(0.1 + 0.9*((double)rand()/INT_MAX));
  702     //p->setMass(0.5);
  703 
  704 
  705     cp = bodies->getbody(k-halfnbody);
  706     cp->PosvecStore(tmpv);
  707     tvecClr(tmpv2);
  708     for(i=0; i < NDIM; i++){
  709       tmpv2[i] = tmpv[i] + offset;
  710       tvecAdd(cmr,tmpv2);
  711     }
  712     p->PosvecAdd(tmpv2);
  713     cp->VelvecStore(tmpv);
  714     tvecClr(tmpv2);
  715     for(i=0; i < NDIM; i++){
  716       tmpv2[i] = tmpv[i];
  717       tvecAdd(cmv,tmpv2);
  718     }
  719     p->VelvecAdd(tmpv2);
  720   }
  721 
  722   tvecScale(cmr, 1.0/nbody);
  723   tvecScale(cmv, 1.0/nbody);
  724 
  725   for(k=0; k < nbody; k++) {
  726     p = bodies->getbody(k);
  727     //p->setMass(0.1 + 0.9*((double)rand()/INT_MAX));
  728     p->setMass(0.5);
  729 
  730     p->PosvecStore(tmpv);
  731     tvecSub(tmpv,cmr);
  732     p->PosvecLoad(tmpv);
  733     p->VelvecStore(tmpv);
  734     tvecSub(tmpv,cmv);
  735     p->VelvecLoad(tmpv);
  736   }
  737 }
  738 
  739 // ----------------------------------------------------------------------------
  740 void NBodySystem::pickshell(double *vec, double rad){
  741   int k;
  742   double rsq, rsc;
  743 
  744   do{
  745     for(k = 0; k < NDIM; k++)
  746       vec[k] = xrand(-1.0, 1.0);
  747     rsq = tvecDotProd(vec,vec);
  748   } while (rsq > 1.0);
  749   rsc = rad / sqrt(rsq);
  750   tvecScale(vec,rsc);
  751 }
  752 
  753 
  754 // ----------------------------------------------------------------------------
  755 void NBodySystem::startrun(parms *p){
  756   int seed;
  757   char *name, headbuf[128];
  758 
  759   p->loadInfile();
  760   name = p->getinfile();
  761   if(strlen(name) > 1) {
  762     sprintf(headbuf, "Hack code: input file %s\n", name);
  763     strcpy(p->getheadline(),headbuf);
  764   } else {
  765     nbody = p->getiparam("nbody");
  766     if(nbody < 1)
  767      cerr << "startrun: absurd nbody\n";
  768     seed = p->getiparam("seed");
  769     strcpy(p->getheadline(),"Hack code: Plummer model");
  770   }
  771 
  772   p->loadParms();
  773 
  774   // allocation of data structures filing is done later
  775   cells = new cellset((int)(p->getfcells()*p->getfleaves()*nbody));
  776   leaves = new leafset((int)(p->getfleaves()*nbody));
  777   bodies = new bodyset(nbody);
  778 
  779   tout = tnow + p->getdtout();
  780   pranset(seed);
  781 
  782   name = p->getinfile();
  783   if(strlen(name) > 1) {
  784     loadData(name);
  785   } else {
  786     loadTestData();
  787   }
  788 
  789 }
  790 
  791 // ----------------------------------------------------------------------------
  792 void NBodySystem::stepsystem(parms *p){
  793   int i;
  794   body *bptr;
  795   double incr, hts, ts, stop, tol, eps;
  796 
  797   int trackstart, trackend;
  798   int treebuildstart, treebuildend;
  799   int forcecalcstart, forcecalcend;
  800   int updatestart, updateend;
  801   int totalstep = 0; double now0 = tnow;
  802 
  803   incr = p->getdtime();
  804   hts = p->getdthf();
  805   ts = p->getdtime();
  806   stop = p->gettstop();
  807   tol = p->gettolsq();
  808   eps = p->getepssq();
  809 
  810   for (int isteps = 0; isteps < (int) ((stop-now0)/incr + 0.1) ; isteps++) {//while (tnow < stop + 0.1 * incr){
  811     get_time(&trackstart);
  812     leaves->LeafSetReset();
  813     cells->CellSetReset();
  814     BH_tree = cells->makecell(NULL);
  815 
  816     get_time(&treebuildstart);
  817 
  818     maketree();
  819     ((node*)BH_tree)->computecofm();
  820     get_time(&treebuildend);
  821     treebuildtime += treebuildend - treebuildstart;
  822 
  823     get_time(&forcecalcstart);
  824     computeForces(tol,hts,eps);
  825     get_time(&forcecalcend);
  826     forcecalctime += forcecalcend - forcecalcstart;
  827 
  828     /* advance bodies */
  829     get_time(&updatestart);
  830     Advance(hts,ts);
  831 /*
  832     for(i = 0; i < nbody; i++) {
  833        bptr = bodies->getbody(i);
  834        bptr->advance(hts, ts);
  835     }
  836 */
  837 
  838     get_time(&updateend);
  839     updatetime += updateend - updatestart;
  840     get_time(&trackend);
  841     tracktime += trackend - trackstart;
  842 
  843     /* compute bounding box */
  844     setbound();
  845 
  846     nstep++;
  847     tnow = tnow + incr;
  848     totalstep++;
  849 
  850     updateEnergy();
  851     outputEnergy();
  852 
  853     if (totalstep == 2) {
  854       reset_stat();
  855       clearTiming();
  856       clearStats();
  857       get_time(&computestart);
  858     }
  859   }
  860   get_time(&computeend);
  861   print_stat();
  862   computeend = computeend - computestart;
  863   computestart = 0;
  864 
  865 }
  866 
  867 
  868 // --- PARMS CODE ------------------------------------------------------------
  869 parms::parms()  {
  870   int i;
  871   // SET DEFAULT PARAMETER VALUES
  872   strcpy(headline,"");
  873   strcpy(infile,"");
  874   strcpy(outfile,"");
  875   dtime = 0.0;
  876   dtout = 0.0;
  877   tstop = 0.0;
  878   fcells = 0.0;
  879   fleaves = 0.0;
  880   tol = 0.0;
  881   tolsq = 0.0;
  882   eps = 0.0;
  883   epssq = 0.0;
  884   dthf = 0.0;
  885 
  886   for(i = 0; i < 12; i++)
  887     defaults[i] = (char*)malloc(MAX_NAME);
  888 
  889    /* file names for input/output                                 */
  890   strcpy(defaults[0],"in=");    /* snapshot of initial conditions  */
  891   strcpy(defaults[1],"out=");   /* stream of output snapshots      */
  892   /* params, used if no input specified, to make a Plummer Model   */
  893   strcpy(defaults[2],"nbody=1024");  /* number of particles to generate */
  894   strcpy(defaults[3],"seed=123");    /* random number generator seed    */
  895   /* params to control N-body integration                          */
  896   strcpy(defaults[4],"dtime=0.025"); /* integration time-step      */
  897   strcpy(defaults[5],"eps=0.05");    /* usual potential softening  */
  898   strcpy(defaults[6],"tol=1.0");     /* cell subdivision tolerence */
  899   strcpy(defaults[7],"fcells=2.0");  /* cell allocation parameter  */
  900   strcpy(defaults[8],"fleaves=0.5"); /* leaf allocation parameter  */
  901   strcpy(defaults[9],"tstop=0.075"); /* time to stop integration   */
  902   strcpy(defaults[10],"dtout=0.25"); /* data-output interval       */
  903   strcpy(defaults[11],"");
  904 }
  905 
  906 void parms::setparms(char **fnames, double *args){
  907   strcpy(headline,fnames[0]);
  908   strcpy(infile,fnames[1]);
  909   strcpy(outfile,fnames[2]);
  910   dtime = args[0];
  911   dtout = args[1];
  912   tstop = args[2];
  913   fcells = args[3];
  914   fleaves = args[4];
  915   tol = args[5];
  916   tolsq = args[6];
  917   eps = args[7];
  918   epssq = args[8];
  919   dthf = args[9];
  920 }
  921 
  922 
  923 int parms::getiparam(char *name){
  924   getparam(name);
  925   return (atoi(buf));
  926 }
  927 
  928 // ---------------------------------------------------------------------------
  929 double parms::getdparam(char *name){
  930    getparam(name);
  931    return (atof(buf));
  932 }
  933 
  934 void parms::getparam(char *name) {
  935   int i, leng;
  936   char *temp;
  937 
  938   buf[0] = '\0';
  939   if(defaults == NULL){
  940     cerr << "getparam: called before initparam\n";
  941   }
  942   i = scanbind(defaults, name);
  943   if(i < 0){
  944      cerr << "getparam: " << name << " unknown\n";
  945   }
  946   temp = extrvalue(defaults[i]);
  947   fgets(buf, 128, stdin);
  948   leng = strlen(buf) + 1;
  949   if(leng <= 1)
  950     strcpy(buf,temp);
  951 }
  952 
  953 // ---------------------------------------------------------------------------
  954 void parms::loadInfile(){
  955   getparam("in");
  956   strcpy(infile,buf);
  957 }
  958 
  959 // ---------------------------------------------------------------------------
  960 void parms::loadParms(){
  961   getparam("out");
  962   strcpy(outfile,buf);
  963   dtime = getdparam("dtime");
  964   dthf = 0.5 * dtime;
  965   eps = getdparam("eps");
  966   epssq = eps*eps;
  967   tol = getdparam("tol");
  968   tolsq = tol*tol;
  969   fcells = getdparam("fcells");
  970   fleaves = getdparam("fleaves");
  971   tstop = getdparam("tstop");
  972   dtout = getdparam("dtout");
  973 }
  974 
  975 
  976 
  977 // --- CELLSET CODE ------------------------------------------------------------
  978 cell *cellset::newcell(){
  979   cell *c;
  980   if(numcell == maxcell){
  981     cerr << "newcell: More than " << maxcell << " cells; increase fcells\n";
  982     return NULL;
  983   }
  984   c = &celltab[numcell];
  985   numcell++;
  986   c->CellClear();
  987   return c;
  988 }
  989 // ---------------------------------------------------------------------------
  990 cell *cellset::makecell(node *parent) {
  991   cell *c;
  992   c = newcell();
  993   c->setParent(parent);
  994   if(c == NULL) return NULL;
  995   c->adjLevel(parent);
  996   return c;
  997 }
  998 
  999 
 1000 // --- LEAFSET CODE ------------------------------------------------------------
 1001 leaf *leafset::newleaf() {
 1002   leaf* p;
 1003   if(numleaf == maxleaf){
 1004     cerr << "newleaf: More than " << maxleaf << " leaves; increase fleaves\n";
 1005     return NULL;
 1006   }
 1007   p = &leaftab[numleaf];
 1008   numleaf++;
 1009   p->LeafClear();
 1010   return p;
 1011 }
 1012 
 1013 // ---------------------------------------------------------------------------
 1014 leaf *leafset::makeleaf(node *parent){
 1015   leaf* l;
 1016   l = newleaf();
 1017   if(l == NULL) return NULL;
 1018   l->setParent(parent);
 1019   l->adjLevel(parent);
 1020   return (l);
 1021 }
 1022 
 1023 
 1024 // --- BODYSET CODE ------------------------------------------------------------
 1025 body* bodyset::newbody(){
 1026   body *p;
 1027   if(numbody == maxbody){
 1028     cerr << "newbody: More than " << maxbody << "body; increase nbody\n";
 1029     return NULL;
 1030   }
 1031   p = &bodytab[numbody];
 1032   numbody++;
 1033   return p;
 1034 }
 1035 
 1036 // --- LEAF CODE ------------------------------------------------------------
 1037 void leaf::LeafAddBody(body *bp){                  // Assumes enough capacity
 1038   bp->setParent(this);
 1039   bp->setLevel(Level());
 1040   bp->setChildNum(num_bodies);
 1041   bodyp[num_bodies] = bp;
 1042   num_bodies++;
 1043 }
 1044 
 1045 void leaf::Subdivide(cell *parent, int idx){
 1046   cell *c;
 1047   leaf *le;
 1048   double tmpv[NDIM];
 1049   int i, index, xp[NDIM], num, Lev;
 1050   body *loctab[LEAFMAXBODIES], *p;
 1051 
 1052   num = num_bodies;
 1053   for(i=0; i < num; i++) {
 1054     loctab[i] = bodyp[i];
 1055     bodyp[i] = NULL;
 1056   }
 1057   num_bodies = 0;
 1058   c = Nbody.getcells()->makecell(parent);
 1059   c->setChildNum(child_num);
 1060   parent->setSubp(idx,c);
 1061   Lev = level;
 1062 
 1063   p = loctab[0];
 1064   p->PosvecStore(tmpv);
 1065   Nbody.intcoord(xp,tmpv);
 1066   index = subindex(xp, Lev);
 1067   c->setSubp(index,this);
 1068   child_num = index;
 1069   parent = c;
 1070   level = (Lev >> 1);
 1071   LeafAddBody(p);
 1072 
 1073   for(i=1; i < num; i++) {
 1074     p = loctab[i];
 1075     p->PosvecStore(tmpv);
 1076     Nbody.intcoord(xp,tmpv);
 1077     index = subindex(xp, Lev);
 1078     le = (leaf*)(c->Subp(index));
 1079     if(le == NULL) {
 1080       le = Nbody.getleaves()->makeleaf(c);
 1081       le->setChildNum(index);
 1082       c->setSubp(index,le);
 1083     }
 1084     le->LeafAddBody(p);
 1085   }
 1086 }
 1087 
 1088 
 1089 // --- CELL CODE ------------------------------------------------------------
 1090 void cell::CellAddBody(body *p, int *coords) {
 1091   int Lev, kidIndex;
 1092   node *ptr;
 1093   leaf *le;
 1094 
 1095   kidIndex = subindex(coords,level);
 1096   ptr = subp[kidIndex];
 1097   Lev = (level >> 1);
 1098   if(Lev != 0){
 1099     if(ptr == NULL) {
 1100       le = Nbody.getleaves()->makeleaf(this);
 1101       le->setChildNum(kidIndex);
 1102       subp[kidIndex] = le;
 1103       le->LeafAddBody(p);
 1104     } else {
 1105       if(ptr->Type() == LEAF){
 1106         le = (leaf*)ptr;
 1107         if(le->full()){
 1108           le->Subdivide(this,kidIndex);
 1109           CellAddBody(p,coords);
 1110         } else {
 1111           le->LeafAddBody(p);
 1112         }
 1113       } else {
 1114         ((cell*)ptr)->CellAddBody(p,coords);
 1115       }
 1116     }
 1117   } else {
 1118     cerr << " *** Error: Not enough levels in tree...(CellAddBody)\n";
 1119   }
 1120 }
 1121 
 1122 // ----------------------------------------------------------------------------
 1123 int node::subindex(int *x, int l){
 1124   int i, k, yes;
 1125 
 1126   i = 0;
 1127   yes = 0;
 1128   for(k = 0; k < NDIM; k++) {
 1129     if(((x[k] & l) && !yes) || (!(x[k] & l) && yes)) {
 1130       i += NSUB >> (k + 1);
 1131       yes = 1;
 1132     }
 1133     else yes = 0;
 1134   }
 1135   return (i);
 1136 }
 1137 
 1138 void NBodySystem::setbound(){
 1139   int i;
 1140   double tmp1[NDIM], tmp2[NDIM], side;
 1141   body *p;
 1142 
 1143   mincorner->vecSetUnit();
 1144   mincorner->vecScale(1E99);
 1145   maxcorner->vecSetUnit();
 1146   maxcorner->vecScale(-1E99);
 1147   for(i=0; i < nbody; i++) {
 1148     p = bodies->getbody(i);
 1149     p->PosvecStore(tmp1);
 1150     mincorner->vecMin(tmp1);
 1151     maxcorner->vecMax(tmp1);
 1152   }
 1153 
 1154   side = 0.0;
 1155   maxcorner->vecStore(tmp1);
 1156   mincorner->vecStore(tmp2);
 1157   tvecSub(tmp1,tmp2);
 1158   for(i=0; i < NDIM; i++)
 1159     if(side < tmp1[i])
 1160       side = tmp1[i];
 1161 
 1162   for(i=0; i < NDIM; i++)
 1163     tmp1[i] = -side/100000;
 1164   mincorner->vecAdd(tmp1);
 1165   size = 1.00002*side;
 1166 }
 1167 
 1168 
 1169 // ----------------------------------------------------------------------------
 1170 // MAKETREE: initialize tree structure for force calculation.
 1171 // ----------------------------------------------------------------------------
 1172 void NBodySystem::maketree(){
 1173   int i, xp[NDIM];
 1174   double v[NDIM];
 1175   body *p;
 1176 
 1177   for(i=0; i < nbody; i++){
 1178     p = bodies->getbody(i);
 1179     p->PosvecStore(v);
 1180     intcoord(xp,v);
 1181     BH_tree->CellAddBody(p,xp);
 1182   }
 1183 }
 1184 
 1185 
 1186 void node::computecofm(){  // <--- Parallelizable...
 1187   int i;
 1188   node *nptr;
 1189   body *p;
 1190   leaf *le;
 1191   cell *ce;
 1192   double m;
 1193   double tmpv[NDIM];
 1194 
 1195   if(type == LEAF){
 1196     le = (leaf*)this;
 1197     for(i=0; i < le->NumBody(); i++) {
 1198       p = le->Bodyp(i);
 1199       m = p->Mass();
 1200       mass += m;
 1201       p->PosvecStore(tmpv);
 1202       tvecScale(tmpv,m);
 1203       pos.vecAdd(tmpv);
 1204     }
 1205     pos.vecStore(tmpv);
 1206     tvecScale(tmpv,1.0/mass);
 1207     pos.vecLoad(tmpv);
 1208   } else {
 1209     ce = (cell*)this;
 1210     for(i=0; i < NSUB; i++){
 1211       nptr = ce->Subp(i);
 1212       if(nptr != NULL){
 1213         nptr->computecofm();
 1214         m = nptr->Mass();
 1215         mass += m;
 1216         nptr->PosvecStore(tmpv);
 1217         tvecScale(tmpv,m);
 1218         pos.vecAdd(tmpv);
 1219       }
 1220     }
 1221     pos.vecStore(tmpv);
 1222     tvecScale(tmpv,1.0/mass);
 1223     pos.vecLoad(tmpv);
 1224   }
 1225 }
 1226 
 1227 double body::subdiv(node *p, double *res) {
 1228   int i;
 1229   double drsq, d;
 1230 
 1231   drsq = 0.0;
 1232   for (i = 0; i < NDIM; i++) {
 1233     d = p->pos.vecField(i)-pos.vecField(i);
 1234     drsq += d*d;
 1235     res[i] = d;
 1236   }
 1237   return(drsq);
 1238 }
 1239 
 1240 inline double body::computeInter(node *p, double epsSq, double *res){
 1241   int i;
 1242   double drabs, inc, mor3;
 1243   double drsq, d;
 1244 
 1245   drsq = epsSq;
 1246   for (i = 0; i < NDIM; i++) {
 1247     d = p->pos.vecField(i)-pos.vecField(i);
 1248     drsq += d*d;
 1249     res[i] = d;
 1250   }
 1251   drabs = sqrt(drsq);
 1252   inc = p->Mass() / drabs;
 1253   mor3 = inc / drsq;
 1254   tvecScale(res,mor3);
 1255   return inc;
 1256 }
 1257 
 1258 void body::updatePhi(double inc) {
 1259   phi -= inc;
 1260 }
 1261 
 1262 void body::fastgravsub(node *p, double *ai, double pdrsq, double epssq) {
 1263   double drsq;
 1264   double drabs, inc, mor3;
 1265   double tmpv[NDIM];
 1266 
 1267   drsq = pdrsq + epssq;
 1268   drabs = sqrt(drsq);
 1269   inc = p->Mass() / drabs;
 1270 //  phi -= inc;
 1271   updatePhi(inc);
 1272   mor3 = inc / drsq;
 1273   tvecLoad(tmpv,ai);
 1274   tvecScale(tmpv,mor3);
 1275   acc.vecAdd(tmpv);
 1276 }
 1277 
 1278 void body::gravsub(node *p, double epsSq){
 1279   double phii;
 1280   double ai[NDIM];
 1281 
 1282   phii = computeInter(p,epsSq,ai);
 1283 //  phi -= phii;
 1284   updatePhi(phii);
 1285   acc.vecAdd(ai);
 1286 }
 1287 
 1288 inline void body::openCell(node *n, double tolsq, double dsq, double epssq){
 1289   int i;
 1290   node* nn;
 1291 
 1292   for(i=0; i < NSUB; i++) {
 1293     nn = ((cell*)n)->Subp(i);
 1294     if(nn != NULL){
 1295       walksub(nn,tolsq,(dsq/4.0),epssq);
 1296     }
 1297   }
 1298 }
 1299 
 1300 inline void body::openLeaf(node *n,double epssq, int iter){
 1301   int i;
 1302   body *bptr;
 1303 
 1304   for(i=0; i < iter; i++){
 1305     bptr = ((leaf*)n)->Bodyp(i);
 1306     if(bptr != this){
 1307       gravsub((node*)bptr, epssq);
 1308     }
 1309   }
 1310 }
 1311 
 1312 void body::walksub(node *n, double tolsq, double dsq, double epssq){
 1313   double drsq;
 1314   double ai[NDIM];
 1315   drsq = subdiv(n,ai);
 1316   if((tolsq * drsq) < dsq){
 1317     if((n->Type()) == CELL){
 1318       openCell(n,tolsq,dsq,epssq);
 1319     } else {
 1320       openLeaf(n,epssq,((leaf*)n)->NumBody());
 1321     }
 1322   } else {
 1323     fastgravsub(n, ai, drsq, epssq);
 1324   }
 1325 }
 1326 
 1327 void body::advance(double hts, double ts){
 1328   double dvel[NDIM], vel1[NDIM], dpos[NDIM];
 1329 
 1330   acc.vecStore(dvel);
 1331   vel.vecStore(vel1);
 1332   tvecScale(dvel,hts);
 1333   tvecAdd(vel1,dvel);
 1334   tvecLoad(dpos,vel1);
 1335   tvecScale(dpos,ts);
 1336   tvecScale(dvel,2.0);
 1337 
 1338   vel.vecAdd(dvel);
 1339   pos.vecAdd(dpos);
 1340 }
 1341 
 1342 void NBodySystem::Advance(double hts, double ts){
 1343   int i;
 1344   body *p;
 1345 
 1346   for(i=0; i < nbody; i++) {
 1347     p = bodies->getbody(i);
 1348     p->advance(hts,ts);
 1349   }
 1350 }
 1351 
 1352 void body::swapAcc() {
 1353   double tmpv[NDIM];
 1354 
 1355   phi = 0.0;
 1356   acc.vecStore(tmpv);
 1357   tmp.vecLoad(tmpv);
 1358   acc.vecClr();
 1359 }
 1360 
 1361 void body::startVel(double f){
 1362    double tmpv1[NDIM], tmpv2[NDIM];
 1363    acc.vecStore(tmpv1);
 1364    tmp.vecStore(tmpv2);
 1365    tvecSub(tmpv1,tmpv2);
 1366    tvecScale(tmpv1,f);
 1367    vel.vecAdd(tmpv1);
 1368 }
 1369 
 1370 void NBodySystem::SwapAccs(){
 1371   int i;
 1372   body *p;
 1373 
 1374   for(i=0; i < nbody; i++) {
 1375     p = bodies->getbody(i);
 1376     p->swapAcc();
 1377   }
 1378 }
 1379 
 1380 void NBodySystem::ComputeAccels(double tol, double eps){
 1381   int i;
 1382   body *p;
 1383   int nbody_local = nbody;
 1384 
 1385    #pragma omp parallel for private (p) schedule(static)
 1386   for(i=0; i < nbody_local; i++) {
 1387     p = bodies->getbody(i);
 1388     p->walksub(BH_tree, tol, size*size, eps);
 1389   }
 1390 }
 1391 
 1392 void NBodySystem::StartVels(double hts){
 1393   int i;
 1394   body *p;
 1395 
 1396   for(i=0; i < nbody; i++) {
 1397     p = bodies->getbody(i);
 1398     p->startVel(hts);
 1399   }
 1400 }
 1401 
 1402 void NBodySystem::computeForces(double tol, double hts, double eps){
 1403   SwapAccs();
 1404   ComputeAccels(tol,eps);
 1405   if(nstep > 0){
 1406     StartVels(hts);
 1407   }
 1408 }
 1409 
 1410 // ----------------------------------------------------------------------------
 1411 int main(int argc, char **argv){
 1412   Nbody.init();
 1413   Nbody.startrun(&simparms);
 1414   Nbody.initOutput(&simparms);
 1415   Nbody.setbound();
 1416   Nbody.clearTiming();
 1417   Nbody.clearStats();
 1418 
 1419   Nbody.stepsystem(&simparms);
 1420 
 1421   //  bodyset * bset = Nbody.getbodies();
 1422   //for (int i = 0; i < bset->size(); i++)
 1423   //  bset->
 1424 
 1425   Nbody.outputTiming();
 1426   //Nbody.outputStats();
 1427   Nbody.output(&simparms);
 1428   Nbody.outputEnergy();
 1429   //Nbody.dump();
 1430 }
 1431