aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
-rwxr-xr-xmass_profile/analyze_entropy_profile.py54
-rwxr-xr-xmass_profile/calc_all_entropy.sh39
-rwxr-xr-xmass_profile/calc_entropy.py104
-rwxr-xr-xmass_profile/fit_beta_entropy.sh150
-rwxr-xr-xmass_profile/fit_dbeta_entropy.sh149
5 files changed, 104 insertions, 392 deletions
diff --git a/mass_profile/analyze_entropy_profile.py b/mass_profile/analyze_entropy_profile.py
deleted file mode 100755
index c123681..0000000
--- a/mass_profile/analyze_entropy_profile.py
+++ /dev/null
@@ -1,54 +0,0 @@
-#!/usr/bin/env python
-
-import sys
-import numpy
-import scipy.interpolate
-
-center_entropy_file=open('entropy_center.qdp')
-entropy_file=open('summary_entropy.qdp')
-confidence_level=.68
-rout=float(sys.argv[1])
-
-center_s=0
-for i in center_entropy_file:
- r,s=i.split()
- r=float(r)
- s=float(s)
- if r>rout:
- center_s=s
- break
-
-new_data=True
-
-
-s_list=[]
-for i in entropy_file:
- if i[0]=='n':
- new_data=True
- continue
- if new_data==False:
- continue
- r,s=i.split()
- r=float(r)
- s=float(s)
- if r>rout:
- new_data=False
- s_list.append(s)
-
-s_idx=-1
-
-s_list.sort()
-for i in range(len(s_list)-1):
- if (center_s-s_list[i])*(center_s-s_list[i+1])<=0:
- m_idx=i
- break
-
-
-slidx=int(s_idx*(1-confidence_level))
-suidx=s_idx-1+int((len(s_list)-s_idx)*confidence_level)
-
-
-serr1=s_list[slidx]-center_s
-serr2=s_list[suidx]-center_s
-
-print("S=\t%e\t %e/+%e keV cm^2 (1 sigma)"%(center_s,serr1,serr2))
diff --git a/mass_profile/calc_all_entropy.sh b/mass_profile/calc_all_entropy.sh
deleted file mode 100755
index 5506ab8..0000000
--- a/mass_profile/calc_all_entropy.sh
+++ /dev/null
@@ -1,39 +0,0 @@
-#!/usr/bin/env bash
-
-if [ $# -eq 1 ]
-then
- :
-else
- echo "Usage: <global.cfg list file>"
- exit
-fi
-
-bdir=`pwd`
-
-for i in `cat $1`
-do
- dname=`dirname $i`
- #dname=`dirname $dname`
- echo $dname
- cd $dname ||continue
-# mkdir -p profile_entropy
-# cp -rf profile/* profile_entropy/
-# cd profile/
-
- [ -e result ] && r200file=result
- [ -e result ] || r200file=result_checked
-
- r200=`grep r200 $r200file|awk -F = '{print $2}' |awk -F '-' '{print $1}'`
- rout=`python -c "print($r200*.1)"`
- echo $rout
- if grep n01 source.cfg >/dev/null
- then
- echo "dbeta"
- fit_dbeta_entropy.sh global.cfg $rout
- else
- echo "beta"
- fit_beta_entropy.sh global.cfg $rout
- fi
-
- cd $bdir
-done
diff --git a/mass_profile/calc_entropy.py b/mass_profile/calc_entropy.py
new file mode 100755
index 0000000..736fe7e
--- /dev/null
+++ b/mass_profile/calc_entropy.py
@@ -0,0 +1,104 @@
+#!/usr/bin/env python3
+#
+# Calculate the entropy within the specified radius.
+#
+# Junhua GU
+# Weitian LI
+# 2016-06-07
+#
+
+import argparse
+import re
+from itertools import groupby
+import numpy as np
+
+
+def isplit(iterable, splitters):
+ """
+ Credit: https://stackoverflow.com/a/4322780/4856091
+ """
+ return [list(g) for k, g in groupby(iterable,
+ lambda x:x in splitters) if not k]
+
+
+def get_entropy(data, r):
+ """
+ Get the entropy *at* the specified radius.
+
+ XXX: whether to interpolate first?
+ """
+ radius = data[:, 0]
+ entropy = data[:, 1]
+ s = np.min(entropy[radius > r])
+ return s
+
+
+def read_merged_qdp(infile):
+ """
+ Read merged QDP with multiple group of data separated by "no no no".
+ """
+ lines = map(lambda line: re.sub(r"^\s*no\s+no\s+no.*$", "X",
+ line.strip(), flags=re.I),
+ open(infile).readlines())
+ lines = isplit(lines, ("X",))
+ data_groups = []
+ for block in lines:
+ data = [list(map(float, l.split())) for l in block]
+ data.append(np.row_stack(data))
+ return data_groups
+
+
+def calc_error(center_value, mc_values, ci=0.683):
+ """
+ Calculate the uncertainties/errors.
+ """
+ data = np.concatenate([[center_value], mc_values])
+ median, q_lower, q_upper = np.percentile(data, q=(50, 50-50*ci, 50+50*ci))
+ mean = np.mean(data)
+ std = np.std(data)
+ return {
+ "mean": mean,
+ "std": std,
+ "median": median,
+ "q_lower": q_lower,
+ "q_upper": q_upper,
+ }
+
+
+def main():
+ parser = argparse.ArgumentParser(
+ description="Calculate the entropy within the given radius")
+ parser.add_argument("-C", "--confidence-level", dest="ci",
+ type=float, default=0.683,
+ help="confidence level to estimate the errors")
+ parser.add_argument("center_data",
+ help="calculate central entropy profile " +
+ "(e.g., entropy_center.qdp)")
+ parser.add_argument("mc_data",
+ help="Merged QDP file of all the Monte Carlo " +
+ "simulated entropy profiles " +
+ "(e.g., summary_entropy.qdp)")
+ parser.add_argument("rout", type=float, help="outer radius (kpc)")
+ args = parser.parse_args()
+
+ center_data = np.loadtxt(args.center_data)
+ center_s = get_entropy(center_data, r=args.rout)
+
+ data_groups = read_merged_qdp(args.mc_data)
+ entropy_list = []
+ for dg in data_groups:
+ s = get_entropy(dg, r=args.rout)
+ entropy_list.append(s)
+ results = calc_error(center_s, entropy_list, ci=args.ci)
+ s_err_lower = results["q_lower"] - center_s
+ s_err_upper = results["q_upper"] - center_s
+
+ print("entropy= %e %+e/%+e keV cm^2 (ci=%.1f%%)" %
+ (center_s, s_err_lower, s_err_upper, args.ci * 100))
+ print("entropy(mean)= %e" % results["mean"])
+ print("entropy(median)= %e" % results["median"])
+ print("entropy(std)= %e" % results["std"])
+
+
+if __name__ == "__main__":
+ main()
diff --git a/mass_profile/fit_beta_entropy.sh b/mass_profile/fit_beta_entropy.sh
deleted file mode 100755
index f52e4a1..0000000
--- a/mass_profile/fit_beta_entropy.sh
+++ /dev/null
@@ -1,150 +0,0 @@
-#!/bin/bash
-
-echo $#
-if [ $# -eq 2 ]
-then
- :
-else
- echo "Usage:$0 <cfg file> <radius in kpc>"
- exit
-fi
-export PGPLOT_FONT=`locate grfont.dat|head -1`
-
-mkdir -p entropy_files
-
-cfg_file=$1
-base_path=`dirname $0`
-echo $base_path
-#initialize profile type name
-t_profile_type=`grep t_profile $cfg_file|awk '{print $2}'`
-#initialize data file name
-t_data_file=`grep t_data_file $cfg_file|awk '{print $2}'`
-#initialize sbp config file
-sbp_cfg=`grep sbp_cfg $cfg_file|awk '{print $2}'`
-#initialize the temperature profile file
-T_file=`grep '^T_file' $sbp_cfg|awk '{print $2}'`
-#echo $t_profile_type
-cm_per_pixel=`grep '^cm_per_pixel' $sbp_cfg|awk '{print $2}'`
-
-#determine which temperature profile to be used, and fit the T profile
-if [ $t_profile_type == zyy ]
-then
- t_param_file=`grep t_param_file $cfg_file|awk '{print $2}'`
- $base_path/fit_zyy_model $t_data_file $t_param_file $cm_per_pixel
- mv -f zyy_dump.qdp ${T_file}
-elif [ $t_profile_type == m0603246 ]
-then
- $base_path/fit_m0603246 $t_data_file $cm_per_pixel
- mv -f m0603246_dump.qdp ${T_file}
-elif [ $t_profile_type == wang2012 ]
-then
- t_param_file=`grep t_param_file $cfg_file|awk '{print $2}'`
- $base_path/fit_wang2012_model $t_data_file $t_param_file $cm_per_pixel
- mv -f wang2012_dump.qdp ${T_file}
-elif [ $t_profile_type == allen ]
-then
- $base_path/fit_allen_model $t_data_file $cm_per_pixel
- mv -f allen_dump.qdp ${T_file}
-elif [ $t_profile_type == zzl ]
-then
- t_param_file=`grep t_param_file $cfg_file|awk '{print $2}'`
- $base_path/fit_zzl_model $t_data_file $t_param_file
- mv -f zzl_dump.qdp ${T_file}
-else
- echo temperature profile name invalid!
- exit
-fi
-
-cfunc_file=`grep '^cfunc_file' ${sbp_cfg} |awk '{print $2}'`
-z=`grep '^z' ${sbp_cfg}|awk '{print $2}'`
-abund=`grep '^abund' ${cfg_file} |awk '{print $2}'`
-nh=`grep '^nh' ${cfg_file} |awk '{print $2}'`
-$base_path/coolfunc_calc.sh ${T_file} $abund $nh $z $cfunc_file
-mv flux_cnt_ratio.txt flux_cnt_ratio_center_entropy.txt
-#fit sbp
-$base_path/fit_beta_sbp $sbp_cfg
-echo $cfunc_file
-#exit
-
-#store central valu
-mv sbp_fit.qdp sbp_fit_center_entropy.qdp
-mv mass_int.qdp mass_int_center_entropy.qdp
-mv overdensity.qdp overdensity_center_entropy.qdp
-mv gas_mass_int.qdp gas_mass_int_center_entropy.qdp
-mv entropy.qdp entropy_center.qdp
-sbp_data_file=`grep sbp_file $sbp_cfg|awk '{print $2}'`
-radius_sbp_file=`grep radius_sbp_file ${cfg_file}|awk '{print $2}'`
-
-if [ x"$radius_sbp_file" == x ]
-then
- echo "Error, must have radius_sbp_file assigned, this file should be a 4-column file, which contains the radius, radius err, sbp, and sbp err"
- exit
-fi
-
-cat ${radius_sbp_file} | sed 's/#.*$//' | grep -Ev '^\s*$' > .tmp.txt
-mv .tmp.txt ${radius_sbp_file}
-
-rm -f summary_entropy.qdp
-#100 times of Monte-carlo simulation to determine error
-#just repeat above steps
-for i in `seq 1 100`
-do
- echo $t_data_file
- $base_path/shuffle_T.py $t_data_file temp_shuffled_t.dat
- $base_path/shuffle_sbp.py $sbp_data_file temp_shuffled_sbp.dat
-
-#exit
- if [ $t_profile_type == zyy ]
- then
- t_param_file=`grep t_param_file $cfg_file|awk '{print $2}'`
- $base_path/fit_zyy_model temp_shuffled_t.dat $t_param_file $cm_per_pixel
- mv -f zyy_dump.qdp ${T_file}
- elif [ $t_profile_type == m0603246 ]
- then
- $base_path/fit_m0603246 temp_shuffled_t.dat $cm_per_pixel
- mv -f m0603246_dump.qdp ${T_file}
- elif [ $t_profile_type == wang2012 ]
- then
- t_param_file=`grep t_param_file $cfg_file|awk '{print $2}'`
- $base_path/fit_wang2012_model temp_shuffled_t.dat $t_param_file $cm_per_pixel
- mv -f wang2012_dump.qdp ${T_file}
- elif [ $t_profile_type == allen ]
- then
- $base_path/fit_allen_model temp_shuffled_t.dat $cm_per_pixel
- mv -f allen_dump.qdp ${T_file}
- elif [ $t_profile_type == zzl ]
- then
- t_param_file=`grep t_param_file $cfg_file|awk '{print $2}'`
- $base_path/fit_zzl_model temp_shuffled_t.dat $t_param_file
- mv -f zzl_dump.qdp ${T_file}
- else
- echo temperature profile name invalid!
- exit
- fi
-
-#exit
- echo >temp_sbp.cfg
-
- cat $sbp_cfg|while read l
-do
- if echo $l|grep sbp_file >/dev/null
- then
- echo sbp_file temp_shuffled_sbp.dat >>temp_sbp.cfg
- elif echo $l|grep T_file >/dev/null
- then
- echo T_file ${T_file} >>temp_sbp.cfg
- else
- echo $l >>temp_sbp.cfg
- fi
-
-done
-
-#$base_path/coolfunc_calc.sh ${T_file} $abund $nh $z $cfunc_file
-
-$base_path/fit_beta_sbp temp_sbp.cfg
-
-cat entropy.qdp >>summary_entropy.qdp
-echo no no no >>summary_entropy.qdp
-done
-
-$base_path/analyze_entropy_profile.py $2 |tee entropy_result.txt
diff --git a/mass_profile/fit_dbeta_entropy.sh b/mass_profile/fit_dbeta_entropy.sh
deleted file mode 100755
index 3f1bbda..0000000
--- a/mass_profile/fit_dbeta_entropy.sh
+++ /dev/null
@@ -1,149 +0,0 @@
-#!/bin/bash
-
-echo $#
-if [ $# -eq 2 ]
-then
- :
-else
- echo "Usage:$0 <cfg file> <radius in kpc>"
- exit
-fi
-export PGPLOT_FONT=`locate grfont.dat|head -1`
-
-cfg_file=$1
-base_path=`dirname $0`
-echo $base_path
-#initialize profile type name
-t_profile_type=`grep t_profile $cfg_file|awk '{print $2}'`
-#initialize data file name
-t_data_file=`grep t_data_file $cfg_file|awk '{print $2}'`
-#initialize sbp config file
-sbp_cfg=`grep sbp_cfg $cfg_file|awk '{print $2}'`
-#initialize the temperature profile file
-T_file=`grep '^T_file' $sbp_cfg|awk '{print $2}'`
-#echo $t_profile_type
-cm_per_pixel=`grep '^cm_per_pixel' $sbp_cfg|awk '{print $2}'`
-#determine which temperature profile to be used, and fit the T profile
-if [ $t_profile_type == zyy ]
-then
- t_param_file=`grep t_param_file $cfg_file|awk '{print $2}'`
- $base_path/fit_zyy_model $t_data_file $t_param_file $cm_per_pixel
- mv -f zyy_dump.qdp ${T_file}
-elif [ $t_profile_type == m0603246 ]
-then
- $base_path/fit_m0603246 $t_data_file $cm_per_pixel
- mv -f m0603246_dump.qdp ${T_file}
-elif [ $t_profile_type == wang2012 ]
-then
- t_param_file=`grep t_param_file $cfg_file|awk '{print $2}'`
- $base_path/fit_wang2012_model $t_data_file $t_param_file $cm_per_pixel
- mv -f wang2012_dump.qdp ${T_file}
-elif [ $t_profile_type == allen ]
-then
- $base_path/fit_allen_model $t_data_file $cm_per_pixel
- mv -f allen_dump.qdp ${T_file}
-elif [ $t_profile_type == zzl ]
-then
- t_param_file=`grep t_param_file $cfg_file|awk '{print $2}'`
- $base_path/fit_zzl_model $t_data_file $t_param_file
- mv -f zzl_dump.qdp ${T_file}
-else
- echo temperature profile name invalid!
- exit
-fi
-
-cfunc_file=`grep '^cfunc_file' ${sbp_cfg} |awk '{print $2}'`
-z=`grep '^z' ${sbp_cfg}|awk '{print $2}'`
-abund=`grep '^abund' ${cfg_file} |awk '{print $2}'`
-nh=`grep '^nh' ${cfg_file} |awk '{print $2}'`
-$base_path/coolfunc_calc.sh ${T_file} $abund $nh $z $cfunc_file
-mv flux_cnt_ratio.txt flux_cnt_ratio_center_entropy.txt
-#fit sbp
-$base_path/fit_dbeta_sbp $sbp_cfg
-echo $cfunc_file
-#exit
-
-#store central valu
-mv sbp_fit.qdp sbp_fit_center_entropy.qdp
-mv mass_int.qdp mass_int_center_entropy.qdp
-mv overdensity.qdp overdensity_center_entropy.qdp
-mv gas_mass_int.qdp gas_mass_int_center_entropy.qdp
-mv entropy.qdp entropy_center.qdp
-sbp_data_file=`grep sbp_file $sbp_cfg|awk '{print $2}'`
-radius_sbp_file=`grep radius_sbp_file ${cfg_file}|awk '{print $2}'`
-
-if [ x"$radius_sbp_file" == x ]
-then
- echo "Error, must have radius_sbp_file assigned, this file should be a 4-column file, which contains the radius, radius err, sbp, and sbp err"
- exit
-fi
-
-cat ${radius_sbp_file} | sed 's/#.*$//' | grep -Ev '^\s*$' > .tmp.txt
-mv .tmp.txt ${radius_sbp_file}
-
-rm -f summary_entropy.qdp
-
-#100 times of Monte-carlo simulation to determine error
-#just repeat above steps
-for i in `seq 1 100`
-do
- echo $t_data_file
- $base_path/shuffle_T.py $t_data_file temp_shuffled_t.dat
- $base_path/shuffle_sbp.py $sbp_data_file temp_shuffled_sbp.dat
-
-#exit
- if [ $t_profile_type == zyy ]
- then
- t_param_file=`grep t_param_file $cfg_file|awk '{print $2}'`
- $base_path/fit_zyy_model temp_shuffled_t.dat $t_param_file $cm_per_pixel
- mv -f zyy_dump.qdp ${T_file}
- elif [ $t_profile_type == m0603246 ]
- then
- $base_path/fit_m0603246 temp_shuffled_t.dat $cm_per_pixel
- mv -f m0603246_dump.qdp ${T_file}
- elif [ $t_profile_type == wang2012 ]
- then
- t_param_file=`grep t_param_file $cfg_file|awk '{print $2}'`
- $base_path/fit_wang2012_model temp_shuffled_t.dat $t_param_file $cm_per_pixel
- mv -f wang2012_dump.qdp ${T_file}
- elif [ $t_profile_type == allen ]
- then
- $base_path/fit_allen_model temp_shuffled_t.dat $cm_per_pixel
- mv -f allen_dump.qdp ${T_file}
- elif [ $t_profile_type == zzl ]
- then
- t_param_file=`grep t_param_file $cfg_file|awk '{print $2}'`
- $base_path/fit_zzl_model temp_shuffled_t.dat $t_param_file
- mv -f zzl_dump.qdp ${T_file}
- else
- echo temperature profile name invalid!
- exit
- fi
-
-#exit
- echo >temp_sbp.cfg
-
- cat $sbp_cfg|while read l
-do
- if echo $l|grep sbp_file >/dev/null
- then
- echo sbp_file temp_shuffled_sbp.dat >>temp_sbp.cfg
- elif echo $l|grep T_file >/dev/null
- then
- echo T_file ${T_file} >>temp_sbp.cfg
- else
- echo $l >>temp_sbp.cfg
- fi
-
-done
-
-#$base_path/coolfunc_calc.sh ${T_file} $abund $nh $z $cfunc_file
-
-$base_path/fit_dbeta_sbp temp_sbp.cfg
-
-cat entropy.qdp >>summary_entropy.qdp
-echo no no no >>summary_entropy.qdp
-done
-
-$base_path/analyze_entropy_profile.py $2 |tee entropy_result.txt
-