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

1 files changed, 0 insertions, 204 deletions

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