diff options
Diffstat (limited to 'src/projector.hpp')
| -rw-r--r-- | src/projector.hpp | 214 |
1 files changed, 214 insertions, 0 deletions
diff --git a/src/projector.hpp b/src/projector.hpp new file mode 100644 index 0000000..4ef4b32 --- /dev/null +++ b/src/projector.hpp @@ -0,0 +1,214 @@ +#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 |