lens.h

/*
 * multiplane.h
 *
 */
 
#ifndef MULTIPLANE_H_
#define MULTIPLANE_H_
 
#include "quadTree.h"
#include "utilities_slsim.h"
#include "planes.h"
#include "geometry.h"
 
#include <map>
 
GLAMER_TEST_USES(LensTest)
 

 
class Lens{
public:
  Lens(long* seed, PosType z_source,CosmoParamSet cosmoset, bool verbose = false);
  //Lens(InputParams& params, long* my_seed, CosmoParamSet cosmoset, bool verbose = false);
  Lens(long* seed, PosType z_source,const COSMOLOGY &cosmo, bool verbose = false);
  //Lens(InputParams& params, long* my_seed, const COSMOLOGY &cosmo, bool verbose = false);
  
 ~Lens();

 bool set;

 int getNplanes() const {return lensing_planes.size();}
  
 PosType getfov() const {return fieldofview;};
 void setfov(PosType fov){fieldofview=fov;};

 void resetFieldNplanes(std::size_t field_Nplanes, bool verbose = false);

 void resetFieldHalos(bool verbose = false);

 void printMultiLens();
 
  
 PosType getZlens() const{
  if(flag_switch_main_halo_on)
   return main_halos[0]->getZlens();
  else{
   ERROR_MESSAGE();
   std::cerr << "error, no main lens present" << std::endl;
   exit(1);
  }
 }

 PosType getAngDistLens() const{
  if(flag_switch_main_halo_on)
   return cosmo.angDist( main_halos[0]->getZlens());
  else{
   ERROR_MESSAGE();
   std::cerr << "error, no main lens present" << std::endl;
   exit(1);
  }
 }
  
  Utilities::Geometry::SphericalPoint<> getCenter() const {return central_point_sphere;}

 void clearMainHalos(bool verbose=false);
  template<typename HaloType>
  void clearMainHalo(bool verbose=false);
 

  //void insertMainHalo(LensHalo* halo,PosType zlens, bool addplanes,bool verbose = false);
  
/*  void insertMainHalo(LensHalo *halo, bool addplanes,bool verbose)
  {
//    LensHaloNFW * halo = new LensHaloNFW(halo_in);
    halo->setCosmology(cosmo);
    main_halos.push_back(halo);
    
    flag_switch_main_halo_on = true;
    
    if(addplanes) addMainHaloToPlane(halo);
    else addMainHaloToNearestPlane(halo);
    
    combinePlanes(verbose);
  }
*/
  
  template <typename T>
  void insertMainHalo(const T &halo_in, bool addplanes,bool verbose=false)
  {
    
    T * halo = new T(halo_in);
    halo->setCosmology(cosmo);
    main_halos.push_back(halo);
    
    flag_switch_main_halo_on = true;
    
    if(addplanes) addMainHaloToPlane(halo);
    else addMainHaloToNearestPlane(halo);
    
    combinePlanes(verbose);
  }
  
  template <typename T>
  void moveinMainHalo(T &halo_in, bool addplanes,bool verbose=false)
  {
    T * halo = new T(std::move(halo_in));
    halo->setCosmology(cosmo);
    main_halos.push_back(halo);
    
    flag_switch_main_halo_on = true;
    
    if(addplanes) addMainHaloToPlane(halo);
    else addMainHaloToNearestPlane(halo);
    
    combinePlanes(verbose);
  }
  
  template <typename T>
  void replaceMainHalo(const T &halo_in,bool addplanes,bool verbose=false)
  {
    main_halos.clear();
    
    T * halo = new T(halo_in);
    halo->setCosmology(cosmo);
    main_halos.push_back(halo);
    
    flag_switch_main_halo_on = true;
    
    Utilities::delete_container(main_planes);
    createMainPlanes();
    combinePlanes(verbose);
  }

  template <typename T>
  void insertMainHalos(std::vector<T> &my_halos,bool addplanes, bool verbose=false)
  {
    T* ptr;
    //for(std::size_t i = 0; i < my_halos.size() ; ++i)
    for(T &h : my_halos){
      ptr = new T(h);
      ptr->setCosmology(cosmo);
      ptr->setDist(cosmo);
      main_halos.push_back( ptr );
      if(addplanes) addMainHaloToPlane( ptr );
      else addMainHaloToNearestPlane( ptr );
    }
    
    flag_switch_main_halo_on = true;
    
    combinePlanes(verbose);
  }

  template <typename T>
  void replaceMainHalos(std::vector<T> &my_halos,bool verbose)
  {
    
    main_halos.clear();
 
    T* ptr;
    //for(std::size_t i = 0; i < my_halos.size() ; ++i)
    for(T &h : my_halos){
      ptr = new T(h);
      ptr->setCosmology(cosmo);
      ptr->setDist(cosmo);
      main_halos.push_back( ptr );
    }
    
    flag_switch_main_halo_on = true;
    
    Utilities::delete_container(main_planes);
    createMainPlanes();
    combinePlanes(verbose);
  }

 //void insertMainHalos(LensHalo** halos, std::size_t Nhalos,bool addplanes,bool verbose = false);

  //void replaceMainHalo(LensHalo* halo,PosType zlens, bool addplanes,bool verbose = false);

 // replaceMainHalos(LensHalo** halos, std::size_t Nhalos,bool verbose = false);
 
 
//  /* \brief Add substructures to the lens.
//
//   This method is meant for inserting substructure to a main lens.  All the substructure will be at
//   one redshift.  The mass function follows a power law.  The density of substructures is constant within
//   a circular region.  The tidal truncation is controlled through the parameter density_contrast which is
//   the average density within the substructures orbit in units of the average density to the universe at
//   the redshift where they are places.  For example density_contrast=200 would give them the truncation radius appropriate at R_200.
//   */
//  void insertSubstructures(
//        PosType Rregion            /// radius of region in which substructures are inserted (radians)
//        ,PosType center[]          /// center of region in which the substructures are inserted (radians)
//        ,PosType NumberDensity     /// number density per radian^2 of all substructures
//        ,PosType Mass_min          /// minimum mass of substructures
//        ,PosType Mass_max          /// maximum mass of substructures
//        ,PosType redshift          /// redshift of substructures
//        ,PosType alpha             /// index of mass function (dN/dm \propto m^alpha)
//        ,PosType density_contrast  ///
//        ,bool verbose
//  );
//
//  /** \brief This function will randomize the substructure without changing the region, mass function, etc.
//
//   The Lens::insertSubstructures() function must have been called on this instance of the Lens before.
//   */
//  void resetSubstructure(bool verbose = false);
//

 std::size_t getNMainHalos() const;
 template<typename HaloType>
 std::size_t getNMainHalos() const;
 
 LensHalo* getMainHalo(std::size_t i);
  
 template<typename HaloType>
 HaloType* getMainHalo(std::size_t i);

  void rayshooter(RAY &ray);

 void rayshooterInternal(unsigned long Npoints    
                          ,Point *i_points         
                          ,bool RSIverbose = false
                          );
  
  
  void rayshooterInternal(unsigned long Npoints   
                          ,LinkedPoint *i_points        
                          ,std::vector<double> &source_zs 
                          ,bool RSIverbose = false
                          );
  
  void rayshooterInternal(unsigned long Npoints   
                         ,RAY *rays            
                         );
  void rayshooterInternal(RAY &ray);  
 
  void info_rayshooter(RAY &i_point
                      ,std::vector<Point_2d> & ang_positions
                      ,std::vector<KappaType> & kappa_on_planes
                      ,std::vector<std::vector<LensHalo*>> & halo_neighbors
                      ,LensHalo &halo_max
                      ,KappaType &kappa_max
                      ,KappaType gamma_max[]
                      ,PosType rmax  
                      ,int tag = 0   
                      ,short mode = 0  
                      ,bool verbose = false
                                       );
  
  void mass_on_planes(const std::vector<RAY> &rays     
                    ,std::vector<double> &masses     
                    ,bool verbose = false
                    );
  
  /*RAY find_image_min(
          Point_2d y_source    /// input position of source (radians)
          ,Point_2d &x_image    /// initial guess for image postion (radians)
          ,PosType z_source     /// redshift of source
          ,PosType ytol2        /// target tolerance in source position squared
          ,PosType &fret        ///
          ,int sign=0             /// sign of magnification, it is found automatically if left out
  );*/
  
 
  RAY find_image_min(
            Point &p              
            ,double zs            
            ,PosType ytol2        
            ,PosType &dy2         
            ,bool use_image_guess 
  );
  RAY find_image_min(const RAY &in_ray    
                     ,PosType ytol2       
  );
  
  RAY find_image_min(
            Point &p              
            ,double zs            
            ,PosType ytol2        
            ,PosType &dy2         
            ,std::vector<Point_2d> &boundary   
  );
  void find_images_min_parallel(std::vector<RAY> &rays
                                        ,double ytol2
                                        ,std::vector<bool> &success
                                      );
  
  std::vector<RAY> find_images(Point_2d y_source
                                    ,double z_source
                                    ,Point_2d center
                                    ,double range  
                                    ,double stop_res
                                    );
  std::vector<RAY> find_images(GridMap &init_grid
                              ,Point_2d y_source
                              ,double z_source
                              ,double stop_res
                              );
  
 
//  void find_point_source_images(
//                              Grid &grid
//                              ,Point_2d y_source
//                              ,PosType r_source
//                              ,PosType z_source
//                              ,std::vector<RAY> &images
//                              ,double dytol2
//                              ,double dxtol
//                              ,bool verbose = false
//                                 );
  
  
 // methods used for use with implanted sources

 short ResetSourcePlane(PosType z,bool nearest=false,bool verbose = false);

  void FindSourcePlane(
                          PosType zs                 
                          ,long &jmax                
                          ,double &Dls               
                          ,double &Ds                
  );

 void RevertSourcePlane(){ toggle_source_plane = false;}
 //void ImplantSource(unsigned long index,CosmoHndl cosmo);
 PosType getSourceZ(){
  if(toggle_source_plane){
   return zs_implant;
  }else{
   return plane_redshifts.back();
  }
 }
 
 PosType getZmax() const{return plane_redshifts.back();}

 void PrintCosmology() { cosmo.PrintCosmology(); }
 
 PosType getSigmaCrit(PosType zsource) const{ return cosmo.SigmaCrit(getZlens(), zsource); }
 

  const COSMOLOGY & getCosmo(){return cosmo;}
  
  void TurnFieldOff() { flag_switch_field_off = true ; }
  void TurnFieldOn() { flag_switch_field_off = false ; }
  
  PosType getFieldMinMass() const { return field_min_mass ; }
 
  // get the field_Off value :
  bool getfieldOff() const {return flag_switch_field_off ;}
  
  void GenerateFieldHalos(double min_mass 
                          ,MassFuncType mass_function 
                          ,double field_of_view  
                          ,int Nplanes           
                          ,LensHaloType halo_type = LensHaloType::nfw_lens  
                          ,GalaxyLensHaloType galaxy_type = GalaxyLensHaloType::null_gal  
                          ,double buffer = 1.0 
                          ,bool verbose = false
                );
  
  Lens & operator=(Lens &&lens){
    
    fieldofview = lens.fieldofview;
    seed = lens.seed;
    init_seed = lens.init_seed;
    cosmo = lens.cosmo;
    toggle_source_plane = lens.toggle_source_plane;
  
    zs_implant = lens.zs_implant;
    zsource = lens.zsource;
 
    NZSamples = lens.NZSamples;
    zbins = lens.zbins;
    NhalosbinZ = lens.NhalosbinZ;
    Nhaloestot_Tab = lens.Nhaloestot_Tab;
    aveNhalosField = lens.aveNhalosField;
    Logm = lens.Logm;
    NhalosbinMass = lens.NhalosbinMass;
    sigma_back_Tab = lens.sigma_back_Tab;
    
    flag_switch_deflection_off = lens.flag_switch_deflection_off;
    flag_switch_lensing_off = lens.flag_switch_lensing_off;
    
    Dl = lens.Dl;
    dDl = lens.dDl;
    dTl = lens.dTl;
    plane_redshifts = lens.plane_redshifts;
    charge = lens.charge;
    flag_switch_field_off = lens.flag_switch_field_off;
    
    field_halos = lens.field_halos;
 
    field_Nplanes_original = lens.field_Nplanes_original;
    field_Nplanes_current = lens.field_Nplanes_current;
    
    field_plane_redshifts = lens.field_plane_redshifts;
    field_plane_redshifts_original = lens.field_plane_redshifts_original;
    field_Dl = lens.field_Dl;
    field_Dl_original = lens.field_Dl_original;
    
    //substructure = lens.substructure;
    
    //WasInsertSubStructuresCalled = lens.WasInsertSubStructuresCalled;
    field_mass_func_type = lens.field_mass_func_type;
    mass_func_PL_slope = lens.mass_func_PL_slope;
    field_min_mass = lens.field_min_mass;
    field_int_prof_type = lens.field_int_prof_type;
    field_prof_internal_slope = lens.field_prof_internal_slope;
    
    flag_field_gal_on = lens.flag_field_gal_on;
    field_int_prof_gal_type = lens.field_int_prof_gal_type;
    field_int_prof_gal_slope = lens.field_int_prof_gal_slope;
    
    redshift_planes_file = lens.redshift_planes_file;
    read_redshift_planes = lens.read_redshift_planes;
    
    field_input_sim_file = lens.field_input_sim_file;
    field_input_sim_format = lens.field_input_sim_format;
    
    sim_input_flag = lens.sim_input_flag;
    read_sim_file = lens.read_sim_file;
    field_buffer = lens.field_buffer;
    
    flag_switch_main_halo_on = lens.flag_switch_main_halo_on;
    
    main_plane_redshifts = lens.main_plane_redshifts;
    main_Dl = lens.main_Dl;
    
    main_halo_type = lens.main_halo_type;
    main_galaxy_halo_type = lens.main_galaxy_halo_type;
 
    pixel_map_input_file = lens.pixel_map_input_file;
    pixel_map_on = lens.pixel_map_on;
    pixel_map_zeropad = lens.pixel_map_zeropad;
    pixel_map_zeromean = lens.pixel_map_zeromean;
    
    central_point_sphere = lens.central_point_sphere;
    sim_angular_radius = lens.sim_angular_radius;
    inv_ang_screening_scale = lens.inv_ang_screening_scale;
    
    std::swap(lensing_planes,lens.lensing_planes);
    std::swap(field_planes,lens.field_planes);
    swap(main_halos,lens.main_halos);  
    std::swap(main_planes,lens.main_planes);
    
    return *this;
  }
  Lens(Lens &&lens){
    *this = std::move(lens);
  }
  
  // get the adress of field_plane_redshifts
  std::vector<PosType>  get_plane_redshifts () { return plane_redshifts; }
 
private:
  Lens & operator=(const Lens &lens);   // block copy
  Lens(const Lens &lens);
  
protected:
  PosType fieldofview;
 
private:
 GLAMER_TEST_FRIEND(LensTest)
  
  
  void _find_images_(std::vector<RAY> &images
                     ,Point_2d &y
                     ,Point_2d ll
                     ,Point_2d ur
                     ,double resolution
                     ,double stop_res
                     );
  
  void _find_images_min_parallel_(RAY *rays
                                  ,size_t begin
                                  ,size_t end
                                  ,double ytol2
                                  ,std::vector<bool> &success
                                  );
  
 // seed for random field generation
 long *seed;
  
  long init_seed;
  //InputParams init_params;
  
  // the cosmology
 COSMOLOGY cosmo;
 
  template<typename P>
  void compute_points_parallel(int start
                             ,int chunk_size
                             ,P *i_point
                             ,double *source_zs
                             ,bool multiZs
                             ,bool verbose = false);
  
  // this is a version that uses RAYs instead of Points or LinkedPoints and must be kept sinked with the above
  void compute_rays_parallel(int start
                                   ,int chunk_size
                                   ,RAY *rays
                                   ,bool verbose = false);
 
 //void readCosmology(InputParams& params);
 //void assignParams(InputParams& params,bool verbose = false);
  void defaultParams(PosType zsource,bool verbose = true);
 
 bool toggle_source_plane;
 
  // redhsift of true source plane if rest from original
  PosType zs_implant;
 
 PosType zsource;
 
 void quicksort(LensHaloHndl *halo,PosType **pos,unsigned long N);
 
 // create the lens planes
 //void buildPlanes(InputParams& params, bool verbose);
 
 void setFieldDist();
 void setFieldDistFromFile();
 void setupFieldPlanes();
  void ComputeHalosDistributionVariables ();
  
  enum DM_Light_Division {All_DM,Moster};
 
 void createFieldHalos(bool verbose,DM_Light_Division division = Moster);
  
 void readInputSimFileMillennium(bool verbose,DM_Light_Division division = Moster);
  
 void readInputSimFileMultiDarkHalos(bool verbose,DM_Light_Division division = Moster);
  
  void readInputSimFileObservedGalaxies(bool verbose);
  
 void createFieldPlanes(bool verbose);
 
 // generate main halo from the parameter file
 //void createMainHalos(InputParams& params);
 void createMainPlanes();
 void addMainHaloToPlane(LensHalo* halo);
  void addMainHaloToNearestPlane(LensHalo* halo);

 void combinePlanes(bool verbose);
 
  /* Variables used by buildPlanes, createFieldHalos, and createFieldPlanes */
  
  const int Nzbins = 64 ;
  const int Nmassbin=64;
  int NZSamples = 50;
  // table for redshift bins for mass function
  std::vector<PosType> zbins ;
  // table of number of halos per redshift bins for mass function
  std::vector<PosType> NhalosbinZ ;
  std::vector<PosType> Nhaloestot_Tab ;
  PosType aveNhalosField ;
  std::vector<PosType> Logm;
  std::vector<std::vector<PosType>> NhalosbinMass;
  std::vector<PosType> sigma_back_Tab;
  
  /* ----- */
  
  size_t getNFieldHalos() const {return field_halos.size();}
  //size_t getNSubHalos() const {return substructure.halos.size();}
  
private: /* force calculation */
 bool flag_switch_deflection_off;
 bool flag_switch_lensing_off;
 
 std::vector<LensPlane *> lensing_planes;
 std::vector<PosType> Dl;
  std::vector<PosType> dDl;
  std::vector<PosType> dTl;
 std::vector<PosType> plane_redshifts;
 PosType charge;
 
private: /* field */
 bool flag_switch_field_off;
 
  std::vector<LensHalo *> field_halos;
  std::size_t field_Nplanes_original;
  std::size_t field_Nplanes_current;
  
 std::vector<LensPlane*> field_planes;
  std::vector<PosType> field_plane_redshifts;
  std::vector<PosType> field_plane_redshifts_original;
  std::vector<PosType> field_Dl;
  std::vector<PosType> field_Dl_original;
  
//  struct SubStructureInfo{
//    // things for substructures
//    /// vector of all substructure halos
//    std::vector<LensHaloNFW *> halos;
//    LensPlane *plane;
//    PosType Rregion = 0;
//    PosType Mmax = 0;
//    PosType Mmin = 0;
//    PosType alpha = 0;
//    PosType Ndensity = 0;
//    Point_2d center;
//    PosType rho_tidal = 0;
//    // Added quantities for the resetting of the substructure
//    // (when WasInsertSubStructuresCalled = MAYBE) :
//    PosType redshift = 0;
//    bool verbose = false;
//  };
//
//  SubStructureInfo substructure;
  
 // Boo WasInsertSubStructuresCalled = NO ;
  
 //PosType **halo_pos;
 
 MassFuncType field_mass_func_type;
 PosType mass_func_PL_slope;
 PosType field_min_mass;
 LensHaloType field_int_prof_type;
 PosType field_prof_internal_slope;
 
 bool flag_field_gal_on;
 GalaxyLensHaloType field_int_prof_gal_type;
  PosType field_int_prof_gal_slope;
 
 // mass fraction in the host galaxy
 //PosType field_galaxy_mass_fraction;
 
 std::string redshift_planes_file;
 bool read_redshift_planes;
 
 std::string field_input_sim_file;
  HaloCatFormats field_input_sim_format;
  
 bool sim_input_flag;
 //std::string input_gal_file;
 //bool gal_input_flag;
 bool read_sim_file;
 
 PosType field_buffer;
 
private: /* main */
 bool flag_switch_main_halo_on;
 
 Utilities::MixedVector<LensHalo*> main_halos;
 std::vector<LensPlane*> main_planes;
 std::vector<PosType> main_plane_redshifts;
 std::vector<PosType> main_Dl;
 
 LensHaloType main_halo_type;
 GalaxyLensHaloType main_galaxy_halo_type;
 
private: /* input */
 std::string pixel_map_input_file;
 
  short pixel_map_on;
  int pixel_map_zeropad;
  bool pixel_map_zeromean;
 void readPixelizedDensity();
  
  Utilities::Geometry::SphericalPoint<PosType> central_point_sphere;
  PosType sim_angular_radius;
  PosType inv_ang_screening_scale;
};
 
inline std::size_t Lens::getNMainHalos() const
{
 return main_halos.size();
}
 
template<typename HaloType>
inline std::size_t Lens::getNMainHalos() const
{
 return main_halos.size<HaloType>();
}
 
inline LensHalo* Lens::getMainHalo(std::size_t i)
{
 return main_halos.at(i);
}
 
template<typename HaloType>
inline HaloType* Lens::getMainHalo(std::size_t i)
{
  if(main_halos.size<HaloType>() == 0 ) return nullptr;
 return main_halos.at<HaloType>(i);
}
 
typedef Lens* LensHndl;
 
// P should be Point or LinkedPoint
template <typename P>
void Lens::compute_points_parallel(int start
                                 ,int chunk_size
                                 ,P *i_points
                                 ,double *source_zs
                                 ,bool multiZs
                                 ,bool verbose)
{
  int end        = start + chunk_size;
 
  int i, j;
  
  PosType xx[2];
  PosType aa,bb;
  PosType alpha[2];
  
  KappaType kappa,gamma[3];
  KappaType phi;
  
  Matrix2x2<PosType> G;
 
  PosType SumPrevAlphas[2];
  Matrix2x2<PosType> SumPrevAG;
  
  PosType *theta;
  
  long jmax = lensing_planes.size();
  double Dls_Ds; // this is the ratio between of the distance between the last lens plane and the source to the distance to the source
  double D_Ds; // this is the ratio between of the distance to the last lens plane and the source to the distance to the source
 
  if(!multiZs){
    if(source_zs[0] == plane_redshifts.back() ){
      Dls_Ds = dDl.back() / Dl.back();
      D_Ds = Dl[Dl.size() - 2] / Dl.back();
    }else{
      PosType Dls,Ds;
      FindSourcePlane(source_zs[0],jmax,Dls,Ds);
      Dls_Ds = Dls / Ds;
      if(jmax > 0) D_Ds = Dl[jmax-1] / Ds;
    }
  }
  
  // Main loop : loop over the points of the image
  for(i = start; i < end; i++)
  {
    // In case e.g. a temporary point is outside of the grid.
    if(i_points[i].in_image == MAYBE) continue;
    
    //theta = i_points[i].image->x;
    theta = i_points[i].ptr_y();
    
    theta[0] = i_points[i].x[0];
    theta[1] = i_points[i].x[1];
 
    // Initializing SumPrevAlphas :
    SumPrevAlphas[0] = theta[0];
    SumPrevAlphas[1] = theta[1];
 
    // Initializing SumPrevAG :
    SumPrevAG.setToI();
    
    // Setting phi on the first plane.
    phi = 0.0;
    
    // Default values :
    i_points[i].A.setToI();
    i_points[i].dt = 0;
    
    // In case we don't want to compute the values :
    if(flag_switch_lensing_off)
    {
      i_points[i].image->A.setToI();
      continue;
    }
    
    // Time delay at first plane : position on the observer plane is (0,0) => no need to take difference of positions.
    i_points[i].dt = 0;
    
    //0.5*( p->i_points[i].image->x[0]*p->i_points[i].image->x[0] + p->i_points[i].image->x[1]*p->i_points[i].image->x[1] )/ p->dDl[0] ;
    
    if(multiZs){
      PosType Dls,Ds;
      FindSourcePlane(source_zs[i],jmax,Dls,Ds);
      Dls_Ds = Dls / Ds;
      if(jmax > 0) D_Ds = Dl[jmax-1]/Ds;
    }
    
    // Begining of the loop through the planes :
    // Each iteration leaves i_point[i].image on plane (j+1)
 
    for(j = 0; j < jmax ; ++j)
      {
      
      double Dphysical = Dl[j]/(1 + plane_redshifts[j]);
      // convert to physical coordinates on the plane j, just for force calculation
      xx[0] = theta[0] *  Dphysical;
      xx[1] = theta[1] *  Dphysical;
      // PhysMpc = ComMpc / (1+z)
      
      assert(xx[0] == xx[0] && xx[1] == xx[1]);
      
      
      lensing_planes[j]->force(alpha,&kappa,gamma,&phi,xx);
      // Computed in physical coordinates, xx is in PhysMpc.
      
      
      // if(alpha[0] != alpha[0] || alpha[1] != alpha[1]){
      //   std::cout << "xx[0] = " << xx[0] << " xx[1] = " << xx[1] << std::endl ;
      //   std::cout << "alpha = " << alpha[0] << " " << alpha[1] << std::endl ;
      //   throw std::runtime_error("alpha is NaN");
      // }
      assert(alpha[0] == alpha[0] && alpha[1] == alpha[1]);
      assert(gamma[0] == gamma[0] && gamma[1] == gamma[1]);
      assert(kappa == kappa);
      if(std::isinf(kappa)) { std::cout << "xx = " << xx[0] << " " << xx[1] << std::endl ;}
      assert(!std::isinf(kappa));
      
      G[0] = kappa + gamma[0];    G[1] = gamma[1];
      G[2] = gamma[1]; G[3] = kappa - gamma[0];
  
      /* multiply by fac to obtain 1/comoving_distance/physical_distance
       * such that a multiplication with the charge (in units of physical distance)
       * will result in a 1/comoving_distance quantity */
      
      G *= charge * Dl[j] / (1 + plane_redshifts[j]);
        
      assert(gamma[0] == gamma[0] && gamma[1] == gamma[1]);
      assert(kappa == kappa);
      //assert(phi == phi);
      
      // This computes \vec{x}^{j+1} in terms of \vec{x}^{j}
      // according to the corrected Eq. (18) of paper GLAMER II ---------------------------------
      
      // Adding the j-plane alpha contribution to the sum \Sum_{k=1}^{j} \vec{alpha_j} :
      SumPrevAlphas[0] -= charge * alpha[0] ;
      SumPrevAlphas[1] -= charge * alpha[1] ;
      
      if(j < jmax-1 ){
        aa = dDl[j+1] / Dl[j+1];
        bb = Dl[j] / Dl[j+1];
      }else{
        aa = Dls_Ds;
        bb = D_Ds;
      }
      
        if(!flag_switch_deflection_off){
          theta[0] = bb * theta[0] + aa * SumPrevAlphas[0];
          theta[1] = bb * theta[1] + aa * SumPrevAlphas[1];
        }
   
      // ----------------------------------------------------------------------------------------
            
      // Sum_{k=1}^{j} Dl[k] A^k.G^k
      SumPrevAG -= (G * (i_points[i].A)) ;
      
      // Computation of the "plus quantities", i.e. the  next plane quantities :
      i_points[i].A = i_points[i].A * bb + SumPrevAG * aa;
      
      // ----------------------------------------------
      
      // Geometric time delay with added potential
      //p->i_points[i].dt += 0.5*( (xplus[0] - xminus[0])*(xplus[0] - xminus[0]) + (xplus[1] - xminus[1])*(xplus[1] - xminus[1]) ) * p->dTl[j+1] /p->dDl[j+1] /p->dDl[j+1] - phi * p->charge ; /// in Mpc  ???
      
      // Check that the 1+z factor must indeed be there (because the x positions have been rescaled, so it may be different compared to the draft).
      // Remark : Here the true lensing potential is not "phi" but "phi * p->charge = phi * 4 pi G".
      
      
    } // End of the loop going through the planes
    
    if(flag_switch_deflection_off){
      i_points[i].A = Matrix2x2<KappaType>::I() - SumPrevAG;
    }
    
    // Subtracting off a term that makes the unperturbed ray to have zero time delay
    //p->i_points[i].dt -= 0.5*( p->i_points[i].image->x[0]*p->i_points[i].image->x[0] + p->i_points[i].image->x[1]*p->i_points[i].image->x[1] ) / p->Dl[NLastPlane];
    
    // Conversion of dt from Mpc (physical Mpc) to Years -----------------
    i_points[i].dt *= MpcToSeconds * SecondToYears ;
    
    // ---------------------------------------------------------------------------------------------
    
    // Putting the final values of the quantities in the image point -----
    i_points[i].image->A = i_points[i].A;
    i_points[i].image->dt = i_points[i].dt;
    // ------------------------------------------------------------------------
    
    
    // TEST : showing final quantities
    // ------------------------------=
    if(verbose)  std::cout << "RSI final : X X | " << i << "  " << source_zs[i] << " | " << i_points[i].kappa() << " " << i_points[i].gamma1() << " " << i_points[i].gamma2() << " " << i_points[i].gamma3() << " " << i_points[i].invmag() << " | " << i_points[i].dt << std::endl ;
    
  } // End of the main loop.
 
}
 
#endif /* MULTIPLANE_H_ */