Add mass_profile tools

* These tools are mainly use to calculate the total gravitational mass
  profile, as well as the intermediate products (e.g., surface
  brightness profile fitting, gas density profile, NFW fitting, etc.)
* There are additional tools for calculating the luminosity and flux.
* These tools mainly developed by Junhua GU, and contributed by
  Weitian (Aaron) LI, and Zhenghao ZHU.

1 files changed, 336 insertions, 0 deletions

diff --git a/mass_profile/fit_lt_bpl.cpp b/mass_profile/fit_lt_bpl.cpp
new file mode 100644
index 0000000..5637dca
--- /dev/null
+++ b/mass_profile/fit_lt_bpl.cpp
@@ -0,0 +1,336 @@
+/*
+  Perform a double-beta density model fitting to the surface brightness data
+  Author: Junhua Gu
+  Last modified: 2011.01.01
+  This code is distributed with no warrant
+*/
+
+//#define HAVE_X_ERROR
+#include <iomanip>
+#include <iostream>
+#include <sstream>
+#include <fstream>
+#include <models/bpl1d.hpp>
+#include <models/lin1d.hpp>
+#include "statistics/chisq.hpp"
+#include "statistics/logchisq.hpp"
+#include "statistics/leastsq.hpp"
+#include <data_sets/default_data_set.hpp>
+#include <methods/powell/powell_method.hpp>
+#include <methods/gsl_simplex/gsl_simplex.hpp>
+#include <core/freeze_param.hpp>
+
+using namespace std;
+using namespace opt_utilities;
+//double s=5.63136645E20;
+const double kpc=3.086E21;//kpc in cm
+const double Mpc=kpc*1000;
+const double pi=4*atan(1);
+double std_norm_rand()
+{
+  double u=0;
+  double v=0;
+  while(u<=0||v<=0)
+    {
+      u=rand()/(double)RAND_MAX;
+      rand();
+      v=rand()/(double)RAND_MAX;
+    }
+  double x=std::sqrt(-log(u))*cos(2*pi*v);
+  return x;
+}
+
+double shuffle_data(double xc,double xl,double xu)
+{
+  double lxc=log(xc);
+  double lxl=log(xc-xl)-log(xc);
+  double lxu=log(xc+xu)-log(xc);
+    
+  if(std_norm_rand()>0)
+    {
+      double result=std::exp(lxc-std::abs(std_norm_rand()*lxl));
+      return result;
+    }
+  else
+    {
+      double result=std::exp(lxc+std::abs(std_norm_rand()*lxu));
+      return result;
+    }
+}
+
+int main(int argc,char* argv[])
+{
+  srand(time(0));
+  if(argc!=4)
+    {
+      cerr<<"Usage:"<<argv[0]<<" <a 5 column file with T -Terr +Terr L Lerr> <initial broken temperature> <T lower limit>"<<endl;
+      return -1;
+    }
+  double T_lower_limit(atof(argv[3]));
+  double Tb=atof(argv[2]);
+  assert(Tb>0);
+  ifstream ifs_data(argv[1]);
+  default_data_set<double,double> ds;
+  ofstream ofs_result("l-t_bpl_result.qdp");
+  ofs_result<<"read terr 1 2"<<endl;
+  ofs_result<<"skip single"<<endl;
+  ofs_result<<"log"<<endl;
+  //ofs_result<<"li on 2"<<endl;
+  ofs_result<<"time off"<<endl;
+  ofs_result<<"la f"<<endl;
+  ofs_result<<"la x temperature (keV)"<<endl;
+  ofs_result<<"la y luminosity (10\\u43\\d erg s\\u-1\\d)"<<endl;
+  double yunit=1E43;
+  double sxxl=0;
+  double s1l=0;
+  double sxl=0;
+  double syl=0;
+  double sxyl=0;
+
+  double sxxu=0;
+  double s1u=0;
+  double sxu=0;
+  double syu=0;
+  double sxyu=0;
+
+  bool is_first_nonono=true;
+
+  for(;;)
+    {
+      double T,Tl,Tu;
+      double L,Ll,Lu;
+      std::string line;
+      getline(ifs_data,line);
+      //ifs_data>>T>>Tl>>Tu>>M>>Ml>>Mu;
+      if(!ifs_data.good())
+	{
+	  break;
+	}
+      line+=" ";
+      istringstream iss(line);
+      
+      if(line[0]=='#')
+	{
+	  if(!is_first_nonono)
+	    {
+	      ofs_result<<"no no no"<<endl;
+	    }
+	  else
+	    {
+	      is_first_nonono=false;
+	    }
+	  continue;
+	}
+      iss>>T>>Tl>>Tu>>L>>Ll;
+      Lu=Ll;
+      //std::cerr<<L<<"\t"<<Lerr<<endl;
+      if(!iss.good())
+	{
+	  continue;
+	}
+
+      if(T<T_lower_limit||L<0)
+	{
+	  continue;
+	}
+      if(std::abs(Lu)+std::abs(Ll)<L*.1)
+	{
+	  double k=L*.1/(std::abs(Lu)+std::abs(Ll));
+	  Lu*=k;
+	  Ll*=k;
+	}
+      if(std::abs(Ll)>std::abs(L))
+	{
+	  continue;
+	}
+      Tl=std::abs(Tl);
+      Tu=std::abs(Tu);
+      Ll=std::abs(Ll);
+      Lu=std::abs(Lu);
+      ofs_result<<T<<"\t"<<-std::abs(Tl)<<"\t"<<+std::abs(Tu)<<"\t"<<L/yunit<<"\t"<<-std::abs(Ll/yunit)<<"\t"<<+std::abs(Lu/yunit)<<endl;
+      double x=(T);
+      double y=(L);
+      double xu=Tu;
+      double xl=Tl;
+
+      double yu=Lu;
+      double yl=Ll;
+      if(T>Tb)
+	{
+	  sxxl+=log(x)*log(x);
+	  sxl+=log(x);
+	  syl+=log(y);
+	  sxyl+=log(y)*log(x);
+	  s1l+=1;
+	}
+      else
+	{
+	  sxxu+=log(x)*log(x);
+	  sxu+=log(x);
+	  syu+=log(y);
+	  sxyu+=log(y)*log(x);
+	  s1u+=1;
+	}
+      data<double,double> d(x,y,std::abs(yl),std::abs(yu),
+			    std::abs(xl),std::abs(xu));
+      ds.add_data(d);
+    }
+
+  double Ml=sxxl*s1l-sxl*sxl;
+  double Mal=sxyl*s1l-syl*sxl;
+  double Mbl=sxxl*syl-sxl*sxyl;
+  double k0l=Mal/Ml;
+  double b0l=Mbl/Ml;
+
+  double Mu=sxxu*s1u-sxu*sxu;
+  double Mau=sxyu*s1u-syu*sxu;
+  double Mbu=sxxu*syu-sxu*sxyu;
+  double k0u=Mau/Mu;
+  double b0u=Mbu/Mu;
+
+  double gamma0l=k0l;
+  double gamma0u=k0u;
+  
+  double ampl0l=exp(b0l)*pow(Tb,gamma0l);
+  double ampl0u=exp(b0u)*pow(Tb,gamma0u);;
+
+  
+
+  ofs_result<<"no no no"<<endl;
+  fitter<double,double,vector<double>,double,std::string> fit;
+  fit.set_opt_method(powell_method<double,vector<double> >());
+  
+  fit.set_statistic(logchisq<double,double,vector<double>,double,std::string>());
+  //fit.set_statistic(leastsq<double,double,vector<double>,double,std::string>());
+  fit.set_model(bpl1d<double>());
+  fit.load_data(ds);
+
+  cerr<<"k0l="<<k0l<<endl;
+  cerr<<"k0u="<<k0u<<endl;
+  cerr<<"Ampl0="<<(ampl0l+ampl0u)/2<<endl;
+  
+  fit.set_param_value("bpx",Tb);
+  fit.set_param_value("bpy",(ampl0l+ampl0u)/2);
+  fit.set_param_value("gamma1",gamma0l);
+  fit.set_param_value("gamma2",gamma0u);
+
+  fit.fit();
+  //fit.set_opt_method(gsl_simplex<double,vector<double> >());
+  fit.fit();
+  std::vector<double> p=fit.fit();
+  Tb=fit.get_param_value("bpx");
+  ///std::cout<<"chi="<<fit.get_statistic().eval(p)<<std::endl;
+  for(double i=.5;i<12;i*=1.01)
+    {
+      ofs_result<<i<<"\t0\t0\t"<<fit.eval_model_raw(i,p)/yunit<<"\t0\t0\n";
+    }
+
+  
+  std::vector<double> mean_p(p.size());
+  std::vector<double> mean_p2(p.size());
+  int cnt=0;
+  for(int n=0;n<100;++n)
+    {
+      ++cnt;
+      cerr<<".";
+      double sxxl=0;
+      double s1l=0;
+      double sxl=0;
+      double syl=0;
+      double sxyl=0;
+      
+      double sxxu=0;
+      double s1u=0;
+      double sxu=0;
+      double syu=0;
+      double sxyu=0;
+
+      opt_utilities::default_data_set<double,double> ds1;
+      for(int i=0;i<ds.size();++i)
+	{
+	  double x=ds.get_data(i).get_x();
+	  double y=ds.get_data(i).get_y();
+	  double xl=ds.get_data(i).get_x_lower_err();
+	  double xu=ds.get_data(i).get_x_upper_err();
+	  double yl=ds.get_data(i).get_y_lower_err();
+	  double yu=ds.get_data(i).get_y_upper_err();
+	  
+	  double new_x=shuffle_data(x,
+				    xl,
+				    xu);
+	  double new_y=shuffle_data(y,
+				    yl,
+				    yu);
+	  
+	  ds1.add_data(data<double,double>(new_x,new_y,
+					   yl/y*new_y,
+					   yu/y*new_y,
+					   xl/x*new_x,
+					   xu/x*new_x));
+
+	  x=new_x;
+	  y=new_y;
+	  if(x>Tb)
+	    {
+	      sxxl+=log(x)*log(x);
+	      sxl+=log(x);
+	      syl+=log(y);
+	      sxyl+=log(y)*log(x);
+	      s1l+=1;
+	    }
+	  else
+	    {
+	      sxxu+=log(x)*log(x);
+	      sxu+=log(x);
+	      syu+=log(y);
+	      sxyu+=log(y)*log(x);
+	      s1u+=1;
+	    }
+	}
+      double Ml=sxxl*s1l-sxl*sxl;
+      double Mal=sxyl*s1l-syl*sxl;
+      double Mbl=sxxl*syl-sxl*sxyl;
+      double k0l=Mal/Ml;
+      double b0l=Mbl/Ml;
+      
+      double Mu=sxxu*s1u-sxu*sxu;
+      double Mau=sxyu*s1u-syu*sxu;
+      double Mbu=sxxu*syu-sxu*sxyu;
+      double k0u=Mau/Mu;
+      double b0u=Mbu/Mu;
+      
+      double gamma0l=k0l;
+      double gamma0u=k0u;
+      
+      double ampl0l=exp(b0l)*pow(Tb,gamma0l);
+      double ampl0u=exp(b0u)*pow(Tb,gamma0u);;
+      
+      fit.set_param_value("bpx",Tb);
+      fit.set_param_value("bpy",(ampl0l+ampl0u)/2);
+      fit.set_param_value("gamma1",gamma0l);
+      fit.set_param_value("gamma2",gamma0u);
+      
+      
+      fit.load_data(ds1);
+      
+      fit.fit();
+      vector<double> p=fit.fit();
+      for(int i=0;i<p.size();++i)
+	{
+	  mean_p[i]+=p[i];
+	  mean_p2[i]+=p[i]*p[i];
+	}
+      //cerr<<fit.get_param_value("gamma1")<<"\t"<<fit.get_param_value("gamma2")<<endl;
+      
+    }
+  vector<double> std_p(p.size());
+  cerr<<endl;
+  for(int i=0;i<mean_p.size();++i)
+    {
+      mean_p[i]/=cnt;
+      mean_p2[i]/=cnt;
+      std_p[i]=std::sqrt(mean_p2[i]-mean_p[i]*mean_p[i]);
+      cout<<fit.get_param_info(i).get_name()<<"= "<<p[i]<<" +/- "<<std_p[i]<<endl;
+    }
+  std::cout<<"Num of sources:"<<ds.size()<<endl;
+}