pytrunc.phase module#

Scattering phase functions and Legendre moments.

This module provides analytic scattering phase functions (Henyey-Greenstein, two-term Henyey-Greenstein and Fournier-Forand) and the calculation of the phase function Legendre moments, either numerically for an arbitrary phase function or analytically for the Henyey-Greenstein family.

pytrunc.phase.calc_hg_moments(g: float, m_max: int) → ndarray[tuple[Any, ...], dtype[float64]][source]#

Compute exact Henyey-Greenstein phase moments

  • see Eq. 8 in Kattawar (1975)

Parameters:
gfloat

The Henyey-Greenstein parameter g (measures the asymmetry of the phase matrix)

m_maxint

The maximum moment number to compute, i.e., compute m[0], …, m[m_max]

Returns:
ndarray

The phase moments, a 1-D ndarray of size m_max + 1

References

Kattawar, G. W. (1975). A three-parameter analytic phase function for multiple scattering calculations. Journal of Quantitative Spectroscopy and Radiative Transfer, 15(9), 839-849.

Examples

>>> from pytrunc.phase import calc_hg_moments
>>> calc_hg_moments(g=0.85, m_max=3)
array([1.      , 0.85    , 0.7225  , 0.614125])
pytrunc.phase.calc_moments(phase: ndarray[tuple[Any, ...], dtype[float64]], theta: ndarray[tuple[Any, ...], dtype[float64]], m_max: int, method: str = 'lobatto', theta_unit: str = 'deg', normalize: bool = False, xk: ndarray[tuple[Any, ...], dtype[float64]] | None = None, wk: ndarray[tuple[Any, ...], dtype[float64]] | None = None, pl_costh: ndarray[tuple[Any, ...], dtype[float64]] | None = None) → ndarray[tuple[Any, ...], dtype[float64]][source]#

Calculate the phase matrix moments until m_max moment

Parameters:
phasendarray

The phase matrix, it must be 1-D

thetandarray

The phase matrix angles, it must be 1-D. See the theta_unit parameter for the unit

m_maxint

The maximum moment number to compute, i.e., compute m[0], …, m[m_max]

methodstr, optional

The method used to calculate the moments. Default is ‘lobatto’ (very efficient with a “gauss kind” theta distribution), other choices are ‘simpson’ and ‘trapezoid’ (efficient with a regular theta distribution, use the scipy simpson and trapezoid functions)

theta_unitstr, optional

The unit of the theta angles. Default is ‘deg’, other choice is ‘rad’

normalizebool, optional

If True, normalize such that the first moment is exactly equal to 1. Default is False

xkndarray or None, optional

Force the Lobatto quadrature abscissas, it must be 1-D. Considered only if wk is also provided

wkndarray or None, optional

Force the Lobatto weights, it must be 1-D. Considered only if xk is also provided

pl_costhndarray or None, optional

Force the Legendre polynomials values of cos(theta). The 2-D ndarray shape must be (m_max+1, len(theta))

Returns:
ndarray

The computed phase moments, a 1-D ndarray of size m_max + 1

Notes

  • See Eq. A2 in Wiscombe (1977) for the moment computation using the Lobatto quadrature in [0,pi]

  • For the Lobatto quadrature abscissas and weights calculation see Michels (1963)

References

Michels, H. (1963). Abscissas and weight coefficients for Lobatto quadrature. Mathematics of Computation, 17(83), 237-244.

Wiscombe, W. J. (1977). The delta-M method: Rapid yet accurate radiative flux calculations for strongly asymmetric phase functions. Journal of Atmospheric Sciences, 34(9), 1408-1422.

Examples

>>> import numpy as np
>>> from pytrunc.phase import calc_moments, henyey_greenstein
>>> theta = np.linspace(0.0, 180.0, 1801)
>>> phase = henyey_greenstein(theta, g=0.85, normalize=2)
>>> calc_moments(phase, theta, m_max=3)
array([0.99998936, 0.84998937, 0.72248935, 0.61411436])
pytrunc.phase.calc_tthg_moments(g1: float, g2: float, f: float, m_max: int) → ndarray[tuple[Any, ...], dtype[float64]][source]#

Compute exact Two-term Henyey-Greenstein phase moments

  • see Eq. 11 in Kattawar (1975)

Parameters:
g1float

The first H-G term parameter g (forward part)

g2float

The second H-G term parameter g (backward part)

ffloat

The fraction parameter between the two H-G terms

m_maxint

The maximum moment number to compute, i.e., compute m[0], …, m[m_max]

Returns:
ndarray

The phase moments, a 1-D ndarray of size m_max + 1

References

Kattawar, G. W. (1975). A three-parameter analytic phase function for multiple scattering calculations. Journal of Quantitative Spectroscopy and Radiative Transfer, 15(9), 839-849.

Examples

>>> from pytrunc.phase import calc_tthg_moments
>>> calc_tthg_moments(g1=0.85, g2=-0.6, f=0.9, m_max=3)
array([1.       , 0.705    , 0.68625  , 0.5311125])
pytrunc.phase.fournier_forand(theta: ArrayLike, n: float, mu: float, theta_unit: str = 'deg', normalize: float | None = None) → ndarray[tuple[Any, ...], dtype[float64]][source]#

Compute the Fournier-Forand phase matrix

Parameters:
thetandarray

The phase matrix angles, it must be 1-D. See the theta_unit parameter for the unit

nfloat

The real index of refraction of the particles

mufloat

The slope parameter of the hyperbolic (Junge) particle size distribution (typically between 3 and 5 for oceanic particles)

theta_unitstr, optional

The unit of the theta angles. Default is ‘deg’, other choice is ‘rad’

normalizefloat or None, optional

The normalization value of the integral of F_FF(θ)dcosθ, where F_FF(θ) is the phase matrix. The scipy simpson function is used for the normalization. Default is None, meaning no normalization

Returns:
ndarray

The phase matrix, a 1-D ndarray with the same shape as theta

Notes

The Fournier-Forand equation:

  • \(F_FF(θ) = [1 / (4*pi*(1-δ)**2*δ**v)] * [v*(1-δ) - (1-δ**v) + (δ*(1-δ**v) - v*(1-δ))*sin(θ/2)**(-2)] + [(1-δ_180**v) / (16*pi*(δ_180-1)*δ_180**v)] * (3*cos(θ)**2 - 1)\)

  • with \(v = (3-μ)/2\) and \(δ = (4 / (3*(n-1)**2))*sin(θ/2)**2\), where δ_180 is δ evaluated at θ = 180°

  • By default the integral of F_FF(θ)dcosθ is equal to 1/(2*pi). The integral value can be different due to a very low discretization of θ and/or a strong forward peak. The use of the normalize parameter can be useful to renormalize the phase function

  • F_FF(θ) diverges as θ → 0: at θ = 0 the result is NaN/inf (with numpy invalid-value warnings). Use a grid starting above 0, or overwrite the forward peak values before integrating

References

Fournier, G. R., & Forand, J. L. (1994). Analytic phase function for ocean water. In Ocean Optics XII (Vol. 2258, pp. 194-201). SPIE.

other reference

Examples

>>> import numpy as np
>>> from pytrunc.phase import fournier_forand
>>> theta = np.linspace(10.0, 170.0, 3)
>>> fournier_forand(theta, n=1.117, mu=3.695)
array([1.24049454, 0.0065938 , 0.00479607])
pytrunc.phase.henyey_greenstein(theta: ArrayLike, g: float, theta_unit: str = 'deg', normalize: float | None = None) → ndarray[tuple[Any, ...], dtype[float64]][source]#

Compute the Henyey-Greenstein phase matrix

Parameters:
thetandarray

The phase matrix angles, it must be 1-D. See the theta_unit parameter for the unit

gfloat

The Henyey-Greenstein parameter g (measures the asymmetry of the phase matrix)

theta_unitstr, optional

The unit of the theta angles. Default is ‘deg’, other choice is ‘rad’

normalizefloat or None, optional

The normalization value of the integral of F_HG(θ)dcosθ, where F_HG(θ) is the phase matrix. The scipy simpson function is used for the normalization. Default is None, meaning no normalization

Returns:
ndarray

The phase matrix, a 1-D ndarray with the same shape as theta

Notes

The Henyey-Greenstein equation:

  • \(F_HG(θ) = (1/(4*pi))*[(1-g**2) / (1+g**2-(2*g*cos(θ)))**(3/2)]\)

  • By default the integral of F_HG(θ)dcosθ is equal to 1/(2*pi). The integral value can be different due to a very low discretization of θ and/or a high g value. The use of the normalize parameter can be useful to renormalize the phase function

References

Henyey, L. G., & Greenstein, J. L. (1941). Diffuse radiation in the galaxy. Astrophysical Journal, 93, 70-83.

other reference

Examples

>>> import numpy as np
>>> from pytrunc.phase import henyey_greenstein
>>> theta = np.linspace(0.0, 180.0, 3)
>>> henyey_greenstein(theta, g=0.85)
array([6.54303655e+00, 9.76819403e-03, 3.48769050e-03])
pytrunc.phase.two_term_henyey_greenstein(theta: ArrayLike, g1: float, g2: float, f: float, theta_unit: str = 'deg', normalize: float | None = None) → ndarray[tuple[Any, ...], dtype[float64]][source]#

Compute the two-term Henyey-Greenstein phase matrix

Parameters:
thetandarray

The phase matrix angles, it must be 1-D. See the theta_unit parameter for the unit

g1float

The first H-G term parameter g (forward part)

g2float

The second H-G term parameter g (backward part)

ffloat

The fraction parameter between the two H-G terms (see notes)

theta_unitstr, optional

The unit of the theta angles. Default is ‘deg’, other choice is ‘rad’

normalizefloat or None, optional

The normalization value of the integral of F_TTHG(θ)dcosθ, where F_TTHG(θ) is the phase matrix. The scipy simpson function is used for the normalization. Default is None, meaning no normalization

Returns:
ndarray

The phase matrix, a 1-D ndarray with the same shape as theta

Notes

The two term Henyey-Greenstein equation:

  • \(F_TTHG(θ) = f*F_HG1(θ) + (1-f)*F_HG2(θ)\)

  • By default the integral of F_TTHG(θ)dcosθ is equal to 1/(2*pi). The integral value can be different due to a very low discretization of θ and/or a high g value. The use of the normalize parameter can be useful to renormalize the phase function

References

Irvine, W. M. (1965). Multiple scattering by large particles (No. NASA-CR-64638).

Examples

>>> import numpy as np
>>> from pytrunc.phase import two_term_henyey_greenstein
>>> theta = np.linspace(0.0, 180.0, 3)
>>> two_term_henyey_greenstein(theta, g1=0.85, g2=-0.6, f=0.9)
array([5.88997629, 0.01200253, 0.08271639])