aboutsummaryrefslogtreecommitdiffstats
path: root/mass_profile/fit_mt_bpl.cpp
diff options
context:
space:
mode:
authorAaron LI <aaronly.me@gmail.com>2016-05-27 22:47:24 +0800
committerAaron LI <aaronly.me@gmail.com>2016-05-27 22:47:24 +0800
commitffd178e0bd72562a3c2cff9747b6e656edc881dc (patch)
tree8800b7b5b2e8bc3df1a6760df5cd54eaaa686702 /mass_profile/fit_mt_bpl.cpp
parent5c35fad9240fb42c1371c721e0b2af7379bd9ea0 (diff)
downloadchandra-acis-analysis-ffd178e0bd72562a3c2cff9747b6e656edc881dc.tar.bz2
Add mass_profile tools
* These tools are mainly use to calculate the total gravitational mass profile, as well as the intermediate products (e.g., surface brightness profile fitting, gas density profile, NFW fitting, etc.) * There are additional tools for calculating the luminosity and flux. * These tools mainly developed by Junhua GU, and contributed by Weitian (Aaron) LI, and Zhenghao ZHU.
Diffstat (limited to 'mass_profile/fit_mt_bpl.cpp')
-rw-r--r--mass_profile/fit_mt_bpl.cpp350
1 files changed, 350 insertions, 0 deletions
diff --git a/mass_profile/fit_mt_bpl.cpp b/mass_profile/fit_mt_bpl.cpp
new file mode 100644
index 0000000..7cf22b2
--- /dev/null
+++ b/mass_profile/fit_mt_bpl.cpp
@@ -0,0 +1,350 @@
+/*
+ Perform a double-beta density model fitting to the surface brightness data
+ Author: Junhua Gu
+ Last modified: 2011.01.01
+ This code is distributed with no warrant
+*/
+
+//#define HAVE_X_ERROR
+#include <iomanip>
+#include <iostream>
+#include <sstream>
+#include <fstream>
+#include <models/bpl1d.hpp>
+#include <models/lin1d.hpp>
+#include "statistics/chisq.hpp"
+#include "statistics/logchisq.hpp"
+#include "statistics/leastsq.hpp"
+#include <data_sets/default_data_set.hpp>
+#include <methods/powell/powell_method.hpp>
+#include <methods/gsl_simplex/gsl_simplex.hpp>
+#include <core/freeze_param.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;
+const double pi=4*atan(1);
+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;
+}
+
+double shuffle_data(double xc,double xl,double xu)
+{
+ double lxc=log(xc);
+ double lxl=log(xc-xl)-log(xc);
+ double lxu=log(xc+xu)-log(xc);
+
+ if(std_norm_rand()>0)
+ {
+ double result=std::exp(lxc-std::abs(std_norm_rand()*lxl));
+ return result;
+ }
+ else
+ {
+ double result=std::exp(lxc+std::abs(std_norm_rand()*lxu));
+ return result;
+ }
+}
+
+int main(int argc,char* argv[])
+{
+ srand(time(0));
+ if(argc!=4)
+ {
+ cerr<<"Usage:"<<argv[0]<<" <a 6 column file with T -Terr +Terr M -Merr +Merr> <initial broken temperature> <T lower limit>"<<endl;
+ return -1;
+ }
+ double T_lower_limit(atof(argv[3]));
+ double Tb=atof(argv[2]);
+ assert(Tb>0);
+ ifstream ifs_data(argv[1]);
+ default_data_set<double,double> ds;
+ ofstream ofs_result("m-t_bpl_result.qdp");
+ ofs_result<<"read terr 1 2"<<endl;
+ ofs_result<<"skip single"<<endl;
+ ofs_result<<"log"<<endl;
+ //ofs_result<<"li on 2"<<endl;
+ ofs_result<<"time off"<<endl;
+ ofs_result<<"la f"<<endl;
+ ofs_result<<"la x temperature (keV)"<<endl;
+ ofs_result<<"la y mass (M\\d\\(2281)\\u)"<<endl;
+ double sxxl=0;
+ double s1l=0;
+ double sxl=0;
+ double syl=0;
+ double sxyl=0;
+
+ double sxxu=0;
+ double s1u=0;
+ double sxu=0;
+ double syu=0;
+ double sxyu=0;
+
+ bool is_first_nonono=true;
+
+ for(;;)
+ {
+ double T,Tl,Tu;
+ double M,Ml,Mu;
+ std::string line;
+ getline(ifs_data,line);
+ //ifs_data>>T>>Tl>>Tu>>M>>Ml>>Mu;
+ if(!ifs_data.good())
+ {
+ break;
+ }
+ line+=" ";
+ istringstream iss(line);
+
+ if(line[0]=='#')
+ {
+ if(!is_first_nonono)
+ {
+ ofs_result<<"no no no"<<endl;
+ }
+ else
+ {
+ is_first_nonono=false;
+ }
+ continue;
+ }
+ iss>>T>>Tl>>Tu>>M>>Ml>>Mu;
+ //std::cerr<<L<<"\t"<<Lerr<<endl;
+ if(!iss.good())
+ {
+ continue;
+ }
+
+ if(T<T_lower_limit||M<0)
+ {
+ continue;
+ }
+
+ if(std::abs(Mu)<M*.1||std::abs(Ml)<M*.1)
+ {
+ cerr<<"mass error less than 10%, skipped"<<endl;
+ cerr<<line<<endl;
+ continue;
+ }
+#if 1
+ if(std::abs(Tu)<.1||std::abs(Tl)<.1)
+ {
+ cerr<<"T error less than 10%, skipped"<<endl;
+ cerr<<line<<endl;
+ continue;
+ }
+#endif
+
+ if(std::abs(Mu)+std::abs(Ml)<M*.1)
+ {
+ double k=M*.1/(std::abs(Mu)+std::abs(Ml));
+ Mu*=k;
+ Ml*=k;
+ }
+ if(std::abs(Ml)>std::abs(M))
+ {
+ continue;
+ }
+ Tl=std::abs(Tl);
+ Tu=std::abs(Tu);
+ Ml=std::abs(Ml);
+ Mu=std::abs(Mu);
+ ofs_result<<T<<"\t"<<-std::abs(Tl)<<"\t"<<+std::abs(Tu)<<"\t"<<M<<"\t"<<-std::abs(Ml)<<"\t"<<+std::abs(Mu)<<endl;
+ double x=(T);
+ double y=(M);
+ double xu=Tu;
+ double xl=Tl;
+
+ double yu=Mu;
+ double yl=Ml;
+ if(T>Tb)
+ {
+ sxxl+=log(x)*log(x);
+ sxl+=log(x);
+ syl+=log(y);
+ sxyl+=log(y)*log(x);
+ s1l+=1;
+ }
+ else
+ {
+ sxxu+=log(x)*log(x);
+ sxu+=log(x);
+ syu+=log(y);
+ sxyu+=log(y)*log(x);
+ s1u+=1;
+ }
+ data<double,double> d(x,y,std::abs(yl),std::abs(yu),
+ std::abs(xl),std::abs(xu));
+ ds.add_data(d);
+ }
+
+ double Ml=sxxl*s1l-sxl*sxl;
+ double Mal=sxyl*s1l-syl*sxl;
+ double Mbl=sxxl*syl-sxl*sxyl;
+ double k0l=Mal/Ml;
+ double b0l=Mbl/Ml;
+
+ double Mu=sxxu*s1u-sxu*sxu;
+ double Mau=sxyu*s1u-syu*sxu;
+ double Mbu=sxxu*syu-sxu*sxyu;
+ double k0u=Mau/Mu;
+ double b0u=Mbu/Mu;
+
+ double gamma0l=k0l;
+ double gamma0u=k0u;
+
+ double ampl0l=exp(b0l)*pow(Tb,gamma0l);
+ double ampl0u=exp(b0u)*pow(Tb,gamma0u);;
+
+
+
+ ofs_result<<"no no no"<<endl;
+ fitter<double,double,vector<double>,double,std::string> fit;
+ fit.set_opt_method(powell_method<double,vector<double> >());
+
+ fit.set_statistic(logchisq<double,double,vector<double>,double,std::string>());
+ //fit.set_statistic(leastsq<double,double,vector<double>,double,std::string>());
+ fit.set_model(bpl1d<double>());
+ fit.load_data(ds);
+
+ cerr<<"k0l="<<k0l<<endl;
+ cerr<<"k0u="<<k0u<<endl;
+ cerr<<"Ampl0="<<(ampl0l+ampl0u)/2<<endl;
+
+ fit.set_param_value("bpx",Tb);
+ fit.set_param_value("bpy",(ampl0l+ampl0u)/2);
+ fit.set_param_value("gamma1",gamma0l);
+ fit.set_param_value("gamma2",gamma0u);
+
+ fit.fit();
+ //fit.set_opt_method(gsl_simplex<double,vector<double> >());
+ fit.fit();
+ std::vector<double> p=fit.fit();
+ Tb=fit.get_param_value("bpx");
+ //std::cout<<"chi="<<fit.get_statistic().eval(p)<<std::endl;
+ for(double i=.5;i<12;i*=1.01)
+ {
+ ofs_result<<i<<"\t0\t0\t"<<fit.eval_model_raw(i,p)<<"\t0\t0\n";
+ }
+
+
+ std::vector<double> mean_p(p.size());
+ std::vector<double> mean_p2(p.size());
+ int cnt=0;
+ for(int n=0;n<100;++n)
+ {
+ ++cnt;
+ cerr<<".";
+ double sxxl=0;
+ double s1l=0;
+ double sxl=0;
+ double syl=0;
+ double sxyl=0;
+
+ double sxxu=0;
+ double s1u=0;
+ double sxu=0;
+ double syu=0;
+ double sxyu=0;
+
+ opt_utilities::default_data_set<double,double> ds1;
+ for(int i=0;i<ds.size();++i)
+ {
+ double x=ds.get_data(i).get_x();
+ double y=ds.get_data(i).get_y();
+ double xl=ds.get_data(i).get_x_lower_err();
+ double xu=ds.get_data(i).get_x_upper_err();
+ double yl=ds.get_data(i).get_y_lower_err();
+ double yu=ds.get_data(i).get_y_upper_err();
+
+ double new_x=shuffle_data(x,
+ xl,
+ xu);
+ double new_y=shuffle_data(y,
+ yl,
+ yu);
+
+ ds1.add_data(data<double,double>(new_x,new_y,
+ yl/y*new_y,
+ yu/y*new_y,
+ xl/x*new_x,
+ xu/x*new_x));
+
+ x=new_x;
+ y=new_y;
+ if(x>Tb)
+ {
+ sxxl+=log(x)*log(x);
+ sxl+=log(x);
+ syl+=log(y);
+ sxyl+=log(y)*log(x);
+ s1l+=1;
+ }
+ else
+ {
+ sxxu+=log(x)*log(x);
+ sxu+=log(x);
+ syu+=log(y);
+ sxyu+=log(y)*log(x);
+ s1u+=1;
+ }
+ }
+ double Ml=sxxl*s1l-sxl*sxl;
+ double Mal=sxyl*s1l-syl*sxl;
+ double Mbl=sxxl*syl-sxl*sxyl;
+ double k0l=Mal/Ml;
+ double b0l=Mbl/Ml;
+
+ double Mu=sxxu*s1u-sxu*sxu;
+ double Mau=sxyu*s1u-syu*sxu;
+ double Mbu=sxxu*syu-sxu*sxyu;
+ double k0u=Mau/Mu;
+ double b0u=Mbu/Mu;
+
+ double gamma0l=k0l;
+ double gamma0u=k0u;
+
+ double ampl0l=exp(b0l)*pow(Tb,gamma0l);
+ double ampl0u=exp(b0u)*pow(Tb,gamma0u);;
+
+ fit.set_param_value("bpx",Tb);
+ fit.set_param_value("bpy",(ampl0l+ampl0u)/2);
+ fit.set_param_value("gamma1",gamma0l);
+ fit.set_param_value("gamma2",gamma0u);
+
+
+ fit.load_data(ds1);
+
+ fit.fit();
+ vector<double> p=fit.fit();
+ for(int i=0;i<p.size();++i)
+ {
+ mean_p[i]+=p[i];
+ mean_p2[i]+=p[i]*p[i];
+ }
+ //cerr<<fit.get_param_value("gamma1")<<"\t"<<fit.get_param_value("gamma2")<<endl;
+
+ }
+ vector<double> std_p(p.size());
+ cerr<<endl;
+ for(int i=0;i<mean_p.size();++i)
+ {
+ mean_p[i]/=cnt;
+ mean_p2[i]/=cnt;
+ std_p[i]=std::sqrt(mean_p2[i]-mean_p[i]*mean_p[i]);
+ cout<<fit.get_param_info(i).get_name()<<"= "<<p[i]<<" +/- "<<std_p[i]<<endl;
+ }
+ std::cout<<"Num of sources:"<<ds.size()<<endl;
+}