GLAMERDOC++
Gravitational Lensing Code Library
grid_maintenance.h
1 /*
2  * grid_maintenance.h
3  *
4  * Created on: Oct 12, 2011
5  * Author: bmetcalf
6  */
7 
8 #ifndef _grid_maintenance_declare_
9 #define _grid_maintenance_declare_
10 
11 #include "lens.h"
12 #include "point.h"
13 #include "Tree.h"
14 #include "source.h"
15 #include <mutex>
16 #include <utilities_slsim.h>
17 #include "concave_hull.h"
18 
19 class LensHaloBaseNSIE;
20 
24 struct Grid{
25 
26  Grid(LensHndl lens,unsigned long N1d,const double center[2],double range);
27  Grid(LensHndl lens ,unsigned long Nx ,const PosType center[2] ,PosType rangeX ,PosType rangeY);
28  ~Grid();
29 
30  Grid ReInitialize(LensHndl lens);
31  //void ReShoot(LensHndl lens);
32  void zoom(LensHndl lens,double *center,double scale,Branch *top = NULL);
33 
34  //unsigned long PruneTrees(double resolution,bool useSB,double fluxlimit);
35  //unsigned long PrunePointsOutside(double resolution,double *y,double r_in ,double r_out);
36 
37  double RefreshSurfaceBrightnesses(Source* source);
38 
39  double AddSurfaceBrightnesses(Source* source);
40 
41 
42  double mark_closest_point_source_images(
43  Point_2d y_source
44  ,PosType r_source_max
45  ,PosType luminosity
46  ,bool verbose=false
47  );
48  void find_point_source_images(
49  Point_2d y_source
50  ,PosType r_source
51  ,PosType z_source
52  ,std::vector<RAY> &images
53  ,bool verbose=false
54  );
55 
56  double ClearSurfaceBrightnesses();
57 
58  /*** Refine the gris based on the smoothness of the surface brightness.
59  Return new total flux.
60 
61  May be slow.
62  */
63  double refine_on_surfacebrightness(Lens &lens,Source &source);
64 
65  unsigned long getNumberOfPoints() const;
67  PosType EinsteinArea() const;
68 
70  TreeHndl i_tree;
72  TreeHndl s_tree;
73 
75  int getInitNgrid(){return Ngrid_init;}
77  int getNgrid_block(){return Ngrid_block;}
79  double getInitRange(){return i_tree->getTop()->boundary_p2[0] - i_tree->getTop()->boundary_p1[0];}
80  Point_2d getInitCenter();
81  Point * RefineLeaf(LensHndl lens,Point *point);
82  Point * RefineLeaves(LensHndl lens,std::vector<Point *>& points);
83  void ClearAllMarks();
84 
85  //void test_mag_matrix();
86  template <typename T>
87  void writeFits(const double center[],size_t Npixels,double resolution,LensingVariable lensvar,std::string filename);
88  template <typename T>
89  void writeFits(const double center[],size_t Nx,size_t Ny,double resolution,LensingVariable lensvar,std::string filename);
90 
92  template <typename T>
93  void writeFits(
94  double strech
95  ,LensingVariable lensvar
96  ,std::string filename
97  );
98  template <typename T>
99  void writePixelFits(size_t Nx
100  ,LensingVariable lensvar
101  ,std::string filename
102  );
103  template <typename T>
104  void writeFitsVector(const double center[],size_t Npixels,double resolution,LensingVariable lensvar,std::string filename);
105  template <typename T>
106  PixelMap<T> writePixelMap(const double center[],size_t Npixels,double resolution,LensingVariable lensvar);
107  template <typename T>
108  PixelMap<T> writePixelMap(const double center[],size_t Nx,size_t Ny,double resolution,LensingVariable lensvar);
109 
111  template <typename T>
112  PixelMap<T> writePixelMap(LensingVariable lensvar);
113 
115  template<typename T>
116  void MapSurfaceBrightness(PixelMap<T> &map){
117  map.Clean();
118  map.AddGridBrightness(*this);
119  }
121  template <typename T>
122  PixelMap<T> MapSurfaceBrightness(double resolution);
123 
124  template <typename T>
125  PixelMap<T> writePixelMapUniform(const PosType center[],size_t Nx,size_t Ny,LensingVariable lensvar);
126  template<typename T>
127  void writePixelMapUniform(PixelMap<T> &map,LensingVariable lensvar);
128 
129  template <typename T>
130  void writeFitsUniform(const PosType center[],size_t Nx,size_t Ny,LensingVariable lensvar,std::string filename);
131 
132  void find_images(
133  PosType *y_source
134  ,PosType r_source
135  ,int &Nimages
136  ,std::vector<ImageInfo> &imageinfo
137  ,unsigned long &Nimagepoints
138  );
139 
140  void map_images(
141  Lens* lens
142  ,Source *source
143  ,int *Nimages
144  ,std::vector<ImageInfo> &imageinfo
145  ,PosType xmax
146  ,PosType xmin
147  ,PosType initial_size
148  ,ExitCriterion criterion
149  ,bool FindCenter
150  ,bool divide_images
151  );
152 
153  Grid(Grid &&grid){
154  *this = std::move(grid);
155  }
156 
157  Grid operator=(Grid &grid) = delete;
158  Grid(Grid &grid) = delete;
159 
160  Grid & operator=(Grid &&grid){
161  assert(&grid != this);
162 
163  i_tree = grid.i_tree;
164  grid.i_tree = nullptr;
165  s_tree = grid.s_tree;
166  grid.s_tree = nullptr;
167  neighbors = grid.neighbors;
168  grid.neighbors = nullptr;
169  //trashkist = grid.trashkist;
170  //grid.trashkist = nullptr;
171 
172  Ngrid_init = grid.Ngrid_init;
173  Ngrid_init2 = grid.Ngrid_init2;
174  Ngrid_block = grid.Ngrid_block;
175  initialized = grid.initialized;
176  maglimit = grid.maglimit;
177  pointID = grid.pointID;
178  axisratio = grid.axisratio;
179 
180  point_factory.clear();
181  point_factory=std::move(grid.point_factory);
182 
183  return *this;
184  }
185 
187  PosType magnification(double sblimit=-1.0e12) const;
188  PosType UnlensedFlux(double sblimit=-1.0e12) const;
189  PosType LensedFlux(double sblimit=-1.0e12) const;
190 
191  //PosType magnification2() const;
192  //PosType magnification3() const;
194  Point_2d centroid() const;
195 
196  private:
197  void xygridpoints(Point *points,double range,const double *center,long Ngrid
198  ,short remove_center);
199 
201  int Ngrid_init;
202  int Ngrid_init2;
203 
205  int Ngrid_block;
206  bool initialized;
207  //Kist<Point> * trashkist;
208 
209  double maglimit;
210  Kist<Point> * neighbors;
211  bool find_mag_matrix(double *a,Point *p0,Point *p1,Point *p2);
212 
213  bool uniform_mag_from_deflect(double *a,Point *point);
214  //bool uniform_mag_from_shooter(double *a,Point *point);
215 
216  double mag_from_deflect(Point *point) const;
217 
218  unsigned long pointID;
219  PosType axisratio;
220 
221  template <typename T>
222  void writePixelMapUniform_(Point *head,size_t N
223  ,PixelMap<T> *map,LensingVariable val);
224 
225  // cluge to make compatible with old method of producing points
226  Point * NewPointArray(size_t N){
227  Point * p = point_factory(N);
228  p[0].head = N;
229  for(size_t i=1; i < N ; ++i) p[i].head = 0;
230  return p;
231  }
232  MemmoryBank<Point> point_factory;
233  static std::mutex grid_mutex;
234 };
235 
236 typedef struct Grid* GridHndl;
237 
239 std::string to_string(CritType crit);
240 
241 // in image_finder_kist.c
242 namespace ImageFinding{
243 
245  struct CriticalCurve{
246 
247  CriticalCurve(){
248  critical_center[0] = critical_center[1] = 0.0;
249  caustic_center[0] = caustic_center[1] = 0.0;
250  critical_area = 0.0;
251  caustic_area = 0.0;
252  contour_ell = 0.0;
253  ellipse_area = 0.0;
254  z_source = 0.0;
255  type = CritType::ND;
256  caustic_intersections = -1;
257  touches_edge = false;
258  crit_id = -1;
259  };
260 
261  CriticalCurve(const CriticalCurve &p){
262  critcurve = p.critcurve;
263  //critical_curve = p.critical_curve;
264  caustic_curve_outline = p.caustic_curve_outline;
265  caustic_curve_intersecting = p.caustic_curve_intersecting;
266  critical_center = p.critical_center;
267  caustic_center = p.caustic_center;
268  critical_area = p.critical_area;
269  caustic_area = p.caustic_area;
270  ellipse_curve = p.ellipse_curve;
271  contour_ell = p.contour_ell;
272  ellipse_area = p.ellipse_area;
273  z_source = p.z_source;
274  type = p.type;
275  caustic_intersections = p.caustic_intersections;
276  touches_edge = p.touches_edge;
277  crit_id = p.crit_id;
278  }
279 
280  CriticalCurve & operator=(const CriticalCurve &p){
281  if(this == &p) return *this;
282 
283  critcurve = p.critcurve;
284  //critical_curve = p.critical_curve;
285  caustic_curve_outline = p.caustic_curve_outline;
286  caustic_curve_intersecting = p.caustic_curve_intersecting;
287  critical_center = p.critical_center;
288  caustic_center = p.caustic_center;
289  critical_area = p.critical_area;
290  caustic_area = p.caustic_area;
291  ellipse_curve = p.ellipse_curve;
292  contour_ell = p.contour_ell;
293  ellipse_area = p.ellipse_area;
294  z_source = p.z_source;
295  type = p.type;
296  caustic_intersections = p.caustic_intersections;
297  touches_edge = p.touches_edge;
298  crit_id = p.crit_id;
299  return *this;
300  }
301 
302  std::vector<RAY> critcurve;
303  std::vector<Point_2d> caustic_curve_outline;
304  std::vector<Point_2d> caustic_curve_intersecting;
305  std::vector<Point_2d> ellipse_curve;
306 
307  PosType z_source;
309  CritType type;
311  int caustic_intersections;
312 
314  Point_2d critical_center;
316  Point_2d caustic_center;
317 
319  PosType critical_area;
321  PosType caustic_area;
322 
324  PosType contour_ell;
326  PosType ellipse_area;
327 
329  bool touches_edge;
330 
332  int crit_id;
333 
335 
336  bool inCausticCurve(Point_2d &x){
337  return Utilities::inhull(x.x,caustic_curve_outline);
338  }
339 
341  bool intersectingCausticCurve(Point_2d &x,double r){
342  return Utilities::circleIntersetsCurve(x, r, caustic_curve_outline);
343  }
344 
345 
347  bool inCriticalCurve(Point_2d &x){
348  return Utilities::inhull<RAY>(x.x,critcurve);
349  }
350 
352  bool EntirelyinCausticCurve(Point_2d &x, PosType sourceRadius)
353  {
354  // Testing if the center of the source is within the caustic :
355  bool IsInCurve = Utilities::inCurve(x,caustic_curve_outline);
356 
357  // Testing now that it is not too close from it (i.e. farther than source radius) :
358  int i=0; // index going over the different points of the caustic curve
359  PosType DistSourceToCautic; // distance between the source center and the considered point of the caustic line.
360 
361  if(IsInCurve == true) // center inside the caustic
362  {
363  while(i<caustic_curve_outline.size()) // testing that the source size does not overlap with the caustic line.
364  {
365  DistSourceToCautic = sqrt((caustic_curve_outline[i].x[0] - x.x[0])*(caustic_curve_outline[i].x[0] - x.x[0]) + (caustic_curve_outline[i].x[1] - x.x[1])*(caustic_curve_outline[i].x[1] - x.x[1]));
366  if (DistSourceToCautic < sourceRadius) return false ; // source too close, we return false and don't consider the point.
367  i++;
368  }
369  return true ; // if not we return true (the source is valid)
370  }
371  else return false ; // center not inside the caustic
372  }
373 
375  bool EntirelyinCriticalCurve(Point_2d x, PosType sourceRadius)
376  {
377  // Testing if the center of the source is within the critical curve :
378  bool IsInCurve = Utilities::inCurve(x,caustic_curve_outline);
379 
380  // Testing now that it is not too close from it (i.e. farther than source radius) :
381  int i=0; // index going over the different points of the critical curve
382  PosType DistSourceToCritCurve; // distance between the source center and the considered point of the critical line.
383 
384  if(IsInCurve == true) // center inside the critical line
385  {
386  while(i<critcurve.size()) // testing that the source size does not overlap with the critical line.
387  {
388  DistSourceToCritCurve = sqrt((critcurve[i].x[0] - x.x[0])*(critcurve[i].x[0] - x.x[0]) + (critcurve[i].x[1] - x.x[1])*(critcurve[i].x[1] - x.x[1]));
389  if (DistSourceToCritCurve < sourceRadius) return false ; // source too close, we return false and don't consider the point.
390  i++;
391  }
392  return true ; // if not we return true (the source is valid)
393  }
394  else return false ; // center not inside the critical line
395  }
396 
397 
400  void RandomSourcesWithinCaustic(
401  int N
402  ,std::vector<Point_2d> &y
403  ,Utilities::RandomNumbers_NR &rng
404  );
405 
406  Point_2d RandomSourceWithinCaustic(
407  Utilities::RandomNumbers_NR &rng
408  );
409 
412  void RandomSourcesNearCaustic(double R
413  ,int N
414  ,std::vector<Point_2d> &y
415  ,Utilities::RandomNumbers_NR &rng
416  );
419  Point_2d RandomSourceNearCaustic(double R
420  ,Utilities::RandomNumbers_NR &rng
421  );
422 
423 
426  void RandomSourceStrictlyWithinCaustic(int N
427  ,std::vector<Point_2d> &y
428  ,Utilities::RandomNumbers_NR &rng
429  ,PosType sourceRadius
430  ,PosType distSourceToCaustic
431  );
433  void CritRange(Point_2d &p1,Point_2d &p2) const;
435  void CausticRange(Point_2d &p1,Point_2d &p2) const;
436 
438  void CriticalRadius(PosType &rmax,PosType &rmin,PosType &rave) const{
439  if(critcurve.size() < 2){
440  rave = rmax = rmin = 0.0;
441  return;
442  }
443 
444  rave = rmin = rmax = (critical_center - critcurve[0].x).length();
445  PosType rad;
446 
447  for(size_t ii=1; ii< critcurve.size(); ++ii){
448  rad = (critical_center - critcurve[ii].x).length();
449 
450  rave += rad;
451  if(rad < rmin) rmin = rad;
452  if(rad > rmax) rmax = rad;
453  }
454 
455  rave /= critcurve.size();
456  }
458  void CausticRadius(PosType &rmax,PosType &rmin,PosType &rave) const{
459  if(caustic_curve_outline.size() < 2){
460  rave = rmax = rmin = 0.0;
461  return;
462  }
463 
464  rave = rmin = rmax = (caustic_center - caustic_curve_outline[0]).length();
465 
466  PosType rad;
467 
468  for(size_t ii=1; ii< caustic_curve_outline.size(); ++ii){
469  rad = (caustic_center - caustic_curve_outline[ii]).length();
470 
471  rave += rad;
472  if(rad < rmin) rmin = rad;
473  if(rad > rmax) rmax = rad;
474  }
475 
476  rave /= caustic_curve_outline.size();
477  }
478 
480  double AreaNearCaustic(double R
481  );
482 
484  static void print(std::string filename,std::vector<ImageFinding::CriticalCurve> &critcurves) {
485 
486  std::ofstream os(filename.c_str());
487 
488  os <<"# type : "
489  << (int)(CritType::tangential) << " - tangential "
490  << (int)(CritType::radial) << " - radial "
491  << (int)(CritType::pseudo) << " - pseudo"
492  << std::endl;
493  os <<"# plane : "
494  << 0 << " - caustic "
495  << 1 << " - critical curve "
496  << std::endl;
497 
498  os <<"caustic_num,plane,x,y,type" << std::endl;
499 
500  int i=0;
501  Point_2d xo(0,0);
502  for(auto &curve : critcurves){
503 
504  if(curve.type == CritType::tangential){
505  for(Point_2d &p : curve.caustic_curve_outline){
506  os << i << ",0," << p.x[0] - xo[0] <<"," << p.x[1] - xo[1] <<","<< (int)(curve.type) << std::endl;
507  }
508  }else{
509  for(Point_2d &p : curve.caustic_curve_intersecting){
510  os << i << ",0," << p.x[0] - xo[0] <<","<< p.x[1] - xo[1] <<","<< (int)(curve.type) << std::endl;
511  }
512  }
513  for(RAY &p : curve.critcurve){
514  os << i << ",1," << p.x[0] - xo[0] <<","<< p.x[1] - xo[1] <<","<< (int)(curve.type) << std::endl;
515  }
516  ++i;
517  }
518  os.close();
519  }
520 
521  private:
522  Point_2d p1,p2;
523  };
524 
525  void find_images_kist(LensHndl lens,PosType *y_source,PosType r_source,GridHndl grid
526  ,int *Nimages,std::vector<ImageInfo> &imageinfo,unsigned long *Nimagepoints
527  ,PosType initial_size,bool splitimages,short edge_refinement
528  ,bool verbose = false);
529 
530  //void find_image_simple(LensHndl lens,Point_2d y_source,PosType z_source,Point_2d &image_x
531  // ,PosType xtol2,PosType &fret);
532 
533  void find_images_microlens(LensHndl lens,double *y_source,double r_source,GridHndl grid
534  ,int *Nimages,std::vector<ImageInfo> &imageinfo,unsigned long *Nimagepoints
535  ,double initial_size,double mu_min,bool splitimages,short edge_refinement
536  ,bool verbose);
537 
538  void find_images_microlens_exper(LensHndl lens,PosType *y_source,PosType r_source
539  ,GridHndl grid,int *Nimages,std::vector<ImageInfo> &imageinfo,unsigned long *Nimagepoints,PosType initial_size ,PosType mu_min
540  ,bool splitimages,short edge_refinement,bool verbose);
541 
542  void image_finder_kist(LensHndl lens, PosType *y_source,PosType r_source,GridHndl grid
543  ,int *Nimages,std::vector<ImageInfo> &imageinfo,unsigned long *Nimagepoints
544  ,short splitparities,short true_images);
545 
546 
547  void find_crit(LensHndl lens,GridHndl grid,std::vector<CriticalCurve> &crtcurve,int *Ncrits
548  ,double resolution,double invmag_min = 0.0,bool verbose = false,bool test=false);
549  void find_crit(Lens &lens,GridMap &gridmap,std::vector<CriticalCurve> &crtcurves,bool verbose = false);
550 
551  // finds the contours of magnification and source plane curve
552  void find_magnification_contour(
553  Lens &lens
554  ,GridMap &gridmap
555  ,double invmag
556  ,std::vector<std::vector<RAY> > &contour
557  ,std::vector<bool> &hits_boundary
558  );
559 
560  //void find_crit2(LensHndl lens,GridHndl grid,std::vector<CriticalCurve> &critcurve,int *Ncrits
561  // ,double resolution,bool *orderingsuccess,bool ordercurve,bool dividecurves,double invmag_min = 0.0,bool verbose = false);
562 
563  CritType find_pseudo(ImageInfo &pseudocurve,ImageInfo &negimage
564  ,PosType pseudolimit,LensHndl lens,GridHndl grid
565  ,PosType resolution,Kist<Point> &paritypoints,bool TEST=false);
566 
567  void find_contour(LensHndl lens,GridHndl grid,std::vector<CriticalCurve> &contour,int *Ncrits,PosType resolution,bool *orderingsuccess,bool ordercurve, bool dividecurves, double contour_value,LensingVariable contour_type,bool verbose = false);
568 
569  namespace Temporary{
570  //PosType *y;
571  //Lens * lens;
572  //Point *point;
573 
574  PosType mindyfunc(PosType *x);
575  }
576 
577  namespace IF_routines{
578  int refine_grid_kist(LensHndl lens,GridHndl grid,ImageInfo *imageinfo
579  ,int Nimages,double res_target,short criterion
580  ,Kist<Point> * newpointkist = NULL,bool batch=true);
581 
582 
583  void refine_crit_in_image(LensHndl lens,GridHndl grid,double r_source,double x_source[],double resolution);
584 
585  int refine_grid(LensHndl lens,GridHndl grid,OldImageInfo *imageinfo
586  ,unsigned long Nimages,double res_target,short criterion,bool batch=true);
587 
588  long refine_edges(LensHndl lens,GridHndl grid,ImageInfo *imageinfo
589  ,int Nimages,double res_target,short criterion
590  ,Kist<Point> * newpointkist = NULL,bool batch=true);
591 
592  long refine_edges2(LensHndl lens,double *y_source,double r_source,GridHndl grid
593  ,ImageInfo *imageinfo,bool *image_overlap,int Nimages,double res_target
594  ,short criterion,bool batch=true);
595 
596  void sort_out_points(Point *i_points,ImageInfo *imageinfo,double r_source,double y_source[]);
597 
598  }
599 
600  void printCriticalCurves(std::string filename
601  ,const std::vector<ImageFinding::CriticalCurve> &critcurves);
602 
606  template<typename T>
607  PixelMap<T> mapCriticalCurves(
609  const std::vector<ImageFinding::CriticalCurve> &critcurves,
611  int Nx
612  );
613 
617  template<typename T>
618  PixelMap<T> mapCausticCurves(const std::vector<ImageFinding::CriticalCurve> &critcurves
619  ,int Nx
620  );
621 }
622 
623 
624 std::ostream &operator<<(std::ostream &os, const ImageFinding::CriticalCurve &p);
625 
627 template <typename T>
628 void Grid::writeFits(
629  const PosType center[]
630  ,size_t Npixels
631  ,PosType resolution
632  ,LensingVariable lensvar
633  ,std::string filename
634  ){
635  writeFits<T>(center,Npixels,Npixels,resolution,lensvar,filename);
636 }
638 template <typename T>
639 void Grid::writeFits(
640  const PosType center[]
641  ,size_t Nx
642  ,size_t Ny
643  ,PosType resolution
644  ,LensingVariable lensvar
645  ,std::string filename
646  ){
647  PixelMap<T> map(center, Nx,Ny, resolution);
648  std::string tag;
649 
650  switch (lensvar) {
651  case LensingVariable::DELAYT:
652  tag = ".dt.fits";
653  break;
654  case LensingVariable::ALPHA1:
655  tag = ".alpha1.fits";
656  break;
657  case LensingVariable::ALPHA2:
658  tag = ".alpha2.fits";
659  break;
660  case LensingVariable::ALPHA:
661  tag = ".alpha.fits";
662  break;
663  case LensingVariable::KAPPA:
664  tag = ".kappa.fits";
665  break;
666  case LensingVariable::GAMMA1:
667  tag = ".gamma1.fits";
668  break;
669  case LensingVariable::GAMMA2:
670  tag = ".gamma2.fits";
671  break;
672  case LensingVariable::GAMMA3:
673  tag = ".gamma3.fits";
674  break;
675  case LensingVariable::GAMMA:
676  tag = ".gamma.fits";
677  break;
678  case LensingVariable::INVMAG:
679  tag = ".invmag.fits";
680  break;
681  case LensingVariable::SurfBrightness:
682  tag = ".surfbright.fits";
683  break;
684  default:
685  break;
686  }
687 
688  map.AddGrid(*this,lensvar);
689  map.printFITS(filename + tag);
690 
691  return;
692 }
693 
695 template <typename T>
696 PixelMap<T> Grid::writePixelMap(
697  const PosType center[]
698  ,size_t Npixels
699  ,PosType resolution
700  ,LensingVariable lensvar
701  ){
702 
703  return writePixelMap<T>(center,Npixels,Npixels,resolution,lensvar);
704 }
706 template <typename T>
707 PixelMap<T> Grid::writePixelMap(
708  const PosType center[]
709  ,size_t Nx
710  ,size_t Ny
711  ,PosType resolution
712  ,LensingVariable lensvar
713 ){
714  PixelMap<T> map(center, Nx, Ny, resolution);
715  map.AddGrid(*this,lensvar);
716 
717  return map;
718 }
720 template <typename T>
721 PixelMap<T> Grid::writePixelMap(
722  LensingVariable lensvar
723 ){
724 
725  Branch *branch = i_tree->getTop();
726  double resolution = (branch->boundary_p2[0] - branch->boundary_p1[0])/Ngrid_init;
727  PixelMap<T> map(branch->center, Ngrid_init, Ngrid_init2, resolution);
728  map.AddGrid(*this,lensvar);
729 
730  return map;
731 }
732 
733 template <typename T>
734 PixelMap<T> Grid::MapSurfaceBrightness(double resolution){
735  Branch *branch = i_tree->getTop();
736  int Nx = (int)( (branch->boundary_p2[0] - branch->boundary_p1[0])/resolution );
737  int Ny = (int)( (branch->boundary_p2[1] - branch->boundary_p1[1])/resolution );
738 
739  PixelMap<T> map(branch->center,Nx,Ny,resolution);
740  map.AddGridBrightness(*this);
741 
742  return map;
743 }
744 
747 template <typename T>
748 void Grid::writePixelFits(
749  size_t Nx
750  ,LensingVariable lensvar
751  ,std::string filename
752 ){
753 
754  Point_2d center = getInitCenter();
755  PosType resolution = getInitRange()/Nx;
756  size_t Ny = (size_t)(Nx*axisratio);
757  writeFits<T>(center.x,Nx,Ny,resolution, lensvar, filename);
758 
759  return;
760 }
761 
769 template <typename T>
770 void Grid::writeFitsUniform(
771  const PosType center[]
772  ,size_t Nx
773  ,size_t Ny
774  ,LensingVariable lensvar
775  ,std::string filename
776  ){
777  std::string tag;
778 
779  switch (lensvar) {
780  case LensingVariable::DELAYT:
781  tag = ".dt.fits";
782  break;
783  case LensingVariable::ALPHA1:
784  tag = ".alpha1.fits";
785  break;
786  case LensingVariable::ALPHA2:
787  tag = ".alpha2.fits";
788  break;
789  case LensingVariable::ALPHA:
790  tag = ".alpha.fits";
791  break;
792  case LensingVariable::KAPPA:
793  tag = ".kappa.fits";
794  break;
795  case LensingVariable::GAMMA1:
796  tag = ".gamma1.fits";
797  break;
798  case LensingVariable::GAMMA2:
799  tag = ".gamma2.fits";
800  break;
801  case LensingVariable::GAMMA3:
802  tag = ".gamma3.fits";
803  break;
804  case LensingVariable::GAMMA:
805  tag = ".gamma.fits";
806  break;
807  case LensingVariable::INVMAG:
808  tag = ".invmag.fits";
809  break;
810  case LensingVariable::SurfBrightness:
811  tag = ".surfbright.fits";
812  break;
813  default:
814  break;
815  }
816 
817  PixelMap<T> map = writePixelMapUniform<T>(center,Nx,Ny,lensvar);
818  map.printFITS(filename + tag);
819 }
820 
828 template <typename T>
829 PixelMap<T> Grid::writePixelMapUniform(
830  const PosType center[]
831  ,size_t Nx
832  ,size_t Ny
833  ,LensingVariable lensvar
834  ){
835 
836  if(getNumberOfPoints() ==0 ) return PixelMap<T>();
837  PixelMap<T> map(center, Nx, Ny,i_tree->pointlist.Top()->gridsize);
838  map.Clean();
839 
840  //int Nblocks = Utilities::GetNThreads();
841  int Nblocks = 16;
842 
843  //std::vector<PointList> lists(Nblocks);
844 
845  std::vector<Point *> heads(Nblocks);
846  std::vector<size_t> sizes(Nblocks,0);
847 
848  bool allowDecent;
849  TreeStruct::iterator i_tree_it(i_tree);
850  int i = 0;
851 
852  do{
853  if((*i_tree_it)->level == 4){
854 
855  heads[i] = (*i_tree_it)->points;
856  sizes[i] = (*i_tree_it)->npoints;
857 
858  //lists[i].setTop( (*i_tree_it)->points );
859  //lists[i].setN( (*i_tree_it)->npoints );
860  ++i;
861  allowDecent = false;
862  }else{
863  allowDecent = true;
864  }
865  }while(i_tree_it.TreeWalkStep(allowDecent) && i < Nblocks);
866 
867  std::vector<std::thread> thrs;
868  for(int ii = 0; ii < i ;++ii){
869  //writePixelMapUniform_(heads[ii],sizes[ii],&map,lensvar);
870  //thrs.push_back(std::thread(&Grid::writePixelMapUniform_,this,lists[ii],&map,lensvar));
871 
872  thrs.push_back(std::thread(&Grid::writePixelMapUniform_<T>,this,heads[ii],sizes[ii],&map,lensvar));
873  }
874  for(int ii = 0; ii < i ;++ii) thrs[ii].join();
875 
876  return map;
877 }
878 
879 template <typename T>
880 void Grid::writePixelMapUniform(
881  PixelMap<T> &map
882  ,LensingVariable lensvar
883  ){
884 
885  if(getNumberOfPoints() ==0 ) return;
886 
887  map.Clean();
888  int Nblocks = 16;
889  //std::vector<PointList> lists(Nblocks);
890  TreeStruct::iterator i_tree_it(i_tree);
891 
892  std::vector<Point *> heads(Nblocks);
893  std::vector<size_t> sizes(Nblocks,0);
894 
895 
896  bool allowDecent;
897  i_tree_it.movetop();
898  int i = 0;
899  do{
900  if((*i_tree_it)->level == 4){
901  assert(i < 16);
902  //lists[i].setTop( (*i_tree_it)->points );
903  //lists[i].setN( (*i_tree_it)->npoints );
904  heads[i] = (*i_tree_it)->points;
905  sizes[i] = (*i_tree_it)->npoints;
906 
907  ++i;
908  allowDecent = false;
909  }else{
910  allowDecent = true;
911  }
912  }while(i_tree_it.TreeWalkStep(allowDecent) && i < Nblocks);
913 
914  std::vector<std::thread> thr;
915  for(int ii = 0; ii < i ;++ii){
916  thr.push_back(std::thread(&Grid::writePixelMapUniform_<T>,this,heads[ii],sizes[ii],&map,lensvar));
917  }
918  for(auto &t : thr) t.join();
919 }
920 
921 template <typename T>
922 void Grid::writePixelMapUniform_(Point *head,size_t N,PixelMap<T> *map,LensingVariable val){
923  double tmp;
924  PosType tmp2[2];
925  long index;
926 
927  Point *ppoint = head;
928 
929  for(size_t i = 0; i< N; ++i){
930 
931  switch (val) {
932  case LensingVariable::ALPHA:
933  tmp2[0] = ppoint->x[0] - ppoint->image->x[0];
934  tmp2[1] = ppoint->x[1] - ppoint->image->x[1];
935  tmp = sqrt(tmp2[0]*tmp2[0] + tmp2[1]*tmp2[1]);
936  break;
937  case LensingVariable::ALPHA1:
938  tmp = (ppoint->x[0] - ppoint->image->x[0]);
939  break;
940  case LensingVariable::ALPHA2:
941  tmp = (ppoint->x[1] - ppoint->image->x[1]);
942  break;
943  case LensingVariable::KAPPA:
944  tmp = ppoint->kappa();
945  break;
946  case LensingVariable::GAMMA:
947  tmp2[0] = ppoint->gamma1();
948  tmp2[1] = ppoint->gamma2();
949  tmp = sqrt(tmp2[0]*tmp2[0] + tmp2[1]*tmp2[1]);
950  break;
951  case LensingVariable::GAMMA1:
952  tmp = ppoint->gamma1();
953  break;
954  case LensingVariable::GAMMA2:
955  tmp = ppoint->gamma2();
956  break;
957  case LensingVariable::GAMMA3:
958  tmp = ppoint->gamma3();
959  break;
960  case LensingVariable::INVMAG:
961  tmp = ppoint->invmag();
962  break;
963  case LensingVariable::DELAYT:
964  tmp = ppoint->dt;
965  break;
966  case LensingVariable::SurfBrightness:
967  tmp = ppoint->surface_brightness;
968  break;
969  default:
970  std::cerr << "PixelMap<T>::AddGrid() does not work for the input LensingVariable" << std::endl;
971  throw std::runtime_error("PixelMap<T>::AddGrid() does not work for the input LensingVariable");
972  break;
973  // If this list is to be expanded to include ALPHA or GAMMA take care to add them as vectors
974  }
975 
976  index = map->find_index(ppoint->x);
977  if(index != -1)(*map)[index] = tmp;
978 
979  ppoint = ppoint->next;
980  }
981 }
982 
984 template <typename T>
985 void Grid::writeFitsVector(
986  const PosType center[]
987  ,size_t Npixels
988  ,PosType resolution
989  ,LensingVariable lensvar
990  ,std::string filename
991  ){
992  //throw std::runtime_error("Not done yet!");
993 
994  PosType range = Npixels*resolution,tmp_x[2];
995  ImageInfo tmp_image,tmp_image_theta;
996  size_t i;
997  std::string tag;
998 
999  i_tree->PointsWithinKist(center,range/sqrt(2.),tmp_image.imagekist,0);
1000  i_tree->PointsWithinKist(center,range/sqrt(2.),tmp_image_theta.imagekist,0);
1001 
1002  std::vector<PosType> tmp_sb_vec(tmp_image.imagekist->Nunits());
1003 
1004  for(tmp_image.imagekist->MoveToTop(),i=0;i<tmp_sb_vec.size();++i,tmp_image.imagekist->Down()){
1005  tmp_sb_vec[i] = tmp_image.imagekist->getCurrent()->surface_brightness;
1006  switch (lensvar) {
1007  case LensingVariable::ALPHA1:
1008  tmp_x[0] = tmp_image.imagekist->getCurrent()->x[0]
1009  - tmp_image.imagekist->getCurrent()->image->x[0];
1010 
1011  tmp_x[1] = tmp_image.imagekist->getCurrent()->x[1]
1012  - tmp_image.imagekist->getCurrent()->image->x[1];
1013 
1014  tmp_image.imagekist->getCurrent()->surface_brightness = sqrt( tmp_x[0]*tmp_x[0] + tmp_x[1]*tmp_x[1]);
1015  tmp_image_theta.imagekist->getCurrent()->surface_brightness = atan2(tmp_x[1],tmp_x[0]);
1016 
1017  tag = ".alphaV.fits";
1018  break;
1019  case LensingVariable::GAMMA:
1020 
1021  tmp_x[0] = tmp_image.imagekist->getCurrent()->gamma1();
1022  tmp_x[1] = tmp_image.imagekist->getCurrent()->gamma2();
1023 
1024  tmp_image.imagekist->getCurrent()->surface_brightness = sqrt( tmp_x[0]*tmp_x[0] + tmp_x[1]*tmp_x[1]);
1025  tmp_image_theta.imagekist->getCurrent()->surface_brightness = atan2(tmp_x[1],tmp_x[0])/2;
1026 
1027  tag = ".gammaV.fits";
1028  break;
1029  default:
1030  std::cout << "Grid::writeFitsVector() does not support the LensVariable you are using." << std::endl;
1031  return;
1032  }
1033  }
1034 
1035  PixelMap<T> map_m(center, Npixels, resolution),map_t(center,Npixels,resolution);
1036 
1037  map_m.Clean();
1038  map_m.AddImages(&tmp_image,1,-1);
1039  map_m = PixelMap<T>(map_m,4);
1040  map_m = PixelMap<T>(map_m,1/4.);
1041 
1042  map_t.Clean();
1043  map_t.AddImages(&tmp_image_theta,1,-1);
1044 
1045  map_m.printFITS(filename + tag);
1046 
1047  for(tmp_image.imagekist->MoveToTop(),i=0;i<tmp_sb_vec.size();++i,tmp_image.imagekist->Down())
1048  tmp_image.imagekist->getCurrent()->surface_brightness = tmp_sb_vec[i];
1049 }
1050 
1051 template <typename T>
1052 void Grid::writeFits(double strech,LensingVariable lensvar ,std::string filename){
1053  Point_2d center = getInitCenter();
1054  size_t N1 = (size_t)(strech*Ngrid_init);
1055  size_t N2 = (size_t)(strech*Ngrid_init2);
1056 
1057  writeFits<T>(center.x,N1,N2,getInitRange()/N1,lensvar,filename);
1058 }
1059 
1060 
1061 
1062 #endif
Grid::RefineLeaf
Point * RefineLeaf(LensHndl lens, Point *point)
Fundamental function used to divide a leaf in the tree into nine subcells.
Definition: grid_maintenance.cpp:543
ImageFinding::find_images_microlens_exper
void find_images_microlens_exper(LensHndl lens, PosType *y_source, PosType r_source, GridHndl grid, int *Nimages, std::vector< ImageInfo > &imageinfo, unsigned long *Nimagepoints, PosType initial_size, PosType mu_min, bool splitimages, short edge_refinement, bool verbose)
experimental version of find_image_microlens()
Definition: image_finder_kist.cpp:1027
to_string
std::string to_string(const Band &band)
convert a Band type to a string name.
Definition: InputParams.cpp:878
Utilities::circleIntersetsCurve
bool circleIntersetsCurve(const Point_2d &x, double r, const std::vector< Point_2d > &v)
returns true if a circle of radius r around the point x intersects with the curve v....
Definition: curve_routines.cpp:2759
Grid::MapSurfaceBrightness
void MapSurfaceBrightness(PixelMap< T > &map)
make image of surface brightness
Definition: grid_maintenance.h:116
Grid::writePixelMap
PixelMap< T > writePixelMap(const double center[], size_t Npixels, double resolution, LensingVariable lensvar)
Outputs a PixelMap of the lensing quantities of a fixed grid.
Definition: grid_maintenance.h:696
ImageFinding::CriticalCurve::crit_id
int crit_id
id of the critical curve
Definition: grid_maintenance.h:332
Grid::getInitRange
double getInitRange()
return initial range of gridded region
Definition: grid_maintenance.h:79
ImageFinding::CriticalCurve::CausticRadius
void CausticRadius(PosType &rmax, PosType &rmin, PosType &rave) const
find 3 measures of the caustic curve radius
Definition: grid_maintenance.h:458
Grid::RefineLeaves
Point * RefineLeaves(LensHndl lens, std::vector< Point * > &points)
Same as RefineLeaf() but multiple points can be passed. The rays are shot all together so that more p...
Definition: grid_maintenance.cpp:681
LensingVariable::GAMMA3
@ GAMMA3
third component of shear
ImageFinding::CriticalCurve::critical_area
PosType critical_area
area of critical curve (radians^2)
Definition: grid_maintenance.h:319
ImageFinding
The ImageFinding namespace is for functions related to finding and mapping images.
Definition: grid_maintenance.h:242
ImageFinding::image_finder_kist
void image_finder_kist(LensHndl lens, PosType *y_source, PosType r_source, GridHndl grid, int *Nimages, std::vector< ImageInfo > &imageinfo, unsigned long *Nimagepoints, short splitparities, short true_images)
Finds images for a given source position and size. Not meant for high level user.
Definition: image_finder_kist.cpp:1530
ImageFinding::CriticalCurve::inCriticalCurve
bool inCriticalCurve(Point_2d &x)
return true if x is inside or on the border of the critical curve
Definition: grid_maintenance.h:347
ImageInfo::imagekist
Kist< Point > * imagekist
Array of points in image, SHOULD NOT BE USED IN FAVOR OF imagekist! Still used by caustic finding rou...
Definition: image_info.h:47
Grid::getInitNgrid
int getInitNgrid()
return initial number of grid points in each direction
Definition: grid_maintenance.h:75
Source
Base class for all sources.
Definition: source.h:44
LensingVariable::ALPHA2
@ ALPHA2
y component of deflection
ImageFinding::find_images_microlens
void find_images_microlens(LensHndl lens, double *y_source, double r_source, GridHndl grid, int *Nimages, std::vector< ImageInfo > &imageinfo, unsigned long *Nimagepoints, double initial_size, double mu_min, bool splitimages, short edge_refinement, bool verbose)
Finds images given a source position and size.
Definition: image_finder_kist.cpp:465
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
Grid
Structure to contain both source and image trees.
Definition: grid_maintenance.h:24
Grid::AddSurfaceBrightnesses
double AddSurfaceBrightnesses(Source *source)
Recalculate surface brightness just like Grid::RefreshSurfaceBrightness but the new source is added t...
Definition: grid_maintenance.cpp:344
Grid::ReInitialize
Grid ReInitialize(LensHndl lens)
Returns a new grid that has not been refined but has the same intial image grid, but calculated with ...
Definition: grid_maintenance.cpp:148
Branch
The box representing a branch of a binary tree structure. Used specifically in TreeStruct for organiz...
Definition: point.h:677
TreeStruct::pointlist
PointList pointlist
list of points
Definition: Tree.h:125
LensingVariable::GAMMA
@ GAMMA
magnitude of shear
ImageFinding::CriticalCurve::touches_edge
bool touches_edge
touches the edge of the gridded region
Definition: grid_maintenance.h:329
Grid::getNumberOfPoints
unsigned long getNumberOfPoints() const
Returns number of points on image plane.
Definition: grid_maintenance.cpp:518
ImageFinding::CriticalCurve::critical_center
Point_2d critical_center
center of critical curve
Definition: grid_maintenance.h:314
ImageFinding::CriticalCurve
Structure to contain information and functions for critical curves and caustics.
Definition: grid_maintenance.h:245
Grid::ClearAllMarks
void ClearAllMarks()
Rest all in_image markers to False.
Definition: grid_maintenance.cpp:984
ImageFinding::CriticalCurve::CritRange
void CritRange(Point_2d &p1, Point_2d &p2) const
find rectangular region enclosing critical curve
Definition: Tree.cpp:972
LensingVariable::KAPPA
@ KAPPA
convergence
ImageFinding::CriticalCurve::RandomSourceNearCaustic
Point_2d RandomSourceNearCaustic(double R, Utilities::RandomNumbers_NR &rng)
Returns a random point within distance R or filly within the caustic.
Definition: Tree.cpp:930
PixelMap::AddGrid
void AddGrid(const Grid &grid, T value=1.0)
Fills in pixels where the image plane points in the grid are located with the value given.
Definition: pixelmap.cpp:2231
ImageFinding::CriticalCurve::RandomSourceWithinCaustic
Point_2d RandomSourceWithinCaustic(Utilities::RandomNumbers_NR &rng)
Definition: Tree.cpp:918
ImageFinding::find_images_kist
void find_images_kist(LensHndl lens, PosType *y_source, PosType r_source, GridHndl grid, int *Nimages, std::vector< ImageInfo > &imageinfo, unsigned long *Nimagepoints, PosType initial_size, bool splitimages, short edge_refinement, bool verbose=false)
Finds images given a source position and size.
Definition: image_finder_kist.cpp:41
Point_2d
Class for representing points or vectors in 2 dimensions. Not that the dereferencing operator is over...
Definition: point.h:48
Grid::zoom
void zoom(LensHndl lens, double *center, double scale, Branch *top=NULL)
Test if point is in a region of uniform magnification using the kappa and gamma calculated from the r...
Definition: grid_maintenance.cpp:1233
Utilities::inhull< RAY >
bool inhull< RAY >(PosType x[], const std::vector< RAY > &H)
finds in x is within the curve discribed by the H[].x points ie image points
Definition: concave_hull.h:1164
ImageFinding::CriticalCurve::RandomSourcesNearCaustic
void RandomSourcesNearCaustic(double R, int N, std::vector< Point_2d > &y, Utilities::RandomNumbers_NR &rng)
Returns a vector of random point within distance R or filly within the caustic.
Definition: Tree.cpp:922
Utilities::RandomNumbers_NR
This is a class for generating random numbers. It simplifies and fool proofs initialization and allow...
Definition: utilities_slsim.h:1068
LensingVariable::ALPHA
@ ALPHA
magnitude of deflection in radians
ImageFinding::CriticalCurve::caustic_intersections
int caustic_intersections
estimated number of intersections of the caustic, -1 if not set
Definition: grid_maintenance.h:311
PixelMap::AddGridBrightness
void AddGridBrightness(Grid &grid)
Add an image from a the surface brightnesses of a Grid to the PixelMap.
Definition: pixelmap.cpp:685
GridMap
A simplified version of the Grid structure for making non-adaptive maps of the lensing quantities (ka...
Definition: gridmap.h:31
Grid::centroid
Point_2d centroid() const
centroid of flux
Definition: grid_maintenance.cpp:456
Grid::find_point_source_images
void find_point_source_images(Point_2d y_source, PosType r_source, PosType z_source, std::vector< RAY > &images, bool verbose=false)
This function finds all the images for a circular source of radius r_source, then finds the points wi...
Definition: image_finder.cpp:1094
LensingVariable::SurfBrightness
@ SurfBrightness
Surface brightness.
ImageFinding::find_contour
void find_contour(LensHndl lens, GridHndl grid, std::vector< CriticalCurve > &contour, int *Ncrits, PosType resolution, bool *orderingsuccess, bool ordercurve, bool dividecurves, double contour_value, LensingVariable contour_type, bool verbose=false)
Finds iso kappa contours.
Definition: find_crit.cpp:1897
TreeStruct::iterator
A iterator class fore TreeStruct that allows for movement through the tree without changing anything ...
Definition: Tree.h:47
TreeStruct::iterator::TreeWalkStep
bool TreeWalkStep(bool allowDescent)
step for walking tree by iteration instead of recursion.
Definition: tree_maintenance.cpp:1610
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
TreeStruct
Tree: Exported struct.
Definition: Tree.h:31
TreeStruct::getTop
Branch * getTop()
root branch
Definition: Tree.h:122
ImageFinding::CriticalCurve::CausticRange
void CausticRange(Point_2d &p1, Point_2d &p2) const
find rectangular region enclosing caustic curve
Definition: Tree.cpp:987
ImageFinding::CriticalCurve::ellipse_area
PosType ellipse_area
area of an ellipse with axis ratio contour_ell and major axis = max distance between center (as given...
Definition: grid_maintenance.h:326
ImageFinding::CriticalCurve::RandomSourcesWithinCaustic
void RandomSourcesWithinCaustic(int N, std::vector< Point_2d > &y, Utilities::RandomNumbers_NR &rng)
Returns a vector of random point within the caustic. It is more efficient to call this once for many ...
Definition: Tree.cpp:904
ImageFinding::CriticalCurve::intersectingCausticCurve
bool intersectingCausticCurve(Point_2d &x, double r)
returns true if a circle of radius r around the point x intersects with the caustic curve
Definition: grid_maintenance.h:341
LensingVariable::GAMMA2
@ GAMMA2
second component of shear
ImageFinding::CriticalCurve::caustic_area
PosType caustic_area
area of caustic curve (radians^2)
Definition: grid_maintenance.h:321
Grid::writeFitsVector
void writeFitsVector(const double center[], size_t Npixels, double resolution, LensingVariable lensvar, std::string filename)
Outputs a fits file for making plots of vector fields.
Definition: grid_maintenance.h:985
RAY
Simple representaion of a light path giving position on the image and source planes and lensing quant...
Definition: point.h:555
PixelMap::AddImages
void AddImages(ImageInfo *imageinfo, int Nimages, float rescale=1.)
Add an image to the map.
Definition: pixelmap.cpp:641
Grid::ClearSurfaceBrightnesses
double ClearSurfaceBrightnesses()
Reset the surface brightness and in_image flag in every point on image and source planes to zero (fal...
Definition: grid_maintenance.cpp:500
Utilities::inCurve
bool inCurve(const Ptype &x, const std::vector< Ptype > &H)
returns true if x is within the polygon H
Definition: concave_hull.h:1054
ImageFinding::CriticalCurve::AreaNearCaustic
double AreaNearCaustic(double R)
returns an estimate of the area inside and within distance R of the caustic
Definition: Tree.cpp:936
LensingVariable::DELAYT
@ DELAYT
time delay
ImageFinding::CriticalCurve::inCausticCurve
bool inCausticCurve(Point_2d &x)
return true if x is inside or on the border of the caustic curve
Definition: grid_maintenance.h:336
Grid::i_tree
TreeHndl i_tree
tree on image plane
Definition: grid_maintenance.h:70
ImageFinding::find_crit
void find_crit(LensHndl lens, GridHndl grid, std::vector< CriticalCurve > &crtcurve, int *Ncrits, double resolution, double invmag_min=0.0, bool verbose=false, bool test=false)
Finds critical curves and caustics.
Definition: find_crit.cpp:37
Grid::writePixelMapUniform
PixelMap< T > writePixelMapUniform(const PosType center[], size_t Nx, size_t Ny, LensingVariable lensvar)
Make a Pixel map of the without distribution the pixels.
Definition: grid_maintenance.h:829
Grid::writeFits
void writeFits(const double center[], size_t Npixels, double resolution, LensingVariable lensvar, std::string filename)
Outputs a fits image of a lensing variable of choice.
Definition: grid_maintenance.h:628
LensingVariable::GAMMA1
@ GAMMA1
first component of shear
ImageFinding::find_pseudo
CritType find_pseudo(ImageInfo &pseudocurve, ImageInfo &negimage, PosType pseudolimit, LensHndl lens, GridHndl grid, PosType resolution, Kist< Point > &paritypoints, bool TEST=false)
Definition: find_crit.cpp:1107
Grid::Grid
Grid(LensHndl lens, unsigned long N1d, const double center[2], double range)
Constructor for initializing square grid.
Definition: grid_maintenance.cpp:22
ImageFinding::mapCriticalCurves
PixelMap< T > mapCriticalCurves(const std::vector< ImageFinding::CriticalCurve > &critcurves, int Nx)
Makes an image of the critical curves. The map will encompose all curves found. The pixel values are ...
Definition: Tree.cpp:1026
LensingVariable::INVMAG
@ INVMAG
inverse of magnification
Grid::getNgrid_block
int getNgrid_block()
return number of cells in each dimension into which each cell is divided when a refinement is made
Definition: grid_maintenance.h:77
ImageFinding::mapCausticCurves
PixelMap< T > mapCausticCurves(const std::vector< ImageFinding::CriticalCurve > &critcurves, int Nx)
Makes an image of the caustic curves. The map will encompose all curves found. The pixel values are t...
Definition: Tree.cpp:1064
ImageFinding::CriticalCurve::type
CritType type
type of caustic, 0 -not defined, 1 -radial, 2 - tangential,3 - pseudo
Definition: grid_maintenance.h:309
OldImageInfo
This is an old version that should not be used anymore in favor of ImageInfo.
Definition: image_info.h:109
ImageFinding::CriticalCurve::EntirelyinCriticalCurve
bool EntirelyinCriticalCurve(Point_2d x, PosType sourceRadius)
return true if x is strictly inside (entirely) the critical curve
Definition: grid_maintenance.h:375
Grid::RefreshSurfaceBrightnesses
double RefreshSurfaceBrightnesses(Source *source)
Reshoot the rays with the same image postions.
Definition: grid_maintenance.cpp:315
Grid::find_images
void find_images(PosType *y_source, PosType r_source, int &Nimages, std::vector< ImageInfo > &imageinfo, unsigned long &Nimagepoints)
Finds images for a given source position and size. It seporates images of different pairities.
Definition: image_finder.cpp:951
ImageFinding::CriticalCurve::EntirelyinCausticCurve
bool EntirelyinCausticCurve(Point_2d &x, PosType sourceRadius)
return true if x is strictly inside (entirely) the caustic curve
Definition: grid_maintenance.h:352
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
TreeStruct::PointsWithinKist
PosType PointsWithinKist(const PosType *center, PosType rmax, Kist< Point > *neighborkist, short markpoints) const
Finds all points in tree that lie within rmax of the point ray[].
Definition: KistDriver.cpp:263
Grid::magnification
PosType magnification(double sblimit=-1.0e12) const
flux weighted local magnification that does not take multiple imaging into effect
Definition: grid_maintenance.cpp:398
LensingVariable
LensingVariable
output lensing variables
Definition: standard.h:89
ImageFinding::CriticalCurve::caustic_center
Point_2d caustic_center
center of caustic curve
Definition: grid_maintenance.h:316
Grid::EinsteinArea
PosType EinsteinArea() const
area of region with negative magnification
Definition: grid_maintenance.cpp:368
ImageFinding::CriticalCurve::RandomSourceStrictlyWithinCaustic
void RandomSourceStrictlyWithinCaustic(int N, std::vector< Point_2d > &y, Utilities::RandomNumbers_NR &rng, PosType sourceRadius, PosType distSourceToCaustic)
Returns a vector of random point strictly within the caustic (i.e. not touching the border)....
Definition: Tree.cpp:958
Grid::writePixelFits
void writePixelFits(size_t Nx, LensingVariable lensvar, std::string filename)
Make a fits map that is automatically centered on the grid and has approximately the same range as th...
Definition: grid_maintenance.h:748
Grid::writeFitsUniform
void writeFitsUniform(const PosType center[], size_t Nx, size_t Ny, LensingVariable lensvar, std::string filename)
Output a fits map of the without distribution the pixels.
Definition: grid_maintenance.h:770
Grid::map_images
void map_images(Lens *lens, Source *source, int *Nimages, std::vector< ImageInfo > &imageinfo, PosType xmax, PosType xmin, PosType initial_size, ExitCriterion criterion, bool FindCenter, bool divide_images)
Find images and refine them based on their surface brightness distribution.
Definition: map_images.cpp:27
ImageInfo
Structure for storing information about images or curves.
Definition: image_info.h:20
ImageFinding::CriticalCurve::CriticalRadius
void CriticalRadius(PosType &rmax, PosType &rmin, PosType &rave) const
find 3 measures of the critical curve radius
Definition: grid_maintenance.h:438
ImageFinding::CriticalCurve::print
static void print(std::string filename, std::vector< ImageFinding::CriticalCurve > &critcurves)
print the critical curves and caustic curves to a csv file
Definition: grid_maintenance.h:484
Grid::mark_closest_point_source_images
double mark_closest_point_source_images(Point_2d y_source, PosType r_source_max, PosType luminosity, bool verbose=false)
This function finds all the images for a circular source of radius r_source, then finds the points wi...
Definition: image_finder.cpp:1061
Grid::~Grid
~Grid()
Destructor for a Grid. Frees all memory.
Definition: grid_maintenance.cpp:136
RAY::x
Point_2d x
image position
Definition: point.h:625
MemmoryBank< Point >
Lens
A class to represents a lens with multiple planes.
Definition: lens.h:71
Grid::s_tree
TreeHndl s_tree
tree on source plane
Definition: grid_maintenance.h:72
ImageFinding::CriticalCurve::contour_ell
PosType contour_ell
axis ratio of a contour defined by the ratio of the max to min distance between center (as given by h...
Definition: grid_maintenance.h:324
Generated on Sun Sep 27 2026 14:12:56 for GLAMERDOC++ by   doxygen 1.8.18