From ffd178e0bd72562a3c2cff9747b6e656edc881dc Mon Sep 17 00:00:00 2001 From: Aaron LI Date: Fri, 27 May 2016 22:47:24 +0800 Subject: Add mass_profile tools * These tools are mainly use to calculate the total gravitational mass profile, as well as the intermediate products (e.g., surface brightness profile fitting, gas density profile, NFW fitting, etc.) * There are additional tools for calculating the luminosity and flux. * These tools mainly developed by Junhua GU, and contributed by Weitian (Aaron) LI, and Zhenghao ZHU. --- mass_profile/fit_dbeta_sbp.cpp | 585 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 585 insertions(+) create mode 100644 mass_profile/fit_dbeta_sbp.cpp (limited to 'mass_profile/fit_dbeta_sbp.cpp') diff --git a/mass_profile/fit_dbeta_sbp.cpp b/mass_profile/fit_dbeta_sbp.cpp new file mode 100644 index 0000000..b74dc0c --- /dev/null +++ b/mass_profile/fit_dbeta_sbp.cpp @@ -0,0 +1,585 @@ +/* + Perform a double-beta density model fitting to the surface brightness data + Author: Junhua Gu + Last modified: 2011.01.01 + This code is distributed with no warrant +*/ + + +#include +#include +#include +using namespace std; +#include "vchisq.hpp" +#include "dbeta.hpp" +#include "beta_cfg.hpp" +#include +#include +#include +#include +#include "spline.h" +using namespace opt_utilities; +//double s=5.63136645E20; +const double kpc=3.086E21;//kpc in cm +const double Mpc=kpc*1000; +double dbeta_func(double r, + double n01,double rc1,double beta1, + double n02,double rc2,double beta2) +{ + + return abs(n01)*pow(1+r*r/rc1/rc1,-3./2.*abs(beta1))+abs(n02)*pow(1+r*r/rc2/rc2,-3./2.*abs(beta2)); +} + + + //calculate critical density from z, under following cosmological constants +static double calc_critical_density(double z, + const double H0=2.3E-18, + const double Omega_m=.27) +{ + const double G=6.673E-8;//cm^3 g^-1 s^2 + const double E=std::sqrt(Omega_m*(1+z)*(1+z)*(1+z)+1-Omega_m); + const double H=H0*E; + return 3*H*H/8/pi/G; +} + + +//A class enclosing the spline interpolation method of cooling function +//check spline.h for more detailed information +//this class is a thin wrapper for the spline class defined in spline.h +class spline_func_obj + :public func_obj +{ + //has an spline object + spline spl; +public: + //This function is used to calculate the intepolated value + double do_eval(const double& x) + { + return spl.get_value(x); + } + + //we need this function, when this object is performing a clone of itself + spline_func_obj* do_clone()const + { + return new spline_func_obj(*this); + } + +public: + //add points to the spline object, after which the spline will be initialized + void add_point(double x,double y) + { + spl.push_point(x,y); + } + + //before getting the intepolated value, the spline should be initialzied by calling this function + void gen_spline() + { + spl.gen_spline(0,0); + } +}; + +int main(int argc,char* argv[]) +{ + if(argc!=2) + { + cerr<"< radii; + std::vector sbps; + std::vector sbpe; + std::vector radii_all; + std::vector sbps_all; + std::vector sbpe_all; + //read sbp and sbp error data + for(ifstream ifs(cfg.sbp_file.c_str());;) + { + assert(ifs.is_open()); + double x,xe; + ifs>>x>>xe; + if(!ifs.good()) + { + break; + } + if(x/xe<2) + { + break; + } + cerr<>x; + if(!ifs.good()) + { + break; + } + cerr<0) + { + rmin=cfg.rmin_pixel; + } + else + { + rmin=cfg.rmin_kpc*kpc/cm_per_pixel; + } + + cerr<<"rmin="< radii_tmp,sbps_tmp,sbpe_tmp; + radii_tmp.resize(radii.size()); + sbps_tmp.resize(sbps.size()); + sbpe_tmp.resize(sbpe.size()); + copy(radii.begin(),radii.end(),radii_tmp.begin()); + copy(sbps.begin(),sbps.end(),sbps_tmp.begin()); + copy(sbpe.begin(),sbpe.end(),sbpe_tmp.begin()); + for(list::iterator i=radii_tmp.begin();i!=radii_tmp.end();) + { + if(*i>x>>y; + if(!ifs.good()) + { + break; + } + cerr<radii.back()) + { + break; + } + //cf.add_point(x,y*2.1249719395939022e-68);//change with source + cf.add_point(x,y);//change with source + } + cf.gen_spline(); + + //read temperature profile data + spline_func_obj Tprof; + int tcnt=0; + for(ifstream ifs1(cfg.T_file.c_str());;++tcnt) + { + assert(ifs1.is_open()); + double x,y; + ifs1>>x>>y; + if(!ifs1.good()) + { + break; + } + cerr<,std::vector > ds; + ds.add_data(data,std::vector >(radii,sbps,sbpe,sbpe,radii,radii)); + + //initial fitter + fitter,vector,vector,double> f; + f.load_data(ds); + //initial the object, which is used to calculate projection effect + projector a; + bool tie_beta=false; + if(cfg.param_map.find("beta")!=cfg.param_map.end() + &&cfg.param_map.find("beta1")==cfg.param_map.end() + &&cfg.param_map.find("beta2")==cfg.param_map.end()) + { + dbeta2 dbetao; + a.attach_model(dbetao); + tie_beta=true; + } + else if((cfg.param_map.find("beta1")!=cfg.param_map.end() + ||cfg.param_map.find("beta2")!=cfg.param_map.end()) + &&cfg.param_map.find("beta")==cfg.param_map.end()) + { + dbeta dbetao; + a.attach_model(dbetao); + tie_beta=false; + } + else + { + cerr<<"Error, cannot decide whether to tie beta together or let them vary freely!"< c; + c.verbose(true); + c.set_limit(); + f.set_statistic(c); + //optimization method + f.set_opt_method(powell_method >()); + //initialize the initial values + double n01=0; + double rc1=0; + double n02=0; + double rc2=0; + double beta=0; + double bkg=0; + if(tie_beta) + { + f.set_param_value("beta",.7); + f.set_param_lower_limit("beta",.3); + f.set_param_upper_limit("beta",1.4); + } + else + { + f.set_param_value("beta1",.7); + f.set_param_lower_limit("beta1",.3); + f.set_param_upper_limit("beta1",1.4); + f.set_param_value("beta2",.7); + f.set_param_lower_limit("beta2",.3); + f.set_param_upper_limit("beta2",1.4); + } + for(std::map >::iterator i=cfg.param_map.begin(); + i!=cfg.param_map.end();++i) + { + std::string pname=i->first; + f.set_param_value(pname,i->second.at(0)); + if(i->second.size()==3) + { + double a1=i->second[1]; + double a2=i->second[2]; + double u=std::max(a1,a2); + double l=std::min(a1,a2); + f.set_param_upper_limit(pname,u); + f.set_param_lower_limit(pname,l); + } + else + { + if(pname=="beta"||pname=="beta1"||pname=="beta2") + { + f.set_param_lower_limit(pname,.3); + f.set_param_upper_limit(pname,1.4); + } + } + } + + + + //perform the fitting, first freeze beta1, beta2, rc1, and rc2 + if(tie_beta) + { + f.set_param_modifier(freeze_param,std::vector,std::vector,std::string>("beta")+ + freeze_param,std::vector,std::vector,std::string>("rc1")+ + freeze_param,std::vector,std::vector,std::string>("rc2") + ); + } + else + { + f.set_param_modifier(freeze_param,std::vector,std::vector,std::string>("beta1")+ + freeze_param,std::vector,std::vector,std::string>("beta2")+ + freeze_param,std::vector,std::vector,std::string>("rc1")+ + freeze_param,std::vector,std::vector,std::string>("rc2") + ); + } + + f.fit(); + + f.clear_param_modifier(); + + //then perform the fitting, freeze beta1 and beta2 + //f.set_param_modifier(freeze_param,std::vector,std::vector,std::string>("beta")); + //f.set_param_modifier(freeze_param,std::vector,std::vector,std::string>("bkg")); + f.fit(); + //f.clear_param_modifier(); + + //finally thaw all parameters + f.fit(); + double beta1=0; + double beta2=0; + + n01=f.get_param_value("n01"); + rc1=f.get_param_value("rc1"); + n02=f.get_param_value("n02"); + rc2=f.get_param_value("rc2"); + if(tie_beta) + { + beta=f.get_param_value("beta"); + beta1=beta; + beta2=beta; + } + else + { + beta1=f.get_param_value("beta1"); + beta2=f.get_param_value("beta2"); + } + //output the params + ofstream param_output("dbeta_param.txt"); + //output the datasets and fitting results + for(int i=0;i