我正在尝试根据以下论文实施圆盘的实体角度的计算：

问题描述 投票：0回答：0
但是我没有得到相同的数值结果，我希望这个问题来自数学符号和Scipy功能之间我的零件的错误解释。

设置

设置如下：观察点P位于磁盘上是否位于磁盘上，具体取决于磁盘半径$ r_m $，以及观察点P的径向距离（投影到磁盘上到O ）R_O。该问题的另一个重要参数是磁盘和观察点之间的高度距离，指出l纸张提供了3个公式此后显示的公式，具体取决于$ r_o $和$ r_m $：

因此，目标是用Pyhon写这三个公式。其他参数（例如RMAX）是基于3个输入RO，RM和L的中间变量，该变量表征了设置。 enter image description here测试案例 一件好事是，可以使用其他计算相同数量的方法，因此很容易生成测试案例 - 此外，本文提供了与其他方法相一致的表格结果。换句话说，我们有一组测试箱，它们使用比率ro/rm和l/rm进行参数化：

enter image description here the是分析数据的Python代码：

import pandas as pd 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)

问题：使用scipy

实施

论文中提供了总体数学推导，我有兴趣仅复制Python函数中的最终方程式 - 不幸的是，我在所有3个情况下都会获得不同的结果（RO

RM），我怀疑我的代码中的某个地方有一个错误。

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) import numpy as np from scipy.special import ellipkinc, ellipk, ellipe, ellipeinc 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) enter image description here最终功能disk_SA_

（SA代表实体角）

depebugging

casince ro = rm具有最简单的方程式，仅依赖于rmax和k（k），我首先检查这些数量。 rmax
此数量是根据基本参数R0，RM和L（在等式（12）之后）给出的：

所以我的代码似乎还可以。
LLE的移动

R_max = np.sqrt(L**2 + (r_o + r_m)**2)



功能k（k）


在等式（12）之后还给出了中间变量k，AS：



enter image description here所以我的代码K(k)
似乎还可以。

现在只有K实现要检查。论文指出：


SoK是Legendre的第一类椭圆形成部分，如果我没记错的话，它将在Scipy中在以下网址实现：
Https：//docs.scipy.org/doc/doc/scipy/scipy/reference/reference/generate/generated/generated/scipy.special.special。 ellipk.html


so使用enter image description here
k = np.sqrt(4 * r_o * r_m / (L**2 + (r_o + r_m)**2))

我希望使用
def K(k):
    """Complete elliptic integral of the first kind."""
    return ellipk(k)

正确的最终结果是正确的。这不是根据测试案例之间的比较：


np.pi - 2 * L / R_max * K(k)

lo/rm = 1的索引= 4的示例（所有结果都是错误的，但是专注于ro = rm允许更容易固定错误，因为方程更简单）：


enter image description here


没有人知道问题的来源？
    
	

文档
r_m = 1
r_os = np.arange(0, 2.2, 0.2)
Ls = [0.5, 1, 1.5, 2]
data = []
for r_o in r_os:
    for L in Ls:
        sa = disk_SA_(L, r_m, r_o)
        sa_integ = disk_SA(0*m, r_o*m, L*m, r_m*m, method='dblquad').value
        sa_integ_trapz = disk_SA(0*m, r_o*m, L*m, r_m*m, method='trapz').value
        data.append([r_o, r_m, r_o/r_m, L, L/r_m, sa, sa_integ, sa_integ_trapz])

df_results = pd.DataFrame(
    data, columns=["r_o", "r_m", "ro_rm", "L", "L_rm", "Omega_implem"]
)
df = pd.merge(df, df_results, on=["ro_rm", 'L_rm'])
df['implem_ok'] = np.isclose(df['Omega'], df['Omega_implem'])
指出存在不同的参数化。看起来Scipy's使用$ m = k^2 $作为参数，因此您要将$ k^2 $传递到其中。
    enter image description here
	




algorithm

math

scipy




    


  

  

    
最新问题


  





  

    © www.soinside.com 2019 - 2025. All rights reserved.