aboutsummaryrefslogtreecommitdiffstats
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/Makefile98
-rw-r--r--src/README.md37
-rw-r--r--src/beta.hpp45
-rw-r--r--src/beta_cfg.cpp85
-rw-r--r--src/beta_cfg.hpp23
-rw-r--r--src/calc_lx_beta.cpp497
-rw-r--r--src/calc_lx_dbeta.cpp548
-rw-r--r--src/chisq.hpp206
-rw-r--r--src/dbeta.hpp108
-rw-r--r--src/dump_fit_qdp.cpp55
-rw-r--r--src/dump_fit_qdp.hpp16
-rw-r--r--src/fit_beta_sbp.cpp534
-rw-r--r--src/fit_dbeta_sbp.cpp586
-rw-r--r--src/fit_nfw_mass.cpp159
-rw-r--r--src/fit_wang2012_model.cpp195
-rw-r--r--src/nfw.hpp56
-rw-r--r--src/projector.hpp214
-rw-r--r--src/report_error.cpp39
-rw-r--r--src/report_error.hpp6
-rw-r--r--src/spline.hpp110
-rw-r--r--src/vchisq.hpp137
-rw-r--r--src/wang2012_model.hpp67
22 files changed, 3821 insertions, 0 deletions
diff --git a/src/Makefile b/src/Makefile
new file mode 100644
index 0000000..76ddeaa
--- /dev/null
+++ b/src/Makefile
@@ -0,0 +1,98 @@
+#
+# 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)
+
+
+install: $(TARGETS)
+ @for f in $(TARGETS); do \
+ (cd ../bin && ln -svf ../src/$$f . ); \
+ done
+
+
+uninstall:
+ @for f in $(TARGETS); do \
+ rm -fv ../bin/$$f; \
+ done
diff --git a/src/README.md b/src/README.md
new file mode 100644
index 0000000..0c61821
--- /dev/null
+++ b/src/README.md
@@ -0,0 +1,37 @@
+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/src/beta.hpp b/src/beta.hpp
new file mode 100644
index 0000000..6ae70ac
--- /dev/null
+++ b/src/beta.hpp
@@ -0,0 +1,45 @@
+#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/src/beta_cfg.cpp b/src/beta_cfg.cpp
new file mode 100644
index 0000000..c373688
--- /dev/null
+++ b/src/beta_cfg.cpp
@@ -0,0 +1,85 @@
+#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/src/beta_cfg.hpp b/src/beta_cfg.hpp
new file mode 100644
index 0000000..27148df
--- /dev/null
+++ b/src/beta_cfg.hpp
@@ -0,0 +1,23 @@
+#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/src/calc_lx_beta.cpp b/src/calc_lx_beta.cpp
new file mode 100644
index 0000000..11d7b04
--- /dev/null
+++ b/src/calc_lx_beta.cpp
@@ -0,0 +1,497 @@
+/**
+ * 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/src/calc_lx_dbeta.cpp b/src/calc_lx_dbeta.cpp
new file mode 100644
index 0000000..31addeb
--- /dev/null
+++ b/src/calc_lx_dbeta.cpp
@@ -0,0 +1,548 @@
+/**
+ * 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/src/chisq.hpp b/src/chisq.hpp
new file mode 100644
index 0000000..0e58200
--- /dev/null
+++ b/src/chisq.hpp
@@ -0,0 +1,206 @@
+/**
+ \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/src/dbeta.hpp b/src/dbeta.hpp
new file mode 100644
index 0000000..7246b44
--- /dev/null
+++ b/src/dbeta.hpp
@@ -0,0 +1,108 @@
+#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/src/dump_fit_qdp.cpp b/src/dump_fit_qdp.cpp
new file mode 100644
index 0000000..85307d1
--- /dev/null
+++ b/src/dump_fit_qdp.cpp
@@ -0,0 +1,55 @@
+#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/src/dump_fit_qdp.hpp b/src/dump_fit_qdp.hpp
new file mode 100644
index 0000000..8364dfd
--- /dev/null
+++ b/src/dump_fit_qdp.hpp
@@ -0,0 +1,16 @@
+#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/src/fit_beta_sbp.cpp b/src/fit_beta_sbp.cpp
new file mode 100644
index 0000000..295fa1e
--- /dev/null
+++ b/src/fit_beta_sbp.cpp
@@ -0,0 +1,534 @@
+/*
+ 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/src/fit_dbeta_sbp.cpp b/src/fit_dbeta_sbp.cpp
new file mode 100644
index 0000000..71b3089
--- /dev/null
+++ b/src/fit_dbeta_sbp.cpp
@@ -0,0 +1,586 @@
+/*
+ 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/src/fit_nfw_mass.cpp b/src/fit_nfw_mass.cpp
new file mode 100644
index 0000000..fae74ef
--- /dev/null
+++ b/src/fit_nfw_mass.cpp
@@ -0,0 +1,159 @@
+/*
+ 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/src/fit_wang2012_model.cpp b/src/fit_wang2012_model.cpp
new file mode 100644
index 0000000..f15a2b1
--- /dev/null
+++ b/src/fit_wang2012_model.cpp
@@ -0,0 +1,195 @@
+/*
+ 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/src/nfw.hpp b/src/nfw.hpp
new file mode 100644
index 0000000..8e3e59e
--- /dev/null
+++ b/src/nfw.hpp
@@ -0,0 +1,56 @@
+/**
+ \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/src/projector.hpp b/src/projector.hpp
new file mode 100644
index 0000000..4ef4b32
--- /dev/null
+++ b/src/projector.hpp
@@ -0,0 +1,214 @@
+#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/src/report_error.cpp b/src/report_error.cpp
new file mode 100644
index 0000000..23e6a26
--- /dev/null
+++ b/src/report_error.cpp
@@ -0,0 +1,39 @@
+#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/src/report_error.hpp b/src/report_error.hpp
new file mode 100644
index 0000000..30607e5
--- /dev/null
+++ b/src/report_error.hpp
@@ -0,0 +1,6 @@
+#ifndef REPORT_ERR
+#define REPORT_ERR
+
+void report_error(const char* message);
+
+#endif
diff --git a/src/spline.hpp b/src/spline.hpp
new file mode 100644
index 0000000..1027b61
--- /dev/null
+++ b/src/spline.hpp
@@ -0,0 +1,110 @@
+#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/src/vchisq.hpp b/src/vchisq.hpp
new file mode 100644
index 0000000..8cd8481
--- /dev/null
+++ b/src/vchisq.hpp
@@ -0,0 +1,137 @@
+/**
+ \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/src/wang2012_model.hpp b/src/wang2012_model.hpp
new file mode 100644
index 0000000..c678970
--- /dev/null
+++ b/src/wang2012_model.hpp
@@ -0,0 +1,67 @@
+/**
+ \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.