Rename mass_profile to src; Add install & uninstall to Makefile

22 files changed, 0 insertions, 3809 deletions

diff --git a/mass_profile/Makefile b/mass_profile/Makefile
deleted file mode 100644
index 17b7a0e..0000000
--- a/mass_profile/Makefile
+++ /dev/null
@@ -1,86 +0,0 @@
-#
-# Makefile for `chandra-acis-analysis/mass_profile` tools
-#
-# Junhua GU
-# Weitian LI
-# 2016-06-07
-#
-
-
-CXX ?= g++
-CXXFLAGS += -Wall
-CXXFLAGS += -Werror
-
-ifdef OPENMP
-	CXXFLAGS += -fopenmp
-endif
-
-ifdef DEBUG
-	CXXFLAGS += -g
-else
-	CXXFLAGS += -O2
-endif
-
-OPT_UTIL_INC ?= -I../opt_utilities
-
-TARGETS= fit_dbeta_sbp fit_beta_sbp fit_wang2012_model \
-		fit_nfw_mass calc_lx_dbeta calc_lx_beta
-HEADERS= projector.hpp spline.hpp vchisq.hpp
-
-all: $(TARGETS)
-
-# NOTE:
-# Object/source files should placed *before* libraries (order matters)
-
-fit_dbeta_sbp: fit_dbeta_sbp.o beta_cfg.o report_error.o
-	$(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC)
-
-fit_beta_sbp: fit_beta_sbp.o beta_cfg.o dump_fit_qdp.o
-	$(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC)
-
-fit_wang2012_model: fit_wang2012_model.o
-	$(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC)
-
-fit_nfw_mass: fit_nfw_mass.o
-	$(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC)
-
-calc_lx_dbeta: calc_lx_dbeta.o beta_cfg.o report_error.o
-	$(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC)
-
-calc_lx_beta: calc_lx_beta.o beta_cfg.o dump_fit_qdp.o
-	$(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC)
-
-
-fit_dbeta_sbp.o: fit_dbeta_sbp.cpp $(HEADERS)
-	$(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC)
-
-fit_beta_sbp.o: fit_beta_sbp.cpp beta.hpp $(HEADERS)
-	$(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC)
-
-fit_wang2012_model.o: fit_wang2012_model.cpp wang2012_model.hpp chisq.hpp
-	$(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC)
-
-fit_nfw_mass.o: fit_nfw_mass.cpp nfw.hpp chisq.hpp
-	$(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC)
-
-calc_lx_dbeta.o: calc_lx_dbeta.cpp $(HEADERS)
-	$(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC)
-
-calc_lx_beta.o: calc_lx_beta.cpp beta.hpp $(HEADERS)
-	$(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC)
-
-beta_cfg.o: beta_cfg.cpp beta_cfg.hpp
-	$(CXX) $(CXXFLAGS) -c $<
-
-report_error.o: report_error.cpp report_error.hpp
-	$(CXX) $(CXXFLAGS) -c $<
-
-dump_fit_qdp.o: dump_fit_qdp.cpp dump_fit_qdp.hpp
-	$(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC)
-
-%.o: %.cpp
-	$(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC)
-
-
-clean:
-	rm -f *.o $(TARGETS)
diff --git a/mass_profile/README.md b/mass_profile/README.md
deleted file mode 100644
index 0c61821..0000000
--- a/mass_profile/README.md
+++ /dev/null
@@ -1,37 +0,0 @@
-Mass/Luminosity/Flux Calculation Tools
-==============================
-
-Junhua GU, Weitian LI, Zhenghao ZHU
-
-
-Introduction
-------------
-This directory contains the tools/scripts used to fit the temperature
-profile and surface brightness file, to calculate the gas density profile
-and gravitational mass profile, to calculate luminosity and flux and
-other related quantities.
-
-
-NOTE
-----
-* Mass calculation references:
-  + Walker et al. 2012, MNRAS, 422, 3503-3515,
-    [ADS:2012MNRAS.422.3503W](http://adsabs.harvard.edu/abs/2012MNRAS.422.3503W)
-  + Ettori et al. 2013, Space Science Reviews, 177, 119-154,
-    [ADS:2013SSRv..177..119E](http://adsabs.harvard.edu/abs/2013SSRv..177..119E)
-* Errors are calculated at 68% confident level.
-* NFW profile is employed to extrapolate the mass profile.
-* We use a self-proposed model (see ``wang2012_model.hpp``) to describe/fit
-  the gas temperature profile.
-* The single-beta and double-beta models (with a constant background) are used
-  to describe the gas density profile.
-
-
-TODO
-----
-* Merge the scripts/tools for single-beta and double-beta SBP models
-  (to reduce code duplications)
-* Uncertainties/errors calculation **maybe** inappropriate/problematic
-  (e.g., ``analyze_mass_profile.py``);
-  why not just use the quantile-based (e.g., Q84.15-Q15.85 ~= 68.3%)
-  uncertainties or standard deviation???
diff --git a/mass_profile/beta.hpp b/mass_profile/beta.hpp
deleted file mode 100644
index 6ae70ac..0000000
--- a/mass_profile/beta.hpp
+++ /dev/null
@@ -1,45 +0,0 @@
-#ifndef BETA
-#define BETA
-#include "projector.hpp"
-
-namespace opt_utilities
-{
-  template <typename T>
-  class beta
-    :public model<std::vector<T>,std::vector<T>,std::vector<T> >
-  {
-  public:
-    beta()
-    {
-      this->push_param_info(param_info<std::vector<T>,std::string>("n0",1,0,1E99));
-      this->push_param_info(param_info<std::vector<T>,std::string>("beta",.66,0,1E99));
-      this->push_param_info(param_info<std::vector<T>,std::string>("rc",100,0,1E99));
-    }
-
-  public:
-    beta<T>* do_clone()const
-    {
-      return new beta<T>(*this);
-    }
-
-    std::vector<T> do_eval(const std::vector<T> & x,
-			   const std::vector<T>& p)
-    {
-      T n0=std::abs(p[0]);
-      T beta=p[1];
-      T rc=p[2];
-
-      std::vector<T> result(x.size()-1);
-      for(size_t i=1;i<x.size();++i)
-	{
-	  T xi=(x[i]+x[i-1])/2;
-	  T yi=0;
-	  yi=n0*pow(1+xi*xi/rc/rc,-3./2.*beta);
-	  result[i-1]=yi;
-	}
-      return result;
-    }
-  };
-}
-
-#endif
diff --git a/mass_profile/beta_cfg.cpp b/mass_profile/beta_cfg.cpp
deleted file mode 100644
index c373688..0000000
--- a/mass_profile/beta_cfg.cpp
+++ /dev/null
@@ -1,85 +0,0 @@
-#include "beta_cfg.hpp"
-#include <sstream>
-#include <cstdlib>
-using namespace std;
-
-cfg_map parse_cfg_file(std::istream& is)
-{
-  cfg_map result;
-  result.rmin_pixel=-1;
-  result.rmin_kpc=-1;
-  for(;;)
-    {
-      std::string line;
-      getline(is,line);
-      line+="\n";
-      if(!is.good())
-	{
-	  break;
-	}
-      string key;
-      istringstream iss(line);
-      iss>>key;
-      if(key=="sbp_data")
-	{
-	  string value;
-	  iss>>value;
-	  result.sbp_data=value;
-	}
-      else if(key=="cfunc_profile")
-	{
-	  string value;
-	  iss>>value;
-	  result.cfunc_profile=value;
-	}
-      else if(key=="tprofile")
-	{
-	  string value;
-	  iss>>value;
-	  result.tprofile=value;
-	}
-      else if(key=="z")
-	{
-	  double z;
-	  iss>>z;
-	  result.z=z;
-	}
-      else if(key=="cm_per_pixel")
-	{
-	  double cm_per_pixel;
-	  iss>>cm_per_pixel;
-	  result.cm_per_pixel=cm_per_pixel;
-	}
-      else if(key=="rmin_pixel")
-	{
-	  double v;
-	  iss>>v;
-	  result.rmin_pixel=v;
-	}
-      else if(key=="rmin_kpc")
-	{
-	  double v;
-	  iss>>v;
-	  result.rmin_kpc=v;
-	}
-      else
-	{
-	  std::vector<double> value;
-	  for(;;)
-	    {
-	      double v;
-	      iss>>v;
-	      if(!iss.good())
-		{
-		  break;
-		}
-	      value.push_back(v);
-	    }
-	  if(!value.empty())
-	    {
-	      result.param_map[key]=value;
-	    }
-	}
-    }
-  return result;
-}
diff --git a/mass_profile/beta_cfg.hpp b/mass_profile/beta_cfg.hpp
deleted file mode 100644
index 27148df..0000000
--- a/mass_profile/beta_cfg.hpp
+++ /dev/null
@@ -1,23 +0,0 @@
-#ifndef BETA_CFG
-#define BETA_CFG
-
-#include <map>
-#include <vector>
-#include <string>
-#include <iostream>
-
-struct cfg_map
-{
-  std::string sbp_data;
-  std::string cfunc_profile;
-  std::string tprofile;
-  double z;
-  double cm_per_pixel;
-  double rmin_kpc;
-  double rmin_pixel;
-  std::map<std::string,std::vector<double> > param_map;
-};
-
-cfg_map parse_cfg_file(std::istream& is);
-
-#endif
diff --git a/mass_profile/calc_lx_beta.cpp b/mass_profile/calc_lx_beta.cpp
deleted file mode 100644
index 11d7b04..0000000
--- a/mass_profile/calc_lx_beta.cpp
+++ /dev/null
@@ -1,497 +0,0 @@
-/**
- * Calculate the total luminosity and flux within the specified radius.
- *
- * Base on 'fit_beta_sbp.cpp' and supersede 'calc_lx.cpp'
- *
- * Author: Junhua Gu
- */
-
-#include <iostream>
-#include <fstream>
-#include <sstream>
-#include <list>
-#include <algorithm>
-#include "beta_cfg.hpp"
-#include "dump_fit_qdp.hpp"
-#include "report_error.hpp"
-#include "vchisq.hpp"
-#include "chisq.hpp"
-#include "beta.hpp"
-#include "models/beta1d.hpp"
-#include <data_sets/default_data_set.hpp>
-#include <methods/powell/powell_method.hpp>
-#include <core/freeze_param.hpp>
-#include <error_estimator/error_estimator.hpp>
-#include "spline.hpp"
-
-using namespace std;
-using namespace opt_utilities;
-//double s=5.63136645E20;
-const double kpc=3.086E21;//kpc in cm
-const double Mpc=kpc*1000;
-
-double beta_func(double r, double n0, double rc, double beta)
-{
-  return abs(n0) * pow(1+r*r/rc/rc, -3./2.*abs(beta));
-}
-
-//A class enclosing the spline interpolation method
-class spline_func_obj
-  :public func_obj<double,double>
-{
-  //has an spline object
-  spline<double> spl;
-public:
-  //This function is used to calculate the intepolated value
-  double do_eval(const double& x)
-  {
-    return spl.get_value(x);
-  }
-
-  //we need this function, when this object is performing a clone of itself
-  spline_func_obj* do_clone()const
-  {
-    return new spline_func_obj(*this);
-  }
-
-public:
-  //add points to the spline object, after which the spline will be initialized
-  void add_point(double x,double y)
-  {
-    spl.push_point(x,y);
-  }
-
-  //before getting the intepolated value, the spline should be initialzied by calling this function
-  void gen_spline()
-  {
-    spl.gen_spline(0,0);
-  }
-};
-
-int main(int argc,char* argv[])
-{
-  if(argc<4)
-    {
-      cerr<<argv[0]<<" <sbp.conf> <rout_kpc> <cfunc_erg> [cfunc2_erg ...]"<<endl;
-      return -1;
-    }
-  //initialize the parameters list
-  ifstream cfg_file(argv[1]);
-  assert(cfg_file.is_open());
-  cfg_map cfg=parse_cfg_file(cfg_file);
-
-  const double z=cfg.z;
-
-  //initialize the radius list, sbp list and sbp error list
-  std::vector<double> radii;
-  std::vector<double> sbps;
-  std::vector<double> sbpe;
-  std::vector<double> radii_all;
-  std::vector<double> sbps_all;
-  std::vector<double> sbpe_all;
-  //prepend the zero point to radius list
-  radii.push_back(0.0);
-  radii_all.push_back(0.0);
-  //read sbp and sbp error data
-  ifstream ifs(cfg.sbp_data.c_str());
-  std::string line;
-  if (ifs.is_open())
-    {
-      while(std::getline(ifs, line))
-        {
-          if (line.empty())
-            continue;
-
-          std::istringstream iss(line);
-          double x, xe, y, ye;
-          if ((iss >> x >> xe >> y >> ye))
-            {
-              std::cerr << "sbprofile data: "
-                        << x << ", " << xe << ", " << y << ", " << ye
-                        << std::endl;
-              radii.push_back(x+xe);  /* NOTE: use outer radii of regions */
-              radii_all.push_back(x+xe);
-              sbps.push_back(y);
-              sbps_all.push_back(y);
-              sbpe.push_back(ye);
-              sbpe_all.push_back(ye);
-            }
-          else
-            {
-              std::cerr << "skipped line: " << line << std::endl;
-            }
-        }
-    }
-  else
-    {
-      std::cerr << "ERROR: cannot open file: " << cfg.sbp_data.c_str()
-                << std::endl;
-      return 1;
-    }
-
-  //initialize the cm/pixel value
-  double cm_per_pixel=cfg.cm_per_pixel;
-  double rmin;
-  if(cfg.rmin_pixel>0)
-    {
-      rmin=cfg.rmin_pixel;
-    }
-  else
-    {
-      rmin=cfg.rmin_kpc*kpc/cm_per_pixel;
-    }
-
-  cerr<<"rmin="<<rmin<<" (pixel)"<<endl;
-  std::list<double> radii_tmp,sbps_tmp,sbpe_tmp;
-  radii_tmp.resize(radii.size());
-  sbps_tmp.resize(sbps.size());
-  sbpe_tmp.resize(sbpe.size());
-  copy(radii.begin(),radii.end(),radii_tmp.begin());
-  copy(sbps.begin(),sbps.end(),sbps_tmp.begin());
-  copy(sbpe.begin(),sbpe.end(),sbpe_tmp.begin());
-  for(list<double>::iterator i=radii_tmp.begin();i!=radii_tmp.end();)
-    {
-      if(*i<rmin)
-	{
-	  radii_tmp.pop_front();
-	  sbps_tmp.pop_front();
-	  sbpe_tmp.pop_front();
-	  i=radii_tmp.begin();
-	  continue;
-	}
-      ++i;
-    }
-  radii.resize(radii_tmp.size());
-  sbps.resize(sbps_tmp.size());
-  sbpe.resize(sbpe_tmp.size());
-  copy(radii_tmp.begin(),radii_tmp.end(),radii.begin());
-  copy(sbps_tmp.begin(),sbps_tmp.end(),sbps.begin());
-  copy(sbpe_tmp.begin(),sbpe_tmp.end(),sbpe.begin());
-
-  //read cooling function data
-  spline_func_obj cf;
-  for(ifstream ifs(cfg.cfunc_profile.c_str());;)
-    {
-      assert(ifs.is_open());
-      double x,y;
-      ifs>>x>>y;
-      if(!ifs.good())
-	{
-	  break;
-	}
-      cerr<<x<<"\t"<<y<<endl;
-      if(x>radii.back())
-	{
-	  break;
-	}
-      cf.add_point(x,y);
-    }
-  cf.gen_spline();
-
-  //read temperature profile data
-  spline_func_obj Tprof;
-  int tcnt=0;
-  for(ifstream ifs1(cfg.tprofile.c_str());;++tcnt)
-    {
-      assert(ifs1.is_open());
-      double x,y;
-      ifs1>>x>>y;
-      if(!ifs1.good())
-      {
-	break;
-      }
-      cerr<<x<<"\t"<<y<<endl;
-#if 0
-      if(tcnt==0)
-	{
-	  Tprof.add_point(0,y);
-	}
-#endif
-      Tprof.add_point(x,y);
-    }
-
-
-  Tprof.gen_spline();
-
-  default_data_set<std::vector<double>,std::vector<double> > ds;
-  ds.add_data(data<std::vector<double>,std::vector<double> >(radii,sbps,sbpe,sbpe,radii,radii));
-
-  //initial fitter
-  fitter<vector<double>,vector<double>,vector<double>,double> f;
-  f.load_data(ds);
-  //initial the object, which is used to calculate projection effect
-  projector<double> a;
-  beta<double> betao;
-  //attach the cooling function
-  a.attach_cfunc(cf);
-  a.set_cm_per_pixel(cm_per_pixel);
-  a.attach_model(betao);
-  f.set_model(a);
-  //chi^2 statistic
-  vchisq<double> c;
-  c.verbose(true);
-  c.set_limit();
-  f.set_statistic(c);
-  //optimization method
-  f.set_opt_method(powell_method<double,std::vector<double> >());
-  //initialize the initial values
-  /*
-  double n0=0;
-  double beta=0;
-  double rc=0;
-  */
-  double bkg_level=0;
-
-  for(std::map<std::string,std::vector<double> >::iterator i=cfg.param_map.begin();
-      i!=cfg.param_map.end();++i)
-    {
-      std::string pname=i->first;
-      f.set_param_value(pname,i->second.at(0));
-      if(i->second.size()==3)
-	{
-	  double a1=i->second[1];
-	  double a2=i->second[2];
-	  double u=std::max(a1,a2);
-	  double l=std::min(a1,a2);
-	  f.set_param_upper_limit(pname,u);
-	  f.set_param_lower_limit(pname,l);
-	}
-      else
-	{
-	  if(pname=="beta")
-	    {
-	      f.set_param_lower_limit(pname,.3);
-	      f.set_param_upper_limit(pname,1.4);
-	    }
-	}
-    }
-
-  f.fit();
-  f.fit();
-  std::vector<double> p=f.get_all_params();
-  /*
-  n0=f.get_param_value("n0");
-  rc=f.get_param_value("rc");
-  beta=f.get_param_value("beta");
-  */
-  //output the datasets and fitting results
-  ofstream param_output("lx_beta_param.txt");
-  for(size_t i=0;i<f.get_num_params();++i)
-    {
-      if(f.get_param_info(i).get_name()=="rc")
-	{
-	  cerr<<"rc_kpc"<<"\t"<<abs(f.get_param_info(i).get_value())*cm_per_pixel/kpc<<endl;
-	  param_output<<"rc_kpc"<<"\t"<<abs(f.get_param_info(i).get_value())*cm_per_pixel/kpc<<endl;
-	}
-      cerr<<f.get_param_info(i).get_name()<<"\t"<<abs(f.get_param_info(i).get_value())<<endl;
-      param_output<<f.get_param_info(i).get_name()<<"\t"<<abs(f.get_param_info(i).get_value())<<endl;
-    }
-  cerr<<"reduced_chi^2="<<f.get_statistic_value()/(radii.size()-f.get_model().get_num_free_params())<<endl;
-  param_output<<"reduced_chi^2="<<f.get_statistic_value()/(radii.size()-f.get_model().get_num_free_params())<<endl;
-
-  std::vector<double> mv=f.eval_model_raw(radii_all,p);
-  int sbps_inner_cut_size=int(sbps_all.size()-sbps.size());
-  ofstream ofs_sbp("lx_sbp_fit.qdp");
-  ofs_sbp<<"read serr 2"<<endl;
-  ofs_sbp<<"skip single"<<endl;
-  ofs_sbp<<"line off "<<endl;
-  if(sbps_inner_cut_size>=1)
-    {
-      ofs_sbp<<"line on 2"<<endl;
-      ofs_sbp<<"line on 3"<<endl;
-      ofs_sbp<<"line on 4"<<endl;
-      ofs_sbp<<"line on 5"<<endl;
-      ofs_sbp<<"ls 2 on 2"<<endl;
-      ofs_sbp<<"ls 2 on 4"<<endl;
-      ofs_sbp<<"ls 2 on 5"<<endl;
-      ofs_sbp<<"line on 7"<<endl;
-      ofs_sbp<<"ls 2 on 7"<<endl;
-
-      ofs_sbp<<"ma 1 on 2"<<endl;
-      ofs_sbp<<"color 1 on 1"<<endl;
-      ofs_sbp<<"color 2 on 2"<<endl;
-      ofs_sbp<<"color 3 on 3"<<endl;
-      ofs_sbp<<"color 4 on 4"<<endl;
-      ofs_sbp<<"color 5 on 5"<<endl;
-
-      ofs_sbp<<"win 1"<<endl;
-      ofs_sbp<<"yplot 1 2 3 4 5"<<endl;
-      ofs_sbp<<"loc 0 0 1 1"<<endl;
-      ofs_sbp<<"vie .1 .4 .9 .9"<<endl;
-      ofs_sbp<<"la y cnt/s/pixel/cm^2"<<endl;
-      ofs_sbp<<"log x"<<endl;
-      ofs_sbp<<"log y"<<endl;
-      ofs_sbp<<"r x "<<(radii[1]+radii[0])/2*cm_per_pixel/kpc<<" "<<(radii[sbps.size()-2]+radii[sbps.size()-1])/2*cm_per_pixel/kpc<<endl;
-      ofs_sbp<<"win 2"<<endl;
-      ofs_sbp<<"yplot 6 7"<<endl;
-      ofs_sbp<<"loc 0 0 1 1"<<endl;
-      ofs_sbp<<"vie .1 .1 .9 .4"<<endl;
-      ofs_sbp<<"la x radius (kpc)"<<endl;
-      ofs_sbp<<"la y chi"<<endl;
-      ofs_sbp<<"log y off"<<endl;
-      ofs_sbp<<"log x"<<endl;
-      ofs_sbp<<"r x "<<(radii[1]+radii[0])/2*cm_per_pixel/kpc<<" "<<(radii[sbps.size()-2]+radii[sbps.size()-1])/2*cm_per_pixel/kpc<<endl;
-    }
-  else
-    {
-      ofs_sbp<<"line on 2"<<endl;
-      ofs_sbp<<"line on 3"<<endl;
-      ofs_sbp<<"line on 4"<<endl;
-      ofs_sbp<<"ls 2 on 3"<<endl;
-      ofs_sbp<<"ls 2 on 4"<<endl;
-      ofs_sbp<<"line on 6"<<endl;
-      ofs_sbp<<"ls 2 on 6"<<endl;
-
-      ofs_sbp<<"color 1 on 1"<<endl;
-      ofs_sbp<<"color 3 on 2"<<endl;
-      ofs_sbp<<"color 4 on 3"<<endl;
-      ofs_sbp<<"color 5 on 4"<<endl;
-      //ofs_sbp<<"ma 1 on 2"<<endl;
-
-      ofs_sbp<<"win 1"<<endl;
-      ofs_sbp<<"yplot 1 2 3 4"<<endl;
-      ofs_sbp<<"loc 0 0 1 1"<<endl;
-      ofs_sbp<<"vie .1 .4 .9 .9"<<endl;
-      ofs_sbp<<"la y cnt/s/pixel/cm^2"<<endl;
-      ofs_sbp<<"log x"<<endl;
-      ofs_sbp<<"log y"<<endl;
-      ofs_sbp<<"r x "<<(radii[1]+radii[0])/2*cm_per_pixel/kpc<<" "<<(radii[radii.size()-2]+radii[radii.size()-1])/2*cm_per_pixel/kpc<<endl;
-      ofs_sbp<<"win 2"<<endl;
-      ofs_sbp<<"yplot 5 6"<<endl;
-      ofs_sbp<<"loc 0 0 1 1"<<endl;
-      ofs_sbp<<"vie .1 .1 .9 .4"<<endl;
-      ofs_sbp<<"la x radius (kpc)"<<endl;
-      ofs_sbp<<"la y chi"<<endl;
-      ofs_sbp<<"log x"<<endl;
-      ofs_sbp<<"log y off"<<endl;
-      ofs_sbp<<"r x "<<(radii[1]+radii[0])/2*cm_per_pixel/kpc<<" "<<(radii[radii.size()-2]+radii[radii.size()-1])/2*cm_per_pixel/kpc<<endl;
-
-    }
-  // cout<<sbps_all.size()<<"\t"<<sbps.size()<<"\t"<<sbps_inner_cut_size<<endl;
-  for(size_t i=1;i<sbps_all.size();++i)
-    {
-      double x=(radii_all[i]+radii_all[i-1])/2;
-      double y=sbps_all[i-1];
-      double ye=sbpe_all[i-1];
-      ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<y<<"\t"<<ye<<endl;
-    }
-  if(sbps_inner_cut_size>=1)
-    {
-      ofs_sbp<<"no no no"<<endl;
-      for(int i=1;i<sbps_inner_cut_size+1;++i)
-	{
-	  double x=(radii_all[i]+radii_all[i-1])/2;
-	  double ym=mv[i-1];
-	  ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<ym<<"\t"<<"0"<<endl;
-	}
-    }
-  ofs_sbp<<"no no no"<<endl;
-  for(size_t i=sbps_inner_cut_size;i<sbps_all.size();++i)
-    {
-      double x=(radii_all[i]+radii_all[i-1])/2;
-      double ym=mv[i-1];
-      ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<ym<<"\t"<<"0"<<endl;
-    }
-  ofs_sbp<<"no no no"<<endl;
-  //bkg level
-  bkg_level=abs(f.get_param_value("bkg"));
-  for(size_t i=0;i<sbps_all.size();++i)
-    {
-      double x=(radii_all[i]+radii_all[i-1])/2;
-      ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<bkg_level<<"\t0"<<endl;
-    }
-  ofs_sbp<<"no no no"<<endl;
-  //rc
-  double rc_kpc=abs(f.get_param_value("rc")*cm_per_pixel/kpc);
-  double max_sbp=*max_element(sbps_all.begin(),sbps_all.end());
-  double min_sbp=*min_element(sbps_all.begin(),sbps_all.end());
-  for(double x=min_sbp;x<=max_sbp;x+=(max_sbp-min_sbp)/100)
-    {
-      ofs_sbp<<rc_kpc<<"\t"<<x<<"\t"<<"0"<<endl;
-    }
-  //resid
-  ofs_sbp<<"no no no"<<endl;
-  for(size_t i=1;i<sbps.size();++i)
-    {
-      double x=(radii[i]+radii[i-1])/2;
-      //double y=sbps[i-1];
-      //double ye=sbpe[i-1];
-      double ym=mv[i-1];
-      ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<(ym-sbps[i-1])/sbpe[i-1]<<"\t"<<1<<endl;
-    }
-
-  //zero level of resid
-  ofs_sbp<<"no no no"<<endl;
-  for(size_t i=1;i<sbps.size();++i)
-    {
-      double x=(radii[i]+radii[i-1])/2;
-      //double y=sbps[i-1];
-      //double ye=sbpe[i-1];
-      //double ym=mv[i-1];
-      ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<0<<"\t"<<0<<endl;
-    }
-
-  mv=betao.eval(radii,p);
-  ofstream ofs_rho("lx_rho_fit.qdp");
-  ofstream ofs_rho_data("lx_rho_fit.dat");
-
-  ofs_rho<<"la x radius (kpc)"<<endl;
-  ofs_rho<<"la y density (cm\\u-3\\d)"<<endl;
-  /*
-  for(int i=1;i<sbps.size();++i)
-    {
-      double x=(radii[i]+radii[i-1])/2;
-      double ym=mv[i-1];
-      ofs_rho<<x*cm_per_pixel/kpc<<"\t"<<ym<<endl;
-    }
-  */
-  p.back()=0;
-  radii.clear();
-  double rout=atof(argv[2])*kpc;
-
-  for(double r=0;r<rout;r+=1*kpc)//step size=1kpc
-    {
-      double r_pix=r/cm_per_pixel;
-      radii.push_back(r_pix);
-    }
-
-  double Da=cm_per_pixel/(.492/3600./180.*pi);
-  double Dl=Da*(1+z)*(1+z);
-  cout<<"dl="<<Dl/kpc<<endl;
-
-  for(int n=3;n<argc;++n)
-    {
-      spline_func_obj cf_erg;
-      for(ifstream ifs(argv[n]);;)
-	{
-	  assert(ifs.is_open());
-	  double x,y;
-	  ifs>>x>>y;
-	  if(!ifs.good())
-	    {
-	      break;
-	    }
-	  //cerr<<x<<"\t"<<y<<endl;
-
-	  cf_erg.add_point(x,y);//change with source
-	}
-      cf_erg.gen_spline();
-
-      projector<double>& pj=dynamic_cast<projector<double>&>(f.get_model());
-      pj.attach_cfunc(cf_erg);
-
-
-
-      mv=f.eval_model_raw(radii,p);
-      double flux_erg=0;
-      for(size_t i=0;i<radii.size()-1;++i)
-	{
-	  double S=pi*(radii[i+1]+radii[i])*(radii[i+1]-radii[i]);
-	  flux_erg+=S*mv[i];
-	}
-      cout<<flux_erg*4*pi*Dl*Dl<<endl;
-      cout<<flux_erg<<endl;
-      param_output<<"Lx"<<n-2<<"\t"<<flux_erg*4*pi*Dl*Dl<<endl;
-      param_output<<"Fx"<<n-2<<"\t"<<flux_erg<<endl;
-    }
-}
diff --git a/mass_profile/calc_lx_dbeta.cpp b/mass_profile/calc_lx_dbeta.cpp
deleted file mode 100644
index 31addeb..0000000
--- a/mass_profile/calc_lx_dbeta.cpp
+++ /dev/null
@@ -1,548 +0,0 @@
-/**
- * Calculate the total luminosity and flux within the specified radius.
- *
- * Base on 'fit_dbeta_sbp.cpp' and supersede 'calc_lx.cpp'
- *
- * Author: Junhua Gu
- */
-
-#include <iostream>
-#include <fstream>
-#include <sstream>
-#include <list>
-#include "vchisq.hpp"
-#include "dbeta.hpp"
-#include "beta_cfg.hpp"
-#include <data_sets/default_data_set.hpp>
-#include <methods/powell/powell_method.hpp>
-#include <core/freeze_param.hpp>
-#include <error_estimator/error_estimator.hpp>
-#include "spline.hpp"
-
-using namespace std;
-using namespace opt_utilities;
-//double s=5.63136645E20;
-const double kpc=3.086E21;//kpc in cm
-const double Mpc=kpc*1000;
-
-double dbeta_func(double r, double n01, double rc1, double beta1,
-                  double n02, double rc2, double beta2)
-{
-  double v1 = abs(n01) * pow(1+r*r/rc1/rc1, -3./2.*abs(beta1));
-  double v2 = abs(n02) * pow(1+r*r/rc2/rc2, -3./2.*abs(beta2));
-  return v1 + v2;
-}
-
-
-//A class enclosing the spline interpolation method
-class spline_func_obj
-  :public func_obj<double,double>
-{
-  //has an spline object
-  spline<double> spl;
-public:
-  //This function is used to calculate the intepolated value
-  double do_eval(const double& x)
-  {
-    return spl.get_value(x);
-  }
-
-  //we need this function, when this object is performing a clone of itself
-  spline_func_obj* do_clone()const
-  {
-    return new spline_func_obj(*this);
-  }
-
-public:
-  //add points to the spline object, after which the spline will be initialized
-  void add_point(double x,double y)
-  {
-    spl.push_point(x,y);
-  }
-
-  //before getting the intepolated value, the spline should be initialzied by calling this function
-  void gen_spline()
-  {
-    spl.gen_spline(0,0);
-  }
-};
-
-int main(int argc,char* argv[])
-{
-  if(argc<4)
-    {
-      cerr<<argv[0]<<" <sbp.conf> <rout_kpc> <cfunc_erg> [cfunc2_erg ...]"<<endl;
-      return -1;
-    }
-  //initialize the parameters list
-  ifstream cfg_file(argv[1]);
-  assert(cfg_file.is_open());
-  cfg_map cfg=parse_cfg_file(cfg_file);
-
-  const double z=cfg.z;
-
-  //initialize the radius list, sbp list and sbp error list
-  std::vector<double> radii;
-  std::vector<double> sbps;
-  std::vector<double> sbpe;
-  std::vector<double> radii_all;
-  std::vector<double> sbps_all;
-  std::vector<double> sbpe_all;
-  //prepend the zero point to radius list
-  radii.push_back(0.0);
-  radii_all.push_back(0.0);
-  //read sbp and sbp error data
-  ifstream ifs(cfg.sbp_data.c_str());
-  std::string line;
-  if (ifs.is_open())
-    {
-      while(std::getline(ifs, line))
-        {
-          if (line.empty())
-            continue;
-
-          std::istringstream iss(line);
-          double x, xe, y, ye;
-          if ((iss >> x >> xe >> y >> ye))
-            {
-              std::cerr << "sbprofile data: "
-                        << x << ", " << xe << ", " << y << ", " << ye
-                        << std::endl;
-              radii.push_back(x+xe);  /* NOTE: use outer radii of regions */
-              radii_all.push_back(x+xe);
-              sbps.push_back(y);
-              sbps_all.push_back(y);
-              sbpe.push_back(ye);
-              sbpe_all.push_back(ye);
-            }
-          else
-            {
-              std::cerr << "skipped line: " << line << std::endl;
-            }
-        }
-    }
-  else
-    {
-      std::cerr << "ERROR: cannot open file: " << cfg.sbp_data.c_str()
-                << std::endl;
-      return 1;
-    }
-
-  //initialize the cm/pixel value
-  double cm_per_pixel=cfg.cm_per_pixel;
-  double rmin;
-  if(cfg.rmin_pixel>0)
-    {
-      rmin=cfg.rmin_pixel;
-    }
-  else
-    {
-      rmin=cfg.rmin_kpc*kpc/cm_per_pixel;
-    }
-
-  cerr<<"rmin="<<rmin<<" (pixel)"<<endl;
-  std::list<double> radii_tmp,sbps_tmp,sbpe_tmp;
-  radii_tmp.resize(radii.size());
-  sbps_tmp.resize(sbps.size());
-  sbpe_tmp.resize(sbpe.size());
-  copy(radii.begin(),radii.end(),radii_tmp.begin());
-  copy(sbps.begin(),sbps.end(),sbps_tmp.begin());
-  copy(sbpe.begin(),sbpe.end(),sbpe_tmp.begin());
-  for(list<double>::iterator i=radii_tmp.begin();i!=radii_tmp.end();)
-    {
-      if(*i<rmin)
-	{
-	  radii_tmp.pop_front();
-	  sbps_tmp.pop_front();
-	  sbpe_tmp.pop_front();
-	  i=radii_tmp.begin();
-	  continue;
-	}
-      ++i;
-    }
-  radii.resize(radii_tmp.size());
-  sbps.resize(sbps_tmp.size());
-  sbpe.resize(sbpe_tmp.size());
-  copy(radii_tmp.begin(),radii_tmp.end(),radii.begin());
-  copy(sbps_tmp.begin(),sbps_tmp.end(),sbps.begin());
-  copy(sbpe_tmp.begin(),sbpe_tmp.end(),sbpe.begin());
-
-  //read cooling function data
-  spline_func_obj cf;
-  for(ifstream ifs(cfg.cfunc_profile.c_str());;)
-    {
-      assert(ifs.is_open());
-      double x,y;
-      ifs>>x>>y;
-      if(!ifs.good())
-	{
-	  break;
-	}
-      cerr<<x<<"\t"<<y<<endl;
-      if(x>radii.back())
-	{
-	  break;
-	}
-      cf.add_point(x,y);
-    }
-  cf.gen_spline();
-
-  //read temperature profile data
-  spline_func_obj Tprof;
-  int tcnt=0;
-  for(ifstream ifs1(cfg.tprofile.c_str());;++tcnt)
-    {
-      assert(ifs1.is_open());
-      double x,y;
-      ifs1>>x>>y;
-      if(!ifs1.good())
-      {
-	break;
-      }
-      cerr<<x<<"\t"<<y<<endl;
-#if 0
-      if(tcnt==0)
-	{
-	  Tprof.add_point(0,y);
-	}
-#endif
-      Tprof.add_point(x,y);
-    }
-
-
-  Tprof.gen_spline();
-
-  default_data_set<std::vector<double>,std::vector<double> > ds;
-  ds.add_data(data<std::vector<double>,std::vector<double> >(radii,sbps,sbpe,sbpe,radii,radii));
-
-  //initial fitter
-  fitter<vector<double>,vector<double>,vector<double>,double> f;
-  f.load_data(ds);
-  //initial the object, which is used to calculate projection effect
-  projector<double> a;
-  bool tie_beta=false;
-  if(cfg.param_map.find("beta")!=cfg.param_map.end()
-     &&cfg.param_map.find("beta1")==cfg.param_map.end()
-     &&cfg.param_map.find("beta2")==cfg.param_map.end())
-    {
-      dbeta2<double> dbetao;
-      a.attach_model(dbetao);
-      tie_beta=true;
-    }
-  else if((cfg.param_map.find("beta1")!=cfg.param_map.end()
-	   ||cfg.param_map.find("beta2")!=cfg.param_map.end())
-	  &&cfg.param_map.find("beta")==cfg.param_map.end())
-    {
-      dbeta<double> dbetao;
-      a.attach_model(dbetao);
-      tie_beta=false;
-    }
-  else
-    {
-      cerr<<"Error, cannot decide whether to tie beta together or let them vary freely!"<<endl;
-      assert(0);
-    }
-
-  //attach the cooling function
-  a.attach_cfunc(cf);
-  a.set_cm_per_pixel(cm_per_pixel);
-
-  f.set_model(a);
-  //chi^2 statistic
-  vchisq<double> c;
-  c.verbose(true);
-  c.set_limit();
-  f.set_statistic(c);
-  //optimization method
-  f.set_opt_method(powell_method<double,std::vector<double> >());
-  //initialize the initial values
-  /*
-  double =0;
-  double rc1=0;
-  double n02=0;
-  double rc2=0;
-  double beta=0;
-  */
-  double bkg_level=0;
-
-  if(tie_beta)
-    {
-      f.set_param_value("beta",.7);
-      f.set_param_lower_limit("beta",.3);
-      f.set_param_upper_limit("beta",1.4);
-    }
-  else
-    {
-      f.set_param_value("beta1",.7);
-      f.set_param_lower_limit("beta1",.3);
-      f.set_param_upper_limit("beta1",1.4);
-      f.set_param_value("beta2",.7);
-      f.set_param_lower_limit("beta2",.3);
-      f.set_param_upper_limit("beta2",1.4);
-    }
-  for(std::map<std::string,std::vector<double> >::iterator i=cfg.param_map.begin();
-      i!=cfg.param_map.end();++i)
-    {
-      std::string pname=i->first;
-      f.set_param_value(pname,i->second.at(0));
-      if(i->second.size()==3)
-	{
-	  double a1=i->second[1];
-	  double a2=i->second[2];
-	  double u=std::max(a1,a2);
-	  double l=std::min(a1,a2);
-	  f.set_param_upper_limit(pname,u);
-	  f.set_param_lower_limit(pname,l);
-	}
-      else
-	{
-	  if(pname=="beta"||pname=="beta1"||pname=="beta2")
-	    {
-	      f.set_param_lower_limit(pname,.3);
-	      f.set_param_upper_limit(pname,1.4);
-	    }
-	}
-    }
-
-
-
-  //perform the fitting, first freeze beta1, beta2, rc1, and rc2
-  if(tie_beta)
-    {
-      f.set_param_modifier(freeze_param<std::vector<double>,std::vector<double>,std::vector<double>,std::string>("beta")+
-			   freeze_param<std::vector<double>,std::vector<double>,std::vector<double>,std::string>("rc1")+
-			   freeze_param<std::vector<double>,std::vector<double>,std::vector<double>,std::string>("rc2")
-			   );
-    }
-  else
-    {
-      f.set_param_modifier(freeze_param<std::vector<double>,std::vector<double>,std::vector<double>,std::string>("beta1")+
-			   freeze_param<std::vector<double>,std::vector<double>,std::vector<double>,std::string>("beta2")+
-			   freeze_param<std::vector<double>,std::vector<double>,std::vector<double>,std::string>("rc1")+
-			   freeze_param<std::vector<double>,std::vector<double>,std::vector<double>,std::string>("rc2")
-			   );
-    }
-
-  f.fit();
-
-  f.clear_param_modifier();
-
-  //then perform the fitting, freeze beta1 and beta2
-  //f.set_param_modifier(freeze_param<std::vector<double>,std::vector<double>,std::vector<double>,std::string>("beta"));
-  //f.set_param_modifier(freeze_param<std::vector<double>,std::vector<double>,std::vector<double>,std::string>("bkg"));
-  f.fit();
-  //f.clear_param_modifier();
-
-  //finally thaw all parameters
-  f.fit();
-
-  /*
-  double beta1=0;
-  double beta2=0;
-
-  n01=f.get_param_value("n01");
-  rc1=f.get_param_value("rc1");
-  n02=f.get_param_value("n02");
-  rc2=f.get_param_value("rc2");
-  if(tie_beta)
-    {
-      beta=f.get_param_value("beta");
-      beta1=beta;
-      beta2=beta;
-    }
-  else
-    {
-      beta1=f.get_param_value("beta1");
-      beta2=f.get_param_value("beta2");
-    }
-  */
-
-  //output the params
-  ofstream param_output("lx_dbeta_param.txt");
-  //output the datasets and fitting results
-  for(size_t i=0;i<f.get_num_params();++i)
-    {
-      if(f.get_param_info(i).get_name()=="rc1")
-	{
-	  cerr<<"rc1_kpc"<<"\t"<<abs(f.get_param_info(i).get_value())*cm_per_pixel/kpc<<endl;
-	  param_output<<"rc1_kpc"<<"\t"<<abs(f.get_param_info(i).get_value())*cm_per_pixel/kpc<<endl;
-	}
-      if(f.get_param_info(i).get_name()=="rc2")
-	{
-	  cerr<<"rc2_kpc"<<"\t"<<abs(f.get_param_info(i).get_value())*cm_per_pixel/kpc<<endl;
-	  param_output<<"rc2_kpc"<<"\t"<<abs(f.get_param_info(i).get_value())*cm_per_pixel/kpc<<endl;
-	}
-      cerr<<f.get_param_info(i).get_name()<<"\t"<<abs(f.get_param_info(i).get_value())<<endl;
-      param_output<<f.get_param_info(i).get_name()<<"\t"<<abs(f.get_param_info(i).get_value())<<endl;
-    }
-  cerr<<"reduced_chi^2="<<f.get_statistic_value()/(radii.size()-f.get_model().get_num_free_params())<<endl;
-  param_output<<"reduced_chi^2="<<f.get_statistic_value()/(radii.size()-f.get_model().get_num_free_params())<<endl;
-
-  //c.verbose(false);
-  //f.set_statistic(c);
-  //f.fit();
-  std::vector<double> p=f.get_all_params();
-  f.clear_param_modifier();
-  std::vector<double> mv=f.eval_model(radii,p);
-
-
-  ofstream ofs_sbp("lx_sbp_fit.qdp");
-  ofs_sbp<<"read serr 2"<<endl;
-  ofs_sbp<<"skip single"<<endl;
-
-  ofs_sbp<<"line on 2"<<endl;
-  ofs_sbp<<"line on 3"<<endl;
-  ofs_sbp<<"line on 4"<<endl;
-  ofs_sbp<<"line on 5"<<endl;
-  ofs_sbp<<"line on 7"<<endl;
-  ofs_sbp<<"ls 2 on 7"<<endl;
-  ofs_sbp<<"ls 2 on 3"<<endl;
-  ofs_sbp<<"ls 2 on 4"<<endl;
-  ofs_sbp<<"ls 2 on 5"<<endl;
-
-
-
-  ofs_sbp<<"!LAB  POS Y  4.00"<<endl;
-  ofs_sbp<<"!LAB  ROT"<<endl;
-  ofs_sbp<<"win 1"<<endl;
-  ofs_sbp<<"yplot 1 2 3 4 5"<<endl;
-  ofs_sbp<<"loc 0 0 1 1"<<endl;
-  ofs_sbp<<"vie .1 .4 .9 .9"<<endl;
-  ofs_sbp<<"la y cnt/s/pixel/cm\\u2\\d"<<endl;
-  ofs_sbp<<"log x"<<endl;
-  ofs_sbp<<"log y"<<endl;
-  ofs_sbp<<"r x "<<(radii[1]+radii[0])/2*cm_per_pixel/kpc<<" "<<(radii[sbps.size()-2]+radii[sbps.size()-1])/2*cm_per_pixel/kpc<<endl;
-  ofs_sbp<<"win 2"<<endl;
-  ofs_sbp<<"yplot 6 7"<<endl;
-  ofs_sbp<<"loc 0 0 1 1"<<endl;
-  ofs_sbp<<"vie .1 .1 .9 .4"<<endl;
-  ofs_sbp<<"la x radius (kpc)"<<endl;
-  ofs_sbp<<"la y chi"<<endl;
-  ofs_sbp<<"log x"<<endl;
-  ofs_sbp<<"log y off"<<endl;
-  ofs_sbp<<"r x "<<(radii[1]+radii[0])/2*cm_per_pixel/kpc<<" "<<(radii[sbps.size()-2]+radii[sbps.size()-1])/2*cm_per_pixel/kpc<<endl;
-  for(size_t i=1;i<sbps.size();++i)
-    {
-      double x=(radii[i]+radii[i-1])/2;
-      double y=sbps[i-1];
-      double ye=sbpe[i-1];
-      //double ym=mv[i-1];
-      ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<y<<"\t"<<ye<<endl;
-    }
-  ofs_sbp<<"no no no"<<endl;
-  for(size_t i=1;i<sbps.size();++i)
-    {
-      double x=(radii[i]+radii[i-1])/2;
-      //double y=sbps[i-1];
-      //double ye=sbpe[i-1];
-      double ym=mv[i-1];
-      ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<ym<<"\t"<<0<<endl;
-    }
-  //bkg
-  ofs_sbp<<"no no no"<<endl;
-  bkg_level=abs(f.get_param_value("bkg"));
-  for(size_t i=0;i<sbps.size();++i)
-    {
-      double x=(radii[i]+radii[i-1])/2;
-      ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<bkg_level<<"\t0"<<endl;
-    }
-  //rc1
-  ofs_sbp<<"no no no"<<endl;
-  double rc1_kpc=abs(f.get_param_value("rc1")*cm_per_pixel/kpc);
-  double max_sbp=*max_element(sbps.begin(),sbps.end());
-  double min_sbp=*min_element(sbps.begin(),sbps.end());
-  for(double x=min_sbp;x<=max_sbp;x+=(max_sbp-min_sbp)/100)
-    {
-      ofs_sbp<<rc1_kpc<<"\t"<<x<<"\t"<<"0"<<endl;
-    }
-  //rc2
-  ofs_sbp<<"no no no"<<endl;
-  double rc2_kpc=abs(f.get_param_value("rc2")*cm_per_pixel/kpc);
-  for(double x=min_sbp;x<=max_sbp;x+=(max_sbp-min_sbp)/100)
-    {
-      ofs_sbp<<rc2_kpc<<"\t"<<x<<"\t"<<"0"<<endl;
-    }
-  //resid
-  ofs_sbp<<"no no no"<<endl;
-  for(size_t i=1;i<sbps.size();++i)
-    {
-      double x=(radii[i]+radii[i-1])/2;
-      //double y=sbps[i-1];
-      //double ye=sbpe[i-1];
-      double ym=mv[i-1];
-      ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<(ym-sbps[i-1])/sbpe[i-1]<<"\t"<<1<<endl;
-    }
-  //zero level in resid map
-  ofs_sbp<<"no no no"<<endl;
-  for(size_t i=1;i<sbps.size();++i)
-    {
-      double x=(radii[i]+radii[i-1])/2;
-      //double y=sbps[i-1];
-      //double ye=sbpe[i-1];
-      //double ym=mv[i-1];
-      ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<0<<"\t"<<0<<endl;
-    }
-
-
-  mv=f.eval_model_raw(radii,p);
-  ofstream ofs_rho("lx_rho_fit.qdp");
-  ofstream ofs_rho_data("lx_rho_fit.dat");
-  //ofstream ofs_entropy("entropy.qdp");
-  ofs_rho<<"la x radius (kpc)"<<endl;
-  ofs_rho<<"la y density (cm\\u-3\\d)"<<endl;
-  /*
-  for(int i=1;i<sbps.size();++i)
-    {
-      double x=(radii[i]+radii[i-1])/2;
-      double ym=mv[i-1];
-      ofs_rho<<x*cm_per_pixel/kpc<<"\t"<<ym<<endl;
-    }
-  */
-
-  p.back()=0;
-
-  radii.clear();
-  double rout=atof(argv[2])*kpc;
-  for(double r=0;r<rout;r+=1*kpc)//step size=1kpc
-    {
-      double r_pix=r/cm_per_pixel;
-      radii.push_back(r_pix);
-    }
-
-  double Da=cm_per_pixel/(.492/3600./180.*pi);
-  double Dl=Da*(1+z)*(1+z);
-  cout<<"dl="<<Dl/kpc<<endl;
-  for(int n=3;n<argc;++n)
-    {
-      spline_func_obj cf_erg;
-      for(ifstream ifs(argv[n]);;)
-	{
-	  assert(ifs.is_open());
-	  double x,y;
-	  ifs>>x>>y;
-	  if(!ifs.good())
-	    {
-	      break;
-	    }
-	  //cerr<<x<<"\t"<<y<<endl;
-
-	  cf_erg.add_point(x,y);//change with source
-	}
-      cf_erg.gen_spline();
-
-      projector<double>& pj=dynamic_cast<projector<double>&>(f.get_model());
-      pj.attach_cfunc(cf_erg);
-
-      mv=f.eval_model_raw(radii,p);
-      double flux_erg=0;
-      for(size_t i=0;i<radii.size()-1;++i)
-	{
-	  double S=pi*(radii[i+1]+radii[i])*(radii[i+1]-radii[i]);
-	  flux_erg+=S*mv[i];
-	}
-      cout<<flux_erg*4*pi*Dl*Dl<<endl;
-      cout<<flux_erg<<endl;
-      param_output<<"Lx"<<n-2<<"\t"<<flux_erg*4*pi*Dl*Dl<<endl;
-      param_output<<"Fx"<<n-2<<"\t"<<flux_erg<<endl;
-    }
-}
diff --git a/mass_profile/chisq.hpp b/mass_profile/chisq.hpp
deleted file mode 100644
index 0e58200..0000000
--- a/mass_profile/chisq.hpp
+++ /dev/null
@@ -1,206 +0,0 @@
-/**

-   \file chisq.hpp

-   \brief chi-square statistic

-   \author Junhua Gu

- */

-

-#ifndef CHI_SQ_HPP

-#define CHI_SQ_HPP

-

-#define OPT_HEADER

-

-#include <core/fitter.hpp>

-#include <iostream>

-#include <vector>

-#include <misc/optvec.hpp>

-#include <cmath>

-

-using std::cerr;

-using std::endl;

-

-namespace opt_utilities

-{

-  static const int display_interval=10;

-  /**

-     \brief chi-square statistic

-     \tparam Ty the return type of model

-     \tparam Tx the type of the self-var

-     \tparam Tp the type of model parameter

-     \tparam Ts the type of the statistic

-     \tparam Tstr the type of the string used

-   */

-  template<typename Ty,typename Tx,typename Tp,typename Ts,typename Tstr>

-  class chisq

-    :public statistic<Ty,Tx,Tp,Ts,Tstr>

-  {

-  };

-  template<>

-  class chisq<double,double,std::vector<double>,double,std::string>

-    :public statistic<double,double,std::vector<double> ,double,std::string>

-  {

-  public:

-    typedef double Ty;

-    typedef double Tx;

-    typedef std::vector<double> Tp;

-    typedef double Ts;

-    typedef std::string Tstr;

-  private:

-    bool verb;

-    bool limit_bound;

-    int n;

-

-    statistic<Ty,Tx,Tp,Ts,Tstr>* do_clone()const

-    {

-      // return const_cast<statistic<Ty,Tx,Tp>*>(this);

-      return new chisq<Ty,Tx,Tp,Ts,Tstr>(*this);

-    }

-

-    const char* do_get_type_name()const

-    {

-      return "chi^2 statistics (specialized for double)";

-    }

-  public:

-    void verbose(bool v)

-    {

-      verb=v;

-    }

-

-    void set_limit()

-    {

-      limit_bound=true;

-    }

-

-    void clear_limit()

-    {

-      limit_bound=false;

-    }

-  public:

-    chisq()

-      :verb(true),limit_bound(false)

-    {}

-

-    Ty do_eval(const Tp& p)

-    {

-      if(limit_bound)

-	{

-	  Tp p1=this->get_fitter().get_model().reform_param(p);

-	  for(size_t i=0;i<p1.size();++i)

-	    {

-	      if(p1[i]>this->get_fitter().get_param_info(i).get_upper_limit()||

-		 p1[i]<this->get_fitter().get_param_info(i).get_lower_limit())

-		{

-		  return 1e99;

-		}

-	    }

-	}

-      Ty result(0);

-      std::vector<float> vx;

-      std::vector<float> vy;

-      std::vector<float> vye1;

-      std::vector<float> vye2;

-      std::vector<float> my;

-      float xmin=1e99,xmax=-1e99,ymin=1e99,ymax=-1e99;

-      if(verb)

-	{

-	  n++;

-	  if(n%display_interval==0)

-	    {

-	      vx.resize(this->get_data_set().size());

-	      vy.resize(this->get_data_set().size());

-	      vye1.resize(this->get_data_set().size());

-	      vye2.resize(this->get_data_set().size());

-	      my.resize(this->get_data_set().size());

-	    }

-	}

-

-      for(int i=(this->get_data_set()).size()-1;i>=0;--i)

-	{

-

-#ifdef HAVE_X_ERROR

-	  Tx x1=this->get_data_set().get_data(i).get_x()-this->get_data_set().get_data(i).get_x_lower_err();

-	  Tx x2=this->get_data_set().get_data(i).get_x()+this->get_data_set().get_data(i).get_x_upper_err();

-	  Tx x=this->get_data_set().get_data(i).get_x();

-	  Ty errx1=(this->eval_model(x1,p)-this->eval_model(x,p));

-	  Ty errx2=(this->eval_model(x2,p)-this->eval_model(x,p));

-	  //Ty errx=0;

-#else

-	  Ty errx1=0;

-	  Ty errx2=0;

-#endif

-

-	  Ty y_model=this->eval_model(this->get_data_set().get_data(i).get_x(),p);

-	  Ty y_obs=this->get_data_set().get_data(i).get_y();

-	  Ty y_err;

-

-	  Ty errx=0;

-	  if(errx1<errx2)

-	    {

-	      if(y_obs<y_model)

-		{

-		  errx=errx1>0?errx1:-errx1;

-		}

-	      else

-		{

-		  errx=errx2>0?errx2:-errx2;

-		}

-	    }

-	  else

-	    {

-	      if(y_obs<y_model)

-		{

-		  errx=errx2>0?errx2:-errx2;

-		}

-	      else

-		{

-		  errx=errx1>0?errx1:-errx1;

-		}

-	    }

-

-	  if(y_model>y_obs)

-	    {

-	      y_err=this->get_data_set().get_data(i).get_y_upper_err();

-	    }

-	  else

-	    {

-	      y_err=this->get_data_set().get_data(i).get_y_lower_err();

-	    }

-

-	  Ty chi=(y_obs-y_model)/std::sqrt(y_err*y_err+errx*errx);

-

-	  result+=chi*chi;

-

-	  if(verb&&n%display_interval==0)

-	    {

-	      vx.at(i)=this->get_data_set().get_data(i).get_x();

-	      vy.at(i)=this->get_data_set().get_data(i).get_y();

-	      vye1.at(i)=std::abs(this->get_data_set().get_data(i).get_y_lower_err());

-	      vye2.at(i)=std::abs(this->get_data_set().get_data(i).get_y_upper_err());

-	      my.at(i)=y_model;

-

-	      xmin=std::min(vx.at(i),xmin);

-	      ymin=std::min(vy.at(i),ymin-vye1[i]);

-	      xmax=std::max(vx.at(i),xmax);

-	      ymax=std::max(vy.at(i),ymax+vye2[i]);

-	    }

-

-	}

-      if(verb)

-	{

-	  if(n%display_interval==0)

-	    {

-	      cerr<<result<<"\t";

-	      for(size_t i=0;i<get_size(p);++i)

-		{

-		  cerr<<get_element(p,i)<<",";

-		}

-	      cerr<<endl;

-	      //cerr<<x1<<"\t"<<x2<<endl;

-	    }

-	}

-

-      return result;

-    }

-  };

-}

-

-#endif /* CHI_SQ_HPP */

diff --git a/mass_profile/dbeta.hpp b/mass_profile/dbeta.hpp
deleted file mode 100644
index 7246b44..0000000
--- a/mass_profile/dbeta.hpp
+++ /dev/null
@@ -1,108 +0,0 @@
-#ifndef DBETA
-#define DBETA
-#include "projector.hpp"
-
-/**
-   dbeta: double beta model for density
-   dbeta2: double beta model for density with only one beta
-*/
-
-
-namespace opt_utilities
-{
-  template <typename T>
-  class dbeta
-    :public model<std::vector<T>,std::vector<T>,std::vector<T> >
-  {
-  public:
-    dbeta()
-    {
-      this->push_param_info(param_info<std::vector<T>,std::string>("n01",1));
-      this->push_param_info(param_info<std::vector<T>,std::string>("beta1",.66));
-      this->push_param_info(param_info<std::vector<T>,std::string>("rc1",100));
-
-      this->push_param_info(param_info<std::vector<T>,std::string>("n02",1));
-      this->push_param_info(param_info<std::vector<T>,std::string>("beta2",.67));
-      this->push_param_info(param_info<std::vector<T>,std::string>("rc2",110));
-
-    }
-
-  public:
-    dbeta<T>* do_clone()const
-    {
-      return new dbeta<T>(*this);
-    }
-
-    std::vector<T> do_eval(const std::vector<T> & x,
-			   const std::vector<T>& p)
-    {
-      T n01=std::abs(p[0]);
-      T beta1=p[1];
-      T rc1=p[2];
-
-      T n02=std::abs(p[3]);
-      T beta2=p[4];
-      T rc2=p[5];
-
-
-
-      std::vector<T> result(x.size()-1);
-      for(size_t i=1;i<x.size();++i)
-	{
-	  T xi=(x[i]+x[i-1])/2;
-	  T yi=0;
-	  yi=n01*pow(1+xi*xi/rc1/rc1,-3./2.*beta1)+n02*pow(1+xi*xi/rc2/rc2,-3./2.*beta2);
-	  result[i-1]=yi;
-	}
-      return result;
-    }
-  };
-
-  template <typename T>
-  class dbeta2
-    :public model<std::vector<T>,std::vector<T>,std::vector<T> >
-  {
-  public:
-    dbeta2()
-    {
-      this->push_param_info(param_info<std::vector<T>,std::string>("n01",1));
-      this->push_param_info(param_info<std::vector<T>,std::string>("rc1",100));
-      this->push_param_info(param_info<std::vector<T>,std::string>("n02",1));
-      this->push_param_info(param_info<std::vector<T>,std::string>("rc2",110));
-      this->push_param_info(param_info<std::vector<T>,std::string>("beta",.67));
-
-    }
-
-  public:
-    dbeta2<T>* do_clone()const
-    {
-      return new dbeta2<T>(*this);
-    }
-
-    std::vector<T> do_eval(const std::vector<T> & x,
-			   const std::vector<T>& p)
-    {
-      T n01=std::abs(p[0]);
-      T rc1=p[1];
-
-      T n02=std::abs(p[2]);
-      T rc2=p[3];
-      T beta=p[4];
-      T beta1=beta;
-      T beta2=beta;
-
-      std::vector<T> result(x.size()-1);
-      for(size_t i=1;i<x.size();++i)
-	{
-	  T xi=(x[i]+x[i-1])/2;
-	  T yi=0;
-	  yi=n01*pow(1+xi*xi/rc1/rc1,-3./2.*beta1)+n02*pow(1+xi*xi/rc2/rc2,-3./2.*beta2);
-	  result[i-1]=yi;
-	}
-      return result;
-    }
-  };
-
-}
-
-#endif
diff --git a/mass_profile/dump_fit_qdp.cpp b/mass_profile/dump_fit_qdp.cpp
deleted file mode 100644
index 85307d1..0000000
--- a/mass_profile/dump_fit_qdp.cpp
+++ /dev/null
@@ -1,55 +0,0 @@
-#include "dump_fit_qdp.hpp"
-
-namespace opt_utilities
-{
-  const static double kpc=3.086E21;
-  void dump_sbp_beta(std::ostream& os,fitter<double,double,std::vector<double>,double,std::string>& f,double cm_per_pixel,const std::vector<double>& r,const std::vector<double>& y,const std::vector<double>& ye)
-  {
-    os<<"read serr 1 2"<<std::endl;
-    os<<"skip single"<<std::endl;
-    os<<"la x \"radius (kpc)\""<<std::endl;
-    os<<"la y \"surface brightness (cts s\\u-1\\d pixel\\u-2\\d)\""<<std::endl;
-    os<<"li on 2"<<std::endl;
-    for(size_t i=1;i<r.size();++i)
-      {
-	os<<(r[i]+r[i-1])/2*cm_per_pixel/kpc<<"\t"<<(r[i]-r[i-1])/2*cm_per_pixel/kpc<<"\t"<<y[i-1]<<"\t"<<ye[i-1]<<std::endl;
-      }
-    os<<"no no no"<<std::endl;
-    std::vector<double> p=f.get_all_params();
-    for(size_t i=1;i<r.size();++i)
-      {
-	double x=(r[i]+r[i-1])/2;
-	os<<x*cm_per_pixel/kpc<<"\t"<<0<<"\t"<<f.eval_model_raw(x,p)<<"\t"<<0<<std::endl;
-      }
-  }
-
-  void dump_rho_beta(std::ostream& os,fitter<std::vector<double>,std::vector<double>,std::vector<double>,double,std::string>& f,double cm_per_pixel,const std::vector<double>& r,const std::vector<double>& sbps,const std::vector<double>& sbpe)
-  {
-    os<<"read serr 1 2"<<std::endl;
-    os<<"skip single"<<std::endl;
-    os<<"la x \"radius (kpc)\""<<std::endl;
-    os<<"la y \"density (cm\\u-3\\d)\""<<std::endl;
-    os<<"li on 2"<<std::endl;
-
-    for(size_t i=1;i<r.size();++i)
-    {
-      double x=(r[i]+r[i-1])/2;
-      double y=sbps[i-1];
-      double ye=sbpe[i-1];
-      os<<x*cm_per_pixel/kpc<<"\t0\t"<<y<<"\t"<<ye<<std::endl;
-    }
-
-    os<<"no no no"<<std::endl;
-    std::vector<double> p=f.get_all_params();
-    std::vector<double> mv=f.eval_model_raw(r,p);
-    for(size_t i=1;i<r.size();++i)
-    {
-      double x=(r[i]+r[i-1])/2;
-      double y=mv[i-1];
-      double ye=0;
-      os<<x*cm_per_pixel/kpc<<"\t0\t"<<y<<"\t"<<ye<<std::endl;
-    }
-
-  }
-
-};
diff --git a/mass_profile/dump_fit_qdp.hpp b/mass_profile/dump_fit_qdp.hpp
deleted file mode 100644
index 8364dfd..0000000
--- a/mass_profile/dump_fit_qdp.hpp
+++ /dev/null
@@ -1,16 +0,0 @@
-#ifndef DUMP_FIT_QDP_HPP
-#define DUMP_FIT_QDP_HPP
-
-#include <core/fitter.hpp>
-#include <vector>
-#include <iostream>
-#include <string>
-
-namespace opt_utilities
-{
-  void dump_sbp_beta(std::ostream& os,fitter<double,double,std::vector<double>,double,std::string>& f,double cm_per_pixel,const std::vector<double>& r,const std::vector<double>& y,const std::vector<double>& ye);
-  void dump_rho_beta(std::ostream& os,fitter<std::vector<double>,std::vector<double>,std::vector<double>,double,std::string>& f,double cm_per_pixel,const std::vector<double>& r,const std::vector<double>& sbps,const std::vector<double>& sbpe);
-  void dump_rho_dbeta(std::ostream& os,fitter<std::vector<double>,std::vector<double>,std::vector<double>,double,std::string>& f,double cm_per_pixel);
-};
-
-#endif
diff --git a/mass_profile/fit_beta_sbp.cpp b/mass_profile/fit_beta_sbp.cpp
deleted file mode 100644
index 295fa1e..0000000
--- a/mass_profile/fit_beta_sbp.cpp
+++ /dev/null
@@ -1,534 +0,0 @@
-/*
-  Perform a double-beta density model fitting to the surface brightness data
-  Author: Junhua Gu
-  Last modified: 2016.06.07
-  This code is distributed with no warrant
-*/
-
-#include <iostream>
-#include <fstream>
-#include <sstream>
-#include <list>
-#include <algorithm>
-#include "beta_cfg.hpp"
-#include "dump_fit_qdp.hpp"
-#include "report_error.hpp"
-#include "vchisq.hpp"
-#include "chisq.hpp"
-#include "beta.hpp"
-#include "models/beta1d.hpp"
-#include <data_sets/default_data_set.hpp>
-#include <methods/powell/powell_method.hpp>
-#include <core/freeze_param.hpp>
-#include <error_estimator/error_estimator.hpp>
-#include "spline.hpp"
-
-using namespace std;
-using namespace opt_utilities;
-//double s=5.63136645E20;
-const double kpc=3.086E21;//kpc in cm
-const double Mpc=kpc*1000;
-
-double beta_func(double r, double n0, double rc, double beta)
-{
-  return abs(n0) * pow(1+r*r/rc/rc, -3./2.*abs(beta));
-}
-
-//calculate critical density from z, under following cosmological constants
-static double calc_critical_density(double z,
-				    const double H0=2.3E-18,
-				      const double Omega_m=.27)
-{
-  const double G=6.673E-8;//cm^3 g^-1 s^2
-  const double E=std::sqrt(Omega_m*(1+z)*(1+z)*(1+z)+1-Omega_m);
-  const double H=H0*E;
-  return 3*H*H/8/pi/G;
-}
-
-
-//A class enclosing the spline interpolation method
-class spline_func_obj
-  :public func_obj<double,double>
-{
-  //has an spline object
-  spline<double> spl;
-public:
-  //This function is used to calculate the intepolated value
-  double do_eval(const double& x)
-  {
-    return spl.get_value(x);
-  }
-
-  //we need this function, when this object is performing a clone of itself
-  spline_func_obj* do_clone()const
-  {
-    return new spline_func_obj(*this);
-  }
-
-public:
-  //add points to the spline object, after which the spline will be initialized
-  void add_point(double x,double y)
-  {
-    spl.push_point(x,y);
-  }
-
-  //before getting the intepolated value, the spline should be initialzied by calling this function
-  void gen_spline()
-  {
-    spl.gen_spline(0,0);
-  }
-};
-
-int main(int argc,char* argv[])
-{
-  if(argc!=2)
-    {
-      cerr<<argv[0]<<" <configure file>"<<endl;
-      return -1;
-    }
-  //initialize the parameters list
-  ifstream cfg_file(argv[1]);
-  assert(cfg_file.is_open());
-  cfg_map cfg=parse_cfg_file(cfg_file);
-
-  const double z=cfg.z;
-
-  //initialize the radius list, sbp list and sbp error list
-  std::vector<double> radii;
-  std::vector<double> sbps;
-  std::vector<double> sbpe;
-  std::vector<double> radii_all;
-  std::vector<double> sbps_all;
-  std::vector<double> sbpe_all;
-  //prepend the zero point to radius list
-  radii.push_back(0.0);
-  radii_all.push_back(0.0);
-  //read sbp and sbp error data
-  cerr << "Read surface brightness profile data ..." << endl;
-  ifstream ifs(cfg.sbp_data.c_str());
-  std::string line;
-  if (ifs.is_open())
-    {
-      while(std::getline(ifs, line))
-        {
-          if (line.empty())
-            continue;
-
-          std::istringstream iss(line);
-          double x, xe, y, ye;
-          if ((iss >> x >> xe >> y >> ye))
-            {
-              std::cerr << "sbprofile data: "
-                        << x << ", " << xe << ", " << y << ", " << ye
-                        << std::endl;
-              radii.push_back(x+xe);  /* NOTE: use outer radii of regions */
-              radii_all.push_back(x+xe);
-              sbps.push_back(y);
-              sbps_all.push_back(y);
-              sbpe.push_back(ye);
-              sbpe_all.push_back(ye);
-            }
-          else
-            {
-              std::cerr << "skipped line: " << line << std::endl;
-            }
-        }
-    }
-  else
-    {
-      std::cerr << "ERROR: cannot open file: " << cfg.sbp_data.c_str()
-                << std::endl;
-      return 1;
-    }
-
-  //initialize the cm/pixel value
-  double cm_per_pixel=cfg.cm_per_pixel;
-  double rmin;
-  if(cfg.rmin_pixel>0)
-    {
-      rmin=cfg.rmin_pixel;
-    }
-  else
-    {
-      rmin=cfg.rmin_kpc*kpc/cm_per_pixel;
-    }
-
-  cerr<<"rmin="<<rmin<<" (pixel)"<<endl;
-  std::list<double> radii_tmp,sbps_tmp,sbpe_tmp;
-  radii_tmp.resize(radii.size());
-  sbps_tmp.resize(sbps.size());
-  sbpe_tmp.resize(sbpe.size());
-  copy(radii.begin(),radii.end(),radii_tmp.begin());
-  copy(sbps.begin(),sbps.end(),sbps_tmp.begin());
-  copy(sbpe.begin(),sbpe.end(),sbpe_tmp.begin());
-  for(list<double>::iterator i=radii_tmp.begin();i!=radii_tmp.end();)
-    {
-      if(*i<rmin)
-	{
-	  radii_tmp.pop_front();
-	  sbps_tmp.pop_front();
-	  sbpe_tmp.pop_front();
-	  i=radii_tmp.begin();
-	  continue;
-	}
-      ++i;
-    }
-  radii.resize(radii_tmp.size());
-  sbps.resize(sbps_tmp.size());
-  sbpe.resize(sbpe_tmp.size());
-  copy(radii_tmp.begin(),radii_tmp.end(),radii.begin());
-  copy(sbps_tmp.begin(),sbps_tmp.end(),sbps.begin());
-  copy(sbpe_tmp.begin(),sbpe_tmp.end(),sbpe.begin());
-
-  //read cooling function data
-  cerr << "Read cooling function profile data ..." << endl;
-  spline_func_obj cf;
-  for(ifstream ifs(cfg.cfunc_profile.c_str());;)
-    {
-      assert(ifs.is_open());
-      double x,y;
-      ifs>>x>>y;
-      if(!ifs.good())
-	{
-	  break;
-	}
-      cerr<<x<<"\t"<<y<<endl;
-      if(x>radii.back())
-	{
-      cerr << "radius_max: " << radii.back() << endl;
-	  break;
-	}
-      cf.add_point(x,y);
-    }
-  cf.gen_spline();
-
-  //read temperature profile data
-  cerr << "Read temperature profile data ..." << endl;
-  spline_func_obj Tprof;
-  int tcnt=0;
-  for(ifstream ifs1(cfg.tprofile.c_str());;++tcnt)
-    {
-      assert(ifs1.is_open());
-      double x,y;
-      ifs1>>x>>y;
-      if(!ifs1.good())
-      {
-	break;
-      }
-      cerr<<x<<"\t"<<y<<endl;
-#if 0
-      if(tcnt==0)
-	{
-	  Tprof.add_point(0,y);
-	}
-#endif
-      Tprof.add_point(x,y);
-    }
-  Tprof.gen_spline();
-
-  default_data_set<std::vector<double>,std::vector<double> > ds;
-  ds.add_data(data<std::vector<double>,std::vector<double> >(radii,sbps,sbpe,sbpe,radii,radii));
-
-  //initial fitter
-  fitter<vector<double>,vector<double>,vector<double>,double> f;
-  f.load_data(ds);
-  //initial the object, which is used to calculate projection effect
-  projector<double> a;
-  beta<double> betao;
-  //attach the cooling function
-  a.attach_cfunc(cf);
-  a.set_cm_per_pixel(cm_per_pixel);
-  a.attach_model(betao);
-  f.set_model(a);
-  //chi^2 statistic
-  vchisq<double> c;
-  c.verbose(true);
-  c.set_limit();
-  f.set_statistic(c);
-  //optimization method
-  f.set_opt_method(powell_method<double,std::vector<double> >());
-  //initialize the initial values
-  double n0=0;
-  //double beta=atof(arg_map["beta"].c_str());
-  double beta=0;
-  double rc=0;
-  double bkg_level=0;
-
-  for(std::map<std::string,std::vector<double> >::iterator i=cfg.param_map.begin();
-      i!=cfg.param_map.end();++i)
-    {
-      std::string pname=i->first;
-      f.set_param_value(pname,i->second.at(0));
-      if(i->second.size()==3)
-	{
-	  double a1=i->second[1];
-	  double a2=i->second[2];
-	  double u=std::max(a1,a2);
-	  double l=std::min(a1,a2);
-	  f.set_param_upper_limit(pname,u);
-	  f.set_param_lower_limit(pname,l);
-	}
-      else
-	{
-	  if(pname=="beta")
-	    {
-	      f.set_param_lower_limit(pname,.3);
-	      f.set_param_upper_limit(pname,1.4);
-	    }
-	}
-    }
-
-  f.fit();
-  f.fit();
-  std::vector<double> p=f.get_all_params();
-  n0=f.get_param_value("n0");
-  rc=f.get_param_value("rc");
-  beta=f.get_param_value("beta");
-  //output the datasets and fitting results
-  ofstream param_output("beta_param.txt");
-  for(size_t i=0;i<f.get_num_params();++i)
-    {
-      if(f.get_param_info(i).get_name()=="rc")
-	{
-	  cerr<<"rc_kpc"<<"\t"<<abs(f.get_param_info(i).get_value())*cm_per_pixel/kpc<<endl;
-	  param_output<<"rc_kpc"<<"\t"<<abs(f.get_param_info(i).get_value())*cm_per_pixel/kpc<<endl;
-	}
-      cerr<<f.get_param_info(i).get_name()<<"\t"<<abs(f.get_param_info(i).get_value())<<endl;
-      param_output<<f.get_param_info(i).get_name()<<"\t"<<abs(f.get_param_info(i).get_value())<<endl;
-    }
-  cerr<<"reduced_chi^2="<<f.get_statistic_value()/(radii.size()-f.get_model().get_num_free_params())<<endl;
-  param_output<<"reduced_chi^2="<<f.get_statistic_value()/(radii.size()-f.get_model().get_num_free_params())<<endl;
-
-  std::vector<double> mv=f.eval_model_raw(radii_all,p);
-  int sbps_inner_cut_size=int(sbps_all.size()-sbps.size());
-  ofstream ofs_sbp("sbp_fit.qdp");
-  ofs_sbp<<"read serr 2"<<endl;
-  ofs_sbp<<"skip single"<<endl;
-  ofs_sbp<<"line off "<<endl;
-  if(sbps_inner_cut_size>=1)
-    {
-      ofs_sbp<<"line on 2"<<endl;
-      ofs_sbp<<"line on 3"<<endl;
-      ofs_sbp<<"line on 4"<<endl;
-      ofs_sbp<<"line on 5"<<endl;
-      ofs_sbp<<"ls 2 on 2"<<endl;
-      ofs_sbp<<"ls 2 on 4"<<endl;
-      ofs_sbp<<"ls 2 on 5"<<endl;
-      ofs_sbp<<"line on 7"<<endl;
-      ofs_sbp<<"ls 2 on 7"<<endl;
-
-      ofs_sbp<<"ma 1 on 2"<<endl;
-      ofs_sbp<<"color 1 on 1"<<endl;
-      ofs_sbp<<"color 2 on 2"<<endl;
-      ofs_sbp<<"color 3 on 3"<<endl;
-      ofs_sbp<<"color 4 on 4"<<endl;
-      ofs_sbp<<"color 5 on 5"<<endl;
-
-      ofs_sbp<<"win 1"<<endl;
-      ofs_sbp<<"yplot 1 2 3 4 5"<<endl;
-      ofs_sbp<<"loc 0 0 1 1"<<endl;
-      ofs_sbp<<"vie .1 .4 .9 .9"<<endl;
-      ofs_sbp<<"la y cnt/s/pixel/cm^2"<<endl;
-      ofs_sbp<<"log x"<<endl;
-      ofs_sbp<<"log y"<<endl;
-      ofs_sbp<<"r x "<<(radii[1]+radii[0])/2*cm_per_pixel/kpc<<" "<<(radii[sbps.size()-2]+radii[sbps.size()-1])/2*cm_per_pixel/kpc<<endl;
-      ofs_sbp<<"win 2"<<endl;
-      ofs_sbp<<"yplot 6 7"<<endl;
-      ofs_sbp<<"loc 0 0 1 1"<<endl;
-      ofs_sbp<<"vie .1 .1 .9 .4"<<endl;
-      ofs_sbp<<"la x radius (kpc)"<<endl;
-      ofs_sbp<<"la y chi"<<endl;
-      ofs_sbp<<"log y off"<<endl;
-      ofs_sbp<<"log x"<<endl;
-      ofs_sbp<<"r x "<<(radii[1]+radii[0])/2*cm_per_pixel/kpc<<" "<<(radii[sbps.size()-2]+radii[sbps.size()-1])/2*cm_per_pixel/kpc<<endl;
-    }
-  else
-    {
-      ofs_sbp<<"line on 2"<<endl;
-      ofs_sbp<<"line on 3"<<endl;
-      ofs_sbp<<"line on 4"<<endl;
-      ofs_sbp<<"ls 2 on 3"<<endl;
-      ofs_sbp<<"ls 2 on 4"<<endl;
-      ofs_sbp<<"line on 6"<<endl;
-      ofs_sbp<<"ls 2 on 6"<<endl;
-
-      ofs_sbp<<"color 1 on 1"<<endl;
-      ofs_sbp<<"color 3 on 2"<<endl;
-      ofs_sbp<<"color 4 on 3"<<endl;
-      ofs_sbp<<"color 5 on 4"<<endl;
-      //ofs_sbp<<"ma 1 on 2"<<endl;
-
-      ofs_sbp<<"win 1"<<endl;
-      ofs_sbp<<"yplot 1 2 3 4"<<endl;
-      ofs_sbp<<"loc 0 0 1 1"<<endl;
-      ofs_sbp<<"vie .1 .4 .9 .9"<<endl;
-      ofs_sbp<<"la y cnt/s/pixel/cm^2"<<endl;
-      ofs_sbp<<"log x"<<endl;
-      ofs_sbp<<"log y"<<endl;
-      ofs_sbp<<"r x "<<(radii[1]+radii[0])/2*cm_per_pixel/kpc<<" "<<(radii[radii.size()-2]+radii[radii.size()-1])/2*cm_per_pixel/kpc<<endl;
-      ofs_sbp<<"win 2"<<endl;
-      ofs_sbp<<"yplot 5 6"<<endl;
-      ofs_sbp<<"loc 0 0 1 1"<<endl;
-      ofs_sbp<<"vie .1 .1 .9 .4"<<endl;
-      ofs_sbp<<"la x radius (kpc)"<<endl;
-      ofs_sbp<<"la y chi"<<endl;
-      ofs_sbp<<"log x"<<endl;
-      ofs_sbp<<"log y off"<<endl;
-      ofs_sbp<<"r x "<<(radii[1]+radii[0])/2*cm_per_pixel/kpc<<" "<<(radii[radii.size()-2]+radii[radii.size()-1])/2*cm_per_pixel/kpc<<endl;
-
-    }
-  // cout<<sbps_all.size()<<"\t"<<sbps.size()<<"\t"<<sbps_inner_cut_size<<endl;
-  for(size_t i=1;i<sbps_all.size();++i)
-    {
-      double x=(radii_all[i]+radii_all[i-1])/2;
-      double y=sbps_all[i-1];
-      double ye=sbpe_all[i-1];
-      ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<y<<"\t"<<ye<<endl;
-    }
-  if(sbps_inner_cut_size>=1)
-    {
-      ofs_sbp<<"no no no"<<endl;
-      for(int i=1;i<sbps_inner_cut_size+1;++i)
-	{
-	  double x=(radii_all[i]+radii_all[i-1])/2;
-	  double ym=mv[i-1];
-	  ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<ym<<"\t"<<"0"<<endl;
-	}
-    }
-  ofs_sbp<<"no no no"<<endl;
-  for(size_t i=sbps_inner_cut_size;i<sbps_all.size();++i)
-    {
-      double x=(radii_all[i]+radii_all[i-1])/2;
-      double ym=mv[i-1];
-      ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<ym<<"\t"<<"0"<<endl;
-    }
-  ofs_sbp<<"no no no"<<endl;
-  //bkg level
-  bkg_level=abs(f.get_param_value("bkg"));
-  for(size_t i=0;i<sbps_all.size();++i)
-    {
-      double x=(radii_all[i]+radii_all[i-1])/2;
-      ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<bkg_level<<"\t0"<<endl;
-    }
-  ofs_sbp<<"no no no"<<endl;
-  //rc
-  double rc_kpc=abs(f.get_param_value("rc")*cm_per_pixel/kpc);
-  double max_sbp=*max_element(sbps_all.begin(),sbps_all.end());
-  double min_sbp=*min_element(sbps_all.begin(),sbps_all.end());
-  for(double x=min_sbp;x<=max_sbp;x+=(max_sbp-min_sbp)/100)
-    {
-      ofs_sbp<<rc_kpc<<"\t"<<x<<"\t"<<"0"<<endl;
-    }
-  //resid
-  ofs_sbp<<"no no no"<<endl;
-  for(size_t i=1;i<sbps.size();++i)
-    {
-      double x=(radii[i]+radii[i-1])/2;
-      //double y=sbps[i-1];
-      //double ye=sbpe[i-1];
-      double ym=mv[i-1];
-      ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<(ym-sbps[i-1])/sbpe[i-1]<<"\t"<<1<<endl;
-    }
-
-  //zero level of resid
-  ofs_sbp<<"no no no"<<endl;
-  for(size_t i=1;i<sbps.size();++i)
-    {
-      double x=(radii[i]+radii[i-1])/2;
-      //double y=sbps[i-1];
-      //double ye=sbpe[i-1];
-      //double ym=mv[i-1];
-      ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<0<<"\t"<<0<<endl;
-    }
-
-  mv=betao.eval(radii,p);
-  ofstream ofs_rho("rho_fit.qdp");
-  ofstream ofs_rho_data("rho_fit.dat");
-  ofstream ofs_entropy("entropy.qdp");
-  ofs_rho<<"la x radius (kpc)"<<endl;
-  ofs_rho<<"la y density (cm\\u-3\\d)"<<endl;
-  /*
-  for(int i=1;i<sbps.size();++i)
-    {
-      double x=(radii[i]+radii[i-1])/2;
-      double ym=mv[i-1];
-      ofs_rho<<x*cm_per_pixel/kpc<<"\t"<<ym<<endl;
-    }
-  */
-
-  //double lower,upper;
-  double dr=1;
-  //calculate the mass profile
-  //const double G=6.673E-8;//cm^3 g^-1 s^-2
-  // Molecular weight per electron
-  // Reference: Ettori et al. 2013, Space Sci. Rev., 177, 119-154; Eq.(9) below
-  static const double mu=1.155;
-  static const double mp=1.67262158E-24;//g
-  static const double M_sun=1.98892E33;//g
-  //static const double k=1.38E-16;
-
-  ofstream ofs_mass("mass_int.qdp");
-  ofstream ofs_mass_dat("mass_int.dat");
-  ofstream ofs_overdensity("overdensity.qdp");
-  ofstream ofs_gas_mass("gas_mass_int.qdp");
-  //ofs_mass<<"la x radius (kpc)"<<endl;
-  //ofs_mass<<"la y mass enclosed (solar mass)"<<endl;
-  //ofs_overdensity<<"la x radius (kpc)"<<endl;
-  //ofs_overdensity<<"la y overdensity"<<endl;
-  double gas_mass=0;
-  for(double r=1;r<200000;r+=dr)
-    {
-      dr=r/100;
-      double r1=r+dr;
-      double r_cm=r*cm_per_pixel;
-      double r1_cm=r1*cm_per_pixel;
-      double dr_cm=dr*cm_per_pixel;
-      double V_cm3=4./3.*pi*(dr_cm*(r1_cm*r1_cm+r_cm*r_cm+r_cm*r1_cm));
-      double ne=beta_func(r,n0,rc,beta);//cm^-3
-
-      double dmgas=V_cm3*ne*mu*mp/M_sun;
-      gas_mass+=dmgas;
-
-      ofs_gas_mass<<r*cm_per_pixel/kpc<<"\t"<<gas_mass<<endl;
-      double ne1=beta_func(r1,n0,rc,beta);//cm^3
-
-      double T_keV=Tprof(r);
-      double T1_keV=Tprof(r1);
-
-      //double T_K=T_keV*11604505.9;
-      //double T1_K=T1_keV*11604505.9;
-
-      double dlnT=log(T1_keV/T_keV);
-      double dlnr=log(r+dr)-log(r);
-      double dlnn=log(ne1/ne);
-
-      //double r_kpc=r_cm/kpc;
-      double r_Mpc=r_cm/Mpc;
-      //double M=-r_cm*T_K*k/G/mu/mp*(dlnT/dlnr+dlnn/dlnr);
-      //ref:http://adsabs.harvard.edu/abs/2012MNRAS.422.3503W
-      //Walker et al. 2012
-      double M=-3.68E13*M_sun*T_keV*r_Mpc*(dlnT/dlnr+dlnn/dlnr);
-      double rho=M/(4./3.*pi*r_cm*r_cm*r_cm);
-
-      double S=T_keV/pow(ne,2./3.);
-      //cout<<r<<"\t"<<M/M_sun<<endl;
-      //cout<<r<<"\t"<<T_keV<<endl;
-
-      ofs_rho<<r*cm_per_pixel/kpc<<"\t"<<ne<<endl;
-      ofs_rho_data<<r*cm_per_pixel/kpc<<"\t"<<ne<<endl;
-      ofs_entropy<<r*cm_per_pixel/kpc<<"\t"<<S<<endl;
-#if 0
-      if(r*cm_per_pixel/kpc<5)
-	{
-	  continue;
-	}
-#endif
-      ofs_mass<<r*cm_per_pixel/kpc<<"\t"<<M/M_sun<<endl;
-      if(r<radii.at(sbps.size()))
-	{
-	  ofs_mass_dat<<r*cm_per_pixel/kpc<<"\t0\t"<<M/M_sun<<"\t"<<M/M_sun*.1<<endl;
-	}
-      ofs_overdensity<<r*cm_per_pixel/kpc<<"\t"<<rho/calc_critical_density(z)<<endl;
-
-    }
-}
diff --git a/mass_profile/fit_dbeta_sbp.cpp b/mass_profile/fit_dbeta_sbp.cpp
deleted file mode 100644
index 71b3089..0000000
--- a/mass_profile/fit_dbeta_sbp.cpp
+++ /dev/null
@@ -1,586 +0,0 @@
-/*
-  Perform a double-beta density model fitting to the surface brightness data
-  Author: Junhua Gu
-  Last modified: 2016.06.07
-  This code is distributed with no warrant
-*/
-
-
-#include <iostream>
-#include <fstream>
-#include <sstream>
-#include <list>
-#include "vchisq.hpp"
-#include "dbeta.hpp"
-#include "beta_cfg.hpp"
-#include <data_sets/default_data_set.hpp>
-#include <methods/powell/powell_method.hpp>
-#include <core/freeze_param.hpp>
-#include <error_estimator/error_estimator.hpp>
-#include "spline.hpp"
-
-using namespace std;
-using namespace opt_utilities;
-//double s=5.63136645E20;
-const double kpc=3.086E21;//kpc in cm
-const double Mpc=kpc*1000;
-
-double dbeta_func(double r, double n01, double rc1, double beta1,
-                  double n02, double rc2, double beta2)
-{
-  double v1 = abs(n01) * pow(1+r*r/rc1/rc1, -3./2.*abs(beta1));
-  double v2 = abs(n02) * pow(1+r*r/rc2/rc2, -3./2.*abs(beta2));
-  return v1 + v2;
-}
-
-//calculate critical density from z, under following cosmological constants
-static double calc_critical_density(double z,
-				    const double H0=2.3E-18,
-				      const double Omega_m=.27)
-{
-  const double G=6.673E-8;//cm^3 g^-1 s^2
-  const double E=std::sqrt(Omega_m*(1+z)*(1+z)*(1+z)+1-Omega_m);
-  const double H=H0*E;
-  return 3*H*H/8/pi/G;
-}
-
-
-//A class enclosing the spline interpolation method
-class spline_func_obj
-  :public func_obj<double,double>
-{
-  //has an spline object
-  spline<double> spl;
-public:
-  //This function is used to calculate the intepolated value
-  double do_eval(const double& x)
-  {
-    return spl.get_value(x);
-  }
-
-  //we need this function, when this object is performing a clone of itself
-  spline_func_obj* do_clone()const
-  {
-    return new spline_func_obj(*this);
-  }
-
-public:
-  //add points to the spline object, after which the spline will be initialized
-  void add_point(double x,double y)
-  {
-    spl.push_point(x,y);
-  }
-
-  //before getting the intepolated value, the spline should be initialzied by calling this function
-  void gen_spline()
-  {
-    spl.gen_spline(0,0);
-  }
-};
-
-int main(int argc,char* argv[])
-{
-  if(argc!=2)
-    {
-      cerr<<argv[0]<<" <configure file>"<<endl;
-      return -1;
-    }
-  //initialize the parameters list
-  ifstream cfg_file(argv[1]);
-  assert(cfg_file.is_open());
-  cfg_map cfg=parse_cfg_file(cfg_file);
-
-  const double z=cfg.z;
-
-  //initialize the radius list, sbp list and sbp error list
-  std::vector<double> radii;
-  std::vector<double> sbps;
-  std::vector<double> sbpe;
-  std::vector<double> radii_all;
-  std::vector<double> sbps_all;
-  std::vector<double> sbpe_all;
-  //prepend the zero point to radius list
-  radii.push_back(0.0);
-  radii_all.push_back(0.0);
-  //read sbp and sbp error data
-  ifstream ifs(cfg.sbp_data.c_str());
-  std::string line;
-  if (ifs.is_open())
-    {
-      while(std::getline(ifs, line))
-        {
-          if (line.empty())
-            continue;
-
-          std::istringstream iss(line);
-          double x, xe, y, ye;
-          if ((iss >> x >> xe >> y >> ye))
-            {
-              std::cerr << "sbprofile data: "
-                        << x << ", " << xe << ", " << y << ", " << ye
-                        << std::endl;
-              radii.push_back(x+xe);  /* NOTE: use outer radii of regions */
-              radii_all.push_back(x+xe);
-              sbps.push_back(y);
-              sbps_all.push_back(y);
-              sbpe.push_back(ye);
-              sbpe_all.push_back(ye);
-            }
-          else
-            {
-              std::cerr << "skipped line: " << line << std::endl;
-            }
-        }
-    }
-  else
-    {
-      std::cerr << "ERROR: cannot open file: " << cfg.sbp_data.c_str()
-                << std::endl;
-      return 1;
-    }
-
-  //initialize the cm/pixel value
-  double cm_per_pixel=cfg.cm_per_pixel;
-  double rmin;
-  if(cfg.rmin_pixel>0)
-    {
-      rmin=cfg.rmin_pixel;
-    }
-  else
-    {
-      rmin=cfg.rmin_kpc*kpc/cm_per_pixel;
-    }
-
-  cerr<<"rmin="<<rmin<<" (pixel)"<<endl;
-  std::list<double> radii_tmp,sbps_tmp,sbpe_tmp;
-  radii_tmp.resize(radii.size());
-  sbps_tmp.resize(sbps.size());
-  sbpe_tmp.resize(sbpe.size());
-  copy(radii.begin(),radii.end(),radii_tmp.begin());
-  copy(sbps.begin(),sbps.end(),sbps_tmp.begin());
-  copy(sbpe.begin(),sbpe.end(),sbpe_tmp.begin());
-  for(list<double>::iterator i=radii_tmp.begin();i!=radii_tmp.end();)
-    {
-      if(*i<rmin)
-	{
-	  radii_tmp.pop_front();
-	  sbps_tmp.pop_front();
-	  sbpe_tmp.pop_front();
-	  i=radii_tmp.begin();
-	  continue;
-	}
-      ++i;
-    }
-  radii.resize(radii_tmp.size());
-  sbps.resize(sbps_tmp.size());
-  sbpe.resize(sbpe_tmp.size());
-  copy(radii_tmp.begin(),radii_tmp.end(),radii.begin());
-  copy(sbps_tmp.begin(),sbps_tmp.end(),sbps.begin());
-  copy(sbpe_tmp.begin(),sbpe_tmp.end(),sbpe.begin());
-
-  //read cooling function data
-  spline_func_obj cf;
-  for(ifstream ifs(cfg.cfunc_profile.c_str());;)
-    {
-      assert(ifs.is_open());
-      double x,y;
-      ifs>>x>>y;
-      if(!ifs.good())
-	{
-	  break;
-	}
-      cerr<<x<<"\t"<<y<<endl;
-      if(x>radii.back())
-	{
-	  break;
-	}
-      cf.add_point(x,y);
-    }
-  cf.gen_spline();
-
-  //read temperature profile data
-  spline_func_obj Tprof;
-  int tcnt=0;
-  for(ifstream ifs1(cfg.tprofile.c_str());;++tcnt)
-    {
-      assert(ifs1.is_open());
-      double x,y;
-      ifs1>>x>>y;
-      if(!ifs1.good())
-      {
-	break;
-      }
-      cerr<<x<<"\t"<<y<<endl;
-#if 0
-      if(tcnt==0)
-	{
-	  Tprof.add_point(0,y);
-	}
-#endif
-      Tprof.add_point(x,y);
-    }
-
-
-  Tprof.gen_spline();
-
-  default_data_set<std::vector<double>,std::vector<double> > ds;
-  ds.add_data(data<std::vector<double>,std::vector<double> >(radii,sbps,sbpe,sbpe,radii,radii));
-
-  //initial fitter
-  fitter<vector<double>,vector<double>,vector<double>,double> f;
-  f.load_data(ds);
-  //initial the object, which is used to calculate projection effect
-  projector<double> a;
-  bool tie_beta=false;
-  if(cfg.param_map.find("beta")!=cfg.param_map.end()
-     &&cfg.param_map.find("beta1")==cfg.param_map.end()
-     &&cfg.param_map.find("beta2")==cfg.param_map.end())
-    {
-      dbeta2<double> dbetao;
-      a.attach_model(dbetao);
-      tie_beta=true;
-    }
-  else if((cfg.param_map.find("beta1")!=cfg.param_map.end()
-	   ||cfg.param_map.find("beta2")!=cfg.param_map.end())
-	  &&cfg.param_map.find("beta")==cfg.param_map.end())
-    {
-      dbeta<double> dbetao;
-      a.attach_model(dbetao);
-      tie_beta=false;
-    }
-  else
-    {
-      cerr<<"Error, cannot decide whether to tie beta together or let them vary freely!"<<endl;
-      assert(0);
-    }
-
-  //attach the cooling function
-  a.attach_cfunc(cf);
-  a.set_cm_per_pixel(cm_per_pixel);
-
-  f.set_model(a);
-  //chi^2 statistic
-  vchisq<double> c;
-  c.verbose(true);
-  c.set_limit();
-  f.set_statistic(c);
-  //optimization method
-  f.set_opt_method(powell_method<double,std::vector<double> >());
-  //initialize the initial values
-  double n01=0;
-  double rc1=0;
-  double n02=0;
-  double rc2=0;
-  double beta=0;
-  double bkg_level=0;
-  if(tie_beta)
-    {
-      f.set_param_value("beta",.7);
-      f.set_param_lower_limit("beta",.3);
-      f.set_param_upper_limit("beta",1.4);
-    }
-  else
-    {
-      f.set_param_value("beta1",.7);
-      f.set_param_lower_limit("beta1",.3);
-      f.set_param_upper_limit("beta1",1.4);
-      f.set_param_value("beta2",.7);
-      f.set_param_lower_limit("beta2",.3);
-      f.set_param_upper_limit("beta2",1.4);
-    }
-  for(std::map<std::string,std::vector<double> >::iterator i=cfg.param_map.begin();
-      i!=cfg.param_map.end();++i)
-    {
-      std::string pname=i->first;
-      f.set_param_value(pname,i->second.at(0));
-      if(i->second.size()==3)
-	{
-	  double a1=i->second[1];
-	  double a2=i->second[2];
-	  double u=std::max(a1,a2);
-	  double l=std::min(a1,a2);
-	  f.set_param_upper_limit(pname,u);
-	  f.set_param_lower_limit(pname,l);
-	}
-      else
-	{
-	  if(pname=="beta"||pname=="beta1"||pname=="beta2")
-	    {
-	      f.set_param_lower_limit(pname,.3);
-	      f.set_param_upper_limit(pname,1.4);
-	    }
-	}
-    }
-
-
-
-  //perform the fitting, first freeze beta1, beta2, rc1, and rc2
-  if(tie_beta)
-    {
-      f.set_param_modifier(freeze_param<std::vector<double>,std::vector<double>,std::vector<double>,std::string>("beta")+
-			   freeze_param<std::vector<double>,std::vector<double>,std::vector<double>,std::string>("rc1")+
-			   freeze_param<std::vector<double>,std::vector<double>,std::vector<double>,std::string>("rc2")
-			   );
-    }
-  else
-    {
-      f.set_param_modifier(freeze_param<std::vector<double>,std::vector<double>,std::vector<double>,std::string>("beta1")+
-			   freeze_param<std::vector<double>,std::vector<double>,std::vector<double>,std::string>("beta2")+
-			   freeze_param<std::vector<double>,std::vector<double>,std::vector<double>,std::string>("rc1")+
-			   freeze_param<std::vector<double>,std::vector<double>,std::vector<double>,std::string>("rc2")
-			   );
-    }
-
-  f.fit();
-
-  f.clear_param_modifier();
-
-  //then perform the fitting, freeze beta1 and beta2
-  //f.set_param_modifier(freeze_param<std::vector<double>,std::vector<double>,std::vector<double>,std::string>("beta"));
-  //f.set_param_modifier(freeze_param<std::vector<double>,std::vector<double>,std::vector<double>,std::string>("bkg"));
-  f.fit();
-  //f.clear_param_modifier();
-
-  //finally thaw all parameters
-  f.fit();
-  double beta1=0;
-  double beta2=0;
-
-  n01=f.get_param_value("n01");
-  rc1=f.get_param_value("rc1");
-  n02=f.get_param_value("n02");
-  rc2=f.get_param_value("rc2");
-  if(tie_beta)
-    {
-      beta=f.get_param_value("beta");
-      beta1=beta;
-      beta2=beta;
-    }
-  else
-    {
-      beta1=f.get_param_value("beta1");
-      beta2=f.get_param_value("beta2");
-    }
-  //output the params
-  ofstream param_output("dbeta_param.txt");
-  //output the datasets and fitting results
-  for(size_t i=0;i<f.get_num_params();++i)
-    {
-      if(f.get_param_info(i).get_name()=="rc1")
-	{
-	  cerr<<"rc1_kpc"<<"\t"<<abs(f.get_param_info(i).get_value())*cm_per_pixel/kpc<<endl;
-	  param_output<<"rc1_kpc"<<"\t"<<abs(f.get_param_info(i).get_value())*cm_per_pixel/kpc<<endl;
-	}
-      if(f.get_param_info(i).get_name()=="rc2")
-	{
-	  cerr<<"rc2_kpc"<<"\t"<<abs(f.get_param_info(i).get_value())*cm_per_pixel/kpc<<endl;
-	  param_output<<"rc2_kpc"<<"\t"<<abs(f.get_param_info(i).get_value())*cm_per_pixel/kpc<<endl;
-	}
-      cerr<<f.get_param_info(i).get_name()<<"\t"<<abs(f.get_param_info(i).get_value())<<endl;
-      param_output<<f.get_param_info(i).get_name()<<"\t"<<abs(f.get_param_info(i).get_value())<<endl;
-    }
-  cerr<<"reduced_chi^2="<<f.get_statistic_value()/(radii.size()-f.get_model().get_num_free_params())<<endl;
-  param_output<<"reduced_chi^2="<<f.get_statistic_value()/(radii.size()-f.get_model().get_num_free_params())<<endl;
-
-  //c.verbose(false);
-  //f.set_statistic(c);
-  //f.fit();
-  std::vector<double> p=f.get_all_params();
-  f.clear_param_modifier();
-  std::vector<double> mv=f.eval_model(radii,p);
-
-
-  ofstream ofs_sbp("sbp_fit.qdp");
-  ofs_sbp<<"read serr 2"<<endl;
-  ofs_sbp<<"skip single"<<endl;
-
-  ofs_sbp<<"line on 2"<<endl;
-  ofs_sbp<<"line on 3"<<endl;
-  ofs_sbp<<"line on 4"<<endl;
-  ofs_sbp<<"line on 5"<<endl;
-  ofs_sbp<<"line on 7"<<endl;
-  ofs_sbp<<"ls 2 on 7"<<endl;
-  ofs_sbp<<"ls 2 on 3"<<endl;
-  ofs_sbp<<"ls 2 on 4"<<endl;
-  ofs_sbp<<"ls 2 on 5"<<endl;
-
-
-
-  ofs_sbp<<"!LAB  POS Y  4.00"<<endl;
-  ofs_sbp<<"!LAB  ROT"<<endl;
-  ofs_sbp<<"win 1"<<endl;
-  ofs_sbp<<"yplot 1 2 3 4 5"<<endl;
-  ofs_sbp<<"loc 0 0 1 1"<<endl;
-  ofs_sbp<<"vie .1 .4 .9 .9"<<endl;
-  ofs_sbp<<"la y cnt/s/pixel/cm^2"<<endl;
-  ofs_sbp<<"log x"<<endl;
-  ofs_sbp<<"log y"<<endl;
-  ofs_sbp<<"r x "<<(radii[1]+radii[0])/2*cm_per_pixel/kpc<<" "<<(radii[sbps.size()-2]+radii[sbps.size()-1])/2*cm_per_pixel/kpc<<endl;
-  ofs_sbp<<"win 2"<<endl;
-  ofs_sbp<<"yplot 6 7"<<endl;
-  ofs_sbp<<"loc 0 0 1 1"<<endl;
-  ofs_sbp<<"vie .1 .1 .9 .4"<<endl;
-  ofs_sbp<<"la x radius (kpc)"<<endl;
-  ofs_sbp<<"la y chi"<<endl;
-  ofs_sbp<<"log x"<<endl;
-  ofs_sbp<<"log y off"<<endl;
-  ofs_sbp<<"r x "<<(radii[1]+radii[0])/2*cm_per_pixel/kpc<<" "<<(radii[sbps.size()-2]+radii[sbps.size()-1])/2*cm_per_pixel/kpc<<endl;
-  for(size_t i=1;i<sbps.size();++i)
-    {
-      double x=(radii[i]+radii[i-1])/2;
-      double y=sbps[i-1];
-      double ye=sbpe[i-1];
-      //double ym=mv[i-1];
-      ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<y<<"\t"<<ye<<endl;
-    }
-  ofs_sbp<<"no no no"<<endl;
-  for(size_t i=1;i<sbps.size();++i)
-    {
-      double x=(radii[i]+radii[i-1])/2;
-      //double y=sbps[i-1];
-      //double ye=sbpe[i-1];
-      double ym=mv[i-1];
-      ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<ym<<"\t"<<0<<endl;
-    }
-  //bkg
-  ofs_sbp<<"no no no"<<endl;
-  bkg_level=abs(f.get_param_value("bkg"));
-  for(size_t i=0;i<sbps.size();++i)
-    {
-      double x=(radii[i]+radii[i-1])/2;
-      ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<bkg_level<<"\t0"<<endl;
-    }
-  //rc1
-  ofs_sbp<<"no no no"<<endl;
-  double rc1_kpc=abs(f.get_param_value("rc1")*cm_per_pixel/kpc);
-  double max_sbp=*max_element(sbps.begin(),sbps.end());
-  double min_sbp=*min_element(sbps.begin(),sbps.end());
-  for(double x=min_sbp;x<=max_sbp;x+=(max_sbp-min_sbp)/100)
-    {
-      ofs_sbp<<rc1_kpc<<"\t"<<x<<"\t"<<"0"<<endl;
-    }
-  //rc2
-  ofs_sbp<<"no no no"<<endl;
-  double rc2_kpc=abs(f.get_param_value("rc2")*cm_per_pixel/kpc);
-  for(double x=min_sbp;x<=max_sbp;x+=(max_sbp-min_sbp)/100)
-    {
-      ofs_sbp<<rc2_kpc<<"\t"<<x<<"\t"<<"0"<<endl;
-    }
-  //resid
-  ofs_sbp<<"no no no"<<endl;
-  for(size_t i=1;i<sbps.size();++i)
-    {
-      double x=(radii[i]+radii[i-1])/2;
-      //double y=sbps[i-1];
-      //double ye=sbpe[i-1];
-      double ym=mv[i-1];
-      ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<(ym-sbps[i-1])/sbpe[i-1]<<"\t"<<1<<endl;
-    }
-  //zero level in resid map
-  ofs_sbp<<"no no no"<<endl;
-  for(size_t i=1;i<sbps.size();++i)
-    {
-      double x=(radii[i]+radii[i-1])/2;
-      //double y=sbps[i-1];
-      //double ye=sbpe[i-1];
-      //double ym=mv[i-1];
-      ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<0<<"\t"<<0<<endl;
-    }
-
-
-  mv=f.eval_model_raw(radii,p);
-  ofstream ofs_rho("rho_fit.qdp");
-  ofstream ofs_rho_data("rho_fit.dat");
-  ofstream ofs_entropy("entropy.qdp");
-  ofs_rho<<"la x radius (kpc)"<<endl;
-  ofs_rho<<"la y density (cm\\u-3\\d)"<<endl;
-  /*
-  for(int i=1;i<sbps.size();++i)
-    {
-      double x=(radii[i]+radii[i-1])/2;
-      double ym=mv[i-1];
-      ofs_rho<<x*cm_per_pixel/kpc<<"\t"<<ym<<endl;
-    }
-  */
-
-  //double lower,upper;
-  double dr=1;
-  //calculate the mass profile
-  //const double G=6.673E-8;//cm^3 g^-1 s^-2
-  // Molecular weight per electron
-  // Reference: Ettori et al. 2013, Space Sci. Rev., 177, 119-154; Eq.(9) below
-  static const double mu=1.155;
-  static const double mp=1.67262158E-24;//g
-  static const double M_sun=1.98892E33;//g
-  //static const double k=1.38E-16;
-
-  ofstream ofs_mass("mass_int.qdp");
-  ofstream ofs_mass_dat("mass_int.dat");
-  ofstream ofs_overdensity("overdensity.qdp");
-  ofstream ofs_gas_mass("gas_mass_int.qdp");
-  //ofs_mass<<"la x radius (kpc)"<<endl;
-  //ofs_mass<<"la y mass enclosed (solar mass)"<<endl;
-  //ofs_overdensity<<"la x radius (kpc)"<<endl;
-  //ofs_overdensity<<"la y overdensity"<<endl;
-  double gas_mass=0;
-  for(double r=1;r<200000;r+=dr)
-    {
-      dr=r/100;
-      double r1=r+dr;
-      double r_cm=r*cm_per_pixel;
-      double r1_cm=r1*cm_per_pixel;
-      double dr_cm=dr*cm_per_pixel;
-      double V_cm3=4./3.*pi*(dr_cm*(r1_cm*r1_cm+r_cm*r_cm+r_cm*r1_cm));
-      double ne=dbeta_func(r,n01,rc1,beta1,
-			   n02,rc2,beta2);//cm^3
-
-      double dmgas=V_cm3*ne*mu*mp/M_sun;
-      gas_mass+=dmgas;
-
-      ofs_gas_mass<<r*cm_per_pixel/kpc<<"\t"<<gas_mass<<endl;
-      double ne_beta1=dbeta_func(r,n01,rc1,beta1, 0,rc2,beta2);
-
-      double ne_beta2=dbeta_func(r,0,rc1,beta1, n02,rc2,beta2);
-
-      double ne1=dbeta_func(r1,n01,rc1,beta1, n02,rc2,beta2);//cm^3
-
-      double T_keV=Tprof(r);
-      double T1_keV=Tprof(r1);
-
-      //double T_K=T_keV*11604505.9;
-      //double T1_K=T1_keV*11604505.9;
-
-      double dlnT=log(T1_keV/T_keV);
-      double dlnr=log(r+dr)-log(r);
-      double dlnn=log(ne1/ne);
-
-      //double r_kpc=r_cm/kpc;
-      double r_Mpc=r_cm/Mpc;
-      //double M=-r_cm*T_K*k/G/mu/mp*(dlnT/dlnr+dlnn/dlnr);
-      //ref:http://adsabs.harvard.edu/abs/2012MNRAS.422.3503W
-      //Walker et al. 2012
-      double M=-3.68E13*M_sun*T_keV*r_Mpc*(dlnT/dlnr+dlnn/dlnr);
-      double rho=M/(4./3.*pi*r_cm*r_cm*r_cm);
-
-      double S=T_keV/pow(ne,2./3.);
-      //cout<<r<<"\t"<<M/M_sun<<endl;
-      //cout<<r<<"\t"<<T_keV<<endl;
-
-      ofs_rho<<r*cm_per_pixel/kpc<<"\t"<<ne<<"\t"<<ne_beta1<<"\t"<<ne_beta2<<endl;
-      ofs_rho_data<<r*cm_per_pixel/kpc<<"\t"<<ne<<endl;
-      ofs_entropy<<r*cm_per_pixel/kpc<<"\t"<<S<<endl;
-#if 0
-      if(r*cm_per_pixel/kpc<5)
-	{
-	  continue;
-	}
-#endif
-      ofs_mass<<r*cm_per_pixel/kpc<<"\t"<<M/M_sun<<endl;
-      if(r<radii.back())
-	{
-	  ofs_mass_dat<<r*cm_per_pixel/kpc<<"\t0\t"<<M/M_sun<<"\t"<<M/M_sun*.1<<endl;
-	}
-      ofs_overdensity<<r*cm_per_pixel/kpc<<"\t"<<rho/calc_critical_density(z)<<endl;
-
-    }
-}
diff --git a/mass_profile/fit_nfw_mass.cpp b/mass_profile/fit_nfw_mass.cpp
deleted file mode 100644
index fae74ef..0000000
--- a/mass_profile/fit_nfw_mass.cpp
+++ /dev/null
@@ -1,159 +0,0 @@
-/*
-  Fitting nfw mass profile model
-  Author: Junhua Gu
-  Last modification 20120721
-*/
-
-#include "nfw.hpp"
-#include <core/optimizer.hpp>
-#include <core/fitter.hpp>
-#include <data_sets/default_data_set.hpp>
-#include "chisq.hpp"
-#include <methods/powell/powell_method.hpp>
-#include <iostream>
-#include <fstream>
-#include <vector>
-#include <string>
-
-using namespace opt_utilities;
-using namespace std;
-const double cm=1;
-const double kpc=3.08568e+21*cm;
-const double pi=4*atan(1);
-static double calc_critical_density(double z,
-				    const double H0=2.3E-18,
-				    const double Omega_m=.27)
-{
-  const double G=6.673E-8;//cm^3 g^-1 s^2
-  const double E=std::sqrt(Omega_m*(1+z)*(1+z)*(1+z)+1-Omega_m);
-  const double H=H0*E;
-  return 3*H*H/8/pi/G;
-}
-
-
-int main(int argc,char* argv[])
-{
-  if(argc<3)
-    {
-      cerr<<"Usage:"<<argv[0]<<" <data file with 4 columns of x, xe, y, ye> <z> [rmin in kpc]"<<endl;
-      return -1;
-    }
-  double rmin_kpc=1;
-  if(argc>=4)
-    {
-      rmin_kpc=atof(argv[3]);
-    }
-  double z=0;
-  z=atof(argv[2]);
-  //define the fitter
-  fitter<double,double,vector<double>,double,std::string> fit;
-  //define the data set
-  default_data_set<double,double> ds;
-  //open the data file
-  ifstream ifs(argv[1]);
-  //cout<<"read serr 2"<<endl;
-  ofstream ofs_fit_result("nfw_fit_result.qdp");
-
-  ofs_fit_result<<"read serr 1 2"<<endl;
-  ofs_fit_result<<"skip single"<<endl;
-  ofs_fit_result<<"line off"<<endl;
-  ofs_fit_result<<"li on 2"<<endl;
-  ofs_fit_result<<"li on 4"<<endl;
-  ofs_fit_result<<"ls 2 on 4"<<endl;
-
-  ofs_fit_result<<"win 1"<<endl;
-  ofs_fit_result<<"yplot 1 2"<<endl;
-  ofs_fit_result<<"loc 0 0 1 1"<<endl;
-  ofs_fit_result<<"vie .1 .4 .9 .9"<<endl;
-  ofs_fit_result<<"la y Mass (solar)"<<endl;
-  ofs_fit_result<<"log x"<<endl;
-  ofs_fit_result<<"log y"<<endl;
-  ofs_fit_result<<"win 2"<<endl;
-  ofs_fit_result<<"yplot 3 4"<<endl;
-  ofs_fit_result<<"loc 0 0 1 1"<<endl;
-  ofs_fit_result<<"vie .1 .1 .9 .4"<<endl;
-  ofs_fit_result<<"la x radius (kpc)"<<endl;
-  ofs_fit_result<<"la y chi"<<endl;
-  ofs_fit_result<<"log x"<<endl;
-  ofs_fit_result<<"log y off"<<endl;
-  for(;;)
-    {
-      //read radius, temperature and error
-      double r,re,m,me;
-      ifs>>r>>re>>m>>me;
-      if(!ifs.good())
-	{
-	  break;
-	}
-      if(r<rmin_kpc)
-	{
-	  continue;
-	}
-      data<double,double> d(r,m,me,me,re,re);
-      ofs_fit_result<<r<<"\t"<<re<<"\t"<<m<<"\t"<<me<<endl;
-      ds.add_data(d);
-    }
-  ofs_fit_result<<"no no no"<<endl;
-  //load data
-  fit.load_data(ds);
-  //define the optimization method
-  fit.set_opt_method(powell_method<double,vector<double> >());
-  //use chi^2 statistic
-  fit.set_statistic(chisq<double,double,vector<double>,double,std::string>());
-  fit.set_model(nfw<double>());
-  //fit.set_param_value("rs",4);
-  //fit.set_param_value("rho0",100);
-  fit.fit();
-  fit.fit();
-  vector<double> p=fit.fit();
-  //output parameters
-  ofstream ofs_param("nfw_param.txt");
-  for(size_t i=0;i<fit.get_num_params();++i)
-    {
-      cout<<fit.get_param_info(i).get_name()<<"\t"<<fit.get_param_info(i).get_value()<<endl;
-      ofs_param<<fit.get_param_info(i).get_name()<<"\t"<<fit.get_param_info(i).get_value()<<endl;
-    }
-  cout<<"reduced chi^2="<<fit.get_statistic_value()<<endl;
-  ofs_param<<"reduced chi^2="<<fit.get_statistic_value()<<endl;
-  ofstream ofs_model("nfw_dump.qdp");
-  ofstream ofs_overdensity("overdensity.qdp");
-  //const double G=6.673E-8;//cm^3 g^-1 s^-2
-  //static const double mu=1.4074;
-  //static const double mp=1.67262158E-24;//g
-  static const double M_sun=1.98892E33;//g
-  //static const double k=1.38E-16;
-
-  double xmax=0;
-  for(double x=std::max(rmin_kpc,ds.get_data(0).get_x());;x+=1)
-    {
-      double model_value=fit.eval_model(x,p);
-      ofs_model<<x<<"\t"<<model_value<<endl;
-      ofs_fit_result<<x<<"\t0\t"<<model_value<<"\t0"<<endl;
-      double V=4./3.*pi*pow(x*kpc,3);
-      double m=model_value*M_sun;
-      double rho=m/V;//g/cm^3
-      double over_density=rho/calc_critical_density(z);
-      ofs_overdensity<<x<<"\t"<<over_density<<endl;
-      xmax=x;
-      if(over_density<100)
-	{
-	  break;
-	}
-    }
-  ofs_fit_result<<"no no no"<<endl;
-  for(size_t i=0;i<ds.size();++i)
-    {
-      data<double,double> d=ds.get_data(i);
-      double x=d.get_x();
-      double y=d.get_y();
-      double ye=d.get_y_lower_err();
-      double ym=fit.eval_model(x,p);
-      ofs_fit_result<<x<<"\t"<<0<<"\t"<<(y-ym)/ye<<"\t"<<1<<endl;
-    }
-  ofs_fit_result<<"no no no"<<endl;
-  for(double x=std::max(rmin_kpc,ds.get_data(0).get_x());x<xmax;x+=1)
-    {
-      ofs_fit_result<<x<<"\t"<<0<<"\t"<<0<<"\t"<<0<<endl;
-    }
-
-}
diff --git a/mass_profile/fit_wang2012_model.cpp b/mass_profile/fit_wang2012_model.cpp
deleted file mode 100644
index f15a2b1..0000000
--- a/mass_profile/fit_wang2012_model.cpp
+++ /dev/null
@@ -1,195 +0,0 @@
-/*
-  Fitting Jy Wang's temperature profile model
-  Author: Jingying Wang
-  Last modification 20120819
-
-*/
-
-#include "wang2012_model.hpp"
-#include <core/optimizer.hpp>
-#include <core/fitter.hpp>
-#include <data_sets/default_data_set.hpp>
-#include "chisq.hpp"
-#include <methods/powell/powell_method.hpp>
-#include <core/freeze_param.hpp>
-#include <iostream>
-#include <fstream>
-#include <vector>
-#include <string>
-
-using namespace opt_utilities;
-using namespace std;
-const double cm=1;
-const double kpc=3.08568e+21*cm;
-
-int main(int argc,char* argv[])
-{
-  if(argc<2)
-    {
-      cerr<<"Usage:"<<argv[0]<<" <data file with 4 columns of x, xe, y, ye> [param file] [cm per pixel]"<<endl;
-      return -1;
-    }
-  double cm_per_pixel=-1;
-  if(argc>=4)
-    {
-      cm_per_pixel=atof(argv[3]);
-    }
-
-  //define the fitter
-  fitter<double,double,vector<double>,double,std::string> fit;
-  //define the data set
-  default_data_set<double,double> ds;
-  //open the data file
-  ifstream ifs(argv[1]);
-  double min_r=1e9;
-  //cout<<"read serr 2"<<endl;
-  ofstream ofs_fit_result("fit_result.qdp");
-  ofs_fit_result<<"read serr 1 2"<<endl;
-  ofs_fit_result<<"skip single"<<endl;
-  if(cm_per_pixel>0)
-    {
-      ofs_fit_result<<"la x radius (kpc)"<<endl;
-    }
-  else
-    {
-      ofs_fit_result<<"la x radius (pixel)"<<endl;
-    }
-  ofs_fit_result<<"la y temperature (keV)"<<endl;
-  for(;;)
-    {
-      //read radius, temperature and error
-      double r,re,t,te;
-      ifs>>r>>re>>t>>te;
-      if(!ifs.good())
-	{
-	  break;
-	}
-      min_r=min(r,min_r);
-      data<double,double> d(r,t,te,te,re,re);
-      //std::cerr<<r<<"\t"<<t<<"\t"<<te<<endl;
-      if(cm_per_pixel>0)
-	{
-	  ofs_fit_result<<r*cm_per_pixel/kpc<<"\t"<<re*cm_per_pixel/kpc<<"\t"<<t<<"\t"<<te<<endl;
-	}
-      else
-	{
-	  ofs_fit_result<<r<<"\t"<<re<<"\t"<<t<<"\t"<<te<<endl;
-	}
-      ds.add_data(d);
-    }
-  ofs_fit_result<<"no no no"<<endl;
-  //load data
-  fit.load_data(ds);
-  //define the optimization method
-  fit.set_opt_method(powell_method<double,vector<double> >());
-  //use chi^2 statistic
-  chisq<double,double,vector<double>,double,std::string> chisq_object;
-  chisq_object.set_limit();
-  fit.set_statistic(chisq_object);
-  //fit.set_statistic(chisq<double,double,vector<double>,double,std::string>());
-  fit.set_model(wang2012_model<double>());
-
-  if(argc>=3&&std::string(argv[2])!="NONE")
-    {
-      std::vector<std::string> freeze_list;
-      ifstream ifs_param(argv[2]);
-      assert(ifs_param.is_open());
-      for(;;)
-	{
-	  string pname;
-	  double pvalue;
-	  double lower,upper;
-	  char param_status;
-	  ifs_param>>pname>>pvalue>>lower>>upper>>param_status;
-	  if(!ifs_param.good())
-	    {
-	      break;
-	    }
-	  if(param_status=='F')
-	    {
-	      freeze_list.push_back(pname);
-	    }
-	  if(pvalue<=lower||pvalue>=upper)
-	    {
-	      cerr<<"Invalid initial value, central value not enclosed by the lower and upper boundaries, adjust automatically"<<endl;
-	      pvalue=std::max(pvalue,lower);
-	      pvalue=std::min(pvalue,upper);
-	    }
-	  fit.set_param_value(pname,pvalue);
-	  fit.set_param_lower_limit(pname,lower);
-	  fit.set_param_upper_limit(pname,upper);
-	}
-      if(!freeze_list.empty())
-	{
-	  freeze_param<double,double,std::vector<double>,std::string> fp(freeze_list[0]);
-	  fit.set_param_modifier(fp);
-	  for(size_t i=1;i<freeze_list.size();++i)
-	    {
-	      dynamic_cast<freeze_param<double,double,std::vector<double>,std::string>&>(fit.get_param_modifier())+=freeze_param<double,double,std::vector<double>,std::string>(freeze_list[i]);
-	    }
-	}
-    }
-
-  for(int i=0;i<100;++i)
-    {
-      fit.fit();
-    }
-  vector<double> p=fit.fit();
-#if 0
-  ofstream output_param;
-  if(argc>=3&&std::string(argv[2])!="NONE")
-    {
-      output_param.open(argv[2]);
-    }
-  else
-    {
-      output_param.open("para0.txt");
-    }
-#endif
-  //output parameters
-  for(size_t i=0;i<fit.get_num_params();++i)
-    {
-      std::string pname=fit.get_param_info(i).get_name();
-      std::string pstatus=fit.get_model().report_param_status(pname);
-
-      if(pstatus==""||pstatus=="thawed")
-	{
-	  pstatus="T";
-	}
-      else
-	{
-	  pstatus="F";
-	}
-      cout<<fit.get_param_info(i).get_name()<<"\t"<<fit.get_param_info(i).get_value()<<"\t"<<fit.get_param_info(i).get_lower_limit()<<"\t"<<fit.get_param_info(i).get_upper_limit()<<"\t"<<pstatus<<endl;
-      //if(argc>=3&&std::string(argv[2])!="NONE")
-#if 0
-	{
-	  output_param<<fit.get_param_info(i).get_name()<<"\t"<<fit.get_param_info(i).get_value()<<"\t"<<fit.get_param_info(i).get_lower_limit()<<"\t"<<fit.get_param_info(i).get_upper_limit()<<"\t"<<pstatus<<endl;
-	}
-#endif
-    }
-
-  //cout<<"T0"<<"\t"<<fit.get_param_value("T0")<<endl;
-  //cout<<"T1"<<"\t"<<fit.get_param_value("T1")<<endl;
-  //cout<<"xt"<<"\t"<<fit.get_param_value("xt")<<endl;
-  //cout<<"eta"<<"\t"<<fit.get_param_value("eta")<<endl;
-  //dump the data for checking
-  ofstream ofs_model("wang2012_dump.qdp");
-
-
-  min_r=0;
-  for(double x=min_r;x<3000;x+=10)
-    {
-      double model_value=fit.eval_model_raw(x,p);
-
-      ofs_model<<x<<"\t"<<model_value<<endl;
-      if(cm_per_pixel>0)
-	{
-	  ofs_fit_result<<x*cm_per_pixel/kpc<<"\t0\t"<<model_value<<"\t0"<<endl;
-	}
-      else
-	{
-	  ofs_fit_result<<x<<"\t0\t"<<model_value<<"\t0"<<endl;
-	}
-    }
-}
diff --git a/mass_profile/nfw.hpp b/mass_profile/nfw.hpp
deleted file mode 100644
index 8e3e59e..0000000
--- a/mass_profile/nfw.hpp
+++ /dev/null
@@ -1,56 +0,0 @@
-/**
-   \file nfw.hpp
-   \brief Jingying Wang's model
-   \author Jingying Wang
- */
-
-
-#ifndef NFW
-#define NFW
-#define OPT_HEADER
-#include <core/fitter.hpp>
-#include <cmath>
-
-namespace opt_utilities
-{
-  template <typename T>
-  class nfw
-    :public model<T,T,std::vector<T>,std::string>
-  {
-  private:
-    model<T,T,std::vector<T> >* do_clone()const
-    {
-      return new nfw<T>(*this);
-    }
-
-    const char* do_get_type_name()const
-    {
-      return "1d power law";
-    }
-  public:
-    nfw()
-    {
-      this->push_param_info(param_info<std::vector<T> >("rho0",1,0,1e99));
-      this->push_param_info(param_info<std::vector<T> >("rs",100,0,1e99));
-    }
-
-    T do_eval(const T& r,const std::vector<T>& param)
-    {
-      T rho0=std::abs(param[0]);
-      T rs=std::abs(param[1]);
-      static const T pi=4*std::atan(1);
-      return 4*pi*rho0*rs*rs*rs*(std::log((r+rs)/rs)-r/(r+rs));
-    }
-
-  private:
-    std::string do_get_information()const
-    {
-      return "";
-    }
-  };
-}
-
-
-
-#endif
-//EOF
diff --git a/mass_profile/projector.hpp b/mass_profile/projector.hpp
deleted file mode 100644
index 4ef4b32..0000000
--- a/mass_profile/projector.hpp
+++ /dev/null
@@ -1,214 +0,0 @@
-#ifndef PROJ_HPP
-#define PROJ_HPP
-/*
-  Defining the class that is used to consider the projection effect
-  Author: Junhua Gu
-  Last modified: 2011.01.01
-*/
-
-
-#include <core/fitter.hpp>
-#include <vector>
-#include <cmath>
-
-static const double pi=4*atan(1);
-// Ratio of the electron density (n_e) to the proton density (n_p)
-//   n_gas = n_e + n_p ~= 1.826 n_e  =>  n_e / n_p ~= 1.21
-// Reference: Ettori et al. 2013, Space Sci. Rev., 177, 119-154; Eq.(9) below
-static const double ne_np_ratio = 1.21;
-
-namespace opt_utilities
-{
-  //This is used to project a 3-D surface brightness model to 2-D profile
-  template <typename T>
-  class projector
-    :public model<std::vector<T>,std::vector<T>,std::vector<T> >
-  {
-  private:
-    //Points to a 3-D model that is to be projected
-    model<std::vector<T>,std::vector<T>,std::vector<T> >* pmodel;
-    func_obj<T,T>* pcfunc;
-    T cm_per_pixel;
-  public:
-    //default cstr
-    projector()
-      :pmodel(NULL_PTR),pcfunc(NULL_PTR),cm_per_pixel(1)
-    {}
-    //copy cstr
-    projector(const projector& rhs)
-      :cm_per_pixel(rhs.cm_per_pixel)
-    {
-      attach_model(*(rhs.pmodel));
-      if(rhs.pcfunc)
-        {
-          pcfunc=rhs.pcfunc->clone();
-        }
-      else
-        {
-          pcfunc=NULL_PTR;
-        }
-    }
-    //assign operator
-    projector& operator=(const projector& rhs)
-    {
-      cm_per_pixel=rhs.cm_per_pixel;
-      if(pmodel)
-        {
-          pmodel->destroy();
-        }
-      if(pcfunc)
-        {
-          pcfunc->destroy();
-        }
-      if(rhs.pcfunc)
-        {
-          pcfunc=rhs.pcfunc->clone();
-        }
-      if(rhs.pmodel)
-        {
-          pmodel=rhs.pmodel->clone();
-        }
-    }
-    //destr
-    ~projector()
-    {
-      if(pmodel)
-        {
-          pmodel->destroy();
-        }
-      if(pcfunc)
-        {
-          pcfunc->destroy();
-        }
-    }
-    //used to clone self
-    model<std::vector<T>,std::vector<T>,std::vector<T> >*
-    do_clone()const
-    {
-      return new projector(*this);
-    }
-
-  public:
-    void set_cm_per_pixel(const T& x)
-    {
-      cm_per_pixel=x;
-    }
-
-    //attach the model that is to be projected
-    void attach_model(const model<std::vector<T>,std::vector<T>,std::vector<T> >& m)
-    {
-      this->clear_param_info();
-      for(size_t i=0;i<m.get_num_params();++i)
-        {
-          this->push_param_info(m.get_param_info(i));
-        }
-      this -> push_param_info(param_info<std::vector<T>,
-                              std::string>("bkg",0,0,1E99));
-      pmodel=m.clone();
-      pmodel->clear_param_modifier();
-    }
-
-    void attach_cfunc(const func_obj<T,T>& cf)
-    {
-      if(pcfunc)
-        {
-          pcfunc->destroy();
-        }
-      pcfunc=cf.clone();
-    }
-
-  public:
-    //calc the volume
-    /*
-      This is a sphere that is subtracted by a cycline.
-       /|     |\
-      / |     | \
-      | |     | |
-      | |     | |
-      \ |     | /
-       \|     |/
-     */
-    T calc_v_ring(T rsph,T rcyc)
-    {
-      if(rcyc<rsph)
-        {
-          double a=rsph*rsph-rcyc*rcyc;
-          return 4.*pi/3.*std::sqrt(a*a*a);
-        }
-      return 0;
-    }
-
-    //calc the No. nsph sphere's projection volume on the No. nrad pie region
-    T calc_v(const std::vector<T>& rlist,int nsph,int nrad)
-    {
-      if(nsph<nrad)
-        {
-          return 0;
-        }
-      else if(nsph==nrad)
-        {
-          return calc_v_ring(rlist[nsph+1], rlist[nrad]);
-        }
-      else {
-        return (calc_v_ring(rlist[nsph+1], rlist[nrad]) -
-                calc_v_ring(rlist[nsph], rlist[nrad]) -
-                calc_v_ring(rlist[nsph+1], rlist[nrad+1]) +
-                calc_v_ring(rlist[nsph], rlist[nrad+1]));
-      }
-    }
-  public:
-    bool do_meets_constraint(const std::vector<T>& p)const
-    {
-      std::vector<T> p1(this->reform_param(p));
-      for(size_t i=0;i!=p1.size();++i)
-        {
-          if(get_element(p1,i)>this->get_param_info(i).get_upper_limit()||
-            get_element(p1,i)<this->get_param_info(i).get_lower_limit())
-            {
-              // std::cerr<<this->get_param_info(i).get_name()<<"\t"<<p1[i]<<std::endl;
-              return false;
-            }
-        }
-      std::vector<T> p2(p1.size()-1);
-      for(size_t i=0;i<p1.size()-1;++i)
-        {
-          p2.at(i)=p1[i];
-        }
-
-      return pmodel->meets_constraint(p2);
-    }
-  public:
-    //Perform the projection
-    std::vector<T> do_eval(const std::vector<T>& x,const std::vector<T>& p)
-    {
-      T bkg=std::abs(p.back());
-      //I think following codes are clear enough :).
-      std::vector<T> unprojected(pmodel->eval(x,p));
-      std::vector<T> projected(unprojected.size());
-
-      for(size_t nrad=0; nrad<x.size()-1; ++nrad)
-        {
-          for(size_t nsph=nrad; nsph<x.size()-1; ++nsph)
-            {
-              double v = calc_v(x, nsph, nrad) * pow(cm_per_pixel, 3);
-              if(pcfunc)
-                {
-                  double cfunc = (*pcfunc)((x[nsph+1] + x[nsph]) / 2.0);
-                  projected[nrad] += (unprojected[nsph] * unprojected[nsph] *
-                                      cfunc * v / ne_np_ratio);
-                }
-              else
-                {
-                  projected[nrad] += unprojected[nsph] * unprojected[nsph] * v;
-                }
-            }
-          double area = pi * (x[nrad+1]*x[nrad+1] - x[nrad]*x[nrad]);
-          projected[nrad] /= area;
-          projected[nrad] += bkg;
-        }
-      return projected;
-    }
-  };
-};
-
-#endif
diff --git a/mass_profile/report_error.cpp b/mass_profile/report_error.cpp
deleted file mode 100644
index 23e6a26..0000000
--- a/mass_profile/report_error.cpp
+++ /dev/null
@@ -1,39 +0,0 @@
-#include <iostream>
-
-using namespace std;
-void report_error(const char* message)
-{
-  cerr<<"MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM"<<endl;
-  cerr<<"MMMMMMMMMMMM        MMMMMMMMMMMM"<<endl;
-  cerr<<"MMMMMMMMMM            MMMMMMMMMM"<<endl;
-  cerr<<"MMMMMMMMM              MMMMMMMMM"<<endl;
-  cerr<<"MMMMMMMM                MMMMMMMM"<<endl;
-  cerr<<"MMMMMMM                 MMMMMMMM"<<endl;
-  cerr<<"MMMMMMM                  MMMMMMM"<<endl;
-  cerr<<"MMMMMMM                  MMMMMMM"<<endl;
-  cerr<<"MMMMMMM    MMM    MMM    MMMMMMM"<<endl;
-  cerr<<"MMMMMMM   MMMMM   MMMM   MMMMMMM"<<endl;
-  cerr<<"MMMMMMM   MMMMM   MMMM   MMMMMMM"<<endl;
-  cerr<<"MMMMMMMM   MMMM M MMMM  MMMMMMMM"<<endl;
-  cerr<<"MMMMMMMM        M        MMMMMMM"<<endl;
-  cerr<<"MMMMMMMM       MMM      MMMMMMMM"<<endl;
-  cerr<<"MMMMMMMMMMMM   MMM  MMMMMMMMMMMM"<<endl;
-  cerr<<"MMMMMMMMMM MM       M  MMMMMMMMM"<<endl;
-  cerr<<"MMMMMMMMMM  M M M M M MMMMMMMMMM"<<endl;
-  cerr<<"MMMM  MMMMM MMMMMMMMM MMMMM   MM"<<endl;
-  cerr<<"MMM    MMMM M MMMMM M MMMM    MM"<<endl;
-  cerr<<"MMM    MMMM   M M M  MMMMM   MMM"<<endl;
-  cerr<<"MMMM    MMMM         MMM      MM"<<endl;
-  cerr<<"MMM       MMMM     MMMM       MM"<<endl;
-  cerr<<"MMM         MMMMMMMM      M  MMM"<<endl;
-  cerr<<"MMMM  MMM      MMM      MMMMMMMM"<<endl;
-  cerr<<"MMMMMMMMMMM  MM       MMMMMMM  M"<<endl;
-  cerr<<"MMM  MMMMMMM       MMMMMMMMM   M"<<endl;
-  cerr<<"MM    MMM        MM            M"<<endl;
-  cerr<<"MM            MMMM            MM"<<endl;
-  cerr<<"MMM        MMMMMMMMMMMMM       M"<<endl;
-  cerr<<"MM      MMMMMMMMMMMMMMMMMMM    M"<<endl;
-  cerr<<"MMM   MMMMMMMMMMMMMMMMMMMMMM   M"<<endl;
-  cerr<<"MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM"<<endl;
-  cerr<<message<<endl;
-}
diff --git a/mass_profile/report_error.hpp b/mass_profile/report_error.hpp
deleted file mode 100644
index 30607e5..0000000
--- a/mass_profile/report_error.hpp
+++ /dev/null
@@ -1,6 +0,0 @@
-#ifndef REPORT_ERR
-#define REPORT_ERR
-
-void report_error(const char* message);
-
-#endif
diff --git a/mass_profile/spline.hpp b/mass_profile/spline.hpp
deleted file mode 100644
index 1027b61..0000000
--- a/mass_profile/spline.hpp
+++ /dev/null
@@ -1,110 +0,0 @@
-#ifndef SPLINE_HPP
-#define SPLINE_HPP
-
-#include <vector>
-#include <cstdlib>
-#include <cassert>
-#include <cmath>
-#include <limits>
-
-template <typename T>
-class spline
-{
-public:
-  std::vector<T> x_list;
-  std::vector<T> y_list;
-  std::vector<T> y2_list;
-
-public:
-  void push_point(T x,T y)
-  {
-    if(!x_list.empty())
-      {
-	assert(x>*(x_list.end()-1));
-      }
-    x_list.push_back(x);
-    y_list.push_back(y);
-  }
-
-  T get_value(T x)
-  {
-    if(x<=x_list[0])
-      {
-	return y_list[0];
-      }
-    if(x>=x_list.back())
-      {
-	return y_list.back();
-      }
-    assert(x_list.size()==y2_list.size());
-    assert(x>x_list[0]);
-    assert(x<x_list.back());
-    int n1,n2;
-    n1=0;
-    n2=x_list.size()-1;
-    while((n2-n1)!=1)
-      {
-	//cerr<<n1<<"\t"<<n2<<endl;
-	if(x_list[n1+1]<=x)
-	  {
-	    n1++;
-	  }
-	if(x_list[n2-1]>x)
-	  {
-	    n2--;
-	  }
-      }
-    T h=x_list[n2]-x_list[n1];
-    double a=(x_list[n2]-x)/h;
-    double b=(x-x_list[n1])/h;
-    return a*y_list[n1]+b*y_list[n2]+((a*a*a-a)*y2_list[n1]+
-				      (b*b*b-b)*y2_list[n2])*(h*h)/6.;
-
-  }
-
-  void gen_spline(T y2_0,T y2_N)
-  {
-    int n=x_list.size();
-    y2_list.resize(0);
-    y2_list.resize(x_list.size());
-    std::vector<T> u(x_list.size());
-    if(std::abs(y2_0)<std::numeric_limits<T>::epsilon())
-      {
-	y2_list[0]=0;
-	u[0]=0;
-      }
-    else
-      {
-	y2_list[0]=-.5;
-	u[0]=(3./(x_list[1]-x_list[0]))*((y_list[1]-y_list[0])/(x_list[1]-x_list[0])-y2_0);
-      }
-    for(int i=1;i<n-1;++i)
-      {
-	double sig=(x_list[i]-x_list[i-1])/(x_list[i+1]-x_list[i-1]);
-	double p=sig*y2_list[i-1]+2.;
-	y2_list[i]=(sig-1.)/p;
-	u[i]=(y_list[i+1]-y_list[i])/(x_list[i+1]-x_list[i])
-	  -(y_list[i]-y_list[i-1])/(x_list[i]-x_list[i-1]);
-	u[i]=(6.*u[i]/(x_list[i+1]-x_list[i-1])-sig*u[i-1])/p;
-      }
-    double qn,un;
-    if(std::abs(y2_N)<std::numeric_limits<T>::epsilon())
-      {
-	qn=un=0;
-      }
-    else
-      {
-	qn=.5;
-	un=(3./(x_list[n-1]-x_list[n-2]))*(y2_N-(y_list[n-1]-y_list[n-2])/(x_list[n-1]-x_list[n-2]));
-
-      }
-    y2_list[n-1]=(un-qn*u[n-2])/(qn*y2_list[n-2]+1.);
-    for(int i=n-2;i>=0;--i)
-      {
-	y2_list[i]=y2_list[i]*y2_list[i+1]+u[i];
-      }
-  }
-
-};
-
-#endif  /* SPLINE_HPP */
diff --git a/mass_profile/vchisq.hpp b/mass_profile/vchisq.hpp
deleted file mode 100644
index 8cd8481..0000000
--- a/mass_profile/vchisq.hpp
+++ /dev/null
@@ -1,137 +0,0 @@
-/**

-   \file vchisq.hpp

-   \brief chi-square statistic

-   \author Junhua Gu

- */

-

-#ifndef VCHI_SQ_HPP

-#define VCHI_SQ_HPP

-

-#define OPT_HEADER

-

-#include <core/fitter.hpp>

-#include <iostream>

-#include <vector>

-#include <misc/optvec.hpp>

-#include <cmath>

-

-using std::cerr;

-using std::endl;

-

-namespace opt_utilities

-{

-  template<typename T>

-  class vchisq

-    :public statistic<std::vector<T>,std::vector<T>,std::vector<T>,T,std::string>

-  {

-  private:

-    bool verb;

-    int n;

-    bool limit_bound;

-    typedef std::vector<T> Tp;

-

-    vchisq<T>* do_clone()const

-    {

-      return new vchisq<T>(*this);

-    }

-

-    const char* do_get_type_name()const

-    {

-      return "chi^2 statistic";

-    }

-

-  public:

-    void verbose(bool v)

-    {

-      verb=v;

-    }

-

-    void set_limit()

-    {

-      limit_bound=true;

-    }

-

-    void clear_limit()

-    {

-      limit_bound=false;

-    }

-  public:

-    vchisq()

-      :verb(false),limit_bound(false)

-    {}

-

-    T do_eval(const std::vector<T>& p)

-    {

-      if(limit_bound)

-	{

-	  if(!this->get_fitter().get_model().meets_constraint(p))

-	    {

-	      return 1e99;

-	    }

-	}

-      T result(0);

-      std::vector<float> vx;

-      std::vector<float> vy;

-      std::vector<float> vye;

-      std::vector<float> my;

-      float x1=1e99,x2=-1e99,y1=1e99,y2=-1e99;

-      if(verb)

-	{

-	  n++;

-	  if(n%100==0)

-	    {

-	      vx.resize(this->get_data_set().get_data(0).get_y().size());

-	      vy.resize(this->get_data_set().get_data(0).get_y().size());

-	      vye.resize(this->get_data_set().get_data(0).get_y().size());

-	      my.resize(this->get_data_set().get_data(0).get_y().size());

-	    }

-	}

-      for(int i=(this->get_data_set()).size()-1;i>=0;--i)

-	{

-	  const std::vector<double> y_model(this->eval_model(this->get_data_set().get_data(i).get_x(),p));

-	  const std::vector<double>& y=this->get_data_set().get_data(i).get_y();

-	  const std::vector<double>& ye=this->get_data_set().get_data(i).get_y_lower_err();

-	  for(size_t j=0;j<y.size();++j)

-	    {

-	      double chi=(y_model[j]-y[j])/ye[j];

-	      result+=chi*chi;

-	    }

-

-	  if(verb&&n%100==0)

-	    {

-	      for(size_t j=0;j<y.size();++j)

-		{

-		  vx.at(j)=((this->get_data_set().get_data(i).get_x().at(j)+this->get_data_set().get_data(i).get_x().at(j+1))/2.);

-		  vy.at(j)=(y[j]);

-		  vye.at(j)=ye[j];

-		  my.at(j)=(y_model[j]);

-		  x1=std::min(vx.at(j),x1);

-		  y1=std::min(vy.at(j),y1);

-		  x2=std::max(vx.at(j),x2);

-		  y2=std::max(vy.at(j),y2);

-		  vye[j]=log10(vy[j]+vye[j])-log10(vy[j]);

-		  vx[j]=log10(vx[j]);

-		  vy[j]=log10(vy[j]);

-		  my[j]=log10(my[j]);

-		}

-	    }

-	}

-      if(verb)

-	{

-	  if(n%100==0)

-	    {

-	      cerr<<result<<"\t";

-	      for(size_t i=0;i<get_size(p);++i)

-		{

-		  cerr<<get_element(p,i)<<",";

-		}

-	      cerr<<endl;

-	    }

-	}

-

-      return result;

-    }

-  };

-

-}

-#endif

diff --git a/mass_profile/wang2012_model.hpp b/mass_profile/wang2012_model.hpp
deleted file mode 100644
index c678970..0000000
--- a/mass_profile/wang2012_model.hpp
+++ /dev/null
@@ -1,67 +0,0 @@
-/**
-   \file wang2012_model.hpp
-   \brief Jingying Wang's model
-   \author Jingying Wang
- */
-
-
-#ifndef WANG2012_MODEL
-#define WANG2012_MODEL
-#define OPT_HEADER
-#include <core/fitter.hpp>
-#include <cmath>
-
-namespace opt_utilities
-{
-  template <typename T>
-  class wang2012_model
-    :public model<T,T,std::vector<T>,std::string>
-  {
-  private:
-    model<T,T,std::vector<T> >* do_clone()const
-    {
-      return new wang2012_model<T>(*this);
-    }
-
-    const char* do_get_type_name()const
-    {
-      return "1d power law";
-    }
-  public:
-    wang2012_model()
-    {
-      this->push_param_info(param_info<std::vector<T> >("A",5,0,500));
-      this->push_param_info(param_info<std::vector<T> >("n",1.66,0,10));
-      this->push_param_info(param_info<std::vector<T> >("xi",0.45,0,1));
-      this->push_param_info(param_info<std::vector<T> >("a2",1500,0,1e8));
-      this->push_param_info(param_info<std::vector<T> >("a3",50,0,1e8));
-      this->push_param_info(param_info<std::vector<T> >("beta",0.49,0.1,0.7));
-      this->push_param_info(param_info<std::vector<T> >("T0",0,0,10));
-      
-    }
-
-    T do_eval(const T& x,const std::vector<T>& param)
-    {
-      T A=param[0];
-      T n=param[1];
-      T xi=param[2];
-      T a2=param[3];
-      T a3=param[4];
-      T beta=param[5];
-      T T0=param[6];
-        return A*(pow(x,n)+xi*a2)/(pow(x,n)+a2)/pow(1+x*x/a3/a3,beta)+T0;
-      //return A*(pow(x,n)+a1)/(pow(x,n)+1)/pow(1+x*x/a3/a3,beta)+T0;
-    }
-
-  private:
-    std::string do_get_information()const
-    {
-      return "";
-    }
-  };
-}
-
-
-
-#endif
-//EOF