GLAMERDOC++
Gravitational Lensing Code Library
simpleTree.h
1 
8 #ifndef SIMP_TREE_H_
9 #define SIMP_TREE_H_
10 
11 #include <queue>
12 
13 #include "standard.h"
14 #include "Tree.h"
15 //#include "lens_halos.h"
16 
20 struct BranchNB{
21 
22  BranchNB(int Ndim){
23  center = new PosType[Ndim*3];
24  boundary_p1 = &center[Ndim];
25  boundary_p2 = &center[2*Ndim];
26  }
27  ~BranchNB(){
28  delete center;
29  }
30 
32  IndexType *particles;
33  IndexType nparticles;
35  IndexType big_particle;
37  PosType *center;
39  int level;
40  unsigned long number;
42  PosType *boundary_p1;
44  PosType *boundary_p2;
45  BranchNB *child1;
46  BranchNB *child2;
48  BranchNB *prev;
51  BranchNB *brother;
52 };
53 
61 template<typename PType>
62 struct TreeNBStruct{
63  BranchNB *top;
64  BranchNB *current;
66  unsigned long Nbranches;
68  short Ndimensions;
70  PType *pp;
71 };
72 
73 //typedef struct TreeNBStruct * TreeNBHndl;
74 
85 template<typename PType>
86 class TreeSimple {
87 public:
88  TreeSimple(PType *xp,IndexType Npoints,int bucket = 5,int dimensions = 2,bool median = true);
89  virtual ~TreeSimple();
90 
92  template<typename T>
93  void PointsWithinCircle(T center[2],float radius,std::list<unsigned long> &neighborkist);
94 
96  template<typename T>
97  void PointsWithinEllipse(T center[2],float a_max,float a_min,float posangle,std::list<unsigned long> &neighborkist);
98 
100  //template<typename T>
101  //void NearestNeighbors(T *ray,int Nneighbors,float *rsph,IndexType *neighbors) const;
102 
103  template<typename T>
104  T NNDistance(T *ray ,int Nneighbors) const;
105 
106  class iterator{
107  public:
108 
109  iterator(BranchNB *branch){
110  current = top = branch;
111  }
112  iterator(const iterator &it){
113  current = it.current;
114  top = it.top;
115  }
116 
117  iterator & operator=(const iterator &it){
118  if(&it == this) return *this;
119  current = it.current;
120  top = it.top;
121 
122  return *this;
123  }
124 
125  bool operator==(const iterator &it){
126  return (current == it.current)*(top == it.top);
127  }
128 
129  BranchNB *operator*(){return current;}
130 
132  bool walk(bool decend){
133 
134  ++count;
135  if(decend && current->child1 != NULL){
136  current = current->child1;
137  return true;
138  }
139  if(decend && current->child2 != NULL){
140  current = current->child2;
141  return true;
142  }
143 
144  if(current->brother == top->brother) return false;
145 
146  current = current->brother;
147  return true;
148  }
149 
150  bool up(){
151  ++count;
152  if(current == top) return false;
153  current = current->prev;
154  return true;
155  }
156 
157  bool down(int child){
158  ++count;
159  if(child == 1){
160  if(current->child1 == NULL) return false;
161  current = current->child1;
162  return true;
163  }
164  if(child == 2){
165  if(current->child2 == NULL) return false;
166  current = current->child2;
167  return true;
168  }
169  return false;
170  }
171 
172  void moveTop(){
173  current = top;
174  }
175 
176  bool atleaf(){
177  return (current->child1 == NULL)*(current->child2 == NULL);
178  }
179  bool atTop(){
180  return (current == top);
181  }
182 
183  size_t getcout(){return count;}
184  private:
185  size_t count = 0;
186  BranchNB *current;
187  BranchNB *top;
188  };
189 protected:
190 
191  int Ndim;//,incell2;
192  int incell;
193  TreeNBStruct<PType>* tree;
194  IndexType *index;
195  IndexType Nparticles;
196  bool median_cut;
197  int Nbucket;
198  PosType realray[3];
199  PType *xxp;
200 
201  TreeNBStruct<PType>* BuildTreeNB(PType *xxp,IndexType Nparticles,IndexType *particles,int Ndimensions,PosType theta);
202  void _BuildTreeNB(TreeNBStruct<PType>* tree,IndexType nparticles,IndexType *particles);
203 
204  template <typename T>
205  void _findleaf(T *ray,TreeSimple::iterator &it) const;
206 
207 
208  template<typename T>
209  void _PointsWithin(T *ray,float *rmax,std::list<unsigned long> &neighborkist);
210 
211  template<typename T>
212  void _NearestNeighbors(T *ray,int Nneighbors,unsigned long *neighbors,PosType *rneighbors);
213 
214  //template<typename T>
215  //void _NearestNeighbors(T *ray,PosType *realray,int Nneighbors,unsigned long *neighbors,PosType *rneighbors
216  // ,TreeSimple<PType>::iterator &it,int &bool) const ;
217 
218  BranchNB *NewBranchNB(IndexType *particles,IndexType nparticles
219  ,PosType boundary_p1[],PosType boundary_p2[]
220  ,PosType center[],int level,unsigned long branchNBnumber);
221  void FreeBranchNB(BranchNB *branchNB);
222  TreeNBStruct<PType> * NewTreeNB(IndexType *particles,IndexType nparticles
223  ,PosType boundary_p1[],PosType boundary_p2[],
224  PosType center[],short Ndimensions);
225  void freeTreeNB(TreeNBStruct<PType> * tree);
226  short emptyTreeNB(TreeNBStruct<PType> * tree);
227  void _freeTreeNB(TreeNBStruct<PType> * tree,short child);
228  bool isEmptyNB(TreeNBStruct<PType> * tree);
229  bool atTopNB(TreeNBStruct<PType> * tree);
230  bool noChildNB(TreeNBStruct<PType> * tree);
231  bool offEndNB(TreeNBStruct<PType> * tree);
232  void getCurrentNB(TreeNBStruct<PType> * tree,IndexType *particles,IndexType *nparticles);
233  unsigned long getNbranchesNB(TreeNBStruct<PType> * tree);
234  void moveTopNB(TreeNBStruct<PType> * tree);
235  void moveUpNB(TreeNBStruct<PType> * tree);
236  void moveToChildNB(TreeNBStruct<PType> * tree,int child);
237  void insertChildToCurrentNB(TreeNBStruct<PType> * tree, IndexType *particles,IndexType nparticles
238  ,PosType boundary_p1[],PosType boundary_p2[]
239  ,PosType center[],int child);
240  void attachChildToCurrentNB(TreeNBStruct<PType> * tree,BranchNB &data,int child);
241  bool TreeNBWalkStep(TreeNBStruct<PType> * tree,bool allowDescent);
242 
243  inline bool atLeaf(){
244  return (tree->current->child1 == NULL)*(tree->current->child2 == NULL);
245  }
246  template<typename T>
247  inline bool inbox(const T* center,PosType *p1,PosType *p2) const{
248  return (center[0]>=p1[0])*(center[0]<=p2[0])*(center[1]>=p1[1])*(center[1]<=p2[1]);
249  }
250 
251  /*TreeNBStruct<PType> * rotate_simulation(PosType **xp,IndexType Nparticles,IndexType *particles
252  ,PosType **coord,PosType theta,float *rsph,float *mass
253  ,bool MultiRadius,bool MultiMass);
254  TreeNBStruct<PType> * rotate_project(PosType **xp,IndexType Nparticles,IndexType *particles
255  ,PosType **coord,PosType theta,float *rsph,float *mass
256  ,bool MultiRadius,bool MultiMass);
257  TreeNBStruct<PType> * spread_particles(PosType **xp,IndexType Nparticles,IndexType *particles
258  ,PosType theta,float *rsph,float *mass,bool MultiRadius,bool MultiMass);
259  void cuttoffscale(TreeNBStruct<PType> * tree,PosType *theta);*/
260 };
261 
262 
263 template<typename PType>
264 TreeSimple<PType>::TreeSimple(PType *xpt,IndexType Npoints,int bucket,int Ndimensions,bool median){
265  index = new IndexType[Npoints];
266  IndexType ii;
267 
268  Nbucket = bucket;
269  Ndim = Ndimensions;
270  median_cut = median;
271  Nparticles = Npoints;
272 
273  xxp = xpt;
274 
275  for(ii=0;ii<Npoints;++ii) index[ii] = ii;
276 
277  tree = TreeSimple::BuildTreeNB(xxp,Npoints,index,Ndimensions,0);
278 
279  return;
280 }
281 
282 template<typename PType>
283 TreeSimple<PType>::~TreeSimple()
284 {
285  freeTreeNB(tree);
286  delete[] index;
287  return;
288 }
289 
290 template<typename PType>
291 template<typename T>
292 void TreeSimple<PType>::PointsWithinEllipse(
293  T *ray
294  ,float rmax
295  ,float rmin
296  ,float posangle
297  ,std::list <unsigned long> &neighborlist
298 ){
299  PosType x,y,cs,sn;
300 
301 
302  if(rmax < rmin){float tmp = rmax; rmax=rmin; rmin=tmp;}
303 
304  if(rmax <=0.0 || rmin <= 0.0){ neighborlist.clear(); return; }
305 
306  // find point within a circle circumscribes the ellipse
307  PointsWithinCircle(ray,rmax,neighborlist);
308 
309  cs = cos(posangle);
310  sn = sin(posangle);
311  // go through points within the circle and reject points outside the ellipse
312  for( std::list<unsigned long>::iterator it = neighborlist.begin();it != neighborlist.end();){
313  x = xxp[*it][0]*cs - xxp[*it][1]*sn;
314  y = xxp[*it][0]*sn + xxp[*it][1]*cs;
315  if( ( pow(x/rmax,2) + pow(y/rmin,2) ) > 1) it = neighborlist.erase(it);
316  else ++it;
317  }
318  return;
319 }
320 
321 template<typename PType>
322 template<typename T>
323 void TreeSimple<PType>::PointsWithinCircle(
324  T *ray
325  ,float rmax
326  ,std::list <unsigned long> &neighborlist
327 ){
328 
329  neighborlist.clear();
330 
331  realray[0]=ray[0];
332  realray[1]=ray[1];
333 
334  //cout << "ray = " << ray[0] << " " << ray[1] << " rmax = " << rmax << endl;
335  moveTopNB(tree);
336  if( inbox(ray,tree->current->boundary_p1,tree->current->boundary_p2) == 0 ){
337 
338  ray[0] = (ray[0] > tree->current->boundary_p1[0]) ? ray[0] : tree->current->boundary_p1[0];
339  ray[0] = (ray[0] < tree->current->boundary_p2[0]) ? ray[0] : tree->current->boundary_p2[0];
340 
341  ray[1] = (ray[1] > tree->current->boundary_p1[1]) ? ray[1] : tree->current->boundary_p1[1];
342  ray[1] = (ray[1] < tree->current->boundary_p2[1]) ? ray[1] : tree->current->boundary_p2[1];
343 
344  //cout << "ray = " << ray[0] << " " << ray[1] << endl;
345  //ray[0]=DMAX(ray[0],tree->current->boundary_p1[0]);
346  //ray[0]=DMIN(ray[0],tree->current->boundary_p2[0]);
347 
348  //ray[1]=DMAX(ray[1],tree->current->boundary_p1[1]);
349  //ray[1]=DMIN(ray[1],tree->current->boundary_p2[1]);
350  }
351  incell=1;
352 
353  _PointsWithin<T>(ray,&rmax,neighborlist);
354 
355  return;
356 }
359 template<typename PType>
360 template<typename T>
361 void TreeSimple<PType>::_PointsWithin(T *ray,float *rmax,std::list <unsigned long> &neighborlist){
362 
363  int j,incell2=1;
364  unsigned long i;
365  PosType radius;
366  short pass;
367 
368  if(incell){ // not found cell yet
369 
370  if( inbox(ray,tree->current->boundary_p1,tree->current->boundary_p2) ){
371  //cout << " In box" << endl;
372 
373  // found the box small enough
374  if( Utilities::cutbox(ray,tree->current->boundary_p1,tree->current->boundary_p2,*rmax)==1
375  || atLeaf() ){
376  //cout << " Found cell" << endl;
377 
378  // whole box in circle or a leaf with ray in it
379 
380  incell=0;
381 
382  ray[0]=realray[0];
383  ray[1]=realray[1];
384 
385  if( atLeaf() ){
386  // if leaf calculate the distance to all the points in cell
387  for(i=0;i<tree->current->nparticles;++i){
388  for(j=0,radius=0.0;j<2;++j) radius+=pow(xxp[tree->current->particles[i]][j]-ray[j],2);
389  if( radius < *rmax**rmax ){
390  neighborlist.push_back(tree->current->particles[i]);
391  //cout << "add point to list" << neighborlist.size() << endl;
392  //InsertAfterCurrentKist(neighborlist,tree->current->particles[i]);
393  }
394  }
395  }else{ // put all of points in box into getCurrentKist(imagelist)
396  for(i=0;i<tree->current->nparticles;++i){
397  neighborlist.push_back(tree->current->particles[i]);
398  //InsertAfterCurrentKist(neighborlist,tree->current->particles[i]);
399  //cout << "add point to list" << neighborlist.size() << endl;
400 
401  }
402  }
403 
404  }else{ // keep going down the tree
405 
406  if(tree->current->child1 !=NULL){
407  moveToChildNB(tree,1);
408  _PointsWithin(ray,rmax,neighborlist);
409  moveUpNB(tree);
410 
411  incell2=incell;
412  }
413 
414  if(tree->current->child2 !=NULL){
415  moveToChildNB(tree,2);
416  _PointsWithin(ray,rmax,neighborlist);
417  moveUpNB(tree);
418  }
419 
420  // if ray found in second child go back to first to search for neighbors
421  if( (incell2==1) && (incell==0) ){
422  if(tree->current->child1 !=NULL){
423  moveToChildNB(tree,1);
424  _PointsWithin(ray,rmax,neighborlist);
425  moveUpNB(tree);
426  }
427  }
428  }
429  } // not in the box
430 
431  }else{ // found cell
432 
433  pass=Utilities::cutbox(ray,tree->current->boundary_p1,tree->current->boundary_p2,*rmax);
434  // does radius cut into the box
435  if( pass ){
436  //cout << " Cell found searching other cells" << endl;
437  if( atLeaf() ){ /* leaf case */
438 
439  for(i=0;i<tree->current->nparticles;++i){
440  for(j=0,radius=0.0;j<2;++j) radius+=pow(xxp[tree->current->particles[i]][j]-ray[j],2);
441  if( radius < *rmax**rmax ){
442  neighborlist.push_back(tree->current->particles[i]);
443  //InsertAfterCurrentKist(neighborlist,tree->current->particles[i]);
444  //cout << "add point to list" << neighborlist.size() << endl;
445 
446  }
447  }
448  }else if(pass==1){ // whole box is inside radius
449 
450  for(i=0;i<tree->current->nparticles;++i){
451  neighborlist.push_back(tree->current->particles[i]);
452  //InsertAfterCurrentKist(neighborlist,tree->current->particles[i]);
453  //cout << "add point to list" << neighborlist.size() << endl;
454 
455  }
456  }else{
457  if(tree->current->child1 !=NULL){
458  moveToChildNB(tree,1);
459  _PointsWithin(ray,rmax,neighborlist);
460  moveUpNB(tree);
461  }
462 
463  if(tree->current->child2 !=NULL){
464  moveToChildNB(tree,2);
465  _PointsWithin(ray,rmax,neighborlist);
466  moveUpNB(tree);
467  }
468  }
469 
470  }
471  }
472 
473  return;
474 }
475 
481 template<typename PType>
482 template<typename T>
483 T TreeSimple<PType>::NNDistance(
484  T * ray
485  ,int Nneighbors
486 ) const{
487 
488  if(tree->top->nparticles <= Nneighbors){
489  ERROR_MESSAGE();
490  printf("ERROR: in TreeSimple::NearestNeighbors, number of neighbors > total number of particles\n");
491  exit(1);
492  }
493 
494  std::vector<PosType> tmp(Ndim);
495  TreeSimple::iterator iter(tree->top);
496 
497  // if the real ray is outside of root box find the closest boundary point that is inside
498  for(int j=0;j<tree->Ndimensions;++j){
499  tmp[j] = (ray[j] > tree->current->boundary_p1[j]) ? ray[j] : tree->current->boundary_p1[j];
500  tmp[j] = (ray[j] < tree->current->boundary_p2[j]) ? tmp[j] : tree->current->boundary_p2[j];
501  }
502 
503  //std::cout << tmp[0] << " " << tmp[1] << " " << tmp[2] << std::endl;
504 
505  // find leaf
506  _findleaf(tmp.data(),iter);
507 
508  //PosType tmp = (10*(tree->top->boundary_p2[0]-tree->top->boundary_p1[0]));
509 
510  int levelup = 0; // ??
511  while((*iter)->nparticles < Nneighbors){ iter.up(); ++levelup;}
512  auto pass = iter;
513 
514  //std::vector<PosType> r2neighbors((*iter)->nparticles);
515  std::priority_queue<PosType> r2heap;
516 
517  int i;
518  for(i=0 ; i < Nneighbors ;++i){
519  PosType r2 = 0;
520  for(int j=0;j<Ndim;++j){
521  //r2neighbors[i] += pow(tree->pp[(*iter)->particles[i]][j]-ray[j],2);
522  r2 += pow(tree->pp[(*iter)->particles[i]][j]-ray[j],2);
523  }
524  r2heap.push(r2);
525  }
526 
527  for(; i < (*iter)->nparticles ;++i){
528  PosType r2 = 0;
529  for(int j=0;j<Ndim;++j){
530  //r2neighbors[i] += pow(tree->pp[(*iter)->particles[i]][j]-ray[j],2);
531  r2 += pow(tree->pp[(*iter)->particles[i]][j]-ray[j],2);
532  }
533  if(r2 < r2heap.top()){
534  r2heap.push(r2);
535  r2heap.pop();
536  }
537  }
538 
539  //std::sort(r2neighbors.begin(),r2neighbors.end());
540  //r2neighbors.resize(Nneighbors);
541 
542  //PosType bestr = sqrt( r2neighbors.back() );
543 
544  PosType bestr = sqrt( r2heap.top() );
545  /*
546  int cutbox = 1;
547  for(int dim = 0 ; dim < Ndim ; ++dim){
548  cutbox *= ( (ray[dim] - bestr) > (*iter)->boundary_p1[dim] )
549  *( (ray[dim] + bestr) < (*iter)->boundary_p2[dim] );
550  }
551 
552  while(!cutbox && !iter.atTop()){
553  iter.up();
554  cutbox = 1;
555  for(int dim = 0 ; dim < Ndim ; ++dim){
556  cutbox *= ( (ray[dim] - bestr) > (*iter)->boundary_p1[dim] )
557  *( (ray[dim] + bestr) < (*iter)->boundary_p2[dim] );
558  }
559  }
560  */
561 
562  iter.moveTop();
563  //auto end = (*iter)->brother;
564  iter.walk(true);
565 
567  bool decend = true;
568  do{
569 
570  if(iter == pass){
571  decend =false;
572  }else{
573 
574  BranchNB * cbranch = *iter;
575 
576  //int cutbox = 1;
577  //for(int dim = 0 ; dim < Ndim ; ++dim){
578  // cutbox *= ( (ray[dim] + bestr) > (*iter)->boundary_p1[dim] )
579  // *( (ray[dim] - bestr) < (*iter)->boundary_p2[dim] );
580  //}
581 
582  short dim = (cbranch->level-1) % Ndim; // this box was cut in this dimension
583 
584  if(( (ray[dim] + bestr) > cbranch->boundary_p1[dim] )
585  *( (ray[dim] - bestr) < cbranch->boundary_p2[dim] ) ){
586  decend = true;
587 
588  if(iter.atleaf()){
589 
590  for(int i=0 ; i< cbranch->nparticles ;++i){
591  PosType r2 = 0;
592  for(int j=0;j<Ndim;++j){
593  r2 += pow(tree->pp[cbranch->particles[i]][j]-ray[j],2);
594  }
595 
596  if(r2 < bestr*bestr){
597 
598  r2heap.push(r2);
599  r2heap.pop();
600  bestr = sqrt( r2heap.top() );
601 
602  //int ii = r2neighbors.size() - 1;
603  //r2neighbors[ii] = r2;
604  //while(r2neighbors[ii] < r2neighbors[ii-1] && ii > 0){
605  // std::swap(r2neighbors[ii],r2neighbors[ii-1]);
606  // --ii;
607  //}
608  //bestr = sqrt( r2neighbors.back() );
609  }
610  }
611  }
612  }else{
613  decend = false;
614  }
615  }
616 // }while(iter.walk(decend) && !(*iter == end) );
617  }while(iter.walk(decend));
618 
619  //std::cout << " count " << iter.getcout() << " branches :" << tree->Nbranches << std::endl;
620  return bestr;
621 }
622 
623 template<typename PType>
624 template<typename T>
625 void TreeSimple<PType>::_findleaf(T *ray,TreeSimple::iterator &it) const {
626 
627  if( it.atleaf() ) return;
628 
629  int d = (*it)->level % Ndim;
630 
631  if( (*it)->child1->boundary_p2[d] > ray[d] ){
632  it.down(1);
633  _findleaf(ray,it);
634  }else{
635  it.down(2);
636  _findleaf(ray,it);
637  }
638 
639  return;
640 }
641 
642 /*
643 template<typename PType>
644 template<typename T>
645 void TreeSimple<PType>::_NearestNeighbors(T *ray,PosType *rlray,int Nneighbors
646  ,unsigned long *neighbors,PosType *rneighbors
647  ,TreeSimple::iterator &it,bool &notfound) const {
648 
649  int incellNB2=1;
650  IndexType i;
651  short j;
652 
653  if(notfound){ // not found cell yet
654 
655  if( inbox(ray,(*it)->boundary_p1,(*it)->boundary_p2) ){
656 
657  // found the box small enough
658  if( (*it)->nparticles <= Nneighbors+Nbucket ){
659  notfound=0;
660  for(j=0;j<tree->Ndimensions;++j) ray[j]=rlray[j];
661 
662  // calculate the distance to all the particles in cell
663  for(i=0;i<(*it)->nparticles;++i){
664  for(j=0,rneighbors[i]=0.0;j<tree->Ndimensions;++j){
665  rneighbors[i] += pow(tree->pp[(*it)->particles[i]][j]-ray[j],2);
666  }
667  rneighbors[i]=sqrt( rneighbors[i] );
668  assert(rneighbors[i] < 10);
669  neighbors[i]=(*it)->particles[i];
670  }
671 
672  Utilities::quicksort(neighbors,rneighbors,(*it)->nparticles);
673 
674 
675  }else{ // keep going down the tree
676 
677  if(it.down(1)){
678  //moveToChildNB(tree,1);
679  _NearestNeighbors(ray,rlray,Nneighbors,neighbors,rneighbors,it,notfound);
680  //printf("moving up from level %i\n",(*it)->level);
681  //moveUpNB(tree);
682  it.up();
683  incellNB2=notfound;
684  }
685 
686  if(it.down(2)){
687  //printf("moving to child2 from level %i\n",(*it)->level);
688  //moveToChildNB(tree,2);
689  _NearestNeighbors(ray,rlray,Nneighbors,neighbors,rneighbors,it,notfound);
690  //printf("moving up from level %i\n",(*it)->level);
691  //moveUpNB(tree);
692  it.up();
693  }
694 
695  // if ray found in second child go back to first to search for neighbors
696  if( (incellNB2==1) && (notfound==0) ){
697  if(it.down(1)){
698  //moveToChildNB(tree,1);
699  _NearestNeighbors(ray,rlray,Nneighbors,neighbors,rneighbors,it,notfound);
700  //moveUpNB(tree);
701  it.up();
702  }
703  }
704  }
705  }
706  }else{ // found cell
707  // does radius cut into the box
708  if( Utilities::cutbox(ray,(*it)->boundary_p1,(*it)->boundary_p2,rneighbors[Nneighbors-1]) ){
709 
710  if( ((*it)->child1 == NULL)*((*it)->child2 == NULL)){ // leaf case
711 
712  // combine found neighbors with particles in box and resort
713  for(i=Nneighbors;i<((*it)->nparticles+Nneighbors);++i){
714  for(j=0,rneighbors[i]=0.0;j<tree->Ndimensions;++j){
715  rneighbors[i]+=pow(tree->pp[(*it)->particles[i-Nneighbors]][j]-ray[j],2);
716  }
717  rneighbors[i]=sqrt( rneighbors[i] );
718  assert(rneighbors[i] < 10);
719  neighbors[i]=(*it)->particles[i-Nneighbors];
720  }
721 
722  Utilities::quicksort(neighbors,rneighbors,Nneighbors+Nbucket);
723 
724  }else{
725 
726  if(it.down(1)){
727  //moveToChildNB(tree,1);
728  _NearestNeighbors(ray,rlray,Nneighbors,neighbors,rneighbors,it,notfound);
729  //moveUpNB(tree);
730  it.up();
731  }
732 
733  if(it.down(2)){
734  //moveToChildNB(tree,2);
735  _NearestNeighbors(ray,rlray,Nneighbors,neighbors,rneighbors,it,notfound);
736  //moveUpNB(tree);
737  it.up();
738  }
739  }
740  }
741  }
742  return;
743 }
744 */
745 
746 template<typename PType>
747 template<typename T>
748 void TreeSimple<PType>::_NearestNeighbors(T *ray,int Nneighbors,unsigned long *neighbors,PosType *rneighbors){
749 
750  int incellNB2=1;
751  IndexType i;
752  short j;
753 
754  if(incell){ /* not found cell yet */
755 
756  if( inbox(ray,tree->current->boundary_p1,tree->current->boundary_p2) ){
757 
758  /* found the box small enough */
759  if( tree->current->nparticles <= Nneighbors+Nbucket ){
760  incell=0;
761  for(j=0;j<tree->Ndimensions;++j) ray[j]=realray[j];
762 
763  /* calculate the distance to all the particles in cell */
764  for(i=0;i<tree->current->nparticles;++i){
765  for(j=0,rneighbors[i]=0.0;j<tree->Ndimensions;++j){
766  rneighbors[i] += pow(tree->xp[tree->current->particles[i]][j]-ray[j],2);
767  }
768  rneighbors[i]=sqrt( rneighbors[i] );
769  assert(rneighbors[i] < 10);
770  neighbors[i]=tree->current->particles[i];
771  }
772 
773  Utilities::quicksort(neighbors,rneighbors,tree->current->nparticles);
774 
775  }else{ /* keep going down the tree */
776 
777  if(tree->current->child1 !=NULL){
778  moveToChildNB(tree,1);
779  _NearestNeighbors(ray,Nneighbors,neighbors,rneighbors);
780  /*printf("moving up from level %i\n",tree->current->level);*/
781  moveUpNB(tree);
782 
783  incellNB2=incell;
784  }
785 
786  if(tree->current->child2 !=NULL){
787  /*printf("moving to child2 from level %i\n",tree->current->level);*/
788  moveToChildNB(tree,2);
789  _NearestNeighbors(ray,Nneighbors,neighbors,rneighbors);
790  /*printf("moving up from level %i\n",tree->current->level);*/
791  moveUpNB(tree);
792  }
793 
795  if( (incellNB2==1) && (incell==0) ){
796  if(tree->current->child1 !=NULL){
797  moveToChildNB(tree,1);
798  _NearestNeighbors(ray,Nneighbors,neighbors,rneighbors);
799  moveUpNB(tree);
800  }
801  }
802  }
803  }
804  }else{ // found cell
805  /* does radius cut into the box */
806  if( Utilities::cutbox(ray,tree->current->boundary_p1,tree->current->boundary_p2,rneighbors[Nneighbors-1]) ){
807 
808  if( (tree->current->child1 == NULL)*(tree->current->child2 == NULL)){ /* leaf case */
809 
810  /* combine found neighbors with particles in box and resort */
811  for(i=Nneighbors;i<(tree->current->nparticles+Nneighbors);++i){
812  for(j=0,rneighbors[i]=0.0;j<tree->Ndimensions;++j){
813  rneighbors[i]+=pow(tree->xp[tree->current->particles[i-Nneighbors]][j]-ray[j],2);
814  }
815  rneighbors[i]=sqrt( rneighbors[i] );
816  assert(rneighbors[i] < 10);
817  neighbors[i]=tree->current->particles[i-Nneighbors];
818  }
819 
820  Utilities::quicksort(neighbors,rneighbors,Nneighbors+Nbucket);
821 
822  }else{
823 
824  if(tree->current->child1 !=NULL){
825  moveToChildNB(tree,1);
826  _NearestNeighbors(ray,Nneighbors,neighbors,rneighbors);
827  moveUpNB(tree);
828  }
829 
830  if(tree->current->child2 !=NULL){
831  moveToChildNB(tree,2);
832  _NearestNeighbors(ray,Nneighbors,neighbors,rneighbors);
833  moveUpNB(tree);
834  }
835  }
836  }
837  }
838  return;
839 }
840 
841 template<typename PType>
842 TreeNBStruct<PType> * TreeSimple<PType>::BuildTreeNB(PType *xp,IndexType Nparticles,IndexType *particles,int Ndims
843  ,PosType theta){
844  TreeNBStruct<PType> * tree;
845  IndexType i;
846  short j;
847  PosType p1[3],p2[3],center[3];
848 
849  assert(Ndim == Ndims);
850 
851  for(j=0;j<Ndim;++j){
852  p1[j]=xp[0][j];
853  p2[j]=xp[0][j];
854  }
855 
856  for(i=0;i<Nparticles;++i){
857  for(j=0;j<Ndim;++j){
858  if(xp[i][j] < p1[j] ) p1[j]=xp[i][j];
859  if(xp[i][j] > p2[j] ) p2[j]=xp[i][j];
860  }
861  }
862 
863  for(j=0;j<Ndim;++j) center[j]=(p1[j]+p2[j])/2;
864 
865  /* Initialize tree root */
866  tree=NewTreeNB(particles,Nparticles,p1,p2,center,Ndim);
867 
868  tree->pp = xp;
869 
870  /* build the tree */
871  _BuildTreeNB(tree,Nparticles,particles);
872 
873  /* visit every branch to find center of mass and cutoff scale */
874  moveTopNB(tree);
875 
876  return tree;
877 }
878 
879 
880 // tree must be created and first branch must be set before start
881 template<typename PType>
882 void TreeSimple<PType>::_BuildTreeNB(TreeNBStruct<PType> * tree,IndexType nparticles,IndexType *particles){
883  IndexType i,cut,dimension;
884  short j;
885  BranchNB *cbranch,branch1(Ndim),branch2(Ndim);
886  PosType xcut;
887  PosType *x;
888 
889  cbranch=tree->current; /* pointer to current branch */
890 
891  cbranch->center[0] = (cbranch->boundary_p1[0] + cbranch->boundary_p2[0])/2;
892  cbranch->center[1] = (cbranch->boundary_p1[1] + cbranch->boundary_p2[1])/2;
893  //cbranch->quad[0]=cbranch->quad[1]=cbranch->quad[2]=0;
894 
895  /* leaf case */
896  if(cbranch->nparticles <= Nbucket){
897  cbranch->big_particle = 0;
898  return;
899  }
900 
901  /* initialize boundaries to old boundaries */
902  for(j=0;j<Ndim;++j){
903  branch1.boundary_p1[j]=cbranch->boundary_p1[j];
904  branch1.boundary_p2[j]=cbranch->boundary_p2[j];
905 
906  branch2.boundary_p1[j]=cbranch->boundary_p1[j];
907  branch2.boundary_p2[j]=cbranch->boundary_p2[j];
908  }
909  cbranch->big_particle = 0;
910 
911  // **** makes sure force does not require nbucket at leaf
912 
913  /* set dimension to cut box */
914  dimension=(cbranch->level % Ndim);
915 
916  x=(PosType *)malloc((cbranch->nparticles-cbranch->big_particle)*sizeof(PosType));
917  for(i=cbranch->big_particle;i<cbranch->nparticles;++i) x[i]=tree->pp[particles[i]][dimension];
918 
919  if(median_cut){
920  Utilities::quicksort(particles,x,cbranch->nparticles-cbranch->big_particle);
921 
922  cut=(cbranch->nparticles-cbranch->big_particle)/2;
923  branch1.boundary_p2[dimension]=x[cut];
924  branch2.boundary_p1[dimension]=x[cut];
925  }else{
926  xcut=(cbranch->boundary_p1[dimension]+cbranch->boundary_p2[dimension])/2;
927  branch1.boundary_p2[dimension]=xcut;
928  branch2.boundary_p1[dimension]=xcut;
929 
930  Utilities::quickPartition(xcut,&cut,particles
931  ,x,cbranch->nparticles-cbranch->big_particle);
932  }
933 
934  /* set particle numbers and pointers to particles */
935  branch1.prev=cbranch;
936  branch1.nparticles=cut;
937  branch1.particles=particles+cbranch->big_particle;
938 
939  branch2.prev=cbranch;
940  branch2.nparticles=cbranch->nparticles-cbranch->big_particle - cut;
941  if(cut < (cbranch->nparticles-cbranch->big_particle) )
942  branch2.particles=&particles[cut+cbranch->big_particle];
943  else branch2.particles=NULL;
944 
945  free(x);
946 
947  if(branch1.nparticles > 0) attachChildToCurrentNB(tree,branch1,1);
948  if(branch2.nparticles > 0) attachChildToCurrentNB(tree,branch2,2);
949 
950  // work out brothers for children
951  if( (cbranch->child1 != NULL) && (cbranch->child2 != NULL) ){
952  cbranch->child1->brother = cbranch->child2;
953  cbranch->child2->brother = cbranch->brother;
954  }
955  if( (cbranch->child1 == NULL) && (cbranch->child2 != NULL) )
956  cbranch->child2->brother = cbranch->brother;
957  if( (cbranch->child1 != NULL) && (cbranch->child2 == NULL) )
958  cbranch->child1->brother = cbranch->brother;
959 
960 
961  if( branch1.nparticles > 0 ){
962  moveToChildNB(tree,1);
963  _BuildTreeNB(tree,branch1.nparticles,branch1.particles);
964  moveUpNB(tree);
965  }
966 
967  if(branch2.nparticles > 0 ){
968  moveToChildNB(tree,2);
969  _BuildTreeNB(tree,branch2.nparticles,branch2.particles);
970  moveUpNB(tree);
971  }
972 
973  return;
974 }
975 
985 /***** Constructors/Destructors *****/
986 
987 /************************************************************************
988  * NewBranchNB
989  * Returns pointer to new BranchNB struct. Initializes children pointers to NULL,
990  * and sets data field to input. Private.
991  ************************************************************************/
992 template<typename PType>
993 BranchNB *TreeSimple<PType>::NewBranchNB(IndexType *particles,IndexType nparticles
994  ,PosType boundary_p1[],PosType boundary_p2[]
995  ,PosType center[],int level,unsigned long branchNBnumber){
996 
997  BranchNB *branchNB = new BranchNB(Ndim);
998  int i;
999 
1000  if (!branchNB){
1001  ERROR_MESSAGE(); fprintf(stderr,"allocation failure in NewBranchNB()\n");
1002  exit(1);
1003  }
1004  branchNB->particles = particles;
1005  branchNB->nparticles = nparticles;
1006 
1007  for(i=0;i<Ndim;++i) branchNB->center[i]=center[i];
1008  branchNB->level=level;
1009 
1010  for(i=0;i<Ndim;++i){
1011  branchNB->boundary_p1[i]= boundary_p1[i];
1012  branchNB->boundary_p2[i]= boundary_p2[i];
1013  }
1014 
1015  branchNB->number=branchNBnumber;
1016 
1017  branchNB->child1 = NULL;
1018  branchNB->child2 = NULL;
1019  branchNB->prev = NULL;
1020  branchNB->brother = NULL;
1021 
1022  return(branchNB);
1023 }
1024 
1025 /************************************************************************
1026  * FreeBranchNB
1027  * Frees memory pointed to by branchNB. Private.
1028  ************************************************************************/
1029 template<typename PType>
1030 void TreeSimple<PType>::FreeBranchNB(BranchNB *branchNB){
1031 
1032  assert( branchNB != NULL);
1033  delete branchNB;
1034 
1035  return;
1036 }
1037 
1038 /************************************************************************
1039  * NewTreeNB
1040  * Returns pointer to new TreeNB struct. Initializes top, last, and
1041  * current pointers to NULL. Sets NbranchNBes field to 0. Exported.
1042  ************************************************************************/
1043 template<typename PType>
1044 TreeNBStruct<PType> * TreeSimple<PType>::NewTreeNB(IndexType *particles,IndexType nparticles
1045  ,PosType boundary_p1[],PosType boundary_p2[],
1046  PosType center[],short Ndimensions){
1047 
1048  TreeNBStruct<PType> * tree;
1049 
1050  tree = new TreeNBStruct<PType>;// *)malloc(sizeof(TreeNBStruct));
1051  if (!tree){
1052  ERROR_MESSAGE(); fprintf(stderr,"allocation failure in NewTreeNB()\n");
1053  exit(1);
1054  }
1055 
1056  tree->top= NewBranchNB(particles,nparticles,boundary_p1,boundary_p2,center,0,0);
1057  if (!(tree->top)){
1058  ERROR_MESSAGE(); fprintf(stderr,"allocation failure in NewTreeNB()\n");
1059  exit(1);
1060  }
1061 
1062  tree->Nbranches = 1;
1063  tree->current = tree->top;
1064  tree->Ndimensions=Ndimensions;
1065 
1066  return(tree);
1067 }
1068 
1069 template<typename PType>
1070 void TreeSimple<PType>::freeTreeNB(TreeNBStruct<PType> * tree){
1071  /* free treeNB
1072  * does not free the particle positions, masses or rsph
1073  */
1074 
1075  if(tree == NULL) return;
1076 
1077  emptyTreeNB(tree);
1078  FreeBranchNB(tree->top);
1079  delete tree;
1080 
1081  return;
1082 }
1083 
1084 template<typename PType>
1085 short TreeSimple<PType>::emptyTreeNB(TreeNBStruct<PType> * tree){
1086 
1087  moveTopNB(tree);
1088  _freeTreeNB(tree,0);
1089 
1090  assert(tree->Nbranches == 1);
1091 
1092  return 1;
1093 }
1094 
1095 template<typename PType>
1096 void TreeSimple<PType>::_freeTreeNB(TreeNBStruct<PType> * tree,short child){
1097  BranchNB *branch;
1098 
1099  assert( tree );
1100  assert( tree->current);
1101 
1102  if(tree->current->child1 != NULL){
1103  moveToChildNB(tree,1);
1104  _freeTreeNB(tree,1);
1105  }
1106 
1107  if(tree->current->child2 != NULL){
1108  moveToChildNB(tree,2);
1109  _freeTreeNB(tree,2);
1110  }
1111 
1112  if( (tree->current->child1 == NULL)*(tree->current->child2 == NULL) ){
1113 
1114  if(atTopNB(tree)) return;
1115 
1116  branch = tree->current;
1117  moveUpNB(tree);
1118  FreeBranchNB(branch);
1119 
1120  /*printf("*** removing branch %i number of branches %i\n",branch->number
1121  ,tree->Nbranches-1);*/
1122 
1123  if(child==1) tree->current->child1 = NULL;
1124  if(child==2) tree->current->child2 = NULL;
1125 
1126  --tree->Nbranches;
1127 
1128  return;
1129  }
1130 
1131  return;
1132 }
1133 
1134 /************************************************************************
1135  * isEmptyNB
1136  * Returns "true" if the TreeNB is empty and "false" otherwise. Exported.
1137  ************************************************************************/
1138 template<typename PType>
1139 bool TreeSimple<PType>::isEmptyNB(TreeNBStruct<PType> * tree){
1140 
1141  assert(tree != NULL);
1142  return(tree->Nbranches == 0);
1143 }
1144 
1145 /************************************************************************
1146  * atTop
1147  * Returns "true" if current is the same as top and "false" otherwise.
1148  * Exported.
1149  * Pre: !isEmptyNB(tree)
1150  ************************************************************************/
1151 template<typename PType>
1152 bool TreeSimple<PType>::atTopNB(TreeNBStruct<PType> * tree){
1153 
1154  assert(tree != NULL);
1155  if( isEmptyNB(tree) ){
1156  ERROR_MESSAGE();
1157  fprintf(stderr, "TreeNB Error: calling atTop() on empty tree\n");
1158  exit(1);
1159  }
1160  return(tree->current == tree->top);
1161 }
1162 
1163 /************************************************************************
1164  * noChild
1165  * Returns "true" if the child of the current branchNB does not exist and "false" otherwise.
1166  * Exported.
1167  * Pre: !isEmptyNB(tree)
1168  ************************************************************************/
1169 template<typename PType>
1170 bool TreeSimple<PType>::noChildNB(TreeNBStruct<PType> * tree){
1171 
1172  assert(tree != NULL);
1173  if( isEmptyNB(tree) ){
1174 
1175  ERROR_MESSAGE(); fprintf(stderr, "TreeNB Error: calling atTop() on empty tree\n");
1176  exit(1);
1177  }
1178 
1179  if( (tree->current->child1 == NULL) || (tree->current->child2 == NULL) ) return true;
1180  return false;
1181 }
1182 
1183 /************************************************************************
1184  * offEndNB
1185  * Returns "true" if current is off end and "false" otherwise. Exported.
1186  ************************************************************************/
1187 template<typename PType>
1188 bool TreeSimple<PType>::offEndNB(TreeNBStruct<PType> * tree){
1189 
1190  assert(tree != NULL);
1191  return(tree->current == NULL);
1192 }
1193 
1194 /************************************************************************
1195  * getCurrentNB
1196  * Returns the particuls of current. Exported.
1197  * Pre: !offEndNB(tree)
1198  ************************************************************************/
1199 template<typename PType>
1200 void TreeSimple<PType>::getCurrentNB(TreeNBStruct<PType> * tree,IndexType *particles,IndexType *nparticles){
1201 
1202  assert(tree != NULL);
1203  if( offEndNB(tree) ){
1204 
1205  ERROR_MESSAGE(); fprintf(stderr, "TreeNB Error: calling getCurrent() when current is off end\n");
1206  exit(1);
1207  }
1208 
1209  *nparticles=tree->current->nparticles;
1210  particles=tree->current->particles;
1211 
1212  return;
1213 }
1214 
1215 /************************************************************************
1216  * getNbranchesNB
1217  * Returns the NbranchNBes of tree. Exported.
1218  ************************************************************************/
1219 template<typename PType>
1220 unsigned long TreeSimple<PType>::getNbranchesNB(TreeNBStruct<PType> * tree){
1221 
1222  assert(tree != NULL);
1223  return(tree->Nbranches);
1224 }
1225 
1226 /***** Manipulation procedures *****/
1227 
1228 /************************************************************************
1229  * moveTopNB
1230  * Moves current to the front of tree. Exported.
1231  * Pre: !isEmptyNB(tree)
1232  ************************************************************************/
1233 template<typename PType>
1234 void TreeSimple<PType>::moveTopNB(TreeNBStruct<PType> * tree){
1235  //std::cout << tree << std::endl;
1236  //std::cout << tree->current << std::endl;
1237  //std::cout << tree->top << std::endl;
1238 
1239  assert(tree != NULL);
1240  if( isEmptyNB(tree) ){
1241 
1242  ERROR_MESSAGE(); fprintf(stderr, "TreeNB Error: calling moveTopNB() on empty tree\n");
1243  exit(1);
1244  }
1245 
1246 
1247  tree->current = tree->top;
1248 
1249  return;
1250 }
1251 
1252 /************************************************************************
1253  * movePrev
1254  * Moves current to the branchNB before it in tree. This can move current
1255  * off end. Exported.
1256  * Pre: !offEndNB(tree)
1257  ************************************************************************/
1258 template<typename PType>
1259 void TreeSimple<PType>::moveUpNB(TreeNBStruct<PType> * tree){
1260 
1261  assert(tree != NULL);
1262  if( offEndNB(tree) ){
1263  ERROR_MESSAGE(); fprintf(stderr, "TreeNB Error: call to moveUpNB() when current is off end\n");
1264  exit(1);
1265  }
1266  if( tree->current == tree->top ){
1267  ERROR_MESSAGE(); fprintf(stderr, "TreeNB Error: call to moveUpNB() tried to move off the top\n");
1268  exit(1);
1269  }
1270 
1271  tree->current = tree->current->prev; /* can move off end */
1272  return;
1273 }
1274 
1275 /************************************************************************
1276  * moveToChildNB
1277  * Moves current to child branchNB after it in tree. This can move current off
1278  * end. Exported.
1279  * Pre: !offEndNB(tree)
1280  ************************************************************************/
1281 template<typename PType>
1282 void TreeSimple<PType>::moveToChildNB(TreeNBStruct<PType> * tree,int child){
1283 
1284  assert(tree != NULL);
1285  if( offEndNB(tree) ){
1286 
1287  ERROR_MESSAGE(); fprintf(stderr, "TreeNB Error: calling moveChildren() when current is off end\n");
1288  exit(1);
1289  }
1290  if(child==1){
1291  if( tree->current->child1 == NULL ){
1292  ERROR_MESSAGE(); fprintf(stderr, "TreeNB Error: moveToChildNB() typing to move to child1 when it doesn't exist\n");
1293  exit(1);
1294  }
1295  tree->current = tree->current->child1;
1296  }
1297  if(child==2){
1298  if( tree->current->child2 == NULL ){
1299  ERROR_MESSAGE(); fprintf(stderr, "TreeNB Error: moveToChildNB() typing to move to child2 when it doesn't exist\n");
1300  exit(1);
1301  }
1302  tree->current = tree->current->child2;
1303  }
1304  return;
1305 }
1306 
1307 /************************************************************************
1308  * insertAfterCurrent
1309  * Inserts a new BranchNB after the current branchNB in the tree and sets the
1310  * data field of the new BranchNB to input. Exported.
1311  * Pre: !offEndNB(tree)
1312  ************************************************************************/
1313 template<typename PType>
1314 void TreeSimple<PType>::insertChildToCurrentNB(TreeNBStruct<PType> * tree, IndexType *particles,IndexType nparticles
1315  ,PosType boundary_p1[],PosType boundary_p2[]
1316  ,PosType center[],int child){
1317 
1318  BranchNB *branchNB;
1319 
1320  /*printf("attaching child%i current paricle number %i\n",child,tree->current->nparticles);*/
1321 
1322  branchNB = NewBranchNB(particles,nparticles,boundary_p1,boundary_p2,center
1323  ,tree->current->level+1,tree->Nbranches);
1324 
1325  assert(tree != NULL);
1326 
1327  if( offEndNB(tree) ){
1328 
1329  ERROR_MESSAGE(); fprintf(stderr, "TreeNB Error: calling insertChildToCurrentNB() when current is off end\n");
1330  exit(1);
1331  }
1332 
1333  branchNB->prev = tree->current;
1334 
1335  if(child==1){
1336  if(tree->current->child1 != NULL){
1337  ERROR_MESSAGE(); fprintf(stderr, "TreeNB Error: calling insertChildToCurrentNB() when child1 alread exists\n");
1338  exit(1);
1339  }
1340  tree->current->child1 = branchNB;
1341  }
1342  if(child==2){
1343  if(tree->current->child2 != NULL){
1344  ERROR_MESSAGE();
1345  fprintf(stderr, "TreeNB Error: calling insertChildToCurrentNB() when child2 alread exists\n current level=%i Nbranches=%li\n"
1346  ,tree->current->level,tree->Nbranches);
1347  exit(1);
1348  }
1349  tree->current->child2 = branchNB;
1350  }
1351 
1352  tree->Nbranches++;
1353 
1354  return;
1355 }
1356 
1357 /* same as above but takes a branchNB structure */
1358 template<typename PType>
1359 void TreeSimple<PType>::attachChildToCurrentNB(TreeNBStruct<PType> * tree,BranchNB &data,int child){
1360 
1361  insertChildToCurrentNB(tree,data.particles,data.nparticles
1362  ,data.boundary_p1,data.boundary_p2,data.center,child);
1363  return;
1364 }
1365 
1366 // step for walking tree by iteration instead of recursion
1367 template<typename PType>
1368 bool TreeSimple<PType>::TreeNBWalkStep(TreeNBStruct<PType> * tree,bool allowDescent){
1369  if(allowDescent && tree->current->child1 != NULL){
1370  moveToChildNB(tree,1);
1371  return true;
1372  }
1373  if(allowDescent && tree->current->child2 != NULL){
1374  moveToChildNB(tree,2);
1375  return true;
1376  }
1377 
1378  if(tree->current->brother != NULL){
1379  tree->current=tree->current->brother;
1380  return true;
1381  }
1382  return false;
1383 }
1384 
1385 #endif /* SIMP_TREE_H_ */
TreeSimple
A C++ class wrapper for the bianary treeNB used in the Nbody force calculation, but also useful for g...
Definition: simpleTree.h:86
TreeSimple::NewBranchNB
BranchNB * NewBranchNB(IndexType *particles, IndexType nparticles, PosType boundary_p1[], PosType boundary_p2[], PosType center[], int level, unsigned long branchNBnumber)
Definition: simpleTree.h:993
BranchNB::boundary_p1
PosType * boundary_p1
bottom, left, back corner of box
Definition: simpleTree.h:42
BranchNB::brother
BranchNB * brother
Definition: simpleTree.h:51
BranchNB::big_particle
IndexType big_particle
the number of particles that aren't in children
Definition: simpleTree.h:35
TreeSimple::_PointsWithin
void _PointsWithin(T *ray, float *rmax, std::list< unsigned long > &neighborkist)
Definition: simpleTree.h:361
BranchNB::prev
BranchNB * prev
father of branch
Definition: simpleTree.h:48
BranchNB::boundary_p2
PosType * boundary_p2
top, right, front corner of box
Definition: simpleTree.h:44
TreeNBStruct::Ndimensions
short Ndimensions
Dimension of tree, 2 or 3. This will dictate how the force is calculated.
Definition: simpleTree.h:68
Utilities::median
T median(std::vector< T > vec)
find median of vector
Definition: utilities.h:243
TreeNBStruct
TreeNBStruct: Tree structure used for force calculation with particles (i.e. stars,...
Definition: simpleTree.h:62
standard.h
TreeSimple::_NearestNeighbors
void _NearestNeighbors(T *ray, int Nneighbors, unsigned long *neighbors, PosType *rneighbors)
Definition: simpleTree.h:748
TreeSimple::PointsWithinCircle
void PointsWithinCircle(T center[2], float radius, std::list< unsigned long > &neighborkist)
Finds the points within a circle around center and puts their index numbers in a list.
BranchNB
Box representing a branch in a tree. It has four children. Used in TreeNBStruct which is used in Tree...
Definition: simpleTree.h:20
TreeSimple::NNDistance
T NNDistance(T *ray, int Nneighbors) const
Finds the nearest N neighbors and puts their index numbers in an array, also returns the distance to ...
Definition: simpleTree.h:483
BranchNB::particles
IndexType * particles
array of particles in BranchNB
Definition: simpleTree.h:32
TreeNBStruct::pp
PType * pp
Array of particle positions.
Definition: simpleTree.h:70
TreeNBStruct::Nbranches
unsigned long Nbranches
number of branches in tree
Definition: simpleTree.h:66
BranchNB::level
int level
level in tree
Definition: simpleTree.h:39
BranchNB::center
PosType * center
center of mass
Definition: simpleTree.h:37
TreeSimple::PointsWithinEllipse
void PointsWithinEllipse(T center[2], float a_max, float a_min, float posangle, std::list< unsigned long > &neighborkist)
Finds the points within an ellipse around center and puts their index numbers in a list.
Generated on Sun Sep 27 2026 14:12:56 for GLAMERDOC++ by   doxygen 1.8.18