sourceAnaGalaxy.h

/*
 * souceAnaLens.h
 *
 *  Created on: Aug 13, 2012
 *      Author: bmetcalf
 */
 
#ifndef SOURCE_ANA_H_
#define SOURCE_ANA_H_
 
#include "source.h"
#include "overzier_source.h"
#include "simpleTreeVec.h"
#include "utilities.h"
 
class SourceMultiAnaGalaxy: public Source{
public:
 SourceMultiAnaGalaxy(PosType mag,PosType mag_bulge,Band band,PosType zero_point, PosType Reff, PosType Rdisk, PosType PA, PosType inclination,PosType my_z,PosType *my_theta
                       ,Utilities::RandomNumbers_NR &ran);
 SourceMultiAnaGalaxy(SourceOverzierPlus *my_galaxy);
 //SourceMultiAnaGalaxy(InputParams& params,Utilities::RandomNumbers_NR &ran);
 ~SourceMultiAnaGalaxy();
 
 PosType SurfaceBrightness(const PosType* x) const {
  PosType sb = galaxies[index].SurfaceBrightness(x);
  if (sb < sb_limit) return 0.;
  return sb;
  }
 
 PosType getTotalFlux() const {return mag_to_flux(galaxies[index].getMag(),mag_zero_point);}
 
 void printSource();
 // Add a pre-constructed galaxy to the source collection
 void AddAGalaxy(SourceOverzierPlus *my_galaxy){galaxies.push_back(*my_galaxy);}
 
 SourceOverzierPlus& setIndex (std::size_t i){
  if(i < galaxies.size())
   index = i;
  return galaxies[index];
 }
 SourceOverzierPlus& operator[] (std::size_t i){
  if(i < galaxies.size())
   return galaxies[i];
  return galaxies[index];
 }
 
 const SourceOverzierPlus& operator[] (std::size_t i) const {
  if(i < galaxies.size())
   return galaxies[i];
  return galaxies[index];
 }
 
 SourceOverzierPlus& CurrentGalaxy(){
  return galaxies[index];
 }
 
 //PosType getZ() const {return galaxies[index].getZ();}
   PosType getZ() const {return galaxies[index].getZ();}
  //PosType getRadius() const {return max(galaxies[index]->Reff,galaxies[index]->Rdisk);}
 PosType getRadius() const {return galaxies[index].getRadius();}
 void setZ(PosType my_z){ galaxies[index].setZ(my_z);}
  void resetBand(Band my_band){
    for(size_t i=0;i<galaxies.size();++i) galaxies[i].changeBand(my_band);
    band = my_band;
  }
 
 unsigned long getID(){return galaxies[index].getID();}
 
 
 Point_2d getTheta() const{return galaxies[index].getTheta();}
 void setTheta(PosType my_theta[2]){galaxies[index].setTheta(my_theta);}
  void setTheta(PosType my_x,PosType my_y){galaxies[index].setTheta(my_x, my_y);}
  void setTheta(const Point_2d &p){galaxies[index].setTheta(p);}
 
 std::size_t getNumberOfGalaxies() const {return galaxies.size();}
 
 void multiplier(PosType z,PosType mag_cut,int Multiplicity,Utilities::RandomNumbers_NR &ran);
  void sortInRedshift();
  void sortInMag(Band tmp_band);
  void sortInID();
  PosType getFOV(){return (rangex[1]-rangex[0])*(rangey[1]-rangey[0])*180*180/PI/PI;}
  
  std::size_t findclosestonsky(PosType theta[],PosType *radius){
    size_t index;
    searchtree->NearestNeighbor(theta,*radius,index);
    return index;
  }
 
  void findclosestonsky(PosType theta[],std::vector<PosType> &radius,std::vector<size_t> &indexes){
    indexes.resize(radius.size());
    searchtree->NearestNeighbors(theta,radius.size(),radius,indexes);
    return;
  }
 
  void findonsky(PosType theta[],float radius,std::list<size_t> &indexes){
    indexes.clear();
    searchtree->PointsWithinCircle(theta,radius,indexes);
    return;
  }
  
  void findnear(PosType theta[],float radius,std::list<size_t> &indexes,PosType z_range[]){
    indexes.clear();
    searchtree->PointsWithinCircle(theta,radius,indexes);
    
    PosType z;
    // make redshift cut for lens sources
    for(std::list<size_t>::iterator itt = indexes.begin();
        itt != indexes.end(); ++itt){
      z = galaxies[*itt].getZ();
      if(z < z_range[0] || z > z_range[1]){
        indexes.erase(itt);
        if(itt != indexes.begin()) --itt;
      }
    }
    
    return;
  }
 
private:
  // make it uncopyable
  //SourceMultiAnaGalaxy(SourceMultiAnaGalaxy &s){};
  SourceMultiAnaGalaxy & operator=(SourceMultiAnaGalaxy &s){return s;}
  
 Band band;
 float mag_limit;
 std::size_t index;
 
 std::vector<SourceOverzierPlus> galaxies;
  TreeSimpleVec<SourceOverzierPlus> *searchtree;
 std::string input_gal_file;
 
 void readDataFileMillenn(std::map<Band,double> zeropoints,Utilities::RandomNumbers_NR &ran);
 void assignParams(InputParams& params);
 
  PosType rangex[2],rangey[2];
};
 
bool redshiftcompare(SourceOverzierPlus s1,SourceOverzierPlus s2);
bool magcompare(SourceOverzierPlus s1,SourceOverzierPlus s2);
bool idcompare(SourceOverzierPlus s1,SourceOverzierPlus s2);
 
class SourceMultiShapelets: public Source{
public:
  
  //SourceMultiShapelets(InputParams& params);
  
  SourceMultiShapelets(double mag_zero_point):Source(0,Point_2d(0,0),0,-1,mag_zero_point){};
  
  SourceMultiShapelets(const std::string &my_shapelets_folder  
                       ,Band my_band  
                       ,double my_max_mag_limit  
                       ,double my_min_mag_limit  
                       ,double my_z_max          
                       ,double my_sb_limit       
                       ,double maximum_radius   
                       ,double zero_point        
                       );
 
  void input(const std::string &my_shapelets_folder  
                       ,Band my_band  
                       ,double my_max_mag_limit  
                       ,double my_min_mag_limit  
                       ,double my_z_max          
                       ,double my_sb_limit     
                       ,double maximum_radius   
                       ,double zero_point        
                     );
 
  ~SourceMultiShapelets();
  void sortInRedshift();
  void sortInMag();
  PosType SurfaceBrightness(const PosType* x)const {
  PosType sb = galaxies[index].SurfaceBrightness(x);
  if (sb < sb_limit) return 0.;
  return sb;
  }
    
 void printSource();
  std::size_t getNumberOfGalaxies() const {return galaxies.size();}
  std::size_t size() const {return galaxies.size();}
 
 PosType getTotalFlux() const {return mag_to_flux(galaxies[index].getMag(),mag_zero_point);}
 
 Point_2d getTheta() const {return galaxies[index].getTheta();}
 void setTheta(PosType my_theta[2]){galaxies[index].setTheta(my_theta);}
 void setTheta(PosType my_x,PosType my_y){galaxies[index].setTheta(my_x, my_y);}
  void setTheta(const Point_2d &p){galaxies[index].setTheta(p);}
 
 PosType getZ() const {return galaxies[index].getZ();}
 //PosType getRadius() const {return max(galaxies[index]->Reff,galaxies[index]->Rdisk);}
 PosType getRadius() const {return galaxies[index].getRadius();}
 
 SourceShapelets& setIndex (std::size_t i){
  if(i < galaxies.size())
   index = i;
  return galaxies[index];
 }
  size_t getIndex() const {return index;}
  
  SourceShapelets& operator[] (std::size_t i){
    if(i < galaxies.size())
      return galaxies[i];
    return galaxies[index];
  }
 
  SourceShapelets& back(){
    return galaxies.back();
  }
 
 const SourceShapelets& operator[] (std::size_t i) const {
  if(i < galaxies.size())
   return galaxies[i];
  return galaxies[index];
 }
    
 SourceShapelets& CurrentGalaxy(){
  return galaxies[index];
 }
    
 SourceShapelets& setBand (Band b,double zeropoint){
        band = b;
        for (int i = 0; i < galaxies.size(); i++)
            galaxies[i].setActiveBand(band);
      return galaxies[index];
  }
  
  int locateZ(PosType z) const {
      return Utilities::locate<SourceShapelets,PosType>(galaxies, z, [](PosType z,const SourceShapelets &s){return (z < s.getZ());});
  }
  
  int getCurrentID(){
    return galaxies[index].id;
  }
 
  std::vector<SourceShapelets> galaxies;
 
private:
 void assignParams(InputParams& params);
  std::size_t index;
  float max_mag_limit;
  float min_mag_limit;
  float z_max;  // maximum redshift allowed
  Band band;
  double radius_max;
 
  void readCatalog();
 std::string shapelets_folder;
};
 
bool redshiftcompare_shap(SourceShapelets s1,SourceShapelets s2);
bool magcompare_shap(SourceShapelets s1,SourceShapelets s2);
 
#endif /* SOURCE_H_ */