aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rw-r--r--mass_profile/Makefile8
-rw-r--r--mass_profile/adapt_trapezoid.h124
-rw-r--r--mass_profile/calc_distance.cc61
-rw-r--r--mass_profile/calc_distance.h8
-rw-r--r--mass_profile/calc_distance_main.cc52
5 files changed, 1 insertions, 252 deletions
diff --git a/mass_profile/Makefile b/mass_profile/Makefile
index 4f037e6..1d7f48b 100644
--- a/mass_profile/Makefile
+++ b/mass_profile/Makefile
@@ -24,7 +24,7 @@ endif
OPT_UTIL_INC ?= -I../opt_utilities
TARGETS= fit_nfw_sbp fit_dbeta_sbp fit_beta_sbp \
- fit_direct_beta calc_distance fit_wang2012_model \
+ fit_direct_beta fit_wang2012_model \
fit_nfw_mass fit_mt_pl fit_lt_pl fit_mt_bpl \
fit_lt_bpl cooling_time calc_lx_dbeta calc_lx_beta
HEADERS= projector.hpp spline.hpp vchisq.hpp
@@ -57,9 +57,6 @@ calc_lx_dbeta: calc_lx_dbeta.o beta_cfg.o report_error.o
calc_lx_beta: calc_lx_beta.o beta_cfg.o dump_fit_qdp.o
$(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC)
-calc_distance: calc_distance_main.o calc_distance.o
- $(CXX) $(CXXFLAGS) -o $@ $^
-
fit_mt_pl: fit_mt_pl.cpp
$(CXX) $(CXXFLAGS) $< -o $@ $(OPT_UTIL_INC)
@@ -106,9 +103,6 @@ report_error.o: report_error.cpp report_error.hpp
dump_fit_qdp.o: dump_fit_qdp.cpp dump_fit_qdp.hpp
$(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC)
-calc_distance.o: calc_distance.cc calc_distance.h adapt_trapezoid.h
- $(CXX) $(CXXFLAGS) -c $<
-
%.o: %.cc
$(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC)
diff --git a/mass_profile/adapt_trapezoid.h b/mass_profile/adapt_trapezoid.h
deleted file mode 100644
index 35d22b4..0000000
--- a/mass_profile/adapt_trapezoid.h
+++ /dev/null
@@ -1,124 +0,0 @@
-#ifndef ADAPT_TRAPEZOID_H
-#define ADAPT_TRAPEZOID_H
-#include <list>
-#include <utility>
-#include <cassert>
-//#include <iostream>
-
-//using namespace std;
-
-template<typename T1,typename T2,typename T3>
-class triple
-{
- public:
- T1 first;
- T2 second;
- T3 third;
- triple(T1 x1,T2 x2,T3 x3)
- :first(x1),second(x2),third(x3)
- {
- }
-};
-
-
-template<typename T1,typename T2,typename T3>
-triple<T1,T2,T3> make_triple(T1 x1,T2 x2,T3 x3)
-{
- return triple<T1,T2,T3>(x1,x2,x3);
-}
-
-template <typename T>
-T trapezoid(T (*fun)(T),T x1,T x2,T err_limit)
-{
- int n=256;
- T result;
- const int max_division=24;
- T old_value=0;
- for(int i=1;i<max_division;i++)
- {
- result=0.;
- n*=2;
- T step=(x2-x1)/n;
- for(int j=0;j<n;j++)
- {
- result+=(fun(x1+(j+1)*step)+fun(x1+j*step))*step/T(2.);
- }
- old_value-=result;
- old_value=old_value<0?-old_value:old_value;
- if(old_value<err_limit)
- {
- return result;
- }
- old_value=result;
- }
-}
-
-
-template <typename T>
-T adapt_trapezoid(T (*fun)(T),T x1,T x2,T err_limit)
-{
- // const err_limit=.001;
- typedef triple<T,T,bool> interval;
- /*T for interval type,
- bool for state trur for still to be updated,
- false for do not need to be updated
- */
- std::list<interval> interval_list;
- T current_sum=((fun(x1)+fun(x2))/2.*(x2-x1));
- interval_list.push_back(make_triple(x1,current_sum,true));
- interval_list.push_back(make_triple(x2,(T)0.,true));
- bool int_state=1;
- int n_intervals=1;
- while(int_state)
- {
- //std::cout<<n_intervals<<std::endl;
- int_state=0;
- typename std::list<interval>::iterator i1=interval_list.begin();
- typename std::list<interval>::iterator i2=interval_list.begin();
- i2++;
- for(;i2!=interval_list.end();i1++,i2=i1,i2++)
- {
- //cout<<i1->first<<"\t"<<i2->first<<endl;
- //assert(i2->first>i1->first);
- if(i1->third)
- {
- interval new_interval((i1->first+i2->first)/2,0,true);
-
- T sum1,sum2;
- sum1=(fun(new_interval.first)+fun(i1->first))/2*(new_interval.first-i1->first);
- sum2=(fun(new_interval.first)+fun(i2->first))/2*(i2->first-new_interval.first);
- new_interval.second=sum2;
- T err;
- err=i1->second-sum1-sum2;
- err=err<0?-err:err;
-
- if(err>err_limit/n_intervals)
- {
- i1->second=sum1;
- interval_list.insert(i2,new_interval);
- n_intervals++;
- if(n_intervals>10e6)
- {
-
- break;
- }
- }
- else
- {
- i1->third=false;
- }
- int_state=1;
- }
- }
-
- }
- T result=0;
- for(typename std::list<interval>::iterator i=interval_list.begin();i!=interval_list.end();i++)
- {
- result+=i->second;
- }
- return result;
-}
-
-#endif
-//end of the file
diff --git a/mass_profile/calc_distance.cc b/mass_profile/calc_distance.cc
deleted file mode 100644
index 91f2ea3..0000000
--- a/mass_profile/calc_distance.cc
+++ /dev/null
@@ -1,61 +0,0 @@
-#include <iostream>
-#include <cmath>
-#include <cstdlib>
-#include <cstddef>
-#include <cassert>
-#include "adapt_trapezoid.h"
-
-//calc_distance
-//usage:
-//calc_distance z
-
-using namespace std;
-
-static double cm=1;
-static double s=1;
-static double km=1000*100;
-static double Mpc=3.08568e+24*cm;
-//static double kpc=3.08568e+21*cm;
-//static double yr=365.*24.*3600.;
-//static double Gyr=1e9*yr;
-static double H=71.*km/s/Mpc;
-static const double c=299792458.*100.*cm;
-//const double c=3e8*100*cm;
-static const double omega_m=0.27;
-static const double omega_l=0.73;
-static const double arcsec2arc_ratio=1./60/60/180*3.1415926;
-
-
-double E(double z)
-{
- double omega_k=1-omega_m-omega_l;
- return sqrt(omega_m*(1+z)*(1+z)*(1+z)+omega_k*(1+z)*(1+z)+omega_l);
-}
-
-double f_dist(double z)
-{
- return 1/E(z);
-}
-
-double f_age(double z)
-{
- return f_dist(1/z)/(z*z);
-}
-
-
-
-double calc_angular_distance(double z)
-{
- //return c/H*integer(f_dist,0,z)/(1+z);
- //return c/H*adapt_trapezoid(f_dist,0.,z,1e-4)/(1+z);
- return adapt_trapezoid(f_dist,0.,z,1e-4)/(1+z);
-}
-
-double calc_luminosity_distance(double z)
-{
- //return c/H*integer(f_dist,0,z)/(1+z);
- return c/H*adapt_trapezoid(f_dist,0.,z,1e-4)*(1+z);
-}
-
-
-//EOF
diff --git a/mass_profile/calc_distance.h b/mass_profile/calc_distance.h
deleted file mode 100644
index ba01316..0000000
--- a/mass_profile/calc_distance.h
+++ /dev/null
@@ -1,8 +0,0 @@
-#ifndef CALC_DISTANCE_H
-#define CALC_DISTANCE_H
-extern double calc_angular_distance(double z);
-extern double calc_luminosity_distance(double z);
-extern double E(double z);
-#endif
-
-//EOF
diff --git a/mass_profile/calc_distance_main.cc b/mass_profile/calc_distance_main.cc
deleted file mode 100644
index 28bdb70..0000000
--- a/mass_profile/calc_distance_main.cc
+++ /dev/null
@@ -1,52 +0,0 @@
-#include "calc_distance.h"
-#include <cstdlib>
-#include <cmath>
-#include <iostream>
-using namespace std;
-
-static double cm=1;
-static double s=1;
-static double km=1000*100;
-static double Mpc=3.08568e+24*cm;
-static double kpc=3.08568e+21*cm;
-//static double yr=365.*24.*3600.;
-//static double Gyr=1e9*yr;
-static double H=71.*km/s/Mpc;
-static const double c=299792458.*100.*cm;
-//const double c=3e8*100*cm;
-static const double pi=4*atan(1);
-static const double omega_m=0.27;
-static const double omega_l=0.73;
-static const double arcsec2arc_ratio=1./60/60/180*pi;
-
-
-
-
-int main(int argc,char* argv[])
-{
- if(argc<2)
- {
- cerr<<"Usage:"<<argv[0]<<" z"<<endl;
- exit(-1);
- }
- if(argc==3)
- {
- H=atof(argv[2])*km/s/Mpc;
- }
- double z=atof(argv[1]);
- double d=c/H*calc_angular_distance(z);
- //double age=calc_age(z);
- //cout<<d<<endl;
- cout<<"d_a_cm= "<<d<<" #angular distance in cm"<<endl;
- cout<<"d_a_mpc= "<<d/Mpc<<"#angular distance in Mpc"<<endl;
- cout<<"d_l_cm= "<<(1+z)*(1+z)*d<<" #luminosity distance in cm"<<endl;
- cout<<"d_l_mpc= "<<(1+z)*(1+z)*d/Mpc<<" #luminosity in Mpc"<<endl;
- cout<<"kpc_per_sec= "<<d/kpc*arcsec2arc_ratio<<" #kpc per arcsec"<<endl;
- cout<<"For Chandra:"<<endl;
- cout<<"kpc_per_pixel= "<<d/kpc*0.492*arcsec2arc_ratio<<" #kpc per chandra pixel"<<endl;
- cout<<"cm_per_pixel= "<<d*.492*arcsec2arc_ratio<<" #cm per chandra pixel"<<endl;
- cout<<"norm= "<<1E-14/(4*pi*pow(d*(1+z),2))<<" #norm used to calc cooling function"<<endl;
- //cout<<ddivid(calc_distance,d,0,1,.0001)<<endl;
- cout<<"E(z)= "<<E(z)<<endl;
-}
-
Copyright © 2017-2018 Aaron LI. All Rights Reserved.