GLAMERDOC++
Gravitational Lensing Code Library
quadTreeHalo.h
1 
7 //**********************************************************************************************************
8 /*
9  * quadTree.cpp
10  *
11  */
12 
13 /*
14  * Programmer: R Ben Metcalf
15  */
16 #ifndef QUAD_TREE_HALO_
17 #define QUAD_TREE_HALO_
18 
19 //#include "utilities_slsim.h"
20 //#include "utilities.h"
21 #include "Tree.h"
22 #include "qTreeNB.h"
23 
43 template<typename PType>
44 class TreeQuadHalos {
45 public:
46  TreeQuadHalos(
47  PType *xpt
48  ,IndexType Npoints
49  ,PosType my_inv_area = 0
50  ,int bucket = 5
51  ,PosType theta_force = 0.1
52  ,PosType my_inv_screening_scale = 0
53  );
54 
55  ~TreeQuadHalos();
56 
57  void force2D(PosType const *ray,PosType *alpha,KappaType *kappa,KappaType *gamma
58  ,KappaType *phi) const;
59 
61  void neighbors(PosType ray[],PosType rmax,std::list<IndexType> &neighbors) const;
62 
63  void printParticlesInBranch(unsigned long number);
64  void printBranchs(int level = -1);
65 
66 protected:
67 
68  PType *xxp;
69  double inv_area;
70 
71  IndexType Nparticles;
72  //PosType sigma_background;
73  int Nbucket;
74 
75  PosType force_theta;
76  PosType max_range;
77 
78  std::unique_ptr<QTreeNB<PType> > tree;
79  std::vector<IndexType> index;
80  std::vector<PosType> workspace;
81 
82  PosType realray[2];
83  int incell,incell2;
84 
86  bool periodic_buffer;
87  PosType inv_screening_scale2;
88  PosType original_xl; // x-axis size of simulation used for peridic buffering. Requrement that it top branch be square my make it differ from the size of top branch.
89  PosType original_yl; // x-axis size of simulation used for peridic buffering.
90 
91  //QTreeNB<PType> * BuildQTreeNB(PType *xp,IndexType Nparticles,IndexType *particles);
92  void BuildQTreeNB(PType *xp,IndexType Nparticles);
93  void _BuildQTreeNB(IndexType nparticles,IndexType *particles);
94 
95  inline short WhichQuad(PosType *x,QBranch &branch);
96 
97  //inline bool atLeaf();
98  inline bool inbox(const PosType *ray,const PosType *p1,const PosType *p2){
99  return (ray[0]>=p1[0])*(ray[0]<=p2[0])*(ray[1]>=p1[1])*(ray[1]<=p2[1]);
100  }
101  //int cutbox(PosType *ray,PosType *p1,PosType *p2,float rmax);
102 
103  void CalcMoments();
104  void rotate_coordinates(PosType **coord);
105 
106  // Internal profiles for a Gaussian particle
107  inline PosType alpha_h(PosType r2s2,PosType sigma) const{
108  return (sigma > 0.0 ) ? ( exp(-0.5*r2s2) - 1.0 ) : -1.0;
109  }
110  inline PosType kappa_h(PosType r2s2,PosType sigma) const{
111  return 0.5*r2s2*exp(-0.5*r2s2);
112  }
113  inline PosType gamma_h(PosType r2s2,PosType sigma) const{
114  return (sigma > 0.0 ) ? (-2.0 + (2.0 + r2s2)*exp(-0.5*r2s2) ) : -2.0;
115  }
116  inline PosType phi_h(PosType r2s2,PosType sigma) const{
117  //ERROR_MESSAGE(); // not yet written
118  //exit(1);
119  return 0;
120  }
121 
122  void walkTree_iter(QBiterator<PType> &treeit, PosType const *ray,PosType *alpha,KappaType *kappa
123  ,KappaType *gamma,KappaType *phi) const;
124 
125  PosType phiintconst;
126 
127  struct QBranch{
128  QBranch(QBranch *parent):prev(parent){
129  static unsigned long n=0;
130  child0 = NULL;
131  child1 = NULL;
132  child2 = NULL;
133  child3 = NULL;
134  brother = NULL;
135  particles = NULL;
136  nparticles = 0;
137  number = n;
138  ++n;
139  };
140 
141  ~QBranch(){ };
142 
144  IndexType *particles;
145  std::unique_ptr<IndexType[]> big_particles;
146  IndexType nparticles;
148  IndexType Nbig_particles;
150  PosType maxrsph;
152  PosType center[2];
153  PosType mass;
155  int level;
156  unsigned long number;
158  PosType boundary_p1[2];
160  PosType boundary_p2[2];
161  PosType boxsize2;
162 
163  QBranch *child0;
164  QBranch *child1;
165  QBranch *child2;
166  QBranch *child3;
167 
169  QBranch *prev;
172  QBranch *brother;
173 
175  PosType quad[3];
177  PosType rmax;
179  PosType r2crit_angle;
180  /* force calculation */
181  PosType rcrit_part;
182  //PosType cm[2]; /* projected center of mass */
183  };
184 };
185 
188 template<typename PType>
189 TreeQuadHalos<PType>::TreeQuadHalos(
190  PType *xpt
191  ,IndexType Npoints
192  ,PosType my_inv_area
193  ,int bucket
194  ,PosType theta_force
195  ,bool my_periodic_buffer
196  ,PosType my_inv_screening_scale
197  ,PosType maximum_range
198 
199 ):
200 xxp(xpt)
201 ,inv_area(my_inv_area)
202 ,Nparticles(Npoints)
203 ,Nbucket(bucket)
204 ,force_theta(theta_force)
205 ,max_range(maximum_range)
206 ,periodic_buffer(my_periodic_buffer)
207 ,inv_screening_scale2(my_inv_screening_scale*my_inv_screening_scale)
208 {
209  index.resize(Npoints);
210  IndexType ii;
211  if(mass_fixed <= 0 ) MultiMass = true; else MultiMass = false;
212  if(size_fixed <= 0 ) MultiRadius = true; else MultiRadius = false;
213 
214  for(ii=0;ii<Npoints;++ii) index[ii] = ii;
215 
216  if(Npoints > 0){
217  BuildQTreeNB(xxp,Npoints);
218  CalcMoments();
219  }
220 
221  phiintconst = (120*log(2.) - 631.)/840 + 19./70;
222 
223  return;
224 }
225 
227 template<typename PType>
228 TreeQuadHalos<PType>::~TreeQuadHalos()
229 {
230  if(Nparticles == 0) return;
231  return;
232 }
233 
234 template<typename PType>
235 //QTreeNB<PType> * TreeQuadHalos<PType>::BuildQTreeNB(PType *xxp,IndexType Nparticles,IndexType *particles){
236 void TreeQuadHalos<PType>::BuildQTreeNB(PType *xxp,IndexType Nparticles){
237  IndexType i;
238  short j;
239  PosType p1[2],p2[2];
240 
241  for(j=0;j<2;++j){
242  p1[j]=xxp[0][j];
243  p2[j]=xxp[0][j];
244  }
245 
246  // Find borders that enclose all particles
247  for(i=0;i<Nparticles;++i){
248  for(j=0;j<2;++j){
249  if(xxp[i][j] < p1[j] ) p1[j]=xxp[i][j];
250  if(xxp[i][j] > p2[j] ) p2[j]=xxp[i][j];
251  }
252  }
253 
254  if(Nparticles == 1){
255  p1[0]=-0.25;
256  p1[1]=-0.25;
257  p2[0]=0.25;
258  p2[1]=0.25;
259  }
260 
261  // store true dimensions of simulation
262  PosType lengths[2] = {p2[0]-p1[0],p2[1]-p1[1]};
263  original_xl = lengths[0];
264  original_yl = lengths[1];
265 
266  if(lengths[0] == 0 || lengths[1] == 0){
267  throw std::invalid_argument("particles in same place.");
268  }
269 
270  // If region is not square, make it square.
271  j = lengths[0] > lengths[1] ? 1 : 0;
272  p2[j] = p1[j] + lengths[!j];
273 
274  /* Initialize tree root */
275  tree.reset( new QTreeNB<PType>(xxp,index.data(),Nparticles,p1,p2) );
276 
277  /* build the tree */
278  workspace.resize(Nparticles);
279  _BuildQTreeNB(Nparticles,index.data());
280  workspace.clear();
281  workspace.shrink_to_fit();
282 
283  /* visit every branch to find center of mass and cutoff scale */
284  tree->moveTop();
285 
286  return;
287 }
288 
290 template<typename PType>
291 inline short TreeQuadHalos<PType>::WhichQuad(PosType *x,QBranch &branch){
292  return (x[0] < branch.center[0]) + 2*(x[1] < branch.center[1]);
293 }
294 
296 template<typename PType>
297 void TreeQuadHalos<PType>::_BuildQTreeNB(IndexType nparticles,IndexType *particles){
298 
299  QBranch *cbranch = tree->current; // pointer to current branch
300  IndexType i,j,cut,cut2;
301 
302  cbranch->center[0] = (cbranch->boundary_p1[0] + cbranch->boundary_p2[0])/2;
303  cbranch->center[1] = (cbranch->boundary_p1[1] + cbranch->boundary_p2[1])/2;
304  cbranch->quad[0] = cbranch->quad[1]=cbranch->quad[2]=0;
305  cbranch->mass = 0.0;
306 
307  double theta_range = 2*force_theta;
308  if(max_range > 0){
309  double boxsize = 1.732*(cbranch->boundary_p2[0] - cbranch->boundary_p1[0]);
310  theta_range = boxsize / MAX(sqrt(cbranch->center[0]*cbranch->center[0]
311  + cbranch->center[1]*cbranch->center[1] )
312  - max_range, boxsize );
313  }
314 
315  // leaf case
316  if(cbranch->nparticles <= Nbucket
317  || force_theta > theta_range
318  ){
319  PosType r;
320  cbranch->Nbig_particles = 0;
321  PosType boxsize = cbranch->boundary_p2[0]-cbranch->boundary_p1[0];
322  for(i=0;i<cbranch->nparticles;++i){
323  if(xxp[particles[i]*MultiRadius].size() < boxsize) ++cbranch->Nbig_particles;
324  }
325  if(cbranch->Nbig_particles){
326  cbranch->big_particles.reset(new IndexType[cbranch->Nbig_particles]);
327 
328  for(i=0,j=0;i<cbranch->nparticles;++i){
329  if(xxp[particles[i]*MultiRadius].size() < boxsize) cbranch->big_particles[j++] = particles[i];
330  }
331  }else{
332  cbranch->big_particles.reset(nullptr);
333  }
334 
335  return;
336  }
337 
338  // find particles too big to be in children
339 
340  //PosType *x = new PosType[cbranch->nparticles];
341  PosType *x = workspace.data();
342 
343  cbranch->Nbig_particles=0;
344 
345  if(MultiRadius){
346  // store the particles that are too large to be in a child at the end of the list of particles in cbranch
347  for(i=0;i<cbranch->nparticles;++i){
348  x[i] = xxp[particles[i]].size();
349  }
350  PosType maxsize = (cbranch->boundary_p2[0]-cbranch->boundary_p1[0])/2;
351 
352  // sort particles in size
353  //quicksort(particles,x,cbranch->nparticles);
354  Utilities::quickPartition(maxsize,&cut,particles,x,cbranch->nparticles);
355 
356  if(cut < cbranch->nparticles){
357  if(tree->atTop()){
358  cbranch->Nbig_particles = cut2 = cbranch->nparticles-cut;
359  }else{
360  Utilities::quickPartition(2*maxsize,&cut2,&particles[cut],&x[cut],cbranch->nparticles-cut);
361  cbranch->Nbig_particles = cut2;
362  }
363 
364  cbranch->big_particles.reset(new IndexType[cbranch->Nbig_particles]);
365  for(i=cut;i<(cut+cut2);++i) cbranch->big_particles[i-cut] = particles[i];
366  }
367 
368  }else{
369  //x[0] = xxp[0].Size();
370  x[0] = size_fixed;
371  if(x[0] > (cbranch->boundary_p2[0]-cbranch->boundary_p1[0])/2
372  && x[0] < (cbranch->boundary_p2[0]-cbranch->boundary_p1[0])){
373  cbranch->Nbig_particles = cbranch->nparticles;
374  //cbranch->big_particles = cbranch->particles;
375  cbranch->big_particles.reset(new IndexType[cbranch->Nbig_particles]);
376  for(i=0;i<cbranch->nparticles;++i) cbranch->big_particles[i] = cbranch->particles[particles[i]];
377  }
378  }
379 
380  assert( cbranch->nparticles >= cbranch->Nbig_particles );
381  IndexType NpInChildren = cbranch->nparticles;// - cbranch->Nbig_particles;
382  assert(NpInChildren >= 0);
383 
384  if(NpInChildren == 0){
385  //delete[] x;
386  return;
387  }
388 
389  IndexType cutx,cuty;
390  PosType xcut,ycut;
391 
392  QBranch *child0 = new QBranch(cbranch);
393  QBranch *child1 = new QBranch(cbranch);
394  QBranch *child2 = new QBranch(cbranch);
395  QBranch *child3 = new QBranch(cbranch);
396 
397  tree->attachChildrenToCurrent(child0,child1,child2,child3);
398 
399  // initialize boundaries of children
400 
401  child0->boundary_p1[0]=cbranch->center[0];
402  child0->boundary_p1[1]=cbranch->center[1];
403  child0->boundary_p2[0]=cbranch->boundary_p2[0];
404  child0->boundary_p2[1]=cbranch->boundary_p2[1];
405 
406  child1->boundary_p1[0]=cbranch->boundary_p1[0];
407  child1->boundary_p1[1]=cbranch->center[1];
408  child1->boundary_p2[0]=cbranch->center[0];
409  child1->boundary_p2[1]=cbranch->boundary_p2[1];
410 
411  child2->boundary_p1[0]=cbranch->center[0];
412  child2->boundary_p1[1]=cbranch->boundary_p1[1];
413  child2->boundary_p2[0]=cbranch->boundary_p2[0];
414  child2->boundary_p2[1]=cbranch->center[1];
415 
416  child3->boundary_p1[0]=cbranch->boundary_p1[0];
417  child3->boundary_p1[1]=cbranch->boundary_p1[1];
418  child3->boundary_p2[0]=cbranch->center[0];
419  child3->boundary_p2[1]=cbranch->center[1];
420 
421 
422  // **** makes sure force does not require nbucket at leaf
423 
424  xcut = cbranch->center[0];
425  ycut = cbranch->center[1];
426 
427  // divide along y direction
428  for(i=0;i<NpInChildren;++i) x[i] = tree->xxp[particles[i]][1];
429  Utilities::quickPartition(ycut,&cuty,particles,x,NpInChildren);
430 
431  if(cuty > 0){
432  // divide first group in the x direction
433  for(i=0;i<cuty;++i) x[i] = tree->xxp[particles[i]][0];
434  Utilities::quickPartition(xcut,&cutx,particles,x,cuty);
435 
436  child3->nparticles=cutx;
437  assert(child3->nparticles >= 0);
438  if(child3->nparticles > 0)
439  child3->particles = particles;
440  else child3->particles = NULL;
441 
442  child2->nparticles = cuty - cutx;
443  assert(child2->nparticles >= 0);
444  if(child2->nparticles > 0)
445  child2->particles = &particles[cutx];
446  else child2->particles = NULL;
447 
448  }else{
449  child3->nparticles = 0;
450  child3->particles = NULL;
451 
452  child2->nparticles = 0;
453  child2->particles = NULL;
454  }
455 
456  if(cuty < NpInChildren){
457  // divide second group in the x direction
458  for(i=cuty;i<NpInChildren;++i) x[i-cuty] = tree->xxp[particles[i]][0];
459  Utilities::quickPartition(xcut,&cutx,&particles[cuty],x,NpInChildren-cuty);
460 
461  child1->nparticles=cutx;
462  assert(child1->nparticles >= 0);
463  if(child1->nparticles > 0)
464  child1->particles = &particles[cuty];
465  else child1->particles = NULL;
466 
467  child0->nparticles=NpInChildren - cuty - cutx;
468  assert(child0->nparticles >= 0);
469  if(child0->nparticles > 0)
470  child0->particles = &particles[cuty + cutx];
471  else child0->particles = NULL;
472 
473  }else{
474  child1->nparticles = 0;
475  child1->particles = NULL;
476 
477  child0->nparticles = 0;
478  child0->particles = NULL;
479  }
480 
481  //delete[] x;
482 
483  tree->moveToChild(0);
484  _BuildQTreeNB(child0->nparticles,child0->particles);
485  tree->moveUp();
486 
487  tree->moveToChild(1);
488  _BuildQTreeNB(child1->nparticles,child1->particles);
489  tree->moveUp();
490 
491  tree->moveToChild(2);
492  _BuildQTreeNB(child2->nparticles,child2->particles);
493  tree->moveUp();
494 
495  tree->moveToChild(3);
496  _BuildQTreeNB(child3->nparticles,child3->particles);
497  tree->moveUp();
498 
499  return;
500 }
501 
502 // calculates moments of the mass and the cutoff scale for each box
503 template<typename PType>
504 void TreeQuadHalos<PType>::CalcMoments(){
505 
506  //*** make compatable
507  IndexType i;
508  PosType rcom,xcm[2],xcut;
509  QBranch *cbranch;
510  PosType tmp;
511  double absmass; // absolute magnitude of mass
512 
513  tree->moveTop();
514  do{
515  cbranch=tree->current; /* pointer to current branch */
516 
517  cbranch->rmax = sqrt( pow(cbranch->boundary_p2[0]-cbranch->boundary_p1[0],2)
518  + pow(cbranch->boundary_p2[1]-cbranch->boundary_p1[1],2) );
519 
520  // calculate mass
521  for(i=0,cbranch->mass=0,absmass=0;i<cbranch->nparticles;++i){
522  cbranch->mass += xxp[cbranch->particles[i]*MultiMass].mass();
523  absmass += fabs(xxp[cbranch->particles[i]*MultiMass].mass());
524  }
525  // calculate center of mass
526  cbranch->center[0]=cbranch->center[1]=0;
527  for(i=0;i<cbranch->nparticles;++i){
528  tmp = fabs(xxp[cbranch->particles[i]*MultiMass].mass());
529 
530  //tmp = haloON ? halos[cbranch->particles[i]*MultiMass].mass : masses[cbranch->particles[i]*MultiMass];
531  cbranch->center[0] += tmp*tree->xxp[cbranch->particles[i]][0]/absmass;
532  cbranch->center[1] += tmp*tree->xxp[cbranch->particles[i]][1]/absmass;
533  }
535  // calculate quadropole moment of branch
537  cbranch->quad[0]=cbranch->quad[1]=cbranch->quad[2]=0;
538  for(i=0;i<cbranch->nparticles;++i){
539  xcm[0]=tree->xxp[cbranch->particles[i]][0]-cbranch->center[0];
540  xcm[1]=tree->xxp[cbranch->particles[i]][1]-cbranch->center[1];
541  xcut = pow(xcm[0],2) + pow(xcm[1],2);
542  tmp = xxp[cbranch->particles[i]*MultiMass].mass();
543 
544  cbranch->quad[0] += (xcut-2*xcm[0]*xcm[0])*tmp;
545  cbranch->quad[1] += (xcut-2*xcm[1]*xcm[1])*tmp;
546  cbranch->quad[2] += -2*xcm[0]*xcm[1]*tmp;
547  }
548 
549  // largest distance from center of mass of cell
550  for(i=0,rcom=0.0;i<2;++i) rcom += ( pow(cbranch->center[i]-cbranch->boundary_p1[i],2) > pow(cbranch->center[i]-cbranch->boundary_p2[i],2) ) ?
551  pow(cbranch->center[i]-cbranch->boundary_p1[i],2) : pow(cbranch->center[i]-cbranch->boundary_p2[i],2);
552 
553  rcom=sqrt(rcom);
554 
555  if(force_theta > 0.0){
556  cbranch->r2crit_angle = 1.15470*rcom/(force_theta);
557  cbranch->r2crit_angle *= cbranch->r2crit_angle;
558  }else cbranch->r2crit_angle=1.0e100;
559 
560  }while(tree->WalkStep(true));
561 
562  return;
563 }
564 
566 template<typename PType>
567 void TreeQuadHalos<PType>::rotate_coordinates(PosType **coord){
568  IndexType i;
569  short j;
570  PosType tmp[3];
571 
572  /* rotate particle positions */
573  for(i=0;i<tree->top->nparticles;++i){
574  for(j=0;j<3;++j) tmp[j]=0.0;
575  for(j=0;j<3;++j){
576  tmp[0]+=coord[0][j]*xxp[i][j];
577  tmp[1]+=coord[1][j]*xxp[i][j];
578  tmp[2]+=coord[2][j]*xxp[i][j];
579  }
580  for(j=0;j<3;++j) xxp[i][j]=tmp[j];
581  }
582  return;
583 }
584 
607 template<typename PType>
608 void TreeQuadHalos<PType>::force2D(const PosType *ray,PosType *alpha,KappaType *kappa,KappaType *gamma,KappaType *phi) const{
609 
610  alpha[0]=alpha[1]=gamma[0]=gamma[1]=gamma[2]=0.0;
611  *kappa=*phi=0.0;
612 
613  if(Nparticles == 0) return;
614 
615  assert(tree);
616  QBiterator<PType> it(tree);
617 
618  if(periodic_buffer){
619  PosType tmp_ray[2];
620 
621  for(int i = -1;i<2;++i){
622  for(int j = -1;j<2;++j){
623  tmp_ray[0] = ray[0] + i*original_xl;
624  tmp_ray[1] = ray[1] + j*original_yl;
625  walkTree_iter(it,tmp_ray,alpha,kappa,gamma,phi);
626  }
627  }
628  }else{
629  walkTree_iter(it,ray,alpha,kappa,gamma,phi);
630  }
631 
632  // Subtract off uniform mass sheet to compensate for the extra mass
633  // added to the universe in the halos.
634  //alpha[0] -= ray[0]*sigma_background*(inv_screening_scale2 == 0);
635  //alpha[1] -= ray[1]*sigma_background*(inv_screening_scale2 == 0);
636  //*kappa -= sigma_background;
637 
638  return;
639 }
640 
643 template<typename PType>
644 void TreeQuadHalos<PType>::neighbors(PosType ray[],PosType rmax,std::list<IndexType> &neighbors) const{
645  QBiterator<PType> it(tree);
646  neighbors.clear();
647 
648  PosType r2 = rmax*rmax;
649  bool decend=true;
650  it.movetop();
651  do{
652  int cut = Utilities::cutbox(ray,(*it)->boundary_p1,(*it)->boundary_p2,rmax);
653  decend = true;
654 
655  if(cut == 0){ // box outside circle
656  decend = false;
657  }else if(cut == 1){ // whole box inside circle
658  for(int i=0;i<(*it)->nparticles;++i) neighbors.push_back((*it)->particles[i]);
659  decend = false;
660  }else if((*it)->child0 == NULL){ // at leaf
661  for(int i=0;i<(*it)->nparticles;++i){
662  if( (tree->xxp[(*it)->particles[i]][0] - ray[0])*(tree->xxp[(*it)->particles[i]][0] - ray[0])
663  + (tree->xxp[(*it)->particles[i]][1] - ray[1])*(tree->xxp[(*it)->particles[i]][1] - ray[1]) < r2) neighbors.push_back((*it)->particles[i]);
664  }
665  decend = false;
666  }
667  }while(it.TreeWalkStep(decend));
668 
669 }
670 
671 template<typename PType>
672 void TreeQuadHalos<PType>::walkTree_iter(
673  QBiterator<PType> &treeit,
674  const PosType *ray
675  ,PosType *alpha
676  ,KappaType *kappa
677  ,KappaType *gamma
678  ,KappaType *phi
679  ) const {
680 
681  PosType xcm[2],rcm2cell,rcm2,tmp,boxarea;
682  IndexType i;
683  bool allowDescent=true;
684  unsigned long count=0,tmp_index;
685  PosType prefac;
686  //bool notscreen = true;
687 
688  assert(tree);
689 
690  treeit.movetop();
691  //tree->moveTop();
692  do{
693  ++count;
694  allowDescent=false;
695 
696  /*
697  if(inv_screening_scale2 != 0) notscreen = BoxIntersectCircle(ray,3*sqrt(1.0/inv_screening_scale2)
698  ,(*treeit)->boundary_p1
699  ,(*treeit)->boundary_p2);
700  else notscreen = true;
701  */
702  //if(tree->current->nparticles > 0)
703  if((*treeit)->nparticles > 0)
704  {
705 
706  xcm[0]=(*treeit)->center[0]-ray[0];
707  xcm[1]=(*treeit)->center[1]-ray[1];
708 
709  rcm2cell = xcm[0]*xcm[0] + xcm[1]*xcm[1];
710 
711  boxarea = ((*treeit)->boundary_p2[0]-(*treeit)->boundary_p1[0])*((*treeit)->boundary_p2[0]-(*treeit)->boundary_p1[0]);
712 
713  if( rcm2cell < ((*treeit)->rcrit_angle)*((*treeit)->rcrit_angle) || rcm2cell < 5.83*boxarea)
714  {
715 
716  // includes rcrit_particle constraint
717  allowDescent=true;
718 
719 
720  // Treat all particles in a leaf as a point particle
721  if(treeit.atLeaf())
722  {
723 
724  for(i = 0 ; i < (*treeit)->nparticles ; ++i)
725  {
726 
727  xcm[0] = tree->xp[(*treeit)->particles[i]][0] - ray[0];
728  xcm[1] = tree->xp[(*treeit)->particles[i]][1] - ray[1];
729 
730  rcm2 = xcm[0]*xcm[0] + xcm[1]*xcm[1];
731  double screening = exp(-rcm2*inv_screening_scale2);
732  if(rcm2 < 1e-20) rcm2 = 1e-20;
733 
734  //tmp_index = MultiMass*(*treeit)->particles[i];
735 
736  //if(haloON){ prefac = halos[tmp_index]->get_mass(); }
737  //else{ prefac = masses[tmp_index]; }
738 
739  double mass;
740  if(MultiMass) mass = xxp[(*treeit)->particles[i]].Mass();
741  else mass = mass_fixed;
742 
743  prefac = mass/rcm2/PI*screening;
744  tmp = -mass*( screening/rcm2/PI - inv_area);
745 
746  alpha[0] += tmp*xcm[0];
747  alpha[1] += tmp*xcm[1];
748 
749  tmp = -2.0*prefac/rcm2;
750 
751  gamma[0] += 0.5*(xcm[0]*xcm[0]-xcm[1]*xcm[1])*tmp;
752  gamma[1] += xcm[0]*xcm[1]*tmp;
753 
754  *kappa -= mass*inv_area;
755  *phi += (prefac*log(rcm2) - mass*inv_area)*rcm2*0.5;
756 
757  } // end of for
758  } // end of if(tree->atLeaf())
759 
760 
761  // Find the particles that intersect with ray and add them individually.
762  if(rcm2cell < 5.83*boxarea)
763  {
764 
765  for(i = 0 ; i < (*treeit)->Nbig_particles ; ++i)
766  {
767 
768  tmp_index = (*treeit)->big_particles[i];
769 
770  xcm[0] = tree->xp[tmp_index][0] - ray[0];
771  xcm[1] = tree->xp[tmp_index][1] - ray[1];
772 
773  rcm2 = xcm[0]*xcm[0] + xcm[1]*xcm[1];
774 
775  if(rcm2 < 1e-20) rcm2 = 1e-20;
776 
777  //PosType size = sizes[tmp_index*MultiRadius];
778  PosType size = xxp[tmp_index*MultiRadius].size();
779 
780  // intersecting, subtract the point particle
781  if(rcm2 < 4*size*size)
782  {
783  b_spline_profile(xcm,sqrt(rcm2),xxp[MultiMass*tmp_index].Mass(),size,alpha,kappa,gamma,phi);
784  }
785 
786  }
787  }
788 
789  }
790  else
791  { // use whole cell
792 
793  allowDescent=false;
794 
795  // monopole
796  double screening = exp(-rcm2cell*inv_screening_scale2);
797  //double tmp = -1.0*(*treeit)->mass/rcm2cell/PI*screening;
798  double tmp = (*treeit)->mass*( inv_area - screening/rcm2cell/PI );
799 
800  alpha[0] += tmp*xcm[0];
801  alpha[1] += tmp*xcm[1];
802 
803  { // taken out to speed up
804  tmp=-2.0*(*treeit)->mass/PI/rcm2cell/rcm2cell*screening;
805  gamma[0] += 0.5*(xcm[0]*xcm[0]-xcm[1]*xcm[1])*tmp;
806  gamma[1] += xcm[0]*xcm[1]*tmp;
807 
808  *phi += 0.5*(*treeit)->mass*log( rcm2cell )/PI*screening;
809  *phi -= 0.5*( (*treeit)->quad[0]*xcm[0]*xcm[0] + (*treeit)->quad[1]*xcm[1]*xcm[1] + 2*(*treeit)->quad[2]*xcm[0]*xcm[1] )/(PI*rcm2cell*rcm2cell)*screening;
810  }
811 
812  // quadrapole contribution
813  // the kappa and gamma are not calculated to this order
814  alpha[0] -= ((*treeit)->quad[0]*xcm[0] + (*treeit)->quad[2]*xcm[1])
815  /(rcm2cell*rcm2cell)/PI*screening;
816  alpha[1] -= ((*treeit)->quad[1]*xcm[1] + (*treeit)->quad[2]*xcm[0])
817  /(rcm2cell*rcm2cell)/PI*screening;
818 
819  tmp = 4*((*treeit)->quad[0]*xcm[0]*xcm[0] + (*treeit)->quad[1]*xcm[1]*xcm[1]
820  + 2*(*treeit)->quad[2]*xcm[0]*xcm[1])/(rcm2cell*rcm2cell*rcm2cell)/PI*screening;
821 
822  alpha[0] += tmp*xcm[0];
823  alpha[1] += tmp*xcm[1];
824 
825  *kappa -= (*treeit)->mass * inv_area;
826  *phi -= 0.5*(*treeit)->mass*inv_area;
827  }
828  }
829  }while(treeit.TreeWalkStep(allowDescent));
830 
831 
832  // Subtract off uniform mass sheet to compensate for the extra mass
833  // added to the universe in the halos.
834  //alpha[0] -= ray[0]*sigma_background*(inv_screening_scale2 == 0.0);
835  //alpha[1] -= ray[1]*sigma_background*(inv_screening_scale2 == 0.0);
836  //{ // taken out to speed up
837  // *kappa -= sigma_background;
838  // // *phi -= sigma_background * sigma_background ;
839  //}
840 
841  return;
842 }
843 
844 
869 template<typename PType>
870 void TreeQuadHalos<PType>::force2D_recur(const PosType *ray,PosType *alpha,KappaType *kappa,KappaType *gamma,KappaType *phi){
871 
872 
873  alpha[0]=alpha[1]=gamma[0]=gamma[1]=gamma[2]=0.0;
874  *kappa=*phi=0.0;
875 
876  if(Nparticles == 0) return;
877  assert(tree);
878 
879  //walkTree_recur(tree->top,ray,&alpha[0],kappa,&gamma[0],phi);
880 
881  if(periodic_buffer){
882  PosType tmp_ray[2],tmp_c[2];
883 
884  // for points outside of primary zone shift to primary zone
885  //if(inbox(ray,tree->top->boundary_p1,tree->top->boundary_p2)){
886  tmp_c[0] = ray[0];
887  tmp_c[1] = ray[1];
888  /*}else{
889  int di[2];
890  PosType center[2] = { (tree->top->boundary_p1[0] + tree->top->boundary_p2[0])/2 ,
891  (tree->top->boundary_p1[1] + tree->top->boundary_p2[1])/2 };
892 
893  di[0] = (int)( 2*(ray[0] - center[0])/lx );
894  di[1] = (int)( 2*(ray[1] - center[1])/ly );
895 
896  di[0] = (di[0] > 0) ? (di[0] + 1)/2 : (di[0] - 1)/2;
897  di[1] = (di[1] > 0) ? (di[1] + 1)/2 : (di[1] - 1)/2;
898 
899  tmp_c[0] = ray[0] - lx * di[0];
900  tmp_c[1] = ray[1] - ly * di[1];
901 
902  assert(inbox(tmp_c,tree->top->boundary_p1,tree->top->boundary_p2));
903  }*/
904 
905  for(int i = -1;i<2;++i){
906  for(int j = -1;j<2;++j){
907  tmp_ray[0] = tmp_c[0] + i*original_xl;
908  tmp_ray[1] = tmp_c[1] + j*original_yl;
909  walkTree_recur(tree->top,tmp_ray,alpha,kappa,gamma,phi);
910  }
911  }
912  }else{
913  walkTree_recur(tree->top,ray,alpha,kappa,gamma,phi);
914  }
915 
916  // Subtract off uniform mass sheet to compensate for the extra mass
917  // added to the universe in the halos.
918  //alpha[0] -= ray[0]*sigma_background*(inv_screening_scale2 == 0);
919  //alpha[1] -= ray[1]*sigma_background*(inv_screening_scale2 == 0);
920 
921  //*kappa -= sigma_background;
922 
923  return;
924 }
925 
926 
927 template<typename PType>
928 void TreeQuadHalos<PType>::walkTree_recur(QBranch *branch,PosType const *ray,PosType *alpha,KappaType *kappa,KappaType *gamma, KappaType *phi){
929 
930  PosType xcm[2],rcm2cell,rcm2,tmp,boxsize2;
931  IndexType i;
932  std::size_t tmp_index;
933  PosType prefac;
934  /*bool notscreen = true;
935 
936  if(inv_screening_scale2 != 0) notscreen = BoxIntersectCircle(ray,3*sqrt(1.0/inv_screening_scale2), branch->boundary_p1, branch->boundary_p2);
937  */
938  if(branch->nparticles > 0)
939  {
940  xcm[0]=branch->center[0]-ray[0];
941  xcm[1]=branch->center[1]-ray[1];
942 
943  rcm2cell = xcm[0]*xcm[0] + xcm[1]*xcm[1];
944 
945  boxsize2 = (branch->boundary_p2[0]-branch->boundary_p1[0])*(branch->boundary_p2[0]-branch->boundary_p1[0]);
946 
947  if( rcm2cell < branch->r2crit_angle || rcm2cell < 5.83*boxsize2)
948  {
949 
950  // Treat all particles in a leaf as a point particle
951  if(tree->atLeaf(branch))
952  {
953 
954  for(i = 0 ; i < branch->nparticles ; ++i){
955 
956  xcm[0] = tree->xxp[branch->particles[i]][0] - ray[0];
957  xcm[1] = tree->xxp[branch->particles[i]][1] - ray[1];
958 
959  rcm2 = xcm[0]*xcm[0] + xcm[1]*xcm[1];
960  double screening = exp(-rcm2*inv_screening_scale2);
961  if(rcm2 < 1e-20) rcm2 = 1e-20;
962 
963  tmp_index = MultiMass*branch->particles[i];
964 
965  //if(haloON ) { prefac = halos[tmp_index]->get_mass(); }
966  //else{ prefac = masses[tmp_index]; }
967  double mass = xxp[tmp_index].mass();
968  prefac = mass*screening/rcm2/PI;
969  //prefac /= rcm2*PI/screening;
970 
971 
972  alpha[0] += (inv_area*mass - prefac)*xcm[0];
973  alpha[1] += (inv_area*mass - prefac)*xcm[1];
974 
975  {
976  tmp = -2.0*prefac/rcm2;
977 
978  gamma[0] += 0.5*(xcm[0]*xcm[0]-xcm[1]*xcm[1])*tmp;
979  gamma[1] += xcm[0]*xcm[1]*tmp;
980 
981  *kappa -= mass*inv_area;
982  *phi += prefac*rcm2*0.5*log(rcm2) - mass*inv_area*rcm2/2;
983  }
984  }
985  }
986 
987  // Find the particles that intersect with ray and add them individually.
988  if(rcm2cell < 5.83*boxsize2)
989  {
990 
991  for(i = 0 ; i < branch->Nbig_particles ; ++i)
992  {
993 
994  tmp_index = branch->big_particles[i];
995 
996  xcm[0] = tree->xxp[tmp_index][0] - ray[0];
997  xcm[1] = tree->xxp[tmp_index][1] - ray[1];
998  rcm2 = xcm[0]*xcm[0] + xcm[1]*xcm[1];
999 
1001  if(rcm2 < 1e-20) rcm2 = 1e-20;
1002 
1003  PosType size = xxp[tmp_index*MultiRadius].size();
1004 
1005  // intersecting, subtract the point particle
1006  if(rcm2 < 4*size*size)
1007  {
1008 
1009  b_spline_profile(xcm,sqrt(rcm2),xxp[MultiMass*tmp_index].mass(),size,alpha,kappa,gamma,phi);
1010 
1011  }
1012 
1013  }
1014  }
1015 
1016  if(branch->child0 != NULL)
1017  walkTree_recur(branch->child0,&ray[0],&alpha[0],kappa,&gamma[0],&phi[0]);
1018  if(branch->child1 != NULL)
1019  walkTree_recur(branch->child1,&ray[0],&alpha[0],kappa,&gamma[0],&phi[0]);
1020  if(branch->child2 != NULL)
1021  walkTree_recur(branch->child2,&ray[0],&alpha[0],kappa,&gamma[0],&phi[0]);
1022  if(branch->child3 != NULL)
1023  walkTree_recur(branch->child3,&ray[0],&alpha[0],kappa,&gamma[0],&phi[0]);
1024 
1025  }
1026  else
1027  { // use whole cell
1028 
1029  double screening = exp(-rcm2cell*inv_screening_scale2);
1030  double tmp = branch->mass*inv_area - branch->mass/rcm2cell/PI*screening;
1031 
1032  alpha[0] += tmp*xcm[0];
1033  alpha[1] += tmp*xcm[1];
1034 
1035  { // taken out to speed up
1036  tmp=-2.0*branch->mass/PI/rcm2cell/rcm2cell*screening;
1037  gamma[0] += 0.5*(xcm[0]*xcm[0]-xcm[1]*xcm[1])*tmp;
1038  gamma[1] += xcm[0]*xcm[1]*tmp;
1039 
1040  *phi += 0.5*tree->current->mass*log( rcm2cell )/PI*screening;
1041  *phi -= 0.5*( tree->current->quad[0]*xcm[0]*xcm[0] + tree->current->quad[1]*xcm[1]*xcm[1] + 2*tree->current->quad[2]*xcm[0]*xcm[1] )/(PI*rcm2cell*rcm2cell)*screening;
1042  }
1043 
1044  // quadrapole contribution
1045  // the kappa and gamma are not calculated to this order
1046  alpha[0] -= (branch->quad[0]*xcm[0] + branch->quad[2]*xcm[1])
1047  /(rcm2cell*rcm2cell)/PI*screening;
1048  alpha[1] -= (branch->quad[1]*xcm[1] + branch->quad[2]*xcm[0])
1049  /(rcm2cell*rcm2cell)/PI*screening;
1050 
1051  tmp = 4*(branch->quad[0]*xcm[0]*xcm[0] + branch->quad[1]*xcm[1]*xcm[1]
1052  + 2*branch->quad[2]*xcm[0]*xcm[1])/(rcm2cell*rcm2cell*rcm2cell)/PI*screening;
1053 
1054  alpha[0] += tmp*xcm[0];
1055  alpha[1] += tmp*xcm[1];
1056 
1057  *kappa -= inv_area * branch->mass;
1058  *phi -= 0.5*branch->mass*rcm2cell*inv_area;
1059 
1060  return;
1061  }
1062 
1063  }
1064  return;
1065 }
1066 
1070 template<typename PType>
1071 void TreeQuadHalos<PType>::printParticlesInBranch(unsigned long number){
1072  unsigned long i;
1073 
1074  tree->moveTop();
1075  do{
1076  if(tree->current->number == number){
1077  std::cout << tree->current->nparticles << std::endl;
1078  for(i=0;i<tree->current->nparticles;++i){
1079  std::cout << xxp[tree->current->particles[i]][0] << " "
1080  << xxp[tree->current->particles[i]][1] << std::endl;
1081  }
1082  return;
1083  }
1084  }while(tree->WalkStep(true));
1085 
1086  return;
1087 }
1092 template<typename PType>
1093 void TreeQuadHalos<PType>::printBranchs(int level){
1094 
1095  bool decend = true;
1096  tree->moveTop();
1097  do{
1098  std::cout << tree->current->boundary_p1[0] << " " << tree->current->boundary_p1[1] << " "
1099  << tree->current->boundary_p2[0] << " " << tree->current->boundary_p2[1] << std::endl;
1100  if(tree->current->level == level) decend = false;
1101  else decend = true;
1102  }while(tree->WalkStep(decend));
1103 
1104  return;
1105 };
1106 
1107 #endif /* QUAD_TREE_H_ */
TreeQuadHalos::TreeQuadHalos
TreeQuadHalos(PType *xpt, IndexType Npoints, PosType my_inv_area=0, int bucket=5, PosType theta_force=0.1, PosType my_inv_screening_scale=0)
Constructor meant for point particles, simulation particles.
Definition: quadTreeHalo.h:189
TreeQuadHalos
TreeQuadHalos is a class for calculating the deflection, kappa and gamma by tree method.
Definition: quadTreeHalo.h:44
TreeQuadHalos::neighbors
void neighbors(PosType ray[], PosType rmax, std::list< IndexType > &neighbors) const
find all points within rmax of ray in 2D
TreeQuadHalos::rotate_coordinates
void rotate_coordinates(PosType **coord)
simple rotates the coordinates in the xp array
Definition: quadTreeHalo.h:567
TreeQuadHalos::periodic_buffer
bool periodic_buffer
if true there is one layer of peridic buffering
Definition: quadTreeHalo.h:86
QBiterator
A iterator class fore TreeStruct that allows for movement through the tree without changing anything ...
Definition: qTreeNB.h:173
TreeQuadHalos::~TreeQuadHalos
~TreeQuadHalos()
Particle positions and other data are not destroyed.
Definition: quadTreeHalo.h:228
TreeQuadHalos::_BuildQTreeNB
void _BuildQTreeNB(IndexType nparticles, IndexType *particles)
tree must be created and first branch must be set before start
Definition: quadTreeHalo.h:297
TreeQuadHalos::WhichQuad
short WhichQuad(PosType *x, QBranch &branch)
returns an index for which of the four quadrangles of the branch the point x[] is in
Definition: quadTreeHalo.h:291
QTreeNB
QTreeNB: Tree structure used for force calculation with particles (i.e. stars, Nbody or halos).
Definition: qTreeNB.h:98
Generated on Sun Sep 27 2026 14:12:56 for GLAMERDOC++ by   doxygen 1.8.18