From 4ea7a05ea9a7352602f1f48a860fd81c36e8bc04 Mon Sep 17 00:00:00 2001 From: Aaron LI Date: Mon, 20 Feb 2017 12:26:17 +0800 Subject: Rename mass_profile to src; Add install & uninstall to Makefile --- src/fit_beta_sbp.cpp | 534 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 534 insertions(+) create mode 100644 src/fit_beta_sbp.cpp (limited to 'src/fit_beta_sbp.cpp') diff --git a/src/fit_beta_sbp.cpp b/src/fit_beta_sbp.cpp new file mode 100644 index 0000000..295fa1e --- /dev/null +++ b/src/fit_beta_sbp.cpp @@ -0,0 +1,534 @@ +/* + Perform a double-beta density model fitting to the surface brightness data + Author: Junhua Gu + Last modified: 2016.06.07 + This code is distributed with no warrant +*/ + +#include +#include +#include +#include +#include +#include "beta_cfg.hpp" +#include "dump_fit_qdp.hpp" +#include "report_error.hpp" +#include "vchisq.hpp" +#include "chisq.hpp" +#include "beta.hpp" +#include "models/beta1d.hpp" +#include +#include +#include +#include +#include "spline.hpp" + +using namespace std; +using namespace opt_utilities; +//double s=5.63136645E20; +const double kpc=3.086E21;//kpc in cm +const double Mpc=kpc*1000; + +double beta_func(double r, double n0, double rc, double beta) +{ + return abs(n0) * pow(1+r*r/rc/rc, -3./2.*abs(beta)); +} + +//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 +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; + //prepend the zero point to radius list + radii.push_back(0.0); + radii_all.push_back(0.0); + //read sbp and sbp error data + cerr << "Read surface brightness profile data ..." << endl; + ifstream ifs(cfg.sbp_data.c_str()); + std::string line; + if (ifs.is_open()) + { + while(std::getline(ifs, line)) + { + if (line.empty()) + continue; + + std::istringstream iss(line); + double x, xe, y, ye; + if ((iss >> x >> xe >> y >> ye)) + { + std::cerr << "sbprofile data: " + << x << ", " << xe << ", " << y << ", " << ye + << std::endl; + radii.push_back(x+xe); /* NOTE: use outer radii of regions */ + radii_all.push_back(x+xe); + sbps.push_back(y); + sbps_all.push_back(y); + sbpe.push_back(ye); + sbpe_all.push_back(ye); + } + else + { + std::cerr << "skipped line: " << line << std::endl; + } + } + } + else + { + std::cerr << "ERROR: cannot open file: " << cfg.sbp_data.c_str() + << std::endl; + return 1; + } + + //initialize the cm/pixel value + double cm_per_pixel=cfg.cm_per_pixel; + double rmin; + if(cfg.rmin_pixel>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()) + { + cerr << "radius_max: " << radii.back() << endl; + break; + } + cf.add_point(x,y); + } + cf.gen_spline(); + + //read temperature profile data + cerr << "Read temperature profile data ..." << endl; + spline_func_obj Tprof; + int tcnt=0; + for(ifstream ifs1(cfg.tprofile.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; + beta betao; + //attach the cooling function + a.attach_cfunc(cf); + a.set_cm_per_pixel(cm_per_pixel); + a.attach_model(betao); + f.set_model(a); + //chi^2 statistic + vchisq c; + c.verbose(true); + c.set_limit(); + f.set_statistic(c); + //optimization method + f.set_opt_method(powell_method >()); + //initialize the initial values + double n0=0; + //double beta=atof(arg_map["beta"].c_str()); + double beta=0; + double rc=0; + double bkg_level=0; + + 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") + { + f.set_param_lower_limit(pname,.3); + f.set_param_upper_limit(pname,1.4); + } + } + } + + f.fit(); + f.fit(); + std::vector p=f.get_all_params(); + n0=f.get_param_value("n0"); + rc=f.get_param_value("rc"); + beta=f.get_param_value("beta"); + //output the datasets and fitting results + ofstream param_output("beta_param.txt"); + for(size_t i=0;i=1) + { + ofs_sbp<<"line on 2"<=1) + { + ofs_sbp<<"no no no"<