Add xspec_avg_tz.tcl to help calc average temperature and abundance

Also update xspec_tprofile.tcl a bit

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 }}}