diff options
Diffstat (limited to 'mass_profile')
| -rw-r--r-- | mass_profile/Makefile | 7 | ||||
| -rw-r--r-- | mass_profile/fit_nfw_sbp.cpp | 416 | ||||
| -rw-r--r-- | mass_profile/nfw_ne.hpp | 204 |
3 files changed, 1 insertions, 626 deletions
diff --git a/mass_profile/Makefile b/mass_profile/Makefile index c13b312..67e9aae 100644 --- a/mass_profile/Makefile +++ b/mass_profile/Makefile @@ -23,7 +23,7 @@ endif OPT_UTIL_INC ?= -I../opt_utilities -TARGETS= fit_nfw_sbp fit_dbeta_sbp fit_beta_sbp \ +TARGETS= fit_dbeta_sbp fit_beta_sbp \ fit_wang2012_model \ fit_nfw_mass fit_mt_pl fit_lt_pl fit_mt_bpl \ fit_lt_bpl calc_lx_dbeta calc_lx_beta @@ -33,8 +33,6 @@ all: $(TARGETS) # NOTE: # Object/source files should placed *before* libraries (order matters) -fit_nfw_sbp: fit_nfw_sbp.o - $(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC) fit_dbeta_sbp: fit_dbeta_sbp.o beta_cfg.o report_error.o $(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC) @@ -67,9 +65,6 @@ fit_lt_bpl: fit_lt_bpl.cpp $(CXX) $(CXXFLAGS) $< -o $@ $(OPT_UTIL_INC) -fit_nfw_sbp.o: fit_nfw_sbp.cpp nfw_ne.hpp $(HEADERS) - $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC) - fit_dbeta_sbp.o: fit_dbeta_sbp.cpp constrained_dbeta.hpp $(HEADERS) $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC) diff --git a/mass_profile/fit_nfw_sbp.cpp b/mass_profile/fit_nfw_sbp.cpp deleted file mode 100644 index 9ac8401..0000000 --- a/mass_profile/fit_nfw_sbp.cpp +++ /dev/null @@ -1,416 +0,0 @@ -/* - Fitting the surface brightness profile with an NFW-based surface brightness model - Author: Junhua Gu - Last modification 20120721 - The temperature is assumed to be an allen model with a minimum temperature assumed -*/ - - -#include <iostream> -#include <fstream> -#include "vchisq.hpp" -#include "nfw_ne.hpp" -#include <data_sets/default_data_set.hpp> -#include <methods/powell/powell_method.hpp> -#include <core/freeze_param.hpp> -#include <error_estimator/error_estimator.hpp> -#include "spline.hpp" - -using namespace std; -using namespace opt_utilities; - -//double s=5.63136645E20; -const double M_sun=1.988E33;//solar mass in g -const double kpc=3.086E21;//kpc in cm - -//A class enclosing the spline interpolation method -class spline_func_obj - :public func_obj<double,double> -{ - //has an spline object - spline<double> spl; -public: - //This function is used to calculate the intepolated value - double do_eval(const double& x) - { - /* - if(x<=spl.x_list[0]) - { - return spl.y_list[0]; - } - if(x>=spl.x_list.back()) - { - return spl.y_list.back(); - } - */ - 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); - } -}; - -//Allen temperature model - -int main(int argc,char* argv[]) -{ - if(argc!=2) - { - cerr<<argv[0]<<" <configure file>"<<endl; - cerr<<"Here is a sample of the configure file"<<endl; - cerr<<"radius_file radius1.dat\n" - "sbp_file\tsbp1.dat\n" - "cfunc_file\tcfunc.dat\n" - "n0\t\t.04\n" - "rs\t\t816\n" - "rho0\t\t.01\n" - "bkg\t\t0\n" - "cm_per_pixel\t1.804E21\n" - "z\t\t0.062476\n" - "T_file\t\tT.dat" - <<endl; - cerr<<"Notes:"<<endl; - cerr<<"n0 is in the unit of cm^-3"<<endl; - cerr<<"rs is in the unit of pixel"<<endl; - cerr<<"rho0 is in the unit of mass of proton per cm^3"<<endl; - - return -1; - } - //define a map to store the parameters - std::map<std::string,std::string> arg_map; - //open the configuration file - ifstream cfg_file(argv[1]); - assert(cfg_file.is_open()); - for(;;) - { - std::string key; - std::string value; - cfg_file>>key>>value; - if(!cfg_file.good()) - { - cfg_file.close(); - break; - } - arg_map[key]=value; - } - //check whether following parameters are defined in the configuration file - assert(arg_map.find("radius_file")!=arg_map.end()); - assert(arg_map.find("sbp_file")!=arg_map.end()); - assert(arg_map.find("cfunc_file")!=arg_map.end()); - assert(arg_map.find("T_file")!=arg_map.end()); - assert(arg_map.find("z")!=arg_map.end()); - const double z=atof(arg_map["z"].c_str()); - double r_min=0; - if(arg_map.find("r_min")!=arg_map.end()) - { - r_min=atof(arg_map["r_min"].c_str()); - cerr<<"r_min presents and its value is "<<r_min<<endl; - } - //note that in this program, the radius are not the central value of each annuli or pie region, but the boundaries. - //for example, if we have a set of radius and surface brightness values as follows: - /* - radius width surface brightness - 1 1 x - 2 1 x - 3 1 x - - then the radius is stored as - 0 1.5 2.5 3.5 - note that there should be 4 radius values, although only to represent 3 annuli. - this will be convenient to calculate the volume of each spherical shell, - and can naturally ensure the annuli are adjacent with each other, with out any gaps. - - - */ - - std::vector<double> radii;//to store radius - std::vector<double> sbps;//to store the surface brightness value - std::vector<double> sbpe;//to store the sbp error - //read in radius file - /* - About the format of the radius file: - the radius file contains only radius, separated by space or line feed (i.e., the <ENTER> key). - the unit should be pixel - - The number of radius can be larger than the number of annuli+1, the exceeded radius can be used - to calculate the influence of outer shells. - */ - int ncut=0; - for(ifstream ifs(arg_map["radius_file"].c_str());;) - { - assert(ifs.is_open()); - double x; - ifs>>x; - if(!ifs.good()) - { - break; - } - if(x<r_min) - { - ++ncut; - continue; - } - cerr<<x<<endl; - radii.push_back(x); - } - //read in surface brightness file - /* - the surface brightness file contains two columns, that are surface brightness and the error, respectively. - */ - for(ifstream ifs(arg_map["sbp_file"].c_str());;) - { - assert(ifs.is_open()); - double x,xe; - ifs>>x>>xe; - if(!ifs.good()) - { - break; - } - if(ncut) - { - --ncut; - continue; - } - cerr<<x<<"\t"<<xe<<endl; - sbps.push_back(x); - sbpe.push_back(xe); - } - //cerr<<radii.size()<<"\t"<<sbps.size()<<endl; - //return 0; - spline_func_obj cf; - - //read in the cooling function file - /* - the cooling function file contains two columns, that are radius (central value, in pixel) and the ``reduced cooling function'' - the reduced cooling function is calculated as follows - by using wabs*apec model, fill the kT, z, nH to the actual value (derived from deproject spectral fitting), - and fill norm as 1E-14/(4*pi*(Da*(1+z))^2). - then use flux e1 e2 (e1 and e2 are the energy limits of the surface brightness profile) to get the cooling function (in photon counts, - not in erg/s) - */ - for(ifstream ifs(arg_map["cfunc_file"].c_str());;) - { - assert(ifs.is_open()); - double x,y; - ifs>>x>>y; - if(!ifs.good()) - { - break; - } - cerr<<x<<"\t"<<y<<endl; - //cf.add_point(x,y*2.1249719395939022e-68);//change with source - cf.add_point(x,y);//change with source - } - cf.gen_spline(); - - for(double x=0;x<1000;x++) - { - //cout<<x<<"\t"<<cf(x)<<endl; - } - //return 0; - //cout<<radii.size()<<endl; - //cout<<sbps.size()<<endl; - - //initial a data set object and put the data together - default_data_set<std::vector<double>,std::vector<double> > ds; - ds.add_data(data<std::vector<double>,std::vector<double> >(radii,sbps,sbpe,sbpe,radii,radii)); - - //initial a fitter object - fitter<vector<double>,vector<double>,vector<double>,double> f; - //load the data set into the fitter object - f.load_data(ds); - //define a projector object - //see projector for more detailed information - projector<double> a; - //define the nfw surface brightness profofile model - nfw_ne<double> nfw; - //attach the cooling function into the projector - a.attach_cfunc(cf); - assert(arg_map.find("cm_per_pixel")!=arg_map.end()); - //set the cm to pixel ratio - double cm_per_pixel=atof(arg_map["cm_per_pixel"].c_str()); - a.set_cm_per_pixel(cm_per_pixel); - nfw.set_cm_per_pixel(cm_per_pixel); - //define the temperature profile model - spline_func_obj tf; - - for(ifstream ifs_tfunc(arg_map["T_file"].c_str());;) - { - assert(ifs_tfunc.is_open()); - double x,y; - ifs_tfunc>>x>>y; - if(!ifs_tfunc.good()) - { - break; - } - if(x<r_min) - { - continue; - } - tf.add_point(x,y); - } - tf.gen_spline(); - - //attach the temperature, surface brightness model and projector together - nfw.attach_Tfunc(tf); - a.attach_model(nfw); - f.set_model(a); - //define the chi-square statistic - vchisq<double> c; - c.verbose(true); - f.set_statistic(c); - //set the optimization method, here we use powell method - f.set_opt_method(powell_method<double,std::vector<double> >()); - //set the initial values - double n0=atof(arg_map["n0"].c_str()); - double rho0=atof(arg_map["rho0"].c_str()); - double rs=atof(arg_map["rs"].c_str()); - double bkg=atof(arg_map["bkg"].c_str()); - f.set_param_value("n0",n0); - f.set_param_value("rho0",rho0); - f.set_param_value("rs",rs); - - - f.set_param_value("bkg",bkg); - - cout<<f.get_data_set().size()<<endl; - cout<<f.get_num_params()<<endl; -#if 1 - //perform the fitting - f.fit(); - f.set_precision(1e-10); - f.fit(); - f.clear_param_modifier(); - f.fit(); -#endif - //output the parameters - ofstream param_output("nfw_param.txt"); - - for(size_t i=0;i<f.get_num_params();++i) - { - cout<<f.get_param_info(i).get_name()<<"\t"<<abs(f.get_param_info(i).get_value())<<endl; - param_output<<f.get_param_info(i).get_name()<<"\t"<<abs(f.get_param_info(i).get_value())<<endl; - } - c.verbose(false); - f.set_statistic(c); -#if 1 - f.fit(); - f.fit(); -#endif - //fetch the fitting result - std::vector<double> p=f.get_all_params(); - f.clear_param_modifier(); - std::vector<double> mv=f.eval_model(radii,p); - cerr<<mv.size()<<endl; - //output the results - ofstream ofs_sbp("sbp_fit.qdp"); - ofstream ofs_resid("resid.qdp"); - ofs_resid<<"read serr 2"<<endl; - //output the surface brightness profile - ofs_sbp<<"read serr 2"<<endl; - ofs_sbp<<"skip single"<<endl; - for(size_t i=1;i<sbps.size();++i) - { - double x=(radii[i]+radii[i-1])/2; - double y=sbps[i-1]; - double ye=sbpe[i-1]; - double ym=mv[i-1]; - ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<y<<"\t"<<ye<<"\t"<<ym<<endl; - ofs_resid<<x*cm_per_pixel/kpc<<"\t"<<(y-ym)/ye<<"\t1\n"; - } - //output the electron density - mv=nfw.eval(radii,p); - ofstream ofs_rho("rho_fit.qdp"); - for(size_t i=1;i<sbps.size();++i) - { - double x=(radii[i]+radii[i-1])/2; - double ym=mv[i-1]; - ofs_rho<<x*cm_per_pixel/kpc<<"\t"<<ym<<endl; - } - //output integral mass profile - rho0=f.get_param_value("rho0")*1.67E-24; - rs=f.get_param_value("rs")*cm_per_pixel; - ofstream ofs_int_mass("mass_int.qdp"); - for(double r=0;r<2000*kpc;r+=kpc) - { - ofs_int_mass<<r/kpc<<"\t"<<nfw_mass_enclosed(r,abs(rho0),abs(rs))/M_sun<<endl; - } - //calculate the overdensity profile - ofstream ofs_overdensity("overdensity.qdp"); - - std::vector<double> radius_list; - std::vector<double> delta_list; - - cerr<<"delta\tr_delta (kpc)\tr_delta (pixel)\tmass_delta (solar mass)\n"; - for(double r=kpc;r<6000*kpc;r+=kpc) - { - double delta=nfw_average_density(r,abs(rho0),abs(rs))/calc_critical_density(z); - radius_list.push_back(r); - delta_list.push_back(delta); - - /* - if(delta<=200&&!hit_200) - { - hit_200=true; - cerr<<200<<"\t"<<r/kpc<<"\t\t"<<r/cm_per_pixel<<"\t\t"<<nfw_mass_enclosed(r,abs(rho0),abs(rs))/M_sun<<endl; - break; - } - if(delta<=500&&!hit_500) - { - hit_500=true; - cerr<<500<<"\t"<<r/kpc<<"\t\t"<<r/cm_per_pixel<<"\t\t"<<nfw_mass_enclosed(r,abs(rho0),abs(rs))/M_sun<<endl; - } - */ - ofs_overdensity<<r/kpc<<"\t"<<delta<<endl; - } - - for(size_t i=0;i<radius_list.size()-1;++i) - { - double r=radius_list[i]; - if(delta_list[i]>=200&&delta_list[i+1]<200) - { - cerr<<200<<"\t"<<r/kpc<<"\t\t"<<r/cm_per_pixel<<"\t\t"<<nfw_mass_enclosed(r,abs(rho0),abs(rs))/M_sun<<endl; - } - - if(delta_list[i]>=500&&delta_list[i+1]<500) - { - cerr<<500<<"\t"<<r/kpc<<"\t\t"<<r/cm_per_pixel<<"\t\t"<<nfw_mass_enclosed(r,abs(rho0),abs(rs))/M_sun<<endl; - } - - if(delta_list[i]>=1500&&delta_list[i+1]<1500) - { - cerr<<1500<<"\t"<<r/kpc<<"\t\t"<<r/cm_per_pixel<<"\t\t"<<nfw_mass_enclosed(r,abs(rho0),abs(rs))/M_sun<<endl; - } - - if(delta_list[i]>=2500&&delta_list[i+1]<2500) - { - cerr<<2500<<"\t"<<r/kpc<<"\t\t"<<r/cm_per_pixel<<"\t\t"<<nfw_mass_enclosed(r,abs(rho0),abs(rs))/M_sun<<endl; - } - - } - - //output the M200 and R200 - //cerr<<"M200="<<nfw_mass_enclosed(r200,abs(rho0),abs(rs))/M_sun<<" solar mass"<<endl; - //cerr<<"R200="<<r200/kpc<<" kpc"<<endl; - //for(int i=0;i<p.size();++i) - //{ - //cerr<<p[i]<<endl; - //} - return 0; -} diff --git a/mass_profile/nfw_ne.hpp b/mass_profile/nfw_ne.hpp deleted file mode 100644 index 1eeb17c..0000000 --- a/mass_profile/nfw_ne.hpp +++ /dev/null @@ -1,204 +0,0 @@ -/* - Gas density profile derived from nfw mass profile and temperature profile - Author: Junhua Gu - Last modification: 20120721 -*/ - -#ifndef NFW_NE -#define NFW_NE -#include "projector.hpp" -#include <algorithm> -#include <functional> -#include <numeric> - -//a series of physical constants -static const double G=6.673E-8;//cm^3 g^-1 s^2 -static const double mu=1.4074; -static const double mp=1.67262158E-24;//g -static const double k=1.60217646E-9;//erg/keV -static const double c=2.99792458E10;//cm/s - - -namespace opt_utilities -{ - //the nfw mass enclosed within a radius r, with parameter rho0 and rs - template <typename T> - T nfw_mass_enclosed(T r,T rho0,T rs) - { - return 4*pi*rho0*rs*rs*rs*(std::log((r+rs)/rs)-r/(r+rs)); - } - - //average mass density - template <typename T> - T nfw_average_density(T r,T rho0,T rs) - { - if(r==0) - { - return rho0; - } - - return nfw_mass_enclosed(r,rho0,rs)/(4.*pi/3*r*r*r); - } - - //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 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 wraps method of calculating gas density from mass profile and temperature profile - template <typename T> - class nfw_ne - :public model<std::vector<T>,std::vector<T>,std::vector<T> > - { - private: - //pointer to temperature profile function - func_obj<T,T>* pTfunc; - //cm per pixel - T cm_per_pixel; - public: - //default constructor - nfw_ne() - :pTfunc(0),cm_per_pixel(1) - { - - this->push_param_info(param_info<std::vector<T>,std::string>("rho0",1));//in mp - this->push_param_info(param_info<std::vector<T>,std::string>("rs",100)); - this->push_param_info(param_info<std::vector<T>,std::string>("n0",.01)); - } - - //copy constructor - nfw_ne(const nfw_ne& rhs) - :cm_per_pixel(rhs.cm_per_pixel) - { - if(rhs.pTfunc) - { - pTfunc=rhs.pTfunc->clone(); - } - else - { - pTfunc=0; - } - //initial parameter list - this->push_param_info(param_info<std::vector<T>,std::string>("rho0",rhs.get_param_info("rho0").get_value())); - this->push_param_info(param_info<std::vector<T>,std::string>("rs",rhs.get_param_info("rs").get_value())); - this->push_param_info(param_info<std::vector<T>,std::string>("n0",rhs.get_param_info("n0").get_value())); - } - - //assignment operator - nfw_ne& operator=(const nfw_ne& rhs) - { - cm_per_pixel=rhs.cm_per_pixel; - if(pTfunc) - { - pTfunc->destroy(); - } - if(rhs.pTfunc) - { - pTfunc=rhs.pTfunc->clone(); - } - } - - //destructor - ~nfw_ne() - { - if(pTfunc) - { - pTfunc->destroy(); - } - } - - public: - //attach the temperature profile function - void attach_Tfunc(const func_obj<T,T>& Tf) - { - if(pTfunc) - { - pTfunc->destroy(); - } - pTfunc=Tf.clone(); - } - - //set the cm per pixel value - void set_cm_per_pixel(const T& x) - { - cm_per_pixel=x; - } - - //clone self - nfw_ne<T>* do_clone()const - { - return new nfw_ne<T>(*this); - } - - - //calculate density under parameters p, at radius r - /* - r is a vector, which stores a series of radius values - the annuli or pie regions are enclosed between any two - adjacent radii. - so the returned value has length smaller than r by 1. - */ - std::vector<T> do_eval(const std::vector<T> & r, - const std::vector<T>& p) - { - assert(pTfunc); - //const T kT_erg=k*5; - T rho0=std::abs(p[0])*mp; - T rs=std::abs(p[1]); - T n0=std::abs(p[2]); - T rs_cm=rs*cm_per_pixel; - - std::vector<T> yvec(r.size()); - const T kT_erg0=pTfunc->eval((r.at(0)+r.at(1))/2)*k; - //calculate the integration -#ifdef _OPENMP -#pragma omp parallel for schedule(dynamic) -#endif - for(size_t i=0;i<r.size();++i) - { - T r_cm=r[i]*cm_per_pixel; - T kT_erg=pTfunc->eval(r[i])*k; - if(abs(r_cm)==0) - { - continue; - } - yvec.at(i)=G*nfw_mass_enclosed(r_cm,rho0,rs_cm)*mu*mp/kT_erg/r_cm/r_cm; - //std::cout<<r_cm/1e20<<"\t"<<nfw_mass_enclosed(r_cm,rho0,rs_cm)/1e45<<std::endl; - //std::cout<<r_cm/1e20<<"\t"<<G*nfw_mass_enclosed(r_cm,rho0,rs_cm)*mu*mp/kT_erg/r_cm/r_cm<<std::endl; - } - - std::vector<T> ydxvec(r.size()-1); -#ifdef _OPENMP -#pragma omp parallel for schedule(dynamic) -#endif - for(size_t i=1;i<r.size();++i) - { - T dr=r[i]-r[i-1]; - T dr_cm=dr*cm_per_pixel; - ydxvec.at(i-1)=(yvec[i]+yvec[i-1])/2*dr_cm; - } - std::partial_sum(ydxvec.begin(),ydxvec.end(),ydxvec.begin()); - //construct the result - std::vector<T> result(r.size()-1); -#ifdef _OPENMP -#pragma omp parallel for schedule(dynamic) -#endif - for(size_t i=0;i<r.size()-1;++i) - { - T y=-ydxvec.at(i); - T kT_erg=pTfunc->eval(r[i])*k; - //std::cout<<y<<std::endl; - result.at(i)=n0*exp(y)*kT_erg0/kT_erg; - } - return result; - } - }; -} - -#endif |