##
## XSPEC Tcl script to extract the temperature profile from the
## fitted *deprojection spectra analysis* `projct*wabs*apec` results.
##
## NOTE:
## * Spectra loaded *in order* (in order to correctly calculate radii)
## * Use XSPEC model `projct*wabs*apec`
## * 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 parameters)
##
## Output:
## * xspec_saveall.xcm : save all for the current state
## * fitted_params.log : log the fitted parameters
## * tprofile.qdp : temp profile QDP file for graphing
## * tprofile.txt : temp profile data file for later fit
##
## Weitian LI <liweitianux@live.com>
## Created: 2012-08-12
##
## Change logs:
## 2017-02-05, Weitian LI
##   * Drop the cooling function calculation
##   * Rename to 'xspec_tprofile.tcl'
##   * Cleanups
## v2.0, 2012-08-14, Weitian LI
##   * add `output a QDP file for temperature profile'
##   * cancel `calc errors of abundance'
##   * improve error treatment (add warning)
## v1.1, 2012-08-12
##   * improve description
##   * modify error treatment


## basic variables {{{
set NAME "xspec_avg_tz.tcl"
set DATE [ exec date ]
set ERR_FLAG "FALSE"
## basic vars }}}

## xspec settings {{{
# 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 "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 "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 "tprofile.qdp"
set tpro_fn "tprofile.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
}
set tpro_qdp_fd [ open $tpro_qdp_fn w ]
set tpro_fd [ open $tpro_fn w ]
## output files }}}

## datasets, number of data group
tclout datasets
set datasets $xspec_tclout

## 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 fitted results
    # 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
    # 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"
}

## close file
close $tpro_qdp_fd
close $tpro_fd

## check `ERR_FLAG', print WARNING
if {[ string equal $ERR_FLAG "TRUE" ]} {
    puts "*** WARNING: there are WRONG error values"
}

# vim: set ts=8 sw=4 tw=0 fenc= ft=tcl: #