Sherpa表面亮度拟合

与致力成为X射线数据分析通用工具的HEASOFT不同,Ciao (Chandra Interactive Analysis of Observations) 是专为Chandra卫星数据开发的。它也包含了一个数据分析工具Sherpa。基于现有的脚本语言,从而可以方便地移植扩充。开始有Python 和S-Lang 两个版本,不过现在已经停止了对S-Lang的支持。下面我们介绍如何用它来拟合星系团的面亮度分布。
为了减少计算量,我们先从去除点源的图像中裁出星系团所在区域.使用 dmstat处理前面生成的图像文件.可以得到图像中的极大极小值坐标,如果目标形状不规则,还可以设置 centroid =yes 计算图像重心.

>dmstat image_tot.fits centroid=yes
SRCIMG(x, y)
    min:	0 	      @:	( 2001 2001 )
    max:	40.420722961 	      @:	( 4387 4121 )
cntrd[log] :	( 1105.1675226 1041.8687731 )
cntrd[phys]:	( 4209.3350451 4082.7375461 )
sigma_cntrd:	( 117.40507997 109.33751444 )
   good:	3062500 
   null:	0 

因为这个星系团形状比较规则, 我们这里直接取极值点.

dmcopy "evt2.fits[sky=circle(4387,4121,300)][bin x=::2,y=::2]" circle.fits clobber=yes

下面就可以进入sherpa环境了

sherpa
from sherpa_contrib import *
#导入模块, 如果在安装Ciao之前系统中已经有了python,要注意将Ciao目录中的python函数路径加入;或者直接运行Ciao目录下的Python。
load_image("circle.fits")
set_coord("physical")
set_stat("cash")
set_method("simplex")
notice2d("circle(4387,4121,100)")
image_data()
# 显示图像
set_source(beta2d.src)
# 设定模型,这里我们采用最简单的β模型 $$S(r)=S_{0}[1+(\frac{r}{r_{c}})^{2}]^{-3\beta+1/2}$$
guess(src)
# 赋初值
print src
# 显示模型

---------> print(src)
beta2d.src
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   src.r0       thawed           20         0.02        20000           
   src.xpos     thawed       4389.5       4305.5       4469.5           
   src.ypos     thawed       4123.5       4039.5       4203.5           
   src.ellip    frozen            0            0        0.999           
   src.theta    frozen            0            0      6.28319    radians
   src.ampl     thawed           14        0.014        14000           
   src.alpha    thawed            1          -10           10           


from sherpa_contrib.profiles import *
prof_data()
# 绘制表面亮度曲线
log_scale()
limits(X_AXIS, 0.5, AUTO)
fit()
# 拟合

Dataset               = 1
Method                = neldermead
Statistic             = cash
Initial fit statistic = 11287
Final fit statistic   = 11287 at function evaluation 343
Data points           = 7857
Degrees of freedom    = 7852
Change in statistic   = 0
   src.r0         18.8162     
   src.xpos       4380.16     
   src.ypos       4124.5      
   src.ampl       3.52568     
   src.alpha      0.837185   

参数误差可以通过covariance()计算. 这样就得到了面亮度分布的β模型拟合参数。可以通过图像直观地检查。


image_fit()
# 调用ds9显示当前模型拟合结果,左上角是原始图片,右上角是拟合结果,左下角为残差

prof_fit_resid()
print_window(“win1.png”)
# 显示 profile 拟合结果,可以用ChIPS内部命令保存。
如果要将拟合得到的图像中心位置转换为WCS坐标, 可用下面的命令:

dmcoords evt2.fits asolfile=asol1.fits opt=sky x=4380.16 y=4124.5 verbose=1

由这个结果我们可以看出,β模型可以较好地拟合星系团外围亮度分布,但无法解释星系团中心的亮度超出。因为这个模型基于等温假设,亮度谱指数中$$\beta = \mu m_{p} \sigma_{r}^{2}/kT_{g}$$(具体含义点此查看)。而在实际星系团内总是存在温度梯度。因此双β模型 set_source(beta2d.src+beta2d.core) 通常能给出更好的拟合结果. 详见 Radial and elliptical profiles of Image Data

订阅评论
提醒
2 评论
最旧
最新 最多投票
内联反馈
查看所有评论
hbcai
2010 年 8 月 20 日 17:42

好。有时间我将TIMING也整理出来。