gridmap.h

//
//  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(SourceHndl source);
  double AdaptiveRefreshSurfaceBrightnesses(Lens &lens,Source &source);
  
  double AddSurfaceBrightnesses(SourceHndl source);
 
  Point_2d image_point(size_t index){return i_points[index];}
  Point_2d source_point(size_t index){return s_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);}
  
  PixelMap writePixelMap(LensingVariable lensvar);
  void writeFits(LensingVariable lensvar,std::string filensame);
 
  void writePixelMapUniform(PixelMap &map,LensingVariable lensvar);
  void writeFitsUniform(const PosType center[],size_t Nx,size_t Ny,LensingVariable lensvar,std::string filename);
  PixelMap writePixelMapUniform(const PosType center[],size_t Nx,size_t Ny,LensingVariable lensvar);
 
  void writeFitsUniform(
                        LensingVariable lensvar    
                        ,std::string filename     
                        ){
    PixelMap map = writePixelMap(lensvar);
    map.printFITS(filename);
  }
  
  PixelMap getPixelMapFlux(int res) const;
  void getPixelMapFlux(PixelMap &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;
    }
    size_t index[3];
    size_t & operator[](int i){return index[i];}
  };
  
    
  /*** \brief Returns a list of  RAYs from a set of source positions.
   
   The image positions are found in parallel by the triangle method.  The order of the
   output rays will be the same as the sources with multiple images consecutive.
   */
  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_crit_boundary(std::vector<std::vector<Point_2d> > &points
//                          ,std::vector<bool> &hits_boundary
//                          ) const;
  
private:
  GridMap & operator=(GridMap &grid);
  
  /* Depricated to Utilities::find_boundaries<>()
   
   finds ordered boundaries to regions where bitmap == true
 
   This can be used to find critical curves or contours.
   `bitmap` should be the same size as the `Gridmap`
   If the boundary curve  touches the edge of the `GridMap` it will be indicated in `hits_boundary` as
   `true`.
   
   Boundaries will never cross or lead off the grid.  On the edges they will leave the edge pixels out even if they should be in.  This is a technical comprimise.
  */
//  void find_boundaries(std::vector<bool> &bitmap  // = true inside
//                       ,std::vector<std::vector<Point_2d> > &points
//                       ,std::vector<bool> &hits_edge
//                       ,bool add_to_vector=false
//                       );
//
  
  // 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;
  void writePixelMapUniform_(Point* points,size_t size,PixelMap *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;
 
  // 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;
};
 
#endif // defined(__GLAMER__gridmap__)