aboutsummaryrefslogtreecommitdiffstats
path: root/bin/calc_lxfx.sh
blob: f1954dbc37b0f6f0165c8a3a18567004ba681505 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
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}