'''Analytic and numeric tools to compute solid angle of a rectangular pyramid
'''
import numpy as np
from numba import njit, float32
from scipy.integrate import dblquad
[docs]def solidangle_rectangle_num(a, b, xp, yp, zp):
"""Compute the solid angle subtended by a rectangle of length a and b
as seen by an observer located at (xp,yp,zp).
Assume that the center of the rectangle is located at the origin and
that the rectangle lie in the (xOy) plane.
The surface integral is performed numerically.
:param a: rectangle side length
:type a: float>0
:param b: rectangle side length (the other one)
:type b: float>0
:param xp: observer position along x axis
:type xp: float
:param yp: observer position along y axis
:type yp: float
:param zp: observer position along z axis
:type zp: float
:return: solid angle real value
:rtype: float>0
"""
xlow = -a / 2
xhigh = a / 2
ylow = -b / 2
yhigh = b / 2
return zp * dblquad(
lambda x, y: ((xp - x) ** 2 + (yp - y) ** 2 + zp ** 2) ** (-3 / 2.),
xlow, xhigh, ylow, yhigh
)[0]
[docs]@njit(float32(float32, float32, float32, float32, float32))
def solidangle_rectangle_ana_numba(a, b, xp, yp, zp): # pragma: no cover
"""Compute the solid angle subtended by a rectangle of length a and b
as seen by an observer located at (xp,yp,zp).
Assume that the center of the rectangle is located at the origin and
that the rectangle lie in the (xOy) plane.
The surface integral is computed analytically.
documentation:
* note_solid_angle_computation_pb_v2.pdf (version v2)
..warning: Some debugging in coordination with AlG reveal some round off
numerical error performed during the JIT compilation step. Hence the
pure python implementation in solidangle_rectangle_ana.
:param a: rectangle side length
:type a: float>0
:param b: rectangle side length (the other one)
:type b: float>0
:param xp: observer position along x axis
:type xp: float
:param yp: observer position along y axis
:type yp: float
:param zp: observer position along z axis
:type zp: float
:return: solid angle real value
:rtype: float>0
"""
c = zp
# K1 integral
d1 = np.sqrt((0.5 * b - yp) ** 2 + c * c)
num_term_part1 = (0.5 * a - xp) * np.sqrt(d1 * d1 - c * c)
den_term_part1 = c * np.sqrt(d1 * d1 + (0.5 * a - xp) ** 2)
num_term_part2 = -(0.5 * a + xp) * np.sqrt(d1 * d1 - c * c)
den_term_part2 = c * np.sqrt(d1 * d1 + (0.5 * a + xp) ** 2)
term_par1 = np.arctan(num_term_part1 / den_term_part1)
term_par2 = np.arctan(num_term_part2 / den_term_part2)
K1 = (1. / (c * np.sqrt(d1 * d1 - c * c))) * (term_par1 - term_par2)
# K2 integral
d2 = np.sqrt((0.5 * b + yp) ** 2 + c * c)
num_term_part1 = (0.5 * a - xp) * np.sqrt(d2 * d2 - c * c)
den_term_part1 = c * np.sqrt(d2 * d2 + (0.5 * a - xp) ** 2)
num_term_part2 = -(0.5 * a + xp) * np.sqrt(d2 * d2 - c * c)
den_term_part2 = c * np.sqrt(d2 * d2 + (0.5 * a + xp) ** 2)
term_par1 = np.arctan(num_term_part1 / den_term_part1)
term_par2 = np.arctan(num_term_part2 / den_term_part2)
K2 = (1. / (c * np.sqrt(d2 * d2 - c * c))) * (term_par1 - term_par2)
return zp * ((0.5 * b - yp) * K1 + (0.5 * b + yp) * K2)