diff options
Diffstat (limited to 'mass_profile/projector.hpp')
| -rw-r--r-- | mass_profile/projector.hpp | 214 |
1 files changed, 0 insertions, 214 deletions
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 |