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

>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

1	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

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

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

1	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