Doc/html/gridmap_8h_source.html

//

//  gridmap.h

//  GLAMER

//

//  Created by bmetcalf on 8/19/14.

//

//


#ifndef __GLAMER__gridmap__

#define __GLAMER__gridmap__


#include <iostream>

#include <memory>

#include "lens.h"

#include "point.h"

#include "concave_hull.h"

#include "Tree.h"

#include "source.h"

#include <mutex>


struct GridMap{


  friend class Lens;


 GridMap(LensHndl lens,unsigned long N1d,const double center[2],double range);

  GridMap(LensHndl lens ,unsigned long Nx ,const PosType center[2] ,PosType rangeX ,PosType rangeY);

  GridMap(unsigned long N1d,const double center[2],double range);

 ~GridMap();


  GridMap ReInitialize(LensHndl lens);


  void deLens();

  double RefreshSurfaceBrightnesses(Source* source);

  double AdaptiveRefreshSurfaceBrightnesses(Lens &lens,Source &source);


  double AddSurfaceBrightnesses(Source* source);


  Point_2d image_point(size_t index){return i_points[index];}

  Point_2d source_point(size_t index){return s_points[index];}

  RAY ray(size_t index){return i_points[index];}


  void ClearSurfaceBrightnesses();

  void assertNAN(); // check for nan in surface prightness

 size_t getNumberOfPoints() const {return Ngrid_init*Ngrid_init2;}


 int getInitNgrid() const {return Ngrid_init;}

 double getXRange() const {return x_range;}

 double getYRange() const {return x_range*axisratio;}

  double getResolution() const {return x_range/(Ngrid_init-1);}


  template<typename T>

  PixelMap<T> writePixelMap(LensingVariable lensvar);

  template <typename T>

  void writeFits(LensingVariable lensvar,std::string filensame);


  template<typename T>

  void writePixelMapUniform(PixelMap<T> &map,LensingVariable lensvar);

  template <typename T>

  void writeFitsUniform(const PosType center[],size_t Nx,size_t Ny,LensingVariable lensvar,std::string filename);

  template<typename T>

  PixelMap<T> writePixelMapUniform(const PosType center[],size_t Nx,size_t Ny,LensingVariable lensvar);


  template<typename T>


  void writeFitsUniform(

                        LensingVariable lensvar

                        ,std::string filename

                        ){

    PixelMap<T> map = writePixelMap<T>(lensvar);

    map.printFITS(filename);

  }


  template<typename T>

  PixelMap<T> getPixelMapFlux() const;


  template<typename T>

  void getPixelMapFlux(PixelMap<T> &map) const;


  PosType EinsteinArea() const;


  //PosType magnification_invmag() const;

  //PosType magnification2() const;


  Point_2d centroid() const;


  Point_2d getCenter(){return center;}


  Point * operator[](size_t i){return i_points.data() + i;};


  GridMap(GridMap &&grid){


    Ngrid_init = grid.Ngrid_init;

    Ngrid_init2 = grid.Ngrid_init2;

    pointID = grid.pointID;

    axisratio = grid.axisratio;

    x_range = grid.x_range;


    std::swap(i_points,grid.i_points);

    std::swap(s_points,grid.s_points);


    center = grid.center;

  }


  GridMap(GridMap &grid){


    Ngrid_init = grid.Ngrid_init;

    Ngrid_init2 = grid.Ngrid_init2;

    pointID = grid.pointID;

    axisratio = grid.axisratio;

    x_range = grid.x_range;


    i_points = grid.i_points;

    s_points = grid.s_points;


    center = grid.center;

  }


  GridMap & operator=(GridMap &&grid){

    Ngrid_init = grid.Ngrid_init;

    Ngrid_init2 = grid.Ngrid_init2;

    pointID = grid.pointID;

    axisratio = grid.axisratio;

    x_range = grid.x_range;


    std::swap(i_points,grid.i_points);

    std::swap(s_points,grid.s_points);


    center = grid.center;


    return *this;

  }


  struct Triangle{

    Triangle(size_t i,size_t j,size_t k){

      index[0] = i;

      index[1] = j;

      index[2] = k;

    }

    Triangle(){

      index[0] = 0;

      index[1] = 0;

      index[2] = 0;

    }

    Triangle(const Triangle &tri){

      index[0] = tri[0];

      index[1] = tri[1];

      index[2] = tri[2];

    }

    Triangle & operator=(const Triangle &tri){

      index[0] = tri[0];

      index[1] = tri[1];

      index[2] = tri[2];

      return *this;

    }

    ~Triangle(){};

    std::vector<long> index = {0,0,0};

    long & operator[](int i){return index[i];}

    long operator[](int i) const {return index[i];}

  };

  struct Rectangle{

    Rectangle(size_t i,size_t j){

      index[0] = i;

      index[1] = j;

    }

    Rectangle(){

      index[0] = 0;

      index[1] = 0;

    }

    Rectangle(const Rectangle &rec){

      index[0] = rec[0];

      index[1] = rec[1];

    }

    Rectangle & operator=(Rectangle &rec){

      index[0] = rec[0];

      index[1] = rec[1];

      return *this;

    }

    ~Rectangle(){};

    std::vector<long> index = {0,0};

    long & operator[](int i){return index[i];}

    long operator[](int i) const {return index[i];}

  };


  // determines if two rectangles touch

  inline bool touch(const Rectangle &tr1,const Rectangle &tr2) const;


  // returns a vector of rectangles that encompase all triangles with one rectangle for each touching group

  std::vector<Rectangle> merge_boxes(

                   std::list<Triangle> &triangles

                   ) const;

  std::vector<Rectangle> merge_boxes(

                   std::vector<Triangle> &triangles

                   ) const;


  std::list<RAY> find_images(std::vector<Point_2d> &ys

                             ,std::vector<int> &multiplicity

                             ) const;


  void find_images(Point_2d y

                   ,std::vector<Point_2d> &image_points

                   ,std::vector<Triangle> &triangles

  ) const ;


  void find_boundaries_of_caustics(std::vector<std::vector<Point_2d> > &boundaries

                           ,std::vector<bool> &hits_edge

                           ){


    size_t N=Ngrid_init*Ngrid_init2;

    std::vector<Point_2d> source_points(N);

    for(size_t i=0 ; i<N ; ++i) source_points[i] = i_points[i];


    std::vector<int> multiplicities(N);

    find_images(source_points,multiplicities);


    std::vector<bool> bitmap(N,false);

    for(size_t i=0 ; i<N ; ++i){

      if(multiplicities[i] > 1) bitmap[i] = true;

    }


    Utilities::find_boundaries<Point_2d>(bitmap,Ngrid_init,boundaries,hits_edge);

  }


  PosType magnificationFlux(Source &source) const ;


  double magnificationTr() const ;


  double magnificationTr(std::vector<size_t> &pixels) const ;


  double AreaCellOnSourcePlane(size_t k) const;


// depricated version of PixalMap::magnificationTr()

//  double magnificationTr2() const {

//

//    size_t k;

//    size_t k1;

//    size_t k2;

//

//    double flux_source=0.0,flux_image=0.0;

//    for(size_t i=1 ; i< Ngrid_init-1 ; i++){

//      for(size_t j=1 ; j< Ngrid_init2-1 ; j++){

//

//        k = i + j*Ngrid_init;

//

//        assert( s_points[k].surface_brightness == i_points[k].surface_brightness);

//        double sb = s_points[k].surface_brightness;

//

//        if(sb !=0.0){

//          flux_source += AreaCellOnSourcePlane(k) *sb;

//          flux_source += AreaCellOnSourcePlane(k-1) *sb;

//          flux_source += AreaCellOnSourcePlane(k-1-Ngrid_init) *sb;

//          flux_source += AreaCellOnSourcePlane(k-Ngrid_init) *sb;

//        }

//

//        flux_image += sb;

//      }

//    }

//

//    return flux_image * getResolution() * getResolution() * 4 / flux_source;

//  }


  double AddPointSource(const Point_2d &y,double flux);


  void find_crit(std::vector<std::vector<Point_2d> > &points

                 ,std::vector<bool> &hits_boundary

                 ,std::vector<CritType> &crit_type

                 );


  void find_magnification_contour(

      std::vector<std::vector<Point_2d> > &curves

      ,std::vector<bool> &hits_boundary

      ,double invmag

  );


//  void find_crit_boundary(std::vector<std::vector<Point_2d> > &points

//                          ,std::vector<bool> &hits_boundary

//                          ) const;

  int getNx(){return Ngrid_init;}

  int getNy(){return Ngrid_init2;}

private:

  GridMap & operator=(GridMap &grid);


  // curve must be in pixel units

  bool  incurve(long k,std::vector<Point_2d> &curve) const;


  // cluge to make compatible with old method of producing points

  std::vector<Point> NewPointArray(size_t N){

    std::vector<Point> p(N);

    p[0].head = N;

    return p;

  }


  void xygridpoints(Point *points,double range,const double *center,long Ngrid

                    ,short remove_center);


 int Ngrid_init;

  int Ngrid_init2;


  unsigned long pointID;

  PosType axisratio;

  PosType x_range;

  template<typename T>

  void writePixelMapUniform_(Point* points,size_t size,PixelMap<T> *map,LensingVariable val);


  std::vector<Point> i_points;

  std::vector<Point> s_points;

  Point_2d center;


  bool to_refine(long i,long j,double total,double f) const ;

  static std::mutex grid_mutex;


  void _find_images_(Point_2d *ys,int *multiplicity,long Nys,std::list<RAY> &rays) const;

  void _find_images2_(size_t j1

                              ,size_t j2

                              ,std::list<Point_2d> &image_points

                              ,std::list<Triangle> &triangles

                              ,Point_2d y

                              ) const;


  // find if there are images of y in specific cells

  void limited_image_search(Point_2d &y

                   ,std::vector<size_t> &cell_numbers

                   ,std::vector<Triangle> &triangles

  ) const;

};


template<typename T>


PixelMap<T> GridMap::getPixelMapFlux() const{


  // The number of pixels on a side of the new map will be

  // N = (Ngrid_init-1)/resf + 1;

  // so that the resolution is resf x the GridMap resolution


  PixelMap<T> map(center.x

                  ,Ngrid_init

                  ,Ngrid_init2

                  ,x_range/(Ngrid_init-1));


  size_t index;

   size_t n = Ngrid_init*Ngrid_init2;

   for(size_t i=0 ; i<n ; ++i){

     index = map.find_index(i_points[i].x);

     map.data()[index] = i_points[i].surface_brightness;

   }


  //for(size_t i = 0 ; i < Ngrid_init ; ++i){

  //  for(size_t j = 0 ; j < Ngrid_init2 ; ++j){

  //    map.data()[i/resf + map.getNx() * (j / resf)] +=

  //    i_points[ i + Ngrid_init * j].surface_brightness/factor;

  //  }

  //}


  map.Renormalize(map.getResolution()*map.getResolution());


  return map;

}


template<typename T>


void GridMap::getPixelMapFlux(PixelMap<T> &map) const{


  //int resf = (Ngrid_init-1)/(map.getNx()-1);


  if(map.getNx() != Ngrid_init) throw std::invalid_argument("PixelMap does not match GripMap! Use the other GridMap::getPixelMapFlux() to contruct a PixelMap.");

  if(map.getNy() != Ngrid_init2) throw std::invalid_argument("PixelMap does not match GripMap! Use the other GridMap::getPixelMapFlux() to contruct a PixelMap.");

  //if(map.getResolution() != x_range*resf/(Ngrid_init-1)) throw std::invalid_argument("PixelMap does not match GripMap resolution! Use the other GridMap::getPixelMapFlux() to contruct a PixelMap.");


  if(map.getCenter()[0] != center[0]) throw std::invalid_argument("PixelMap does not match GripMap!");

  if(map.getCenter()[1] != center[1]) throw std::invalid_argument("PixelMap does not match GripMap!");


  double res = getResolution();

  if((map.getResolution()-res) < 1.0e-6*res ) throw std::invalid_argument("PixelMap resolution does not match GripMap!");


  map.Clean();


  //long nx = map.getNx();

  //for(size_t i = 0 ; i < Ngrid_init ; ++i){

  //  for(size_t j = 0 ; j < Ngrid_init2 ; ++j){

  //    size_t k = i + nx * j;

  //    map.data()[k] += i_points[k].surface_brightness;

  //  }

  //}


   size_t index;

   size_t n = Ngrid_init*Ngrid_init2;

   for(size_t i=0 ; i<n ; ++i){

     index = map.find_index(i_points[i].x);

     map.data()[index] = i_points[i].surface_brightness;

   }


  map.Renormalize(map.getResolution()*map.getResolution());

}


template<typename T>


PixelMap<T> GridMap::writePixelMapUniform(

                                       const PosType center[]

                                       ,size_t Nx

                                       ,size_t Ny

                                       ,LensingVariable lensvar

){


  if(getNumberOfPoints() == 0 ) return PixelMap<T>();

  PixelMap<T> map(center, Nx, Ny,x_range/(Nx-1));


  map.Clean();


  writePixelMapUniform(map,lensvar);


  return map;

}


template<typename T>


void GridMap::writeFits(

                        LensingVariable lensvar

                        ,std::string filename

                        ){

                          PixelMap<T> map = writePixelMap<T>(lensvar);

                          map.printFITS(filename);

}


template<typename T>


PixelMap<T> GridMap::writePixelMap(

              LensingVariable lensvar

){

  size_t Nx =  Ngrid_init;

  size_t Ny = Ngrid_init2;


  PixelMap<T> map( center.x, Nx, Ny,x_range/(Nx-1) );


  size_t N = map.size();

  assert(N == Nx*Ny);


  double tmp2[2];

  switch (lensvar) {

    case LensingVariable::ALPHA:

      for(size_t i=0 ; i<N ; ++i){

        tmp2[0] = i_points[i].x[0] - i_points[i].image->x[0];

        tmp2[1] = i_points[i].x[1] - i_points[i].image->x[1];

        map[i] = sqrt(tmp2[0]*tmp2[0] + tmp2[1]*tmp2[1]);

      }

      break;

    case LensingVariable::ALPHA1:

      for(size_t i=0 ; i<N ; ++i)

        map[i] = (i_points[i].x[0] - i_points[i].image->x[0]);

      break;

    case LensingVariable::ALPHA2:

      for(size_t i=0 ; i<N ; ++i)

        map[i] = (i_points[i].x[1] - i_points[i].image->x[1]);

      break;

    case LensingVariable::KAPPA:

      for(size_t i=0 ; i<N ; ++i)

        map[i] = i_points[i].kappa();

      break;

    case LensingVariable::GAMMA:

       for(size_t i=0 ; i<N ; ++i){

         tmp2[0] = i_points[i].gamma1();

         tmp2[1] = i_points[i].gamma2();

         map[i] = sqrt(tmp2[0]*tmp2[0] + tmp2[1]*tmp2[1]);

      }

      break;

    case LensingVariable::GAMMA1:

      for(size_t i=0 ; i<N ; ++i)

        map[i] = i_points[i].gamma1();

      break;

    case LensingVariable::GAMMA2:

      for(size_t i=0 ; i<N ; ++i)

        map[i] = i_points[i].gamma2();

      break;

    case LensingVariable::GAMMA3:

      for(size_t i=0 ; i<N ; ++i)

        map[i] = i_points[i].gamma3();

      break;

    case LensingVariable::INVMAG:

      for(size_t i=0 ; i<N ; ++i)

        map[i] = i_points[i].invmag();

      break;

    case LensingVariable::DELAYT:

      for(size_t i=0 ; i<N ; ++i)

        map[i] = i_points[i].dt;

      break;

    case LensingVariable::SurfBrightness:

      for(size_t i=0 ; i<N ; ++i)

        map[i] = i_points[i].surface_brightness;

      break;

    default:

      std::cerr << "GridMap::writePixelMapUniform() does not work for the input LensingVariable" << std::endl;

      throw std::runtime_error("GridMap::writePixelMapUniform() does not work for the input LensingVariable");

      break;

      // If this list is to be expanded to include ALPHA or GAMMA take care to add them as vectors

  }

  return map;

}


template <typename T>


void GridMap::writePixelMapUniform(

                                   PixelMap<T> &map

                                   ,LensingVariable lensvar

){


  if(getNumberOfPoints() ==0 ) return;


  map.Clean();


  //writePixelMapUniform_(i_points,getNumberOfPoints(),&map,lensvar);

  //return;


  std::vector<std::thread> thr;

  int nthreads = Utilities::GetNThreads();


  int chunk_size;

  do{

    chunk_size =  getNumberOfPoints()/nthreads;

    if(chunk_size == 0) nthreads /= 2;

  }while(chunk_size == 0);


  size_t size = chunk_size;

  for(int ii = 0; ii < nthreads ;++ii){

    if(ii == nthreads-1)

    size = getNumberOfPoints() - (nthreads-1)*chunk_size;

    thr.push_back(std::thread(&GridMap::writePixelMapUniform_<T>,this,&(i_points[ii*chunk_size]),size,&map,lensvar));

  }

  for(auto &t : thr) t.join();

}


template <typename T>

void GridMap::writePixelMapUniform_(Point* points,size_t size,PixelMap<T> *map,LensingVariable val){

  double tmp;

  PosType tmp2[2];

  long index;


  for(size_t i = 0; i< size; ++i){

    switch (val) {

      case LensingVariable::ALPHA:

        tmp2[0] = points[i].x[0] - points[i].image->x[0];

        tmp2[1] = points[i].x[1] - points[i].image->x[1];

        tmp = sqrt(tmp2[0]*tmp2[0] + tmp2[1]*tmp2[1]);

        break;

      case LensingVariable::ALPHA1:

        tmp = (points[i].x[0] - points[i].image->x[0]);

        break;

      case LensingVariable::ALPHA2:

        tmp = (points[i].x[1] - points[i].image->x[1]);

        break;

      case LensingVariable::KAPPA:

        tmp = points[i].kappa();

        break;

      case LensingVariable::GAMMA:

        tmp2[0] = points[i].gamma1();

        tmp2[1] = points[i].gamma2();

        tmp = sqrt(tmp2[0]*tmp2[0] + tmp2[1]*tmp2[1]);

        break;

      case LensingVariable::GAMMA1:

        tmp = points[i].gamma1();

        break;

      case LensingVariable::GAMMA2:

        tmp = points[i].gamma2();

        break;

      case LensingVariable::GAMMA3:

        tmp = points[i].gamma3();

        break;

      case LensingVariable::INVMAG:

        tmp = points[i].invmag();

        break;

      case LensingVariable::DELAYT:

        tmp = points[i].dt;

        break;

      case LensingVariable::SurfBrightness:

        tmp = points[i].surface_brightness;

        break;

      default:

        std::cerr << "PixelMap<T>::AddGrid() does not work for the input LensingVariable" << std::endl;

        throw std::runtime_error("PixelMap<T>::AddGrid() does not work for the input LensingVariable");

        break;

        // If this list is to be expanded to include ALPHA or GAMMA take care to add them as vectors

    }


    index = map->find_index(points[i].x);

    if(index != -1)(*map)[index] = tmp;

  }

}


template <typename T>


void GridMap::writeFitsUniform(

                               const PosType center[]

                               ,size_t Nx

                               ,size_t Ny

                               ,LensingVariable lensvar

                               ,std::string filename

){

  std::string tag;


  switch (lensvar) {

    case LensingVariable::DELAYT:

      tag = ".dt.fits";

      break;

    case LensingVariable::ALPHA1:

      tag = ".alpha1.fits";

      break;

    case LensingVariable::ALPHA2:

      tag = ".alpha2.fits";

      break;

    case LensingVariable::ALPHA:

      tag = ".alpha.fits";

      break;

    case LensingVariable::KAPPA:

      tag = ".kappa.fits";

      break;

    case LensingVariable::GAMMA1:

      tag = ".gamma1.fits";

      break;

    case LensingVariable::GAMMA2:

      tag = ".gamma2.fits";

      break;

    case LensingVariable::GAMMA3:

      tag = ".gamma3.fits";

      break;

    case LensingVariable::GAMMA:

      tag = ".gamma.fits";

      break;

    case LensingVariable::INVMAG:

      tag = ".invmag.fits";

      break;

    case LensingVariable::SurfBrightness:

      tag = ".surfbright.fits";

      break;

 default:

      break;

  }


  PixelMap<T> map = writePixelMapUniform<T>(center,Nx,Ny,lensvar);

  map.printFITS(filename + tag);

}


#endif // defined(__GLAMER__gridmap__)

PixelMap
Image structure that can be manipulated and exported to/from fits files.
Definition pixelmap.h:42

PixelMap::Renormalize
void Renormalize(T factor)
Multiplies the whole map by a scalar factor.
Definition pixelmap.cpp:509

PixelMap::find_index
long find_index(PosType const x[], long &ix, long &iy) const
get the index for a position, returns -1 if out of map, this version returns the 2D grid coordinates
Definition pixelmap.cpp:2227

PixelMap::printFITS
void printFITS(std::string filename, bool Xflip=false, bool ctype=false, bool verbose=false)
Output the pixel map as a fits file.
Definition pixelmap.cpp:1318

Source
Base class for all sources.
Definition source.h:44

Utilities::find_boundaries
void find_boundaries(std::vector< bool > &bitmap, long nx, std::vector< std::vector< P > > &points, std::vector< bool > &hits_edge, bool add_to_vector=false, bool outer_only=false)
finds ordered boundaries to regions where bitmap == true
Definition concave_hull.h:201

Utilities::GetNThreads
int GetNThreads()
returns the compiler variable N_THREADS that is maximum number of threads to be used.
Definition utilities.cpp:553

LensingVariable
LensingVariable
output lensing variables
Definition standard.h:89

LensingVariable::ALPHA
@ ALPHA
magnitude of deflection in radians
Definition standard.h:91

LensingVariable::GAMMA2
@ GAMMA2
second component of shear
Definition standard.h:97

LensingVariable::DELAYT
@ DELAYT
time delay
Definition standard.h:90

LensingVariable::GAMMA3
@ GAMMA3
third component of shear
Definition standard.h:98

LensingVariable::ALPHA2
@ ALPHA2
y component of deflection
Definition standard.h:93

LensingVariable::SurfBrightness
@ SurfBrightness
Surface brightness.
Definition standard.h:101

LensingVariable::GAMMA1
@ GAMMA1
first component of shear
Definition standard.h:96

LensingVariable::GAMMA
@ GAMMA
magnitude of shear
Definition standard.h:95

LensingVariable::KAPPA
@ KAPPA
convergence
Definition standard.h:94

LensingVariable::ALPHA1
@ ALPHA1
x component of deflection
Definition standard.h:92

LensingVariable::INVMAG
@ INVMAG
inverse of magnification
Definition standard.h:99

GridMap::getPixelMapFlux
PixelMap< T > getPixelMapFlux() const
returns a PixelMap with the flux in pixels at a resolution of res times the original resolution
Definition gridmap.h:431

GridMap::AddSurfaceBrightnesses
double AddSurfaceBrightnesses(Source *source)
Recalculate surface brightness just like GridMap::RefreshSurfaceBrightness but the new source is adde...
Definition gridmap.cpp:260

GridMap::image_point
Point_2d image_point(size_t index)
get the image point for a index number
Definition gridmap.h:72

GridMap::find_images
std::list< RAY > find_images(std::vector< Point_2d > &ys, std::vector< int > &multiplicity) const
Returns a list of RAYs from a set of source positions.
Definition gridmap.cpp:927

GridMap::EinsteinArea
PosType EinsteinArea() const
returns the area (radians^2) of the region with negative magnification at resolution of fixed grid
Definition gridmap.cpp:323

GridMap::ray
RAY ray(size_t index)
get the image point for a index number
Definition gridmap.h:76

GridMap::getXRange
double getXRange() const
return initial range of gridded region. This is the distance from the first ray in a row to the last ...
Definition gridmap.h:85

GridMap::RefreshSurfaceBrightnesses
double RefreshSurfaceBrightnesses(Source *source)
Recalculate surface brightness at every point without changing the positions of the gridmap or any le...
Definition gridmap.cpp:174

GridMap::AddPointSource
double AddPointSource(const Point_2d &y, double flux)
add flux to the rays that are nearest to the source on the source plane for each image
Definition gridmap.cpp:466

GridMap::source_point
Point_2d source_point(size_t index)
get the image point for a index number
Definition gridmap.h:74

GridMap::find_magnification_contour
void find_magnification_contour(std::vector< std::vector< Point_2d > > &curves, std::vector< bool > &hits_boundary, double invmag)
Find image-plane contours of magnification.   This is usually only used within ImageFinding:: functio...
Definition gridmap.cpp:494

GridMap::writePixelMap
PixelMap< T > writePixelMap(LensingVariable lensvar)
make pixel map of lensing quantities at the resolution of the GridMap
Definition gridmap.h:533

GridMap::writeFitsUniform
void writeFitsUniform(const PosType center[], size_t Nx, size_t Ny, LensingVariable lensvar, std::string filename)
Definition gridmap.h:693

GridMap::writeFitsUniform
void writeFitsUniform(LensingVariable lensvar, std::string filename)
this will make a fits map of the grid as is.
Definition gridmap.h:106

GridMap::getResolution
double getResolution() const
resolution in radians, this is range / (N-1)
Definition gridmap.h:88

GridMap::ReInitialize
GridMap ReInitialize(LensHndl lens)
reshoot the rays for example when the source plane has been changed
Definition gridmap.cpp:148

GridMap::AdaptiveRefreshSurfaceBrightnesses
double AdaptiveRefreshSurfaceBrightnesses(Lens &lens, Source &source)
Definition gridmap.cpp:189

GridMap::magnificationTr
double magnificationTr() const
calculate the LOCAL magnification by triangel method weighted by interpolated surface brightness
Definition gridmap.cpp:355

GridMap::writePixelMapUniform
void writePixelMapUniform(PixelMap< T > &map, LensingVariable lensvar)
Definition gridmap.h:605

GridMap::magnificationFlux
PosType magnificationFlux(Source &source) const
Definition gridmap.cpp:344

GridMap::getInitNgrid
int getInitNgrid() const
return initial number of grid points in each direction
Definition gridmap.h:83

GridMap::GridMap
GridMap(LensHndl lens, unsigned long N1d, const double center[2], double range)
Constructor for initializing square grid.
Definition gridmap.cpp:93

GridMap::find_boundaries_of_caustics
void find_boundaries_of_caustics(std::vector< std::vector< Point_2d > > &boundaries, std::vector< bool > &hits_edge)
finds the boundary of the region on the source plane where there are more than one image
Definition gridmap.h:278

GridMap::AreaCellOnSourcePlane
double AreaCellOnSourcePlane(size_t k) const
area of a cell (pixel size region with its lower left at point k) on source plane - calculated by tri...
Definition gridmap.cpp:452

GridMap::writeFits
void writeFits(LensingVariable lensvar, std::string filensame)
fits output of lensing quantities at the resolution of the GridMap
Definition gridmap.h:524

GridMap::centroid
Point_2d centroid() const
returns centroid of flux on the grid
Definition gridmap.cpp:915

GridMap::find_crit
void find_crit(std::vector< std::vector< Point_2d > > &points, std::vector< bool > &hits_boundary, std::vector< CritType > &crit_type)
Find critical curves. This is usually not used outside of ImageFinding::find_crit()
Definition gridmap.cpp:527

GridMap::deLens
void deLens()
resets to state without lensing
Definition gridmap.cpp:162

Point_2d
Class for representing points or vectors in 2 dimensions. Not that the dereferencing operator is over...
Definition point.h:48

Point
A point on the source or image plane that contains a position and the lensing quantities.
Definition point.h:421

RAY
Simple representaion of a light path giving position on the image and source planes and lensing quant...
Definition point.h:517