Download
1 / 24
240 likes | 367 Views
NGC3603 中的尘埃颗粒. 新模型. 内容. 拟和 : 检验 模型 方法 结果 : 温度 尘埃密度 Gas-to-dust ratio TODO. 拟和检验. 卡方检验 : 测量次数 >50, 样本 ? 个体 ?(BingJiang) Cash 检验 (cstat): 中心极限定理之一 , 大数近似 , 复合泊松 -> 高斯分布 . 测量次数不确定 (MSX,helpdesk) Ref: Ciao(Yangchen),xspec(Yangchen),cash et al Xspec manual 及附录 ,
E N D
NGC3603中的尘埃颗粒 新模型
内容 拟和: • 检验 • 模型 • 方法 结果: • 温度 • 尘埃密度 • Gas-to-dust ratio TODO
拟和检验 • 卡方检验: 测量次数>50,样本?个体?(BingJiang) • Cash检验(cstat): • 中心极限定理之一,大数近似,复合泊松->高斯分布. • 测量次数不确定(MSX,helpdesk) • Ref: • Ciao(Yangchen),xspec(Yangchen),cash et al • Xspec manual及附录, • Xspec作者主页http://lheawww.gsfc.nasa.gov/~kaa/head_stats/title.html • Least-square. Goodness. – 没有好的goodness评价方案,拟和时舍弃sigma数据(errorbar).
拟和模型 • 4PAH+blackbody • 考虑到21um附近没有PAH辐射,令c=0,得到3PAH+blackbody 结果都不好,怀疑mathematica的结果,局部最小?且有警告. 换其他软件拟和?编程找到合适的a,b,c值?
拟和软件 • Mathematica: NonlinearFit – ?分段函数 • Matlab: cftool, Optimization Toolbox - ? .. • GSL:http://www.gnu.org/software/gsl • 1stOpt : “世界上最好的拟和软件” 算法最多(UGO)!支持全局优化!不受初值影响!可设定参数范围!支持代码分段!
GSL代码 • #include "t.c" int main(void) { size_t np = 3; double par[2] = {1.39, 0.8}; const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex; gsl_multimin_fminimizer *s = NULL; gsl_vector *ss, *x; gsl_multimin_function minex_func; size_t iter = 0, i; int status; double size; /* Initial vertex size vector */ ss = gsl_vector_alloc (np); /* Set all step sizes to 1 */ gsl_vector_set_all (ss, 1.0); /* Starting point */ x = gsl_vector_alloc (np); gsl_vector_set (x, 0, 1.0);/* U */ gsl_vector_set (x, 1, 100);/* v */ gsl_vector_set (x, 2, 20);/* t */ /* Initialize method and iterate */ minex_func.f = &my_f; minex_func.n = np; minex_func.params = (void *)∥ s = gsl_multimin_fminimizer_alloc (T, np); gsl_multimin_fminimizer_set (s, &minex_func, x, ss); do { iter++; status = gsl_multimin_fminimizer_iterate(s); if (status) break; size = gsl_multimin_fminimizer_size (s); status = gsl_multimin_test_size (size, 1e-10); if (status == GSL_SUCCESS) { printf ("converged to minimum at\n"); } printf ("%5d ", iter); for (i = 0; i < np; i++) { printf ("%G ", gsl_vector_get (s->x, i)); } printf ("f() = %G size = %G\n", s->fval, size); } while (status == GSL_CONTINUE && iter < 100); //while (status == GSL_CONTINUE ); gsl_vector_free(x); gsl_vector_free(ss); gsl_multimin_fminimizer_free (s); return status; }
GSL的结果 • LM方法非线性拟和:T=11K,目标函数值偏大,明显不是全局最小! • 最小值方法: 直接退出 . 无法编程寻找a,b 于是决定用土方法: 1stOpt穷举a,b对.
1stOpt: 限制nd,寻找合理的a,b对 a 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 -ring 0.8 0.746 0.699 0.660 0.631 0.605 0.584 0.566 0.552 0.539 0.529 -A 0.915 0.845 0.790 0.745 0.708 0.678 0.653 0.633 0.616 0.6 0.588 -B 0.746 0.692 0.649 0.615 0.587 0.563 0.544 0.528 0.515 0.504 0.495 -C 0.563 0.527 0.499 0.476 0.458 0.443 0.432 0.422 0.414 0.408 0.403 -D 0.725 0.671 0.627 0.592 0.564 0.541 0.521 0.505 0.492 0.480 0.471 可以看到各个区域的b值单调减小，在a[1,2]内,b永不相交!
新的参数方案 • 理由: 所有区域使用一致的PAH比例是不可能的. • 做法: exp[500/30*0.6]=exp[10]>>1,考虑到黑体谱只对E波段有影响,所以对于ACD,PAH比例从数据直接估算.
估算方案 • 对原始数据做除法,求出对应区域每个点的流量比再对区域平均. Bin->div->平均. Msx+ds9 Region文件的bin: ciao格式(BingJiang). • 直接对最终数据(Fnu)做除法.
结果分析 • C区域nd偏大3个量级,T偏小(20K) • BD区域偏小. • A~10ring. 仍然无法得到一致的尘埃密度 只选A,ring区.
1stOpt代码 Title "ring"; Parameter u = 1 [0.1,10],v=100000000[10000,100000000000] ,t = 30[20,300]; Constant h=479.92156638356494; Constant a=4.86/3.72,b=3.20/4.86; DataSet; x,y,s= 3.62319 3.72061 0.771501 2.47321 4.85386 1.44736 2.04778 3.20737 0.763918 1.40581 6.71629 2.05759 EndDataSet; ConstStr FA=u+v*(x1^4)/(exp(h/t*x1)-1); ConstStr FC=a*u+v*(x2^4)/(exp(h/t*x2)-1); ConstStr FD=b*a*u+v*(x3^4)/(exp(h/t*x3)-1); ConstStr FE=v*(x4^4)/(exp(h/t*x4)-1); MinFunction (FA-y1)^2+(FC-y2)^2+(FD-y3)^2+(FE-y4)^2;
尘埃密度和尘气比 • Q取ss433文:1.45*10^-16nu,得到尘埃密度:10 ^-9 • (model2使用draine 1993之SiC数据) Vring=6.401E+08 Va=1.10011E+07 Lunit=0.000203622 Unit=8.44253E-12 (kpc)(kpc)(kpc) region ring: model1: nd=1.49419E-10 g/d=4427.27 model2: nd=3.46971E-10 g/d=1906.55 region A: model1: nd=2.99752E-09 g/d=220.689 model2: nd=6.93232E-09 g/d=95.4254
Q,Lambd-2 • Draine 1993: :Q(SiC)~Q(ss433)~0.1Q(Si) • 选sic还是si?why?or 其他颗粒? • Q的选取.只照顾到E就ok?(>21,q~lambd-2,draine 1984)
TODO • Draine 1984 和 Takeuchi文比较,后者做作lambda-2 fit(Yangchen). • 寻找NGC3603温度依据30kk? (22kk now). • 1000颗恒星的依据 • Q值.
