求助：Sheth-Tormen质量函数代码空图、零值及红移结果一致问题

你的代码出现空图、M>1e13时n值为0、红移z不影响结果的问题，核心是公式实现错误、红移依赖缺失、积分计算不准确这几个原因，以下是具体修正步骤：

1. 核心公式错误修正 #

（1）Sheth-Tormen质量函数公式写错 #

原代码中n_H函数的公式存在两处关键错误：

• 项(1+(a*delta_c**2/(s**2)))应该是(1+(s**2/(a*delta_c**2)))
• 多余了*(s/a/delta_c**2)这一项，这会导致大质量时s很小，该项趋近于0，最终n值为0

正确的Sheth-Tormen公式形式为：
$$
n(M,z) = A \sqrt{\frac{2a}{\pi}} \frac{\rho_m}{M} \left| \frac{d\ln\sigma}{d\ln M} \right| \left[ 1 + \left( \frac{\sigma^2}{a \delta_c^2} \right) \right] \exp\left( -\frac{a \delta_c^2}{2\sigma2} \right)
$$

（2）红移依赖缺失 #

原代码中delta_c固定为z=0时的1.686，但临界密度阈值随红移变化，正确的拟合公式为：
$$
\delta_c(z) \approx 1.686 \left( 1 + 0.0123 \ln\left( \frac{1+z}{1} \right) \right)
$$
需要在n_H函数中根据红移z计算对应的delta_c。

（3）dlnsigma_dlnM的括号错误 #

原代码分母2*np.log(M+eps/M-eps)少了括号，应该是2*np.log((M+eps)/(M-eps))，否则会导致斜率计算完全错误。

2. Sigma(M)计算修正 #

（1）功率谱与积分范围错误 #

• 原Eisenstein-Hu转移函数的实现过于简化，这里改用更准确的近似形式，同时积分的k范围应该覆盖从1e-4/h到1e2/h（远大于原代码的10k），确保积分收敛。
• R的计算中，平均物质密度的数值错误，应该使用rho_crit = 2.775e11 * h**2（单位：Msun/Mpc³），而非直接用1e11。
• Sigma(M)的积分逻辑错误，正确的计算是对x²Pk(x)exp(-x²R²)积分后除以2π²，再开根号，原代码错误地乘以了k³。

修正后的完整代码 #

import numpy as np
import matplotlib.pyplot as plt

# 宇宙学参数（Planck 2015）
h = 0.678
Omega_m = 0.308
A = 0.3222 / h**3
a = 0.707
rho_crit = 2.775e11 * h**2  # 临界密度，单位：Msun/Mpc³

# 定义质量范围
logMmin = 12  # 最小log10(M/Msun)
logMmax = 15  # 最大log10(M/Msun)
nM = 100
logM = np.linspace(logMmin, logMmax, nM)
M = 10**logM / h  # 转换为h^-1 Msun


def sigma(M):
    # 计算线性密度场的均方根涨落sigma(M)
    # 功率谱近似（Eisenstein-Hu简化版）
    def Pk(k):
        k_eq = 0.015 * h  # 等价尺度
        alpha = 1.0 - 0.328 * np.log(431*h**2) * Omega_m + 0.38*np.log(22.3*h**2)*Omega_m**2
        T = np.where(k < k_eq, 
                     np.log(1+2.34*alpha*k/k_eq) / (2.34*alpha*k/k_eq) * (1+3.89*k/k_eq + (16.1*k/k_eq)**2 + (5.46*k/k_eq)**3 + (6.71*k/k_eq)**4)**(-0.25),
                     (k_eq/k)**1.5 * np.log(1+2.34*alpha*k/k_eq)/(2.34*alpha*k/k_eq))
        return 2 * np.pi**2 / k**3 * T**2 * (k/h)**0.96  # 归一化到k=0.1h/Mpc处的sigma=1（近似）

    # 质量对应的共动半径（Mpc/h）
    R = (3 * M / (4 * np.pi * Omega_m * rho_crit)) ** (1/3)
    # 积分范围：k从1e-4/h到1e2/h
    k_arr = np.geomspace(1e-4/h, 1e2/h, 200)
    integrand = k_arr**2 * Pk(k_arr) * np.exp(-(k_arr * R)**2)
    # 计算积分并得到sigma(M)
    s2 = np.trapz(integrand, k_arr) / (2 * np.pi**2)
    return np.sqrt(s2)


def dlnsigma_dlnM(M):
    # 计算dlnσ/dlnM，用中心差分
    eps = 1e-5 * M
    sigma_plus = sigma(M + eps)
    sigma_minus = sigma(M - eps)
    return (np.log(sigma_plus) - np.log(sigma_minus)) / (np.log((M + eps)/(M - eps)))


def n_H(M, z):
    # 随红移变化的临界密度阈值delta_c(z)
    delta_c = 1.686 * (1 + 0.0123 * np.log((1 + z)/1))
    rho_m = Omega_m * rho_crit  # 平均物质密度，单位：Msun/Mpc³
    s = sigma(M)
    ds_dlnM = np.abs(dlnsigma_dlnM(M))  # 取绝对值，确保n为正
    
    # 正确的Sheth-Tormen公式
    term1 = A * np.sqrt(2 * a / np.pi)
    term2 = rho_m / M
    term3 = ds_dlnM
    term4 = 1 + (s**2 / (a * delta_c**2))
    term5 = np.exp(-a * delta_c**2 / (2 * s**2))
    return term1 * term2 * term3 * term4 * term5 / h**4  # 转换为h^4 Mpc^-3 Msun^-1的单位


# 红移数组
z_arr = [0, 0.1, 0.5, 2, 4, 5]

# 绘图
fig, ax = plt.subplots(figsize=(8,6))
for z in z_arr:
    n = np.array([n_H(m, z) for m in M])
    ax.plot(10**logM, n, label=f"z = {z}")
    # 保存数据
    data = np.column_stack((10**logM, n))
    np.savetxt(f'mass_function_z{z:.1f}.txt', data, header='M[Msun] n[h^4 Mpc^-3 Msun^-1]', fmt='%.6e')

ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlabel('$M\ [M_\odot]$')
ax.set_ylabel('$n(M,z)\ [h^4\ \mathrm{Mpc}^{-3}\ M_\odot^{-1}]$')
ax.set_title('Sheth-Tormen 质量函数')
ax.legend()
ax.set_xlim([10**12, 10**15])
ax.set_ylim([1e-8, 1e-2])
plt.show()

修正后的效果 #

• 红移z的影响会正确体现：高红移下，质量函数会向小质量端偏移，且整体数值降低
• M>1e13时不再出现n值为0的情况
• 对数刻度下会显示正常的质量函数曲线

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