lens_halos.h

/*
 * lens_halos.h
 *
 *  Created on: 06.05.2013
 */
 
#ifndef LENS_HALOS_H_
#define LENS_HALOS_H_
 
//#include "standard.h"
#include "InputParams.h"
//#include "source.h"
//#include "point.h"
#include "quadTree.h"
#include "particle_types.h"
#include "image_processing.h"
 
#include <complex>
#include <complex.h>
#ifdef ENABLE_CERF
#include <cerf.h>
#endif
#ifdef ENABLE_EIGEN
//#include </usr/local/include/eigen3/Eigen/Dense>
#include <eigen3/Eigen/Dense>
#include <eigen3/Eigen/StdVector>
#include <boost/math/special_functions/gamma.hpp>
#endif

 
class LensHalo{
public:
 LensHalo();
  LensHalo(PosType z,const COSMOLOGY &cosmo);
  //LensHalo(InputParams& params,COSMOLOGY &cosmo,bool needRsize = true);
  //LensHalo(InputParams& params,bool needRsize = true);
  LensHalo(const LensHalo &h);
  
  LensHalo(LensHalo &&h){
//    LensHalo(LensHalo &&h):star_tree(nullptr){
    *this = std::move(h);
  }
 
 virtual ~LensHalo();
  
  LensHalo & operator=(const LensHalo &h);
  LensHalo & operator=(LensHalo &&h);
 
  int tag=0; 
  
  inline float get_Rmax() const { return Rmax; }
  inline float getRsize() const { return Rsize; }
  inline float get_mass() const { return mass; }
 inline float get_rscale() const { return rscale; }
 inline PosType getZlens() const { return zlens; }
    
  // set the position of the Halo in physical Mpc on the lens plane
  //void setX(PosType PosX, PosType PosY) { posHalo[0] = PosX ; posHalo[1] = PosY ; }
  // set the position of the Halo in physical Mpc on the lens plane
  //void setX(PosType *PosXY) { posHalo[0] = PosXY[0] ; posHalo[1] = PosXY[1] ; }
  
  void getX(PosType * MyPosHalo) const {
    assert(Dist != -1);
    MyPosHalo[0] = posHalo[0]*Dist;
    MyPosHalo[1] = posHalo[1]*Dist;
  }

  PosType operator[](int i) const{return posHalo[i]*Dist;}
  
  void setTheta(PosType PosX, PosType PosY) { posHalo[0] = PosX ; posHalo[1] = PosY ; }
  void setTheta(PosType *PosXY) { posHalo[0] = PosXY[0] ; posHalo[1] = PosXY[1] ; }
  void setTheta(const Point_2d &p) { posHalo[0] = p[0] ; posHalo[1] = p[1] ; }
  void getTheta(PosType * MyPosHalo) const { MyPosHalo[0] = posHalo[0] ; MyPosHalo[1] = posHalo[1]; }
  
  void setDist(COSMOLOGY &co){Dist = co.angDist(zlens);}

  PosType getDist() const {return Dist;}

  void displayPos() { std::cout << "Halo PosX = " << posHalo[0] << " ; Halo PosY = " << posHalo[1] << std::endl; }
  
 virtual void initFromFile(float my_mass, long *seed, float vmax, float r_halfmass){};
 virtual void initFromMassFunc(float my_mass, float my_Rsize, float my_rscale, PosType my_slope, long *seed);
  
  virtual void set_RsizeRmax(float my_Rsize){Rmax = Rmax*my_Rsize/Rsize; Rsize = my_Rsize; xmax = Rsize/rscale;};
 void set_mass(float my_mass){mass=my_mass;};
 virtual void set_rscale(float my_rscale){rscale=my_rscale; xmax = Rsize/rscale;};
  //void setZlens(PosType my_zlens){zlens=my_zlens; Dist=-1;}
  void setZlens(PosType my_zlens,const COSMOLOGY &cosmo){
    zlens=my_zlens;
    Dist=cosmo.angDist(my_zlens);
  }
  void setRsize(PosType R){Rsize = R;}
  
  // ste redshift and distance
  void setZlensDist(PosType my_zlens,const COSMOLOGY &cos){
    zlens=my_zlens;
    Dist = cos.angDist(zlens);
  }
  void setMass(PosType m){mass = m;}
  
  
 virtual void set_slope(PosType my_slope){beta=my_slope;};
  virtual PosType get_slope(){return beta;};
  bool get_flag_elliptical(){return elliptical_flag;};
  void set_flag_elliptical(bool ell){elliptical_flag=ell;};
  bool get_switch_flag(){return switch_flag;}; 
  void set_switch_flag(bool swt){switch_flag=swt;}; 
  
  
 virtual void setCosmology(const COSMOLOGY& cosmo) {}
 
 /* calculate the lensing properties -- deflection, convergence, and shear
   if not overwritten by derived class it uses alpha_h(), gamma_h(), etc. of the
   derived case or for a point source in this class
   
   xcm - the physical position on the lens plane relative to the center of the LensHalo in Mpc
  Units :
   
  ALPHA    -    mass/PhysMpc - ALPHA / Sig_crit / Dl is the deflection in radians
  KAPPA    -    surface mass density , mass / /PhysMpc/PhysMpc - KAPPA / Sig_crit is the convergence
  GAMMA    -    mass / /PhysMpc/PhysMpc - GAMMA / Sig_crit is the shear
  PHI      -    mass - PHI / Sig_crit is the lensing potential
   
* returns the lensing quantities of a ray in center of mass coordinates.
   *
   *  Warning: This adds to input value of alpha, kappa, gamma, and phi.  They need
   *  to be zeroed out if the contribution of just this halo is wanted.
   */
  virtual void force_halo(
      PosType *alpha          
      ,KappaType *kappa     
      ,KappaType *gamma     
      ,KappaType *phi       
      ,PosType const *xcm   
      ,bool subtract_point=false 
      ,PosType screening=1.0   
  );

 //void force_stars(PosType *alpha,KappaType *kappa,KappaType *gamma,PosType const *xcm);
  
 bool compareZ(PosType z){return z > zlens;};
  
//  bool AreStarsImplanted() const {return stars_implanted;}
//  void implant_stars(PosType **centers,int Nregions,long *seed, IMFtype type=One);
//  void implant_stars(PosType *center,long *seed,IMFtype type = One);
//  //void implant_stars(PosType *x,PosType *y,int Nregions,long *seed,IMFtype type=One);
//  void remove_stars();
//  IMFtype getStarIMF_type() const {return main_stars_imf_type;}
//  /// Fraction of surface density in stars
//  PosType getFstars() const {return star_fstars;}
//  /// The mass of the stars if they are all the same mass
//  PosType getStarMass() const {if(stars_implanted)return stars_xp[0].mass(); else return 0.0;}
  
  EllipMethod getEllipMethod() const {return main_ellip_method;}
  std::vector<double> get_mod() { std::vector<double> fmodes(Nmod); for(int i=0;i<Nmod;i++){fmodes[i]= mod[i] ;}  ;return fmodes;}
  const static int get_Nmod() {return Nmod;}
  
 virtual std::size_t Nparams() const;
 virtual PosType getParam(std::size_t p) const;
 virtual PosType setParam(std::size_t p, PosType value);
 
 virtual void printCSV(std::ostream&, bool header = false) const;
  
 //void PrintStars(bool show_stars);
  
  PosType MassBy2DIntegation(PosType R);
  PosType MassBy1DIntegation(PosType R);
  PosType test_average_gt(PosType R);
  PosType test_average_kappa(PosType R);
  
  // renomalize to make mass match
  
  void set_norm_factor(){mass_norm_factor=1;mass_norm_factor=mass/MassBy1DIntegation(Rsize);}
  
  void set_rsize(float my_rsize){ Rsize = my_rsize;};
 float get_rsize(){return Rsize;};
 
  // all of the following functions were used for Ansatz III w derivatives of the Fourier modes
  
  bool test();
  
  size_t getID() const {return idnumber;}
  void setID(size_t id){idnumber = id;}
  
  PosType renormalization(PosType r_max);
 
 
  PixelMap<double> map_variables(
                         LensingVariable lensvar 
                         ,size_t Nx
                         ,size_t Ny
                         ,double res             
  );
  
private:
  size_t idnumber; 
  PosType posHalo[2];
  PosType zlens;
 
protected:
  // This is the size of the halo beyond which it does not have the expected profile.
  float Rsize = 0;
 
  // total mass in Msun
  float mass;
  // angular size distance to this lens
  PosType Dist;
  PosType mnorm;
 
  // Beyond Rmax the halo will be treated as a point mass.  Between Rsize and Rmax the deflection
  // and shear are interpolated.  For circularly symmetric lenses Rmax should be equal to Rsize
  float Rmax;
 
 
  PosType alpha_int(PosType x) const;
  PosType norm_int(PosType r_max);
 
  void force_halo_sym(PosType *alpha,KappaType *kappa,KappaType *gamma,KappaType *phi,PosType const *xcm,bool subtract_point=false,PosType screening = 1.0);
  void force_halo_asym(PosType *alpha,KappaType *kappa,KappaType *gamma,KappaType *phi,PosType const *xcm,bool subtract_point=false,PosType screening = 1.0);
  
  bool force_point(PosType *alpha,KappaType *kappa,KappaType *gamma
                   ,KappaType *phi,PosType const *xcm,PosType rcm2
                   ,bool subtract_point,PosType screening);
  
  struct norm_func{
    norm_func(LensHalo& halo, PosType my_r_max): halo(halo), r_max(my_r_max){};
    LensHalo& halo;
    PosType r_max;
    //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);}
    PosType operator ()(PosType theta) {return halo.alpha_ell(r_max, theta);}
  };
  
  struct Ialpha_func{
    Ialpha_func(LensHalo& halo): halo(halo){};
    LensHalo& halo;
    PosType operator ()(PosType x) {return halo.alpha_h(x)/x ;}
  };
  
  struct Ig_func{
    Ig_func(const LensHalo& halo): halo(halo){};
    const LensHalo& halo;
    PosType operator ()(PosType x) {return halo.gfunction(x)/x ;}
  };
  
//  std::vector<IndexType> stars_index;
//  std::vector<StarType> stars_xp;
//  TreeQuadParticles<StarType> *star_tree;
//
//  int stars_N;
//  PosType star_massscale;
//  /// star masses relative to star_massscles
//  //float *star_masses;
//  PosType star_fstars;
//  PosType star_theta_force;
//  int star_Nregions;
//  std::vector<PosType> star_region;
  PosType beta;
//  void substract_stars_disks(PosType const *ray,PosType *alpha
//                             ,KappaType *kappa,KappaType *gamma);
//  std::vector<float> stellar_mass_function(IMFtype type, unsigned long Nstars, long *seed, PosType minmass=0.0, PosType maxmass=0.0
//                               ,PosType bendmass=0.0, PosType powerlo=0.0, PosType powerhi=0.0);
  
  
  void assignParams(InputParams& params,bool needRsize);
  //void assignParams_stars(InputParams& params);
  
  void error_message1(std::string name,std::string filename);
  
  
  float Rmax_to_Rsize_ratio = 1.2;
  
  float rscale;
  // redshift
  //PosType zlens;
 
  EllipMethod main_ellip_method;
 
//  bool stars_implanted;
//  /// Number of regions to be subtracted to compensate for the mass in stars
//  IMFtype main_stars_imf_type;
//  PosType main_stars_min_mass;
//  PosType main_stars_max_mass;
//  PosType bend_mstar;
//  PosType lo_mass_slope;
//  PosType hi_mass_slope;
//  /// parameters for stellar mass function: minimal and maximal stellar mass, bending point for a broken power law IMF
//  std::vector<PosType> star_Sigma;
//  std::vector<Point_2d> star_xdisk;
//
  
  virtual PosType inline alpha_h(PosType x) const {return -1;};
  virtual KappaType inline kappa_h(PosType x) const {return 0;};
  virtual KappaType inline gamma_h(PosType x) const {return -2;};
  virtual KappaType inline phi_h(PosType x) const {return 1;};
  virtual KappaType inline phi_int(PosType x) const {return 1;};
  virtual PosType inline ffunction(PosType x)const {return 0;};
  virtual PosType inline gfunction(PosType x) const {return -1;};
  virtual PosType inline dgfunctiondx(PosType x){return 0;};
  virtual PosType inline bfunction(PosType x){return -1;};
  virtual PosType inline dhfunction(PosType x) const {return 1;};
  virtual PosType inline ddhfunction(PosType x, bool numerical){return 0;};
  virtual PosType inline dddhfunction(PosType x, bool numerical){return 0;};
  virtual PosType inline bnumfunction(PosType x){return -1;};
  virtual PosType inline dbfunction(PosType x){return 0;};
  virtual PosType inline ddbfunction(PosType x){return 0;};
  virtual PosType inline dmoddb(int whichmod, PosType q, PosType b){return 0;};
  virtual PosType inline ddmoddb(int whichmod, PosType q, PosType b){return 0;};
  virtual PosType inline dmoddq(int whichmod, PosType q, PosType b){return 0;};
  virtual PosType inline ddmoddq(int whichmod, PosType q, PosType b){return 0;};
  
  PosType xmax;  
  
  PosType mass_norm_factor=1;
  
  // Functions for calculating axial dependence
  float pa;
  float fratio=1.0;
  bool elliptical_flag = false;
  bool switch_flag = false; 
  
  
  void faxial(PosType x,PosType theta,PosType f[]);
  void faxial0(PosType theta,PosType f0[]);
  void faxial1(PosType theta,PosType f1[]);
  void faxial2(PosType theta,PosType f2[]);
  void gradial(PosType r,PosType g[]);
  void gradial2(PosType r,PosType mu, PosType sigma,PosType g[]);
  
  void felliptical(PosType x, PosType q, PosType theta, PosType f[], PosType g[]);
  
  virtual void gamma_asym(PosType x,PosType theta, PosType gamma[]);
  virtual PosType kappa_asym(PosType x,PosType theta);
  virtual void alphakappagamma_asym(PosType x,PosType theta, PosType alpha[]
                                    ,PosType *kappa,PosType gamma[],PosType *phi);
  virtual void alphakappagamma1asym(PosType x,PosType theta, PosType alpha[2]
                                    ,PosType *kappa,PosType gamma[],PosType *phi);
  virtual void alphakappagamma2asym(PosType x,PosType theta, PosType alpha[2]
                                    ,PosType *kappa,PosType gamma[],PosType *phi);
  virtual void alphakappagamma3asym(PosType x,PosType theta, PosType alpha[2]
                                    ,PosType *kappa,PosType gamma[],PosType *phi);
  
  virtual PosType alpha_ell(PosType x,PosType theta);
  
  double fourier_coeff(double n, double q, double beta);
  double IDAXDM(double lambda, double a2, double b2, double x[], double rmax, double mo);
  double IDAYDM(double lambda, double a2, double b2, double x[], double rmax, double mo);
  double SCHRAMMKN(double n, double x[], double rmax);
  double SCHRAMMJN(double n, double x[], double rmax);
  double SCHRAMMI(double x[], double rmax);
 
  
  void calcModes(double q, double beta, double rottheta, PosType newmod[]);
  void calcModesB(PosType x, double q, double beta, double rottheta, PosType newmod[]);
  void calcModesC(PosType beta_r, double q, double rottheta, PosType newmod[]);
  
  virtual PosType inline InterpolateModes(int whichmod, PosType q, PosType b){return 0;};
  
  void analModes(int modnumber, PosType my_beta, PosType q, PosType amod[3]);
  
  struct fourier_func{
    fourier_func(double my_n, double my_q, double my_beta): n(my_n),q(my_q),beta(my_beta){};
    double n;
    double q;
    double beta;
    double operator ()(double theta) {return cos(n*theta)/pow(cos(theta)*cos(theta) + 1/q/q*sin(theta)*sin(theta),beta/2) ;}
  };
 
  struct SCHRAMMJ_func{
    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){};
    double n, *x, rmx;
    double operator ()(double u) {
      
      // PosType xisq=sqrt(u*(x[0]*x[0]+x[1]*x[1]/(1-(1-halo->fratio*halo->fratio)*u)));
      PosType xisq=sqrt((x[0]*x[0]+x[1]*x[1]/(1-u*(1-halo->fratio*halo->fratio))));
      
      if(xisq*xisq < 1e-20) xisq = 1.0e-10;
      KappaType kappa=halo->kappa_h(xisq)/PI/xisq/xisq;
      //std::cout << "Schramm: n=" << n << " " << u << " " << xisq << " " << halo->fratio << " " << rmx << " " << halo->Rmax << " " << halo->rscale << std::endl;
      return kappa*halo->mass/pow((1.-(1.-halo->fratio*halo->fratio)*u),n+0.5);
    }
  private:
    LensHalo *halo;
  };
  
  struct SCHRAMMK_func{
    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){};
    double n, *x, rmx;
    double operator ()(double u) {
      
     // PosType xisq=sqrt(u*(x[0]*x[0]+x[1]*x[1]/(1-(1-halo->fratio*halo->fratio)*u)));
      PosType xisq=sqrt((x[0]*x[0]+x[1]*x[1]/(1-u*(1-halo->fratio*halo->fratio))));
      if(xisq*xisq < 1e-20) xisq = 1.0e-10;
      PosType h=0.0001;
      PosType kp1=halo->kappa_h(xisq+h)/((xisq+h)*(xisq+h));
      PosType km1=halo->kappa_h(xisq-h)/((xisq-h)*(xisq-h));
      PosType kp2=halo->kappa_h(xisq+2*h)/((xisq+2*h)*(xisq+2*h));
      PosType km2=halo->kappa_h(xisq-2*h)/((xisq-2*h)*(xisq-2*h));
      PosType dkdxi=(8*kp1-8*km1-kp2+km2)/12/h/PI;
      
 
 
      //std::cout << "Schramm: n=" << n << " " << u << " " << xisq << " " << halo->fratio << " " << rmx << " " << halo->Rmax << " " << halo->rscale << std::endl;
      return u*dkdxi/pow((1.-(1.-halo->fratio*halo->fratio)*u),n+0.5);
      }
    private:
      LensHalo *halo;
    };
 
  
  struct SCHRAMMI_func{
    SCHRAMMI_func(double my_x[], double my_rmax, LensHalo *my_halo): x(my_x), rmx(my_rmax), halo(my_halo){};
    double *x, rmx;
    double operator ()(double u) {
      
      // PosType xisq=sqrt(u*(x[0]*x[0]+x[1]*x[1]/(1-(1-halo->fratio*halo->fratio)*u)));
      PosType xisq=sqrt((x[0]*x[0]+x[1]*x[1]/(1-u*(1-halo->fratio*halo->fratio))));
      
      if(xisq*xisq < 1e-20) xisq = 1.0e-10;
      PosType alpha=halo->alpha_h(xisq)/PI/xisq;
      //std::cout << "Schramm: n=" << n << " " << u << " " << xisq << " " << halo->fratio << " " << rmx << " " << halo->Rmax << " " << halo->rscale << std::endl;
      return xisq*alpha*halo->mass/pow((1.-(1.-halo->fratio*halo->fratio)*u),0.5);
    }
  private:
    LensHalo *halo;
  };
 
  
  struct IDAXDM_func{
    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){};
    double lambda, a2,b2, *x, rmx;
    double operator ()(double m) {
      double ap = m*m*a2 + lambda,bp = m*m*b2 + lambda;
      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
      PosType tmp = m*rmx;
      if(tmp*tmp < 1e-20) tmp = 1.0e-10;
      PosType xiso=tmp/halo->rscale;
      KappaType kappa=halo->kappa_h(xiso)/PI/xiso/xiso;
      
      PosType alpha[2]={0,0},tm[2] = {m*rmx,0};
      KappaType kappam=0,gamma[2]={0,0},phi=0;
      halo->force_halo_sym(alpha,&kappam,gamma,&phi,tm);
      
      
      //std::cout << "output x: " << m << " " << xiso << " " << m*kappa/(ap*ap*ap*bp*p2)*halo->mass << " " << m*kappam/(ap*ap*ap*bp*p2)<< std::endl;
      double integrandA=m*kappa/(ap*ap*ap*bp*p2)*halo->mass;
      double integrandB=m*kappam/(ap*ap*ap*bp*p2);
      //std::cout << integrandA-integrandB << std::endl;
      assert( std::abs((integrandA - integrandB)/integrandA)<1E-5);
 
      assert(kappa >= 0.0);
 
      //return m*kappa/(ap*ap*ap*bp*p2)*halo->mass; // integrand of equation (28) in Schramm 1990
      return integrandB;
    }
  private:
    LensHalo *halo;
  };
  
  struct IDAYDM_func{
    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){};
    double lambda, a2,b2, *x, rmx;
    double operator ()(double m) {
      double ap = m*m*a2 + lambda,bp = m*m*b2 + lambda;
      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
      PosType tmp = m*rmx;
      if(tmp*tmp < 1e-20) tmp = 1.0e-10;
      PosType xiso=tmp/halo->rscale;
      KappaType kappa=halo->kappa_h(xiso)/PI/xiso/xiso;
      
      PosType alpha[2]={0,0},tm[2] = {m*rmx,0};
      KappaType kappam=0,gamma[2]={0,0},phi=0;
      halo->force_halo_sym(alpha,&kappam,gamma,&phi,tm);
      
      double integrandA=m*kappa/(ap*ap*ap*bp*p2)*halo->mass;
      double integrandB=m*kappam/(ap*bp*bp*bp*p2);
      assert( std::abs((integrandA - integrandB)/integrandA)<1E-5);
 
      assert(kappa >= 0.0);
      //return m*kappa/(ap*bp*bp*bp*p2)*halo->mass; // integrand of equation (28) in Schramm 1990
      return integrandB;
    }
  private:
    LensHalo *halo;
  };
 
  
  const static int Nmod = 32;
  
  // Analytic description of Fourier modes
  
  PosType mod[Nmod];
  PosType mod1[Nmod];
  PosType mod2[Nmod];
  PosType r_eps;
  
  
  // These are stucts used in doing tests
  
  /*struct test_gt_func{
    test_gt_func(LensHalo& halo,PosType my_r): halo(halo),r(my_r){};
    LensHalo& halo;
    PosType r;
    PosType a[2] = {0,0},x[2] = {0,0};
    KappaType k = 0,g[3] = {0,0,0} ,p=0;
    //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));}
    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));}
    
  };*/
  
  
  struct test_gt_func{
    test_gt_func(PosType my_r,LensHalo *halo): r(my_r),halo(halo){};
    double operator ()(PosType t) {
      PosType alpha[2] = {0,0},x[2] = {0,0};
      KappaType kappa = 0,gamma[3] = {0,0,0},phi =0 ;
      x[0]=r*cos(t);
      x[1]=r*sin(t);
      halo->force_halo(alpha,&kappa,gamma,&phi,x);
      assert(gamma[0]==gamma[0]);
      assert(gamma[1]==gamma[1]);
       return (gamma[0]*cos(2*t)+gamma[1]*sin(2*t));}
    //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));}
  private:
    PosType r;
    LensHalo *halo;
  };
 
  
 /* struct test_kappa_func{
    test_kappa_func(LensHalo& halo,PosType my_r): halo(halo),r(my_r){};
    LensHalo& halo;
    PosType r;
    PosType a[2] = {0,0},x[2] = {0,0};
    KappaType k = 0,g[3] = {0,0,0} ,p=0;
    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; }
  };
*/
  
  struct test_kappa_func{
    test_kappa_func(PosType my_r,LensHalo *halo): r(my_r),halo(halo){};
    double operator ()(PosType t) {
      PosType alpha[2] = {0,0},x[2] = {0,0};
      KappaType kappa = 0,gamma[3] = {0,0,0},phi =0 ;
      x[0]=r*cos(t);
      x[1]=r*sin(t);
      halo->force_halo(alpha,&kappa,gamma,&phi,x);
      assert(kappa==kappa);
      return kappa; }
  private:
    PosType r;
    LensHalo *halo;
  };
  
  struct DMDTHETA{
    DMDTHETA(PosType R,LensHalo *halo): R(R),halo(halo){};
    PosType operator()(PosType theta){
      PosType alpha[2] = {0,0},x[2] = {0,0};
      KappaType kappa = 0,gamma[3] = {0,0,0},phi =0 ;
      
      x[0] = R*cos(theta);
      x[1] = R*sin(theta);
      
      halo->force_halo(alpha,&kappa,gamma,&phi,x);
      
      PosType alpha_r = -alpha[0]*cos(theta) - alpha[1]*sin(theta);
      //std::cout << alpha[0] << "  " << alpha[1] << std::endl;
      assert( alpha_r == alpha_r );
      //std::cout << theta << "  " << alpha_r << std::endl;
      return alpha_r;
    }
  private:
    PosType R;
    LensHalo *halo;
  };
  
  struct DMDRDTHETA{
    DMDRDTHETA(PosType R,LensHalo *halo): R(R),halo(halo){};
    PosType operator()(PosType theta){
      
      PosType alpha[2] = {0,0},x[2] = {0,0};
      KappaType kappa = 0,gamma[3] = {0,0,0} ,phi=0;
      
      x[0] = R*cos(theta);
      x[1] = R*sin(theta);
      
      halo->force_halo(alpha,&kappa,gamma,&phi,x);
      assert(kappa == kappa);
      assert(kappa != INFINITY);
      return kappa;
    }
  protected:
    PosType R;
    LensHalo *halo;
  };
  
  struct DMDR{
    DMDR(LensHalo *halo): halo(halo){};
    PosType operator()(PosType logR){
       //if(halo->get_flag_elliptical()){
        LensHalo::DMDRDTHETA dmdrdtheta(exp(logR),halo);
        //std::cout << " R = " << exp(logR) << std::endl;
        
        if(exp(2*logR) == 0.0) return 0.0;
         return Utilities::nintegrate<LensHalo::DMDRDTHETA,PosType>(dmdrdtheta,0,2*PI,1.0e-7)
         *exp(2*logR);
//      }else{
//        PosType alpha[2] = {0,0},x[2] = {0,0};
//        KappaType kappa = 0,gamma[3] = {0,0,0} ,phi=0;
//
//        x[0] = exp(logR);
//        x[1] = 0;
//
//        halo->force_halo(alpha,&kappa,gamma,&phi,x);
//        return 2*PI*kappa*exp(2*logR);
//      }
    }
  protected:
    LensHalo *halo;
    
  };
 
  struct DALPHAXDM{
    DALPHAXDM(PosType lambda,PosType a2,PosType b2,PosType X[],LensHalo* myisohalo)
    :lambda(lambda),a2(a2),b2(b2),x(X),isohalo(myisohalo){};
    
    PosType operator()(PosType m);
  private:
    PosType lambda;
    PosType a2;
    PosType b2;
    PosType *x;
    LensHalo* isohalo;
  };
  struct DALPHAYDM{
    DALPHAYDM(PosType lambda,PosType a2,PosType b2,PosType X[],LensHalo* isohalo)
    :lambda(lambda),a2(a2),b2(b2),x(X),isohalo(isohalo){};
    
    PosType operator()(PosType m);
  private:
    PosType lambda;
    PosType a2;
    PosType b2;
    PosType *x;
    LensHalo* isohalo;
  };
};

class LensHaloNFW: public LensHalo{
public:
 LensHaloNFW();
  LensHaloNFW(float my_mass   
              ,float my_Rsize  
              ,PosType my_zlens   
              ,float my_concentration
              ,float my_fratio    
              ,float my_pa        
              ,const COSMOLOGY &cosmo
              ,EllipMethod my_ellip_method=EllipMethod::Pseudo
              );
 
  LensHaloNFW(const LensHaloNFW &h):LensHalo(h){
    ++LensHaloNFW::count;
    gmax = h.gmax;
  }
  LensHaloNFW(const LensHaloNFW &&h):LensHalo(std::move(h)){
    ++LensHaloNFW::count;
    gmax = h.gmax;
  }
 
  LensHaloNFW & operator=(const LensHaloNFW &h){
    if(this != &h){
      LensHalo::operator=(h);
      gmax = h.gmax;
    }
    return *this;
  }
  LensHaloNFW & operator=(const LensHaloNFW &&h){
    if(this != &h){
      LensHalo::operator=(std::move(h));
      gmax = h.gmax;
    }
    return *this;
  }
  
 virtual ~LensHaloNFW();
  
 PosType ffunction(PosType x) const;
 PosType gfunction(PosType x) const;
  PosType dgfunctiondx(PosType x);
 PosType g2function(PosType x) const;
 PosType hfunction(PosType x) const;
  PosType dhfunction(PosType x) const;
  PosType ddhfunction(PosType x, bool numerical);
  PosType dddhfunction(PosType x, bool numerical);
  PosType bfunction(PosType x);
  PosType dbfunction(PosType x);
  PosType ddbfunction(PosType x);
  PosType dmoddb(int whichmod, PosType q, PosType b);
  PosType ddmoddb(int whichmod, PosType q, PosType b);
  PosType dmoddq(int whichmod, PosType q, PosType b);
  PosType ddmoddq(int whichmod, PosType q, PosType b);
  
  //PosType dmod(PosType x, int modnumber, PosType my_slope, PosType my_fratio);    // was used for Ansatz III w derivatives of the Fourier modes
  //PosType ddmod(PosType x, int modnumber, PosType my_slope, PosType my_fratio);   // was used for Ansatz III w derivatives of the Fourier modes
  
  
 // TODO: BEN: the below functions alphaNFW, kappaNFW and gammaNFW are obsolete and better to be deleted to avoid confusion
 void alphaNFW(PosType *alpha,PosType *x,PosType Rtrunc,PosType mass,PosType r_scale
                ,PosType *center,PosType Sigma_crit);
 KappaType kappaNFW(PosType *x,PosType Rtrunc,PosType mass,PosType r_scale
                     ,PosType *center,PosType Sigma_crit);
 void gammaNFW(KappaType *gamma,PosType *x,PosType Rtrunc,PosType mass,PosType r_scale
                ,PosType *center,PosType Sigma_crit);
  
 void initFromFile(float my_mass, long *seed, float vmax, float r_halfmass);
 void initFromMassFunc(float my_mass, float my_Rsize, float my_rscale, PosType my_slope, long *seed);
  void set_RsizeXmax(float my_Rsize){LensHalo::setRsize(my_Rsize); xmax = LensHalo::getRsize()/rscale; gmax = InterpolateFromTable(gtable,xmax);};
 void set_rscaleXmax(float my_rscale){rscale=my_rscale; xmax = LensHalo::getRsize()/rscale; gmax = InterpolateFromTable(gtable,xmax);};
  
protected:
 static const long NTABLE;
 static const PosType maxrm;
 static int count;
  
 static PosType *ftable,*gtable,*g2table,*htable,*xtable,*xgtable,***modtable; // modtable was used for Ansatz IV and worked well
 void make_tables();
 PosType InterpolateFromTable(PosType *table, PosType y) const;
  PosType InterpolateModes(int whichmod, PosType q, PosType b);
  
 void assignParams(InputParams& params);
  
 // Override internal structure of halos
 inline PosType alpha_h(PosType x) const {
  //return -1.0*InterpolateFromTable(gtable,x)/InterpolateFromTable(gtable,xmax);
  return -1.0*InterpolateFromTable(gtable,x)/gmax;
    // return -0.5/x*InterpolateFromTable(gtable,x)/gmax;
 }
 inline KappaType kappa_h(PosType x) const{
  return 0.5*x*x*InterpolateFromTable(ftable,x)/gmax;
 }
 inline KappaType gamma_h(PosType x) const{
  //return -0.25*x*x*InterpolateFromTable(g2table,x)/gmax;
  return -0.5*x*x*InterpolateFromTable(g2table,x)/gmax;
 }
 inline KappaType phi_h(PosType x) const{
    return 0.25*(InterpolateFromTable(htable,x) - InterpolateFromTable(htable,LensHalo::getRsize()/rscale))/gmax + log(LensHalo::getRsize()) ;
    // The constant contribution is made to match with the point mass at x = Rsize/rscale.
 }
  inline KappaType phi_int(PosType x) const{
    return -1.0*InterpolateFromTable(xgtable,x)/gmax; //alpha_int(x);
  }
  
private:
  PosType gmax;
};
// ********************
 

class LensHaloPseudoNFW: public LensHalo{
public:
 LensHaloPseudoNFW();
  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);
 //LensHaloPseudoNFW(InputParams& params);
 ~LensHaloPseudoNFW();
  
 PosType mhat(PosType y, PosType beta) const;
  PosType gfunction(PosType x) const;
  
 void set_slope(PosType my_slope){beta=my_slope; make_tables();};
  PosType get_slope(){return beta;};
  //void set_Rsize(float my_Rsize){Rsize = my_Rsize; xmax = Rsize/rscale;};
 
 void initFromMassFunc(float my_mass, float my_Rsize, float my_rscale, PosType my_slope, long *seed);
  
private:
 static const long NTABLE;
 static const PosType maxrm;
 static int count;
  
 static PosType *mhattable,*xtable;
 void make_tables();
 PosType InterpolateFromTable(PosType y) const;
  
 void assignParams(InputParams& params);
  
 PosType beta;
  
 // Override internal structure of halos
 inline PosType alpha_h(PosType x) const {
  return -1.0*InterpolateFromTable(x)/InterpolateFromTable(xmax);
 }
 inline KappaType kappa_h(PosType x) const {
  return 0.5*x*x/InterpolateFromTable(xmax)/pow(1+x,beta);
 }
 inline KappaType gamma_h(PosType x) const {
  //return (0.5*x*x/pow(1+x,beta) - InterpolateFromTable(x))/InterpolateFromTable(xmax);
  return (x*x/pow(1+x,beta) - 2*InterpolateFromTable(x))/InterpolateFromTable(xmax);
 }
 inline KappaType phi_h(PosType x) const {
  return -1.0*alpha_int(x)/InterpolateFromTable(xmax) ;
    //ERROR_MESSAGE();
  //std::cout << "time delay has not been fixed for PseudoNFW profile yet." << std::endl;
  //exit(1);
  //return 0.0;
 }
    inline KappaType phi_int(PosType x) const{
        return -1.0*alpha_int(x)/InterpolateFromTable(xmax) ;
    }
};
 

class LensHaloPowerLaw: public LensHalo{
public:
 LensHaloPowerLaw();
  LensHaloPowerLaw(float my_mass,float my_Rsize,PosType my_zlens,PosType my_beta
                   ,float my_fratio,float my_pa,const COSMOLOGY &cosmo
                   ,EllipMethod my_ellip_method=EllipMethod::Pseudo);
 //LensHaloPowerLaw(InputParams& params);
 ~LensHaloPowerLaw();
  
 //void set_slope(PosType my_slope){beta=my_slope;};
  
  //PosType get_slope(){return beta;};
  
 void initFromMassFunc(float my_mass, float my_Rsize, float my_rscale, PosType my_slope, long *seed);
  
private:
 void assignParams(InputParams& params);
  
 // PosType beta;
  
 // Override internal structure of halos
 inline PosType alpha_h(
                         PosType x  
                         ) const{
  if(x<xmax*1.0e-6) x=1e-6*xmax;
  return -1.0*pow(x/xmax,-beta+2);
 }
 inline KappaType kappa_h(
                           PosType x   
                           ) const {
  if(x==0) x=1e-6*xmax;
    double tmp = x/xmax;
  return 0.5*(-beta+2)*pow(tmp,2-beta);
 }
 inline KappaType gamma_h(PosType x) const {
  if(x==0) x=1e-6*xmax;
    return -beta*pow(x/xmax,-beta+2);
 }
  inline KappaType phi_h(PosType x) const {
  if(x==0) x=1e-6*xmax;
    return ( pow(x/xmax,2-beta) - 1 )/(2-beta) + log(LensHalo::getRsize()) ;
 }
  inline KappaType phi_int(PosType x) const{
  //return alpha_int(x);
    return -1.0*pow(x/xmax,-beta+2)/(2-beta);
  }
};
 

class LensHaloRealNSIE : public LensHalo{
public:

  LensHaloRealNSIE(float my_mass,PosType my_zlens,float my_sigma
                   ,float my_rcore,float my_fratio,float my_pa,const COSMOLOGY &cosmo);
  
  // Warning: If my_rcore > 0.0 and my_fratio < 1 then the mass will be somewhat less than my_mass.
 //LensHaloRealNSIE(InputParams& params);
  
  LensHaloRealNSIE(const LensHaloRealNSIE &h):
  LensHalo(h)
  {
    ++objectCount;
    sigma = h.sigma;
    fratio = h.fratio;
    pa = h.pa;
    rcore = h.rcore;
    units = h.units;
  }
  
  LensHaloRealNSIE &operator=(const LensHaloRealNSIE &h){
    if(&h == this) return *this;
    LensHalo::operator=(h);
    sigma = h.sigma;
    fratio = h.fratio;
    pa = h.pa;
    rcore = h.rcore;
    units = h.units;
    return *this;
  }
  
 ~LensHaloRealNSIE();
  
 void force_halo(PosType *alpha,KappaType *kappa,KappaType *gamma,KappaType *phi,PosType const *xcm,bool subtract_point=false,PosType screening = 1.0);
 
 float get_sigma(){return sigma;};
 // get the NSIE radius
 //float get_Rsize(){return Rsize;};
 float get_fratio(){return fratio;};
 float get_pa(){return pa;};
 float get_rcore(){return rcore;};
  
 void set_sigma(float my_sigma){sigma=my_sigma;};
 // set the NSIE radius
 //void set_Rsize(float my_Rsize){Rsize=my_Rsize;};
 void set_fratio(float my_fratio){fratio=my_fratio;};
 void set_pa(float my_pa){pa=my_pa;};
 void set_rcore(float my_rcore){rcore=my_rcore;};
  
  void setZlens(PosType my_zlens,const COSMOLOGY &cosmo){
    LensHalo::setZlens(my_zlens,cosmo);
  }
 
  
protected:
  
  float units;
  
  static size_t objectCount;
  static std::vector<double> q_table;
  static std::vector<double> Fofq_table;
  
 //void initFromFile(float my_mass, long *seed, float vmax, float r_halfmass);
 //void initFromMassFunc(float my_mass, float my_Rsize, float my_rscale, PosType my_slope, long *seed);
 //void initFromMass(float my_mass, long *seed);
  
 void assignParams(InputParams& params);
  
 float sigma;
 float fratio;
 float pa;
 float rcore;
  
  PosType rmax_calc();
  void construct_ellip_tables();
  
  struct NormFuncer{
    NormFuncer(double my_q):q(my_q){};
    
    double operator()(double t){
      return 1.0/sqrt( 1 + (q*q-1)*sin(t)*sin(t));
    };
    
  private:
    double q;
  };
  
};

class LensHaloTNSIE : public LensHalo{
public:
 
  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 
                );
  
  LensHaloTNSIE(const LensHaloTNSIE &h):
  LensHalo(h)
  {
    sigma = h.sigma;
    fratio = h.fratio;
    pa = h.pa;
    rcore = h.rcore;
    rtrunc = h.rtrunc;
    units = h.units;
  }
  
  LensHaloTNSIE &operator=(const LensHaloTNSIE &h){
    if(&h == this) return *this;
    LensHalo::operator=(h);
    sigma = h.sigma;
    fratio = h.fratio;
    pa = h.pa;
    rcore = h.rcore;
    rtrunc = h.rtrunc;
    units = h.units;
    return *this;
  }
  
  ~LensHaloTNSIE(){};
  
  static double calc_sigma(float mass,float Rtrunc,float Rcore,float fratio){
    return lightspeed*sqrt( Grav*mass*sqrt(fratio) / (Rtrunc-Rcore) / PI);
  }

  void force_halo(PosType *alpha,KappaType *kappa,KappaType *gamma,KappaType *phi,PosType const *xcm,bool subtract_point=false,PosType screening = 1.0);
  
  float get_sigma(){return sigma;};
  // get the NSIE radius
  //float get_Rsize(){return Rsize;};
  float get_fratio(){return fratio;};
  float get_pa(){return pa;};
  float get_rcore(){return rcore;};
  float get_rtrunc(){return rtrunc;};
  
  void setZlens(PosType my_zlens,const COSMOLOGY &cosmo){
    LensHalo::setZlens(my_zlens,cosmo);
  }

  void set_pa(float my_pa){pa=my_pa;};
  
protected:
  
  float units;
  
  void assignParams(InputParams& params);
  float sigma;
  float fratio;
  float pa;
  float rcore;
  float rtrunc;
};

 
class LensHaloTEPL : public LensHalo{
public:
 
  LensHaloTEPL(float my_mass  
                ,PosType my_zlens 
                ,PosType r_trunc  
                ,PosType gamma    
                ,float my_fratio 
                ,float my_pa     
                ,const COSMOLOGY &cosmo  
                ,float f=10 
  );
  
  LensHaloTEPL(const LensHaloTEPL &h):
  LensHalo(h)
  {
    tt = h.tt;
    x_T = h.x_T;
    q = h.q;
    q_prime = h.q_prime;
    SigmaT = h.SigmaT;
    pa = h.pa;
    mass_pi = h.mass_pi;
    R = h.R;
  }
  
  LensHaloTEPL &operator=(const LensHaloTEPL &h){
    if(&h == this) return *this;
    LensHalo::operator=(h);
 
    tt = h.tt;
    x_T = h.x_T;
    q = h.q;
    q_prime = h.q_prime;
    SigmaT = h.SigmaT;
    pa = h.pa;
    mass_pi = h.mass_pi;
    R = h.R;
 
    return *this;
  }
  
  ~LensHaloTEPL(){};
 

  void force_halo(PosType *alpha,KappaType *kappa,KappaType *gamma,KappaType *phi,PosType const *xcm,bool subtract_point=false,PosType screening = 1.0);
  
  float get_fratio(){return q;};
  float get_pa(){return pa;};
  float get_rtrunc(){return x_T;};
  float get_t(){return tt;}
  
  void set_pa(double p){pa = p;}
  
  void deflection(std::complex<double> &z
                  ,std::complex<double> &a
                  ,std::complex<double> &g
                  ,KappaType &sigma) const;
  
  
protected:
  
  //float units;
  double tt;
  double x_T;       // truncation radius
  double q;
  double q_prime;
  double SigmaT;   // SigmaT - surface density at truncation radius
  double pa;
  double mass_pi;
  std::complex<double> R; // rotation
  
  std::complex<double> F(double r_e,double t,std::complex<double> z) const;
};
 
/*** \brief Truncated eliptical double power-law
 
 */
class LensHaloTEBPL : public LensHalo{
public:
  
  LensHaloTEBPL(float my_mass  
                ,PosType my_zlens 
                ,PosType r_break  
                ,PosType r_trunc  
                ,PosType t1    
                ,PosType t2    
                ,float my_fratio 
                ,float my_pa     
                ,const COSMOLOGY &cosmo  
                ,float f=10 
  );
 
  ~LensHaloTEBPL(){}
  
  LensHaloTEBPL(const LensHaloTEBPL &h):
  LensHalo(h),h1(h.h1),h2(h.h2),h3(h.h3)
  {
 
    m2 = h.m2;
    m3 = h.m3;
    m1 = h.m1;
    
    rb = h.rb;
    rt = h.rt;
    q = h.q;
 
    R = h.R;
  }
  
  LensHaloTEBPL &operator=(const LensHaloTEBPL &h){
    if(&h == this) return *this;
    LensHalo::operator=(h);
 
    h1 = h.h1;
    h2 = h.h2;
    h3 = h.h3;
    
    m2 = h.m2;
    m3 = h.m3;
    m1 = h.m1;
    
    rb = h.rb;
    rt = h.rt;
    q = h.q;
 
    R = h.R;
 
    return *this;
  }
  
  float get_fratio(){return q;};
  float get_pa(){return pa;};
  float get_rtrunc(){return rt;};
  float get_rbreak(){return rb;};
  float get_t1(){return h1.get_t();}
  float get_t2(){return h2.get_t();}
 
  void set_pa(double p){pa = p;}
 
  void force_halo(PosType *alpha,KappaType *kappa,KappaType *gamma,KappaType *phi,PosType const *xcm,bool subtract_point=false,PosType screening = 1.0);
  
private:
 
  double q;
  double rb;
  double rt;
 
  double m2 ;
  double m3 ;
  double m1 ;
 
  LensHaloTEPL h1;
  LensHaloTEPL h2;
  LensHaloTEPL h3;
  
  std::complex<PosType> R;
  
};
 
#ifdef ENABLE_CERF
/***
\brief  A class for making elliptical Gaussian lenses.
 
 This class uses the libcerf library (https://jugit.fz-juelich.de/mlz/libcerf).
 It can be installed with homebrew and it must be linked by setting the cmake variable ENABLE_CERF = ON
 */
class LensHaloGaussian : public LensHalo{
public:
 
  LensHaloGaussian(float my_mass  
                ,PosType my_zlens 
                ,PosType r_scale  
                ,float my_fratio 
                ,float my_pa     
                ,const COSMOLOGY &cosmo  
                ,float f=100 
  );
  
  LensHaloGaussian(const LensHaloGaussian &h):
  LensHalo(h)
  {
    Rhight = h.Rhight;
    q = h.q;
    q_prime = h.q_prime;
    SigmaO = h.SigmaO;
    pa = h.pa;
    R = h.R;
    
    norm = h.norm;
    norm_g = h.norm_g;
    ss = h.ss;
    I = h.I;
    I_sqpi = h.I_sqpi;
    one_sqpi = h.one_sqpi;
  }
  
  LensHaloGaussian &operator=(const LensHaloGaussian &h){
    if(&h == this) return *this;
    LensHalo::operator=(h);
    
    Rhight = h.Rhight;
    q = h.q;
    q_prime = h.q_prime;
    SigmaO = h.SigmaO;
    pa = h.pa;
    R = h.R;
    
    norm = h.norm;
    norm_g = h.norm_g;
    ss = h.ss;
    I = h.I;
    I_sqpi = h.I_sqpi;
    one_sqpi = h.one_sqpi;
    
    return *this;
  }
  
  ~LensHaloGaussian(){};
 

  void force_halo(PosType *alpha,KappaType *kappa,KappaType *gamma,KappaType *phi,PosType const *xcm,bool subtract_point=false,PosType screening = 1.0);
  
  float get_fratio(){return q;};
  float get_pa(){return pa;};
  float get_scalehight(){return Rhight;};

  float get_SigmaCentral(){return SigmaO;};
 
  void set_pa(double p){pa = p;}
  
  void deflection(std::complex<double> &z
                  ,std::complex<double> &a
                  ,std::complex<double> &g
                  ,KappaType &sigma) const;
  
  
protected:
  
  std::complex<double> dwdz(std::complex<double> z) const{
    return 2.0*(I_sqpi - z*wF(z));
  };
  
  std::complex<double> wF(std::complex<double> z) const;
  std::complex<double> my_erfc(std::complex<double> z) const;
 
  // https://arxiv.org/pdf/1407.0748.pdf
//  std::complex<double> erfc(std::complex<double> z) const{
//    if(z.real() >= 0){
//    std::complex<double> ans = wf(I*z);
//    ans *= exp(-z*z);
    //return ans;
//      return exp(-z*z)*wf(I*z);
//    }else{
//      return 2. - exp(-z*z)*wf(-I*z);
//    }
//  }
  
  double Rhight;  // scale hight in larges dimension
  double q;
  double q_prime;
  double SigmaO;   // SigmaO - central surface density
  double pa;
  std::complex<double> R; // rotation
  
  double norm;
  double norm_g;
  double ss;
  
  //std::complex<double> F(double r_e,double t,std::complex<double> z) const;
  std::complex<double> erfi(std::complex<double> z) const{
    return I*(1. - exp(z*z)*wF(z));
  }
  
//  std::complex<double> rho(std::complex<double> z) const{
//
//    //std::complex<double> y=(q*q*z.real() + I*z.imag())/ss/sqrt(1-q*q) ;
//    //return -I * ( wF( z/ss/sqrt(q_prime) )
//    //                        -  exp( y*y -z*z/ss/ss/q_prime ) * wF(y));
//
//    return wbar(z/ss/ss/sqrt(q_prime),1) - wbar(z/ss/ss/sqrt(q_prime),q);
//  }
//
//  std::complex<double> wbar(std::complex<double> z, double p) const{
//    double x = z.real();
//    double y = z.imag();
//    std::cout << exp(-z*z) << "," <<
//    exp(-x*x*(1-p*p) - y*y*(1./p/p -1)) << "," << -I* wF(p*x + I*y/p) << std::endl;
//
//    double tmp = exp(-x*x*(1-p*p) - y*y*(1./p/p -1));
//    if(tmp==0) return 0;
//
//    return  -(exp(-z*z)*I*tmp)*wF(p*x + I*y/p);
//  }
 
  // Faddeeva function according to Zaghloul (2017)
  // It oesn't seem right near the real axis for small |z|!
  // not used
//  std::complex<double> wf(std::complex<double> z) const{
//    //std::cout << z << std::endl;
//    double y2 = z.imag()*z.imag();
//    double z2 = std::norm(z);
//    if( z2 >= 3.8e4 ){
//      return I_sqpi / z;
//    }
//    if(z2 >= 256){
//      //std::cout << z << "," << (z*z -0.5) << "," << I * PI << std::endl;
//      return z/(z*z -0.5) * I_sqpi;
//    }
//    else if(z2 >= 62 && z2 < 256){
//      return (z*z - 1.) / (z * (z*z-1.5)) * I_sqpi;
//    }
//    else if(z2 < 62 && z2 >= 30 && y2 >= 1.0e-13){
//      return z*(z*z -2.5) / (z*z*(z*z-3.) + 0.75) * I_sqpi;
//    }
//    else if ( (z2 < 62 && z2 >= 30 && y2 < 1.0e-13) ||  (z2 < 30 && z2 >= 2.5 && y2 < 0.072) ){
//      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) ))))))
//      /(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)))))));
//    }
//
//    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) ))))))
//      /(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)))))));
//  }
  
  std::complex<double> I;
  std::complex<double> I_sqpi;
  double one_sqpi;
//  double U1[6] = {1.320522,35.7668,219.031,1540.787,3321.990,36183.31};
//  double V1[7] = {1.841439,61.57037,364.2191,2186.181,9022.228,24322.84,32066.6};
//
//  double U2[6] = {5.9126262,30.180142,93.15558,181.92853,214.38239,122.60793};
//  double V2[7] = {10.479857,53.992907,170.35400,348.70392,457.33448,352.73063,122.60793};
};
 
#ifdef ENABLE_EIGEN

namespace MultiGauss{
 
struct PROFILE{
  virtual double operator()(double r) = 0;
  virtual double cum(double r) = 0;
};
 
struct POWERLAW : public PROFILE
{
  POWERLAW(
           float g  // power-law index of surfcae density
  ){
    if(g <= -2){
      std::cerr << "power-law index <= 2 give unbounded mass!" << std::endl;
      throw std::runtime_error("bad mass");
    }
    gamma = g;
  }
  double operator()(double r){return pow(r,gamma);}
  double cum(double r){return 2*PI*pow(r,gamma + 2)/(gamma+2);}
  
  void reset(float g){gamma = g;}
  float gamma;
};
 
struct SERSIC : public PROFILE
{
  SERSIC(){reset(0,0);}
  SERSIC(float nn,float Re){
    reset(nn,Re);
  }
  double cum(double r){
    return cum_norm * boost::math::tgamma_lower(2*n, bn*pow(r/re,t) );
  }
  double operator()(double r){ //check that these are consistant
    return exp(-bn*pow(r/re,t) );
  }
  
  void reset(float nn,float Re){
    re = Re;
    n = nn;
    bn = 1.9992*n - 0.3271;  // approximation valid for 0.5 < n < 8
    t = 1.0/n;
    cum_norm = 2*PI*n/pow(bn,2*n)*re*re;
  }
  
private:
  float t;
  float re;
  float n;
  float bn;
  float cum_norm;
};
 
struct NFW: public PROFILE
{
  NFW(float rs):rs(rs){
  }
  double cum(double r){
    double x =r/rs;
    
    double ans=log(x/2);
    if(x==1.0) ans += 1.0;
    if(x>1.0) ans += 2*atan(sqrt((x-1)/(x+1)))/sqrt(x*x-1);
    if(x<1.0) ans += 2*atanh(sqrt((1-x)/(x+1)))/sqrt(1-x*x);
    
    return 2*PI*rs*rs*ans;
  }
  double operator()(double r){ //check that these are consistant
    double x =r/rs;
    
    if(x==1.0){return 1.0/3.0;}
    if(x>1.0) return (1-2*atan(sqrt((x-1)/(x+1)))/sqrt(x*x-1))/(x*x-1);
    return (1-2*atanh(sqrt((1-x)/(x+1)))/sqrt(1-x*x))/(x*x-1);
  }
  
  void reset(double Rs){rs = Rs;}
  
  float rs;
};
}
 
/***
\brief  A class for constructing and approximation to any elliptical profile out of a series of elliptical gaussians.
 
 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
 
 The profile is fit Nradii points logarithmicly distributed between r_min and r_max using Ngaussians Gaussians in that range.
 Typically Nradii ~ 2 * Ngaussians.
 
 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.
 
 This class uses both the libcerf library (https://jugit.fz-juelich.de/mlz/libcerf).
 and the eigen library.  They must be installed and linked using the cmake variable ENABLE_CERF = ON and ENABLE_EIGEN=ON
 */
 
class LensHaloMultiGauss: public LensHalo{
  
private:
  
  int nn; // number of gaussians
  int mm; // number of fit radii
  double q; // axis ratio
  double pa; // position angle
  double mass_norm;
  double r_norm;
  std::vector<double> sigmas;
  std::vector<double> radii;
  std::complex<double> Rotation;
  std::vector<LensHaloGaussian> gaussians;
  std::vector<double> A;
  float rms_error;
  //static int count;
 
public:
  
  LensHaloMultiGauss(
                     double mass_norm  
                     ,double Rnorm       
                     ,MultiGauss::PROFILE &profile  
                     ,int Ngaussians    
                     ,int Nradii    
                     ,PosType r_min
                     ,PosType r_max
                     ,PosType my_zlens 
                     ,float my_fratio 
                     ,float my_pa     
                     ,const COSMOLOGY &cosmo  
                     ,float f=10 
                     ,bool verbose = false
                     );
  
  LensHaloMultiGauss(
                     double my_mass_norm
                     ,double Rnorm
                     ,double my_scale   // radial scale in units of the scale that was used to produce relative_scales
                     ,const std::vector<double> &relative_scales
                     ,const std::vector<double> &relative_masses
                     ,PosType my_zlens 
                     ,float my_fratio 
                     ,float my_pa     
                     ,const COSMOLOGY &cosmo  
                     ,float f=100 
                     ,bool verbose = false
                     );
  
  
 
  LensHaloMultiGauss(LensHaloMultiGauss &&halo);
  LensHaloMultiGauss(const LensHaloMultiGauss &halo);
 
  LensHaloMultiGauss & operator=(const LensHaloMultiGauss &&halo);
  LensHaloMultiGauss & operator=(const LensHaloMultiGauss &halo);
  
  ~LensHaloMultiGauss(){};//--LensHaloMultiGauss::count;}
  
  void set_pa(double my_pa);
  
  /***
    This is a static function that can be used to find the masses and scales for the gaussians which can
   then be fed into the scond constructor with a rescaling.  This avoids having to recalculate them when the same profile is used multiple times.
   
   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.
   
   <p>
   MultiGauss::sersic profile(1,1);
   double mass_out;
   float rms_error;
   std::vector<double> scales,masses;
   LensHaloMultiGauss::calc_masses_scales(profile, 13, 25, 0.01,3,0.5,1,mass_out,scales,masses,rms_error);
   LensHaloMultiGauss halo1(mass,rscale,scales,masses,zl,q,-pa,cosmo);
   LensHaloMultiGauss halo2(mass/2,0.1*rscale,scales,masses,zl2,q2,-pa2,cosmo);
   <\p>
   */
  
  static void calc_masses_scales(MultiGauss::PROFILE &profile  
                                 ,int Ngaussians    
                                 ,int Nradii    
                                 ,PosType r_min
                                 ,PosType r_max
                                 ,double mass_norm
                                 ,double r_norm
                                 ,double &totalmass_out
                                 ,std::vector<double> &scales
                                 ,std::vector<double> &masses
                                 ,float &rms_error
                                 ,bool verbose=false
                                 ){
    // make logorithmic
    scales.resize(Ngaussians);
    {
      double tmp = pow(r_max/r_min,1.0/(Ngaussians-2));
      scales[0] = r_min/5;
      scales[1] = r_min;
      for(int i=2 ; i < Ngaussians ; ++i){
        scales[i] = tmp*scales[i-1];
      }
    }
    
    std::vector<double> radii(Nradii);
    {
      double tmp = pow(r_max/r_min,1.0/(Nradii-1));
      radii[0] = r_min;
      for(int i=1 ; i < Nradii ; ++i){
        radii[i] = tmp*radii[i-1];
      }
    }
    
    Eigen::MatrixXd M(Nradii,Ngaussians);
    
    // first radius is fit to mass within radius
    {
      double cp = profile.cum(radii[0]);
      for(int n=0 ; n < Ngaussians ; ++n){
        M(0,n) = 2*PI*scales[n] * scales[n]
        * ( 1 -exp(-radii[0]*radii[0] / 2 / scales[n] /scales[n]) )
        / cp ;
      }
    }
    
    // all others are fit to relative density
    for(int m=1 ; m < Nradii ; ++m){
      double p = profile(radii[m]);
      for(int n=0 ; n < Ngaussians ; ++n){
        M(m,n) = exp( - radii[m]*radii[m] / 2 / scales[n] /scales[n])
        / p ;
      }
    }
    
    // invert M
    
    masses.resize(Ngaussians);
    {
      std::vector<double> b(Nradii,1);
      Eigen::Map<Eigen::VectorXd> v(b.data(),Nradii);
      Eigen::Map<Eigen::VectorXd> a(masses.data(),Ngaussians);
      
      //M = M.inverse();
      a = M.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(v);
      
      if(verbose) std::cout << "test in LensHaloMultiGauss::calc_masses_scales()" << std::endl << M*a << std::endl;
    }
    
    // tests
    {
      rms_error=0;
      for(int m=0 ; m < 1 ; ++m){
        double surf=0;
        for(int n=0 ; n < Ngaussians ; ++n){
          surf += masses[n] * scales[n] * scales[n]
          * 2*PI*( 1 - exp(-radii[m]*radii[m] / 2 / scales[n] /scales[n]) );
        }
        double tmp = profile.cum(radii[m]);
        if(verbose) std::cout << tmp << " "  << surf << std::endl;
        
        rms_error += (surf - tmp )*(surf - tmp) / tmp / tmp;
      }
      
      for(int m=1 ; m < Nradii ; ++m){
        double surf=0;
        for(int n=0 ; n < Ngaussians ; ++n){
          surf += masses[n] * exp(-radii[m]*radii[m] / 2 / scales[n] /scales[n]) ;
        }
        double tmp = profile(radii[m]);
        if(verbose) std::cout << tmp << " "  << surf << std::endl;
        
        rms_error += (surf - tmp )*(surf - tmp ) / tmp /tmp;
      }
      rms_error = sqrt(rms_error)/Ngaussians;
      if(verbose) std::cout << " MSE error " << rms_error << std::endl;
    }
    
    // find masses of Gaussian components
   {
      double mass_tmp = 0;
      for(int n=0 ; n<Ngaussians ; ++n){
        mass_tmp += scales[n]*scales[n]*masses[n]
        *( 1 - exp(-r_norm*r_norm / 2 / scales[n] /scales[n]) );
      }
      
      {  // convert to mass units of each component
        double tmp = (mass_norm/mass_tmp);
        // rescale so total mass is correct
        for(int n=0 ; n<Ngaussians ; ++n) masses[n] *= tmp*scales[n]*scales[n];
      }
     
   }
    totalmass_out = 0;
    for(int n=0 ; n<Ngaussians ; ++n) totalmass_out += masses[n];
  }
  
  double profile(double r);
  
  double mass_cum(double r);
  
  float error();
  
  void force_halo(PosType *alpha,KappaType *kappa,KappaType *gamma,KappaType *phi,PosType const *xcm,bool subtract_point=false,PosType screening = 1);
  
 
};
 
#endif // eigen
#endif // libcerf
 
/*
 *
 * \brief A class for calculating the deflection, kappa and gamma caused by a halos
 * with truncated Hernquist mass profiles.
 *
 * The profile is \f$ \rho \propto \left( \frac{r}{r_s} \right)^{-1} \left( 1 + \frac{r}{r_s} \right)^{-3} \f$.
 *
 *
 */
 
class LensHaloHernquist: public LensHalo{
public:
 //LensHaloHernquist();
  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);
 //LensHaloHernquist(InputParams& params);
 virtual ~LensHaloHernquist();
  
  PosType ffunction(PosType x) const;
 PosType gfunction(PosType x) const;
 PosType hfunction(PosType x) const;
 PosType g2function(PosType x) const;
  
 /* the below functions alphaHern, kappaHern and gammaHern are obsolete and better to be deleted to avoid confusion
   void alphaHern(PosType *alpha,PosType *x,PosType Rtrunc,PosType mass,PosType r_scale
   ,PosType *center,PosType Sigma_crit);
   KappaType kappaHern(PosType *x,PosType Rtrunc,PosType mass,PosType r_scale
   ,PosType *center,PosType Sigma_crit);
   void gammaHern(KappaType *gamma,PosType *x,PosType Rtrunc,PosType mass,PosType r_scale
   ,PosType *center,PosType Sigma_crit);
   */
 //void initFromFile(float my_mass, long *seed, float vmax, float r_halfmass);
  
  void set_rscale(float my_rscale){rscale=my_rscale; xmax = LensHalo::getRsize()/rscale; gmax = InterpolateFromTable(gtable,xmax);};
  // friend struct Ialpha_func;
  
protected:
 static const long NTABLE;
 static const PosType maxrm;
 static int count;
  
 static PosType *ftable,*gtable,*g2table,*htable,*xtable,*xgtable;
 void make_tables();
 PosType InterpolateFromTable(PosType *table, PosType y) const;
  
 void assignParams(InputParams& params);
  
 // Override internal structure of halos
 inline PosType alpha_h(PosType x) const {
  return -0.25*InterpolateFromTable(gtable,x)/gmax;
 }
 inline KappaType kappa_h(PosType x) const {
  return 0.25*x*x*InterpolateFromTable(ftable,x)/gmax; // 0.5*
 }
 inline KappaType gamma_h(PosType x) const {
  return -0.5*x*x*InterpolateFromTable(g2table,x)/gmax;
 }
 inline KappaType phi_h(PosType x) const {
  return -0.25*InterpolateFromTable(htable,x)/gmax;
  //return -1.0*InterpolateFromTable(htable,x)/gmax;
 }
    inline KappaType phi_int(PosType x) const{
  return -0.25*InterpolateFromTable(xgtable,x)/gmax;
  }
  
  
private:
  PosType gmax;
};
 
/*
 *
 * \brief A class for calculating the deflection, kappa and gamma caused by a collection of halos
 * with truncated Jaffe mass profiles.
 *
 * 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.
 *
 *
 */
 
class LensHaloJaffe: public LensHalo{
public:
 //LensHaloJaffe();
  LensHaloJaffe(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);
 //LensHaloJaffe(InputParams& params);
 virtual ~LensHaloJaffe();
  
  void set_rscale(float my_rscale){rscale=my_rscale; xmax = LensHalo::getRsize()/rscale; gmax = InterpolateFromTable(gtable,xmax);};
  
  PosType ffunction(PosType x) const;
 PosType gfunction(PosType x) const;
 PosType hfunction(PosType x) const;
 PosType g2function(PosType x) const;
  PosType bfunction(PosType x);
  PosType dbfunction(PosType x);
  PosType ddbfunction(PosType x);
  
protected:
  
 static const long NTABLE;
 static const PosType maxrm;
 static int count;
  
 static PosType *ftable,*gtable,*g2table,*htable,*xtable,*xgtable;
 void make_tables();
 PosType InterpolateFromTable(PosType *table, PosType y) const;
  
 void assignParams(InputParams& params);
  
 // Override internal structure of halos
 inline PosType alpha_h(PosType x) const{
  return -0.25*InterpolateFromTable(gtable,x)/gmax;
 }
 inline KappaType kappa_h(PosType x) const {
  return 0.125*x*x*InterpolateFromTable(ftable,x)/gmax;
 }
 inline KappaType gamma_h(PosType x) const {
  //return -0.125*x*x*InterpolateFromTable(g2table,x)/gmax;
  return -0.25*x*x*InterpolateFromTable(g2table,x)/gmax;
 }
 inline KappaType phi_h(PosType x) const {
        return -0.25*InterpolateFromTable(xgtable,x)/gmax;
    }
    inline KappaType phi_int(PosType x) const{
  return -0.25*InterpolateFromTable(xgtable,x)/gmax;
  }
  
private:
  PosType gmax;
  
  // I have temporarily set these functions to 0 to make the code compile, Ben
  //  PosType ffunction(PosType x){throw std::runtime_error("Set to temporary invalid value"); return 0;}
  // PosType gfunction(PosType x){throw std::runtime_error("Set to temporary invalid value"); return 0;}
  // PosType hfunction(PosType x){throw std::runtime_error("Set to temporary invalid value"); return 0;}
  // PosType g2function(PosType x){throw std::runtime_error("Set to temporary invalid value"); return 0;}
  
};
 
 
 

class LensHaloDummy: public LensHalo{
public:
 LensHaloDummy();
  LensHaloDummy(float my_mass,float my_Rsize,PosType my_zlens,float my_rscale,const COSMOLOGY &cosmo);
 //LensHaloDummy(InputParams& params);
 ~LensHaloDummy(){};
 
 void force_halo(PosType *alpha,KappaType *kappa,KappaType *gamma,KappaType *phi,PosType const *xcm,bool subtract_point=false,PosType screening = 1.0);
  // void force_halo(PosType *alpha,KappaType *kappa,KappaType *gamma,PosType *xcm,bool subtract_point=false,PosType screening = 1.0);
  
 void initFromMassFunc(float my_mass, float my_Rsize, float my_rscale, PosType my_slope, long *seed);
  
 
private:
 void assignParams(InputParams& params);
 inline PosType alpha_h(PosType x) const {return  0.;}
 inline KappaType kappa_h(PosType x) const {return  0.;}
 inline KappaType gamma_h(PosType x) const {return  0.;}
};
 
 
typedef LensHalo* LensHaloHndl;
 
//bool LensHaloZcompare(LensHalo *lh1,LensHalo *lh2);//{return (lh1->getZlens() < lh2->getZlens());}
//bool LensHaloZcompare(LensHalo lh1,LensHalo lh2){return (lh1.getZlens() < lh2.getZlens());}
 
#endif /* LENS_HALOS_H_ */