Remove unused 'fit_{lt,mt}_{bpl,pl}.cpp' tools

1 files changed, 0 insertions, 350 deletions

diff --git a/mass_profile/fit_mt_bpl.cpp b/mass_profile/fit_mt_bpl.cpp
deleted file mode 100644
index e956f4f..0000000
--- a/mass_profile/fit_mt_bpl.cpp
+++ /dev/null
@@ -1,350 +0,0 @@
-/*
-  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 6 column file with T -Terr +Terr M -Merr +Merr> <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("m-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 mass (M\\d\\(2281)\\u)"<<endl;
-  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 M,Ml,Mu;
-      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>>M>>Ml>>Mu;
-      //std::cerr<<L<<"\t"<<Lerr<<endl;
-      if(!iss.good())
-	{
-	  continue;
-	}
-
-      if(T<T_lower_limit||M<0)
-	{
-	  continue;
-	}
-
-      if(std::abs(Mu)<M*.1||std::abs(Ml)<M*.1)
-	{
-	  cerr<<"mass error less than 10%, skipped"<<endl;
-	  cerr<<line<<endl;
-	  continue;
-	}
-#if 1
-      if(std::abs(Tu)<.1||std::abs(Tl)<.1)
-	{
-	  cerr<<"T error less than 10%, skipped"<<endl;
-	  cerr<<line<<endl;
-	  continue;
-	}
-#endif
-
-      if(std::abs(Mu)+std::abs(Ml)<M*.1)
-	{
-	  double k=M*.1/(std::abs(Mu)+std::abs(Ml));
-	  Mu*=k;
-	  Ml*=k;
-	}
-      if(std::abs(Ml)>std::abs(M))
-	{
-	  continue;
-	}
-      Tl=std::abs(Tl);
-      Tu=std::abs(Tu);
-      Ml=std::abs(Ml);
-      Mu=std::abs(Mu);
-      ofs_result<<T<<"\t"<<-std::abs(Tl)<<"\t"<<+std::abs(Tu)<<"\t"<<M<<"\t"<<-std::abs(Ml)<<"\t"<<+std::abs(Mu)<<endl;
-      double x=(T);
-      double y=(M);
-      double xu=Tu;
-      double xl=Tl;
-
-      double yu=Mu;
-      double yl=Ml;
-      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)<<"\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(size_t 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(size_t 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(size_t 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;
-}