#!/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'.
##
## Weitian LI
## Created: 2012-08-17
## Updated: 2016-06-08
##

## cmdline arguments {{{
if [ $# -eq 5 ]; then
    :
elif [ $# -eq 6 ]; then
    COOLFUNC_BOLO="$6"
    [ -e "${COOLFUNC_BOLO}" ] && rm -f ${COOLFUNC_BOLO}
else
    printf "usage:\n"
    printf "    `basename $0` <tprofile> <avg_abund> <nH> <redshift> <coolfunc_outfile> [coolfunc_bolo]\n"
    exit 1
fi
TPROFILE=$1
ABUNDANCE=$2
N_H=$3
REDSHIFT=$4
COOLFUNC_DAT=$5
COOLFUNC_DAT_RATIO="flux_cnt_ratio.txt"
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}
[ -e "${COOLFUNC_DAT_RATIO}" ] && rm -f ${COOLFUNC_DAT_RATIO}
## arguments }}}

## specify variable name outside while loop
## otherwise the inside vars invisible
XSPEC_CF_XCM="_coolfunc_calc.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}"
set cff_fn "${COOLFUNC_DAT_RATIO}"
if { [ file exists \${cf_fn} ] } {
    exec rm -fv \${cf_fn}
}
if { [ file exists \${cff_fn} ] } {
    exec rm -fv \${cff_fn}
}

## open files
set tpro_fd [ open \${tpro_fn} r ]
set cf_fd [ open \${cf_fn} w ]
set cff_fd [ open \${cff_fn} w ]

_EOF_

if  [ ! -z "${COOLFUNC_BOLO}" ]; then
    cat >> ${XSPEC_CF_XCM} << _EOF_
# coolfunc bolometric
set cfbolo_fn "${COOLFUNC_BOLO}"
if { [ file exists \${cfbolo_fn} ] } {
    exec rm -fv \${cfbolo_fn}
}
set cfbolo_fd [ open \${cfbolo_fn} w ]

_EOF_
fi

cat >> ${XSPEC_CF_XCM} << _EOF_
## 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}"
    flux 0.01 100.0
    tclout flux 1
    scan \${xspec_tclout} "%f" cff_erg
    puts \${cff_fd} "\${radius}   [expr \${cff_erg}/\${cf_photon}]"
_EOF_
if  [ ! -z "${COOLFUNC_BOLO}" ]; then
    cat >> ${XSPEC_CF_XCM} << _EOF_
    # coolfunc bolometric
    puts \${cfbolo_fd} "\${radius}    \${cff_erg}"
_EOF_
fi
cat >> ${XSPEC_CF_XCM} << _EOF_
}

## close opened files
close \${tpro_fd}
close \${cf_fd}
_EOF_

if  [ ! -z "${COOLFUNC_BOLO}" ]; then
    cat >> ${XSPEC_CF_XCM} << _EOF_
# coolfunc bolometric
close \${cfbolo_fd}

_EOF_
fi

cat >> ${XSPEC_CF_XCM} << _EOF_
## exit
tclexit
_EOF_
## xcm generation }}}

## invoke xspec to calc
printf "invoking XSPEC to calculate cooling function profile ...\n"
xspec - ${XSPEC_CF_XCM} > /dev/null