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 --- mass_profile/Makefile | 86 ------ mass_profile/README.md | 37 --- mass_profile/beta.hpp | 45 --- mass_profile/beta_cfg.cpp | 85 ------ mass_profile/beta_cfg.hpp | 23 -- mass_profile/calc_lx_beta.cpp | 497 ------------------------------ mass_profile/calc_lx_dbeta.cpp | 548 --------------------------------- mass_profile/chisq.hpp | 206 ------------- mass_profile/dbeta.hpp | 108 ------- mass_profile/dump_fit_qdp.cpp | 55 ---- mass_profile/dump_fit_qdp.hpp | 16 - mass_profile/fit_beta_sbp.cpp | 534 -------------------------------- mass_profile/fit_dbeta_sbp.cpp | 586 ------------------------------------ mass_profile/fit_nfw_mass.cpp | 159 ---------- mass_profile/fit_wang2012_model.cpp | 195 ------------ mass_profile/nfw.hpp | 56 ---- mass_profile/projector.hpp | 214 ------------- mass_profile/report_error.cpp | 39 --- mass_profile/report_error.hpp | 6 - mass_profile/spline.hpp | 110 ------- mass_profile/vchisq.hpp | 137 --------- mass_profile/wang2012_model.hpp | 67 ----- src/Makefile | 98 ++++++ src/README.md | 37 +++ src/beta.hpp | 45 +++ src/beta_cfg.cpp | 85 ++++++ src/beta_cfg.hpp | 23 ++ src/calc_lx_beta.cpp | 497 ++++++++++++++++++++++++++++++ src/calc_lx_dbeta.cpp | 548 +++++++++++++++++++++++++++++++++ src/chisq.hpp | 206 +++++++++++++ src/dbeta.hpp | 108 +++++++ src/dump_fit_qdp.cpp | 55 ++++ src/dump_fit_qdp.hpp | 16 + src/fit_beta_sbp.cpp | 534 ++++++++++++++++++++++++++++++++ src/fit_dbeta_sbp.cpp | 586 ++++++++++++++++++++++++++++++++++++ src/fit_nfw_mass.cpp | 159 ++++++++++ src/fit_wang2012_model.cpp | 195 ++++++++++++ src/nfw.hpp | 56 ++++ src/projector.hpp | 214 +++++++++++++ src/report_error.cpp | 39 +++ src/report_error.hpp | 6 + src/spline.hpp | 110 +++++++ src/vchisq.hpp | 137 +++++++++ src/wang2012_model.hpp | 67 +++++ 44 files changed, 3821 insertions(+), 3809 deletions(-) delete mode 100644 mass_profile/Makefile delete mode 100644 mass_profile/README.md delete mode 100644 mass_profile/beta.hpp delete mode 100644 mass_profile/beta_cfg.cpp delete mode 100644 mass_profile/beta_cfg.hpp delete mode 100644 mass_profile/calc_lx_beta.cpp delete mode 100644 mass_profile/calc_lx_dbeta.cpp delete mode 100644 mass_profile/chisq.hpp delete mode 100644 mass_profile/dbeta.hpp delete mode 100644 mass_profile/dump_fit_qdp.cpp delete mode 100644 mass_profile/dump_fit_qdp.hpp delete mode 100644 mass_profile/fit_beta_sbp.cpp delete mode 100644 mass_profile/fit_dbeta_sbp.cpp delete mode 100644 mass_profile/fit_nfw_mass.cpp delete mode 100644 mass_profile/fit_wang2012_model.cpp delete mode 100644 mass_profile/nfw.hpp delete mode 100644 mass_profile/projector.hpp delete mode 100644 mass_profile/report_error.cpp delete mode 100644 mass_profile/report_error.hpp delete mode 100644 mass_profile/spline.hpp delete mode 100644 mass_profile/vchisq.hpp delete mode 100644 mass_profile/wang2012_model.hpp create mode 100644 src/Makefile create mode 100644 src/README.md create mode 100644 src/beta.hpp create mode 100644 src/beta_cfg.cpp create mode 100644 src/beta_cfg.hpp create mode 100644 src/calc_lx_beta.cpp create mode 100644 src/calc_lx_dbeta.cpp create mode 100644 src/chisq.hpp create mode 100644 src/dbeta.hpp create mode 100644 src/dump_fit_qdp.cpp create mode 100644 src/dump_fit_qdp.hpp create mode 100644 src/fit_beta_sbp.cpp create mode 100644 src/fit_dbeta_sbp.cpp create mode 100644 src/fit_nfw_mass.cpp create mode 100644 src/fit_wang2012_model.cpp create mode 100644 src/nfw.hpp create mode 100644 src/projector.hpp create mode 100644 src/report_error.cpp create mode 100644 src/report_error.hpp create mode 100644 src/spline.hpp create mode 100644 src/vchisq.hpp create mode 100644 src/wang2012_model.hpp diff --git a/mass_profile/Makefile b/mass_profile/Makefile deleted file mode 100644 index 17b7a0e..0000000 --- a/mass_profile/Makefile +++ /dev/null @@ -1,86 +0,0 @@ -# -# Makefile for `chandra-acis-analysis/mass_profile` tools -# -# Junhua GU -# Weitian LI -# 2016-06-07 -# - - -CXX ?= g++ -CXXFLAGS += -Wall -CXXFLAGS += -Werror - -ifdef OPENMP - CXXFLAGS += -fopenmp -endif - -ifdef DEBUG - CXXFLAGS += -g -else - CXXFLAGS += -O2 -endif - -OPT_UTIL_INC ?= -I../opt_utilities - -TARGETS= fit_dbeta_sbp fit_beta_sbp fit_wang2012_model \ - fit_nfw_mass calc_lx_dbeta calc_lx_beta -HEADERS= projector.hpp spline.hpp vchisq.hpp - -all: $(TARGETS) - -# NOTE: -# Object/source files should placed *before* libraries (order matters) - -fit_dbeta_sbp: fit_dbeta_sbp.o beta_cfg.o report_error.o - $(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC) - -fit_beta_sbp: fit_beta_sbp.o beta_cfg.o dump_fit_qdp.o - $(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC) - -fit_wang2012_model: fit_wang2012_model.o - $(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC) - -fit_nfw_mass: fit_nfw_mass.o - $(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC) - -calc_lx_dbeta: calc_lx_dbeta.o beta_cfg.o report_error.o - $(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC) - -calc_lx_beta: calc_lx_beta.o beta_cfg.o dump_fit_qdp.o - $(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC) - - -fit_dbeta_sbp.o: fit_dbeta_sbp.cpp $(HEADERS) - $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC) - -fit_beta_sbp.o: fit_beta_sbp.cpp beta.hpp $(HEADERS) - $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC) - -fit_wang2012_model.o: fit_wang2012_model.cpp wang2012_model.hpp chisq.hpp - $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC) - -fit_nfw_mass.o: fit_nfw_mass.cpp nfw.hpp chisq.hpp - $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC) - -calc_lx_dbeta.o: calc_lx_dbeta.cpp $(HEADERS) - $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC) - -calc_lx_beta.o: calc_lx_beta.cpp beta.hpp $(HEADERS) - $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC) - -beta_cfg.o: beta_cfg.cpp beta_cfg.hpp - $(CXX) $(CXXFLAGS) -c $< - -report_error.o: report_error.cpp report_error.hpp - $(CXX) $(CXXFLAGS) -c $< - -dump_fit_qdp.o: dump_fit_qdp.cpp dump_fit_qdp.hpp - $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC) - -%.o: %.cpp - $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC) - - -clean: - rm -f *.o $(TARGETS) diff --git a/mass_profile/README.md b/mass_profile/README.md deleted file mode 100644 index 0c61821..0000000 --- a/mass_profile/README.md +++ /dev/null @@ -1,37 +0,0 @@ -Mass/Luminosity/Flux Calculation Tools -============================== - -Junhua GU, Weitian LI, Zhenghao ZHU - - -Introduction ------------- -This directory contains the tools/scripts used to fit the temperature -profile and surface brightness file, to calculate the gas density profile -and gravitational mass profile, to calculate luminosity and flux and -other related quantities. - - -NOTE ----- -* Mass calculation references: - + Walker et al. 2012, MNRAS, 422, 3503-3515, - [ADS:2012MNRAS.422.3503W](http://adsabs.harvard.edu/abs/2012MNRAS.422.3503W) - + Ettori et al. 2013, Space Science Reviews, 177, 119-154, - [ADS:2013SSRv..177..119E](http://adsabs.harvard.edu/abs/2013SSRv..177..119E) -* Errors are calculated at 68% confident level. -* NFW profile is employed to extrapolate the mass profile. -* We use a self-proposed model (see ``wang2012_model.hpp``) to describe/fit - the gas temperature profile. -* The single-beta and double-beta models (with a constant background) are used - to describe the gas density profile. - - -TODO ----- -* Merge the scripts/tools for single-beta and double-beta SBP models - (to reduce code duplications) -* Uncertainties/errors calculation **maybe** inappropriate/problematic - (e.g., ``analyze_mass_profile.py``); - why not just use the quantile-based (e.g., Q84.15-Q15.85 ~= 68.3%) - uncertainties or standard deviation??? diff --git a/mass_profile/beta.hpp b/mass_profile/beta.hpp deleted file mode 100644 index 6ae70ac..0000000 --- a/mass_profile/beta.hpp +++ /dev/null @@ -1,45 +0,0 @@ -#ifndef BETA -#define BETA -#include "projector.hpp" - -namespace opt_utilities -{ - template - class beta - :public model,std::vector,std::vector > - { - public: - beta() - { - this->push_param_info(param_info,std::string>("n0",1,0,1E99)); - this->push_param_info(param_info,std::string>("beta",.66,0,1E99)); - this->push_param_info(param_info,std::string>("rc",100,0,1E99)); - } - - public: - beta* do_clone()const - { - return new beta(*this); - } - - std::vector do_eval(const std::vector & x, - const std::vector& p) - { - T n0=std::abs(p[0]); - T beta=p[1]; - T rc=p[2]; - - std::vector result(x.size()-1); - for(size_t i=1;i -#include -using namespace std; - -cfg_map parse_cfg_file(std::istream& is) -{ - cfg_map result; - result.rmin_pixel=-1; - result.rmin_kpc=-1; - for(;;) - { - std::string line; - getline(is,line); - line+="\n"; - if(!is.good()) - { - break; - } - string key; - istringstream iss(line); - iss>>key; - if(key=="sbp_data") - { - string value; - iss>>value; - result.sbp_data=value; - } - else if(key=="cfunc_profile") - { - string value; - iss>>value; - result.cfunc_profile=value; - } - else if(key=="tprofile") - { - string value; - iss>>value; - result.tprofile=value; - } - else if(key=="z") - { - double z; - iss>>z; - result.z=z; - } - else if(key=="cm_per_pixel") - { - double cm_per_pixel; - iss>>cm_per_pixel; - result.cm_per_pixel=cm_per_pixel; - } - else if(key=="rmin_pixel") - { - double v; - iss>>v; - result.rmin_pixel=v; - } - else if(key=="rmin_kpc") - { - double v; - iss>>v; - result.rmin_kpc=v; - } - else - { - std::vector value; - for(;;) - { - double v; - iss>>v; - if(!iss.good()) - { - break; - } - value.push_back(v); - } - if(!value.empty()) - { - result.param_map[key]=value; - } - } - } - return result; -} diff --git a/mass_profile/beta_cfg.hpp b/mass_profile/beta_cfg.hpp deleted file mode 100644 index 27148df..0000000 --- a/mass_profile/beta_cfg.hpp +++ /dev/null @@ -1,23 +0,0 @@ -#ifndef BETA_CFG -#define BETA_CFG - -#include -#include -#include -#include - -struct cfg_map -{ - std::string sbp_data; - std::string cfunc_profile; - std::string tprofile; - double z; - double cm_per_pixel; - double rmin_kpc; - double rmin_pixel; - std::map > param_map; -}; - -cfg_map parse_cfg_file(std::istream& is); - -#endif diff --git a/mass_profile/calc_lx_beta.cpp b/mass_profile/calc_lx_beta.cpp deleted file mode 100644 index 11d7b04..0000000 --- a/mass_profile/calc_lx_beta.cpp +++ /dev/null @@ -1,497 +0,0 @@ -/** - * Calculate the total luminosity and flux within the specified radius. - * - * Base on 'fit_beta_sbp.cpp' and supersede 'calc_lx.cpp' - * - * Author: Junhua Gu - */ - -#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)); -} - -//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<4) - { - cerr< [cfunc2_erg ...]"< 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 - 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()) - { - break; - } - cf.add_point(x,y); - } - cf.gen_spline(); - - //read temperature profile data - 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=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("lx_beta_param.txt"); - for(size_t i=0;i=1) - { - ofs_sbp<<"line on 2"<=1) - { - ofs_sbp<<"no no no"<& pj=dynamic_cast&>(f.get_model()); - pj.attach_cfunc(cf_erg); - - - - mv=f.eval_model_raw(radii,p); - double flux_erg=0; - for(size_t i=0;i -#include -#include -#include -#include "vchisq.hpp" -#include "dbeta.hpp" -#include "beta_cfg.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 dbeta_func(double r, double n01, double rc1, double beta1, - double n02, double rc2, double beta2) -{ - double v1 = abs(n01) * pow(1+r*r/rc1/rc1, -3./2.*abs(beta1)); - double v2 = abs(n02) * pow(1+r*r/rc2/rc2, -3./2.*abs(beta2)); - return v1 + v2; -} - - -//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<4) - { - cerr< [cfunc2_erg ...]"< 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 - 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()) - { - break; - } - cf.add_point(x,y); - } - cf.gen_spline(); - - //read temperature profile data - 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; - 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 =0; - double rc1=0; - double n02=0; - double rc2=0; - double beta=0; - */ - double bkg_level=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("lx_dbeta_param.txt"); - //output the datasets and fitting results - for(size_t i=0;i& pj=dynamic_cast&>(f.get_model()); - pj.attach_cfunc(cf_erg); - - mv=f.eval_model_raw(radii,p); - double flux_erg=0; - for(size_t i=0;i -#include -#include -#include -#include - -using std::cerr; -using std::endl; - -namespace opt_utilities -{ - static const int display_interval=10; - /** - \brief chi-square statistic - \tparam Ty the return type of model - \tparam Tx the type of the self-var - \tparam Tp the type of model parameter - \tparam Ts the type of the statistic - \tparam Tstr the type of the string used - */ - template - class chisq - :public statistic - { - }; - template<> - class chisq,double,std::string> - :public statistic ,double,std::string> - { - public: - typedef double Ty; - typedef double Tx; - typedef std::vector Tp; - typedef double Ts; - typedef std::string Tstr; - private: - bool verb; - bool limit_bound; - int n; - - statistic* do_clone()const - { - // return const_cast*>(this); - return new chisq(*this); - } - - const char* do_get_type_name()const - { - return "chi^2 statistics (specialized for double)"; - } - public: - void verbose(bool v) - { - verb=v; - } - - void set_limit() - { - limit_bound=true; - } - - void clear_limit() - { - limit_bound=false; - } - public: - chisq() - :verb(true),limit_bound(false) - {} - - Ty do_eval(const Tp& p) - { - if(limit_bound) - { - Tp p1=this->get_fitter().get_model().reform_param(p); - for(size_t i=0;ithis->get_fitter().get_param_info(i).get_upper_limit()|| - p1[i]get_fitter().get_param_info(i).get_lower_limit()) - { - return 1e99; - } - } - } - Ty result(0); - std::vector vx; - std::vector vy; - std::vector vye1; - std::vector vye2; - std::vector my; - float xmin=1e99,xmax=-1e99,ymin=1e99,ymax=-1e99; - if(verb) - { - n++; - if(n%display_interval==0) - { - vx.resize(this->get_data_set().size()); - vy.resize(this->get_data_set().size()); - vye1.resize(this->get_data_set().size()); - vye2.resize(this->get_data_set().size()); - my.resize(this->get_data_set().size()); - } - } - - for(int i=(this->get_data_set()).size()-1;i>=0;--i) - { - -#ifdef HAVE_X_ERROR - Tx x1=this->get_data_set().get_data(i).get_x()-this->get_data_set().get_data(i).get_x_lower_err(); - Tx x2=this->get_data_set().get_data(i).get_x()+this->get_data_set().get_data(i).get_x_upper_err(); - Tx x=this->get_data_set().get_data(i).get_x(); - Ty errx1=(this->eval_model(x1,p)-this->eval_model(x,p)); - Ty errx2=(this->eval_model(x2,p)-this->eval_model(x,p)); - //Ty errx=0; -#else - Ty errx1=0; - Ty errx2=0; -#endif - - Ty y_model=this->eval_model(this->get_data_set().get_data(i).get_x(),p); - Ty y_obs=this->get_data_set().get_data(i).get_y(); - Ty y_err; - - Ty errx=0; - if(errx10?errx1:-errx1; - } - else - { - errx=errx2>0?errx2:-errx2; - } - } - else - { - if(y_obs0?errx2:-errx2; - } - else - { - errx=errx1>0?errx1:-errx1; - } - } - - if(y_model>y_obs) - { - y_err=this->get_data_set().get_data(i).get_y_upper_err(); - } - else - { - y_err=this->get_data_set().get_data(i).get_y_lower_err(); - } - - Ty chi=(y_obs-y_model)/std::sqrt(y_err*y_err+errx*errx); - - result+=chi*chi; - - if(verb&&n%display_interval==0) - { - vx.at(i)=this->get_data_set().get_data(i).get_x(); - vy.at(i)=this->get_data_set().get_data(i).get_y(); - vye1.at(i)=std::abs(this->get_data_set().get_data(i).get_y_lower_err()); - vye2.at(i)=std::abs(this->get_data_set().get_data(i).get_y_upper_err()); - my.at(i)=y_model; - - xmin=std::min(vx.at(i),xmin); - ymin=std::min(vy.at(i),ymin-vye1[i]); - xmax=std::max(vx.at(i),xmax); - ymax=std::max(vy.at(i),ymax+vye2[i]); - } - - } - if(verb) - { - if(n%display_interval==0) - { - cerr< - class dbeta - :public model,std::vector,std::vector > - { - public: - dbeta() - { - this->push_param_info(param_info,std::string>("n01",1)); - this->push_param_info(param_info,std::string>("beta1",.66)); - this->push_param_info(param_info,std::string>("rc1",100)); - - this->push_param_info(param_info,std::string>("n02",1)); - this->push_param_info(param_info,std::string>("beta2",.67)); - this->push_param_info(param_info,std::string>("rc2",110)); - - } - - public: - dbeta* do_clone()const - { - return new dbeta(*this); - } - - std::vector do_eval(const std::vector & x, - const std::vector& p) - { - T n01=std::abs(p[0]); - T beta1=p[1]; - T rc1=p[2]; - - T n02=std::abs(p[3]); - T beta2=p[4]; - T rc2=p[5]; - - - - std::vector result(x.size()-1); - for(size_t i=1;i - class dbeta2 - :public model,std::vector,std::vector > - { - public: - dbeta2() - { - this->push_param_info(param_info,std::string>("n01",1)); - this->push_param_info(param_info,std::string>("rc1",100)); - this->push_param_info(param_info,std::string>("n02",1)); - this->push_param_info(param_info,std::string>("rc2",110)); - this->push_param_info(param_info,std::string>("beta",.67)); - - } - - public: - dbeta2* do_clone()const - { - return new dbeta2(*this); - } - - std::vector do_eval(const std::vector & x, - const std::vector& p) - { - T n01=std::abs(p[0]); - T rc1=p[1]; - - T n02=std::abs(p[2]); - T rc2=p[3]; - T beta=p[4]; - T beta1=beta; - T beta2=beta; - - std::vector result(x.size()-1); - for(size_t i=1;i,double,std::string>& f,double cm_per_pixel,const std::vector& r,const std::vector& y,const std::vector& ye) - { - os<<"read serr 1 2"< p=f.get_all_params(); - for(size_t i=1;i,std::vector,std::vector,double,std::string>& f,double cm_per_pixel,const std::vector& r,const std::vector& sbps,const std::vector& sbpe) - { - os<<"read serr 1 2"< p=f.get_all_params(); - std::vector mv=f.eval_model_raw(r,p); - for(size_t i=1;i -#include -#include -#include - -namespace opt_utilities -{ - void dump_sbp_beta(std::ostream& os,fitter,double,std::string>& f,double cm_per_pixel,const std::vector& r,const std::vector& y,const std::vector& ye); - void dump_rho_beta(std::ostream& os,fitter,std::vector,std::vector,double,std::string>& f,double cm_per_pixel,const std::vector& r,const std::vector& sbps,const std::vector& sbpe); - void dump_rho_dbeta(std::ostream& os,fitter,std::vector,std::vector,double,std::string>& f,double cm_per_pixel); -}; - -#endif diff --git a/mass_profile/fit_beta_sbp.cpp b/mass_profile/fit_beta_sbp.cpp deleted file mode 100644 index 295fa1e..0000000 --- a/mass_profile/fit_beta_sbp.cpp +++ /dev/null @@ -1,534 +0,0 @@ -/* - 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"< -#include -#include -#include -#include "vchisq.hpp" -#include "dbeta.hpp" -#include "beta_cfg.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 dbeta_func(double r, double n01, double rc1, double beta1, - double n02, double rc2, double beta2) -{ - double v1 = abs(n01) * pow(1+r*r/rc1/rc1, -3./2.*abs(beta1)); - double v2 = abs(n02) * pow(1+r*r/rc2/rc2, -3./2.*abs(beta2)); - return v1 + v2; -} - -//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 - 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()) - { - break; - } - cf.add_point(x,y); - } - cf.gen_spline(); - - //read temperature profile data - 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; - 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_level=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(size_t i=0;i -#include -#include -#include "chisq.hpp" -#include -#include -#include -#include -#include - -using namespace opt_utilities; -using namespace std; -const double cm=1; -const double kpc=3.08568e+21*cm; -const double pi=4*atan(1); -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; -} - - -int main(int argc,char* argv[]) -{ - if(argc<3) - { - cerr<<"Usage:"< [rmin in kpc]"<=4) - { - rmin_kpc=atof(argv[3]); - } - double z=0; - z=atof(argv[2]); - //define the fitter - fitter,double,std::string> fit; - //define the data set - default_data_set ds; - //open the data file - ifstream ifs(argv[1]); - //cout<<"read serr 2"<>r>>re>>m>>me; - if(!ifs.good()) - { - break; - } - if(r d(r,m,me,me,re,re); - ofs_fit_result< >()); - //use chi^2 statistic - fit.set_statistic(chisq,double,std::string>()); - fit.set_model(nfw()); - //fit.set_param_value("rs",4); - //fit.set_param_value("rho0",100); - fit.fit(); - fit.fit(); - vector p=fit.fit(); - //output parameters - ofstream ofs_param("nfw_param.txt"); - for(size_t i=0;i d=ds.get_data(i); - double x=d.get_x(); - double y=d.get_y(); - double ye=d.get_y_lower_err(); - double ym=fit.eval_model(x,p); - ofs_fit_result< -#include -#include -#include "chisq.hpp" -#include -#include -#include -#include -#include -#include - -using namespace opt_utilities; -using namespace std; -const double cm=1; -const double kpc=3.08568e+21*cm; - -int main(int argc,char* argv[]) -{ - if(argc<2) - { - cerr<<"Usage:"< [param file] [cm per pixel]"<=4) - { - cm_per_pixel=atof(argv[3]); - } - - //define the fitter - fitter,double,std::string> fit; - //define the data set - default_data_set ds; - //open the data file - ifstream ifs(argv[1]); - double min_r=1e9; - //cout<<"read serr 2"<0) - { - ofs_fit_result<<"la x radius (kpc)"<>r>>re>>t>>te; - if(!ifs.good()) - { - break; - } - min_r=min(r,min_r); - data d(r,t,te,te,re,re); - //std::cerr<0) - { - ofs_fit_result< >()); - //use chi^2 statistic - chisq,double,std::string> chisq_object; - chisq_object.set_limit(); - fit.set_statistic(chisq_object); - //fit.set_statistic(chisq,double,std::string>()); - fit.set_model(wang2012_model()); - - if(argc>=3&&std::string(argv[2])!="NONE") - { - std::vector freeze_list; - ifstream ifs_param(argv[2]); - assert(ifs_param.is_open()); - for(;;) - { - string pname; - double pvalue; - double lower,upper; - char param_status; - ifs_param>>pname>>pvalue>>lower>>upper>>param_status; - if(!ifs_param.good()) - { - break; - } - if(param_status=='F') - { - freeze_list.push_back(pname); - } - if(pvalue<=lower||pvalue>=upper) - { - cerr<<"Invalid initial value, central value not enclosed by the lower and upper boundaries, adjust automatically"<,std::string> fp(freeze_list[0]); - fit.set_param_modifier(fp); - for(size_t i=1;i,std::string>&>(fit.get_param_modifier())+=freeze_param,std::string>(freeze_list[i]); - } - } - } - - for(int i=0;i<100;++i) - { - fit.fit(); - } - vector p=fit.fit(); -#if 0 - ofstream output_param; - if(argc>=3&&std::string(argv[2])!="NONE") - { - output_param.open(argv[2]); - } - else - { - output_param.open("para0.txt"); - } -#endif - //output parameters - for(size_t i=0;i=3&&std::string(argv[2])!="NONE") -#if 0 - { - output_param<0) - { - ofs_fit_result< -#include - -namespace opt_utilities -{ - template - class nfw - :public model,std::string> - { - private: - model >* do_clone()const - { - return new nfw(*this); - } - - const char* do_get_type_name()const - { - return "1d power law"; - } - public: - nfw() - { - this->push_param_info(param_info >("rho0",1,0,1e99)); - this->push_param_info(param_info >("rs",100,0,1e99)); - } - - T do_eval(const T& r,const std::vector& param) - { - T rho0=std::abs(param[0]); - T rs=std::abs(param[1]); - static const T pi=4*std::atan(1); - return 4*pi*rho0*rs*rs*rs*(std::log((r+rs)/rs)-r/(r+rs)); - } - - private: - std::string do_get_information()const - { - return ""; - } - }; -} - - - -#endif -//EOF diff --git a/mass_profile/projector.hpp b/mass_profile/projector.hpp deleted file mode 100644 index 4ef4b32..0000000 --- a/mass_profile/projector.hpp +++ /dev/null @@ -1,214 +0,0 @@ -#ifndef PROJ_HPP -#define PROJ_HPP -/* - Defining the class that is used to consider the projection effect - Author: Junhua Gu - Last modified: 2011.01.01 -*/ - - -#include -#include -#include - -static const double pi=4*atan(1); -// Ratio of the electron density (n_e) to the proton density (n_p) -// n_gas = n_e + n_p ~= 1.826 n_e => n_e / n_p ~= 1.21 -// Reference: Ettori et al. 2013, Space Sci. Rev., 177, 119-154; Eq.(9) below -static const double ne_np_ratio = 1.21; - -namespace opt_utilities -{ - //This is used to project a 3-D surface brightness model to 2-D profile - template - class projector - :public model,std::vector,std::vector > - { - private: - //Points to a 3-D model that is to be projected - model,std::vector,std::vector >* pmodel; - func_obj* pcfunc; - T cm_per_pixel; - public: - //default cstr - projector() - :pmodel(NULL_PTR),pcfunc(NULL_PTR),cm_per_pixel(1) - {} - //copy cstr - projector(const projector& rhs) - :cm_per_pixel(rhs.cm_per_pixel) - { - attach_model(*(rhs.pmodel)); - if(rhs.pcfunc) - { - pcfunc=rhs.pcfunc->clone(); - } - else - { - pcfunc=NULL_PTR; - } - } - //assign operator - projector& operator=(const projector& rhs) - { - cm_per_pixel=rhs.cm_per_pixel; - if(pmodel) - { - pmodel->destroy(); - } - if(pcfunc) - { - pcfunc->destroy(); - } - if(rhs.pcfunc) - { - pcfunc=rhs.pcfunc->clone(); - } - if(rhs.pmodel) - { - pmodel=rhs.pmodel->clone(); - } - } - //destr - ~projector() - { - if(pmodel) - { - pmodel->destroy(); - } - if(pcfunc) - { - pcfunc->destroy(); - } - } - //used to clone self - model,std::vector,std::vector >* - do_clone()const - { - return new projector(*this); - } - - public: - void set_cm_per_pixel(const T& x) - { - cm_per_pixel=x; - } - - //attach the model that is to be projected - void attach_model(const model,std::vector,std::vector >& m) - { - this->clear_param_info(); - for(size_t i=0;ipush_param_info(m.get_param_info(i)); - } - this -> push_param_info(param_info, - std::string>("bkg",0,0,1E99)); - pmodel=m.clone(); - pmodel->clear_param_modifier(); - } - - void attach_cfunc(const func_obj& cf) - { - if(pcfunc) - { - pcfunc->destroy(); - } - pcfunc=cf.clone(); - } - - public: - //calc the volume - /* - This is a sphere that is subtracted by a cycline. - /| |\ - / | | \ - | | | | - | | | | - \ | | / - \| |/ - */ - T calc_v_ring(T rsph,T rcyc) - { - if(rcyc& rlist,int nsph,int nrad) - { - if(nsph& p)const - { - std::vector p1(this->reform_param(p)); - for(size_t i=0;i!=p1.size();++i) - { - if(get_element(p1,i)>this->get_param_info(i).get_upper_limit()|| - get_element(p1,i)get_param_info(i).get_lower_limit()) - { - // std::cerr<get_param_info(i).get_name()<<"\t"< p2(p1.size()-1); - for(size_t i=0;imeets_constraint(p2); - } - public: - //Perform the projection - std::vector do_eval(const std::vector& x,const std::vector& p) - { - T bkg=std::abs(p.back()); - //I think following codes are clear enough :). - std::vector unprojected(pmodel->eval(x,p)); - std::vector projected(unprojected.size()); - - for(size_t nrad=0; nrad - -using namespace std; -void report_error(const char* message) -{ - cerr<<"MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM"< -#include -#include -#include -#include - -template -class spline -{ -public: - std::vector x_list; - std::vector y_list; - std::vector y2_list; - -public: - void push_point(T x,T y) - { - if(!x_list.empty()) - { - assert(x>*(x_list.end()-1)); - } - x_list.push_back(x); - y_list.push_back(y); - } - - T get_value(T x) - { - if(x<=x_list[0]) - { - return y_list[0]; - } - if(x>=x_list.back()) - { - return y_list.back(); - } - assert(x_list.size()==y2_list.size()); - assert(x>x_list[0]); - assert(xx) - { - n2--; - } - } - T h=x_list[n2]-x_list[n1]; - double a=(x_list[n2]-x)/h; - double b=(x-x_list[n1])/h; - return a*y_list[n1]+b*y_list[n2]+((a*a*a-a)*y2_list[n1]+ - (b*b*b-b)*y2_list[n2])*(h*h)/6.; - - } - - void gen_spline(T y2_0,T y2_N) - { - int n=x_list.size(); - y2_list.resize(0); - y2_list.resize(x_list.size()); - std::vector u(x_list.size()); - if(std::abs(y2_0)::epsilon()) - { - y2_list[0]=0; - u[0]=0; - } - else - { - y2_list[0]=-.5; - u[0]=(3./(x_list[1]-x_list[0]))*((y_list[1]-y_list[0])/(x_list[1]-x_list[0])-y2_0); - } - for(int i=1;i::epsilon()) - { - qn=un=0; - } - else - { - qn=.5; - un=(3./(x_list[n-1]-x_list[n-2]))*(y2_N-(y_list[n-1]-y_list[n-2])/(x_list[n-1]-x_list[n-2])); - - } - y2_list[n-1]=(un-qn*u[n-2])/(qn*y2_list[n-2]+1.); - for(int i=n-2;i>=0;--i) - { - y2_list[i]=y2_list[i]*y2_list[i+1]+u[i]; - } - } - -}; - -#endif /* SPLINE_HPP */ diff --git a/mass_profile/vchisq.hpp b/mass_profile/vchisq.hpp deleted file mode 100644 index 8cd8481..0000000 --- a/mass_profile/vchisq.hpp +++ /dev/null @@ -1,137 +0,0 @@ -/** - \file vchisq.hpp - \brief chi-square statistic - \author Junhua Gu - */ - -#ifndef VCHI_SQ_HPP -#define VCHI_SQ_HPP - -#define OPT_HEADER - -#include -#include -#include -#include -#include - -using std::cerr; -using std::endl; - -namespace opt_utilities -{ - template - class vchisq - :public statistic,std::vector,std::vector,T,std::string> - { - private: - bool verb; - int n; - bool limit_bound; - typedef std::vector Tp; - - vchisq* do_clone()const - { - return new vchisq(*this); - } - - const char* do_get_type_name()const - { - return "chi^2 statistic"; - } - - public: - void verbose(bool v) - { - verb=v; - } - - void set_limit() - { - limit_bound=true; - } - - void clear_limit() - { - limit_bound=false; - } - public: - vchisq() - :verb(false),limit_bound(false) - {} - - T do_eval(const std::vector& p) - { - if(limit_bound) - { - if(!this->get_fitter().get_model().meets_constraint(p)) - { - return 1e99; - } - } - T result(0); - std::vector vx; - std::vector vy; - std::vector vye; - std::vector my; - float x1=1e99,x2=-1e99,y1=1e99,y2=-1e99; - if(verb) - { - n++; - if(n%100==0) - { - vx.resize(this->get_data_set().get_data(0).get_y().size()); - vy.resize(this->get_data_set().get_data(0).get_y().size()); - vye.resize(this->get_data_set().get_data(0).get_y().size()); - my.resize(this->get_data_set().get_data(0).get_y().size()); - } - } - for(int i=(this->get_data_set()).size()-1;i>=0;--i) - { - const std::vector y_model(this->eval_model(this->get_data_set().get_data(i).get_x(),p)); - const std::vector& y=this->get_data_set().get_data(i).get_y(); - const std::vector& ye=this->get_data_set().get_data(i).get_y_lower_err(); - for(size_t j=0;jget_data_set().get_data(i).get_x().at(j)+this->get_data_set().get_data(i).get_x().at(j+1))/2.); - vy.at(j)=(y[j]); - vye.at(j)=ye[j]; - my.at(j)=(y_model[j]); - x1=std::min(vx.at(j),x1); - y1=std::min(vy.at(j),y1); - x2=std::max(vx.at(j),x2); - y2=std::max(vy.at(j),y2); - vye[j]=log10(vy[j]+vye[j])-log10(vy[j]); - vx[j]=log10(vx[j]); - vy[j]=log10(vy[j]); - my[j]=log10(my[j]); - } - } - } - if(verb) - { - if(n%100==0) - { - cerr< -#include - -namespace opt_utilities -{ - template - class wang2012_model - :public model,std::string> - { - private: - model >* do_clone()const - { - return new wang2012_model(*this); - } - - const char* do_get_type_name()const - { - return "1d power law"; - } - public: - wang2012_model() - { - this->push_param_info(param_info >("A",5,0,500)); - this->push_param_info(param_info >("n",1.66,0,10)); - this->push_param_info(param_info >("xi",0.45,0,1)); - this->push_param_info(param_info >("a2",1500,0,1e8)); - this->push_param_info(param_info >("a3",50,0,1e8)); - this->push_param_info(param_info >("beta",0.49,0.1,0.7)); - this->push_param_info(param_info >("T0",0,0,10)); - - } - - T do_eval(const T& x,const std::vector& param) - { - T A=param[0]; - T n=param[1]; - T xi=param[2]; - T a2=param[3]; - T a3=param[4]; - T beta=param[5]; - T T0=param[6]; - return A*(pow(x,n)+xi*a2)/(pow(x,n)+a2)/pow(1+x*x/a3/a3,beta)+T0; - //return A*(pow(x,n)+a1)/(pow(x,n)+1)/pow(1+x*x/a3/a3,beta)+T0; - } - - private: - std::string do_get_information()const - { - return ""; - } - }; -} - - - -#endif -//EOF diff --git a/src/Makefile b/src/Makefile new file mode 100644 index 0000000..76ddeaa --- /dev/null +++ b/src/Makefile @@ -0,0 +1,98 @@ +# +# Makefile for `chandra-acis-analysis/mass_profile` tools +# +# Junhua GU +# Weitian LI +# 2016-06-07 +# + + +CXX ?= g++ +CXXFLAGS += -Wall +CXXFLAGS += -Werror + +ifdef OPENMP + CXXFLAGS += -fopenmp +endif + +ifdef DEBUG + CXXFLAGS += -g +else + CXXFLAGS += -O2 +endif + +OPT_UTIL_INC ?= -I../opt_utilities + +TARGETS= fit_dbeta_sbp fit_beta_sbp fit_wang2012_model \ + fit_nfw_mass calc_lx_dbeta calc_lx_beta +HEADERS= projector.hpp spline.hpp vchisq.hpp + +all: $(TARGETS) + +# NOTE: +# Object/source files should placed *before* libraries (order matters) + +fit_dbeta_sbp: fit_dbeta_sbp.o beta_cfg.o report_error.o + $(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC) + +fit_beta_sbp: fit_beta_sbp.o beta_cfg.o dump_fit_qdp.o + $(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC) + +fit_wang2012_model: fit_wang2012_model.o + $(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC) + +fit_nfw_mass: fit_nfw_mass.o + $(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC) + +calc_lx_dbeta: calc_lx_dbeta.o beta_cfg.o report_error.o + $(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC) + +calc_lx_beta: calc_lx_beta.o beta_cfg.o dump_fit_qdp.o + $(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC) + + +fit_dbeta_sbp.o: fit_dbeta_sbp.cpp $(HEADERS) + $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC) + +fit_beta_sbp.o: fit_beta_sbp.cpp beta.hpp $(HEADERS) + $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC) + +fit_wang2012_model.o: fit_wang2012_model.cpp wang2012_model.hpp chisq.hpp + $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC) + +fit_nfw_mass.o: fit_nfw_mass.cpp nfw.hpp chisq.hpp + $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC) + +calc_lx_dbeta.o: calc_lx_dbeta.cpp $(HEADERS) + $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC) + +calc_lx_beta.o: calc_lx_beta.cpp beta.hpp $(HEADERS) + $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC) + +beta_cfg.o: beta_cfg.cpp beta_cfg.hpp + $(CXX) $(CXXFLAGS) -c $< + +report_error.o: report_error.cpp report_error.hpp + $(CXX) $(CXXFLAGS) -c $< + +dump_fit_qdp.o: dump_fit_qdp.cpp dump_fit_qdp.hpp + $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC) + +%.o: %.cpp + $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC) + + +clean: + rm -f *.o $(TARGETS) + + +install: $(TARGETS) + @for f in $(TARGETS); do \ + (cd ../bin && ln -svf ../src/$$f . ); \ + done + + +uninstall: + @for f in $(TARGETS); do \ + rm -fv ../bin/$$f; \ + done diff --git a/src/README.md b/src/README.md new file mode 100644 index 0000000..0c61821 --- /dev/null +++ b/src/README.md @@ -0,0 +1,37 @@ +Mass/Luminosity/Flux Calculation Tools +============================== + +Junhua GU, Weitian LI, Zhenghao ZHU + + +Introduction +------------ +This directory contains the tools/scripts used to fit the temperature +profile and surface brightness file, to calculate the gas density profile +and gravitational mass profile, to calculate luminosity and flux and +other related quantities. + + +NOTE +---- +* Mass calculation references: + + Walker et al. 2012, MNRAS, 422, 3503-3515, + [ADS:2012MNRAS.422.3503W](http://adsabs.harvard.edu/abs/2012MNRAS.422.3503W) + + Ettori et al. 2013, Space Science Reviews, 177, 119-154, + [ADS:2013SSRv..177..119E](http://adsabs.harvard.edu/abs/2013SSRv..177..119E) +* Errors are calculated at 68% confident level. +* NFW profile is employed to extrapolate the mass profile. +* We use a self-proposed model (see ``wang2012_model.hpp``) to describe/fit + the gas temperature profile. +* The single-beta and double-beta models (with a constant background) are used + to describe the gas density profile. + + +TODO +---- +* Merge the scripts/tools for single-beta and double-beta SBP models + (to reduce code duplications) +* Uncertainties/errors calculation **maybe** inappropriate/problematic + (e.g., ``analyze_mass_profile.py``); + why not just use the quantile-based (e.g., Q84.15-Q15.85 ~= 68.3%) + uncertainties or standard deviation??? diff --git a/src/beta.hpp b/src/beta.hpp new file mode 100644 index 0000000..6ae70ac --- /dev/null +++ b/src/beta.hpp @@ -0,0 +1,45 @@ +#ifndef BETA +#define BETA +#include "projector.hpp" + +namespace opt_utilities +{ + template + class beta + :public model,std::vector,std::vector > + { + public: + beta() + { + this->push_param_info(param_info,std::string>("n0",1,0,1E99)); + this->push_param_info(param_info,std::string>("beta",.66,0,1E99)); + this->push_param_info(param_info,std::string>("rc",100,0,1E99)); + } + + public: + beta* do_clone()const + { + return new beta(*this); + } + + std::vector do_eval(const std::vector & x, + const std::vector& p) + { + T n0=std::abs(p[0]); + T beta=p[1]; + T rc=p[2]; + + std::vector result(x.size()-1); + for(size_t i=1;i +#include +using namespace std; + +cfg_map parse_cfg_file(std::istream& is) +{ + cfg_map result; + result.rmin_pixel=-1; + result.rmin_kpc=-1; + for(;;) + { + std::string line; + getline(is,line); + line+="\n"; + if(!is.good()) + { + break; + } + string key; + istringstream iss(line); + iss>>key; + if(key=="sbp_data") + { + string value; + iss>>value; + result.sbp_data=value; + } + else if(key=="cfunc_profile") + { + string value; + iss>>value; + result.cfunc_profile=value; + } + else if(key=="tprofile") + { + string value; + iss>>value; + result.tprofile=value; + } + else if(key=="z") + { + double z; + iss>>z; + result.z=z; + } + else if(key=="cm_per_pixel") + { + double cm_per_pixel; + iss>>cm_per_pixel; + result.cm_per_pixel=cm_per_pixel; + } + else if(key=="rmin_pixel") + { + double v; + iss>>v; + result.rmin_pixel=v; + } + else if(key=="rmin_kpc") + { + double v; + iss>>v; + result.rmin_kpc=v; + } + else + { + std::vector value; + for(;;) + { + double v; + iss>>v; + if(!iss.good()) + { + break; + } + value.push_back(v); + } + if(!value.empty()) + { + result.param_map[key]=value; + } + } + } + return result; +} diff --git a/src/beta_cfg.hpp b/src/beta_cfg.hpp new file mode 100644 index 0000000..27148df --- /dev/null +++ b/src/beta_cfg.hpp @@ -0,0 +1,23 @@ +#ifndef BETA_CFG +#define BETA_CFG + +#include +#include +#include +#include + +struct cfg_map +{ + std::string sbp_data; + std::string cfunc_profile; + std::string tprofile; + double z; + double cm_per_pixel; + double rmin_kpc; + double rmin_pixel; + std::map > param_map; +}; + +cfg_map parse_cfg_file(std::istream& is); + +#endif diff --git a/src/calc_lx_beta.cpp b/src/calc_lx_beta.cpp new file mode 100644 index 0000000..11d7b04 --- /dev/null +++ b/src/calc_lx_beta.cpp @@ -0,0 +1,497 @@ +/** + * Calculate the total luminosity and flux within the specified radius. + * + * Base on 'fit_beta_sbp.cpp' and supersede 'calc_lx.cpp' + * + * Author: Junhua Gu + */ + +#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)); +} + +//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<4) + { + cerr< [cfunc2_erg ...]"< 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 + 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()) + { + break; + } + cf.add_point(x,y); + } + cf.gen_spline(); + + //read temperature profile data + 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=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("lx_beta_param.txt"); + for(size_t i=0;i=1) + { + ofs_sbp<<"line on 2"<=1) + { + ofs_sbp<<"no no no"<& pj=dynamic_cast&>(f.get_model()); + pj.attach_cfunc(cf_erg); + + + + mv=f.eval_model_raw(radii,p); + double flux_erg=0; + for(size_t i=0;i +#include +#include +#include +#include "vchisq.hpp" +#include "dbeta.hpp" +#include "beta_cfg.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 dbeta_func(double r, double n01, double rc1, double beta1, + double n02, double rc2, double beta2) +{ + double v1 = abs(n01) * pow(1+r*r/rc1/rc1, -3./2.*abs(beta1)); + double v2 = abs(n02) * pow(1+r*r/rc2/rc2, -3./2.*abs(beta2)); + return v1 + v2; +} + + +//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<4) + { + cerr< [cfunc2_erg ...]"< 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 + 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()) + { + break; + } + cf.add_point(x,y); + } + cf.gen_spline(); + + //read temperature profile data + 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; + 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 =0; + double rc1=0; + double n02=0; + double rc2=0; + double beta=0; + */ + double bkg_level=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("lx_dbeta_param.txt"); + //output the datasets and fitting results + for(size_t i=0;i& pj=dynamic_cast&>(f.get_model()); + pj.attach_cfunc(cf_erg); + + mv=f.eval_model_raw(radii,p); + double flux_erg=0; + for(size_t i=0;i +#include +#include +#include +#include + +using std::cerr; +using std::endl; + +namespace opt_utilities +{ + static const int display_interval=10; + /** + \brief chi-square statistic + \tparam Ty the return type of model + \tparam Tx the type of the self-var + \tparam Tp the type of model parameter + \tparam Ts the type of the statistic + \tparam Tstr the type of the string used + */ + template + class chisq + :public statistic + { + }; + template<> + class chisq,double,std::string> + :public statistic ,double,std::string> + { + public: + typedef double Ty; + typedef double Tx; + typedef std::vector Tp; + typedef double Ts; + typedef std::string Tstr; + private: + bool verb; + bool limit_bound; + int n; + + statistic* do_clone()const + { + // return const_cast*>(this); + return new chisq(*this); + } + + const char* do_get_type_name()const + { + return "chi^2 statistics (specialized for double)"; + } + public: + void verbose(bool v) + { + verb=v; + } + + void set_limit() + { + limit_bound=true; + } + + void clear_limit() + { + limit_bound=false; + } + public: + chisq() + :verb(true),limit_bound(false) + {} + + Ty do_eval(const Tp& p) + { + if(limit_bound) + { + Tp p1=this->get_fitter().get_model().reform_param(p); + for(size_t i=0;ithis->get_fitter().get_param_info(i).get_upper_limit()|| + p1[i]get_fitter().get_param_info(i).get_lower_limit()) + { + return 1e99; + } + } + } + Ty result(0); + std::vector vx; + std::vector vy; + std::vector vye1; + std::vector vye2; + std::vector my; + float xmin=1e99,xmax=-1e99,ymin=1e99,ymax=-1e99; + if(verb) + { + n++; + if(n%display_interval==0) + { + vx.resize(this->get_data_set().size()); + vy.resize(this->get_data_set().size()); + vye1.resize(this->get_data_set().size()); + vye2.resize(this->get_data_set().size()); + my.resize(this->get_data_set().size()); + } + } + + for(int i=(this->get_data_set()).size()-1;i>=0;--i) + { + +#ifdef HAVE_X_ERROR + Tx x1=this->get_data_set().get_data(i).get_x()-this->get_data_set().get_data(i).get_x_lower_err(); + Tx x2=this->get_data_set().get_data(i).get_x()+this->get_data_set().get_data(i).get_x_upper_err(); + Tx x=this->get_data_set().get_data(i).get_x(); + Ty errx1=(this->eval_model(x1,p)-this->eval_model(x,p)); + Ty errx2=(this->eval_model(x2,p)-this->eval_model(x,p)); + //Ty errx=0; +#else + Ty errx1=0; + Ty errx2=0; +#endif + + Ty y_model=this->eval_model(this->get_data_set().get_data(i).get_x(),p); + Ty y_obs=this->get_data_set().get_data(i).get_y(); + Ty y_err; + + Ty errx=0; + if(errx10?errx1:-errx1; + } + else + { + errx=errx2>0?errx2:-errx2; + } + } + else + { + if(y_obs0?errx2:-errx2; + } + else + { + errx=errx1>0?errx1:-errx1; + } + } + + if(y_model>y_obs) + { + y_err=this->get_data_set().get_data(i).get_y_upper_err(); + } + else + { + y_err=this->get_data_set().get_data(i).get_y_lower_err(); + } + + Ty chi=(y_obs-y_model)/std::sqrt(y_err*y_err+errx*errx); + + result+=chi*chi; + + if(verb&&n%display_interval==0) + { + vx.at(i)=this->get_data_set().get_data(i).get_x(); + vy.at(i)=this->get_data_set().get_data(i).get_y(); + vye1.at(i)=std::abs(this->get_data_set().get_data(i).get_y_lower_err()); + vye2.at(i)=std::abs(this->get_data_set().get_data(i).get_y_upper_err()); + my.at(i)=y_model; + + xmin=std::min(vx.at(i),xmin); + ymin=std::min(vy.at(i),ymin-vye1[i]); + xmax=std::max(vx.at(i),xmax); + ymax=std::max(vy.at(i),ymax+vye2[i]); + } + + } + if(verb) + { + if(n%display_interval==0) + { + cerr< + class dbeta + :public model,std::vector,std::vector > + { + public: + dbeta() + { + this->push_param_info(param_info,std::string>("n01",1)); + this->push_param_info(param_info,std::string>("beta1",.66)); + this->push_param_info(param_info,std::string>("rc1",100)); + + this->push_param_info(param_info,std::string>("n02",1)); + this->push_param_info(param_info,std::string>("beta2",.67)); + this->push_param_info(param_info,std::string>("rc2",110)); + + } + + public: + dbeta* do_clone()const + { + return new dbeta(*this); + } + + std::vector do_eval(const std::vector & x, + const std::vector& p) + { + T n01=std::abs(p[0]); + T beta1=p[1]; + T rc1=p[2]; + + T n02=std::abs(p[3]); + T beta2=p[4]; + T rc2=p[5]; + + + + std::vector result(x.size()-1); + for(size_t i=1;i + class dbeta2 + :public model,std::vector,std::vector > + { + public: + dbeta2() + { + this->push_param_info(param_info,std::string>("n01",1)); + this->push_param_info(param_info,std::string>("rc1",100)); + this->push_param_info(param_info,std::string>("n02",1)); + this->push_param_info(param_info,std::string>("rc2",110)); + this->push_param_info(param_info,std::string>("beta",.67)); + + } + + public: + dbeta2* do_clone()const + { + return new dbeta2(*this); + } + + std::vector do_eval(const std::vector & x, + const std::vector& p) + { + T n01=std::abs(p[0]); + T rc1=p[1]; + + T n02=std::abs(p[2]); + T rc2=p[3]; + T beta=p[4]; + T beta1=beta; + T beta2=beta; + + std::vector result(x.size()-1); + for(size_t i=1;i,double,std::string>& f,double cm_per_pixel,const std::vector& r,const std::vector& y,const std::vector& ye) + { + os<<"read serr 1 2"< p=f.get_all_params(); + for(size_t i=1;i,std::vector,std::vector,double,std::string>& f,double cm_per_pixel,const std::vector& r,const std::vector& sbps,const std::vector& sbpe) + { + os<<"read serr 1 2"< p=f.get_all_params(); + std::vector mv=f.eval_model_raw(r,p); + for(size_t i=1;i +#include +#include +#include + +namespace opt_utilities +{ + void dump_sbp_beta(std::ostream& os,fitter,double,std::string>& f,double cm_per_pixel,const std::vector& r,const std::vector& y,const std::vector& ye); + void dump_rho_beta(std::ostream& os,fitter,std::vector,std::vector,double,std::string>& f,double cm_per_pixel,const std::vector& r,const std::vector& sbps,const std::vector& sbpe); + void dump_rho_dbeta(std::ostream& os,fitter,std::vector,std::vector,double,std::string>& f,double cm_per_pixel); +}; + +#endif 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"< +#include +#include +#include +#include "vchisq.hpp" +#include "dbeta.hpp" +#include "beta_cfg.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 dbeta_func(double r, double n01, double rc1, double beta1, + double n02, double rc2, double beta2) +{ + double v1 = abs(n01) * pow(1+r*r/rc1/rc1, -3./2.*abs(beta1)); + double v2 = abs(n02) * pow(1+r*r/rc2/rc2, -3./2.*abs(beta2)); + return v1 + v2; +} + +//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 + 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()) + { + break; + } + cf.add_point(x,y); + } + cf.gen_spline(); + + //read temperature profile data + 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; + 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_level=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(size_t i=0;i +#include +#include +#include "chisq.hpp" +#include +#include +#include +#include +#include + +using namespace opt_utilities; +using namespace std; +const double cm=1; +const double kpc=3.08568e+21*cm; +const double pi=4*atan(1); +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; +} + + +int main(int argc,char* argv[]) +{ + if(argc<3) + { + cerr<<"Usage:"< [rmin in kpc]"<=4) + { + rmin_kpc=atof(argv[3]); + } + double z=0; + z=atof(argv[2]); + //define the fitter + fitter,double,std::string> fit; + //define the data set + default_data_set ds; + //open the data file + ifstream ifs(argv[1]); + //cout<<"read serr 2"<>r>>re>>m>>me; + if(!ifs.good()) + { + break; + } + if(r d(r,m,me,me,re,re); + ofs_fit_result< >()); + //use chi^2 statistic + fit.set_statistic(chisq,double,std::string>()); + fit.set_model(nfw()); + //fit.set_param_value("rs",4); + //fit.set_param_value("rho0",100); + fit.fit(); + fit.fit(); + vector p=fit.fit(); + //output parameters + ofstream ofs_param("nfw_param.txt"); + for(size_t i=0;i d=ds.get_data(i); + double x=d.get_x(); + double y=d.get_y(); + double ye=d.get_y_lower_err(); + double ym=fit.eval_model(x,p); + ofs_fit_result< +#include +#include +#include "chisq.hpp" +#include +#include +#include +#include +#include +#include + +using namespace opt_utilities; +using namespace std; +const double cm=1; +const double kpc=3.08568e+21*cm; + +int main(int argc,char* argv[]) +{ + if(argc<2) + { + cerr<<"Usage:"< [param file] [cm per pixel]"<=4) + { + cm_per_pixel=atof(argv[3]); + } + + //define the fitter + fitter,double,std::string> fit; + //define the data set + default_data_set ds; + //open the data file + ifstream ifs(argv[1]); + double min_r=1e9; + //cout<<"read serr 2"<0) + { + ofs_fit_result<<"la x radius (kpc)"<>r>>re>>t>>te; + if(!ifs.good()) + { + break; + } + min_r=min(r,min_r); + data d(r,t,te,te,re,re); + //std::cerr<0) + { + ofs_fit_result< >()); + //use chi^2 statistic + chisq,double,std::string> chisq_object; + chisq_object.set_limit(); + fit.set_statistic(chisq_object); + //fit.set_statistic(chisq,double,std::string>()); + fit.set_model(wang2012_model()); + + if(argc>=3&&std::string(argv[2])!="NONE") + { + std::vector freeze_list; + ifstream ifs_param(argv[2]); + assert(ifs_param.is_open()); + for(;;) + { + string pname; + double pvalue; + double lower,upper; + char param_status; + ifs_param>>pname>>pvalue>>lower>>upper>>param_status; + if(!ifs_param.good()) + { + break; + } + if(param_status=='F') + { + freeze_list.push_back(pname); + } + if(pvalue<=lower||pvalue>=upper) + { + cerr<<"Invalid initial value, central value not enclosed by the lower and upper boundaries, adjust automatically"<,std::string> fp(freeze_list[0]); + fit.set_param_modifier(fp); + for(size_t i=1;i,std::string>&>(fit.get_param_modifier())+=freeze_param,std::string>(freeze_list[i]); + } + } + } + + for(int i=0;i<100;++i) + { + fit.fit(); + } + vector p=fit.fit(); +#if 0 + ofstream output_param; + if(argc>=3&&std::string(argv[2])!="NONE") + { + output_param.open(argv[2]); + } + else + { + output_param.open("para0.txt"); + } +#endif + //output parameters + for(size_t i=0;i=3&&std::string(argv[2])!="NONE") +#if 0 + { + output_param<0) + { + ofs_fit_result< +#include + +namespace opt_utilities +{ + template + class nfw + :public model,std::string> + { + private: + model >* do_clone()const + { + return new nfw(*this); + } + + const char* do_get_type_name()const + { + return "1d power law"; + } + public: + nfw() + { + this->push_param_info(param_info >("rho0",1,0,1e99)); + this->push_param_info(param_info >("rs",100,0,1e99)); + } + + T do_eval(const T& r,const std::vector& param) + { + T rho0=std::abs(param[0]); + T rs=std::abs(param[1]); + static const T pi=4*std::atan(1); + return 4*pi*rho0*rs*rs*rs*(std::log((r+rs)/rs)-r/(r+rs)); + } + + private: + std::string do_get_information()const + { + return ""; + } + }; +} + + + +#endif +//EOF diff --git a/src/projector.hpp b/src/projector.hpp new file mode 100644 index 0000000..4ef4b32 --- /dev/null +++ b/src/projector.hpp @@ -0,0 +1,214 @@ +#ifndef PROJ_HPP +#define PROJ_HPP +/* + Defining the class that is used to consider the projection effect + Author: Junhua Gu + Last modified: 2011.01.01 +*/ + + +#include +#include +#include + +static const double pi=4*atan(1); +// Ratio of the electron density (n_e) to the proton density (n_p) +// n_gas = n_e + n_p ~= 1.826 n_e => n_e / n_p ~= 1.21 +// Reference: Ettori et al. 2013, Space Sci. Rev., 177, 119-154; Eq.(9) below +static const double ne_np_ratio = 1.21; + +namespace opt_utilities +{ + //This is used to project a 3-D surface brightness model to 2-D profile + template + class projector + :public model,std::vector,std::vector > + { + private: + //Points to a 3-D model that is to be projected + model,std::vector,std::vector >* pmodel; + func_obj* pcfunc; + T cm_per_pixel; + public: + //default cstr + projector() + :pmodel(NULL_PTR),pcfunc(NULL_PTR),cm_per_pixel(1) + {} + //copy cstr + projector(const projector& rhs) + :cm_per_pixel(rhs.cm_per_pixel) + { + attach_model(*(rhs.pmodel)); + if(rhs.pcfunc) + { + pcfunc=rhs.pcfunc->clone(); + } + else + { + pcfunc=NULL_PTR; + } + } + //assign operator + projector& operator=(const projector& rhs) + { + cm_per_pixel=rhs.cm_per_pixel; + if(pmodel) + { + pmodel->destroy(); + } + if(pcfunc) + { + pcfunc->destroy(); + } + if(rhs.pcfunc) + { + pcfunc=rhs.pcfunc->clone(); + } + if(rhs.pmodel) + { + pmodel=rhs.pmodel->clone(); + } + } + //destr + ~projector() + { + if(pmodel) + { + pmodel->destroy(); + } + if(pcfunc) + { + pcfunc->destroy(); + } + } + //used to clone self + model,std::vector,std::vector >* + do_clone()const + { + return new projector(*this); + } + + public: + void set_cm_per_pixel(const T& x) + { + cm_per_pixel=x; + } + + //attach the model that is to be projected + void attach_model(const model,std::vector,std::vector >& m) + { + this->clear_param_info(); + for(size_t i=0;ipush_param_info(m.get_param_info(i)); + } + this -> push_param_info(param_info, + std::string>("bkg",0,0,1E99)); + pmodel=m.clone(); + pmodel->clear_param_modifier(); + } + + void attach_cfunc(const func_obj& cf) + { + if(pcfunc) + { + pcfunc->destroy(); + } + pcfunc=cf.clone(); + } + + public: + //calc the volume + /* + This is a sphere that is subtracted by a cycline. + /| |\ + / | | \ + | | | | + | | | | + \ | | / + \| |/ + */ + T calc_v_ring(T rsph,T rcyc) + { + if(rcyc& rlist,int nsph,int nrad) + { + if(nsph& p)const + { + std::vector p1(this->reform_param(p)); + for(size_t i=0;i!=p1.size();++i) + { + if(get_element(p1,i)>this->get_param_info(i).get_upper_limit()|| + get_element(p1,i)get_param_info(i).get_lower_limit()) + { + // std::cerr<get_param_info(i).get_name()<<"\t"< p2(p1.size()-1); + for(size_t i=0;imeets_constraint(p2); + } + public: + //Perform the projection + std::vector do_eval(const std::vector& x,const std::vector& p) + { + T bkg=std::abs(p.back()); + //I think following codes are clear enough :). + std::vector unprojected(pmodel->eval(x,p)); + std::vector projected(unprojected.size()); + + for(size_t nrad=0; nrad + +using namespace std; +void report_error(const char* message) +{ + cerr<<"MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM"< +#include +#include +#include +#include + +template +class spline +{ +public: + std::vector x_list; + std::vector y_list; + std::vector y2_list; + +public: + void push_point(T x,T y) + { + if(!x_list.empty()) + { + assert(x>*(x_list.end()-1)); + } + x_list.push_back(x); + y_list.push_back(y); + } + + T get_value(T x) + { + if(x<=x_list[0]) + { + return y_list[0]; + } + if(x>=x_list.back()) + { + return y_list.back(); + } + assert(x_list.size()==y2_list.size()); + assert(x>x_list[0]); + assert(xx) + { + n2--; + } + } + T h=x_list[n2]-x_list[n1]; + double a=(x_list[n2]-x)/h; + double b=(x-x_list[n1])/h; + return a*y_list[n1]+b*y_list[n2]+((a*a*a-a)*y2_list[n1]+ + (b*b*b-b)*y2_list[n2])*(h*h)/6.; + + } + + void gen_spline(T y2_0,T y2_N) + { + int n=x_list.size(); + y2_list.resize(0); + y2_list.resize(x_list.size()); + std::vector u(x_list.size()); + if(std::abs(y2_0)::epsilon()) + { + y2_list[0]=0; + u[0]=0; + } + else + { + y2_list[0]=-.5; + u[0]=(3./(x_list[1]-x_list[0]))*((y_list[1]-y_list[0])/(x_list[1]-x_list[0])-y2_0); + } + for(int i=1;i::epsilon()) + { + qn=un=0; + } + else + { + qn=.5; + un=(3./(x_list[n-1]-x_list[n-2]))*(y2_N-(y_list[n-1]-y_list[n-2])/(x_list[n-1]-x_list[n-2])); + + } + y2_list[n-1]=(un-qn*u[n-2])/(qn*y2_list[n-2]+1.); + for(int i=n-2;i>=0;--i) + { + y2_list[i]=y2_list[i]*y2_list[i+1]+u[i]; + } + } + +}; + +#endif /* SPLINE_HPP */ diff --git a/src/vchisq.hpp b/src/vchisq.hpp new file mode 100644 index 0000000..8cd8481 --- /dev/null +++ b/src/vchisq.hpp @@ -0,0 +1,137 @@ +/** + \file vchisq.hpp + \brief chi-square statistic + \author Junhua Gu + */ + +#ifndef VCHI_SQ_HPP +#define VCHI_SQ_HPP + +#define OPT_HEADER + +#include +#include +#include +#include +#include + +using std::cerr; +using std::endl; + +namespace opt_utilities +{ + template + class vchisq + :public statistic,std::vector,std::vector,T,std::string> + { + private: + bool verb; + int n; + bool limit_bound; + typedef std::vector Tp; + + vchisq* do_clone()const + { + return new vchisq(*this); + } + + const char* do_get_type_name()const + { + return "chi^2 statistic"; + } + + public: + void verbose(bool v) + { + verb=v; + } + + void set_limit() + { + limit_bound=true; + } + + void clear_limit() + { + limit_bound=false; + } + public: + vchisq() + :verb(false),limit_bound(false) + {} + + T do_eval(const std::vector& p) + { + if(limit_bound) + { + if(!this->get_fitter().get_model().meets_constraint(p)) + { + return 1e99; + } + } + T result(0); + std::vector vx; + std::vector vy; + std::vector vye; + std::vector my; + float x1=1e99,x2=-1e99,y1=1e99,y2=-1e99; + if(verb) + { + n++; + if(n%100==0) + { + vx.resize(this->get_data_set().get_data(0).get_y().size()); + vy.resize(this->get_data_set().get_data(0).get_y().size()); + vye.resize(this->get_data_set().get_data(0).get_y().size()); + my.resize(this->get_data_set().get_data(0).get_y().size()); + } + } + for(int i=(this->get_data_set()).size()-1;i>=0;--i) + { + const std::vector y_model(this->eval_model(this->get_data_set().get_data(i).get_x(),p)); + const std::vector& y=this->get_data_set().get_data(i).get_y(); + const std::vector& ye=this->get_data_set().get_data(i).get_y_lower_err(); + for(size_t j=0;jget_data_set().get_data(i).get_x().at(j)+this->get_data_set().get_data(i).get_x().at(j+1))/2.); + vy.at(j)=(y[j]); + vye.at(j)=ye[j]; + my.at(j)=(y_model[j]); + x1=std::min(vx.at(j),x1); + y1=std::min(vy.at(j),y1); + x2=std::max(vx.at(j),x2); + y2=std::max(vy.at(j),y2); + vye[j]=log10(vy[j]+vye[j])-log10(vy[j]); + vx[j]=log10(vx[j]); + vy[j]=log10(vy[j]); + my[j]=log10(my[j]); + } + } + } + if(verb) + { + if(n%100==0) + { + cerr< +#include + +namespace opt_utilities +{ + template + class wang2012_model + :public model,std::string> + { + private: + model >* do_clone()const + { + return new wang2012_model(*this); + } + + const char* do_get_type_name()const + { + return "1d power law"; + } + public: + wang2012_model() + { + this->push_param_info(param_info >("A",5,0,500)); + this->push_param_info(param_info >("n",1.66,0,10)); + this->push_param_info(param_info >("xi",0.45,0,1)); + this->push_param_info(param_info >("a2",1500,0,1e8)); + this->push_param_info(param_info >("a3",50,0,1e8)); + this->push_param_info(param_info >("beta",0.49,0.1,0.7)); + this->push_param_info(param_info >("T0",0,0,10)); + + } + + T do_eval(const T& x,const std::vector& param) + { + T A=param[0]; + T n=param[1]; + T xi=param[2]; + T a2=param[3]; + T a3=param[4]; + T beta=param[5]; + T T0=param[6]; + return A*(pow(x,n)+xi*a2)/(pow(x,n)+a2)/pow(1+x*x/a3/a3,beta)+T0; + //return A*(pow(x,n)+a1)/(pow(x,n)+1)/pow(1+x*x/a3/a3,beta)+T0; + } + + private: + std::string do_get_information()const + { + return ""; + } + }; +} + + + +#endif +//EOF -- cgit v1.2.2