diff options
| -rw-r--r-- | xspec/xspec_avg_tz.tcl | 109 | ||||
| -rw-r--r-- | xspec/xspec_tprofile.tcl | 14 |
2 files changed, 115 insertions, 8 deletions
diff --git a/xspec/xspec_avg_tz.tcl b/xspec/xspec_avg_tz.tcl new file mode 100644 index 0000000..7720f64 --- /dev/null +++ b/xspec/xspec_avg_tz.tcl @@ -0,0 +1,109 @@ +# Copyright (c) 2017 Weitian LI <liweitianux@live.com> +# MIT license +# +# XSPEC Tcl script to calculate the average temperature and abundance +# by ting the temperatures and abundances of all regions. +# +# NOTE: +# * 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: +# * tz_average.txt : average temperature and abundance (and errors) +# +# Change logs: +# 2017-02-06, Weitian LI +# * Initial version based on `xspec_tprofile.tcl` + + +set NAME "xspec_avg_tz.tcl" +set DATE [ exec date ] +# Set the error flag +set ERR_FLAG "FALSE" + +# Return TCL results for XSPEC commands. +set xs_return_result 1 + +# Keep going until fit converges. +query yes + +# Get the value of the specified parameter +proc get_value {pn} { + tclout param $pn + scan $xspec_tclout "%f" value + return $value +} + +# Calculate the error of the specified parameter +proc get_error {pn} { + global ERR_FLAG + error 1.0 $pn + tclout error $pn + scan $xspec_tclout "%f %f" val_l val_u + set err [ expr {($val_u - $val_l) / 2.0} ] + # error treatment + if {abs($err) < 1.0e-7} { + puts "*** WARNING: invalid error value ***" + set err "NULL" + set ERR_FLAG "TRUE" + } elseif {$err < 0.01} { + set err 0.01 + } + return $err +} + +# Output file +set tz_fn "tz_average.txt" +if {[ file exists $tz_fn ]} { + exec mv -fv $tz_fn ${tz_fn}_bak +} +set tz_fd [ open $tz_fn w ] + +# Header +puts $tz_fd "# Average temperature and abundance" +puts $tz_fd "#" +puts $tz_fd "# Generated by: ${NAME}" +puts $tz_fd "# Created: ${DATE}" +puts $tz_fd "#" +puts $tz_fd "# Item Value Error" + +# Number of data group +tclout datasets +set datasets $xspec_tclout + +# Untie and thaw the temperature and abundance of the first region +untie 5 6 +thaw 5 6 + +# Tie all other temperature and abundance to that of the first region +for {set i 2} {$i <= ${datasets}} {incr i} { + # Parameter number of temperature and abundance + set temp_pn [ expr {8 * $i - 3} ] + set abund_pn [ expr {8 * $i - 2} ] + newpar $temp_pn = 5 + newpar $abund_pn = 6 +} + +# Fit the spectra again +fit + +# Get the values and errors of average temperature and abundance +set temp_val [ get_value 5 ] +set temp_err [ get_error 5 ] +set abund_val [ get_value 6 ] +set abund_err [ get_error 6 ] + +# Output the average temperature and abundance +puts "Average temperature: $temp_val, $temp_err" +puts "Average abundance: $abund_val, $abund_err" +puts $tz_fd "temp $temp_val $temp_err" +puts $tz_fd "abund $abund_val $abund_err" + +close $tz_fd + +# Print a WARNING if there are any errors +if {[ string equal $ERR_FLAG "TRUE" ]} { + puts "*** WARNING: there are invalid error values ***" +} diff --git a/xspec/xspec_tprofile.tcl b/xspec/xspec_tprofile.tcl index e792d3c..ec7263c 100644 --- a/xspec/xspec_tprofile.tcl +++ b/xspec/xspec_tprofile.tcl @@ -2,20 +2,18 @@ ## XSPEC Tcl script to extract the temperature profile from the ## fitted *deprojection spectra analysis* `projct*wabs*apec` results. ## -## Requirements: +## NOTE: ## * Spectra loaded *in order* (in order to correctly calculate radii) ## * Use XSPEC model `projct*wabs*apec` -## -## NOTES: ## * 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) +## * 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 @@ -35,7 +33,7 @@ ## basic variables {{{ -## record process date +set NAME "xspec_avg_tz.tcl" set DATE [ exec date ] set ERR_FLAG "FALSE" ## basic vars }}} |