overzier_source.h

/*
 * overzier_source.h
 *
 *  Created on: Mar 6, 2010
 *      Author: R.B. Metcalf
 */
#ifndef OVERZIER_SOURCE_H_
#define OVERZIER_SOURCE_H_
 
#include "source.h"
#include "sersic_source.h"
 
// define PI here if not done via include
#ifndef PI
#define PI 3.141592653589793238462643383279502884
#endif
 
class SourceOverzier : public Source
{
public:
 //SourceOverzier();
 SourceOverzier(PosType mag,PosType mag_bulge,Band band,double zero_point,PosType Reff,PosType Rdisk,PosType PA,PosType inclination,unsigned long my_id,PosType my_z,const PosType *theta);
  
  SourceOverzier(const SourceOverzier &s);
  SourceOverzier& operator=(const SourceOverzier &s);
 virtual ~SourceOverzier();
 
 void setInternals(PosType mag,PosType mag_bulge,Band band,PosType Reff,PosType Rdisk,PosType PA,PosType inclination,unsigned long my_id,PosType my_z,const PosType *my_theta);
  virtual PosType SurfaceBrightness(const PosType *x) const;
 PosType getTotalFlux() const;
 void printSource();
  
  PosType getMag() const { return current.mag; }
  PosType getMag(Band band) const ;
  PosType getMagBulge() const { return current.mag_bulge; }
  PosType getMagBulge(Band band) const;
  
  /*
 void setUMag(PosType m) { mag_u = m; }
 void setGMag(PosType m) { mag_g = m; }
 void setRMag(PosType m) { mag_r = m; }
 void setIMag(PosType m) { mag_i = m; }
 void setZMag(PosType m) { mag_z = m; }
 void setJMag(PosType m) { mag_J = m; }
 void setHMag(PosType m) { mag_H = m; }
 void setKMag(PosType m) { mag_Ks = m; }
  */
  
  virtual void setMag(PosType my_mag,Band band,PosType zero_point){
    current.mag_map[band] = my_mag;
    zeropoints[band] = zero_point;
  }
  virtual void setMagBulge(PosType my_mag,Band band,PosType zero_point){
    current.bulge_mag_map[band] = my_mag;
    zeropoints[band] = zero_point;
  }
  
 PosType getReff() const { return current.Reff/arcsecTOradians; }
 PosType getRdisk() const { return current.Rdisk/arcsecTOradians; }
 
 PosType getBtoT() const { return pow(10,(-current.mag_bulge + current.mag)/2.5); }
 PosType getPA() const { return current.PA; }
 PosType getInclination() const { return current.inclination;}
  float getSEDtype() const {return sedtype;}
  void setSEDtype(float s){ sedtype = s;}
 
  virtual void changeBand(Band band);
 
 inline PosType getMinSize(PosType f) {return std::min(1.678*current.Reff*fabs(cos(current.inclination))*pow(-log (f)/7.67,4),current.Rdisk*(-log (f)/1.67));}
 
  static PosType *getx(SourceOverzier &sourceo){return sourceo.source_x.x;}
 
protected:
  
  SourceSersic spheroid;
  static std::map<Band,PosType> zeropoints;
  
  float sedtype = -1;
  // renormalize the disk and bulge to agree with current mag and mag_bulge
  void renormalize_current();
 void assignParams(InputParams& params);
 
  struct Params{
 
    PosType Reff=0;
    PosType Rdisk=0;
 
    //PosType BtoT;
    PosType PA=0;
    PosType inclination=0;
    
    double cosPA;
    double sinPA;
 
    PosType cxx=0,cyy=0;
    PosType sbDo=0;
    PosType mag=0;
    PosType mag_bulge=0;
    
    // colors
    std::map<Band,double> mag_map;
    std::map<Band,double> bulge_mag_map;
    Band band;
    
    void print(){
      std::cout << "Reff :" << Reff/arcsecTOradians << " arcsec ";
      std::cout << "Rdisk :" << Rdisk/arcsecTOradians << " arcsec ";
      std::cout << "PA :" << PA << " ";
      std::cout << "inclination :" << inclination << " radians";
      std::cout << "sbDo :" << sbDo << " ";
      std::cout << "mag :" << mag << " ";
      std::cout << "mag_bulge :" << mag_bulge << " ";
      
      std::cout << "BtoT :" << pow(10,(-mag_bulge + mag)/2.5) << std::endl;
    }
  };
 
  Params current;
 // optional position variables
};
 
class SourceOverzierPlus : public SourceOverzier
{
public:
 
  SourceOverzierPlus(
                     PosType my_mag          
                     ,PosType my_mag_bulge   
                     ,Band band              
                     ,PosType mag_zero_point       
                     ,PosType my_Reff        
                     ,PosType my_bulge_q     
                     ,PosType my_bulge_index 
                     ,PosType my_Rdisk       
                     ,PosType my_PA          
                     ,PosType inclination    
                     ,unsigned long my_id
                     ,PosType my_z           
                     ,const PosType *theta   
                     ,Utilities::RandomNumbers_NR &ran
                     );
  
  SourceOverzierPlus(
                     PosType my_mag         
                     ,PosType my_mag_bulge  
                     ,Band band              
                     ,PosType mag_zero_point       
                     ,PosType my_Reff       
                     ,PosType my_Rdisk      
                     ,PosType my_PA         
                     ,PosType inclination   
                     ,unsigned long my_id
                     ,PosType my_z           
                     ,const PosType *theta   
                     ,PosType my_sersic_index_max
                     ,PosType my_sersic_index_min
                     ,Utilities::RandomNumbers_NR &ran
);
 
 
  ~SourceOverzierPlus();
  
  SourceOverzierPlus(const SourceOverzierPlus &p);
  SourceOverzierPlus & operator=(const SourceOverzierPlus &p);
  
  //*** meed to be able to change band
  //*** put modes and phases into surface brightness
  //*** possible put gaussian texture on disk
 
  PosType SurfaceBrightness(const PosType *y) const;
 
  void setMag(PosType my_mag,Band band,double zeropoint){
    current.mag_map[band] = my_mag;
    original.mag_map[band] = my_mag;
    zeropoints[band] = zeropoint;
  }
  void setMagBulge(PosType my_mag,Band band,double zeropoint){
    current.bulge_mag_map[band] = my_mag;
    original.bulge_mag_map[band] = my_mag;
    zeropoints[band] = zeropoint;
    spheroid.setMag(my_mag,band,zeropoint);
  }
  
  PosType getPA() const { return PA; }
  void setPA(double pa); // in radians
 
  int getNarms() const {return Narms;}
  PosType getArmAmplitude() const {return Ad;}
  PosType getArmAlpha() const {return arm_alpha;}
  PosType getSphIndex() const {return spheroid.getSersicIndex();}
  PosType getSphAxisRatio() const {return spheroid.getAxesRatio();}
  //PosType getSphPA() const {return spheroid.getPA();}
  
  void changeBand(Band band);
  static PosType* getx(SourceOverzierPlus &sourceo){return sourceo.source_x.x;}
 
  virtual inline void setTheta(PosType *xx){
    source_x[0] = xx[0];
    source_x[1] = xx[1];
  }
  virtual void setTheta(PosType my_x,PosType my_y){
    source_x[0] = my_x;
    source_x[1] = my_y;
  }
  virtual void setTheta(const Point_2d &p){
    source_x = p;
  }
  void setBulgeAxisRatio(PosType q){
    spheroid.setAxesRatio(q);
  }
  void setBulgeSersicIndex(double index){
    spheroid.setSersicIndex(index);
  }
  void randomize(Utilities::RandomNumbers_NR &ran);
private:
  int Narms;
  PosType Ad,mctalpha,arm_alpha;
  SourceSersic spheroid;
  std::vector<PosType> modes;
  PosType disk_phase;
  PosType cosPA,sinPA,cosi,PA;
  
  SourceOverzier::Params original;  // original parameters
 
  PosType sersic_index_max;
  PosType sersic_index_min;
};
#endif /* GALAXIES_OVERZIER_H_ */