Move shell/python scripts from 'mass_profile/' to 'bin/'

12 files changed, 1406 insertions, 0 deletions

diff --git a/bin/analyze_lxfx.py b/bin/analyze_lxfx.py
new file mode 100755
index 0000000..62859ff
--- /dev/null
+++ b/bin/analyze_lxfx.py
@@ -0,0 +1,57 @@
+#!/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/bin/analyze_mass_profile.py b/bin/analyze_mass_profile.py
new file mode 100755
index 0000000..3ee128c
--- /dev/null
+++ b/bin/analyze_mass_profile.py
@@ -0,0 +1,195 @@
+#!/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/bin/calc_coolfunc.sh b/bin/calc_coolfunc.sh
new file mode 100755
index 0000000..d1d0b35
--- /dev/null
+++ b/bin/calc_coolfunc.sh
@@ -0,0 +1,113 @@
+#!/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/bin/calc_coolfunc_bands.sh b/bin/calc_coolfunc_bands.sh
new file mode 100755
index 0000000..bebdce2
--- /dev/null
+++ b/bin/calc_coolfunc_bands.sh
@@ -0,0 +1,119 @@
+#!/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/bin/calc_entropy.py b/bin/calc_entropy.py
new file mode 100755
index 0000000..736fe7e
--- /dev/null
+++ b/bin/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/bin/calc_lxfx.sh b/bin/calc_lxfx.sh
new file mode 100755
index 0000000..f1954db
--- /dev/null
+++ b/bin/calc_lxfx.sh
@@ -0,0 +1,162 @@
+#!/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/bin/calc_lxfx_wrapper.sh b/bin/calc_lxfx_wrapper.sh
new file mode 100755
index 0000000..7dea0d0
--- /dev/null
+++ b/bin/calc_lxfx_wrapper.sh
@@ -0,0 +1,76 @@
+#!/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/bin/fg_2500_500.py b/bin/fg_2500_500.py
new file mode 100755
index 0000000..67a1a11
--- /dev/null
+++ b/bin/fg_2500_500.py
@@ -0,0 +1,153 @@
+#!/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/bin/fit_mass.sh b/bin/fit_mass.sh
new file mode 100755
index 0000000..f349f24
--- /dev/null
+++ b/bin/fit_mass.sh
@@ -0,0 +1,240 @@
+#!/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/bin/fit_sbp.sh b/bin/fit_sbp.sh
new file mode 100755
index 0000000..5989ff1
--- /dev/null
+++ b/bin/fit_sbp.sh
@@ -0,0 +1,61 @@
+#!/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/bin/get_lxfx_data.sh b/bin/get_lxfx_data.sh
new file mode 100755
index 0000000..dbb07dd
--- /dev/null
+++ b/bin/get_lxfx_data.sh
@@ -0,0 +1,89 @@
+#!/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/bin/shuffle_profile.py b/bin/shuffle_profile.py
new file mode 100755
index 0000000..124f505
--- /dev/null
+++ b/bin/shuffle_profile.py
@@ -0,0 +1,37 @@
+#!/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)