Rename mass_profile to src; Add install & uninstall to Makefile

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