GLAMERDOC++
Gravitational Lensing Code Library
oTreeNB.h
1 //
2 // OTreeNB.h
3 // 3d octal tree
4 // GLAMER
5 //
6 // Created by Ben Metcalf on 5/2025.
7 //
8 
9 #ifndef OTreeNB_h
10 #define OTreeNB_h
11 
12 #include "utilities.h"
13 #include "cosmo.h"
14 #include "point.h"
15 
16 #ifndef PI
17 #define PI 3.141593
18 #endif
19 
22 #ifndef IndexType_declare
23 #define IndexType_declare
24 typedef unsigned long IndexType;
25 #endif
26 
27 
40 template <typename PType = Point_3d<PosType>,typename DType = PosType >
41 struct OTreeNB
42 {
43 
45 
47 struct Branch
48 {
49  Branch()
50  {
51  prev = nullptr;
52  child = nullptr;
53  brother = nullptr;
54  particles = nullptr;
55  nparticles = 0;
56  level = -1;
57  boxsize = 0.0;
58  root_inv_density = 0.0;
59  //++Nbranches;
60  };
61 
62  // static int Nbranches; ///< number of branches in tree, used for debugging
63  //Branch(Branch *parent) : prev(parent)
64  //{
65  // child = nullptr;
66  // brother = nullptr;
67  // particles = nullptr;
68  // nparticles = 0;
69  // level = parent->level + 1;
70  //};
71 
72  ~Branch() {
73  //std::cout << " deleting branch " << Nbranches << " ";
74  //--Nbranches;
75  if(children != nullptr) delete[] children;
76  };//{std::cout << "Branch destructor called" << std::endl; };
77 
79  IndexType *particles;
80  IndexType nparticles;
81 
83  PType boundary_p1;
85  PType boundary_p2;
87  PType center;
88 
89  PosType boxsize; // one dimensional size of box
90  PosType root_inv_density; // cube root of number density in branch
92  int level;
93  //unsigned long number;
94 
95  // std::vector<Branch> children;
96  Branch *children = nullptr;
97 
98  Branch *child;
100  Branch *prev;
103  Branch *brother;
104 
105  // projected quantities
106  // center of mass
107  PosType xcm[3] = {0.0, 0.0, 0.0};
108  // quadropole moment of branch in particle mass units
109  PosType quad[3] = {0.0, 0.0, 0.0};
110  DType mass = 0;
111 };
112 
117  class iterator
118  {
119 
120  private:
121  Branch *current;
122 
123  public:
125  iterator(OTreeNB<PType> *tree) { current = tree->top; }
127  iterator(Branch *branch) { current = branch; }
128 
129  iterator& operator=(iterator &it)
130  {
131  current = it.current;
132  return *this;
133  }
134 
135  iterator& operator=(const Branch *branch)
136  {
137  current = branch;
138  return *this;
139  }
140 
142  Branch *operator*() { return current; }
143 
144  bool up()
145  {
146  if (current->prev == nullptr)
147  {
148  current = nullptr;
149  return false; // at top
150  }
151  current = current->prev;
152  return true;
153  }
154 
156  bool down(short child)
157  {
158  assert(child >= 0 && child < 8);
159  current = current->child;
160  }
161 
162  bool atTop() const
163  {
164  if (current->prev == nullptr)
165  return true;
166  return false;
167  }
168 
169  bool atLeaf() const
170  {
171  for (int i = 0; i < 8; ++i)
172  {
173  if (current->child != nullptr)
174  return false;
175  }
176  return true;
177  }
178 
186  bool walk(bool allowDescent,iterator end)
187  {
188  if (current == nullptr)
189  return false; // at top
190 
191  // if there is a valid chaild, go down
192  if (allowDescent && current->child != nullptr){
193  current = current->child;
194  return true;
195  }
196 
197  if (current->brother != end.current->brother)
198  {
199  current = current->brother;
200  return true;
201  }
202  current = end.current;
203  return false; // at end
204  }
205 
206  };
207 
208 
209  OTreeNB(PType *xp,IndexType nparticles):
210  xxp(xp)
211  {
212  index.resize(nparticles);
213  for(unsigned long i=0;i<nparticles;++i) index[i] = i;
214  Nbranches = 1;
215  top.particles = index.data();
216  //std::cout << top.particles[10] << std::endl;
217  top.nparticles = nparticles;
218 
219  top.level = 0;
220  top.boundary_p1[0] = top.boundary_p2[0] = xp[0][0];
221  top.boundary_p1[1] = top.boundary_p2[1] = xp[0][1];
222  top.boundary_p1[2] = top.boundary_p2[2] = xp[0][2];
223 
224  for(IndexType i = 1 ; i < nparticles ; ++i)
225  {
226  for(int j=0;j<3;++j){
227  top.boundary_p1[j] = top.boundary_p1[j] < xp[i][j] ? top.boundary_p1[j] : xp[i][j];
228  top.boundary_p2[j] = top.boundary_p2[j] > xp[i][j] ? top.boundary_p2[j] : xp[i][j];
229 
230  //top.boundary_p1[j] = min<PType>(top.boundary_p1[j], xp[i][j]);
231  //top.boundary_p2[j] = max<PType>(top.boundary_p2[j], xp[i][j]);
232  }
233  }
234 
235 
236  // make it into a cube
237  PType length = top.boundary_p2 - top.boundary_p1;
238  PosType max_length = length[0];
239  int dim = 0;
240  for(int j=1;j<3;++j){
241  if(length[j] > max_length){
242  max_length = length[j];
243  dim = j;
244  }
245  }
246  for(int j=0;j<3;++j){
247  if(j != dim){
248  top.boundary_p2[j] = top.boundary_p1[j] + max_length;
249  }
250  }
251 
252  top.center = (top.boundary_p1 + top.boundary_p2) * 0.5;
253  top.boxsize = max_length;
254  top.root_inv_density = max_length/pow(top.nparticles,1.0/3.);
255  }
256  ~OTreeNB(){ }; // all the Branches are deleted in the destructor of Branch
257 
258 
259  // returns number of branches in tree
260  unsigned long size() { return Nbranches; };
261 
262  // build the tree down to bucket size, does not require the particles to have sizes
263  void build(int Nbucket){
264  auto it = begin();
265  if((*it)->child != nullptr){
266  throw std::runtime_error("OTreeNB Error: calling build() on already built tree");
267  }
268  span8(*it);
269 
270  while(it.walk(true,begin())){
271 
272  //std::cout << "level : " << (*it)->level << std::endl;
273  //std::cout << (*it)->boundary_p1 << std::endl;
274  //std::cout << (*it)->boundary_p2 << std::endl;
275  //std::cout << "N : " << (*it)->nparticles << std::endl << std::endl;
276 
277  if((*it)->nparticles > Nbucket) span8(*it);
278  }
279  calcMoments(); // assumes all particles are same mass
280  }
281 
282  // build the tree down to bucket size, does not require the particles to have sizes
283  void build_size(int Nbucket){
284  auto it = begin();
285  if((*it)->child != nullptr){
286  throw std::runtime_error("OTreeNB Error: calling build() on already built tree");
287  }
288  span8(*it);
289  while(it.walk(true,begin())){
290 
291  // find largest particle in branch
292  double max_size = 0.0;
293  for(IndexType i=0;i<(*it)->nparticles;++i){
294  max_size = max(max_size, xxp[(*it)->particles[i]].size);
295  }
296  //std::cout << "level : " << (*it)->level << std::endl;
297  //std::cout << (*it)->boundary_p1 << std::endl;
298  //std::cout << (*it)->boundary_p2 << std::endl;
299  //std::cout << "N : " << (*it)->nparticles << std::endl << std::endl;
300  if((*it)->nparticles > Nbucket && max_size < (*it)->boxsize ) span8(*it);
301  }
302  }
303 
304  // calculate the moments assuming the particles are equal mass and symmetric
305  void calcMoments(){
306  auto it = begin();
307  while(it.walk(true,begin())){
308 
309  Branch *branch = *it;
310  // calculate center of mass
311  branch->xcm[0] = 0.0;
312  branch->xcm[1] = 0.0;
313  branch->xcm[2] = 0.0;
314  if(branch->nparticles == 0){
315  branch->mass = 0.0;
316  }else{
317 
318  for(IndexType i=0;i<branch->nparticles;++i){
319  PType &x = xxp[branch->particles[i]];
320  branch->xcm[0] += x[0];
321  branch->xcm[1] += x[1];
322  branch->xcm[2] += x[2];
323  }
324  branch->xcm[0] /= branch->nparticles;
325  branch->xcm[1] /= branch->nparticles;
326  branch->xcm[2] /= branch->nparticles;
327  branch->mass = branch->nparticles;
328  }
329 
343  }
344  }
345 
346  // calculate the moments assuming the particles have different masses
347  void calcMoments(DType *masses){
348  moments_set = true;
349  auto it = begin();
350  while(it.walk(true,begin())){
351 
352  Branch *branch = *it;
353  if(branch->nparticles == 0) continue;
354 
355  // calculate center of mass
356  branch->xcm[0] = 0.0;
357  branch->xcm[1] = 0.0;
358  branch->xcm[2] = 0.0;
359  branch->mass = 0;
360 
361  for(IndexType i=0;i<branch->nparticles;++i){
362  IndexType j = branch->particles[i];
363  PType &x = xxp[j];
364  double m = masses[j];
365  branch->xcm[0] += x[0]*m;
366  branch->xcm[1] += x[1]*m;
367  branch->xcm[2] += x[2]*m;
368  branch->mass += m;
369  }
370  branch->xcm[0] /= branch->mass;
371  branch->xcm[1] /= branch->mass;
372  branch->xcm[2] /= branch->mass;
373 
389  }
390  }
391 // calculate the moments from halos
392  template<typename HType>
393  void calcMoments(std::vector<HType> &halos){
394  moments_set = true;
395  auto it = begin();
396  while(it.walk(true,begin())){
397 
398  Branch *branch = *it;
399  if(branch->nparticles == 0) continue;
400 
401  // calculate center of mass
402  branch->xcm[0] = 0.0;
403  branch->xcm[1] = 0.0;
404  branch->xcm[2] = 0.0;
405  branch->mass = 0;
406 
407  for(IndexType i=0;i<branch->nparticles;++i){
408  IndexType j = branch->particles[i];
409  PType &x = xxp[j];
410  double m = halos[j].get_mass();
411 
412  branch->xcm[0] += x[0]*m;
413  branch->xcm[1] += x[1]*m;
414  branch->xcm[2] += x[2]*m;
415  branch->mass += m;
416  }
417  branch->xcm[0] /= branch->mass;
418  branch->xcm[1] /= branch->mass;
419  branch->xcm[2] /= branch->mass;
420 
437  }
438  }
439 
440  void force2D(const PosType *ray
441  ,PosType particle_mass
442  ,PosType smooth_factor
443  ,PosType theta2 // opening angle in radians
444  ,PosType inv_area // compensating negative mass area
445  ,PosType *alpha // not zeroed
446  ,KappaType *kappa // not zeroed
447  ,KappaType *gamma // not zeroed
448  ,KappaType *phi // not zeroed
449  ) {
450 
451  //alpha[0]=alpha[1]=gamma[0]=gamma[1]=gamma[2]=0.0;
452  //*kappa=*phi=0.0;
453 
454  PosType inv_area_pi = inv_area*PI;
455  PosType xcm[2],r2cm;
456  auto it = begin();
457  bool decend = true;
458  while(it.walk(decend,begin())){
459 
460  xcm[0] = (*it)->xcm[0] - ray[0];
461  xcm[1] = (*it)->xcm[1] - ray[1];
462  r2cm = xcm[0]*xcm[0] + xcm[1]*xcm[1];
463 
464  if( r2cm*theta2 > ((*it)->boxsize)*((*it)->boxsize) ){
465 
466  if(r2cm*inv_area_pi < 1){
467  //????? use moments
468  double mass = (*it)->mass * particle_mass;
469  double prefac = mass/r2cm/PI;
470  double tmp = -( prefac - mass*inv_area);
471 
472  alpha[0] -= tmp*xcm[0];
473  alpha[1] -= tmp*xcm[1];
474 
475  tmp = -2.0*prefac/r2cm;
476 
477  gamma[0] += 0.5*(xcm[0]*xcm[0]-xcm[1]*xcm[1])*tmp;
478  gamma[1] += xcm[0]*xcm[1]*tmp;
479 
480  *kappa -= mass*inv_area;
481  *phi += (prefac*log(r2cm) - mass*inv_area)*r2cm*0.5;
482  }
483  decend = false;
484  }else{
485 
486  if((*it)->child == nullptr){ // leaf
487 
488  for(IndexType i=0;i<(*it)->nparticles;++i){
489 
490  PType &x = xxp[(*it)->particles[i]];
491  xcm[0] = x[0] - ray[0];
492  xcm[1] = x[1] - ray[1];
493  r2cm = xcm[0]*xcm[0] + xcm[1]*xcm[1];
494  if(r2cm*inv_area_pi < 1){
495  double prefac = particle_mass /r2cm/PI;
496  double tmp = -( prefac - particle_mass*inv_area);
497 
498  alpha[0] -= tmp*xcm[0];
499  alpha[1] -= tmp*xcm[1];
500 
501  tmp = -2.0*prefac/r2cm;
502 
503  gamma[0] += 0.5*(xcm[0]*xcm[0]-xcm[1]*xcm[1])*tmp;
504  gamma[1] += xcm[0]*xcm[1]*tmp;
505 
506  *kappa -= particle_mass*inv_area;
507  *phi += (prefac*log(r2cm)- particle_mass*inv_area)*r2cm*0.5;
508 
509  double scale = smooth_factor*(*it)->root_inv_density;
510  if(r2cm < 4*scale*scale){
511  // cubic B-spline profile
512  b_spline_profile(
513  xcm
514  ,sqrt(r2cm)
515  ,particle_mass
516  ,scale
517  ,alpha,kappa,gamma,phi
518  );
519  }
520  }
521  }
522  decend = false;
523  }else{
524  // decend
525  decend = true;
526  }
527 
528  }
529  }
530  }
531  void force2D(const PosType *ray
532  ,DType *masses
533  ,PosType smooth_factor
534  ,PosType theta2 // opening angle in radians
535  ,PosType inv_area // compensating negative mass area
536  ,PosType *alpha // not zeroed
537  ,KappaType *kappa // not zeroed
538  ,KappaType *gamma // not zeroed
539  ,KappaType *phi // not zeroed
540  ) {
541 
542  assert( moments_set );
543  //alpha[0]=alpha[1]=gamma[0]=gamma[1]=gamma[2]=0.0;
544  //*kappa=*phi=0.0;
545  PosType inv_area_pi = inv_area*PI;
546  PosType xcm[2],r2cm;
547  auto it = begin();
548  bool decend = true;
549  while(it.walk(decend,begin())){
550 
551  xcm[0] = (*it)->xcm[0] - ray[0];
552  xcm[1] = (*it)->xcm[1] - ray[1];
553  r2cm = xcm[0]*xcm[0] + xcm[1]*xcm[1];
554 
555  if( r2cm*theta2 > ((*it)->boxsize)*((*it)->boxsize) ){
556 
557  if(r2cm*inv_area_pi < 1){
558  //????? use moments
559  DType mass = (*it)->mass;
560  double prefac = mass/r2cm/PI;
561  double tmp = -( prefac - mass*inv_area);
562 
563  alpha[0] -= tmp*xcm[0];
564  alpha[1] -= tmp*xcm[1];
565 
566  tmp = -2.0*prefac/r2cm;
567 
568  gamma[0] += 0.5*(xcm[0]*xcm[0]-xcm[1]*xcm[1])*tmp;
569  gamma[1] += xcm[0]*xcm[1]*tmp;
570 
571  *kappa -= mass*inv_area;
572  *phi += (prefac*log(r2cm) - mass*inv_area)*r2cm*0.5;
573  }
574  decend = false;
575  }else{
576 
577  if((*it)->child == nullptr){ // leaf
578 
579  for(IndexType i=0;i<(*it)->nparticles;++i){
580  IndexType j = (*it)->particles[i];
581  PType &x = xxp[j];
582  DType mass = masses[j];
583 
584  xcm[0] = x[0] - ray[0];
585  xcm[1] = x[1] - ray[1];
586  r2cm = xcm[0]*xcm[0] + xcm[1]*xcm[1];
587  if(r2cm*inv_area_pi < 1){
588  double prefac = mass /r2cm/PI;
589  double tmp = -( prefac - mass * inv_area);
590 
591  alpha[0] -= tmp*xcm[0];
592  alpha[1] -= tmp*xcm[1];
593 
594  tmp = -2.0*prefac/r2cm;
595 
596  gamma[0] += 0.5*(xcm[0]*xcm[0]-xcm[1]*xcm[1])*tmp;
597  gamma[1] += xcm[0]*xcm[1]*tmp;
598 
599  *kappa -= mass * inv_area;
600  *phi += (prefac*log(r2cm)- mass * inv_area)*r2cm*0.5;
601 
602  double scale = smooth_factor*(*it)->root_inv_density;
603  if(r2cm < 4*scale*scale){
604  // cubic B-spline profile
605  b_spline_profile(
606  xcm
607  ,sqrt(r2cm)
608  ,mass
609  ,scale
610  ,alpha,kappa,gamma,phi
611  );
612  }
613  }
614  }
615  decend = false;
616  }else{
617  // decend
618  decend = true;
619  }
620 
621  }
622  }
623  }
624 
625  void force2D_hole(const PosType *ray
626  ,PosType particle_mass
627  ,PosType smooth_factor
628  ,PosType theta2 // opening angle in radians
629  ,PosType inv_area // compensating negative mass area
630  ,PosType rmin // minimum radius for hole
631  ,PosType *xo // center of hole in Mpc
632  ,PosType *alpha // not zeroed
633  ,KappaType *kappa // not zeroed
634  ,KappaType *gamma // not zeroed
635  ,KappaType *phi // not zeroed
636  ) {
637 
638  //alpha[0]=alpha[1]=gamma[0]=gamma[1]=gamma[2]=0.0;
639  //*kappa=*phi=0.0;
640  PosType inv_area_pi = inv_area*PI;
641  PosType xcm[2],r2cm,dx[2];
642 
643  auto it = begin();
644  bool decend = true;
645  while(it.walk(decend,begin())){
646 
647  dx[0] = (*it)->center[0] - xo[0];
648  dx[1] = (*it)->center[1] - xo[1];
649  PosType d = sqrt(dx[0]*dx[0] + dx[1]*dx[1]); // distance to center of branch
650 
651  if(d < rmin - (*it)->boxsize) {
652  // branch is inside hole, skip it
653  decend = false;
654  continue;
655  }
656 
657  xcm[0] = (*it)->xcm[0] - ray[0];
658  xcm[1] = (*it)->xcm[1] - ray[1];
659  r2cm = xcm[0]*xcm[0] + xcm[1]*xcm[1];
660 
661  if(d > rmin + (*it)->boxsize){
662  if( r2cm*theta2 > ((*it)->boxsize)*((*it)->boxsize) ){
663 
664  if(r2cm*inv_area_pi < 1){
665  //????? use moments
666  double mass = (*it)->nparticles * particle_mass;
667  double prefac = mass/r2cm/PI;
668  double tmp = ( prefac - mass*inv_area);
669 
670  alpha[0] += tmp*xcm[0];
671  alpha[1] += tmp*xcm[1];
672 
673  tmp = -2.0*prefac/r2cm;
674 
675  gamma[0] += 0.5*(xcm[0]*xcm[0]-xcm[1]*xcm[1])*tmp;
676  gamma[1] += xcm[0]*xcm[1]*tmp;
677 
678  *kappa -= mass*inv_area;
679  *phi += (prefac*log(r2cm) - mass*inv_area)*r2cm*0.5;
680  }
681  decend = false;
682  }else{
683 
684  if((*it)->child == nullptr){ // leaf
685 
686  for(IndexType i=0;i<(*it)->nparticles;++i){
687 
688  PType &x = xxp[(*it)->particles[i]];
689  xcm[0] = x[0] - ray[0];
690  xcm[1] = x[1] - ray[1];
691  r2cm = xcm[0]*xcm[0] + xcm[1]*xcm[1];
692  if(r2cm*inv_area_pi < 1){
693  double prefac = particle_mass /r2cm/PI;
694  double tmp = -( prefac - particle_mass*inv_area);
695 
696  alpha[0] -= tmp*xcm[0];
697  alpha[1] -= tmp*xcm[1];
698 
699  tmp = -2.0*prefac/r2cm;
700 
701  gamma[0] += 0.5*(xcm[0]*xcm[0]-xcm[1]*xcm[1])*tmp;
702  gamma[1] += xcm[0]*xcm[1]*tmp;
703 
704  *kappa -= particle_mass*inv_area;
705  *phi += (prefac*log(r2cm)- particle_mass*inv_area)*r2cm*0.5;
706 
707  double scale = smooth_factor*(*it)->root_inv_density;
708  if(r2cm < 4*scale*scale){
709  // cubic B-spline profile
710  b_spline_profile(
711  xcm
712  ,sqrt(r2cm)
713  ,particle_mass
714  ,scale
715  ,alpha,kappa,gamma,phi
716  );
717  }
718  }
719  }
720  decend = false;
721  }else{
722  // decend
723  decend = true;
724  }
725  }
726  }else{
727  // possible particle in and outside of hole
728  if((*it)->child == nullptr){ // leaf
729 
730  for(IndexType i=0;i<(*it)->nparticles;++i){
731 
732  PType &x = xxp[(*it)->particles[i]];
733  xcm[0] = x[0] - ray[0];
734  xcm[1] = x[1] - ray[1];
735  r2cm = xcm[0]*xcm[0] + xcm[1]*xcm[1];
736 
737  dx[0] = x[0] - xo[0];
738  dx[1] = x[1] - xo[1];
739 
740  if( dx[0]*dx[0] + dx[1]*dx[1] > rmin*rmin &&
741  r2cm*inv_area_pi < 1){
742  double prefac = particle_mass /r2cm/PI;
743  double tmp = -( prefac - particle_mass*inv_area);
744 
745  alpha[0] -= tmp*xcm[0];
746  alpha[1] -= tmp*xcm[1];
747 
748  tmp = -2.0*prefac/r2cm;
749 
750  gamma[0] += 0.5*(xcm[0]*xcm[0]-xcm[1]*xcm[1])*tmp;
751  gamma[1] += xcm[0]*xcm[1]*tmp;
752 
753  *kappa -= particle_mass*inv_area;
754  *phi += (prefac*log(r2cm)- particle_mass*inv_area)*r2cm*0.5;
755 
756  double scale = smooth_factor*(*it)->root_inv_density;
757  if(r2cm < 4*scale*scale){
758  // cubic B-spline profile
759  b_spline_profile(
760  xcm
761  ,sqrt(r2cm)
762  ,particle_mass
763  ,scale
764  ,alpha,kappa,gamma,phi
765  );
766  }
767  }
768  }
769  decend = false;
770  }else{
771  // decend
772  decend = true;
773  }
774  }
775  }
776  }
777 
778  void force2D_hole(
779  const PosType *ray
780  ,PType *masses
781  ,PosType smooth_factor
782  ,PosType theta2 // opening angle in radians
783  ,PosType inv_area // compensating negative mass area
784  ,PosType rmin // minimum radius for hole
785  ,PosType *xo // center of hole in Mpc
786  ,PosType *alpha // not zeroed
787  ,KappaType *kappa // not zeroed
788  ,KappaType *gamma // not zeroed
789  ,KappaType *phi // not zeroed
790  ) {
791 
792  assert( moments_set );
793  //alpha[0]=alpha[1]=gamma[0]=gamma[1]=gamma[2]=0.0;
794  //*kappa=*phi=0.0;
795 
796  PosType inv_area_pi = inv_area*PI;
797  PosType xcm[2],r2cm,dx[2];
798 
799  auto it = begin();
800  bool decend = true;
801  while(it.walk(decend,begin())){
802 
803  dx[0] = (*it)->center[0] - xo[0];
804  dx[1] = (*it)->center[1] - xo[1];
805  PosType d = sqrt(dx[0]*dx[0] + dx[1]*dx[1]); // distance to center of branch
806 
807  if(d < rmin - (*it)->boxsize) {
808  // branch is inside hole, skip it
809  decend = false;
810  continue;
811  }
812 
813  xcm[0] = (*it)->xcm[0] - ray[0];
814  xcm[1] = (*it)->xcm[1] - ray[1];
815  r2cm = xcm[0]*xcm[0] + xcm[1]*xcm[1];
816 
817  if(d > rmin + (*it)->boxsize){
818  if( r2cm*theta2 > ((*it)->boxsize)*((*it)->boxsize) ){
819 
820  if(r2cm*inv_area_pi < 1){
821  // use moments of the branch
822  DType mass = (*it)->mass;
823  double prefac = mass/r2cm/PI;
824  double tmp = -( prefac - mass*inv_area);
825 
826  alpha[0] -= tmp*xcm[0];
827  alpha[1] -= tmp*xcm[1];
828 
829  tmp = -2.0*prefac/r2cm;
830 
831  gamma[0] += 0.5*(xcm[0]*xcm[0]-xcm[1]*xcm[1])*tmp;
832  gamma[1] += xcm[0]*xcm[1]*tmp;
833 
834  *kappa -= mass*inv_area;
835  *phi += (prefac*log(r2cm) - mass*inv_area)*r2cm*0.5;
836  }
837  decend = false;
838  }else{
839 
840  if((*it)->child == nullptr){ // leaf
841 
842  for(IndexType i=0;i<(*it)->nparticles;++i){
843  IndexType j = (*it)->particles[i];
844  PType &x = xxp[j];
845  xcm[0] = x[0] - ray[0];
846  xcm[1] = x[1] - ray[1];
847  r2cm = xcm[0]*xcm[0] + xcm[1]*xcm[1];
848  if(r2cm*inv_area_pi < 1){
849  DType mass = masses[j];
850  double prefac = mass /r2cm/PI;
851  double tmp = -( prefac - mass * inv_area);
852 
853  alpha[0] -= tmp*xcm[0];
854  alpha[1] -= tmp*xcm[1];
855 
856  tmp = -2.0*prefac/r2cm;
857 
858  gamma[0] += 0.5*(xcm[0]*xcm[0]-xcm[1]*xcm[1])*tmp;
859  gamma[1] += xcm[0]*xcm[1]*tmp;
860 
861  *kappa -= mass * inv_area;
862  *phi += (prefac*log(r2cm) - mass * inv_area)*r2cm*0.5;
863 
864  double scale = smooth_factor*(*it)->root_inv_density;
865  if(r2cm < 4*scale*scale){
866  // cubic B-spline profile
867  b_spline_profile(
868  xcm
869  ,sqrt(r2cm)
870  ,mass
871  ,scale
872  ,alpha,kappa,gamma,phi
873  );
874  }
875  }
876  }
877  decend = false;
878  }else{
879  // decend
880  decend = true;
881  }
882  }
883  }else{
884  // possible particle in and outside of hole
885  if((*it)->child == nullptr){ // leaf
886 
887  for(IndexType i=0;i<(*it)->nparticles;++i){
888  IndexType j = (*it)->particles[i];
889  PType &x = xxp[j];
890  xcm[0] = x[0] - ray[0];
891  xcm[1] = x[1] - ray[1];
892  r2cm = xcm[0]*xcm[0] + xcm[1]*xcm[1];
893 
894  dx[0] = x[0] - xo[0];
895  dx[1] = x[1] - xo[1];
896 
897  if( dx[0]*dx[0] + dx[1]*dx[1] > rmin*rmin &&
898  r2cm*inv_area_pi < 1){
899  DType mass = masses[j];
900  double prefac = mass /r2cm/PI;
901  double tmp = -( prefac - mass * inv_area);
902 
903  alpha[0] -= tmp*xcm[0];
904  alpha[1] -= tmp*xcm[1];
905 
906  tmp = -2.0*prefac/r2cm;
907 
908  gamma[0] += 0.5*(xcm[0]*xcm[0]-xcm[1]*xcm[1])*tmp;
909  gamma[1] += xcm[0]*xcm[1]*tmp;
910 
911  *kappa -= mass * inv_area;
912  *phi += (prefac*log(r2cm)- mass * inv_area)*r2cm*0.5;
913 
914  double scale = smooth_factor*(*it)->root_inv_density;
915  if(r2cm < 4*scale*scale){
916  // cubic B-spline profile
917  b_spline_profile(
918  xcm
919  ,sqrt(r2cm)
920  ,mass
921  ,scale
922  ,alpha,kappa,gamma,phi
923  );
924  }
925  }
926  }
927  decend = false;
928  }else{
929  // decend
930  decend = true;
931  }
932  }
933  }
934  }
935 
936  template<typename HType>
937  void force2D_halo(
938  const PosType *ray
939  ,HType *halos
940  //,PosType smooth_factor
941  ,PosType theta2 // opening angle in radians
942  ,PosType inv_area // compensating negative mass area
943  ,PosType rmin // minimum radius for hole
944  ,const PosType *xo // center of hole in Mpc
945  ,PosType *alpha // not zeroed
946  ,KappaType *kappa // not zeroed
947  ,KappaType *gamma // not zeroed
948  ,KappaType *phi // not zeroed
949  ) {
950 
951  assert( moments_set );
952  //alpha[0]=alpha[1]=gamma[0]=gamma[1]=gamma[2]=0.0;
953  //*kappa=*phi=0.0;
954  PosType inv_area_pi = inv_area*PI;
955  PosType xcm[2],r2cm,dx[2], tmp_v[2];
956 
957  auto it = begin();
958  bool decend = true;
959  while(it.walk(decend,begin())){
960 
961  assert( (*it)->nparticles > 0 );
962 
963  dx[0] = (*it)->center[0] - xo[0];
964  dx[1] = (*it)->center[1] - xo[1];
965  PosType d = sqrt(dx[0]*dx[0] + dx[1]*dx[1]); // distance to center of branch
966 
967  if(d < rmin - (*it)->boxsize) {
968  // branch is all inside hole, use halos
969  for(IndexType i=0;i<(*it)->nparticles;++i){
970  IndexType j = (*it)->particles[i];
971  PType &x = xxp[j];
972  xcm[0] = x[0] - ray[0];
973  xcm[1] = x[1] - ray[1];
974  r2cm = xcm[0]*xcm[0] + xcm[1]*xcm[1];
975  if(r2cm*inv_area_pi < 1){
976  // subtract off the negative mass sheet
977  DType mass = halos[j].get_mass();
978  double prefac = mass /r2cm/PI;
979  double tmp = - mass * inv_area;
980 
981  alpha[0] += tmp*xcm[0];
982  alpha[1] += tmp*xcm[1];
983 
984  *kappa -= mass * inv_area;
985  *phi += - mass * inv_area*r2cm*0.5;
986 
987  /* point mass for testing
988  prefac = mass/r2cm/PI;
989  tmp = prefac;
990 
991  alpha[0] += tmp*xcm[0];
992  alpha[1] += tmp*xcm[1];
993 
994  tmp = -2.0*prefac/r2cm;
995 
996  gamma[0] += 0.5*(xcm[0]*xcm[0]-xcm[1]*xcm[1])*tmp;
997  gamma[1] += xcm[0]*xcm[1]*tmp;
998 
999  *phi += prefac*log(r2cm) * r2cm*0.5;
1000  */
1001  // inside hole
1002  tmp_v[0] = 0;
1003  tmp_v[1] = 0;
1004  halos[j].force_halo(tmp_v,kappa,gamma,phi,xcm);
1005  alpha[0] -= tmp_v[0];
1006  alpha[1] -= tmp_v[1];
1007  }
1008  }
1009  decend = false;
1010  }else{
1011  // box outside or partly outside hole
1012  xcm[0] = (*it)->xcm[0] - ray[0];
1013  xcm[1] = (*it)->xcm[1] - ray[1];
1014  r2cm = xcm[0]*xcm[0] + xcm[1]*xcm[1];
1015 
1016  if(d > rmin + (*it)->boxsize){
1017  // box all outside hole
1018  if( r2cm*theta2 > ((*it)->boxsize)*((*it)->boxsize) ){
1019  // use moments
1020  if(r2cm*inv_area_pi < 1){
1021  // all the way out, use moments of the branch
1022  DType mass = (*it)->mass;
1023  double prefac = mass/r2cm/PI;
1024  double tmp = ( prefac - mass*inv_area);
1025 
1026  alpha[0] += tmp*xcm[0];
1027  alpha[1] += tmp*xcm[1];
1028 
1029  tmp = -2.0*prefac/r2cm;
1030 
1031  gamma[0] += 0.5*(xcm[0]*xcm[0]-xcm[1]*xcm[1])*tmp;
1032  gamma[1] += xcm[0]*xcm[1]*tmp;
1033 
1034  *kappa -= mass*inv_area;
1035  *phi += (prefac*log(r2cm) - mass*inv_area)*r2cm*0.5;
1036  }
1037  decend = false;
1038  }else{
1039  // box is too large to use moments
1040  if((*it)->child == nullptr){ // leaf
1041 
1042  for(IndexType i=0;i<(*it)->nparticles;++i){
1043  IndexType j = (*it)->particles[i];
1044  PType &x = xxp[j];
1045  xcm[0] = x[0] - ray[0];
1046  xcm[1] = x[1] - ray[1];
1047  r2cm = xcm[0]*xcm[0] + xcm[1]*xcm[1];
1048 
1049  if(r2cm*inv_area_pi < 1){
1050  DType mass = halos[j].get_mass();
1051  double prefac = mass /r2cm/PI;
1052  double tmp = ( prefac - mass * inv_area);
1053 
1054  alpha[0] += tmp*xcm[0];
1055  alpha[1] += tmp*xcm[1];
1056 
1057  tmp = -2.0*prefac/r2cm;
1058 
1059  gamma[0] += 0.5*(xcm[0]*xcm[0]-xcm[1]*xcm[1])*tmp;
1060  gamma[1] += xcm[0]*xcm[1]*tmp;
1061 
1062  *kappa -= mass * inv_area;
1063  *phi += (prefac*log(r2cm) - mass * inv_area)*r2cm*0.5;
1064 
1065  }
1066  }
1067 
1068  decend = false;
1069  }else{
1070  // decend not a leaf
1071  decend = true;
1072  }
1073  }
1074  }else{
1075  // branch is split by edge of hole
1076  if((*it)->child == nullptr){ // leaf
1077 
1078  for(IndexType i=0 ; i<(*it)->nparticles ; ++i){
1079  IndexType j = (*it)->particles[i];
1080  PType &x = xxp[j];
1081  xcm[0] = x[0] - ray[0];
1082  xcm[1] = x[1] - ray[1];
1083  r2cm = xcm[0]*xcm[0] + xcm[1]*xcm[1];
1084  if(r2cm*inv_area_pi < 1){
1085 
1086  dx[0] = x[0] - xo[0];
1087  dx[1] = x[1] - xo[1];
1088 
1089  if( dx[0]*dx[0] + dx[1]*dx[1] < rmin*rmin){
1090  DType mass = halos[j].get_mass();
1091  double prefac = mass /r2cm/PI;
1092  double tmp = - mass * inv_area;
1093 
1094  alpha[0] += tmp*xcm[0];
1095  alpha[1] += tmp*xcm[1];
1096 
1097  *kappa -= mass * inv_area;
1098  *phi += - mass * inv_area*r2cm*0.5;
1099  // inside hole
1100  tmp_v[0] = 0;
1101  tmp_v[1] = 0;
1102  halos[j].force_halo(tmp_v,kappa,gamma,phi,xcm);
1103  alpha[0] -= tmp_v[0];
1104  alpha[1] -= tmp_v[1];
1105 
1106  }else{
1107  //assert( dx[0]*dx[0] + dx[1]*dx[1] < (*begin())->boxsize*(*begin())->boxsize/4);
1108  // outside hole point mass
1109  DType mass = halos[j].get_mass();
1110  double prefac = mass/r2cm/PI;
1111  double tmp = ( prefac - mass*inv_area);
1112 
1113  alpha[0] += tmp*xcm[0];
1114  alpha[1] += tmp*xcm[1];
1115 
1116  tmp = -2.0*prefac/r2cm;
1117 
1118  gamma[0] += 0.5*(xcm[0]*xcm[0]-xcm[1]*xcm[1])*tmp;
1119  gamma[1] += xcm[0]*xcm[1]*tmp;
1120 
1121  *kappa -= mass*inv_area;
1122  *phi += (prefac*log(r2cm) - mass*inv_area)*r2cm*0.5;
1123 
1124  }
1125  }
1126  }
1127  decend = false;
1128  }else{
1129  // decend not a leaf
1130  decend = true;
1131  }
1132  }
1133  }
1134 
1135  }
1136  }
1137 
1138 
1139  size_t getUsedBranches() { return Nbranches; }
1140  size_t getTotalBranches() { return total_branches; }
1141 
1142  iterator begin()
1143  {
1144  return iterator(&top);
1145  }
1146 
1147  // maximum depth of tree
1148  int getDepth() const
1149  {
1150  return depth;
1151  }
1152 
1153  void getBoundingBox(
1154  PType &p1,
1155  PType &p2
1156  ) const {
1157  p1 = top.boundary_p1;
1158  p2 = top.boundary_p2;
1159  }
1160 private:
1161 
1162  void span8(Branch *current);
1163 
1164  Branch top;
1166  PType *xxp;
1167  std::vector<unsigned long> index;
1168 
1170  unsigned long Nbranches;
1171  unsigned long total_branches = 1;
1172  int depth = 1; // maximum depth of tree, used for debugging
1173  bool moments_set = false;
1174 
1175  PosType phiintconst = (120*log(2.) - 631.)/840 + 19./70; // ????
1176 
1177  /* cubic B-spline kernel for particle profile
1178 
1179  The lensing quantities are added to and a point mass is subtracted
1180  */
1181  inline void b_spline_profile(
1182  PosType *xcm // vector in Mpc connecting ray to center of particle
1183  ,PosType r // distance from center in Mpc
1184  ,PosType Mass // mass in solar masses
1185  ,PosType size // size scale in Mpc
1186  ,PosType *alpha // deflection angle times Sigma_crit
1187  ,KappaType *kappa // surface density
1188  ,KappaType *gamma // shear times Sigma_crit
1189  ,KappaType *phi
1190  ) const {
1191 
1192 
1193  PosType q = r/size;
1194  PosType M,sigma;
1195  if(q > 2){
1196  sigma = 0;
1197  M = 1;
1198  }else{
1199  PosType q2=q*q,q3=q2*q,q4=q2*q2,q5=q4*q;
1200 
1201  sigma = (8 - 12*q + 6*q2 - q3)/4;
1202  if(q > 1){
1203  sigma *= 10/size/size/7/PI;
1204  M = (-1 + 20*q2*(1 - q + 3*q2/8 - q3/20) )/7;
1205  *phi += Mass*(-1232. + 1200*q2 - 800.*q3 + 225.*q4 - 24*q5 + 120*log(2./q) )/840/PI;
1206  }else{
1207  sigma = 10*( sigma - 1 + 3*q - 3*q2 + q3)/size/size/7/PI;
1208  M = 10*q2*(1 - 3*q/4 + 3*q3/10)/7;
1209 
1210  *phi += Mass*( phiintconst + 10*(q2/2 - 3*q4/4 + 3*q5/50)/7
1211  )/PI;
1212  }
1213  }
1214 
1215  PosType alpha_r,gt; // deflection * Sig_crit / Mass
1216  alpha_r = (M-1)/PI/r;
1217  gt = alpha_r/r - sigma;
1218 
1219  alpha[0] += Mass*alpha_r*xcm[0]/r;
1220  alpha[1] += Mass*alpha_r*xcm[1]/r;
1221  gamma[0] -= gt*Mass*(xcm[0]*xcm[0]-xcm[1]*xcm[1])/r/r;
1222  gamma[1] -= gt*Mass*2*xcm[0]*xcm[1]/r/r;
1223  *kappa += Mass*sigma;
1224  *phi -= Mass*log(r)/PI;
1225  }
1226 
1227 };
1228 
1229 // creates 8 children to current with links to their parent and brother
1230 // and sets their center and boxsize
1231 template <typename PType,typename DType>
1232 void OTreeNB<PType,DType>::span8(Branch *current)
1233 {
1234  if (current == nullptr)
1235  {
1236  ERROR_MESSAGE();
1237  std::cerr << "OTreeNB Error: calling spon() on empty tree" << std::endl;
1238  exit(1);
1239  }
1240 
1241  current->children = new Branch[8];
1242  Branch *children = current->children;
1243 
1244  for(int i=0; i<8; ++i)
1245  {
1246  children[i].prev = current;
1247  children[i].level = current->level + 1;
1248  }
1249  total_branches += 8;
1250  if(depth < current->level + 1) depth = current->level + 1;
1251 
1252  PType &center = current->center;
1253  //center[0] = (current->boundary_p1[0] + current->boundary_p2[0]) / 2.0;
1254  //center[1] = (current->boundary_p1[1] + current->boundary_p2[1]) / 2.0;
1255  //center[2] = (current->boundary_p1[2] + current->boundary_p2[2]) / 2.0;
1256  PosType boxsize = current->boxsize / 2.0;
1257 
1258  size_t m = 0;
1259  for (int i = -1; i < 2; i += 2)
1260  {
1261  for (int j = -1; j < 2; j += 2)
1262  {
1263  for (int k = -1; k < 2; k += 2)
1264  {
1265 
1266  children[m].boundary_p1[0] = ((1 - i) / 2) * current->boundary_p1[0]
1267  + ((1 + i) / 2) * center[0];
1268  children[m].boundary_p1[1] = ((1 - j) / 2) * current->boundary_p1[1]
1269  + ((1 + j) / 2) * center[1];
1270  children[m].boundary_p1[2] = ((1 - k) / 2) * current->boundary_p1[2]
1271  + ((1 + k) / 2) * center[2];
1272 
1273  children[m].boundary_p2[0] = ((1 + i) / 2) * current->boundary_p2[0]
1274  + ((1 - i) / 2) * center[0];
1275  children[m].boundary_p2[1] = ((1 + j) / 2) * current->boundary_p2[1]
1276  + ((1 - j) / 2) * center[1];
1277  children[m].boundary_p2[2] = ((1 + k) / 2) * current->boundary_p2[2]
1278  + ((1 - k) / 2) * center[2];
1279 
1280  children[m].center = (children[m].boundary_p1 + children[m].boundary_p2) * 0.5;
1281 
1282  children[m].boxsize = boxsize;
1283  ++m;
1284  }
1285  }
1286  }
1287 
1288  // sort particles into children
1289  IndexType *p = current->particles;
1290  //std::cout << p[0] << " " << p[1] << " " << p[2] << std::endl;
1291  IndexType *p_end = p + current->nparticles;
1292  IndexType np = 0;
1293  for (int m = 0; m < 8; ++m)
1294  {
1295  children[m].particles = p;
1296  IndexType *pp = p;
1297  //std::cout << pp[0] << " " << pp[1] << " " << pp[2] << std::endl;
1298  //std::cout << xxp[pp[0]][1] << " " << xxp[*pp][1] << " " << pp[2] << std::endl;
1299 
1300  while (pp != p_end)
1301  {
1302  PType &x = xxp[*pp];
1303  if (x[0] >= children[m].boundary_p1[0] &&
1304  x[0] <= children[m].boundary_p2[0] &&
1305  x[1] >= children[m].boundary_p1[1] &&
1306  x[1] <= children[m].boundary_p2[1] &&
1307  x[2] >= children[m].boundary_p1[2] &&
1308  x[2] <= children[m].boundary_p2[2])
1309  {
1310  std::swap(*p, *pp);
1311  children[m].nparticles++;
1312  ++p;
1313  }
1314  ++pp;
1315  }
1316  np += children[m].nparticles;
1317 
1318  children[m].root_inv_density = children[m].boxsize/pow(children[m].nparticles,1.0/3.); // reset density
1319  }
1320  assert( current->nparticles == np );
1321 
1322  // remove empty branches from the brotherhood
1323  int i=0;
1324  while(children[i].nparticles==0 && i<8) ++i;
1325  assert(i != 8);
1326  current->child = children+i;
1327  ++Nbranches;
1328  int j=i+1;
1329  while(j<8){
1330  if(children[j].nparticles > 0){
1331  children[i].brother = children + j;
1332  i=j;
1333  ++Nbranches;
1334  }
1335  ++j;
1336  }
1337  children[i].brother = current->brother;
1338 }
1339 
1340 #endif /* OTreeNB_h */
OTreeNB::Branch::particles
IndexType * particles
array of particles in Branch
Definition: oTreeNB.h:76
OTreeNB::iterator::walk
bool walk(bool allowDescent, iterator end)
Definition: oTreeNB.h:186
Branch
The box representing a branch of a binary tree structure. Used specifically in TreeStruct for organiz...
Definition: point.h:677
OTreeNB::Branch::boundary_p1
PType boundary_p1
bottom, left, back corner of box
Definition: oTreeNB.h:83
OTreeNB::iterator::down
bool down(short child)
Move to child.
Definition: oTreeNB.h:156
OTreeNB::getBoundingBox
void getBoundingBox(PType &p1, PType &p2) const
Definition: oTreeNB.h:1153
OTreeNB::Branch::boundary_p2
PType boundary_p2
top, right, front corner of box
Definition: oTreeNB.h:85
OTreeNB::Branch::level
int level
level in tree
Definition: oTreeNB.h:92
OTreeNB::iterator::operator*
Branch * operator*()
Returns a pointer to the current Branch.
Definition: oTreeNB.h:142
OTreeNB::Branch::prev
Branch * prev
father of branch
Definition: oTreeNB.h:100
OTreeNB::Branch::brother
Branch * brother
Definition: oTreeNB.h:103
OTreeNB::calcMoments
void calcMoments()
Definition: oTreeNB.h:305
OTreeNB
OTreeNB: Tree structure used for force calculation with particles (i.e. stars, Nbody or halos).
Definition: oTreeNB.h:42
OTreeNB::calcMoments
void calcMoments(DType *masses)
Definition: oTreeNB.h:347
OTreeNB::iterator::iterator
iterator(OTreeNB< PType > *tree)
Sets the top or root to the top of "tree".
Definition: oTreeNB.h:125
OTreeNB::Branch
classes
Definition: oTreeNB.h:48
OTreeNB::iterator::iterator
iterator(Branch *branch)
Sets the root to the input branch so that this will be a subtree in branch is not the real root.
Definition: oTreeNB.h:127
OTreeNB::calcMoments
void calcMoments(std::vector< HType > &halos)
Definition: oTreeNB.h:393
OTreeNB::iterator
A iterator class that allows for movement through the tree without changing anything in the tree itse...
Definition: oTreeNB.h:118
Generated on Sun Sep 27 2026 14:12:56 for GLAMERDOC++ by   doxygen 1.8.18