r""" Emission functions.
An emission function, :math:`p_\Omega(\theta, \phi)`, are
determined as a probability function
in terms of the solid angle (steradian, :math:`\Omega`).
.. math::
\int p_\Omega(\theta, \phi) d\Omega = 1
It is *not* interms of :math`\theta, \phi`.
In this module, the surface-emission is considered, so that the angle
:math:`\theta` (polar angle) shall be 0 to 90 degrees. Beyond 90 degrees, the emission should be zero.
The angle :math:`\theta` is defined from the local zenith (normal vector of surface). The angle :math:`\phi` is measured from the opposite direction of the reference vector. (Forward scattering shall be defined as zero phi). Increasing phi corresponds to the anti-clockwize looking from zenith.
- For the uniform (isotropic) case, the emission function is :math:`1/2\pi`.
- For the cosine-normal case, the emission function is :math:`(\gamma+1)/2\pi cos^\gamma \theta`.
- For the forward-specular model, the cosine-normal function is tilted by the given angle.
- For the backward-specular model, the cosine-nomal function is tilted (backward) by the given angle.
"""
import numpy as np
# bfb34c1
[docs]def emission_cosnorm(theta, phi, gamma=1):
""" Cosine-normal emission function.
:param theta: The polar angle in deg. [0, 180]. (>90, the returned value is nan).
:param phi: The azimuth angle in deg. (-180, 180]
>>> theta = [30, 50, 120]
>>> phi = [50, 20, -30]
>>> print(emission_uniform(theta, phi)) # doctest: +NORMALIZE_WHITESPACE
[0.15915494 0.15915494 nan]
"""
u = np.power(np.cos(np.radians(theta)), gamma)
emfunc = u / 2 / np.pi * (1 + gamma)
emfunc = np.where(np.array(theta) <= 90, emfunc, np.nan)
return emfunc
[docs]class POMSP:
_CACHED_POMSP = {}
[docs] @classmethod
def getObject(cls, tilt_angle, gamma, normalize_accuracy=360):
key = (tilt_angle, gamma, normalize_accuracy)
if not key in cls._CACHED_POMSP:
cls._CACHED_POMSP[key] = _POMSP(tilt_angle, gamma, normalize_accuracy)
return cls._CACHED_POMSP[key]
[docs] @classmethod
def clear_cache(cls):
cls._CACHED_POMSP = {}
[docs]def clear_cache():
POMSP.clear_cache()
class _POMSP: # I forgot what the name stand for. Specular class
def __init__(self, incom_deg, gamma, normalize_accuracy=360):
self.theta0 = np.radians(incom_deg)
self.gamma = gamma
theta = np.linspace(0, 90, normalize_accuracy)
phi = np.linspace(-180, 180, normalize_accuracy * 4)
dtheta = theta[1] - theta[0]
dphi = phi[1] - phi[0]
theta, phi = np.meshgrid(theta, phi, indexing='ij')
p = self.getP(theta, phi, normalize=False)
steradian = np.sin(np.radians(theta)) * np.radians(dtheta) * np.radians(dphi)
self.normalize_factor = (p * steradian).sum()
def getP(self, theta, phi, normalize=True):
theta_0 = self.theta0
t = np.radians(theta)
p = np.radians(phi)
u = np.sin(theta_0) * np.sin(t) * np.cos(p) + np.cos(theta_0) * np.cos(t)
u = np.where(u > 0, u, 0)
u = np.power(u, self.gamma)
if normalize:
u = u / self.normalize_factor
return u
[docs]def emission_specular(theta, phi, tilt_angle=0, gamma=1, normalize_accuracy=360):
""" Specular reflection
:param theta: The polar angle in deg. [0, 180]. (>90, the returned value is nan).
:param phi: The azimuth angle in deg. (-180, 180]
:keyword tilt_angle: The tilt angle for specular reflection for the forward scattering. [0, 90]
Usually, the angle is the same as the incoming angle of the projectile.
If you set it as 0 (default), the result is the same as cosing_normal.
:keyword gamma: The gamma paramter for the cosine normal function.
:keyword normalize_accuracy: To normalize the results, this function integrate numerically over the half-sphere.
This parameter defines how accurate the normalization should be.
>>> print('{:.3g}'.format(emission_specular(30, 0, tilt_angle=0, gamma=1)))
0.275
>>> print('{:.3g}'.format(emission_specular(30, 0, tilt_angle=30, gamma=1)))
0.341
>>> print('{:.3g}'.format(emission_specular(30, 20, tilt_angle=30, gamma=1)))
0.336
>>> print('{:.3g}'.format(emission_specular(30, -20, tilt_angle=30, gamma=1)))
0.336
>>> print('{:.3g}'.format(emission_specular(30, -30, tilt_angle=30, gamma=1)))
0.329
>>> print('{:.3g}'.format(emission_specular(80, 20, tilt_angle=70, gamma=3)))
0.682
>>> print('{:.3g}'.format(emission_specular(80, -20, tilt_angle=70, gamma=3)))
0.682
>>> print('{:.3g}'.format(emission_specular(100, -20, tilt_angle=70, gamma=3)))
nan
"""
pom = POMSP.getObject(tilt_angle, gamma, normalize_accuracy=normalize_accuracy)
emfunc = pom.getP(theta, phi)
emfunc = np.where(np.array(theta) <= 90, emfunc, np.nan)
return emfunc
[docs]def emission_backward(theta, phi, tilt_angle=0, gamma=1, normalize_accuracy=360):
""" Emission function of backscattering
:param theta: The polar angle in deg. [0, 180]. (>90, the returned value is nan).
:param phi: The azimuth angle in deg. (-180, 180]
:keyword tilt_angle: The tilt angle for specular reflection for the backward scattering. [0, 90]
Usually, the angle is the same as the incoming angle of the projectile.
If you set it as 0 (default), the result is the same as cosing_normal.
:keyword gamma: The gamma paramter for the cosine normal function.
:keyword normalize_accuracy: To normalize the results, this function integrate numerically over the half-sphere.
This parameter defines how accurate the normalization should be.
>>> print('{:.3g}'.format(emission_backward(30, 180, tilt_angle=0, gamma=1)))
0.275
>>> print('{:.3g}'.format(emission_backward(30, 180, tilt_angle=30, gamma=1)))
0.341
>>> print('{:.3g}'.format(emission_backward(30, 160, tilt_angle=30, gamma=1)))
0.336
>>> print('{:.3g}'.format(emission_backward(30, -160, tilt_angle=30, gamma=1)))
0.336
>>> print('{:.3g}'.format(emission_backward(30, -150, tilt_angle=30, gamma=1)))
0.329
>>> print('{:.3g}'.format(emission_backward(80, 160, tilt_angle=70, gamma=3)))
0.682
>>> print('{:.3g}'.format(emission_backward(80, -160, tilt_angle=70, gamma=3)))
0.682
>>> print('{:.3g}'.format(emission_backward(100, -160, tilt_angle=70, gamma=3)))
nan
"""
phibs = np.where(phi < 0, -phi - 180, -phi + 180)
return emission_specular(theta, phibs,
tilt_angle=tilt_angle, gamma=gamma,
normalize_accuracy=normalize_accuracy)