// // treecode.C #define PR(x) cerr << #x << " = " << x << " " #define PRC(x) cerr << #x << " = " << x << ", " #define PRL(x) cerr << #x << " = " << x << "\n" #include "cpgplot.h" #include #include #include #define real double #include "BHtree.h" void initgraph() { if(cpgopen("") !=1 ) exit(-1); cpgask(0); } void plot_particles(particle * p, int n, double xmax, int first) { static float * xp; static float * yp; static int nbuf = -1; cpgbbuf(); if (first == 1){ cpgpage(); cpgenv(-xmax, xmax, -xmax, xmax,1,-2); } if ( nbuf < n){ if (nbuf > 0){ delete [] xp; delete [] yp; } xp = new float[n]; yp = new float[n]; nbuf = n; } for(int i = 0;irsize) { rsize *= 2; } } } return rsize; } void calculate_gravity(bhnode* bn, int nnodes, particle * pp, int n, real eps2, real theta) { real rsize = calculate_size(pp,n); bn->assign_root(myvector(0.0), rsize*2, pp, n); int heap_remainder = nnodes-1; bhnode * btmp = bn + 1; bn->create_tree_recursive(btmp, heap_remainder); // bn->dump(); // PRL(bn->sanity_check()); bn->set_cm_quantities(); clear_acc_and_phi(pp, n); particle * p = pp; for(int i = 0; ipos);PR(p->phi); PRL(p->acc); p++; } } void integrate(bhnode * bn, int nnodes, particle * pp, int n, real eps2, real theta, real dt) { for(int i = 0;i> n ; pp = new particle[n]; double rsize = 1.0; cerr << "Enter power index:"; real power_index; cin >> power_index ; cerr << "power index = " << power_index <> eps2 >>theta >>dt >>tend >> iout; cerr << "eps2=" << eps2 << " theta=" <> xmax_plot; cerr << "xmax for plot = " << xmax_plot <