光谱立方体建模#
用例： 从数据立方体中提取感兴趣的空间-光谱特征并测量其属性。根据模型参数图创建谱线和多环芳烃（PAH）图。
数据： 使用CUBISM生成的SL1和SL2光谱立方体，来源于斯皮策（Spitzer）红外光谱（IRS）对附近星系梅西耶58中心的光谱映射观测。
工具： astropy, specutils, dust-extinction, matplotlib。
跨仪器： NIRSpec, MIRI。
文档： 本笔记本是STScI更大后处理数据分析工具生态系统的一部分。
引言#
JWST光谱立方体的分析需要提取感兴趣的空间-光谱特征并测量其属性。在这里，我们演示了应用于单个spaxel的1-D光谱建模，立方体中的900个spaxels，以及在空间区域上求和的光谱。分析大型JWST光谱数据立方体可能会计算上非常昂贵，这取决于建模的spaxel数量和自由模型参数的数量。通过关注感兴趣的空间子区域或最小化模型复杂性，可以使这项任务变得可管理。立方体中的子区域通过位置和谱线比率进行选择。然后对这些区域中的spaxels进行求和，并使用连续体、发射线、多环芳烃（PAH）和硅酸盐尘埃特征的组合进行建模。最佳拟合模型光谱使用lmfit包进行拟合，采用Levenberg-Marquardt最小二乘法。
导入#
time 用于计时
multiprocessing.Pool 用于对立方体像素进行并行处理
pickle 用于在多进程池中进行序列化和反序列化输入输出
numpy 用于数组处理和数学运算
matplotlib.pyplot 用于绘制图像和光谱
astropy.io 用于读取和写入FITS立方体和图像
astropy.visualization 用于构建RGB图像
lmfit Levenberg-Marquart拟合包，用于模型拟合
import time # 导入时间模块
import multiprocess as mp # 导入多进程模块
from multiprocess import Pool # 从多进程模块导入进程池
import pickle # 导入序列化模块
import numpy as np # 导入NumPy库
from numpy import exp, loadtxt, pi, sqrt, arange, interp # 从NumPy导入常用函数
import matplotlib.pyplot as plt # 导入Matplotlib库用于绘图
from matplotlib.colors import LogNorm # 从Matplotlib导入对数归一化
from astropy.io import fits # 从Astropy导入FITS文件处理模块
from astropy.nddata import StdDevUncertainty # 从Astropy导入标准偏差不确定性类
import astropy.units as u # 导入Astropy单位模块
from astropy.visualization import make_lupton_rgb # 从Astropy导入Lupton RGB图像生成函数
from astropy.visualization import quantity_support # 从Astropy导入量支持模块
from astropy.modeling import models # 从Astropy导入模型模块
from specutils import Spectrum1D # 从specutils导入一维光谱类
import lmfit # 导入lmfit模块用于模型拟合
from lmfit import Model, Minimizer, Parameters # 从lmfit导入模型、最小化器和参数类
from lmfit.model import save_modelresult # 从lmfit导入保存模型结果的函数
函数的定义#
滚动#
限制单元输出中的滚动。
%%javascript
// 设置IPython输出区域的自动滚动阈值为70
IPython.OutputArea.auto_scroll_threshold=70;
模型函数#
用于光谱拟合的多种发射线、连续谱和消光模型。
def Modelcero(x):
"""简单模型0，用于初始化一些模型"""
return 0*x # 返回一个与x相同形状的零数组
# powerlaw, zblackbody, gaussian 和 drude 函数包装了现有的 astropy 模型。
# lmfit 需要将波长 (x) 作为函数参数传递。
def powerlaw(x, c1, c2):
"""幂律函数 (c1*nu^c2)，在1微米处归一化，采用正指数约定"""
x0 = 1.0 # 在1微米处归一化
c2 = -c2 # 反转指数
return models.PowerLaw1D.evaluate(x, c1, x0, c2) # 计算并返回幂律模型的值
def zblackbody(x, OnePlusZ, T, scale):
"""黑体函数"""
xs = x * u.um / OnePlusZ # 将波长转换为红移后的波长
bb = models.BlackBody(T * u.K, scale) # 创建黑体模型
return bb(xs) # 返回黑体辐射值
def gaussian(x, amplitude, xcen, std):
"""一维高斯函数"""
return models.Gaussian1D.evaluate(x, amplitude, xcen, std) # 计算并返回高斯模型的值
def drude(x, amplitude, peakx, frac_FWHM):
"""尘埃发射模型"""
FWHM = peakx * frac_FWHM # 计算全宽半高
return models.Drude1D.evaluate(x, amplitude, peakx, FWHM) # 计算并返回Drude模型的值
# pahdust 和 sidust 函数在 astropy 中不存在
def pahdust(x, OnePlusZ, amplitude_76, amplitude_113):
"""PAH尘埃发射模型"""
# 离子化PAH特征
PAH_peakx = [5.27, 5.70, 6.22, 6.69, 7.42, 7.60, 7.85, 8.33, 8.61] # PAH峰值波长
PAH_frac_FWHM = [0.034, 0.035, 0.030, 0.07, 0.126, 0.044, 0.053, 0.050, 0.039] # PAH全宽半高比例
PAH_rel_amplitude = [0.0, 0.0, 0.8, 0.0, 0.0, 1.0, 0.6, 0.0, 0.5] # PAH相对振幅
PAH_amplitude = [] # 初始化PAH振幅列表
# 计算PAH振幅
for ampl in PAH_rel_amplitude:
PAH_amplitude.append(amplitude_76 * (ampl / PAH_rel_amplitude[5])) # 归一化振幅
# 中性PAH特征
PAH_peakx += [10.68, 11.23, 11.33] # 添加中性PAH峰值
PAH_peakx += [11.99, 12.62, 12.69, 13.48, 14.04, 14.19, 15.90, 16.45, 17.04, 17.375, 17.87, 18.92] # 添加更多峰值
PAH_frac_FWHM += [0.020, 0.012, 0.022] # 添加中性PAH全宽半高比例
PAH_frac_FWHM += [0.045, 0.042, 0.013, 0.040, 0.016, 0.025, 0.020, 0.014, 0.065, 0.012, 0.016, 0.019] # 添加更多比例
PAH_rel_amplitude = [0.0, 0.6, 1.0] # 中性PAH相对振幅
PAH_rel_amplitude += [0.2, 0.3, 0.1, 0.0, 0.1, 0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] # 添加更多相对振幅
# 计算中性PAH振幅
for ampl in PAH_rel_amplitude:
PAH_amplitude.append(amplitude_113 * (ampl / PAH_rel_amplitude[2])) # 归一化振幅
pahflux = x - x # 初始化PAH通量
# 计算PAH通量
for peakx, frac_FWHM, ampl in zip(PAH_peakx, PAH_frac_FWHM, PAH_amplitude):
pahflux = pahflux + drude(x, ampl, peakx * OnePlusZ, frac_FWHM) # 累加每个Drude模型的贡献
return pahflux # 返回PAH通量（振幅单位）
def sidust(x, T, amplitude):
"""硅酸盐尘埃发射模型"""
# 自定义消光曲线，由三种成分的加权和构成：
# 两个Drude函数和一个指数1.7的幂律。
d1 = drude(x, 0.80, 10.0, 0.25) # 第一个Drude模型
d2 = drude(x, 0.25, 11.1, 0.25) # 第二个Drude模型
# d3 = drude(x, 0.40, 17.0, 0.40) # 超出波长范围。
ext = d1 + d2 # 计算总的消光
# 形成修改后的硅酸盐和幂律的线性组合。
beta = 0.1 # 线性组合的权重
ext = (1. - beta) * ext + beta * (9.7 / x) ** 1.7 # 计算加权消光
# 计算硅酸盐发射
si_em = amplitude * 1.0E-6 * 3.97289E13 / x ** 3 / (exp(1.4387752E4 / x / T) - 1.) * ext
return si_em # 返回硅酸盐发射值
# 启用表格消光函数。现有的dust_extinction选项在中红外波长下不足
def extinction_a(x, tau, a):
"""来自表格数据的消光"""
ext = np.interp(x, wave, a) # 插值计算消光
return exp(-tau * ext) # 返回消光值
# 从Chiar & Tielens (2005)读取消光数据
# 注意，这个数据是归一化到K带的消光（在2.14微米处为1.004）
# agal是银河中心的消光，alocal是局部ISM消光
Boxdata = "https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/cube_fitting/"
wave, agal, alocal = np.loadtxt(Boxdata + 'chiar+tielens_2005.dat', skiprows=14, usecols=(0, 1, 2), unpack=True) # 读取数据
# 与dust_extinction进行比较
# P92是唯一一个扩展到中红外范围的模型，但它比较粗糙
# 归一化到V带的消光
from dust_extinction.parameter_averages import F99 # 导入F99模型
from dust_extinction.shapes import P92 # 导入P92模型
a92 = P92() # 创建P92模型实例
# 消光比较
tau = 1. # 消光的光学深度
xa = arange(5., 15., 0.01) # 定义波长范围
extin1 = extinction_a(xa, tau, alocal) # 计算局部ISM消光
extin2 = extinction_a(xa, tau, agal) # 计算银河中心消光
lam = xa * u.um # 将波长转换为单位微米
scale = 21.9 # 经验缩放因子，用于调整P2曲线以匹配Chiar & Tielens局部ISM曲线
aa92 = scale * a92(lam) # 计算缩放后的P92消光
extin3 = exp(-tau * aa92) # 计算缩放后的消光值
with quantity_support(): # 启用量支持
plt.plot(lam, extin1, label='CT+05 Local ISM') # 绘制局部ISM消光
plt.plot(lam, extin2, label='CT+05 Gal. Cen.') # 绘制银河中心消光
plt.plot(lam, extin3, label='P92 scaled') # 绘制缩放后的P92消光
plt.legend() # 显示图例
plt.show() # 显示图形
pl1 = powerlaw(xa, 1.0, -1.0) # 计算幂律模型
模型参数#
声明模型函数的参数。通过在lmfit的setpar命令中设置vary=’False’来冻结参数。
目前我们只是在拟合发射特征的幅度。
具有更多可变参数的拟合需要更长的时间。
def setpar(pars, name, value, vary, minus):
"""设置任意参数"""
pars[name].set(value, vary=vary, min=minus) # 设置参数的值、是否可变及最小值
def gaussian_defpar(name, amplitude, mean, stddev, pars):
"""设置高斯模型的参数"""
std = stddev # 标准差
name = str(name) # 将名称转换为字符串
xcen = mean * OneZ # 中心位置乘以红移因子
setpar(pars, name + '_std', std, False, None) # 设置标准差参数
# setpar(pars, name + '_std', std, True, 0) # 可变标准差参数（注释掉）
setpar(pars, name + '_xcen', xcen, False, None) # 设置中心位置参数
setpar(pars, name + '_amplitude', amplitude, True, 0) # 设置幅度参数
return # 返回时不需要任何值
def drude_defpar(name, amplitude, peakx, frac_FWHM, pars):
"""设置Drude模型的参数"""
peakx = peakx * OneZ # 峰值乘以红移因子
setpar(pars, name + '_frac_FWHM', frac_FWHM, False, None) # 设置FWHM参数
# setpar(pars, name + '_frac_FWHM', frac_FWHM, True, 0) # 可变FWHM参数（注释掉）
setpar(pars, name + '_peakx', peakx, False, None) # 设置峰值参数
setpar(pars, name + '_amplitude', amplitude, True, 0) # 设置幅度参数
return # 返回时不需要任何值
def pahdust_defpar(name, OnePlusZ, amplitude_76, amplitude_113, pars):
"""设置PDR尘埃模型的参数"""
setpar(pars, name + '_amplitude_76', amplitude_76, True, 0) # 设置76微米幅度参数
setpar(pars, name + '_amplitude_113', amplitude_113, True, 0) # 设置113微米幅度参数
setpar(pars, name + '_OnePlusZ', OnePlusZ, False, None) # 设置红移参数
return # 返回时不需要任何值
def sidust_defpar(name, T, amplitude, pars):
"""设置SIDUST模型的参数"""
setpar(pars, name + '_amplitude', amplitude, True, 0) # 设置幅度参数
setpar(pars, name + '_T', T, False, 0) # 设置温度参数
return # 返回时不需要任何值
def gaussian_extractpars(prefix, result):
"""提取高斯模型的参数"""
amp = result.params[prefix + 'amplitude'].value # 提取幅度
cen = result.params[prefix + 'xcen'].value # 提取中心位置
std = result.params[prefix + 'std'].value # 提取标准差
return amp, cen, std # 返回幅度、中心位置和标准差
def drude_extractpars(prefix, result):
"""提取Drude模型的参数"""
ampl = result.params[prefix + 'amplitude'].value # 提取幅度
peak = result.params[prefix + 'peakx'].value # 提取峰值
frac = result.params[prefix + 'frac_FWHM'].value # 提取FWHM
return ampl, peak, frac # 返回幅度、峰值和FWHM
def pahdust_extractpars(prefix, result):
"""提取PDR尘埃模型的参数"""
ampl_76 = result.params[prefix + 'amplitude_76'].value # 提取76微米幅度
ampl_113 = result.params[prefix + 'amplitude_113'].value # 提取113微米幅度
return ampl_76, ampl_113 # 返回76和113微米的幅度
def sidust_extractpars(prefix, result):
"""提取SIDUST模型的参数"""
T = result.params[prefix + 'T'].value # 提取温度
ampl = result.params[prefix + 'amplitude'].value # 提取幅度
return T, ampl # 返回温度和幅度
模型参数的线通量#
用于根据模型参数计算发射特征的通量和宽度的辅助函数。
# 注意：这些是线通量的解析估计，
# 不可在astropy建模或specutils中获得
def gaussianline_flux(amp, std, cen):
"""计算高斯线的积分通量"""
# 单位：amp: Jy， cen: 微米
c = 29979.2458 * 10**10 # 光速，单位为微米/s
gaussianfactor = sqrt(2 * pi) # 高斯因子
return amp * gaussianfactor * c * std / (cen**2) * 10**(-23) # erg s^{-1} cm^{-2}
def drudeline_flux(amp, frac, peak):
"""计算Drude线的积分通量"""
# 单位：amp: Jy， peak: 微米
c = 29979.2458 * 10**10 # 光速，单位为微米/s
drudefactor = pi # Drude因子
return (amp * frac * drudefactor * c / (2 * peak)) * 10**(-23) # erg s^{-1} cm^{-2}
提取并绘制一维光谱及模型组件#
def extract_spec(a, b):
"""从a,b坐标提取光谱和误差"""
spec_pix = data_cube[:, a, b] # 从数据立方体中提取指定坐标的光谱数据
err_pix = error_cube[:, a, b] # 从误差立方体中提取指定坐标的误差数据
return spec_pix, err_pix # 返回光谱数据和误差数据
def model_comps(model_result):
"""评估模型组件"""
comps = model_result.eval_components() # 评估模型的各个组件
plcomp = [] # 用于存储功率律组件
h2comp = [] # 用于存储氢分子组件
pahcomp = [] # 用于存储PAH组件
sidustcomp = [] # 用于存储硅尘埃组件
# print(comps.keys()) # 打印组件的键（可选）
for key in comps.keys(): # 遍历所有组件的键
keyl = key.lower() # 将键转换为小写
if keyl[0:3] == 'pwl':
plcomp.append(key) # 如果是功率律组件，添加到plcomp列表
if keyl[0] == 'h':
h2comp.append(key) # 如果是氢分子组件，添加到h2comp列表
if keyl[0:3] == 'pah':
pahcomp.append(key) # 如果是PAH组件，添加到pahcomp列表
if keyl[0:3] == 'sid':
sidustcomp.append(key) # 如果是硅尘埃组件，添加到sidustcomp列表
# 初始化模型
plaw_model = x - x # 初始化功率律模型
h2_model = x - x # 初始化氢分子模型
ion_model = x - x # 初始化离子模型
pah_model = x - x # 初始化PAH模型
sidust_model = x - x # 初始化硅尘埃模型
# 累加各个组件的贡献
for comp in plcomp:
plaw_model += comps[comp] # 累加功率律组件
for comp in h2comp:
h2_model += comps[comp] # 累加氢分子组件
for comp in pahcomp:
pah_model += comps[comp] # 累加PAH组件
for comp in sidustcomp:
sidust_model += comps[comp] # 累加硅尘埃组件
atomiclines = ['ArII', 'SIV', 'NeII'] # 定义原子线组件
for comp in atomiclines:
ion_model += comps[comp + '_'] # 累加原子线组件
return ([plaw_model, h2_model, ion_model, pah_model, sidust_model], ['PL', 'H2', 'Ions', 'PAHs', 'Si Dust']) # 返回模型和标签
def plot_fit(x, spec, specerr, model_result):
"""绘制光谱、模型组件和残差"""
# 模型结果
fit_model = model_result.best_fit # 获取最佳拟合模型
fit_residual = spec - fit_model # 计算残差
# 评估模型组件
mod_comps, mod_labels = model_comps(model_result) # 获取模型组件和标签
# 创建光谱1D对象
spec1 = Spectrum1D(spectral_axis=x * u.um, flux=spec * u.Jy, uncertainty=StdDevUncertainty(specerr)) # 创建光谱对象
with quantity_support():
f, ax = plt.subplots() # 创建子图
ax.plot(spec1.spectral_axis, spec1.flux, label='Data', color='k') # 绘制数据光谱
ax.grid(linestyle=':') # 设置网格样式
plt.show() # 显示图形
# 绘制结果
fig1 = plt.figure(figsize=(15, 6), dpi=150) # 创建图形
frame1 = fig1.add_axes((.1, .3, .8, .6)) # 添加第一个坐标轴
frame1.set_xticklabels([]) # 隐藏x轴标签
plt.ylabel(r"$Flux\ (Jy)$") # 设置y轴标签
plt.grid(linestyle=':') # 设置网格样式
plt.errorbar(x, spec, yerr=specerr, label='Data', color='k') # 绘制带误差条的数据
plt.plot(x, fit_model, label='Model', c='red') # 绘制模型
mod_colors = ['brown', 'g', 'b', 'magenta', 'orange'] # 定义模型组件颜色
for mcomp, mlabel, mcolr in zip(mod_comps, mod_labels, mod_colors):
plt.plot(x, mcomp, label=mlabel, c=mcolr) # 绘制模型组件
plt.legend() # 显示图例
frame2 = fig1.add_axes((.1, .1, .8, .2)) # 添加第二个坐标轴
plt.plot(x, 0. * x, '-', c='r') # 绘制零线
plt.xlabel(r"$Wavelength\ (\mu m)$") # 设置x轴标签
plt.grid(linestyle=':') # 设置网格样式
plt.errorbar(x, fit_residual, yerr=specerr, c='k') # 绘制残差带误差条
plt.show() # 显示图形
return # 返回
输入数据#
由于目前缺乏JWST飞行数据，我们利用使用CUBISM生成的SL1和SL2光谱立方体，这些立方体来自对附近星系梅西耶58中心的斯皮策IRS光谱映射观测。该星系具有一个发射强烈射电辐射的塞弗特核，正在通过冲击激发强烈的H2发射。此外，PAH发射特征揭示了显著的恒星形成活动。
# Spitzer IRS (CUBISM) cube loader does not exist in specutils
# 目标名称
targname='M58'
# 红移
z=0.005060
OneZ=1.+z # 计算红移后的单位
# 下载并打开数据立方体及其不确定性
BoxPath="https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/cube_fitting/"
cubeSL1 = fits.open(BoxPath+targname+'_SL1_cube.fits') # 打开SL1立方体
errorsSL1 = fits.open(BoxPath+targname+'_SL1_cube_unc.fits') # 打开SL1不确定性立方体
cubeSL2 = fits.open(BoxPath+targname+'_SL2_cube.fits') # 打开SL2立方体
errorsSL2 = fits.open(BoxPath+targname+'_SL2_cube_unc.fits') # 打开SL2不确定性立方体
# 立方体信息和头文件
cubeSL1.info() # 输出SL1立方体信息
cubeSL2.info() # 输出SL2立方体信息
hdr1 = cubeSL1[0].header # 获取SL1头文件
er_hdr1=errorsSL1[0].header # 获取SL1不确定性头文件
hdr2 = cubeSL2[0].header # 获取SL2头文件
er_hdr2=errorsSL2[0].header # 获取SL2不确定性头文件
# Flux 数据
data_cube1 = cubeSL1[0].data # 获取SL1数据
error_cube1 = errorsSL1[0].data # 获取SL1不确定性数据
data_cube2 = cubeSL2[0].data # 获取SL2数据
error_cube2 = errorsSL2[0].data # 获取SL2不确定性数据
# 波长数据
xwave1 = cubeSL1[1].data # 获取SL1波长数据
xwave1 = xwave1.field(0)[0] # 提取SL1波长
xwave2 = cubeSL2[1].data # 获取SL2波长数据
xwave2 = xwave2.field(0)[0] # 提取SL2波长
x=[] # 初始化波长列表
for line in xwave2: # 遍历SL2波长
x.append(float(line)) # 将SL2波长添加到列表
for line in xwave1: # 遍历SL1波长
x.append(float(line)) # 将SL1波长添加到列表
x=np.array(x) # 转换为NumPy数组
# 将单位从MJy/sr转换为Jy/pix
correct_unit1=abs(hdr1['CDELT1'])*abs(hdr1['CDELT2'])*0.0003046118*(10**6) # 计算SL1单位转换因子
data_cube1=data_cube1*correct_unit1 # 应用单位转换
error_cube1=error_cube1*correct_unit1 # 应用单位转换
correct_unit2=abs(hdr2['CDELT1'])*abs(hdr2['CDELT2'])*0.0003046118*(10**6) # 计算SL2单位转换因子
data_cube2=data_cube2*correct_unit2 # 应用单位转换
error_cube2=error_cube2*correct_unit2 # 应用单位转换
# 连接SL1和SL2数据
data_cube=np.concatenate((data_cube2,data_cube1),axis=0) # 沿第一个轴连接数据
error_cube=np.concatenate((error_cube2,error_cube1),axis=0) # 沿第一个轴连接不确定性数据
# 重新排序x和data_cube
xcor = x.argsort() # 获取排序索引
data_cube = data_cube[xcor] # 根据索引重新排序数据
error_cube = error_cube[xcor] # 根据索引重新排序不确定性数据
x=np.sort(x) # 对波长进行排序
# 立方体维度和修剪
xsize, ysize, zsize = data_cube.shape # 获取立方体的维度
ytrim=0; ysize=ysize-ytrim # 修剪y维度
ztrim=3; zsize=zsize-ztrim # 修剪z维度
print('Trimmed cube dimensions:', xsize, ysize, zsize) # 输出修剪后的维度
# 压缩的2D图像
cube_2dflux=np.sum(data_cube,axis=0)[0:ysize,0:zsize] # 沿z轴求和生成2D图像
# 压缩的1D光谱
cube_1dflux=np.zeros(xsize,) # 初始化1D光谱数组
for a in arange(0,ysize): # 遍历y维度
for b in arange(0,zsize): # 遍历z维度
spec_pix,err_pix=extract_spec(a,b) # 提取光谱和不确定性
cube_1dflux=cube_1dflux+spec_pix # 累加光谱数据
# 绘图
f,(ax1,ax2)=plt.subplots(1,2, figsize=(10,5)) # 创建子图
ax1.set_title('log Flux (Sum over Wavelength)') # 设置标题
ax1.imshow(cube_2dflux, origin='lower', cmap='gray', norm=LogNorm()) # 显示2D图像
ax2.set_title('Flux (Sum over Spaxels)') # 设置标题
ax2.plot(x,cube_1dflux) # 绘制1D光谱
ax2.set_xlabel(r"Wavelength $(\mu m)$") # 设置x轴标签
ax2.set_ylabel('Flux (Jy)') # 设置y轴标签
plt.show() # 显示图形
模型参数起始值#
设置光谱模型组件的起始参数。使用仅有两个自由参数的PAH尘埃模型（7.6微米和11.3微米的PAH通量）。其余的PAH比率由这个“pahdust”模型固定，以便整个数据立方体的拟合不会花费太长时间。稍后，我们将对所有PAH幅度自由的区域光谱进行拟合。
# 功率律模型
p_law = Model(powerlaw, prefix='pwl_') # 创建功率律模型，前缀为'pwl_'
generic_pars = p_law.make_params() # 生成模型的参数
setpar(generic_pars, 'pwl_c1', 0.0002, True, None) # 设置参数pwl_c1的初始值
setpar(generic_pars, 'pwl_c2', 0.0015, True, None) # 设置参数pwl_c2的初始值
# PAH尘埃模型，固定PAH特征比率
prefix = 'pahdust1_' # 设置前缀为'pahdust1_'
pahs = Model(pahdust, prefix=prefix) # 创建PAH尘埃模型
generic_pars.update(pahs.make_params()) # 更新参数
pahdust_defpar('pahdust1', 1.00506, 0.003, 0.005, generic_pars) # 定义PAH尘埃模型的参数
# 硅酸盐尘埃发射模型
prefix = 'sidust1_' # 设置前缀为'sidust1_'
silicate = Model(sidust, prefix=prefix) # 创建硅酸盐尘埃模型
generic_pars.update(silicate.make_params()) # 更新参数
sidust_defpar('sidust1', 190., 1.E-4, generic_pars) # 定义硅酸盐尘埃模型的参数
# 使用Smith et al 2007的波长的高斯模型
# H2纯旋转线
h2lines = ['55', '61', '69', '80', '96', '122'] # H2线的编号
h2wavel = [5.511, 6.109, 6.909, 8.026, 9.665, 12.278] # H2线对应的波长
amp_guess = 0.05 # 初始幅度猜测
std_guess = 0.04 # 初始标准差猜测
gaussH2 = Model(Modelcero) # 创建高斯模型的基础
for line in h2lines: # 遍历每一条H2线
prefix = 'h' + line + '_' # 设置前缀
gaussH2 += Model(gaussian, prefix=prefix) # 将高斯模型添加到gaussH2中
generic_pars.update(gaussH2.make_params()) # 更新参数
for h2lin, h2w in zip(h2lines, h2wavel): # 遍历H2线和波长
gaussian_defpar('h' + h2lin, amp_guess, h2w, std_guess, generic_pars) # 定义高斯参数
# 原子线模型
atomiclines = ['ArII', 'SIV', 'NeII'] # 原子线的名称
atomicwavel = [6.985, 10.511, 12.813] # 原子线对应的波长
amp_guess = 0.05 # 初始幅度猜测
std_guess = 0.05 # 初始标准差猜测
gaussAt = Model(Modelcero) # 创建原子线高斯模型的基础
for line in atomiclines: # 遍历每一条原子线
prefix = line + '_' # 设置前缀
gaussAt += Model(gaussian, prefix=prefix) # 将高斯模型添加到gaussAt中
generic_pars.update(gaussAt.make_params()) # 更新参数
for atlin, atw in zip(atomiclines, atomicwavel): # 遍历原子线和波长
gaussian_defpar(atlin, amp_guess, atw, std_guess, generic_pars) # 定义高斯参数
拟合一个代表性光谱像素#
在立方体中拟合一个单一的光谱像素，并使用拟合参数来调整通用模型的起始参数。
# 选择特定的spaxel（空间像素）
spec = data_cube[:, 8, 16] # 从数据立方体中提取第8行第16列的光谱数据
specerr = error_cube[:, 8, 16] # 从误差立方体中提取对应的误差数据
# 复合模型（包含功率律 + 硅酸盐 + PAHs + H2 + 原子线，不考虑消光）
generic_mod = p_law + silicate + pahs + gaussH2 + gaussAt # 定义复合模型
print(generic_mod.param_names) # 打印模型参数名称
# 拟合并绘制图形
spax_result = generic_mod.fit(spec, generic_pars, weights=1. / specerr, x=x) # 拟合光谱数据
plot_fit(x, spec, specerr, spax_result) # 绘制拟合结果图
# 保存结果
save_modelresult(spax_result, 'OnePointResult.sav') # 将拟合结果保存到文件中
# 更新通用模型的起始参数
for par in spax_result.params: # 遍历拟合结果中的每个参数
value = spax_result.params[par].value # 获取参数的值
vary = spax_result.params[par].vary # 获取参数是否可变的标志
minus = spax_result.params[par].min # 获取参数的最小值
if minus == float("-inf"): # 如果最小值为负无穷
minus = None # 将最小值设为None
generic_pars[par].set(value, vary=vary, min=minus) # 更新通用模型的参数
观察到的通量密度图（黑色），模型拟合（红色），以及残差（底部面板）。该像素在9.6微米处有强烈的H2发射，以及来自星形成区域的8微米PAH和12.8微米[Ne II]发射。
拟合整个数据立方体#
识别NaN值，然后创建模型列表（基于通用模型），每个无NaN的spaxel对应一个条目。接下来启动多进程池，限制为可用核心数减去一个。拟合整个数据立方体（900个spaxel）大约需要2分钟，使用7个进程并行运行。
cube_res = mp.Manager().dict() # 创建一个共享的字典，用于存储结果
# The fit method is redefined as a top-level function to make it pickle-able for multiprocessing.Pool
def modfit(spaxel, generic_pars, weights, x):
# 使用给定的参数拟合spaxel数据，并返回序列化的模型
modf = generic_mod.fit(spaxel, generic_pars, weights=1./spaxerr, x=x)
return modf.dumps() # 返回序列化的模型
# Helper to either assign a NaN result or fit a single non-NaN spaxel
def fit_point(args):
a, b, spaxel, spaxerr, x, generic_pars, nancheck = args # 解包参数
# print(a, b, len(spaxel), mp.current_process().name)
if nancheck == False: # 如果没有NaN值
res_point = modfit(spaxel, generic_pars, weights=1./spaxerr, x=x) # 拟合数据
cube_res[(a, b)] = res_point # 将结果存入共享字典
elif nancheck == True: # 如果存在NaN值
cube_res[(a, b)] = float('nan') # 将结果设为NaN
return # 返回
# Extraction region dimensions (trimmed cube)
astart = 0; aend = ysize # 定义提取区域的维度
bstart = 0; bend = zsize # 定义提取区域的维度
# Check for NaNs, select data, and set model start parameters for each spaxel
pooldata = [] # 初始化用于存储待处理数据的列表
for a in np.arange(astart, aend): # 遍历a维度
for b in np.arange(bstart, bend): # 遍历b维度
nancheck = False # 初始化NaN检查标志
for point in data_cube[:, a, b]: # 检查每个spaxel中的数据点
if np.isnan(point) == True or point == 0: # 如果发现NaN或零
nancheck = True # 设置NaN检查标志为True
# Copy spaxel data vectors
spaxel = 1.0 * np.array(data_cube[:, a, b]) # 复制spaxel数据
spaxerr = 1.0 * np.array(error_cube[:, a, b]) # 复制误差数据
wavelen = 1.0 * np.array(x) # 复制波长数据
pooldata.append([a, b, spaxel, spaxerr, wavelen, generic_pars, nancheck]) # 将数据添加到待处理列表
# Launch Multiprocessing Pool
start_time = time.time() # 记录开始时间
if __name__ == '__main__': # 确保在主程序中运行
with Pool(mp.cpu_count() - 1) as pool: # 创建进程池
pool.map(fit_point, pooldata) # 并行处理数据
print('Time count') # 打印时间计数
print("--- %s seconds ---" % (time.time() - start_time)) # 打印处理时间
# Save model fit
with open('MPF_cube_result_model.sav', 'wb') as fp: # 打开文件以保存模型
pickle.dump(list(cube_res.items()), fp) # 将结果序列化并保存
# Load model fit to full cube
cube_res = mp.Manager().dict() # 创建一个新的共享字典
with open('MPF_cube_result_model.sav', 'rb') as fp: # 打开文件以加载模型
restore_cube = pickle.load(fp) # 反序列化加载结果
for line in restore_cube: # 遍历加载的结果
pos = line[0] # 获取位置
cube_res[pos] = line[1] # 将结果存入共享字典
总通量和简化卡方（Reduced Chi^2）图#
显示总通量图像（沿波长方向压缩的观测通量立方体）和模型拟合的简化卡方图像，位于方形子区域内。
# 计算每个spaxel的观测通量总和
cube_tflux = np.sum(data_cube, axis=0)[0:ysize, 0:zsize] # 在第0轴上求和，得到的结果是一个二维数组
# 计算每个spaxel的模型通量总和
tcube_chi = np.zeros((ysize, zsize)) # 初始化一个与观测通量相同大小的数组，用于存储卡方值
tcube_modflux = np.zeros((ysize, zsize)) # 初始化一个与观测通量相同大小的数组，用于存储模型通量
funcdefs = {} # 定义一个空字典，用于存储函数定义
modres = spax_result # 将spax_result赋值给modres，方便后续使用
for pos, dumpval in cube_res.items(): # 遍历cube_res字典中的每个位置和对应值
# for pos, dumpval in cube_res: # 这行代码被注释掉，表示可以直接使用上面的for循环
if str(cube_res[pos]) == 'nan': # 检查当前值是否为NaN
tcube_modflux[pos] = float('nan') # 如果是NaN，则模型通量也设为NaN
tcube_chi[pos] = float('nan') # 如果是NaN，则卡方值也设为NaN
else:
result = modres.loads(dumpval, funcdefs=funcdefs) # 加载当前位置的模型结果
tcube_modflux[pos] = np.sum(result.best_fit) # 计算模型通量的总和并存储
tcube_chi[pos] = result.redchi # 存储卡方值
# 创建一个包含3个子图的图形
f, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(10, 15)) # 创建1行3列的子图
ax1.set_title('log Observed Flux') # 设置第一个子图的标题
ax1.imshow(cube_tflux, origin='lower', cmap='gray', norm=LogNorm()) # 显示观测通量的对数图像
ax2.set_title('log Model Flux') # 设置第二个子图的标题
ax2.imshow(tcube_modflux, origin='lower', cmap='gray', norm=LogNorm()) # 显示模型通量的对数图像
ax3.set_title('Reduced Chi Square') # 设置第三个子图的标题
ax3.imshow(cube_tflux - tcube_modflux, origin='lower', cmap='hot') # 显示观测通量与模型通量之差的热图
plt.show() # 显示所有子图
折叠立方体图，显示观测通量（observed flux）、模型通量（model flux）和卡方/自由度（Chi^2/DF）。 输入数据立方体中存在NaN的白色区域，在输出模型立方体的左侧和顶部边缘以及第26列中也存在NaN。 核心区域的残差可能是由于PSF采样不足或AGN模型拟合不佳所导致的。
线和PAH特征通量图。#
从模型参数立方体中提取谱线和PAH图的函数#
def line_map(line, cube_res, modres):
"""创建高斯发射线模型分量的总发射线通量图"""
line_flux = np.zeros((35, 35)) # 初始化一个35x35的零数组用于存储发射线通量
for pos, dumpval in cube_res.items(): # 遍历立方体结果中的每个位置和对应值
if str(cube_res[pos]) == 'nan': # 检查当前值是否为NaN
line_flux[pos] = float('nan') # 如果是NaN，则在通量图中设置为NaN
else:
result = modres.loads(dumpval, funcdefs=funcdefs) # 加载模型结果
prefix = line + '_' # 创建前缀以便提取参数
amp, cen, std = gaussian_extractpars(prefix, result) # 提取高斯模型参数
line_flux[pos] = gaussianline_flux(amp, std, cen) # 计算并存储发射线通量
return line_flux # 返回发射线通量图
def pahdust_maps(pahdust_comp, cube_res, modres):
"""创建高斯发射线模型分量的总发射线通量图"""
pah_ion_flux = np.zeros((35, 35)) # 初始化一个35x35的零数组用于存储离子通量
pah_neutral_flux = np.zeros((35, 35)) # 初始化一个35x35的零数组用于存储中性通量
for pos, dumpval in cube_res.items(): # 遍历立方体结果中的每个位置和对应值
if str(cube_res[pos]) == 'nan': # 检查当前值是否为NaN
pah_ion_flux[pos] = float('nan') # 如果是NaN，则在离子通量图中设置为NaN
pah_neutral_flux[pos] = float('nan') # 如果是NaN，则在中性通量图中设置为NaN
else:
result = modres.loads(dumpval, funcdefs=funcdefs) # 加载模型结果
prefix = pahdust_comp + '_' # 创建前缀以便提取参数
ampl_76, ampl_113 = pahdust_extractpars(prefix, result) # 提取PAH尘埃的参数
pah_76_flux = drudeline_flux(ampl_76, 0.044, 7.60) # 计算76微米的通量
pah_113_flux = drudeline_flux(ampl_113, 0.032, 11.33) # 计算113微米的通量
pah_ion_flux[pos] = pah_76_flux # 存储76微米的离子通量
pah_neutral_flux[pos] = pah_113_flux # 存储113微米的中性通量
return pah_ion_flux, pah_neutral_flux # 返回离子和中性通量图
def sidust_map(sidust_comp, cube_res, modres):
"""创建硅酸盐尘埃发射图"""
sidust_flux = np.zeros((35, 35)) # 初始化一个35x35的零数组用于存储硅酸盐尘埃通量
for pos, dumpval in cube_res.items(): # 遍历立方体结果中的每个位置和对应值
if str(cube_res[pos]) == 'nan': # 检查当前值是否为NaN
sidust_flux[pos] = float('nan') # 如果是NaN，则在通量图中设置为NaN
else:
result = modres.loads(dumpval, funcdefs=funcdefs) # 加载模型结果
prefix = sidust_comp + '_' # 创建前缀以便提取参数
t, si_ampl = sidust_extractpars(prefix, result) # 提取硅酸盐尘埃的参数
si_norm = 1.0 # 设置标准化因子
sidust_flux[pos] = si_ampl * si_norm # 计算并存储硅酸盐尘埃通量
return sidust_flux # 返回硅酸盐尘埃通量图
H2 9.6 微米, [Ne II] 12.8 微米, PAH 7.6, PAH 11.3, 以及来自模型的硅酸盐尘埃发射图#
# 计算H2 9.6微米的光谱图
h2_flux = line_map('h96', cube_res, modres)
# 计算Ne II 12.8微米的光谱图
ne_flux = line_map('NeII', cube_res, modres)
# 计算PAH离子和中性光谱图
pah_ion_flux, pah_neutral_flux = pahdust_maps('pahdust1', cube_res, modres)
# 计算硅尘埃的光谱图
sidust_flux = sidust_map('sidust1', cube_res, modres)
# 创建一个1行5列的子图，设置图像大小
f, (ax1, ax2, ax3, ax4, ax5) = plt.subplots(1, 5, figsize=(10, 20))
# 设置Ne II 12.8微米的标题
ax1.set_title('[Ne II] 12.8 um')
# 显示Ne II光谱图
ax1.imshow(ne_flux, origin='lower')
# 设置H2 9.6微米的标题
ax2.set_title('H2 9.6 um')
# 显示H2光谱图
ax2.imshow(h2_flux, origin='lower')
# 设置PAH 7.6微米的标题
ax3.set_title('PAH 7.6 um')
# 显示PAH离子光谱图
ax3.imshow(pah_ion_flux, origin='lower')
# 设置PAH 11.3微米的标题
ax4.set_title('PAH 11.3 um')
# 显示PAH中性光谱图
ax4.imshow(pah_neutral_flux, origin='lower')
# 设置硅尘埃发射的标题
ax5.set_title('Si dust em.')
# 显示硅尘埃光谱图
ax5.imshow(sidust_flux, origin='lower')
# 显示所有子图
plt.show()
# 创建RGB图像
# 设置拉伸因子
stretch1 = 0.000000000000004
stretch2 = 0.000000000000005
# 生成RGB图像，包含Ne II、H2和PAH中性光谱
m58_neii_pah0_h2 = make_lupton_rgb(2.5 * ne_flux, h2_flux, pah_neutral_flux, minimum=0, Q=5, stretch=stretch1, filename="m58_NeII_H96_PAH113.png")
# 生成另一个RGB图像，包含PAH中性、H2和PAH离子光谱
m58_pah_h2 = make_lupton_rgb(pah_neutral_flux, h2_flux, 1.2 * pah_ion_flux, minimum=0, Q=5, stretch=stretch2, filename="m58_PAH113_H96_PAH76.png")
# 创建一个2面板的图形
f, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 15))
# 关闭第一个面板的坐标轴
ax1.axis('off')
# 显示第一个RGB图像
ax1.imshow(m58_neii_pah0_h2, origin='lower')
# 关闭第二个面板的坐标轴
ax2.axis('off')
# 显示第二个RGB图像
ax2.imshow(m58_pah_h2, origin='lower')
# 显示所有面板
plt.show()
3色RGB图像。顶部：单个特征图。左下：红色 = [Ne II] 12.8微米，绿色 = H2 S(3) 9.6微米，蓝色 = PAH 11.3微米。活动星系核中的电离原子气体显示为红色，冲击区域为绿色，尘埃发射为蓝色。右下：红色 = PAH 11.3微米，绿色 = H2 S(3) 9.6微米，蓝色 = PAH 7.6微米。中性PAH发射显示为红色，来自星形成区的光致电离区（PDRs）的电离PAH发射为蓝色，冲击区域为绿色。
提取和建模立方体子区域#
制作光谱区域掩膜#
根据空间位置、线通量和线或特征比率定义光谱提取区域。
# 定义一个5x5的核提取区域
nuc_mask=np.zeros((35, 35)) # 创建一个35x35的零矩阵作为核掩膜
for a in arange(-2,3): # 遍历-2到2的范围
for b in arange(-2,3): # 遍历-2到2的范围
nuc_mask[14-a,17-b]=int(1) # 在指定位置设置掩膜值为1
# 螺旋臂提取区域
arm_mask=np.where(((pah_neutral_flux>=0.155e-14) &((pah_ion_flux/pah_neutral_flux)>0.8) & (nuc_mask==0)),1,0) # 根据条件生成螺旋臂掩膜
# 内盘提取区域
pah_mask=np.where(((pah_neutral_flux>=0.155e-14) & (h2_flux<1e-15) & (nuc_mask==0)),1,0) - arm_mask # 根据条件生成内盘掩膜
pah_mask[:,0:3]=int(0) # 修剪立方体的边缘
pah_mask[28:,:]=int(0) # 修剪立方体的底部
# 最强H2发射区域
h2s_mask=np.where(((h2_flux>=1e-15)&((h2_flux/pah_neutral_flux)>0.7) & (nuc_mask==0)),1,0) # 根据条件生成强H2掩膜
# H2提取区域
h2_mask=np.where(((h2_flux>=1e-15) & (nuc_mask==0)),1,0)-h2s_mask # 根据条件生成H2掩膜
# M58组合中心区域掩膜（排除螺旋臂）
m58_mask=nuc_mask + h2_mask*2 + h2s_mask*3 + pah_mask*4 +arm_mask*5 # 组合所有掩膜
# 创建子图
f,((ax1,ax2,ax3,ax4,ax5,ax6))=plt.subplots(1,6,figsize=(15,10)) # 创建1行6列的子图
ax1.set_title('Nucleus') # 设置标题为“核”
ax1.imshow(nuc_mask,origin='lower') # 显示核掩膜
ax2.set_title('Spiral Arm') # 设置标题为“螺旋臂”
ax2.imshow(arm_mask,origin='lower') # 显示螺旋臂掩膜
ax3.set_title('Inner Disk') # 设置标题为“内盘”
ax3.imshow(pah_mask,origin='lower') # 显示内盘掩膜
ax4.set_title('H2') # 设置标题为“H2”
ax4.imshow(h2_mask,origin='lower') # 显示H2掩膜
ax5.set_title('H2-strongest') # 设置标题为“最强H2”
ax5.imshow(h2s_mask,origin='lower') # 显示强H2掩膜
ax6.set_title('All regions') # 设置标题为“所有区域”
ax6.imshow(m58_mask,origin='lower') # 显示所有区域掩膜
plt.show() # 显示图形
基于空间和发射线或PAH特征强度的提取区域掩膜。
从掩膜区域提取光谱#
五个不同区域的光谱汇总：螺旋臂 (Spiral Arm)、核 (Nucleus)、多环芳烃 (PAH)、氢分子 (H2) 和最强氢分子 (Strongest H2)
nuc_spec = np.zeros(xsize,) # 初始化核区域光谱数组
nuc_err = np.zeros(xsize,) # 初始化核区域误差数组
pah_spec = np.zeros(xsize,) # 初始化PAH区域光谱数组
pah_err = np.zeros(xsize,) # 初始化PAH区域误差数组
h2_spec = np.zeros(xsize,) # 初始化H2区域光谱数组
h2_err = np.zeros(xsize,) # 初始化H2区域误差数组
arm_spec = np.zeros(xsize,) # 初始化螺旋臂区域光谱数组
arm_err = np.zeros(xsize,) # 初始化螺旋臂区域误差数组
h2s_spec = np.zeros(xsize,) # 初始化H2强区域光谱数组
h2s_err = np.zeros(xsize,) # 初始化H2强区域误差数组
for a in arange(0, ysize): # 遍历y轴
for b in arange(0, zsize): # 遍历z轴
if nuc_mask[a, b] == 1: # 如果核区域掩膜为1
spec_pix, err_pix = extract_spec(a, b) # 提取光谱和误差
nuc_spec = nuc_spec + spec_pix # 累加光谱
nuc_err = nuc_err + err_pix**2 # 累加误差的平方
if pah_mask[a, b] == 1: # 如果PAH区域掩膜为1
spec_pix, err_pix = extract_spec(a, b) # 提取光谱和误差
pah_spec = pah_spec + spec_pix # 累加光谱
pah_err = pah_err + err_pix**2 # 累加误差的平方
if h2_mask[a, b] == 1: # 如果H2区域掩膜为1
spec_pix, err_pix = extract_spec(a, b) # 提取光谱和误差
h2_spec = h2_spec + spec_pix # 累加光谱
h2_err = h2_err + err_pix**2 # 累加误差的平方
if arm_mask[a, b] == 1: # 如果螺旋臂区域掩膜为1
spec_pix, err_pix = extract_spec(a, b) # 提取光谱和误差
arm_spec = arm_spec + spec_pix # 累加光谱
arm_err = arm_err + err_pix**2 # 累加误差的平方
if h2s_mask[a, b] == 1: # 如果H2强区域掩膜为1
spec_pix, err_pix = extract_spec(a, b) # 提取光谱和误差
h2s_spec = h2s_spec + spec_pix # 累加光谱
h2s_err = h2s_err + err_pix**2 # 累加误差的平方
# 计算误差的平方根
nuc_err = sqrt(nuc_err)
pah_err = sqrt(pah_err)
h2_err = sqrt(h2_err)
arm_err = sqrt(arm_err)
h2s_err = sqrt(h2s_err)
# 创建子图
f, ((ax1, ax2, ax3), (ax4, ax5, ax6)) = plt.subplots(2, 3, figsize=(15, 10))
ax1.set_title('Spiral Arm') # 设置螺旋臂子图标题
ax1.errorbar(x, arm_spec, yerr=arm_err, c='gold') # 绘制螺旋臂光谱及误差
ax1.set_ylabel('Flux (Jy)') # 设置y轴标签
ax2.set_title('Nucleus') # 设置核区域子图标题
ax2.errorbar(x, nuc_spec, yerr=nuc_err, c='navy') # 绘制核区域光谱及误差
ax3.set_title('Inner Disk') # 设置内盘子图标题
ax3.errorbar(x, pah_spec, yerr=pah_err, c='limegreen') # 绘制内盘光谱及误差
ax4.set_title('H2') # 设置H2子图标题
ax4.errorbar(x, h2_spec, yerr=h2_err, c='steelblue') # 绘制H2光谱及误差
ax4.set_xlabel(r"Wavelength $(\mu m)$") # 设置x轴标签
ax4.set_ylabel('Flux (Jy)') # 设置y轴标签
ax5.set_title('All regions') # 设置所有区域子图标题
ax5.imshow(m58_mask, origin='lower') # 显示掩膜图像
ax5.set_xlabel('Pixel') # 设置x轴标签
ax6.set_title('H2-Strongest') # 设置H2强区域子图标题
ax6.errorbar(x, h2s_spec, yerr=h2s_err, c='teal') # 绘制H2强区域光谱及误差
ax6.set_xlabel(r"Wavelength $(\mu m)$") # 设置x轴标签
plt.show() # 显示图形
将光谱保存到文件中#
# 将核区域的光谱数据保存到文件，文件名为targname+'_spec_nuc.dat'
# np.savetxt(targname+'_spec_nuc.dat', np.column_stack((x, nuc_spec, nuc_err)), header='Wavelength[microns] Flux_Nucleus[Jy] Err_Flux_Nucleus')
# 将H2区域的光谱数据保存到文件，文件名为targname+'_spec_h2.dat'
# np.savetxt(targname+'_spec_h2.dat', np.column_stack((x, h2_spec, h2_err)), header='Wavelength[microns] Flux_H2[Jy] Err_H2')
# 将H2S区域的光谱数据保存到文件，文件名为targname+'_spec_h2s.dat'
# np.savetxt(targname+'_spec_h2s.dat', np.column_stack((x, h2_spec, h2_err)), header='Wavelength[microns] Flux_H2s[Jy] Err_H2s')
# 将PAH区域的光谱数据保存到文件，文件名为targname+'_spec_pah.dat'
# np.savetxt(targname+'_spec_pah.dat', np.column_stack((x, pah_spec, pah_err)), header='Wavelength[microns] Flux_PAH[Jy] Err_PAH')
# 将Arm区域的光谱数据保存到文件，文件名为targname+'_spec_arm.dat'
# np.savetxt(targname+'_spec_arm.dat', np.column_stack((x, arm_spec, arm_err)), header='Wavelength[microns] Flux_Arm[Jy] Err_Arm')
区域光谱拟合#
通用模型组件#
对14个独立的PAH（多环芳香烃）组件进行建模，所有幅度均为自由变量。对整个立方体进行这样的计算将会非常耗费计算资源，但对于汇总区域的光谱，仅需几秒钟。
# 功率律模型
p_law = Model(powerlaw, prefix='pwl_') # 创建功率律模型
generic_pars = p_law.make_params() # 生成模型参数
setpar(generic_pars, 'pwl_c1', 0.005, True, None) # 设置参数pwl_c1
setpar(generic_pars, 'pwl_c2', 0.0015, True, None) # 设置参数pwl_c2
# 单个PAH（多环芳烃）模型
drudes = Model(Modelcero) # 创建一个空模型
allpahlines = ['52', '62', '74', '76', '78', '83', '86', '112', '113', '119', '126', '134', '140', '141'] # PAH线列表
for line in allpahlines: # 遍历所有PAH线
prefix = 'PAH' + line + '_' # 为每个PAH线创建前缀
drudes += Model(drude, prefix=prefix) # 将每个PAH线模型添加到drudes模型中
generic_pars.update(drudes.make_params()) # 更新参数
# 定义每个PAH的参数
drude_defpar('PAH52', 0.05, 5.27, 0.034, generic_pars) # PAH52的参数
drude_defpar('PAH62', 0.05, 6.22, 0.030, generic_pars) # PAH62的参数
drude_defpar('PAH74', 0.05, 7.42, 0.126, generic_pars) # PAH74的参数
drude_defpar('PAH76', 0.05, 7.60, 0.044, generic_pars) # PAH76的参数
drude_defpar('PAH78', 0.5, 7.85, 0.053, generic_pars) # PAH78的参数
drude_defpar('PAH83', 0.5, 8.33, 0.05, generic_pars) # PAH83的参数
drude_defpar('PAH86', 0.5, 8.61, 0.039, generic_pars) # PAH86的参数
drude_defpar('PAH112', 0.5, 11.23, 0.012, generic_pars) # PAH112的参数
drude_defpar('PAH113', 0.5, 11.33, 0.032, generic_pars) # PAH113的参数
drude_defpar('PAH119', 0.5, 11.99, 0.045, generic_pars) # PAH119的参数
drude_defpar('PAH126', 0.5, 12.62, 0.042, generic_pars) # PAH126的参数
drude_defpar('PAH134', 0.5, 13.48, 0.04, generic_pars) # PAH134的参数
drude_defpar('PAH140', 0.5, 14.04, 0.016, generic_pars) # PAH140的参数
drude_defpar('PAH141', 0.5, 14.19, 0.025, generic_pars) # PAH141的参数
# 硅酸盐尘埃发射模型
prefix = 'sidust1_' # 设置前缀
silicate = Model(sidust, prefix=prefix) # 创建硅酸盐模型
generic_pars.update(silicate.make_params()) # 更新参数
sidust_defpar('sidust1', 190., 1.E-4, generic_pars) # 定义硅酸盐的参数
# 采用Smith等人2007年的FWHM和X中心的高斯模型
# H2分子模型
h2lines = ['55', '61', '69', '80', '96', '122'] # H2线列表
gaussH2 = Model(Modelcero) # 创建一个空模型
for line in h2lines: # 遍历所有H2线
prefix = 'h' + line + '_' # 为每个H2线创建前缀
gaussH2 += Model(gaussian, prefix=prefix) # 将每个H2线模型添加到gaussH2模型中
generic_pars.update(gaussH2.make_params()) # 更新参数
# 定义每个H2线的参数
gaussian_defpar('h55', 0.05, 5.511, 0.02, generic_pars) # H2线55的参数
gaussian_defpar('h61', 0.05, 6.109, 0.02, generic_pars) # H2线61的参数
gaussian_defpar('h69', 0.05, 6.909, 0.02, generic_pars) # H2线69的参数
gaussian_defpar('h80', 0.05, 8.026, 0.04, generic_pars) # H2线80的参数
gaussian_defpar('h96', 0.05, 9.665, 0.04, generic_pars) # H2线96的参数
gaussian_defpar('h122', 0.05, 12.278, 0.04, generic_pars) # H2线122的参数
# 原子线模型
atomiclines = ['ArII', 'SIV', 'NeII'] # 原子线列表
gaussAt = Model(Modelcero) # 创建一个空模型
for line in atomiclines: # 遍历所有原子线
prefix = line + '_' # 为每个原子线创建前缀
gaussAt += Model(gaussian, prefix=prefix) # 将每个原子线模型添加到gaussAt模型中
generic_pars.update(gaussAt.make_params()) # 更新参数
# 定义每个原子线的参数
gaussian_defpar('ArII', 0.05, 6.985, 0.02, generic_pars) # ArII的参数
gaussian_defpar('SIV', 0.05, 10.511, 0.04, generic_pars) # SIV的参数
gaussian_defpar('NeII', 0.05, 12.813, 0.04, generic_pars) # NeII的参数
螺旋臂光谱拟合#
# 螺旋臂光谱
spec = arm_spec # 将光谱数据赋值给变量spec
specerr = arm_err # 将光谱误差数据赋值给变量specerr
# 模型（幂律 + 各个PAH + H2和原子线，无消光）
model = p_law + drudes + gaussAt + gaussH2 # 定义模型，包括幂律、德鲁德模型、Gauss模型和H2模型
# 调整幂律起始参数
setpar(generic_pars, 'pwl_c1', 0.01, True, None) # 设置幂律模型的第一个参数
setpar(generic_pars, 'pwl_c2', 0.001, True, None) # 设置幂律模型的第二个参数
# 模型结果
arm_result = model.fit(spec, generic_pars, weights=1./specerr, x=x) # 拟合模型到光谱数据
arm_model = arm_result.best_fit # 获取最佳拟合模型
arm_residual = spec - arm_model # 计算残差，即原始光谱与拟合模型的差
arm_modelcomps, arm_modlabels = model_comps(arm_result) # 获取模型组成部分及其标签
# 绘图
# plot_fit(x, spec, specerr, arm_result) # 绘制拟合结果（注释掉的代码）
模型拟合螺旋臂区域的综合光谱。
H2-最强区域光谱拟合#
# H2区域光谱
spec = h2s_spec # 定义光谱数据
specerr = h2s_err # 定义光谱误差
# 模型（功率律 + 单独的PAHs + H2和原子线，无消光）
model = p_law + drudes + gaussAt + gaussH2 # 定义模型，包括功率律、德鲁德模型和高斯模型
# 调整功率律起始参数
setpar(generic_pars, 'pwl_c1', 0.002, True, None) # 设置功率律参数c1
setpar(generic_pars, 'pwl_c2', 0.0015, True, None) # 设置功率律参数c2
# 模型结果
h2s_result = model.fit(spec, generic_pars, weights=1./specerr, x=x) # 拟合模型到光谱数据
h2s_model = h2s_result.best_fit # 获取最佳拟合模型
h2s_residual = spec - h2s_model # 计算残差
h2s_modelcomps, h2s_modlabels = model_comps(h2s_result) # 获取模型组成部分及其标签
# 绘图
# plot_fit(x, spec, specerr, h2s_result) # 绘制拟合结果（注释掉的代码）
模型拟合H2-最强区域的综合光谱。
核心光谱拟合#
# 核心光谱
spec = nuc_spec # 将核光谱赋值给变量spec
specerr = nuc_err # 将核光谱误差赋值给变量specerr
# 模型（功率律 + AGN硅酸盐发射 + 各个PAH + H2和原子线，无消光）
model = p_law + silicate + drudes + gaussH2 + gaussAt # 定义模型为多个成分的组合
# 调整功率律和硅酸盐尘埃的起始参数
setpar(generic_pars, 'pwl_c1', 0.005, True, None) # 设置功率律参数c1
setpar(generic_pars, 'pwl_c2', 0.01, True, None) # 设置功率律参数c2
sidust_defpar('sidust1', 190., 2.E-4, generic_pars) # 定义硅酸盐尘埃的参数
# 模型结果
nuc_result = model.fit(spec, generic_pars, weights=1./specerr, x=x) # 拟合模型并获取结果
nuc_model = nuc_result.best_fit # 获取最佳拟合模型
nuc_residual = spec - nuc_model # 计算残差
nuc_modelcomps, nuc_modlabels = model_comps(nuc_result) # 获取模型成分和标签
# 绘图
plot_fit(x, spec, specerr, nuc_result) # 绘制拟合结果
模型拟合到核，包含固定温度的活动星系核（AGN）硅酸盐发射。
nuc_result # 变量nuc_result存储拟合结果
# result.fit_report(show_correl=False) # 显示拟合报告，不显示相关性
核谱的拟合参数
概述图#
# 创建一个2行3列的子图，设置图形大小为15x10
f,((ax1,ax2,ax3),(ax4,ax5,ax6))=plt.subplots(2,3,figsize=(15,10))
# 设置ax1的标题为'Spiral Arm'
ax1.set_title('Spiral Arm')
# 设置y轴标签为'Flux (Jy)'
ax1.set_ylabel('Flux (Jy)')
# 在ax1上绘制带误差条的光谱数据，颜色为金色
ax1.errorbar(x,arm_spec,yerr=arm_err, c='gold')
# 定义模型颜色列表
mod_colors=['brown','g','b','magenta','orange']
# 遍历每个模型组件，标签和颜色，并在ax1上绘制
for mcomp,mlabl,mcolr in zip(arm_modelcomps,arm_modlabels,mod_colors):
ax1.plot(x,mcomp,label=mlabl,c=mcolr)
# 设置ax2的标题为'Nucleus'
ax2.set_title('Nucleus')
# 在ax2上绘制带误差条的核光谱数据，颜色为海军蓝
ax2.errorbar(x,nuc_spec,yerr=nuc_err, c='navy')
# 在ax2上绘制模型（已注释掉）
#ax2.plot(x, fit_model,label='Model',c='red')
# 遍历每个模型组件，标签和颜色，并在ax2上绘制
for mcomp,mlabl,mcolr in zip(nuc_modelcomps,nuc_modlabels,mod_colors):
ax2.plot(x,mcomp,label=mlabl,c=mcolr)
# 设置ax3的标题为'Inner Disk'
ax3.set_title('Inner Disk')
# 在ax3上绘制带误差条的PAH光谱数据，颜色为青柠绿
ax3.errorbar(x,pah_spec,yerr=pah_err, c='limegreen')
# 设置ax4的标题为'H2'
ax4.set_title('H2')
# 设置ax4的x轴标签为'Wavelength $(\mu m)$'
ax4.set_xlabel(r"Wavelength $(\mu m)$")
# 设置ax4的y轴标签为'Flux (Jy)'
ax4.set_ylabel('Flux (Jy)')
# 在ax4上绘制带误差条的H2光谱数据，颜色为钢蓝色
ax4.errorbar(x,h2_spec,yerr=h2_err, c='steelblue')
# 设置ax5的标题为'All regions'
ax5.set_title('All regions')
# 设置ax5的x轴标签为'Pixel'
ax5.set_xlabel('Pixel')
# 在ax5上显示m58_mask图像，原点在下方
ax5.imshow(m58_mask,origin='lower')
# 设置ax6的标题为'H2-Strongest'
ax6.set_title('H2-Strongest')
# 设置ax6的x轴标签为'Wavelength $(\mu m)$'
ax6.set_xlabel(r"Wavelength $(\mu m)$")
# 在ax6上绘制带误差条的H2强光谱数据，颜色为水鸭色
ax6.errorbar(x,h2s_spec,yerr=h2s_err,c='teal')
# 遍历每个模型组件，标签和颜色，并在ax6上绘制
for mcomp,mlabl,mcolr in zip(h2s_modelcomps,arm_modlabels,mod_colors):
ax6.plot(x,mcomp,label=mlabl,c=mcolr)
# 显示所有绘制的图形
plt.show()
由Patrick Ogle和Ivan Lopez创建的笔记本。