GLAMERDOC++
Gravitational Lensing Code Library
gridmap.h
1 //
2 // gridmap.h
3 // GLAMER
4 //
5 // Created by bmetcalf on 8/19/14.
6 //
7 //
8 
9 #ifndef __GLAMER__gridmap__
10 #define __GLAMER__gridmap__
11 
12 #include <iostream>
13 #include <memory>
14 #include "lens.h"
15 #include "point.h"
16 #include "concave_hull.h"
17 #include "Tree.h"
18 #include "source.h"
19 #include <mutex>
20 
31 struct GridMap{
32 
33  friend class Lens;
34 
35  GridMap(LensHndl lens,unsigned long N1d,const double center[2],double range);
36  GridMap(LensHndl lens ,unsigned long Nx ,const PosType center[2] ,PosType rangeX ,PosType rangeY);
38  GridMap(unsigned long N1d,const double center[2],double range);
39  ~GridMap();
40 
42  GridMap ReInitialize(LensHndl lens);
43 
45  void deLens();
55  double RefreshSurfaceBrightnesses(Source* source);
61  double AdaptiveRefreshSurfaceBrightnesses(Lens &lens,Source &source);
62 
69  double AddSurfaceBrightnesses(Source* source);
70 
72  Point_2d image_point(size_t index){return i_points[index];}
74  Point_2d source_point(size_t index){return s_points[index];}
76  RAY ray(size_t index){return i_points[index];}
77 
78  void ClearSurfaceBrightnesses();
79  void assertNAN(); // check for nan in surface prightness
80  size_t getNumberOfPoints() const {return Ngrid_init*Ngrid_init2;}
81 
83  int getInitNgrid() const {return Ngrid_init;}
85  double getXRange() const {return x_range;}
86  double getYRange() const {return x_range*axisratio;}
88  double getResolution() const {return x_range/(Ngrid_init-1);}
89 
91  template<typename T>
92  PixelMap<T> writePixelMap(LensingVariable lensvar);
94  template <typename T>
95  void writeFits(LensingVariable lensvar,std::string filensame,bool flipX = false);
96 
97  template<typename T>
98  void writePixelMapUniform(PixelMap<T> &map,LensingVariable lensvar,bool flipX = false);
99  template <typename T>
100  void writeFitsUniform(const PosType center[],size_t Nx,size_t Ny,LensingVariable lensvar
101  ,std::string filename,bool flipX = false);
102  template<typename T>
103  PixelMap<T> writePixelMapUniform(const PosType center[],size_t Nx,size_t Ny,LensingVariable lensvar);
104 
106  template<typename T>
107  void writeFitsUniform(
108  LensingVariable lensvar
109  ,std::string filename
110  ,bool xflip = false
111  ){
112  PixelMap<T> map = writePixelMap<T>(lensvar);
113  map.printFITS(filename,xflip);
114  }
115 
117  template<typename T>
118  PixelMap<T> getPixelMapFlux() const;
119 
123  template<typename T>
124  void getPixelMapFlux(PixelMap<T> &map) const;
125 
127  PosType EinsteinArea() const;
128 
131  //PosType magnification_invmag() const;
132  //PosType magnification2() const;
133 
135  Point_2d centroid() const;
136 
137  Point_2d getCenter() const {return center;}
138 
139  Point * operator[](size_t i){return i_points.data() + i;};
140 
141  GridMap(GridMap &&grid){
142 
143  Ngrid_init = grid.Ngrid_init;
144  Ngrid_init2 = grid.Ngrid_init2;
145  pointID = grid.pointID;
146  axisratio = grid.axisratio;
147  x_range = grid.x_range;
148 
149  std::swap(i_points,grid.i_points);
150  std::swap(s_points,grid.s_points);
151 
152  center = grid.center;
153  }
154 
155  GridMap(GridMap &grid){
156 
157  Ngrid_init = grid.Ngrid_init;
158  Ngrid_init2 = grid.Ngrid_init2;
159  pointID = grid.pointID;
160  axisratio = grid.axisratio;
161  x_range = grid.x_range;
162 
163  i_points = grid.i_points;
164  s_points = grid.s_points;
165 
166  center = grid.center;
167  }
168 
169  GridMap & operator=(GridMap &&grid){
170  Ngrid_init = grid.Ngrid_init;
171  Ngrid_init2 = grid.Ngrid_init2;
172  pointID = grid.pointID;
173  axisratio = grid.axisratio;
174  x_range = grid.x_range;
175 
176  std::swap(i_points,grid.i_points);
177  std::swap(s_points,grid.s_points);
178 
179  center = grid.center;
180 
181  return *this;
182  }
183 
184  struct Triangle{
185  Triangle(size_t i,size_t j,size_t k){
186  index[0] = i;
187  index[1] = j;
188  index[2] = k;
189  }
190  Triangle(){
191  index[0] = 0;
192  index[1] = 0;
193  index[2] = 0;
194  }
195  Triangle(const Triangle &tri){
196  index[0] = tri[0];
197  index[1] = tri[1];
198  index[2] = tri[2];
199  }
200  Triangle & operator=(const Triangle &tri){
201  index[0] = tri[0];
202  index[1] = tri[1];
203  index[2] = tri[2];
204  return *this;
205  }
206  ~Triangle(){};
207  std::vector<long> index = {0,0,0};
208  long & operator[](int i){return index[i];}
209  long operator[](int i) const {return index[i];}
210  };
211  struct Rectangle{
212  Rectangle(size_t i,size_t j){
213  index[0] = i;
214  index[1] = j;
215  }
216  Rectangle(){
217  index[0] = 0;
218  index[1] = 0;
219  }
220  Rectangle(const Rectangle &rec){
221  index[0] = rec[0];
222  index[1] = rec[1];
223  }
224  Rectangle & operator=(Rectangle &rec){
225  index[0] = rec[0];
226  index[1] = rec[1];
227  return *this;
228  }
229  ~Rectangle(){};
230  std::vector<long> index = {0,0};
231  long & operator[](int i){return index[i];}
232  long operator[](int i) const {return index[i];}
233  };
234 
235  // determines if two rectangles touch
236  inline bool touch(const Rectangle &tr1,const Rectangle &tr2) const;
237 
238  // returns a vector of rectangles that encompase all triangles with one rectangle for each touching group
239  std::vector<Rectangle> merge_boxes(
240  std::list<Triangle> &triangles
241  ) const;
242  std::vector<Rectangle> merge_boxes(
243  std::vector<Triangle> &triangles
244  ) const;
245 
256  std::list<RAY> find_images(std::vector<Point_2d> &ys
257  ,std::vector<int> &multiplicity
258  ) const;
259 
266  void find_images(Point_2d y
267  ,std::vector<Point_2d> &image_points
268  ,std::vector<Triangle> &triangles
269  ) const ;
270 
271 
280  void find_boundaries_of_caustics(std::vector<std::vector<Point_2d> > &boundaries
281  ,std::vector<bool> &hits_edge
282  ){
283 
284  size_t N=Ngrid_init*Ngrid_init2;
285  std::vector<Point_2d> source_points(N);
286  for(size_t i=0 ; i<N ; ++i) source_points[i] = i_points[i];
287 
288  std::vector<int> multiplicities(N);
289  find_images(source_points,multiplicities);
290 
291  std::vector<bool> bitmap(N,false);
292  for(size_t i=0 ; i<N ; ++i){
293  if(multiplicities[i] > 1) bitmap[i] = true;
294  }
295 
296  Utilities::find_boundaries<Point_2d>(bitmap,Ngrid_init,boundaries,hits_edge);
297  }
298 
312  Point_2d magnificationFlux(Source &source) const ;
313 
321  double magnificationTr() const ;
322 
325  double magnificationTr(std::vector<size_t> &pixels) const ;
326 
328  double AreaCellOnSourcePlane(size_t k) const;
329 
330 
331 // depricated version of PixalMap::magnificationTr()
332 // double magnificationTr2() const {
333 //
334 // size_t k;
335 // size_t k1;
336 // size_t k2;
337 //
338 // double flux_source=0.0,flux_image=0.0;
339 // for(size_t i=1 ; i< Ngrid_init-1 ; i++){
340 // for(size_t j=1 ; j< Ngrid_init2-1 ; j++){
341 //
342 // k = i + j*Ngrid_init;
343 //
344 // assert( s_points[k].surface_brightness == i_points[k].surface_brightness);
345 // double sb = s_points[k].surface_brightness;
346 //
347 // if(sb !=0.0){
348 // flux_source += AreaCellOnSourcePlane(k) *sb;
349 // flux_source += AreaCellOnSourcePlane(k-1) *sb;
350 // flux_source += AreaCellOnSourcePlane(k-1-Ngrid_init) *sb;
351 // flux_source += AreaCellOnSourcePlane(k-Ngrid_init) *sb;
352 // }
353 //
354 // flux_image += sb;
355 // }
356 // }
357 //
358 // return flux_image * getResolution() * getResolution() * 4 / flux_source;
359 // }
360 
367  double AddPointSource(const Point_2d &y,double flux);
368 
374  void find_crit(std::vector<std::vector<Point_2d> > &points
375  ,std::vector<bool> &hits_boundary
376  ,std::vector<CritType> &crit_type
377  );
378 
382  void find_magnification_contour(
383  std::vector<std::vector<Point_2d> > &curves
384  ,std::vector<bool> &hits_boundary
385  ,double invmag
386  );
395  void border_area(double &source_area,double &image_area) const;
396 
397  int getNx(){return Ngrid_init;}
398  int getNy(){return Ngrid_init2;}
399 private:
400  GridMap & operator=(GridMap &grid);
401 
402  // curve must be in pixel units
403  bool incurve(long k,std::vector<Point_2d> &curve) const;
404 
405  // cluge to make compatible with old method of producing points
406  std::vector<Point> NewPointArray(size_t N){
407  std::vector<Point> p(N);
408  p[0].head = N;
409  return p;
410  }
411 
412  void xygridpoints(Point *points,double range,const double *center,long Ngrid
413  ,short remove_center);
414 
416  int Ngrid_init;
417  int Ngrid_init2;
418 
419  unsigned long pointID;
420  PosType axisratio;
421  PosType x_range;
422  template<typename T>
423  void writePixelMapUniform_(Point* points,size_t size,PixelMap<T> *map,LensingVariable val);
424 
425  std::vector<Point> i_points;
426  std::vector<Point> s_points;
427  Point_2d center;
428 
429  bool to_refine(long i,long j,double total,double f) const ;
430  static std::mutex grid_mutex;
431 
432  void _find_images_(Point_2d *ys,int *multiplicity,long Nys,std::list<RAY> &rays) const;
433  void _find_images2_(size_t j1
434  ,size_t j2
435  ,std::list<Point_2d> &image_points
436  ,std::list<Triangle> &triangles
437  ,Point_2d y
438  ) const;
439 
440  void _magnificationFlux_parallel(
441  long start
442  ,long end
443  ,double &magnified_flux
444  ,double &unmagnified_flux
445  ,Source &source
446  //,Point_2d &recenter
447  ) const;
448 
449  void _AddSurfaceBrightnesses_parallel(
450  long start
451  ,long end
452  ,double &flux
453  ,Source *source
454  );
455 
456  // find if there are images of y in specific cells
457  void limited_image_search(Point_2d &y
458  ,std::vector<size_t> &cell_numbers
459  ,std::vector<Triangle> &triangles
460  ) const;
461 };
462 
464 template<typename T>
465 PixelMap<T> GridMap::getPixelMapFlux() const{
466 
467  // The number of pixels on a side of the new map will be
468  // N = (Ngrid_init-1)/resf + 1;
469  // so that the resolution is resf x the GridMap resolution
470 
471  PixelMap<T> map(center.x
472  ,Ngrid_init
473  ,Ngrid_init2
474  ,x_range/(Ngrid_init-1));
475 
476  size_t n = Ngrid_init*Ngrid_init2;
477  for(size_t i=0 ; i<n ; ++i){
478  map.data()[map.find_index(i_points[i].x)] = i_points[i].surface_brightness;
479  }
480 
481  //for(size_t i = 0 ; i < Ngrid_init ; ++i){
482  // for(size_t j = 0 ; j < Ngrid_init2 ; ++j){
483  // map.data()[i/resf + map.getNx() * (j / resf)] +=
484  // i_points[ i + Ngrid_init * j].surface_brightness/factor;
485  // }
486  //}
487 
488  map.Renormalize(map.getResolution()*map.getResolution());
489 
490  return map;
491 }
492 
494 template<typename T>
495 void GridMap::getPixelMapFlux(PixelMap<T> &map) const{
496 
497  //int resf = (Ngrid_init-1)/(map.getNx()-1);
498 
499  if(map.getNx() != Ngrid_init) throw std::invalid_argument("PixelMap does not match GripMap! Use the other GridMap::getPixelMapFlux() to contruct a PixelMap.");
500  if(map.getNy() != Ngrid_init2) throw std::invalid_argument("PixelMap does not match GripMap! Use the other GridMap::getPixelMapFlux() to contruct a PixelMap.");
501  //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.");
502 
503  if(map.getCenter()[0] != center[0]) throw std::invalid_argument("PixelMap does not match GripMap!");
504  if(map.getCenter()[1] != center[1]) throw std::invalid_argument("PixelMap does not match GripMap!");
505 
506  double res = getResolution();
507  if((map.getResolution()-res) < 1.0e-6*res ) throw std::invalid_argument("PixelMap resolution does not match GripMap!");
508 
509  map.Clean();
510 
511  //long nx = map.getNx();
512  //for(size_t i = 0 ; i < Ngrid_init ; ++i){
513  // for(size_t j = 0 ; j < Ngrid_init2 ; ++j){
514  // size_t k = i + nx * j;
515  // map.data()[k] += i_points[k].surface_brightness;
516  // }
517  //}
518 
519  size_t index;
520  size_t n = Ngrid_init*Ngrid_init2;
521  for(size_t i=0 ; i<n ; ++i){
522  index = map.find_index(i_points[i].x);
523  map.data()[index] = i_points[i].surface_brightness;
524  }
525 
526 
527  map.Renormalize(map.getResolution()*map.getResolution());
528 }
529 
537 template<typename T>
538 PixelMap<T> GridMap::writePixelMapUniform(
539  const PosType center[]
540  ,size_t Nx
541  ,size_t Ny
542  ,LensingVariable lensvar
543 ){
544 
545  if(getNumberOfPoints() == 0 ) return PixelMap<T>();
546  PixelMap<T> map(center, Nx, Ny,x_range/(Nx-1));
547 
548  map.Clean();
549 
550  writePixelMapUniform(map,lensvar);
551 
552  return map;
553 }
554 
555 template<typename T>
556 void GridMap::writeFits(
557  LensingVariable lensvar
558  ,std::string filename
559  ,bool xflip
560  ){
561  PixelMap<T> map = writePixelMap<T>(lensvar);
562  map.printFITS(filename,xflip);
563 }
564 
565 template<typename T>
566 PixelMap<T> GridMap::writePixelMap(
567  LensingVariable lensvar
568 ){
569  size_t Nx = Ngrid_init;
570  size_t Ny = Ngrid_init2;
571 
572  PixelMap<T> map( center.x, Nx, Ny,x_range/(Nx-1) );
573 
574  size_t N = map.size();
575  assert(N == Nx*Ny);
576  Point_2d p1,p2;
577 
578  double tmp2[2];
579  switch (lensvar) {
580  case LensingVariable::ALPHA:
581  for(size_t i=0 ; i<N ; ++i){
582  tmp2[0] = i_points[i].x[0] - i_points[i].image->x[0];
583  tmp2[1] = i_points[i].x[1] - i_points[i].image->x[1];
584  map[i] = sqrt(tmp2[0]*tmp2[0] + tmp2[1]*tmp2[1]);
585  }
586  break;
587  case LensingVariable::ALPHA1:
588  for(size_t i=0 ; i<N ; ++i)
589  map[i] = (i_points[i].x[0] - i_points[i].image->x[0]);
590  break;
591  case LensingVariable::ALPHA2:
592  for(size_t i=0 ; i<N ; ++i)
593  map[i] = (i_points[i].x[1] - i_points[i].image->x[1]);
594  break;
595  case LensingVariable::KAPPA:
596  for(size_t i=0 ; i<N ; ++i)
597  map[i] = i_points[i].kappa();
598  break;
599  case LensingVariable::GAMMA:
600  for(size_t i=0 ; i<N ; ++i){
601  tmp2[0] = i_points[i].gamma1();
602  tmp2[1] = i_points[i].gamma2();
603  map[i] = sqrt(tmp2[0]*tmp2[0] + tmp2[1]*tmp2[1]);
604  }
605  break;
606  case LensingVariable::GAMMA1:
607  for(size_t i=0 ; i<N ; ++i)
608  map[i] = i_points[i].gamma1();
609  break;
610  case LensingVariable::GAMMA2:
611  for(size_t i=0 ; i<N ; ++i)
612  map[i] = i_points[i].gamma2();
613  break;
614  case LensingVariable::GAMMA3:
615  for(size_t i=0 ; i<N ; ++i)
616  map[i] = i_points[i].gamma3();
617  break;
618  case LensingVariable::INVMAG:
619  for(size_t i=0 ; i<N ; ++i)
620  map[i] = i_points[i].invmag();
621  break;
622  case LensingVariable::DELAYT:
623  for(size_t i=0 ; i<N ; ++i)
624  map[i] = i_points[i].dt;
625  break;
626  case LensingVariable::SurfBrightness:
627  for(size_t i=0 ; i<N ; ++i)
628  map[i] = i_points[i].surface_brightness;
629  break;
630  case LensingVariable::EigenV:
631  for(size_t i=0 ; i<N ; ++i){
632  i_points[i].A.eigen_vec(p1,p2,tmp2);
633  if(abs(tmp2[0]) > abs(tmp2[1])){
634  map[i] = atan(p2[1]/p2[0]);
635  }else{
636  map[i] = atan(p1[1]/p1[0]);
637  }
638  }
639  break;
640 
641  default:
642  std::cerr << "GridMap::writePixelMapUniform() does not work for the input LensingVariable" << std::endl;
643  throw std::runtime_error("GridMap::writePixelMapUniform() does not work for the input LensingVariable");
644  break;
645  // If this list is to be expanded to include ALPHA or GAMMA take care to add them as vectors
646  }
647  return map;
648 }
649 template <typename T>
650 void GridMap::writePixelMapUniform(
651  PixelMap<T> &map
652  ,LensingVariable lensvar
653  ,bool flipX
654 ){
655 
656  if(getNumberOfPoints() ==0 ) return;
657 
658  map.Clean();
659 
660  //writePixelMapUniform_(i_points,getNumberOfPoints(),&map,lensvar);
661  //return;
662 
663  std::vector<std::thread> thr;
664  int nthreads = Utilities::GetNThreads();
665 
666  size_t Npoints = getNumberOfPoints();
667  if(nthreads >= Npoints) nthreads = 1;
668  size_t chunk_size = (Npoints + nthreads - 1) / nthreads; // divide by threads rounded up
669 
670  size_t size = chunk_size;
671  for(int ii = 0; ii < nthreads ;++ii){
672  long start = ii*chunk_size;
673  long end = start + chunk_size;
674  if(end > Npoints) end = Npoints;
675  thr.push_back(std::thread(&GridMap::writePixelMapUniform_<T>,this,&(i_points[start])
676  ,end-start,&map,lensvar));
677  }
678  for(auto &t : thr) t.join();
679 }
680 
681 template <typename T>
682 void GridMap::writePixelMapUniform_(Point* points,size_t size,PixelMap<T> *map,LensingVariable val){
683  double tmp;
684  PosType tmp2[2];
685  Point_2d p1,p2;
686  long index;
687 
688  for(size_t i = 0; i< size; ++i){
689  switch (val) {
690  case LensingVariable::ALPHA:
691  tmp2[0] = points[i].x[0] - points[i].image->x[0];
692  tmp2[1] = points[i].x[1] - points[i].image->x[1];
693  tmp = sqrt(tmp2[0]*tmp2[0] + tmp2[1]*tmp2[1]);
694  break;
695  case LensingVariable::ALPHA1:
696  tmp = (points[i].x[0] - points[i].image->x[0]);
697  break;
698  case LensingVariable::ALPHA2:
699  tmp = (points[i].x[1] - points[i].image->x[1]);
700  break;
701  case LensingVariable::KAPPA:
702  tmp = points[i].kappa();
703  break;
704  case LensingVariable::GAMMA:
705  tmp2[0] = points[i].gamma1();
706  tmp2[1] = points[i].gamma2();
707  tmp = sqrt(tmp2[0]*tmp2[0] + tmp2[1]*tmp2[1]);
708  break;
709  case LensingVariable::GAMMA1:
710  tmp = points[i].gamma1();
711  break;
712  case LensingVariable::GAMMA2:
713  tmp = points[i].gamma2();
714  break;
715  case LensingVariable::GAMMA3:
716  tmp = points[i].gamma3();
717  break;
718  case LensingVariable::INVMAG:
719  tmp = points[i].invmag();
720  break;
721  case LensingVariable::DELAYT:
722  tmp = points[i].dt;
723  break;
724  case LensingVariable::SurfBrightness:
725  tmp = points[i].surface_brightness;
726  break;
727  case LensingVariable::EigenV:
728  points[i].A.eigen_vec(p1,p2,tmp2);
729  if(abs(tmp2[0]) > abs(tmp2[1])){
730  tmp = atan(p2[1]/p2[0]);
731  }else{
732  tmp = atan(p1[1]/p1[0]);
733  }
734  break;
735 
736  default:
737  std::cerr << "PixelMap<T>::AddGrid() does not work for the input LensingVariable" << std::endl;
738  throw std::runtime_error("PixelMap<T>::AddGrid() does not work for the input LensingVariable");
739  break;
740  // If this list is to be expanded to include ALPHA or GAMMA take care to add them as vectors
741  }
742 
743  index = map->find_index(points[i].x);
744  if(index != -1)(*map)[index] = tmp;
745  }
746 }
747 
748 template <typename T>
749 void GridMap::writeFitsUniform(
750  const PosType center[]
751  ,size_t Nx
752  ,size_t Ny
753  ,LensingVariable lensvar
754  ,std::string filename
755  ,bool xflip
756 ){
757  std::string tag;
758 
759  switch (lensvar) {
760  case LensingVariable::DELAYT:
761  tag = ".dt.fits";
762  break;
763  case LensingVariable::ALPHA1:
764  tag = ".alpha1.fits";
765  break;
766  case LensingVariable::ALPHA2:
767  tag = ".alpha2.fits";
768  break;
769  case LensingVariable::ALPHA:
770  tag = ".alpha.fits";
771  break;
772  case LensingVariable::KAPPA:
773  tag = ".kappa.fits";
774  break;
775  case LensingVariable::GAMMA1:
776  tag = ".gamma1.fits";
777  break;
778  case LensingVariable::GAMMA2:
779  tag = ".gamma2.fits";
780  break;
781  case LensingVariable::GAMMA3:
782  tag = ".gamma3.fits";
783  break;
784  case LensingVariable::GAMMA:
785  tag = ".gamma.fits";
786  break;
787  case LensingVariable::INVMAG:
788  tag = ".invmag.fits";
789  break;
790  case LensingVariable::SurfBrightness:
791  tag = ".surfbright.fits";
792  break;
793  case LensingVariable::EigenV:
794  tag = ".eignangle.fits";
795  break;
796  default:
797  break;
798  }
799 
800  PixelMap<T> map = writePixelMapUniform<T>(center,Nx,Ny,lensvar);
801  map.printFITS(filename + tag,xflip);
802 }
803 
804 #endif // defined(__GLAMER__gridmap__)
LensingVariable::GAMMA3
@ GAMMA3
third component of shear
GridMap::writeFits
void writeFits(LensingVariable lensvar, std::string filensame, bool flipX=false)
fits output of lensing quantities at the resolution of the GridMap
Definition: gridmap.h:556
GridMap::EinsteinArea
PosType EinsteinArea() const
returns the area (radians^2) of the region with negative magnification at resolution of fixed grid
Definition: gridmap.cpp:394
Source
Base class for all sources.
Definition: source.h:44
LensingVariable::ALPHA2
@ ALPHA2
y component of deflection
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:2522
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::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
LensingVariable::GAMMA
@ GAMMA
magnitude of shear
GridMap::magnificationFlux
Point_2d magnificationFlux(Source &source) const
Definition: gridmap.cpp:433
GridMap::centroid
Point_2d centroid() const
returns centroid of flux on the grid
Definition: gridmap.cpp:1052
GridMap::writePixelMap
PixelMap< T > writePixelMap(LensingVariable lensvar)
make pixel map of lensing quantities at the resolution of the GridMap
Definition: gridmap.h:566
LensingVariable::KAPPA
@ KAPPA
convergence
GridMap::writeFitsUniform
void writeFitsUniform(const PosType center[], size_t Nx, size_t Ny, LensingVariable lensvar, std::string filename, bool flipX=false)
Definition: gridmap.h:749
GridMap::source_point
Point_2d source_point(size_t index)
get the image point for a index number
Definition: gridmap.h:74
Point_2d
Class for representing points or vectors in 2 dimensions. Not that the dereferencing operator is over...
Definition: point.h:48
GridMap::deLens
void deLens()
resets to state without lensing
Definition: gridmap.cpp:162
GridMap::AddSurfaceBrightnesses
double AddSurfaceBrightnesses(Source *source)
Recalculate surface brightness just like GridMap::RefreshSurfaceBrightness but the new source is adde...
Definition: gridmap.cpp:300
LensingVariable::ALPHA
@ ALPHA
magnitude of deflection in radians
GridMap
A simplified version of the Grid structure for making non-adaptive maps of the lensing quantities (ka...
Definition: gridmap.h:31
GridMap::magnificationTr
double magnificationTr() const
calculate the LOCAL magnification by triangel method weighted by interpolated surface brightness
Definition: gridmap.cpp:492
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:280
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:1064
LensingVariable::SurfBrightness
@ SurfBrightness
Surface brightness.
LensingVariable::EigenV
@ EigenV
angle of the eigenvector of the magnification matrix with the smallest absolute eigenvalue
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:1613
LensingVariable::ALPHA1
@ ALPHA1
x component of deflection
GridMap::border_area
void border_area(double &source_area, double &image_area) const
Returns the area in square radians within the border of the GridMap on both the source and image plan...
Definition: gridmap.cpp:189
LensingVariable::GAMMA2
@ GAMMA2
second component of shear
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:603
GridMap::ReInitialize
GridMap ReInitialize(LensHndl lens)
reshoot the rays for example when the source plane has been changed
Definition: gridmap.cpp:148
GridMap::image_point
Point_2d image_point(size_t index)
get the image point for a index number
Definition: gridmap.h:72
RAY
Simple representaion of a light path giving position on the image and source planes and lensing quant...
Definition: point.h:555
GridMap::getInitNgrid
int getInitNgrid() const
return initial number of grid points in each direction
Definition: gridmap.h:83
GridMap::getResolution
double getResolution() const
resolution in radians, this is range / (N-1)
Definition: gridmap.h:88
LensingVariable::DELAYT
@ DELAYT
time delay
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:631
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:465
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:664
LensingVariable::GAMMA1
@ GAMMA1
first component of shear
GridMap::AdaptiveRefreshSurfaceBrightnesses
double AdaptiveRefreshSurfaceBrightnesses(Lens &lens, Source &source)
Definition: gridmap.cpp:211
GridMap::GridMap
GridMap(LensHndl lens, unsigned long N1d, const double center[2], double range)
Constructor for initializing square grid.
Definition: gridmap.cpp:93
LensingVariable::INVMAG
@ INVMAG
inverse of magnification
GridMap::writeFitsUniform
void writeFitsUniform(LensingVariable lensvar, std::string filename, bool xflip=false)
this will make a fits map of the grid as is.
Definition: gridmap.h:107
Point
A point on the source or image plane that contains a position and the lensing quantities.
Definition: point.h:459
PixelMap
Image structure that can be manipulated and exported to/from fits files.
Definition: pixelmap.h:42
LensingVariable
LensingVariable
output lensing variables
Definition: standard.h:89
Lens::operator=
Lens & operator=(Lens &&lens)
Definition: lens.h:475
Utilities::GetNThreads
int GetNThreads()
returns the compiler variable N_THREADS that is maximum number of threads to be used.
Definition: utilities.cpp:553
GridMap::writePixelMapUniform
void writePixelMapUniform(PixelMap< T > &map, LensingVariable lensvar, bool flipX=false)
Definition: gridmap.h:650
PixelMap::Renormalize
void Renormalize(T factor)
Multiplies the whole map by a scalar factor.
Definition: pixelmap.cpp:509
Lens
A class to represents a lens with multiple planes.
Definition: lens.h:71
GridMap::ray
RAY ray(size_t index)
get the image point for a index number
Definition: gridmap.h:76
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:589
Generated on Sun Sep 27 2026 14:12:56 for GLAMERDOC++ by   doxygen 1.8.18