aboutsummaryrefslogtreecommitdiffstats
path: root/mass_profile
diff options
context:
space:
mode:
Diffstat (limited to 'mass_profile')
-rw-r--r--mass_profile/Makefile86
-rw-r--r--mass_profile/README.md37
-rw-r--r--mass_profile/beta.hpp45
-rw-r--r--mass_profile/beta_cfg.cpp85
-rw-r--r--mass_profile/beta_cfg.hpp23
-rw-r--r--mass_profile/calc_lx_beta.cpp497
-rw-r--r--mass_profile/calc_lx_dbeta.cpp548
-rw-r--r--mass_profile/chisq.hpp206
-rw-r--r--mass_profile/dbeta.hpp108
-rw-r--r--mass_profile/dump_fit_qdp.cpp55
-rw-r--r--mass_profile/dump_fit_qdp.hpp16
-rw-r--r--mass_profile/fit_beta_sbp.cpp534
-rw-r--r--mass_profile/fit_dbeta_sbp.cpp586
-rw-r--r--mass_profile/fit_nfw_mass.cpp159
-rw-r--r--mass_profile/fit_wang2012_model.cpp195
-rw-r--r--mass_profile/nfw.hpp56
-rw-r--r--mass_profile/projector.hpp214
-rw-r--r--mass_profile/report_error.cpp39
-rw-r--r--mass_profile/report_error.hpp6
-rw-r--r--mass_profile/spline.hpp110
-rw-r--r--mass_profile/vchisq.hpp137
-rw-r--r--mass_profile/wang2012_model.hpp67
22 files changed, 0 insertions, 3809 deletions
diff --git a/mass_profile/Makefile b/mass_profile/Makefile
deleted file mode 100644
index 17b7a0e..0000000
--- a/mass_profile/Makefile
+++ /dev/null
@@ -1,86 +0,0 @@
-#
-# Makefile for `chandra-acis-analysis/mass_profile` tools
-#
-# Junhua GU
-# Weitian LI
-# 2016-06-07
-#
-
-
-CXX ?= g++
-CXXFLAGS += -Wall
-CXXFLAGS += -Werror
-
-ifdef OPENMP
- CXXFLAGS += -fopenmp
-endif
-
-ifdef DEBUG
- CXXFLAGS += -g
-else
- CXXFLAGS += -O2
-endif
-
-OPT_UTIL_INC ?= -I../opt_utilities
-
-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)
-
-# NOTE:
-# Object/source files should placed *before* libraries (order matters)
-
-fit_dbeta_sbp: fit_dbeta_sbp.o beta_cfg.o report_error.o
- $(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC)
-
-fit_beta_sbp: fit_beta_sbp.o beta_cfg.o dump_fit_qdp.o
- $(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC)
-
-fit_wang2012_model: fit_wang2012_model.o
- $(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC)
-
-fit_nfw_mass: fit_nfw_mass.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)
-
-calc_lx_beta: calc_lx_beta.o beta_cfg.o dump_fit_qdp.o
- $(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC)
-
-
-fit_dbeta_sbp.o: fit_dbeta_sbp.cpp $(HEADERS)
- $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC)
-
-fit_beta_sbp.o: fit_beta_sbp.cpp beta.hpp $(HEADERS)
- $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC)
-
-fit_wang2012_model.o: fit_wang2012_model.cpp wang2012_model.hpp chisq.hpp
- $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC)
-
-fit_nfw_mass.o: fit_nfw_mass.cpp nfw.hpp chisq.hpp
- $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC)
-
-calc_lx_dbeta.o: calc_lx_dbeta.cpp $(HEADERS)
- $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC)
-
-calc_lx_beta.o: calc_lx_beta.cpp beta.hpp $(HEADERS)
- $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC)
-
-beta_cfg.o: beta_cfg.cpp beta_cfg.hpp
- $(CXX) $(CXXFLAGS) -c $<
-
-report_error.o: report_error.cpp report_error.hpp
- $(CXX) $(CXXFLAGS) -c $<
-
-dump_fit_qdp.o: dump_fit_qdp.cpp dump_fit_qdp.hpp
- $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC)
-
-%.o: %.cpp
- $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC)
-
-
-clean:
- rm -f *.o $(TARGETS)
diff --git a/mass_profile/README.md b/mass_profile/README.md
deleted file mode 100644
index 0c61821..0000000
--- a/mass_profile/README.md
+++ /dev/null
@@ -1,37 +0,0 @@
-Mass/Luminosity/Flux Calculation Tools
-==============================
-
-Junhua GU, Weitian LI, Zhenghao ZHU
-
-
-Introduction
-------------
-This directory contains the tools/scripts used to fit the temperature
-profile and surface brightness file, to calculate the gas density profile
-and gravitational mass profile, to calculate luminosity and flux and
-other related quantities.
-
-
-NOTE
-----
-* Mass calculation references:
- + Walker et al. 2012, MNRAS, 422, 3503-3515,
- [ADS:2012MNRAS.422.3503W](http://adsabs.harvard.edu/abs/2012MNRAS.422.3503W)
- + Ettori et al. 2013, Space Science Reviews, 177, 119-154,
- [ADS:2013SSRv..177..119E](http://adsabs.harvard.edu/abs/2013SSRv..177..119E)
-* Errors are calculated at 68% confident level.
-* NFW profile is employed to extrapolate the mass profile.
-* We use a self-proposed model (see ``wang2012_model.hpp``) to describe/fit
- the gas temperature profile.
-* The single-beta and double-beta models (with a constant background) are used
- to describe the gas density profile.
-
-
-TODO
-----
-* Merge the scripts/tools for single-beta and double-beta SBP models
- (to reduce code duplications)
-* Uncertainties/errors calculation **maybe** inappropriate/problematic
- (e.g., ``analyze_mass_profile.py``);
- why not just use the quantile-based (e.g., Q84.15-Q15.85 ~= 68.3%)
- uncertainties or standard deviation???
diff --git a/mass_profile/beta.hpp b/mass_profile/beta.hpp
deleted file mode 100644
index 6ae70ac..0000000
--- a/mass_profile/beta.hpp
+++ /dev/null
@@ -1,45 +0,0 @@
-#ifndef BETA
-#define BETA
-#include "projector.hpp"
-
-namespace opt_utilities
-{
- template <typename T>
- class beta
- :public model<std::vector<T>,std::vector<T>,std::vector<T> >
- {
- public:
- beta()
- {
- this->push_param_info(param_info<std::vector<T>,std::string>("n0",1,0,1E99));
- this->push_param_info(param_info<std::vector<T>,std::string>("beta",.66,0,1E99));
- this->push_param_info(param_info<std::vector<T>,std::string>("rc",100,0,1E99));
- }
-
- public:
- beta<T>* do_clone()const
- {
- return new beta<T>(*this);
- }
-
- std::vector<T> do_eval(const std::vector<T> & x,
- const std::vector<T>& p)
- {
- T n0=std::abs(p[0]);
- T beta=p[1];
- T rc=p[2];
-
- std::vector<T> result(x.size()-1);
- for(size_t i=1;i<x.size();++i)
- {
- T xi=(x[i]+x[i-1])/2;
- T yi=0;
- yi=n0*pow(1+xi*xi/rc/rc,-3./2.*beta);
- result[i-1]=yi;
- }
- return result;
- }
- };
-}
-
-#endif
diff --git a/mass_profile/beta_cfg.cpp b/mass_profile/beta_cfg.cpp
deleted file mode 100644
index c373688..0000000
--- a/mass_profile/beta_cfg.cpp
+++ /dev/null
@@ -1,85 +0,0 @@
-#include "beta_cfg.hpp"
-#include <sstream>
-#include <cstdlib>
-using namespace std;
-
-cfg_map parse_cfg_file(std::istream& is)
-{
- cfg_map result;
- result.rmin_pixel=-1;
- result.rmin_kpc=-1;
- for(;;)
- {
- std::string line;
- getline(is,line);
- line+="\n";
- if(!is.good())
- {
- break;
- }
- string key;
- istringstream iss(line);
- iss>>key;
- if(key=="sbp_data")
- {
- string value;
- iss>>value;
- result.sbp_data=value;
- }
- else if(key=="cfunc_profile")
- {
- string value;
- iss>>value;
- result.cfunc_profile=value;
- }
- else if(key=="tprofile")
- {
- string value;
- iss>>value;
- result.tprofile=value;
- }
- else if(key=="z")
- {
- double z;
- iss>>z;
- result.z=z;
- }
- else if(key=="cm_per_pixel")
- {
- double cm_per_pixel;
- iss>>cm_per_pixel;
- result.cm_per_pixel=cm_per_pixel;
- }
- else if(key=="rmin_pixel")
- {
- double v;
- iss>>v;
- result.rmin_pixel=v;
- }
- else if(key=="rmin_kpc")
- {
- double v;
- iss>>v;
- result.rmin_kpc=v;
- }
- else
- {
- std::vector<double> value;
- for(;;)
- {
- double v;
- iss>>v;
- if(!iss.good())
- {
- break;
- }
- value.push_back(v);
- }
- if(!value.empty())
- {
- result.param_map[key]=value;
- }
- }
- }
- return result;
-}
diff --git a/mass_profile/beta_cfg.hpp b/mass_profile/beta_cfg.hpp
deleted file mode 100644
index 27148df..0000000
--- a/mass_profile/beta_cfg.hpp
+++ /dev/null
@@ -1,23 +0,0 @@
-#ifndef BETA_CFG
-#define BETA_CFG
-
-#include <map>
-#include <vector>
-#include <string>
-#include <iostream>
-
-struct cfg_map
-{
- std::string sbp_data;
- std::string cfunc_profile;
- std::string tprofile;
- double z;
- double cm_per_pixel;
- double rmin_kpc;
- double rmin_pixel;
- std::map<std::string,std::vector<double> > param_map;
-};
-
-cfg_map parse_cfg_file(std::istream& is);
-
-#endif
diff --git a/mass_profile/calc_lx_beta.cpp b/mass_profile/calc_lx_beta.cpp
deleted file mode 100644
index 11d7b04..0000000
--- a/mass_profile/calc_lx_beta.cpp
+++ /dev/null
@@ -1,497 +0,0 @@
-/**
- * Calculate the total luminosity and flux within the specified radius.
- *
- * Base on 'fit_beta_sbp.cpp' and supersede 'calc_lx.cpp'
- *
- * Author: Junhua Gu
- */
-
-#include <iostream>
-#include <fstream>
-#include <sstream>
-#include <list>
-#include <algorithm>
-#include "beta_cfg.hpp"
-#include "dump_fit_qdp.hpp"
-#include "report_error.hpp"
-#include "vchisq.hpp"
-#include "chisq.hpp"
-#include "beta.hpp"
-#include "models/beta1d.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 kpc=3.086E21;//kpc in cm
-const double Mpc=kpc*1000;
-
-double beta_func(double r, double n0, double rc, double beta)
-{
- return abs(n0) * pow(1+r*r/rc/rc, -3./2.*abs(beta));
-}
-
-//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)
- {
- 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);
- }
-};
-
-int main(int argc,char* argv[])
-{
- if(argc<4)
- {
- cerr<<argv[0]<<" <sbp.conf> <rout_kpc> <cfunc_erg> [cfunc2_erg ...]"<<endl;
- return -1;
- }
- //initialize the parameters list
- ifstream cfg_file(argv[1]);
- assert(cfg_file.is_open());
- cfg_map cfg=parse_cfg_file(cfg_file);
-
- const double z=cfg.z;
-
- //initialize the radius list, sbp list and sbp error list
- std::vector<double> radii;
- std::vector<double> sbps;
- std::vector<double> sbpe;
- std::vector<double> radii_all;
- std::vector<double> sbps_all;
- std::vector<double> sbpe_all;
- //prepend the zero point to radius list
- radii.push_back(0.0);
- radii_all.push_back(0.0);
- //read sbp and sbp error data
- ifstream ifs(cfg.sbp_data.c_str());
- std::string line;
- if (ifs.is_open())
- {
- while(std::getline(ifs, line))
- {
- if (line.empty())
- continue;
-
- std::istringstream iss(line);
- double x, xe, y, ye;
- if ((iss >> x >> xe >> y >> ye))
- {
- std::cerr << "sbprofile data: "
- << x << ", " << xe << ", " << y << ", " << ye
- << std::endl;
- radii.push_back(x+xe); /* NOTE: use outer radii of regions */
- radii_all.push_back(x+xe);
- sbps.push_back(y);
- sbps_all.push_back(y);
- sbpe.push_back(ye);
- sbpe_all.push_back(ye);
- }
- else
- {
- std::cerr << "skipped line: " << line << std::endl;
- }
- }
- }
- else
- {
- std::cerr << "ERROR: cannot open file: " << cfg.sbp_data.c_str()
- << std::endl;
- return 1;
- }
-
- //initialize the cm/pixel value
- double cm_per_pixel=cfg.cm_per_pixel;
- double rmin;
- if(cfg.rmin_pixel>0)
- {
- rmin=cfg.rmin_pixel;
- }
- else
- {
- rmin=cfg.rmin_kpc*kpc/cm_per_pixel;
- }
-
- cerr<<"rmin="<<rmin<<" (pixel)"<<endl;
- std::list<double> radii_tmp,sbps_tmp,sbpe_tmp;
- radii_tmp.resize(radii.size());
- sbps_tmp.resize(sbps.size());
- sbpe_tmp.resize(sbpe.size());
- copy(radii.begin(),radii.end(),radii_tmp.begin());
- copy(sbps.begin(),sbps.end(),sbps_tmp.begin());
- copy(sbpe.begin(),sbpe.end(),sbpe_tmp.begin());
- for(list<double>::iterator i=radii_tmp.begin();i!=radii_tmp.end();)
- {
- if(*i<rmin)
- {
- radii_tmp.pop_front();
- sbps_tmp.pop_front();
- sbpe_tmp.pop_front();
- i=radii_tmp.begin();
- continue;
- }
- ++i;
- }
- radii.resize(radii_tmp.size());
- sbps.resize(sbps_tmp.size());
- sbpe.resize(sbpe_tmp.size());
- copy(radii_tmp.begin(),radii_tmp.end(),radii.begin());
- copy(sbps_tmp.begin(),sbps_tmp.end(),sbps.begin());
- copy(sbpe_tmp.begin(),sbpe_tmp.end(),sbpe.begin());
-
- //read cooling function data
- spline_func_obj cf;
- for(ifstream ifs(cfg.cfunc_profile.c_str());;)
- {
- assert(ifs.is_open());
- double x,y;
- ifs>>x>>y;
- if(!ifs.good())
- {
- break;
- }
- cerr<<x<<"\t"<<y<<endl;
- if(x>radii.back())
- {
- break;
- }
- cf.add_point(x,y);
- }
- cf.gen_spline();
-
- //read temperature profile data
- spline_func_obj Tprof;
- int tcnt=0;
- for(ifstream ifs1(cfg.tprofile.c_str());;++tcnt)
- {
- assert(ifs1.is_open());
- double x,y;
- ifs1>>x>>y;
- if(!ifs1.good())
- {
- break;
- }
- cerr<<x<<"\t"<<y<<endl;
-#if 0
- if(tcnt==0)
- {
- Tprof.add_point(0,y);
- }
-#endif
- Tprof.add_point(x,y);
- }
-
-
- Tprof.gen_spline();
-
- 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 fitter
- fitter<vector<double>,vector<double>,vector<double>,double> f;
- f.load_data(ds);
- //initial the object, which is used to calculate projection effect
- projector<double> a;
- beta<double> betao;
- //attach the cooling function
- a.attach_cfunc(cf);
- a.set_cm_per_pixel(cm_per_pixel);
- a.attach_model(betao);
- f.set_model(a);
- //chi^2 statistic
- vchisq<double> c;
- c.verbose(true);
- c.set_limit();
- f.set_statistic(c);
- //optimization method
- f.set_opt_method(powell_method<double,std::vector<double> >());
- //initialize the initial values
- /*
- double n0=0;
- double beta=0;
- double rc=0;
- */
- double bkg_level=0;
-
- for(std::map<std::string,std::vector<double> >::iterator i=cfg.param_map.begin();
- i!=cfg.param_map.end();++i)
- {
- std::string pname=i->first;
- f.set_param_value(pname,i->second.at(0));
- if(i->second.size()==3)
- {
- double a1=i->second[1];
- double a2=i->second[2];
- double u=std::max(a1,a2);
- double l=std::min(a1,a2);
- f.set_param_upper_limit(pname,u);
- f.set_param_lower_limit(pname,l);
- }
- else
- {
- if(pname=="beta")
- {
- f.set_param_lower_limit(pname,.3);
- f.set_param_upper_limit(pname,1.4);
- }
- }
- }
-
- f.fit();
- f.fit();
- std::vector<double> p=f.get_all_params();
- /*
- n0=f.get_param_value("n0");
- rc=f.get_param_value("rc");
- beta=f.get_param_value("beta");
- */
- //output the datasets and fitting results
- ofstream param_output("lx_beta_param.txt");
- for(size_t i=0;i<f.get_num_params();++i)
- {
- if(f.get_param_info(i).get_name()=="rc")
- {
- cerr<<"rc_kpc"<<"\t"<<abs(f.get_param_info(i).get_value())*cm_per_pixel/kpc<<endl;
- param_output<<"rc_kpc"<<"\t"<<abs(f.get_param_info(i).get_value())*cm_per_pixel/kpc<<endl;
- }
- cerr<<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;
- }
- cerr<<"reduced_chi^2="<<f.get_statistic_value()/(radii.size()-f.get_model().get_num_free_params())<<endl;
- param_output<<"reduced_chi^2="<<f.get_statistic_value()/(radii.size()-f.get_model().get_num_free_params())<<endl;
-
- std::vector<double> mv=f.eval_model_raw(radii_all,p);
- int sbps_inner_cut_size=int(sbps_all.size()-sbps.size());
- ofstream ofs_sbp("lx_sbp_fit.qdp");
- ofs_sbp<<"read serr 2"<<endl;
- ofs_sbp<<"skip single"<<endl;
- ofs_sbp<<"line off "<<endl;
- if(sbps_inner_cut_size>=1)
- {
- ofs_sbp<<"line on 2"<<endl;
- ofs_sbp<<"line on 3"<<endl;
- ofs_sbp<<"line on 4"<<endl;
- ofs_sbp<<"line on 5"<<endl;
- ofs_sbp<<"ls 2 on 2"<<endl;
- ofs_sbp<<"ls 2 on 4"<<endl;
- ofs_sbp<<"ls 2 on 5"<<endl;
- ofs_sbp<<"line on 7"<<endl;
- ofs_sbp<<"ls 2 on 7"<<endl;
-
- ofs_sbp<<"ma 1 on 2"<<endl;
- ofs_sbp<<"color 1 on 1"<<endl;
- ofs_sbp<<"color 2 on 2"<<endl;
- ofs_sbp<<"color 3 on 3"<<endl;
- ofs_sbp<<"color 4 on 4"<<endl;
- ofs_sbp<<"color 5 on 5"<<endl;
-
- ofs_sbp<<"win 1"<<endl;
- ofs_sbp<<"yplot 1 2 3 4 5"<<endl;
- ofs_sbp<<"loc 0 0 1 1"<<endl;
- ofs_sbp<<"vie .1 .4 .9 .9"<<endl;
- ofs_sbp<<"la y cnt/s/pixel/cm^2"<<endl;
- ofs_sbp<<"log x"<<endl;
- ofs_sbp<<"log y"<<endl;
- ofs_sbp<<"r x "<<(radii[1]+radii[0])/2*cm_per_pixel/kpc<<" "<<(radii[sbps.size()-2]+radii[sbps.size()-1])/2*cm_per_pixel/kpc<<endl;
- ofs_sbp<<"win 2"<<endl;
- ofs_sbp<<"yplot 6 7"<<endl;
- ofs_sbp<<"loc 0 0 1 1"<<endl;
- ofs_sbp<<"vie .1 .1 .9 .4"<<endl;
- ofs_sbp<<"la x radius (kpc)"<<endl;
- ofs_sbp<<"la y chi"<<endl;
- ofs_sbp<<"log y off"<<endl;
- ofs_sbp<<"log x"<<endl;
- ofs_sbp<<"r x "<<(radii[1]+radii[0])/2*cm_per_pixel/kpc<<" "<<(radii[sbps.size()-2]+radii[sbps.size()-1])/2*cm_per_pixel/kpc<<endl;
- }
- else
- {
- ofs_sbp<<"line on 2"<<endl;
- ofs_sbp<<"line on 3"<<endl;
- ofs_sbp<<"line on 4"<<endl;
- ofs_sbp<<"ls 2 on 3"<<endl;
- ofs_sbp<<"ls 2 on 4"<<endl;
- ofs_sbp<<"line on 6"<<endl;
- ofs_sbp<<"ls 2 on 6"<<endl;
-
- ofs_sbp<<"color 1 on 1"<<endl;
- ofs_sbp<<"color 3 on 2"<<endl;
- ofs_sbp<<"color 4 on 3"<<endl;
- ofs_sbp<<"color 5 on 4"<<endl;
- //ofs_sbp<<"ma 1 on 2"<<endl;
-
- ofs_sbp<<"win 1"<<endl;
- ofs_sbp<<"yplot 1 2 3 4"<<endl;
- ofs_sbp<<"loc 0 0 1 1"<<endl;
- ofs_sbp<<"vie .1 .4 .9 .9"<<endl;
- ofs_sbp<<"la y cnt/s/pixel/cm^2"<<endl;
- ofs_sbp<<"log x"<<endl;
- ofs_sbp<<"log y"<<endl;
- ofs_sbp<<"r x "<<(radii[1]+radii[0])/2*cm_per_pixel/kpc<<" "<<(radii[radii.size()-2]+radii[radii.size()-1])/2*cm_per_pixel/kpc<<endl;
- ofs_sbp<<"win 2"<<endl;
- ofs_sbp<<"yplot 5 6"<<endl;
- ofs_sbp<<"loc 0 0 1 1"<<endl;
- ofs_sbp<<"vie .1 .1 .9 .4"<<endl;
- ofs_sbp<<"la x radius (kpc)"<<endl;
- ofs_sbp<<"la y chi"<<endl;
- ofs_sbp<<"log x"<<endl;
- ofs_sbp<<"log y off"<<endl;
- ofs_sbp<<"r x "<<(radii[1]+radii[0])/2*cm_per_pixel/kpc<<" "<<(radii[radii.size()-2]+radii[radii.size()-1])/2*cm_per_pixel/kpc<<endl;
-
- }
- // cout<<sbps_all.size()<<"\t"<<sbps.size()<<"\t"<<sbps_inner_cut_size<<endl;
- for(size_t i=1;i<sbps_all.size();++i)
- {
- double x=(radii_all[i]+radii_all[i-1])/2;
- double y=sbps_all[i-1];
- double ye=sbpe_all[i-1];
- ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<y<<"\t"<<ye<<endl;
- }
- if(sbps_inner_cut_size>=1)
- {
- ofs_sbp<<"no no no"<<endl;
- for(int i=1;i<sbps_inner_cut_size+1;++i)
- {
- double x=(radii_all[i]+radii_all[i-1])/2;
- double ym=mv[i-1];
- ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<ym<<"\t"<<"0"<<endl;
- }
- }
- ofs_sbp<<"no no no"<<endl;
- for(size_t i=sbps_inner_cut_size;i<sbps_all.size();++i)
- {
- double x=(radii_all[i]+radii_all[i-1])/2;
- double ym=mv[i-1];
- ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<ym<<"\t"<<"0"<<endl;
- }
- ofs_sbp<<"no no no"<<endl;
- //bkg level
- bkg_level=abs(f.get_param_value("bkg"));
- for(size_t i=0;i<sbps_all.size();++i)
- {
- double x=(radii_all[i]+radii_all[i-1])/2;
- ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<bkg_level<<"\t0"<<endl;
- }
- ofs_sbp<<"no no no"<<endl;
- //rc
- double rc_kpc=abs(f.get_param_value("rc")*cm_per_pixel/kpc);
- double max_sbp=*max_element(sbps_all.begin(),sbps_all.end());
- double min_sbp=*min_element(sbps_all.begin(),sbps_all.end());
- for(double x=min_sbp;x<=max_sbp;x+=(max_sbp-min_sbp)/100)
- {
- ofs_sbp<<rc_kpc<<"\t"<<x<<"\t"<<"0"<<endl;
- }
- //resid
- ofs_sbp<<"no no no"<<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"<<(ym-sbps[i-1])/sbpe[i-1]<<"\t"<<1<<endl;
- }
-
- //zero level of resid
- ofs_sbp<<"no no no"<<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"<<0<<"\t"<<0<<endl;
- }
-
- mv=betao.eval(radii,p);
- ofstream ofs_rho("lx_rho_fit.qdp");
- ofstream ofs_rho_data("lx_rho_fit.dat");
-
- ofs_rho<<"la x radius (kpc)"<<endl;
- ofs_rho<<"la y density (cm\\u-3\\d)"<<endl;
- /*
- for(int 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;
- }
- */
- p.back()=0;
- radii.clear();
- double rout=atof(argv[2])*kpc;
-
- for(double r=0;r<rout;r+=1*kpc)//step size=1kpc
- {
- double r_pix=r/cm_per_pixel;
- radii.push_back(r_pix);
- }
-
- double Da=cm_per_pixel/(.492/3600./180.*pi);
- double Dl=Da*(1+z)*(1+z);
- cout<<"dl="<<Dl/kpc<<endl;
-
- for(int n=3;n<argc;++n)
- {
- spline_func_obj cf_erg;
- for(ifstream ifs(argv[n]);;)
- {
- assert(ifs.is_open());
- double x,y;
- ifs>>x>>y;
- if(!ifs.good())
- {
- break;
- }
- //cerr<<x<<"\t"<<y<<endl;
-
- cf_erg.add_point(x,y);//change with source
- }
- cf_erg.gen_spline();
-
- projector<double>& pj=dynamic_cast<projector<double>&>(f.get_model());
- pj.attach_cfunc(cf_erg);
-
-
-
- mv=f.eval_model_raw(radii,p);
- double flux_erg=0;
- for(size_t i=0;i<radii.size()-1;++i)
- {
- double S=pi*(radii[i+1]+radii[i])*(radii[i+1]-radii[i]);
- flux_erg+=S*mv[i];
- }
- cout<<flux_erg*4*pi*Dl*Dl<<endl;
- cout<<flux_erg<<endl;
- param_output<<"Lx"<<n-2<<"\t"<<flux_erg*4*pi*Dl*Dl<<endl;
- param_output<<"Fx"<<n-2<<"\t"<<flux_erg<<endl;
- }
-}
diff --git a/mass_profile/calc_lx_dbeta.cpp b/mass_profile/calc_lx_dbeta.cpp
deleted file mode 100644
index 31addeb..0000000
--- a/mass_profile/calc_lx_dbeta.cpp
+++ /dev/null
@@ -1,548 +0,0 @@
-/**
- * Calculate the total luminosity and flux within the specified radius.
- *
- * Base on 'fit_dbeta_sbp.cpp' and supersede 'calc_lx.cpp'
- *
- * Author: Junhua Gu
- */
-
-#include <iostream>
-#include <fstream>
-#include <sstream>
-#include <list>
-#include "vchisq.hpp"
-#include "dbeta.hpp"
-#include "beta_cfg.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 kpc=3.086E21;//kpc in cm
-const double Mpc=kpc*1000;
-
-double dbeta_func(double r, double n01, double rc1, double beta1,
- double n02, double rc2, double beta2)
-{
- double v1 = abs(n01) * pow(1+r*r/rc1/rc1, -3./2.*abs(beta1));
- double v2 = abs(n02) * pow(1+r*r/rc2/rc2, -3./2.*abs(beta2));
- return v1 + v2;
-}
-
-
-//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)
- {
- 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);
- }
-};
-
-int main(int argc,char* argv[])
-{
- if(argc<4)
- {
- cerr<<argv[0]<<" <sbp.conf> <rout_kpc> <cfunc_erg> [cfunc2_erg ...]"<<endl;
- return -1;
- }
- //initialize the parameters list
- ifstream cfg_file(argv[1]);
- assert(cfg_file.is_open());
- cfg_map cfg=parse_cfg_file(cfg_file);
-
- const double z=cfg.z;
-
- //initialize the radius list, sbp list and sbp error list
- std::vector<double> radii;
- std::vector<double> sbps;
- std::vector<double> sbpe;
- std::vector<double> radii_all;
- std::vector<double> sbps_all;
- std::vector<double> sbpe_all;
- //prepend the zero point to radius list
- radii.push_back(0.0);
- radii_all.push_back(0.0);
- //read sbp and sbp error data
- ifstream ifs(cfg.sbp_data.c_str());
- std::string line;
- if (ifs.is_open())
- {
- while(std::getline(ifs, line))
- {
- if (line.empty())
- continue;
-
- std::istringstream iss(line);
- double x, xe, y, ye;
- if ((iss >> x >> xe >> y >> ye))
- {
- std::cerr << "sbprofile data: "
- << x << ", " << xe << ", " << y << ", " << ye
- << std::endl;
- radii.push_back(x+xe); /* NOTE: use outer radii of regions */
- radii_all.push_back(x+xe);
- sbps.push_back(y);
- sbps_all.push_back(y);
- sbpe.push_back(ye);
- sbpe_all.push_back(ye);
- }
- else
- {
- std::cerr << "skipped line: " << line << std::endl;
- }
- }
- }
- else
- {
- std::cerr << "ERROR: cannot open file: " << cfg.sbp_data.c_str()
- << std::endl;
- return 1;
- }
-
- //initialize the cm/pixel value
- double cm_per_pixel=cfg.cm_per_pixel;
- double rmin;
- if(cfg.rmin_pixel>0)
- {
- rmin=cfg.rmin_pixel;
- }
- else
- {
- rmin=cfg.rmin_kpc*kpc/cm_per_pixel;
- }
-
- cerr<<"rmin="<<rmin<<" (pixel)"<<endl;
- std::list<double> radii_tmp,sbps_tmp,sbpe_tmp;
- radii_tmp.resize(radii.size());
- sbps_tmp.resize(sbps.size());
- sbpe_tmp.resize(sbpe.size());
- copy(radii.begin(),radii.end(),radii_tmp.begin());
- copy(sbps.begin(),sbps.end(),sbps_tmp.begin());
- copy(sbpe.begin(),sbpe.end(),sbpe_tmp.begin());
- for(list<double>::iterator i=radii_tmp.begin();i!=radii_tmp.end();)
- {
- if(*i<rmin)
- {
- radii_tmp.pop_front();
- sbps_tmp.pop_front();
- sbpe_tmp.pop_front();
- i=radii_tmp.begin();
- continue;
- }
- ++i;
- }
- radii.resize(radii_tmp.size());
- sbps.resize(sbps_tmp.size());
- sbpe.resize(sbpe_tmp.size());
- copy(radii_tmp.begin(),radii_tmp.end(),radii.begin());
- copy(sbps_tmp.begin(),sbps_tmp.end(),sbps.begin());
- copy(sbpe_tmp.begin(),sbpe_tmp.end(),sbpe.begin());
-
- //read cooling function data
- spline_func_obj cf;
- for(ifstream ifs(cfg.cfunc_profile.c_str());;)
- {
- assert(ifs.is_open());
- double x,y;
- ifs>>x>>y;
- if(!ifs.good())
- {
- break;
- }
- cerr<<x<<"\t"<<y<<endl;
- if(x>radii.back())
- {
- break;
- }
- cf.add_point(x,y);
- }
- cf.gen_spline();
-
- //read temperature profile data
- spline_func_obj Tprof;
- int tcnt=0;
- for(ifstream ifs1(cfg.tprofile.c_str());;++tcnt)
- {
- assert(ifs1.is_open());
- double x,y;
- ifs1>>x>>y;
- if(!ifs1.good())
- {
- break;
- }
- cerr<<x<<"\t"<<y<<endl;
-#if 0
- if(tcnt==0)
- {
- Tprof.add_point(0,y);
- }
-#endif
- Tprof.add_point(x,y);
- }
-
-
- Tprof.gen_spline();
-
- 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 fitter
- fitter<vector<double>,vector<double>,vector<double>,double> f;
- f.load_data(ds);
- //initial the object, which is used to calculate projection effect
- projector<double> a;
- bool tie_beta=false;
- if(cfg.param_map.find("beta")!=cfg.param_map.end()
- &&cfg.param_map.find("beta1")==cfg.param_map.end()
- &&cfg.param_map.find("beta2")==cfg.param_map.end())
- {
- dbeta2<double> dbetao;
- a.attach_model(dbetao);
- tie_beta=true;
- }
- else if((cfg.param_map.find("beta1")!=cfg.param_map.end()
- ||cfg.param_map.find("beta2")!=cfg.param_map.end())
- &&cfg.param_map.find("beta")==cfg.param_map.end())
- {
- dbeta<double> dbetao;
- a.attach_model(dbetao);
- tie_beta=false;
- }
- else
- {
- cerr<<"Error, cannot decide whether to tie beta together or let them vary freely!"<<endl;
- assert(0);
- }
-
- //attach the cooling function
- a.attach_cfunc(cf);
- a.set_cm_per_pixel(cm_per_pixel);
-
- f.set_model(a);
- //chi^2 statistic
- vchisq<double> c;
- c.verbose(true);
- c.set_limit();
- f.set_statistic(c);
- //optimization method
- f.set_opt_method(powell_method<double,std::vector<double> >());
- //initialize the initial values
- /*
- double =0;
- double rc1=0;
- double n02=0;
- double rc2=0;
- double beta=0;
- */
- double bkg_level=0;
-
- if(tie_beta)
- {
- f.set_param_value("beta",.7);
- f.set_param_lower_limit("beta",.3);
- f.set_param_upper_limit("beta",1.4);
- }
- else
- {
- f.set_param_value("beta1",.7);
- f.set_param_lower_limit("beta1",.3);
- f.set_param_upper_limit("beta1",1.4);
- f.set_param_value("beta2",.7);
- f.set_param_lower_limit("beta2",.3);
- f.set_param_upper_limit("beta2",1.4);
- }
- for(std::map<std::string,std::vector<double> >::iterator i=cfg.param_map.begin();
- i!=cfg.param_map.end();++i)
- {
- std::string pname=i->first;
- f.set_param_value(pname,i->second.at(0));
- if(i->second.size()==3)
- {
- double a1=i->second[1];
- double a2=i->second[2];
- double u=std::max(a1,a2);
- double l=std::min(a1,a2);
- f.set_param_upper_limit(pname,u);
- f.set_param_lower_limit(pname,l);
- }
- else
- {
- if(pname=="beta"||pname=="beta1"||pname=="beta2")
- {
- f.set_param_lower_limit(pname,.3);
- f.set_param_upper_limit(pname,1.4);
- }
- }
- }
-
-
-
- //perform the fitting, first freeze beta1, beta2, rc1, and rc2
- if(tie_beta)
- {
- f.set_param_modifier(freeze_param<std::vector<double>,std::vector<double>,std::vector<double>,std::string>("beta")+
- freeze_param<std::vector<double>,std::vector<double>,std::vector<double>,std::string>("rc1")+
- freeze_param<std::vector<double>,std::vector<double>,std::vector<double>,std::string>("rc2")
- );
- }
- else
- {
- f.set_param_modifier(freeze_param<std::vector<double>,std::vector<double>,std::vector<double>,std::string>("beta1")+
- freeze_param<std::vector<double>,std::vector<double>,std::vector<double>,std::string>("beta2")+
- freeze_param<std::vector<double>,std::vector<double>,std::vector<double>,std::string>("rc1")+
- freeze_param<std::vector<double>,std::vector<double>,std::vector<double>,std::string>("rc2")
- );
- }
-
- f.fit();
-
- f.clear_param_modifier();
-
- //then perform the fitting, freeze beta1 and beta2
- //f.set_param_modifier(freeze_param<std::vector<double>,std::vector<double>,std::vector<double>,std::string>("beta"));
- //f.set_param_modifier(freeze_param<std::vector<double>,std::vector<double>,std::vector<double>,std::string>("bkg"));
- f.fit();
- //f.clear_param_modifier();
-
- //finally thaw all parameters
- f.fit();
-
- /*
- double beta1=0;
- double beta2=0;
-
- n01=f.get_param_value("n01");
- rc1=f.get_param_value("rc1");
- n02=f.get_param_value("n02");
- rc2=f.get_param_value("rc2");
- if(tie_beta)
- {
- beta=f.get_param_value("beta");
- beta1=beta;
- beta2=beta;
- }
- else
- {
- beta1=f.get_param_value("beta1");
- beta2=f.get_param_value("beta2");
- }
- */
-
- //output the params
- ofstream param_output("lx_dbeta_param.txt");
- //output the datasets and fitting results
- for(size_t i=0;i<f.get_num_params();++i)
- {
- if(f.get_param_info(i).get_name()=="rc1")
- {
- cerr<<"rc1_kpc"<<"\t"<<abs(f.get_param_info(i).get_value())*cm_per_pixel/kpc<<endl;
- param_output<<"rc1_kpc"<<"\t"<<abs(f.get_param_info(i).get_value())*cm_per_pixel/kpc<<endl;
- }
- if(f.get_param_info(i).get_name()=="rc2")
- {
- cerr<<"rc2_kpc"<<"\t"<<abs(f.get_param_info(i).get_value())*cm_per_pixel/kpc<<endl;
- param_output<<"rc2_kpc"<<"\t"<<abs(f.get_param_info(i).get_value())*cm_per_pixel/kpc<<endl;
- }
- cerr<<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;
- }
- cerr<<"reduced_chi^2="<<f.get_statistic_value()/(radii.size()-f.get_model().get_num_free_params())<<endl;
- param_output<<"reduced_chi^2="<<f.get_statistic_value()/(radii.size()-f.get_model().get_num_free_params())<<endl;
-
- //c.verbose(false);
- //f.set_statistic(c);
- //f.fit();
- std::vector<double> p=f.get_all_params();
- f.clear_param_modifier();
- std::vector<double> mv=f.eval_model(radii,p);
-
-
- ofstream ofs_sbp("lx_sbp_fit.qdp");
- ofs_sbp<<"read serr 2"<<endl;
- ofs_sbp<<"skip single"<<endl;
-
- ofs_sbp<<"line on 2"<<endl;
- ofs_sbp<<"line on 3"<<endl;
- ofs_sbp<<"line on 4"<<endl;
- ofs_sbp<<"line on 5"<<endl;
- ofs_sbp<<"line on 7"<<endl;
- ofs_sbp<<"ls 2 on 7"<<endl;
- ofs_sbp<<"ls 2 on 3"<<endl;
- ofs_sbp<<"ls 2 on 4"<<endl;
- ofs_sbp<<"ls 2 on 5"<<endl;
-
-
-
- ofs_sbp<<"!LAB POS Y 4.00"<<endl;
- ofs_sbp<<"!LAB ROT"<<endl;
- ofs_sbp<<"win 1"<<endl;
- ofs_sbp<<"yplot 1 2 3 4 5"<<endl;
- ofs_sbp<<"loc 0 0 1 1"<<endl;
- ofs_sbp<<"vie .1 .4 .9 .9"<<endl;
- ofs_sbp<<"la y cnt/s/pixel/cm\\u2\\d"<<endl;
- ofs_sbp<<"log x"<<endl;
- ofs_sbp<<"log y"<<endl;
- ofs_sbp<<"r x "<<(radii[1]+radii[0])/2*cm_per_pixel/kpc<<" "<<(radii[sbps.size()-2]+radii[sbps.size()-1])/2*cm_per_pixel/kpc<<endl;
- ofs_sbp<<"win 2"<<endl;
- ofs_sbp<<"yplot 6 7"<<endl;
- ofs_sbp<<"loc 0 0 1 1"<<endl;
- ofs_sbp<<"vie .1 .1 .9 .4"<<endl;
- ofs_sbp<<"la x radius (kpc)"<<endl;
- ofs_sbp<<"la y chi"<<endl;
- ofs_sbp<<"log x"<<endl;
- ofs_sbp<<"log y off"<<endl;
- ofs_sbp<<"r x "<<(radii[1]+radii[0])/2*cm_per_pixel/kpc<<" "<<(radii[sbps.size()-2]+radii[sbps.size()-1])/2*cm_per_pixel/kpc<<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<<endl;
- }
- ofs_sbp<<"no no no"<<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"<<ym<<"\t"<<0<<endl;
- }
- //bkg
- ofs_sbp<<"no no no"<<endl;
- bkg_level=abs(f.get_param_value("bkg"));
- for(size_t i=0;i<sbps.size();++i)
- {
- double x=(radii[i]+radii[i-1])/2;
- ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<bkg_level<<"\t0"<<endl;
- }
- //rc1
- ofs_sbp<<"no no no"<<endl;
- double rc1_kpc=abs(f.get_param_value("rc1")*cm_per_pixel/kpc);
- double max_sbp=*max_element(sbps.begin(),sbps.end());
- double min_sbp=*min_element(sbps.begin(),sbps.end());
- for(double x=min_sbp;x<=max_sbp;x+=(max_sbp-min_sbp)/100)
- {
- ofs_sbp<<rc1_kpc<<"\t"<<x<<"\t"<<"0"<<endl;
- }
- //rc2
- ofs_sbp<<"no no no"<<endl;
- double rc2_kpc=abs(f.get_param_value("rc2")*cm_per_pixel/kpc);
- for(double x=min_sbp;x<=max_sbp;x+=(max_sbp-min_sbp)/100)
- {
- ofs_sbp<<rc2_kpc<<"\t"<<x<<"\t"<<"0"<<endl;
- }
- //resid
- ofs_sbp<<"no no no"<<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"<<(ym-sbps[i-1])/sbpe[i-1]<<"\t"<<1<<endl;
- }
- //zero level in resid map
- ofs_sbp<<"no no no"<<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"<<0<<"\t"<<0<<endl;
- }
-
-
- mv=f.eval_model_raw(radii,p);
- ofstream ofs_rho("lx_rho_fit.qdp");
- ofstream ofs_rho_data("lx_rho_fit.dat");
- //ofstream ofs_entropy("entropy.qdp");
- ofs_rho<<"la x radius (kpc)"<<endl;
- ofs_rho<<"la y density (cm\\u-3\\d)"<<endl;
- /*
- for(int 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;
- }
- */
-
- p.back()=0;
-
- radii.clear();
- double rout=atof(argv[2])*kpc;
- for(double r=0;r<rout;r+=1*kpc)//step size=1kpc
- {
- double r_pix=r/cm_per_pixel;
- radii.push_back(r_pix);
- }
-
- double Da=cm_per_pixel/(.492/3600./180.*pi);
- double Dl=Da*(1+z)*(1+z);
- cout<<"dl="<<Dl/kpc<<endl;
- for(int n=3;n<argc;++n)
- {
- spline_func_obj cf_erg;
- for(ifstream ifs(argv[n]);;)
- {
- assert(ifs.is_open());
- double x,y;
- ifs>>x>>y;
- if(!ifs.good())
- {
- break;
- }
- //cerr<<x<<"\t"<<y<<endl;
-
- cf_erg.add_point(x,y);//change with source
- }
- cf_erg.gen_spline();
-
- projector<double>& pj=dynamic_cast<projector<double>&>(f.get_model());
- pj.attach_cfunc(cf_erg);
-
- mv=f.eval_model_raw(radii,p);
- double flux_erg=0;
- for(size_t i=0;i<radii.size()-1;++i)
- {
- double S=pi*(radii[i+1]+radii[i])*(radii[i+1]-radii[i]);
- flux_erg+=S*mv[i];
- }
- cout<<flux_erg*4*pi*Dl*Dl<<endl;
- cout<<flux_erg<<endl;
- param_output<<"Lx"<<n-2<<"\t"<<flux_erg*4*pi*Dl*Dl<<endl;
- param_output<<"Fx"<<n-2<<"\t"<<flux_erg<<endl;
- }
-}
diff --git a/mass_profile/chisq.hpp b/mass_profile/chisq.hpp
deleted file mode 100644
index 0e58200..0000000
--- a/mass_profile/chisq.hpp
+++ /dev/null
@@ -1,206 +0,0 @@
-/**
- \file chisq.hpp
- \brief chi-square statistic
- \author Junhua Gu
- */
-
-#ifndef CHI_SQ_HPP
-#define CHI_SQ_HPP
-
-#define OPT_HEADER
-
-#include <core/fitter.hpp>
-#include <iostream>
-#include <vector>
-#include <misc/optvec.hpp>
-#include <cmath>
-
-using std::cerr;
-using std::endl;
-
-namespace opt_utilities
-{
- static const int display_interval=10;
- /**
- \brief chi-square statistic
- \tparam Ty the return type of model
- \tparam Tx the type of the self-var
- \tparam Tp the type of model parameter
- \tparam Ts the type of the statistic
- \tparam Tstr the type of the string used
- */
- template<typename Ty,typename Tx,typename Tp,typename Ts,typename Tstr>
- class chisq
- :public statistic<Ty,Tx,Tp,Ts,Tstr>
- {
- };
- template<>
- class chisq<double,double,std::vector<double>,double,std::string>
- :public statistic<double,double,std::vector<double> ,double,std::string>
- {
- public:
- typedef double Ty;
- typedef double Tx;
- typedef std::vector<double> Tp;
- typedef double Ts;
- typedef std::string Tstr;
- private:
- bool verb;
- bool limit_bound;
- int n;
-
- statistic<Ty,Tx,Tp,Ts,Tstr>* do_clone()const
- {
- // return const_cast<statistic<Ty,Tx,Tp>*>(this);
- return new chisq<Ty,Tx,Tp,Ts,Tstr>(*this);
- }
-
- const char* do_get_type_name()const
- {
- return "chi^2 statistics (specialized for double)";
- }
- public:
- void verbose(bool v)
- {
- verb=v;
- }
-
- void set_limit()
- {
- limit_bound=true;
- }
-
- void clear_limit()
- {
- limit_bound=false;
- }
- public:
- chisq()
- :verb(true),limit_bound(false)
- {}
-
- Ty do_eval(const Tp& p)
- {
- if(limit_bound)
- {
- Tp p1=this->get_fitter().get_model().reform_param(p);
- for(size_t i=0;i<p1.size();++i)
- {
- if(p1[i]>this->get_fitter().get_param_info(i).get_upper_limit()||
- p1[i]<this->get_fitter().get_param_info(i).get_lower_limit())
- {
- return 1e99;
- }
- }
- }
- Ty result(0);
- std::vector<float> vx;
- std::vector<float> vy;
- std::vector<float> vye1;
- std::vector<float> vye2;
- std::vector<float> my;
- float xmin=1e99,xmax=-1e99,ymin=1e99,ymax=-1e99;
- if(verb)
- {
- n++;
- if(n%display_interval==0)
- {
- vx.resize(this->get_data_set().size());
- vy.resize(this->get_data_set().size());
- vye1.resize(this->get_data_set().size());
- vye2.resize(this->get_data_set().size());
- my.resize(this->get_data_set().size());
- }
- }
-
- for(int i=(this->get_data_set()).size()-1;i>=0;--i)
- {
-
-#ifdef HAVE_X_ERROR
- Tx x1=this->get_data_set().get_data(i).get_x()-this->get_data_set().get_data(i).get_x_lower_err();
- Tx x2=this->get_data_set().get_data(i).get_x()+this->get_data_set().get_data(i).get_x_upper_err();
- Tx x=this->get_data_set().get_data(i).get_x();
- Ty errx1=(this->eval_model(x1,p)-this->eval_model(x,p));
- Ty errx2=(this->eval_model(x2,p)-this->eval_model(x,p));
- //Ty errx=0;
-#else
- Ty errx1=0;
- Ty errx2=0;
-#endif
-
- Ty y_model=this->eval_model(this->get_data_set().get_data(i).get_x(),p);
- Ty y_obs=this->get_data_set().get_data(i).get_y();
- Ty y_err;
-
- Ty errx=0;
- if(errx1<errx2)
- {
- if(y_obs<y_model)
- {
- errx=errx1>0?errx1:-errx1;
- }
- else
- {
- errx=errx2>0?errx2:-errx2;
- }
- }
- else
- {
- if(y_obs<y_model)
- {
- errx=errx2>0?errx2:-errx2;
- }
- else
- {
- errx=errx1>0?errx1:-errx1;
- }
- }
-
- if(y_model>y_obs)
- {
- y_err=this->get_data_set().get_data(i).get_y_upper_err();
- }
- else
- {
- y_err=this->get_data_set().get_data(i).get_y_lower_err();
- }
-
- Ty chi=(y_obs-y_model)/std::sqrt(y_err*y_err+errx*errx);
-
- result+=chi*chi;
-
- if(verb&&n%display_interval==0)
- {
- vx.at(i)=this->get_data_set().get_data(i).get_x();
- vy.at(i)=this->get_data_set().get_data(i).get_y();
- vye1.at(i)=std::abs(this->get_data_set().get_data(i).get_y_lower_err());
- vye2.at(i)=std::abs(this->get_data_set().get_data(i).get_y_upper_err());
- my.at(i)=y_model;
-
- xmin=std::min(vx.at(i),xmin);
- ymin=std::min(vy.at(i),ymin-vye1[i]);
- xmax=std::max(vx.at(i),xmax);
- ymax=std::max(vy.at(i),ymax+vye2[i]);
- }
-
- }
- if(verb)
- {
- if(n%display_interval==0)
- {
- cerr<<result<<"\t";
- for(size_t i=0;i<get_size(p);++i)
- {
- cerr<<get_element(p,i)<<",";
- }
- cerr<<endl;
- //cerr<<x1<<"\t"<<x2<<endl;
- }
- }
-
- return result;
- }
- };
-}
-
-#endif /* CHI_SQ_HPP */
diff --git a/mass_profile/dbeta.hpp b/mass_profile/dbeta.hpp
deleted file mode 100644
index 7246b44..0000000
--- a/mass_profile/dbeta.hpp
+++ /dev/null
@@ -1,108 +0,0 @@
-#ifndef DBETA
-#define DBETA
-#include "projector.hpp"
-
-/**
- dbeta: double beta model for density
- dbeta2: double beta model for density with only one beta
-*/
-
-
-namespace opt_utilities
-{
- template <typename T>
- class dbeta
- :public model<std::vector<T>,std::vector<T>,std::vector<T> >
- {
- public:
- dbeta()
- {
- this->push_param_info(param_info<std::vector<T>,std::string>("n01",1));
- this->push_param_info(param_info<std::vector<T>,std::string>("beta1",.66));
- this->push_param_info(param_info<std::vector<T>,std::string>("rc1",100));
-
- this->push_param_info(param_info<std::vector<T>,std::string>("n02",1));
- this->push_param_info(param_info<std::vector<T>,std::string>("beta2",.67));
- this->push_param_info(param_info<std::vector<T>,std::string>("rc2",110));
-
- }
-
- public:
- dbeta<T>* do_clone()const
- {
- return new dbeta<T>(*this);
- }
-
- std::vector<T> do_eval(const std::vector<T> & x,
- const std::vector<T>& p)
- {
- T n01=std::abs(p[0]);
- T beta1=p[1];
- T rc1=p[2];
-
- T n02=std::abs(p[3]);
- T beta2=p[4];
- T rc2=p[5];
-
-
-
- std::vector<T> result(x.size()-1);
- for(size_t i=1;i<x.size();++i)
- {
- T xi=(x[i]+x[i-1])/2;
- T yi=0;
- yi=n01*pow(1+xi*xi/rc1/rc1,-3./2.*beta1)+n02*pow(1+xi*xi/rc2/rc2,-3./2.*beta2);
- result[i-1]=yi;
- }
- return result;
- }
- };
-
- template <typename T>
- class dbeta2
- :public model<std::vector<T>,std::vector<T>,std::vector<T> >
- {
- public:
- dbeta2()
- {
- this->push_param_info(param_info<std::vector<T>,std::string>("n01",1));
- this->push_param_info(param_info<std::vector<T>,std::string>("rc1",100));
- this->push_param_info(param_info<std::vector<T>,std::string>("n02",1));
- this->push_param_info(param_info<std::vector<T>,std::string>("rc2",110));
- this->push_param_info(param_info<std::vector<T>,std::string>("beta",.67));
-
- }
-
- public:
- dbeta2<T>* do_clone()const
- {
- return new dbeta2<T>(*this);
- }
-
- std::vector<T> do_eval(const std::vector<T> & x,
- const std::vector<T>& p)
- {
- T n01=std::abs(p[0]);
- T rc1=p[1];
-
- T n02=std::abs(p[2]);
- T rc2=p[3];
- T beta=p[4];
- T beta1=beta;
- T beta2=beta;
-
- std::vector<T> result(x.size()-1);
- for(size_t i=1;i<x.size();++i)
- {
- T xi=(x[i]+x[i-1])/2;
- T yi=0;
- yi=n01*pow(1+xi*xi/rc1/rc1,-3./2.*beta1)+n02*pow(1+xi*xi/rc2/rc2,-3./2.*beta2);
- result[i-1]=yi;
- }
- return result;
- }
- };
-
-}
-
-#endif
diff --git a/mass_profile/dump_fit_qdp.cpp b/mass_profile/dump_fit_qdp.cpp
deleted file mode 100644
index 85307d1..0000000
--- a/mass_profile/dump_fit_qdp.cpp
+++ /dev/null
@@ -1,55 +0,0 @@
-#include "dump_fit_qdp.hpp"
-
-namespace opt_utilities
-{
- const static double kpc=3.086E21;
- void dump_sbp_beta(std::ostream& os,fitter<double,double,std::vector<double>,double,std::string>& f,double cm_per_pixel,const std::vector<double>& r,const std::vector<double>& y,const std::vector<double>& ye)
- {
- os<<"read serr 1 2"<<std::endl;
- os<<"skip single"<<std::endl;
- os<<"la x \"radius (kpc)\""<<std::endl;
- os<<"la y \"surface brightness (cts s\\u-1\\d pixel\\u-2\\d)\""<<std::endl;
- os<<"li on 2"<<std::endl;
- for(size_t i=1;i<r.size();++i)
- {
- os<<(r[i]+r[i-1])/2*cm_per_pixel/kpc<<"\t"<<(r[i]-r[i-1])/2*cm_per_pixel/kpc<<"\t"<<y[i-1]<<"\t"<<ye[i-1]<<std::endl;
- }
- os<<"no no no"<<std::endl;
- std::vector<double> p=f.get_all_params();
- for(size_t i=1;i<r.size();++i)
- {
- double x=(r[i]+r[i-1])/2;
- os<<x*cm_per_pixel/kpc<<"\t"<<0<<"\t"<<f.eval_model_raw(x,p)<<"\t"<<0<<std::endl;
- }
- }
-
- void dump_rho_beta(std::ostream& os,fitter<std::vector<double>,std::vector<double>,std::vector<double>,double,std::string>& f,double cm_per_pixel,const std::vector<double>& r,const std::vector<double>& sbps,const std::vector<double>& sbpe)
- {
- os<<"read serr 1 2"<<std::endl;
- os<<"skip single"<<std::endl;
- os<<"la x \"radius (kpc)\""<<std::endl;
- os<<"la y \"density (cm\\u-3\\d)\""<<std::endl;
- os<<"li on 2"<<std::endl;
-
- for(size_t i=1;i<r.size();++i)
- {
- double x=(r[i]+r[i-1])/2;
- double y=sbps[i-1];
- double ye=sbpe[i-1];
- os<<x*cm_per_pixel/kpc<<"\t0\t"<<y<<"\t"<<ye<<std::endl;
- }
-
- os<<"no no no"<<std::endl;
- std::vector<double> p=f.get_all_params();
- std::vector<double> mv=f.eval_model_raw(r,p);
- for(size_t i=1;i<r.size();++i)
- {
- double x=(r[i]+r[i-1])/2;
- double y=mv[i-1];
- double ye=0;
- os<<x*cm_per_pixel/kpc<<"\t0\t"<<y<<"\t"<<ye<<std::endl;
- }
-
- }
-
-};
diff --git a/mass_profile/dump_fit_qdp.hpp b/mass_profile/dump_fit_qdp.hpp
deleted file mode 100644
index 8364dfd..0000000
--- a/mass_profile/dump_fit_qdp.hpp
+++ /dev/null
@@ -1,16 +0,0 @@
-#ifndef DUMP_FIT_QDP_HPP
-#define DUMP_FIT_QDP_HPP
-
-#include <core/fitter.hpp>
-#include <vector>
-#include <iostream>
-#include <string>
-
-namespace opt_utilities
-{
- void dump_sbp_beta(std::ostream& os,fitter<double,double,std::vector<double>,double,std::string>& f,double cm_per_pixel,const std::vector<double>& r,const std::vector<double>& y,const std::vector<double>& ye);
- void dump_rho_beta(std::ostream& os,fitter<std::vector<double>,std::vector<double>,std::vector<double>,double,std::string>& f,double cm_per_pixel,const std::vector<double>& r,const std::vector<double>& sbps,const std::vector<double>& sbpe);
- void dump_rho_dbeta(std::ostream& os,fitter<std::vector<double>,std::vector<double>,std::vector<double>,double,std::string>& f,double cm_per_pixel);
-};
-
-#endif
diff --git a/mass_profile/fit_beta_sbp.cpp b/mass_profile/fit_beta_sbp.cpp
deleted file mode 100644
index 295fa1e..0000000
--- a/mass_profile/fit_beta_sbp.cpp
+++ /dev/null
@@ -1,534 +0,0 @@
-/*
- Perform a double-beta density model fitting to the surface brightness data
- Author: Junhua Gu
- Last modified: 2016.06.07
- This code is distributed with no warrant
-*/
-
-#include <iostream>
-#include <fstream>
-#include <sstream>
-#include <list>
-#include <algorithm>
-#include "beta_cfg.hpp"
-#include "dump_fit_qdp.hpp"
-#include "report_error.hpp"
-#include "vchisq.hpp"
-#include "chisq.hpp"
-#include "beta.hpp"
-#include "models/beta1d.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 kpc=3.086E21;//kpc in cm
-const double Mpc=kpc*1000;
-
-double beta_func(double r, double n0, double rc, double beta)
-{
- return abs(n0) * pow(1+r*r/rc/rc, -3./2.*abs(beta));
-}
-
-//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 G=6.673E-8;//cm^3 g^-1 s^2
- 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 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)
- {
- 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);
- }
-};
-
-int main(int argc,char* argv[])
-{
- if(argc!=2)
- {
- cerr<<argv[0]<<" <configure file>"<<endl;
- return -1;
- }
- //initialize the parameters list
- ifstream cfg_file(argv[1]);
- assert(cfg_file.is_open());
- cfg_map cfg=parse_cfg_file(cfg_file);
-
- const double z=cfg.z;
-
- //initialize the radius list, sbp list and sbp error list
- std::vector<double> radii;
- std::vector<double> sbps;
- std::vector<double> sbpe;
- std::vector<double> radii_all;
- std::vector<double> sbps_all;
- std::vector<double> sbpe_all;
- //prepend the zero point to radius list
- radii.push_back(0.0);
- radii_all.push_back(0.0);
- //read sbp and sbp error data
- cerr << "Read surface brightness profile data ..." << endl;
- ifstream ifs(cfg.sbp_data.c_str());
- std::string line;
- if (ifs.is_open())
- {
- while(std::getline(ifs, line))
- {
- if (line.empty())
- continue;
-
- std::istringstream iss(line);
- double x, xe, y, ye;
- if ((iss >> x >> xe >> y >> ye))
- {
- std::cerr << "sbprofile data: "
- << x << ", " << xe << ", " << y << ", " << ye
- << std::endl;
- radii.push_back(x+xe); /* NOTE: use outer radii of regions */
- radii_all.push_back(x+xe);
- sbps.push_back(y);
- sbps_all.push_back(y);
- sbpe.push_back(ye);
- sbpe_all.push_back(ye);
- }
- else
- {
- std::cerr << "skipped line: " << line << std::endl;
- }
- }
- }
- else
- {
- std::cerr << "ERROR: cannot open file: " << cfg.sbp_data.c_str()
- << std::endl;
- return 1;
- }
-
- //initialize the cm/pixel value
- double cm_per_pixel=cfg.cm_per_pixel;
- double rmin;
- if(cfg.rmin_pixel>0)
- {
- rmin=cfg.rmin_pixel;
- }
- else
- {
- rmin=cfg.rmin_kpc*kpc/cm_per_pixel;
- }
-
- cerr<<"rmin="<<rmin<<" (pixel)"<<endl;
- std::list<double> radii_tmp,sbps_tmp,sbpe_tmp;
- radii_tmp.resize(radii.size());
- sbps_tmp.resize(sbps.size());
- sbpe_tmp.resize(sbpe.size());
- copy(radii.begin(),radii.end(),radii_tmp.begin());
- copy(sbps.begin(),sbps.end(),sbps_tmp.begin());
- copy(sbpe.begin(),sbpe.end(),sbpe_tmp.begin());
- for(list<double>::iterator i=radii_tmp.begin();i!=radii_tmp.end();)
- {
- if(*i<rmin)
- {
- radii_tmp.pop_front();
- sbps_tmp.pop_front();
- sbpe_tmp.pop_front();
- i=radii_tmp.begin();
- continue;
- }
- ++i;
- }
- radii.resize(radii_tmp.size());
- sbps.resize(sbps_tmp.size());
- sbpe.resize(sbpe_tmp.size());
- copy(radii_tmp.begin(),radii_tmp.end(),radii.begin());
- copy(sbps_tmp.begin(),sbps_tmp.end(),sbps.begin());
- copy(sbpe_tmp.begin(),sbpe_tmp.end(),sbpe.begin());
-
- //read cooling function data
- cerr << "Read cooling function profile data ..." << endl;
- spline_func_obj cf;
- for(ifstream ifs(cfg.cfunc_profile.c_str());;)
- {
- assert(ifs.is_open());
- double x,y;
- ifs>>x>>y;
- if(!ifs.good())
- {
- break;
- }
- cerr<<x<<"\t"<<y<<endl;
- if(x>radii.back())
- {
- cerr << "radius_max: " << radii.back() << endl;
- break;
- }
- cf.add_point(x,y);
- }
- cf.gen_spline();
-
- //read temperature profile data
- cerr << "Read temperature profile data ..." << endl;
- spline_func_obj Tprof;
- int tcnt=0;
- for(ifstream ifs1(cfg.tprofile.c_str());;++tcnt)
- {
- assert(ifs1.is_open());
- double x,y;
- ifs1>>x>>y;
- if(!ifs1.good())
- {
- break;
- }
- cerr<<x<<"\t"<<y<<endl;
-#if 0
- if(tcnt==0)
- {
- Tprof.add_point(0,y);
- }
-#endif
- Tprof.add_point(x,y);
- }
- Tprof.gen_spline();
-
- 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 fitter
- fitter<vector<double>,vector<double>,vector<double>,double> f;
- f.load_data(ds);
- //initial the object, which is used to calculate projection effect
- projector<double> a;
- beta<double> betao;
- //attach the cooling function
- a.attach_cfunc(cf);
- a.set_cm_per_pixel(cm_per_pixel);
- a.attach_model(betao);
- f.set_model(a);
- //chi^2 statistic
- vchisq<double> c;
- c.verbose(true);
- c.set_limit();
- f.set_statistic(c);
- //optimization method
- f.set_opt_method(powell_method<double,std::vector<double> >());
- //initialize the initial values
- double n0=0;
- //double beta=atof(arg_map["beta"].c_str());
- double beta=0;
- double rc=0;
- double bkg_level=0;
-
- for(std::map<std::string,std::vector<double> >::iterator i=cfg.param_map.begin();
- i!=cfg.param_map.end();++i)
- {
- std::string pname=i->first;
- f.set_param_value(pname,i->second.at(0));
- if(i->second.size()==3)
- {
- double a1=i->second[1];
- double a2=i->second[2];
- double u=std::max(a1,a2);
- double l=std::min(a1,a2);
- f.set_param_upper_limit(pname,u);
- f.set_param_lower_limit(pname,l);
- }
- else
- {
- if(pname=="beta")
- {
- f.set_param_lower_limit(pname,.3);
- f.set_param_upper_limit(pname,1.4);
- }
- }
- }
-
- f.fit();
- f.fit();
- std::vector<double> p=f.get_all_params();
- n0=f.get_param_value("n0");
- rc=f.get_param_value("rc");
- beta=f.get_param_value("beta");
- //output the datasets and fitting results
- ofstream param_output("beta_param.txt");
- for(size_t i=0;i<f.get_num_params();++i)
- {
- if(f.get_param_info(i).get_name()=="rc")
- {
- cerr<<"rc_kpc"<<"\t"<<abs(f.get_param_info(i).get_value())*cm_per_pixel/kpc<<endl;
- param_output<<"rc_kpc"<<"\t"<<abs(f.get_param_info(i).get_value())*cm_per_pixel/kpc<<endl;
- }
- cerr<<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;
- }
- cerr<<"reduced_chi^2="<<f.get_statistic_value()/(radii.size()-f.get_model().get_num_free_params())<<endl;
- param_output<<"reduced_chi^2="<<f.get_statistic_value()/(radii.size()-f.get_model().get_num_free_params())<<endl;
-
- std::vector<double> mv=f.eval_model_raw(radii_all,p);
- int sbps_inner_cut_size=int(sbps_all.size()-sbps.size());
- ofstream ofs_sbp("sbp_fit.qdp");
- ofs_sbp<<"read serr 2"<<endl;
- ofs_sbp<<"skip single"<<endl;
- ofs_sbp<<"line off "<<endl;
- if(sbps_inner_cut_size>=1)
- {
- ofs_sbp<<"line on 2"<<endl;
- ofs_sbp<<"line on 3"<<endl;
- ofs_sbp<<"line on 4"<<endl;
- ofs_sbp<<"line on 5"<<endl;
- ofs_sbp<<"ls 2 on 2"<<endl;
- ofs_sbp<<"ls 2 on 4"<<endl;
- ofs_sbp<<"ls 2 on 5"<<endl;
- ofs_sbp<<"line on 7"<<endl;
- ofs_sbp<<"ls 2 on 7"<<endl;
-
- ofs_sbp<<"ma 1 on 2"<<endl;
- ofs_sbp<<"color 1 on 1"<<endl;
- ofs_sbp<<"color 2 on 2"<<endl;
- ofs_sbp<<"color 3 on 3"<<endl;
- ofs_sbp<<"color 4 on 4"<<endl;
- ofs_sbp<<"color 5 on 5"<<endl;
-
- ofs_sbp<<"win 1"<<endl;
- ofs_sbp<<"yplot 1 2 3 4 5"<<endl;
- ofs_sbp<<"loc 0 0 1 1"<<endl;
- ofs_sbp<<"vie .1 .4 .9 .9"<<endl;
- ofs_sbp<<"la y cnt/s/pixel/cm^2"<<endl;
- ofs_sbp<<"log x"<<endl;
- ofs_sbp<<"log y"<<endl;
- ofs_sbp<<"r x "<<(radii[1]+radii[0])/2*cm_per_pixel/kpc<<" "<<(radii[sbps.size()-2]+radii[sbps.size()-1])/2*cm_per_pixel/kpc<<endl;
- ofs_sbp<<"win 2"<<endl;
- ofs_sbp<<"yplot 6 7"<<endl;
- ofs_sbp<<"loc 0 0 1 1"<<endl;
- ofs_sbp<<"vie .1 .1 .9 .4"<<endl;
- ofs_sbp<<"la x radius (kpc)"<<endl;
- ofs_sbp<<"la y chi"<<endl;
- ofs_sbp<<"log y off"<<endl;
- ofs_sbp<<"log x"<<endl;
- ofs_sbp<<"r x "<<(radii[1]+radii[0])/2*cm_per_pixel/kpc<<" "<<(radii[sbps.size()-2]+radii[sbps.size()-1])/2*cm_per_pixel/kpc<<endl;
- }
- else
- {
- ofs_sbp<<"line on 2"<<endl;
- ofs_sbp<<"line on 3"<<endl;
- ofs_sbp<<"line on 4"<<endl;
- ofs_sbp<<"ls 2 on 3"<<endl;
- ofs_sbp<<"ls 2 on 4"<<endl;
- ofs_sbp<<"line on 6"<<endl;
- ofs_sbp<<"ls 2 on 6"<<endl;
-
- ofs_sbp<<"color 1 on 1"<<endl;
- ofs_sbp<<"color 3 on 2"<<endl;
- ofs_sbp<<"color 4 on 3"<<endl;
- ofs_sbp<<"color 5 on 4"<<endl;
- //ofs_sbp<<"ma 1 on 2"<<endl;
-
- ofs_sbp<<"win 1"<<endl;
- ofs_sbp<<"yplot 1 2 3 4"<<endl;
- ofs_sbp<<"loc 0 0 1 1"<<endl;
- ofs_sbp<<"vie .1 .4 .9 .9"<<endl;
- ofs_sbp<<"la y cnt/s/pixel/cm^2"<<endl;
- ofs_sbp<<"log x"<<endl;
- ofs_sbp<<"log y"<<endl;
- ofs_sbp<<"r x "<<(radii[1]+radii[0])/2*cm_per_pixel/kpc<<" "<<(radii[radii.size()-2]+radii[radii.size()-1])/2*cm_per_pixel/kpc<<endl;
- ofs_sbp<<"win 2"<<endl;
- ofs_sbp<<"yplot 5 6"<<endl;
- ofs_sbp<<"loc 0 0 1 1"<<endl;
- ofs_sbp<<"vie .1 .1 .9 .4"<<endl;
- ofs_sbp<<"la x radius (kpc)"<<endl;
- ofs_sbp<<"la y chi"<<endl;
- ofs_sbp<<"log x"<<endl;
- ofs_sbp<<"log y off"<<endl;
- ofs_sbp<<"r x "<<(radii[1]+radii[0])/2*cm_per_pixel/kpc<<" "<<(radii[radii.size()-2]+radii[radii.size()-1])/2*cm_per_pixel/kpc<<endl;
-
- }
- // cout<<sbps_all.size()<<"\t"<<sbps.size()<<"\t"<<sbps_inner_cut_size<<endl;
- for(size_t i=1;i<sbps_all.size();++i)
- {
- double x=(radii_all[i]+radii_all[i-1])/2;
- double y=sbps_all[i-1];
- double ye=sbpe_all[i-1];
- ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<y<<"\t"<<ye<<endl;
- }
- if(sbps_inner_cut_size>=1)
- {
- ofs_sbp<<"no no no"<<endl;
- for(int i=1;i<sbps_inner_cut_size+1;++i)
- {
- double x=(radii_all[i]+radii_all[i-1])/2;
- double ym=mv[i-1];
- ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<ym<<"\t"<<"0"<<endl;
- }
- }
- ofs_sbp<<"no no no"<<endl;
- for(size_t i=sbps_inner_cut_size;i<sbps_all.size();++i)
- {
- double x=(radii_all[i]+radii_all[i-1])/2;
- double ym=mv[i-1];
- ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<ym<<"\t"<<"0"<<endl;
- }
- ofs_sbp<<"no no no"<<endl;
- //bkg level
- bkg_level=abs(f.get_param_value("bkg"));
- for(size_t i=0;i<sbps_all.size();++i)
- {
- double x=(radii_all[i]+radii_all[i-1])/2;
- ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<bkg_level<<"\t0"<<endl;
- }
- ofs_sbp<<"no no no"<<endl;
- //rc
- double rc_kpc=abs(f.get_param_value("rc")*cm_per_pixel/kpc);
- double max_sbp=*max_element(sbps_all.begin(),sbps_all.end());
- double min_sbp=*min_element(sbps_all.begin(),sbps_all.end());
- for(double x=min_sbp;x<=max_sbp;x+=(max_sbp-min_sbp)/100)
- {
- ofs_sbp<<rc_kpc<<"\t"<<x<<"\t"<<"0"<<endl;
- }
- //resid
- ofs_sbp<<"no no no"<<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"<<(ym-sbps[i-1])/sbpe[i-1]<<"\t"<<1<<endl;
- }
-
- //zero level of resid
- ofs_sbp<<"no no no"<<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"<<0<<"\t"<<0<<endl;
- }
-
- mv=betao.eval(radii,p);
- ofstream ofs_rho("rho_fit.qdp");
- ofstream ofs_rho_data("rho_fit.dat");
- ofstream ofs_entropy("entropy.qdp");
- ofs_rho<<"la x radius (kpc)"<<endl;
- ofs_rho<<"la y density (cm\\u-3\\d)"<<endl;
- /*
- for(int 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;
- }
- */
-
- //double lower,upper;
- double dr=1;
- //calculate the mass profile
- //const double G=6.673E-8;//cm^3 g^-1 s^-2
- // Molecular weight per electron
- // Reference: Ettori et al. 2013, Space Sci. Rev., 177, 119-154; Eq.(9) below
- static const double mu=1.155;
- static const double mp=1.67262158E-24;//g
- static const double M_sun=1.98892E33;//g
- //static const double k=1.38E-16;
-
- ofstream ofs_mass("mass_int.qdp");
- ofstream ofs_mass_dat("mass_int.dat");
- ofstream ofs_overdensity("overdensity.qdp");
- ofstream ofs_gas_mass("gas_mass_int.qdp");
- //ofs_mass<<"la x radius (kpc)"<<endl;
- //ofs_mass<<"la y mass enclosed (solar mass)"<<endl;
- //ofs_overdensity<<"la x radius (kpc)"<<endl;
- //ofs_overdensity<<"la y overdensity"<<endl;
- double gas_mass=0;
- for(double r=1;r<200000;r+=dr)
- {
- dr=r/100;
- double r1=r+dr;
- double r_cm=r*cm_per_pixel;
- double r1_cm=r1*cm_per_pixel;
- double dr_cm=dr*cm_per_pixel;
- double V_cm3=4./3.*pi*(dr_cm*(r1_cm*r1_cm+r_cm*r_cm+r_cm*r1_cm));
- double ne=beta_func(r,n0,rc,beta);//cm^-3
-
- double dmgas=V_cm3*ne*mu*mp/M_sun;
- gas_mass+=dmgas;
-
- ofs_gas_mass<<r*cm_per_pixel/kpc<<"\t"<<gas_mass<<endl;
- double ne1=beta_func(r1,n0,rc,beta);//cm^3
-
- double T_keV=Tprof(r);
- double T1_keV=Tprof(r1);
-
- //double T_K=T_keV*11604505.9;
- //double T1_K=T1_keV*11604505.9;
-
- double dlnT=log(T1_keV/T_keV);
- double dlnr=log(r+dr)-log(r);
- double dlnn=log(ne1/ne);
-
- //double r_kpc=r_cm/kpc;
- double r_Mpc=r_cm/Mpc;
- //double M=-r_cm*T_K*k/G/mu/mp*(dlnT/dlnr+dlnn/dlnr);
- //ref:http://adsabs.harvard.edu/abs/2012MNRAS.422.3503W
- //Walker et al. 2012
- double M=-3.68E13*M_sun*T_keV*r_Mpc*(dlnT/dlnr+dlnn/dlnr);
- double rho=M/(4./3.*pi*r_cm*r_cm*r_cm);
-
- double S=T_keV/pow(ne,2./3.);
- //cout<<r<<"\t"<<M/M_sun<<endl;
- //cout<<r<<"\t"<<T_keV<<endl;
-
- ofs_rho<<r*cm_per_pixel/kpc<<"\t"<<ne<<endl;
- ofs_rho_data<<r*cm_per_pixel/kpc<<"\t"<<ne<<endl;
- ofs_entropy<<r*cm_per_pixel/kpc<<"\t"<<S<<endl;
-#if 0
- if(r*cm_per_pixel/kpc<5)
- {
- continue;
- }
-#endif
- ofs_mass<<r*cm_per_pixel/kpc<<"\t"<<M/M_sun<<endl;
- if(r<radii.at(sbps.size()))
- {
- ofs_mass_dat<<r*cm_per_pixel/kpc<<"\t0\t"<<M/M_sun<<"\t"<<M/M_sun*.1<<endl;
- }
- ofs_overdensity<<r*cm_per_pixel/kpc<<"\t"<<rho/calc_critical_density(z)<<endl;
-
- }
-}
diff --git a/mass_profile/fit_dbeta_sbp.cpp b/mass_profile/fit_dbeta_sbp.cpp
deleted file mode 100644
index 71b3089..0000000
--- a/mass_profile/fit_dbeta_sbp.cpp
+++ /dev/null
@@ -1,586 +0,0 @@
-/*
- Perform a double-beta density model fitting to the surface brightness data
- Author: Junhua Gu
- Last modified: 2016.06.07
- This code is distributed with no warrant
-*/
-
-
-#include <iostream>
-#include <fstream>
-#include <sstream>
-#include <list>
-#include "vchisq.hpp"
-#include "dbeta.hpp"
-#include "beta_cfg.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 kpc=3.086E21;//kpc in cm
-const double Mpc=kpc*1000;
-
-double dbeta_func(double r, double n01, double rc1, double beta1,
- double n02, double rc2, double beta2)
-{
- double v1 = abs(n01) * pow(1+r*r/rc1/rc1, -3./2.*abs(beta1));
- double v2 = abs(n02) * pow(1+r*r/rc2/rc2, -3./2.*abs(beta2));
- return v1 + v2;
-}
-
-//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 G=6.673E-8;//cm^3 g^-1 s^2
- 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 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)
- {
- 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);
- }
-};
-
-int main(int argc,char* argv[])
-{
- if(argc!=2)
- {
- cerr<<argv[0]<<" <configure file>"<<endl;
- return -1;
- }
- //initialize the parameters list
- ifstream cfg_file(argv[1]);
- assert(cfg_file.is_open());
- cfg_map cfg=parse_cfg_file(cfg_file);
-
- const double z=cfg.z;
-
- //initialize the radius list, sbp list and sbp error list
- std::vector<double> radii;
- std::vector<double> sbps;
- std::vector<double> sbpe;
- std::vector<double> radii_all;
- std::vector<double> sbps_all;
- std::vector<double> sbpe_all;
- //prepend the zero point to radius list
- radii.push_back(0.0);
- radii_all.push_back(0.0);
- //read sbp and sbp error data
- ifstream ifs(cfg.sbp_data.c_str());
- std::string line;
- if (ifs.is_open())
- {
- while(std::getline(ifs, line))
- {
- if (line.empty())
- continue;
-
- std::istringstream iss(line);
- double x, xe, y, ye;
- if ((iss >> x >> xe >> y >> ye))
- {
- std::cerr << "sbprofile data: "
- << x << ", " << xe << ", " << y << ", " << ye
- << std::endl;
- radii.push_back(x+xe); /* NOTE: use outer radii of regions */
- radii_all.push_back(x+xe);
- sbps.push_back(y);
- sbps_all.push_back(y);
- sbpe.push_back(ye);
- sbpe_all.push_back(ye);
- }
- else
- {
- std::cerr << "skipped line: " << line << std::endl;
- }
- }
- }
- else
- {
- std::cerr << "ERROR: cannot open file: " << cfg.sbp_data.c_str()
- << std::endl;
- return 1;
- }
-
- //initialize the cm/pixel value
- double cm_per_pixel=cfg.cm_per_pixel;
- double rmin;
- if(cfg.rmin_pixel>0)
- {
- rmin=cfg.rmin_pixel;
- }
- else
- {
- rmin=cfg.rmin_kpc*kpc/cm_per_pixel;
- }
-
- cerr<<"rmin="<<rmin<<" (pixel)"<<endl;
- std::list<double> radii_tmp,sbps_tmp,sbpe_tmp;
- radii_tmp.resize(radii.size());
- sbps_tmp.resize(sbps.size());
- sbpe_tmp.resize(sbpe.size());
- copy(radii.begin(),radii.end(),radii_tmp.begin());
- copy(sbps.begin(),sbps.end(),sbps_tmp.begin());
- copy(sbpe.begin(),sbpe.end(),sbpe_tmp.begin());
- for(list<double>::iterator i=radii_tmp.begin();i!=radii_tmp.end();)
- {
- if(*i<rmin)
- {
- radii_tmp.pop_front();
- sbps_tmp.pop_front();
- sbpe_tmp.pop_front();
- i=radii_tmp.begin();
- continue;
- }
- ++i;
- }
- radii.resize(radii_tmp.size());
- sbps.resize(sbps_tmp.size());
- sbpe.resize(sbpe_tmp.size());
- copy(radii_tmp.begin(),radii_tmp.end(),radii.begin());
- copy(sbps_tmp.begin(),sbps_tmp.end(),sbps.begin());
- copy(sbpe_tmp.begin(),sbpe_tmp.end(),sbpe.begin());
-
- //read cooling function data
- spline_func_obj cf;
- for(ifstream ifs(cfg.cfunc_profile.c_str());;)
- {
- assert(ifs.is_open());
- double x,y;
- ifs>>x>>y;
- if(!ifs.good())
- {
- break;
- }
- cerr<<x<<"\t"<<y<<endl;
- if(x>radii.back())
- {
- break;
- }
- cf.add_point(x,y);
- }
- cf.gen_spline();
-
- //read temperature profile data
- spline_func_obj Tprof;
- int tcnt=0;
- for(ifstream ifs1(cfg.tprofile.c_str());;++tcnt)
- {
- assert(ifs1.is_open());
- double x,y;
- ifs1>>x>>y;
- if(!ifs1.good())
- {
- break;
- }
- cerr<<x<<"\t"<<y<<endl;
-#if 0
- if(tcnt==0)
- {
- Tprof.add_point(0,y);
- }
-#endif
- Tprof.add_point(x,y);
- }
-
-
- Tprof.gen_spline();
-
- 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 fitter
- fitter<vector<double>,vector<double>,vector<double>,double> f;
- f.load_data(ds);
- //initial the object, which is used to calculate projection effect
- projector<double> a;
- bool tie_beta=false;
- if(cfg.param_map.find("beta")!=cfg.param_map.end()
- &&cfg.param_map.find("beta1")==cfg.param_map.end()
- &&cfg.param_map.find("beta2")==cfg.param_map.end())
- {
- dbeta2<double> dbetao;
- a.attach_model(dbetao);
- tie_beta=true;
- }
- else if((cfg.param_map.find("beta1")!=cfg.param_map.end()
- ||cfg.param_map.find("beta2")!=cfg.param_map.end())
- &&cfg.param_map.find("beta")==cfg.param_map.end())
- {
- dbeta<double> dbetao;
- a.attach_model(dbetao);
- tie_beta=false;
- }
- else
- {
- cerr<<"Error, cannot decide whether to tie beta together or let them vary freely!"<<endl;
- assert(0);
- }
-
- //attach the cooling function
- a.attach_cfunc(cf);
- a.set_cm_per_pixel(cm_per_pixel);
-
- f.set_model(a);
- //chi^2 statistic
- vchisq<double> c;
- c.verbose(true);
- c.set_limit();
- f.set_statistic(c);
- //optimization method
- f.set_opt_method(powell_method<double,std::vector<double> >());
- //initialize the initial values
- double n01=0;
- double rc1=0;
- double n02=0;
- double rc2=0;
- double beta=0;
- double bkg_level=0;
- if(tie_beta)
- {
- f.set_param_value("beta",.7);
- f.set_param_lower_limit("beta",.3);
- f.set_param_upper_limit("beta",1.4);
- }
- else
- {
- f.set_param_value("beta1",.7);
- f.set_param_lower_limit("beta1",.3);
- f.set_param_upper_limit("beta1",1.4);
- f.set_param_value("beta2",.7);
- f.set_param_lower_limit("beta2",.3);
- f.set_param_upper_limit("beta2",1.4);
- }
- for(std::map<std::string,std::vector<double> >::iterator i=cfg.param_map.begin();
- i!=cfg.param_map.end();++i)
- {
- std::string pname=i->first;
- f.set_param_value(pname,i->second.at(0));
- if(i->second.size()==3)
- {
- double a1=i->second[1];
- double a2=i->second[2];
- double u=std::max(a1,a2);
- double l=std::min(a1,a2);
- f.set_param_upper_limit(pname,u);
- f.set_param_lower_limit(pname,l);
- }
- else
- {
- if(pname=="beta"||pname=="beta1"||pname=="beta2")
- {
- f.set_param_lower_limit(pname,.3);
- f.set_param_upper_limit(pname,1.4);
- }
- }
- }
-
-
-
- //perform the fitting, first freeze beta1, beta2, rc1, and rc2
- if(tie_beta)
- {
- f.set_param_modifier(freeze_param<std::vector<double>,std::vector<double>,std::vector<double>,std::string>("beta")+
- freeze_param<std::vector<double>,std::vector<double>,std::vector<double>,std::string>("rc1")+
- freeze_param<std::vector<double>,std::vector<double>,std::vector<double>,std::string>("rc2")
- );
- }
- else
- {
- f.set_param_modifier(freeze_param<std::vector<double>,std::vector<double>,std::vector<double>,std::string>("beta1")+
- freeze_param<std::vector<double>,std::vector<double>,std::vector<double>,std::string>("beta2")+
- freeze_param<std::vector<double>,std::vector<double>,std::vector<double>,std::string>("rc1")+
- freeze_param<std::vector<double>,std::vector<double>,std::vector<double>,std::string>("rc2")
- );
- }
-
- f.fit();
-
- f.clear_param_modifier();
-
- //then perform the fitting, freeze beta1 and beta2
- //f.set_param_modifier(freeze_param<std::vector<double>,std::vector<double>,std::vector<double>,std::string>("beta"));
- //f.set_param_modifier(freeze_param<std::vector<double>,std::vector<double>,std::vector<double>,std::string>("bkg"));
- f.fit();
- //f.clear_param_modifier();
-
- //finally thaw all parameters
- f.fit();
- double beta1=0;
- double beta2=0;
-
- n01=f.get_param_value("n01");
- rc1=f.get_param_value("rc1");
- n02=f.get_param_value("n02");
- rc2=f.get_param_value("rc2");
- if(tie_beta)
- {
- beta=f.get_param_value("beta");
- beta1=beta;
- beta2=beta;
- }
- else
- {
- beta1=f.get_param_value("beta1");
- beta2=f.get_param_value("beta2");
- }
- //output the params
- ofstream param_output("dbeta_param.txt");
- //output the datasets and fitting results
- for(size_t i=0;i<f.get_num_params();++i)
- {
- if(f.get_param_info(i).get_name()=="rc1")
- {
- cerr<<"rc1_kpc"<<"\t"<<abs(f.get_param_info(i).get_value())*cm_per_pixel/kpc<<endl;
- param_output<<"rc1_kpc"<<"\t"<<abs(f.get_param_info(i).get_value())*cm_per_pixel/kpc<<endl;
- }
- if(f.get_param_info(i).get_name()=="rc2")
- {
- cerr<<"rc2_kpc"<<"\t"<<abs(f.get_param_info(i).get_value())*cm_per_pixel/kpc<<endl;
- param_output<<"rc2_kpc"<<"\t"<<abs(f.get_param_info(i).get_value())*cm_per_pixel/kpc<<endl;
- }
- cerr<<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;
- }
- cerr<<"reduced_chi^2="<<f.get_statistic_value()/(radii.size()-f.get_model().get_num_free_params())<<endl;
- param_output<<"reduced_chi^2="<<f.get_statistic_value()/(radii.size()-f.get_model().get_num_free_params())<<endl;
-
- //c.verbose(false);
- //f.set_statistic(c);
- //f.fit();
- std::vector<double> p=f.get_all_params();
- f.clear_param_modifier();
- std::vector<double> mv=f.eval_model(radii,p);
-
-
- ofstream ofs_sbp("sbp_fit.qdp");
- ofs_sbp<<"read serr 2"<<endl;
- ofs_sbp<<"skip single"<<endl;
-
- ofs_sbp<<"line on 2"<<endl;
- ofs_sbp<<"line on 3"<<endl;
- ofs_sbp<<"line on 4"<<endl;
- ofs_sbp<<"line on 5"<<endl;
- ofs_sbp<<"line on 7"<<endl;
- ofs_sbp<<"ls 2 on 7"<<endl;
- ofs_sbp<<"ls 2 on 3"<<endl;
- ofs_sbp<<"ls 2 on 4"<<endl;
- ofs_sbp<<"ls 2 on 5"<<endl;
-
-
-
- ofs_sbp<<"!LAB POS Y 4.00"<<endl;
- ofs_sbp<<"!LAB ROT"<<endl;
- ofs_sbp<<"win 1"<<endl;
- ofs_sbp<<"yplot 1 2 3 4 5"<<endl;
- ofs_sbp<<"loc 0 0 1 1"<<endl;
- ofs_sbp<<"vie .1 .4 .9 .9"<<endl;
- ofs_sbp<<"la y cnt/s/pixel/cm^2"<<endl;
- ofs_sbp<<"log x"<<endl;
- ofs_sbp<<"log y"<<endl;
- ofs_sbp<<"r x "<<(radii[1]+radii[0])/2*cm_per_pixel/kpc<<" "<<(radii[sbps.size()-2]+radii[sbps.size()-1])/2*cm_per_pixel/kpc<<endl;
- ofs_sbp<<"win 2"<<endl;
- ofs_sbp<<"yplot 6 7"<<endl;
- ofs_sbp<<"loc 0 0 1 1"<<endl;
- ofs_sbp<<"vie .1 .1 .9 .4"<<endl;
- ofs_sbp<<"la x radius (kpc)"<<endl;
- ofs_sbp<<"la y chi"<<endl;
- ofs_sbp<<"log x"<<endl;
- ofs_sbp<<"log y off"<<endl;
- ofs_sbp<<"r x "<<(radii[1]+radii[0])/2*cm_per_pixel/kpc<<" "<<(radii[sbps.size()-2]+radii[sbps.size()-1])/2*cm_per_pixel/kpc<<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<<endl;
- }
- ofs_sbp<<"no no no"<<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"<<ym<<"\t"<<0<<endl;
- }
- //bkg
- ofs_sbp<<"no no no"<<endl;
- bkg_level=abs(f.get_param_value("bkg"));
- for(size_t i=0;i<sbps.size();++i)
- {
- double x=(radii[i]+radii[i-1])/2;
- ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<bkg_level<<"\t0"<<endl;
- }
- //rc1
- ofs_sbp<<"no no no"<<endl;
- double rc1_kpc=abs(f.get_param_value("rc1")*cm_per_pixel/kpc);
- double max_sbp=*max_element(sbps.begin(),sbps.end());
- double min_sbp=*min_element(sbps.begin(),sbps.end());
- for(double x=min_sbp;x<=max_sbp;x+=(max_sbp-min_sbp)/100)
- {
- ofs_sbp<<rc1_kpc<<"\t"<<x<<"\t"<<"0"<<endl;
- }
- //rc2
- ofs_sbp<<"no no no"<<endl;
- double rc2_kpc=abs(f.get_param_value("rc2")*cm_per_pixel/kpc);
- for(double x=min_sbp;x<=max_sbp;x+=(max_sbp-min_sbp)/100)
- {
- ofs_sbp<<rc2_kpc<<"\t"<<x<<"\t"<<"0"<<endl;
- }
- //resid
- ofs_sbp<<"no no no"<<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"<<(ym-sbps[i-1])/sbpe[i-1]<<"\t"<<1<<endl;
- }
- //zero level in resid map
- ofs_sbp<<"no no no"<<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"<<0<<"\t"<<0<<endl;
- }
-
-
- mv=f.eval_model_raw(radii,p);
- ofstream ofs_rho("rho_fit.qdp");
- ofstream ofs_rho_data("rho_fit.dat");
- ofstream ofs_entropy("entropy.qdp");
- ofs_rho<<"la x radius (kpc)"<<endl;
- ofs_rho<<"la y density (cm\\u-3\\d)"<<endl;
- /*
- for(int 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;
- }
- */
-
- //double lower,upper;
- double dr=1;
- //calculate the mass profile
- //const double G=6.673E-8;//cm^3 g^-1 s^-2
- // Molecular weight per electron
- // Reference: Ettori et al. 2013, Space Sci. Rev., 177, 119-154; Eq.(9) below
- static const double mu=1.155;
- static const double mp=1.67262158E-24;//g
- static const double M_sun=1.98892E33;//g
- //static const double k=1.38E-16;
-
- ofstream ofs_mass("mass_int.qdp");
- ofstream ofs_mass_dat("mass_int.dat");
- ofstream ofs_overdensity("overdensity.qdp");
- ofstream ofs_gas_mass("gas_mass_int.qdp");
- //ofs_mass<<"la x radius (kpc)"<<endl;
- //ofs_mass<<"la y mass enclosed (solar mass)"<<endl;
- //ofs_overdensity<<"la x radius (kpc)"<<endl;
- //ofs_overdensity<<"la y overdensity"<<endl;
- double gas_mass=0;
- for(double r=1;r<200000;r+=dr)
- {
- dr=r/100;
- double r1=r+dr;
- double r_cm=r*cm_per_pixel;
- double r1_cm=r1*cm_per_pixel;
- double dr_cm=dr*cm_per_pixel;
- double V_cm3=4./3.*pi*(dr_cm*(r1_cm*r1_cm+r_cm*r_cm+r_cm*r1_cm));
- double ne=dbeta_func(r,n01,rc1,beta1,
- n02,rc2,beta2);//cm^3
-
- double dmgas=V_cm3*ne*mu*mp/M_sun;
- gas_mass+=dmgas;
-
- ofs_gas_mass<<r*cm_per_pixel/kpc<<"\t"<<gas_mass<<endl;
- double ne_beta1=dbeta_func(r,n01,rc1,beta1, 0,rc2,beta2);
-
- double ne_beta2=dbeta_func(r,0,rc1,beta1, n02,rc2,beta2);
-
- double ne1=dbeta_func(r1,n01,rc1,beta1, n02,rc2,beta2);//cm^3
-
- double T_keV=Tprof(r);
- double T1_keV=Tprof(r1);
-
- //double T_K=T_keV*11604505.9;
- //double T1_K=T1_keV*11604505.9;
-
- double dlnT=log(T1_keV/T_keV);
- double dlnr=log(r+dr)-log(r);
- double dlnn=log(ne1/ne);
-
- //double r_kpc=r_cm/kpc;
- double r_Mpc=r_cm/Mpc;
- //double M=-r_cm*T_K*k/G/mu/mp*(dlnT/dlnr+dlnn/dlnr);
- //ref:http://adsabs.harvard.edu/abs/2012MNRAS.422.3503W
- //Walker et al. 2012
- double M=-3.68E13*M_sun*T_keV*r_Mpc*(dlnT/dlnr+dlnn/dlnr);
- double rho=M/(4./3.*pi*r_cm*r_cm*r_cm);
-
- double S=T_keV/pow(ne,2./3.);
- //cout<<r<<"\t"<<M/M_sun<<endl;
- //cout<<r<<"\t"<<T_keV<<endl;
-
- ofs_rho<<r*cm_per_pixel/kpc<<"\t"<<ne<<"\t"<<ne_beta1<<"\t"<<ne_beta2<<endl;
- ofs_rho_data<<r*cm_per_pixel/kpc<<"\t"<<ne<<endl;
- ofs_entropy<<r*cm_per_pixel/kpc<<"\t"<<S<<endl;
-#if 0
- if(r*cm_per_pixel/kpc<5)
- {
- continue;
- }
-#endif
- ofs_mass<<r*cm_per_pixel/kpc<<"\t"<<M/M_sun<<endl;
- if(r<radii.back())
- {
- ofs_mass_dat<<r*cm_per_pixel/kpc<<"\t0\t"<<M/M_sun<<"\t"<<M/M_sun*.1<<endl;
- }
- ofs_overdensity<<r*cm_per_pixel/kpc<<"\t"<<rho/calc_critical_density(z)<<endl;
-
- }
-}
diff --git a/mass_profile/fit_nfw_mass.cpp b/mass_profile/fit_nfw_mass.cpp
deleted file mode 100644
index fae74ef..0000000
--- a/mass_profile/fit_nfw_mass.cpp
+++ /dev/null
@@ -1,159 +0,0 @@
-/*
- Fitting nfw mass profile model
- Author: Junhua Gu
- Last modification 20120721
-*/
-
-#include "nfw.hpp"
-#include <core/optimizer.hpp>
-#include <core/fitter.hpp>
-#include <data_sets/default_data_set.hpp>
-#include "chisq.hpp"
-#include <methods/powell/powell_method.hpp>
-#include <iostream>
-#include <fstream>
-#include <vector>
-#include <string>
-
-using namespace opt_utilities;
-using namespace std;
-const double cm=1;
-const double kpc=3.08568e+21*cm;
-const double pi=4*atan(1);
-static double calc_critical_density(double z,
- const double H0=2.3E-18,
- const double Omega_m=.27)
-{
- const double G=6.673E-8;//cm^3 g^-1 s^2
- 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;
-}
-
-
-int main(int argc,char* argv[])
-{
- if(argc<3)
- {
- cerr<<"Usage:"<<argv[0]<<" <data file with 4 columns of x, xe, y, ye> <z> [rmin in kpc]"<<endl;
- return -1;
- }
- double rmin_kpc=1;
- if(argc>=4)
- {
- rmin_kpc=atof(argv[3]);
- }
- double z=0;
- z=atof(argv[2]);
- //define the fitter
- fitter<double,double,vector<double>,double,std::string> fit;
- //define the data set
- default_data_set<double,double> ds;
- //open the data file
- ifstream ifs(argv[1]);
- //cout<<"read serr 2"<<endl;
- ofstream ofs_fit_result("nfw_fit_result.qdp");
-
- ofs_fit_result<<"read serr 1 2"<<endl;
- ofs_fit_result<<"skip single"<<endl;
- ofs_fit_result<<"line off"<<endl;
- ofs_fit_result<<"li on 2"<<endl;
- ofs_fit_result<<"li on 4"<<endl;
- ofs_fit_result<<"ls 2 on 4"<<endl;
-
- ofs_fit_result<<"win 1"<<endl;
- ofs_fit_result<<"yplot 1 2"<<endl;
- ofs_fit_result<<"loc 0 0 1 1"<<endl;
- ofs_fit_result<<"vie .1 .4 .9 .9"<<endl;
- ofs_fit_result<<"la y Mass (solar)"<<endl;
- ofs_fit_result<<"log x"<<endl;
- ofs_fit_result<<"log y"<<endl;
- ofs_fit_result<<"win 2"<<endl;
- ofs_fit_result<<"yplot 3 4"<<endl;
- ofs_fit_result<<"loc 0 0 1 1"<<endl;
- ofs_fit_result<<"vie .1 .1 .9 .4"<<endl;
- ofs_fit_result<<"la x radius (kpc)"<<endl;
- ofs_fit_result<<"la y chi"<<endl;
- ofs_fit_result<<"log x"<<endl;
- ofs_fit_result<<"log y off"<<endl;
- for(;;)
- {
- //read radius, temperature and error
- double r,re,m,me;
- ifs>>r>>re>>m>>me;
- if(!ifs.good())
- {
- break;
- }
- if(r<rmin_kpc)
- {
- continue;
- }
- data<double,double> d(r,m,me,me,re,re);
- ofs_fit_result<<r<<"\t"<<re<<"\t"<<m<<"\t"<<me<<endl;
- ds.add_data(d);
- }
- ofs_fit_result<<"no no no"<<endl;
- //load data
- fit.load_data(ds);
- //define the optimization method
- fit.set_opt_method(powell_method<double,vector<double> >());
- //use chi^2 statistic
- fit.set_statistic(chisq<double,double,vector<double>,double,std::string>());
- fit.set_model(nfw<double>());
- //fit.set_param_value("rs",4);
- //fit.set_param_value("rho0",100);
- fit.fit();
- fit.fit();
- vector<double> p=fit.fit();
- //output parameters
- ofstream ofs_param("nfw_param.txt");
- for(size_t i=0;i<fit.get_num_params();++i)
- {
- cout<<fit.get_param_info(i).get_name()<<"\t"<<fit.get_param_info(i).get_value()<<endl;
- ofs_param<<fit.get_param_info(i).get_name()<<"\t"<<fit.get_param_info(i).get_value()<<endl;
- }
- cout<<"reduced chi^2="<<fit.get_statistic_value()<<endl;
- ofs_param<<"reduced chi^2="<<fit.get_statistic_value()<<endl;
- ofstream ofs_model("nfw_dump.qdp");
- ofstream ofs_overdensity("overdensity.qdp");
- //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 M_sun=1.98892E33;//g
- //static const double k=1.38E-16;
-
- double xmax=0;
- for(double x=std::max(rmin_kpc,ds.get_data(0).get_x());;x+=1)
- {
- double model_value=fit.eval_model(x,p);
- ofs_model<<x<<"\t"<<model_value<<endl;
- ofs_fit_result<<x<<"\t0\t"<<model_value<<"\t0"<<endl;
- double V=4./3.*pi*pow(x*kpc,3);
- double m=model_value*M_sun;
- double rho=m/V;//g/cm^3
- double over_density=rho/calc_critical_density(z);
- ofs_overdensity<<x<<"\t"<<over_density<<endl;
- xmax=x;
- if(over_density<100)
- {
- break;
- }
- }
- ofs_fit_result<<"no no no"<<endl;
- for(size_t i=0;i<ds.size();++i)
- {
- data<double,double> d=ds.get_data(i);
- double x=d.get_x();
- double y=d.get_y();
- double ye=d.get_y_lower_err();
- double ym=fit.eval_model(x,p);
- ofs_fit_result<<x<<"\t"<<0<<"\t"<<(y-ym)/ye<<"\t"<<1<<endl;
- }
- ofs_fit_result<<"no no no"<<endl;
- for(double x=std::max(rmin_kpc,ds.get_data(0).get_x());x<xmax;x+=1)
- {
- ofs_fit_result<<x<<"\t"<<0<<"\t"<<0<<"\t"<<0<<endl;
- }
-
-}
diff --git a/mass_profile/fit_wang2012_model.cpp b/mass_profile/fit_wang2012_model.cpp
deleted file mode 100644
index f15a2b1..0000000
--- a/mass_profile/fit_wang2012_model.cpp
+++ /dev/null
@@ -1,195 +0,0 @@
-/*
- Fitting Jy Wang's temperature profile model
- Author: Jingying Wang
- Last modification 20120819
-
-*/
-
-#include "wang2012_model.hpp"
-#include <core/optimizer.hpp>
-#include <core/fitter.hpp>
-#include <data_sets/default_data_set.hpp>
-#include "chisq.hpp"
-#include <methods/powell/powell_method.hpp>
-#include <core/freeze_param.hpp>
-#include <iostream>
-#include <fstream>
-#include <vector>
-#include <string>
-
-using namespace opt_utilities;
-using namespace std;
-const double cm=1;
-const double kpc=3.08568e+21*cm;
-
-int main(int argc,char* argv[])
-{
- if(argc<2)
- {
- cerr<<"Usage:"<<argv[0]<<" <data file with 4 columns of x, xe, y, ye> [param file] [cm per pixel]"<<endl;
- return -1;
- }
- double cm_per_pixel=-1;
- if(argc>=4)
- {
- cm_per_pixel=atof(argv[3]);
- }
-
- //define the fitter
- fitter<double,double,vector<double>,double,std::string> fit;
- //define the data set
- default_data_set<double,double> ds;
- //open the data file
- ifstream ifs(argv[1]);
- double min_r=1e9;
- //cout<<"read serr 2"<<endl;
- ofstream ofs_fit_result("fit_result.qdp");
- ofs_fit_result<<"read serr 1 2"<<endl;
- ofs_fit_result<<"skip single"<<endl;
- if(cm_per_pixel>0)
- {
- ofs_fit_result<<"la x radius (kpc)"<<endl;
- }
- else
- {
- ofs_fit_result<<"la x radius (pixel)"<<endl;
- }
- ofs_fit_result<<"la y temperature (keV)"<<endl;
- for(;;)
- {
- //read radius, temperature and error
- double r,re,t,te;
- ifs>>r>>re>>t>>te;
- if(!ifs.good())
- {
- break;
- }
- min_r=min(r,min_r);
- data<double,double> d(r,t,te,te,re,re);
- //std::cerr<<r<<"\t"<<t<<"\t"<<te<<endl;
- if(cm_per_pixel>0)
- {
- ofs_fit_result<<r*cm_per_pixel/kpc<<"\t"<<re*cm_per_pixel/kpc<<"\t"<<t<<"\t"<<te<<endl;
- }
- else
- {
- ofs_fit_result<<r<<"\t"<<re<<"\t"<<t<<"\t"<<te<<endl;
- }
- ds.add_data(d);
- }
- ofs_fit_result<<"no no no"<<endl;
- //load data
- fit.load_data(ds);
- //define the optimization method
- fit.set_opt_method(powell_method<double,vector<double> >());
- //use chi^2 statistic
- chisq<double,double,vector<double>,double,std::string> chisq_object;
- chisq_object.set_limit();
- fit.set_statistic(chisq_object);
- //fit.set_statistic(chisq<double,double,vector<double>,double,std::string>());
- fit.set_model(wang2012_model<double>());
-
- if(argc>=3&&std::string(argv[2])!="NONE")
- {
- std::vector<std::string> freeze_list;
- ifstream ifs_param(argv[2]);
- assert(ifs_param.is_open());
- for(;;)
- {
- string pname;
- double pvalue;
- double lower,upper;
- char param_status;
- ifs_param>>pname>>pvalue>>lower>>upper>>param_status;
- if(!ifs_param.good())
- {
- break;
- }
- if(param_status=='F')
- {
- freeze_list.push_back(pname);
- }
- if(pvalue<=lower||pvalue>=upper)
- {
- cerr<<"Invalid initial value, central value not enclosed by the lower and upper boundaries, adjust automatically"<<endl;
- pvalue=std::max(pvalue,lower);
- pvalue=std::min(pvalue,upper);
- }
- fit.set_param_value(pname,pvalue);
- fit.set_param_lower_limit(pname,lower);
- fit.set_param_upper_limit(pname,upper);
- }
- if(!freeze_list.empty())
- {
- freeze_param<double,double,std::vector<double>,std::string> fp(freeze_list[0]);
- fit.set_param_modifier(fp);
- for(size_t i=1;i<freeze_list.size();++i)
- {
- dynamic_cast<freeze_param<double,double,std::vector<double>,std::string>&>(fit.get_param_modifier())+=freeze_param<double,double,std::vector<double>,std::string>(freeze_list[i]);
- }
- }
- }
-
- for(int i=0;i<100;++i)
- {
- fit.fit();
- }
- vector<double> p=fit.fit();
-#if 0
- ofstream output_param;
- if(argc>=3&&std::string(argv[2])!="NONE")
- {
- output_param.open(argv[2]);
- }
- else
- {
- output_param.open("para0.txt");
- }
-#endif
- //output parameters
- for(size_t i=0;i<fit.get_num_params();++i)
- {
- std::string pname=fit.get_param_info(i).get_name();
- std::string pstatus=fit.get_model().report_param_status(pname);
-
- if(pstatus==""||pstatus=="thawed")
- {
- pstatus="T";
- }
- else
- {
- pstatus="F";
- }
- cout<<fit.get_param_info(i).get_name()<<"\t"<<fit.get_param_info(i).get_value()<<"\t"<<fit.get_param_info(i).get_lower_limit()<<"\t"<<fit.get_param_info(i).get_upper_limit()<<"\t"<<pstatus<<endl;
- //if(argc>=3&&std::string(argv[2])!="NONE")
-#if 0
- {
- output_param<<fit.get_param_info(i).get_name()<<"\t"<<fit.get_param_info(i).get_value()<<"\t"<<fit.get_param_info(i).get_lower_limit()<<"\t"<<fit.get_param_info(i).get_upper_limit()<<"\t"<<pstatus<<endl;
- }
-#endif
- }
-
- //cout<<"T0"<<"\t"<<fit.get_param_value("T0")<<endl;
- //cout<<"T1"<<"\t"<<fit.get_param_value("T1")<<endl;
- //cout<<"xt"<<"\t"<<fit.get_param_value("xt")<<endl;
- //cout<<"eta"<<"\t"<<fit.get_param_value("eta")<<endl;
- //dump the data for checking
- ofstream ofs_model("wang2012_dump.qdp");
-
-
- min_r=0;
- for(double x=min_r;x<3000;x+=10)
- {
- double model_value=fit.eval_model_raw(x,p);
-
- ofs_model<<x<<"\t"<<model_value<<endl;
- if(cm_per_pixel>0)
- {
- ofs_fit_result<<x*cm_per_pixel/kpc<<"\t0\t"<<model_value<<"\t0"<<endl;
- }
- else
- {
- ofs_fit_result<<x<<"\t0\t"<<model_value<<"\t0"<<endl;
- }
- }
-}
diff --git a/mass_profile/nfw.hpp b/mass_profile/nfw.hpp
deleted file mode 100644
index 8e3e59e..0000000
--- a/mass_profile/nfw.hpp
+++ /dev/null
@@ -1,56 +0,0 @@
-/**
- \file nfw.hpp
- \brief Jingying Wang's model
- \author Jingying Wang
- */
-
-
-#ifndef NFW
-#define NFW
-#define OPT_HEADER
-#include <core/fitter.hpp>
-#include <cmath>
-
-namespace opt_utilities
-{
- template <typename T>
- class nfw
- :public model<T,T,std::vector<T>,std::string>
- {
- private:
- model<T,T,std::vector<T> >* do_clone()const
- {
- return new nfw<T>(*this);
- }
-
- const char* do_get_type_name()const
- {
- return "1d power law";
- }
- public:
- nfw()
- {
- this->push_param_info(param_info<std::vector<T> >("rho0",1,0,1e99));
- this->push_param_info(param_info<std::vector<T> >("rs",100,0,1e99));
- }
-
- T do_eval(const T& r,const std::vector<T>& param)
- {
- T rho0=std::abs(param[0]);
- T rs=std::abs(param[1]);
- static const T pi=4*std::atan(1);
- return 4*pi*rho0*rs*rs*rs*(std::log((r+rs)/rs)-r/(r+rs));
- }
-
- private:
- std::string do_get_information()const
- {
- return "";
- }
- };
-}
-
-
-
-#endif
-//EOF
diff --git a/mass_profile/projector.hpp b/mass_profile/projector.hpp
deleted file mode 100644
index 4ef4b32..0000000
--- a/mass_profile/projector.hpp
+++ /dev/null
@@ -1,214 +0,0 @@
-#ifndef PROJ_HPP
-#define PROJ_HPP
-/*
- Defining the class that is used to consider the projection effect
- Author: Junhua Gu
- Last modified: 2011.01.01
-*/
-
-
-#include <core/fitter.hpp>
-#include <vector>
-#include <cmath>
-
-static const double pi=4*atan(1);
-// Ratio of the electron density (n_e) to the proton density (n_p)
-// n_gas = n_e + n_p ~= 1.826 n_e => n_e / n_p ~= 1.21
-// Reference: Ettori et al. 2013, Space Sci. Rev., 177, 119-154; Eq.(9) below
-static const double ne_np_ratio = 1.21;
-
-namespace opt_utilities
-{
- //This is used to project a 3-D surface brightness model to 2-D profile
- template <typename T>
- class projector
- :public model<std::vector<T>,std::vector<T>,std::vector<T> >
- {
- private:
- //Points to a 3-D model that is to be projected
- model<std::vector<T>,std::vector<T>,std::vector<T> >* pmodel;
- func_obj<T,T>* pcfunc;
- T cm_per_pixel;
- public:
- //default cstr
- projector()
- :pmodel(NULL_PTR),pcfunc(NULL_PTR),cm_per_pixel(1)
- {}
- //copy cstr
- projector(const projector& rhs)
- :cm_per_pixel(rhs.cm_per_pixel)
- {
- attach_model(*(rhs.pmodel));
- if(rhs.pcfunc)
- {
- pcfunc=rhs.pcfunc->clone();
- }
- else
- {
- pcfunc=NULL_PTR;
- }
- }
- //assign operator
- projector& operator=(const projector& rhs)
- {
- cm_per_pixel=rhs.cm_per_pixel;
- if(pmodel)
- {
- pmodel->destroy();
- }
- if(pcfunc)
- {
- pcfunc->destroy();
- }
- if(rhs.pcfunc)
- {
- pcfunc=rhs.pcfunc->clone();
- }
- if(rhs.pmodel)
- {
- pmodel=rhs.pmodel->clone();
- }
- }
- //destr
- ~projector()
- {
- if(pmodel)
- {
- pmodel->destroy();
- }
- if(pcfunc)
- {
- pcfunc->destroy();
- }
- }
- //used to clone self
- model<std::vector<T>,std::vector<T>,std::vector<T> >*
- do_clone()const
- {
- return new projector(*this);
- }
-
- public:
- void set_cm_per_pixel(const T& x)
- {
- cm_per_pixel=x;
- }
-
- //attach the model that is to be projected
- void attach_model(const model<std::vector<T>,std::vector<T>,std::vector<T> >& m)
- {
- this->clear_param_info();
- for(size_t i=0;i<m.get_num_params();++i)
- {
- this->push_param_info(m.get_param_info(i));
- }
- this -> push_param_info(param_info<std::vector<T>,
- std::string>("bkg",0,0,1E99));
- pmodel=m.clone();
- pmodel->clear_param_modifier();
- }
-
- void attach_cfunc(const func_obj<T,T>& cf)
- {
- if(pcfunc)
- {
- pcfunc->destroy();
- }
- pcfunc=cf.clone();
- }
-
- public:
- //calc the volume
- /*
- This is a sphere that is subtracted by a cycline.
- /| |\
- / | | \
- | | | |
- | | | |
- \ | | /
- \| |/
- */
- T calc_v_ring(T rsph,T rcyc)
- {
- if(rcyc<rsph)
- {
- double a=rsph*rsph-rcyc*rcyc;
- return 4.*pi/3.*std::sqrt(a*a*a);
- }
- return 0;
- }
-
- //calc the No. nsph sphere's projection volume on the No. nrad pie region
- T calc_v(const std::vector<T>& rlist,int nsph,int nrad)
- {
- if(nsph<nrad)
- {
- return 0;
- }
- else if(nsph==nrad)
- {
- return calc_v_ring(rlist[nsph+1], rlist[nrad]);
- }
- else {
- return (calc_v_ring(rlist[nsph+1], rlist[nrad]) -
- calc_v_ring(rlist[nsph], rlist[nrad]) -
- calc_v_ring(rlist[nsph+1], rlist[nrad+1]) +
- calc_v_ring(rlist[nsph], rlist[nrad+1]));
- }
- }
- public:
- bool do_meets_constraint(const std::vector<T>& p)const
- {
- std::vector<T> p1(this->reform_param(p));
- for(size_t i=0;i!=p1.size();++i)
- {
- if(get_element(p1,i)>this->get_param_info(i).get_upper_limit()||
- get_element(p1,i)<this->get_param_info(i).get_lower_limit())
- {
- // std::cerr<<this->get_param_info(i).get_name()<<"\t"<<p1[i]<<std::endl;
- return false;
- }
- }
- std::vector<T> p2(p1.size()-1);
- for(size_t i=0;i<p1.size()-1;++i)
- {
- p2.at(i)=p1[i];
- }
-
- return pmodel->meets_constraint(p2);
- }
- public:
- //Perform the projection
- std::vector<T> do_eval(const std::vector<T>& x,const std::vector<T>& p)
- {
- T bkg=std::abs(p.back());
- //I think following codes are clear enough :).
- std::vector<T> unprojected(pmodel->eval(x,p));
- std::vector<T> projected(unprojected.size());
-
- for(size_t nrad=0; nrad<x.size()-1; ++nrad)
- {
- for(size_t nsph=nrad; nsph<x.size()-1; ++nsph)
- {
- double v = calc_v(x, nsph, nrad) * pow(cm_per_pixel, 3);
- if(pcfunc)
- {
- double cfunc = (*pcfunc)((x[nsph+1] + x[nsph]) / 2.0);
- projected[nrad] += (unprojected[nsph] * unprojected[nsph] *
- cfunc * v / ne_np_ratio);
- }
- else
- {
- projected[nrad] += unprojected[nsph] * unprojected[nsph] * v;
- }
- }
- double area = pi * (x[nrad+1]*x[nrad+1] - x[nrad]*x[nrad]);
- projected[nrad] /= area;
- projected[nrad] += bkg;
- }
- return projected;
- }
- };
-};
-
-#endif
diff --git a/mass_profile/report_error.cpp b/mass_profile/report_error.cpp
deleted file mode 100644
index 23e6a26..0000000
--- a/mass_profile/report_error.cpp
+++ /dev/null
@@ -1,39 +0,0 @@
-#include <iostream>
-
-using namespace std;
-void report_error(const char* message)
-{
- cerr<<"MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM"<<endl;
- cerr<<"MMMMMMMMMMMM MMMMMMMMMMMM"<<endl;
- cerr<<"MMMMMMMMMM MMMMMMMMMM"<<endl;
- cerr<<"MMMMMMMMM MMMMMMMMM"<<endl;
- cerr<<"MMMMMMMM MMMMMMMM"<<endl;
- cerr<<"MMMMMMM MMMMMMMM"<<endl;
- cerr<<"MMMMMMM MMMMMMM"<<endl;
- cerr<<"MMMMMMM MMMMMMM"<<endl;
- cerr<<"MMMMMMM MMM MMM MMMMMMM"<<endl;
- cerr<<"MMMMMMM MMMMM MMMM MMMMMMM"<<endl;
- cerr<<"MMMMMMM MMMMM MMMM MMMMMMM"<<endl;
- cerr<<"MMMMMMMM MMMM M MMMM MMMMMMMM"<<endl;
- cerr<<"MMMMMMMM M MMMMMMM"<<endl;
- cerr<<"MMMMMMMM MMM MMMMMMMM"<<endl;
- cerr<<"MMMMMMMMMMMM MMM MMMMMMMMMMMM"<<endl;
- cerr<<"MMMMMMMMMM MM M MMMMMMMMM"<<endl;
- cerr<<"MMMMMMMMMM M M M M M MMMMMMMMMM"<<endl;
- cerr<<"MMMM MMMMM MMMMMMMMM MMMMM MM"<<endl;
- cerr<<"MMM MMMM M MMMMM M MMMM MM"<<endl;
- cerr<<"MMM MMMM M M M MMMMM MMM"<<endl;
- cerr<<"MMMM MMMM MMM MM"<<endl;
- cerr<<"MMM MMMM MMMM MM"<<endl;
- cerr<<"MMM MMMMMMMM M MMM"<<endl;
- cerr<<"MMMM MMM MMM MMMMMMMM"<<endl;
- cerr<<"MMMMMMMMMMM MM MMMMMMM M"<<endl;
- cerr<<"MMM MMMMMMM MMMMMMMMM M"<<endl;
- cerr<<"MM MMM MM M"<<endl;
- cerr<<"MM MMMM MM"<<endl;
- cerr<<"MMM MMMMMMMMMMMMM M"<<endl;
- cerr<<"MM MMMMMMMMMMMMMMMMMMM M"<<endl;
- cerr<<"MMM MMMMMMMMMMMMMMMMMMMMMM M"<<endl;
- cerr<<"MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM"<<endl;
- cerr<<message<<endl;
-}
diff --git a/mass_profile/report_error.hpp b/mass_profile/report_error.hpp
deleted file mode 100644
index 30607e5..0000000
--- a/mass_profile/report_error.hpp
+++ /dev/null
@@ -1,6 +0,0 @@
-#ifndef REPORT_ERR
-#define REPORT_ERR
-
-void report_error(const char* message);
-
-#endif
diff --git a/mass_profile/spline.hpp b/mass_profile/spline.hpp
deleted file mode 100644
index 1027b61..0000000
--- a/mass_profile/spline.hpp
+++ /dev/null
@@ -1,110 +0,0 @@
-#ifndef SPLINE_HPP
-#define SPLINE_HPP
-
-#include <vector>
-#include <cstdlib>
-#include <cassert>
-#include <cmath>
-#include <limits>
-
-template <typename T>
-class spline
-{
-public:
- std::vector<T> x_list;
- std::vector<T> y_list;
- std::vector<T> y2_list;
-
-public:
- void push_point(T x,T y)
- {
- if(!x_list.empty())
- {
- assert(x>*(x_list.end()-1));
- }
- x_list.push_back(x);
- y_list.push_back(y);
- }
-
- T get_value(T x)
- {
- if(x<=x_list[0])
- {
- return y_list[0];
- }
- if(x>=x_list.back())
- {
- return y_list.back();
- }
- assert(x_list.size()==y2_list.size());
- assert(x>x_list[0]);
- assert(x<x_list.back());
- int n1,n2;
- n1=0;
- n2=x_list.size()-1;
- while((n2-n1)!=1)
- {
- //cerr<<n1<<"\t"<<n2<<endl;
- if(x_list[n1+1]<=x)
- {
- n1++;
- }
- if(x_list[n2-1]>x)
- {
- n2--;
- }
- }
- T h=x_list[n2]-x_list[n1];
- double a=(x_list[n2]-x)/h;
- double b=(x-x_list[n1])/h;
- return a*y_list[n1]+b*y_list[n2]+((a*a*a-a)*y2_list[n1]+
- (b*b*b-b)*y2_list[n2])*(h*h)/6.;
-
- }
-
- void gen_spline(T y2_0,T y2_N)
- {
- int n=x_list.size();
- y2_list.resize(0);
- y2_list.resize(x_list.size());
- std::vector<T> u(x_list.size());
- if(std::abs(y2_0)<std::numeric_limits<T>::epsilon())
- {
- y2_list[0]=0;
- u[0]=0;
- }
- else
- {
- y2_list[0]=-.5;
- u[0]=(3./(x_list[1]-x_list[0]))*((y_list[1]-y_list[0])/(x_list[1]-x_list[0])-y2_0);
- }
- for(int i=1;i<n-1;++i)
- {
- double sig=(x_list[i]-x_list[i-1])/(x_list[i+1]-x_list[i-1]);
- double p=sig*y2_list[i-1]+2.;
- y2_list[i]=(sig-1.)/p;
- u[i]=(y_list[i+1]-y_list[i])/(x_list[i+1]-x_list[i])
- -(y_list[i]-y_list[i-1])/(x_list[i]-x_list[i-1]);
- u[i]=(6.*u[i]/(x_list[i+1]-x_list[i-1])-sig*u[i-1])/p;
- }
- double qn,un;
- if(std::abs(y2_N)<std::numeric_limits<T>::epsilon())
- {
- qn=un=0;
- }
- else
- {
- qn=.5;
- un=(3./(x_list[n-1]-x_list[n-2]))*(y2_N-(y_list[n-1]-y_list[n-2])/(x_list[n-1]-x_list[n-2]));
-
- }
- y2_list[n-1]=(un-qn*u[n-2])/(qn*y2_list[n-2]+1.);
- for(int i=n-2;i>=0;--i)
- {
- y2_list[i]=y2_list[i]*y2_list[i+1]+u[i];
- }
- }
-
-};
-
-#endif /* SPLINE_HPP */
diff --git a/mass_profile/vchisq.hpp b/mass_profile/vchisq.hpp
deleted file mode 100644
index 8cd8481..0000000
--- a/mass_profile/vchisq.hpp
+++ /dev/null
@@ -1,137 +0,0 @@
-/**
- \file vchisq.hpp
- \brief chi-square statistic
- \author Junhua Gu
- */
-
-#ifndef VCHI_SQ_HPP
-#define VCHI_SQ_HPP
-
-#define OPT_HEADER
-
-#include <core/fitter.hpp>
-#include <iostream>
-#include <vector>
-#include <misc/optvec.hpp>
-#include <cmath>
-
-using std::cerr;
-using std::endl;
-
-namespace opt_utilities
-{
- template<typename T>
- class vchisq
- :public statistic<std::vector<T>,std::vector<T>,std::vector<T>,T,std::string>
- {
- private:
- bool verb;
- int n;
- bool limit_bound;
- typedef std::vector<T> Tp;
-
- vchisq<T>* do_clone()const
- {
- return new vchisq<T>(*this);
- }
-
- const char* do_get_type_name()const
- {
- return "chi^2 statistic";
- }
-
- public:
- void verbose(bool v)
- {
- verb=v;
- }
-
- void set_limit()
- {
- limit_bound=true;
- }
-
- void clear_limit()
- {
- limit_bound=false;
- }
- public:
- vchisq()
- :verb(false),limit_bound(false)
- {}
-
- T do_eval(const std::vector<T>& p)
- {
- if(limit_bound)
- {
- if(!this->get_fitter().get_model().meets_constraint(p))
- {
- return 1e99;
- }
- }
- T result(0);
- std::vector<float> vx;
- std::vector<float> vy;
- std::vector<float> vye;
- std::vector<float> my;
- float x1=1e99,x2=-1e99,y1=1e99,y2=-1e99;
- if(verb)
- {
- n++;
- if(n%100==0)
- {
- vx.resize(this->get_data_set().get_data(0).get_y().size());
- vy.resize(this->get_data_set().get_data(0).get_y().size());
- vye.resize(this->get_data_set().get_data(0).get_y().size());
- my.resize(this->get_data_set().get_data(0).get_y().size());
- }
- }
- for(int i=(this->get_data_set()).size()-1;i>=0;--i)
- {
- const std::vector<double> y_model(this->eval_model(this->get_data_set().get_data(i).get_x(),p));
- const std::vector<double>& y=this->get_data_set().get_data(i).get_y();
- const std::vector<double>& ye=this->get_data_set().get_data(i).get_y_lower_err();
- for(size_t j=0;j<y.size();++j)
- {
- double chi=(y_model[j]-y[j])/ye[j];
- result+=chi*chi;
- }
-
- if(verb&&n%100==0)
- {
- for(size_t j=0;j<y.size();++j)
- {
- vx.at(j)=((this->get_data_set().get_data(i).get_x().at(j)+this->get_data_set().get_data(i).get_x().at(j+1))/2.);
- vy.at(j)=(y[j]);
- vye.at(j)=ye[j];
- my.at(j)=(y_model[j]);
- x1=std::min(vx.at(j),x1);
- y1=std::min(vy.at(j),y1);
- x2=std::max(vx.at(j),x2);
- y2=std::max(vy.at(j),y2);
- vye[j]=log10(vy[j]+vye[j])-log10(vy[j]);
- vx[j]=log10(vx[j]);
- vy[j]=log10(vy[j]);
- my[j]=log10(my[j]);
- }
- }
- }
- if(verb)
- {
- if(n%100==0)
- {
- cerr<<result<<"\t";
- for(size_t i=0;i<get_size(p);++i)
- {
- cerr<<get_element(p,i)<<",";
- }
- cerr<<endl;
- }
- }
-
- return result;
- }
- };
-
-}
-#endif
diff --git a/mass_profile/wang2012_model.hpp b/mass_profile/wang2012_model.hpp
deleted file mode 100644
index c678970..0000000
--- a/mass_profile/wang2012_model.hpp
+++ /dev/null
@@ -1,67 +0,0 @@
-/**
- \file wang2012_model.hpp
- \brief Jingying Wang's model
- \author Jingying Wang
- */
-
-
-#ifndef WANG2012_MODEL
-#define WANG2012_MODEL
-#define OPT_HEADER
-#include <core/fitter.hpp>
-#include <cmath>
-
-namespace opt_utilities
-{
- template <typename T>
- class wang2012_model
- :public model<T,T,std::vector<T>,std::string>
- {
- private:
- model<T,T,std::vector<T> >* do_clone()const
- {
- return new wang2012_model<T>(*this);
- }
-
- const char* do_get_type_name()const
- {
- return "1d power law";
- }
- public:
- wang2012_model()
- {
- this->push_param_info(param_info<std::vector<T> >("A",5,0,500));
- this->push_param_info(param_info<std::vector<T> >("n",1.66,0,10));
- this->push_param_info(param_info<std::vector<T> >("xi",0.45,0,1));
- this->push_param_info(param_info<std::vector<T> >("a2",1500,0,1e8));
- this->push_param_info(param_info<std::vector<T> >("a3",50,0,1e8));
- this->push_param_info(param_info<std::vector<T> >("beta",0.49,0.1,0.7));
- this->push_param_info(param_info<std::vector<T> >("T0",0,0,10));
-
- }
-
- T do_eval(const T& x,const std::vector<T>& param)
- {
- T A=param[0];
- T n=param[1];
- T xi=param[2];
- T a2=param[3];
- T a3=param[4];
- T beta=param[5];
- T T0=param[6];
- return A*(pow(x,n)+xi*a2)/(pow(x,n)+a2)/pow(1+x*x/a3/a3,beta)+T0;
- //return A*(pow(x,n)+a1)/(pow(x,n)+1)/pow(1+x*x/a3/a3,beta)+T0;
- }
-
- private:
- std::string do_get_information()const
- {
- return "";
- }
- };
-}
-
-
-
-#endif
-//EOF
Copyright © 2017-2018 Aaron LI. All Rights Reserved.