Solid angle of a disk computation with scipy
I am trying to implement the computation of the solid angle of a circular disk based on the following paper:
https://public.websites.umich.edu/~ners311/CourseLibrary/SolidAngleOfADiskOffAxis.pdf
but I am not getting the same numerical results, I expect the problem to come from mis-interpretation of my part between math notations and scipy's functions.
The setup
The setup is the following : the point of observation P is located either over the disk or not over the disk, depending on the disk radius $r_m$, and the radial distance of the point of observation P (projected onto the disk to O) r_O. Another important parameter of the problem is the height distance between the disk and the point of observation, noted L The paper provides 3 formulas that are shown here after, depending on $r_O$ and $r_m$:
So the goal is to write in pyhon those 3 formulas. Other paramters like Rmax are intermediate variables based on the 3 input rO, rm and L that characterize the setup.
Test-cases
A good thing is that other methods to compute the same quantity are available so it is easy to generate test-cases - in addition, the paper provides tabulated results which are coherent from those other methods. In other words, we have a set of test-cases, that are parametrized using the ratio rO/rm and L/rm:
Here is the python code from the parsed data:
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)
The problem: implementation using scipy
The overall math derivation is provided in the paper and I am interested in reproducing just the final equations in python functions - unfortunately I get different results for all 3 cases (rO<rm, rO=rm, rO>rm), and I suspect a mistake somewhere in my code.
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)
The final function disk_SA_ (SA stands for Solid Angle)
Debugging
Since the case rO=rm has the simplest equation and only relies on Rmax and K(k), I first check those quantities.
Rmax This quantity is given in terms of the base parameters r0, rm, and L (right after eq (12)):
so my code R_max = np.sqrt(L**2 + (r_o + r_m)**2) seems alright.
Let's move onto K(k)
Function K(k)
The intermediate variable k is also given after eq (12), as:
so my code k = np.sqrt(4 * r_o * r_m / (L**2 + (r_o + r_m)**2)) also seems ok.
Now there's just K implementation to check. The paper states :
so K is Legendre's complete elliptic integral of the first kind, which if I'm not mistaken, is implemented in scipy at : https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.ellipk.html
so using
def K(k):
"""Complete elliptic integral of the first kind."""
return ellipk(k)
I'd expect the final result using np.pi - 2 * L / R_max * K(k) to be correct. Which is not according to the comparison between the test-cases:
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'])
for example at index=4 for rO/rm=1 (all results are wrong, but focusing on rO=rm allows to pinpout error easier since the equation is simpler):
Does anyone has any idea where the problem comes from ?
The documentation of ellipk mentions that there are different parametrizations.
Note the correspondence between $k^2$ in equation 13 of the paper and $m$ in SciPy's version of the integral.
It looks like SciPy's uses "modulus" m = k**2$ as the argument, so you would want to pass k**2 into it. Passing k**2 instead of k into the elliptic functions everywhere, I'm seeing agreement with the tabulated results.