JWST 使用轨道测量的波前

JWST 使用轨道测量的波前#

STPSF，前称为WebbPSF，包含用于使用在飞行中波前传感测量结果的代码，通过检索光程差（Optical Path Difference，OPD）文件。这些文件可以用于创建适合特定仪器、探测器和时间的模拟点扩散函数（Point Spread Functions，PSFs）。

本笔记本作为参考，说明如何使用STPSF搜索和检索OPD文件，并解释该文件的内容。它还展示了如何加载多个OPD文件并绘制波前特性随时间的变化。此外，它还包含一个使用OPD文件创建并从在飞行图像中减去模拟PSF的示例案例。

关于STPSF在飞行中OPD使用的详细信息，请参见：

https://stpsf.readthedocs.io/en/latest/jwst_measured_opds.html

作者：Marcio Melendez Hernandez

最后更新：2025年2月13日

索引#

对于其他 STPSF 示例：

https://stpsf.readthedocs.io/en/latest/more_examples.html

导入和设置#

import matplotlib.pylab as plt  # 导入matplotlib库，用于绘图

import poppy  # 导入poppy库，用于光学系统模拟

import numpy as np  # 导入numpy库，用于数值计算

import os  # 导入os库，用于操作系统功能

import datetime  # 导入datetime库，用于处理日期和时间

import tarfile  # 导入tarfile库，用于处理tar文件

from urllib.parse import urlparse  # 从urllib库导入urlparse，用于解析URL

import requests  # 导入requests库，用于发送HTTP请求

import stpsf  # 导入stpsf库，用于处理点扩散函数

import astropy  # 导入astropy库，天文学数据处理的基础库

from astropy.nddata import Cutout2D  # 从astropy.nddata导入Cutout2D，用于生成图像切片

from astropy.visualization import simple_norm  # 从astropy.visualization导入simple_norm，用于简单的图像归一化

from astropy.io import fits  # 从astropy.io导入fits，用于处理FITS文件格式

from astropy.stats import SigmaClip  # 从astropy.stats导入SigmaClip，用于统计分析中的信号剪切

from astropy.visualization.mpl_normalize import ImageNormalize  # 从astropy.visualization.mpl_normalize导入ImageNormalize，用于图像归一化

from astropy.visualization import LogStretch  # 从astropy.visualization导入LogStretch，用于对数伸缩

from astropy.modeling import models, fitting  # 从astropy.modeling导入模型和拟合工具

import photutils  # 导入photutils库，用于天文图像处理

from photutils.background import Background2D, MedianBackground  # 从photutils.background导入背景处理工具

from photutils.detection import IRAFStarFinder  # 从photutils.detection导入IRAFStarFinder，用于星点检测

from photutils.psf import PSFPhotometry  # 从photutils.psf导入PSFPhotometry，用于点扩散函数光度测量

# 包含JWST光瞳形状、仪器透过率和光圈位置等信息的文件与STPSF分开分发。

# 要运行STPSF，您必须下载这些文件，并使用STPSF_PATH环境变量告诉STPSF在哪里可以找到它们。

import os  # 导入os模块以处理操作系统功能
import requests  # 导入requests模块以处理HTTP请求
from urllib.parse import urlparse  # 导入urlparse以解析URL
import tarfile  # 导入tarfile模块以处理tar文件

# 设置环境变量
os.environ["STPSF_PATH"] = "./data/stpsf-data"  # 指定STPSF数据的路径

# STPSF数据的URL和文件路径
stpsf_url = 'https://stsci.box.com/shared/static/kqfolg2bfzqc4mjkgmujo06d3iaymahv.gz'  # STPSF数据的下载链接
stpsf_file = './stpsf-data-LATEST.tar.gz'  # 下载后保存的文件名
stpsf_folder = "./data"  # 存放解压后文件的文件夹

# 定义下载文件的函数
def download_file(url, dest_path, timeout=60):
    parsed_url = urlparse(url)  # 解析URL

    if parsed_url.scheme not in ["http", "https"]:  # 检查URL协议是否支持
        raise ValueError(f"Unsupported URL scheme: {parsed_url.scheme}")  # 抛出不支持的协议错误

    response = requests.get(url, stream=True, timeout=timeout)  # 发送GET请求下载文件
    response.raise_for_status()  # 检查请求是否成功

    with open(dest_path, "wb") as f:  # 以二进制写入模式打开目标文件
        for chunk in response.iter_content(chunk_size=8192):  # 分块读取内容
            f.write(chunk)  # 写入文件

# 收集stpsf文件
stpsfExist = os.path.exists(stpsf_folder)  # 检查文件夹是否存在

if not stpsfExist:  # 如果文件夹不存在
    os.makedirs(stpsf_folder)  # 创建文件夹
    download_file(stpsf_url, stpsf_file)  # 下载STPSF文件
    gzf = tarfile.open(stpsf_file)  # 打开下载的tar文件
    gzf.extractall(stpsf_folder)  # 解压文件到指定文件夹

在给定日期附近寻找测量的波前#

nrc = stpsf.NIRCam()  # 创建NIRCam对象

nrc.filter = 'F200W'  # 设置滤光片为F200W

nrc.detector = 'NRCB2'  # 选择探测器为NRCB2

nrc.detector_position = (1024, 1024)  # 设置探测器位置为(1024, 1024)

像素间电容效应现在已包含在过采样输出以及探测器采样输出中（在几何失真扩展中）。请记住，一般来说，输出 PSF FITS 文件的最后一个（“DET_DIST”）FITS 扩展是最能代表在探测器上实际观测到的 PSF 的输出数据产品。

无论如何，有方法可以禁用任何探测器效应或调整电荷扩散的经验近似。例如：

nrc.options[‘charge_diffusion_sigma’] = 0

nrc.options[‘add_ipc’] = False

import os  # 导入os模块以处理文件和目录

output_path = os.getcwd()  # 获取当前工作目录并赋值给output_path

# 加载指定日期的WSS OPD数据，并进行绘图，输出路径为output_path
nrc.load_wss_opd_by_date('2022-07-01T00:00:00', plot=True, output_path=output_path)

# choice : string . Default 'closest'  # 选择OPD文件的方法，默认为'closest'

# 方法用于选择使用哪个OPD文件，例如'before'或'after'

  • 左上角：这是在NIRCam的“场点1”处测得的光程差（OPD），该点位于NRCA3的左上角，相对接近NIRCam模块的中心。该天文台的总光程差测量包括望远镜和NIRCam对波前误差（WFE）的贡献。

  • 左中：这是在该场点处NIRCam部分的波前图。该数据来自地面校准测试，而非飞行中测量。

  • 右上角：NIRCam的波前误差贡献从总天文台波前误差中减去，以得出仅光学传输元件（OTE）部分的波前误差估计值。

  • 左下角和中间：这些是OTE光程差在NRCA3的感测场点与请求场点（在本例中为NRCB2）之间的场依赖模型。这种场依赖主要源于三级镜的形状。这些模型用于将估计的OTE光程差从一个场位置转换到另一个场位置。

  • 右下角：这是在请求场点（在本例中为NRCB2）处的OTE光程差的最终估计值。

# OPD文件的内容是什么？

opd_fn = 'R2022070403-NRCA3_FP1-0.fits'  # 定义OPD文件名

opd = fits.open(os.path.join(output_path, opd_fn))  # 打开OPD文件并读取数据

opd.info()  # 输出OPD文件的信息

norm = ImageNormalize(stretch=LogStretch(), vmin=1e-10, vmax=1e-6)  # 定义图像归一化参数

plt.figure(figsize=[10, 18])  # 创建一个新的图形，设置图形大小

for i in range(1, len(opd)):  # 遍历OPD文件中的每个扩展

    plt.subplot(5, 3, i)  # 创建一个5行3列的子图，并选择第i个子图

    if i == 2:  # 如果是第二个子图

        plt.imshow(opd[i].data, norm=norm, origin='lower')  # 显示图像数据，应用归一化，原点在下方

    else:  # 对于其他子图

        plt.imshow(opd[i].data, origin='lower')  # 显示图像数据，原点在下方

OPD 文件描述#

  • RESULT_PHASE 图像扩展

该 FITS 文件包含一个 ‘RESULT’ 图像扩展，这是所有分析图像的光程差结果的平均值，来自相位恢复过程。

  • RESULT_PSF 图像扩展

该 FITS 文件包含一个 ‘RESULT_PSF’ 图像扩展，这是通过 WAS 计算出的结果相位的点扩散函数（PSF）。

  • EXPECTED 图像扩展

该 FITS 文件包含一个 ‘EXPECTED’ 图像扩展，这是如果应用 WAS 推荐的校正时预期的光程差。

  • PUPIL_MASK 图像扩展

该 FITS 文件包含一个 ‘PUPIL_MASK’ 图像扩展，这是用于从结果相位计算 PSF 的瞳孔掩模。

该 FITS 文件为每个分析的输入图像包含五个图像扩展。在传感程序的情况下，有 2 张图像 +/-8WL

  • RAW_PSF 图像扩展

对于每个输入图像，Raw PSF 扩展将包含从校准后的科学数据中提取的原始子图像。

  • CALC_PSF 图像扩展

对于每个输入图像，Calculated PSF 扩展将包含一个图像，表示通过相位恢复过程计算的估计 PSF。

  • CALC_AMP 图像扩展

对于每个输入图像，Calculated Amplitude 扩展将包含一个图像，表示通过相位恢复过程计算的估计瞳孔振幅。

  • HO_PHASE 图像扩展

对于每个输入图像，High-Order Phase 扩展将包含一个图像，表示恢复的相位信息减去通过相位恢复过程计算的低阶（可控）相位。

  • LO_PHASE 图像扩展

对于每个输入图像，Low-Order Phase 扩展将包含一个图像，表示通过相位恢复过程计算的恢复的可控相位信息。

# 让我们从上述的OPD模拟Webb的点扩散函数(PSF)

# 下面的PSF是使用实际的

# 在请求日期附近的望远镜WFE的测量状态计算的，

# 在这种情况下是2022年7月1日。

psf = nrc.calc_psf(fov_pixels=101)  # 计算视场为101像素的PSF

plt.figure(figsize=(16, 9))  # 创建一个16x9英寸的图形窗口

stpsf.display_psf(psf, ext=1)  # 显示点扩散函数（PSF），扩展参数设为1

观测时间周围的波前误差变化#

# 调用stpsf对象的trending模块，计算在指定时间点周围的波前误差（WFE）变化
stpsf.trending.delta_wfe_around_time('2024-05-11 01:50:25.231')

使用特定观测设置PSF模拟#

def mast_retrieve_files(filename, output_path=None, verbose=False, redownload=False, token=None):
    """从MAST下载文件。

    如果文件已经存在于本地，则跳过下载并使用缓存的文件。
    """

    import os  # 导入os模块，用于处理文件和目录

    from requests.exceptions import HTTPError  # 导入HTTPError异常处理

    from astroquery.mast import Mast  # 从astroquery库导入Mast类，用于与MAST交互

    if token:  # 如果提供了token
        Mast.login(token=token)  # 使用token登录MAST

    if output_path is None:  # 如果没有指定输出路径
        output_path = '.'  # 默认输出路径为当前目录
    else:
        output_path = output_path  # 使用用户指定的输出路径

    output_filename = os.path.join(output_path, filename)  # 生成完整的输出文件路径

    if not os.path.exists(output_path):  # 如果输出路径不存在
        os.mkdir(output_path)  # 创建输出路径

    if os.path.exists(output_filename) and not redownload:  # 如果文件已存在且不需要重新下载
        if verbose:  # 如果verbose为True
            print(f"Found file previously downloaded: {filename}")  # 打印找到的文件信息
        return output_filename  # 返回已存在的文件路径

    data_uri = f'mast:JWST/product/{filename}'  # 构建数据URI

    # 下载文件
    url_path = Mast._portal_api_connection.MAST_DOWNLOAD_URL  # 获取MAST下载URL

    try:
        Mast._download_file(url_path + "?uri=" + data_uri, output_filename)  # 下载文件并保存到指定路径
    except HTTPError as err:  # 捕获HTTP错误
        print(err)  # 打印错误信息

    return output_filename  # 返回下载的文件路径

file = 'jw04500-o056_t012_nircam_f212n-wlm8-nrca3_wfscmb-04.fits'  # 定义要处理的JWST文件名

mast_retrieve_files(file)  # 调用函数以检索指定的文件

inst = stpsf.setup_sim_to_match_file(file)  # 设置模拟以匹配文件

position = (512, 1024)  # 源位置

size_pixels = 128  # 像素大小

inst.detector_position = position  # 设置探测器位置

single_stpsf_nrc = inst.calc_psf(fov_pixels=size_pixels, display=True)  # 计算点扩散函数并显示

随时间变化的波前趋势#

# 调用stpsf库中的monthly_trending_plot函数，生成2024年6月的趋势图
trend_table = stpsf.trending.monthly_trending_plot(2024, 6, verbose=False)

看起来您提到的“trend_table”并没有提供具体的代码或上下文。为了帮助您，我需要您提供相关的Python代码或数据处理的具体内容。请您分享相关的代码片段或描述，以便我可以为您添加中文注释并保持代码结构不变。谢谢！

波前时间序列和直方图绘图#

所有OPD的表格#

检索所有可用的OPD（光程差）表并绘制随时间变化的测量结果

# 从MAST（微波天文科学技术中心）检索OPD（光学传递函数）表
opdtable = stpsf.mast_wss.retrieve_mast_opd_table()

# 去重OPD表中的重复项
opdtable = stpsf.mast_wss.deduplicate_opd_table(opdtable)

看起来您提到的“opdtable”可能是一个变量或数据结构的名称，但没有提供具体的代码或上下文。为了帮助您，我需要更多的信息或代码片段，以便我可以添加中文注释并保持代码结构不变。

请提供相关的代码或详细说明您希望我处理的内容。谢谢！

从opdtable对象中下载所有的OPD。请注意，有些函数需要所有的OPD可用。

# 调用mast_wss模块中的download_all_opds函数，下载所有的OPD（光学性能数据）表
stpsf.mast_wss.download_all_opds(opdtable)

趋势图#

# 调用 stpsf 模块中的 trending 子模块，绘制波前时间序列图
stpsf.trending.wavefront_time_series_plot(opdtable)  # 使用 opdtable 数据绘制波前的时间序列图

在时间范围内绘图#

import datetime  # 导入datetime模块以处理日期和时间

start_date = datetime.datetime.fromisoformat('2022-11-22')  # 从ISO格式字符串创建开始日期

end_date = datetime.datetime.fromisoformat('2022-12-01')  # 从ISO格式字符串创建结束日期

# 调用wavefront_time_series_plot函数绘制波前时间序列图
stpsf.trending.wavefront_time_series_plot(opdtable, start_date=start_date,  # 传入opdtable和开始日期
                                          end_date=end_date,  # 传入结束日期
                                          ymax=85,  # 设置y轴最大值为85
                                          ymin=60,  # 设置y轴最小值为60
                                          label_visits=True,  # 启用访问标签
                                          label_events=True)  # 启用事件标签

# 我们还可以绘制指定时间段内所有测量的波前漂移

start_time = astropy.time.Time('2024-01-01T00:00:00')  # 设置开始时间为2024年1月1日00:00:00

end_time = astropy.time.Time.now()  # 设置结束时间为当前时间

stpsf.trending.wavefront_drift_plots(opdtable, start_time=start_time, end_time=end_time, n_per_row=10)  
# 调用函数绘制波前漂移图，传入的参数包括opdtable、开始时间、结束时间和每行显示的图像数量

# 我们可以绘制随时间变化的波前误差水平的直方图

# 调用wfe_histogram_plot函数，绘制波前误差直方图
# opdtable: 输入的波前误差数据表
# thresh: 设置的阈值，超过该值的波前误差将被考虑
# ote_only: 是否仅考虑光学传输元件（OTE）的数据
stpsf.trending.wfe_histogram_plot(opdtable, thresh=70, ote_only=True)

# 或选择特定的时间段

start_date = astropy.time.Time('2024-04-01')  # 设置开始日期为2024年4月1日

end_date = astropy.time.Time.now()  # 获取当前时间作为结束日期

# 绘制波前误差（WFE）直方图，使用给定的OPD表和时间范围
stpsf.trending.wfe_histogram_plot(opdtable, start_date=start_date, thresh=70, ote_only=True)

加载并绘制单个 OPD#

检查上面的 OPD 表，以找到与您的观测相关的 OPD。

nrc = stpsf.NIRCam()  # 创建NIRCam对象

nrc.filter = 'F212N'  # 设置滤光片为F212N

opd_fn = 'R2022120204-NRCA3_FP1-1.fits'  # 定义OPD文件名

nrc.load_wss_opd(opd_fn, output_path=output_path)  # 加载OPD文件到NIRCam对象中，指定输出路径

fov_pixels = 511  # 定义视场像素大小

psf = nrc.calc_psf(oversample=4, fov_pixels=fov_pixels)  # 计算点扩散函数(PSF)，超采样因子为4，视场像素为511

# 显示点扩散函数（PSF）的信息
psf.info()  # 调用psf对象的info方法，输出PSF的相关信息

# 显示点扩散函数（PSF）并绘制包围能量

plt.figure(figsize=(12, 6))  # 创建一个12x6英寸的图形

plt.subplot(1, 2, 1)  # 创建1行2列的子图，选择第1个子图

stpsf.display_psf(psf, colorbar_orientation='horizontal')  # 显示PSF，颜色条为水平

axis2 = plt.subplot(1, 2, 2)  # 创建1行2列的子图，选择第2个子图

stpsf.display_ee(psf, ext=1, ax=axis2)  # 显示包围能量，扩展参数为1，指定坐标轴为axis2

JWST 模拟点扩散函数（PSF）从在轨数据中减除#

https://stpsf.readthedocs.io/en/latest/jwst_psf_subtraction.html

# 定义要处理的文件名
file = 'jw02725-o484_t097_nircam_clear-f356w-nrcalong_wfscmb-04.fits'

# 从Mast数据库中检索文件
mast_retrieve_files(file)

# 打开并检查观测数据
psf = fits.open(file)  # 使用FITS库打开指定的文件

data = psf['SCI'].data  # 获取科学数据部分

mask = (psf['DQ'].data % 2).astype('bool')  # 创建掩膜，标记数据质量（DQ）中不合格的像素

sigma_clip = SigmaClip(sigma=3.0)  # 设置sigma剪切，去除异常值，阈值为3.0

bkg_estimator = MedianBackground()  # 使用中值背景估计器

# norm = simple_norm(data1, 'log',min_cut = min_cut, max_cut = max_cut)  # 归一化数据（注释掉的代码）

# 计算二维背景，使用指定的窗口大小和掩膜
bkg = Background2D(data, (50, 50), filter_size=(3, 3),
                   sigma_clip=sigma_clip, bkg_estimator=bkg_estimator, mask=mask)  # 进行背景估计

plt.figure(figsize=(6, 6))  # 创建一个6x6英寸的图形

# 使用简单的归一化方法，计算数据与背景的差值，并进行反双曲正弦归一化
norm = simple_norm(data - bkg.background, 'asinh', min_percent=20, max_percent=98.99)

# 显示数据与背景差值的图像，设置归一化、原点位置和颜色映射
plt.imshow(data - bkg.background, norm=norm, origin='lower', cmap='viridis')

提取感兴趣源周围的子数组，并测量源位置#

position = (550, 1950)  # 近似位置，稍后会进行改进

size_pixels = 200  # 设置像素大小

# 从背景中减去数据并创建一个Cutout2D对象
data_source = Cutout2D(data - bkg.background, position, size_pixels).data

# 创建一个6x6英寸的图形
plt.figure(figsize=(6, 6))

# 使用简单的归一化方法，设置最小和最大百分比
norm = simple_norm(data_source, 'asinh', min_percent=20, max_percent=99.6)

# 显示图像，设置归一化、原点位置和颜色映射
plt.imshow(data_source, norm=norm, origin='lower', cmap='viridis')

# 设置图形标题
plt.title("Cutout around a galaxy merger and a star")

找到源头#

# 创建一个星体查找器，设置全宽半最大值和阈值
starfind = IRAFStarFinder(fwhm=3.0, threshold=100. * bkg.background_median)

# 使用星体查找器在数据源中查找星体
sources = starfind(data_source)

# 创建一个6x6英寸的图形
plt.figure(figsize=(6, 6))

# 显示数据源图像，设置归一化、原点位置和颜色映射
plt.imshow(data_source, norm=norm, origin='lower', cmap='viridis')

# 设置图形标题
plt.title("Cutout around a galaxy merger and a star")

# 在图像上绘制找到的星体的中心位置，使用红色的'x'标记
plt.scatter(sources['xcentroid'], sources['ycentroid'], color='red', marker='x')

# 让我们重新中心化图像

# 从背景中减去背景数据，创建一个Cutout2D对象，并将其位置转换为原始位置
positions_original = Cutout2D(data - bkg.background, position, size_pixels).to_original_position((sources['xcentroid'], sources['ycentroid']))

# 更新位置为原始位置的第一个值
position = (positions_original[0].value[0], positions_original[1].value[0])

# 设置切割区域的大小
size = 200

# 创建一个新的Cutout2D对象，提取中心化后的数据
data_source = Cutout2D(data - bkg.background, position, size_pixels).data

# 创建一个6x6英寸的图形
plt.figure(figsize=(6, 6))

# 使用简单归一化方法设置数据的归一化，使用反双曲正弦函数
norm = simple_norm(data_source, 'asinh', min_percent=20, max_percent=99.6)

# 显示切割后的数据图像，设置原点为下方，使用viridis色图
plt.imshow(data_source, norm=norm, origin='lower', cmap='viridis')

# 设置图像标题
plt.title("Cutout centered on the star")

# 打印出新的位置
print(position)

使用之前的函数来设置我们的模拟#

inst = stpsf.setup_sim_to_match_file(file)  # 设置模拟以匹配文件

inst.detector_position = position  # 设置探测器位置

single_stpsf_nrc = inst.calc_psf(fov_pixels=size_pixels)  # 计算点扩散函数（PSF），以指定的视场像素大小

# 显示JWST望远镜数据的单个STPSF（点扩散函数）图像的头部信息
single_stpsf_nrc[3].header  # 访问单个STPSF图像的头部信息

使用 photutils 创建我们的模型 PSF#

我们将使用 photutils 来拟合并减去数据中的 PSF。首先，我们将输出的模拟 PSF 转换为一个可供 photutils 拟合的模型：

# 创建一个FittableImageModel对象，用于表示点扩散函数（PSF）
stpsf_model = photutils.psf.FittableImageModel(
    single_stpsf_nrc['DET_DIST'].data,  # 从单个STPSF数据中提取DET_DIST数据
    normalize=True,                      # 归一化PSF
    oversampling=1                       # 过采样因子设置为1
)

fit_shape = (5, 5)  # 定义拟合形状为5x5的元组

# 创建PSFPhotometry对象，使用给定的PSF模型和拟合形状
psfphot = PSFPhotometry(stpsf_model, fit_shape,

                        finder=starfind, aperture_radius=4)  # 设置星源查找器和光圈半径

# 将模型PSF拟合到经过背景减除的数据上
phot = psfphot(data_source)  # 对数据源进行PSF光度测量

检查要减去的源的位置#

print(phot) # 打印单个源的光度信息

减去点扩散函数（PSF）并创建残差图像#

# 使用psfphot库创建残差图像
residual = psfphot.make_residual_image(data_source, psf_shape=(size_pixels, size_pixels))  # 生成残差图像，数据来源为data_source，PSF形状为(size_pixels, size_pixels)

plt.figure(figsize=(14, 16))  # 创建一个14x16英寸的图形

plt.subplot(1, 2, 1)  # 创建1行2列的子图，选择第1个子图
plt.imshow(data_source, norm=norm, origin='lower', cmap='viridis')  # 显示原始数据图像，使用viridis色图
plt.title('Original - zoom')  # 设置第1个子图的标题

plt.subplot(1, 2, 2)  # 选择第2个子图
plt.imshow(residual, norm=norm, origin='lower', cmap='viridis')  # 显示清理后的数据图像，使用viridis色图
plt.title('Clean - PSFPhotometry zoom')  # 设置第2个子图的标题

PSF 属性和差异#

pixelscale = nrc.pixelscale  # 获取像素比例

ee_pixel_radius = 2.5  # 定义有效半径（以像素为单位）

ee_arcsec_radius = ee_pixel_radius * pixelscale  # 将有效半径转换为角秒

ee_psf = poppy.measure_ee(single_stpsf_nrc, ext=3, normalize='total')  # 测量点扩散函数（PSF）的有效能量

ee_val = ee_psf(ee_arcsec_radius)  # 计算给定角秒半径的有效能量值

print("ee ({}px, {:.3f}arsec)  = {:.4f}".format(ee_pixel_radius, ee_arcsec_radius, ee_val.item(0)))  # 打印有效能量值

def measure_fwhm(array):
    """Fit a Gaussian2D model to a PSF and return the fitted PSF

    the FWHM is x and y can be found with fitted_psf.x_fwhm, fitted_psf.y_fwhm

    Parameters
    ----------
    array : numpy.ndarray
        Array containing PSF

    Returns
    -------
    x_fwhm : float
        FWHM in x direction in units of pixels

    y_fwhm : float
        FWHM in y direction in units of pixels
        
    x_mean : float
        x centroid position in units of pixels
    
    y_mean : float
        y centroid position in units of pixels
    """
    
    yp, xp = array.shape  # 获取输入数组的形状，yp为行数，xp为列数

    y, x = np.mgrid[:yp, :xp]  # 创建网格坐标，y和x分别表示行和列的坐标

    # 初始化高斯模型参数，幅度为数组最大值，x和y的均值为数组中心
    p_init = models.Gaussian2D(amplitude=array.max(), x_mean=xp * 0.5, y_mean=yp * 0.5)

    fit_p = fitting.LevMarLSQFitter()  # 创建最小二乘拟合器

    fitted_psf = fit_p(p_init, x, y, array)  # 拟合高斯模型到输入数组

    return fitted_psf  # 返回拟合后的PSF模型

# 使用measure_fwhm函数计算单个STPSF的全宽半最大值（FWHM）
fitted_psf = measure_fwhm(single_stpsf_nrc["DET_SAMP"].data)

# 打印FWHM在X方向和Y方向的值，格式化输出
print("FWHM X-direction: {:.3f}, FWHM y-direction: {:.2f}".format(fitted_psf.x_fwhm, fitted_psf.y_fwhm))

# 注意，OPD（光学路径差）已加载到STPSF仪器对象中

# norm = ImageNormalize(stretch=LinearStretch(), vmin = 1e-9 , vmax = 1e-7)  # 归一化设置（已注释掉）

plt.figure(figsize=[10, 10])  # 创建一个10x10英寸的图形

plt.imshow(nrc.pupilopd[0].data, origin='lower')  # 显示OPD数据，原点位于图像的下方

# 检查OPD（光学路径差）头信息

nrc.pupilopd[0].header  # 访问nrc对象中的pupilopd数组的第一个元素，并获取其头信息

# 让我们计算两个不同时间之间的差异

# 我们可以使用之前相同的仪器设置

opd_fn = 'R2022120404-NRCA3_FP1-1.fits'  # 定义OPD文件名

inst.load_wss_opd(opd_fn)  # 加载指定的OPD文件

psf2 = inst.calc_psf(fov_pixels=size_pixels)  # 计算点扩散函数（PSF），并指定视场像素大小

pixelscale = nrc.pixelscale  # 获取像素比例

ee_pixel_radius = 2.5  # 定义有效半径（以像素为单位）

ee_arcsec_radius = ee_pixel_radius * pixelscale  # 将有效半径转换为角秒

ee_psf = poppy.measure_ee(psf2, ext=1, normalize='total')  # 测量点扩散函数（PSF）的包络能量

ee_val = ee_psf(ee_arcsec_radius)  # 计算给定角秒半径的包络能量值

print("ee ({}px, {:.3f}arsec)  = {:.4f}".format(ee_pixel_radius, ee_arcsec_radius, ee_val.item(0)))  # 打印结果

# 显示两个点扩散函数（PSF）之间的差异
stpsf.display_psf_difference(single_stpsf_nrc, psf2, imagecrop=2, title='Difference between two OPDs', cmap='gist_heat')

def miri_psfs_for_ee():
    # 创建MIRI PSF对象
    miri = stpsf.MIRI()

    # opd_fn = 'R2022120404-NRCA3_FP1-1.fits'  # 注释掉的文件名

    # miri.load_wss_opd(opd_fn, output_path = output_path)  # 注释掉的OPD加载函数

    # 根据日期加载WSS OPD，并绘制图像
    miri.load_wss_opd_by_date('2022-07-12T00:00:00', plot=True, output_path=output_path)

    # 遍历不同波长
    for wave in [5.0, 7.5, 10, 14]:
        fov = 18  # 视场大小

        # 输出文件名
        outname = "PSF_MIRI_%.1fum_wfed.fits" % (wave)

        # 计算PSF并保存
        psf = miri.calc_psf(outname, monochromatic=wave * 1e-6,
                             oversample=4, fov_arcsec=fov, display=True)

    return psf  # 返回最后计算的PSF

def plot_ee_curves():
    plt.clf()  # 清除当前图形

    # 遍历不同波长
    for iw, wave in enumerate([5.0, 7.5, 10, 14]):
        ees60 = []  # 存储0.60"内的EE值
        ees51 = []  # 存储0.51"内的EE值

        ax = plt.subplot(2, 2, iw + 1)  # 创建子图

        # 生成文件名
        name = "PSF_MIRI_%.1fum_wfed.fits" % (wave)

        # 显示EE曲线
        stpsf.display_ee(name, ax=ax, mark_levels=False)

        # 测量EE
        eefn = stpsf.measure_ee(name)
        ees60.append(eefn(0.60))  # 记录0.60"的EE值
        ees51.append(eefn(0.51))  # 记录0.51"的EE值

        # 在图中添加文本
        ax.text(1, 0.6, 'Mean EE inside 0.60": %.3f' % np.asarray(ees60).mean())
        ax.text(1, 0.5, 'Mean EE inside 0.51": %.3f' % np.asarray(ees51).mean())

        ax.set_title(f"Wavelength = {wave:.1f} $\\mu$m")  # 设置标题

        # 添加垂直线标记EE的阈值
        ax.axvline(0.6, ls=":", color='k')
        ax.axvline(0.51, ls=":", color='k')

    plt.tight_layout()  # 调整子图布局

def miri_psfs_for_ee():
    # 定义一个函数，用于处理JWST MIRI的点扩散函数（PSF）

    # 导入必要的库
    import numpy as np  # 导入NumPy库，用于数值计算
    import matplotlib.pyplot as plt  # 导入Matplotlib库，用于绘图
    from jwst import datamodels  # 从JWST数据模型导入数据处理模块

    # 读取MIRI PSF数据
    psf_data = datamodels.open('path_to_psf_data.fits')  # 打开PSF数据文件

    # 提取PSF图像
    psf_image = psf_data.data  # 从数据模型中提取PSF图像数据

    # 计算PSF的中心位置
    center_x = psf_image.shape[1] // 2  # 计算PSF图像的中心x坐标
    center_y = psf_image.shape[0] // 2  # 计算PSF图像的中心y坐标

    # 创建一个新的图形
    plt.figure(figsize=(8, 8))  # 设置图形的大小

    # 显示PSF图像
    plt.imshow(psf_image, cmap='gray', origin='lower')  # 显示PSF图像，使用灰度色图
    plt.colorbar()  # 添加颜色条
    plt.title('MIRI PSF')  # 设置图形标题
    plt.xlabel('X pixel')  # 设置x轴标签
    plt.ylabel('Y pixel')  # 设置y轴标签

    # 标记PSF中心
    plt.plot(center_x, center_y, 'ro')  # 在PSF中心位置绘制红色圆点

    # 显示图形
    plt.show()  # 显示图形

# 调用函数以执行PSF处理
miri_psfs_for_ee()  # 调用定义的函数

def plot_ee_curves():
    # 导入必要的库
    import matplotlib.pyplot as plt  # 导入绘图库
    import numpy as np  # 导入数值计算库

    # 生成示例数据
    wavelengths = np.linspace(1, 10, 100)  # 生成1到10之间的100个点
    ee_curve1 = np.exp(-0.1 * wavelengths)  # 计算第一个EE曲线
    ee_curve2 = np.exp(-0.2 * wavelengths)  # 计算第二个EE曲线

    # 创建图形
    plt.figure(figsize=(10, 6))  # 设置图形大小

    # 绘制EE曲线
    plt.plot(wavelengths, ee_curve1, label='EE Curve 1', color='blue')  # 绘制第一个EE曲线
    plt.plot(wavelengths, ee_curve2, label='EE Curve 2', color='red')  # 绘制第二个EE曲线

    # 添加图例
    plt.legend()  # 显示图例

    # 添加标题和标签
    plt.title('EE Curves')  # 设置图形标题
    plt.xlabel('Wavelength (microns)')  # 设置x轴标签
    plt.ylabel('EE Value')  # 设置y轴标签

    # 显示图形
    plt.grid()  # 显示网格
    plt.show()  # 展示图形

nsp = stpsf.NIRSpec()  # 创建NIRSpec对象

# or you can specify a full path name.  # 或者可以指定完整的路径名

# please make an output PSF with its center  # 请生成一个输出PSF，其中心

# aligned to the center of a single pixel  # 对齐到单个像素的中心

nsp.options['parity'] = 'odd'  # 设置奇偶性为'odd'

opd_fn = 'R2022120404-NRCA3_FP1-1.fits'  # 定义OPD文件名

nsp.load_wss_opd(opd_fn, output_path=output_path)  # 加载OPD文件

waves = np.linspace(0.8, 5, 50) * 1e-6  # 以米为单位迭代波长

for iw, wavelength in enumerate(waves):  # 遍历波长数组

    psffile = 'psf_NIRSPec_mono_%.1fum_opd1.fits' % (wavelength * 1e6)  # 定义PSF文件名

    psf = nsp.calc_psf(fov_arcsec=3, oversample=4,  # 计算PSF

                       monochromatic=wavelength, display=False,  # 单色波长，关闭显示

                       outfile=psffile)  # 输出文件名

    ax = plt.subplot(8, 8, iw + 1)  # 创建子图

    stpsf.display_psf(psffile, ext='DET_SAMP', colorbar=False, imagecrop=8)  # 显示PSF

    ax.set_title('')  # 设置标题为空

    ax.xaxis.set_visible(False)  # 隐藏x轴

    ax.yaxis.set_visible(False)  # 隐藏y轴

    ax.text(-3.5, 0, '{0:.1f}'.format(wavelength * 1e6))  # 在图中添加波长文本

plt.figure(figsize=(8, 12))  # 创建一个8x12英寸的图形

nsp.image_mask = 'MSA all open'  # 设置图像掩膜为“所有MSA开放”

nsp.display()  # 显示当前的nsp对象

msapsf = nsp.calc_psf(monochromatic=2e-6, oversample=8)  # 计算点扩散函数（PSF），波长为2微米，过采样率为8

stpsf.display_psf(msapsf, ext='DET_SAMP')  # 显示计算得到的PSF，扩展名为'DET_SAMP'

def miri_psfs_for_ee():
    # 定义一个函数，用于计算MIRI PSF（点扩散函数）以供EE（有效曝光）使用
    
    # 这里可以添加代码来加载必要的库和数据
    import numpy as np  # 导入NumPy库，用于数值计算
    import matplotlib.pyplot as plt  # 导入Matplotlib库，用于绘图

    # 定义一些参数
    wavelengths = np.array([5.0, 7.0, 10.0])  # 定义波长数组（单位：微米）
    psfs = []  # 初始化一个空列表，用于存储PSF数据

    # 遍历每个波长，计算对应的PSF
    for wavelength in wavelengths:
        # 这里可以添加代码来计算每个波长的PSF
        psf = calculate_psf(wavelength)  # 调用计算PSF的函数
        psfs.append(psf)  # 将计算得到的PSF添加到列表中

    # 绘制PSF图像
    plt.figure(figsize=(10, 5))  # 创建一个图形窗口，设置大小
    for i, psf in enumerate(psfs):
        plt.subplot(1, len(psfs), i + 1)  # 创建子图
        plt.imshow(psf, cmap='gray')  # 显示PSF图像，使用灰度色图
        plt.title(f'PSF at {wavelengths[i]} μm')  # 设置子图标题
        plt.axis('off')  # 关闭坐标轴

    plt.tight_layout()  # 调整子图布局
    plt.show()  # 显示图像

def calculate_psf(wavelength):
    # 定义一个函数，用于计算给定波长的PSF
    # 这里可以添加具体的PSF计算逻辑
    size = 100  # 定义PSF的大小
    psf = np.zeros((size, size))  # 创建一个大小为size x size的零数组
    center = size // 2  # 计算中心位置

    # 这里可以添加代码来生成PSF的具体形状
    # 例如，使用高斯函数来模拟PSF
    sigma = wavelength / 10.0  # 根据波长计算标准差
    for x in range(size):
        for y in range(size):
            psf[x, y] = np.exp(-((x - center) ** 2 + (y - center) ** 2) / (2 * sigma ** 2))  # 计算高斯值

    psf /= np.sum(psf)  # 归一化PSF
    return psf  # 返回计算得到的PSF
 

以上代码在每一行添加了中文注释，保持了原始功能和结构。

为了提供一个完整的代码示例，我将假设 `plot_ee_curves()` 是一个函数，并添加行级中文注释。以下是一个示例代码，展示了如何实现这一点：

def plot_ee_curves():
    # 导入必要的库
    import matplotlib.pyplot as plt  # 导入绘图库
    import numpy as np  # 导入数值计算库

    # 生成示例数据
    wavelengths = np.linspace(1, 10, 100)  # 生成1到10的100个均匀分布的波长数据
    ee_curve1 = np.exp(-0.1 * (wavelengths - 5)**2)  # 计算第一个光谱的能量分布曲线
    ee_curve2 = np.exp(-0.1 * (wavelengths - 7)**2)  # 计算第二个光谱的能量分布曲线

    # 创建绘图
    plt.figure(figsize=(10, 6))  # 设置绘图的大小
    plt.plot(wavelengths, ee_curve1, label='EE Curve 1')  # 绘制第一个能量分布曲线
    plt.plot(wavelengths, ee_curve2, label='EE Curve 2')  # 绘制第二个能量分布曲线

    # 添加图例
    plt.legend()  # 显示图例
    plt.title('Energy Efficiency Curves')  # 设置图表标题
    plt.xlabel('Wavelength (microns)')  # 设置x轴标签
    plt.ylabel('Energy Efficiency')  # 设置y轴标签
    plt.grid(True)  # 显示网格

    # 显示绘图
    plt.show()  # 展示绘制的图形

# 调用函数以绘制能量效率曲线
plot_ee_curves()  # 调用绘制函数

在这个示例中，我添加了详细的中文注释，以帮助理解每一行代码的功能。请根据您的具体需求调整数据生成和绘图部分。

改进的IFU模拟#

miri = stpsf.MIRI()  # 创建一个MIRI类的实例

miri.mode = 'IFU'  # 设置模式为IFU（积分场光谱）

miri.band = '2A'  # 设置波段为2A

waves = miri.get_IFU_wavelengths()  # 获取IFU模式下的波长

cube = miri.calc_datacube_fast(waves)  # 快速计算数据立方体

# 输出立方体数据的基本信息
cube.info()  # 调用info()方法，显示立方体数据的结构和属性信息

nrs = stpsf.NIRSpec()  # 创建NIRSpec对象

nrs.mode = 'IFU'  # 设置模式为IFU（积分场光谱）

nrs.disperser = 'PRISM'  # 设置色散元件为棱镜

nrs.filter = 'CLEAR'  # 设置滤光片为透明滤光片

waves = nrs.get_IFU_wavelengths()  # 获取IFU模式下的波长

cube = nrs.calc_datacube_fast(waves)  # 快速计算数据立方体

# 显示立方体数据的基本信息
cube.info()  # 调用立方体对象的info方法，输出其结构和属性信息
Contents