path: root/xspec/xspec_coolfunc_v2.tcl

#####################################################################
## 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: #