星系O2/O3表面亮度-半径数据适配函数及绘图实现问询

推荐的拟合函数 #

• 双Sersic函数：适配同时包含核球与盘结构的星系，比单Sersic函数更灵活，能覆盖更复杂的亮度分布
  公式：$I(r) = I_{e1} \exp\left(-b_{n1}\left(\left(\frac{r}{r_{e1}}\right)^{1/n1} - 1\right)\right) + I_{e2} \exp\left(-b_{n2}\left(\left(\frac{r}{r_{e2}}\right)^{1/n2} - 1\right)\right)$
  其中$b_n$是满足$\gamma(n, b_n) = \gamma(n, \infty)/2$的参数（$\gamma$为不完全伽马函数），可通过近似公式快速计算
• Nuker Law：针对具有陡峭核结构的椭圆星系设计，对中心亮度骤变的情况适配性更强
  公式：$I(r) = I_b \left(\frac{r}{r_b}\right)^{-\gamma} \left(1 + \left(\frac{r}{r_b}\right)^{\alpha\right)}{(\gamma-\beta)/\alpha}$
  参数说明：$r_b$为转折半径，$\gamma$是内区斜率，$\beta$是外区斜率，$\alpha$控制转折的尖锐程度

Python实现代码 #

依赖库导入 #

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.special import gammainc, gamma

定义拟合函数 #

1. 双Sersic函数

def b_n(n):
    # 快速计算Sersic公式的b_n参数（近似拟合公式）
    return 2*n - 1/3 + 4/(405*n) + 46/(25515*n**2) + 131/(1148175*n**3) - 2194697/(30690717750*n**4)

def double_sersic(r, Ie1, re1, n1, Ie2, re2, n2):
    term1 = Ie1 * np.exp(-b_n(n1) * ((r/re1)**(1/n1) - 1))
    term2 = Ie2 * np.exp(-b_n(n2) * ((r/re2)**(1/n2) - 1))
    return term1 + term2

2. Nuker Law函数

def nuker_law(r, Ib, rb, gamma, alpha, beta):
    return Ib * (r/rb)**(-gamma) * (1 + (r/rb)**alpha)**((gamma - beta)/alpha)

数据加载与拟合示例 #

将示例数据替换为你的真实星系O2、O3的半径、表面亮度及误差数据：

# 替换为你的真实数据
radii_o2 = np.array([0.1, 0.5, 1.0, 2.0, 3.0, 5.0, 10.0])
sb_o2 = np.array([1000, 500, 200, 80, 30, 8, 1])
sb_err_o2 = np.array([50, 30, 15, 8, 5, 2, 0.5])

# 双Sersic拟合：初始参数需根据数据估算（中心亮度、特征半径、Sersic指数）
p0_double = [500, 0.5, 1.0, 100, 2.0, 4.0]
params_double_o2, cov_double_o2 = curve_fit(double_sersic, radii_o2, sb_o2, p0=p0_double, sigma=sb_err_o2)

# Nuker Law拟合：初始参数需根据数据估算（转折亮度、转折半径、内外斜率）
p0_nuker = [200, 1.0, 1.5, 2.0, 4.0]
params_nuker_o2, cov_nuker_o2 = curve_fit(nuker_law, radii_o2, sb_o2, p0=p0_nuker, sigma=sb_err_o2)

绘图（x线性、y对数刻度+误差棒） #

# 生成拟合曲线的连续x点
r_fit = np.linspace(min(radii_o2), max(radii_o2), 100)
sb_fit_double = double_sersic(r_fit, *params_double_o2)
sb_fit_nuker = nuker_law(r_fit, *params_nuker_o2)

plt.figure(figsize=(8, 6))
# 绘制原始数据与误差棒
plt.errorbar(radii_o2, sb_o2, yerr=sb_err_o2, fmt='o', capsize=5, label='O2星系观测数据')
# 绘制拟合曲线
plt.plot(r_fit, sb_fit_double, label='双Sersic拟合', linestyle='--')
plt.plot(r_fit, sb_fit_nuker, label='Nuker Law拟合', linestyle=':')

# 坐标轴设置
plt.xlabel('半径')
plt.ylabel('表面亮度')
plt.yscale('log')
plt.xscale('linear')
plt.legend()
plt.grid(alpha=0.3)
plt.show()

实用提示 #

• 初始参数的猜测直接影响拟合收敛性，建议先通过观测数据的分布大致估算（比如中心亮度值、特征半径范围）
• 若拟合效果仍不理想，可尝试指数函数+Sersic函数的组合（盘结构用指数函数，核球用Sersic），或加入常数项处理背景噪声
• 针对O3星系的拟合，只需复用上述代码替换对应数据即可

内容的提问来源于stack exchange，提问作者Chinmaya