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   for(k = 0; k < nbody; k++) {
  501     p = bodies->getbody(k);
  502     matrix m1;
  503     double v0[3], v1[3], v2[3], v3[3], v4[3];
  504     m = p->Mass();
  505     mtot.addval(m);
  506     p->VelvecStore(v1);
  507     velsq = tvecDotProd(v1,v1);
  508     v4[1] = 0.5 * m * velsq;
  509     v4[2] = 0.5 * m * p->Phi();
  510     v4[0] = v4[1] + v4[2];
  511     etot->vecAdd(v4);
  512     p->PosvecStore(v1);
  513     p->VelvecStore(v2);
  514     p->AccvecStore(v3);
  515     tvecLoad(v0,v2);
  516     tvecScale(v0,0.5*m);
  517     m1.matOutProd(v0,v2);
  518     keten->matAdd(&m1);
  519     tvecLoad(v0,v1);
  520     tvecScale(v0,m);
  521     m1.matOutProd(v0,v3);
  522     peten->matAdd(&m1);
  523     tvecLoad(v0,v1);
  524     tvecScale(v0,m);
  525     cmphasePos->vecAdd(v0);
  526     tvecLoad(v0,v2);
  527     tvecScale(v0,m);
  528     cmphaseVel->vecAdd(v0);
  529     tvecCrossProd(v0,v1,v2);
  530     tvecScale(v0,m);
  531     amvec->vecAdd(v0);
  532   }
  533 
  534    if (mtot.readval() != 0.0){
  535       cmphasePos->vecScale(1.0/mtot.readval());
  536       cmphaseVel->vecScale(1.0/mtot.readval());
  537    }
  538 }
  539 
  540 
  541 // ----------------------------------------------------------------------------
  542 void NBodySystem::outputEnergy(){
  543   cout << "\n\nENERGY TOTALS\n";
  544   cout << " MASS:        "<< mtot.readval() << "\n";
  545   cout << " ENERGY:     "; etot->vecPrint(); cout << "\n";
  546   cout << " KIN ENERGY: "; keten->matPrint(); cout << "\n";
  547   cout << " POT ENERGY: "; peten->matPrint(); cout << "\n";
  548   cout << " AGGREGATE CM POS: "; cmphasePos->vecPrint(); cout << "\n";
  549   cout << " AGGREGATE CM VEL: "; cmphaseVel->vecPrint(); cout << "\n";
  550   cout << " ANGULAR MOMENTUM: "; amvec->vecPrint(); cout << "\n\n";
  551 }
  552 
  553 // ----------------------------------------------------------------------------
  554 void NBodySystem::output(parms *sim) {
  555   if( (tout - 0.01 * sim->getdtime()) <= tnow)
  556     tout += sim->getdtout();
  557 }
  558 
  559 // ----------------------------------------------------------------------------
  560 void NBodySystem::clearTiming(){
  561    tracktime = 0;
  562    updatetime = 0;
  563    partitiontime = 0;
  564    treebuildtime = 0;
  565    forcecalctime = 0;
  566 }
  567 
  568 // ----------------------------------------------------------------------------
  569 void NBodySystem::initOutput(parms *p){
  570   cout << "\n\t\tOutput: " << p->getheadline() << "\n\n";
  571   cout << "nbody   dtime   eps   tol   dtout   tstop   fcells\n";
  572   cout << nbody << "     " << p->getdtime() << "   " << p->geteps()
  573       << "   " << p->gettol() << "   " << p->getdtout() << "    "
  574       << p->gettstop() << "   " << p->getfcells() << "\n\n";
  575 }
  576 
  577 // ----------------------------------------------------------------------------
  578 void NBodySystem::outputTiming(){
  579   int ticks;
  580   double dticks;
  581 
  582   get_ticks_time(&ticks);
  583   dticks = ticks;
  584   cout << "COMPUTESTART  = " <<
  585     (computestart / dticks) << "\n";
  586   cout << "COMPUTEEND    = " <<
  587     (computeend / dticks) << "\n";
  588   cout << "COMPUTETIME   = " <<
  589     ((computeend - computestart) / dticks) << "\n";
  590   cout << "TRACKTIME     = " <<
  591     (tracktime / dticks) << "\n";
  592   cout << "PARTITIONTIME = " <<
  593     (partitiontime / dticks) << "\t"
  594         << ((float)partitiontime)/tracktime << "\n";
  595   cout << "TREEBUILDTIME = " <<
  596     (treebuildtime / dticks) << "\t"
  597         << ((float)treebuildtime)/tracktime << "\n";
  598   cout << "FORCECALCTIME = " <<
  599     (forcecalctime / dticks) << "\t"
  600         << ((float)forcecalctime)/tracktime << "\n";
  601   cout << "UPDATEPOSTIME = " <<
  602     (updatetime / dticks) << "\t"
  603         << ((float)updatetime)/tracktime << "\n";
  604   cout << "RESTTIME      = " <<
  605     ((tracktime - partitiontime -
  606      updatetime - treebuildtime - forcecalctime)
  607     / dticks) << "\t"
  608         << ((float)(tracktime-partitiontime- updatetime - treebuildtime-
  609         forcecalctime))/tracktime<<"\n";
  610   cout.flush();
  611 }
  612 
  613 // ----------------------------------------------------------------------------
  614 void NBodySystem::loadData(char *name){
  615   int i, ndim;
  616   body *p;
  617   double m, tmp[3];
  618   ifstream instr(name);
  619 
  620   cerr << "Reading input file : " << name << '\n';
  621   cerr.flush();
  622   if (!instr)
  623      cerr << "loadData: cannot find file " << name << '\n';
  624   instr >> nbody;
  625   if (nbody < 1)
  626      cerr << "loadData: nbody = " << nbody << " is absurd\n";
  627   instr >> ndim;
  628   if(ndim != NDIM)
  629     cerr << "inputdata: NDIM = " << NDIM << "ndim = " << ndim << " is absurd\n";
  630   instr >> tnow;
  631 
  632   for(i = 0; i < nbody; i++){
  633     p = bodies->getbody(i);
  634     instr >> m;
  635     p->setMass(m);
  636   }
  637   for(i = 0; i < nbody; i++){
  638     p = bodies->getbody(i);
  639     in_vector(instr,tmp);
  640     p->PosvecLoad(tmp);
  641   }
  642   for(i = 0; i < nbody; i++){
  643     p = bodies->getbody(i);
  644     in_vector(instr,tmp);
  645     p->VelvecLoad(tmp);
  646   }
  647   instr.close();
  648 }
  649 
  650 // ----------------------------------------------------------------------------
  651 void NBodySystem::loadTestData(){
  652   double rsc, vsc, r, v, x, y, offset, tmp;
  653   double cmr[NDIM], cmv[NDIM], tmpv[NDIM], tmpv2[NDIM];
  654   body *p, *cp;
  655   int halfnbody, i, k;
  656 
  657   tnow = 0.0;
  658   rsc = 9 * PI / 16;
  659   vsc = sqrt(1.0 / rsc);
  660 
  661   tvecClr(cmr);
  662   tvecClr(cmv);
  663 
  664   srand(123);
  665 
  666 
  667 
  668   halfnbody = nbody / 2;
  669   if (nbody % 2 != 0) halfnbody++;
  670 
  671   for(k = 0; k < halfnbody; k++) {
  672     p = bodies->newbody();
  673     //p->setMass((double)(1.0/nbody));
  674     //p->setMass(0.1 + 0.9*((double)rand()/INT_MAX));
  675     //p->setMass(0.5);
  676 
  677     /* reject radii greater than 10 */
  678     do {
  679       tmp = xrand(0.0, 0.999);
  680       r = 1.0 / sqrt(pow(tmp, -2.0/3.0) - 1);
  681     } while (r > 9.0);
  682 
  683     pickshell(tmpv, rsc * r);
  684     p->PosvecAdd(tmpv);
  685     tvecAdd(cmr,tmpv);
  686     do{
  687       x = xrand(0.0, 1.0);
  688       y = xrand(0.0, 0.1);
  689     } while (y > x*x * pow(1 - x*x, 3.5));
  690     v = sqrt(2.0) * x / pow(1 + r*r, 0.25);
  691     pickshell(tmpv, vsc * v);
  692     p->VelvecAdd(tmpv);
  693     tvecAdd(cmv,tmpv);
  694   }
  695 
  696   offset = 4.0;
  697   for(k = halfnbody; k < nbody; k++) {
  698     p = bodies->newbody();
  699     //p->setMass((double)(1.0/nbody));
  700     //p->setMass(0.1 + 0.9*((double)rand()/INT_MAX));
  701     //p->setMass(0.5);
  702 
  703 
  704     cp = bodies->getbody(k-halfnbody);
  705     cp->PosvecStore(tmpv);
  706     tvecClr(tmpv2);
  707     for(i=0; i < NDIM; i++){
  708       tmpv2[i] = tmpv[i] + offset;
  709       tvecAdd(cmr,tmpv2);
  710     }
  711     p->PosvecAdd(tmpv2);
  712     cp->VelvecStore(tmpv);
  713     tvecClr(tmpv2);
  714     for(i=0; i < NDIM; i++){
  715       tmpv2[i] = tmpv[i];
  716       tvecAdd(cmv,tmpv2);
  717     }
  718     p->VelvecAdd(tmpv2);
  719   }
  720 
  721   tvecScale(cmr, 1.0/nbody);
  722   tvecScale(cmv, 1.0/nbody);
  723 
  724   for(k=0; k < nbody; k++) {
  725     p = bodies->getbody(k);
  726     //p->setMass(0.1 + 0.9*((double)rand()/INT_MAX));
  727     p->setMass(0.5);
  728 
  729     p->PosvecStore(tmpv);
  730     tvecSub(tmpv,cmr);
  731     p->PosvecLoad(tmpv);
  732     p->VelvecStore(tmpv);
  733     tvecSub(tmpv,cmv);
  734     p->VelvecLoad(tmpv);
  735   }
  736 }
  737 
  738 // ----------------------------------------------------------------------------
  739 void NBodySystem::pickshell(double *vec, double rad){
  740   int k;
  741   double rsq, rsc;
  742 
  743   do{
  744     for(k = 0; k < NDIM; k++)
  745       vec[k] = xrand(-1.0, 1.0);
  746     rsq = tvecDotProd(vec,vec);
  747   } while (rsq > 1.0);
  748   rsc = rad / sqrt(rsq);
  749   tvecScale(vec,rsc);
  750 }
  751 
  752 
  753 // ----------------------------------------------------------------------------
  754 void NBodySystem::startrun(parms *p){
  755   int seed;
  756   char *name, headbuf[128];
  757 
  758   p->loadInfile();
  759   name = p->getinfile();
  760   if(strlen(name) > 1) {
  761     sprintf(headbuf, "Hack code: input file %s\n", name);
  762     strcpy(p->getheadline(),headbuf);
  763   } else {
  764     nbody = p->getiparam("nbody");
  765     if(nbody < 1)
  766      cerr << "startrun: absurd nbody\n";
  767     seed = p->getiparam("seed");
  768     strcpy(p->getheadline(),"Hack code: Plummer model");
  769   }
  770 
  771   p->loadParms();
  772 
  773   // allocation of data structures filing is done later
  774   cells = new cellset((int)(p->getfcells()*p->getfleaves()*nbody));
  775   leaves = new leafset((int)(p->getfleaves()*nbody));
  776   bodies = new bodyset(nbody);
  777 
  778   tout = tnow + p->getdtout();
  779   pranset(seed);
  780 
  781   name = p->getinfile();
  782   if(strlen(name) > 1) {
  783     loadData(name);
  784   } else {
  785     loadTestData();
  786   }
  787 
  788 }
  789 
  790 // ----------------------------------------------------------------------------
  791 void NBodySystem::stepsystem(parms *p){
  792   int i;
  793   body *bptr;
  794   double incr, hts, ts, stop, tol, eps;
  795 
  796   int trackstart, trackend;
  797   int treebuildstart, treebuildend;
  798   int forcecalcstart, forcecalcend;
  799   int updatestart, updateend;
  800   int totalstep = 0; double now0 = tnow;
  801 
  802   incr = p->getdtime();
  803   hts = p->getdthf();
  804   ts = p->getdtime();
  805   stop = p->gettstop();
  806   tol = p->gettolsq();
  807   eps = p->getepssq();
  808 
  809   for (int isteps = 0; isteps < (int) ((stop-now0)/incr + 0.1) ; isteps++) {//while (tnow < stop + 0.1 * incr){
  810     get_time(&trackstart);
  811     leaves->LeafSetReset();
  812     cells->CellSetReset();
  813     BH_tree = cells->makecell(NULL);
  814 
  815     get_time(&treebuildstart);
  816 
  817     maketree();
  818     ((node*)BH_tree)->computecofm();
  819     get_time(&treebuildend);
  820     treebuildtime += treebuildend - treebuildstart;
  821 
  822     get_time(&forcecalcstart);
  823     computeForces(tol,hts,eps);
  824     get_time(&forcecalcend);
  825     forcecalctime += forcecalcend - forcecalcstart;
  826 
  827     /* advance bodies */
  828     get_time(&updatestart);
  829     Advance(hts,ts);
  830 /*
  831     for(i = 0; i < nbody; i++) {
  832        bptr = bodies->getbody(i);
  833        bptr->advance(hts, ts);
  834     }
  835 */
  836 
  837     get_time(&updateend);
  838     updatetime += updateend - updatestart;
  839     get_time(&trackend);
  840     tracktime += trackend - trackstart;
  841 
  842     /* compute bounding box */
  843     setbound();
  844 
  845     nstep++;
  846     tnow = tnow + incr;
  847     totalstep++;
  848 
  849     updateEnergy();
  850     outputEnergy();
  851 
  852     if (totalstep == 2) {
  853       reset_stat();
  854       clearTiming();
  855       clearStats();
  856       get_time(&computestart);
  857     }
  858   }
  859   get_time(&computeend);
  860   print_stat();
  861   computeend = computeend - computestart;
  862   computestart = 0;
  863 
  864 }
  865 
  866 
  867 // --- PARMS CODE ------------------------------------------------------------
  868 parms::parms()  {
  869   int i;
  870   // SET DEFAULT PARAMETER VALUES
  871   strcpy(headline,"");
  872   strcpy(infile,"");
  873   strcpy(outfile,"");
  874   dtime = 0.0;
  875   dtout = 0.0;
  876   tstop = 0.0;
  877   fcells = 0.0;
  878   fleaves = 0.0;
  879   tol = 0.0;
  880   tolsq = 0.0;
  881   eps = 0.0;
  882   epssq = 0.0;
  883   dthf = 0.0;
  884 
  885   for(i = 0; i < 12; i++)
  886     defaults[i] = (char*)malloc(MAX_NAME);
  887 
  888    /* file names for input/output                                 */
  889   strcpy(defaults[0],"in=");    /* snapshot of initial conditions  */
  890   strcpy(defaults[1],"out=");   /* stream of output snapshots      */
  891   /* params, used if no input specified, to make a Plummer Model   */
  892   strcpy(defaults[2],"nbody=1024");  /* number of particles to generate */
  893   strcpy(defaults[3],"seed=123");    /* random number generator seed    */
  894   /* params to control N-body integration                          */
  895   strcpy(defaults[4],"dtime=0.025"); /* integration time-step      */
  896   strcpy(defaults[5],"eps=0.05");    /* usual potential softening  */
  897   strcpy(defaults[6],"tol=1.0");     /* cell subdivision tolerence */
  898   strcpy(defaults[7],"fcells=2.0");  /* cell allocation parameter  */
  899   strcpy(defaults[8],"fleaves=0.5"); /* leaf allocation parameter  */
  900   strcpy(defaults[9],"tstop=0.075"); /* time to stop integration   */
  901   strcpy(defaults[10],"dtout=0.25"); /* data-output interval       */
  902   strcpy(defaults[11],"");
  903 }
  904 
  905 void parms::setparms(char **fnames, double *args){
  906   strcpy(headline,fnames[0]);
  907   strcpy(infile,fnames[1]);
  908   strcpy(outfile,fnames[2]);
  909   dtime = args[0];
  910   dtout = args[1];
  911   tstop = args[2];
  912   fcells = args[3];
  913   fleaves = args[4];
  914   tol = args[5];
  915   tolsq = args[6];
  916   eps = args[7];
  917   epssq = args[8];
  918   dthf = args[9];
  919 }
  920 
  921 
  922 int parms::getiparam(char *name){
  923   getparam(name);
  924   return (atoi(buf));
  925 }
  926 
  927 // ---------------------------------------------------------------------------
  928 double parms::getdparam(char *name){
  929    getparam(name);
  930    return (atof(buf));
  931 }
  932 
  933 void parms::getparam(char *name) {
  934   int i, leng;
  935   char *temp;
  936 
  937   buf[0] = '\0';
  938   if(defaults == NULL){
  939     cerr << "getparam: called before initparam\n";
  940   }
  941   i = scanbind(defaults, name);
  942   if(i < 0){
  943      cerr << "getparam: " << name << " unknown\n";
  944   }
  945   temp = extrvalue(defaults[i]);
  946   fgets(buf, 128, stdin);
  947   leng = strlen(buf) + 1;
  948   if(leng <= 1)
  949     strcpy(buf,temp);
  950 }
  951 
  952 // ---------------------------------------------------------------------------
  953 void parms::loadInfile(){
  954   getparam("in");
  955   strcpy(infile,buf);
  956 }
  957 
  958 // ---------------------------------------------------------------------------
  959 void parms::loadParms(){
  960   getparam("out");
  961   strcpy(outfile,buf);
  962   dtime = getdparam("dtime");
  963   dthf = 0.5 * dtime;
  964   eps = getdparam("eps");
  965   epssq = eps*eps;
  966   tol = getdparam("tol");
  967   tolsq = tol*tol;
  968   fcells = getdparam("fcells");
  969   fleaves = getdparam("fleaves");
  970   tstop = getdparam("tstop");
  971   dtout = getdparam("dtout");
  972 }
  973 
  974 
  975 
  976 // --- CELLSET CODE ------------------------------------------------------------
  977 cell *cellset::newcell(){
  978   cell *c;
  979   if(numcell == maxcell){
  980     cerr << "newcell: More than " << maxcell << " cells; increase fcells\n";
  981     return NULL;
  982   }
  983   c = &celltab[numcell];
  984   numcell++;
  985   c->CellClear();
  986   return c;
  987 }
  988 // ---------------------------------------------------------------------------
  989 cell *cellset::makecell(node *parent) {
  990   cell *c;
  991   c = newcell();
  992   c->setParent(parent);
  993   if(c == NULL) return NULL;
  994   c->adjLevel(parent);
  995   return c;
  996 }
  997 
  998 
  999 // --- LEAFSET CODE ------------------------------------------------------------
 1000 leaf *leafset::newleaf() {
 1001   leaf* p;
 1002   if(numleaf == maxleaf){
 1003     cerr << "newleaf: More than " << maxleaf << " leaves; increase fleaves\n";
 1004     return NULL;
 1005   }
 1006   p = &leaftab[numleaf];
 1007   numleaf++;
 1008   p->LeafClear();
 1009   return p;
 1010 }
 1011 
 1012 // ---------------------------------------------------------------------------
 1013 leaf *leafset::makeleaf(node *parent){
 1014   leaf* l;
 1015   l = newleaf();
 1016   if(l == NULL) return NULL;
 1017   l->setParent(parent);
 1018   l->adjLevel(parent);
 1019   return (l);
 1020 }
 1021 
 1022 
 1023 // --- BODYSET CODE ------------------------------------------------------------
 1024 body* bodyset::newbody(){
 1025   body *p;
 1026   if(numbody == maxbody){
 1027     cerr << "newbody: More than " << maxbody << "body; increase nbody\n";
 1028     return NULL;
 1029   }
 1030   p = &bodytab[numbody];
 1031   numbody++;
 1032   return p;
 1033 }
 1034 
 1035 // --- LEAF CODE ------------------------------------------------------------
 1036 void leaf::LeafAddBody(body *bp){                  // Assumes enough capacity
 1037   bp->setParent(this);
 1038   bp->setLevel(Level());
 1039   bp->setChildNum(num_bodies);
 1040   bodyp[num_bodies] = bp;
 1041   num_bodies++;
 1042 }
 1043 
 1044 void leaf::Subdivide(cell *parent, int idx){
 1045   cell *c;
 1046   leaf *le;
 1047   double tmpv[NDIM];
 1048   int i, index, xp[NDIM], num, Lev;
 1049   body *loctab[LEAFMAXBODIES], *p;
 1050 
 1051   num = num_bodies;
 1052   for(i=0; i < num; i++) {
 1053     loctab[i] = bodyp[i];
 1054     bodyp[i] = NULL;
 1055   }
 1056   num_bodies = 0;
 1057   c = Nbody.getcells()->makecell(parent);
 1058   c->setChildNum(child_num);
 1059   parent->setSubp(idx,c);
 1060   Lev = level;
 1061 
 1062   p = loctab[0];
 1063   p->PosvecStore(tmpv);
 1064   Nbody.intcoord(xp,tmpv);
 1065   index = subindex(xp, Lev);
 1066   c->setSubp(index,this);
 1067   child_num = index;
 1068   parent = c;
 1069   level = (Lev >> 1);
 1070   LeafAddBody(p);
 1071 
 1072   for(i=1; i < num; i++) {
 1073     p = loctab[i];
 1074     p->PosvecStore(tmpv);
 1075     Nbody.intcoord(xp,tmpv);
 1076     index = subindex(xp, Lev);
 1077     le = (leaf*)(c->Subp(index));
 1078     if(le == NULL) {
 1079       le = Nbody.getleaves()->makeleaf(c);
 1080       le->setChildNum(index);
 1081       c->setSubp(index,le);
 1082     }
 1083     le->LeafAddBody(p);
 1084   }
 1085 }
 1086 
 1087 
 1088 // --- CELL CODE ------------------------------------------------------------
 1089 void cell::CellAddBody(body *p, int *coords) {
 1090   int Lev, kidIndex;
 1091   node *ptr;
 1092   leaf *le;
 1093 
 1094   kidIndex = subindex(coords,level);
 1095   ptr = subp[kidIndex];
 1096   Lev = (level >> 1);
 1097   if(Lev != 0){
 1098     if(ptr == NULL) {
 1099       le = Nbody.getleaves()->makeleaf(this);
 1100       le->setChildNum(kidIndex);
 1101       subp[kidIndex] = le;
 1102       le->LeafAddBody(p);
 1103     } else {
 1104       if(ptr->Type() == LEAF){
 1105         le = (leaf*)ptr;
 1106         if(le->full()){
 1107           le->Subdivide(this,kidIndex);
 1108           CellAddBody(p,coords);
 1109         } else {
 1110           le->LeafAddBody(p);
 1111         }
 1112       } else {
 1113         ((cell*)ptr)->CellAddBody(p,coords);
 1114       }
 1115     }
 1116   } else {
 1117     cerr << " *** Error: Not enough levels in tree...(CellAddBody)\n";
 1118   }
 1119 }
 1120 
 1121 // ----------------------------------------------------------------------------
 1122 int node::subindex(int *x, int l){
 1123   int i, k, yes;
 1124 
 1125   i = 0;
 1126   yes = 0;
 1127   for(k = 0; k < NDIM; k++) {
 1128     if(((x[k] & l) && !yes) || (!(x[k] & l) && yes)) {
 1129       i += NSUB >> (k + 1);
 1130       yes = 1;
 1131     }
 1132     else yes = 0;
 1133   }
 1134   return (i);
 1135 }
 1136 
 1137 void NBodySystem::setbound(){
 1138   int i;
 1139   double tmp1[NDIM], tmp2[NDIM], side;
 1140   body *p;
 1141 
 1142   mincorner->vecSetUnit();
 1143   mincorner->vecScale(1E99);
 1144   maxcorner->vecSetUnit();
 1145   maxcorner->vecScale(-1E99);
 1146   for(i=0; i < nbody; i++) {
 1147     p = bodies->getbody(i);
 1148     p->PosvecStore(tmp1);
 1149     mincorner->vecMin(tmp1);
 1150     maxcorner->vecMax(tmp1);
 1151   }
 1152 
 1153   side = 0.0;
 1154   maxcorner->vecStore(tmp1);
 1155   mincorner->vecStore(tmp2);
 1156   tvecSub(tmp1,tmp2);
 1157   for(i=0; i < NDIM; i++)
 1158     if(side < tmp1[i])
 1159       side = tmp1[i];
 1160 
 1161   for(i=0; i < NDIM; i++)
 1162     tmp1[i] = -side/100000;
 1163   mincorner->vecAdd(tmp1);
 1164   size = 1.00002*side;
 1165 }
 1166 
 1167 
 1168 // ----------------------------------------------------------------------------
 1169 // MAKETREE: initialize tree structure for force calculation.
 1170 // ----------------------------------------------------------------------------
 1171 void NBodySystem::maketree(){
 1172   int i, xp[NDIM];
 1173   double v[NDIM];
 1174   body *p;
 1175 
 1176   for(i=0; i < nbody; i++){
 1177     p = bodies->getbody(i);
 1178     p->PosvecStore(v);
 1179     intcoord(xp,v);
 1180     BH_tree->CellAddBody(p,xp);
 1181   }
 1182 }
 1183 
 1184 
 1185 void node::computecofm(){  // <--- Parallelizable...
 1186   int i;
 1187   node *nptr;
 1188   body *p;
 1189   leaf *le;
 1190   cell *ce;
 1191   double m;
 1192   double tmpv[NDIM];
 1193 
 1194   if(type == LEAF){
 1195     le = (leaf*)this;
 1196     for(i=0; i < le->NumBody(); i++) {
 1197       p = le->Bodyp(i);
 1198       m = p->Mass();
 1199       mass += m;
 1200       p->PosvecStore(tmpv);
 1201       tvecScale(tmpv,m);
 1202       pos.vecAdd(tmpv);
 1203     }
 1204     pos.vecStore(tmpv);
 1205     tvecScale(tmpv,1.0/mass);
 1206     pos.vecLoad(tmpv);
 1207   } else {
 1208     ce = (cell*)this;
 1209     for(i=0; i < NSUB; i++){
 1210       nptr = ce->Subp(i);
 1211       if(nptr != NULL){
 1212         nptr->computecofm();
 1213         m = nptr->Mass();
 1214         mass += m;
 1215         nptr->PosvecStore(tmpv);
 1216         tvecScale(tmpv,m);
 1217         pos.vecAdd(tmpv);
 1218       }
 1219     }
 1220     pos.vecStore(tmpv);
 1221     tvecScale(tmpv,1.0/mass);
 1222     pos.vecLoad(tmpv);
 1223   }
 1224 }
 1225 
 1226 double body::subdiv(node *p, double *res) {
 1227   int i;
 1228   double drsq, d;
 1229 
 1230   drsq = 0.0;
 1231   for (i = 0; i < NDIM; i++) {
 1232     d = p->pos.vecField(i)-pos.vecField(i);
 1233     drsq += d*d;
 1234     res[i] = d;
 1235   }
 1236   return(drsq);
 1237 }
 1238 
 1239 inline double body::computeInter(node *p, double epsSq, double *res){
 1240   int i;
 1241   double drabs, inc, mor3;
 1242   double drsq, d;
 1243 
 1244   drsq = epsSq;
 1245   for (i = 0; i < NDIM; i++) {
 1246     d = p->pos.vecField(i)-pos.vecField(i);
 1247     drsq += d*d;
 1248     res[i] = d;
 1249   }
 1250   drabs = sqrt(drsq);
 1251   inc = p->Mass() / drabs;
 1252   mor3 = inc / drsq;
 1253   tvecScale(res,mor3);
 1254   return inc;
 1255 }
 1256 
 1257 void body::updatePhi(double inc) {
 1258   phi -= inc;
 1259 }
 1260 
 1261 void body::fastgravsub(node *p, double *ai, double pdrsq, double epssq) {
 1262   double drsq;
 1263   double drabs, inc, mor3;
 1264   double tmpv[NDIM];
 1265 
 1266   drsq = pdrsq + epssq;
 1267   drabs = sqrt(drsq);
 1268   inc = p->Mass() / drabs;
 1269 //  phi -= inc;
 1270   updatePhi(inc);
 1271   mor3 = inc / drsq;
 1272   tvecLoad(tmpv,ai);
 1273   tvecScale(tmpv,mor3);
 1274   acc.vecAdd(tmpv);
 1275 }
 1276 
 1277 void body::gravsub(node *p, double epsSq){
 1278   double phii;
 1279   double ai[NDIM];
 1280 
 1281   phii = computeInter(p,epsSq,ai);
 1282 //  phi -= phii;
 1283   updatePhi(phii);
 1284   acc.vecAdd(ai);
 1285 }
 1286 
 1287 inline void body::openCell(node *n, double tolsq, double dsq, double epssq){
 1288   int i;
 1289   node* nn;
 1290 
 1291   for(i=0; i < NSUB; i++) {
 1292     nn = ((cell*)n)->Subp(i);
 1293     if(nn != NULL){
 1294       walksub(nn,tolsq,(dsq/4.0),epssq);
 1295     }
 1296   }
 1297 }
 1298 
 1299 inline void body::openLeaf(node *n,double epssq, int iter){
 1300   int i;
 1301   body *bptr;
 1302 
 1303   for(i=0; i < iter; i++){
 1304     bptr = ((leaf*)n)->Bodyp(i);
 1305     if(bptr != this){
 1306       gravsub((node*)bptr, epssq);
 1307     }
 1308   }
 1309 }
 1310 
 1311 void body::walksub(node *n, double tolsq, double dsq, double epssq){
 1312   double drsq;
 1313   double ai[NDIM];
 1314   drsq = subdiv(n,ai);
 1315   if((tolsq * drsq) < dsq){
 1316     if((n->Type()) == CELL){
 1317       openCell(n,tolsq,dsq,epssq);
 1318     } else {
 1319       openLeaf(n,epssq,((leaf*)n)->NumBody());
 1320     }
 1321   } else {
 1322     fastgravsub(n, ai, drsq, epssq);
 1323   }
 1324 }
 1325 
 1326 void body::advance(double hts, double ts){
 1327   double dvel[NDIM], vel1[NDIM], dpos[NDIM];
 1328 
 1329   acc.vecStore(dvel);
 1330   vel.vecStore(vel1);
 1331   tvecScale(dvel,hts);
 1332   tvecAdd(vel1,dvel);
 1333   tvecLoad(dpos,vel1);
 1334   tvecScale(dpos,ts);
 1335   tvecScale(dvel,2.0);
 1336 
 1337   vel.vecAdd(dvel);
 1338   pos.vecAdd(dpos);
 1339 }
 1340 
 1341 void NBodySystem::Advance(double hts, double ts){
 1342   int i;
 1343   body *p;
 1344 
 1345   for(i=0; i < nbody; i++) {
 1346     p = bodies->getbody(i);
 1347     p->advance(hts,ts);
 1348   }
 1349 }
 1350 
 1351 void body::swapAcc() {
 1352   double tmpv[NDIM];
 1353 
 1354   phi = 0.0;
 1355   acc.vecStore(tmpv);
 1356   tmp.vecLoad(tmpv);
 1357   acc.vecClr();
 1358 }
 1359 
 1360 void body::startVel(double f){
 1361    double tmpv1[NDIM], tmpv2[NDIM];
 1362    acc.vecStore(tmpv1);
 1363    tmp.vecStore(tmpv2);
 1364    tvecSub(tmpv1,tmpv2);
 1365    tvecScale(tmpv1,f);
 1366    vel.vecAdd(tmpv1);
 1367 }
 1368 
 1369 void NBodySystem::SwapAccs(){
 1370   int i;
 1371   body *p;
 1372 
 1373   for(i=0; i < nbody; i++) {
 1374     p = bodies->getbody(i);
 1375     p->swapAcc();
 1376   }
 1377 }
 1378 
 1379 void NBodySystem::ComputeAccels(double tol, double eps){
 1380   int i;
 1381   body *p;
 1382   int nbody_local = nbody;
 1383 
 1384    #pragma omp parallel for private (p) schedule(static)
 1385   for(i=0; i < nbody_local; i++) {
 1386     p = bodies->getbody(i);
 1387     p->walksub(BH_tree, tol, size*size, eps);
 1388   }
 1389 }
 1390 
 1391 void NBodySystem::StartVels(double hts){
 1392   int i;
 1393   body *p;
 1394 
 1395   for(i=0; i < nbody; i++) {
 1396     p = bodies->getbody(i);
 1397     p->startVel(hts);
 1398   }
 1399 }
 1400 
 1401 void NBodySystem::computeForces(double tol, double hts, double eps){
 1402   SwapAccs();
 1403   ComputeAccels(tol,eps);
 1404   if(nstep > 0){
 1405     StartVels(hts);
 1406   }
 1407 }
 1408 
 1409 // ----------------------------------------------------------------------------
 1410 int main(int argc, char **argv){
 1411   Nbody.init();
 1412   Nbody.startrun(&simparms);
 1413   Nbody.initOutput(&simparms);
 1414   Nbody.setbound();
 1415   Nbody.clearTiming();
 1416   Nbody.clearStats();
 1417 
 1418   Nbody.stepsystem(&simparms);
 1419 
 1420   //  bodyset * bset = Nbody.getbodies();
 1421   //for (int i = 0; i < bset->size(); i++)
 1422   //  bset->
 1423 
 1424   Nbody.outputTiming();
 1425   //Nbody.outputStats();
 1426   Nbody.output(&simparms);
 1427   Nbody.outputEnergy();
 1428   //Nbody.dump();
 1429 }
 1430