aboutsummaryrefslogtreecommitdiffstats
path: root/mass_profile
diff options
context:
space:
mode:
Diffstat (limited to 'mass_profile')
-rwxr-xr-xmass_profile/analyze_lxfx.py57
-rwxr-xr-xmass_profile/analyze_mass_profile.py195
-rwxr-xr-xmass_profile/calc_coolfunc.sh113
-rwxr-xr-xmass_profile/calc_coolfunc_bands.sh119
-rwxr-xr-xmass_profile/calc_entropy.py104
-rwxr-xr-xmass_profile/calc_lxfx.sh162
-rwxr-xr-xmass_profile/calc_lxfx_wrapper.sh76
-rwxr-xr-xmass_profile/fg_2500_500.py153
-rwxr-xr-xmass_profile/fit_mass.sh240
-rwxr-xr-xmass_profile/fit_sbp.sh61
-rwxr-xr-xmass_profile/get_lxfx_data.sh89
-rwxr-xr-xmass_profile/shuffle_profile.py37
12 files changed, 0 insertions, 1406 deletions
diff --git a/mass_profile/analyze_lxfx.py b/mass_profile/analyze_lxfx.py
deleted file mode 100755
index 62859ff..0000000
--- a/mass_profile/analyze_lxfx.py
+++ /dev/null
@@ -1,57 +0,0 @@
-#!/usr/bin/env python3
-#
-# Calculate the mean and standard deviation of the (Monte Carlo)
-# Lx and Fx data
-#
-# Junhua GU
-# Weitian LI
-# 2016-06-07
-#
-
-import sys
-import argparse
-import numpy as np
-
-
-def read_bands(bands):
- """
- Read energy bands list, each band per line.
- """
- bands = map(str.split, open(bands).readlines())
- bands = ["-".join(b) for b in bands]
- return bands
-
-
-def output(name, bands, means, sigmas, outfile=sys.stdout):
- if outfile is not sys.stdout:
- outfile = open(outfile, "w")
- for b, m, s in zip(bands, means, sigmas):
- print("%s(%s)= %4.2E +/- %4.2E erg/s" % (name, b, m, s),
- file=outfile)
- if outfile is not sys.stdout:
- outfile.close()
-
-
-def main():
- parser = argparse.ArgumentParser(
- description="Analyze Lx/Fx results")
- parser.add_argument("name", help="Lx or Fx")
- parser.add_argument("infile", help="input data file")
- parser.add_argument("outfile", help="output results file")
- parser.add_argument("bands", help="energy bands of the input data columns")
- args = parser.parse_args()
-
- data = np.loadtxt(args.infile)
- bands = read_bands(args.bands)
- if len(bands) != data.shape[1]:
- raise ValueError("number of energy bands != number of data columns")
-
- means = np.mean(data, axis=0)
- sigmas = np.std(data, axis=0)
- output(name=args.name, bands=bands, means=means, sigmas=sigmas)
- output(name=args.name, bands=bands, means=means, sigmas=sigmas,
- outfile=args.outfile)
-
-
-if __name__ == "__main__":
- main()
diff --git a/mass_profile/analyze_mass_profile.py b/mass_profile/analyze_mass_profile.py
deleted file mode 100755
index 3ee128c..0000000
--- a/mass_profile/analyze_mass_profile.py
+++ /dev/null
@@ -1,195 +0,0 @@
-#!/usr/bin/env python
-
-import sys
-import numpy
-import scipy.interpolate
-
-confidence_level=.68
-def read_file(param):
- delta=float(param[0])
-
- file_mass_center=open("mass_int_center.qdp").readlines();
- file_delta_center=open("overdensity_center.qdp").readlines();
-
- center_r=0
- center_m=0
- center_gm=0
- center_gf=0
-
-
- for i in range(0,len(file_mass_center)):
- lm=file_mass_center[i].strip();
- ld=file_delta_center[i].strip();
- r,m=lm.split()
- r,d=ld.split()
- r=float(r)
- d=float(d)
- m=float(m)
- if m<1e11:
- continue
- if d<delta:
- center_r=r
- center_m=m
- for j in open("gas_mass_int_center.qdp"):
- rgm,gm=j.strip().split()
- rgm=float(rgm)
- gm=float(gm)
- if rgm>r:
-
- center_gm=gm
- center_gf=gm/m
- break
- break
- if len(param)>1 and param[1]=='c':
- #print("%s(<r%d)=%E solar mass"%("mass",delta,center_m))
- #print("%s%d=%E kpc"%("r",delta,center_r))
- #print("%s(<r%d)=%E solar mass"%("gas mass",delta,center_gm))
- #print("%s(<r%d)=%E"%("gas fraction",delta,center_gf))
- return center_m,center_r,center_gm,center_gf,None,None,None,None
-
-
-#print(center_gm,center_gf)
- file_mass=open('summary_mass_profile.qdp').readlines()
- file_delta=open('summary_overdensity.qdp').readlines()
- file_gm=open('summary_gas_mass_profile.qdp')
-
-
- flag=True
- rlist=[]
- mlist=[]
- gmlist=[]
- gflist=[]
- old_m=0
- invalid_count=0
- for i in range(0,len(file_mass)):
- lm=file_mass[i].strip()
- ld=file_delta[i].strip()
- if lm[0]=='n':
- flag=True
- old_m=0
- continue
- if not flag:
- continue
- r,m=lm.split()
- m=float(m)
- if m<1e12:
- continue
- if m<old_m:
- invalid_count+=1
- flag=False
- continue
- r,d=ld.split()
- d=float(d)
- r=float(r)
-
- if d<delta:
- #print("%s %e"%(d,m))
- mlist.append(m)
- rlist.append(r)
- flag1=True
- while True:
- lgm=file_gm.readline().strip()
- if lgm[0]=='n':
- break
- rgm,gm=lgm.split()
- rgm=float(rgm)
- gm=float(gm)
- if rgm>r and flag1:
- gmlist.append(gm)
-
- flag1=False
- gflist.append(gm/mlist[-1])
- #print(gm,gflist[-1])
- flag=False
- old_m=m
- print("%d abnormal data dropped"%(invalid_count))
-
-
- return center_m,center_r,center_gm,center_gf,mlist,rlist,gmlist,gflist
-#center_m=numpy.mean(mlist)
-#center_r=numpy.mean(rlist)
-
-center_m,center_r,center_gm,center_gf,mlist,rlist,gmlist,gflist=read_file(sys.argv[1:])
-delta=float(sys.argv[1])
-
-if len(sys.argv)>2 and sys.argv[2]=='c':
- print("%s(<r%d)=%E solar mass"%("mass",delta,center_m))
- print("%s%d=%E kpc"%("r",delta,center_r))
- print("%s(<r%d)=%E solar mass"%("gas mass",delta,center_gm))
- print("%s(<r%d)=%E"%("gas fraction",delta,center_gf))
- sys.exit(0)
-
-
-mlist.sort()
-rlist.sort()
-gflist.sort()
-gmlist.sort()
-
-m_idx=-1
-r_idx=-1
-gm_idx=-1
-gf_idx=-1
-delta=float(sys.argv[1])
-for i in range(len(mlist)-1):
- if (center_m-mlist[i])*(center_m-mlist[i+1])<=0:
- m_idx=i
- break
-
-for i in range(len(rlist)-1):
- if (center_r-rlist[i])*(center_r-rlist[i+1])<=0:
- r_idx=i
- break
-
-for i in range(len(gmlist)-1):
- if (center_gm-gmlist[i])*(center_gm-gmlist[i+1])<=0:
- gm_idx=i
- break
-
-for i in range(len(gflist)-1):
- if (center_gf-gflist[i])*(center_gf-gflist[i+1])<=0:
- gf_idx=i
- break
-
-
-if m_idx==-1 or r_idx==-1 or gf_idx==-1 or gm_idx==-1:
- print("Error, the center value is not enclosed by the Monte-Carlo realizations, please check the result!")
- print("m:%E %E %E"%(center_m,mlist[0],mlist[-1]))
- print("gm:%E %E %E"%(center_gm,gmlist[0],gmlist[-1]))
- print("gf:%E %E %E"%(center_gf,gflist[0],gflist[-1]))
- print("r:%E %E %E"%(center_r,rlist[0],rlist[-1]))
- sys.exit(1)
-
-
-mlidx=int(m_idx*(1-confidence_level))
-muidx=m_idx-1+int((len(mlist)-m_idx)*confidence_level)
-
-
-rlidx=int(r_idx*(1-confidence_level))
-ruidx=r_idx-1+int((len(rlist)-r_idx)*confidence_level)
-
-gmlidx=int(gm_idx*(1-confidence_level))
-gmuidx=gm_idx-1+int((len(gmlist)-gm_idx)*confidence_level)
-
-gflidx=int(gf_idx*(1-confidence_level))
-gfuidx=gf_idx-1+int((len(gflist)-gf_idx)*confidence_level)
-
-
-merr1=mlist[mlidx]-center_m
-merr2=mlist[muidx]-center_m
-
-rerr1=rlist[rlidx]-center_r
-rerr2=rlist[ruidx]-center_r
-
-gmerr1=gmlist[gmlidx]-center_gm
-gmerr2=gmlist[gmuidx]-center_gm
-
-gferr1=gflist[gflidx]-center_gf
-gferr2=gflist[gfuidx]-center_gf
-
-#print("%d %d %d"%(mlidx,m_idx,muidx))
-#print("%d %d %d"%(rlidx,r_idx,ruidx))
-
-print("m%d=\t%e\t %e/+%e solar mass (1 sigma)"%(delta,center_m,merr1,merr2))
-print("gas_m%d=\t%e\t %e/+%e solar mass (1 sigma)"%(delta,center_gm,gmerr1,gmerr2))
-print("gas_fraction%d=\t%e\t %e/+%e (1 sigma)"%(delta,center_gf,gferr1,gferr2))
-print("r%d=\t%d\t %d/+%d kpc (1 sigma)"%(delta,center_r,rerr1,rerr2))
diff --git a/mass_profile/calc_coolfunc.sh b/mass_profile/calc_coolfunc.sh
deleted file mode 100755
index d1d0b35..0000000
--- a/mass_profile/calc_coolfunc.sh
+++ /dev/null
@@ -1,113 +0,0 @@
-#!/bin/sh
-##
-## Calculate the 'cooling function' profile with respect to the
-## given 'temperature profile' and the average abundance, redshift,
-## and column density nH, using the XSPEC model 'wabs*apec'.
-##
-## NOTE:
-## The output cooling function values should be the 'flux' values
-## with unit 'photon/s/cm^2' (different to 'calc_coolfunc_bands.sh').
-## These results will be used by 'fit_{beta,dbeta}_sbp' to derive the
-## (3D) gas density profile from (2D) surface brightness profile,
-## whose values have unit 'photon/cm^2/pixel/s'.
-##
-## Weitian LI
-## Created: 2012-08-17
-##
-## Change logs:
-## 2017-02-17, Weitian LI
-## * Rename from 'coolfunc_calc.sh' to 'calc_coolfunc.sh'
-## * Clean up that do not calculate and output <coolfunc_bolo>
-## 2017-02-16, Weitian LI
-## * Do not calculate and output 'flux_cnt_ratio.txt'
-##
-
-## cmdline arguments {{{
-if [ $# -ne 5 ]; then
- printf "usage:\n"
- printf " `basename $0` <tprofile> <avg_abund> <nH> <redshift> <coolfunc_outfile>\n"
- exit 1
-fi
-TPROFILE=$1
-ABUNDANCE=$2
-N_H=$3
-REDSHIFT=$4
-COOLFUNC_DAT=$5
-NORM=`cosmo_calc ${REDSHIFT} | grep 'norm.*cooling_function' | awk -F':' '{ print $2 }'`
-
-if [ ! -r "${TPROFILE}" ]; then
- printf "ERROR: given tprofile '${TPROFILE}' NOT accessiable\n"
- exit 2
-fi
-[ -e "${COOLFUNC_DAT}" ] && rm -f ${COOLFUNC_DAT}
-## arguments }}}
-
-XSPEC_CF_XCM="_calc_coolfunc.xcm"
-[ -e "${XSPEC_CF_XCM}" ] && rm -f ${XSPEC_CF_XCM}
-
-## generate xspec script {{{
-cat > ${XSPEC_CF_XCM} << _EOF_
-## XSPEC Tcl script
-## Calculate the cooling function profile w.r.t the temperature profile.
-##
-## Generated by: `basename $0`
-## Date: `date`
-
-set xs_return_results 1
-set xs_echo_script 0
-# set tcl_precision 12
-## set basic data {{{
-set nh ${N_H}
-set redshift ${REDSHIFT}
-set abundance ${ABUNDANCE}
-set norm ${NORM}
-## basic }}}
-
-## xspec related {{{
-# debug settings {{{
-chatter 0
-# debug }}}
-query yes
-abund grsa
-dummyrsp 0.01 100.0 4096 linear
-# load model 'wabs*apec' to calc cooling function
-model wabs*apec & \${nh} & 1.0 & \${abundance} & \${redshift} & \${norm} & /*
-## xspec }}}
-
-## set input and output filename & open files
-set tpro_fn "${TPROFILE}"
-set cf_fn "${COOLFUNC_DAT}"
-if { [ file exists \${cf_fn} ] } {
- exec rm -fv \${cf_fn}
-}
-
-## open files
-set tpro_fd [ open \${tpro_fn} r ]
-set cf_fd [ open \${cf_fn} w ]
-
-## read data from tprofile line by line
-while { [ gets \${tpro_fd} tpro_line ] != -1 } {
- scan \${tpro_line} "%f %f" radius temperature
- #puts "radius: \${radius}, temperature: \${temperature}"
- # set temperature value
- newpar 2 \${temperature}
- # calc flux & tclout
- flux 0.7 7.0
- tclout flux 1
- scan \${xspec_tclout} "%f %f %f %f" _ _ _ cf_photon
- #puts "cf: \${cf_photon}"
- puts \${cf_fd} "\${radius} \${cf_photon}"
-}
-
-## close opened files
-close \${tpro_fd}
-close \${cf_fd}
-
-## exit
-tclexit
-_EOF_
-## xcm generation }}}
-
-## invoke xspec to calc
-printf "invoking XSPEC to calculate cooling function profile ...\n"
-xspec - ${XSPEC_CF_XCM} > /dev/null
diff --git a/mass_profile/calc_coolfunc_bands.sh b/mass_profile/calc_coolfunc_bands.sh
deleted file mode 100755
index bebdce2..0000000
--- a/mass_profile/calc_coolfunc_bands.sh
+++ /dev/null
@@ -1,119 +0,0 @@
-#!/bin/sh
-##
-## Calculate the 'cooling function' profile for each of the energy band
-## specified in a bands file, with respect to the given 'temperature profile'
-## and the average abundance, redshift, and column density nH, using the
-## XSPEC model 'wabs*apec'.
-##
-## NOTE:
-## To calculate the luminosity and flux from the source using the
-## 'calc_lx_{beta,dbeta}', set 'nH=0'.
-## Also the output cooling function values should be the 'flux' values
-## with unit 'erg/s/cm^2'.
-##
-## Weitian LI
-## Updated: 2017-02-17
-##
-
-## cmdline arguments {{{
-if [ $# -ne 6 ]; then
- printf "usage:\n"
- printf " `basename $0` <tprofile> <avg_abund> <nH> <redshift> <coolfunc_prefix> <band_list>\n"
- exit 1
-fi
-TPROFILE=$1
-ABUNDANCE=$2
-N_H=$3
-REDSHIFT=$4
-COOLFUNC_PREFIX=$5
-BLIST=$6
-NORM=`cosmo_calc ${REDSHIFT} | grep 'norm.*cooling_function' | awk -F':' '{ print $2 }'`
-
-if [ ! -r "${TPROFILE}" ]; then
- printf "ERROR: given tprofile '${TPROFILE}' NOT accessiable\n"
- exit 2
-fi
-## arguments }}}
-
-XSPEC_CF_XCM="_calc_coolfunc_bands.xcm"
-[ -e "${XSPEC_CF_XCM}" ] && rm -f ${XSPEC_CF_XCM}
-
-## generate xspec script {{{
-cat > ${XSPEC_CF_XCM} << _EOF_
-## XSPEC Tcl script
-## Calculate the cooling function profile w.r.t the temperature profile,
-## for each specified energy band.
-##
-## Generated by: `basename $0`
-## Date: `date`
-
-set xs_return_results 1
-set xs_echo_script 0
-# set tcl_precision 12
-## set basic data {{{
-set nh ${N_H}
-set redshift ${REDSHIFT}
-set abundance ${ABUNDANCE}
-set norm ${NORM}
-## basic }}}
-
-## xspec related {{{
-# debug settings {{{
-chatter 0
-# debug }}}
-query yes
-abund grsa
-dummyrsp 0.01 100.0 4096 linear
-# load model 'wabs*apec' to calc cooling function
-model wabs*apec & \${nh} & 1.0 & \${abundance} & \${redshift} & \${norm} &
-## xspec }}}
-
-## set input and output filename
-set tpro_fn "${TPROFILE}"
-set blist_fn "${BLIST}"
-set cf_prefix "${COOLFUNC_PREFIX}"
-set blist_fd [ open \${blist_fn} r ]
-
-## loop over each energy band
-while { [ gets \${blist_fd} blist_line ] != -1 } {
- if { "\${blist_line}" == "bolo" } {
- set e1 0.01
- set e2 100.0
- set name_suffix bolo
- } else {
- set e1 [ lindex \${blist_line} 0 ]
- set e2 [ lindex \${blist_line} 1 ]
- set name_suffix "\${e1}-\${e2}"
- }
- set cf_fn "\${cf_prefix}\${name_suffix}.dat"
- if { [ file exists \${cf_fn} ] } {
- exec rm -fv \${cf_fn}
- }
- set cf_fd [ open \${cf_fn} w ]
- set tpro_fd [ open \${tpro_fn} r ]
-
- ## read data from tprofile line by line
- while { [ gets \${tpro_fd} tpro_line ] != -1 } {
- scan \${tpro_line} "%f %f" radius temperature
- #puts "radius: \${radius}, temperature: \${temperature}"
- # set temperature value
- newpar 2 \${temperature}
- # calc flux & tclout
- flux \${e1} \${e2}
- tclout flux 1
- scan \${xspec_tclout} "%f" cf_erg
- #puts "cf: \${cf_erg}"
- puts \${cf_fd} "\${radius} \${cf_erg}"
- }
- close \${tpro_fd}
- close \${cf_fd}
-}
-
-## exit
-tclexit
-_EOF_
-## xcm generation }}}
-
-## invoke xspec to calc
-printf "invoking XSPEC to calculate cooling function profile ...\n"
-xspec - ${XSPEC_CF_XCM} > /dev/null
diff --git a/mass_profile/calc_entropy.py b/mass_profile/calc_entropy.py
deleted file mode 100755
index 736fe7e..0000000
--- a/mass_profile/calc_entropy.py
+++ /dev/null
@@ -1,104 +0,0 @@
-#!/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/calc_lxfx.sh b/mass_profile/calc_lxfx.sh
deleted file mode 100755
index f1954db..0000000
--- a/mass_profile/calc_lxfx.sh
+++ /dev/null
@@ -1,162 +0,0 @@
-#!/bin/sh
-#
-# Wrapper script used to calculate the luminosity (Lx) and flux (Fx) data,
-# which invokes the programming 'calc_lx_beta' (single-beta SBP) or
-# 'calc_lx_dbeta' (double-beta SBP).
-#
-# Output:
-# * lx_result.txt
-# * fx_result.txt
-# * summary_lx.dat
-# * summary_fx.dat
-# * lx_beta_param.txt / lx_dbeta_param.txt
-#
-# Author: Junhua GU
-# Created: 2013-06-24
-#
-# Weitian LI
-# 2016-06-07
-#
-
-if [ $# -eq 2 ] || [ $# -eq 3 ]; then
- :
-else
- echo "usage:"
- echo " `basename $0` <mass.conf> <rout_kpc> [c]"
- echo ""
- echo "arguments:"
- echo " <mass.conf>: config file for mass profile calculation"
- echo " <rout_kpc>: outer/cut radius within which to calculate Lx & Fx"
- echo " e.g., r500, r200 (unit: kpc)"
- echo " [c]: optional; if specified, do not calculate the errors"
- exit 1
-fi
-export PATH="/usr/local/bin:/usr/bin:/bin:$PATH"
-
-mass_cfg="$1"
-rout="$2"
-case "$3" in
- [cC])
- F_C="YES"
- ;;
- *)
- F_C="NO"
- ;;
-esac
-
-base_path=$(dirname $(realpath $0))
-
-## Extract settings/values from the config file
-abund=`grep '^abund' ${mass_cfg} | awk '{ print $2 }'`
-tprofile_data=`grep '^tprofile_data' ${mass_cfg} | awk '{ print $2 }'`
-tprofile_cfg=`grep '^tprofile_cfg' ${mass_cfg} | awk '{ print $2 }'`
-sbp_cfg=`grep '^sbp_cfg' ${mass_cfg} | awk '{ print $2 }'`
-
-sbp_data=`grep '^sbp_data' ${sbp_cfg} | awk '{ print $2 }'`
-tprofile=`grep '^tprofile' ${sbp_cfg} | awk '{ print $2 }'`
-z=`grep '^z' ${sbp_cfg} | awk '{ print $2 }'`
-cm_per_pixel=`grep '^cm_per_pixel' ${sbp_cfg} | awk '{ print $2 }'`
-
-if grep -q '^beta2' $sbp_cfg; then
- MODEL="dbeta"
-else
- MODEL="beta"
-fi
-
-PROG_TPROFILE="fit_wang2012_model"
-tprofile_dump="wang2012_dump.qdp"
-${base_path}/${PROG_TPROFILE} ${tprofile_data} ${tprofile_cfg} \
- ${cm_per_pixel} 2> /dev/null
-mv -fv ${tprofile_dump} ${tprofile}
-
-# energy bands for which the cooling function data will be calculated
-BLIST="blist.txt"
-[ -e "${BLIST}" ] && mv -f ${BLIST} ${BLIST}_bak
-cat > ${BLIST} << _EOF_
-bolo
-0.7 7
-0.1 2.4
-_EOF_
-
-# NOTE:
-# Set 'nh=0' when calculating the cooling function values, and use the
-# value given by 'flux' with unit 'erg/s/cm^2'.
-${base_path}/calc_coolfunc_bands.sh ${tprofile} ${abund} \
- 0 ${z} "cfunc_" ${BLIST}
-
-PROG="calc_lx_${MODEL}"
-LXF_RES="lx_${MODEL}_param.txt"
-${base_path}/${PROG} ${sbp_cfg} ${rout} \
- cfunc_bolo.dat \
- cfunc_0.7-7.dat \
- cfunc_0.1-2.4.dat 2> /dev/null
-LX1=`grep '^Lx1' ${LX_RES} | awk '{ print $2 }'`
-LX2=`grep '^Lx2' ${LX_RES} | awk '{ print $2 }'`
-LX3=`grep '^Lx3' ${LX_RES} | awk '{ print $2 }'`
-FX1=`grep '^Fx1' ${LX_RES} | awk '{ print $2 }'`
-FX2=`grep '^Fx2' ${LX_RES} | awk '{ print $2 }'`
-FX3=`grep '^Fx3' ${LX_RES} | awk '{ print $2 }'`
-
-echo "${LX1} ${LX2} ${LX3}" >summary_lx.dat
-echo "${FX1} ${FX2} ${FX3}" >summary_fx.dat
-
-# save the calculated central values
-mv ${LX_RES} ${LX_RES%.txt}_center.txt
-mv lx_sbp_fit.qdp lx_sbp_fit_center.qdp
-mv lx_rho_fit.dat lx_rho_fit_center.dat
-
-# only calculate the central values
-if [ "${F_C}" = "YES" ]; then
- echo "Calculate the central values only ..."
- ${base_path}/analyze_lxfx.py "Lx" summary_lx.dat lx_result.txt ${BLIST}
- ${base_path}/analyze_lxfx.py "Fx" summary_fx.dat fx_result.txt ${BLIST}
- exit 0
-fi
-
-
-###########################################################
-# Estimate the errors of Lx and Fx by Monte Carlo simulation
-MC_TIMES=100
-for i in `seq 1 ${MC_TIMES}`; do
- ${base_path}/shuffle_profile.py ${tprofile_data} tmp_tprofile.txt
- ${base_path}/shuffle_profile.py ${sbp_data} tmp_sbprofile.txt
-
- # temperature profile
- ${base_path}/${PROG_TPROFILE} tmp_tprofile.txt ${tprofile_cfg} \
- ${cm_per_pixel} 2> /dev/null
- mv -f ${tprofile_dump} ${tprofile}
-
- TMP_SBP_CFG="tmp_sbp.cfg"
- [ -e "${TMP_SBP_CFG}" ] && rm -f ${TMP_SBP_CFG}
- cat ${sbp_cfg} | while read l; do
- if echo "${l}" | grep -q '^sbp_data' >/dev/null; then
- echo "sbp_data tmp_sbprofile.txt" >> ${TMP_SBP_CFG}
- elif echo "${l}" | grep -q '^tprofile' >/dev/null; then
- echo "tprofile ${tprofile}" >> ${TMP_SBP_CFG}
- else
- echo "${l}" >> ${TMP_SBP_CFG}
- fi
- done
-
- echo "### `pwd -P`"
- echo "### ${i} / ${MC_TIMES} ###"
- ${base_path}/calc_coolfunc_bands.sh ${tprofile} ${abund} \
- 0 ${z} "cfunc_" ${BLIST}
- ${base_path}/${PROG} ${TMP_SBP_CFG} ${rout} \
- cfunc_bolo.dat \
- cfunc_0.7-7.dat \
- cfunc_0.1-2.4.dat 2> /dev/null
- LX1=`grep '^Lx1' ${LX_RES} | awk '{ print $2 }'`
- LX2=`grep '^Lx2' ${LX_RES} | awk '{ print $2 }'`
- LX3=`grep '^Lx3' ${LX_RES} | awk '{ print $2 }'`
- FX1=`grep '^Fx1' ${LX_RES} | awk '{ print $2 }'`
- FX2=`grep '^Fx2' ${LX_RES} | awk '{ print $2 }'`
- FX3=`grep '^Fx3' ${LX_RES} | awk '{ print $2 }'`
-
- echo "${LX1} ${LX2} ${LX3}" >>summary_lx.dat
- echo "${FX1} ${FX2} ${FX3}" >>summary_fx.dat
-done # end of 'for'
-
-# analyze Lx & Fx Monte Carlo results
-${base_path}/analyze_lxfx.py "Lx" summary_lx.dat lx_result.txt ${BLIST}
-${base_path}/analyze_lxfx.py "Fx" summary_fx.dat fx_result.txt ${BLIST}
diff --git a/mass_profile/calc_lxfx_wrapper.sh b/mass_profile/calc_lxfx_wrapper.sh
deleted file mode 100755
index 7dea0d0..0000000
--- a/mass_profile/calc_lxfx_wrapper.sh
+++ /dev/null
@@ -1,76 +0,0 @@
-#!/bin/sh
-#
-# Calculate the Lx & Fx data.
-#
-# Based on 'loop_lx.sh', but only process one source
-#
-# Weitian LI
-# 2013-10-30
-#
-
-
-full_path=`readlink -f $0`
-base_dir=`dirname $full_path`
-
-if [ $# -lt 2 ]; then
- printf "usage:\n"
- printf " `basename $0` <global.cfg> [c] < 500 | 200 > ...\n"
- exit 1
-fi
-
-cfg_file="$1"
-pre_results="final_result.txt"
-
-case "$2" in
- [cC]*)
- F_C="YES"
- shift
- ;;
- *)
- F_C="NO"
- ;;
-esac
-
-shift
-echo "delta: $@" # 'printf' not work
-
-if [ "${F_C}" = "YES" ]; then
- printf "MODE: center\n"
-fi
-
-# exit
-
-if [ ! -r "${cfg_file}" ]; then
- printf "ERROR: global cfg not accessible\n"
- exit 11
-elif [ ! -r "${pre_results}" ]; then
- printf "ERROR: previous '${pre_results}' not accessible\n"
- exit 12
-else
- sbp_cfg=`grep '^sbp_cfg' $cfg_file | awk '{ print $2 }'`
- ##
- for delta in $@; do
- if grep -q '^beta2' $sbp_cfg; then
- MODEL="dbeta"
- else
- MODEL="beta"
- fi
- rout=`grep "^r${delta}" ${pre_results} | sed -e 's/=/ /' | awk '{ print $2 }'`
- if [ "${F_C}" = "YES" ]; then
- lx_res="lx_result_${delta}_c.txt"
- fx_res="fx_result_${delta}_c.txt"
- CMD="$base_dir/calc_lxfx.sh $cfg_file $rout c"
- else
- lx_res="lx_result_${delta}.txt"
- fx_res="fx_result_${delta}.txt"
- CMD="$base_dir/calc_lxfx.sh $cfg_file $rout"
- fi
- [ -e "${lx_res}" ] && mv -f ${lx_res} ${lx_res}_bak
- [ -e "${fx_res}" ] && mv -f ${fx_res} ${fx_res}_bak
- ${CMD}
- mv -f lx_result.txt ${lx_res}
- mv -f fx_result.txt ${fx_res}
- done
-fi
-
-exit 0
diff --git a/mass_profile/fg_2500_500.py b/mass_profile/fg_2500_500.py
deleted file mode 100755
index 67a1a11..0000000
--- a/mass_profile/fg_2500_500.py
+++ /dev/null
@@ -1,153 +0,0 @@
-#!/usr/bin/env python
-
-import sys
-import numpy
-import scipy.interpolate
-
-confidence_level=.68
-def read_file(param):
- delta=float(param[0])
-
- file_mass_center=open("mass_int_center.qdp").readlines();
- file_delta_center=open("overdensity_center.qdp").readlines();
-
- center_r=0
- center_m=0
- center_gm=0
- center_gf=0
-
-
- for i in range(0,len(file_mass_center)):
- lm=file_mass_center[i].strip();
- ld=file_delta_center[i].strip();
- r,m=lm.split()
- r,d=ld.split()
- r=float(r)
- d=float(d)
- m=float(m)
- if m<1e11:
- continue
- if d<delta:
- center_r=r
- center_m=m
- for j in open("gas_mass_int_center.qdp"):
- rgm,gm=j.strip().split()
- rgm=float(rgm)
- gm=float(gm)
- if rgm>r:
-
- center_gm=gm
- center_gf=gm/m
- break
- break
- if len(param)>1 and param[1]=='c':
- #print("%s(<r%d)=%E solar mass"%("mass",delta,center_m))
- #print("%s%d=%E kpc"%("r",delta,center_r))
- #print("%s(<r%d)=%E solar mass"%("gas mass",delta,center_gm))
- #print("%s(<r%d)=%E"%("gas fraction",delta,center_gf))
- return center_m,center_r,center_gm,center_gf,None,None,None,None
-
-
-#print(center_gm,center_gf)
- file_mass=open('summary_mass_profile.qdp').readlines()
- file_delta=open('summary_overdensity.qdp').readlines()
- file_gm=open('summary_gas_mass_profile.qdp')
-
-
- flag=True
- rlist=[]
- mlist=[]
- gmlist=[]
- gflist=[]
- old_m=0
- invalid_count=0
- for i in range(0,len(file_mass)):
- lm=file_mass[i].strip()
- ld=file_delta[i].strip()
- if lm[0]=='n':
- flag=True
- old_m=0
- continue
- if not flag:
- continue
- r,m=lm.split()
- m=float(m)
- if m<1e12:
- continue
- if m<old_m:
- invalid_count+=1
- flag=False
- continue
- r,d=ld.split()
- d=float(d)
- r=float(r)
-
- if d<delta:
- #print("%s %e"%(d,m))
- mlist.append(m)
- rlist.append(r)
- flag1=True
- while True:
- lgm=file_gm.readline().strip()
- if lgm[0]=='n':
- break
- rgm,gm=lgm.split()
- rgm=float(rgm)
- gm=float(gm)
- if rgm>r and flag1:
- gmlist.append(gm)
-
- flag1=False
- gflist.append(gm/mlist[-1])
- #print(gm,gflist[-1])
- flag=False
- old_m=m
- print("%d abnormal data dropped"%(invalid_count))
-
-
- return center_m,center_r,center_gm,center_gf,mlist,rlist,gmlist,gflist
-#center_m=numpy.mean(mlist)
-#center_r=numpy.mean(rlist)
-
-if len(sys.argv)>1:
- center_m2500,center_r2500,center_gm2500,center_gf2500,mlist2500,rlist2500,gmlist2500,gflist2500=read_file([2500,sys.argv[1]])
- center_m500,center_r500,center_gm500,center_gf500,mlist500,rlist500,gmlist500,gflist500=read_file([500,sys.argv[1]])
-else:
- center_m2500,center_r2500,center_gm2500,center_gf2500,mlist2500,rlist2500,gmlist2500,gflist2500=read_file([2500])
- center_m500,center_r500,center_gm500,center_gf500,mlist500,rlist500,gmlist500,gflist500=read_file([500])
-
-if mlist2500!=None and len(mlist2500)!=len(mlist500):
- raise Exception("Something wrong, the number of 2500 and 500 data are different")
-
-
-if mlist2500==None:
- print("gas fraction between r2500 and r500 is %E"%((center_gm500-center_gm2500)/(center_m500-center_m2500)))
- sys.exit(0)
-
-gf_2500_500=[]
-
-for i in range(0,len(mlist500)):
- if mlist500[i]-mlist2500[i]<=0:
- continue
- gf_2500_500.append((gmlist500[i]-gmlist2500[i])/(mlist500[i]-mlist2500[i]))
-
-gf_2500_500.sort();
-
-
-center_gf_2500_500=(center_gm500-center_gm2500)/(center_m500-center_m2500)
-gf_idx=-1
-
-for i in range(len(gf_2500_500)-1):
- if (center_gf_2500_500-gf_2500_500[i])*(center_gf_2500_500-gf_2500_500[i+1])<=0:
- gf_idx=i
- break
-if gf_idx==-1:
- raise Exception("Something wrong!")
-
-gflidx=int(gf_idx*(1-confidence_level))
-gfuidx=gf_idx-1+int((len(gf_2500_500)-gf_idx)*confidence_level)
-
-gferr1=gf_2500_500[gflidx]-center_gf_2500_500
-gferr2=gf_2500_500[gfuidx]-center_gf_2500_500
-
-print("gas_fraction between r2500 and r500=\t%e\t %e/+%e (1 sigma)"%(center_gf_2500_500,gferr1,gferr2))
diff --git a/mass_profile/fit_mass.sh b/mass_profile/fit_mass.sh
deleted file mode 100755
index f349f24..0000000
--- a/mass_profile/fit_mass.sh
+++ /dev/null
@@ -1,240 +0,0 @@
-#!/bin/sh
-#
-# Front-end script used to calculate the mass and related values.
-#
-# Output:
-# * final_result.txt / center_only_results.txt
-# * beta_param.txt / dbeta_param_center.txt
-# * gas_mass_int_center.qdp
-# * mass_int_center.qdp
-# * nfw_fit_center.qdp
-# * nfw_param_center.txt
-# * overdensity_center.qdp
-# * rho_fit_center.dat
-# * rho_fit_center.qdp
-# * sbp_fit_center.qdp
-# * entropy_center.qdp
-# * wang2012_param_center.txt
-# * tprofile_dump_center.qdp
-# * tprofile_fit_center.qdp
-# * summary_mass_profile.qdp
-# * summary_overdensity.qdp
-# * summary_gas_mass_profile.qdp
-# * summary_entropy.qdp
-#
-# Junhua Gu
-# Weitian LI
-# 2016-06-07
-#
-
-if [ $# -eq 1 ] || [ $# -eq 2 ]; then
- :
-else
- echo "usage:"
- echo " `basename $0` <mass.conf> [c]"
- echo ""
- echo "arguments:"
- echo " <mass.conf>: config file for mass profile calculation"
- echo " [c]: optional; if specified, do not calculate the errors"
- exit 1
-fi
-
-if ! which xspec > /dev/null 2>&1; then
- printf "*** ERROR: please initialize HEASOFT first\n"
- exit 2
-fi
-
-export PATH="/usr/local/bin:/usr/bin:/bin:$PATH"
-
-base_path=$(dirname $(realpath $0))
-printf "## base_path: \`${base_path}'\n"
-mass_cfg="$1"
-printf "## use configuration file: \`${mass_cfg}'\n"
-case "$2" in
- [cC])
- F_C="YES"
- ;;
- *)
- F_C="NO"
- ;;
-esac
-
-nh=`grep '^nh' ${mass_cfg} | awk '{ print $2 }'`
-abund=`grep '^abund' ${mass_cfg} | awk '{ print $2 }'`
-nfw_rmin_kpc=`grep '^nfw_rmin_kpc' ${mass_cfg} | awk '{ print $2 }'`
-tprofile_data=`grep '^tprofile_data' ${mass_cfg} | awk '{ print $2 }'`
-tprofile_cfg=`grep '^tprofile_cfg' ${mass_cfg} | awk '{ print $2 }'`
-sbp_cfg=`grep '^sbp_cfg' ${mass_cfg} | awk '{ print $2 }'`
-
-# sbp config file
-sbp_data=`grep '^sbp_data' ${sbp_cfg} | awk '{ print $2 }'`
-tprofile=`grep '^tprofile' ${sbp_cfg} | awk '{ print $2 }'`
-cfunc_profile=`grep '^cfunc_profile' ${sbp_cfg} | awk '{ print $2 }'`
-z=`grep '^z' ${sbp_cfg} | awk '{ print $2 }'`
-cm_per_pixel=`cosmo_calc ${z} | grep 'cm/pixel' | awk -F':' '{ print $2 }'`
-sed -i'' "s/^cm_per_pixel.*$/cm_per_pixel ${cm_per_pixel}/" ${sbp_cfg}
-
-if grep -q '^beta2' $sbp_cfg; then
- MODEL="dbeta"
- MODEL_NAME="double-beta"
-else
- MODEL="beta"
- MODEL_NAME="single-beta"
-fi
-
-PROG_TPROFILE="fit_wang2012_model"
-tprofile_dump="wang2012_dump.qdp"
-tprofile_param_center="wang2012_param_center.txt"
-tprofile_fit_center="tprofile_fit_center.qdp"
-tprofile_center="tprofile_dump_center.qdp"
-
-${base_path}/${PROG_TPROFILE} ${tprofile_data} ${tprofile_cfg} \
- ${cm_per_pixel} 2> /dev/null | tee ${tprofile_param_center}
-cp -fv ${tprofile_dump} ${tprofile}
-mv -fv ${tprofile_dump} ${tprofile_center}
-mv -fv fit_result.qdp ${tprofile_fit_center}
-
-${base_path}/calc_coolfunc.sh ${tprofile_center} \
- ${abund} ${nh} ${z} ${cfunc_profile}
-cfunc_profile_center="coolfunc_profile_center.txt"
-cp -f ${cfunc_profile} ${cfunc_profile_center}
-
-PROG_SBPFIT="fit_${MODEL}_sbp"
-RES_SBPFIT="${MODEL}_param.txt"
-RES_SBPFIT_CENTER="${MODEL}_param_center.txt"
-${base_path}/${PROG_SBPFIT} ${sbp_cfg} 2> /dev/null
-mv -fv ${RES_SBPFIT} ${RES_SBPFIT_CENTER}
-cat ${RES_SBPFIT_CENTER}
-mv -fv sbp_fit.qdp sbp_fit_center.qdp
-mv -fv rho_fit.qdp rho_fit_center.qdp
-mv -fv rho_fit.dat rho_fit_center.dat
-mv -fv entropy.qdp entropy_center.qdp
-${base_path}/fit_nfw_mass mass_int.dat ${z} ${nfw_rmin_kpc} 2> /dev/null
-mv -fv nfw_param.txt nfw_param_center.txt
-mv -fv nfw_fit_result.qdp nfw_fit_center.qdp
-mv -fv nfw_dump.qdp mass_int_center.qdp
-mv -fv overdensity.qdp overdensity_center.qdp
-mv -fv gas_mass_int.qdp gas_mass_int_center.qdp
-
-## only calculate central value {{{
-if [ "${F_C}" = "YES" ]; then
- RES_CENTER="center_only_results.txt"
- [ -e "${RES_CENTER}" ] && mv -f ${RES_CENTER} ${RES_CENTER}_bak
- ${base_path}/analyze_mass_profile.py 200 c | tee -a ${RES_CENTER}
- ${base_path}/analyze_mass_profile.py 500 c | tee -a ${RES_CENTER}
- ${base_path}/analyze_mass_profile.py 1500 c | tee -a ${RES_CENTER}
- ${base_path}/analyze_mass_profile.py 2500 c | tee -a ${RES_CENTER}
- ${base_path}/fg_2500_500.py c | tee -a ${RES_CENTER}
- exit 0
-fi
-## central value }}}
-
-## ------------------------------------------------------------------
-
-# clean previous files
-rm -f summary_overdensity.qdp
-rm -f summary_mass_profile.qdp
-rm -f summary_gas_mass_profile.qdp
-rm -f summary_entropy.qdp
-
-# Estimate the errors of Lx and Fx by Monte Carlo simulation
-printf "\n+++++++++++++++++++ Monte Carlo +++++++++++++++++++++\n"
-MC_TIMES=100
-for i in `seq 1 ${MC_TIMES}`; do
- ${base_path}/shuffle_profile.py ${tprofile_data} tmp_tprofile.txt
- ${base_path}/shuffle_profile.py ${sbp_data} tmp_sbprofile.txt
-
- # temperature profile
- ${base_path}/${PROG_TPROFILE} tmp_tprofile.txt ${tprofile_cfg} \
- ${cm_per_pixel} 2> /dev/null
- mv -f ${tprofile_dump} ${tprofile}
-
- TMP_SBP_CFG="tmp_sbp.cfg"
- [ -e "${TMP_SBP_CFG}" ] && rm -f ${TMP_SBP_CFG}
- cat ${sbp_cfg} | while read l; do
- if echo "${l}" | grep -q '^sbp_data' >/dev/null; then
- echo "sbp_data tmp_sbprofile.txt" >> ${TMP_SBP_CFG}
- elif echo "${l}" | grep -q '^tprofile' >/dev/null; then
- echo "tprofile ${tprofile}" >> ${TMP_SBP_CFG}
- else
- echo "${l}" >> ${TMP_SBP_CFG}
- fi
- done
-
- printf "## ${i} / ${MC_TIMES} ##\n"
- printf "## `pwd -P` ##\n"
- ${base_path}/calc_coolfunc.sh ${tprofile} ${abund} ${nh} ${z} ${cfunc_profile}
- ${base_path}/${PROG_SBPFIT} ${TMP_SBP_CFG} 2> /dev/null
- cat ${RES_SBPFIT}
- ${base_path}/fit_nfw_mass mass_int.dat ${z} ${nfw_rmin_kpc} 2> /dev/null
- cat nfw_dump.qdp >> summary_mass_profile.qdp
- echo "no no no" >> summary_mass_profile.qdp
- cat overdensity.qdp >> summary_overdensity.qdp
- echo "no no no" >> summary_overdensity.qdp
- cat gas_mass_int.qdp >> summary_gas_mass_profile.qdp
- echo "no no no" >> summary_gas_mass_profile.qdp
- cat entropy.qdp >> summary_entropy.qdp
- echo "no no no" >> summary_entropy.qdp
-done # end of `for'
-
-# recover the files of original center values
-cp -f ${cfunc_profile_center} ${cfunc_profile}
-cp -f ${tprofile_center} ${tprofile}
-printf "\n+++++++++++++++++ MONTE CARLO END +++++++++++++++++++\n"
-
-## analyze results
-RES_TMP="_tmp_result_mrl.txt"
-RES_FINAL="final_result.txt"
-[ -e "${RES_TMP}" ] && mv -fv ${RES_TMP} ${RES_TMP}_bak
-[ -e "${RES_FINAL}" ] && mv -fv ${RES_FINAL} ${RES_FINAL}_bak
-
-${base_path}/analyze_mass_profile.py 200 | tee -a ${RES_TMP}
-${base_path}/analyze_mass_profile.py 500 | tee -a ${RES_TMP}
-${base_path}/analyze_mass_profile.py 1500 | tee -a ${RES_TMP}
-${base_path}/analyze_mass_profile.py 2500 | tee -a ${RES_TMP}
-
-R200_VAL=`grep '^r200' ${RES_TMP} | awk '{ print $2 }'`
-R500_VAL=`grep '^r500' ${RES_TMP} | awk '{ print $2 }'`
-R1500_VAL=`grep '^r1500' ${RES_TMP} | awk '{ print $2 }'`
-R2500_VAL=`grep '^r2500' ${RES_TMP} | awk '{ print $2 }'`
-
-R200E=`grep '^r200' ${RES_TMP} | tail -n 1 | awk '{ print $2,$3 }'`
-R500E=`grep '^r500' ${RES_TMP} | tail -n 1 | awk '{ print $2,$3 }'`
-R1500E=`grep '^r1500' ${RES_TMP} | tail -n 1 | awk '{ print $2,$3 }'`
-R2500E=`grep '^r2500' ${RES_TMP} | tail -n 1 | awk '{ print $2,$3 }'`
-M200E=`grep '^m200' ${RES_TMP} | tail -n 1 | awk '{ print $2,$3 }'`
-M500E=`grep '^m500' ${RES_TMP} | tail -n 1 | awk '{ print $2,$3 }'`
-M1500E=`grep '^m1500' ${RES_TMP} | tail -n 1 | awk '{ print $2,$3 }'`
-M2500E=`grep '^m2500' ${RES_TMP} | tail -n 1 | awk '{ print $2,$3 }'`
-MG200E=`grep '^gas_m200' ${RES_TMP} | tail -n 1 | awk '{ print $2,$3 }'`
-MG500E=`grep '^gas_m500' ${RES_TMP} | tail -n 1 | awk '{ print $2,$3 }'`
-MG1500E=`grep '^gas_m1500' ${RES_TMP} | tail -n 1 | awk '{ print $2,$3 }'`
-MG2500E=`grep '^gas_m2500' ${RES_TMP} | tail -n 1 | awk '{ print $2,$3 }'`
-FG200E=`grep '^gas_fraction200' ${RES_TMP} | tail -n 1 | awk '{ print $2,$3 }'`
-FG500E=`grep '^gas_fraction500' ${RES_TMP} | tail -n 1 | awk '{ print $2,$3 }'`
-FG1500E=`grep '^gas_fraction1500' ${RES_TMP} | tail -n 1 | awk '{ print $2,$3 }'`
-FG2500E=`grep '^gas_fraction2500' ${RES_TMP} | tail -n 1 | awk '{ print $2,$3 }'`
-
-printf "\n+++++++++++++++ RESULTS (${MODEL_NAME}) +++++++++++++++\n"
-printf "model: ${MODEL_NAME}\n" | tee -a ${RES_FINAL}
-cat ${RES_SBPFIT_CENTER} | tee -a ${RES_FINAL}
-printf "\n" | tee -a ${RES_FINAL}
-printf "r200= ${R200E} kpc\n" | tee -a ${RES_FINAL}
-printf "m200= ${M200E} Msun\n" | tee -a ${RES_FINAL}
-printf "gas_m200= ${MG200E} Msun\n" | tee -a ${RES_FINAL}
-printf "gas_fraction200= ${FG200E} x100%%\n" | tee -a ${RES_FINAL}
-printf "r500= ${R500E} kpc\n" | tee -a ${RES_FINAL}
-printf "m500= ${M500E} Msun\n" | tee -a ${RES_FINAL}
-printf "gas_m500= ${MG500E} Msun\n" | tee -a ${RES_FINAL}
-printf "gas_fraction500= ${FG500E} x100%%\n" | tee -a ${RES_FINAL}
-printf "r1500= ${R1500E} kpc\n" | tee -a ${RES_FINAL}
-printf "m1500= ${M1500E} Msun\n" | tee -a ${RES_FINAL}
-printf "gas_m1500= ${MG1500E} Msun\n" | tee -a ${RES_FINAL}
-printf "gas_fraction1500= ${FG1500E} x100%%\n" | tee -a ${RES_FINAL}
-printf "r2500= ${R2500E} kpc\n" | tee -a ${RES_FINAL}
-printf "m2500= ${M2500E} Msun\n" | tee -a ${RES_FINAL}
-printf "gas_m2500= ${MG2500E} Msun\n" | tee -a ${RES_FINAL}
-printf "gas_fraction2500= ${FG2500E} x100%%\n" | tee -a ${RES_FINAL}
-printf "\n" | tee -a ${RES_FINAL}
-${base_path}/fg_2500_500.py | tee -a ${RES_FINAL}
-printf "\n+++++++++++++++++++++++++++++++++++++++++++++++++++++\n"
diff --git a/mass_profile/fit_sbp.sh b/mass_profile/fit_sbp.sh
deleted file mode 100755
index 5989ff1..0000000
--- a/mass_profile/fit_sbp.sh
+++ /dev/null
@@ -1,61 +0,0 @@
-#!/bin/sh
-#
-# Handy script for SBP fitting.
-# This script wraps the 'fit_beta_sbp' and 'fit_dbeta_sbp',
-# and automatically determine the sbp model according to the config file.
-#
-# Weitian LI
-# 2013-02-20
-#
-# Change logs:
-# 2017-02-07, Weitian LI
-# * Use `sbp_cfg` from command line argument, instead of the one specified
-# in the `mass.conf`
-# * Update the variable names according to the updated config files
-# * Some cleanups
-#
-
-if [ $# -ne 2 ]; then
- echo "usage: $0 <sbp.conf> <mass.conf>"
- exit 1
-fi
-
-sbp_cfg="$1"
-mass_cfg="$2"
-
-if [ "$0" = `basename $0` ]; then
- script_path=`which $0`
- base_path=`dirname ${script_path}`
-else
- base_path=`dirname $0`
-fi
-
-nh=`grep '^nh' ${mass_cfg} | awk '{ print $2 }'`
-abund=`grep '^abund' ${mass_cfg} | awk '{ print $2 }'`
-tprofile_data=`grep '^tprofile_data' ${mass_cfg} | awk '{ print $2 }'`
-tprofile_cfg=`grep '^tprofile_cfg' ${mass_cfg} | awk '{ print $2 }'`
-
-z=`grep '^z' ${sbp_cfg} | awk '{ print $2 }'`
-cm_per_pixel=`cosmo_calc ${z} | grep 'cm/pixel' | awk -F':' '{ print $2 }'`
-sed -i'' "s/^cm_per_pixel.*$/cm_per_pixel ${cm_per_pixel}/" ${sbp_cfg}
-cfunc_profile=`grep '^cfunc_profile' ${sbp_cfg} | awk '{ print $2 }'`
-tprofile=`grep '^tprofile' ${sbp_cfg} | awk '{ print $2 }'`
-
-if grep -q '^beta2' ${sbp_cfg}; then
- MODEL="double-beta"
- PROG=fit_dbeta_sbp
-else
- MODEL="single-beta"
- PROG=fit_beta_sbp
-fi
-
-${base_path}/fit_wang2012_model ${tprofile_data} ${tprofile_cfg} \
- ${cm_per_pixel} 2> /dev/null
-cp wang2012_dump.qdp ${tprofile}
-if [ ! -f ${cfunc_profile} ]; then
- ${base_path}/calc_coolfunc.sh ${tprofile} ${abund} ${nh} ${z} \
- ${cfunc_profile}
-fi
-${base_path}/${PROG} ${sbp_cfg}
-echo "## MODEL: ${MODEL}"
-echo "## z: ${z}"
diff --git a/mass_profile/get_lxfx_data.sh b/mass_profile/get_lxfx_data.sh
deleted file mode 100755
index dbb07dd..0000000
--- a/mass_profile/get_lxfx_data.sh
+++ /dev/null
@@ -1,89 +0,0 @@
-#!/bin/sh
-#
-# collect data for 'loop_lx.sh' & 'calc_lxfx_simple.sh'
-#
-
-if [ $# -lt 2 ]; then
- printf "usage:\n"
- printf " `basename $0` <dir> [c] <delta ...>\n"
- exit 1
-fi
-
-DIR="$1"
-shift
-case "$1" in
- [cC]*)
- F_C="YES"
- shift
- ;;
- *)
- F_C="NO"
- ;;
-esac
-printf "CENTER_MODE: $F_C\n"
-echo "DELTA: $@"
-
-cd $DIR
-pwd -P
-
-INFO=`ls ../*_INFO.json 2> /dev/null`
-if [ ! -z "$INFO" ]; then
- OI=`grep '"Obs\.\ ID' ${INFO} | sed 's/.*"Obs.*":\ //' | sed 's/\ *,$//'`
- NAME=`grep '"Source\ Name' ${INFO} | sed 's/.*"Source.*":\ //' | sed 's/^"//' | sed 's/"\ *,$//'`
- UNAME=`grep '"Unified\ Name' ${INFO} | sed 's/.*"Unified.*":\ //' | sed 's/^"//' | sed 's/"\ *,$//'`
- Z=`grep '"redshift' ${INFO} | sed 's/.*"redshift.*":\ //' | sed 's/\ *,$//'`
-fi
-
-printf "# OI,NAME,UNAME,Z"
-if [ "${F_C}" = "YES" ]; then
- for DELTA in $@; do
- printf ",L${DELTA}(bolo),L${DELTA}(0.7-7),L${DELTA}(0.1-2.4),F${DELTA}(bolo),F${DELTA}(0.7-7),F${DELTA}(0.1-2.4)"
- done
- printf "\n"
-else
- for DELTA in $@; do
- printf ",L${DELTA}(bolo),L${DELTA}ERR(bolo),L${DELTA}(0.7-7),L${DELTA}ERR(0.7-7),L${DELTA}(0.1-2.4),L${DELTA}ERR(0.1-2.4),F${DELTA}(bolo),F${DELTA}ERR(bolo),F${DELTA}(0.7-7),F${DELTA}ERR(0.7-7),F${DELTA}(0.1-2.4),F${DELTA}ERR(0.1-2.4)"
- done
- printf "\n"
-fi
-
-printf "# $OI,$NAME,$UNAME,$Z"
-
-if [ "${F_C}" = "YES" ]; then
- for DELTA in $@; do
- LX_RES="lx_result_${DELTA}_c.txt"
- FX_RES="fx_result_${DELTA}_c.txt"
- if [ -r ${LX_RES} ] && [ -r ${FX_RES} ]; then
- Lbolo=`grep '^Lx(bolo' ${LX_RES} | awk '{ print $2 }'`
- L077=`grep '^Lx(0\.7-7' ${LX_RES} | awk '{ print $2 }'`
- L0124=`grep '^Lx(0\.1-2\.4' ${LX_RES} | awk '{ print $2 }'`
- Fbolo=`grep '^Fx(bolo' ${FX_RES} | awk '{ print $2 }'`
- F077=`grep '^Fx(0\.7-7' ${FX_RES} | awk '{ print $2 }'`
- F0124=`grep '^Fx(0\.1-2\.4' ${FX_RES} | awk '{ print $2 }'`
- printf ",$Lbolo,$L077,$L0124,$Fbolo,$F077,$F0124"
- fi
- done
- printf "\n"
-else
- for DELTA in $@; do
- LX_RES="lx_result_${DELTA}.txt"
- FX_RES="fx_result_${DELTA}.txt"
- if [ -r ${LX_RES} ] && [ -r ${FX_RES} ]; then
- Lbolo=`grep '^Lx(bolo' ${LX_RES} | awk '{ print $2 }'`
- LboloERR=`grep '^Lx(bolo' ${LX_RES} | awk '{ print $4 }'`
- L077=`grep '^Lx(0\.7-7' ${LX_RES} | awk '{ print $2 }'`
- L077ERR=`grep '^Lx(0\.7-7' ${LX_RES} | awk '{ print $4 }'`
- L0124=`grep '^Lx(0\.1-2\.4' ${LX_RES} | awk '{ print $2 }'`
- L0124ERR=`grep '^Lx(0\.1-2\.4' ${LX_RES} | awk '{ print $4 }'`
- Fbolo=`grep '^Fx(bolo' ${FX_RES} | awk '{ print $2 }'`
- FboloERR=`grep '^Fx(bolo' ${FX_RES} | awk '{ print $4 }'`
- F077=`grep '^Fx(0\.7-7' ${FX_RES} | awk '{ print $2 }'`
- F077ERR=`grep '^Fx(0\.7-7' ${FX_RES} | awk '{ print $4 }'`
- F0124=`grep '^Fx(0\.1-2\.4' ${FX_RES} | awk '{ print $2 }'`
- F0124ERR=`grep '^Fx(0\.1-2\.4' ${FX_RES} | awk '{ print $4 }'`
- printf ",$Lbolo,$LboloERR,$L077,$L077ERR,$L0124,$L0124ERR,$Fbolo,$FboloERR,$F077,$F077ERR,$F0124,$F0124ERR"
- fi
- done
- printf "\n"
-fi
-
diff --git a/mass_profile/shuffle_profile.py b/mass_profile/shuffle_profile.py
deleted file mode 100755
index 124f505..0000000
--- a/mass_profile/shuffle_profile.py
+++ /dev/null
@@ -1,37 +0,0 @@
-#!/usr/bin/env python3
-#
-# Shuffle the profile data point values according to their errors.
-#
-# Weitian LI
-# 2017-02-07
-
-import sys
-import numpy as np
-
-
-if len(sys.argv) != 3:
- print("Usage: %s <input_profile> <shuffled_profile>")
- sys.exit(1)
-
-
-# 4-column data: radius, err, temperature/brightness, err
-data = np.loadtxt(sys.argv[1])
-
-x1 = data[:, 2]
-xe = data[:, 3]
-x2 = np.zeros(shape=x1.shape)
-
-for i in range(len(x2)):
- if x1[i] <= 0 or xe[i] <= 0:
- # Skip shuffle
- x2[i] = x1[i]
-
- v = -1.0
- while v <= 0:
- v = np.random.normal(0.0, 1.0) * xe[i] + x1[i]
- x2[i] = v
-
-# Replace original values
-data[:, 2] = x2
-
-np.savetxt(sys.argv[2], data)
Copyright © 2017-2018 Aaron LI. All Rights Reserved.