Rename mass_profile to src; Add install & uninstall to Makefile

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