Tag: Matlab

Astronomy · Note

关于冷尘埃SED程序的一些笔记

Posted on by starwave

(这里有一些理论公式推导：http://eagle.lamost.org/?p=37798 ，比此文更加详细，可以参看)

前一阵被X.L. Fan师兄问到关于尘埃能谱 (Dust SED) 做图的问题。然后仔细检查了一遍自己本科论文时候编写的数据，发现有错误，于是花了些时间重新编写了一下。虽然拟合尘埃在亚毫米波段的灰体谱是一件很简单的事情，但在我看来，其中涉及到的单位的问题比较典型。我们平常在做SED的时候，常常会搞混 Flux Density，Intensity，Flux 等概念，他们的单位等等。并且平常老师在课程中讲到各种SED的时候，也没有仔细的讲讲观测中这些东西要注意的地方。这点比较遗憾。

我把程序贴出来跟大家共享一下(点击下载压缩包)。如果你也是在做冷尘埃灰体谱(following the gray-body function or the modified Planck distribution function)辐射的，可以试试这个程序。程序是用 Matlab 写的，一共包含一个主程序 dustsed.m 和两个调用的函数: errorbare.m 和 fitblackbody.m, 其中 errorbare 是由Gerry提供的 ( 详见：Gerry-云居 )。拟合的公式是：

(1)   \begin{equation*} F_{\nu}=\Omega_s Q_{em}B(\nu,T), Q_{em} \propto (\lambda_0 /\lambda)^\beta \end{equation*}

下面测试一下，输入的波长单位是 micron，我取了60, 100, 122, 158, 170, 450, 850 micron 这几个波段的数据，其流量按照Jy的单位分别是65.52, 114.74, 92.0, 91.5, 75.60, 8.832, 1.393 Jy，然后输入测量的误差，调用函数。拟合出单一成份模型中的谱指数、\beta和温度T以及\Omega_{em}（实际上，我在拟合的过程中用的这个参量包含了尘埃的投影面积和距离的信息，但是把含有β的那个因子独立出去了，单独作为一个参量拟合，所以程序里面的 emissivity 不是严格的）。图如下:

给出的结果是T_d=35.581167 K，\beta=1.160865 ，给出的相对误差为0.208285。可以看出，这个Beta明显偏小了。所以应该尝试用更多的温度来拟合，如果用两个温度成份拟合的话会得到更好的结果T_{d}(C)=15.624843 K， T_{d}(W)=32.219577 K， 相对误差为0.205005，这里基于 Dunne L. et al., 2001 MNRAS 的工作，取 \beta=2(有点跑题了……)。

目前 Matlab程序里面拟合使用的函数是 fminsearch ，因为 Q_{em} 很小的缘故，必须用相对误差才能有效的找到比较好的拟合参数（Q_{em}通常在1E-10以下）。下一步要尝试其他的拟合办法，然后在程序中加入 Chi-square 检验的部分。