GLAMERDOC++
Gravitational Lensing Code Library
Loading...
Searching...
No Matches
quadTreeHalos.h
1//
2// quadTreeHalos.h
3// GLAMER
4//
5// Created by Robert Benton Metcalf on 30/10/2020.
6//
7
8#ifndef quadTreeHalos_h
9#define quadTreeHalos_h
10
11#include "standard.h"
12#include "utilities_slsim.h"
13#include "Tree.h"
14#include "qTreeNB.h"
15
34template <typename LensHaloType>
35class TreeQuadHalos {
36public:
37 TreeQuadHalos(
38 LensHaloType **my_halos
39 ,IndexType Npoints
40 ,PosType my_inv_area = 0
41 ,int bucket = 5
42 ,PosType theta_force = 0.5
43 ,bool my_periodic_buffer = false
44 ,PosType my_inv_screening_scale = 0
45 ,PosType maximum_range = -1
46 );
47 ~TreeQuadHalos();
48
49 friend class LensPlaneTree;
50 void force2D(PosType const *ray,PosType *alpha,KappaType *kappa,KappaType *gamma
51 ,KappaType *phi) const;
52
53 //void force2D_recur(const PosType *ray,PosType *alpha,KappaType *kappa
54 // ,KappaType *gamma,KappaType *phi);
55
57 void neighbors(PosType ray[],PosType rmax,std::list<IndexType> &neighbors) const;
58 void neighbors(PosType ray[],PosType rmax,std::vector<LensHaloType *> &neighbors) const;
59
60 void printParticlesInBranch(unsigned long number);
61
62 void printBranchs(int level = -1);
63
64protected:
65
66 std::vector<PosType> workspace;
67 PosType **xp;
68 bool MultiMass;
69 bool MultiRadius;
70
71 IndexType Nparticles;
72 PosType inv_area;
73 int Nbucket;
74
75 PosType force_theta;
76 PosType max_range;
77
78 QTreeNB<PosType *> * tree;
79 std::vector<IndexType> index;
80
81 //bool haloON;
82 LensHaloType **halos;
83
84 PosType realray[2];
85 int incell,incell2;
86
88 bool periodic_buffer;
89 PosType inv_screening_scale2;
90 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.
91 PosType original_yl; // x-axis size of simulation used for peridic buffering.
92
93 QTreeNB<PosType *> * BuildQTreeNB(PosType **xp,IndexType Nparticles,IndexType *particles);
94 void _BuildQTreeNB(IndexType nparticles,IndexType *particles);
95
96 inline short WhichQuad(PosType *x,QBranchNB &branch);
97
98 //inline bool atLeaf();
99 inline bool inbox(const PosType *ray,const PosType *p1,const PosType *p2){
100 return (ray[0]>=p1[0])*(ray[0]<=p2[0])*(ray[1]>=p1[1])*(ray[1]<=p2[1]);
101 }
102 //int cutbox(PosType *ray,PosType *p1,PosType *p2,float rmax);
103
104 void CalcMoments();
105 void rotate_coordinates(PosType **coord);
106
107 //QTreeNBHndl rotate_simulation(PosType **xp,IndexType Nparticles,IndexType *particles
108 // ,PosType **coord,PosType theta,float *rsph,float *mass
109 // ,bool MultiRadius,bool MultiMass);
110 //QTreeNBHndl rotate_project(PosType **xp,IndexType Nparticles,IndexType *particles
111 // ,PosType **coord,PosType theta,float *rsph,float *mass
112 // ,bool MultiRadius,bool MultiMass);
113 void cuttoffscale(QTreeNB<PosType *> * tree,PosType *theta);
114
115 //void walkTree_recur(QBranchNB *branch,PosType const *ray,PosType *alpha,KappaType *kappa,KappaType *gamma,KappaType *phi);
116 void walkTree_iter(QBiterator<PosType *> &treeit, PosType const *ray,PosType *alpha,KappaType *kappa
117 ,KappaType *gamma,KappaType *phi) const;
118
119 PosType phiintconst;
120
121};
122
123
128template <typename LensHaloType>
129TreeQuadHalos<LensHaloType>::TreeQuadHalos(
130 LensHaloType **my_halos
131 ,IndexType Npoints
132 ,PosType my_inv_area
133 ,int bucket
134 ,PosType my_theta_force
135 ,bool periodic_buffer
136 ,PosType my_inv_screening_scale
137 ,PosType maximum_range
138):
139MultiMass(true),MultiRadius(true)//,masses(NULL),sizes(NULL)
140,Nparticles(Npoints),inv_area(my_inv_area),Nbucket(bucket)
141,force_theta(my_theta_force)
142,max_range(maximum_range),halos(my_halos),periodic_buffer(periodic_buffer)
143,inv_screening_scale2(my_inv_screening_scale*my_inv_screening_scale)
144{
145 index.resize(Npoints);
146 IndexType ii;
147
148 for(ii=0;ii<Npoints;++ii) index[ii] = ii;
149
150 // copy halo positions into xp array to make compatable with QTreeNB
151 xp = Utilities::PosTypeMatrix(Npoints,2);
152 for(ii=0;ii<Npoints;++ii) halos[ii]->getX(xp[ii]);
153
154 if(Npoints > 0){
155 tree = BuildQTreeNB(xp,Npoints,index.data());
156 CalcMoments();
157 }
158
159 phiintconst = (120*log(2.) - 631.)/840 + 19./70;
160 return;
161}
162
164template <typename LensHaloType>
165TreeQuadHalos<LensHaloType>::~TreeQuadHalos()
166{
167 if(Nparticles == 0) return;
168 delete tree;
169 Utilities::free_PosTypeMatrix(xp,Nparticles,2);
170 return;
171}
172
173template <typename LensHaloType>
174QTreeNB<PosType *> * TreeQuadHalos<LensHaloType>::BuildQTreeNB(PosType **xp,IndexType Nparticles,IndexType *particles){
175 IndexType i;
176 short j;
177 PosType p1[2],p2[2];
178
179 for(j=0;j<2;++j){
180 p1[j]=xp[0][j];
181 p2[j]=xp[0][j];
182 }
183
184 // Find borders that enclose all particles
185 for(i=0;i<Nparticles;++i){
186 for(j=0;j<2;++j){
187 if(xp[i][j] < p1[j] ) p1[j]=xp[i][j];
188 if(xp[i][j] > p2[j] ) p2[j]=xp[i][j];
189 }
190 }
191
192 if(Nparticles == 1){
193 p1[0]=-0.25;
194 p1[1]=-0.25;
195 p2[0]=0.25;
196 p2[1]=0.25;
197 }
198
199 // store true dimensions of simulation
200 PosType lengths[2] = {p2[0]-p1[0],p2[1]-p1[1]};
201 original_xl = lengths[0];
202 original_yl = lengths[1];
203
204 if(lengths[0] == 0 || lengths[1] == 0){
205 throw std::invalid_argument("particles in same place.");
206 }
207
208 // If region is not square, make it square.
209 j = lengths[0] > lengths[1] ? 1 : 0;
210 p2[j] = p1[j] + lengths[!j];
211
212 /* Initialize tree root */
213 tree = new QTreeNB<PosType *>(xp,particles,Nparticles,p1,p2);
214
215 /* build the tree */
216 workspace.resize(Nparticles);
217 _BuildQTreeNB(Nparticles,particles);
218 workspace.clear();
219 workspace.shrink_to_fit();
220
221 /* visit every branch to find center of mass and cutoff scale */
222 tree->moveTop();
223
224 return tree;
225}
226
228template <typename LensHaloType>
229inline short TreeQuadHalos<LensHaloType>::WhichQuad(PosType *x,QBranchNB &branch){
230 return (x[0] < branch.center[0]) + 2*(x[1] < branch.center[1]);
231}
232
234template <typename LensHaloType>
235void TreeQuadHalos<LensHaloType>::_BuildQTreeNB(IndexType nparticles,IndexType *particles){
236
237 QBranchNB *cbranch = tree->current; /* pointer to current branch */
238 IndexType i,j,cut,cut2;
239
240 cbranch->center[0] = (cbranch->boundary_p1[0] + cbranch->boundary_p2[0])/2;
241 cbranch->center[1] = (cbranch->boundary_p1[1] + cbranch->boundary_p2[1])/2;
242 cbranch->quad[0]=cbranch->quad[1]=cbranch->quad[2]=0;
243 cbranch->mass = 0.0;
244
245 double theta_range = 2*force_theta;
246 if(max_range > 0){
247 double boxsize = 1.732*(cbranch->boundary_p2[0] - cbranch->boundary_p1[0]);
248 theta_range = boxsize / MAX(sqrt(cbranch->center[0]*cbranch->center[0]
249 + cbranch->center[1]*cbranch->center[1] )
250 - max_range, boxsize );
251 }
252
253 // leaf case
254 if(cbranch->nparticles <= Nbucket
255 || force_theta > theta_range
256 ){
257
258 cbranch->Nbig_particles = cbranch->nparticles;
259 cbranch->big_particles.reset( new IndexType[cbranch->Nbig_particles] );
260 for(i=0,j=0;i<cbranch->nparticles;++i){
261 cbranch->big_particles[i] = particles[i];
262 }
263
264// PosType r;
265// cbranch->Nbig_particles = 0;
266// for(i=0;i<cbranch->nparticles;++i){
267// jt = particles[i]*MultiRadius;
268// r = halos[jt]->get_Rmax();
269// //r = haloON ? halos[particles[i]*MultiRadius].Rmax : sizes[particles[i]*MultiRadius];
270// if(r < (cbranch->boundary_p2[0]-cbranch->boundary_p1[0])) ++(cbranch->Nbig_particles);
271// }
272// if(cbranch->Nbig_particles){
273// //cbranch->big_particles = new IndexType[cbranch->Nbig_particles];
274// cbranch->big_particles.reset( new IndexType[cbranch->Nbig_particles] );
275// for(i=0,j=0;i<cbranch->nparticles;++i){
276// jt = particles[i]*MultiRadius;
277// r = halos[jt]->get_Rmax();
278// if(r < (cbranch->boundary_p2[0]-cbranch->boundary_p1[0])) cbranch->big_particles[j++] = particles[i];
279// }
280// }
281// else{
282// cbranch->big_particles.reset(nullptr);
283// }
284//
285 return;
286 }
287
288 // find particles too big to be in children
289
290 //PosType *x = new PosType[cbranch->nparticles];
291 PosType *x = workspace.data();
292
293 cbranch->Nbig_particles=0;
294
295 if(MultiRadius){
296 // store the particles that are too large to be in a child at the end of the list of particles in cbranch
297 for(i=0;i<cbranch->nparticles;++i){
298 x[i] = halos[particles[i]]->get_Rmax();
299 }
300 PosType half_box = (cbranch->boundary_p2[0]-cbranch->boundary_p1[0])/2;
301
302 // sort particles in size
303 Utilities::quickPartition(half_box,&cut,particles,x,cbranch->nparticles);
304
305 if(cut < cbranch->nparticles){
306 if(tree->atTop()){
307 cbranch->Nbig_particles = cut2 = cbranch->nparticles-cut;
308 }else{
309 Utilities::quickPartition(2*half_box,&cut2,&particles[cut],&x[cut],cbranch->nparticles-cut);
310 cbranch->Nbig_particles = cut2;
311 }
312
313 //cbranch->big_particles = new IndexType[cbranch->Nbig_particles];
314 cbranch->big_particles.reset( new IndexType[cbranch->Nbig_particles] );
315 for(i=cut;i<(cut+cut2);++i) cbranch->big_particles[i-cut] = particles[i];
316 }
317
318 }else{
319 x[0] = halos[0]->get_Rmax();
320 if(x[0] > (cbranch->boundary_p2[0]-cbranch->boundary_p1[0])/2
321 && x[0] < (cbranch->boundary_p2[0]-cbranch->boundary_p1[0])){
322 cbranch->Nbig_particles = cbranch->nparticles;
323 cbranch->big_particles.reset( new IndexType[cbranch->Nbig_particles] );
324 for(i=0;i<cbranch->Nbig_particles;++i) cbranch->big_particles[i] = particles[i];
325
326 }
327 }
328
329 assert( cbranch->nparticles >= cbranch->Nbig_particles );
330 IndexType NpInChildren = cbranch->nparticles;// - cbranch->Nbig_particles;
331 assert(NpInChildren >= 0);
332
333 if(NpInChildren == 0){
334 //delete[] x;
335 return;
336 }
337
338 IndexType cutx,cuty;
339 PosType xcut,ycut;
340
341 QBranchNB *child0 = new QBranchNB(cbranch);
342 QBranchNB *child1 = new QBranchNB(cbranch);
343 QBranchNB *child2 = new QBranchNB(cbranch);
344 QBranchNB *child3 = new QBranchNB(cbranch);
345
346 tree->attachChildrenToCurrent(child0,child1,child2,child3);
347
348 // initialize boundaries of children
349
350 child0->boundary_p1[0]=cbranch->center[0];
351 child0->boundary_p1[1]=cbranch->center[1];
352 child0->boundary_p2[0]=cbranch->boundary_p2[0];
353 child0->boundary_p2[1]=cbranch->boundary_p2[1];
354 child0->boxsize2 = (child0->boundary_p2[0] - child0->boundary_p1[0])*(child0->boundary_p2[0] - child0->boundary_p1[0]);
355
356 child1->boundary_p1[0]=cbranch->boundary_p1[0];
357 child1->boundary_p1[1]=cbranch->center[1];
358 child1->boundary_p2[0]=cbranch->center[0];
359 child1->boundary_p2[1]=cbranch->boundary_p2[1];
360 child1->boxsize2 = (child1->boundary_p2[0] - child1->boundary_p1[0])*(child1->boundary_p2[0] - child1->boundary_p1[0]);
361
362 child2->boundary_p1[0]=cbranch->center[0];
363 child2->boundary_p1[1]=cbranch->boundary_p1[1];
364 child2->boundary_p2[0]=cbranch->boundary_p2[0];
365 child2->boundary_p2[1]=cbranch->center[1];
366 child2->boxsize2 = (child2->boundary_p2[0] - child2->boundary_p1[0])*(child2->boundary_p2[0] - child2->boundary_p1[0]);
367
368 child3->boundary_p1[0]=cbranch->boundary_p1[0];
369 child3->boundary_p1[1]=cbranch->boundary_p1[1];
370 child3->boundary_p2[0]=cbranch->center[0];
371 child3->boundary_p2[1]=cbranch->center[1];
372 child3->boxsize2 = (child3->boundary_p2[0] - child3->boundary_p1[0])*(child3->boundary_p2[0] - child3->boundary_p1[0]);
373
374
375 // **** makes sure force does not require nbucket at leaf
376 xcut = cbranch->center[0];
377 ycut = cbranch->center[1];
378
379 // divide along y direction
380 for(i=0;i<NpInChildren;++i) x[i] = tree->xxp[particles[i]][1];
381 Utilities::quickPartition(ycut,&cuty,particles,x,NpInChildren);
382
383 if(cuty > 0){
384 // divide first group in the x direction
385 for(i=0;i<cuty;++i) x[i] = tree->xxp[particles[i]][0];
386 Utilities::quickPartition(xcut,&cutx,particles,x,cuty);
387
388 child3->nparticles=cutx;
389 assert(child3->nparticles >= 0);
390 if(child3->nparticles > 0)
391 child3->particles = particles;
392 else child3->particles = NULL;
393
394 child2->nparticles = cuty - cutx;
395 assert(child2->nparticles >= 0);
396 if(child2->nparticles > 0)
397 child2->particles = &particles[cutx];
398 else child2->particles = NULL;
399
400 }else{
401 child3->nparticles = 0;
402 child3->particles = NULL;
403
404 child2->nparticles = 0;
405 child2->particles = NULL;
406 }
407
408 if(cuty < NpInChildren){
409 // divide second group in the x direction
410 for(i=cuty;i<NpInChildren;++i) x[i-cuty] = tree->xxp[particles[i]][0];
411 Utilities::quickPartition(xcut,&cutx,&particles[cuty],x,NpInChildren-cuty);
412
413 child1->nparticles=cutx;
414 assert(child1->nparticles >= 0);
415 if(child1->nparticles > 0)
416 child1->particles = &particles[cuty];
417 else child1->particles = NULL;
418
419 child0->nparticles=NpInChildren - cuty - cutx;
420 assert(child0->nparticles >= 0);
421 if(child0->nparticles > 0)
422 child0->particles = &particles[cuty + cutx];
423 else child0->particles = NULL;
424
425 }else{
426 child1->nparticles = 0;
427 child1->particles = NULL;
428
429 child0->nparticles = 0;
430 child0->particles = NULL;
431 }
432
433 //delete[] x;
434
435 tree->moveToChild(0);
436 _BuildQTreeNB(child0->nparticles,child0->particles);
437 tree->moveUp();
438
439 tree->moveToChild(1);
440 _BuildQTreeNB(child1->nparticles,child1->particles);
441 tree->moveUp();
442
443 tree->moveToChild(2);
444 _BuildQTreeNB(child2->nparticles,child2->particles);
445 tree->moveUp();
446
447 tree->moveToChild(3);
448 _BuildQTreeNB(child3->nparticles,child3->particles);
449 tree->moveUp();
450
451 return;
452}
453
454// calculates moments of the mass and the cutoff scale for each box
455template <typename LensHaloType>
456void TreeQuadHalos<LensHaloType>::CalcMoments(){
457
458 //*** make compatable
459 IndexType i;
460 PosType rcom,xcm[2],xcut;
461 QBranchNB *cbranch;
462 PosType tmp;
463
464 tree->moveTop();
465 do{
466 cbranch=tree->current; /* pointer to current branch */
467
468 cbranch->rmax = sqrt( pow(cbranch->boundary_p2[0]-cbranch->boundary_p1[0],2)
469 + pow(cbranch->boundary_p2[1]-cbranch->boundary_p1[1],2) );
470
471 // calculate mass
472 for(i=0,cbranch->mass=0;i<cbranch->nparticles;++i)
473 cbranch->mass += halos[cbranch->particles[i]]->get_mass();
474 //cbranch->mass += haloON ? halos[cbranch->particles[i]*MultiMass].mass : masses[cbranch->particles[i]*MultiMass];
475
476 // calculate center of mass
477 cbranch->center[0]=cbranch->center[1]=0;
478 for(i=0;i<cbranch->nparticles;++i){
479 tmp = halos[cbranch->particles[i]]->get_mass();
480 //tmp = haloON ? halos[cbranch->particles[i]*MultiMass].mass : masses[cbranch->particles[i]*MultiMass];
481 cbranch->center[0] += tmp*tree->xxp[cbranch->particles[i]][0]/cbranch->mass;
482 cbranch->center[1] += tmp*tree->xxp[cbranch->particles[i]][1]/cbranch->mass;
483 }
485 // calculate quadropole moment of branch assuming point masses
487 cbranch->quad[0]=cbranch->quad[1]=cbranch->quad[2]=0;
488 for(i=0;i<cbranch->nparticles;++i){
489 xcm[0]=tree->xxp[cbranch->particles[i]][0]-cbranch->center[0];
490 xcm[1]=tree->xxp[cbranch->particles[i]][1]-cbranch->center[1];
491 xcut = pow(xcm[0],2) + pow(xcm[1],2);
492
493 tmp = halos[cbranch->particles[i]]->get_mass();
494
495 cbranch->quad[0] += (xcut-2*xcm[0]*xcm[0])*tmp;
496 cbranch->quad[1] += (xcut-2*xcm[1]*xcm[1])*tmp;
497 cbranch->quad[2] += -2*xcm[0]*xcm[1]*tmp;
498 }
499
500 // largest distance from center of mass of cell
501 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) ) ?
502 pow(cbranch->center[i]-cbranch->boundary_p1[i],2) : pow(cbranch->center[i]-cbranch->boundary_p2[i],2);
503
504 rcom=sqrt(rcom);
505
506 if(force_theta > 0.0){
507 cbranch->r2crit_angle = 1.15470*rcom/(force_theta);
508 cbranch->r2crit_angle *= cbranch->r2crit_angle;
509 }else cbranch->r2crit_angle=1.0e100;
510
511 }while(tree->WalkStep(true));
512
513 return;
514}
515
517template <typename LensHaloType>
518void TreeQuadHalos<LensHaloType>::rotate_coordinates(PosType **coord){
519 IndexType i;
520 short j;
521 PosType tmp[3];
522
523 /* rotate particle positions */
524 for(i=0;i<tree->top->nparticles;++i){
525 for(j=0;j<3;++j) tmp[j]=0.0;
526 for(j=0;j<3;++j){
527 tmp[0]+=coord[0][j]*xp[i][j];
528 tmp[1]+=coord[1][j]*xp[i][j];
529 tmp[2]+=coord[2][j]*xp[i][j];
530 }
531 for(j=0;j<3;++j) xp[i][j]=tmp[j];
532 }
533
534 return;
535}
536
556template <typename LensHaloType>
557void TreeQuadHalos<LensHaloType>::force2D(const PosType *ray,PosType *alpha,KappaType *kappa,KappaType *gamma,KappaType *phi) const{
558
559 alpha[0]=alpha[1]=gamma[0]=gamma[1]=gamma[2]=0.0;
560 *kappa=*phi=0.0;
561
562 if(Nparticles == 0) return;
563
564 assert(tree);
565 //QTreeNB<PosType *>::iterator it(tree);
566 QBiterator<PosType *> it(tree);
567
568
569 if(periodic_buffer){
570 PosType tmp_ray[2];
571
572 for(int i = -1;i<2;++i){
573 for(int j = -1;j<2;++j){
574 tmp_ray[0] = ray[0] + i*original_xl;
575 tmp_ray[1] = ray[1] + j*original_yl;
576 walkTree_iter(it,tmp_ray,alpha,kappa,gamma,phi);
577 }
578 }
579 }else{
580 walkTree_iter(it,ray,alpha,kappa,gamma,phi);
581 }
582
583 // Subtract off uniform mass sheet to compensate for the extra mass
584 // added to the universe in the halos.
585 //alpha[0] -= ray[0]*sigma_background*(inv_screening_scale2 == 0);
586 //alpha[1] -= ray[1]*sigma_background*(inv_screening_scale2 == 0);
587
588 //*kappa -= sigma_background;
589
590 return;
591}
592
595template <typename LensHaloType>
596void TreeQuadHalos<LensHaloType>::neighbors(PosType ray[],PosType rmax,std::list<IndexType> &neighbors) const{
597 //QTreeNB::iterator it(tree);
598 QBiterator<PosType *> it(tree);
599
600 neighbors.clear();
601
602 PosType r2 = rmax*rmax;
603 bool decend=true;
604 it.movetop();
605 do{
606 int cut = Utilities::cutbox(ray,(*it)->boundary_p1,(*it)->boundary_p2,rmax);
607 decend = true;
608
609 if(cut == 0){ // box outside circle
610 decend = false;
611 }else if(cut == 1){ // whole box inside circle
612 for(int i=0;i<(*it)->nparticles;++i) neighbors.push_back((*it)->particles[i]);
613 decend = false;
614 }else if((*it)->child0 == NULL){ // at leaf
615 for(int i=0;i<(*it)->nparticles;++i){
616 if( (tree->xxp[(*it)->particles[i]][0] - ray[0])*(tree->xxp[(*it)->particles[i]][0] - ray[0])
617 + (tree->xxp[(*it)->particles[i]][1] - ray[1])*(tree->xxp[(*it)->particles[i]][1] - ray[1]) < r2) neighbors.push_back((*it)->particles[i]);
618 }
619 decend = false;
620 }
621 }while(it.TreeWalkStep(decend));
622
623}
626template <typename LensHaloType>
627void TreeQuadHalos<LensHaloType>::neighbors(PosType ray[],PosType rmax,std::vector<LensHaloType *> &neighbors) const{
628 neighbors.clear();
629
630 if(halos == NULL){
631 std::cerr << "TreeQuadHalos<LensHaloType>::neighbors - The are no halos in this tree use other version of this function" << std::endl;
632 throw std::runtime_error("no halos");
633 return;
634 }
635
636 std::list<IndexType> neighbor_indexes;
637 TreeQuadHalos<LensHaloType>::neighbors(ray,rmax,neighbor_indexes);
638
639 neighbors.resize(neighbor_indexes.size());
640 std::list<IndexType>::iterator it = neighbor_indexes.begin();
641 for(size_t i=0;i<neighbors.size();++i,++it){
642 neighbors[i] = halos[*it];
643 }
644
645 return;
646}
647
648template <typename LensHaloType>
649void TreeQuadHalos<LensHaloType>::walkTree_iter(
650 QBiterator<PosType *> &treeit,
651 const PosType *ray
652 ,PosType *alpha
653 ,KappaType *kappa
654 ,KappaType *gamma
655 ,KappaType *phi
656 ) const
657{
658
659 bool allowDescent;
660 unsigned long count=0;//,tmp_index;
661
662 treeit.movetop();
663 do{
664 ++count;
665 allowDescent=true;
666
667 // add big paritcles, these are all the particles in tha case of a leaf
668
669 for(size_t i = 0 ; i < (*treeit)->Nbig_particles ; ++i){
670
671 size_t tmp_index = (*treeit)->big_particles[i];
672
673 PosType xcm[2];
674 xcm[0] = tree->xxp[tmp_index][0] - ray[0];
675 xcm[1] = tree->xxp[tmp_index][1] - ray[1];
676
677 PosType rcm2 = xcm[0]*xcm[0] + xcm[1]*xcm[1];
678
679 double screening = exp(-rcm2*inv_screening_scale2);
680 halos[tmp_index]->force_halo(alpha,kappa,gamma,phi,xcm,false,screening);
681
682 // background part
683 double mass = halos[tmp_index]->get_mass();
684 double prefac = mass * inv_area ;
685
686 alpha[0] += prefac*xcm[0];
687 alpha[1] += prefac*xcm[1];
688
689 *phi += - inv_area*rcm2 * mass*0.5;
690 *kappa -= mass * inv_area;
691
692 }
693
694 if( !(treeit.atLeaf()) ){
695 PosType xcm_cell[2];
696
697 xcm_cell[0]=(*treeit)->center[0]-ray[0];
698 xcm_cell[1]=(*treeit)->center[1]-ray[1];
699
700 PosType rcm2cell = xcm_cell[0]*xcm_cell[0] + xcm_cell[1]*xcm_cell[1];
701
702 if(rcm2cell > (*treeit)->r2crit_angle ){ // use whole cell
703
704 allowDescent=false;
705
706 double screening = exp(-rcm2cell*inv_screening_scale2);
707 double tmp = (*treeit)->mass*(inv_area - screening/rcm2cell/PI);
708
709 alpha[0] += tmp*xcm_cell[0];
710 alpha[1] += tmp*xcm_cell[1];
711
712 { // taken out to speed up
713 tmp=-2.0*(*treeit)->mass/PI/rcm2cell/rcm2cell*screening;
714 gamma[0] += 0.5*(xcm_cell[0]*xcm_cell[0]-xcm_cell[1]*xcm_cell[1])*tmp;
715 gamma[1] += xcm_cell[0]*xcm_cell[1]*tmp;
716
717 *phi += 0.5*(*treeit)->mass*log( rcm2cell )/PI*screening;
718 *phi -= 0.5*( (*treeit)->quad[0]*xcm_cell[0]*xcm_cell[0] + (*treeit)->quad[1]*xcm_cell[1]*xcm_cell[1] + 2*(*treeit)->quad[2]*xcm_cell[0]*xcm_cell[1] )/(PI*rcm2cell*rcm2cell)*screening;
719 }
720
721 // quadrapole contribution
722 // the kappa and gamma are not calculated to this order
723 alpha[0] -= ((*treeit)->quad[0]*xcm_cell[0] + (*treeit)->quad[2]*xcm_cell[1])
724 /(rcm2cell*rcm2cell)/PI*screening;
725 alpha[1] -= ((*treeit)->quad[1]*xcm_cell[1] + (*treeit)->quad[2]*xcm_cell[0])
726 /(rcm2cell*rcm2cell)/PI*screening;
727
728 tmp = 4*((*treeit)->quad[0]*xcm_cell[0]*xcm_cell[0] + (*treeit)->quad[1]*xcm_cell[1]*xcm_cell[1]
729 + 2*(*treeit)->quad[2]*xcm_cell[0]*xcm_cell[1])/(rcm2cell*rcm2cell*rcm2cell)/PI*screening;
730
731 alpha[0] += tmp*xcm_cell[0];
732 alpha[1] += tmp*xcm_cell[1];
733
734 *kappa -= (*treeit)->mass*inv_area;
735 *phi -= (*treeit)->mass*inv_area*rcm2cell;
736
737 }
738 }
739 assert(count <= treeit.getNbranches());
740 }while(treeit.TreeWalkStep(allowDescent));
741
742 // if((*treeit)->nparticles > 0)
743 // {
744 //
745 //
746 //
747 // // Find the particles that intersect with ray and add them individually.
748 // //if(rcm2cell < 5.83*((*treeit)->boxsize2))
749 //
750 // //boxsize2 = ((*treeit)->boundary_p2[0]-(*treeit)->boundary_p1[0])*((*treeit)->boundary_p2[0]-(*treeit)->boundary_p1[0]);
751 //
752 // if( rcm2cell < ((*treeit)->rcrit_angle)*((*treeit)->rcrit_angle) )
754 // {
755 //
756 // // includes rcrit_particle constraint
757 // allowDescent=true;
758 //
759 //
760 // // Treat all particles in a leaf as a point particle
761 // if(treeit.atLeaf()){
762 //
763 // for(i = 0 ; i < (*treeit)->nparticles ; ++i)
764 // {
765 //
766 // xcm_cell[0] = tree->xxp[(*treeit)->particles[i]][0] - ray[0];
767 // xcm_cell[1] = tree->xxp[(*treeit)->particles[i]][1] - ray[1];
768 //
769 //
770 // rcm2 = xcm_cell[0]*xcm_cell[0] + xcm_cell[1]*xcm_cell[1];
771 // double screening = exp(-rcm2*inv_screening_scale2);
772 // if(rcm2 < 1e-20) rcm2 = 1e-20;
773 //
774 // size_t tmp_index = MultiMass*(*treeit)->particles[i];
775 //
776 // double mass = halos[tmp_index]->get_mass();
777 // prefac = mass * (inv_area - screening/rcm2/PI) ;
778 // //prefac /= rcm2*PI/screening;
779 // //tmp = -1.0*prefac;
780 //
781 // alpha[0] += prefac*xcm_cell[0];
782 // alpha[1] += prefac*xcm_cell[1];
783 //
784 // tmp = -2.0*mass*screening/PI/rcm2/rcm2;
785 //
786 // gamma[0] += 0.5*(xcm_cell[0]*xcm_cell[0]-xcm_cell[1]*xcm_cell[1])*tmp;
787 // gamma[1] += xcm_cell[0]*xcm_cell[1]*tmp;
788 //
789 // *phi += (screening*log(rcm2)/PI - inv_area*rcm2) *mass*0.5;
790 // *kappa -= mass * inv_area;
791 //
792 // } // end of for
793 // } // end of if(tree->atLeaf())
794 //
795 //
796 //
797 // }else
798 // }
799 // }while(treeit.TreeWalkStep(allowDescent));
800 //
801
802 // Subtract off uniform mass sheet to compensate for the extra mass
803 // added to the universe in the halos.
804 //alpha[0] -= ray[0]*sigma_background*(inv_screening_scale2 == 0.0);
805 //alpha[1] -= ray[1]*sigma_background*(inv_screening_scale2 == 0.0);
806 //{ // taken out to speed up
807 // *kappa -= sigma_background;
808 // *phi -= sigma_background * sigma_background ;
809 //}
810
811 return;
812}
813
838//template <typename LensHaloType>
839//void TreeQuadHalos<LensHaloType>::force2D_recur(const PosType *ray,PosType *alpha,KappaType *kappa,KappaType *gamma,KappaType *phi){
840//
841//
842// alpha[0]=alpha[1]=gamma[0]=gamma[1]=gamma[2]=0.0;
843// *kappa=*phi=0.0;
844//
845// if(Nparticles == 0) return;
846// assert(tree);
847//
848// //walkTree_recur(tree->top,ray,&alpha[0],kappa,&gamma[0],phi);
849//
850// if(periodic_buffer){
851// PosType tmp_ray[2],tmp_c[2];
852//
853// // for points outside of primary zone shift to primary zone
854// //if(inbox(ray,tree->top->boundary_p1,tree->top->boundary_p2)){
855// tmp_c[0] = ray[0];
856// tmp_c[1] = ray[1];
857// /*}else{
858// int di[2];
859// PosType center[2] = { (tree->top->boundary_p1[0] + tree->top->boundary_p2[0])/2 ,
860// (tree->top->boundary_p1[1] + tree->top->boundary_p2[1])/2 };
861//
862// di[0] = (int)( 2*(ray[0] - center[0])/lx );
863// di[1] = (int)( 2*(ray[1] - center[1])/ly );
864//
865// di[0] = (di[0] > 0) ? (di[0] + 1)/2 : (di[0] - 1)/2;
866// di[1] = (di[1] > 0) ? (di[1] + 1)/2 : (di[1] - 1)/2;
867//
868// tmp_c[0] = ray[0] - lx * di[0];
869// tmp_c[1] = ray[1] - ly * di[1];
870//
871// assert(inbox(tmp_c,tree->top->boundary_p1,tree->top->boundary_p2));
872// }*/
873//
874// for(int i = -1;i<2;++i){
875// for(int j = -1;j<2;++j){
876// tmp_ray[0] = tmp_c[0] + i*original_xl;
877// tmp_ray[1] = tmp_c[1] + j*original_yl;
878// walkTree_recur(tree->top,tmp_ray,alpha,kappa,gamma,phi);
879// }
880// }
881// }else{
882// walkTree_recur(tree->top,ray,alpha,kappa,gamma,phi);
883// }
884//
885// // Subtract off uniform mass sheet to compensate for the extra mass
886// // added to the universe in the halos.
887// //alpha[0] -= ray[0]*sigma_background*(inv_screening_scale2 == 0);
888// //alpha[1] -= ray[1]*sigma_background*(inv_screening_scale2 == 0);
889//
890// //*kappa -= sigma_background;
891//
892// return;
893//}
894
895//template <typename LensHaloType>
896//void TreeQuadHalos<LensHaloType>::walkTree_recur(QBranchNB *branch,PosType const *ray,PosType *alpha,KappaType *kappa,KappaType *gamma, KappaType *phi){
897//
898// PosType xcm[2],rcm2cell,rcm2,tmp,boxsize2;
899// IndexType i;
900// std::size_t tmp_index;
901// PosType prefac;
902// /*bool notscreen = true;
903//
904// if(inv_screening_scale2 != 0) notscreen = BoxIntersectCircle(ray,3*sqrt(1.0/inv_screening_scale2), branch->boundary_p1, branch->boundary_p2);
905// */
906// if(branch->nparticles > 0)
907// {
908// xcm[0]=branch->center[0]-ray[0];
909// xcm[1]=branch->center[1]-ray[1];
910//
911// rcm2cell = xcm[0]*xcm[0] + xcm[1]*xcm[1];
912//
913// boxsize2 = (branch->boundary_p2[0]-branch->boundary_p1[0])*(branch->boundary_p2[0]-branch->boundary_p1[0]);
914//
915// if( rcm2cell < (branch->rcrit_angle)*(branch->rcrit_angle) || rcm2cell < 5.83*boxsize2)
916// {
917//
918// // Treat all particles in a leaf as a point particle
919// if(tree->atLeaf(branch))
920// {
921// Point_2d initgamma(gamma[0],gamma[1]);
922// for(i = 0 ; i < branch->nparticles ; ++i){
923//
924// xcm[0] = tree->xxp[branch->particles[i]][0] - ray[0];
925// xcm[1] = tree->xxp[branch->particles[i]][1] - ray[1];
926//
927// rcm2 = xcm[0]*xcm[0] + xcm[1]*xcm[1];
928// double screening = exp(-rcm2*inv_screening_scale2);
929// if(rcm2 < 1e-20) rcm2 = 1e-20;
930//
931// tmp_index = MultiMass*branch->particles[i];
932//
933// double mass = halos[tmp_index]->get_mass();
934// prefac = mass * (inv_area - screening/rcm2/PI);
935//
936// alpha[0] += prefac*xcm[0];
937// alpha[1] += prefac*xcm[1];
938//
939// //assert(abs(gamma[0]/2.74889e+15) < 10);
940// //assert(abs(gamma[1]/2.74889e+15) < 10);
941//
942// tmp = -2.0 * mass * screening/rcm2/rcm2/PI;
943// gamma[0] += 0.5*(xcm[0]*xcm[0]-xcm[1]*xcm[1])*tmp;
944// gamma[1] += xcm[0]*xcm[1]*tmp;
945//
946// //assert(abs(gamma[0]/2.74889e+15) < 10);
947// //assert(abs(gamma[1]/2.74889e+15) < 10);
948//
949// *phi += (log(rcm2)/PI - inv_area*rcm2)*0.5*mass;
950// *kappa -= mass*inv_area;
951// }
952// }
953//
954// // Find the particles that intersect with ray and add them individually.
955// if(rcm2cell < 5.83*boxsize2)
956// {
957//
958// for(i = 0 ; i < branch->Nbig_particles ; ++i)
959// {
960//
961// tmp_index = branch->big_particles[i];
962//
963// xcm[0] = tree->xxp[tmp_index][0] - ray[0];
964// xcm[1] = tree->xxp[tmp_index][1] - ray[1];
965// rcm2 = xcm[0]*xcm[0] + xcm[1]*xcm[1];
966// //assert(abs(gamma[0]/2.74889e+15) < 10);
967// //assert(abs(gamma[1]/2.74889e+15) < 10);
968//
969// /////////////////////////////////////////
970// double screening = exp(-rcm2*inv_screening_scale2);
971// halos[tmp_index]->force_halo(alpha,kappa,gamma,phi,xcm,true,screening);
972// //assert(abs(gamma[0]/2.74889e+15) < 10);
973// //assert(abs(gamma[1]/2.74889e+15) < 10);
974//
975// }
976// }
977//
978// if(branch->child0 != NULL)
979// walkTree_recur(branch->child0,&ray[0],&alpha[0],kappa,&gamma[0],&phi[0]);
980// if(branch->child1 != NULL)
981// walkTree_recur(branch->child1,&ray[0],&alpha[0],kappa,&gamma[0],&phi[0]);
982// if(branch->child2 != NULL)
983// walkTree_recur(branch->child2,&ray[0],&alpha[0],kappa,&gamma[0],&phi[0]);
984// if(branch->child3 != NULL)
985// walkTree_recur(branch->child3,&ray[0],&alpha[0],kappa,&gamma[0],&phi[0]);
986//
987// }
988// else
989// { // use whole cell
990//
991// double screening = exp(-rcm2cell*inv_screening_scale2);
992// //double tmp = -1.0*branch->mass/rcm2cell/PI*screening;
993// double tmp = branch->mass*(inv_area - screening/rcm2cell/PI);
994//
995// alpha[0] += tmp*xcm[0];
996// alpha[1] += tmp*xcm[1];
997//
998// { // taken out to speed up
999// assert(abs(gamma[0]/2.74889e+15) < 10);
1000// assert(abs(gamma[1]/2.74889e+15) < 10);
1001//
1002// tmp=-2.0*branch->mass/PI/rcm2cell/rcm2cell*screening;
1003// gamma[0] += 0.5*(xcm[0]*xcm[0]-xcm[1]*xcm[1])*tmp;
1004// gamma[1] += xcm[0]*xcm[1]*tmp;
1005//
1006// assert(abs(gamma[0]/2.74889e+15) < 10);
1007// assert(abs(gamma[1]/2.74889e+15) < 10);
1008//
1009// *phi += 0.5*tree->current->mass*log( rcm2cell )/PI*screening;
1010// *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;
1011// }
1012//
1013// // quadrapole contribution
1014// // the kappa and gamma are not calculated to this order
1015// alpha[0] -= (branch->quad[0]*xcm[0] + branch->quad[2]*xcm[1])
1016// /(rcm2cell*rcm2cell)/PI*screening;
1017// alpha[1] -= (branch->quad[1]*xcm[1] + branch->quad[2]*xcm[0])
1018// /(rcm2cell*rcm2cell)/PI*screening;
1019//
1020// tmp = 4*(branch->quad[0]*xcm[0]*xcm[0] + branch->quad[1]*xcm[1]*xcm[1]
1021// + 2*branch->quad[2]*xcm[0]*xcm[1])/(rcm2cell*rcm2cell*rcm2cell)/PI*screening;
1022//
1023// alpha[0] += tmp*xcm[0];
1024// alpha[1] += tmp*xcm[1];
1025//
1026// *kappa -= branch->mass * inv_area;
1027// *phi -= 0.5 * rcm2cell * inv_area * branch->mass;
1028// return;
1029// }
1030//
1031// }
1032// return;
1033//}
1034
1038template <typename LensHaloType>
1039void TreeQuadHalos<LensHaloType>::printParticlesInBranch(unsigned long number){
1040 unsigned long i;
1041
1042 tree->moveTop();
1043 do{
1044 if(tree->current->number == number){
1045 std::cout << tree->current->nparticles << std::endl;
1046 for(i=0;i<tree->current->nparticles;++i){
1047 std::cout << xp[tree->current->particles[i]][0] << " " << xp[tree->current->particles[i]][1] << std::endl;
1048 }
1049 return;
1050 }
1051 }while(tree->WalkStep(true));
1052
1053 return;
1054}
1060template <typename LensHaloType>
1061void TreeQuadHalos<LensHaloType>::printBranchs(int level){
1062
1063 bool decend = true;
1064 tree->moveTop();
1065 do{
1066 std::cout << tree->current->boundary_p1[0] << " " << tree->current->boundary_p1[1] << " "
1067 << tree->current->boundary_p2[0] << " " << tree->current->boundary_p2[1] << std::endl;
1068 if(tree->current->level == level) decend = false;
1069 else decend = true;
1070 }while(tree->WalkStep(decend));
1071
1072 return;
1073}
1074
1075#endif /* quadTreeHalos_h */
LensHaloType
LensHaloType
Type of halo profile.
Definition InputParams.h:44
LensPlaneTree
A LensPlane with a TreeQuad on it to calculate the deflection caused by field lenses.
Definition planes.h:35
QBiterator
A iterator class fore TreeStruct that allows for movement through the tree without changing anything ...
Definition qTreeNB.h:173
QBranchNB
Box representing a branch in a tree. It has four children. Used in QTreeNB which is used in TreeQuad.
Definition qTreeNB.h:31
QBranchNB::r2crit_angle
PosType r2crit_angle
the critical distance below which a branch is opened in the
Definition qTreeNB.h:84
QBranchNB::center
PosType center[2]
center of mass
Definition qTreeNB.h:56
QBranchNB::quad
PosType quad[3]
quadropole moment of branch
Definition qTreeNB.h:80
QBranchNB::level
int level
level in tree
Definition qTreeNB.h:59
QBranchNB::rmax
PosType rmax
largest dimension of box
Definition qTreeNB.h:82
QBranchNB::boundary_p1
PosType boundary_p1[2]
bottom, left, back corner of box
Definition qTreeNB.h:62
QBranchNB::particles
IndexType * particles
array of particles in QBranchNB
Definition qTreeNB.h:48
QBranchNB::Nbig_particles
IndexType Nbig_particles
the number of particles that aren't in children
Definition qTreeNB.h:52
QBranchNB::boundary_p2
PosType boundary_p2[2]
top, right, front corner of box
Definition qTreeNB.h:64
QTreeNB
QTreeNB: Tree structure used for force calculation with particles (i.e. stars, Nbody or halos).
Definition qTreeNB.h:98
QTreeNB::xxp
PType * xxp
Array of particle positions.
Definition qTreeNB.h:127
Generated on Sun Sep 8 2024 19:31:03 for GLAMERDOC++ by doxygen 1.12.0