WFSS光谱 - 第2部分：交叉相关以确定红移

WFSS光谱 - 第2部分：交叉相关以确定红移#

用例： 最优提取光谱；红移测量；发射线图。简化版的 JDox 科学用例 # 33。

数据： JWST 模拟的 NIRISS 图像来自 MIRAGE，经过 JWST 校准管道 处理；星系团。

工具： specutils, astropy, pandas, emcee, lmfit, corner, h5py。

跨仪器： NIRSpec

文档： 本笔记本是 STScI 更大 后处理数据分析工具生态系统 的一部分。

引言#

本笔记本是关于 NIRISS WFSS 数据的一组中的第 3 个，共 4 个：

  1. 1D 最优提取，因为 JWST 管道仅提供框提取。最优提取提高了微弱源光谱的信噪比（S/N）。

  2. 合并并归一化 1D 光谱。

  3. 与模板进行交叉相关以获取红移。

  4. 空间分辨的发射线图。

本笔记本通过使用 specutils 模板相关性推导出具有多个发射线的星系的红移。该笔记本从第 2 个笔记本中获得的最优提取的 1D 光谱开始。

注意： 没有发射线的光谱（例如，前一个笔记本中的 ID 00003）可能会在此函数中失败。

%matplotlib inline  # 在Jupyter Notebook中内联显示Matplotlib图形

import os  # 导入操作系统模块，用于文件和目录操作

import numpy as np  # 导入NumPy库，用于数值计算和数组操作

import h5py  # 导入h5py库，用于处理HDF5文件格式

from astropy.io import fits  # 从Astropy库导入FITS模块，用于处理FITS文件

from astropy.table import QTable  # 从Astropy库导入QTable，用于处理表格数据

import astropy.units as u  # 导入Astropy单位模块，用于物理量的单位处理

from astropy.nddata import StdDevUncertainty  # 导入标准偏差不确定性类，用于处理不确定性数据

from astropy.modeling.polynomial import Chebyshev1D  # 导入Chebyshev多项式模型，用于数据拟合

from astropy import constants as const  # 导入Astropy常数模块，用于物理常数的使用

import astropy  # 导入Astropy库

print('astropy', astropy.__version__)  # 打印Astropy库的版本信息

import matplotlib.pyplot as plt  # 导入绘图库 matplotlib 的 pyplot 模块

import matplotlib as mpl  # 导入 matplotlib 库

# 设置保存图像时的分辨率为 80 dpi
mpl.rcParams['savefig.dpi'] = 80

# 设置图像显示时的分辨率为 80 dpi
mpl.rcParams['figure.dpi'] = 80

import specutils  # 导入specutils库

from specutils.fitting import continuum  # 从specutils库中导入连续谱拟合模块

from specutils.spectra.spectrum1d import Spectrum1D  # 从specutils库中导入一维光谱类

from specutils.analysis import correlation  # 从specutils库中导入相关性分析模块

from specutils.spectra.spectral_region import SpectralRegion  # 从specutils库中导入光谱区域类

print("Specutils: ", specutils.__version__)  # 打印specutils库的版本

0. 下载笔记本 01 产品#

这些也可以通过运行笔记本获得。

import os  # 导入os模块，用于文件和目录操作

# 检查输出目录是否存在
if not os.path.exists('./output'):
    
    import zipfile  # 导入zipfile模块，用于处理zip文件
    import urllib.request  # 导入urllib.request模块，用于下载文件

    # 定义要下载的zip文件的链接
    boxlink = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/NIRISS_lensing_cluster/output.zip'
    
    # 定义下载后保存的zip文件路径
    boxfile = './output.zip'
    
    # 下载zip文件
    urllib.request.urlretrieve(boxlink, boxfile)
    
    # 打开下载的zip文件
    zf = zipfile.ZipFile(boxfile, 'r')
    
    # 获取zip文件内所有文件的名称列表
    list_names = zf.namelist()
    
    # 遍历zip文件内的每个文件
    for item in list_names:
        # 提取当前文件，使用extract比extractall更安全，避免绝对路径或相对路径问题
        zf.extract(member=item) 
        
else:
    # 如果输出目录已存在，打印提示信息
    print('Already exists')

1. 打开最佳提取的1D光谱文本文件；#

这些是来自笔记本01a和01b的最佳提取、归一化的光谱。

DIR_OUT = './output/'  # 输出目录

filt = 'f200w'  # 选择的滤光片

grism = 'G150C'  # 选择的光栅

# grism = 'G150R'  # 另一种光栅（被注释掉）

id = '00004'  # 文件标识符

# 构建1D光谱文件的路径
file_1d = f'{DIR_OUT}l3_nis_{filt}_{grism}_s{id}_combine_1d_opt.fits'

# 打开FITS文件并读取数据部分
fd = fits.open(file_1d)[1].data

# 打开FITS文件并读取头部信息
hd = fits.open(file_1d)[1].header

fd  # 输出数据

# 归一化观测光谱；仅用于视觉目的

flux_normalize = 500.  # 设置归一化的光通量值为500

wave_200 = fd['wavelength']  # 从数据集中提取波长信息

flux_200 = fd['flux'] / flux_normalize  # 归一化光谱数据

flux_err_200 = fd['uncertainty'] / flux_normalize  # 归一化光谱不确定性数据

# 打开其他滤光片的文件。

filt = 'f150w'  # 设置滤光片为f150w

# 构建文件名，包含输出目录、滤光片、光栅和ID
file_1d = f'{DIR_OUT}l3_nis_{filt}_{grism}_s{id}_combine_1d_opt.fits'

# 打开FITS文件并读取数据
fd = fits.open(file_1d)[1].data

# 提取波长数据
wave_150 = fd['wavelength']

# 提取并归一化光谱数据
flux_150 = fd['flux'] / flux_normalize

# 提取并归一化光谱不确定性
flux_err_150 = fd['uncertainty'] / flux_normalize

# 

filt = 'f115w'  # 设置滤光片为f115w

# 构建文件名，包含输出目录、滤光片、光栅和ID
file_1d = f'{DIR_OUT}l3_nis_{filt}_{grism}_s{id}_combine_1d_opt.fits'

# 打开FITS文件并读取数据
fd = fits.open(file_1d)[1].data

# 提取波长数据
wave_115 = fd['wavelength']

# 提取并归一化光谱数据
flux_115 = fd['flux'] / flux_normalize

# 提取并归一化光谱不确定性
flux_err_115 = fd['uncertainty'] / flux_normalize

# 绘制所有数据

plt.figure()  # 创建一个新的图形

# 绘制F115W波段的光谱数据，带误差条
plt.errorbar(wave_115, flux_115, ls='-', color='lightblue', label='F115W')

# 绘制F150W波段的光谱数据，带误差条
plt.errorbar(wave_150, flux_150, ls='-', color='orange', label='F150W')

# 绘制F200W波段的光谱数据，带误差条
plt.errorbar(wave_200, flux_200, ls='-', color='purple', label='F200W')

# 填充F115W波段的误差区域
plt.fill_between(wave_115, -flux_err_115, flux_err_115, ls='-', color='gray', label='', alpha=0.3)

# 填充F150W波段的误差区域
plt.fill_between(wave_150, -flux_err_150, flux_err_150, ls='-', color='gray', label='', alpha=0.3)

# 填充F200W波段的误差区域
plt.fill_between(wave_200, -flux_err_200, flux_err_200, ls='-', color='gray', label='', alpha=0.3)

# 设置x轴标签
plt.xlabel('Wavelength / um')

# 设置y轴的范围
plt.ylim(0., 0.2)

# 显示图例，位置在左上角
plt.legend(loc=2)

连续谱是可见的。在模板相关之前，我们需要将其减去。#

  • 观察到强发射线。（Hb+Oiii 双线在 ~1.55μm，混合的 Ha+Nii 在 ~2.1μm）

  • 请注意，每个滤光片边缘的通量过剩并不真实，应予以屏蔽。

# F200W

spec_unit = u.MJy  # 定义光谱单位为MJy

# 创建掩膜，选择波长在1.75到1.97微米和2.08到2.23微米之间的区域
mask = ((wave_200 > 1.75) & (wave_200 < 1.97)) | ((wave_200 > 2.08) & (wave_200 < 2.23))

# 使用掩膜选择的波长和通量创建光谱对象
obs_200 = Spectrum1D(spectral_axis=wave_200[mask]*u.um, flux=flux_200[mask]*spec_unit)

# 拟合光谱的连续背景，使用7阶的切比雪夫多项式
continuum_200 = continuum.fit_generic_continuum(obs_200, model=Chebyshev1D(7))

plt.figure()  # 创建一个新的图形

# 绘制原始光谱数据的误差条，线型为实线，颜色为紫色，标签为'F200W'
plt.errorbar(wave_200, flux_200, ls='-', color='purple', label='F200W')

# 绘制拟合的连续背景，颜色为红色，标签为'Fitted continuum'
plt.plot(wave_200[mask], continuum_200(obs_200.spectral_axis), color='r', label='Fitted continuum')

plt.legend(loc=0)  # 显示图例，位置自动选择

# 计算减去连续背景后的通量
flux_200_sub = flux_200 * spec_unit - continuum_200(wave_200*u.um)

# 创建新的掩膜，选择波长在1.75到2.23微米之间的区域
mask_200 = (wave_200 > 1.75) & (wave_200 < 2.23)

# F150W

spec_unit = u.MJy  # 定义光谱单位为MJy

# 创建掩膜，选择波长在1.35到1.48微米和1.6到1.65微米之间的范围
mask = ((wave_150 > 1.35) & (wave_150 < 1.48)) | ((wave_150 > 1.6) & (wave_150 < 1.65))

# 使用掩膜选择的波长和流量创建Spectrum1D对象
obs_150 = Spectrum1D(spectral_axis=wave_150[mask]*u.um, flux=flux_150[mask]*spec_unit)

# 拟合通量的连续谱，使用7阶的切比雪夫多项式
continuum_150 = continuum.fit_generic_continuum(obs_150, model=Chebyshev1D(7))

# 从原始流量中减去拟合的连续谱，得到减去连续谱后的流量
flux_150_sub = flux_150 * spec_unit - continuum_150(wave_150 * u.um)

# 创建一个新的图形
plt.figure()

# 绘制原始流量的误差条图，线型为实线，颜色为橙色，标签为'F150W'
plt.errorbar(wave_150, flux_150, ls='-', color='orange', label='F150W')

# 绘制拟合的连续谱，颜色为红色，标签为'Fitted continuum'
plt.plot(wave_150[mask], continuum_150(obs_150.spectral_axis), color='r', label='Fitted continuum')

# 添加图例，位置为自动选择
plt.legend(loc=0)

# 创建掩膜，选择波长在1.35到1.65微米之间的范围
mask_150 = (wave_150 > 1.35) & (wave_150 < 1.65)

# F150W

spec_unit = u.MJy  # 定义光谱单位为MJy

# 创建掩膜，选择特定波长范围内的数据
mask = ((wave_115 > 1.05) & (wave_115 < 1.13)) | ((wave_115 > 1.17) & (wave_115 < 1.25))

# 使用掩膜选择波长和flux数据，创建Spectrum1D对象
obs_115 = Spectrum1D(spectral_axis=wave_115[mask]*u.um, flux=flux_115[mask]*spec_unit)

# 拟合连续谱，使用Chebyshev多项式模型
continuum_115 = continuum.fit_generic_continuum(obs_115, model=Chebyshev1D(7))

# 创建图形
plt.figure()

# 绘制F115W的光谱数据，带误差条
plt.errorbar(wave_115, flux_115, ls='-', color='lightblue', label='F115W')

# 绘制拟合的连续谱
plt.plot(wave_115[mask], continuum_115(obs_115.spectral_axis), color='r', label='Fitted continuum')

plt.ylim(0, 0.5)  # 设置y轴范围

plt.legend(loc=0)  # 显示图例

# 计算减去连续谱后的flux
flux_115_sub = flux_115 * spec_unit - continuum_115(wave_115 * u.um)

# 创建掩膜，选择特定波长范围内的数据
mask_115 = (wave_115 > 1.02) & (wave_115 < 1.25)

# 绘制所有数据

plt.figure()  # 创建一个新的图形

# 绘制F115W波段的光谱数据，带误差条
plt.errorbar(wave_115[mask_115], flux_115_sub[mask_115], ls='-', color='lightblue', label='F115W')

# 绘制F150W波段的光谱数据，带误差条
plt.errorbar(wave_150[mask_150], flux_150_sub[mask_150], ls='-', color='orange', label='F150W')

# 绘制F200W波段的光谱数据，带误差条
plt.errorbar(wave_200[mask_200], flux_200_sub[mask_200], ls='-', color='purple', label='F200W')

# 填充F115W波段的误差区域
plt.fill_between(wave_115, -flux_err_115, flux_err_115, ls='-', color='gray', label='', alpha=0.3)

# 填充F150W波段的误差区域
plt.fill_between(wave_150, -flux_err_150, flux_err_150, ls='-', color='gray', label='', alpha=0.3)

# 填充F200W波段的误差区域
plt.fill_between(wave_200, -flux_err_200, flux_err_200, ls='-', color='gray', label='', alpha=0.3)

# 设置x轴标签
plt.xlabel('Wavelength / um')

# 显示图例，位置在左上角
plt.legend(loc=2)

# 连接数组

# 使用掩码将波长数组连接起来
wave_obs = np.concatenate([wave_115[mask_115], wave_150[mask_150], wave_200[mask_200]])

# 使用掩码将流量数组连接起来
flux_obs = np.concatenate([flux_115_sub[mask_115], flux_150_sub[mask_150], flux_200_sub[mask_200]])

# 使用掩码将流量误差数组连接起来
flux_err_obs = np.concatenate([flux_err_115[mask_115], flux_err_150[mask_150], flux_err_200[mask_200]])

# 创建一个新的图形
plt.figure()

# 绘制波长与流量的关系图
plt.plot(wave_obs, flux_obs)

# 在流量误差范围内填充区域，使用灰色表示
plt.fill_between(wave_obs, -flux_err_obs, flux_err_obs, ls='-', color='gray', label='', alpha=0.3)

wht_obs = 1 / flux_err_obs**2  # 计算权重，权重为观测误差的倒数平方

spec_unit = u.MJy  # 定义光谱单位为MJy（兆焦耳每平方弧秒）

# 创建一个QTable，包含波长、光谱、权重和不确定性
dataspec = QTable([wave_obs*u.um, flux_obs*spec_unit, wht_obs, flux_err_obs],
                  names=('wavelength', 'flux', 'weight', 'uncertainty'))

dataspec_sub = dataspec[dataspec['weight'] > 0.]  # 筛选出权重大于0的观测数据

# 现在将其转换为Spectrum1D实例
p_obs = Spectrum1D(spectral_axis=dataspec_sub['wavelength'],  # 设置光谱轴为波长
                   flux=dataspec_sub['flux'],  # 设置光谱值为光谱
                   uncertainty=StdDevUncertainty(dataspec_sub['uncertainty']),  # 设置不确定性为标准偏差不确定性
                   unit='MJy')  # 设置单位为MJy

使用输入模板进行交叉相关；#

(另一个关于Specutils交叉相关的笔记本，由Ivo Busko编写： https://spacetelescope.github.io/jdat_notebooks/pages/redshift_crosscorr/redshift_crosscorr.html)

加载合成模板；#

  • 目前为了测试功能，我们使用用于模拟的输入光谱模板。

def read_hdf5(filename):
    # 定义一个函数，用于读取HDF5文件并提取数据

    contents = {}  # 创建一个空字典，用于存储提取的数据

    with h5py.File(filename, 'r') as file_obj:  # 以只读模式打开HDF5文件
        for key in file_obj.keys():  # 遍历文件中的所有键
            dataset = file_obj[key]  # 获取当前键对应的数据集

            try:
                wave_units_string = dataset.attrs['wavelength_units']  # 尝试获取波长单位属性
            except KeyError:
                wave_units_string = 'micron'  # 如果没有找到，默认设置为'micron'

            # 捕获常见错误
            if wave_units_string.lower() in ['microns', 'angstroms', 'nanometers']:  # 检查单位是否在常见单位中
                wave_units_string = wave_units_string[0:-1]  # 去掉单位的复数形式

            # 获取数据
            waves = dataset[0]  # 提取波长数据
            fluxes = dataset[1]  # 提取通量数据

            contents[int(key)] = {'wavelengths': waves, 'fluxes': fluxes}  # 将提取的数据存入字典

    return contents  # 返回包含所有数据的字典

# 下载文件，如果尚不存在。

if not os.path.exists('./pipeline_products'):  # 检查'./pipeline_products'目录是否存在

    import zipfile  # 导入zipfile模块，用于处理zip文件
    import urllib.request  # 导入urllib.request模块，用于下载文件

    # 定义zip文件的下载链接
    boxlink = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/NIRISS_lensing_cluster/pipeline_products.zip'

    boxfile = './pipeline_products.zip'  # 定义下载后保存的zip文件路径

    # 下载zip文件
    urllib.request.urlretrieve(boxlink, boxfile)

    zf = zipfile.ZipFile(boxfile, 'r')  # 打开下载的zip文件

    list_names = zf.namelist()  # 获取zip文件中所有文件的名称列表

    for item in list_names:  # 遍历zip文件中的每个文件
        # 使用extract提取单个文件，而不是extractall，因为在文件具有绝对（/）或相对（..）路径时，extract更安全
        zf.extract(member=item)

# 模板文件路径
file_temp = './pipeline_products/source_sed_file_from_sources_extend_01_and_sed_file.hdf5'

# 读取HDF5文件的内容
content = read_hdf5(file_temp)

# 输出读取的内容
content

# 模板已经为模拟进行了红移处理。

# 现在将其蓝移到 z=0;

z_input = 2.1  # 输入的红移值

iix = 10709  # 特定模板的索引

flux_all = content[iix]['fluxes']  # 获取特定模板的光谱通量

wave_all = content[iix]['wavelengths'] / (1. + z_input) * 1e4  # 计算蓝移后的波长并转换为合适的单位

plt.plot(wave_all, flux_all, color='orange')  # 绘制蓝移后的光谱

plt.xlim(2000, 8000)  # 设置x轴的范围

# 切割模板以获得更好的拟合

con_tmp = (wave_all > 3600) & (wave_all < 13000)  # 选择波长在3600到13000之间的索引

flux = flux_all[con_tmp]  # 根据选择的索引提取flux数据

wave = wave_all[con_tmp]  # 根据选择的索引提取wave数据

flux /= flux.max()  # 将flux归一化到最大值为1

wave *= 1.e-4  # 将波长转换为微米（um），with_spectral_unit不进行转换

factor = p_obs.flux.unit  # 获取观测数据的flux单位以进行归一化

template = Spectrum1D(spectral_axis=wave*u.um, flux=flux*factor)  # 创建一个Spectrum1D对象，包含波长和flux

plt.figure()  # 创建一个新的图形

plt.plot(template.spectral_axis, template.flux, color='orange')  # 绘制模板的光谱轴和flux，颜色为橙色

plt.xlabel(template.spectral_axis.unit)  # 设置x轴标签为波长单位

将模板分箱到粗数组中，如果您有太多元素并导致下面的问题；#

-> 在这个例子中，这是不必要的。#

# 如果你想要进行分箱处理，将此设置为True;

if False:  # 检查是否启用分箱处理

    wave_bin = np.arange(wave.min(), wave.max(), 20.0)  # 创建波长分箱数组，从最小波长到最大波长，步长为20.0

    flux_bin = np.interp(wave_bin, wave, flux)  # 使用插值方法计算分箱后的光谱数据

else:  # 如果不进行分箱处理

    wave_bin = wave  # 直接使用原始波长数据

    flux_bin = flux  # 直接使用原始光谱数据

factor = p_obs.flux.unit  # 获取观测数据的单位，以便将模板归一化到合理范围

template = Spectrum1D(spectral_axis=wave_bin*u.um, flux=flux_bin*factor)  # 创建一个Spectrum1D对象，包含波长和光谱数据

plt.figure()  # 创建一个新的图形

plt.plot(template.spectral_axis, template.flux, color='orange')  # 绘制光谱数据，颜色为橙色

plt.xlabel(template.spectral_axis.unit)  # 设置x轴标签为波长单位

# 对拟合温度进行连续谱拟合

# 定义排除的光谱区域
regions = [SpectralRegion(0.37*u.um, 0.452*u.um),  # 第一段光谱区域
           SpectralRegion(0.476*u.um, 0.53*u.um),  # 第二段光谱区域
           SpectralRegion(0.645*u.um, 0.67*u.um)]   # 第三段光谱区域

# 使用Chebyshev多项式模型拟合连续谱，排除指定的光谱区域
continuum_model = continuum.fit_generic_continuum(template, model=Chebyshev1D(3), exclude_regions=regions)

# 从模板中减去拟合的连续谱，得到去除连续谱后的模板
p_template = template - continuum_model(template.spectral_axis)

# 创建一个新的图形
plt.figure()

# 绘制去除连续谱后的模板光谱
plt.plot(p_template.spectral_axis, p_template.flux, color='orange')

# 设置x轴标签为光谱轴的单位
plt.xlabel(p_template.spectral_axis.unit)

# 设置图形标题
plt.title('Continuum subtracted template')

# 将数据数组转换为Specutils格式

sflux = p_template.flux  # 获取模板的光谱通量

# 创建一个Spectrum1D对象，包含光谱轴和归一化后的光谱通量
p_smooth_template = Spectrum1D(spectral_axis=p_template.spectral_axis,  # 设置光谱轴
                               flux=sflux/sflux.max())  # 归一化光谱通量

plt.figure()  # 创建一个新的图形

plt.plot(p_smooth_template.spectral_axis, p_smooth_template.flux, color='b', label='Smoothed template')  # 绘制平滑模板的光谱，颜色为蓝色

plt.plot(p_obs.spectral_axis, p_obs.flux, color='red', label='Spectrum')  # 绘制观测光谱，颜色为红色

plt.xlabel(p_smooth_template.spectral_axis.unit)  # 设置x轴标签为光谱轴的单位

plt.legend(loc=0)  # 显示图例，位置自动选择

注意：上述两个输入光谱必须有重叠部分，即使只有一小部分，以确保Specutils正常工作。#

使用specutils功能进行交叉相关；#

# 在没有额外规格的情况下，整个模板和整个光谱

# 将被包含在相关性计算中。这通常会导致

# 执行时间显著增加。建议将模板裁剪

# 以仅在有用区域内工作。

# 进行模板相关性计算，返回相关系数和滞后值
corr, lag = correlation.template_correlate(p_obs, p_smooth_template)

plt.figure()  # 创建一个新的图形窗口

plt.gcf().set_size_inches((8., 4.))  # 设置图形的尺寸为8x4英寸

plt.plot(lag, corr, linewidth=0.5)  # 绘制lag与corr的关系图，线宽为0.5

plt.xlabel(lag.unit)  # 设置x轴标签为lag的单位

# 基于最大值计算红移

index_peak = np.argmax(corr)  # 找到相关性数组中的最大值索引

v = lag[index_peak]  # 根据最大值索引获取对应的速度

z = v / const.c.to('km/s')  # 计算红移，使用光速进行单位转换

print("Peak maximum at: ", v)  # 输出峰值最大值
print("Redshift from peak maximum: ", z)  # 输出从峰值最大值计算得到的红移

# 基于抛物线拟合的红移计算

n = 8  # 定义相关最大值左右各取8个点

peak_lags = lag[index_peak-n:index_peak+n+1].value  # 获取最大值附近的延迟值

peak_vals = corr[index_peak-n:index_peak+n+1].value  # 获取最大值附近的相关值

p = np.polyfit(peak_lags, peak_vals, deg=2)  # 对延迟值和相关值进行二次多项式拟合

roots = np.roots(p)  # 计算拟合多项式的根

v_fit = np.mean(roots) * u.km/u.s  # 计算最大值对应的速度，取根的平均值

z = v_fit / const.c.to('km/s')  # 计算红移，速度与光速的比值

print("")  # 打印空行

print("Parabolic fit with maximum at: ", v_fit)  # 输出抛物线拟合的最大值

print("Redshift from parabolic fit: ", z)  # 输出通过抛物线拟合得到的红移

它看起来怎么样？#

print('z =', z)  # 打印红移值 z

plt.figure()  # 创建一个新的图形

plt.plot(  # 绘制模板光谱
    template.spectral_axis * (1. + z),  # 将模板的光谱轴乘以 (1 + z) 以进行红移校正
    p_smooth_template.flux / np.max(p_smooth_template.flux),  # 归一化模板的光谱通量
    color='b', label='Template'  # 设置颜色为蓝色，并标记为 'Template'
)

plt.plot(p_obs.spectral_axis, p_obs.flux / np.max(p_obs.flux), label='Observed')  # 绘制观测光谱并归一化

plt.legend(loc=2)  # 显示图例，位置在左上角

plt.xlabel(p_obs.spectral_axis.unit)  # 设置 x 轴标签为观测光谱的单位
Contents