Remove the unused 'fit_nfw_sbp.cpp' with 'nfw_ne.hpp'

3 files changed, 1 insertions, 626 deletions

diff --git a/mass_profile/Makefile b/mass_profile/Makefile
index c13b312..67e9aae 100644
--- a/mass_profile/Makefile
+++ b/mass_profile/Makefile
@@ -23,7 +23,7 @@ endif
 
 OPT_UTIL_INC ?= -I../opt_utilities
 
-TARGETS= fit_nfw_sbp fit_dbeta_sbp fit_beta_sbp \
+TARGETS= fit_dbeta_sbp fit_beta_sbp \
 		fit_wang2012_model \
 		fit_nfw_mass fit_mt_pl fit_lt_pl fit_mt_bpl \
 		fit_lt_bpl calc_lx_dbeta calc_lx_beta
@@ -33,8 +33,6 @@ all: $(TARGETS)
 
 # NOTE:
 # Object/source files should placed *before* libraries (order matters)
-fit_nfw_sbp: fit_nfw_sbp.o
-	$(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC)
 
 fit_dbeta_sbp: fit_dbeta_sbp.o beta_cfg.o report_error.o
 	$(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC)
@@ -67,9 +65,6 @@ fit_lt_bpl: fit_lt_bpl.cpp
 	$(CXX) $(CXXFLAGS) $< -o $@ $(OPT_UTIL_INC)
 
 
-fit_nfw_sbp.o: fit_nfw_sbp.cpp nfw_ne.hpp $(HEADERS)
-	$(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC)
-
 fit_dbeta_sbp.o: fit_dbeta_sbp.cpp constrained_dbeta.hpp $(HEADERS)
 	$(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC)
 
diff --git a/mass_profile/fit_nfw_sbp.cpp b/mass_profile/fit_nfw_sbp.cpp
deleted file mode 100644
index 9ac8401..0000000
--- a/mass_profile/fit_nfw_sbp.cpp
+++ /dev/null
@@ -1,416 +0,0 @@
-/*
-  Fitting the surface brightness profile with an NFW-based surface brightness model
-  Author: Junhua Gu
-  Last modification 20120721
-  The temperature is assumed to be an allen model with a minimum temperature assumed
-*/
-
-
-#include <iostream>
-#include <fstream>
-#include "vchisq.hpp"
-#include "nfw_ne.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 M_sun=1.988E33;//solar mass in g
-const double kpc=3.086E21;//kpc in cm
-
-//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)
-  {
-    /*
-    if(x<=spl.x_list[0])
-      {
-	return spl.y_list[0];
-      }
-    if(x>=spl.x_list.back())
-      {
-	return spl.y_list.back();
-      }
-    */
-    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);
-  }
-};
-
-//Allen temperature model
-
-int main(int argc,char* argv[])
-{
-  if(argc!=2)
-    {
-      cerr<<argv[0]<<" <configure file>"<<endl;
-      cerr<<"Here is a sample of the configure file"<<endl;
-      cerr<<"radius_file	radius1.dat\n"
-	"sbp_file\tsbp1.dat\n"
-	"cfunc_file\tcfunc.dat\n"
-	"n0\t\t.04\n"
-	"rs\t\t816\n"
-	"rho0\t\t.01\n"
-	"bkg\t\t0\n"
-	"cm_per_pixel\t1.804E21\n"
-	"z\t\t0.062476\n"
-	"T_file\t\tT.dat"
-	<<endl;
-      cerr<<"Notes:"<<endl;
-      cerr<<"n0 is in the unit of cm^-3"<<endl;
-      cerr<<"rs is in the unit of pixel"<<endl;
-      cerr<<"rho0 is in the unit of mass of proton per cm^3"<<endl;
-
-      return -1;
-    }
-  //define a map to store the parameters
-  std::map<std::string,std::string> arg_map;
-  //open the configuration file
-  ifstream cfg_file(argv[1]);
-  assert(cfg_file.is_open());
-  for(;;)
-    {
-      std::string key;
-      std::string value;
-      cfg_file>>key>>value;
-      if(!cfg_file.good())
-	{
-	  cfg_file.close();
-	  break;
-	}
-      arg_map[key]=value;
-    }
-  //check whether following parameters are defined in the configuration file
-  assert(arg_map.find("radius_file")!=arg_map.end());
-  assert(arg_map.find("sbp_file")!=arg_map.end());
-  assert(arg_map.find("cfunc_file")!=arg_map.end());
-  assert(arg_map.find("T_file")!=arg_map.end());
-  assert(arg_map.find("z")!=arg_map.end());
-  const double z=atof(arg_map["z"].c_str());
-  double r_min=0;
-  if(arg_map.find("r_min")!=arg_map.end())
-    {
-      r_min=atof(arg_map["r_min"].c_str());
-      cerr<<"r_min presents and its value is "<<r_min<<endl;
-    }
-  //note that in this program, the radius are not the central value of each annuli or pie region, but the boundaries.
-  //for example, if we have a set of radius and surface brightness values as follows:
-  /*
-    radius      width       surface brightness
-    1           1           x
-    2           1           x
-    3           1           x
-
-    then the radius is stored as
-    0 1.5 2.5 3.5
-    note that there should be 4 radius values, although only to represent 3 annuli.
-    this will be convenient to calculate the volume of each spherical shell,
-    and can naturally ensure the annuli are adjacent with each other, with out any gaps.
-
-
-   */
-
-  std::vector<double> radii;//to store radius
-  std::vector<double> sbps;//to store the surface brightness value
-  std::vector<double> sbpe;//to store the sbp error
-  //read in radius file
-  /*
-    About the format of the radius file:
-    the radius file contains only radius, separated by space or line feed (i.e., the <ENTER> key).
-    the unit should be pixel
-
-    The number of radius can be larger than the number of annuli+1, the exceeded radius can be used
-    to calculate the influence of outer shells.
-   */
-  int ncut=0;
-  for(ifstream ifs(arg_map["radius_file"].c_str());;)
-    {
-      assert(ifs.is_open());
-      double x;
-      ifs>>x;
-      if(!ifs.good())
-	{
-	  break;
-	}
-      if(x<r_min)
-	{
-	  ++ncut;
-	  continue;
-	}
-      cerr<<x<<endl;
-      radii.push_back(x);
-    }
-  //read in surface brightness file
-  /*
-    the surface brightness file contains two columns, that are surface brightness and the error, respectively.
-   */
-  for(ifstream ifs(arg_map["sbp_file"].c_str());;)
-    {
-      assert(ifs.is_open());
-      double x,xe;
-      ifs>>x>>xe;
-      if(!ifs.good())
-	{
-	  break;
-	}
-      if(ncut)
-	{
-	  --ncut;
-	  continue;
-	}
-      cerr<<x<<"\t"<<xe<<endl;
-      sbps.push_back(x);
-      sbpe.push_back(xe);
-    }
-  //cerr<<radii.size()<<"\t"<<sbps.size()<<endl;
-  //return 0;
-  spline_func_obj cf;
-
-  //read in the cooling function file
-  /*
-    the cooling function file contains two columns, that are radius (central value, in pixel) and the ``reduced cooling function''
-    the reduced cooling function is calculated as follows
-    by using wabs*apec model, fill the kT, z, nH to the actual value (derived from deproject spectral fitting),
-    and fill norm as 1E-14/(4*pi*(Da*(1+z))^2).
-    then use flux e1 e2 (e1 and e2 are the energy limits of the surface brightness profile) to get the cooling function (in photon counts,
-    not in erg/s)
-   */
-  for(ifstream ifs(arg_map["cfunc_file"].c_str());;)
-    {
-      assert(ifs.is_open());
-      double x,y;
-      ifs>>x>>y;
-      if(!ifs.good())
-	{
-	  break;
-	}
-      cerr<<x<<"\t"<<y<<endl;
-      //cf.add_point(x,y*2.1249719395939022e-68);//change with source
-      cf.add_point(x,y);//change with source
-    }
-  cf.gen_spline();
-
-  for(double x=0;x<1000;x++)
-    {
-      //cout<<x<<"\t"<<cf(x)<<endl;
-    }
-  //return 0;
-  //cout<<radii.size()<<endl;
-  //cout<<sbps.size()<<endl;
-
-  //initial a data set object and put the data together
-  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 a fitter object
-  fitter<vector<double>,vector<double>,vector<double>,double> f;
-  //load the data set into the fitter object
-  f.load_data(ds);
-  //define a projector object
-  //see projector for more detailed information
-  projector<double> a;
-  //define the nfw surface brightness profofile model
-  nfw_ne<double> nfw;
-  //attach the cooling function into the projector
-  a.attach_cfunc(cf);
-  assert(arg_map.find("cm_per_pixel")!=arg_map.end());
-  //set the cm to pixel ratio
-  double cm_per_pixel=atof(arg_map["cm_per_pixel"].c_str());
-  a.set_cm_per_pixel(cm_per_pixel);
-  nfw.set_cm_per_pixel(cm_per_pixel);
-  //define the temperature profile model
-  spline_func_obj tf;
-
-  for(ifstream ifs_tfunc(arg_map["T_file"].c_str());;)
-    {
-      assert(ifs_tfunc.is_open());
-      double x,y;
-      ifs_tfunc>>x>>y;
-      if(!ifs_tfunc.good())
-	{
-	  break;
-	}
-      if(x<r_min)
-	{
-	  continue;
-	}
-      tf.add_point(x,y);
-    }
-  tf.gen_spline();
-
-  //attach the temperature, surface brightness model and projector together
-  nfw.attach_Tfunc(tf);
-  a.attach_model(nfw);
-  f.set_model(a);
-  //define the chi-square statistic
-  vchisq<double> c;
-  c.verbose(true);
-  f.set_statistic(c);
-  //set the optimization method, here we use powell method
-  f.set_opt_method(powell_method<double,std::vector<double> >());
-  //set the initial values
-  double n0=atof(arg_map["n0"].c_str());
-  double rho0=atof(arg_map["rho0"].c_str());
-  double rs=atof(arg_map["rs"].c_str());
-  double bkg=atof(arg_map["bkg"].c_str());
-  f.set_param_value("n0",n0);
-  f.set_param_value("rho0",rho0);
-  f.set_param_value("rs",rs);
-
-
-  f.set_param_value("bkg",bkg);
-
-  cout<<f.get_data_set().size()<<endl;
-  cout<<f.get_num_params()<<endl;
-#if 1
-  //perform the fitting
-  f.fit();
-  f.set_precision(1e-10);
-  f.fit();
-  f.clear_param_modifier();
-  f.fit();
-#endif
-  //output the parameters
-  ofstream param_output("nfw_param.txt");
-
-  for(size_t i=0;i<f.get_num_params();++i)
-    {
-      cout<<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;
-    }
-  c.verbose(false);
-  f.set_statistic(c);
-#if 1
-  f.fit();
-  f.fit();
-#endif
-  //fetch the fitting result
-  std::vector<double> p=f.get_all_params();
-  f.clear_param_modifier();
-  std::vector<double> mv=f.eval_model(radii,p);
-  cerr<<mv.size()<<endl;
-  //output the results
-  ofstream ofs_sbp("sbp_fit.qdp");
-  ofstream ofs_resid("resid.qdp");
-  ofs_resid<<"read serr 2"<<endl;
-  //output the surface brightness profile
-  ofs_sbp<<"read serr 2"<<endl;
-  ofs_sbp<<"skip single"<<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<<"\t"<<ym<<endl;
-      ofs_resid<<x*cm_per_pixel/kpc<<"\t"<<(y-ym)/ye<<"\t1\n";
-    }
-  //output the electron density
-  mv=nfw.eval(radii,p);
-  ofstream ofs_rho("rho_fit.qdp");
-  for(size_t 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;
-    }
-  //output integral mass profile
-  rho0=f.get_param_value("rho0")*1.67E-24;
-  rs=f.get_param_value("rs")*cm_per_pixel;
-  ofstream ofs_int_mass("mass_int.qdp");
-  for(double r=0;r<2000*kpc;r+=kpc)
-    {
-      ofs_int_mass<<r/kpc<<"\t"<<nfw_mass_enclosed(r,abs(rho0),abs(rs))/M_sun<<endl;
-    }
-  //calculate the overdensity profile
-  ofstream ofs_overdensity("overdensity.qdp");
-
-  std::vector<double> radius_list;
-  std::vector<double> delta_list;
-
-  cerr<<"delta\tr_delta (kpc)\tr_delta (pixel)\tmass_delta (solar mass)\n";
-  for(double r=kpc;r<6000*kpc;r+=kpc)
-    {
-      double delta=nfw_average_density(r,abs(rho0),abs(rs))/calc_critical_density(z);
-      radius_list.push_back(r);
-      delta_list.push_back(delta);
-
-      /*
-      if(delta<=200&&!hit_200)
-	{
-	  hit_200=true;
-	  cerr<<200<<"\t"<<r/kpc<<"\t\t"<<r/cm_per_pixel<<"\t\t"<<nfw_mass_enclosed(r,abs(rho0),abs(rs))/M_sun<<endl;
-	  break;
-	}
-      if(delta<=500&&!hit_500)
-	{
-	  hit_500=true;
-	  cerr<<500<<"\t"<<r/kpc<<"\t\t"<<r/cm_per_pixel<<"\t\t"<<nfw_mass_enclosed(r,abs(rho0),abs(rs))/M_sun<<endl;
-	}
-      */
-      ofs_overdensity<<r/kpc<<"\t"<<delta<<endl;
-    }
-
-  for(size_t i=0;i<radius_list.size()-1;++i)
-    {
-      double r=radius_list[i];
-      if(delta_list[i]>=200&&delta_list[i+1]<200)
-	{
-	  cerr<<200<<"\t"<<r/kpc<<"\t\t"<<r/cm_per_pixel<<"\t\t"<<nfw_mass_enclosed(r,abs(rho0),abs(rs))/M_sun<<endl;
-	}
-
-      if(delta_list[i]>=500&&delta_list[i+1]<500)
-	{
-	  cerr<<500<<"\t"<<r/kpc<<"\t\t"<<r/cm_per_pixel<<"\t\t"<<nfw_mass_enclosed(r,abs(rho0),abs(rs))/M_sun<<endl;
-	}
-
-      if(delta_list[i]>=1500&&delta_list[i+1]<1500)
-	{
-	  cerr<<1500<<"\t"<<r/kpc<<"\t\t"<<r/cm_per_pixel<<"\t\t"<<nfw_mass_enclosed(r,abs(rho0),abs(rs))/M_sun<<endl;
-	}
-
-      if(delta_list[i]>=2500&&delta_list[i+1]<2500)
-	{
-	  cerr<<2500<<"\t"<<r/kpc<<"\t\t"<<r/cm_per_pixel<<"\t\t"<<nfw_mass_enclosed(r,abs(rho0),abs(rs))/M_sun<<endl;
-	}
-
-    }
-
-  //output the M200 and R200
-  //cerr<<"M200="<<nfw_mass_enclosed(r200,abs(rho0),abs(rs))/M_sun<<" solar mass"<<endl;
-  //cerr<<"R200="<<r200/kpc<<" kpc"<<endl;
-  //for(int i=0;i<p.size();++i)
-  //{
-  //cerr<<p[i]<<endl;
-  //}
-  return 0;
-}
diff --git a/mass_profile/nfw_ne.hpp b/mass_profile/nfw_ne.hpp
deleted file mode 100644
index 1eeb17c..0000000
--- a/mass_profile/nfw_ne.hpp
+++ /dev/null
@@ -1,204 +0,0 @@
-/*
-  Gas density profile derived from nfw mass profile and temperature profile
-  Author: Junhua Gu
-  Last modification: 20120721
-*/
-
-#ifndef NFW_NE
-#define NFW_NE
-#include "projector.hpp"
-#include <algorithm>
-#include <functional>
-#include <numeric>
-
-//a series of physical constants
-static 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 k=1.60217646E-9;//erg/keV
-static const double c=2.99792458E10;//cm/s
-
-
-namespace opt_utilities
-{
-  //the nfw mass enclosed within a radius r, with parameter rho0 and rs
-  template <typename T>
-  T nfw_mass_enclosed(T r,T rho0,T rs)
-  {
-    return 4*pi*rho0*rs*rs*rs*(std::log((r+rs)/rs)-r/(r+rs));
-  }
-
-  //average mass density
-  template <typename T>
-  T nfw_average_density(T r,T rho0,T rs)
-  {
-    if(r==0)
-      {
-	return rho0;
-      }
-
-    return nfw_mass_enclosed(r,rho0,rs)/(4.*pi/3*r*r*r);
-  }
-
-  //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 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 wraps method of calculating gas density from mass profile and temperature profile
-  template <typename T>
-  class nfw_ne
-    :public model<std::vector<T>,std::vector<T>,std::vector<T> >
-  {
-  private:
-    //pointer to temperature profile function
-    func_obj<T,T>* pTfunc;
-    //cm per pixel
-    T cm_per_pixel;
-  public:
-    //default constructor
-    nfw_ne()
-      :pTfunc(0),cm_per_pixel(1)
-    {
-
-      this->push_param_info(param_info<std::vector<T>,std::string>("rho0",1));//in mp
-      this->push_param_info(param_info<std::vector<T>,std::string>("rs",100));
-      this->push_param_info(param_info<std::vector<T>,std::string>("n0",.01));
-    }
-
-    //copy constructor
-    nfw_ne(const nfw_ne& rhs)
-      :cm_per_pixel(rhs.cm_per_pixel)
-    {
-      if(rhs.pTfunc)
-	{
-	  pTfunc=rhs.pTfunc->clone();
-	}
-      else
-	{
-	  pTfunc=0;
-	}
-      //initial parameter list
-      this->push_param_info(param_info<std::vector<T>,std::string>("rho0",rhs.get_param_info("rho0").get_value()));
-      this->push_param_info(param_info<std::vector<T>,std::string>("rs",rhs.get_param_info("rs").get_value()));
-      this->push_param_info(param_info<std::vector<T>,std::string>("n0",rhs.get_param_info("n0").get_value()));
-    }
-
-    //assignment operator
-    nfw_ne& operator=(const nfw_ne& rhs)
-    {
-      cm_per_pixel=rhs.cm_per_pixel;
-      if(pTfunc)
-	{
-	  pTfunc->destroy();
-	}
-      if(rhs.pTfunc)
-	{
-	  pTfunc=rhs.pTfunc->clone();
-	}
-    }
-
-    //destructor
-    ~nfw_ne()
-    {
-      if(pTfunc)
-	{
-	  pTfunc->destroy();
-	}
-    }
-
-  public:
-    //attach the temperature profile function
-    void attach_Tfunc(const func_obj<T,T>& Tf)
-    {
-      if(pTfunc)
-	{
-	  pTfunc->destroy();
-	}
-      pTfunc=Tf.clone();
-    }
-
-    //set the cm per pixel value
-    void set_cm_per_pixel(const T& x)
-    {
-      cm_per_pixel=x;
-    }
-
-    //clone self
-    nfw_ne<T>* do_clone()const
-    {
-      return new nfw_ne<T>(*this);
-    }
-
-
-    //calculate density under parameters p, at radius r
-    /*
-      r is a vector, which stores a series of radius values
-      the annuli or pie regions are enclosed between any two
-      adjacent radii.
-      so the returned value has length smaller than r by 1.
-     */
-    std::vector<T> do_eval(const std::vector<T> & r,
-			   const std::vector<T>& p)
-    {
-      assert(pTfunc);
-      //const T kT_erg=k*5;
-      T rho0=std::abs(p[0])*mp;
-      T rs=std::abs(p[1]);
-      T n0=std::abs(p[2]);
-      T rs_cm=rs*cm_per_pixel;
-
-      std::vector<T> yvec(r.size());
-      const T kT_erg0=pTfunc->eval((r.at(0)+r.at(1))/2)*k;
-      //calculate the integration
-#ifdef _OPENMP
-#pragma omp parallel for schedule(dynamic)
-#endif
-      for(size_t i=0;i<r.size();++i)
-	{
-	  T r_cm=r[i]*cm_per_pixel;
-	  T kT_erg=pTfunc->eval(r[i])*k;
-	  if(abs(r_cm)==0)
-	    {
-	      continue;
-	    }
-	  yvec.at(i)=G*nfw_mass_enclosed(r_cm,rho0,rs_cm)*mu*mp/kT_erg/r_cm/r_cm;
-	  //std::cout<<r_cm/1e20<<"\t"<<nfw_mass_enclosed(r_cm,rho0,rs_cm)/1e45<<std::endl;
-	  //std::cout<<r_cm/1e20<<"\t"<<G*nfw_mass_enclosed(r_cm,rho0,rs_cm)*mu*mp/kT_erg/r_cm/r_cm<<std::endl;
-	}
-
-      std::vector<T> ydxvec(r.size()-1);
-#ifdef _OPENMP
-#pragma omp parallel for schedule(dynamic)
-#endif
-      for(size_t i=1;i<r.size();++i)
-	{
-	  T dr=r[i]-r[i-1];
-	  T dr_cm=dr*cm_per_pixel;
-	  ydxvec.at(i-1)=(yvec[i]+yvec[i-1])/2*dr_cm;
-	}
-      std::partial_sum(ydxvec.begin(),ydxvec.end(),ydxvec.begin());
-      //construct the result
-      std::vector<T> result(r.size()-1);
-#ifdef _OPENMP
-#pragma omp parallel for schedule(dynamic)
-#endif
-      for(size_t i=0;i<r.size()-1;++i)
-	{
-	  T y=-ydxvec.at(i);
-	  T kT_erg=pTfunc->eval(r[i])*k;
-	  //std::cout<<y<<std::endl;
-	  result.at(i)=n0*exp(y)*kT_erg0/kT_erg;
-	}
-      return result;
-    }
-  };
-}
-
-#endif