diff options
Diffstat (limited to 'mass_profile/nfw_ne.hpp')
| -rw-r--r-- | mass_profile/nfw_ne.hpp | 198 |
1 files changed, 198 insertions, 0 deletions
diff --git a/mass_profile/nfw_ne.hpp b/mass_profile/nfw_ne.hpp new file mode 100644 index 0000000..5842bb8 --- /dev/null +++ b/mass_profile/nfw_ne.hpp @@ -0,0 +1,198 @@ +/* + 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 +#pragma omp parallel for schedule(dynamic) + for(int 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); +#pragma omp parallel for schedule(dynamic) + for(int 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); +#pragma omp parallel for schedule(dynamic) + for(int 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 |