analytic_lens.h

/*
 * nsie.h
 *
 *  Created on: Dec 8, 2009
 *      Author: R.B. Metcalf
 */
 
#ifndef analens_declare
#define analens_declare
 
#include "Tree.h"
#include "source.h"
#include "base_analens.h"
 
// TODO: BEN finish this documentation for perturbation parameters.
class LensHaloFit : public LensHaloBaseNSIE{
  
  void abstractfunction(){}; // pure virtual
  
public:
  //LensHaloAnaNSIE(InputParams& params,bool verbose = false);
  LensHaloFit(const COSMOLOGY& cosmo, int MyNmodes, PosType beta,PosType zlensref, PosType zsourceref);
 ~LensHaloFit();
 
  void PrintLens(bool show_substruct,bool show_stars);
 
  void FindLensSimple(int Nimages,Point *image_positions,double *y,double **dx_sub);
  bool SafeFindLensSimple(int Nimages,Point *image_positions,double *y,double **dx_sub,int SafetyNum,PosType PixelSizeRad,std::vector<std::vector<PosType>> & PrecisionBackTracedPos,std::vector<std::vector<PosType>> & alphaTab,bool verbose);
  void FindLensSimple(ImageInfo *imageinfo ,int Nimages ,double *y,double **dx_sub);
  
  // these need to be written so that they translate between modes and these quantities
  virtual PosType get_sigma(){return 0.0;};
  virtual PosType get_fratio(){return 0.0;};
  virtual PosType get_pa(){return 0.0;};
  virtual PosType get_rcore(){return 0.0;};
  
 
  void setNmodes(int my_Nmodes){perturb_Nmodes = my_Nmodes;};
 
  int getNmodes(){return perturb_Nmodes;};
  void set_perturbmodes(PosType * ListModes, const int Nmodes);
  void get_perturbmodes(PosType * ListModes, const int Nmodes);
  void get_perturbmodes(std::vector<PosType> & ListModes);
  double * getq() { return qpriv; };
 
  void set_sigma(PosType my_sigma){sigma = my_sigma; };
  
private:
 
   // Things added to manipulate and fit lenses.
   int check_model(int Nimages,int Nsources,int Nlenses,int *pairing,double **xob,double *x_center
     ,int Nmod,double *mod,double **xg,double Re2,double **dx_sub,double **Amag,double ytol);
   double find_axis(double *mod,int Nmod);
   double deflect_translated(double beta,double *mod,double *x,double *y,double *mag,int N
      ,int Nlenses,double Re2,double *x2);
   double ElliptisizeLens(int Nimages,int Nsources,int Nlenses,int *pairing,double **xob
            ,double *xc,double **xg,double sigG,double beta,int Nmod
            ,double *mod,double **dx,double *re2,double *q);
  void setCosmology(const COSMOLOGY& cosmo);
  // output of ElliptisizeLens
  double qpriv[7];
 
  //void find_lens(int Nimages,int Nsources,int *pairing,double **xob,double *xg,double beta
  //               ,int N,int *degen,double *mod,double **v,double **dx);
 
  //void assignParams(InputParams& params);
  
  // Copied from lens_expand.c
  double lens_expand(double beta,double *mod,int Nmodes,PosType const *x,PosType *alpha,KappaType *gamma,KappaType *phi);
};
 
// TODO: BEN finish this documentation for perturbation parameters.
class LensHaloAnaNSIE : public LensHaloBaseNSIE{
  
  void abstractfunction(){}; // pure virtual
 
public:
  //LensHaloAnaNSIE(InputParams& params,const COSMOLOGY &cosmo,bool verbose = false);
  //LensHaloAnaNSIE(const COSMOLOGY& cosmo);
  ~LensHaloAnaNSIE();
  
  void assignParams(InputParams& params);
  //double FractionWithinRe(double rangeInRei);
  void PrintLens(bool show_substruct,bool show_stars);
  
  void setCosmology(const COSMOLOGY& cosmo);
  
  // in randoimize_lens.c
  void RandomizeHost(long *seed,bool tables);
  void RandomizeSigma(long *seed,bool tables);
  void RandomlyDistortLens(long *seed,int Nmodes);
  void AlignedRandomlyDistortLens(long *seed,double theta,int n);
  //void RandomizeSubstructure(double rangeInRei,long *seed);
  //void RandomizeSubstructure2(double rangeInRei,long *seed);
  //void RandomizeSubstructure3(double rangeInRei,long *seed);
  
  virtual PosType get_sigma(){return sigma;};
  //PosType getRsize(){return Rsize;};
  virtual PosType get_fratio(){return fratio;};
  virtual PosType get_pa(){return pa;};
  virtual PosType get_rcore(){return rcore;};
  
 
private:
  
  // Things added to manipulate and fit lenses.
  int check_model(int Nimages,int Nsources,int Nlenses,int *pairing,double **xob,double *x_center
                  ,int Nmod,double *mod,double **xg,double Re2,double **dx_sub,double **Amag,double ytol);
  double find_axis(double *mod,int Nmod);
  double deflect_translated(double beta,double *mod,double *x,double *y,double *mag,int N
                            ,int Nlenses,double Re2,double *x2);
  
};
 
 
// in mark_points.c
void MarkPoints(TreeHndl s_tree,LensHaloAnaNSIE *lens,bool sb_cut,short invert);
void _MarkPoints(TreeHndl s_tree,LensHaloAnaNSIE *lens,bool *sbcut);
bool InSource(double *ray,LensHaloAnaNSIE *lens,bool surfacebright);
 
 
#endif