aboutsummaryrefslogtreecommitdiffstats
path: root/mass_profile
diff options
context:
space:
mode:
Diffstat (limited to 'mass_profile')
-rw-r--r--mass_profile/Makefile18
-rw-r--r--mass_profile/fit_lt_bpl.cpp336
-rw-r--r--mass_profile/fit_lt_pl.cpp239
-rw-r--r--mass_profile/fit_mt_bpl.cpp350
-rw-r--r--mass_profile/fit_mt_pl.cpp265
5 files changed, 2 insertions, 1206 deletions
diff --git a/mass_profile/Makefile b/mass_profile/Makefile
index 67e9aae..a2067ad 100644
--- a/mass_profile/Makefile
+++ b/mass_profile/Makefile
@@ -23,10 +23,8 @@ endif
OPT_UTIL_INC ?= -I../opt_utilities
-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
+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)
@@ -52,18 +50,6 @@ calc_lx_dbeta: calc_lx_dbeta.o beta_cfg.o report_error.o
calc_lx_beta: calc_lx_beta.o beta_cfg.o dump_fit_qdp.o
$(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC)
-fit_mt_pl: fit_mt_pl.cpp
- $(CXX) $(CXXFLAGS) $< -o $@ $(OPT_UTIL_INC)
-
-fit_lt_pl: fit_lt_pl.cpp
- $(CXX) $(CXXFLAGS) $< -o $@ $(OPT_UTIL_INC)
-
-fit_mt_bpl: fit_mt_bpl.cpp
- $(CXX) $(CXXFLAGS) $< -o $@ $(OPT_UTIL_INC)
-
-fit_lt_bpl: fit_lt_bpl.cpp
- $(CXX) $(CXXFLAGS) $< -o $@ $(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_lt_bpl.cpp b/mass_profile/fit_lt_bpl.cpp
deleted file mode 100644
index 0c06b22..0000000
--- a/mass_profile/fit_lt_bpl.cpp
+++ /dev/null
@@ -1,336 +0,0 @@
-/*
- 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 5 column file with T -Terr +Terr L Lerr> <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("l-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 luminosity (10\\u43\\d erg s\\u-1\\d)"<<endl;
- double yunit=1E43;
- 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 L,Ll,Lu;
- 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>>L>>Ll;
- Lu=Ll;
- //std::cerr<<L<<"\t"<<Lerr<<endl;
- if(!iss.good())
- {
- continue;
- }
-
- if(T<T_lower_limit||L<0)
- {
- continue;
- }
- if(std::abs(Lu)+std::abs(Ll)<L*.1)
- {
- double k=L*.1/(std::abs(Lu)+std::abs(Ll));
- Lu*=k;
- Ll*=k;
- }
- if(std::abs(Ll)>std::abs(L))
- {
- continue;
- }
- Tl=std::abs(Tl);
- Tu=std::abs(Tu);
- Ll=std::abs(Ll);
- Lu=std::abs(Lu);
- ofs_result<<T<<"\t"<<-std::abs(Tl)<<"\t"<<+std::abs(Tu)<<"\t"<<L/yunit<<"\t"<<-std::abs(Ll/yunit)<<"\t"<<+std::abs(Lu/yunit)<<endl;
- double x=(T);
- double y=(L);
- double xu=Tu;
- double xl=Tl;
-
- double yu=Lu;
- double yl=Ll;
- 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)/yunit<<"\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(size_t 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(size_t 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(size_t 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;
-}
diff --git a/mass_profile/fit_lt_pl.cpp b/mass_profile/fit_lt_pl.cpp
deleted file mode 100644
index fa37249..0000000
--- a/mass_profile/fit_lt_pl.cpp
+++ /dev/null
@@ -1,239 +0,0 @@
-/*
- 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 <iostream>
-#include <sstream>
-#include <fstream>
-#include <models/pl1d.hpp>
-#include <models/lin1d.hpp>
-#include "statistics/chisq.hpp"
-#include "statistics/leastsq.hpp"
-#include <data_sets/default_data_set.hpp>
-#include <methods/powell/powell_method.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 x=0;
- double u=0;
- double v=0;
-
- do
- {
- u=rand()/(double)RAND_MAX;
- rand();
- v=rand()/(double)RAND_MAX;
-
- x=std::sqrt(-log(u))*cos(2*pi*v);
- }while(isnan(x));
- return x;
-}
-
-double shuffle_data(double xc,double xl,double xu)
-{
- double result=0;
- assert(!isnan(xc));
- assert(!isnan(xl));
- assert(!isnan(xu));
- if(std_norm_rand()>0)
- {
- result=xc-std::abs(std_norm_rand()*xl);
- }
- else
- {
- result=xc+std::abs(std_norm_rand()*xu);
- }
- assert(!isnan(result));
- return result;
-}
-
-int main(int argc,char* argv[])
-{
- if(argc!=3)
- {
- cerr<<"Usage:"<<argv[0]<<" <a 5 column file with T -Terr +Terr L Lerr> <T lower limit>"<<endl;
- return -1;
- }
- double T_lower_limit(atof(argv[2]));
- ifstream ifs_data(argv[1]);
- default_data_set<double,double> ds;
- ofstream ofs_result("l-t_result.qdp");
- ofs_result<<"read terr 1"<<endl;
- ofs_result<<"read serr 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 Luminosity (10\\u43\\d erg s\\u-1\\d)"<<endl;
- double yunit=1E43;
- double sxx=0;
- double s1=0;
- double sx=0;
- double sy=0;
- double sxy=0;
- bool is_first_nonono=true;
- for(;;)
- {
- double T,Tl,Tu;
- double L,Lerr;
- std::string line;
- getline(ifs_data,line);
-
-
- 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>>L>>Lerr;
- //std::cerr<<L<<"\t"<<Lerr<<endl;
- if(!iss.good())
- {
- continue;
- }
- if(T<T_lower_limit||L<0)
- {
- continue;
- }
- if(Lerr<L*.1)
- {
- Lerr=L*.1;
- }
- double Ll=Lerr;
- double Lu=Lerr;
- Tl=std::abs(Tl);
- Tu=std::abs(Tu);
- Ll=std::abs(Ll);
- Lu=std::abs(Lu);
- ofs_result<<T<<"\t"<<-std::abs(Tl)<<"\t"<<+std::abs(Tu)<<"\t"<<L/yunit<<"\t"<<std::abs(Lerr)/yunit<<endl;
- double x=log(T);
- double y=log(L);
- double xu=log(T+Tu)-log(T);
- double xl=log(T-Tl)-log(T);
-
- double yu=log(L+Lu)-log(L);
- double yl=log(L-Ll)-log(L);
- if(isnan(x)||isnan(y)||isnan(yl)||isnan(yu)||
- isnan(xl)||isnan(xu))
- {
- std::cerr<<"one data with error > data, skipped"<<endl;
- std::cerr<<line<<endl;
- continue;
- }
- sxx+=x*x;
- sx+=x;
- sy+=y;
- sxy+=y*x;
- s1+=1;
- data<double,double> d(x,y,std::abs(yl),std::abs(yu),
- std::abs(xl),std::abs(xu));
- ds.add_data(d);
- }
-
- double M=sxx*s1-sx*sx;
- double Ma=sxy*s1-sy*sx;
- double Mb=sxx*sy-sx*sxy;
- double k0=Ma/M;
- double b0=Mb/M;
-
- 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(chisq<double,double,vector<double>,double,std::string>());
- //fit.set_statistic(leastsq<double,double,vector<double>,double,std::string>());
- fit.set_model(lin1d<double>());
- fit.load_data(ds);
-
- cerr<<"k0="<<k0<<endl;
- cerr<<"b0="<<b0<<endl;
- cerr<<"Ampl0="<<exp(b0)<<endl;
- cerr<<"gamma0="<<k0<<endl;
- fit.set_param_value("k",k0);
- fit.set_param_value("b",b0);
- fit.fit();
- std::vector<double> p=fit.fit();
- for(double i=.5;i<12;i*=1.01)
- {
- ofs_result<<i<<"\t0\t0\t"<<exp(fit.eval_model_raw(log(i),p))/yunit<<"\t0\n";
- }
-
-
- double mean_A=0;
- double mean_A2=0;
- double mean_g=0;
- double mean_g2=0;
- int cnt=0;
- for(int n=0;n<100;++n)
- {
- ++cnt;
- cerr<<".";
- opt_utilities::default_data_set<double,double> ds1;
- for(size_t i=0;i<ds.size();++i)
- {
- double new_x=shuffle_data(ds.get_data(i).get_x(),
- ds.get_data(i).get_x_lower_err(),
- ds.get_data(i).get_x_upper_err());
- double new_y=shuffle_data(ds.get_data(i).get_y(),
- ds.get_data(i).get_y_lower_err(),
- ds.get_data(i).get_y_upper_err());
- ds1.add_data(data<double,double>(new_x,new_y,
- ds.get_data(i).get_y_lower_err(),
- ds.get_data(i).get_y_upper_err(),
- ds.get_data(i).get_y_lower_err(),
- ds.get_data(i).get_y_upper_err()));
- //cerr<<new_x<<"\t"<<new_y<<endl;
- }
- fit.load_data(ds1);
-
- fit.fit();
- double k=fit.get_param_value("k");
- double b=fit.get_param_value("b");
- double A=exp(b);
- double g=k;
- mean_A+=A;
- mean_A2+=A*A;
- mean_g+=g;
- mean_g2+=g*g;
- }
- std::cerr<<endl;
- mean_A/=cnt;
- mean_A2/=cnt;
- mean_g/=cnt;
- mean_g2/=cnt;
- double std_A=std::sqrt(mean_A2-mean_A*mean_A);
- double std_g=std::sqrt(mean_g2-mean_g*mean_g);
-
- std::cerr<<"L=L0*T^gamma"<<endl;
- std::cout<<"L0= "<<exp(p[1])<<"+/-"<<std_A<<endl;
- std::cout<<"gamma= "<<p[0]<<"+/-"<<std_g<<endl;
- std::cout<<"Num of sources:"<<ds.size()<<endl;
-
-}
diff --git a/mass_profile/fit_mt_bpl.cpp b/mass_profile/fit_mt_bpl.cpp
deleted file mode 100644
index e956f4f..0000000
--- a/mass_profile/fit_mt_bpl.cpp
+++ /dev/null
@@ -1,350 +0,0 @@
-/*
- 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(size_t 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(size_t 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(size_t 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;
-}
diff --git a/mass_profile/fit_mt_pl.cpp b/mass_profile/fit_mt_pl.cpp
deleted file mode 100644
index bfccaf4..0000000
--- a/mass_profile/fit_mt_pl.cpp
+++ /dev/null
@@ -1,265 +0,0 @@
-/*
- 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/pl1d.hpp>
-#include <models/lin1d.hpp>
-#include "statistics/chisq.hpp"
-#include "statistics/leastsq.hpp"
-#include "statistics/robust_chisq.hpp"
-#include <data_sets/default_data_set.hpp>
-#include <methods/powell/powell_method.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)
-{
- if(std_norm_rand()>0)
- {
- double result=xc-std::abs(std_norm_rand()*xl);
- return result;
- }
- else
- {
- double result= xc+std::abs(std_norm_rand()*xu);
- return result;
- }
-}
-
-int main(int argc,char* argv[])
-{
- if(argc!=3)
- {
- cerr<<"Usage:"<<argv[0]<<" <a 6 column file with T -Terr +Terr M -Merr +Merr> <lower T limit>"<<endl;
- return -1;
- }
- double T_lower_limit(atof(argv[2]));
- ifstream ifs_data(argv[1]);
- default_data_set<double,double> ds;
- ofstream ofs_result("m-t_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\\dsun\\u)"<<endl;
- double sxx=0;
- double s1=0;
- double sx=0;
- double sy=0;
- double sxy=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;
- }
- 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=log(T);
- double y=log(M);
- double xu=log(T+Tu)-log(T);
- double xl=log(T-Tl)-log(T);
-
- double yu=log(M+Mu)-log(M);
- double yl=log(M-Ml)-log(M);
- if(isnan(x)||isnan(y)||isnan(yl)||isnan(yu)||
- isnan(xl)||isnan(xu))
- {
- std::cerr<<"one data with error > data, skipped"<<endl;
- std::cerr<<line<<endl;
- continue;
- }
- sxx+=x*x;
- sx+=x;
- sy+=y;
- sxy+=y*x;
- s1+=1;
- data<double,double> d(x,y,std::abs(yl),std::abs(yu),
- std::abs(xl),std::abs(xu));
- ds.add_data(d);
- }
-
- double M=sxx*s1-sx*sx;
- double Ma=sxy*s1-sy*sx;
- double Mb=sxx*sy-sx*sxy;
- double k0=Ma/M;
- double b0=Mb/M;
-
- 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(chisq<double,double,vector<double>,double,std::string>());
- //fit.set_statistic(robust_chisq<double,double,vector<double>,double,std::string>());
- //fit.set_statistic(leastsq<double,double,vector<double>,double,std::string>());
- fit.set_model(lin1d<double>());
- fit.load_data(ds);
-
- cerr<<"k0="<<k0<<endl;
- cerr<<"b0="<<b0<<endl;
- cerr<<"Ampl0="<<exp(b0)<<endl;
- cerr<<"gamma0="<<k0<<endl;
- fit.set_param_value("k",k0);
- fit.set_param_value("b",b0);
- std::vector<double> p=fit.get_all_params();
- std::cout<<"chi="<<fit.get_statistic().eval(p)<<std::endl;
- fit.fit();
- fit.fit();
- p=fit.fit();
-
- 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"<<exp(fit.eval_model_raw(log(i),p))<<"\t0\t0\n";
- }
-
- ofstream ofs_resid("resid.qdp");
- ofs_resid<<"read terr 1 2 3"<<endl;
- ofs_resid<<"skip single"<<endl;
- ofs_resid<<"ma 3 on 1"<<endl;
- ofs_resid<<"log x"<<endl;
- for(size_t i=0;i<ds.size();++i)
- {
- double x=ds.get_data(i).get_x();
- double y=ds.get_data(i).get_y();
- //double xe1=-ds.get_data(i).get_x_lower_err()*0;
- //double xe2=ds.get_data(i).get_x_upper_err()*0;
- double ye1=-ds.get_data(i).get_y_lower_err();
- double ye2=ds.get_data(i).get_y_upper_err();
- ofs_resid<<exp(x)<<"\t"<<0<<"\t"<<0<<"\t"<<y-fit.eval_model_raw(x,p)<<"\t"<<ye1<<"\t"<<ye2<<"\t"<<"0\t0\t0"<<endl;
- }
- double mean_A=0;
- double mean_A2=0;
- double mean_g=0;
- double mean_g2=0;
- int cnt=0;
- for(int n=0;n<100;++n)
- {
- ++cnt;
- cerr<<".";
- opt_utilities::default_data_set<double,double> ds1;
- for(size_t i=0;i<ds.size();++i)
- {
- double new_x=shuffle_data(ds.get_data(i).get_x(),
- ds.get_data(i).get_x_lower_err(),
- ds.get_data(i).get_x_upper_err());
- double new_y=shuffle_data(ds.get_data(i).get_y(),
- ds.get_data(i).get_y_lower_err(),
- ds.get_data(i).get_y_upper_err());
- ds1.add_data(data<double,double>(new_x,new_y,
- ds.get_data(i).get_y_lower_err(),
- ds.get_data(i).get_y_upper_err(),
- ds.get_data(i).get_y_lower_err(),
- ds.get_data(i).get_y_upper_err()));
- }
- fit.load_data(ds1);
-
- fit.fit();
- double k=fit.get_param_value("k");
- double b=fit.get_param_value("b");
- double A=exp(b);
- double g=k;
- mean_A+=A;
- mean_A2+=A*A;
- mean_g+=g;
- mean_g2+=g*g;
- }
- std::cerr<<endl;
- mean_A/=cnt;
- mean_A2/=cnt;
- mean_g/=cnt;
- mean_g2/=cnt;
- double std_A=std::sqrt(mean_A2-mean_A*mean_A);
- double std_g=std::sqrt(mean_g2-mean_g*mean_g);
-
- std::cerr<<"M=M0*T^gamma"<<endl;
- std::cout<<"M0= "<<exp(p[1])<<"+/-"<<std_A<<endl;
- std::cout<<"gamma= "<<p[0]<<"+/-"<<std_g<<endl;
- std::cout<<"Num of sources:"<<ds.size()<<endl;
-
-}