From cba36462e9e70f45341e432074c726dda5e31092 Mon Sep 17 00:00:00 2001 From: Aaron LI Date: Mon, 20 Feb 2017 12:11:25 +0800 Subject: Move shell/python scripts from 'mass_profile/' to 'bin/' --- mass_profile/analyze_lxfx.py | 57 --------- mass_profile/analyze_mass_profile.py | 195 ---------------------------- mass_profile/calc_coolfunc.sh | 113 ----------------- mass_profile/calc_coolfunc_bands.sh | 119 ----------------- mass_profile/calc_entropy.py | 104 --------------- mass_profile/calc_lxfx.sh | 162 ----------------------- mass_profile/calc_lxfx_wrapper.sh | 76 ----------- mass_profile/fg_2500_500.py | 153 ---------------------- mass_profile/fit_mass.sh | 240 ----------------------------------- mass_profile/fit_sbp.sh | 61 --------- mass_profile/get_lxfx_data.sh | 89 ------------- mass_profile/shuffle_profile.py | 37 ------ 12 files changed, 1406 deletions(-) delete mode 100755 mass_profile/analyze_lxfx.py delete mode 100755 mass_profile/analyze_mass_profile.py delete mode 100755 mass_profile/calc_coolfunc.sh delete mode 100755 mass_profile/calc_coolfunc_bands.sh delete mode 100755 mass_profile/calc_entropy.py delete mode 100755 mass_profile/calc_lxfx.sh delete mode 100755 mass_profile/calc_lxfx_wrapper.sh delete mode 100755 mass_profile/fg_2500_500.py delete mode 100755 mass_profile/fit_mass.sh delete mode 100755 mass_profile/fit_sbp.sh delete mode 100755 mass_profile/get_lxfx_data.sh delete mode 100755 mass_profile/shuffle_profile.py (limited to 'mass_profile') 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 dr: - - center_gm=gm - center_gf=gm/m - break - break - if len(param)>1 and param[1]=='c': - #print("%s(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( -## 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` \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` \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` [c]" - echo "" - echo "arguments:" - echo " : config file for mass profile calculation" - echo " : 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` [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 dr: - - center_gm=gm - center_gf=gm/m - break - break - if len(param)>1 and param[1]=='c': - #print("%s(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` [c]" - echo "" - echo "arguments:" - echo " : 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 " - 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` [c] \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 ") - 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) -- cgit v1.2.2