GLAMERDOC++
Gravitational Lensing Code Library
lens_halos.h
1 /*
2  * lens_halos.h
3  *
4  * Created on: 06.05.2013
5  */
6 
7 #ifndef LENS_HALOS_H_
8 #define LENS_HALOS_H_
9 
10 //#include "standard.h"
11 #include "InputParams.h"
12 //#include "source.h"
13 //#include "point.h"
14 #include "quadTree.h"
15 #include "particle_types.h"
16 #include "image_processing.h"
17 #include "fftw3.h"
18 
19 #include <complex>
20 #include <complex.h>
21 #ifdef ENABLE_CERF
22 #include <cerf.h>
23 #endif
24 #ifdef ENABLE_EIGEN
25 //#include </usr/local/include/eigen3/Eigen/Dense>
26 #include <eigen3/Eigen/Dense>
27 #include <eigen3/Eigen/StdVector>
28 #include <boost/math/special_functions/gamma.hpp>
29 #endif
30 
57 class LensHalo{
58 public:
59  LensHalo();
60  LensHalo(PosType z,const COSMOLOGY &cosmo);
61  //LensHalo(InputParams& params,COSMOLOGY &cosmo,bool needRsize = true);
62  //LensHalo(InputParams& params,bool needRsize = true);
63  LensHalo(const LensHalo &h);
64 
65  LensHalo(LensHalo &&h){
66 // LensHalo(LensHalo &&h):star_tree(nullptr){
67  *this = std::move(h);
68  }
69 
70  virtual ~LensHalo();
71 
72  LensHalo & operator=(const LensHalo &h);
73  LensHalo & operator=(LensHalo &&h);
74 
75  int tag=0;
76 
80  inline float get_Rmax() const { return Rmax; }
82  inline float getRsize() const { return Rsize; }
84  inline float get_mass() const { return mass; }
86  inline float get_rscale() const { return rscale; }
88  inline PosType getZlens() const { return zlens; }
89 
90  // set the position of the Halo in physical Mpc on the lens plane
91  //void setX(PosType PosX, PosType PosY) { posHalo[0] = PosX ; posHalo[1] = PosY ; }
92  // set the position of the Halo in physical Mpc on the lens plane
93  //void setX(PosType *PosXY) { posHalo[0] = PosXY[0] ; posHalo[1] = PosXY[1] ; }
94 
96  void getX(PosType * MyPosHalo) const {
97  assert(Dist != -1);
98  MyPosHalo[0] = posHalo[0]*Dist;
99  MyPosHalo[1] = posHalo[1]*Dist;
100  }
101  Point_2d getX() const { return Point_2d(posHalo[0]*Dist , posHalo[1]*Dist); }
102 
104  PosType operator[](int i) const{return posHalo[i]*Dist;}
105 
107  void setTheta(PosType PosX, PosType PosY) { posHalo[0] = PosX ; posHalo[1] = PosY ; }
109  void setTheta(PosType *PosXY) { posHalo[0] = PosXY[0] ; posHalo[1] = PosXY[1] ; }
111  void setTheta(const Point_2d &p) { posHalo[0] = p[0] ; posHalo[1] = p[1] ; }
113  void getTheta(PosType * MyPosHalo) const { MyPosHalo[0] = posHalo[0] ; MyPosHalo[1] = posHalo[1]; }
114 
115  Point_2d getTheta() const { return Point_2d(posHalo[0] , posHalo[1]); }
116 
118  void setDist(COSMOLOGY &co){Dist = co.angDist(zlens);}
119 
122  PosType getDist() const {return Dist;}
123 
126  void displayPos() { std::cout << "Halo PosX = " << posHalo[0] << " ; Halo PosY = " << posHalo[1] << std::endl; }
127 
129  virtual void initFromFile(float my_mass, long *seed, float vmax, float r_halfmass){};
131  virtual void initFromMassFunc(float my_mass, float my_Rsize, float my_rscale, PosType my_slope, long *seed);
132 
134  virtual void set_RsizeRmax(float my_Rsize){Rmax = Rmax*my_Rsize/Rsize; Rsize = my_Rsize; xmax = Rsize/rscale;};
136  void set_mass(float my_mass){mass=my_mass;};
138  virtual void set_rscale(float my_rscale){rscale=my_rscale; xmax = Rsize/rscale;};
140  //void setZlens(PosType my_zlens){zlens=my_zlens; Dist=-1;}
141  void setZlens(PosType my_zlens,const COSMOLOGY &cosmo){
142  zlens=my_zlens;
143  Dist=cosmo.angDist(my_zlens);
144  }
145  void setRsize(PosType R){Rsize = R;}
146 
147  // set redshift and distance
148  void setZlensDist(PosType my_zlens,const COSMOLOGY &cos){
149  zlens=my_zlens;
150  Dist = cos.angDist(zlens);
151  }
152  void setMass(PosType m){mass = m;}
153 
154 
156  virtual void set_slope(PosType my_slope){beta=my_slope;};
158  virtual PosType get_slope(){return beta;};
160  bool get_flag_elliptical(){return elliptical_flag;};
161  void set_flag_elliptical(bool ell){elliptical_flag=ell;};
162  bool get_switch_flag(){return switch_flag;};
163  void set_switch_flag(bool swt){switch_flag=swt;};
164 
165 
167  virtual void setCosmology(const COSMOLOGY& cosmo) {}
168 
169  /* calculate the lensing properties -- deflection, convergence, and shear
170  if not overwritten by derived class it uses alpha_h(), gamma_h(), etc. of the
171  derived case or for a point source in this class
172 
173  xcm - the physical position on the lens plane relative to the center of the LensHalo in Mpc
174  Units :
175 
176  ALPHA - mass/PhysMpc - ALPHA / Sig_crit / Dl is the deflection in radians
177  KAPPA - surface mass density , mass / /PhysMpc/PhysMpc - KAPPA / Sig_crit is the convergence
178  GAMMA - mass / /PhysMpc/PhysMpc - GAMMA / Sig_crit is the shear
179  PHI - mass - PHI / Sig_crit is the lensing potential
180 
181 * returns the lensing quantities of a ray in center of mass coordinates.
182  *
183  * Warning: This adds to input value of alpha, kappa, gamma, and phi. They need
184  * to be zeroed out if the contribution of just this halo is wanted.
185  */
186  virtual void force_halo(
187  PosType *alpha
188  ,KappaType *kappa
189  ,KappaType *gamma
190  ,KappaType *phi
191  ,PosType const *xcm
192  ,bool subtract_point=false
193  ,PosType screening=1.0
194  );
195 
197  //void force_stars(PosType *alpha,KappaType *kappa,KappaType *gamma,PosType const *xcm);
198 
200  bool compareZ(PosType z){return z > zlens;};
201 
203 // bool AreStarsImplanted() const {return stars_implanted;}
204 // void implant_stars(PosType **centers,int Nregions,long *seed, IMFtype type=One);
205 // void implant_stars(PosType *center,long *seed,IMFtype type = One);
206 // //void implant_stars(PosType *x,PosType *y,int Nregions,long *seed,IMFtype type=One);
207 // void remove_stars();
208 // IMFtype getStarIMF_type() const {return main_stars_imf_type;}
209 // /// Fraction of surface density in stars
210 // PosType getFstars() const {return star_fstars;}
211 // /// The mass of the stars if they are all the same mass
212 // PosType getStarMass() const {if(stars_implanted)return stars_xp[0].mass(); else return 0.0;}
213 
215  EllipMethod getEllipMethod() const {return main_ellip_method;}
217  std::vector<double> get_mod() { std::vector<double> fmodes(Nmod); for(int i=0;i<Nmod;i++){fmodes[i]= mod[i] ;} ;return fmodes;}
219  const static int get_Nmod() {return Nmod;}
220 
222  virtual std::size_t Nparams() const;
224  virtual PosType getParam(std::size_t p) const;
226  virtual PosType setParam(std::size_t p, PosType value);
227 
229  virtual void printCSV(std::ostream&, bool header = false) const;
230 
232  //void PrintStars(bool show_stars);
233 
234  PosType MassBy2DIntegation(PosType R);
235  PosType MassBy1DIntegation(PosType R);
236  PosType test_average_gt(PosType R);
237  PosType test_average_kappa(PosType R);
238 
239  // renomalize to make mass match
240 
241  void set_norm_factor(){mass_norm_factor=1;mass_norm_factor=mass/MassBy1DIntegation(Rsize);}
242 
244  void set_rsize(float my_rsize){ Rsize = my_rsize;};
245  float get_rsize(){return Rsize;};
246 
247  // all of the following functions were used for Ansatz III w derivatives of the Fourier modes
248 
250  bool test();
251 
252  size_t getID() const {return idnumber;}
253  void setID(size_t id){idnumber = id;}
254 
255  PosType renormalization(PosType r_max);
256 
270  PixelMap<double> map_variables(
271  LensingVariable lensvar
272  ,size_t Nx
273  ,size_t Ny
274  ,double res
275  );
276 
277 private:
278  size_t idnumber;
279  PosType posHalo[2];
281  PosType zlens;
282 
283 protected:
284  // This is the size of the halo beyond which it does not have the expected profile.
285  float Rsize = 0;
286 
287  // total mass in Msun
288  float mass;
289  // angular size distance to this lens
290  PosType Dist;
291  PosType mnorm;
292 
293  // Beyond Rmax the halo will be treated as a point mass. Between Rsize and Rmax the deflection
294  // and shear are interpolated. For circularly symmetric lenses Rmax should be equal to Rsize
295  float Rmax;
296 
297 
298  PosType alpha_int(PosType x) const;
299  PosType norm_int(PosType r_max);
300 
301  void force_halo_sym(PosType *alpha,KappaType *kappa,KappaType *gamma,KappaType *phi,PosType const *xcm,bool subtract_point=false,PosType screening = 1.0);
302  void force_halo_asym(PosType *alpha,KappaType *kappa,KappaType *gamma,KappaType *phi,PosType const *xcm,bool subtract_point=false,PosType screening = 1.0);
303 
304  bool force_point(PosType *alpha,KappaType *kappa,KappaType *gamma
305  ,KappaType *phi,PosType const *xcm,PosType rcm2
306  ,bool subtract_point,PosType screening);
307 
308  struct norm_func{
309  norm_func(LensHalo& halo, PosType my_r_max): halo(halo), r_max(my_r_max){};
310  LensHalo& halo;
311  PosType r_max;
312  //PosType operator ()(PosType theta) {halo.alpha_asym(r_max, theta, alpha_arr); return alpha_arr[0]*cos(theta)*cos(theta)+alpha_arr[1]*sin(theta)*sin(theta);}
313  PosType operator ()(PosType theta) {return halo.alpha_ell(r_max, theta);}
314  };
315 
316  struct Ialpha_func{
317  Ialpha_func(LensHalo& halo): halo(halo){};
318  LensHalo& halo;
319  PosType operator ()(PosType x) {return halo.alpha_h(x)/x ;}
320  };
321 
322  struct Ig_func{
323  Ig_func(const LensHalo& halo): halo(halo){};
324  const LensHalo& halo;
325  PosType operator ()(PosType x) {return halo.gfunction(x)/x ;}
326  };
327 
328 // std::vector<IndexType> stars_index;
329 // std::vector<StarType> stars_xp;
330 // TreeQuadParticles<StarType> *star_tree;
331 //
332 // int stars_N;
333 // PosType star_massscale;
334 // /// star masses relative to star_massscles
335 // //float *star_masses;
336 // PosType star_fstars;
337 // PosType star_theta_force;
338 // int star_Nregions;
339 // std::vector<PosType> star_region;
340  PosType beta;
341 // void substract_stars_disks(PosType const *ray,PosType *alpha
342 // ,KappaType *kappa,KappaType *gamma);
343 // std::vector<float> stellar_mass_function(IMFtype type, unsigned long Nstars, long *seed, PosType minmass=0.0, PosType maxmass=0.0
344 // ,PosType bendmass=0.0, PosType powerlo=0.0, PosType powerhi=0.0);
345 
346 
348  void assignParams(InputParams& params,bool needRsize);
350  //void assignParams_stars(InputParams& params);
351 
353  void error_message1(std::string name,std::string filename);
354 
355 
357  float Rmax_to_Rsize_ratio = 1.2;
358 
360  float rscale;
361  // redshift
362  //PosType zlens;
363 
364  EllipMethod main_ellip_method;
365 
366 // bool stars_implanted;
367 // /// Number of regions to be subtracted to compensate for the mass in stars
368 // IMFtype main_stars_imf_type;
369 // PosType main_stars_min_mass;
370 // PosType main_stars_max_mass;
371 // PosType bend_mstar;
372 // PosType lo_mass_slope;
373 // PosType hi_mass_slope;
374 // /// parameters for stellar mass function: minimal and maximal stellar mass, bending point for a broken power law IMF
375 // std::vector<PosType> star_Sigma;
376 // std::vector<Point_2d> star_xdisk;
377 //
378 
381  virtual PosType inline alpha_h(PosType x) const {return -1;};
382  virtual KappaType inline kappa_h(PosType x) const {return 0;};
383  virtual KappaType inline gamma_h(PosType x) const {return -2;};
384  virtual KappaType inline phi_h(PosType x) const {return 1;};
385  virtual KappaType inline phi_int(PosType x) const {return 1;};
386  virtual PosType inline ffunction(PosType x)const {return 0;};
387  virtual PosType inline gfunction(PosType x) const {return -1;};
388  virtual PosType inline dgfunctiondx(PosType x){return 0;};
389  virtual PosType inline bfunction(PosType x){return -1;};
390  virtual PosType inline dhfunction(PosType x) const {return 1;};
391  virtual PosType inline ddhfunction(PosType x, bool numerical){return 0;};
392  virtual PosType inline dddhfunction(PosType x, bool numerical){return 0;};
393  virtual PosType inline bnumfunction(PosType x){return -1;};
394  virtual PosType inline dbfunction(PosType x){return 0;};
395  virtual PosType inline ddbfunction(PosType x){return 0;};
396  virtual PosType inline dmoddb(int whichmod, PosType q, PosType b){return 0;};
397  virtual PosType inline ddmoddb(int whichmod, PosType q, PosType b){return 0;};
398  virtual PosType inline dmoddq(int whichmod, PosType q, PosType b){return 0;};
399  virtual PosType inline ddmoddq(int whichmod, PosType q, PosType b){return 0;};
400 
401  PosType xmax;
402 
403  PosType mass_norm_factor=1;
404 
405  // Functions for calculating axial dependence
406  float pa;
407  float fratio=1.0;
408  bool elliptical_flag = false;
409  bool switch_flag = false;
410 
411 
412  void faxial(PosType x,PosType theta,PosType f[]);
413  void faxial0(PosType theta,PosType f0[]);
414  void faxial1(PosType theta,PosType f1[]);
415  void faxial2(PosType theta,PosType f2[]);
416  void gradial(PosType r,PosType g[]);
417  void gradial2(PosType r,PosType mu, PosType sigma,PosType g[]);
418 
419  void felliptical(PosType x, PosType q, PosType theta, PosType f[], PosType g[]);
420 
421  virtual void gamma_asym(PosType x,PosType theta, PosType gamma[]);
422  virtual PosType kappa_asym(PosType x,PosType theta);
423  virtual void alphakappagamma_asym(PosType x,PosType theta, PosType alpha[]
424  ,PosType *kappa,PosType gamma[],PosType *phi);
425  virtual void alphakappagamma1asym(PosType x,PosType theta, PosType alpha[2]
426  ,PosType *kappa,PosType gamma[],PosType *phi);
427  virtual void alphakappagamma2asym(PosType x,PosType theta, PosType alpha[2]
428  ,PosType *kappa,PosType gamma[],PosType *phi);
429  virtual void alphakappagamma3asym(PosType x,PosType theta, PosType alpha[2]
430  ,PosType *kappa,PosType gamma[],PosType *phi);
431 
432  virtual PosType alpha_ell(PosType x,PosType theta);
433 
434  double fourier_coeff(double n, double q, double beta);
435  double IDAXDM(double lambda, double a2, double b2, double x[], double rmax, double mo);
436  double IDAYDM(double lambda, double a2, double b2, double x[], double rmax, double mo);
437  double SCHRAMMKN(double n, double x[], double rmax);
438  double SCHRAMMJN(double n, double x[], double rmax);
439  double SCHRAMMI(double x[], double rmax);
440 
441 
442  void calcModes(double q, double beta, double rottheta, PosType newmod[]);
443  void calcModesB(PosType x, double q, double beta, double rottheta, PosType newmod[]);
444  void calcModesC(PosType beta_r, double q, double rottheta, PosType newmod[]);
445 
446  virtual PosType inline InterpolateModes(int whichmod, PosType q, PosType b){return 0;};
447 
448  void analModes(int modnumber, PosType my_beta, PosType q, PosType amod[3]);
449 
450  struct fourier_func{
451  fourier_func(double my_n, double my_q, double my_beta): n(my_n),q(my_q),beta(my_beta){};
452  double n;
453  double q;
454  double beta;
455  double operator ()(double theta) {return cos(n*theta)/pow(cos(theta)*cos(theta) + 1/q/q*sin(theta)*sin(theta),beta/2) ;}
456  };
457 
458  struct SCHRAMMJ_func{
459  SCHRAMMJ_func(double my_n, double my_x[], double my_rmax, LensHalo *my_halo): n(my_n), x(my_x), rmx(my_rmax), halo(my_halo){};
460  double n, *x, rmx;
461  double operator ()(double u) {
462 
463  // PosType xisq=sqrt(u*(x[0]*x[0]+x[1]*x[1]/(1-(1-halo->fratio*halo->fratio)*u)));
464  PosType xisq=sqrt((x[0]*x[0]+x[1]*x[1]/(1-u*(1-halo->fratio*halo->fratio))));
465 
466  if(xisq*xisq < 1e-20) xisq = 1.0e-10;
467  KappaType kappa=halo->kappa_h(xisq)/PI/xisq/xisq;
468  //std::cout << "Schramm: n=" << n << " " << u << " " << xisq << " " << halo->fratio << " " << rmx << " " << halo->Rmax << " " << halo->rscale << std::endl;
469  return kappa*halo->mass/pow((1.-(1.-halo->fratio*halo->fratio)*u),n+0.5);
470  }
471  private:
472  LensHalo *halo;
473  };
474 
475  struct SCHRAMMK_func{
476  SCHRAMMK_func(double my_n, double my_x[], double my_rmax, LensHalo *my_halo): n(my_n), x(my_x), rmx(my_rmax), halo(my_halo){};
477  double n, *x, rmx;
478  double operator ()(double u) {
479 
480  // PosType xisq=sqrt(u*(x[0]*x[0]+x[1]*x[1]/(1-(1-halo->fratio*halo->fratio)*u)));
481  PosType xisq=sqrt((x[0]*x[0]+x[1]*x[1]/(1-u*(1-halo->fratio*halo->fratio))));
482  if(xisq*xisq < 1e-20) xisq = 1.0e-10;
483  PosType h=0.0001;
484  PosType kp1=halo->kappa_h(xisq+h)/((xisq+h)*(xisq+h));
485  PosType km1=halo->kappa_h(xisq-h)/((xisq-h)*(xisq-h));
486  PosType kp2=halo->kappa_h(xisq+2*h)/((xisq+2*h)*(xisq+2*h));
487  PosType km2=halo->kappa_h(xisq-2*h)/((xisq-2*h)*(xisq-2*h));
488  PosType dkdxi=(8*kp1-8*km1-kp2+km2)/12/h/PI;
489 
490 
491 
492  //std::cout << "Schramm: n=" << n << " " << u << " " << xisq << " " << halo->fratio << " " << rmx << " " << halo->Rmax << " " << halo->rscale << std::endl;
493  return u*dkdxi/pow((1.-(1.-halo->fratio*halo->fratio)*u),n+0.5);
494  }
495  private:
496  LensHalo *halo;
497  };
498 
499 
500  struct SCHRAMMI_func{
501  SCHRAMMI_func(double my_x[], double my_rmax, LensHalo *my_halo): x(my_x), rmx(my_rmax), halo(my_halo){};
502  double *x, rmx;
503  double operator ()(double u) {
504 
505  // PosType xisq=sqrt(u*(x[0]*x[0]+x[1]*x[1]/(1-(1-halo->fratio*halo->fratio)*u)));
506  PosType xisq=sqrt((x[0]*x[0]+x[1]*x[1]/(1-u*(1-halo->fratio*halo->fratio))));
507 
508  if(xisq*xisq < 1e-20) xisq = 1.0e-10;
509  PosType alpha=halo->alpha_h(xisq)/PI/xisq;
510  //std::cout << "Schramm: n=" << n << " " << u << " " << xisq << " " << halo->fratio << " " << rmx << " " << halo->Rmax << " " << halo->rscale << std::endl;
511  return xisq*alpha*halo->mass/pow((1.-(1.-halo->fratio*halo->fratio)*u),0.5);
512  }
513  private:
514  LensHalo *halo;
515  };
516 
517 
518  struct IDAXDM_func{
519  IDAXDM_func(double my_lambda,double my_a2, double my_b2, double my_x[], double my_rmax, LensHalo *my_halo): lambda(my_lambda), a2(my_a2), b2(my_b2), x(my_x), rmx(my_rmax), halo(my_halo){};
520  double lambda, a2,b2, *x, rmx;
521  double operator ()(double m) {
522  double ap = m*m*a2 + lambda,bp = m*m*b2 + lambda;
523  double p2 = x[0]*x[0]/ap/ap/ap/ap + x[1]*x[1]/bp/bp/bp/bp; // actually the inverse of equation (5) in Schramm 1990
524  PosType tmp = m*rmx;
525  if(tmp*tmp < 1e-20) tmp = 1.0e-10;
526  PosType xiso=tmp/halo->rscale;
527  KappaType kappa=halo->kappa_h(xiso)/PI/xiso/xiso;
528 
529  PosType alpha[2]={0,0},tm[2] = {m*rmx,0};
530  KappaType kappam=0,gamma[2]={0,0},phi=0;
531  halo->force_halo_sym(alpha,&kappam,gamma,&phi,tm);
532 
533 
534  //std::cout << "output x: " << m << " " << xiso << " " << m*kappa/(ap*ap*ap*bp*p2)*halo->mass << " " << m*kappam/(ap*ap*ap*bp*p2)<< std::endl;
535  double integrandA=m*kappa/(ap*ap*ap*bp*p2)*halo->mass;
536  double integrandB=m*kappam/(ap*ap*ap*bp*p2);
537  //std::cout << integrandA-integrandB << std::endl;
538  assert( std::abs((integrandA - integrandB)/integrandA)<1E-5);
539 
540  assert(kappa >= 0.0);
541 
542  //return m*kappa/(ap*ap*ap*bp*p2)*halo->mass; // integrand of equation (28) in Schramm 1990
543  return integrandB;
544  }
545  private:
546  LensHalo *halo;
547  };
548 
549  struct IDAYDM_func{
550  IDAYDM_func(double my_lambda,double my_a2, double my_b2, double my_x[], double my_rmax, LensHalo *my_halo): lambda(my_lambda), a2(my_a2), b2(my_b2), x(my_x), rmx(my_rmax), halo(my_halo){};
551  double lambda, a2,b2, *x, rmx;
552  double operator ()(double m) {
553  double ap = m*m*a2 + lambda,bp = m*m*b2 + lambda;
554  double p2 = x[0]*x[0]/ap/ap/ap/ap + x[1]*x[1]/bp/bp/bp/bp; // actually the inverse of equation (5) in Schramm 1990
555  PosType tmp = m*rmx;
556  if(tmp*tmp < 1e-20) tmp = 1.0e-10;
557  PosType xiso=tmp/halo->rscale;
558  KappaType kappa=halo->kappa_h(xiso)/PI/xiso/xiso;
559 
560  PosType alpha[2]={0,0},tm[2] = {m*rmx,0};
561  KappaType kappam=0,gamma[2]={0,0},phi=0;
562  halo->force_halo_sym(alpha,&kappam,gamma,&phi,tm);
563 
564  double integrandA=m*kappa/(ap*ap*ap*bp*p2)*halo->mass;
565  double integrandB=m*kappam/(ap*bp*bp*bp*p2);
566  assert( std::abs((integrandA - integrandB)/integrandA)<1E-5);
567 
568  assert(kappa >= 0.0);
569  //return m*kappa/(ap*bp*bp*bp*p2)*halo->mass; // integrand of equation (28) in Schramm 1990
570  return integrandB;
571  }
572  private:
573  LensHalo *halo;
574  };
575 
576 
577  const static int Nmod = 32;
578 
579  // Analytic description of Fourier modes
580 
581  PosType mod[Nmod];
582  PosType mod1[Nmod];
583  PosType mod2[Nmod];
584  PosType r_eps;
585 
586 
587  // These are stucts used in doing tests
588 
589  /*struct test_gt_func{
590  test_gt_func(LensHalo& halo,PosType my_r): halo(halo),r(my_r){};
591  LensHalo& halo;
592  PosType r;
593  PosType a[2] = {0,0},x[2] = {0,0};
594  KappaType k = 0,g[3] = {0,0,0} ,p=0;
595  //double operator ()(PosType t) {x[0]=r*cos(t); x[1]=r*sin(t); halo.force_halo(a,&k,g,&p,x); if(r>1){std::cout << r << " " << g[0]*cos(2*t)+g[1]*sin(2*t) << std::endl;} return (g[0]*cos(2*t)+g[1]*sin(2*t));}
596  double operator ()(PosType t) {x[0]=r*cos(t); x[1]=r*sin(t); halo.force_halo(a,&k,g,&p,x); return (g[0]*cos(2*t)+g[1]*sin(2*t));}
597 
598  };*/
599 
600 
601  struct test_gt_func{
602  test_gt_func(PosType my_r,LensHalo *halo): r(my_r),halo(halo){};
603  double operator ()(PosType t) {
604  PosType alpha[2] = {0,0},x[2] = {0,0};
605  KappaType kappa = 0,gamma[3] = {0,0,0},phi =0 ;
606  x[0]=r*cos(t);
607  x[1]=r*sin(t);
608  halo->force_halo(alpha,&kappa,gamma,&phi,x);
609  assert(gamma[0]==gamma[0]);
610  assert(gamma[1]==gamma[1]);
611  return (gamma[0]*cos(2*t)+gamma[1]*sin(2*t));}
612  //double operator ()(PosType t) {x[0]=r*cos(t); x[1]=r*sin(t); halo.force_halo(a,&k,g,&p,x); if(r>1){std::cout << r << " " << g[0]*cos(2*t)+g[1]*sin(2*t) << std::endl;} return (g[0]*cos(2*t)+g[1]*sin(2*t));}
613  private:
614  PosType r;
615  LensHalo *halo;
616  };
617 
618 
619  /* struct test_kappa_func{
620  test_kappa_func(LensHalo& halo,PosType my_r): halo(halo),r(my_r){};
621  LensHalo& halo;
622  PosType r;
623  PosType a[2] = {0,0},x[2] = {0,0};
624  KappaType k = 0,g[3] = {0,0,0} ,p=0;
625  double operator ()(PosType t) {x[0]=r*cos(t); x[1]=r*sin(t); halo.force_halo(a,&k,g,&p,x);return 2*PI*k*r*r; }
626  };
627 */
628 
629  struct test_kappa_func{
630  test_kappa_func(PosType my_r,LensHalo *halo): r(my_r),halo(halo){};
631  double operator ()(PosType t) {
632  PosType alpha[2] = {0,0},x[2] = {0,0};
633  KappaType kappa = 0,gamma[3] = {0,0,0},phi =0 ;
634  x[0]=r*cos(t);
635  x[1]=r*sin(t);
636  halo->force_halo(alpha,&kappa,gamma,&phi,x);
637  assert(kappa==kappa);
638  return kappa; }
639  private:
640  PosType r;
641  LensHalo *halo;
642  };
643 
644  struct DMDTHETA{
645  DMDTHETA(PosType R,LensHalo *halo): R(R),halo(halo){};
646  PosType operator()(PosType theta){
647  PosType alpha[2] = {0,0},x[2] = {0,0};
648  KappaType kappa = 0,gamma[3] = {0,0,0},phi =0 ;
649 
650  x[0] = R*cos(theta);
651  x[1] = R*sin(theta);
652 
653  halo->force_halo(alpha,&kappa,gamma,&phi,x);
654 
655  PosType alpha_r = -alpha[0]*cos(theta) - alpha[1]*sin(theta);
656  //std::cout << alpha[0] << " " << alpha[1] << std::endl;
657  assert( alpha_r == alpha_r );
658  //std::cout << theta << " " << alpha_r << std::endl;
659  return alpha_r;
660  }
661  private:
662  PosType R;
663  LensHalo *halo;
664  };
665 
666  struct DMDRDTHETA{
667  DMDRDTHETA(PosType R,LensHalo *halo): R(R),halo(halo){};
668  PosType operator()(PosType theta){
669 
670  PosType alpha[2] = {0,0},x[2] = {0,0};
671  KappaType kappa = 0,gamma[3] = {0,0,0} ,phi=0;
672 
673  x[0] = R*cos(theta);
674  x[1] = R*sin(theta);
675 
676  halo->force_halo(alpha,&kappa,gamma,&phi,x);
677  assert(kappa == kappa);
678  assert(kappa != INFINITY);
679  return kappa;
680  }
681  protected:
682  PosType R;
683  LensHalo *halo;
684  };
685 
686  struct DMDR{
687  DMDR(LensHalo *halo): halo(halo){};
688  PosType operator()(PosType logR){
689  //if(halo->get_flag_elliptical()){
690  LensHalo::DMDRDTHETA dmdrdtheta(exp(logR),halo);
691  //std::cout << " R = " << exp(logR) << std::endl;
692 
693  if(exp(2*logR) == 0.0) return 0.0;
694  return Utilities::nintegrate<LensHalo::DMDRDTHETA,PosType>(dmdrdtheta,0,2*PI,1.0e-7)
695  *exp(2*logR);
696 // }else{
697 // PosType alpha[2] = {0,0},x[2] = {0,0};
698 // KappaType kappa = 0,gamma[3] = {0,0,0} ,phi=0;
699 //
700 // x[0] = exp(logR);
701 // x[1] = 0;
702 //
703 // halo->force_halo(alpha,&kappa,gamma,&phi,x);
704 // return 2*PI*kappa*exp(2*logR);
705 // }
706  }
707  protected:
708  LensHalo *halo;
709 
710  };
711 
712  struct DALPHAXDM{
713  DALPHAXDM(PosType lambda,PosType a2,PosType b2,PosType X[],LensHalo* myisohalo)
714  :lambda(lambda),a2(a2),b2(b2),x(X),isohalo(myisohalo){};
715 
716  PosType operator()(PosType m);
717  private:
718  PosType lambda;
719  PosType a2;
720  PosType b2;
721  PosType *x;
722  LensHalo* isohalo;
723  };
724  struct DALPHAYDM{
725  DALPHAYDM(PosType lambda,PosType a2,PosType b2,PosType X[],LensHalo* isohalo)
726  :lambda(lambda),a2(a2),b2(b2),x(X),isohalo(isohalo){};
727 
728  PosType operator()(PosType m);
729  private:
730  PosType lambda;
731  PosType a2;
732  PosType b2;
733  PosType *x;
734  LensHalo* isohalo;
735  };
736 };
737 
749 class LensHaloNFW: public LensHalo{
750 public:
752  LensHaloNFW();
753  LensHaloNFW(float my_mass
754  ,float my_Rsize
755  ,PosType my_zlens
756  ,float my_concentration
757  ,float my_fratio
758  ,float my_pa
759  ,const COSMOLOGY &cosmo
760  ,EllipMethod my_ellip_method=EllipMethod::Pseudo
761  );
762 
763  LensHaloNFW(const LensHaloNFW &h):LensHalo(h){
764  ++LensHaloNFW::count;
765  gmax = h.gmax;
766  }
767  LensHaloNFW(const LensHaloNFW &&h):LensHalo(std::move(h)){
768  ++LensHaloNFW::count;
769  gmax = h.gmax;
770  }
771 
772  LensHaloNFW & operator=(const LensHaloNFW &h){
773  if(this != &h){
774  LensHalo::operator=(h);
775  gmax = h.gmax;
776  }
777  return *this;
778  }
779  LensHaloNFW & operator=(const LensHaloNFW &&h){
780  if(this != &h){
781  LensHalo::operator=(std::move(h));
782  gmax = h.gmax;
783  }
784  return *this;
785  }
786 
787  virtual ~LensHaloNFW();
788 
789  PosType ffunction(PosType x) const;
790  PosType gfunction(PosType x) const;
791  PosType dgfunctiondx(PosType x);
792  PosType g2function(PosType x) const;
793  PosType hfunction(PosType x) const;
794  PosType dhfunction(PosType x) const;
795  PosType ddhfunction(PosType x, bool numerical);
796  PosType dddhfunction(PosType x, bool numerical);
797  PosType bfunction(PosType x);
798  PosType dbfunction(PosType x);
799  PosType ddbfunction(PosType x);
800  PosType dmoddb(int whichmod, PosType q, PosType b);
801  PosType ddmoddb(int whichmod, PosType q, PosType b);
802  PosType dmoddq(int whichmod, PosType q, PosType b);
803  PosType ddmoddq(int whichmod, PosType q, PosType b);
804 
805  //PosType dmod(PosType x, int modnumber, PosType my_slope, PosType my_fratio); // was used for Ansatz III w derivatives of the Fourier modes
806  //PosType ddmod(PosType x, int modnumber, PosType my_slope, PosType my_fratio); // was used for Ansatz III w derivatives of the Fourier modes
807 
808 
809  // TODO: BEN: the below functions alphaNFW, kappaNFW and gammaNFW are obsolete and better to be deleted to avoid confusion
810  void alphaNFW(PosType *alpha,PosType *x,PosType Rtrunc,PosType mass,PosType r_scale
811  ,PosType *center,PosType Sigma_crit);
812  KappaType kappaNFW(PosType *x,PosType Rtrunc,PosType mass,PosType r_scale
813  ,PosType *center,PosType Sigma_crit);
814  void gammaNFW(KappaType *gamma,PosType *x,PosType Rtrunc,PosType mass,PosType r_scale
815  ,PosType *center,PosType Sigma_crit);
816 
817  void initFromFile(float my_mass, long *seed, float vmax, float r_halfmass);
818  void initFromMassFunc(float my_mass, float my_Rsize, float my_rscale, PosType my_slope, long *seed);
820  void set_RsizeXmax(float my_Rsize){LensHalo::setRsize(my_Rsize); xmax = LensHalo::getRsize()/rscale; gmax = InterpolateFromTable(gtable,xmax);};
823  void set_rscaleXmax(float my_rscale){rscale=my_rscale; xmax = LensHalo::getRsize()/rscale; gmax = InterpolateFromTable(gtable,xmax);};
824 
825 protected:
827  static const long NTABLE;
829  static const PosType maxrm;
831  static int count;
832 
834  static PosType *ftable,*gtable,*g2table,*htable,*xtable,*xgtable,***modtable; // modtable was used for Ansatz IV and worked well
836  void make_tables();
838  PosType InterpolateFromTable(PosType *table, PosType y) const;
839  PosType InterpolateModes(int whichmod, PosType q, PosType b);
840 
842  void assignParams(InputParams& params);
843 
844  // Override internal structure of halos
846  inline PosType alpha_h(PosType x) const {
847  //return -1.0*InterpolateFromTable(gtable,x)/InterpolateFromTable(gtable,xmax);
848  return -1.0*InterpolateFromTable(gtable,x)/gmax;
849  // return -0.5/x*InterpolateFromTable(gtable,x)/gmax;
850  }
851  inline KappaType kappa_h(PosType x) const{
852  return 0.5*x*x*InterpolateFromTable(ftable,x)/gmax;
853  }
854  inline KappaType gamma_h(PosType x) const{
855  //return -0.25*x*x*InterpolateFromTable(g2table,x)/gmax;
856  return -0.5*x*x*InterpolateFromTable(g2table,x)/gmax;
857  }
858  inline KappaType phi_h(PosType x) const{
859  return 0.0;
860  return 0.25*(InterpolateFromTable(htable,x) - InterpolateFromTable(htable,LensHalo::getRsize()/rscale))/gmax + log(LensHalo::getRsize()) ;
861  // The constant contribution is made to match with the point mass at x = Rsize/rscale.
862  }
863  inline KappaType phi_int(PosType x) const{
864  return -1.0*InterpolateFromTable(xgtable,x)/gmax; //alpha_int(x);
865  }
866 
867 private:
868  PosType gmax;
869 };
870 // ********************
871 
872 
885 class LensHaloPseudoNFW: public LensHalo{
886 public:
888  LensHaloPseudoNFW();
889  LensHaloPseudoNFW(float my_mass,float my_Rsize,PosType my_zlens,float my_rscale,PosType my_beta,float my_fratio,float my_pa,const COSMOLOGY &cosmo, EllipMethod my_ellip_method=EllipMethod::Pseudo);
890  //LensHaloPseudoNFW(InputParams& params);
891  ~LensHaloPseudoNFW();
892 
893  PosType mhat(PosType y, PosType beta) const;
894  PosType gfunction(PosType x) const;
895 
897  void set_slope(PosType my_slope){beta=my_slope; make_tables();};
899  PosType get_slope(){return beta;};
901  //void set_Rsize(float my_Rsize){Rsize = my_Rsize; xmax = Rsize/rscale;};
902 
903  void initFromMassFunc(float my_mass, float my_Rsize, float my_rscale, PosType my_slope, long *seed);
904 
905 private:
907  static const long NTABLE;
909  static const PosType maxrm;
911  static int count;
912 
914  static PosType *mhattable,*xtable;
916  void make_tables();
918  PosType InterpolateFromTable(PosType y) const;
919 
921  void assignParams(InputParams& params);
922 
924  PosType beta;
925 
926  // Override internal structure of halos
928  inline PosType alpha_h(PosType x) const {
929  return -1.0*InterpolateFromTable(x)/InterpolateFromTable(xmax);
930  }
931  inline KappaType kappa_h(PosType x) const {
932  return 0.5*x*x/InterpolateFromTable(xmax)/pow(1+x,beta);
933  }
934  inline KappaType gamma_h(PosType x) const {
935  //return (0.5*x*x/pow(1+x,beta) - InterpolateFromTable(x))/InterpolateFromTable(xmax);
936  return (x*x/pow(1+x,beta) - 2*InterpolateFromTable(x))/InterpolateFromTable(xmax);
937  }
938  inline KappaType phi_h(PosType x) const {
939  return -1.0*alpha_int(x)/InterpolateFromTable(xmax) ;
940  //ERROR_MESSAGE();
941  //std::cout << "time delay has not been fixed for PseudoNFW profile yet." << std::endl;
942  //exit(1);
943  //return 0.0;
944  }
945  inline KappaType phi_int(PosType x) const{
946  return -1.0*alpha_int(x)/InterpolateFromTable(xmax) ;
947  }
948 };
949 
950 
957 class LensHaloPowerLaw: public LensHalo{
958 public:
959  LensHaloPowerLaw();
960  LensHaloPowerLaw(float my_mass,float my_Rsize,PosType my_zlens,PosType my_beta
961  ,float my_fratio,float my_pa,const COSMOLOGY &cosmo
962  ,EllipMethod my_ellip_method=EllipMethod::Pseudo);
963  //LensHaloPowerLaw(InputParams& params);
964  ~LensHaloPowerLaw();
965 
967  //void set_slope(PosType my_slope){beta=my_slope;};
968 
970  //PosType get_slope(){return beta;};
971 
973  void initFromMassFunc(float my_mass, float my_Rsize, float my_rscale, PosType my_slope, long *seed);
974 
975 private:
977  void assignParams(InputParams& params);
978 
980  // PosType beta;
981 
982  // Override internal structure of halos
984  inline PosType alpha_h(
985  PosType x
986  ) const{
987  if(x<xmax*1.0e-6) x=1e-6*xmax;
988  return -1.0*pow(x/xmax,-beta+2);
989  }
991  inline KappaType kappa_h(
992  PosType x
993  ) const {
994  if(x==0) x=1e-6*xmax;
995  double tmp = x/xmax;
996  return 0.5*(-beta+2)*pow(tmp,2-beta);
997  }
999  inline KappaType gamma_h(PosType x) const {
1000  if(x==0) x=1e-6*xmax;
1001  return -beta*pow(x/xmax,-beta+2);
1002  }
1004  inline KappaType phi_h(PosType x) const {
1005  if(x==0) x=1e-6*xmax;
1006  return ( pow(x/xmax,2-beta) - 1 )/(2-beta) + log(LensHalo::getRsize()) ;
1007  }
1008  inline KappaType phi_int(PosType x) const{
1009  //return alpha_int(x);
1010  return -1.0*pow(x/xmax,-beta+2)/(2-beta);
1011  }
1012 };
1013 
1014 
1028 class LensHaloRealNSIE : public LensHalo{
1029 public:
1030 
1032  LensHaloRealNSIE(float my_mass,PosType my_zlens,float my_sigma
1033  ,float my_rcore,float my_fratio,float my_pa,const COSMOLOGY &cosmo);
1034 
1035  // Warning: If my_rcore > 0.0 and my_fratio < 1 then the mass will be somewhat less than my_mass.
1036  //LensHaloRealNSIE(InputParams& params);
1037 
1038  LensHaloRealNSIE(const LensHaloRealNSIE &h):
1039  LensHalo(h)
1040  {
1041  ++objectCount;
1042  sigma = h.sigma;
1043  fratio = h.fratio;
1044  pa = h.pa;
1045  rcore = h.rcore;
1046  units = h.units;
1047  }
1048 
1049  LensHaloRealNSIE &operator=(const LensHaloRealNSIE &h){
1050  if(&h == this) return *this;
1051  LensHalo::operator=(h);
1052  sigma = h.sigma;
1053  fratio = h.fratio;
1054  pa = h.pa;
1055  rcore = h.rcore;
1056  units = h.units;
1057  return *this;
1058  }
1059 
1060  ~LensHaloRealNSIE();
1061 
1063  void force_halo(PosType *alpha,KappaType *kappa,KappaType *gamma,KappaType *phi,PosType const *xcm,bool subtract_point=false,PosType screening = 1.0);
1064 
1066  float get_sigma(){return sigma;};
1067  // get the NSIE radius
1068  //float get_Rsize(){return Rsize;};
1070  float get_fratio(){return fratio;};
1072  float get_pa(){return pa;};
1074  float get_rcore(){return rcore;};
1075 
1077  void set_sigma(float my_sigma){sigma=my_sigma;};
1078  // set the NSIE radius
1079  //void set_Rsize(float my_Rsize){Rsize=my_Rsize;};
1081  void set_fratio(float my_fratio){fratio=my_fratio;};
1083  void set_pa(float my_pa){pa=my_pa;};
1085  void set_rcore(float my_rcore){rcore=my_rcore;};
1086 
1087  void setZlens(PosType my_zlens,const COSMOLOGY &cosmo){
1088  LensHalo::setZlens(my_zlens,cosmo);
1089  }
1090 
1091 
1092 protected:
1093 
1094  float units;
1095 
1096  static size_t objectCount;
1097  static std::vector<double> q_table;
1098  static std::vector<double> Fofq_table;
1099 
1101  //void initFromFile(float my_mass, long *seed, float vmax, float r_halfmass);
1103  //void initFromMassFunc(float my_mass, float my_Rsize, float my_rscale, PosType my_slope, long *seed);
1105  //void initFromMass(float my_mass, long *seed);
1106 
1108  void assignParams(InputParams& params);
1109 
1111  float sigma;
1113  float fratio;
1115  float pa;
1117  float rcore;
1118 
1120  PosType rmax_calc();
1121  void construct_ellip_tables();
1122 
1123  struct NormFuncer{
1124  NormFuncer(double my_q):q(my_q){};
1125 
1126  double operator()(double t){
1127  return 1.0/sqrt( 1 + (q*q-1)*sin(t)*sin(t));
1128  };
1129 
1130  private:
1131  double q;
1132  };
1133 
1134 };
1135 
1140 class LensHaloTNSIE : public LensHalo{
1141 public:
1142 
1143  LensHaloTNSIE(float my_mass
1144  ,PosType my_zlens
1145  ,float my_sigma
1146  ,float my_rcore
1147  ,float my_fratio
1148  ,float my_pa
1149  ,const COSMOLOGY &cosmo
1150  ,float f=20
1151  );
1152 
1153  LensHaloTNSIE(const LensHaloTNSIE &h):
1154  LensHalo(h)
1155  {
1156  sigma = h.sigma;
1157  fratio = h.fratio;
1158  pa = h.pa;
1159  rcore = h.rcore;
1160  rtrunc = h.rtrunc;
1161  units = h.units;
1162  }
1163 
1164  LensHaloTNSIE &operator=(const LensHaloTNSIE &h){
1165  if(&h == this) return *this;
1166  LensHalo::operator=(h);
1167  sigma = h.sigma;
1168  fratio = h.fratio;
1169  pa = h.pa;
1170  rcore = h.rcore;
1171  rtrunc = h.rtrunc;
1172  units = h.units;
1173  return *this;
1174  }
1175 
1176  ~LensHaloTNSIE(){};
1177 
1179  static double calc_sigma(float mass,float Rtrunc,float Rcore,float fratio){
1180  return lightspeed*sqrt( Grav*mass*sqrt(fratio) / (Rtrunc-Rcore) / PI);
1181  }
1182 
1184  void force_halo(PosType *alpha,KappaType *kappa,KappaType *gamma,KappaType *phi,PosType const *xcm,bool subtract_point=false,PosType screening = 1.0);
1185 
1187  float get_sigma(){return sigma;};
1188  // get the NSIE radius
1189  //float get_Rsize(){return Rsize;};
1191  float get_fratio(){return fratio;};
1193  float get_pa(){return pa;};
1195  float get_rcore(){return rcore;};
1197  float get_rtrunc(){return rtrunc;};
1198 
1199  void setZlens(PosType my_zlens,const COSMOLOGY &cosmo){
1200  LensHalo::setZlens(my_zlens,cosmo);
1201  }
1202 
1204  void set_pa(float my_pa){pa=my_pa;};
1205 
1206 protected:
1207 
1208  float units;
1209 
1211  void assignParams(InputParams& params);
1213  float sigma;
1215  float fratio;
1217  float pa;
1219  float rcore;
1221  float rtrunc;
1222 };
1223 
1229 class LensHaloTEPL : public LensHalo{
1230 public:
1231 
1232  LensHaloTEPL(float my_mass
1233  ,PosType my_zlens
1234  ,PosType r_trunc
1235  ,PosType gamma
1236  ,float my_fratio
1237  ,float my_pa
1238  ,const COSMOLOGY &cosmo
1239  ,float f=10
1240  );
1241 
1242  LensHaloTEPL(const LensHaloTEPL &h):
1243  LensHalo(h)
1244  {
1245  tt = h.tt;
1246  x_T = h.x_T;
1247  q = h.q;
1248  q_prime = h.q_prime;
1249  SigmaT = h.SigmaT;
1250  pa = h.pa;
1251  mass_pi = h.mass_pi;
1252  R = h.R;
1253  }
1254 
1255  LensHaloTEPL &operator=(const LensHaloTEPL &h){
1256  if(&h == this) return *this;
1257  LensHalo::operator=(h);
1258 
1259  tt = h.tt;
1260  x_T = h.x_T;
1261  q = h.q;
1262  q_prime = h.q_prime;
1263  SigmaT = h.SigmaT;
1264  pa = h.pa;
1265  mass_pi = h.mass_pi;
1266  R = h.R;
1267 
1268  return *this;
1269  }
1270 
1271  ~LensHaloTEPL(){};
1272 
1273 
1275  void force_halo(PosType *alpha,KappaType *kappa,KappaType *gamma,KappaType *phi,PosType const *xcm,bool subtract_point=false,PosType screening = 1.0);
1276 
1278  float get_fratio(){return q;};
1280  float get_pa(){return pa;};
1282  float get_rtrunc(){return x_T;};
1284  float get_t(){return tt;}
1285 
1286  void set_pa(double p){pa = p;}
1287 
1288  void deflection(std::complex<double> &z
1289  ,std::complex<double> &a
1290  ,std::complex<double> &g
1291  ,KappaType &sigma) const;
1292 
1293 
1294 protected:
1295 
1296  //float units;
1297  double tt;
1298  double x_T; // truncation radius
1299  double q;
1300  double q_prime;
1301  double SigmaT; // SigmaT - surface density at truncation radius
1302  double pa;
1303  double mass_pi;
1304  std::complex<double> R; // rotation
1305 
1306  std::complex<double> F(double r_e,double t,std::complex<double> z) const;
1307 };
1308 
1309 /*** \brief Truncated eliptical double power-law
1310 
1311  */
1312 class LensHaloTEBPL : public LensHalo{
1313 public:
1314 
1315  LensHaloTEBPL(float my_mass
1316  ,PosType my_zlens
1317  ,PosType r_break
1318  ,PosType r_trunc
1319  ,PosType t1
1320  ,PosType t2
1321  ,float my_fratio
1322  ,float my_pa
1323  ,const COSMOLOGY &cosmo
1324  ,float f=10
1325  );
1326 
1327  ~LensHaloTEBPL(){}
1328 
1329  LensHaloTEBPL(const LensHaloTEBPL &h):
1330  LensHalo(h),h1(h.h1),h2(h.h2),h3(h.h3)
1331  {
1332 
1333  m2 = h.m2;
1334  m3 = h.m3;
1335  m1 = h.m1;
1336 
1337  rb = h.rb;
1338  rt = h.rt;
1339  q = h.q;
1340 
1341  R = h.R;
1342  }
1343 
1344  LensHaloTEBPL &operator=(const LensHaloTEBPL &h){
1345  if(&h == this) return *this;
1346  LensHalo::operator=(h);
1347 
1348  h1 = h.h1;
1349  h2 = h.h2;
1350  h3 = h.h3;
1351 
1352  m2 = h.m2;
1353  m3 = h.m3;
1354  m1 = h.m1;
1355 
1356  rb = h.rb;
1357  rt = h.rt;
1358  q = h.q;
1359 
1360  R = h.R;
1361 
1362  return *this;
1363  }
1364 
1366  float get_fratio(){return q;};
1368  float get_pa(){return pa;};
1370  float get_rtrunc(){return rt;};
1372  float get_rbreak(){return rb;};
1374  float get_t1(){return h1.get_t();}
1376  float get_t2(){return h2.get_t();}
1377 
1378  void set_pa(double p){pa = p;}
1379 
1380  void force_halo(PosType *alpha,KappaType *kappa,KappaType *gamma,KappaType *phi,PosType const *xcm,bool subtract_point=false,PosType screening = 1.0);
1381 
1382 private:
1383 
1384  double q;
1385  double rb;
1386  double rt;
1387 
1388  double m2 ;
1389  double m3 ;
1390  double m1 ;
1391 
1392  LensHaloTEPL h1;
1393  LensHaloTEPL h2;
1394  LensHaloTEPL h3;
1395 
1396  std::complex<PosType> R;
1397 
1398 };
1399 
1400 #ifdef ENABLE_CERF
1401 /***
1402 \brief A class for making elliptical Gaussian lenses.
1403 
1404  This class uses the libcerf library (https://jugit.fz-juelich.de/mlz/libcerf).
1405  It can be installed with homebrew and it must be linked by setting the cmake variable ENABLE_CERF = ON
1406  */
1407 class LensHaloGaussian : public LensHalo{
1408 public:
1409 
1410  LensHaloGaussian(float my_mass
1411  ,PosType my_zlens
1412  ,PosType r_scale
1413  ,float my_fratio
1414  ,float my_pa
1415  ,const COSMOLOGY &cosmo
1416  ,float f=100
1417  );
1418 
1419  LensHaloGaussian(const LensHaloGaussian &h):
1420  LensHalo(h)
1421  {
1422  Rhight = h.Rhight;
1423  q = h.q;
1424  q_prime = h.q_prime;
1425  SigmaO = h.SigmaO;
1426  pa = h.pa;
1427  R = h.R;
1428 
1429  norm = h.norm;
1430  norm_g = h.norm_g;
1431  ss = h.ss;
1432  I = h.I;
1433  I_sqpi = h.I_sqpi;
1434  one_sqpi = h.one_sqpi;
1435  }
1436 
1437  LensHaloGaussian &operator=(const LensHaloGaussian &h){
1438  if(&h == this) return *this;
1439  LensHalo::operator=(h);
1440 
1441  Rhight = h.Rhight;
1442  q = h.q;
1443  q_prime = h.q_prime;
1444  SigmaO = h.SigmaO;
1445  pa = h.pa;
1446  R = h.R;
1447 
1448  norm = h.norm;
1449  norm_g = h.norm_g;
1450  ss = h.ss;
1451  I = h.I;
1452  I_sqpi = h.I_sqpi;
1453  one_sqpi = h.one_sqpi;
1454 
1455  return *this;
1456  }
1457 
1458  ~LensHaloGaussian(){};
1459 
1460 
1462  void force_halo(PosType *alpha,KappaType *kappa,KappaType *gamma,KappaType *phi,PosType const *xcm,bool subtract_point=false,PosType screening = 1.0);
1463 
1465  float get_fratio(){return q;};
1467  float get_pa(){return pa;};
1469  float get_scalehight(){return Rhight;};
1470 
1472  float get_SigmaCentral(){return SigmaO;};
1473 
1474  void set_pa(double p){pa = p;}
1475 
1476  void deflection(std::complex<double> &z
1477  ,std::complex<double> &a
1478  ,std::complex<double> &g
1479  ,KappaType &sigma) const;
1480 
1481 
1482 protected:
1483 
1484  std::complex<double> dwdz(std::complex<double> z) const{
1485  return 2.0*(I_sqpi - z*wF(z));
1486  };
1487 
1488  std::complex<double> wF(std::complex<double> z) const;
1489  std::complex<double> my_erfc(std::complex<double> z) const;
1490 
1491  // https://arxiv.org/pdf/1407.0748.pdf
1492 // std::complex<double> erfc(std::complex<double> z) const{
1493 // if(z.real() >= 0){
1494 // std::complex<double> ans = wf(I*z);
1495 // ans *= exp(-z*z);
1496  //return ans;
1497 // return exp(-z*z)*wf(I*z);
1498 // }else{
1499 // return 2. - exp(-z*z)*wf(-I*z);
1500 // }
1501 // }
1502 
1503  double Rhight; // scale hight in larges dimension
1504  double q;
1505  double q_prime;
1506  double SigmaO; // SigmaO - central surface density
1507  double pa;
1508  std::complex<double> R; // rotation
1509 
1510  double norm;
1511  double norm_g;
1512  double ss;
1513 
1514  //std::complex<double> F(double r_e,double t,std::complex<double> z) const;
1515  std::complex<double> erfi(std::complex<double> z) const{
1516  return I*(1. - exp(z*z)*wF(z));
1517  }
1518 
1519 // std::complex<double> rho(std::complex<double> z) const{
1520 //
1521 // //std::complex<double> y=(q*q*z.real() + I*z.imag())/ss/sqrt(1-q*q) ;
1522 // //return -I * ( wF( z/ss/sqrt(q_prime) )
1523 // // - exp( y*y -z*z/ss/ss/q_prime ) * wF(y));
1524 //
1525 // return wbar(z/ss/ss/sqrt(q_prime),1) - wbar(z/ss/ss/sqrt(q_prime),q);
1526 // }
1527 //
1528 // std::complex<double> wbar(std::complex<double> z, double p) const{
1529 // double x = z.real();
1530 // double y = z.imag();
1531 // std::cout << exp(-z*z) << "," <<
1532 // exp(-x*x*(1-p*p) - y*y*(1./p/p -1)) << "," << -I* wF(p*x + I*y/p) << std::endl;
1533 //
1534 // double tmp = exp(-x*x*(1-p*p) - y*y*(1./p/p -1));
1535 // if(tmp==0) return 0;
1536 //
1537 // return -(exp(-z*z)*I*tmp)*wF(p*x + I*y/p);
1538 // }
1539 
1540  // Faddeeva function according to Zaghloul (2017)
1541  // It oesn't seem right near the real axis for small |z|!
1542  // not used
1543 // std::complex<double> wf(std::complex<double> z) const{
1544 // //std::cout << z << std::endl;
1545 // double y2 = z.imag()*z.imag();
1546 // double z2 = std::norm(z);
1547 // if( z2 >= 3.8e4 ){
1548 // return I_sqpi / z;
1549 // }
1550 // if(z2 >= 256){
1551 // //std::cout << z << "," << (z*z -0.5) << "," << I * PI << std::endl;
1552 // return z/(z*z -0.5) * I_sqpi;
1553 // }
1554 // else if(z2 >= 62 && z2 < 256){
1555 // return (z*z - 1.) / (z * (z*z-1.5)) * I_sqpi;
1556 // }
1557 // else if(z2 < 62 && z2 >= 30 && y2 >= 1.0e-13){
1558 // return z*(z*z -2.5) / (z*z*(z*z-3.) + 0.75) * I_sqpi;
1559 // }
1560 // else if ( (z2 < 62 && z2 >= 30 && y2 < 1.0e-13) || (z2 < 30 && z2 >= 2.5 && y2 < 0.072) ){
1561 // return exp(-z*z) + I * z * (U1[5] + z*z*(U1[4] + z*z*(U1[3] + z*z*(U1[2]+ z*z*(U1[1] + z*z*(U1[0] + z*z*sqrt(PI) ))))))
1562 // /(V1[6] + z*z*(V1[5]+ z*z*(V1[4] + z*z*(V1[3] + z*z*(V1[2] + z*z*(V1[1] + z*z*(V1[0]+z*z)))))));
1563 // }
1564 //
1565 // return (U2[5]-I*z*(U2[4]-I*z*(U2[3]-I*z*(U2[2]-I*z*(U2[1]-I*z*(U2[0]-I*z*sqrt(PI) ))))))
1566 // /(V2[6] - I*z*(V2[5] - I*z*(V2[4] - I*z*(V2[3] - I*z*(V2[2]-I*z*(V2[1] - I*z*(V2[0] - I*z)))))));
1567 // }
1568 
1569  std::complex<double> I;
1570  std::complex<double> I_sqpi;
1571  double one_sqpi;
1572 // double U1[6] = {1.320522,35.7668,219.031,1540.787,3321.990,36183.31};
1573 // double V1[7] = {1.841439,61.57037,364.2191,2186.181,9022.228,24322.84,32066.6};
1574 //
1575 // double U2[6] = {5.9126262,30.180142,93.15558,181.92853,214.38239,122.60793};
1576 // double V2[7] = {10.479857,53.992907,170.35400,348.70392,457.33448,352.73063,122.60793};
1577 };
1578 
1579 #ifdef ENABLE_EIGEN
1580 
1582 namespace MultiGauss{
1583 
1584 struct PROFILE{
1585  virtual double operator()(double r) = 0;
1586  virtual double cum(double r) = 0;
1587 };
1588 
1589 struct POWERLAW : public PROFILE
1590 {
1591  POWERLAW(
1592  float g // power-law index of surfcae density
1593  ){
1594  if(g <= -2){
1595  std::cerr << "power-law index <= 2 give unbounded mass!" << std::endl;
1596  throw std::runtime_error("bad mass");
1597  }
1598  gamma = g;
1599  }
1600  double operator()(double r){return pow(r,gamma);}
1601  double cum(double r){return 2*PI*pow(r,gamma + 2)/(gamma+2);}
1602 
1603  void reset(float g){gamma = g;}
1604  float gamma;
1605 };
1606 
1607 struct SERSIC : public PROFILE
1608 {
1609  SERSIC(){reset(0,0);}
1610  SERSIC(float nn,float Re){
1611  reset(nn,Re);
1612  }
1613  double cum(double r){
1614  return cum_norm * boost::math::tgamma_lower(2*n, bn*pow(r/re,t) );
1615  }
1616  double operator()(double r){ //check that these are consistant
1617  return exp(-bn*pow(r/re,t) );
1618  }
1619 
1620  void reset(float nn,float Re){
1621  re = Re;
1622  n = nn;
1623  bn = 1.9992*n - 0.3271; // approximation valid for 0.5 < n < 8
1624  t = 1.0/n;
1625  cum_norm = 2*PI*n/pow(bn,2*n)*re*re;
1626  }
1627 
1628 private:
1629  float t;
1630  float re;
1631  float n;
1632  float bn;
1633  float cum_norm;
1634 };
1635 
1636 struct NFW: public PROFILE
1637 {
1638  NFW(float rs):rs(rs){
1639  }
1640  double cum(double r){
1641  double x =r/rs;
1642 
1643  double ans=log(x/2);
1644  if(x==1.0) ans += 1.0;
1645  if(x>1.0) ans += 2*atan(sqrt((x-1)/(x+1)))/sqrt(x*x-1);
1646  if(x<1.0) ans += 2*atanh(sqrt((1-x)/(x+1)))/sqrt(1-x*x);
1647 
1648  return 2*PI*rs*rs*ans;
1649  }
1650  double operator()(double r){ //check that these are consistant
1651  double x =r/rs;
1652 
1653  if(x==1.0){return 1.0/3.0;}
1654  if(x>1.0) return (1-2*atan(sqrt((x-1)/(x+1)))/sqrt(x*x-1))/(x*x-1);
1655  return (1-2*atanh(sqrt((1-x)/(x+1)))/sqrt(1-x*x))/(x*x-1);
1656  }
1657 
1658  void reset(double Rs){rs = Rs;}
1659 
1660  float rs;
1661 };
1662 }
1663 
1664 /***
1665 \brief A class for constructing and approximation to any elliptical profile out of a series of elliptical gaussians.
1666 
1667  The profile class must have two functions. The profile(r) must returns the surface density and profile.cum(r) must return the mass within the radius.
1668  Thier units are unimportant, but they must be consistant with eachother. Some implemented cases are MultiGauss::sersic, MultiGauss::powerlaw and MultiGauss::nfw
1669 
1670  The profile is fit Nradii points logarithmicly distributed between r_min and r_max using Ngaussians Gaussians in that range.
1671  Typically Nradii ~ 2 * Ngaussians.
1672 
1673  The mass is normalized so that mass_norm is within the elliptical distance Rnorm. The total mass with be calculated and can be recovered after construction.
1674 
1675  This class uses both the libcerf library (https://jugit.fz-juelich.de/mlz/libcerf).
1676  and the eigen library. They must be installed and linked using the cmake variable ENABLE_CERF = ON and ENABLE_EIGEN=ON
1677  */
1678 
1679 class LensHaloMultiGauss: public LensHalo{
1680 
1681 private:
1682 
1683  int nn; // number of gaussians
1684  int mm; // number of fit radii
1685  double q; // axis ratio
1686  double pa; // position angle
1687  double mass_norm;
1688  double r_norm;
1689  std::vector<double> sigmas;
1690  std::vector<double> radii;
1691  std::complex<double> Rotation;
1692  std::vector<LensHaloGaussian> gaussians;
1693  std::vector<double> A;
1694  float rms_error;
1695  //static int count;
1696 
1697 public:
1698 
1700  LensHaloMultiGauss(
1701  double mass_norm
1702  ,double Rnorm
1703  ,MultiGauss::PROFILE &profile
1704  ,int Ngaussians
1705  ,int Nradii
1706  ,PosType r_min
1707  ,PosType r_max
1708  ,PosType my_zlens
1709  ,float my_fratio
1710  ,float my_pa
1711  ,const COSMOLOGY &cosmo
1712  ,float f=10
1713  ,bool verbose = false
1714  );
1715 
1717  LensHaloMultiGauss(
1718  double my_mass_norm
1719  ,double Rnorm
1720  ,double my_scale // radial scale in units of the scale that was used to produce relative_scales
1721  ,const std::vector<double> &relative_scales
1722  ,const std::vector<double> &relative_masses
1723  ,PosType my_zlens
1724  ,float my_fratio
1725  ,float my_pa
1726  ,const COSMOLOGY &cosmo
1727  ,float f=100
1728  ,bool verbose = false
1729  );
1730 
1731 
1732 
1733  LensHaloMultiGauss(LensHaloMultiGauss &&halo);
1734  LensHaloMultiGauss(const LensHaloMultiGauss &halo);
1735 
1736  LensHaloMultiGauss & operator=(const LensHaloMultiGauss &&halo);
1737  LensHaloMultiGauss & operator=(const LensHaloMultiGauss &halo);
1738 
1739  ~LensHaloMultiGauss(){};//--LensHaloMultiGauss::count;}
1740 
1742  void set_pa(double my_pa);
1743 
1744  /***
1745  This is a static function that can be used to find the masses and scales for the gaussians which can
1746  then be fed into the second constructor with a rescaling. This avoids having to recalculate them when
1747  the same profile is used multiple times.
1748 
1749  For example, an NFW can be calculated here with a scale size 1 and an arbitrary mass.
1750  Then it can be used in the constructor with different scale sizes and masses.
1751 
1752  <p>
1753  MultiGauss::sersic profile(1,1);
1754  double mass_out;
1755  float rms_error;
1756  std::vector<double> scales,masses;
1757  LensHaloMultiGauss::calc_masses_scales(profile, 13, 25, 0.01,3,0.5,1,mass_out,scales,masses,rms_error);
1758  LensHaloMultiGauss halo1(mass,rscale,scales,masses,zl,q,-pa,cosmo);
1759  LensHaloMultiGauss halo2(mass/2,0.1*rscale,scales,masses,zl2,q2,-pa2,cosmo);
1760  <\p>
1761  */
1762 
1763  static void calc_masses_scales(MultiGauss::PROFILE &profile
1764  ,int Ngaussians
1765  ,int Nradii
1766  ,PosType r_min
1767  ,PosType r_max
1768  ,double mass_norm
1769  ,double r_norm
1770  ,double &totalmass_out
1771  ,std::vector<double> &scales
1772  ,std::vector<double> &masses
1773  ,float &rms_error
1774  ,bool verbose=false
1775  ){
1776  // make logorithmic
1777  scales.resize(Ngaussians);
1778  {
1779  double tmp = pow(r_max/r_min,1.0/(Ngaussians-2));
1780  scales[0] = r_min/5;
1781  scales[1] = r_min;
1782  for(int i=2 ; i < Ngaussians ; ++i){
1783  scales[i] = tmp*scales[i-1];
1784  }
1785  }
1786 
1787  std::vector<double> radii(Nradii);
1788  {
1789  double tmp = pow(r_max/r_min,1.0/(Nradii-1));
1790  radii[0] = r_min;
1791  for(int i=1 ; i < Nradii ; ++i){
1792  radii[i] = tmp*radii[i-1];
1793  }
1794  }
1795 
1796  Eigen::MatrixXd M(Nradii,Ngaussians);
1797 
1798  // first radius is fit to mass within radius
1799  {
1800  double cp = profile.cum(radii[0]);
1801  for(int n=0 ; n < Ngaussians ; ++n){
1802  M(0,n) = 2*PI*scales[n] * scales[n]
1803  * ( 1 -exp(-radii[0]*radii[0] / 2 / scales[n] /scales[n]) )
1804  / cp ;
1805  }
1806  }
1807 
1808  // all others are fit to relative density
1809  for(int m=1 ; m < Nradii ; ++m){
1810  double p = profile(radii[m]);
1811  for(int n=0 ; n < Ngaussians ; ++n){
1812  M(m,n) = exp( - radii[m]*radii[m] / 2 / scales[n] /scales[n])
1813  / p ;
1814  }
1815  }
1816 
1817  // invert M
1818 
1819  masses.resize(Ngaussians);
1820  {
1821  std::vector<double> b(Nradii,1);
1822  Eigen::Map<Eigen::VectorXd> v(b.data(),Nradii);
1823  Eigen::Map<Eigen::VectorXd> a(masses.data(),Ngaussians);
1824 
1825  //M = M.inverse();
1826  a = M.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(v);
1827 
1828  if(verbose) std::cout << "test in LensHaloMultiGauss::calc_masses_scales()" << std::endl << M*a << std::endl;
1829  }
1830 
1831  // tests
1832  {
1833  rms_error=0;
1834  for(int m=0 ; m < 1 ; ++m){
1835  double surf=0;
1836  for(int n=0 ; n < Ngaussians ; ++n){
1837  surf += masses[n] * scales[n] * scales[n]
1838  * 2*PI*( 1 - exp(-radii[m]*radii[m] / 2 / scales[n] /scales[n]) );
1839  }
1840  double tmp = profile.cum(radii[m]);
1841  if(verbose) std::cout << tmp << " " << surf << std::endl;
1842 
1843  rms_error += (surf - tmp )*(surf - tmp) / tmp / tmp;
1844  }
1845 
1846  for(int m=1 ; m < Nradii ; ++m){
1847  double surf=0;
1848  for(int n=0 ; n < Ngaussians ; ++n){
1849  surf += masses[n] * exp(-radii[m]*radii[m] / 2 / scales[n] /scales[n]) ;
1850  }
1851  double tmp = profile(radii[m]);
1852  if(verbose) std::cout << tmp << " " << surf << std::endl;
1853 
1854  rms_error += (surf - tmp )*(surf - tmp ) / tmp /tmp;
1855  }
1856  rms_error = sqrt(rms_error)/Ngaussians;
1857  if(verbose) std::cout << " MSE error " << rms_error << std::endl;
1858  }
1859 
1860  // find masses of Gaussian components
1861  {
1862  double mass_tmp = 0;
1863  for(int n=0 ; n<Ngaussians ; ++n){
1864  mass_tmp += scales[n]*scales[n]*masses[n]
1865  *( 1 - exp(-r_norm*r_norm / 2 / scales[n] /scales[n]) );
1866  }
1867 
1868  { // convert to mass units of each component
1869  double tmp = (mass_norm/mass_tmp);
1870  // rescale so total mass is correct
1871  for(int n=0 ; n<Ngaussians ; ++n) masses[n] *= tmp*scales[n]*scales[n];
1872  }
1873 
1874  }
1875  totalmass_out = 0;
1876  for(int n=0 ; n<Ngaussians ; ++n) totalmass_out += masses[n];
1877  }
1878 
1879  double profile(double r);
1880 
1882  double mass_cum(double r);
1883 
1885  float error();
1886 
1887  void force_halo(PosType *alpha,KappaType *kappa,KappaType *gamma,KappaType *phi,PosType const *xcm,bool subtract_point=false,PosType screening = 1);
1888 
1889 
1890 };
1891 
1892 #endif // eigen
1893 #endif // libcerf
1894 
1895 /*
1896  *
1897  * \brief A class for calculating the deflection, kappa and gamma caused by a halos
1898  * with truncated Hernquist mass profiles.
1899  *
1900  * The profile is \f$ \rho \propto \left( \frac{r}{r_s} \right)^{-1} \left( 1 + \frac{r}{r_s} \right)^{-3} \f$.
1901  *
1902  *
1903  */
1904 
1905 class LensHaloHernquist: public LensHalo{
1906 public:
1907  //LensHaloHernquist();
1908  LensHaloHernquist(float my_mass,float my_Rsize,PosType my_zlens,float my_rscale,float my_fratio,float my_pa,const COSMOLOGY &cosmo, EllipMethod my_ellip_method=EllipMethod::Pseudo);
1909  //LensHaloHernquist(InputParams& params);
1910  virtual ~LensHaloHernquist();
1911 
1912  PosType ffunction(PosType x) const;
1913  PosType gfunction(PosType x) const;
1914  PosType hfunction(PosType x) const;
1915  PosType g2function(PosType x) const;
1916 
1917  /* the below functions alphaHern, kappaHern and gammaHern are obsolete and better to be deleted to avoid confusion
1918  void alphaHern(PosType *alpha,PosType *x,PosType Rtrunc,PosType mass,PosType r_scale
1919  ,PosType *center,PosType Sigma_crit);
1920  KappaType kappaHern(PosType *x,PosType Rtrunc,PosType mass,PosType r_scale
1921  ,PosType *center,PosType Sigma_crit);
1922  void gammaHern(KappaType *gamma,PosType *x,PosType Rtrunc,PosType mass,PosType r_scale
1923  ,PosType *center,PosType Sigma_crit);
1924  */
1925  //void initFromFile(float my_mass, long *seed, float vmax, float r_halfmass);
1926 
1928  void set_rscale(float my_rscale){rscale=my_rscale; xmax = LensHalo::getRsize()/rscale; gmax = InterpolateFromTable(gtable,xmax);};
1929  // friend struct Ialpha_func;
1930 
1931 protected:
1933  static const long NTABLE;
1935  static const PosType maxrm;
1937  static int count;
1938 
1940  static PosType *ftable,*gtable,*g2table,*htable,*xtable,*xgtable;
1942  void make_tables();
1944  PosType InterpolateFromTable(PosType *table, PosType y) const;
1945 
1947  void assignParams(InputParams& params);
1948 
1949  // Override internal structure of halos
1951  inline PosType alpha_h(PosType x) const {
1952  return -0.25*InterpolateFromTable(gtable,x)/gmax;
1953  }
1954  inline KappaType kappa_h(PosType x) const {
1955  return 0.25*x*x*InterpolateFromTable(ftable,x)/gmax; // 0.5*
1956  }
1957  inline KappaType gamma_h(PosType x) const {
1958  return -0.5*x*x*InterpolateFromTable(g2table,x)/gmax;
1959  }
1960  inline KappaType phi_h(PosType x) const {
1961  return -0.25*InterpolateFromTable(htable,x)/gmax;
1962  //return -1.0*InterpolateFromTable(htable,x)/gmax;
1963  }
1964  inline KappaType phi_int(PosType x) const{
1965  return -0.25*InterpolateFromTable(xgtable,x)/gmax;
1966  }
1967 
1968 
1969 private:
1970  PosType gmax;
1971 };
1972 
1973 /*
1974  *
1975  * \brief A class for calculating the deflection, kappa and gamma caused by a collection of halos
1976  * with truncated Jaffe mass profiles.
1977  *
1978  * The profile is \f$ \rho \propto \left( \frac{r}{r_s} \right)^{-2} \left( 1 + \frac{r}{r_s} \right)^{-2} \f$ so beta would usually be negative.
1979  *
1980  *
1981  */
1982 
1983 class LensHaloJaffe: public LensHalo{
1984 public:
1985  //LensHaloJaffe();
1986  LensHaloJaffe(float my_mass,float my_Rsize,PosType my_zlens,float my_rscale,float my_fratio
1987  ,float my_pa,const COSMOLOGY &cosmo
1988  , EllipMethod my_ellip_method=EllipMethod::Pseudo);
1989  //LensHaloJaffe(InputParams& params);
1990  virtual ~LensHaloJaffe();
1991 
1993  void set_rscale(float my_rscale){rscale=my_rscale; xmax = LensHalo::getRsize()/rscale; gmax = InterpolateFromTable(gtable,xmax);};
1994 
1995  PosType ffunction(PosType x) const;
1996  PosType gfunction(PosType x) const;
1997  PosType hfunction(PosType x) const;
1998  PosType g2function(PosType x) const;
1999  PosType bfunction(PosType x);
2000  PosType dbfunction(PosType x);
2001  PosType ddbfunction(PosType x);
2002 
2003 protected:
2004 
2006  static const long NTABLE;
2008  static const PosType maxrm;
2010  static int count;
2011 
2013  static PosType *ftable,*gtable,*g2table,*htable,*xtable,*xgtable;
2015  void make_tables();
2017  PosType InterpolateFromTable(PosType *table, PosType y) const;
2018 
2020  void assignParams(InputParams& params);
2021 
2022  // Override internal structure of halos
2024  inline PosType alpha_h(PosType x) const{
2025  return -0.25*InterpolateFromTable(gtable,x)/gmax;
2026  }
2027  inline KappaType kappa_h(PosType x) const {
2028  return 0.125*x*x*InterpolateFromTable(ftable,x)/gmax;
2029  }
2030  inline KappaType gamma_h(PosType x) const {
2031  //return -0.125*x*x*InterpolateFromTable(g2table,x)/gmax;
2032  return -0.25*x*x*InterpolateFromTable(g2table,x)/gmax;
2033  }
2035  inline KappaType phi_h(PosType x) const {
2036  return -0.25*InterpolateFromTable(xgtable,x)/gmax;
2037  }
2038  inline KappaType phi_int(PosType x) const{
2039  return -0.25*InterpolateFromTable(xgtable,x)/gmax;
2040  }
2041 
2042 private:
2043  PosType gmax;
2044 
2045  // I have temporarily set these functions to 0 to make the code compile, Ben
2046  // PosType ffunction(PosType x){throw std::runtime_error("Set to temporary invalid value"); return 0;}
2047  // PosType gfunction(PosType x){throw std::runtime_error("Set to temporary invalid value"); return 0;}
2048  // PosType hfunction(PosType x){throw std::runtime_error("Set to temporary invalid value"); return 0;}
2049  // PosType g2function(PosType x){throw std::runtime_error("Set to temporary invalid value"); return 0;}
2050 
2051 };
2052 
2053 
2054 
2055 
2059 class LensHaloDummy: public LensHalo{
2060 public:
2061  LensHaloDummy();
2062  LensHaloDummy(float my_mass,float my_Rsize,PosType my_zlens,float my_rscale,const COSMOLOGY &cosmo);
2063  //LensHaloDummy(InputParams& params);
2064  ~LensHaloDummy(){};
2065 
2067  void force_halo(PosType *alpha,KappaType *kappa,KappaType *gamma,KappaType *phi,PosType const *xcm,bool subtract_point=false,PosType screening = 1.0);
2068  // void force_halo(PosType *alpha,KappaType *kappa,KappaType *gamma,PosType *xcm,bool subtract_point=false,PosType screening = 1.0);
2069 
2071  void initFromMassFunc(float my_mass, float my_Rsize, float my_rscale, PosType my_slope, long *seed);
2072 
2073 
2074 private:
2076  void assignParams(InputParams& params);
2077  inline PosType alpha_h(PosType x) const {return 0.;}
2078  inline KappaType kappa_h(PosType x) const {return 0.;}
2079  inline KappaType gamma_h(PosType x) const {return 0.;}
2080 };
2081 
2082 
2083 typedef LensHalo* LensHaloHndl;
2084 
2085 //bool LensHaloZcompare(LensHalo *lh1,LensHalo *lh2);//{return (lh1->getZlens() < lh2->getZlens());}
2086 //bool LensHaloZcompare(LensHalo lh1,LensHalo lh2){return (lh1.getZlens() < lh2.getZlens());}
2087 
2088 #endif /* LENS_HALOS_H_ */
LensHaloTNSIE::force_halo
void force_halo(PosType *alpha, KappaType *kappa, KappaType *gamma, KappaType *phi, PosType const *xcm, bool subtract_point=false, PosType screening=1.0)
overridden function to calculate the lensing properties
Definition: lens_halos.cpp:956
LensHaloRealNSIE::get_fratio
float get_fratio()
get the axis ratio
Definition: lens_halos.h:1070
LensHalo::operator=
LensHalo & operator=(const LensHalo &h)
Definition: lens_halos.cpp:169
LensHalo::test_average_gt
PosType test_average_gt(PosType R)
calculates the average gamma_t for LensHalo::test()
Definition: lens_halos.cpp:2450
LensHalo::kappa_asym
virtual PosType kappa_asym(PosType x, PosType theta)
Definition: base_analens.cpp:1033
LensHaloTNSIE::LensHaloTNSIE
LensHaloTNSIE(float my_mass, PosType my_zlens, float my_sigma, float my_rcore, float my_fratio, float my_pa, const COSMOLOGY &cosmo, float f=20)
Definition: lens_halos.cpp:927
LensHaloRealNSIE::rcore
float rcore
core size of NSIE
Definition: lens_halos.h:1117
LensHaloTNSIE::rcore
float rcore
core size of NSIE
Definition: lens_halos.h:1219
LensHaloRealNSIE::sigma
float sigma
velocity dispersion of NSIE
Definition: lens_halos.h:1111
LensHalo::get_mass
float get_mass() const
get the mass solar units
Definition: lens_halos.h:84
LensHaloTNSIE::set_pa
void set_pa(float my_pa)
set the position angle
Definition: lens_halos.h:1204
LensHaloRealNSIE::rmax_calc
PosType rmax_calc()
for the set fratio, sigma and rcore calculate the radius that contains the correct mass
Definition: lens_halos.cpp:1169
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:1520
LensHalo::setZlens
void setZlens(PosType my_zlens, const COSMOLOGY &cosmo)
set redshift
Definition: lens_halos.h:141
LensHaloTNSIE::sigma
float sigma
velocity dispersion of TNSIE
Definition: lens_halos.h:1213
LensHaloNFW::dmoddb
PosType dmoddb(int whichmod, PosType q, PosType b)
dbfunction and ddbfunction are approximation formulae for dbeta/dr and d^2beta/dr^2 whereas dbnum and...
Definition: nfw_lens.cpp:205
LensHaloRealNSIE::LensHaloRealNSIE
LensHaloRealNSIE(float my_mass, PosType my_zlens, float my_sigma, float my_rcore, float my_fratio, float my_pa, const COSMOLOGY &cosmo)
explicit constructor, Warning: If my_rcore > 0.0 and my_fratio < 1 then the mass will be somewhat les...
Definition: lens_halos.cpp:1058
LensHaloNFW::LensHaloNFW
LensHaloNFW()
Shell constructor. Sets the halo to zero mass.
Definition: lens_halos.cpp:365
LensHalo::get_Rmax
float get_Rmax() const
this can be used to tag types of LensHalos
Definition: lens_halos.h:80
EllipMethod
EllipMethod
Methods to make a previously isotropic halo elliptical.
Definition: InputParams.h:86
LensHalo::get_Nmod
static const int get_Nmod()
get length of mod array, which is Nmod. Not to be confused with getNmodes in the class LensHaloFit
Definition: lens_halos.h:219
LensHaloNFW::assignParams
void assignParams(InputParams &params)
read in parameters from a parameterfile in InputParams params
Definition: lens_halos.cpp:583
LensHaloNFW::kappaNFW
KappaType kappaNFW(PosType *x, PosType Rtrunc, PosType mass, PosType r_scale, PosType *center, PosType Sigma_crit)
Convergence for an NFW halo.
Definition: nfw_lens.cpp:36
LensHaloTEPL::get_rtrunc
float get_rtrunc()
get the truncation radius
Definition: lens_halos.h:1282
LensHaloRealNSIE::force_halo
void force_halo(PosType *alpha, KappaType *kappa, KappaType *gamma, KappaType *phi, PosType const *xcm, bool subtract_point=false, PosType screening=1.0)
overridden function to calculate the lensing properties
Definition: lens_halos.cpp:1912
LensHaloDummy::force_halo
void force_halo(PosType *alpha, KappaType *kappa, KappaType *gamma, KappaType *phi, PosType const *xcm, bool subtract_point=false, PosType screening=1.0)
overridden function to calculate the lensing properties
Definition: lens_halos.cpp:2364
LensHalo::force_halo_sym
void force_halo_sym(PosType *alpha, KappaType *kappa, KappaType *gamma, KappaType *phi, PosType const *xcm, bool subtract_point=false, PosType screening=1.0)
returns the lensing quantities of a ray in center of mass coordinates for a symmetric halo
Definition: lens_halos.cpp:1594
LensHalo::mass_norm_factor
PosType mass_norm_factor
This is Rsize/rscale !!
Definition: lens_halos.h:403
LensHaloDummy::initFromMassFunc
void initFromMassFunc(float my_mass, float my_Rsize, float my_rscale, PosType my_slope, long *seed)
initialize from a mass function
Definition: lens_halos.cpp:2354
LensHaloNFW::count
static int count
keeps track of how many time the tables are created, default is just once
Definition: lens_halos.h:831
LensHalo::setTheta
void setTheta(const Point_2d &p)
set the position of the Halo in radians
Definition: lens_halos.h:111
LensHalo::get_flag_elliptical
bool get_flag_elliptical()
flag=True if halo elliptical
Definition: lens_halos.h:160
LensHaloTNSIE::rtrunc
float rtrunc
core size of NSIE
Definition: lens_halos.h:1221
LensHalo::getTheta
void getTheta(PosType *MyPosHalo) const
get the position of the Halo in radians
Definition: lens_halos.h:113
LensHalo
A base class for all types of lensing "halos" which are any mass distribution that cause lensing.
Definition: lens_halos.h:57
LensHalo::set_switch_flag
void set_switch_flag(bool swt)
flag permits case distinction in force_halo_asym for elliptical NFWs only (get_switch_flag==true),...
Definition: lens_halos.h:163
LensHalo::set_RsizeRmax
virtual void set_RsizeRmax(float my_Rsize)
set Rsize (in Mpc) and reset Rmax
Definition: lens_halos.h:134
LensHaloNFW::InterpolateFromTable
PosType InterpolateFromTable(PosType *table, PosType y) const
interpolates from the specific tables
Definition: lens_halos.cpp:566
LensHalo::rscale
float rscale
scale length or core size. Different meaning in different cases. Not used in NSIE case.
Definition: lens_halos.h:360
LensHaloTEPL::LensHaloTEPL
LensHaloTEPL(float my_mass, PosType my_zlens, PosType r_trunc, PosType gamma, float my_fratio, float my_pa, const COSMOLOGY &cosmo, float f=10)
Definition: lens_halos.cpp:1188
LensHalo::getX
void getX(PosType *MyPosHalo) const
get the position of the Halo in physical Mpc on the lens plane
Definition: lens_halos.h:96
LensHaloPowerLaw
A class for calculating the deflection, kappa and gamma caused by a collection of halos with truncate...
Definition: lens_halos.h:957
LensHaloPseudoNFW::get_slope
PosType get_slope()
initialize from a mass function
Definition: lens_halos.h:899
LensHaloPseudoNFW::LensHaloPseudoNFW
LensHaloPseudoNFW()
shell constructor, should be avoided
Definition: lens_halos.cpp:646
LensHaloTEPL::get_pa
float get_pa()
get the position angle
Definition: lens_halos.h:1280
LensHalo::alpha_h
virtual PosType alpha_h(PosType x) const
Definition: lens_halos.h:381
Point_2d
Class for representing points or vectors in 2 dimensions. Not that the dereferencing operator is over...
Definition: point.h:48
LensHalo::analModes
void analModes(int modnumber, PosType my_beta, PosType q, PosType amod[3])
Definition: modes_ana.cpp:156
LensHalo::Nparams
virtual std::size_t Nparams() const
get the number of halo parameters
Definition: lens_halos.cpp:2406
LensHaloNFW::gammaNFW
void gammaNFW(KappaType *gamma, PosType *x, PosType Rtrunc, PosType mass, PosType r_scale, PosType *center, PosType Sigma_crit)
Shear for and NFW halo. this might have a flaw in it.
Definition: nfw_lens.cpp:56
LensHalo::getEllipMethod
EllipMethod getEllipMethod() const
stars
Definition: lens_halos.h:215
LensHaloRealNSIE::pa
float pa
position angle on sky, radians
Definition: lens_halos.h:1115
LensHaloTNSIE
Truncated non-singular isothermal ellipsoid.
Definition: lens_halos.h:1140
LensHaloPowerLaw::initFromMassFunc
void initFromMassFunc(float my_mass, float my_Rsize, float my_rscale, PosType my_slope, long *seed)
set the slope of the surface density profile
Definition: lens_halos.cpp:903
LensHalo::LensHalo
LensHalo()
Shell constructor.
Definition: lens_halos.cpp:14
LensHaloPseudoNFW
A class for calculating the deflection, kappa and gamma caused by a collection of halos with a double...
Definition: lens_halos.h:885
LensHaloNFW::maxrm
static const PosType maxrm
maximum Rsize/rscale
Definition: lens_halos.h:829
LensHaloRealNSIE::set_sigma
void set_sigma(float my_sigma)
set the velocity dispersion
Definition: lens_halos.h:1077
LensHaloRealNSIE::get_sigma
float get_sigma()
get the velocity dispersion
Definition: lens_halos.h:1066
LensHalo::alpha_int
PosType alpha_int(PosType x) const
Calculates potential (phi_int) from alpha_h. If flag is_alphah_a_table is True it takes and integrate...
Definition: base_analens.cpp:589
LensHalo::getParam
virtual PosType getParam(std::size_t p) const
get the value of a scaled halo parameter by index
Definition: lens_halos.cpp:2411
LensHaloTEPL
A truncated elliptical power-law profile.
Definition: lens_halos.h:1229
LensHalo::displayPos
void displayPos()
Definition: lens_halos.h:126
LensHaloTEPL::get_t
float get_t()
get the power-law index
Definition: lens_halos.h:1284
LensHalo::alphakappagamma_asym
virtual void alphakappagamma_asym(PosType x, PosType theta, PosType alpha[], PosType *kappa, PosType gamma[], PosType *phi)
Pseudo-elliptical profiles by Phi(G)-Ansatz.
Definition: base_analens.cpp:752
LensHalo::set_slope
virtual void set_slope(PosType my_slope)
set slope
Definition: lens_halos.h:156
LensHaloTNSIE::get_rcore
float get_rcore()
get the core radius
Definition: lens_halos.h:1195
LensHalo::compareZ
bool compareZ(PosType z)
force tree calculation for stars
Definition: lens_halos.h:200
LensHalo::set_rscale
virtual void set_rscale(float my_rscale)
set scale radius (in Mpc)
Definition: lens_halos.h:138
LensHalo::MassBy1DIntegation
PosType MassBy1DIntegation(PosType R)
calculates the mass within radius R by integating alpha on a ring and using Gauss' law,...
Definition: lens_halos.cpp:2443
LensHaloNFW::NTABLE
static const long NTABLE
table size
Definition: lens_halos.h:823
LensHalo::getRsize
float getRsize() const
get the Rsize which is the size of the halo in Mpc
Definition: lens_halos.h:82
LensHalo::fourier_coeff
double fourier_coeff(double n, double q, double beta)
Calculates fourier-coefficients for power law halo.
Definition: base_analens.cpp:583
LensHalo::set_rsize
void set_rsize(float my_rsize)
set radius rsize beyond which interpolation values between alpha_ellip and alpha_iso are computed
Definition: lens_halos.h:244
LensHaloPseudoNFW::set_slope
void set_slope(PosType my_slope)
set the slope of the surface density profile
Definition: lens_halos.h:897
LensHaloRealNSIE::get_rcore
float get_rcore()
get the core radius
Definition: lens_halos.h:1074
LensHaloRealNSIE::assignParams
void assignParams(InputParams &params)
initialize from a simulation file
Definition: lens_halos.cpp:1123
LensHaloTNSIE::get_pa
float get_pa()
get the position angle
Definition: lens_halos.h:1193
LensHalo::setTheta
void setTheta(PosType *PosXY)
set the position of the Halo in radians
Definition: lens_halos.h:109
LensHalo::setTheta
void setTheta(PosType PosX, PosType PosY)
set the position of the Halo in radians
Definition: lens_halos.h:107
LensHaloRealNSIE::get_pa
float get_pa()
get the position angle
Definition: lens_halos.h:1072
LensHaloTEPL::force_halo
void force_halo(PosType *alpha, KappaType *kappa, KappaType *gamma, KappaType *phi, PosType const *xcm, bool subtract_point=false, PosType screening=1.0)
overridden function to calculate the lensing properties
Definition: lens_halos.cpp:1224
LensHalo::map_variables
PixelMap< double > map_variables(LensingVariable lensvar, size_t Nx, size_t Ny, double res)
Map a PixelMap of the surface, density, potential and potential gradient centred on (0,...
Definition: lens_halos.cpp:269
LensHalo::printCSV
virtual void printCSV(std::ostream &, bool header=false) const
print the halo parameters in CSV format
Definition: lens_halos.cpp:2429
LensHaloTNSIE::get_fratio
float get_fratio()
get the axis ratio
Definition: lens_halos.h:1191
LensHalo::test
bool test()
perform some basic consistancy checks for halo
Definition: lens_halos.cpp:2462
LensHaloTNSIE::get_rtrunc
float get_rtrunc()
get the truncation radius
Definition: lens_halos.h:1197
LensHaloTNSIE::fratio
float fratio
axis ratio of surface mass distribution
Definition: lens_halos.h:1215
LensHaloNFW::set_RsizeXmax
void set_RsizeXmax(float my_Rsize)
set Rsize, xmax and gmax
Definition: lens_halos.h:820
LensHalo::setCosmology
virtual void setCosmology(const COSMOLOGY &cosmo)
used for elliptical NFWs only, in that case get_switch_flag==true
Definition: lens_halos.h:167
LensHalo::calcModes
void calcModes(double q, double beta, double rottheta, PosType newmod[])
Calculates the modes for fourier expansion of power law halo. All the modes are relative to the zero ...
Definition: base_analens.cpp:605
LensHalo::MassBy2DIntegation
PosType MassBy2DIntegation(PosType R)
Prints star parameters; if show_stars is true, prints data for single stars.
Definition: lens_halos.cpp:2437
LensHalo::alphakappagamma1asym
virtual void alphakappagamma1asym(PosType x, PosType theta, PosType alpha[2], PosType *kappa, PosType gamma[], PosType *phi)
Elliptical profiles by Fourier-Ansatz.
Definition: base_analens.cpp:798
LensHaloRealNSIE::fratio
float fratio
axis ratio of surface mass distribution
Definition: lens_halos.h:1113
LensHalo::assignParams
void assignParams(InputParams &params, bool needRsize)
read in parameters from a parameterfile in InputParams params
Definition: lens_halos.cpp:243
LensHalo::getZlens
PosType getZlens() const
get the redshift
Definition: lens_halos.h:88
LensHaloRealNSIE::set_fratio
void set_fratio(float my_fratio)
set the axis ratio
Definition: lens_halos.h:1081
LensHalo::get_slope
virtual PosType get_slope()
get slope
Definition: lens_halos.h:158
LensHaloDummy
This is a lens that does no lensing. It is useful for testing and for running refinement code on sour...
Definition: lens_halos.h:2059
LensHalo::getDist
PosType getDist() const
Definition: lens_halos.h:122
LensHaloTNSIE::assignParams
void assignParams(InputParams &params)
read-in parameters from a parameter file
LensHaloTNSIE::calc_sigma
static double calc_sigma(float mass, float Rtrunc, float Rcore, float fratio)
returns the sigma in km / s for the other parameters fixed
Definition: lens_halos.h:1179
InputParams.h
LensHalo::error_message1
void error_message1(std::string name, std::string filename)
read in star parameters. This is valid for all halos and not overloaded.
Definition: lens_halos.cpp:237
LensHalo::faxial
void faxial(PosType x, PosType theta, PosType f[])
If set to true the correct normalization is applied for asymmetric NFW profiles, the mass_norm_factor...
Definition: base_analens.cpp:470
LensHaloNFW::initFromFile
void initFromFile(float my_mass, long *seed, float vmax, float r_halfmass)
Sets the profile to match the mass, Vmax and R_halfmass.
Definition: lens_halos.cpp:618
LensHalo::set_mass
void set_mass(float my_mass)
set mass (in solar masses)
Definition: lens_halos.h:136
LensHaloTNSIE::get_sigma
float get_sigma()
get the velocity dispersion
Definition: lens_halos.h:1187
LensHalo::force_halo_asym
void force_halo_asym(PosType *alpha, KappaType *kappa, KappaType *gamma, KappaType *phi, PosType const *xcm, bool subtract_point=false, PosType screening=1.0)
Definition: lens_halos.cpp:1648
LensHalo::Rmax_to_Rsize_ratio
float Rmax_to_Rsize_ratio
The factor by which Rmax is larger than Rsize.
Definition: lens_halos.h:357
LensHalo::gradial
void gradial(PosType r, PosType g[])
Derivatives of the potential damping factor with respect to r ... TODO: come up with a better damping...
Definition: base_analens.cpp:554
LensHalo::get_rscale
float get_rscale() const
get the scale radius in Mpc
Definition: lens_halos.h:86
LensHalo::initFromMassFunc
virtual void initFromMassFunc(float my_mass, float my_Rsize, float my_rscale, PosType my_slope, long *seed)
initialize from a mass function
Definition: lens_halos.cpp:229
LensHaloPseudoNFW::mhat
PosType mhat(PosType y, PosType beta) const
Auxiliary function for PseudoNFW profile.
Definition: lens_halos.cpp:724
LensHaloNFW::alphaNFW
void alphaNFW(PosType *alpha, PosType *x, PosType Rtrunc, PosType mass, PosType r_scale, PosType *center, PosType Sigma_crit)
deflection caused by NFW halo
Definition: nfw_lens.cpp:9
LensHaloNFW::set_rscaleXmax
void set_rscaleXmax(float my_rscale)
Definition: lens_halos.h:823
LensHaloTEPL::get_fratio
float get_fratio()
get the axis ratio
Definition: lens_halos.h:1278
PixelMap< double >
LensHaloRealNSIE
Represents a non-singular isothermal elliptical lens.
Definition: lens_halos.h:1028
particle_types.h
LensingVariable
LensingVariable
output lensing variables
Definition: standard.h:89
LensHalo::initFromFile
virtual void initFromFile(float my_mass, long *seed, float vmax, float r_halfmass)
initialize from a simulation file
Definition: lens_halos.h:129
LensHalo::felliptical
void felliptical(PosType x, PosType q, PosType theta, PosType f[], PosType g[])
Calculate the derivatives of the G function = r*sqrt(cos(theta)^2 + q(r)^2 sin(theta))
Definition: base_analens.cpp:982
LensHaloNFW::make_tables
void make_tables()
make the specific tables
Definition: lens_halos.cpp:482
LensHaloNFW::ftable
static PosType * ftable
tables for lensing properties specific functions
Definition: lens_halos.h:834
LensHalo::get_mod
std::vector< double > get_mod()
get vector of Fourier modes, which are calculated in the constructors of the LensHaloes when main_ell...
Definition: lens_halos.h:217
LensHaloNFW
A class for calculating the deflection, kappa and gamma caused by an NFW halos.
Definition: lens_halos.h:749
LensHalo::setDist
void setDist(COSMOLOGY &co)
Set the angular size distance to the halo. This should be the distance to the lens plane.
Definition: lens_halos.h:118
LensHaloRealNSIE::set_rcore
void set_rcore(float my_rcore)
set the core radius Einstein radius
Definition: lens_halos.h:1085
InputParams
Structure for reading and writing parameters to and from a parameter file as well as a container for ...
Definition: InputParams.h:99
LensHaloTNSIE::pa
float pa
position angle on sky, radians
Definition: lens_halos.h:1217
LensHalo::operator[]
PosType operator[](int i) const
returns position of the Halo in physical Mpc on the lens plane
Definition: lens_halos.h:104
LensHalo::setParam
virtual PosType setParam(std::size_t p, PosType value)
set the value of a scaled halo parameter by index
Definition: lens_halos.cpp:2420
LensHaloPseudoNFW::initFromMassFunc
void initFromMassFunc(float my_mass, float my_Rsize, float my_rscale, PosType my_slope, long *seed)
set Rsize
Definition: lens_halos.cpp:770
LensHaloNFW::alpha_h
PosType alpha_h(PosType x) const
r |alpha(r = x*rscale)| PI Sigma_crit / Mass
Definition: lens_halos.h:846
LensHaloRealNSIE::set_pa
void set_pa(float my_pa)
set the position angle
Definition: lens_halos.h:1083
LensHaloNFW::initFromMassFunc
void initFromMassFunc(float my_mass, float my_Rsize, float my_rscale, PosType my_slope, long *seed)
initialize from a mass function
Definition: lens_halos.cpp:633
Generated on Sun Sep 27 2026 14:12:56 for GLAMERDOC++ by   doxygen 1.8.18