GLAMERDOC++
Gravitational Lensing Code Library
particle_halo.h
1 //
2 // particle_halo.h
3 // GLAMER
4 //
5 // Created by bmetcalf on 16/06/15.
6 //
7 //
8 
9 #ifndef GLAMER_particle_halo_h
10 #define GLAMER_particle_halo_h
11 
12 #include "geometry.h"
13 #include "quadTree.h"
14 #include "simpleTree.h"
15 #include "particle_types.h"
16 #include "utilities_slsim.h"
17 #include "lens_halos.h"
18 #include "quadTreeHalos.h"
19 
39 template<typename PType>
40 class LensHaloParticles : public LensHalo
41 {
42 public:
43 
44  LensHaloParticles(const std::string& simulation_filename
45  ,SimFileFormat format
46  ,PosType redshift
47  ,int Nsmooth
48  ,const COSMOLOGY& cosmo
49  ,Point_2d theta_rotate
50  ,bool recenter
51  ,bool my_multimass
52  ,double inv_area
53  ,PosType MinPSize=0
54  ,PosType rescale_mass = 1.0
55  ,bool verbose=false
56  );
57 
58  LensHaloParticles(std::vector<PType> &pvector
59  ,float redshift
60  ,const COSMOLOGY& cosmo
61  ,Point_2d theta_rotate
62  ,bool recenter
63  ,double my_inv_area
64  ,float MinPSize = 0
65  ,double max_range = -1
66  ,bool verbose=false
67  ):LensHalo(redshift,cosmo),min_size(MinPSize),multimass(true),inv_area(my_inv_area)
68  {
69  Npoints = pvector.size();
70 
71  std::swap(pvector,trash_collector);
72  pp = trash_collector.data();
73 
74  set_up(redshift,cosmo,theta_rotate,max_range,recenter,verbose);
75  }
76 
78  LensHaloParticles(
79  PType *begin
80  ,size_t n
81  ,float redshift
82  ,const COSMOLOGY& cosmo
83  ,Point_2d theta_rotate
84  ,bool recenter
85  ,double my_inv_area
86  ,float MinPSize = 0
87  ,double max_range = -1
88  ,bool verbose=false
89  ):LensHalo(redshift,cosmo),min_size(MinPSize),multimass(true),inv_area(my_inv_area),Npoints(n)
90  {
91  pp = begin;
92  set_up(redshift,cosmo,theta_rotate,max_range,recenter,verbose);
93  }
94  ~LensHaloParticles();
95 
96  LensHaloParticles(LensHaloParticles &&h):LensHalo(std::move(h)){
97  mcenter = h.mcenter;
98  densest_point = h.densest_point;
99  trash_collector =std::move(h.trash_collector);
100  pp = trash_collector.data();
101  h.pp = nullptr;
102 
103  min_size = h.min_size;
104  multimass = h.multimass;
105  center = h.center;
106  Npoints = h.Npoints;
107  simfile = h.simfile;
108  sizefile = h.sizefile;
109  inv_area = h.inv_area;
110 
111  qtree = h.qtree;
112  h.qtree = nullptr;
113  }
114  LensHaloParticles<PType> & operator=(LensHaloParticles<PType> &&h);
115 
116 
117  void force_halo(double *alpha,KappaType *kappa,KappaType *gamma,KappaType *phi,double const *xcm
118  ,bool subtract_point=false,PosType screening = 1.0);
119 
120  size_t getN() const { return Npoints; };
121 
123  void rotate(Point_2d theta);
124 
126  Point_3d<> CenterOfMass(){return mcenter;}
128  Point_3d<> DensestPoint(){return densest_point;}
129 
134  static void makeSIE(
135  std::string new_filename
136  ,PosType redshift
137  ,double particle_mass
138  ,double total_mass
139  ,double sigma
140  ,double q
141  ,Utilities::RandomNumbers_NR &ran
142  );
143 
148  static LensHaloParticles<ParticleTypeSimple> SIE(
149  PosType redshift
150  ,double particle_mass
151  ,double total_mass
152  ,double sigma
153  ,double q
154  ,int Nneighbors
155  ,COSMOLOGY &cosmo
156  ,Utilities::RandomNumbers_NR &ran);
157 
158 
159  static void calculate_smoothing(int Nsmooth
160  ,PType *pp
161  ,size_t Npoints
162  ,bool verbose = false);
163 
164  void readPositionFileASCII(const std::string &filename);
165 
166  static void writeSizes(const std::string &filename
167  ,int Nsmooth
168  ,const PType *pp
169  ,size_t Npoints
170  );
171 
172  static bool readSizesFile(const std::string &filename
173  ,PType * pp
174  ,size_t Npoints
175  ,int Nsmooth
176  ,PosType min_size);
177 
178  friend class MakeParticleLenses;
179 
180 protected:
181  // constructure for derived classes
182  LensHaloParticles(float redshift
183  ,const COSMOLOGY& cosmo
184  ): LensHalo(redshift,cosmo){}
185 
186  // This constructor is really only for use by MakeParticleLenses. It does not take
187  // possession of the data and so they will not be deleted on destruction
188  LensHaloParticles(PType *pdata
189  ,size_t Nparticles
190  ,float redshift
191  ,const COSMOLOGY& cosmo
192  ,Point_2d theta_rotate
193  ,bool recenter
194  ,double my_inv_area
195  ,bool verbose = false
196  ):LensHalo(redshift,cosmo),pp(pdata),min_size(0),multimass(true),Npoints(Nparticles),inv_area(my_inv_area)
197  {
198  set_up(redshift,cosmo,theta_rotate,-1,recenter,verbose);
199  }
200 
201  void rotate_particles(PosType theta_x,PosType theta_y);
202  static void smooth_(TreeSimple<PType> *tree3d,PType *xp,size_t N,int Nsmooth);
203  void assignParams(InputParams& params);
204 
205  Point_3d<> mcenter;
206  Point_3d<> densest_point;
207 
208  PType *pp;
209  std::vector<PType> trash_collector;
210 
211  PosType min_size; // minimum allowed size
212  bool multimass;
213 
214  Utilities::Geometry::SphericalPoint<> center;
215 
216  size_t Npoints;
217  PosType inv_area;
218 
219  std::string simfile;
220  std::string sizefile;
221 
222 
223  TreeQuadParticles<PType> * qtree;
224  void set_up(float redshift,const COSMOLOGY& cosmo,Point_2d theta_rotate,double max_range,bool recenter,bool verbose);
225 };
226 
227 template<typename PType>
228 LensHaloParticles<PType>::LensHaloParticles(const std::string& simulation_filename
229  ,SimFileFormat format
230  ,PosType redshift
231  ,int Nsmooth
232  ,const COSMOLOGY& cosmo
233  ,Point_2d theta_rotate
234  ,bool recenter
235  ,bool my_multimass
236  ,double my_inv_area
237  ,PosType MinPSize
238  ,PosType massscaling
239  ,bool verbose
240  )
241 :LensHalo(redshift,cosmo),min_size(MinPSize),multimass(my_multimass),simfile(simulation_filename),inv_area(my_inv_area)
242 {
243 
244  LensHalo::setZlens(redshift,cosmo);
245  LensHalo::setCosmology(cosmo);
246  LensHalo::set_flag_elliptical(false);
247 
248  //stars_N = 0;
249  //stars_implanted = false;
250 
251  Rmax = 1.0e3;
252  LensHalo::setRsize(Rmax);
253 
254  switch (format) {
255  case SimFileFormat::ascii:
256  readPositionFileASCII(simulation_filename);
257  break;
258  default:
259  std::cerr << "LensHaloParticles does not accept format of particle data file." << std::endl;
260  throw std::invalid_argument("bad format");
261  }
262 
263  sizefile = simfile + "." + std::to_string(Nsmooth) + "sizes";
264 
265  if(!readSizesFile(sizefile,pp,Npoints,Nsmooth,min_size)){
266 
267  // calculate sizes
268  calculate_smoothing(Nsmooth,pp,Npoints);
269 
270  // save result to a file for future use
271  writeSizes(sizefile,Nsmooth,pp,Npoints);
272  //for(size_t i=0; i<Npoints ; ++i) if(sizes[i] < min_size) sizes[i] = min_size;
273  for(size_t i=0; i<Npoints ; ++i) if(pp[i].size() < min_size) pp[i].Size = min_size;
274  }
275 
276  double min_s = pp[0].Size;
277  size_t smallest_part = 0;
278  for(size_t i = 1 ; i<Npoints ; ++i){
279  if(pp[i].Size < min_s){
280  min_s = pp[i].Size;
281  smallest_part = i;
282  }
283  pp[i].Mass *= massscaling;
284  }
285 
286  densest_point[0] = pp[smallest_part].x[0];
287  densest_point[1] = pp[smallest_part].x[1];
288  densest_point[2] = pp[smallest_part].x[2];
289 
290 
291  // convert from comoving to physical coordinates
292  PosType scale_factor = 1/(1+redshift);
293  mcenter *= 0.0;
294  PosType max_mass = 0.0,min_mass = HUGE_VALF,mass=0;
295  for(size_t i=0;i<Npoints;++i){
296  pp[i][0] *= scale_factor;
297  pp[i][1] *= scale_factor;
298  pp[i][2] *= scale_factor;
299 
300  mcenter[0] += pp[i][0]*pp[multimass*i].mass();
301  mcenter[1] += pp[i][1]*pp[multimass*i].mass();
302  mcenter[2] += pp[i][2]*pp[multimass*i].mass();
303 
304  mass += pp[multimass*i].mass();
305 
306  max_mass = (pp[multimass*i].mass() > max_mass) ? pp[multimass*i].mass() : max_mass;
307  min_mass = (pp[multimass*i].mass() < min_mass) ? pp[multimass*i].mass() : min_mass;
308  }
309  LensHalo::setMass(mass);
310 
311  mcenter /= mass;
312 
313  if(verbose) std::cout << " Particle mass range : " << min_mass << " to " << max_mass << " ratio of : " << max_mass/min_mass << std::endl;
314 
315 
316  if(recenter){
317  PosType r2,r2max=0;
318  for(size_t i=0;i<Npoints;++i){
319  pp[i][0] -= mcenter[0];
320  pp[i][1] -= mcenter[1];
321  pp[i][2] -= mcenter[2];
322 
323  r2 = pp[i][0]*pp[i][0] + pp[i][1]*pp[i][1] + pp[i][2]*pp[i][2];
324  if(r2 > r2max) r2max = r2;
325  }
326 
327  LensHalo::setRsize( sqrt(r2max) );
328 
329  densest_point -= mcenter;
330 
331  mcenter *= 0;
332  }
333 
334  // rotate positions
335  rotate_particles(theta_rotate[0],theta_rotate[1]);
336 
337  qtree = new TreeQuadParticles<PType>(pp,Npoints,-1,-1,inv_area,20);
338 
339 }
340 
341 template<typename PType>
342 void LensHaloParticles<PType>::set_up(
343  float redshift
344  ,const COSMOLOGY& cosmo
345  ,Point_2d theta_rotate
346  ,double max_range
347  ,bool recenter
348  ,bool verbose
349 ){
350 
351  //LensHalo::setZlens(redshift);
352  //LensHalo::setCosmology(cosmo);
353  LensHalo::set_flag_elliptical(false);
354 
355  //stars_N = 0;
356  //stars_implanted = false;
357 
358  Rmax = 1.0e100;
359  LensHalo::setRsize(Rmax);
360 
361  double min_s = pp[0].Size;
362  size_t smallest_part = 0;
363  for(size_t i =0 ; i<Npoints ; ++i){
364  if(pp[i].Size < min_s){
365  min_s = pp[i].Size;
366  smallest_part = i;
367  }
368  }
369 
370  densest_point[0] = pp[smallest_part].x[0];
371  densest_point[1] = pp[smallest_part].x[1];
372  densest_point[2] = pp[smallest_part].x[2];
373 
374  // convert from comoving to physical coordinates
375  //PosType scale_factor = 1/(1+redshift);
376  mcenter *= 0.0;
377  PosType max_mass = 0.0,min_mass = HUGE_VALF,mass=0;
378  for(size_t i=0;i<Npoints;++i){
379 
380  mcenter[0] += pp[i][0]*pp[i].mass();
381  mcenter[1] += pp[i][1]*pp[i].mass();
382  mcenter[2] += pp[i][2]*pp[i].mass();
383 
384  mass += pp[i].mass();
385 
386  max_mass = (pp[i].mass() > max_mass) ? pp[i].mass() : max_mass;
387  min_mass = (pp[i].mass() < min_mass) ? pp[i].mass() : min_mass;
388  }
389  LensHalo::setMass(mass);
390 
391  mcenter /= mass;
392 
393  if(verbose) std::cout << " Particle mass range : " << min_mass << " to " << max_mass << " ratio of : " << max_mass/min_mass << std::endl;
394 
395  if(recenter){
396  PosType r2,r2max=0;
397  for(size_t i=0;i<Npoints;++i){
398  pp[i][0] -= mcenter[0];
399  pp[i][1] -= mcenter[1];
400  pp[i][2] -= mcenter[2];
401 
402  r2 = pp[i][0]*pp[i][0] + pp[i][1]*pp[i][1] + pp[i][2]*pp[i][2];
403  if(r2 > r2max) r2max = r2;
404 
405  densest_point -= mcenter;
406 
407  mcenter *= 0;
408  }
409 
410  LensHalo::setRsize( sqrt(r2max) );
411  }
412 
413  // rotate positions
414  rotate_particles(theta_rotate[0],theta_rotate[1]);
415 
416  qtree = new TreeQuadParticles<PType>(pp,Npoints,-1,-1,inv_area,20,0.1,false,max_range);
417 }
418 
419 
420 
421 template<typename PType>
422 LensHaloParticles<PType>::~LensHaloParticles(){
423  delete qtree;
424 }
425 
426 template<typename PType>
427 void LensHaloParticles<PType>::force_halo(double *alpha,KappaType *kappa,KappaType *gamma,KappaType *phi
428  ,double const *xcm
429  ,bool subtract_point,PosType screening){
430  qtree->force2D_recur(xcm,alpha,kappa,gamma,phi);
431 
432  alpha[0] *= -1;
433  alpha[1] *= -1;
434 }
435 
437 template<typename PType>
438 void LensHaloParticles<PType>::rotate(Point_2d theta){
439  rotate_particles(theta[0],theta[1]);
440  delete qtree;
441 
442  qtree = new TreeQuadParticles<ParticleType<float> >(pp,Npoints,-1,-1,inv_area,20);
443 }
444 
452 template<typename PType>
453 void LensHaloParticles<PType>::readPositionFileASCII(const std::string &filename
454  ){
455 
456  int ncoll = Utilities::IO::CountColumns(filename);
457  if(!multimass && ncoll != 3 ){
458  std::cerr << filename << " should have three columns!" << std::endl;
459  }
460  if(multimass && ncoll != 4 ){
461  std::cerr << filename << " should have four columns!" << std::endl;
462  }
463 
464  std::ifstream myfile(filename);
465 
466  //size_t Npoints = 0;
467 
468  // find number of particles
469 
470  if (myfile.is_open()){
471 
472  float tmp_mass = 0.0;
473  std::string str,label;
474  int count =0;
475  while(std::getline(myfile, str)){
476  std::stringstream ss(str);
477  ss >> label;
478  if(label == "#"){
479  ss >> label;
480  if(label == "nparticles"){
481  ss >> Npoints;
482  ++count;
483  }
484  if(!multimass){
485  if(label == "mass"){
486  ss >> tmp_mass;
487  ++count;
488  }
489  }
490  }else break;
491  if(multimass && count == 1 ) break;
492  if(!multimass && count == 2 ) break;
493  }
494 
495  if(count == 0){
496  if(multimass) std::cerr << "File " << filename << " must have the header lines: " << std::endl
497  << "# nparticles ****" << std::endl << "# mass ****" << std::endl;
498  if(!multimass) std::cerr << "File " << filename << " must have the header lines: " << std::endl
499  << "# nparticles ****" << std::endl;
500  throw std::runtime_error("file reading error");
501  }
502 
503  trash_collector.resize(Npoints); // this is here just to make sure this is deleted
504  pp = trash_collector.data();
505  //xp = Utilities::PosTypeMatrix(Npoints,3);
506  //if(multimass) masses.resize(Npoints);
507  //else masses.push_back(tmp_mass);
508 
509  size_t row = 0;
510 
511  // read in particle positions
512  if(!multimass){
513  while(std::getline(myfile, str) && row < Npoints){
514  if(str[0] == '#') continue; //for comments
515  std::stringstream ss(str);
516 
517  ss >> pp[row][0];
518  if(!(ss >> pp[row][1])) std::cerr << "3 columns are expected in line " << row
519  << " of " << filename << std::endl;
520  if(!(ss >> pp[row][2])) std::cerr << "3 columns are expected in line " << row
521  << " of " << filename << std::endl;
522 
523  row++;
524  }
525  }else{
526  while(std::getline(myfile, str) && row < Npoints){
527  if(str[0] == '#') continue; //for comments
528  std::stringstream ss(str);
529 
530  ss >> pp[row][0];
531  if(!(ss >> pp[row][1])) std::cerr << "4 columns are expected in line " << row
532  << " of " << filename << std::endl;
533  if(!(ss >> pp[row][2])) std::cerr << "4 columns are expected in line " << row
534  << " of " << filename << std::endl;
535  if(!(ss >> pp[row].Mass)) std::cerr << "4 columns are expected in line " << row
536  << " of " << filename << std::endl;
537 
538  row++;
539  }
540  }
541 
542  if(row != Npoints){
543  std::cerr << "Number of data rows in " << filename << " does not match expected number of particles."
544  << std::endl;
545  throw std::runtime_error("file reading error");
546  }
547  }else{
548  std::cerr << "Unable to open file " << filename << std::endl;
549  throw std::runtime_error("file reading error");
550  }
551 
552  std::cout << Npoints << " particle positions read from file " << filename << std::endl;
553 
554 }
555 
556 template<typename PType>
557 bool LensHaloParticles<PType>::readSizesFile(const std::string &filename
558  ,PType * pp
559  ,size_t Npoints
560  ,int Nsmooth
561  ,PosType min_size){
562 
563  std::ifstream myfile(filename);
564 
565  // find number of particles
566 
567  PosType min=HUGE_VALF,max=0.0;
568  if (myfile.is_open()){
569 
570  std::string str,label;
571  int count =0;
572  size_t Ntmp;
573  int NStmp;
574  while(std::getline(myfile, str)){
575  std::stringstream ss(str);
576  ss >> label;
577  if(label == "#"){
578  ss >> label;
579  if(label == "nparticles"){
580  ss >> Ntmp;
581  if(Ntmp != Npoints){
582  std::cerr << "Number of particles in " << filename << " does not match expected number" << std::endl;
583  throw std::runtime_error("file reading error");
584  }
585  ++count;
586  }
587  if(label == "nsmooth"){
588  ss >> NStmp;
589  if(NStmp != Nsmooth) return false;
590  ++count;
591  }
592 
593  }else break;
594  if(count == 2) break;
595  }
596 
597  if(count != 2){
598  std::cerr << "File " << filename << " must have the header lines: " << std::endl
599  << "# nparticles ****" << std::endl;
600  throw std::runtime_error("file reading error");
601  }
602 
603  size_t row = 0;
604 
605  std::cout << "reading in particle sizes from " << filename << "..." << std::endl;
606 
607  // read in particle sizes
608  while(std::getline(myfile, str)){
609  if(str[0] == '#') continue; //for comments
610  std::stringstream ss(str);
611 
612  ss >> pp[row].Size;
613  if(min_size > pp[row].size() ) pp[row].Size = min_size;
614  min = min < pp[row].size() ? min : pp[row].size();
615  max = max > pp[row].size() ? max : pp[row].size();
616  row++;
617  }
618 
619  if(row != Npoints){
620  std::cerr << "Number of data rows in " << filename << " does not match expected number of particles."
621  << std::endl;
622  throw std::runtime_error("file reading error");
623  }
624  }else{
625  return false;
626  }
627 
628  std::cout << Npoints << " particle sizes read from file " << filename << std::endl;
629  std::cout << " maximun particle sizes " << max << " minimum " << min << " Mpc" << std::endl;
630 
631  return true;
632 }
633 
634 template<typename PType>
635 void LensHaloParticles<PType>::rotate_particles(PosType theta_x,PosType theta_y){
636 
637  if(theta_x == 0.0 && theta_y == 0.0) return;
638 
639  PosType coord[3][3];
640  PosType cx,cy,sx,sy;
641 
642  cx = cos(theta_x); sx = sin(theta_x);
643  cy = cos(theta_y); sy = sin(theta_y);
644 
645  coord[0][0] = cy; coord[1][0] = -sy*sx; coord[2][0] = cx;
646  coord[0][1] = 0; coord[1][1] = cx; coord[2][1] = sx;
647  coord[0][2] = -sy; coord[1][2] = -cy*sx; coord[2][2] = cy*cx;
648 
649  //PosType tmp[3];
650  Point_3d<> tmp;
651  int j;
652  // rotate particle positions */
653  for(size_t i=0;i<Npoints;++i){
654  tmp *= 0;
655  for(j=0;j<3;++j){
656  tmp[0] += coord[0][j]*pp[i][j];
657  tmp[1] += coord[1][j]*pp[i][j];
658  tmp[2] += coord[2][j]*pp[i][j];
659  }
660  for(j=0;j<3;++j) pp[i][j]=tmp[j];
661  }
662 
663  // rotate center of mass
664  tmp *= 0.0;
665  for(j=0;j<3;++j){
666  tmp[0] += coord[0][j]*mcenter[j];
667  tmp[1] += coord[1][j]*mcenter[j];
668  tmp[2] += coord[2][j]*mcenter[j];
669  }
670  mcenter = tmp;
671 
672  // rotate center of mass
673  tmp *= 0.0;
674  for(j=0;j<3;++j){
675  tmp[0] += coord[0][j]*densest_point[j];
676  tmp[1] += coord[1][j]*densest_point[j];
677  tmp[2] += coord[2][j]*densest_point[j];
678  }
679  densest_point = tmp;
680 }
681 
682 template<typename PType>
683 void LensHaloParticles<PType>::calculate_smoothing(int Nsmooth,PType *pp
684  ,size_t Npoints
685  ,bool verbose
686  ){
687 
688  int nthreads = Utilities::GetNThreads();
689 
690  if(verbose) std::cout << "Calculating smoothing of particles ..." << std::endl
691  << Nsmooth << " neighbors. If there are a lot of particles this could take a while." << std::endl;
692 
693  time_t to,t;
694  time(&to);
695 
696  // make 3d tree of particle postions
697  TreeSimple<PType> tree3d(pp,Npoints,2*Nsmooth,3,true);
698  // find distance to nth neighbour for every particle
699  if(Npoints < 1000){
700  //IndexType neighbors[Nsmooth];
701  for(size_t i=0;i<Npoints;++i){
702  pp[i].Size = tree3d.NNDistance(&pp[i][0],Nsmooth + 1);
703  }
704  }else{
705  size_t chunksize = Npoints/nthreads;
706  std::vector<std::thread> thr(nthreads);
707 
708  size_t N;
709  for(int ii = 0; ii < nthreads ;++ii){
710  if(ii == nthreads - 1){
711  N = Npoints - ii*chunksize;
712  }else N = chunksize;
713 
714  thr[ii] = std::thread(LensHaloParticles<PType>::smooth_,&tree3d
715  ,&(pp[ii*chunksize]),N,Nsmooth);
716  }
717  for(int ii = 0; ii < nthreads ;++ii) thr[ii].join();
718  }
719  time(&t);
720  if(verbose) std::cout << "done in " << difftime(t,to) << " secs" << std::endl;
721 }
722 
723 template<typename PType>
724 void LensHaloParticles<PType>::smooth_(TreeSimple<PType> *tree3d,PType *pp,size_t N,int Nsmooth){
725 
726  //IndexType neighbors[Nsmooth];
727  for(size_t i=0;i<N;++i){
728  pp[i].Size = tree3d->NNDistance(&(pp[i][0]),Nsmooth + 1);
729  }
730 }
731 
732 
733 template<typename PType>
734 void LensHaloParticles<PType>::writeSizes(const std::string &filename,int Nsmooth
735  ,const PType *pp,size_t Npoints
736  ){
737 
738  std::ofstream myfile(filename);
739 
740  // find number of particles
741 
742  if (myfile.is_open()){
743 
744  std::cout << "Writing particle size information to file " << filename << " ...." << std::endl;
745 
746  myfile << "# nparticles " << Npoints << std::endl;
747  myfile << "# nsmooth " << Nsmooth << std::endl;
748  for(size_t i=0;i<Npoints;++i){
749  myfile << pp[i].Size << std::endl;
750  if(!myfile){
751  std::cerr << "Unable to write to file " << filename << std::endl;
752  throw std::runtime_error("file writing error");
753  }
754  }
755 
756  std::cout << "done" << std::endl;
757 
758  }else{
759  std::cerr << "Unable to write to file " << filename << std::endl;
760  throw std::runtime_error("file writing error");
761  }
762 }
763 
764 template<typename PType>
765 void LensHaloParticles<PType>::makeSIE(
766  std::string new_filename
767  ,PosType redshift
768  ,double particle_mass
769  ,double total_mass
770  ,double sigma
771  ,double q
772  ,Utilities::RandomNumbers_NR &ran
773  ){
774 
775  size_t Npoints = total_mass/particle_mass;
776  PosType Rmax = (1+redshift)*total_mass*Grav*lightspeed*lightspeed/sigma/sigma/2;
777  Point_3d<> point;
778  double qq = sqrt(q);
779 
780  std::ofstream datafile;
781  datafile.open(new_filename);
782 
783  datafile << "# nparticles " << Npoints << std::endl;
784  datafile << "# mass " << particle_mass << std::endl;
785  // create particles
786  for(size_t i=0; i< Npoints ;++i){
787  point[0] = ran.gauss();
788  point[1] = ran.gauss();
789  point[2] = ran.gauss();
790 
791  point *= Rmax*ran()/point.length();
792 
793  point[0] *= qq;
794  point[1] /= qq;
795 
796  datafile << point[0] << " " << point[1] << " "
797  << point[2] << " " << std::endl;
798 
799  }
800 
801  datafile.close();
802 }
803 
804  template<typename PType>
805  LensHaloParticles<ParticleTypeSimple> LensHaloParticles<PType>::SIE(
806  PosType redshift
807  ,double particle_mass
808  ,double total_mass
809  ,double sigma
810  ,double q
811  ,int Nneighbors
812  ,COSMOLOGY &cosmo
813  ,Utilities::RandomNumbers_NR &ran
814  ){
815 
816  size_t Npoints = total_mass/particle_mass;
817  PosType Rmax = total_mass*Grav*lightspeed*lightspeed/sigma/sigma/2;
818 
819  double qq = sqrt(q);
820 
821  std::cout << "making SIE LensHaloParticles..." << std::endl;
822  std::cout << "# nparticles " << Npoints << std::endl;
823  std::cout << "# mass " << particle_mass << std::endl;
824 
825  std::vector<ParticleTypeSimple> parts(Npoints);
826  // create particles
827  for(auto &p : parts){
828  p[0] = ran.gauss();
829  p[1] = ran.gauss();
830  p[2] = ran.gauss();
831 
832  double s = Rmax*ran()/sqrt(p[0]*p[0] + p[1]*p[1] + p[2]*p[2] );
833  p[0] *= s;
834  p[1] *= s;
835  p[2] *= s;
836 
837  p[0] *= qq;
838  p[1] /= qq;
839 
840  p.Mass = particle_mass;
841  }
842 
843  LensHaloParticles<ParticleTypeSimple> ::calculate_smoothing(Nneighbors,parts.data(),Npoints);
844 
845  return LensHaloParticles<ParticleTypeSimple>(parts,redshift
846  ,cosmo,Point_2d(0,0),false,0);
847  }
848 
849 template<typename PType>
850 LensHaloParticles<PType> & LensHaloParticles<PType>::operator=(LensHaloParticles<PType> &&h){
851  if(this == &h) return *this;
852  LensHalo::operator=(std::move(h));
853  mcenter = h.mcenter;
854  densest_point = h.densest_point;
855  trash_collector =std::move(h.trash_collector);
856  pp = h.pp; // note: this depends on std:move keeping the pointer valid if it is constructed with the public constructor
857  h.pp = nullptr;
858 
859  min_size = h.min_size;
860  multimass = h.multimass;
861  center = h.center;
862  Npoints = h.Npoints;
863  simfile = h.simfile;
864  sizefile = h.sizefile;
865  inv_area = h.inv_area;
866 
867  qtree = h.qtree;
868  h.qtree = nullptr;
869 
870  return *this;
871 }
872 
873 
964 class MakeParticleLenses{
965 
966 public:
967 
969  std::vector<LensHaloParticles<ParticleType<float> > *> halos;
970 
972  std::vector<size_t> getnp(){return nparticles;}
973 
975  std::vector<float> getmp(){return masses;}
976 
978  double sim_redshift(){return z_original;}
979 
980  MakeParticleLenses(const std::string &filename
981  ,SimFileFormat format
982  ,int Nsmooth
983  ,bool recenter
984  //,bool compensate=false /// if true a negative mass will be included so that the region wil have zero massl
985  ,bool ignore_type_in_smoothing = false
986  );
987 
988  MakeParticleLenses(const std::string &filename
989  ,bool recenter
990  //,bool compensate=false /// if true a negative mass will be included so that the region wil have zero massl
991  );
992 
993  ~MakeParticleLenses(){
994  for(auto p : halos) delete p;
995  }
996 
998 
999  void Recenter(Point_3d<> x);
1000 
1001  void CreateHalos(const COSMOLOGY &cosmo,double redshift,double inv_area);
1002 
1004  void radialCut(Point_3d<> center,double radius);
1006  void cylindricalCut(Point_2d center,double radius);
1007 
1009  Point_3d<> getCenterOfMass() const{return cm;}
1010 
1012  Point_3d<> densest_particle() const;
1013 
1015  void getBoundingBox(Point_3d<> &Xmin,Point_3d<> &Xmax) const{
1016  Xmin = bbox_ll;
1017  Xmax = bbox_ur;
1018  }
1019 
1020  double getZoriginal(){return z_original;}
1021 
1022  std::vector<ParticleType<float> > data;
1023 private:
1024  const std::string filename;
1025  int Nsmooth;
1026  //bool compensate;
1027 
1028  Point_3d<> bbox_ll; // minumum coordinate values of particles
1029  Point_3d<> bbox_ur; // maximim coordinate values of particles
1030 
1031  double z_original = -1;
1032  std::vector<size_t> nparticles;
1033  std::vector<float> masses;
1034  Point_3d<> cm;
1035 
1036  // write a glamB format file with all required particle data
1037  static void writeSizesB(const std::string &filename
1038  ,std::vector<ParticleType<float> > &pv
1039  ,int Nsmooth,std::vector<size_t> numbytype,double redshift){
1040 
1041  assert(numbytype.size() == 6);
1042  size_t ntot = pv.size();
1043  std::ofstream myfile(filename, std::ios::out | std::ios::binary);
1044  if(!myfile.write((char*)&Nsmooth,sizeof(int))){
1045  std::cerr << "Unable to write to file " << filename << std::endl;
1046  throw std::runtime_error("file writing error");
1047  }else{
1048  myfile.write((char*)&ntot,sizeof(size_t));
1049  myfile.write((char*)(numbytype.data()),sizeof(size_t)*6);
1050  myfile.write((char*)&redshift,sizeof(double));
1051  myfile.write((char*)(pv.data()),sizeof(ParticleType<float>)*ntot);
1052  }
1053  }
1054 
1055  // read a glamB format file
1056  static bool readSizesB(const std::string &filename
1057  ,std::vector<ParticleType<float> > &pv
1058  ,int &Nsmooth,std::vector<size_t> &numbytype,double &redshift){
1059 
1060  numbytype.resize(6);
1061  size_t ntot = pv.size();
1062  std::ifstream myfile(filename, std::ios::in | std::ios::binary);
1063  if(!myfile.read((char*)&Nsmooth,sizeof(int))){
1064  return false;
1065  }else{
1066  myfile.read((char*)&ntot,sizeof(size_t));
1067  myfile.read((char*)(numbytype.data()),sizeof(size_t)*6);
1068  myfile.read((char*)&redshift,sizeof(double));
1069  pv.resize(ntot);
1070  myfile.read((char*)(pv.data()),sizeof(ParticleType<float>)*ntot);
1071  }
1072  return true;
1073  }
1074 
1075  // Read particle from Gadget-2 format file
1076  bool readGadget2(bool ignore_type);
1077 
1078  // Reads particles from first 4 columns of csv file
1079  bool readCSV(int columns_used);
1080 
1081 #ifdef ENABLE_HDF5
1082  bool readHDF5();
1083 #endif
1084 
1085 };
1086 
1091 template<typename HType>
1092 class LensHaloHalos : public LensHalo
1093 {
1094 public:
1095 
1096  LensHaloHalos(std::vector<HType> &pvector
1097  ,float redshift
1098  ,const COSMOLOGY& cosmo
1099  ,float inv_area
1100  ,float my_theta_force=0.1
1101  ,bool mu_periodic_buffer = false
1102  ,bool verbose=false
1103  ):LensHalo(redshift,cosmo),inv_area(inv_area),theta_force(my_theta_force),periodic_buffer(mu_periodic_buffer)
1104  {
1105  std::swap(pvector,trash_collector);
1106  Nhalos = trash_collector.size();
1107  vpp.resize(Nhalos);
1108 
1109  for(size_t ii = 0 ; ii < Nhalos ; ++ii) vpp[ii] = &trash_collector[ii];
1110  set_up(redshift,cosmo,verbose);
1111  }
1112  ~LensHaloHalos();
1113 
1114  LensHaloHalos(LensHaloHalos &&h):LensHalo(std::move(h)){
1115  mcenter = h.mcenter;
1116  trash_collector =std::move(h.trash_collector);
1117  vpp = std::move(h.vpp);
1118 
1119  center = h.center;
1120  Nhalos = h.Nhalos;
1121 
1122  qtree = h.qtree;
1123  h.qtree = nullptr;
1124 
1125  inv_area = h.inv_area;
1126  }
1127  LensHaloHalos & operator=(LensHaloHalos &&h);
1128 
1129 
1130  void force_halo(double *alpha,KappaType *kappa,KappaType *gamma,KappaType *phi
1131  ,double const *xcm
1132  ,bool subtract_point=false,PosType screening = 1.0);
1133 
1134  size_t getN() const { return Nhalos; };
1135 
1137  Point_2d CenterOfMass(){return mcenter;}
1138 
1139 protected:
1140  Point_2d mcenter;
1141  double inv_area;
1142 
1143  std::vector<HType *> vpp;
1144  std::vector<HType> trash_collector;
1145 
1146  Utilities::Geometry::SphericalPoint<> center;
1147 
1148  size_t Nhalos;
1149 
1150  float theta_force;
1151  bool periodic_buffer;
1152 
1153  TreeQuadHalos<HType> * qtree;
1154 
1155  void set_up(float redshift,const COSMOLOGY& cosmo,bool verbose);
1156 };
1157 
1158 template<typename HType>
1159 void LensHaloHalos<HType>::set_up(
1160  float redshift
1161  ,const COSMOLOGY& cosmo
1162  ,bool verbose
1163 ){
1164 
1165  LensHalo::set_flag_elliptical(false);
1166 
1167  Rmax = 1.0e100; // ????
1168  LensHalo::setRsize(Rmax);
1169 
1170  // convert from comoving to physical coordinates
1171  //PosType scale_factor = 1/(1+redshift);
1172  mcenter *= 0.0;
1173  PosType max_mass = 0.0,min_mass = HUGE_VALF,mass=0;
1174 
1175  for(HType &h : trash_collector){
1176 
1177  mcenter[0] += h[0]*h.get_mass();
1178  mcenter[1] += h[1]*h.get_mass();
1179 
1180  mass += h.get_mass();
1181 
1182  max_mass = (h.get_mass() > max_mass) ? h.get_mass() : max_mass;
1183  min_mass = (h.get_mass() < min_mass) ? h.get_mass() : min_mass;
1184  }
1185  LensHalo::setMass(mass);
1186 
1187  mcenter /= mass;
1188 
1189  if(verbose) std::cout << " Particle mass range : " << min_mass << " to " << max_mass << " ratio of : " << max_mass/min_mass << std::endl;
1190 
1191  qtree = new TreeQuadHalos<HType>(vpp.data(),Nhalos,inv_area,5,theta_force,periodic_buffer);
1192 }
1193 
1194 template<typename HType>
1195 LensHaloHalos<HType>::~LensHaloHalos(){
1196  delete qtree;
1197 }
1198 
1199 template<typename HType>
1200 void LensHaloHalos<HType>::force_halo(double *alpha,KappaType *kappa
1201  ,KappaType *gamma,KappaType *phi
1202  ,double const *xcm
1203  ,bool subtract_point,PosType screening){
1204  qtree->force2D(xcm,alpha,kappa,gamma,phi);
1205 
1206  // ?????
1207  alpha[0] *= -1;
1208  alpha[1] *= -1;
1209 }
1210 
1211 template<typename HType>
1212 LensHaloHalos<HType> & LensHaloHalos<HType>::operator=(LensHaloHalos<HType> &&h){
1213  if(this == &h) return *this;
1214  LensHalo::operator=(std::move(h));
1215  mcenter = h.mcenter;
1216 
1217  trash_collector =std::move(h.trash_collector);
1218  vpp = std::move(h.vpp); // note: this depends on std:move keeping the pointer valid if it is constructed with the public constructor
1219 
1220  center = h.center;
1221  Nhalos = h.Nhalos;
1222 
1223  qtree = h.qtree;
1224  h.qtree = nullptr;
1225 
1226  return *this;
1227 }
1228 
1229 #endif
MakeParticleLenses::cylindricalCut
void cylindricalCut(Point_2d center, double radius)
remove particles that are beyond cylindrical radius (Mpc/h) of center
Definition: particle_lens.cpp:497
LensHalo::operator=
LensHalo & operator=(const LensHalo &h)
Definition: lens_halos.cpp:169
to_string
std::string to_string(const Band &band)
convert a Band type to a string name.
Definition: InputParams.cpp:878
TreeSimple
A C++ class wrapper for the bianary treeNB used in the Nbody force calculation, but also useful for g...
Definition: simpleTree.h:86
MakeParticleLenses::radialCut
void radialCut(Point_3d<> center, double radius)
remove particles that are beyond radius (Mpc/h) from center
Definition: particle_lens.cpp:447
LensHalo::force_halo
virtual void force_halo(PosType *alpha, KappaType *kappa, KappaType *gamma, KappaType *phi, PosType const *xcm, bool subtract_point=false, PosType screening=1.0)
Definition: lens_halos.cpp:1528
MakeParticleLenses::halos
std::vector< LensHaloParticles< ParticleType< float > > * > halos
vector of LensHalos, one for each type of particle type
Definition: particle_halo.h:969
LensHalo::setZlens
void setZlens(PosType my_zlens, const COSMOLOGY &cosmo)
set redshift
Definition: lens_halos.h:137
MakeParticleLenses::MakeParticleLenses
MakeParticleLenses(const std::string &filename, SimFileFormat format, int Nsmooth, bool recenter, bool ignore_type_in_smoothing=false)
Definition: particle_lens.cpp:22
MakeParticleLenses::getCenterOfMass
Point_3d getCenterOfMass() const
returns the original center of mass of all the particles
Definition: particle_halo.h:1009
LensHaloParticles::LensHaloParticles
LensHaloParticles(std::vector< PType > &pvector, float redshift, const COSMOLOGY &cosmo, Point_2d theta_rotate, bool recenter, double my_inv_area, float MinPSize=0, double max_range=-1, bool verbose=false)
Definition: particle_halo.h:58
SimFileFormat::ascii
@ ascii
the original ascii GLAMER format.
LensHaloParticles::LensHaloParticles
LensHaloParticles(const std::string &simulation_filename, SimFileFormat format, PosType redshift, int Nsmooth, const COSMOLOGY &cosmo, Point_2d theta_rotate, bool recenter, bool my_multimass, double inv_area, PosType MinPSize=0, PosType rescale_mass=1.0, bool verbose=false)
Definition: particle_halo.h:228
MakeParticleLenses::sim_redshift
double sim_redshift()
returns original redshift of snapshot, redshifts of the halos can be changed
Definition: particle_halo.h:978
LensHaloParticles::set_up
void set_up(float redshift, const COSMOLOGY &cosmo, Point_2d theta_rotate, double max_range, bool recenter, bool verbose)
Definition: particle_halo.h:342
LensHalo
A base class for all types of lensing "halos" which are any mass distribution that cause lensing.
Definition: lens_halos.h:56
LensHaloParticles::LensHaloParticles
LensHaloParticles(PType *begin, size_t n, float redshift, const COSMOLOGY &cosmo, Point_2d theta_rotate, bool recenter, double my_inv_area, float MinPSize=0, double max_range=-1, bool verbose=false)
this constructor does not take possession of the particles
Definition: particle_halo.h:78
TreeQuadParticles
TreeQuadParticles is a class for calculating the deflection, kappa and gamma by tree method.
Definition: quadTree.h:47
Point_2d
Class for representing points or vectors in 2 dimensions. Not that the dereferencing operator is over...
Definition: point.h:48
Utilities::RandomNumbers_NR
This is a class for generating random numbers. It simplifies and fool proofs initialization and allow...
Definition: utilities_slsim.h:1041
MakeParticleLenses::getnp
std::vector< size_t > getnp()
returns number of particles of each type
Definition: particle_halo.h:972
Point_3d::length
T length() const
length
Definition: point.h:973
LensHaloParticles::DensestPoint
Point_3d DensestPoint()
get the densistt point in input coordinates
Definition: particle_halo.h:128
MakeParticleLenses::densest_particle
Point_3d densest_particle() const
returns the location of the densest particle in (Mpc/h)
Definition: particle_lens.cpp:480
MakeParticleLenses::Recenter
void Recenter(Point_3d<> x)
recenter the particles to 3d point in physical Mpc/h units If the halos have already been created the...
Definition: particle_lens.cpp:139
MakeParticleLenses
A class for constructing LensHalos from particles in a data file.
Definition: particle_halo.h:964
LensHaloParticles::LensHaloParticles
LensHaloParticles(float redshift, const COSMOLOGY &cosmo)
Definition: particle_halo.h:182
LensHaloParticles::readPositionFileASCII
void readPositionFileASCII(const std::string &filename)
Reads number of particle and particle positons into Npoint and xp from a ASCII file.
Definition: particle_halo.h:453
Utilities::Geometry::SphericalPoint
represents a point in spherical coordinates, theta = 0 is equator
Definition: geometry.h:30
Utilities::RandomNumbers_NR::gauss
PosType gauss()
generates a Gaussian distributed number with unit variance by polar Box-Muller transform
Definition: utilities_slsim.h:1050
LensHaloParticles::SIE
static LensHaloParticles< ParticleTypeSimple > SIE(PosType redshift, double particle_mass, double total_mass, double sigma, double q, int Nneighbors, COSMOLOGY &cosmo, Utilities::RandomNumbers_NR &ran)
This is a test static function that makes a truncated SIE out of particles and puts it into a file in...
Definition: particle_halo.h:805
MakeParticleLenses::getmp
std::vector< float > getmp()
returns mass of particles of each type, if 0 they can have different masses
Definition: particle_halo.h:975
TreeSimple::NNDistance
T NNDistance(T *ray, int Nneighbors) const
Finds the nearest N neighbors and puts their index numbers in an array, also returns the distance to ...
Definition: simpleTree.h:483
LensHaloParticles::LensHaloParticles
LensHaloParticles(PType *pdata, size_t Nparticles, float redshift, const COSMOLOGY &cosmo, Point_2d theta_rotate, bool recenter, double my_inv_area, bool verbose=false)
Definition: particle_halo.h:188
LensHaloParticles
A class that represents the lensing by a collection of simulation particles.
Definition: particle_halo.h:41
MakeParticleLenses::getBoundingBox
void getBoundingBox(Point_3d<> &Xmin, Point_3d<> &Xmax) const
return the maximum and minimum coordinates of the particles in each dimension in for the original sim...
Definition: particle_halo.h:1015
LensHalo::setCosmology
virtual void setCosmology(const COSMOLOGY &cosmo)
used for elliptical NFWs only, in that case get_switch_flag==true
Definition: lens_halos.h:163
Utilities::IO::CountColumns
int CountColumns(std::string filename, char comment_char='#', char deliniator=' ')
Count the number of columns in a ASCII data file.
Definition: utilities.cpp:907
LensHaloParticles::makeSIE
static void makeSIE(std::string new_filename, PosType redshift, double particle_mass, double total_mass, double sigma, double q, Utilities::RandomNumbers_NR &ran)
This is a test static function that makes a truncated SIE out of particles and puts it into a file in...
Definition: particle_halo.h:765
Point_3d<>
particle_types.h
LensHaloParticles::CenterOfMass
Point_3d CenterOfMass()
get current center of mass in input coordinates
Definition: particle_halo.h:126
Utilities::GetNThreads
int GetNThreads()
returns the compiler variable N_THREADS that is maximum number of threads to be used.
Definition: utilities.cpp:846
LensHaloParticles::rotate
void rotate(Point_2d theta)
rotate the simulation around the origin of the simulation coordinates, (radians)
Definition: particle_halo.h:438
InputParams
Structure for reading and writing parameters to and from a parameter file as well as a container for ...
Definition: InputParams.h:99
MakeParticleLenses::CreateHalos
void CreateHalos(const COSMOLOGY &cosmo, double redshift, double inv_area)
Definition: particle_lens.cpp:153
SimFileFormat
SimFileFormat
Definition: particle_types.h:32
Generated on Sat Mar 16 2024 11:07:14 for GLAMERDOC++ by   doxygen 1.8.18