diff options
Diffstat (limited to 'xspec/xspec_coolfunc_v2.tcl')
| -rw-r--r-- | xspec/xspec_coolfunc_v2.tcl | 246 |
1 files changed, 246 insertions, 0 deletions
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: # |