BOTS 时间序列观测#
用例： 明亮天体时间序列；提取系外行星光谱。
数据： 来自地面观测活动的 JWST 模拟 NIRSpec 数据；GJ436b 光谱来自 Goyal 等人 (2018)。
工具： scikit, lmfit, scipy, matplotlib, astropy, pandas。
跨仪器： .
文档： 本笔记本是 STScI 更大 后处理数据分析工具生态系统 的一部分。
作者： David K. Sing (dsing@jhu.edu)
最后更新： 2020年7月2日
引言#
本笔记本使用在地面观测活动中获取的时间序列 JWST NIRSpec 数据，演示如何从时间序列观测中提取系外行星光谱。
这些数据源自 ISIM-CV3，即 JWST 综合科学仪器模块 (ISIM) 的低温真空测试活动，该活动于 2015-2016 年冬季在戈达德太空飞行中心进行（Kimble 等人，2016）。数据可以在 https://www.cosmos.esa.int/web/jwst-nirspec/test-data 找到，G. Giardino, S. Birkmann, P. Ferruit, B. Dorner, B. Rauscher 对数据的详细和深刻的报告可以在此处找到：ftp://ftp.cosmos.esa.int/jwstlib/ReleasedCV3dataTimeSeries/CV3_TimeSeries_PRM.tgz
该 NIRSpec 时间序列数据集在像素级别注入了过境光曲线，紧密模拟了过境系外行星的明亮天体时间序列 (BOTS) 观测。在这种情况下，选择了来自
的系外行星网格中的 GJ436b 光谱（太阳金属丰度下的清晰大气）。使用实际的 NIRSpec 数据集，可以更准确地模拟探测器的噪声特性、抖动以及从时间序列观测中提取系外行星光谱的影响。
总体而言，本笔记本的目的是处理这些时间序列观测，以：
从 2D 光谱图像中提取 1D 光谱。
定义一个时间序列模型，以拟合波长依赖的过境光曲线。
拟合 1D 光谱的每个时间序列波长区间，测量所需的量 \(R_{pl}(\lambda)/R_{star}\)。
生成可与模型进行比较的测量透射光谱。
该示例输出每个光谱区间的拟合光曲线及其拟合统计信息。
加载包#
本笔记本使用的包（matplotlib、astropy、scipy、glob、lmfit、pickle、os、sklearn）均可通过pip以标准方式安装。
从ExoTiC-ISM中提取了几种计算边缘变暗（limb-darkening）和过境模型的例程（__ Laginja & Wakeford 2020__；https://github.com/hrwakeford/ExoTiC-ISM），并进行了稍微的调整。用于边缘变暗计算的完整恒星模型集也可以从ExoTiC-ISM下载，因为本笔记本仅下载并加载了用于生成边缘变暗的单个恒星模型。
import numpy as np # 导入NumPy库，用于数值计算
import matplotlib.pyplot as plt # 导入Matplotlib库，用于绘图
import matplotlib.ticker as ticker # 导入Matplotlib的ticker模块，用于控制坐标轴刻度
from matplotlib.backends.backend_pdf import PdfPages # 导入PdfPages，用于将图形保存为PDF文件
from astropy.utils.data import download_file # 从Astropy库导入下载文件的功能
from astropy.table import Table # 从Astropy库导入Table类，用于处理表格数据
from astropy.io import fits, ascii # 导入FITS和ASCII文件的读写功能
from astropy.modeling.models import custom_model # 导入自定义模型功能
from astropy.modeling.fitting import LevMarLSQFitter # 导入Levenberg-Marquardt最小二乘拟合器
from scipy.interpolate import interp1d, splev, splrep # 导入SciPy库中的插值功能
from scipy.io import readsav # 导入读取SAV文件的功能
from scipy import stats # 导入SciPy库中的统计功能
import glob # 导入glob模块，用于文件路径匹配
import lmfit # 导入lmfit库，用于非线性最小二乘拟合
import pickle # 导入pickle模块，用于对象序列化和反序列化
from os import path, mkdir # 从os模块导入路径操作和创建目录的功能
from sklearn.linear_model import LinearRegression # 导入线性回归模型
import pandas as pd # 导入Pandas库，用于数据处理和分析
import os # 导入os模块，用于操作系统相关功能
import shutil # 导入shutil模块，用于文件和目录的高级操作
设置参数#
拟合的参数包括存放数据和边缘变暗（limb darkening）恒星模型的目录，以及行星和恒星的属性。这里输入的恒星和行星值（以GJ436为模型）与用于模拟注入过境的值相同。请注意，用于注入过境的4500K恒星模型比GJ436A更热。
# SETUP ----------------------------------------------
# 设置目录
save_directory = './notebookrun2/' # 本地目录，用于保存文件
data_directory = './' # 本地数据目录，用于处理fits文件（如果需要）
# 设置探测器属性和红噪声测量时间尺度
gain = 1.0 # 2D光谱已转换为计数，探测器增益为1.0
binmeasure = 256 # 测量红噪声的分箱技术，选择评估sigma_r的箱大小
number_of_images = 8192 # 数据集中图像的数量
# 设置行星属性
grating = 'NIRSpecPrism' # 使用的光栅类型
ld_model = '3D' # 选择3D/1D恒星模型（过境是用3D模型注入的）
# 设置恒星属性以进行光晕变暗计算
Teff = 4500 # 有效温度（K）
logg = 4.5 # 表面重力
M_H = 0.0 # 恒星金属丰度 log_10[M/H]
Rstar = 0.455 # 行星半径（以太阳半径为单位）
# 设置过境参数（可以从NASA系外行星档案中获取）
t0 = 2454865.084034 # 最低交点的BJD时间
per = 2.64389803 # 轨道周期（天）BJD_TDB
rp = 0.0804 # 行星半径（以恒星半径为单位）
a_Rs = 14.54 # 半长轴（输入a/Rstar，因此单位为恒星半径）
inc = 86.858 * (2*np.pi/360) # 轨道倾角（从度转换为弧度）
ecc = 0.0 # 离心率
omega = 0.0 * (2*np.pi/360) # 近日点经度（从度转换为弧度）
# 计算恒星密度（g/cm^3），基于a/Rs
rho_star = (3*np.pi)/(6.67259E-8*(per*86400)**2)*(a_Rs)**3 # 从a/Rs计算恒星密度
# a_Rs=(rho_star*6.67259E-8*per_sec*per_sec/(3*np.pi))**(1/3) # 从恒星密度计算a/Rs（g/cm^3）
# 创建本地目录
if not path.exists(save_directory): # 检查保存目录是否存在
mkdir(save_directory) # 如果不存在，则创建一个新的目录以保存输出
if not path.exists(save_directory+'3DGrid'): # 检查3DGrid目录是否存在
mkdir(save_directory+'3DGrid') # 如果不存在，则创建新的3DGrid目录以保存数据
limb_dark_directory = save_directory # 指向包含3DGrid/目录的光晕暗化目录
下载和加载NIRSpec数据#
FITS图像被加载，包含图像日期和科学光谱的信息被保存。
一个默认的通量偏移值BZERO也从头文件中提取，并从每个科学帧中减去。
读取2^13个FITS文件的速度较慢。为了加快速度，我们创建了一个pickle文件，用于首次加载FITS图像时。这个1GB的pickle文件会被加载，而不是读取FITS文件（如果找到的话）。
另外，FITS文件可以在这里下载：
https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/transit_spectroscopy_notebook/Archive.Trace_SLIT_A_1600_SRAD-PRM-PS-6007102143_37803_JLAB88_injected.tar.gz。图像位于一个tar.gz归档文件中，需要解压缩，并在上面的SETUP单元中设置data_directory变量为该目录。
下面的单元将下载1GB的JWST数据pickle文件，以及其他几个所需的文件。
# 下载1GB NIRSpec数据
fn_jw = save_directory + 'jwst_data.pickle' # 定义JWST数据文件的路径
if not path.exists(fn_jw): # 检查文件是否已存在
fn = download_file('https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/transit_spectroscopy_notebook/jwst_data.pickle') # 下载JWST数据文件
dest = shutil.move(fn, save_directory + 'jwst_data.pickle') # 移动文件到指定的保存目录
print('JWST Data Download Complete') # 打印下载完成的消息
# 下载进一步需要的文件，并移动到本地目录以便于重复访问
fn_sens = download_file('https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/transit_spectroscopy_notebook/NIRSpec.prism.sensitivity.sav') # 下载灵敏度文件
dest = shutil.move(fn_sens, save_directory + 'NIRSpec.prism.sensitivity.sav') # 移动灵敏度文件到保存目录
fn_ld = download_file('https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/transit_spectroscopy_notebook/3DGrid/mmu_t45g45m00v05.flx') # 下载3D网格文件
destld = shutil.move(fn_ld, save_directory + '3DGrid/mmu_t45g45m00v05.flx') # 移动3D网格文件到保存目录
加载Pickle文件数据。或者，数据也可以从原始的FITS文件中读取。
if path.exists(fn_jw): # 检查JWST数据文件是否存在
dbfile = open(fn_jw, 'rb') # 以二进制模式打开文件进行读取
jwst_data = pickle.load(dbfile) # 从pickle文件加载JWST数据
print('Loading JWST data from Pickle File') # 打印加载数据的提示信息
bjd = jwst_data['bjd'] # 提取BJD时间数据
wsdata_all = jwst_data['wsdata_all'] # 提取波长数据
shx = jwst_data['shx'] # 提取X方向的偏移数据
shy = jwst_data['shy'] # 提取Y方向的偏移数据
common_mode = jwst_data['common_mode'] # 提取公共模式数据
all_spec = jwst_data['all_spec'] # 提取所有光谱数据
exposure_length = jwst_data['exposure_length'] # 提取曝光时间数据
dbfile.close() # 关闭文件
print('Done') # 打印完成的提示信息
elif not path.exists(fn_jw): # 如果JWST数据文件不存在
# 加载所有的fits图像
# 为BJD时间和白光曲线总计数创建数组
list = glob.glob(data_directory + "*.fits") # 获取所有fits文件的列表
index_of_images = np.arange(number_of_images) # 创建图像索引数组
bjd = np.zeros((number_of_images)) # 初始化BJD数组
exposure_length = np.zeros((number_of_images)) # 初始化曝光时间数组
all_spec = np.zeros((32, 512, number_of_images)) # 初始化光谱数组
for i in index_of_images: # 遍历每个图像索引
img = list[i] # 获取当前图像文件名
print(img) # 打印当前图像文件名
hdul = fits.open(img) # 打开当前fits图像
# hdul.info() # 可选：打印图像信息
bjd_image = hdul[0].header['BJD_TDB'] # 从头文件中提取BJD时间
BZERO = hdul[0].header['BZERO'] # 提取flux值偏移量
bjd[i] = bjd_image # 将BJD时间存入数组
expleng = hdul[0].header['INTTIME'] # 提取单次MULTIACCUM的总曝光时间（秒）
exposure_length[i] = expleng / 86400. # 将曝光时间转换为天并存入数组
print(bjd[i]) # 打印当前BJD时间
data = hdul[0].data # 提取图像数据
# total counts in image
# total_counts[i]=gain*np.sum(data[11:18,170:200]-BZERO) # 计算在特定区域的总计数（注释掉）
all_spec[:, :, i] = gain * (data - BZERO) # 将所有光谱数据加载到数组中并减去flux值偏移量
hdul.close() # 关闭当前图像文件
# 对数据进行排序
srt = np.argsort(bjd) # 获取排序索引
bjd = bjd[srt] # 根据索引排序BJD时间
# total_counts=total_counts[srt] # （注释掉）根据索引排序总计数
exposure_length = exposure_length[srt] # 根据索引排序曝光时间
all_spec[:, :, :] = all_spec[:, :, srt] # 根据索引排序所有光谱数据
# 获取数据的波长
file_wave = download_file('https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/transit_spectroscopy_notebook/JWST_NIRSpec_wavelength_microns.txt') # 下载波长数据文件
f = open(file_wave, 'r') # 打开波长数据文件
wsdata_all = np.genfromtxt(f) # 从文件中读取波长数据
print('wsdata size :', wsdata_all.shape) # 打印波长数据的大小
print('Data wavelength Loaded :', wsdata_all) # 打印加载的波长数据
print('wsdata new size :', wsdata_all.shape) # 打印新的波长数据大小
# 读取去趋势参数
# 参数的均值必须为0.0以便正确归一化
# 理想情况下参数的标准差应为1.0
file_xy = download_file('https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/transit_spectroscopy_notebook/JWST_NIRSpec_Xposs_Yposs_CM_detrending.txt') # 下载去趋势参数文件
f = open(file_xy, 'r') # 打开去趋势参数文件
data = np.genfromtxt(f, delimiter=',') # 从文件中读取去趋势参数数据
shx = data[:, 0] # 提取X方向偏移数据
shy = data[:, 1] # 提取Y方向偏移数据
common_mode = data[:, 2] # 提取公共模式数据
# 将数据存储到pickle文件中
jwst_data = {'bjd': bjd, 'wsdata_all': wsdata_all, 'shx': shx, 'shy': shy, 'common_mode': common_mode, 'all_spec': all_spec, 'exposure_length': exposure_length} # 创建数据字典
dbfile = open('jwst_data.pickle', 'ab') # 以二进制模式打开pickle文件进行追加
pickle.dump(jwst_data, dbfile) # 将数据字典写入pickle文件
dbfile.close() # 关闭pickle文件
可视化二维光谱数据#
expnum = 2 # 选择要查看的曝光编号
plt.rcParams['figure.figsize'] = [10.0, 3.0] # 设置图形尺寸
plt.rcParams['figure.dpi'] = 200 # 设置分辨率
plt.rcParams['savefig.dpi'] = 200 # 设置保存图像的分辨率
plt.rcParams['image.aspect'] = 5 # 设置长宽比（CCD非常长！）
plt.cmap = plt.cm.magma # 使用magma色图
plt.cmap.set_bad('k', 1.) # 设置无效值的颜色为黑色
plt.rcParams['image.cmap'] = 'magma' # 设置图像的色图为magma
plt.rcParams['image.interpolation'] = 'none' # 设置图像插值方式为无
plt.rcParams['image.origin'] = 'lower' # 设置图像原点在下方
plt.rcParams['font.family'] = "monospace" # 设置字体为等宽字体
plt.rcParams['font.monospace'] = 'DejaVu Sans Mono' # 设置等宽字体为DejaVu Sans Mono
img = all_spec[:, :, expnum] # 获取指定曝光编号的图像数据
zeros = np.where(img <= 0) # 找到图像中小于等于0的值
img[zeros] = 1E-10 # 将小于等于0的值设为一个小数以避免对数计算错误
fig, axs = plt.subplots() # 创建子图
f = axs.imshow(np.log10(img), vmin=0) # 绘制图像，使用对数尺度
plt.xlabel('x-pixel') # 设置x轴标签
plt.ylabel('y-pixel') # 设置y轴标签
axs.yaxis.set_major_locator(ticker.MultipleLocator(5)) # 设置y轴主刻度间隔为5
axs.yaxis.set_minor_locator(ticker.MultipleLocator(1)) # 设置y轴次刻度间隔为1
axs.xaxis.set_major_locator(ticker.MultipleLocator(50)) # 设置x轴主刻度间隔为50
axs.xaxis.set_minor_locator(ticker.MultipleLocator(10)) # 设置x轴次刻度间隔为10
plt.title('2D NIRSpec Image of Exposure ' + str(expnum)) # 设置图像标题
fig.colorbar(f, label='Log$_{10}$ Electron counts', ax=axs) # 添加颜色条，标记为对数电子计数
plt.show() # 显示图像
从2D图像数组中提取1D光谱#
理想情况下，从2D图像中提取1D光谱应使用沿拟合轨迹的最佳光圈提取，使用与IRAF/apall等效的例程。此功能在astro-py中尚不可用。
已经应用了几个处理步骤。这里的2D光谱已经进行了平场校正，并且通过从每列未照明像素中减去中位数计数率，去除了每个像素的1/f噪声（有关1/f噪声的更多信息，请参见 Giardino et al.）。每个2D图像在X和Y方向上也进行了对齐，以确保每个像素对应相同的波长。由于CV3测试对光通量稳定性没有要求，因此LED的~1%光通量变化也已被去除。
对于光谱提取，此处的示例简单地使用了一个简单的求和框。8192个2D光谱已预加载到一个numpy数组中。光谱在像素Y=16处达到峰值。对于每一列，光圈总和是在Y轴像素11到18之间进行的，这包含了大部分光谱计数。更宽的光圈会增加更多计数，但也会引入更多噪声。
此处未进行进一步清理步骤
理想情况下，因各种原因被标记为坏的像素应进行清理。
应识别并去除宇宙射线。
all_spec.shape # 查看 all_spec 数组的形状
y_lower = 11 # 下提取孔径的像素行数
y_upper = 18 # 上提取孔径的像素行数
all_spec_1D = np.sum(all_spec[y_lower:y_upper, :, :], axis=0) # 在 Y 轴上对像素 11 到 18 进行求和
# 绘图设置
plt.rcParams['figure.figsize'] = [10.0, 3.0] # 设置图形的尺寸
plt.rcParams['figure.dpi'] = 200 # 设置分辨率
plt.rcParams['savefig.dpi'] = 200 # 设置保存图形的分辨率
plt.rcParams['image.aspect'] = 5 # 设置图像的长宽比（CCD 很长！）
plt.cmap = plt.cm.magma # 设置颜色映射为 magma
plt.cmap.set_bad('k', 1.) # 设置无效值的颜色为黑色，透明度为 1
plt.rcParams['image.cmap'] = 'magma' # 设置图像的颜色映射
plt.rcParams['image.interpolation'] = 'none' # 设置图像插值方式为无
plt.rcParams['image.origin'] = 'lower' # 设置图像原点在下方
plt.rcParams['font.family'] = "monospace" # 设置字体为等宽字体
plt.rcParams['font.monospace'] = 'DejaVu Sans Mono' # 设置等宽字体的具体类型
img = all_spec[:, :, expnum] # 获取特定曝光编号的图像数据
zeros = np.where(img <= 0) # 找到图像中小于等于零的值的索引
img[zeros] = 1E-10 # 将零或负值设置为一个小数，以便在对数尺度上绘图
fig, axs = plt.subplots() # 创建一个图形和一组子图
f = axs.imshow(np.log10(img), vmin=0) # 绘制图像，使用对数尺度
plt.xlabel('x-pixel') # 设置 x 轴标签
plt.ylabel('y-pixel') # 设置 y 轴标签
axs.yaxis.set_major_locator(ticker.MultipleLocator(5)) # 设置 y 轴主刻度间隔为 5
axs.yaxis.set_minor_locator(ticker.MultipleLocator(1)) # 设置 y 轴次刻度间隔为 1
axs.xaxis.set_major_locator(ticker.MultipleLocator(50)) # 设置 x 轴主刻度间隔为 50
axs.xaxis.set_minor_locator(ticker.MultipleLocator(10)) # 设置 x 轴次刻度间隔为 10
plt.axhline(y_lower, color='w', ls='dashed') # 绘制下提取孔径的虚线
plt.axhline(y_upper, color='w', ls='dashed') # 绘制上提取孔径的虚线
plt.title('2D NIRSpec Image of Exposure ' + str(expnum)) # 设置图形标题
fig.colorbar(f, label='Log$_{10}$ Electron counts', ax=axs) # 添加颜色条，标注为对数电子计数
plt.show() # 显示图形
可视化一维光谱数据#
fig, axs = plt.subplots() # 创建一个子图
f = plt.plot(wsdata_all, all_spec_1D[:, 0], linewidth=2, zorder=0) # 在数据上绘制过渡模型
plt.xlabel(r'Wavelength ($\mu$m)') # 设置x轴标签为波长
plt.ylabel('Flux (e-)') # 设置y轴标签为通量
axs.xaxis.set_major_locator(ticker.MultipleLocator(0.5)) # 设置x轴主刻度间隔为0.5
axs.xaxis.set_minor_locator(ticker.MultipleLocator(0.1)) # 设置x轴次刻度间隔为0.1
plt.annotate('H$_2$O', xy=(3.0, 42000)) # 在图中添加'H2O'注释，位置为(3.0, 42000)
plt.annotate('CO$_2$', xy=(4.2, 42000)) # 在图中添加'CO2'注释，位置为(4.2, 42000)
plt.show() # 显示图形
CV3测试观察了一盏灯，其点扩散函数（PSF）与JWST将拥有的相似，并且在大约1.5到4.5 \(\mu\)m的范围内有显著的计数。
低温测试舱的窗口上积累了CO\(_2\)和H\(_2\)O冰，这可以在二维光谱中看到作为光谱吸收特征。
计算轨道相位及用于绘图目的的单独细网模型#
# 计算轨道相位
phase = (bjd-t0) / (per) # 计算相位，单位为天，相对于T0历元
phase = phase - np.fix(phase[number_of_images-1]) # 将当前相位调整为0.0
t_fine = np.linspace(np.min(bjd), np.max(bjd), 1000) # 生成用于计算光曲线的时间点，共1000个点
phase_fine = (t_fine-t0)/(per) # 计算相位，单位为天，相对于T0历元
phase_fine = phase_fine - np.fix(phase[number_of_images-1]) # 将当前相位调整为0.0
# 计算b0，表示光度变化的幅度
b0 = a_Rs * np.sqrt((np.sin(phase * 2 * np.pi)) ** 2 + (np.cos(inc) * np.cos(phase * 2 * np.pi)) ** 2)
intransit = (b0 - rp < 1.0E0).nonzero() # 选择在首次接触和第四次接触之间的索引
outtransit = (b0 - rp > 1.0E0).nonzero() # 选择在过境之外的索引
处理探测器上的系统漂移#
CV3测试通过引入较大的空间抖动和漂移来评估仪器的稳定性。这导致光谱在2D探测器上的X,Y位置发生显著移动。虽然这种大规模的位移已经被移除，从而对齐了光谱，但像素内和像素间的灵敏度引入了通量变化，这些变化需要被去除。CV3测试中的抖动超过了30毫弧秒（mas），这大约是JWST稳定性要求的4倍。因此，在轨道上，这些探测器效应预计会显著减小，但仍然会存在，并需要在时间序列观测中进行建模和去除。
这里的探测器X,Y位置是通过对2D图像进行交叉相关测量得到的（首先沿一个维度压缩光谱），并保存在数组\(shx\)和\(shy\)中。这些去趋势向量理想情况下应该使用每次积分的光谱提取中的轨迹位置值进行测量，因为这也可以准确测量光谱在探测器上如何在每次积分中发生空间变化。
探测器的位移原始幅度接近0.2像素，尽管这些向量已经进行了初步归一化。为了去趋势，这些数组的均值应该为0，标准差为1.0。
在CV3数据中，还可以看到与LED灯相关的残余颜色依赖趋势，这可以通过缩放原始共模灯趋势部分去除，该趋势是使用CV3白光曲线测量的。
shx_tmp = shx / np.mean(shx) - 1.0E0 # 将shx的均值调整为0.0
shx_detrend = shx_tmp / np.std(shx_tmp) # 将shx的标准差调整为1.0
shy_tmp = shy / np.mean(shy) - 1.0E0 # 将shy的均值调整为0.0
shy_detrend = shy_tmp / np.std(shy_tmp) # 将shy的标准差调整为1.0
cm = common_mode / np.mean(common_mode) - 1.0E0 # 将common_mode的均值调整为0.0
cm_detrend = cm / np.std(cm) # 将common_mode的标准差调整为1.0
fig, axs = plt.subplots() # 创建一个子图
plt.plot(shx_detrend, label='X-possition') # 绘制X位置的去趋势数据
plt.plot(shy_detrend, label='Y-possition') # 绘制Y位置的去趋势数据
plt.xlabel('Image Sequence Number') # 设置X轴标签
plt.ylabel('Relative Detector Possition') # 设置Y轴标签
plt.title('Time-series Detrending Vectors') # 设置图表标题
axs.xaxis.set_major_locator(ticker.MultipleLocator(1000)) # 设置X轴主刻度间隔为1000
axs.xaxis.set_minor_locator(ticker.MultipleLocator(100)) # 设置X轴次刻度间隔为100
axs.yaxis.set_major_locator(ticker.MultipleLocator(0.5)) # 设置Y轴主刻度间隔为0.5
axs.yaxis.set_minor_locator(ticker.MultipleLocator(0.1)) # 设置Y轴次刻度间隔为0.1
plt.legend() # 显示图例
plt.show() # 显示图表
创建用于去趋势的向量数组。#
来自Sing等人（2019年）的研究：系统误差通常通过参数化的确定性模型来去除，其中非过境光度趋势与多个外部参数（或光学状态向量，\(X\)）相关。这些参数描述了在观测过程中，仪器或其他外部因素随时间的变化，并为每个光学状态参数拟合一个系数\(p_n\)，以建模并去除（或去趋势）光度光曲线。
在包含系统趋势时，随时间变化的通用参数化模型的通量测量\(f(t)\)可以建模为理论过境模型\(T(t,\theta)\)（依赖于过境参数\(\theta\)）、从恒星检测到的总基线通量\(F_0\)以及系统误差模型\(S(x)\)的组合，给出如下公式：
\[f(t) = T(t,\theta) \times F_0 \times S(x)\]
我们将使用线性模型来处理仪器的系统效应。
\[S(x) = p_1 x + p_2 y + p_3 x^2 + p_4 y^2 + p_5 x y + p_6 cm + p_7 \phi\]
其中\(cm\)是共模趋势，而\(\phi\)是线性时间趋势，有助于去除H\(_2\)O光谱特征中变化的H\(_2\)O冰。
shx = shx_detrend # 将去趋势后的shx赋值给shx
shy = shy_detrend # 将去趋势后的shy赋值给shy
common_mode = cm_detrend # 将去趋势后的公共模式赋值给common_mode
# 创建一个去趋势数组，不包含线性时间趋势
XX = np.array([shx, shy, shx**2, shy**2, shx*shy, common_mode, np.ones(number_of_images)])
XX = np.transpose(XX) # 转置数组XX
# 创建一个去趋势数组，包含线性时间趋势
XXX = np.array([shx, shy, shx**2, shy**2, shx*shy, common_mode, phase, np.ones(number_of_images)])
XXX = np.transpose(XXX) # 转置数组XXX
线性回归可以用来快速确定参数 \(p_n\)，使用的是非过境数据。
在这里，我们选择数据的一个波长区间（像素 170 到 200）来制作时间序列。选择非过境点，并对 \(S(x)\) 进行线性回归，以确定光学状态参数 \(p_n\)。
pix1 = 170 # 波长区间下限
pix2 = 200 # 波长区间上限
y = np.sum(all_spec_1D[pix1:pix2, :], axis=0) # 在选定波长区间内的总通量
msize = plt.rcParams['lines.markersize'] ** 2. # 默认标记大小
plt.rcParams['figure.figsize'] = [10.0, 3.0] # 图形尺寸
fig, axs = plt.subplots() # 创建子图
f = plt.plot(wsdata_all, all_spec_1D[:, 0], linewidth=2, zorder=0) # 绘制波长区间的光谱
plt.fill_between(wsdata_all[pix1:pix2], 0, all_spec_1D[pix1:pix2, 0], alpha=0.5) # 填充波长区间下方的区域
plt.xlabel(r'Wavelength ($\mu$m)') # 设置x轴标签
plt.ylabel('Flux (e-)') # 设置y轴标签
plt.title('1D Extracted Spectrum') # 设置图形标题
axs.xaxis.set_major_locator(ticker.MultipleLocator(0.5)) # 设置x轴主刻度
axs.xaxis.set_minor_locator(ticker.MultipleLocator(0.1)) # 设置x轴次刻度
plt.annotate('H$_2$O', xy=(3.0, 42000)) # 注释H2O的位置
plt.annotate('CO$_2$', xy=(4.2, 42000)) # 注释CO2的位置
plt.show() # 显示图形
fig, axs = plt.subplots() # 创建新的子图
plt.scatter(bjd, y/np.mean(y[outtransit]), label='$f(t)$ Data', zorder=1, s=msize*0.75, linewidth=1, alpha=0.4, marker='+', edgecolors='blue') # 绘制散点图
plt.xlabel('Barycentric Julian Date (days)') # 设置x轴标签
plt.ylabel('Relative Flux') # 设置y轴标签
plt.title(r'Time-series Transit Light Curve $\lambda=$['+str(wsdata_all[pix1])+':'+str(wsdata_all[pix2]) + r'] $\mu$m') # 设置图形标题
plt.legend() # 显示图例
plt.show() # 显示图形
regressor = LinearRegression() # 创建线性回归模型
regressor.fit(XX[outtransit], y[outtransit]/np.mean(y[outtransit])) # 拟合模型
print('Linear Regression Coefficients:') # 打印线性回归系数的提示
print(regressor.coef_) # 打印线性回归系数
系数的数量级约为 ~10\(^{-4}\)，因此趋势的幅度在 100 的 ppm 级别。
可视化拟合结果
yfit = regressor.predict(XX) # 在整个时间序列上进行拟合预测
plt.rcParams['figure.figsize'] = [10.0, 7.0] # 设置图形的尺寸
msize = plt.rcParams['lines.markersize'] ** 2. # 获取默认的标记大小并平方
plt.scatter(bjd, y/np.mean(y[outtransit]), label='$f(t)$ Data', zorder=1, s=msize*0.75, linewidth=1, alpha=0.5, marker='+', edgecolors='blue') # 绘制散点图，显示数据点
f = plt.plot(bjd, yfit, label='$S(x)$ Regression fit ', linewidth=2, color='orange', zorder=2, alpha=0.85) # 绘制回归拟合曲线
plt.xlabel('Barycentric Julian Date (days)') # 设置x轴标签为“日心朱利安日期（天）”
plt.ylabel('Relative Flux') # 设置y轴标签为“相对通量”
plt.title(r'Time-series Transit Light Curve $\lambda=$['+str(wsdata_all[pix1])+':'+str(wsdata_all[pix2]) + r'] $\mu$m') # 设置图表标题，包含波长信息
axs.xaxis.set_major_locator(ticker.MultipleLocator(0.01)) # 设置x轴主刻度为0.01的间隔
axs.xaxis.set_minor_locator(ticker.MultipleLocator(0.005)) # 设置x轴次刻度为0.005的间隔
axs.yaxis.set_major_locator(ticker.MultipleLocator(0.002)) # 设置y轴主刻度为0.002的间隔
axs.yaxis.set_minor_locator(ticker.MultipleLocator(0.001)) # 设置y轴次刻度为0.001的间隔
yplot = y / np.mean(y[outtransit]) # 计算归一化的y值
plt.ylim(yplot.min() * 0.999, yplot.max()*1.001) # 设置y轴范围
plt.xlim(bjd.min()-0.001, bjd.max()+0.001) # 设置x轴范围
plt.legend(loc='lower left') # 在左下角显示图例
plt.show() # 显示图形
过境和边缘暗化模型函数#
定义用于拟合例程的函数。
这些函数将接受过境和系统参数，并创建我们的完整过境光曲线模型。
\(model = T(t,\theta)\times F_0 \times S(x)\)
通过返回残差与数据进行比较。
\(y = f(t)\)
残差为：
\((y-model)/(\sigma_y)\)
为了计算过境模型，我们使用 Mandel and Agol (2002)，该模型由 H. Wakeford 用 Python 编写 (ExoTiC-ISM)。
为了计算恒星的边缘暗化，我们使用 Sing 等人（2010）的方法，该方法利用恒星模型并拟合非线性边缘暗化系数，由 H. Wakeford 用 Python 编写的模块 (ExoTiC-ISM)。
首先，根据系统参数 \(a/R_{star}\)、倾角的余弦 \(cos(i)\) 和轨道相位 \(\phi\) 计算新的轨道。
输入参数包括每个相位的行星-恒星中心之间的轨道距离 \(b\)、边缘暗化参数 (\(c_1,c_2,c_3,c_4\)) 和行星与恒星半径比 \(R_p/R_{star}\)。
@custom_model
def nonlinear_limb_darkening(x, c0=0.0, c1=0.0, c2=0.0, c3=0.0):
"""
定义具有四个参数c0, c1, c2, c3的非线性边缘变暗模型。
"""
model = (1. - (c0 * (1. - x ** (1. / 2)) + c1 * (1. - x ** (2. / 2)) + c2 * (1. - x ** (3. / 2)) + c3 *
(1. - x ** (4. / 2))))
return model
@custom_model
def quadratic_limb_darkening(x, aLD=0.0, bLD=0.0):
"""
定义具有参数aLD和bLD的二次边缘变暗模型。
"""
model = 1. - aLD * (1. - x) - bLD * (1. - x) ** (4. / 2.)
return model
def limb_dark_fit(grating, wsdata, M_H, Teff, logg, dirsen, ld_model='1D'):
"""
计算给定波长区间的恒星边缘变暗系数。
当前支持：
HST STIS G750L, G750M, G430L光栅
HST WFC3 UVIS/G280, IR/G102, IR/G141光栅
使用1D模型 - Kurucz (?)
程序来自Sing等人（2010, A&A, 510, A21）。
使用Magic等人（2015, A&A, 573, 90）的3D边缘变暗。
使用光子通量对（lambda*dlamba）进行求和。
:param grating: 字符串; 要使用的光栅（'G430L','G750L','G750M', 'G280', 'G102', 'G141'）
:param wsdata: 数组; 数据波长解决方案
:param M_H: 浮点数; 恒星金属丰度
:param Teff: 浮点数; 恒星有效温度（K）
:param logg: 浮点数; 恒星重力
:param dirsen: 字符串; 主边缘变暗目录的路径
:param ld_model: 字符串; '1D'或'3D'，在边缘变暗模型之间进行选择；默认是1D
:return: uLD: 浮点数; 线性边缘变暗系数
aLD, bLD: 浮点数; 二次边缘变暗系数
cp1, cp2, cp3, cp4: 浮点数; 三参数边缘变暗系数
c1, c2, c3, c4: 浮点数; 非线性边缘变暗系数
"""
print('您正在使用', str(ld_model), '边缘变暗模型。')
if ld_model == '1D':
direc = os.path.join(dirsen, 'Kurucz') # 设置目录为Kurucz
print('当前输入的目录:')
print(' ' + dirsen)
print(' ' + direc)
# 选择金属丰度
M_H_Grid = np.array([-0.1, -0.2, -0.3, -0.5, -1.0, -1.5, -2.0, -2.5, -3.0, -3.5, -4.0, -4.5, -5.0, 0.0, 0.1, 0.2, 0.3, 0.5, 1.0])
M_H_Grid_load = np.array([0, 1, 2, 3, 5, 7, 8, 9, 10, 11, 12, 13, 14, 17, 20, 21, 22, 23, 24])
optM = (abs(M_H - M_H_Grid)).argmin() # 找到与输入金属丰度最接近的索引
MH_ind = M_H_Grid_load[optM] # 获取对应的加载索引
# 确定使用哪个模型，通过输入的金属丰度M_H来确定所需的文件名
direc = 'Kurucz'
file_list = 'kuruczlist.sav'
sav1 = readsav(os.path.join(dirsen, file_list)) # 读取文件列表
model = bytes.decode(sav1['li'][MH_ind]) # 将字节对象转换为字符串
# 选择Teff并随后选择logg
Teff_Grid = np.array([3500, 3750, 4000, 4250, 4500, 4750, 5000, 5250, 5500, 5750, 6000, 6250, 6500])
optT = (abs(Teff - Teff_Grid)).argmin() # 找到与输入Teff最接近的索引
logg_Grid = np.array([4.0, 4.5, 5.0])
optG = (abs(logg - logg_Grid)).argmin() # 找到与输入logg最接近的索引
# 根据logg选择Teff的加载索引
if logg_Grid[optG] == 4.0:
Teff_Grid_load = np.array([8, 19, 30, 41, 52, 63, 74, 85, 96, 107, 118, 129, 138])
elif logg_Grid[optG] == 4.5:
Teff_Grid_load = np.array([9, 20, 31, 42, 53, 64, 75, 86, 97, 108, 119, 129, 139])
elif logg_Grid[optG] == 5.0:
Teff_Grid_load = np.array([10, 21, 32, 43, 54, 65, 76, 87, 98, 109, 120, 130, 140])
# 在模型文件中找到所需Teff的部分，索引T_ind告诉我们这一点。
T_ind = Teff_Grid_load[optT]
header_rows = 3 # 每个部分忽略的行数
data_rows = 1221 # 读取的数据行数
line_skip_data = (T_ind + 1) * header_rows + T_ind * data_rows # 计算需要跳过的行数
line_skip_header = T_ind * (data_rows + header_rows) # 计算跳过的头部行数
# 读取头部信息，以便获取实际的Teff、logg和M_H信息。
headerinfo = pd.read_csv(os.path.join(dirsen, direc, model), delim_whitespace=True, header=None,
skiprows=line_skip_header, nrows=1)
Teff_model = headerinfo[1].values[0] # 获取模型Teff
logg_model = headerinfo[3].values[0] # 获取模型logg
MH_model = headerinfo[6].values[0] # 获取模型M_H
MH_model = float(MH_model[1:-1]) # 转换为浮点数
print('\n与您的输入最接近的值:')
print('Teff: ', Teff_model)
print('M_H: ', MH_model)
print('log(g): ', logg_model)
# 读取数据; 数据是一个pandas对象。
data = pd.read_csv(os.path.join(dirsen, direc, model), delim_whitespace=True, header=None,
skiprows=line_skip_data, nrows=data_rows)
# 解包数据
ws = data[0].values * 10 # 导入波长数据
f0 = data[1].values / (ws * ws) # 计算f0
f1 = data[2].values * f0 / 100000. # 计算f1
f2 = data[3].values * f0 / 100000. # 计算f2
f3 = data[4].values * f0 / 100000. # 计算f3
f4 = data[5].values * f0 / 100000. # 计算f4
f5 = data[6].values * f0 / 100000. # 计算f5
f6 = data[7].values * f0 / 100000. # 计算f6
f7 = data[8].values * f0 / 100000. # 计算f7
f8 = data[9].values * f0 / 100000. # 计算f8
f9 = data[10].values * f0 / 100000. # 计算f9
f10 = data[11].values * f0 / 100000. # 计算f10
f11 = data[12].values * f0 / 100000. # 计算f11
f12 = data[13].values * f0 / 100000. # 计算f12
f13 = data[14].values * f0 / 100000. # 计算f13
f14 = data[15].values * f0 / 100000. # 计算f14
f15 = data[16].values * f0 / 100000. # 计算f15
f16 = data[17].values * f0 / 100000. # 计算f16
# 将它们合并为一个大的数组
fcalc = np.array([f0, f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11, f12, f13, f14, f15, f16])
phot1 = np.zeros(fcalc.shape[0]) # 初始化phot1数组
# 定义mu
mu = np.array([1.000, .900, .800, .700, .600, .500, .400, .300, .250, .200, .150, .125, .100, .075, .050, .025, .010])
# 传递给函数主体的参数有：ws, fcalc, phot1, mu
elif ld_model == '3D':
direc = os.path.join(dirsen, '3DGrid') # 设置目录为3DGrid
print('当前输入的目录:')
print(' ' + dirsen)
print(' ' + direc)
# 选择金属丰度
M_H_Grid = np.array([-3.0, -2.0, -1.0, 0.0]) # 3D模型中可用的金属丰度值
M_H_Grid_load = ['30', '20', '10', '00'] # 对应于各个可用M_H值的标识符
optM = (abs(M_H - M_H_Grid)).argmin() # 找到与输入M_H最接近的索引
# 选择Teff
Teff_Grid = np.array([4000, 4500, 5000, 5500, 5777, 6000, 6500, 7000]) # 3D模型中可用的Teff值
optT = (abs(Teff - Teff_Grid)).argmin() # 找到与输入Teff最接近的索引
# 根据Teff选择logg。如果对于一个Teff给出了多个logg可能性，选择与用户输入最接近的。
if Teff_Grid[optT] == 4000:
logg_Grid = np.array([1.5, 2.0, 2.5])
optG = (abs(logg - logg_Grid)).argmin()
elif Teff_Grid[optT] == 4500:
logg_Grid = np.array([2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0])
optG = (abs(logg - logg_Grid)).argmin()
elif Teff_Grid[optT] == 5000:
logg_Grid = np.array([2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0])
optG = (abs(logg - logg_Grid)).argmin()
elif Teff_Grid[optT] == 5500:
logg_Grid = np.array([3.0, 3.5, 4.0, 4.5, 5.0])
optG = (abs(logg - logg_Grid)).argmin()
elif Teff_Grid[optT] == 5777:
logg_Grid = np.array([4.4])
optG = 0
elif Teff_Grid[optT] == 6000:
logg_Grid = np.array([3.5, 4.0, 4.5])
optG = (abs(logg - logg_Grid)).argmin()
elif Teff_Grid[optT] == 6500:
logg_Grid = np.array([4.0, 4.5])
optG = (abs(logg - logg_Grid)).argmin()
elif Teff_Grid[optT] == 7000:
logg_Grid = np.array([4.5])
optG = 0
# 选择Teff和Log g。Mtxt, Ttxt和Gtxt将作为字符串组合以加载正确的文件。
Mtxt = M_H_Grid_load[optM]
Ttxt = "{:2.0f}".format(Teff_Grid[optT] / 100) # 格式化Teff
if Teff_Grid[optT] == 5777:
Ttxt = "{:4.0f}".format(Teff_Grid[optT]) # 特殊处理Teff
Gtxt = "{:2.0f}".format(logg_Grid[optG] * 10) # 格式化logg
file = 'mmu_t' + Ttxt + 'g' + Gtxt + 'm' + Mtxt + 'v05.flx' # 生成文件名
print('文件名:', file)
# 从IDL .sav文件中读取数据
sav = readsav(os.path.join(direc, file)) # 读取IDL .sav文件
ws = sav['mmd'].lam[0] # 读取波长
flux = sav['mmd'].flx # 读取通量
Teff_model = Teff_Grid[optT] # 获取模型Teff
logg_model = logg_Grid[optG] # 获取模型logg
MH_model = str(M_H_Grid[optM]) # 获取模型M_H
print('\n与您的输入最接近的值:')
print('Teff : ', Teff_model)
print('M_H : ', MH_model)
print('log(g): ', logg_model)
# 解包flux数据
f0 = flux[0]
f1 = flux[1]
f2 = flux[2]
f3 = flux[3]
f4 = flux[4]
f5 = flux[5]
f6 = flux[6]
f7 = flux[7]
f8 = flux[8]
f9 = flux[9]
f10 = flux[10]
# 将它们合并为一个大的数组
fcalc = np.array([f0, f1, f2, f3, f4, f5, f6, f7, f8, f9, f10])
phot1 = np.zeros(fcalc.shape[0]) # 初始化phot1数组
# 从网格中获取Mu
mu = sav['mmd'].mu # 读取mu
# 传递给函数主体的参数有：ws, fcalc, phot1, mu
# 加载响应函数并插值到kurucz模型网格上
# 对于STIS
if grating == 'G430L':
sav = readsav(os.path.join(dirsen, 'G430L.STIS.sensitivity.sav')) # 读取灵敏度文件
wssens = sav['wssens'] # 获取波长灵敏度
sensitivity = sav['sensitivity'] # 获取灵敏度
wdel = 3 # 设置波长间隔
if grating == 'G750M':
sav = readsav(os.path.join(dirsen, 'G750M.STIS.sensitivity.sav')) # 读取灵敏度文件
wssens = sav['wssens'] # 获取波长灵敏度
sensitivity = sav['sensitivity'] # 获取灵敏度
wdel = 0.554 # 设置波长间隔
if grating == 'G750L':
sav = readsav(os.path.join(dirsen, 'G750L.STIS.sensitivity.sav')) # 读取灵敏度文件
wssens = sav['wssens'] # 获取波长灵敏度
sensitivity = sav['sensitivity'] # 获取灵敏度
wdel = 4.882 # 设置波长间隔
# 对于WFC3
if grating == 'G141': # http://www.stsci.edu/hst/acs/analysis/reference_files/synphot_tables.html
sav = readsav(os.path.join(dirsen, 'G141.WFC3.sensitivity.sav')) # 读取灵敏度文件
wssens = sav['wssens'] # 获取波长灵敏度
sensitivity = sav['sensitivity'] # 获取灵敏度
wdel = 1 # 设置波长间隔
if grating == 'G102': # http://www.stsci.edu/hst/acs/analysis/reference_files/synphot_tables.html
sav = readsav(os.path.join(dirsen, 'G141.WFC3.sensitivity.sav')) # 读取灵敏度文件
wssens = sav['wssens'] # 获取波长灵敏度
sensitivity = sav['sensitivity'] # 获取灵敏度
wdel = 1 # 设置波长间隔
if grating == 'G280': # http://www.stsci.edu/hst/acs/analysis/reference_files/synphot_tables.html
sav = readsav(os.path.join(dirsen, 'G280.WFC3.sensitivity.sav')) # 读取灵敏度文件
wssens = sav['wssens'] # 获取波长灵敏度
sensitivity = sav['sensitivity'] # 获取灵敏度
wdel = 1 # 设置波长间隔
# 对于JWST
if grating == 'NIRSpecPrism': # http://www.stsci.edu/hst/acs/analysis/reference_files/synphot_tables.html
sav = readsav(os.path.join(dirsen, 'NIRSpec.prism.sensitivity.sav')) # 读取灵敏度文件
wssens = sav['wssens'] # 获取波长灵敏度
sensitivity = sav['sensitivity'] # 获取灵敏度
wdel = 12 # 设置波长间隔
widek = np.arange(len(wsdata)) # 创建宽度数组
wsHST = wssens # 赋值波长灵敏度
wsHST = np.concatenate((np.array([wsHST[0] - wdel - wdel, wsHST[0] - wdel]),
wsHST,
np.array([wsHST[len(wsHST) - 1] + wdel,
wsHST[len(wsHST) - 1] + wdel + wdel]))) # 扩展波长灵敏度数组
respoutHST = sensitivity / np.max(sensitivity) # 归一化灵敏度
respoutHST = np.concatenate((np.zeros(2), respoutHST, np.zeros(2))) # 在灵敏度前后添加零
inter_resp = interp1d(wsHST, respoutHST, bounds_error=False, fill_value=0) # 插值灵敏度
respout = inter_resp(ws) # 将灵敏度插值到模型波长网格上
wsdata = np.concatenate((np.array([wsdata[0] - wdel - wdel, wsdata[0] - wdel]), wsdata,
np.array([wsdata[len(wsdata) - 1] + wdel, wsdata[len(wsdata) - 1] + wdel + wdel]))) # 扩展数据波长数组
respwavebin = wsdata / wsdata * 0.0 # 初始化响应波长数组
widek = widek + 2 # 需要添加两个索引以补偿用两个零填充
respwavebin[widek] = 1.0 # 设置响应波长数组的值
data_resp = interp1d(wsdata, respwavebin, bounds_error=False, fill_value=0) # 插值数据
reswavebinout = data_resp(ws) # 将数据插值到模型波长网格上
# 对光谱进行积分以生成合成光度点。
for i in range(fcalc.shape[0]): # 遍历不同角度的光谱
fcal = fcalc[i, :] # 获取当前光谱
Tot = int_tabulated(ws, ws * respout * reswavebinout) # 计算总积分
phot1[i] = (int_tabulated(ws, ws * respout * reswavebinout * fcal, sort=True)) / Tot # 计算光度
if ld_model == '1D':
yall = phot1 / phot1[0] # 归一化光度
elif ld_model == '3D':
yall = phot1 / phot1[10] # 归一化光度
x = mu[1:] # 获取波长
y = yall[1:] # 获取光度
# 开始拟合不同的模型
fitter = LevMarLSQFitter() # 创建拟合器
# 拟合四参数非线性边缘变暗模型并获取拟合变量c1, c2, c3, c4。
corot_4_param = nonlinear_limb_darkening() # 创建模型
corot_4_param = fitter(corot_4_param, x, y) # 拟合模型
c1, c2, c3, c4 = corot_4_param.parameters # 获取拟合参数
# 拟合三参数非线性边缘变暗模型并获取拟合变量cp2, cp3, cp4（cp1 = 0）。
corot_3_param = nonlinear_limb_darkening() # 创建模型
corot_3_param.c0.fixed = True # 3参数模型的c0固定为0.0
corot_3_param = fitter(corot_3_param, x, y) # 拟合模型
cp1, cp2, cp3, cp4 = corot_3_param.parameters # 获取拟合参数
# 拟合二次边缘变暗模型并获取拟合参数aLD和bLD。
quadratic = quadratic_limb_darkening() # 创建模型
quadratic = fitter(quadratic, x, y) # 拟合模型
aLD, bLD = quadratic.parameters # 获取拟合参数
# 拟合线性边缘变暗模型并获取拟合变量uLD。
linear = nonlinear_limb_darkening() # 创建模型
linear.c0.fixed = True # 固定c0
linear.c2.fixed = True # 固定c2
linear.c3.fixed = True # 固定c3
linear = fitter(linear, x, y) # 拟合模型
uLD = linear.c1.value # 获取拟合参数
print('\n边缘变暗参数:')
print("4param \t{:0.8f}\t{:0.8f}\t{:0.8f}\t{:0.8f}".format(c1, c2, c3, c4))
print("3param \t{:0.8f}\t{:0.8f}\t{:0.8f}".format(cp2, cp3, cp4))
print("Quad \t{:0.8f}\t{:0.8f}".format(aLD, bLD))
print("Linear \t{:0.8f}".format(uLD))
return uLD, c1, c2, c3, c4, cp1, cp2, cp3, cp4, aLD, bLD
def int_tabulated(X, F, sort=False):
Xsegments = len(X) - 1 # 计算X的段数
# 将向量按升序排序。
if not sort:
ii = np.argsort(X) # 获取排序索引
X = X[ii] # 排序X
F = F[ii] # 排序F
while (Xsegments % 4) != 0: # 确保段数是4的倍数
Xsegments = Xsegments + 1
Xmin = np.min(X) # 获取X的最小值
Xmax = np.max(X) # 获取X的最大值
# 统一步长。
h = (Xmax + 0.0 - Xmin) / Xsegments # 计算步长
# 在Xgrid上计算插值。
z = splev(h * np.arange(Xsegments + 1) + Xmin, splrep(X, F)) # 计算插值
# 使用5点牛顿-科特斯公式计算积分。
ii = (np.arange((len(z) - 1) / 4, dtype=int) + 1) * 4 # 计算索引
return np.sum(2.0 * h * (7.0 * (z[ii - 4] + z[ii]) + 32.0 * (z[ii - 3] + z[ii - 1]) + 12.0 * z[ii - 2]) / 45.0) # 返回积分结果
现在定义过境模型函数#
def occultnl(rl, c1, c2, c3, c4, b0):
"""
MANDEL & AGOL (2002) transit model.
:param rl: float, transit depth (Rp/R*)
:param c1: float, limb darkening parameter 1
:param c2: float, limb darkening parameter 2
:param c3: float, limb darkening parameter 3
:param c4: float, limb darkening parameter 4
:param b0: impact parameter in stellar radii
:return: mulimb0: limb-darkened transit model, mulimbf: lightcurves for each component that you put in the model
"""
mulimb0 = occultuniform(b0, rl) # 计算均匀源的光曲线
bt0 = b0 # 保存初始的影响参数
fac = np.max(np.abs(mulimb0 - 1)) # 计算光度因子的最大绝对值
if fac == 0:
fac = 1e-6 # 如果fac为0，设置一个小的正数以避免除零错误
omega = 4 * ((1 - c1 - c2 - c3 - c4) / 4 + c1 / 5 + c2 / 6 + c3 / 7 + c4 / 8) # 计算权重因子
nb = len(b0) # 获取影响参数数组的长度
indx = np.where(mulimb0 != 1.0)[0] # 找到光度因子不等于1的索引
if len(indx) == 0:
indx = -1 # 如果没有找到，设置为-1
mulimb = mulimb0[indx] # 获取对应的光度因子
mulimbf = np.zeros((5, nb)) # 初始化光曲线数组
mulimbf[0, :] = mulimbf[0, :] + 1. # 设置第一个组件的光度
mulimbf[1, :] = mulimbf[1, :] + 0.8 # 设置第二个组件的光度
mulimbf[2, :] = mulimbf[2, :] + 2 / 3 # 设置第三个组件的光度
mulimbf[3, :] = mulimbf[3, :] + 4 / 7 # 设置第四个组件的光度
mulimbf[4, :] = mulimbf[4, :] + 0.5 # 设置第五个组件的光度
nr = np.int64(2) # 初始化迭代次数
dmumax = 1.0 # 初始化最大变化量
while (dmumax > fac * 1.e-3) and (nr <= 131072): # 迭代直到满足收敛条件
mulimbp = mulimb # 保存当前光度因子
nr = nr * 2 # 每次迭代翻倍
dt = 0.5 * np.pi / nr # 计算时间步长
t = dt * np.arange(nr + 1) # 生成时间数组
th = t + 0.5 * dt # 计算中间时间
r = np.sin(t) # 计算正弦值
sig = np.sqrt(np.cos(th[nr - 1])) # 计算信号强度
# 计算不同光度因子的中间值
mulimbhalf = sig ** 3 * mulimb0[indx] / (1 - r[nr - 1])
mulimb1 = sig ** 4 * mulimb0[indx] / (1 - r[nr - 1])
mulimb3half = sig ** 5 * mulimb0[indx] / (1 - r[nr - 1])
mulimb2 = sig ** 6 * mulimb0[indx] / (1 - r[nr - 1])
for i in range(1, nr): # 遍历每个时间步
mu = occultuniform(b0[indx] / r[i], rl / r[i]) # 计算当前光度因子
sig1 = np.sqrt(np.cos(th[i - 1])) # 计算前一个时间步的信号强度
sig2 = np.sqrt(np.cos(th[i])) # 计算当前时间步的信号强度
# 更新光度因子
mulimbhalf = mulimbhalf + r[i] ** 2 * mu * (sig1 ** 3 / (r[i] - r[i - 1]) - sig2 ** 3 / (r[i + 1] - r[i]))
mulimb1 = mulimb1 + r[i] ** 2 * mu * (sig1 ** 4 / (r[i] - r[i - 1]) - sig2 ** 4 / (r[i + 1] - r[i]))
mulimb3half = mulimb3half + r[i] ** 2 * mu * (sig1 ** 5 / (r[i] - r[i - 1]) - sig2 ** 5 / (r[i + 1] - r[i]))
mulimb2 = mulimb2 + r[i] ** 2 * mu * (sig1 ** 6 / (r[i] - r[i - 1]) - sig2 ** 6 / (r[i + 1] - r[i]))
# 计算新的光度因子
mulimb = ((1 - c1 - c2 - c3 - c4) * mulimb0[indx] +
c1 * mulimbhalf * dt +
c2 * mulimb1 * dt +
c3 * mulimb3half * dt +
c4 * mulimb2 * dt) / omega
ix1 = np.where(mulimb + mulimbp != 0.)[0] # 找到光度因子不为零的索引
if len(ix1) == 0:
ix1 = -1 # 如果没有找到，设置为-1
# 计算最大变化量
dmumax = np.max(np.abs(np.atleast_1d(mulimb)[ix1] - np.atleast_1d(mulimbp)[ix1]) / (
np.atleast_1d(mulimb)[ix1] + np.atleast_1d(mulimbp)[ix1]))
# 更新光度曲线
mulimbf[0, indx] = np.atleast_1d(mulimb0)[indx]
mulimbf[1, indx] = mulimbhalf * dt
mulimbf[2, indx] = mulimb1 * dt
mulimbf[3, indx] = mulimb3half * dt
mulimbf[4, indx] = mulimb2 * dt
np.atleast_1d(mulimb0)[indx] = mulimb # 更新光度因子
b0 = bt0 # 恢复初始影响参数
return mulimb0, mulimbf # 返回光度因子和光曲线
def occultuniform(b0, w):
"""
Compute the lightcurve for occultation of a uniform source without microlensing (Mandel & Agol 2002).
:param b0: array; impact parameter in units of stellar radii
:param w: array; occulting star size in units of stellar radius
:return: muo1: float; fraction of flux at each b0 for a uniform source
"""
if np.abs(w - 0.5) < 1.0e-3: # 如果w接近0.5，设置为0.5
w = 0.5
nb = len(np.atleast_1d(b0)) # 获取影响参数数组的长度
muo1 = np.zeros(nb) # 初始化光度因子数组
for i in range(nb): # 遍历每个影响参数
z = np.atleast_1d(b0)[i] # 获取当前影响参数
if z >= 1 + w: # 如果影响参数大于等于1+w
muo1[i] = 1.0 # 光度因子为1
continue
if w >= 1 and z <= w - 1: # 如果w大于等于1且影响参数小于等于w-1
muo1[i] = 0.0 # 光度因子为0
continue
if z >= np.abs(1 - w) and z <= 1 + w: # 如果影响参数在[|1-w|, 1+w]之间
kap1 = np.arccos(np.min(np.append((1 - w ** 2 + z ** 2) / 2 / z, 1.))) # 计算角度
kap0 = np.arccos(np.min(np.append((w ** 2 + z ** 2 - 1) / 2 / w / z, 1.))) # 计算角度
lambdae = w ** 2 * kap0 + kap1 # 计算光度因子
lambdae = (lambdae - 0.5 * np.sqrt(np.max(np.append(4. * z ** 2 - (1 + z ** 2 - w ** 2) ** 2, 0.)))) / np.pi # 归一化
muo1[i] = 1 - lambdae # 更新光度因子
if z <= 1 - w: # 如果影响参数小于等于1-w
muo1[i] = 1 - w ** 2 # 更新光度因子
continue
return muo1 # 返回光度因子数组
现在定义生成过境光曲线的函数并将其与数据进行比较#
# 函数用于调用和计算模型
def residual(p, phase, x, y, err, c1, c2, c3, c4):
# 计算新的轨道
b0 = p['a_Rs'].value * np.sqrt((np.sin(phase * 2 * np.pi)) ** 2 + (p['cosinc'].value * np.cos(phase * 2 * np.pi)) ** 2)
# 选择在第一次和第四次接触之间的索引
intransit = (b0 - p['rprs'].value < 1.0E0).nonzero()
# 创建光曲线模型，初始所有值设为1.0
light_curve = b0 / b0
# 计算光度减弱，使用Madel和Agol的方法
mulimb0, mulimbf = occultnl(p['rprs'].value, c1, c2, c3, c4, b0[intransit])
# 在过境期间更新光曲线
light_curve[intransit] = mulimb0
# 计算模型：过境模型是基线通量 X 过境模型 X 系统模型
model = (light_curve) * p['f0'].value * (p['Fslope'].value * phase +
p['xsh'].value * shx +
p['x2sh'].value * shx**2. +
p['ysh'].value * shy +
p['y2sh'].value * shy**2. +
p['xysh'].value * shy * shx +
p['comm'].value * common_mode + 1.0)
# 计算当前的卡方值
chi2now = np.sum((y - model)**2 / err**2)
# 计算残差的标准差
res = np.std((y - model) / p['f0'].value)
# 打印当前的rprs值、卡方值和散射值
print("rprs: ", p['rprs'].value, "current chi^2=", chi2now, ' scatter ', res, end="\r")
# 返回标准化的残差
return (y - model) / err
# return np.sum((y-model)**2/err**2) # 这行代码被注释掉了
一个函数也被定义为返回仅仅是过境模型 \(T(t,\theta)\)。
def model_fine(p): # 创建用于绘图目的的精细网格的过境模型
# 计算b0，表示行星相对于恒星的距离
b0 = p['a_Rs'].value * np.sqrt((np.sin(phase_fine * 2 * np.pi)) ** 2 + (p['cosinc'].value * np.cos(phase_fine * 2 * np.pi)) ** 2)
# 使用occultnl函数计算模型的光度，mulimb0和mulimbf分别表示不同的光度
mulimb0, mulimbf = occultnl(p['rprs'].value, c1, c2, c3, c4, b0) # Madel 和 Agol
# 将模型光度赋值给model_fine
model_fine = mulimb0
return model_fine # 返回计算得到的过境模型
现在在示例光曲线中添加一个过境模型。在这里，我们已经计算了边缘变暗系数（limb darkening coefficients），然后在过境光曲线中使用它们。
wave1 = wsdata_all[pix1] # 从wsdata_all中提取第pix1个像素的波长值，赋值给wave1
wave2 = wsdata_all[pix2] # 从wsdata_all中提取第pix2个像素的波长值，赋值给wave2
bin_wave_index = ((wsdata_all > wave1) & (wsdata_all <= wave2)).nonzero() # 找到wsdata_all中在wave1和wave2之间的波长索引
wsdata = wsdata_all[bin_wave_index]*1E4 # 选择波长区间内的值并将单位从微米转换为埃（angstroms）
_uLD, c1, c2, c3, c4, _cp1, _cp2, _cp3, _cp4, aLD, bLD = limb_dark_fit(grating, wsdata, M_H, Teff, logg, limb_dark_directory, ld_model)
# 调用limb_dark_fit函数进行光度暗化拟合，返回多个参数
现在运行过境模型。
过境参数，例如倾角（inclination）和 \(a/R_{star}\)，已在笔记本的开头设置。
# 运行过境模型
rl = 0.0825 # 行星与恒星半径比
# 计算接触点的投影参数 b0
b0 = a_Rs * np.sqrt((np.sin(phase * 2 * np.pi)) ** 2 + (np.cos(inc) * np.cos(phase * 2 * np.pi)) ** 2)
# 选择在第一次和第四次接触之间的索引
intransit = (b0-rl < 1.0E0).nonzero()
# 使用 Mandel & Agol 的非线性暗化过境模型
mulimb0, mulimbf = occultnl(rl, c1, c2, c3, c4, b0)
# 计算模型
model = mulimb0 * yfit
# 绘图设置
plt.rcParams['figure.figsize'] = [10.0, 7.0] # 图形尺寸
msize = plt.rcParams['lines.markersize'] ** 2. # 默认标记大小
# 创建图形
fig = plt.figure(constrained_layout=True)
# 创建网格布局
gs = fig.add_gridspec(3, 1, hspace=0.00, wspace=0.00)
# 添加第一个子图
ax1 = fig.add_subplot(gs[0:2, :])
# 绘制数据散点图
ax1.scatter(bjd, y/np.mean(y[outtransit]), label='$f(t)$ 数据', zorder=1, s=msize*0.75, linewidth=1, alpha=0.5, marker='+', edgecolors='blue')
# 绘制模型拟合曲线
ax1.plot(bjd, model, label='$S(x)$ 回归拟合 ', linewidth=2, color='orange', zorder=2, alpha=0.85)
# 设置 x 轴标签
ax1.xaxis.set_ticklabels([])
# 设置 y 轴标签
plt.ylabel('相对通量')
# 设置图形标题
plt.title(r'时间序列过境光曲线 $\lambda=$['+str(wsdata_all[pix1])+':'+str(wsdata_all[pix2]) + r'] $\mu$m')
# 设置 x 轴主刻度和次刻度
ax1.xaxis.set_major_locator(ticker.MultipleLocator(0.01))
ax1.xaxis.set_minor_locator(ticker.MultipleLocator(0.005))
# 设置 y 轴主刻度和次刻度
ax1.yaxis.set_major_locator(ticker.MultipleLocator(0.002))
ax1.yaxis.set_minor_locator(ticker.MultipleLocator(0.001))
# 计算标准化的 y 值
yplot = y/np.mean(y[outtransit])
# 设置 y 轴范围
plt.ylim(yplot.min()*0.999, yplot.max()*1.001)
# 设置 x 轴范围
plt.xlim(bjd.min()-0.001, bjd.max()+0.001)
# 添加图例
plt.legend()
# 将 ax1 添加到图形中
fig.add_subplot(ax1)
# 残差图
# 添加第二个子图
ax2 = fig.add_subplot(gs[2, :])
# 绘制残差散点图
ax2.scatter(bjd, 1E6*(y/np.mean(y[outtransit])-model), label='$f(t)$ 数据', zorder=1, s=msize*0.75, linewidth=1, alpha=0.5, marker='+', edgecolors='blue')
# 计算并绘制残差的统计信息
wsb, wsb_bin_edges, binnumber = stats.binned_statistic(bjd, 1E6*(y/np.mean(y[outtransit])-model), bins=256)
plt.scatter(wsb_bin_edges[1:], wsb, linewidth=2, alpha=0.75, facecolors='orange', edgecolors='none', marker='o', zorder=25)
# 设置 x 轴和 y 轴标签
plt.xlabel('日心朱利安日期 (天)')
plt.ylabel('残差 (ppm)')
# 设置 x 轴主刻度和次刻度
ax2.xaxis.set_major_locator(ticker.MultipleLocator(0.01))
ax2.xaxis.set_minor_locator(ticker.MultipleLocator(0.005))
# 计算标准化的 y 值
yplot = y / np.mean(y[outtransit])
# 设置 x 轴范围
plt.xlim(bjd.min()-0.001, bjd.max()+0.001)
# 将 ax2 添加到图形中
fig.add_subplot(ax2)
# 显示图形
plt.show()
# 打印卡方值
err = np.sqrt(y) / np.mean(y[outtransit])
print('Chi^2 = '+str(np.sum((y/np.mean(y[outtransit])-model)**2/err**2)))
# 打印残差标准差
print('残差标准差 : '+str(1E6*np.std((y/np.mean(y[outtransit])-model)))+' ppm')
print('256 Bin 标准差 :'+str(np.std(wsb))+' ppm')
请注意，模型的过境深度与数据相比略显过深。行星半径需要更小，参数 \(rl\) 更接近 0.08。作为练习，您可以重新运行上述单元，将行星半径更改为 \(rl\)=0.0805，并将 \(\chi^2\) 值与之前的默认值进行比较（在 \(rl\) = 0.0825 时，\(\chi^2\)=9265.4）。
FIT过渡光曲线#
现在我们可以对每个光曲线进行拟合，通过最小二乘法优化拟合参数。在这里，使用Levenberg-Marquardt拟合来寻找\(\chi^2\)最小值，并使用lmfit包（https://lmfit.github.io/lmfit-py/fitting.html）估计不确定性。
实际上，首先我们会拟合白光光曲线，该光曲线是对整个1D光谱中所有波长的总和。这可以用来拟合系统参数，例如倾角和过渡时间，然后光谱通道被固定为这些值，因为它们与波长无关。然而，CV3数据需要去除灯光的整体变化，这使得我们无法使用这些数据进行白光光曲线分析。在这里，我们继续对光谱光曲线的光谱箱进行拟合。
步骤如下：
选择波长箱
从恒星模型中计算每个箱的光晕暗化系数。
在过渡外数据上进行初始线性回归，以开始系统拟合参数，这大大加快了拟合速度，因为这些参数接近其全局最小值。
开始拟合，并在最小化过程中输出一些统计信息。
一旦找到最佳拟合，将显示多个统计信息。
最后，生成几个图表并将其存储为PDF，然后开始下一个箱。
这些步骤对每个光谱箱进行。
在这个例子中，行星半径被设置为在拟合中变化，同时基线通量和仪器系统参数也会变化。
设置待拟合的波长#
光谱必须在波长上进行分箱，以获得足够的计数，以达到所需的约100 ppm水平。
光谱在大约像素100到400之间有显著的计数，我们从像素 \(k0\) 开始，并以 \(wk\) 像素进行分箱。
还定义了几个数组。
k0 = 113 # 起始波长索引
kend = 392 # 结束波长索引
wk = 15 # 每个bin的宽度
# 计算bin的数量
number_of_bins = int((kend - k0) / wk)
# 初始化各个数组
wsd = np.zeros((number_of_bins)) # 存储波长数据
werr = np.zeros((number_of_bins)) # 存储波长误差
rprs = np.zeros((number_of_bins)) # 存储半径比
rerr = np.zeros((number_of_bins)) # 存储半径比误差
sig_r = np.zeros((number_of_bins)) # 存储半径信号
sig_w = np.zeros((number_of_bins)) # 存储波长信号
beta = np.zeros((number_of_bins)) # 存储beta参数
depth = np.zeros((number_of_bins)) # 存储深度
depth_err = np.zeros((number_of_bins)) # 存储深度误差
循环遍历波长区间，拟合每个光曲线#
注意：此步骤需要相当长的时间来完成（约20分钟，每个波长区间几分钟）#
每个波长区间都将拟合过渡+系统误差模型。各种输出，包括图表，将被保存。
可以跳到下一个单元以加载预先计算的结果。
k = k0 # 设置起始波长
# --------------------------------------------------------------------------
# 循环遍历波长区间并为每个区间进行拟合
for bin in range(0, number_of_bins):
# 选择波长区间
wave1 = wsdata_all[k] # 当前波长区间的下限
wave2 = wsdata_all[k+wk] # 当前波长区间的上限
# 选择波长区间的索引
bin_wave_index = ((wsdata_all > wave1) & (wsdata_all <= wave2)).nonzero() # 获取在当前波长区间内的索引
# 制作光曲线区间
wave_bin_counts = np.sum(all_spec_1D[k+1:k+wk, :], axis=0) # 对波长像素进行求和
wave_bin_counts_err = np.sqrt(wave_bin_counts) # 采用光子噪声作为误差
# 计算边缘暗化
wsdata = wsdata_all[bin_wave_index] * 1E4 # 选择波长区间的值（从微米转换为埃）
_uLD, c1, c2, c3, c4, _cp1, _cp2, _cp3, _cp4, aLD, bLD = limb_dark_fit(grating, wsdata, M_H, Teff, logg, limb_dark_directory, ld_model) # 进行边缘暗化拟合
print('\nc1 = {}'.format(c1)) # 打印边缘暗化系数c1
print('c2 = {}'.format(c2)) # 打印边缘暗化系数c2
print('c3 = {}'.format(c3)) # 打印边缘暗化系数c3
print('c4 = {}'.format(c4)) # 打印边缘暗化系数c4
print('')
# u = [c1,c2,c3,c4] # 边缘暗化系数
# u = [aLD, bLD]
# 制作初始模型
# 设置LMFIT
x = bjd # X数据
y = wave_bin_counts # Y数据
err = wave_bin_counts_err # Y误差
# 在无过境数据上执行快速线性回归以获得准确的起始探测器拟合值
if wave1 > 2.7 and wave1 < 3.45: # 根据波长选择不同的回归器
regressor.fit(XXX[outtransit], y[outtransit] / np.mean(y[outtransit])) # 拟合无过境数据
else:
regressor.fit(XX[outtransit], y[outtransit] / np.mean(y[outtransit])) # 拟合无过境数据
# 创建LMFIT参数集 https://lmfit.github.io/lmfit-py/parameters.html
# class Parameter(name, value=None, vary=True, min=- inf, max=inf, expr=None, brute_step=None, user_data=None)¶
# 设置vary=0以固定
# 设置vary=1以进行拟合
p = lmfit.Parameters() # 存储L-M拟合参数的对象
p.add('cosinc', value=np.cos(inc), vary=0) # 倾角，固定cos(倾角)
p.add('rho_star', value=rho_star, vary=0) # 恒星密度
p.add('a_Rs', value=a_Rs, vary=0) # a/Rstar
p.add('rprs', value=rp, vary=1, min=0, max=1) # 行星与恒星半径比
p.add('t0', value=t0, vary=0) # 过境T0
p.add('f0', value=np.mean(y[outtransit]), vary=1, min=0) # 基线通量
p.add('ecc', value=ecc, vary=0, min=0, max=1) # 偏心率
p.add('omega', value=omega, vary=0) # 近日点参数
# 在水特征中开启线性斜率以考虑在低温测试期间窗户上H2O冰的变化
if wave1 > 2.7 and wave1 < 3.45:
p.add('Fslope', value=regressor.coef_[6], vary=1) # 轨道相位
else:
p.add('Fslope', value=0, vary=0) # 轨道相位
p.add('xsh', value=regressor.coef_[0], vary=1) # 探测器X偏移去趋势
p.add('ysh', value=regressor.coef_[1], vary=1) # 探测器Y偏移去趋势
p.add('x2sh', value=regressor.coef_[2], vary=1) # 探测器X^2偏移去趋势
p.add('y2sh', value=regressor.coef_[3], vary=1) # 探测器Y^2偏移去趋势
p.add('xysh', value=regressor.coef_[4], vary=1) # 探测器XY去趋势
p.add('comm', value=regressor.coef_[5], vary=1) # 共模去趋势
# 执行最小化 https://lmfit.github.io/lmfit-py/fitting.html
# 创建最小化器
# mini = lmfit.Minimizer(residual, p, nan_policy='omit',fcn_args=(phase,x,y,err)
print('')
print('Fitting Bin', bin, ' Wavelength =', np.mean(wsdata) / 1E4, ' Range= [', wave1, ':', wave2, ']') # 打印当前拟合的波长区间
# 使用Levenberg-Marquardt求解
result = lmfit.minimize(residual, params=p, args=(phase, x, y, err, c1, c2, c3, c4)) # 最小化残差
# result = mini.minimize(method='emcee')
print('')
print('Re-Fitting Bin', bin, ' Wavelength =', np.mean(wsdata) / 1E4, ' Range= [', wave1, ':', wave2, ']') # 打印重新拟合的波长区间
# --------------------------------------------------------------------------
print("")
print("redchi", result.redchi) # 打印红卡方
print("chi2", result.chisqr) # 打印卡方
print("nfree", result.nfree) # 打印自由度
print("bic", result.bic) # 打印贝叶斯信息准则
print("aic", result.aic) # 打印赤池信息准则
print("L-M FIT Variable") # 打印L-M拟合变量
print(lmfit.fit_report(result.params)) # 打印拟合报告
text_file = open(save_directory + 'JWST_NIRSpec_Prism_fit_light_curve_bin' + str(bin) + '_statistics.txt', "w") # 打开文件以保存统计信息
n = text_file.write("\nredchi " + str(result.redchi)) # 写入红卡方
n = text_file.write("\nchi2 " + str(result.chisqr)) # 写入卡方
n = text_file.write("\nnfree " + str(result.nfree)) # 写入自由度
n = text_file.write("\nbic " + str(result.bic)) # 写入贝叶斯信息准则
n = text_file.write("\naic " + str(result.aic)) # 写入赤池信息准则
n = text_file.write(lmfit.fit_report(result.params)) # 写入拟合报告
# file-output.py
# 更新最佳拟合参数
p['rho_star'].value = result.params['rho_star'].value # 更新恒星密度
p['cosinc'].value = result.params['cosinc'].value # 更新cos(倾角)
p['rprs'].value = result.params['rprs'].value # 更新行星与恒星半径比
p['t0'].value = result.params['t0'].value # 更新过境T0
p['f0'].value = result.params['f0'].value # 更新基线通量
p['Fslope'].value = result.params['Fslope'].value # 更新斜率
p['xsh'].value = result.params['xsh'].value # 更新X偏移
p['ysh'].value = result.params['ysh'].value # 更新Y偏移
p['x2sh'].value = result.params['x2sh'].value # 更新X^2偏移
p['y2sh'].value = result.params['y2sh'].value # 更新Y^2偏移
p['xysh'].value = result.params['xysh'].value # 更新XY偏移
p['comm'].value = result.params['comm'].value # 更新共模偏移
# 更新拟合光谱数组
wsd[bin] = np.mean(wsdata) / 1E4 # 更新波长中心
werr[bin] = (wsdata.max() - wsdata.min()) / 2E4 # 更新波长半宽
rprs[bin] = result.params['rprs'].value # 更新行星与恒星半径比
rerr[bin] = result.params['rprs'].stderr # 更新行星与恒星半径比的误差
# 计算最佳拟合模型
final_model = y - result.residual * err # 计算最终模型
final_model_fine = model_fine(p) # 计算细化模型
# 更多统计信息
resid = (y - final_model) / p['f0'].value # 计算残差
residppm = 1E6 * (y - final_model) / p['f0'].value # 计算残差的ppm
residerr = err / p['f0'].value # 计算残差误差
sigma = np.std((y - final_model) / p['f0'].value) * 1E6 # 计算残差的标准差
print("Residual standard deviation (ppm) : ", 1E6 * np.std((y - final_model) / p['f0'].value)) # 打印残差标准差
print("Photon noise (ppm) : ", (1 / np.sqrt(p['f0'].value)) * 1E6) # 打印光子噪声
print("Photon noise performance (%) : ", (1 / np.sqrt(p['f0'].value)) * 1E6 / (sigma) * 100) # 打印光子噪声性能
n = text_file.write("\nResidual standard deviation (ppm) : " + str(1E6 * np.std((y - final_model) / p['f0'].value))) # 写入残差标准差
n = text_file.write("\nPhoton noise (ppm) : " + str((1 / np.sqrt(p['f0'].value)) * 1E6)) # 写入光子噪声
n = text_file.write("\nPhoton noise performance (%) : " + str((1 / np.sqrt(p['f0'].value)) * 1E6 / (sigma) * 100)) # 写入光子噪声性能
# 使用分箱技术测量红噪声
sig0 = np.std(resid) # 计算残差的标准差
bins = number_of_images / binmeasure # 计算分箱数
wsb, wsb_bin_edges, binnumber = stats.binned_statistic(bjd, resid, bins=bins) # 进行分箱统计
sig_binned = np.std(wsb) # 计算分箱后的标准差
sigrednoise = np.sqrt(sig_binned**2 - sig0**2 / binmeasure) # 计算红噪声
if np.isnan(sigrednoise): # 如果没有检测到红噪声，则设置为零
sigrednoise = 0
sigwhite = np.sqrt(sig0**2 - sigrednoise**2) # 计算白噪声
sigrednoise = np.sqrt(sig_binned**2 - sigwhite**2 / binmeasure) # 重新计算红噪声
if np.isnan(sigrednoise): # 如果没有检测到红噪声，则设置为零
sigrednoise = 0
beta[bin] = np.sqrt(sig0**2 + binmeasure * sigrednoise**2) / sig0 # 计算红噪声膨胀因子
print("White noise (ppm) : ", 1E6 * sigwhite) # 打印白噪声
print("Red noise (ppm) : ", 1E6 * sigrednoise) # 打印红噪声
print("Transit depth measured error (ppm) : ", 2E6 * result.params['rprs'].value * result.params['rprs'].stderr) # 打印过境深度测量误差
n = text_file.write("\nWhite noise (ppm) : " + str(1E6 * sigwhite)) # 写入白噪声
n = text_file.write("\nRed noise (ppm) : " + str(1E6 * sigrednoise)) # 写入红噪声
n = text_file.write("\nTransit depth measured error (ppm) : " + str(2E6 * result.params['rprs'].value * result.params['rprs'].stderr)) # 写入过境深度测量误差
text_file.close() # 关闭文件
depth[bin] = 1E6 * result.params['rprs'].value**2 # 计算过境深度
depth_err[bin] = 2E6 * result.params['rprs'].value * result.params['rprs'].stderr # 计算过境深度误差
sig_r[bin] = sigrednoise * 1E6 # 更新红噪声
sig_w[bin] = sigwhite * 1E6 # 更新白噪声
# --------------------------------------------------------------------------
# ---------------------------------------------------------
# 将拟合光谱写入ascii文件
ascii_data = Table([wsd, werr, rprs, rerr, depth, depth_err, sig_w, sig_r, beta], names=['Wavelength Center (um)', 'Wavelength half-width (um)', 'Rp/Rs', 'Rp/Rs 1-sigma error', 'Transit Depth (ppm)', 'Transit Depth error', 'Sigma_white (ppm)', 'Sigma_red (ppm)', 'Beta Rednoise Inflation factor']) # 创建表格数据
ascii.write(ascii_data, save_directory + 'JWST_NIRSpec_Prism_fit_transmission_spectra.csv', format='csv', overwrite=True) # 写入CSV文件
# ---------------------------------------------------------
msize = plt.rcParams['lines.markersize'] ** 2. # 默认标记大小
# 绘制数据模型
# 绘图
plt.rcParams['figure.figsize'] = [10.0, 7.0] # 设置图形尺寸
msize = plt.rcParams['lines.markersize'] ** 2. # 默认标记大小
fig = plt.figure(constrained_layout=True) # 创建图形
gs = fig.add_gridspec(3, 1, hspace=0.00, wspace=0.00) # 创建网格规格
ax1 = fig.add_subplot(gs[0:2, :]) # 添加子图
ax1.scatter(x, y / p['f0'].value, s=msize * 0.75, linewidth=1, zorder=0, alpha=0.5, marker='+', edgecolors='blue') # 绘制散点图
ax1.plot(x, final_model / p['f0'].value, linewidth=1, color='orange', alpha=0.8, zorder=15) # 在数据上叠加过境模型
ax1.xaxis.set_ticklabels([]) # 隐藏X轴标签
plt.ylabel('Relative Flux') # 设置Y轴标签
plt.title(r'Time-series Transit Light Curve $\lambda=$[' + str(wave1) + ':' + str(wave2) + r'] $\mu$m') # 设置标题
ax1.xaxis.set_major_locator(ticker.MultipleLocator(0.01)) # 设置X轴主刻度
ax1.xaxis.set_minor_locator(ticker.MultipleLocator(0.005)) # 设置X轴次刻度
ax1.yaxis.set_major_locator(ticker.MultipleLocator(0.002)) # 设置Y轴主刻度
ax1.yaxis.set_minor_locator(ticker.MultipleLocator(0.001)) # 设置Y轴次刻度
yplot = y / np.mean(y[outtransit]) # 归一化Y数据
plt.ylim(yplot.min() * 0.999, yplot.max() * 1.001) # 设置Y轴范围
plt.xlim(bjd.min() - 0.001, bjd.max() + 0.001) # 设置X轴范围
fig.add_subplot(ax1) # 添加子图
# 残差
ax2 = fig.add_subplot(gs[2, :]) # 添加残差子图
ax2.scatter(x, residppm, s=msize * 0.75, linewidth=1, alpha=0.5, marker='+', edgecolors='blue', zorder=0) # 绘制残差散点图
wsb, wsb_bin_edges, binnumber = stats.binned_statistic(bjd, residppm, bins=256) # 进行分箱统计
plt.scatter(wsb_bin_edges[1:], wsb, linewidth=2, alpha=0.75, facecolors='orange', edgecolors='none', marker='o', zorder=25) # 绘制分箱后的残差
plt.xlabel('Barycentric Julian Date (days)') # 设置X轴标签
plt.ylabel('Residual (ppm)') # 设置Y轴标签
plt.plot([bjd.min(), bjd.max()], [0, 0], color='black', zorder=10) # 绘制零线
plt.plot([bjd.min(), bjd.max()], [sigma, sigma], linestyle='--', color='black', zorder=15) # 绘制标准差上限
plt.plot([bjd.min(), bjd.max()], [-sigma, -sigma], linestyle='--', color='black', zorder=20) # 绘制标准差下限
ax2.xaxis.set_major_locator(ticker.MultipleLocator(0.01)) # 设置X轴主刻度
ax2.xaxis.set_minor_locator(ticker.MultipleLocator(0.005)) # 设置X轴次刻度
yplot = y / np.mean(y[outtransit]) # 归一化Y数据
plt.xlim(bjd.min() - 0.001, bjd.max() + 0.001) # 设置X轴范围
fig.add_subplot(ax2) # 添加子图
# 保存
pp = PdfPages(save_directory + 'JWST_NIRSpec_Prism_fit_light_curve_bin' + str(bin) + '_lightcurve.pdf') # 创建PDF文件
plt.savefig(pp, format='pdf') # 保存图形为PDF
pp.close() # 关闭PDF文件
plt.clf() # 清空当前图形
# 绘制系统校正后的光曲线
b0 = p['a_Rs'].value * np.sqrt((np.sin(phase * 2 * np.pi)) ** 2 + (p['cosinc'].value * np.cos(phase * 2 * np.pi)) ** 2) # 计算b0
intransit = (b0 - p['rprs'].value < 1.0E0).nonzero() # 找到过境点
light_curve = b0 / b0 # 初始化光曲线
mulimb0, mulimbf = occultnl(p['rprs'].value, c1, c2, c3, c4, b0[intransit]) # 计算边缘暗化
light_curve[intransit] = mulimb0 # 更新过境点的光曲线
fig, axs = plt.subplots() # 创建新的图形
plt.scatter(x, light_curve + resid, s=msize * 0.75, linewidth=1, zorder=0, alpha=0.5, marker='+', edgecolors='blue') # 绘制光曲线
plt.xlabel('BJD') # 设置X轴标签
plt.ylabel('Relative Flux') # 设置Y轴标签
plt.plot(x, light_curve, linewidth=2, color='orange', alpha=0.8, zorder=15) # 在数据上叠加光曲线
pp = PdfPages(save_directory + 'JWST_NIRSpec_Prism_fit_light_curve_bin' + str(bin) + '_corrected.pdf') # 创建PDF文件
plt.savefig(pp, format='pdf') # 保存图形为PDF
pp.close() # 关闭PDF文件
plt.clf() # 清空当前图形
plt.close('all') # 关闭所有图形
# --------------------------------------------------------------------------
k = k + wk # 步进到下一个波长区间
print('** 现在可以查看输出PDF文件在 ', save_directory) # 打印输出目录
绘制测量的系外行星传输光谱与注入光谱对比#
# --------------------------------------------------------------------------
# 加载注入的传输光谱以与恢复值进行比较
# 下载注入的光谱数据
fn_tm = download_file(r'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/transit_spectroscopy_notebook/trans-iso_GJ-0436_0669.0_+0.0_0.56_0001_0.00_model.NIRSpec_PRISM.txt')
# 移动下载的文件到指定的保存目录
destld = shutil.move(fn_tm, save_directory+'trans-iso_GJ-0436_0669.0_+0.0_0.56_0001_0.00_model.NIRSpec_PRISM.txt')
# 打开保存的光谱文件
f = open(save_directory+'trans-iso_GJ-0436_0669.0_+0.0_0.56_0001_0.00_model.NIRSpec_PRISM.txt', 'r')
# 从文件中读取数据
data = np.genfromtxt(f, delimiter=' ')
# 提取波长和光谱数据
model_ws = data[:, 0] # 波长
model_spec = data[:, 1] # 光谱强度
# 读取拟合的过境深度数据
data = ascii.read(save_directory+'JWST_NIRSpec_Prism_fit_transmission_spectra.csv', format='csv')
# 提取相关数据列
wsd = data['Wavelength Center (um)'] # 波长中心
werr = data['Wavelength half-width (um)'] # 波长半宽
rprs = data['Rp/Rs'] # 行星半径与恒星半径之比
rerr = data['Rp/Rs 1-sigma error'] # 行星半径与恒星半径之比的1σ误差
beta = data['Beta Rednoise Inflation factor'] # 红噪声膨胀因子
# 绘图
fig, axs = plt.subplots() # 创建子图
# 绘制注入的光谱
plt.plot(model_ws, model_spec**2*1E6, linewidth=2, zorder=0, color='blue', label='Injected Spectra') # 在数据上叠加过境模型
# 绘制带误差条的恢复光谱（考虑红噪声）
plt.errorbar(wsd, rprs**2*1E6, xerr=werr, yerr=2*rerr*rprs*1E6*beta, fmt='o', zorder=5, alpha=0.4, color='orange', label=r'Recovered Spectra with $\sigma_r$')
# 绘制不考虑红噪声的恢复光谱
plt.errorbar(wsd, rprs**2*1E6, xerr=werr, yerr=2*rerr*rprs*1E6, fmt='o', zorder=10, color='orange', label='Recovered Spectra')
# 设置坐标轴标签
plt.xlabel(r'Wavelength ($\mu$m)') # 波长
plt.ylabel('Transit Depth ($R_p/R_s$)$^2$ (ppm)') # 过境深度
# 设置y轴主刻度和次刻度
axs.yaxis.set_major_locator(ticker.MultipleLocator(200)) # 主刻度每200
axs.yaxis.set_minor_locator(ticker.MultipleLocator(100)) # 次刻度每100
# 设置x轴主刻度和次刻度
axs.xaxis.set_major_locator(ticker.MultipleLocator(0.5)) # 主刻度每0.5
axs.xaxis.set_minor_locator(ticker.MultipleLocator(0.1)) # 次刻度每0.1
# 添加注释文本
axs.text(3.3, 6850, 'CH$_4$') # CH4注释
axs.text(4.25, 6700, 'CO') # CO注释
axs.text(2.3, 6750, 'CH$_4$') # CH4注释
axs.text(2.75, 6550, 'H$_2$O') # H2O注释
# 设置y轴和x轴的显示范围
plt.ylim(5700, 7000) # y轴范围
plt.xlim(0.9, 5.25) # x轴范围
# 显示图例
plt.legend(loc='lower right')
# 显示图形
plt.show()
# 清空当前图形
plt.clf()
总体而言，注入的过境深度在光谱中得到了很好的恢复，H\(_2\)O 和 CH\(_4\) 等特征易于被检测到。在 3.5 \(\mu\)m 以上的波长数据点存在一些偏移，这可能是由于在 CV3 测试期间，低温窗口上积累的冰导致 CO\(_2\) 或 H\(_2\)O 吸收特征的变化。这些波长显示出一些与时间相关的噪声 (\(\sigma_r\)) 增加，这里已进行了测量，图中的误差也显示了包含该误差的过境深度。
来自地面测试的精度非常令人鼓舞，最佳测量的 bin（发生在光谱中计数率高且干净的部分）在仅 2 小时的数据中达到了近光子极限的过境深度，测量精度约为 30 ppm，并且时间相关噪声 (\(\sigma_r\)) 最小。
为了获得更稳健的误差估计，实际上这里执行的最小二乘法将被 MCMC 例程替代。此外，使用实际的过境数据，过境拟合参数（例如 \(i\)、\(a/R_{star}\)、T\(_0\)）也必须首先进行拟合，因为它们可能会与文献中高精度过境光曲线的估计值有所不同，而 JWST 将提供这些数据。