1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
|
##
## 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: #
|
