r"""
.. autosummary::
henyey_greenstein_1
henyey_greenstein_2
hapke_shoe
h_function_o2
pg_legendre_an
pg_legendre_bn_hg1
pg_legendre_bn_hg2
pg_legendre_pn
hapke_roughness
.. autofunction:: henyey_greenstein_1
.. autofunction:: henyey_greenstein_2
.. autofunction:: henyey_greenstein_3
.. autofunction:: hapke_shoe
.. autofunction:: hapke_shoe_erf
.. autofunction:: hapke_shoe_retro
.. autofunction:: h_function_o1
.. autofunction:: h_function_o2
.. autofunction:: pg_legendre_an
.. autofunction:: pg_legendre_bn_isotropic
.. autofunction:: pg_legendre_bn_hg1
.. autofunction:: pg_legendre_bn_hg2
.. autofunction:: pg_legendre_pn
.. autofunction:: roughness_dxs
.. autofunction:: hapke_roughness
"""
import logging
import numpy as np
from scipy.special import erf
logger = logging.getLogger(__name__)
FLT_EPS = 2.3e-16
FLT_XMAX = 1.0e256
IPI = 1.0e0 / np.pi
[docs]
def henyey_greenstein_1(phase_angle: float | np.ndarray, b: float) -> np.ndarray:
r"""
Henyey-Greenstein single lobe
Single particle phase function as defined in :cite:t:`Hapke_1993`
.. math::
P_{HG,1}(g, b) = \frac{(1 - b^2)}{(1 + 2\cdot b\cdot \cos(g) + b^2)^{3/2}}
:name: eq_hg1
Parameters
----------
array[float] phase_angle:
Phase angle array in radians, noted as g
float b:
Asymmetry parameter (-1 :math:`\leq` b :math:`\leq` 1).
Returns
-------
array[float]
Henyey-Greenstein phase function values.
"""
return np.asarray(
[
(1.0e0 - b * b)
* np.exp(-1.5e0 * np.log(1.0e0 + 2.0e0 * b * np.cos(phase_angle) + b * b))
]
).flatten()
[docs]
def henyey_greenstein_2(phase_angle: np.ndarray, b: float, c: float) -> np.ndarray:
r"""
Henyey-Greenstein double conjugated lobes
Single particle phase function, with conjugated forward
and backward lobes, as defined in :cite:t:`Hapke_1993` and in :cite:t:`McGuire_1995`
.. math::
P_{HG,2}(g, b, c) = \frac{(1 + c)}{2}\cdot P_{HG,1}(g, -b)
+ \frac{(1 - c)}{2}\cdot P_{HG,1}(g, b)
:name: eq_hg2
Parameters
----------
array[float] phase_angle:
Phase angle array in radians, noted as g
float b:
Asymmetry parameter (0 :math:`\leq` b :math:`\leq` 1).
float c:
Backscattering fraction (-1 :math:`\leq` c :math:`\leq` 1).
Returns
-------
array[float]
Two-term Henyey-Greenstein phase function values.
"""
return 0.5 * (1.0e0 + c) * henyey_greenstein_1(phase_angle, -b) + 0.5 * (
1.0e0 - c
) * henyey_greenstein_1(phase_angle, b)
[docs]
def henyey_greenstein_3(
phase_angle: np.ndarray, b1: float, b2: float, c: float
) -> np.ndarray:
r"""
Henyey-Greenstein with two independant lobes.
Single particle phase function, with two independant forward and backward
lobes, as defined in :cite:t:`Hapke_1993` and in :cite:t:`McGuire_1995`.
.. math::
P_{HG,3}(g, b_1, b_2, c) = \frac{(1 + c)}{2}\cdot P_{HG,1}(g, b_1)
+ \frac{(1 - c)}{2}\cdot P_{HG,1}(g, -b_2)
:name: eq_hg3
Parameters
----------
array[float] phase_angle:
Phase angle array in radians, noted as g
float b1:
Asymmetry parameter for backward-scattering lobe (0 :math:`\leq` bm :math:`\leq` 1).
float b2:
Asymmetry parameter for forward-scattering lobe (0 :math:`\leq` bp :math:`\leq` 1).
float c:
Relative lobe strength (-1 :math:`\leq` c :math:`\leq` 1).
Returns
-------
array[float]
Three-term Henyey-Greenstein phase function values.
"""
return 0.5 * (1.0e0 + c) * henyey_greenstein_1(phase_angle, b1) + 0.5 * (
1.0e0 - c
) * henyey_greenstein_1(phase_angle, -b2)
[docs]
def hapke_shoe(phase_angle: np.ndarray, bsh: float, hsh: float) -> np.ndarray:
r"""
Hapke shadow hiding opposition effect (SHOE) function. This function
is an approximation of the exact Hapke SHOE integral, corresponding
to the case of a particulate medium with either a step distribution of
matter or a hyperbolic distribution.
.. math::
:name: eq_hapke_shoe
B_{SH}(g) = \frac{1}{1 + \frac{1}{h_{sh}} \tan(\frac{g}{2})}
Parameters
----------
array[float] phase_angle:
Phase angle array in radians, noted as g
float bsh:
Amplitude of the shadow hiding opposition effect (0 :math:`\leq` bsh :math:`\leq` 1).
float hsh:
Angular width of the shadow hiding opposition effect (hsh > 0).
Returns
-------
array[float]
Hapke shadow hiding opposition effect function values.
.. note::
Using this function implies a *porous* particulate medium, thus any
model using this function should also include the porosity parameter
(K) as defined in :cite:t:`Hapke_1986`.
"""
return np.asarray(
1.0e0 + bsh * np.exp(-1.0e0 * np.log(1.0e0 + np.tan(0.5 * phase_angle) / hsh))
)
[docs]
def hapke_shoe_erf(phase_angle: np.ndarray, bsh: float, hsh: float) -> np.ndarray:
r"""
Hapke shadow hiding opposition effect (SHOE) function using error function.
This function is the exact expression of the SHOE when considering a slab
of particles, which a step distribution of matter.
.. math::
:name: eq_hapke_shoe_erf
B_{SH}(g) = 1 + b_{sh}
\biggl[ \exp{\left(-3 y^2\right)} - 1 +
2 \sqrt{\pi} \cdot y \cdot \exp{(y^2)} \cdot
\bigl[\erf{(2y)} - \erf{(y)}\bigr]
\biggr]
where :math:`y = \sqrt{\frac{h_{sh}}{\tan{\left(\frac{g}{2}\right)}}}`
Parameters
----------
array[float] phase_angle:
Phase angle array in radians, noted as g
float bsh:
Amplitude of the shadow hiding opposition effect (0 :math:`\leq` bsh :math:`\leq` 1).
float hsh:
Angular width of the shadow hiding opposition effect (hsh > 0).
Returns
-------
array[float]
Hapke shadow hiding opposition effect function values.
"""
y = np.sqrt(hsh / np.tan(0.5e0 * phase_angle))
return np.asarray(
1.0e0
+ bsh
* (
np.exp(-3.0e0 * y * y)
- 1.0e0
+ 2.0e0
* np.sqrt(np.pi)
* (erf(2.0e0 * y) - erf(y))
* ((y * np.exp(y * y)) < FLT_XMAX)
)
)
[docs]
def hapke_shoe_retro(phase_angle: np.ndarray, bsh: float, hsh: float) -> np.ndarray:
r"""
Hapke's initial analytic model for the opposition effect. This function
is referred to as the *retrodirective function* in the 1963, 1981 papers.
.. math::
:name: eq_hapke_shoe_retro
B_{SH}(g, h) = 1+B_{0} \cdot
\Biggl[
1-\frac{ \tan{\left(g\right)} }{ 2\cdot h }\cdot
\left(3-e^{ -\frac{h}{ \tan{\left(g\right)} } } \right)\cdot
\left(1-e^{ -\frac{h}{ \tan{\left(g\right)} } } \right)
\Biggr]
where :math:`B_{0}` is the amplitude of the opposition effect, and
:math:`h` is the angular width of the opposition effect.
This expression is valid for phase angles \|g\| <= 90 degrees, beyond
which :math:`B_{SH}(g)` = 1.
Parameters
----------
array[float] phase_angle:
Phase angle array in radians, noted as g
float bsh:
Amplitude of the shadow hiding opposition effect (0 :math:`\leq` bsh :math:`\leq` 1).
float hsh:
Angular width of the shadow hiding opposition effect (hsh > 0).
Returns
-------
array[float]
Evaluation of the retrodirective function.
.. note::
This function is deprecated and should not be used in new models.
"""
# One everywhere...
f_arr = np.ones_like(phase_angle)
# Except where |g| <= 90 degrees
mask = np.abs(phase_angle) < 0.5 * np.pi
f_arr[mask] += bsh * (
1.0e0
- 0.5e0
* np.tan(np.abs(phase_angle[mask]))
/ hsh
* (3.0e0 - np.exp(-hsh / np.tan(np.abs(phase_angle[mask]))))
* (1.0e0 - np.exp(-hsh / np.tan(np.abs(phase_angle[mask]))))
)
return f_arr
# ------ Functions related to the multiple scattering component ------ #
# -------------------------------------------------------------------- #
# This contribution is approached either assuming an isotropic multiple
# scattering through the medium, in which case this contribution is evaluated
# with the approximation of the Ambartsumian-Chandrasekhar H function, and
# leads to the IMSA model, or assuming an anisotropic multiple scattering,
# in which case the single particle phase function is expanded in Legendre
# polynomials, thus leading to the MIMSA model.
[docs]
def h_function_o1(w: float, x: np.ndarray) -> np.ndarray:
r"""
1st order approximation of the Ambartzumian-Chandrasekhar H function.
.. math::
:name: eq_h_function_o1
H(x) = \frac{1+2\cdot x}{1 - 2\cdot x\cdot \sqrt{1 - w}}
Parameters
----------
float w:
Single scattering albedo (0 :math:`\leq` w :math:`\leq` 1).
array[float] x:
Cosine of incidence or emergence angle (mu0 or mu).
Returns
-------
array[float]
1st order approximation of Ambartsumian-Chandrasekhar H function values.
.. note::
This function is an approximation of the exact H function, supposedly
with relative errors less than 4% over the single scattering albedo
range 0 < w < 1, and for all x values between 0 and 1.
The use of the 2nd order approximation should preferred in all cases,
in due of its better accuracy.
"""
return np.asarray((1.0e0 + 2.0e0 * x) / (1.0e0 + 2.0e0 * x * np.sqrt(1.0e0 - w)))
[docs]
def h_function_o2(w: float, x: np.ndarray) -> np.ndarray:
r"""
2nd order approximation of the Ambartzumian-Chandrasekhar H function.
.. math::
:name: eq_h_function_o2
H(x) = \biggl[1 - w \cdot x \cdot
\left(r_{0} +
\frac{1 - 2\cdot x\cdot r_{0}}{2} \cdot
\log\left(1 + \frac{1}{x}\right)
\right)\biggr]^{-1}
Parameters
----------
float w:
Single scattering albedo (0 :math:`\leq` w :math:`\leq` 1).
array[float] x:
Cosine of incidence or emergence angle (mu0 or mu).
Returns
-------
array[float]
2nd order approximation of Ambartsumian-Chandrasekhar H function values.
.. note::
This function is an approximation of the exact H function, supposedly
with relative errors less than 1% up to w ~ 1.
"""
r0 = (1.0e0 - np.sqrt(1.0e0 - w)) / (1.0e0 + np.sqrt(1.0e0 - w))
return np.asarray(
(
1.0e0
- w * x * (r0 + 0.5 * (1.0e0 - 2.0e0 * r0 * x) * np.log((1.0e0 + x) / x))
)
** (-1.0)
)
def imsa_ho1(w: float, mi: np.ndarray, me: np.ndarray) -> np.ndarray:
r"""
Evaluate the multiple scattering contribution using the first order
approximation of the H function.
Parameters
----------
float w:
Single scattering albedo (0 :math:`\leq` w :math:`\leq` 1).
array[float] mi:
cosine of incidence angles (noted µ_{0} or mu0).
array[float] me:
Cosine of emergence angles (noted µ or mue).
Returns
-------
array[float]
1st order approximation of the multiple scattering term
"""
return h_function_o1(w, mi) * h_function_o1(w, me) - 1.0e0
def imsa_ho2(w: float, mi: np.ndarray, me: np.ndarray) -> np.ndarray:
r"""
Evaluate the multiple scattering contribution using the seconc order
approximation of the H function.
Parameters
----------
float w:
Single scattering albedo (0 :math:`\leq` w :math:`\leq` 1).
array[float] mi:
cosine of incidence angles (noted µ_{0} or mu0).
array[float] me:
Cosine of emergence angles (noted µ or mue).
Returns
-------
array[float]
2nd order approximation of the multiple scattering term
"""
return h_function_o2(w, mi) * h_function_o2(w, me) - 1.0e0
# --- MIMSA related function: Legendre polynomial expansions --- #
[docs]
def pg_legendre_an(n: int):
r"""
Compute the A$_{n}$ coefficients for the Legendre polynomial expansion
in the case of an anisotropic multiple scattering.
Using the expression from :cite:t:`Hapke_2002`
.. math::
:name: eq_pg_legendre_an
A_{1} = -0.5 \\
A_{2p} = 0, p in \mathbb{N} \\
A_{2p+1} = A_{2p-1} \cdot \frac{1-2\cdot p}}{2\cdot p+3)}, p in \mathbb{N}^*
Parameters
----------
int n:
Number of coefficients to compute.
Returns
-------
array[float]
An coefficients array.
"""
an = np.zeros(n, dtype=float)
an[1] = -0.5e0
for i in range(3, n, 2):
an[i] = an[i - 2] * (2.0e0 - 1.0e0 * i) / (1.0e0 * i + 1.0e0)
return an
[docs]
def pg_legendre_bn_isotropic(n: int):
r"""
Sets to one the B$_{n}$ coefficients for the Legendre polynomial expansion
in the case of an isotropic phase function.
.. math::
:name: eq_pg_legendre_bn_isotropic
B_{n} = 1, \forall n \in \mathbb{N}
Parameters
----------
int n:
Number of coefficients to compute.
Returns
-------
array[float]
Array set to one.
"""
# TODO: modify this function to become an alias of np.ones
return np.ones(n, dtype=float)
[docs]
def pg_legendre_bn_hg1(n: int, b: float):
r"""
Compute the B$_{n}$ coefficients for the Legendre polynomial expansion
of the single lobe Henyey-Greenstein phase function.
Using the expression from :cite:t:`Hapke_2002`:
.. math::
:name: eq_pg_legendre_bn_hg1
B_{n} = (2\cdot n + 1) \cdot (-b)^{n}
Parameters
----------
int n:
Number of coefficients to compute.
float b:
Asymmetry parameter (-1 :math:`\leq` b :math:`\leq` 1).
Returns
-------
array[float]
Bn coefficients array.
"""
bn = np.zeros(n, dtype=float)
for i in range(0, n):
bn[i] = (2 * i + 1) * ((-b) ** i)
return bn
[docs]
def pg_legendre_bn_hg2(n: int, b: float, c: float):
r"""
Compute the B$_{n}$ coefficients for the Legendre polynomial expansion
of the two-term Henyey-Greenstein phase function.
Using the expression from :cite:t:`Hapke_2002`:
.. math::
:name: eq_pg_legendre_bn_hg2
B_{n} = (2\cdot n + 1) \cdot c \cdot (b)^{n}
Parameters
----------
int n:
Number of coefficients to compute.
float b:
Asymmetry parameter (0 :math:`\leq` b :math:`\leq` 1).
float c:
Backscattering fraction (0 :math:`\leq` c :math:`\leq` 1).
Returns
-------
array[float]
Bn coefficients array.
"""
bn = np.zeros(n, dtype=float)
for i in range(0, n):
bn[i] = c * (2 * i + 1) * (b**i)
return bn
[docs]
def pg_legendre_pn(n: int, x: np.ndarray):
r"""
Compute the first n Legendre polynomials Pn(x) using the Bonnet’s
recursion formula.
.. math::
:name: eq_pg_legendre_pn
P_{0}(x) = 1 \\
P_{1}(x) = x \\
P_{n}(x) = \frac{(2n - 1) \cdot x \cdot P_{n-1}(x) - (n - 1) \cdot P_{n-2}(x)}{n}, n \geq 2
Parameters
----------
int n:
Number of polynomials to compute.
array[float] x:
Input x values array (in the case of the Hapke model x = cos(incidence/emergence angle)).
Returns
-------
array[float]:
2D array containing the first n Legendre polynomials evaluated at each x value.
"""
pn = np.zeros((n, len(x)), dtype=float)
pn[0, :] = 1.0e0
if n >= 1:
pn[1, :] = x
if n >= 2:
for i in range(2, n):
pn[i, :] = (2.0e0 - 1.0e0 / i) * x * pn[i - 1, :] + (
1.0e0 / i - 1.0e0
) * pn[i - 2, :]
return pn
# ------------ Roughness correction related functions ---------------- #
[docs]
def roughness_dxs(
*, inc: np.ndarray, emi: np.ndarray, phi: np.ndarray, theta: float
) -> np.ndarray:
r"""Computes the disk function modulated by the roughness function,
as defined in :cite:t:`Hapke_1984` and :cite:t:`Hapke_1993`.
Parameters
----------
array[float] inc:
the array of incidence angle values (in radians)
array[float] emi:
the array of emergence angle values (in radians)
array[float] azi:
the array of azimuth angle values (in radians)
float theta:
The photometric roughness value (in radians)
Returns
--------
array[float]
the evaluated product of the roughness function
per the roughness-modulated disk function.
Notes
-------
- The function is valid only for azimuth angles between
0 and pi, in due of the 'f' function, approximation chosen
by Hapke (see :cite:t:`Hapke_1993`, chap 12).
- Benchmarking against PH's implementation for 10k runs on an i9-9900K
via jupyter NB and VScode.
old implementation : 0.0422771990 ± 0.0033080535 (s)
this implementation : 0.0191985337 ± 0.0014033128 (s)
"""
msk = inc <= emi
ksm = ~msk
# floats
tt = np.tan(theta)
# arrays
titt = np.tan(inc) * tt
tett = np.tan(emi) * tt
# cos phi and sin(0.5*phi)^2.
cphi = np.cos(phi)
s2p2 = 0.5e0 * (1.0e0 - cphi)
# note to self: tan(0.5 * x) = tan(x) / ( 1 + sqrt(1 + tan(x)^2.) )
# alt f = np.exp(-2.0e0 * np.tan(phi) / (1.e0 + np.sqrt(1.e0+np.tan(phi)*np.tan(phi))))
f = np.exp(-2.0e0 * np.tan(0.5 * phi))
ei = 1.0e0 / titt
ee = 1.0e0 / tett
e1_i = np.exp(-2.0e0 * ei * IPI)
e2_i = np.exp(-1.0e0 * ei * ei * IPI)
e1_e = np.exp(-2.0e0 * ee * IPI)
e2_e = np.exp(-1.0e0 * ee * ee * IPI)
etai = 1.0e0 + titt * e2_i / (2 - e1_i)
etae = 1.0e0 + tett * e2_e / (2 - e1_e)
# Case 1: i <= e
den = 1.0e0 / (2.0e0 - e1_e - phi * e1_i * IPI)
e1ie = (e2_e * cphi + e2_i * s2p2) * den
e2ie = (e2_e - e2_i * s2p2) * den
mieie = 1.0e0 + titt * e1ie
mieee = 1.0e0 + tett * e2ie
s1 = mieee / (etae * (etai * (1.0e0 - f) + f))
d1 = mieie / (mieie + mieee)
s1[ksm] = 0.0e0
d1[ksm] = 0.0e0
# Case 2: i > e
den = 1.0e0 / (2.0e0 - e1_i - phi * e1_e * IPI)
e1ei = (e2_i - e2_e * s2p2) / (2.0e0 - e1_i - phi * e1_e * IPI) * den
e2ei = (e2_i * cphi + e2_e * s2p2) * den
meiie = 1.0e0 + titt * e1ei
meiee = 1.0e0 + tett * e2ei
s2 = meiee / (etai * (etae * (1.0e0 - f) + f))
d2 = meiie / (meiie + meiee)
s2[msk] = 0.0e0
d2[msk] = 0.0e0
return (d1 + d2) * (s1 + s2)
[docs]
def hapke_roughness(
inc: np.ndarray, emi: np.ndarray, azi: np.ndarray, theta: float
) -> np.ndarray:
r"""
Computes the roughness function as defined in :cite:t:`Hapke_1986`. It
returns the ratio of the reflectance masked by mutual shadows.
Parameters
----------
array[float] inc:
the array of incidence angle values (in radians)
array[float] emi:
the array of emergence angle values (in radians)
array[float] azi:
the array of azimuth angle values (in radians)
float theta:
The photometric roughness value (in radians)
Returns
--------
array[float]
the evaluated roughness function
"""
azi_p = azi * np.pi
cos_a = np.cos(azi)
sn2a2 = np.sin(0.5 * azi) ** 2
tan_t = np.tan(theta)
cos_i = np.cos(inc)
tan_i = np.tan(inc)
ti_tt = tan_i * tan_t
si_tt = np.sin(inc) * tan_t
cos_e = np.cos(emi)
te_tt = np.tan(emi) * tan_t
f = np.exp(-2.0e0 * np.tan(0.5e0 * azi))
f2 = 1.0e0 - f
ei = 1.0e0 / ti_tt
e1_i = np.exp(-2.0e0 * ei / np.pi)
e1bi = 2.0e0 - e1_i
e2_i = np.exp(-ei * ei / np.pi)
ee = 1.0e0 / te_tt
e1_e = np.exp(-2.0e0 * ee / np.pi)
e1be = 2.0e0 - e1_e
e2_e = np.exp(-ee * ee / np.pi)
# defining different eta than hapke -> eta(x) = eta_hpk(x)/(\chi{\theta}*cos(x))
eta_i = 1.0e0 + ti_tt * e2_i / e1bi
eta_e = 1.0e0 + te_tt * e2_e / e1be
# Case e greater than i
# -> Hapke 1993 equations: 12.46 -> 12.51
# -> Sofie Spujth equations: 2.36->2.38 + 2.42->2.45
d = e1be - azi_p * e1_i
x = sn2a2 * e2_i
a = cos_e / (cos_i * d + si_tt * (cos_a * e2_e + x))
b = 1.0e0 / (d + te_tt * (e2_e - x))
c = d * eta_e * (f + f2 * eta_i)
dxs_e_gt_i = 1.0e0 / ((a + b) * c)
# Case i greater than e
# -> Hapke 1993 equations:
# -> Sofie Spujth equations:
# -> Shkuratov equations: 12
d = e1bi - azi_p * e1_e
x = sn2a2 * e2_e
a = cos_e / (cos_i * d + si_tt * (e2_i - x))
b = 1.0e0 / (d + te_tt * (cos_a * e2_i + x))
c = d * eta_i * (f + f2 * eta_e)
dxs_i_gt_e = 1.0e0 / ((a + b) * c)
dxs = np.zeros_like(inc)
dxs[emi > inc] = dxs_e_gt_i[emi > inc]
dxs[inc > emi] = dxs_i_gt_e[inc > emi]
return dxs
# -------------------------------------------------------------------- #