基于Scipy实现圆盘立体角计算的结果不符问题排查

问题背景 #

参考论文《Solid Angle Of A Disk OffAxis》实现圆盘立体角计算，将论文公式转为Python代码（使用Scipy椭圆积分函数），但计算结果与论文测试用例不一致，推测是数学符号与Scipy函数的映射存在误解。

观测点P的位置由圆盘半径$r_m$、观测点投影到圆盘的径向距离$r_O$、观测点与圆盘的高度距离$L$决定，论文根据$r_O$与$r_m$的关系提供了3种计算公式。

测试用例 #

论文提供基于$r_O/r_m$和$L/r_m$比值的测试数据，已解析为Python DataFrame，代码如下：

import pandas as pd
import numpy as np

data = [
    0,   3.4732594, np.nan     , 0.5,   
    0.2, 3.4184435, 3.41844,   0.5,   
    0.4, 3.2435434, 3.24354,   0.5,   
    0.6, 2.9185178, 2.91852,   0.5,   
    0.8, 2.4122535, 2.41225,   0.5,   
    1.0, 1.7687239, 1.76872,   0.5,   
    1.2, 1.1661307, 1.16614,   0.5,   
    1.4, 0.7428889, 0.742893,  0.5,   
    1.6, 0.4841273, 0.484130,  0.5,   
    1.8, 0.3287007, 0.328702,  0.5,   
    2.0, 0.2324189, 0.232420,  0.5,   
    0,   1.8403024, np.nan   , 1,
    0.2, 1.8070933, 1.80709, 1,
    0.4, 1.7089486, 1.70895, 1,
    0.6, 1.5517370, 1.55174, 1,
    0.8, 1.3488367, 1.34883, 1,
    1.0, 1.1226876, 1.12269, 1,
    1.2, 0.9003572, 0.900369, 1,
    1.4, 0.7039130, 0.703917, 1,
    1.6, 0.5436956, 0.543705, 1,
    1.8, 0.4195415, 0.419543, 1,
    2.0, 0.3257993, 0.325801, 1,
    0,   1.0552591, np.nan, 1.5,
    0.2, 1.0405177, 1.04052,  1.5,
    0.4, 0.9975504, 0.997549, 1.5,
    0.6, 0.9301028, 0.930101, 1.5,
    0.8, 0.8441578, 0.844152, 1.5,
    1.0, 0.7472299, 0.747229, 1.5,
    1.2, 0.6472056, 0.647217, 1.5,
    1.4, 0.5509617, 0.550965, 1.5,
    1.6, 0.4632819, 0.463285, 1.5,
    1.8, 0.3866757, 0.386678, 1.5,
    2.0, 0.3217142, 0.321716, 1.5,
    0,   0.6633335, np.nan    , 2,
    0.2, 0.6566352, 0.656633, 2,
    0.4, 0.6370508, 0.637049, 2,
    0.6, 0.6060694, 0.606068, 2,
    0.8, 0.5659755, 0.565969, 2,
    1.0, 0.5195359, 0.519535, 2,
    1.2, 0.4696858, 0.469697, 2,
    1.4, 0.4191714, 0.419175, 2,
    1.6, 0.3702014, 0.370204, 2,
    1.8, 0.3243908, 0.324392, 2,
    2.0, 0.282707, 0.282709,  2,
]

columns = ['ro_rm', 'Omega', 'Omega_a', 'L_rm']
reshaped_data = [data[i:i + 4] for i in range(0, len(data), 4)]
df = pd.DataFrame(reshaped_data, columns=columns)

实现代码 #

已实现包含Heuman Lambda函数、椭圆积分函数及最终立体角计算函数disk_SA_的代码：

import numpy as np
from scipy.special import ellipkinc, ellipk, ellipe, ellipeinc

def E(k):
    """Complete elliptic integral of the second kind."""
    return ellipe(k)

def E_(ksi, k):
    """Incomplete elliptic integral of the second kind."""
    return ellipeinc(ksi, k)

def K(k):
    """Complete elliptic integral of the first kind."""
    return ellipk(k)

def F(ksi, k):
    """Incomplete elliptic integral of the first kind."""
    return ellipkinc(ksi, k)

def GAMMA0(ksi, k):
    """
    Heuman's Lambda function (Λ₀(ξ, k)).
    - ksi: ξ, calculated from the formula in the image.
    - k: Modulus of the elliptic integrals.
    """
    kp = np.sqrt(1 - k**2)  # Complementary modulus
    # Complete and incomplete elliptic integrals
    K_k = ellipk(k)  # K(k): Complete elliptic integral of the first kind
    E_k = ellipe(k)  # E(k): Complete elliptic integral of the second kind
    F_ksi_kp = ellipkinc(ksi, kp)  # F(ξ, k'): Incomplete elliptic integral of the first kind
    E_ksi_kp = ellipeinc(ksi, kp)  # E(ξ, k'): Incomplete elliptic integral of the second kind
    
    # Formula for Λ₀(ξ, k)
    return (2 / np.pi) * (E_k * F_ksi_kp + K_k * E_ksi_kp - K_k * F_ksi_kp)

def compute_xi(alpha, k):
    """
    Compute ξ (ksi) based on the formula:
    ξ = sin⁻¹(sqrt((α² - k²) / (α² * k'²))).
    - alpha: α, related to the geometry of the system.
    - k: Modulus of the elliptic integrals.
    """
    kp = np.sqrt(1 - k**2)  # Complementary modulus
    numerator = alpha**2 - k**2
    denominator = alpha**2 * kp**2
    sqrt_argument = numerator / denominator
    return np.arcsin(np.sqrt(sqrt_argument))


def disk_SA_(L, r_m, r_o):
    """Solid angle of a disk viewed from an off-axis point.

    - L: Height above the disk (z coordinate).
    - r_m: Radius of the disk.
    - r_o: Radial position of the observation point.
    """
    R_max = np.sqrt(L**2 + (r_o + r_m)**2)
    R1 = np.sqrt(L**2 + (r_o - r_m)**2)
    
    k = np.sqrt(4 * r_o * r_m / (L**2 + (r_o + r_m)**2))
    alpha = np.sqrt(4 * r_o * r_m / ((r_o + r_m)**2))

    ksi = compute_xi(alpha, k)

    if r_o <= r_m:
        if r_o == r_m:
            return np.pi - 2 * L / R_max * K(k)
        elif r_o == 0:
            return 2 * np.pi * (1 - L / R_max)
        else:
            return 2 * np.pi - 2 * L / R_max * K(k) - np.pi * GAMMA0(ksi, k)
    else:
        return -2 * L / R_max * K(k) + np.pi * GAMMA0(ksi, k)

调试情况 #

针对$r_O=r_m$的简单情况，已验证$R_{max}$、参数$k$的计算逻辑正确，且确认使用Scipy的ellipk实现Legendre第一类完全椭圆积分$K(k)$，但计算结果仍与测试用例不符，请求定位错误。

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