# Source code for colour.phenomenons.rayleigh

#!/usr/bin/env python
# -*- coding: utf-8 -*-

"""
Rayleigh Optical Depth - Scattering in the Atmosphere
=====================================================

Implements rayleigh scattering / optical depth in the atmosphere computation:

-   :func:scattering_cross_section
-   :func:rayleigh_optical_depth
-   :func:rayleigh_scattering

--------
Rayleigh Optical Depth - Scattering in the Atmosphere IPython Notebook
<http://nbviewer.ipython.org/github/colour-science/colour-ipython/blob/master/notebooks/phenomenons/rayleigh.ipynb>_  # noqa

References
----------
.. [1]  On Rayleigh Optical Depth Calculations
<http://journals.ametsoc.org/doi/pdf/10.1175/1520-0426(1999)016%3C1854%3AORODC%3E2.0.CO%3B2>_,  # noqa
DOI: https://doi.org/10.1175/1520-0426(1999)016%3C1854:ORODC%3E2.0.CO;2  # noqa
.. [2]  http://en.wikipedia.org/wiki/Rayleigh_scattering
"""

from __future__ import division, unicode_literals

import math

from colour.colorimetry import (
DEFAULT_SPECTRAL_SHAPE,
SpectralPowerDistribution)

__author__ = 'Colour Developers'
__maintainer__ = 'Colour Developers'
__status__ = 'Production'

__all__ = ['STANDARD_AIR_TEMPERATURE',
'STANDARD_CO2_CONCENTRATION',
'AVERAGE_PRESSURE_MEAN_SEA_LEVEL',
'DEFAULT_LATITUDE',
'DEFAULT_ALTITUDE',
'air_refraction_index_penndorf1957',
'air_refraction_index_edlen1966',
'air_refraction_index_peck1972',
'air_refraction_index_bodhaine1999',
'N2_depolarisation',
'O2_depolarisation',
'F_air_penndorf1957',
'F_air_young1981',
'F_air_bates1984',
'F_air_bodhaine1999',
'molecular_density',
'mean_molecular_weights',
'gravity_list1968',
'scattering_cross_section',
'rayleigh_optical_depth',
'rayleigh_scattering']

STANDARD_AIR_TEMPERATURE = 288.15
"""
*Standard air* temperature :math:T[K] in kelvin degrees (:math:15\circ C).

STANDARD_AIR_TEMPERATURE : numeric
"""

STANDARD_CO2_CONCENTRATION = 300
"""
*Standard air* :math:CO_2 concentration in parts per million (ppm).

STANDARD_CO2_CONCENTRATION : numeric
"""

AVERAGE_PRESSURE_MEAN_SEA_LEVEL = 101325
"""
*Standard air* average pressure :math:Hg at mean sea-level in pascal (Pa).

AVERAGE_PRESSURE_MEAN_SEA_LEVEL : numeric
"""

DEFAULT_LATITUDE = 0
"""
Default latitude in degrees (equator).

DEFAULT_LATITUDE : numeric
"""

DEFAULT_ALTITUDE = 0
"""
Default altitude in meters (sea level).

DEFAULT_ALTITUDE : numeric
"""

[docs]def air_refraction_index_penndorf1957(wavelength, *args):
"""
Returns the air refraction index :math:n_s from given wavelength
:math:\lambda in  micrometers (:math:\mu m) using *Penndorf (1957)*
method.

Parameters
----------
wavelength : numeric
Wavelength :math:\lambda in micrometers (:math:\mu m).
\*args : \*
Arguments.

Returns
-------
numeric
Air refraction index :math:n_s.

--------
air_refraction_index_edlen1966, air_refraction_index_peck1972,
air_refraction_index_bodhaine1999

Examples
--------
>>> air_refraction_index_penndorf1957(0.555)  # doctest: +ELLIPSIS
1.0002777...
"""

wl = wavelength
n = 6432.8 + 2949810 / (146 - wl ** (-2)) + 25540 / (41 - wl ** (-2))
n = n / 1.0e8 + 1
return n

[docs]def air_refraction_index_edlen1966(wavelength, *args):
"""
Returns the air refraction index :math:n_s from given wavelength
:math:\lambda in micrometers (:math:\mu m) using *Edlen (1966)* method.

Parameters
----------
wavelength : numeric
Wavelength :math:\lambda in micrometers (:math:\mu m).
\*args : \*
Arguments.

Returns
-------
numeric
Air refraction index :math:n_s.

--------
air_refraction_index_penndorf1957, air_refraction_index_peck1972,
air_refraction_index_bodhaine1999

Examples
--------
>>> air_refraction_index_edlen1966(0.555)  # doctest: +ELLIPSIS
1.0002777...
"""

wl = wavelength
n = 8342.13 + 2406030 / (130 - wl ** (-2)) + 15997 / (38.9 - wl ** (-2))
n = n / 1.0e8 + 1
return n

[docs]def air_refraction_index_peck1972(wavelength, *args):
"""
Returns the air refraction index :math:n_s from given wavelength
:math:\lambda in micrometers (:math:\mu m) using
*Peck and Reeder (1972)* method.

Parameters
----------
wavelength : numeric
Wavelength :math:\lambda in micrometers (:math:\mu m).
\*args : \*
Arguments.

Returns
-------
numeric
Air refraction index :math:n_s.

--------
air_refraction_index_penndorf1957, air_refraction_index_edlen1966,
air_refraction_index_bodhaine1999

Examples
--------
>>> air_refraction_index_peck1972(0.555)  # doctest: +ELLIPSIS
1.0002777...
"""

wl = wavelength
n = (8060.51 + 2480990 / (132.274 - wl ** (-2)) + 17455.7 /
(39.32957 - wl ** (-2)))
n = n / 1.0e8 + 1
return n

[docs]def air_refraction_index_bodhaine1999(
wavelength,
CO2_concentration=STANDARD_CO2_CONCENTRATION):
"""
Returns the air refraction index :math:n_s from given wavelength
:math:\lambda in micrometers (:math:\mu m) using
*Bodhaine, Wood, Dutton and Slusser (1999)* method.

Parameters
----------
wavelength : numeric
Wavelength :math:\lambda in micrometers (:math:\mu m).
CO2_concentration : numeric
:math:CO_2 concentration in parts per million (ppm).

Returns
-------
numeric
Air refraction index :math:n_s.

--------
air_refraction_index_penndorf1957, air_refraction_index_edlen1966,
air_refraction_index_peck1972

Examples
--------
>>> air_refraction_index_bodhaine1999(0.555)  # doctest: +ELLIPSIS
1.0002777...
"""

wl = wavelength
CCO2 = CO2_concentration
n = ((1 + 0.54 * ((CCO2 * 1e-6) - 300e-6)) *
(air_refraction_index_peck1972(wl) - 1) + 1)
return n

[docs]def N2_depolarisation(wavelength):
"""
Returns the depolarisation of nitrogen :math:N_2 as function of
wavelength :math:\lambda in micrometers (:math:\mu m).

Parameters
----------
wavelength : numeric
Wavelength :math:\lambda in micrometers (:math:\mu m).

Returns
-------
numeric
Nitrogen :math:N_2 depolarisation.

Examples
--------
>>> N2_depolarisation(0.555)  # doctest: +ELLIPSIS
1.0350291...
"""

wl = wavelength
N2 = 1.034 + 3.17 * 1.0e-4 * (1 / wl ** 2)

return N2

[docs]def O2_depolarisation(wavelength):
"""
Returns the depolarisation of oxygen :math:O_2 as function of
wavelength :math:\lambda in micrometers (:math:\mu m).

Parameters
----------
wavelength : numeric
Wavelength :math:\lambda in micrometers (:math:\mu m).

Returns
-------
numeric
Oxygen :math:O_2 depolarisation.

Examples
--------
>>> O2_depolarisation(0.555)  # doctest: +ELLIPSIS
1.1020225...
"""

wl = wavelength
O2 = (1.096 + 1.385 * 1.0e-3 * (1 / wl ** 2) +
1.448 * 1.0e-4 * (1 / wl ** 4))

return O2

[docs]def F_air_penndorf1957(*args):
"""
Returns :math:(6+3_p)/(6-7_p), the depolarisation term :math:F(air) or
*King Factor* using *Penndorf (1957)* method.

Parameters
----------
\*args : \*
Arguments.

Returns
-------
numeric
Air depolarisation.

--------
F_air_young1981, F_air_bates1984, F_air_bodhaine1999

Notes
-----
-   The argument *wavelength* is only provided for consistency with the
other air depolarisation methods but is actually not used as this
definition is essentially a constant in its current implementation.

Examples
--------
>>> F_air_penndorf1957(0.555)
1.0608
"""

return 1.0608

[docs]def F_air_young1981(*args):
"""
Returns :math:(6+3_p)/(6-7_p), the depolarisation term :math:F(air) or
*King Factor* using *Young (1981)* method.

Parameters
----------
\*args : \*
Arguments.

Returns
-------
numeric
Air depolarisation.

--------
F_air_penndorf1957, F_air_bates1984, F_air_bodhaine1999

Notes
-----
-   The argument *wavelength* is only provided for consistency with the
other air depolarisation methods but is actually not used as this
definition is essentially a constant in its current implementation.

Examples
--------
>>> F_air_young1981(0.555)
1.048
"""

return 1.0480

[docs]def F_air_bates1984(wavelength):
"""
Returns :math:(6+3_p)/(6-7_p), the depolarisation term :math:F(air) or
*King Factor* as function of wavelength :math:\lambda in micrometers
(:math:\mu m) using *Bates (1984)* method.

Parameters
----------
wavelength : numeric
Wavelength :math:\lambda in micrometers (:math:\mu m).

Returns
-------
numeric
Air depolarisation.

--------
F_air_penndorf1957, F_air_young1981, F_air_bodhaine1999

Examples
--------
>>> F_air_bates1984(0.555)  # doctest: +ELLIPSIS
1.0481535...
"""

wl = wavelength
O2 = O2_depolarisation(wl)
N2 = N2_depolarisation(wl)
Ar = 1.00
CO2 = 1.15

F_air = ((78.084 * N2 + 20.946 * O2 + CO2 + Ar) /
(78.084 + 20.946 + Ar + CO2))

return F_air

[docs]def F_air_bodhaine1999(wavelength,
CO2_concentration=STANDARD_CO2_CONCENTRATION):
"""
Returns :math:(6+3_p)/(6-7_p), the depolarisation term :math:F(air) or
*King Factor* as function of wavelength :math:\lambda in micrometers
(:math:\mu m) and :math:CO_2 concentration in parts per million (ppm)
using *Bodhaine, Wood, Dutton and Slusser (1999)* method.

Parameters
----------
wavelength : numeric
Wavelength :math:\lambda in micrometers (:math:\mu m).
CO2_concentration : numeric, optional
:math:CO_2 concentration in parts per million (ppm).

Returns
-------
numeric
Air depolarisation.

--------
F_air_penndorf1957, F_air_young1981, F_air_bates1984

Examples
--------
>>> F_air_bodhaine1999(0.555)  # doctest: +ELLIPSIS
1.1246916...
"""

wl = wavelength
O2 = O2_depolarisation(wl)
N2 = N2_depolarisation(wl)
CO2_c = CO2_concentration

F_air = ((78.084 * N2 + 20.946 * O2 + 0.934 * 1 + CO2_c * 1.15) /
(78.084 + 20.946 + 0.934 + CO2_c))

return F_air

[docs]def molecular_density(temperature=STANDARD_AIR_TEMPERATURE,
"""
Returns the molecular density :math:N_s (molecules :math:cm^{-3})
as function of air temperature :math:T[K] in kelvin degrees.

Parameters
----------
temperature : numeric, optional
Air temperature :math:T[K] in kelvin degrees.
*Avogadro*'s number (molecules :math:mol^{-1}).

Returns
-------
numeric
Molecular density :math:N_s (molecules :math:cm^{-3}).

Notes
-----
-   The *Avogadro*'s number used in this implementation is the one given by
by the Committee on Data for Science and Technology (CODATA):
:math:6.02214179x10^{23}, which is different from the reference [1]_
value :math:6.0221367x10^{23}.

Examples
--------
>>> molecular_density(288.15)  # doctest: +ELLIPSIS
2.5469021...e+19
>>> molecular_density(288.15, 6.0221367e23)  # doctest: +ELLIPSIS
2.5468999...e+19
"""

T = temperature
N_s = (avogadro_constant / 22.4141) * (273.15 / T) * (1 / 1000)
return N_s

[docs]def mean_molecular_weights(CO2_concentration=STANDARD_CO2_CONCENTRATION):
"""
Returns the mean molecular weights :math:m_a for dry air as function of
:math:CO_2 concentration in parts per million (ppm).

Parameters
----------
CO2_concentration : numeric, optional
:math:CO_2 concentration in parts per million (ppm).

Returns
-------
numeric
Mean molecular weights :math:m_a for dry air.

Examples
--------
>>> mean_molecular_weights()  # doctest: +ELLIPSIS
28.9640166...
"""

CO2_c = CO2_concentration * 1.0e-6

m_a = 15.0556 * CO2_c + 28.9595
return m_a

[docs]def gravity_list1968(latitude=DEFAULT_LATITUDE, altitude=DEFAULT_ALTITUDE):
"""
Returns the gravity :math:g in :math:cm/s_2 (gal) representative of the
mass-weighted column of air molecules above the site of given latitude and
altitude using *List (1968)* method.

Parameters
----------
latitude : numeric, optional
Latitude of the site in degrees.
altitude : numeric, optional
Altitude of the site in meters.

Returns
-------
numeric
Gravity :math:g in :math:cm/s_2 (gal).

Examples
--------
>>> gravity_list1968()  # doctest: +ELLIPSIS
978.0356070...
>>> gravity_list1968(0, 1500)  # doctest: +ELLIPSIS
977.5726106...

Gravity :math:g for Paris:

>>> gravity_list1968(48.8567, 35)  # doctest: +ELLIPSIS
980.9524178...
"""

# Sea level acceleration of gravity.
g0 = 980.6160 * (1 - 0.0026373 * cos2phi + 0.0000059 * cos2phi ** 2)

g = (g0 - (3.085462e-4 + 2.27e-7 * cos2phi) * altitude +
(7.254e-11 + 1.0e-13 * cos2phi) * altitude ** 2 -
(1.517e-17 + 6e-20 * cos2phi) * altitude ** 3)
return g

[docs]def scattering_cross_section(wavelength,
CO2_concentration=STANDARD_CO2_CONCENTRATION,
temperature=STANDARD_AIR_TEMPERATURE,
n_s=air_refraction_index_bodhaine1999,
F_air=F_air_bodhaine1999):
"""
Returns the scattering cross section per molecule :math:\sigma of dry air
as function of wavelength :math:\lambda in centimeters (cm) using given
:math:CO_2 concentration in parts per million (ppm) and temperature
:math:T[K] in kelvin degrees following *Van de Hulst (1957)* method.

Parameters
----------
wavelength : numeric
Wavelength :math:\lambda in centimeters (cm).
CO2_concentration : numeric, optional
:math:CO_2 concentration in parts per million (ppm).
temperature : numeric, optional
Air temperature :math:T[K] in kelvin degrees.
*Avogadro*'s number (molecules :math:mol^{-1}).
n_s : object
Air refraction index :math:n_s computation method.
F_air : object
:math:(6+3_p)/(6-7_p), the depolarisation term :math:F(air) or
*King Factor* computation method.

Returns
-------
numeric
Scattering cross section per molecule :math:\sigma of dry air.

Warning
-------
Unlike most objects of :mod:colour.phenomenons.rayleigh module,
:func:colour.phenomenons.rayleigh.scattering_cross_section expects
wavelength :math:\lambda to be expressed in centimeters (cm).

Examples
--------
>>> scattering_cross_section(555 * 10e-8)  # doctest: +ELLIPSIS
4.6613309...e-27
"""

wl = wavelength
wl_micrometers = wl * 10e3
n_s = n_s(wl_micrometers)
F_air = F_air(wl_micrometers, CO2_concentration)

sigma = (24 * math.pi ** 3 * (n_s ** 2 - 1) ** 2 /
(wl ** 4 * N_s ** 2 * (n_s ** 2 + 2) ** 2))
sigma *= F_air

return sigma

[docs]def rayleigh_optical_depth(wavelength,
CO2_concentration=STANDARD_CO2_CONCENTRATION,
temperature=STANDARD_AIR_TEMPERATURE,
pressure=AVERAGE_PRESSURE_MEAN_SEA_LEVEL,
latitude=DEFAULT_LATITUDE,
altitude=DEFAULT_ALTITUDE,
n_s=air_refraction_index_bodhaine1999,
F_air=F_air_bodhaine1999):
"""
Returns the rayleigh optical depth :math:T_r(\lambda) as function of
wavelength :math:\lambda in centimeters (cm).

Parameters
----------
wavelength : numeric
Wavelength :math:\lambda in centimeters (cm).
CO2_concentration : numeric, optional
:math:CO_2 concentration in parts per million (ppm).
temperature : numeric, optional
Air temperature :math:T[K] in kelvin degrees.
pressure : numeric
Surface pressure :math:P of the measurement site.
latitude : numeric, optional
Latitude of the site in degrees.
altitude : numeric, optional
Altitude of the site in meters.
*Avogadro*'s number (molecules :math:mol^{-1}).
n_s : object
Air refraction index :math:n_s computation method.
F_air : object
:math:(6+3_p)/(6-7_p), the depolarisation term :math:F(air) or
*King Factor* computation method.

Returns
-------
numeric
Rayleigh optical depth :math:T_r(\lambda).

Warning
-------
Unlike most objects of :mod:colour.phenomenons.rayleigh module,
:func:colour.phenomenons.rayleigh.rayleigh_optical_depth expects
wavelength :math:\lambda to be expressed in centimeters (cm).

Examples
--------
>>> rayleigh_optical_depth(555 * 10e-8)  # doctest: +ELLIPSIS
0.1004070...
"""

CO2_c = CO2_concentration
sigma = scattering_cross_section(wavelength,
CO2_c,
temperature,
n_s,
F_air)

# Conversion from pascal to dyne/cm2.
P = pressure * 10
m_a = mean_molecular_weights(CO2_c)
g = gravity_list1968(latitude, altitude)

T_R = sigma * (P * avogadro_constant) / (m_a * g)

return T_R

rayleigh_scattering = rayleigh_optical_depth

def rayleigh_scattering_spd(shape=DEFAULT_SPECTRAL_SHAPE,
CO2_concentration=STANDARD_CO2_CONCENTRATION,
temperature=STANDARD_AIR_TEMPERATURE,
pressure=AVERAGE_PRESSURE_MEAN_SEA_LEVEL,
latitude=DEFAULT_LATITUDE,
altitude=DEFAULT_ALTITUDE,
n_s=air_refraction_index_bodhaine1999,
F_air=F_air_bodhaine1999):
"""
Returns the rayleigh spectral power distribution for given spectral shape.

Parameters
----------
shape : SpectralShape, optional
Spectral shape used to create the rayleigh scattering spectral power
distribution.
CO2_concentration : numeric, optional
:math:CO_2 concentration in parts per million (ppm).
temperature : numeric, optional
Air temperature :math:T[K] in kelvin degrees.
pressure : numeric
Surface pressure :math:P of the measurement site.
latitude : numeric, optional
Latitude of the site in degrees.
altitude : numeric, optional
Altitude of the site in meters.
*Avogadro*'s number (molecules :math:mol^{-1}).
n_s : object
Air refraction index :math:n_s computation method.
F_air : object
:math:(6+3_p)/(6-7_p), the depolarisation term :math:F(air) or
*King Factor* computation method.

Returns
-------
SpectralPowerDistribution
Rayleigh optical depth spectral power distribution.

Examples
--------
>>> rayleigh_scattering_spd()  # doctest: +ELLIPSIS
<colour.colorimetry.spectrum.SpectralPowerDistribution object at 0x...>
"""

return SpectralPowerDistribution(
name=('Rayleigh Scattering - {0} ppm, {1} K, {2} Pa, {3} Degrees, '
'{4} m').format(CO2_concentration,
temperature,
pressure,
latitude,
altitude),
data=dict((wavelength, rayleigh_scattering(wavelength * 10e-8,
CO2_concentration,
temperature,
pressure,
latitude,
altitude,