added xspec scripts: xspec_bkgcorr_v2.tcl, xspec_coolfunc_v2.tcl

2 files changed, 616 insertions, 0 deletions

diff --git a/xspec/xspec_bkgcorr_v2.tcl b/xspec/xspec_bkgcorr_v2.tcl
new file mode 100644
index 0000000..61a09f1
--- /dev/null
+++ b/xspec/xspec_bkgcorr_v2.tcl
@@ -0,0 +1,370 @@
+#####################################################################
+## XSPEC Tcl script
+##
+## Task:
+## analysis chandra acis background components
+## output model fitting results, then generate
+## fake spectra to correct the blanksky background
+##
+## LIweitiaNux <liweitianux@gmail.com>
+## August 1, 2012
+##
+## NOTES: needs XSPEC v12.x
+##
+## ChangeLogs:
+## v2, 2012/08/14, LIweitiaNux
+##   skip calc `error' of params
+##   change variable names
+##   improve the final comparison graph
+##   replace `trim' with `regexp'
+#####################################################################
+
+## about {{{
+set NAME "xspec_bkgcorr_v2.tcl"
+set VERSION "v2, 2012-08-14"
+## about }}}
+
+## set basic variables {{{
+## tclout: create tcl var from current state
+set tcl_date [ exec date ]
+tclout filename 1
+set spec_name [ exec basename $xspec_tclout ]
+# `trim' is not reliable: `lbkg_grp.pi' --> `lbkg_gr'
+# set spec_rootname "[ string trimright $spec_name {.pi} ]"
+regexp -- {(.+)\.(pi)} $spec_name match spec_rootname fext
+# get `rootname' from `spec_name'
+regexp -- {(.+)_(bin|grp)([0-9]*)\.(pi|fits)} $spec_name match rootname grptype grpval fext
+puts "## rootname: $rootname"
+puts "## grptype: $grptype"
+puts "## grpval: $grpval"
+puts "## file extension: $fext"
+# backgrnd name
+tclout backgrnd 1
+set back_name $xspec_tclout
+# set back_rootname "[ string trimright $back_name {.pi} ]"
+regexp -- {(.+)\.(pi)} $back_name match back_rootname fext
+# regexp -- {(.+)\.(pi)} $spec_name match spec_rootname fext
+# `specscal': scale factor to create fake spectrum,
+# use a large factor (i.e. 10) to get a good statistic for final spectrum
+set specscal 10.0
+## basic variables }}}
+
+## save current xspec fitting results {{{
+set xspec_outroot "xspec_${rootname}_fit"
+# `save all', to save a xcm script
+if {[ file exists ${xspec_outroot}.xcm ]} {
+    exec mv -fv ${xspec_outroot}.xcm ${xspec_outroot}.xcm_bak
+}
+save all "${xspec_outroot}.xcm"
+# writefits: (tcl scripts, v12.x)
+#if {[ file exists ${xspec_outroot}.fits ]} {
+#    exec mv -fv ${xspec_outroot}.fits ${xspec_outroot}.fits_bak
+#}
+#writefits "${xspec_outroot}.fits"
+## save & writefits }}}
+
+## set output file {{{
+set out_fn "bkgcorr_${rootname}.log"
+if {[ file exists ${out_fn} ]} {
+    exec mv -fv ${out_fn} ${out_fn}_bak
+}
+set out_fd [ open ${out_fn} w ]
+## output file }}}
+
+## header, output basic process info {{{
+puts $out_fd "# process by script: ${NAME}"
+puts $out_fd "# version: ${VERSION}"
+puts $out_fd "#"
+puts $out_fd "# process date: $tcl_date"
+puts $out_fd "#"
+puts $out_fd ""
+## header }}}
+
+## save basic info about current fitting {{{
+# files in use
+puts $out_fd "data $spec_name"
+tclout response 1
+set rmf $xspec_tclout
+puts $out_fd "response $rmf"
+tclout arf 1
+set arf $xspec_tclout
+puts $out_fd "arf $arf"
+puts $out_fd "backgrnd $back_name"
+puts $out_fd ""
+# exposure time, backscale
+tclout expos 1 s
+scan $xspec_tclout "%f" expos_src
+tclout expos 1 b
+scan $xspec_tclout "%f" expos_bkg
+puts $out_fd "# src/bkg exptime: $expos_src/$expos_bkg"
+tclout backscal 1 s
+scan $xspec_tclout "%f" backscal_src
+tclout backscal 1 b
+scan $xspec_tclout "%f" backscal_bkg
+puts $out_fd "# src/bkg backscal: $backscal_src/$backscal_bkg"
+puts $out_fd ""
+# model, and basic fitting results
+tclout model
+puts $out_fd "model $xspec_tclout"
+tclout weight
+puts $out_fd "# weight: $xspec_tclout"
+tclout stat
+scan $xspec_tclout "%f" chisq
+tclout dof
+scan $xspec_tclout "%d" dof
+set rechisq_cmd { format "%.5f" [ expr { $chisq / $dof } ] }
+set rechisq [ eval $rechisq_cmd ]
+puts $out_fd "# chisq/dof = $chisq/$dof = $rechisq"
+tclout noticed 1
+puts $out_fd "# noticed channel: $xspec_tclout"
+tclout noticed energy 1
+puts $out_fd "# noticed energy (keV): $xspec_tclout"
+puts $out_fd ""
+## basic fitting info }}}
+
+## model fitting results {{{
+puts $out_fd "# fitting results of each components"
+puts $out_fd "# errors are emitted for simplicity"
+puts $out_fd "#"
+puts $out_fd "# name   value   sigma_err"
+puts $out_fd ""
+# apec_1, Galactic X-ray background, soft-soft
+tclout param 1
+scan $xspec_tclout "%f" temp_a1
+puts $out_fd "temp_a1: $temp_a1"
+tclout param 2
+scan $xspec_tclout "%f" abund_a1
+puts $out_fd "abund_a1: $abund_a1"
+tclout param 3
+scan $xspec_tclout "%f" redshift_a1
+puts $out_fd "redshift_a1: $redshift_a1"
+tclout param 4
+scan $xspec_tclout "%f" norm_a1
+tclout sigma 4
+set norm_sigma_a1 $xspec_tclout
+puts $out_fd "norm_a1: $norm_a1  $norm_sigma_a1"
+puts $out_fd ""
+# apec_2, Galactic X-ray Background, soft
+tclout param 5
+scan $xspec_tclout "%f" temp_a2
+puts $out_fd "temp_a2: $temp_a2"
+tclout param 6
+scan $xspec_tclout "%f" abund_a2
+puts $out_fd "abund_a2: $abund_a2"
+tclout param 7
+scan $xspec_tclout "%f" redshift_a2
+puts $out_fd "redshift_a2: $redshift_a2"
+tclout param 8
+scan $xspec_tclout "%f" norm_a2
+tclout sigma 8
+set norm_sigma_a2 $xspec_tclout
+puts $out_fd "norm_a2: $norm_a2  $norm_sigma_a2"
+puts $out_fd ""
+# wabs
+tclout param 9
+scan $xspec_tclout "%f" wabs
+puts $out_fd "wabs: $wabs"
+puts $out_fd ""
+# powerlaw, Cosmological X-ray Background, hard
+tclout param 10
+scan $xspec_tclout "%f" gamma_po
+puts $out_fd "gamma_po: $gamma_po"
+tclout param 11
+scan $xspec_tclout "%f" norm_po
+tclout sigma 11
+set norm_sigma_po $xspec_tclout
+puts $out_fd "norm_po: $norm_po  $norm_sigma_po"
+puts $out_fd ""
+# apec_3, gas emission from object source
+tclout param 12
+scan $xspec_tclout "%f" temp_gas
+tclout sigma 12
+set temp_sigma_gas $xspec_tclout
+puts $out_fd "temp_gas: $temp_gas $temp_sigma_gas"
+tclout param 13
+scan $xspec_tclout "%f" abund_gas
+tclout sigma 13
+set abund_sigma_gas $xspec_tclout
+puts $out_fd "abund_gas: $abund_gas $abund_sigma_gas"
+tclout param 14
+scan $xspec_tclout "%f" redshift_gas
+puts $out_fd "redshift_gas: $redshift_gas"
+tclout param 15
+scan $xspec_tclout "%f" norm_gas
+tclout sigma 15
+set norm_sigma_gas $xspec_tclout
+puts $out_fd "norm_gas: $norm_gas  $norm_sigma_gas"
+puts $out_fd ""
+## fittins results }}}
+# end output results
+close $out_fd
+
+## prepare data for fake spectrum {{{
+# see `specscal' in top `basic variables' section
+# apec_1
+if {$norm_a1 < 0} {
+    set pm_a1 "-"
+} else {
+    set pm_a1 "+"
+}
+set pars_a1 "$temp_a1 & $abund_a1 & $redshift_a1 & [ expr abs($norm_a1)*$specscal ]"
+set model_a1 "apec"
+# apec_2
+if {$norm_a2 < 0} {
+    set pm_a2 "-"
+} else {
+    set pm_a2 "+"
+}
+set pars_a2 "$temp_a2 & $abund_a2 & $redshift_a2 & [ expr abs($norm_a2)*$specscal ]"
+set model_a2 "apec"
+# powerlaw
+if {$norm_po < 0} {
+    set pm_po "-"
+} else {
+    set pm_po "+"
+}
+set pars_po "$wabs & $gamma_po & [ expr abs($norm_po)*$specscal ]"
+set model_po "wabs*powerlaw"
+# gas emssion
+# norm of `gas component' cannot be negative
+set pars_gas "$wabs & $temp_gas & $abund_gas & $redshift_gas & [ expr abs($norm_gas)*$specscal ]"
+set model_gas "wabs*apec"
+## prepare data }}}
+
+## functions to load model, fake spectrum {{{
+proc tcl_model {model_str pars args} {
+    model none
+    model $model_str & $pars & /*
+}
+proc tcl_fakeit {rmf arf fakedata exptime args} {
+    data none
+    fakeit none & $rmf & $arf & y & & $fakedata & $exptime & /*
+}
+## functions }}}
+
+## set fake spectrum and check previous files {{{
+set fake_a1 "fake_${rootname}_a1.pi"
+set fake_a2 "fake_${rootname}_a2.pi"
+set fake_po "fake_${rootname}_po.pi"
+set fake_gas "fake_${rootname}_gas.pi"
+if {[ file exists $fake_a1 ]} {
+    exec mv -fv $fake_a1 ${fake_a1}_bak
+}
+if {[ file exists $fake_a2 ]} {
+    exec mv -fv $fake_a2 ${fake_a2}_bak
+}
+if {[ file exists $fake_po ]} {
+    exec mv -fv $fake_po ${fake_po}_bak
+}
+if {[ file exists $fake_gas ]} {
+    exec mv -fv $fake_gas ${fake_gas}_bak
+}
+## fake spectrum }}}
+
+## generate fake spectrum {{{
+# blanksky, apec_1
+tcl_model $model_a1 $pars_a1
+tcl_fakeit $rmf $arf $fake_a1 $expos_bkg
+# blanksky, apec_2
+tcl_model $model_a2 $pars_a2
+tcl_fakeit $rmf $arf $fake_a2 $expos_bkg
+# blanksky, wabs*powerlaw
+tcl_model $model_po $pars_po
+tcl_fakeit $rmf $arf $fake_po $expos_bkg
+# src gas, wabs*apec
+tcl_model $model_gas $pars_gas
+tcl_fakeit $rmf $arf $fake_gas $expos_src
+## fake spectrum }}}
+
+## background correction {{{
+## blanksky {{{
+set bbkg_expr "${specscal}*${back_name} ${pm_a1}${fake_a1} ${pm_a2}${fake_a2} ${pm_po}${fake_po}"
+set bbkg_outf "bkgcorr_blanksky_${rootname}.pi"
+set bbkg_back [ expr { $backscal_bkg * $specscal } ]
+set bbkg_cmm1 "corrected background spectrum, blanksky based"
+set bbkg_cmm2 "norm*${specscal}, backscal*${specscal}, properr=no"
+set bbkg_mathpha "mathpha expr=\"${bbkg_expr}\" units=C outfil=${bbkg_outf} exposure=${expos_bkg} areascal=% ncomments=2 comment1=\"${bbkg_cmm1}\" comment2=\"${bbkg_cmm2}\" properr=no backscal=${bbkg_back}"
+## blanksky }}}
+
+## localbkg {{{
+set lbkg_expr "${specscal}*${spec_name} -${fake_gas}"
+set lbkg_outf "bkgcorr_localbkg_${rootname}.pi"
+set lbkg_back [ expr { $backscal_src * $specscal } ]
+set lbkg_cmm1 "corrected background spectrum, local background based"
+set lbkg_cmm2 "norm*${specscal}, backscal*${specscal}, properr=no"
+set lbkg_mathpha "mathpha expr=\"${lbkg_expr}\" units=C outfil=${lbkg_outf} exposure=${expos_src} areascal=% ncomments=2 comment1=\"${lbkg_cmm1}\" comment2=\"${lbkg_cmm2}\" properr=no backscal=${lbkg_back}"
+## local }}}
+
+## XXX: cannot figure out how to `eval' these cmd in Tcl. !!!
+## generate a shell script to exec mathpha {{{
+if {[ file exists ${bbkg_outf} ]} {
+    exec mv -fv ${bbkg_outf} ${bbkg_outf}_bak
+}
+if {[ file exists ${lbkg_outf} ]} {
+    exec mv -fv ${lbkg_outf} ${lbkg_outf}_bak
+}
+set mathpha_fn "_mathpha.sh"
+if {[ file exists ${mathpha_fn} ]} {
+    exec rm -fv ${mathpha_fn}
+}
+set mathpha_fd [ open ${mathpha_fn} w ]
+puts $mathpha_fd $bbkg_mathpha
+puts $mathpha_fd $lbkg_mathpha
+close $mathpha_fd
+# exec
+exec sh ${mathpha_fn}
+#exec rm -rf ${mathpha_sh}
+## mathpha }}}
+## background correction }}}
+
+## load the corrected spectra, save a picture {{{
+set bkgcorr_img "bkgcorr_${rootname}.ps"
+if {[ file exists ${bkgcorr_img} ]} {
+    exec mv -fv ${bkgcorr_img} ${bkgcorr_img}_bak
+}
+
+model none
+data none
+
+# corrected spectrum (blanksky based)
+data 1:1 ${bbkg_outf}
+resp 1:1 ${rmf}
+arf 1:1 ${arf}
+# corrected spectrum (local bkg based)
+data 2:2 ${lbkg_outf}
+resp 1:2 ${rmf}
+arf 1:2 ${arf}
+# original blanksky spectrum
+data 3:3 ${back_name}
+resp 1:3 ${rmf}
+arf 1:3 ${arf}
+# original local spectrum
+data 4:4 ${spec_name}
+resp 1:4 ${rmf}
+arf 1:4 ${arf}
+
+ignore bad
+ignore **:0.0-0.3,11.0-**
+
+setplot energy
+
+# rebin for plot
+setplot rebin 20 20 1
+setplot rebin 20 20 2
+setplot rebin 20 20 3
+setplot rebin 10 20 4
+
+# label
+setplot command LABEL T "Corrected & Original Background Spectrum Comparison"
+setplot command LABEL F "corr_blank(BLACK), corr_local(RED), orig_blank(GREEN), orig_local(BLUE)"
+
+setplot device ${bkgcorr_img}/cps
+plot ldata
+setplot device /xw
+plot
+## bkgcorr picture }}}
+
+## end analysis
+# tclexit
+
+# vim: set ts=8 sw=4 tw=0 fenc=utf-8 ft=tcl: #
diff --git a/xspec/xspec_coolfunc_v2.tcl b/xspec/xspec_coolfunc_v2.tcl
new file mode 100644
index 0000000..0bfcbf5
--- /dev/null
+++ b/xspec/xspec_coolfunc_v2.tcl
@@ -0,0 +1,246 @@
+#####################################################################
+## XSPEC Tcl script
+##
+## Task:
+## To generate `Temperature profile'
+## and `cooling function data file'
+## according to `deprojection spectra analysis' results
+##
+## Requirements:
+## deprojection spectra analysis,
+## spectra loaded *in order* (in order to correctly calc *radius*)
+## use XSPEC model `projct*wabs*apec' (AS IT IS)
+## because this script cannot handle other *parameter numbers*
+##
+## NOTES:
+## *solar abundance std* assumed to be `GRSA' (1998)
+## if error < 0.01, then assume error = 0.01
+## if fabs(error) < 1e-7, then set error as `NULL' in output
+## (this may be caused by frozen or tied params)
+##
+## Output:
+## tcl_xspec_saveall.xcm    (save all for the current state)
+## tcl_fitted_params.log    (log the fitted parameters)
+## tcl_temp_profile.qdp     (temp profile QDP file for graphing)
+## tcl_temp_profile.txt     (temp profile data file for later fit)
+## tcl_coolfunc_dat.txt     (coolfunc data file for later fit)
+##
+## LIweitiaNux <liweitianux@gmail.com>
+## August 12, 2012
+##
+## v1.1, 2012/08/12
+##   improve description
+##   modify error treatment
+## v2.0, 2012/08/14, LIweitiaNux
+##   add `output a QDP file for temperature profile'
+##   cancel `calc errors of abundance'
+##   improve error treatment (add warning)
+#####################################################################
+
+## about {{{
+set NAME "xspec_coolfunc_v2.tcl"
+set VERSION "v2, 2012-08-14"
+## about }}}
+
+## ask `normalization' used to calculate cooling function
+puts -nonewline "LY> Input the NORMALIZATION used to calc cooling function: "
+set norm [ gets stdin ]
+
+## basic variables {{{
+## record process date
+set DATE [ exec date ]
+set ERR_FLAG "FALSE"
+## basic vars }}}
+
+## xspec settings {{{
+# solar abundance std
+# use `grsa' (1998)
+set abund_std "grsa"
+
+# Return TCL results for XSPEC commands.
+set xs_return_result 1
+
+# Keep going until fit converges.
+query yes
+## xspec settings }}}
+
+## xspec, save current results {{{
+set fn_save "tcl_xspec_saveall.xcm"
+if {[ file exists $fn_save ]} {
+    exec mv -fv $fn_save ${fn_save}_bak
+}
+save all $fn_save
+## xspec save all }}}
+
+## xspec log, log current fitted params {{{
+set fn_log "tcl_fitted_params.log"
+if {[ file exists $fn_log ]} {
+    exec mv -fv $fn_log ${fn_log}_bak
+}
+log $fn_log
+show param
+log none
+## xspec log }}}
+
+## set output file {{{
+set tpro_qdp_fn "tcl_temp_profile.qdp"
+set tpro_fn "tcl_temp_profile.txt"
+set cool_fn "tcl_coolfunc_dat.txt"
+# tmp file used to calc `cooling function'
+set tmpxspec_fn "_tcl_xspec_tmp.xcm"
+set tmpres_fn "_tcl_coolres_tmp.txt"
+# check file status
+if {[ file exists $tpro_qdp_fn ]} {
+    exec mv -fv $tpro_qdp_fn ${tpro_qdp_fn}_bak
+}
+if {[ file exists $tpro_fn ]} {
+    exec mv -fv $tpro_fn ${tpro_fn}_bak
+}
+if {[ file exists $cool_fn ]} {
+    exec mv -fv $cool_fn ${cool_fn}_bak
+}
+set tpro_qdp_fd [ open $tpro_qdp_fn w ]
+set tpro_fd [ open $tpro_fn w ]
+set cool_fd [ open $cool_fn w ]
+## output files }}}
+
+## datasets, number of data group
+tclout datasets
+set datasets $xspec_tclout
+
+## get `nH' and `redshift' at first
+tclout param 4
+scan $xspec_tclout "%f" nh
+tclout param 7
+scan $xspec_tclout "%f" redshift
+
+## QDP file header {{{
+puts $tpro_qdp_fd "! generated by: ${NAME}"
+puts $tpro_qdp_fd "! created date: ${DATE}"
+puts $tpro_qdp_fd "SKIP SINGLE"
+puts $tpro_qdp_fd "READ SERR 1 2"
+puts $tpro_qdp_fd ""
+puts $tpro_qdp_fd "LABEL T Temperature Profile"
+puts $tpro_qdp_fd "LABEL F temperature error at 68% confidence level"
+puts $tpro_qdp_fd "LABEL X Radius (pixel)"
+puts $tpro_qdp_fd "LABEL Y Temperature (keV)"
+puts $tpro_qdp_fd ""
+puts $tpro_qdp_fd "! radius(pixel) radius_err temperature(keV) temp_err(60%)"
+## QDP }}}
+
+## output temperature profile and abundance
+for {set i 1} {$i <= ${datasets}} {incr i} {
+    # get the MID-radius of current spectra
+    tclout xflt $i
+    scan $xspec_tclout "%f %f" holder r_out
+    if {$i == 1} {
+        # the innermost region
+        set r_in 0.0
+    } else {
+        set j [ expr {$i - 1} ]
+        tclout xflt $j
+        scan $xspec_tclout "%f %f" holder r_in
+    }
+    set r [ expr {($r_in + $r_out) / 2.0} ]
+    set r_err [ expr {($r_out - $r_in) / 2.0} ]
+    puts "radius: $r, radius_err: $r_err ($r_in, $r_out)"
+    # output values: temperature and abundance
+    # determine the param number of temperature and abundance
+    # temperature
+    set temp_pn [ expr {8 * $i - 3} ]
+    tclout param $temp_pn
+    scan $xspec_tclout "%f" temp_val
+    error 1.0 $temp_pn
+    tclout error $temp_pn
+    scan $xspec_tclout "%f %f" temp_val_l temp_val_u
+    set temp_err [ expr {($temp_val_u - $temp_val_l) / 2.0} ]
+    # error treatment
+    if {abs($temp_err) < 1.0e-7} {
+        puts "*** WARNING: WRONG error values"
+        set temp_err "NULL"
+        set ERR_FLAG "TRUE"
+    } elseif {$temp_err < 0.01} {
+        set temp_err 0.01
+    }
+    puts "temperature (p${temp_pn}): $temp_val, $temp_err"
+    # abundance
+    set abund_pn [ expr {8 * $i - 2} ]
+    tclout param $abund_pn
+    scan $xspec_tclout "%f" abund_val
+    #error 1.0 $abund_pn
+    #tclout error $abund_pn
+    #scan $xspec_tclout "%f %f" abund_val_l abund_val_u
+    #set abund_err [ expr {($abund_val_u - $abund_val_l) / 2.0} ]
+    #if {abs($abund_err) < 1.0e-7} {
+    #    puts "*** WARNING: WRONG error values"
+    #    set abund_err "NULL"
+    #    set ERR_FLAG "TRUE"
+    #} elseif {$abund_err < 0.01} {
+    #    set abund_err 0.01
+    #}
+    #puts "abundance (p${abund_pn}): $abund_val $abund_err"
+    # output `temp_profile' related results
+    puts $tpro_qdp_fd "$r    $r_err    $temp_val    $temp_err"
+    puts $tpro_fd "$r    $r_err    $temp_val    $temp_err"
+
+    ## generate a xspec script to calc cooling function {{{
+    # remove previous files
+    if {[ file exists $tmpxspec_fn ]} {
+        exec rm -fv $tmpxspec_fn
+    }
+    if {[ file exists $tmpres_fn ]} {
+        exec rm -fv $tmpres_fn
+    }
+    # open script file
+    set tmpxspec_fd [ open $tmpxspec_fn w ]
+    puts $tmpxspec_fd "set tmpres_fn \"$tmpres_fn\""
+    puts $tmpxspec_fd "set tmpres_fd \[ open \$tmpres_fn w \]"
+    puts $tmpxspec_fd "query yes"
+    puts $tmpxspec_fd "statistic chi"
+    puts $tmpxspec_fd "method leven 10 0.01"
+    puts $tmpxspec_fd "abund ${abund_std}"
+    puts $tmpxspec_fd "xsect bcmc"
+    puts $tmpxspec_fd "cosmo 70 0 0.73"
+    puts $tmpxspec_fd "xset delta 0.01"
+    puts $tmpxspec_fd "systematic 0"
+    puts $tmpxspec_fd "dummyrsp 0.3 11.0 1024"
+    puts $tmpxspec_fd "model wabs*apec & ${nh} & ${temp_val} & ${abund_val} & ${redshift} & $norm & /*"
+    puts $tmpxspec_fd "flux 0.7 7.0"
+    puts $tmpxspec_fd "tclout flux 1"
+    puts $tmpxspec_fd "scan \$xspec_tclout \"%f %f %f %f\" holder holder holder res_flux"
+    puts $tmpxspec_fd "puts \$res_flux"
+    puts $tmpxspec_fd "puts \$tmpres_fd \$res_flux"
+    puts $tmpxspec_fd "close \$tmpres_fd"
+    puts $tmpxspec_fd "tclexit"
+    close $tmpxspec_fd
+    ## xspec script }}}
+    ## calc cooling function and output results {{{
+    # exec XSPEC to calculate
+    exec xspec - $tmpxspec_fn
+    # read results from files
+    set tmpres_fd [ open $tmpres_fn r ]
+    set tmpline [ gets $tmpres_fd ]
+    scan $tmpline "%f" coolf_val
+    puts "cooling function: $coolf_val"
+    close $tmpres_fd
+    puts $cool_fd "$r    $coolf_val"
+    ## cooling function }}}
+}
+
+## remove tmp files
+if {[ file exists $tmpres_fn ]} {
+    exec rm -fv $tmpres_fn
+}
+
+## close file
+close $tpro_qdp_fd
+close $tpro_fd
+close $cool_fd
+
+## check `ERR_FLAG', print WARNING
+if {[ string equal $ERR_FLAG "TRUE" ]} {
+    puts "*** WARNING: there are WRONG error values"
+}
+
+# EOF
+# vim: set ts=8 sw=4 tw=0 fenc= ft=tcl: #