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