aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rw-r--r--README.md7
-rw-r--r--mass_profile/Makefile76
-rw-r--r--mass_profile/chisq.hpp23
-rw-r--r--mass_profile/fit_nfw_sbp.cpp3
-rw-r--r--mass_profile/plot_reporter.cpp70
-rw-r--r--mass_profile/plot_reporter.hpp23
-rw-r--r--mass_profile/vchisq.hpp30
7 files changed, 42 insertions, 190 deletions
diff --git a/README.md b/README.md
index f17d72f..ad6ce69 100644
--- a/README.md
+++ b/README.md
@@ -50,14 +50,13 @@ Installation
------------
Dependencies:
+ GSL
-+ X11 libraries
-+ a working HEASoft installation (for ``libpgplot.a`` and ``libcpgplot.a``)
1. ``mass_profile``
```
$ cd mass_profile
-$ make cleanall
-$ heainit # initilize HEASoft, in order to find libpgplot.a and libcpgplot.a
+$ make clean
+$ make
+(or use this to enable OpenMP)
$ make OPENMP=yes
```
diff --git a/mass_profile/Makefile b/mass_profile/Makefile
index b6418ab..3bf9aa0 100644
--- a/mass_profile/Makefile
+++ b/mass_profile/Makefile
@@ -22,8 +22,6 @@ else
endif
OPT_UTIL_INC ?= -I../opt_utilities
-PGPLOT_INC= -I$(HEADAS)/include
-PGPLOT_LIB= -L. -lcpgplot -lpgplot -lX11 -lgfortran
TARGETS= fit_nfw_sbp fit_dbeta_sbp fit_beta_sbp \
fit_direct_beta calc_distance fit_wang2012_model \
@@ -31,52 +29,36 @@ TARGETS= fit_nfw_sbp fit_dbeta_sbp fit_beta_sbp \
fit_lt_bpl cooling_time calc_lx_dbeta calc_lx_beta
HEADERS= projector.hpp spline.hpp vchisq.hpp
-all: pgplot $(TARGETS)
-
-pgplot: libcpgplot.a libpgplot.a
-
-libcpgplot.a:
-ifdef HEADAS
- ln -s $(HEADAS)/lib/libcpgplot*.a ./$@
-else
- $(error HEASoft not initialized)
-endif
-
-libpgplot.a:
-ifdef HEADAS
- ln -s $(HEADAS)/lib/libpgplot*.a ./$@
-else
- $(error HEASoft not initialized)
-endif
+all: $(TARGETS)
# NOTE:
# Object/source files should placed *before* libraries (order matters)
-fit_nfw_sbp: fit_nfw_sbp.o plot_reporter.o
- $(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC) $(PGPLOT_INC) $(PGPLOT_LIB)
+fit_nfw_sbp: fit_nfw_sbp.o
+ $(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC)
-fit_dbeta_sbp: fit_dbeta_sbp.o beta_cfg.o report_error.o plot_reporter.o
- $(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC) $(PGPLOT_INC) $(PGPLOT_LIB)
+fit_dbeta_sbp: fit_dbeta_sbp.o beta_cfg.o report_error.o
+ $(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC)
-fit_beta_sbp: fit_beta_sbp.o beta_cfg.o dump_fit_qdp.o plot_reporter.o
- $(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC) $(PGPLOT_INC) $(PGPLOT_LIB)
+fit_beta_sbp: fit_beta_sbp.o beta_cfg.o dump_fit_qdp.o
+ $(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC)
-fit_wang2012_model: fit_wang2012_model.o plot_reporter.o
- $(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC) $(PGPLOT_INC) $(PGPLOT_LIB)
+fit_wang2012_model: fit_wang2012_model.o
+ $(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC)
-fit_nfw_mass: fit_nfw_mass.o plot_reporter.o
- $(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC) $(PGPLOT_INC) $(PGPLOT_LIB)
+fit_nfw_mass: fit_nfw_mass.o
+ $(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC)
-fit_direct_beta: fit_direct_beta.o plot_reporter.o
- $(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC) $(PGPLOT_INC) $(PGPLOT_LIB)
+fit_direct_beta: fit_direct_beta.o
+ $(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC)
-calc_lx: calc_lx.o calc_distance.o plot_reporter.o
- $(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC) $(PGPLOT_INC) $(PGPLOT_LIB)
+calc_lx: calc_lx.o calc_distance.o
+ $(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC)
-calc_lx_dbeta: calc_lx_dbeta.o beta_cfg.o report_error.o plot_reporter.o
- $(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC) $(PGPLOT_INC) $(PGPLOT_LIB)
+calc_lx_dbeta: calc_lx_dbeta.o beta_cfg.o report_error.o
+ $(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC)
-calc_lx_beta: calc_lx_beta.o beta_cfg.o dump_fit_qdp.o plot_reporter.o
- $(CXX) $(CXXFLAGS) $^ -o $@ $(OPT_UTIL_INC) $(PGPLOT_INC) $(PGPLOT_LIB)
+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 $@ $^
@@ -98,28 +80,25 @@ cooling_time: cooling_time.cpp
fit_nfw_sbp.o: fit_nfw_sbp.cpp nfw_ne.hpp $(HEADERS)
- $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC) $(PGPLOT_INC)
+ $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC)
fit_dbeta_sbp.o: fit_dbeta_sbp.cpp constrained_dbeta.hpp $(HEADERS)
- $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC) $(PGPLOT_INC)
+ $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC)
fit_beta_sbp.o: fit_beta_sbp.cpp beta.hpp $(HEADERS)
- $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC) $(PGPLOT_INC)
+ $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC)
fit_wang2012_model.o: fit_wang2012_model.cpp wang2012_model.hpp chisq.hpp
- $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC) $(PGPLOT_INC)
+ $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC)
fit_nfw_mass.o: fit_nfw_mass.cpp nfw.hpp chisq.hpp
- $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC) $(PGPLOT_INC)
+ $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC)
calc_lx_dbeta.o: calc_lx_dbeta.cpp constrained_dbeta.hpp $(HEADERS)
- $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC) $(PGPLOT_INC)
+ $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC)
calc_lx_beta.o: calc_lx_beta.cpp beta.hpp $(HEADERS)
- $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC) $(PGPLOT_INC)
-
-plot_reporter.o: plot_reporter.cpp plot_reporter.hpp
- $(CXX) $(CXXFLAGS) -c $< $(PGPLOT_INC)
+ $(CXX) $(CXXFLAGS) -c $< $(OPT_UTIL_INC)
beta_cfg.o: beta_cfg.cpp beta_cfg.hpp
$(CXX) $(CXXFLAGS) -c $<
@@ -142,6 +121,3 @@ calc_distance.o: calc_distance.cc calc_distance.h adapt_trapezoid.h
clean:
rm -f *.o $(TARGETS)
-
-cleanall: clean
- rm -f libpgplot.a libcpgplot.a
diff --git a/mass_profile/chisq.hpp b/mass_profile/chisq.hpp
index bdfadcb..0e58200 100644
--- a/mass_profile/chisq.hpp
+++ b/mass_profile/chisq.hpp
@@ -6,15 +6,17 @@
#ifndef CHI_SQ_HPP
#define CHI_SQ_HPP
+
#define OPT_HEADER
+
#include <core/fitter.hpp>
#include <iostream>
#include <vector>
#include <misc/optvec.hpp>
#include <cmath>
-#include "plot_reporter.hpp"
-#include <cpgplot.h>
-using std::cerr;using std::endl;
+
+using std::cerr;
+using std::endl;
namespace opt_utilities
{
@@ -77,8 +79,6 @@ namespace opt_utilities
:verb(true),limit_bound(false)
{}
-
-
Ty do_eval(const Tp& p)
{
if(limit_bound)
@@ -111,7 +111,6 @@ namespace opt_utilities
vye2.resize(this->get_data_set().size());
my.resize(this->get_data_set().size());
}
-
}
for(int i=(this->get_data_set()).size()-1;i>=0;--i)
@@ -157,7 +156,6 @@ namespace opt_utilities
}
}
-
if(y_model>y_obs)
{
y_err=this->get_data_set().get_data(i).get_y_upper_err();
@@ -185,7 +183,6 @@ namespace opt_utilities
ymax=std::max(vy.at(i),ymax+vye2[i]);
}
-
}
if(verb)
{
@@ -198,20 +195,12 @@ namespace opt_utilities
}
cerr<<endl;
//cerr<<x1<<"\t"<<x2<<endl;
- pr.init_xyrange(xmin,xmax,ymin,ymax,0);
- pr.plot_err2_dot(vx,vy,vye1,vye2);
- pr.plot_line(vx,my);
- cpgask(0);
}
-
}
return result;
}
};
-
-
}
-#endif
-//EOF
+#endif /* CHI_SQ_HPP */
diff --git a/mass_profile/fit_nfw_sbp.cpp b/mass_profile/fit_nfw_sbp.cpp
index 35e42a6..9ac8401 100644
--- a/mass_profile/fit_nfw_sbp.cpp
+++ b/mass_profile/fit_nfw_sbp.cpp
@@ -15,9 +15,10 @@
#include <core/freeze_param.hpp>
#include <error_estimator/error_estimator.hpp>
#include "spline.hpp"
-#include <cpgplot.h>
+
using namespace std;
using namespace opt_utilities;
+
//double s=5.63136645E20;
const double M_sun=1.988E33;//solar mass in g
const double kpc=3.086E21;//kpc in cm
diff --git a/mass_profile/plot_reporter.cpp b/mass_profile/plot_reporter.cpp
deleted file mode 100644
index b5fc9ff..0000000
--- a/mass_profile/plot_reporter.cpp
+++ /dev/null
@@ -1,70 +0,0 @@
-#include "plot_reporter.hpp"
-#include <cpgplot.h>
-#include <cassert>
-#include <cstdlib>
-
-plot_reporter::plot_reporter()
-{
- const char* pgplot_device=getenv("PGPLOT_DEVICE");
- if(pgplot_device==NULL)
- {
- if (cpgopen("/null") < 1)
- {
- assert(0);
- }
- }
- else
- {
- if (cpgopen(pgplot_device) < 1)
- {
- assert(0);
- }
- }
- cpgask(0);
-}
-
-
-plot_reporter::~plot_reporter()
-{
- cpgclos();
-}
-
-
-void plot_reporter::init_xyrange(float x1,
- float x2,
- float y1,
- float y2,
- int axis_flag)
-{
- cpgenv(x1, x2, y1, y2, 0, axis_flag);
-}
-
-
-void plot_reporter::plot_line(std::vector<float>& x,std::vector<float>& y)
-{
- cpgbbuf();
- cpgline(x.size(),x.data(),y.data());
- cpgebuf();
-}
-
-void plot_reporter::plot_err1_dot(std::vector<float>& x,std::vector<float>& y,
- std::vector<float>& e)
-{
- cpgbbuf();
- cpgpt(x.size(),x.data(),y.data(),1);
- cpgerrb(6,x.size(),x.data(),y.data(),e.data(),0);
- cpgebuf();
-}
-
-
-void plot_reporter::plot_err2_dot(std::vector<float>& x,std::vector<float>& y,
- std::vector<float>& e1,std::vector<float>& e2)
-{
- cpgbbuf();
- cpgpt(x.size(),x.data(),y.data(),1);
- cpgerrb(2,x.size(),x.data(),y.data(),e1.data(),0);
- cpgerrb(4,x.size(),x.data(),y.data(),e2.data(),0);
- cpgebuf();
-}
-
-plot_reporter pr;
diff --git a/mass_profile/plot_reporter.hpp b/mass_profile/plot_reporter.hpp
deleted file mode 100644
index 5625f47..0000000
--- a/mass_profile/plot_reporter.hpp
+++ /dev/null
@@ -1,23 +0,0 @@
-#ifndef PLOT_REPORTER_HPP
-#define PLOT_REPORTER_HPP
-#include <vector>
-class plot_reporter
-{
-private:
- plot_reporter(const plot_reporter&);
- plot_reporter& operator=(const plot_reporter&);
-public:
- plot_reporter();
- ~plot_reporter();
- void init_xyrange(float x1,float x2,float y1,float y2,int axis_flag);
- void plot_line(std::vector<float>& x,std::vector<float>& y);
- void plot_err1_dot(std::vector<float>& x,std::vector<float>& y,
- std::vector<float>& e);
- void plot_err2_dot(std::vector<float>& x,std::vector<float>& y,
- std::vector<float>& e1,std::vector<float>& e2);
-
-};
-
-extern plot_reporter pr;
-
-#endif
diff --git a/mass_profile/vchisq.hpp b/mass_profile/vchisq.hpp
index c3d38b8..8cd8481 100644
--- a/mass_profile/vchisq.hpp
+++ b/mass_profile/vchisq.hpp
@@ -6,19 +6,20 @@
#ifndef VCHI_SQ_HPP
#define VCHI_SQ_HPP
+
#define OPT_HEADER
+
#include <core/fitter.hpp>
#include <iostream>
#include <vector>
#include <misc/optvec.hpp>
#include <cmath>
-#include "plot_reporter.hpp"
-#include <cpgplot.h>
-using std::cerr;using std::endl;
+
+using std::cerr;
+using std::endl;
namespace opt_utilities
{
-
template<typename T>
class vchisq
:public statistic<std::vector<T>,std::vector<T>,std::vector<T>,T,std::string>
@@ -59,8 +60,6 @@ namespace opt_utilities
:verb(false),limit_bound(false)
{}
-
-
T do_eval(const std::vector<T>& p)
{
if(limit_bound)
@@ -71,7 +70,6 @@ namespace opt_utilities
}
}
T result(0);
-
std::vector<float> vx;
std::vector<float> vy;
std::vector<float> vye;
@@ -80,7 +78,6 @@ namespace opt_utilities
if(verb)
{
n++;
-
if(n%100==0)
{
vx.resize(this->get_data_set().get_data(0).get_y().size());
@@ -88,7 +85,6 @@ namespace opt_utilities
vye.resize(this->get_data_set().get_data(0).get_y().size());
my.resize(this->get_data_set().get_data(0).get_y().size());
}
-
}
for(int i=(this->get_data_set()).size()-1;i>=0;--i)
{
@@ -101,10 +97,8 @@ namespace opt_utilities
result+=chi*chi;
}
-
if(verb&&n%100==0)
{
-
for(size_t j=0;j<y.size();++j)
{
vx.at(j)=((this->get_data_set().get_data(i).get_x().at(j)+this->get_data_set().get_data(i).get_x().at(j+1))/2.);
@@ -119,18 +113,13 @@ namespace opt_utilities
vx[j]=log10(vx[j]);
vy[j]=log10(vy[j]);
my[j]=log10(my[j]);
-
}
-
}
-
}
if(verb)
{
-
if(n%100==0)
{
-
cerr<<result<<"\t";
for(size_t i=0;i<get_size(p);++i)
{
@@ -138,15 +127,6 @@ namespace opt_utilities
}
cerr<<endl;
}
- if(n%100==0)
- {
- //cpgask(1);
- //std::cerr<<x1<<"\t"<<y1<<std::endl;
- pr.init_xyrange(log10(x1/1.5),log10(x2*1.5),log10(y1/1.5),log10(y2*1.5),30);
- //pr.init_xyrange((x1),(x2),(y1),(y2));
- pr.plot_err1_dot(vx,vy,vye);
- pr.plot_line(vx,my);
- }
}
return result;
Copyright © 2017-2018 Aaron LI. All Rights Reserved.