diff options
Diffstat (limited to 'mass_profile')
| -rw-r--r-- | mass_profile/Makefile | 5 | ||||
| -rw-r--r-- | mass_profile/calc_lx.cpp | 218 |
2 files changed, 1 insertions, 222 deletions
diff --git a/mass_profile/Makefile b/mass_profile/Makefile index 3bf9aa0..4f037e6 100644 --- a/mass_profile/Makefile +++ b/mass_profile/Makefile @@ -25,7 +25,7 @@ OPT_UTIL_INC ?= -I../opt_utilities TARGETS= fit_nfw_sbp fit_dbeta_sbp fit_beta_sbp \ fit_direct_beta calc_distance fit_wang2012_model \ - fit_nfw_mass calc_lx fit_mt_pl fit_lt_pl fit_mt_bpl \ + fit_nfw_mass fit_mt_pl fit_lt_pl fit_mt_bpl \ fit_lt_bpl cooling_time calc_lx_dbeta calc_lx_beta HEADERS= projector.hpp spline.hpp vchisq.hpp @@ -51,9 +51,6 @@ fit_nfw_mass: fit_nfw_mass.o fit_direct_beta: fit_direct_beta.o $(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC) -calc_lx: calc_lx.o calc_distance.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) diff --git a/mass_profile/calc_lx.cpp b/mass_profile/calc_lx.cpp deleted file mode 100644 index 21a508c..0000000 --- a/mass_profile/calc_lx.cpp +++ /dev/null @@ -1,218 +0,0 @@ -#include "calc_distance.h" -#include "spline.hpp" -#include <iostream> -#include <string> -#include <vector> -#include <statistics/chisq.hpp> -#include <methods/powell/powell_method.hpp> -#include <data_sets/default_data_set.hpp> -#include <misc/data_loaders.hpp> -#include "methods/aga/aga.hpp" -#include <models/beta1d.hpp> -#include <cstdlib> -using namespace std; -using namespace opt_utilities; - - -static const double cm=1; -static const double s=1; -static const double km=1000*100; -static const double Mpc=3.08568e+24*cm; -static const double kpc=3.08568e+21*cm; -static const double yr=365.*24.*3600.; -static const double Gyr=1e9*yr; -static const double H=71.*km/s/Mpc; -static const double c=299792458.*100.*cm; -//const double c=3e8*100*cm; -static const double pi=4*atan(1); -static const double omega_m=0.27; -static const double omega_l=0.73; -static const double arcsec2arc_ratio=1./60/60/180*pi; - -double std_norm_rand() -{ - double u=0; - double v=0; - while(u<=0||v<=0) - { - u=rand()/(double)RAND_MAX; - rand(); - v=rand()/(double)RAND_MAX; - } - double x=std::sqrt(-log(u))*cos(2*pi*v); - return x; -} - - -int main(int argc,char* argv[]) -{ - srand(time(0)); - if(argc<5) - { - cerr<<"Usage:"<<argv[0]<<" <sbp data> <ratio file> <z> <r500 in kpc> [Tprofile.dat]"<<endl; - return -1; - } - double z=atof(argv[3]); - assert(z>0); - double d_a_cm=c/H*calc_angular_distance(z); - double d_l_cm=(1+z)*(1+z)*d_a_cm; - double cm_per_pixel=d_a_cm*.492*arcsec2arc_ratio; - ////////////////////////////// - //perform a 1-beta fitting//// - ////////////////////////////// - fitter<double,double,vector<double>,double,string> f; - - f.set_statistic(chisq<double,double,vector<double>,double,string>()); - f.set_opt_method(powell_method<double,vector<double> >()); - f.set_model(beta1d<double>()); - dl_x_xe_y_ye<double,double> dl; - ifstream ifs(argv[1]); - ifs>>dl; - f.load_data(dl.get_data_set()); - f.fit(); - double rmin=f.get_data_set().get_data(0).get_x(); - double rmax=f.get_data_set().get_data(f.get_data_set().size()-1).get_x(); - ofstream lx_fit_result("lx_fit_result.qdp"); - lx_fit_result<<"read terr 1 2\nskip single\n"; - for(size_t i=0;i<f.get_data_set().size();++i) - { - lx_fit_result<<f.get_data_set().get_data(i).get_x()<<"\t"<< - -abs(f.get_data_set().get_data(i).get_x_lower_err())<<"\t"<< - abs(f.get_data_set().get_data(i).get_x_upper_err())<<"\t"<< - f.get_data_set().get_data(i).get_y()<<"\t"<< - -abs(f.get_data_set().get_data(i).get_y_lower_err())<<"\t"<< - abs(f.get_data_set().get_data(i).get_y_upper_err())<<endl; - } - lx_fit_result<<"no no no\n"; - - for(double i=rmin;i<rmax;i+=1) - { - lx_fit_result<<i<<"\t0\t0\t"<<f.eval_model(i,f.get_all_params())<<"\t0\t0"<<endl; - } - - for(size_t i=0;i<f.get_num_params();++i) - { - cerr<<f.get_param_info(i).get_name()<<"="<< - f.get_param_info(i).get_value()<<endl; - } - - vector<double> p=f.get_all_params(); - - double r500kpc=atof(argv[4]); - assert(r500kpc>0); - double r500pixel=r500kpc*kpc/cm_per_pixel; - - const double dr=1; - double sum_cnt=0; - double sum_flux=0; - spline<double> spl; - ifstream ifs_ratio(argv[2]); - assert(ifs_ratio.is_open()); - for(;;) - { - double x,y; - ifs_ratio>>x>>y; - if(!ifs_ratio.good()) - { - break; - } - spl.push_point(x,y); - } - spl.gen_spline(0,0); - - for(double r=0;r<r500pixel;r+=dr) - { - double r1=r<.2*r500pixel?.2*r500pixel:r; - double sbp=f.eval_model_raw(r1,p); - double cnt=sbp*pi*(2*r+dr)*dr; - double ratio=spl.get_value(r); - //cerr<<sbp<<endl; - sum_cnt+=cnt; - sum_flux+=cnt*ratio; - } - double lx=sum_flux*4*pi*d_l_cm*d_l_cm; - - double l_mean=0; - double l2_mean=0; - double cnt=0; - for(int n=0;n<100;++n) - { - cerr<<"."; - opt_utilities::default_data_set<double,double> ds(dl.get_data_set()); - opt_utilities::default_data_set<double,double> ds1; - for(size_t i=0;i<ds.size();++i) - { - double yc=ds.get_data(i).get_y(); - double ye=(std::abs(ds.get_data(i).get_y_lower_err())+std::abs(ds.get_data(i).get_y_lower_err()))/2; - double xc=ds.get_data(i).get_x(); - double xe=(std::abs(ds.get_data(i).get_x_lower_err())+std::abs(ds.get_data(i).get_x_lower_err()))/2; - double newy=std_norm_rand()*ye+yc; - //cout<<yc<<"\t"<<newy<<endl; - ds1.add_data(data<double,double>(xc,newy,ye,ye,xe,xe)); - } - f.load_data(ds1); - chisq<double,double,vector<double>,double,string> c; - //c.verbose(true); - f.set_statistic(c); - f.fit(); - vector<double> p=f.get_all_params(); - //cout<<f.get_param_value("beta")<<endl; - double sum_cnt=0; - double sum_flux=0; - - for(double r=0;r<r500pixel;r+=dr) - { - double r1=r<.2*r500pixel?.2*r500pixel:r; - double sbp=f.eval_model_raw(r1,p); - double cnt=sbp*pi*(2*r+dr)*dr; - double ratio=spl.get_value(r); - sum_cnt+=cnt; - sum_flux+=cnt*ratio; - } - //cout<<sum_cnt<<endl; - double lx=sum_flux*4*pi*d_l_cm*d_l_cm; - l_mean+=lx; - l2_mean+=lx*lx; - cnt+=1; - //std::cerr<<lx<<endl; - //std::cerr<<f.get_param_value( - } - l_mean/=cnt; - l2_mean/=cnt; - cerr<<endl; - double std_l=std::sqrt(l2_mean-l_mean*l_mean); - - double std_l2=0; - if(argc==6) - { - ifstream ifs_tfile(argv[5]); - assert(ifs_tfile.is_open()); - double mean_T=0; - int cnt=0; - double mean_Te=0; - for(;;) - { - double r,re,t,te; - ifs_tfile>>r>>re>>t>>te; - if(!ifs_tfile.good()) - { - break; - } - cnt+=1; - mean_T+=t; - mean_Te+=te*te; - } - mean_T/=cnt; - mean_Te/=cnt; - mean_Te=std::sqrt(mean_Te); - if(mean_Te>mean_T) - { - mean_Te=mean_T; - } - std_l2=mean_Te/mean_T*lx; - //cout<<mean_Te<<"\t"<<mean_T<<endl; - } - std_l=std::sqrt(std_l*std_l+std_l2*std_l2); - std::cout<<"Lx(bol): "<<lx<<" +/- "<<std_l<<" erg/s #Lx(bol) within r500"<<std::endl; - //std::cout<<sum_cnt<<std::endl; -} |