10 #include "utilities_slsim.h"
16 GLAMER_TEST_USES(LensTest)
73 Lens(
long* seed, PosType z_source,CosmoParamSet cosmoset,
bool verbose =
false);
75 Lens(
long* seed, PosType z_source,
const COSMOLOGY &cosmo,
bool verbose =
false);
87 PosType
getfov()
const {
return fieldofview;};
88 void setfov(PosType fov){fieldofview=fov;};
91 void resetFieldNplanes(std::size_t field_Nplanes,
bool verbose =
false);
95 void resetFieldHalos(
bool verbose =
false);
98 void printMultiLens();
103 if(flag_switch_main_halo_on)
104 return main_halos[0]->getZlens();
107 std::cerr <<
"error, no main lens present" << std::endl;
113 if(flag_switch_main_halo_on)
114 return cosmo.angDist( main_halos[0]->getZlens());
117 std::cerr <<
"error, no main lens present" << std::endl;
125 void clearMainHalos(
bool verbose=
false);
127 template<
typename HaloType>
149 template <
typename T>
153 T * halo =
new T(halo_in);
154 halo->setCosmology(cosmo);
155 main_halos.push_back(halo);
157 flag_switch_main_halo_on =
true;
159 if(addplanes) addMainHaloToPlane(halo);
160 else addMainHaloToNearestPlane(halo);
162 combinePlanes(verbose);
172 template <
typename T>
175 T * halo =
new T(std::move(halo_in));
176 halo->setCosmology(cosmo);
177 main_halos.push_back(halo);
179 flag_switch_main_halo_on =
true;
181 if(addplanes) addMainHaloToPlane(halo);
182 else addMainHaloToNearestPlane(halo);
184 combinePlanes(verbose);
193 template <
typename T>
198 T * halo =
new T(halo_in);
199 halo->setCosmology(cosmo);
200 main_halos.push_back(halo);
202 flag_switch_main_halo_on =
true;
206 combinePlanes(verbose);
219 template <
typename T>
224 for(T &h : my_halos){
226 ptr->setCosmology(cosmo);
228 main_halos.push_back( ptr );
229 if(addplanes) addMainHaloToPlane( ptr );
230 else addMainHaloToNearestPlane( ptr );
233 flag_switch_main_halo_on =
true;
235 combinePlanes(verbose);
243 template <
typename T>
251 for(T &h : my_halos){
253 ptr->setCosmology(cosmo);
255 main_halos.push_back( ptr );
258 flag_switch_main_halo_on =
true;
262 combinePlanes(verbose);
303 std::size_t getNMainHalos()
const;
305 template<
typename HaloType>
306 std::size_t getNMainHalos()
const;
309 LensHalo* getMainHalo(std::size_t i);
312 template<
typename HaloType>
313 HaloType* getMainHalo(std::size_t i);
322 void rayshooter(
RAY &ray);
330 void rayshooterInternal(
unsigned long Npoints
332 ,
bool RSIverbose =
false
340 void rayshooterInternal(
unsigned long Npoints
342 ,std::vector<double> &source_zs
343 ,
bool RSIverbose =
false
346 void rayshooterInternal(
unsigned long Npoints
349 void rayshooterInternal(
RAY &ray);
351 void info_rayshooter(
RAY &i_point
352 ,std::vector<Point_2d> & ang_positions
353 ,std::vector<KappaType> & kappa_on_planes
354 ,std::vector<std::vector<LensHalo*>> & halo_neighbors
356 ,KappaType &kappa_max
357 ,KappaType gamma_max[]
361 ,
bool verbose =
false
365 void mass_on_planes(
const std::vector<RAY> &rays
366 ,std::vector<double> &masses
367 ,
bool verbose =
false
398 ,
bool use_image_guess
404 RAY find_image_min(
const RAY &in_ray
420 ,std::vector<Point_2d> &boundary
424 std::vector<RAY> find_images(
Point_2d y_source
431 std::vector<RAY> find_images(
GridMap &init_grid
437 void find_images_min_parallel(std::vector<RAY> &rays
439 ,std::vector<bool> &success
456 short ResetSourcePlane(PosType z,
bool nearest=
false,
bool verbose =
false);
459 void FindSourcePlane(
469 PosType getSourceZ(){
470 if(toggle_source_plane){
473 return plane_redshifts.back();
477 PosType getZmax()
const{
return plane_redshifts.back();}
483 PosType
getSigmaCrit(PosType zsource)
const{
return cosmo.SigmaCrit(getZlens(), zsource); }
491 void TurnFieldOn() { flag_switch_field_off = false ; }
497 bool getfieldOff()
const {
return flag_switch_field_off ;}
503 void GenerateFieldHalos(
double min_mass
505 ,
double field_of_view
510 ,
bool verbose =
false
517 init_seed = lens.init_seed;
519 toggle_source_plane = lens.toggle_source_plane;
521 zs_implant = lens.zs_implant;
522 zsource = lens.zsource;
524 NZSamples = lens.NZSamples;
526 NhalosbinZ = lens.NhalosbinZ;
527 Nhaloestot_Tab = lens.Nhaloestot_Tab;
528 aveNhalosField = lens.aveNhalosField;
530 NhalosbinMass = lens.NhalosbinMass;
531 sigma_back_Tab = lens.sigma_back_Tab;
533 flag_switch_deflection_off = lens.flag_switch_deflection_off;
534 flag_switch_lensing_off = lens.flag_switch_lensing_off;
539 plane_redshifts = lens.plane_redshifts;
540 charge = lens.charge;
541 flag_switch_field_off = lens.flag_switch_field_off;
543 field_halos = lens.field_halos;
545 field_Nplanes_original = lens.field_Nplanes_original;
546 field_Nplanes_current = lens.field_Nplanes_current;
548 field_plane_redshifts = lens.field_plane_redshifts;
549 field_plane_redshifts_original = lens.field_plane_redshifts_original;
550 field_Dl = lens.field_Dl;
551 field_Dl_original = lens.field_Dl_original;
556 field_mass_func_type = lens.field_mass_func_type;
557 mass_func_PL_slope = lens.mass_func_PL_slope;
558 field_min_mass = lens.field_min_mass;
559 field_int_prof_type = lens.field_int_prof_type;
560 field_prof_internal_slope = lens.field_prof_internal_slope;
562 flag_field_gal_on = lens.flag_field_gal_on;
563 field_int_prof_gal_type = lens.field_int_prof_gal_type;
564 field_int_prof_gal_slope = lens.field_int_prof_gal_slope;
566 redshift_planes_file = lens.redshift_planes_file;
567 read_redshift_planes = lens.read_redshift_planes;
569 field_input_sim_file = lens.field_input_sim_file;
570 field_input_sim_format = lens.field_input_sim_format;
572 sim_input_flag = lens.sim_input_flag;
573 read_sim_file = lens.read_sim_file;
574 field_buffer = lens.field_buffer;
576 flag_switch_main_halo_on = lens.flag_switch_main_halo_on;
578 main_plane_redshifts = lens.main_plane_redshifts;
579 main_Dl = lens.main_Dl;
581 main_halo_type = lens.main_halo_type;
582 main_galaxy_halo_type = lens.main_galaxy_halo_type;
584 pixel_map_input_file = lens.pixel_map_input_file;
585 pixel_map_on = lens.pixel_map_on;
586 pixel_map_zeropad = lens.pixel_map_zeropad;
587 pixel_map_zeromean = lens.pixel_map_zeromean;
589 central_point_sphere = lens.central_point_sphere;
590 sim_angular_radius = lens.sim_angular_radius;
591 inv_ang_screening_scale = lens.inv_ang_screening_scale;
593 std::swap(lensing_planes,lens.lensing_planes);
594 std::swap(field_planes,lens.field_planes);
595 swap(main_halos,lens.main_halos);
596 std::swap(main_planes,lens.main_planes);
601 *
this = std::move(lens);
605 std::vector<PosType> get_plane_redshifts () {
return plane_redshifts; }
616 GLAMER_TEST_FRIEND(LensTest)
619 void _find_images_(std::vector<RAY> &images
625 void _find_images_min_parallel_(
RAY *rays
629 ,std::vector<bool> &success
642 void compute_points_parallel(
int start
647 ,
bool verbose =
false);
650 void compute_rays_parallel(
int start
653 ,
bool verbose =
false);
657 void defaultParams(PosType zsource,
bool verbose =
true);
660 bool toggle_source_plane;
668 void quicksort(
LensHaloHndl *halo,PosType **pos,
unsigned long N);
676 void setFieldDistFromFile();
678 void setupFieldPlanes();
681 void ComputeHalosDistributionVariables ();
683 enum DM_Light_Division {All_DM,Moster};
685 void createFieldHalos(
bool verbose,DM_Light_Division division = Moster);
688 void readInputSimFileMillennium(
bool verbose,DM_Light_Division division = Moster);
691 void readInputSimFileMultiDarkHalos(
bool verbose,DM_Light_Division division = Moster);
694 void readInputSimFileObservedGalaxies(
bool verbose);
697 void createFieldPlanes(
bool verbose);
702 void createMainPlanes();
704 void addMainHaloToPlane(
LensHalo* halo);
705 void addMainHaloToNearestPlane(
LensHalo* halo);
708 void combinePlanes(
bool verbose);
714 const int Nzbins = 64 ;
716 const int Nmassbin=64;
720 std::vector<PosType> zbins ;
722 std::vector<PosType> NhalosbinZ ;
724 std::vector<PosType> Nhaloestot_Tab ;
726 PosType aveNhalosField ;
728 std::vector<PosType> Logm;
730 std::vector<std::vector<PosType>> NhalosbinMass;
732 std::vector<PosType> sigma_back_Tab;
737 size_t getNFieldHalos()
const {
return field_halos.size();}
743 bool flag_switch_deflection_off;
744 bool flag_switch_lensing_off;
747 std::vector<LensPlane *> lensing_planes;
749 std::vector<PosType> Dl;
751 std::vector<PosType> dDl;
753 std::vector<PosType> dTl;
755 std::vector<PosType> plane_redshifts;
761 bool flag_switch_field_off;
764 std::vector<LensHalo *> field_halos;
766 std::size_t field_Nplanes_original;
767 std::size_t field_Nplanes_current;
770 std::vector<LensPlane*> field_planes;
772 std::vector<PosType> field_plane_redshifts;
774 std::vector<PosType> field_plane_redshifts_original;
776 std::vector<PosType> field_Dl;
778 std::vector<PosType> field_Dl_original;
809 PosType mass_func_PL_slope;
811 PosType field_min_mass;
815 PosType field_prof_internal_slope;
818 bool flag_field_gal_on;
822 PosType field_int_prof_gal_slope;
827 std::string redshift_planes_file;
828 bool read_redshift_planes;
830 std::string field_input_sim_file;
831 HaloCatFormats field_input_sim_format;
839 PosType field_buffer;
843 bool flag_switch_main_halo_on;
849 std::vector<LensPlane*> main_planes;
851 std::vector<PosType> main_plane_redshifts;
853 std::vector<PosType> main_Dl;
862 std::string pixel_map_input_file;
867 int pixel_map_zeropad;
868 bool pixel_map_zeromean;
869 void readPixelizedDensity();
874 PosType sim_angular_radius;
876 PosType inv_ang_screening_scale;
879 MINyFunction(
Lens &mylens,
Point_2d y,
int sign):lens(mylens),y(y),sign(sign),r2max(0){}
881 double operator()(
double *x){
884 lens.rayshooterInternal(1,&point);
885 double r2 = (y[0]-point.image->x[0])*(y[0]-point.image->x[0])
886 + (y[1]-point.image->x[1])*(y[1]-point.image->x[1]);
888 r2max = MAX(r2,r2max);
889 return r2 + r2max*abs(sign - sgn(point.invmag()));
902 return main_halos.
size();
905 template<
typename HaloType>
908 return main_halos.
size<HaloType>();
913 return main_halos.
at(i);
916 template<
typename HaloType>
919 if(main_halos.
size<HaloType>() == 0 )
return nullptr;
920 return main_halos.
at<HaloType>(i);
926 template <
typename P>
927 void Lens::compute_points_parallel(
int start
934 int end = start + chunk_size;
942 KappaType kappa,gamma[3];
945 Matrix2x2<PosType> G;
947 PosType SumPrevAlphas[2];
948 Matrix2x2<PosType> SumPrevAG;
952 long jmax = lensing_planes.size();
957 if(source_zs[0] == plane_redshifts.back() ){
958 Dls_Ds = dDl.back() / Dl.back();
959 D_Ds = Dl[Dl.size() - 2] / Dl.back();
964 if(jmax > 0) D_Ds = Dl[jmax-1] / Ds;
969 for(i = start; i < end; i++)
972 if(i_points[i].in_image == MAYBE)
continue;
975 theta = i_points[i].ptr_y();
977 theta[0] = i_points[i].x[0];
978 theta[1] = i_points[i].x[1];
981 SumPrevAlphas[0] = theta[0];
982 SumPrevAlphas[1] = theta[1];
991 i_points[i].A.setToI();
995 if(flag_switch_lensing_off)
997 i_points[i].image->A.setToI();
1010 if(jmax > 0) D_Ds = Dl[jmax-1]/Ds;
1016 for(j = 0; j < jmax ; ++j)
1019 double Dphysical = Dl[j]/(1 + plane_redshifts[j]);
1021 xx[0] = theta[0] * Dphysical;
1022 xx[1] = theta[1] * Dphysical;
1025 assert(xx[0] == xx[0] && xx[1] == xx[1]);
1029 lensing_planes[j]->force(alpha,&kappa,gamma,&phi,xx);
1034 assert(alpha[0] == alpha[0] && alpha[1] == alpha[1]);
1035 assert(gamma[0] == gamma[0] && gamma[1] == gamma[1]);
1036 assert(kappa == kappa);
1037 if(std::isinf(kappa)) { std::cout <<
"xx = " << xx[0] <<
" " << xx[1] << std::endl ;}
1038 assert(!std::isinf(kappa));
1040 G[0] = kappa + gamma[0]; G[1] = gamma[1];
1041 G[2] = gamma[1]; G[3] = kappa - gamma[0];
1047 G *= charge * Dl[j] / (1 + plane_redshifts[j]);
1049 assert(gamma[0] == gamma[0] && gamma[1] == gamma[1]);
1050 assert(kappa == kappa);
1057 SumPrevAlphas[0] -= charge * alpha[0] ;
1058 SumPrevAlphas[1] -= charge * alpha[1] ;
1061 aa = dDl[j+1] / Dl[j+1];
1062 bb = Dl[j] / Dl[j+1];
1068 if(!flag_switch_deflection_off){
1069 theta[0] = bb * theta[0] + aa * SumPrevAlphas[0];
1070 theta[1] = bb * theta[1] + aa * SumPrevAlphas[1];
1076 SumPrevAG -= (G * (i_points[i].A)) ;
1079 i_points[i].A = i_points[i].A * bb + SumPrevAG * aa;
1092 if(flag_switch_deflection_off){
1093 i_points[i].A = Matrix2x2<KappaType>::I() - SumPrevAG;
1100 i_points[i].dt *= MpcToSeconds * SecondToYears ;
1105 i_points[i].image->A = i_points[i].A;
1106 i_points[i].image->dt = i_points[i].dt;
1112 if(verbose) std::cout <<
"RSI final : X X | " << i <<
" " << source_zs[i] <<
" | " << i_points[i].kappa() <<
" " << i_points[i].gamma1() <<
" " << i_points[i].gamma2() <<
" " << i_points[i].gamma3() <<
" " << i_points[i].invmag() <<
" | " << i_points[i].dt << std::endl ;