diff options
Diffstat (limited to 'mass_profile/calc_lx.cpp')
| -rw-r--r-- | mass_profile/calc_lx.cpp | 218 |
1 files changed, 218 insertions, 0 deletions
diff --git a/mass_profile/calc_lx.cpp b/mass_profile/calc_lx.cpp new file mode 100644 index 0000000..e5cff62 --- /dev/null +++ b/mass_profile/calc_lx.cpp @@ -0,0 +1,218 @@ +#include "calc_distance.h" +#include "spline.h" +#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(int 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(int 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(int 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; +} |