# Source code for fluids.atmosphere

    """Chemical Engineering Design Library (ChEDL). Utilities for process modeling.
    Copyright (C) 2016, 2017, 2018 Caleb Bell <Caleb.Andrew.Bell@gmail.com>

    Permission is hereby granted, free of charge, to any person obtaining a copy
    of this software and associated documentation files (the "Software"), to deal
    in the Software without restriction, including without limitation the rights
    to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
    copies of the Software, and to permit persons to whom the Software is
    furnished to do so, subject to the following conditions:

    The above copyright notice and this permission notice shall be included in all
    copies or substantial portions of the Software.

    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
    IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
    FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
    AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
    LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
    OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
    SOFTWARE.

    This module contains models of earth's atmosphere. Models are empirical and
    based on extensive research, primarily by NASA.

    For reporting bugs, adding feature requests, or submitting pull requests,
    please use the `GitHub issue tracker <https://github.com/CalebBell/fluids/>`_
    or contact the author at Caleb.Andrew.Bell@gmail.com.

    .. contents:: :local:

    Atmospheres
    -----------
    .. autoclass:: ATMOSPHERE_1976
        :members:
    .. autoclass:: ATMOSPHERE_NRLMSISE00
        :members:
    .. autofunction:: airmass

    Solar Radiation and Position
    ----------------------------
    .. autofunction:: solar_position
    .. autofunction:: solar_irradiation
    .. autofunction:: sunrise_sunset
    .. autofunction:: earthsun_distance

    """

    import os
    from math import cos, exp, pi, radians, sin, sqrt

    from fluids.constants import N_A, R, au
    from fluids.numerics import numpy as np
    from fluids.numerics import quad, secant

    try:
        from datetime import datetime
    except:
        pass

    __all__ = ['ATMOSPHERE_1976', 'ATMOSPHERE_NRLMSISE00',
               'earthsun_distance', 'solar_position', 'solar_irradiation',
               'sunrise_sunset']

    H_std = [0.0, 11E3, 20E3, 32E3, 47E3, 51E3, 71E3, 84852.0]
    T_grad = [-6.5E-3, 0.0, 1E-3, 2.8E-3, 0.0, -2.8E-3, -2E-3, 0.0]
    T_std = [288.15, 216.65, 216.65, 228.65, 270.65, 270.65, 214.65, 186.946]
    P_std = [101325, 22632.06397346291, 5474.8886696777745, 868.0186847552279,
            110.90630555496608, 66.93887311868738, 3.956420428040732,
            0.3733835899762159]

    r0 = 6356766.0
    P0 = 101325.0
    M0 = 28.9644
    g0 = 9.80665
    gamma = 1.400

    def H_for_P_ATMOSPHERE_1976_err(H, P1):
        return ATMOSPHERE_1976(H, 0.0).P - P1

    def to_int_dP_ATMOSPHERE_1976(Z, dT):
        atm = ATMOSPHERE_1976(Z, dT)
        return atm.g*atm.rho

    class ATMOSPHERE_1976:
        r'''US Standard Atmosphere 1976 class, which calculates `T`, `P`,
        `rho`, `v_sonic`, `mu`, `k`, and `g` as a function of altitude above
        sea level. Designed to provide reasonable results up to an elevation
        of 86,000 m (0.4 Pa). The model is also valid under sea level, to
        -610 meters.

        Parameters
        ----------
        Z : float
            Elevation, [m]
        dT : float, optional
            Temperature difference from standard conditions used in determining
            the properties of the atmosphere, [K]

        Attributes
        ----------
        T : float
            Temperature of atmosphere at specified conditions, [K]
        P : float
            Pressure of atmosphere at specified conditions, [Pa]
        rho : float
            Mass density of atmosphere at specified conditions [kg/m^3]
        H : float
            Geopotential height, [m]
        g : float
            Acceleration due to gravity, [m/s^2]
        mu : float
            Viscosity of atmosphere at specified conditions, [Pa*s]
        k : float
            Thermal conductivity of atmosphere at specified conditions, [W/m/K]
        v_sonic : float
            Speed of sound of atmosphere at specified conditions, [m/s]

        Examples
        --------
        >>> five_km = ATMOSPHERE_1976(5000)
        >>> five_km.P, five_km.rho, five_km.mu
        (54048.28614576141, 0.7364284207799743, 1.628248135362207e-05)
        >>> five_km.k, five_km.g, five_km.v_sonic
        (0.02273190295142526, 9.791241076982665, 320.5455196704035)

        Notes
        -----
        Up to 32 km, the International Standard Atmosphere (ISA) and World
        Meteorological Organization (WMO) standard atmosphere are identical.

        This is a revision of the US 1962 atmosphere.

        References
        ----------
        .. [1] NOAA, NASA, and USAF. "U.S. Standard Atmosphere, 1976" October 15,
           1976. http://ntrs.nasa.gov/search.jsp?R=19770009539.
        .. [2] "ISO 2533:1975 - Standard Atmosphere." ISO.
           http://www.iso.org/iso/catalogue_detail.htm?csnumber=7472.
        .. [3] Yager, Robert J. "Calculating Atmospheric Conditions (Temperature,
           Pressure, Air Density, and Speed of Sound) Using C++," June 2013.
           http://www.dtic.mil/cgi-bin/GetTRDoc?AD=ADA588839
        '''

        def __init__(self, Z, dT=0.0):
            self.Z = Z
            self.dT = dT
            self.H = r0*Z/(r0+Z)

            i = self._get_ind_from_H(self.H)
            self.T_layer = T_std[i]
            self.T_increase = T_grad[i]
            self.P_layer = P_std[i]
            self.H_layer = H_std[i]

            self.H_above_layer = self.H - self.H_layer
            self.T = self.T_layer + self.T_increase*self.H_above_layer

            R = 8314.32
            if self.T_increase == 0.0:
                self.P = self.P_layer*exp(-g0*M0*(self.H_above_layer)/(R*self.T_layer))
            else:
                self.P = self.P_layer*(self.T_layer/self.T)**(g0*M0/(R*self.T_increase))

            # Affects only the following properties
            self.T += dT

            self.rho = self.density(self.T, self.P)
            self.v_sonic = self.sonic_velocity(self.T)
            self.mu = self.viscosity(self.T)
            self.k = self.thermal_conductivity(self.T)
            self.g = self.gravity(self.Z)

        @staticmethod
        def _get_ind_from_H(H):
            r'''Method defined in the US Standard Atmosphere 1976 for determining
            the index of the layer a specified elevation is above. Levels are
            0, 11E3, 20E3, 32E3, 47E3, 51E3, 71E3, 84852 meters respectively.
            '''
            if H <= 0.0:
                return 0
            for ind, Hi in enumerate(H_std):
                if Hi >= H :
                    return ind - 1
            return 7 # case for > 84852 m.

        @staticmethod
        def thermal_conductivity(T):
            r'''Method defined in the US Standard Atmosphere 1976 for calculating
            thermal conductivity of air as a function of `T` only.

            .. math::
                k_g = \frac{2.64638\times10^{-3}T^{1.5}}
                {T + 245.4\cdot 10^{-12./T}}

            Parameters
            ----------
            T : float
                Temperature, [K]

            Returns
            -------
            kg : float
                Thermal conductivity, [W/m/K]
            '''
            # 10**(-12./T) = exp(-12*log(10)/T) = -27.63102111...
            return 2.64638E-3*T*sqrt(T)/(T + 245.4*exp(-27.63102111592855/T))

        @staticmethod
        def viscosity(T):
            r'''Method defined in the US Standard Atmosphere 1976 for calculating
            viscosity of air as a function of `T` only.

            .. math::
                \mu_g = \frac{1.458\times10^{-6}T^{1.5}}{T+110.4}

            Parameters
            ----------
            T : float
                Temperature, [K]

            Returns
            -------
            mug : float
                Viscosity, [Pa*s]
            '''
            return 1.458E-6*T*sqrt(T)/(T + 110.4)

        @staticmethod
        def density(T, P):
            r'''Method defined in the US Standard Atmosphere 1976 for calculating
            density of air as a function of `T` and `P`. MW is defined as 28.9644
            g/mol, and R as 8314.32 J/kmol/K

            .. math::
                \rho_g = \frac{P\cdot MW}{T\cdot R\cdot 1000}

            Parameters
            ----------
            T : float
                Temperature, [K]
            P : float
                Pressure, [Pa]

            Returns
            -------
            rho : float
                Mass density, [kg/m^3]
            '''
            # 0.00348367635597379 = M0/R
            return P*0.00348367635597379/T

        @staticmethod
        def sonic_velocity(T):
            r'''Method defined in the US Standard Atmosphere 1976 for calculating
            the speed of sound in air as a function of `T` only.

            .. math::
                c = \left(\frac{\gamma R T}{MW}\right)^{0.5}

            Parameters
            ----------
            T : float
                Temperature, [K]

            Returns
            -------
            c : float
                Speed of sound, [m/s]
            '''
            # 401.87... = gamma*R/MO
            return sqrt(401.87430086589046*T)

        @staticmethod
        def gravity(Z):
            r'''Method defined in the US Standard Atmosphere 1976 for calculating
            the gravitational acceleration above earth as a function of elevation
            only.

            .. math::
                g = g_0\left(\frac{r_0}{r_0+Z}\right)^2

            Parameters
            ----------
            Z : float
                Elevation above sea level, [m]

            Returns
            -------
            g : float
                Acceleration due to gravity, [m/s^2]
            '''
            x0 = (r0/(r0+Z))
            return g0*x0*x0

        @staticmethod
        def pressure_integral(T1, P1, dH):
            r'''Method to compute an integral of the pressure differential of an
            elevation difference with a base elevation defined by temperature `T1`
            and pressure `P1`. This is
            similar to subtracting the pressures at two different elevations,
            except it allows for local conditions (temperature and pressure) to be
            taken into account. This is useful for e.x. evaluating the pressure
            difference between the top and bottom of a natural draft cooling tower.

            Parameters
            ----------
            T1 : float
                Temperature at the lower elevation condition, [K]
            P1 : float
                Pressure at the lower elevation condition, [Pa]
            dH : float
                Elevation difference for which to evaluate the pressure difference,
                [m]

            Returns
            -------
            delta_P : float
                Pressure difference between the elevations, [Pa]
            '''
            # Compute the elevation to obtain the pressure specified
            H_ref = secant(H_for_P_ATMOSPHERE_1976_err, x0=10.0, low=-610.0, high=86000.0, bisection=True, args=(P1,))

            # Compute the temperature delta
            dT = T1 - ATMOSPHERE_1976(H_ref, 0.0).T

            return quad(to_int_dP_ATMOSPHERE_1976, H_ref, H_ref+dH, args=(dT,))[0]

    class ATMOSPHERE_NRLMSISE00:
        r'''NRLMSISE 00 model for calculating temperature and density of gases in
        the atmosphere, from ground level to 1000 km, as a function of time of year,
        longitude and latitude, solar activity and earth's geomagnetic disturbance.

        NRLMSISE stands for the `US Naval Research Laboratory Mass Spectrometer and
        Incoherent Scatter Radar Exosphere` model, released in 2001; see [1]_ for
        details.

        Parameters
        ----------
        Z : float
            Elevation, [m]
        latitude : float, optional
            Latitude, between -90 and 90 [degrees]
        longitude : float, optional
            Longitude, between -180 and 180 or 0 and 360, [degrees]
        day : float, optional
            Day of year, 0-366 [day]
        seconds : float, optional
            Seconds since start of day, in UT1 time; using UTC provides no loss in
            accuracy [s]
        f107 : float, optional
            Daily average 10.7 cm solar flux measurement of the strength of solar
            emissions on the 100 MHz band centered on 2800 MHz, averaged hourly; in
            sfu units, which are multiples of 10^-22 W/m^2/Hz; use 150 as a default
            [10^-22 W/m^2/Hz]
        f107_avg : float, optional
            81-day sfu average; centered on specified day if possible, otherwise
            use the previous days [10^-22 W/m^2/Hz]
        geomagnetic_disturbance_indices : list of float, optional
            List of the 7 following `Ap` indexes also known as planetary magnetic
            indexes. Has a negligible effect on the calculation. 4 is the default
            value often used for each of these values, [-]

            * Average daily `Ap`.
            * 3-hour average `Ap` centered on the current time.
            * 3-hour average `Ap` before the current time.
            * 6-hour average `Ap` before the current time.
            * 9-hour average `Ap` before the current time.
            * Average `Ap` from 12 to 33 hours before the current time, based on
              eight 3-hour average `Ap` values.
            * Average `Ap` from 36 to 57 hours before the current time, based on
              eight 3-hour average `Ap` values.

        Attributes
        ----------
        rho : float
            Mass density [kg/m^3]
        T : float
            Temperature, [K]
        P : float
            Pressure, calculated with ideal gas law [Pa]
        He_density : float
            Density of helium atoms [count/m^3]
        O_density : float
            Density of monatomic oxygen [count/m^3]
        N2_density : float
            Density of nitrogen molecules [count/m^3]
        O2_density : float
            Density of oxygen molecules [count/m^3]
        Ar_density : float
            Density of Argon atoms [count/m^3]
        H_density : float
            Density of hydrogen atoms [count/m^3]
        N_density : float
            Density of monatomic nitrogen [count/m^3]
        O_anomalous_density : float
            Density of `anomalous` oxygen; see [1]_ for details [count/m^3]
        particle_density : float
            Total density of molecules [count/m^3]
        components : list[str]
            List of species making up the atmosphere [-]
        zs : list[float]
            Mole fractions of each molecule in the atmosphere, in order of
            `components` [-]

        Examples
        --------
        >>> atmosphere = ATMOSPHERE_NRLMSISE00(1E3, 45, 45, 150)
        >>> atmosphere.T, atmosphere.rho
        (285.5440860623, 1.10190620264)

        Notes
        -----
        No full description has been published of this model; it has been defined by
        its implementation only. It was written in FORTRAN, and is accessible
        at ftp://hanna.ccmc.gsfc.nasa.gov/pub/modelweb/atmospheric/msis/nrlmsise00/

        A C port of the model by Dominik Brodowskihas become popular, and is
        available on his website: http://www.brodo.de/space/nrlmsise/.

        In 2013 Joshua Milas ported the C port to Python. This is an interface to
        his excellent port. It is a 1000-sloc model, and has
        been rigorously tested against the C version, and the online calculation
        tool available at [3]_ for parametric inputs of latitude, longitude,
        altitude, time of day and day of year.

        This model is based on measurements other than gravity; it does not provide
        a calculation method for `g`. It does not provide transport properties.

        This model takes on the order of ~2 ms.

        References
        ----------
        .. [1] Picone, J. M., A. E. Hedin, D. P. Drob, and A. C. Aikin.
           "NRLMSISE-00 Empirical Model of the Atmosphere: Statistical Comparisons
           and Scientific Issues." Journal of Geophysical Research: Space Physics
           107, no. A12 (December 1, 2002): 1468. doi:10.1029/2002JA009430.
        .. [2] Tapping, K. F. "The 10.7 Cm Solar Radio Flux (F10.7)." Space Weather
           11, no. 7 (July 1, 2013): 394-406. doi:10.1002/swe.20064.
        .. [3] Natalia Papitashvili. "NRLMSISE-00 Atmosphere Model." Accessed
           November 27, 2016. http://ccmc.gsfc.nasa.gov/modelweb/models/nrlmsise00.php.
        '''

        components = ['N2', 'O2', 'Ar', 'He', 'O', 'H', 'N']
        atrrs = ['N2_density', 'O2_density', 'Ar_density', 'He_density',
                 'O_density', 'H_density', 'N_density']
        MWs = [28.0134, 31.9988, 39.948, 4.002602, 15.9994, 1.00794, 14.0067]

        def __init__(self, Z, latitude=0.0, longitude=0.0, day=0, seconds=0.0,
                     f107=150., f107_avg=150., geomagnetic_disturbance_indices=None):
            self.Z = Z
            self.latitude = latitude
            self.longitude = longitude
            self.day = day
            self.seconds = seconds
            self.f107 = f107
            self.f107_avg = f107_avg
            self.geomagnetic_disturbance_indices = geomagnetic_disturbance_indices

            from fluids.nrlmsise00 import ap_array, gtd7, nrlmsise_flags, nrlmsise_input, nrlmsise_output
            alt = Z*1e-3
            output_obj = nrlmsise_output()
            input_obj = nrlmsise_input()
            flags = nrlmsise_flags()

            flags.switches = [0] + [1]*23

            if geomagnetic_disturbance_indices:
                aph = ap_array()
                aph.a = geomagnetic_disturbance_indices
                flags.switches[9] = -1
                input_obj.ap = geomagnetic_disturbance_indices[0]
                input_obj.ap_a = aph

            input_obj.doy = day
            input_obj.year = 0
            input_obj.sec = seconds
            input_obj.alt = alt
            input_obj.g_lat = latitude
            input_obj.g_long = longitude
            input_obj.lst = seconds/3600. + longitude/15.
            input_obj.f107A = f107_avg
            input_obj.f107 = f107
            gtd7(input_obj, flags, output_obj)

            self.He_density = output_obj.d[0]*1E6 # 1/cm^3 to 1/m^3
            self.O_density = output_obj.d[1]*1E6 # 1/cm^3 to 1/m^3
            self.N2_density = output_obj.d[2]*1E6 # 1/cm^3 to 1/m^3
            self.O2_density = output_obj.d[3]*1E6 # 1/cm^3 to 1/m^3
            self.Ar_density = output_obj.d[4]*1E6 # 1/cm^3 to 1/m^3
            self.rho = output_obj.d[5]*1000 # gram/cm^3 to kg/m^3
            self.H_density = output_obj.d[6]*1E6 # 1/cm^3 to 1/m^3
            self.N_density = output_obj.d[7]*1E6 # 1/cm^3 to 1/m^3
            self.O_anomalous_density = output_obj.d[8]*1E6 # 1/cm^3 to 1/m^3
            self.T_exospheric = output_obj.t[0]
            self.T = output_obj.t[1]

            # Calculate pressure with the ideal gas law PV = nRT with V = 1 m^3
            self.P = sum([getattr(self, a) for a in self.atrrs])*self.T*R/N_A
            # Calculate mass density with known MWs
            self.rho_calculated = sum([getattr(self, a)*MW for c, a, MW in
                                       zip(self.components, self.atrrs, self.MWs)])/(1000.*N_A)

            self.particle_density = sum(getattr(self, a) for a in self.atrrs)
            self.zs = [getattr(self, a)/self.particle_density for a in self.atrrs]

    def to_int_airmass(Z, c1, c2, angle_term, R_planet_inv, func):
        rho = func(Z)
        t1 = c2 - rho*c1
        x0 = angle_term/(1.0 + Z*R_planet_inv)
        t2 = x0*x0
        t3 = 1.0/sqrt(1.0 - t1*t2)
        return rho*t3

    def airmass(func, angle, H_max=86400.0, R_planet=6.371229E6, RI=1.000276):
        r'''Calculates mass of air per square meter in the atmosphere using a
        provided atmospheric model. The lowest air mass is calculated straight up;
        as the angle is lowered to nearer and nearer the horizon, the air mass
        increases, and can approach 40x or more the minimum airmass.

        .. math::
            m(\gamma) = \int_0^\infty \rho \left\{1 - \left[1 + 2(\text{RI}-1)
            (1-\rho/\rho_0)\right]
            \left[\frac{\cos \gamma}{(1+h/R)}\right]^2\right\}^{-1/2} dH

        Parameters
        ----------
        func : float
            Function which returns the density of the atmosphere as a function of
            elevation
        angle : float
            Degrees above the horizon (90 = straight up), [degrees]
        H_max : float, optional
            Maximum height to compute the integration up to before the contribution
            of density becomes negligible, [m]
        R_planet : float, optional
            The radius of the planet for which the integration is being performed,
            [m]
        RI : float, optional
            The refractive index of the atmosphere (air on earth at 0.7 um as
            default) assumed a constant, [-]

        Returns
        -------
        m : float
            Mass of air per square meter in the atmosphere, [kg/m^2]

        Notes
        -----
        Numerical integration via SciPy's `quad` is used to perform the
        calculation.

        Examples
        --------
        >>> airmass(lambda Z : ATMOSPHERE_1976(Z).rho, 90)
        10356.12

        References
        ----------
        .. [1] Kasten, Fritz, and Andrew T. Young. "Revised Optical Air Mass Tables
           and Approximation Formula." Applied Optics 28, no. 22 (November 15,
           1989): 4735-38. https://doi.org/10.1364/AO.28.004735.
        '''
        delta0 = RI - 1.0
        rho0_inv = 1.0/func(0.0)
        angle_term = cos(radians(angle))
        R_planet_inv = 1.0/R_planet

        c0 = delta0 + delta0
        c1 = c0*rho0_inv
        c2 = 1.0 + c0
        return quad(to_int_airmass, 0.0, 86400.0, args=(c1, c2, angle_term, R_planet_inv, func))[0]

    PVLIB_MISSING_MSG = 'The module pvlib is required for this function; install it first'

    def earthsun_distance(moment):
        r'''Calculates the distance between the earth and the sun as a function
        of date and time. Uses the Reda and Andreas (2004) model described in [1]_,
        originally incorporated into the excellent
        `pvlib library <https://github.com/pvlib/pvlib-python>`_

        Parameters
        ----------
        moment : datetime
            Time and date for the calculation, in UTC time (or GMT, which is
            almost the same thing); OR a timezone-aware datetime instance
            which will be internally converted to UTC, [-]

        Returns
        -------
        distance : float
            Distance between the center of the earth and the center of the sun,
            [m]

        Examples
        --------
        >>> from datetime import datetime, timedelta
        >>> earthsun_distance(datetime(2003, 10, 17, 13, 30, 30))
        149090925951.18338

        The distance at perihelion, which occurs at 4:21 according to this
        algorithm. The real value is 04:38 (January 2nd).

        >>> earthsun_distance(datetime(2013, 1, 2, 4, 21, 50))
        147098089490.67123

        The distance at aphelion, which occurs at 14:44 according to this
        algorithm. The real value is dead on - 14:44 (July 5).

        >>> earthsun_distance(datetime(2013, 7, 5, 14, 44, 51, 0))
        152097354414.36044

        Using a timezone-aware date:

        >>> import pytz
        >>> earthsun_distance(pytz.timezone('America/Edmonton').localize(datetime(2020, 6, 6, 10, 0, 0, 0)))
        151817805599.67142

        This has a slightly different value than the value without a timezone;
        almost 5000 km further away!

        >>> earthsun_distance(datetime(2020, 6, 6, 10, 0, 0, 0))
        151812898579.44104

        Notes
        -----
        This function is quite accurate. The difference comes from the impact of
        the moon.

        Note this function is not continuous; the sun-earth distance is not
        sufficiently accurately modeled for the change to be continuous throughout
        each day.

        References
        ----------
        .. [1] Reda, Ibrahim, and Afshin Andreas. "Solar Position Algorithm for
           Solar Radiation Applications." Solar Energy 76, no. 5 (January 1, 2004):
           577-89. https://doi.org/10.1016/j.solener.2003.12.003.
        '''
        from fluids.optional import spa
        delta_t = spa.calculate_deltat(moment.year, moment.month)
        import calendar
        unixtime = calendar.timegm(moment.utctimetuple())
        # Convert datetime object to unixtime
        return spa.earthsun_distance(unixtime, delta_t=delta_t)*au

    def solar_position(moment, latitude, longitude, Z=0.0, T=298.15, P=101325.0,
                       atmos_refract=0.5667):
        r'''Calculate the position of the sun in the sky. It is defined in terms of
        two angles - the zenith and the azimith. The azimuth tells where a sundial
        would see the sun as coming from; the zenith tells how high in the sky it
        is. The solar elevation angle is returned for convenience; it is the
        complimentary angle of the zenith.

        The sun's refraction changes how high it appears as though the sun is;
        so values are returned with an optional conversion to the apparent angle.
        This impacts only the zenith/elevation.

        Uses the Reda and Andreas (2004) model described in [1]_,
        originally incorporated into the excellent
        `pvlib library <https://github.com/pvlib/pvlib-python>`_

        Parameters
        ----------
        moment : datetime, optionally with pytz info
            Time and date for the calculation, in UTC time OR in the time zone
            of the latitude/longitude specified BUT WITH A TZINFO ATTACHED!
            Please be careful with this argument, time zones are confusing. [-]
        latitude : float
            Latitude, between -90 and 90 [degrees]
        longitude : float
            Longitude, between -180 and 180, [degrees]
        Z : float, optional
            Elevation above sea level for the solar position calculation, [m]
        T : float, optional
            Temperature of atmosphere at ground level, [K]
        P : float, optional
            Pressure of atmosphere at ground level, [Pa]
        atmos_refract : float, optional
            Atmospheric refractivity, [degrees]

        Returns
        -------
        apparent_zenith : float
            Zenith of the sun as observed from the ground based after accounting
            for atmospheric refraction, [degrees]
        zenith : float
            Actual zenith of the sun (ignores atmospheric refraction), [degrees]
        apparent_altitude : float
            Altitude of the sun as observed from the ground based after accounting
            for atmospheric refraction, [degrees]
        altitude : float
            Actual altitude of the sun (ignores atmospheric refraction), [degrees]
        azimuth : float
            The azimuth of the sun, [degrees]
        equation_of_time : float
            Equation of time - the number of seconds to be added to the day's
            mean solar time to obtain the apparent solar noon time, [seconds]

        Examples
        --------
        >>> import pytz
        >>> from datetime import datetime, timedelta

        Perth, Australia - sunrise

        >>> solar_position(pytz.timezone('Australia/Perth').localize(datetime(2020, 6, 6, 7, 10, 57)), -31.95265, 115.85742)
        [90.89617025931, 90.89617025931, -0.896170259317, -0.896170259317, 63.6016017691, 79.0711232143]

        Perth, Australia - Comparing against an online source
        https://www.suncalc.org/#/-31.9526,115.8574,9/2020.06.06/14:30/1/0

        >>> solar_position(pytz.timezone('Australia/Perth').localize(datetime(2020, 6, 6, 14, 30, 0)), -31.95265, 115.85742)
        [63.4080568623, 63.4400018158, 26.59194313766, 26.55999818417, 325.121376246, 75.7467475485]

        Perth, Australia - time input without timezone; must be converted by user to UTC!

        >>> solar_position(datetime(2020, 6, 6, 14, 30, 0) - timedelta(hours=8), -31.95265, 115.85742)
        [63.4080568623, 63.4400018158, 26.59194313766, 26.55999818417, 325.121376246, 75.7467475485]

        Sunrise occurs when the zenith is 90 degrees (Calgary, AB):

        >>> local_time = datetime(2018, 4, 15, 6, 43, 5)
        >>> local_time = pytz.timezone('America/Edmonton').localize(local_time)
        >>> solar_position(local_time, 51.0486, -114.07)[0]
        90.0005468548

        Sunset occurs when the zenith is 90 degrees (13.5 hours later in this case):

        >>> solar_position(pytz.timezone('America/Edmonton').localize(datetime(2018, 4, 15, 20, 30, 28)), 51.0486, -114.07)
        [89.999569566, 90.5410381216, 0.000430433876, -0.541038121618, 286.831378190, 6.63142952587]

        Notes
        -----
        If you were standing at the same longitude of the sun such that it was no
        further east or west than you were, the amount of angle it was south or
        north of you is the *zenith*. If it were directly overhead it would be 0°;
        a little north or south and it would be a little positive;
        near sunset or sunrise, near 90°; and at night, between 90° and 180°.

        The *solar altitude angle* is defined as 90° -`zenith`.
        Note the *elevation* angle is just another name for the *altitude* angle.

        The *azimuth* the angle in degrees that the sun is East of the North angle.
        It is positive North eastwards 0° to 360°. Other conventions may be used.

        Note that due to differences in atmospheric refractivity, estimation of
        sunset and sunrise are accuract to no more than one minute. Refraction
        conditions truly vary across the atmosphere; so characterizing it by an
        average value is limiting as well.

        References
        ----------
        .. [1] Reda, Ibrahim, and Afshin Andreas. "Solar Position Algorithm for
           Solar Radiation Applications." Solar Energy 76, no. 5 (January 1, 2004):
           577-89. https://doi.org/10.1016/j.solener.2003.12.003.
        .. [2] "Navigation - What Azimuth Description Systems Are in Use? -
           Astronomy Stack Exchange."
           https://astronomy.stackexchange.com/questions/237/what-azimuth-description-systems-are-in-use?rq=1.
        '''
        import calendar

        from fluids.optional import spa
        tt = moment.utctimetuple()
        delta_t = spa.calculate_deltat(tt.tm_year, tt.tm_mon)
        unixtime = calendar.timegm(tt)
        # Input pressure in milibar; input temperature in deg C
    #    print(dict(unixtime=unixtime, lat=latitude, lon=longitude, elev=Z,
    #                          pressure=P*1E-2, temp=T-273.15, delta_t=delta_t,
    #                          atmos_refract=atmos_refract, sst=False))
        result = spa.solar_position(unixtime, lat=latitude, lon=longitude, elev=Z,
                              pressure=P*1E-2, temp=T-273.15, delta_t=delta_t,
                              atmos_refract=atmos_refract, sst=False)
        # confirmed equation of time https://www.minasi.com/figeot.asp
        # Convert minutes to seconds; sometimes negative, sometimes positive

        result[-1] = result[-1]*60.0
        return result

    def sunrise_sunset(moment, latitude, longitude):
        r'''Calculates the times at which the sun is at sunset; sunrise; and
        halfway between sunrise and sunset (transit).

        Uses the Reda and Andreas (2004) model described in [1]_,
        originally incorporated into the excellent
        `pvlib library <https://github.com/pvlib/pvlib-python>`_

        Parameters
        ----------
        moment : datetime
            Date for the calculation; needs to contain only the year, month, and
            day; if it is timezone-aware, the return values will be localized to
            this timezone [-]
        latitude : float
            Latitude, between -90 and 90 [degrees]
        longitude : float
            Longitude, between -180 and 180, [degrees]

        Returns
        -------
        sunrise : datetime
            The time at the specified day when the sun rises **IN UTC IF MOMENT
            DOES NOT HAVE A TIMEZONE, OTHERWISE THE TIMEZONE GIVEN WITH IT**, [-]
        sunset : datetime
            The time at the specified day when the sun sets **IN UTC IF MOMENT
            DOES NOT HAVE A TIMEZONE, OTHERWISE THE TIMEZONE GIVEN WITH IT**, [-]
        transit : datetime
            The time at the specified day when the sun is at solar noon - halfway
            between sunrise and sunset **IN UTC IF MOMENT
            DOES NOT HAVE A TIMEZONE, OTHERWISE THE TIMEZONE GIVEN WITH IT**, [-]

        Examples
        --------
        >>> sunrise, sunset, transit = sunrise_sunset(datetime(2018, 4, 17),
        ... 51.0486, -114.07)
        >>> sunrise
        datetime.datetime(2018, 4, 17, 12, 36, 55, 782660)
        >>> sunset
        datetime.datetime(2018, 4, 18, 2, 34, 4, 249326)
        >>> transit
        datetime.datetime(2018, 4, 17, 19, 35, 46, 686265)

        Example with time zone:

        >>> import pytz
        >>> sunrise_sunset(pytz.timezone('America/Edmonton').localize(datetime(2018, 4, 17)), 51.0486, -114.07)
        (datetime.datetime(2018, 4, 16, 6, 39, 1, 570479, tzinfo=<DstTzInfo 'America/Edmonton' MDT-1 day, 18:00:00 DST>), datetime.datetime(2018, 4, 16, 20, 32, 25, 778162, tzinfo=<DstTzInfo 'America/Edmonton' MDT-1 day, 18:00:00 DST>), datetime.datetime(2018, 4, 16, 13, 36, 0, 386341, tzinfo=<DstTzInfo 'America/Edmonton' MDT-1 day, 18:00:00 DST>))

        Note that the year/month/day as input with a timezone, is converted to UTC
        time as well.

        Notes
        -----
        This functions takes on the order of 2 ms per calculation.

        References
        ----------
        .. [1] Reda, Ibrahim, and Afshin Andreas. "Solar Position Algorithm for
           Solar Radiation Applications." Solar Energy 76, no. 5 (January 1, 2004):
           577-89. https://doi.org/10.1016/j.solener.2003.12.003.
        '''
        import calendar

        from fluids.optional import spa
        if moment.utcoffset() is not None:
            moment_utc = moment + moment.utcoffset()
        else:
            moment_utc = moment

        delta_t = spa.calculate_deltat(moment_utc.year, moment_utc.month)
        # Strip the part of the day
        ymd_moment_utc = datetime(moment_utc.year, moment_utc.month, moment_utc.day)
        unixtime = calendar.timegm(ymd_moment_utc.utctimetuple())

        unixtime = unixtime - unixtime % (86400) # Remove the remainder of the value, rounding it to the day it is
        transit, sunrise, sunset = spa.transit_sunrise_sunset(unixtime, lat=latitude, lon=longitude, delta_t=delta_t)

        transit = datetime.utcfromtimestamp(transit)
        sunrise = datetime.utcfromtimestamp(sunrise)
        sunset = datetime.utcfromtimestamp(sunset)

        if moment.tzinfo is not None:
            sunrise = moment.tzinfo.fromutc(sunrise)
            sunset = moment.tzinfo.fromutc(sunset)
            transit = moment.tzinfo.fromutc(transit)
        return sunrise, sunset, transit

    apparent_zenith_airmass_models = {'simple', 'kasten1966', 'kastenyoung1989',
                                       'gueymard1993', 'pickering2002'}
    true_zenith_airmass_models = {'youngirvine1967', 'young1994'}

    def _get_extra_radiation_shim(datetime_or_doy, solar_constant=1366.1,
        method='spencer', epoch_year=2014, **kwargs):
        if method == 'spencer':
            if not isinstance(datetime_or_doy, (float, int)):
                dayofyear = datetime_or_doy.timetuple().tm_yday
            else:
                dayofyear = datetime_or_doy
            B = (2.*pi/365.)*(dayofyear - 1)
            RoverR0sqrd = (1.00011 + 0.034221*cos(B) + 0.00128*sin(B) +
            0.000719*cos(2.0*B) + 7.7e-05*sin(2.0*B))

            Ea = solar_constant * RoverR0sqrd
            return Ea
        from pvlib.irradiance import get_extra_radiation
        return get_extra_radiation(datetime_or_doy=datetime_or_doy,
                                  solar_constant=solar_constant,
                                  method=method,
                                  epoch_year=epoch_year,
                                  **kwargs)

    def solar_irradiation(latitude, longitude, Z, moment, surface_tilt,
                          surface_azimuth, T=None, P=None, solar_constant=1366.1,
                          atmos_refract=0.5667, albedo=0.25, linke_turbidity=None,
                          extraradiation_method='spencer',
                          airmass_model='kastenyoung1989',
                          cache=None):
        r'''Calculates the amount of solar radiation and radiation reflected back
        the atmosphere which hits a surface at a specified tilt, and facing a
        specified azimuth.

        This functions is a wrapper for the incredibly
        comprehensive `pvlib library <https://github.com/pvlib/pvlib-python>`_,
        and requires it to be installed.

        Parameters
        ----------
        latitude : float
            Latitude, between -90 and 90 [degrees]
        longitude : float
            Longitude, between -180 and 180, [degrees]
        Z : float, optional
            Elevation above sea level for the position, [m]
        moment : datetime, optionally with pytz info
            Time and date for the calculation, in UTC time OR in the time zone
            of the latitude/longitude specified BUT WITH A TZINFO ATTACHED!
            Please be careful with this argument, time zones are confusing. [-]
        surface_tilt : float
            The angle above the horizontal of the object being hit by radiation,
            [degrees]
        surface_azimuth : float
            The angle the object is facing (positive, North eastwards 0° to 360°),
            [degrees]
        T : float, optional
            Temperature of atmosphere at ground level, [K]
        P : float, optional
            Pressure of atmosphere at ground level, [Pa]
        solar_constant : float, optional
            The amount of solar radiation which reaches earth's disk (at a
            standardized distance of 1 AU); this constant is independent of
            activity or conditions on earth, but will vary throughout the sun's
            lifetime and may increase or decrease slightly due to solar activity,
            [W/m^2]
        atmos_refract : float, optional
            Atmospheric refractivity at sunrise/sunset (0.5667 deg is an often used
            value; this varies substantially and has an impact of a few minutes on
            when sunrise and sunset is), [degrees]
        albedo : float, optional
            The average amount of reflection of the terrain surrounding the object
            at quite a distance; this impacts how much sunlight reflected off the
            ground, gets reflected back off clouds, [-]
        linke_turbidity : float, optional
            The amount of pollution/water in the sky versus a perfect clear sky;
            If not specified, this will be retrieved from a historical grid;
            typical values are 3 for cloudy, and 7 for severe pollution around a
            city, [-]
        extraradiation_method : str, optional
            The specified method to calculate the effect of earth's position on the
            amount of radiation which reaches earth according to the methods
            available in the `pvlib` library, [-]
        airmass_model : str, optional
            The specified method to calculate the amount of air the sunlight
            needs to travel through to reach the earth according to the methods
            available in the `pvlib` library, [-]
        cache : dict, optional
            Dictionary to to check for values to use to skip some calculations;
            `apparent_zenith`, `zenith`, `azimuth` supported, [-]

        Returns
        -------
        poa_global : float
            The total irradiance in the plane of the surface, [W/m^2]
        poa_direct : float
            The total beam irradiance in the plane of the surface, [W/m^2]
        poa_diffuse : float
            The total diffuse irradiance in the plane of the surface, [W/m^2]
        poa_sky_diffuse : float
            The sky component of the diffuse irradiance, excluding the impact
            from the ground, [W/m^2]
        poa_ground_diffuse : float
            The ground-sky diffuse irradiance component, [W/m^2]

        Examples
        --------
        >>> import pytz
        >>> solar_irradiation(Z=1100.0, latitude=51.0486, longitude=-114.07, linke_turbidity=3,
        ... moment=pytz.timezone('America/Edmonton').localize(datetime(2018, 4, 15, 13, 43, 5)), surface_tilt=41.0,
        ... surface_azimuth=180.0)
        (1065.7621896280, 945.2656564506, 120.49653317744, 95.31535344213, 25.181179735317)

        >>> cache = {'apparent_zenith': 41.099082295767545, 'zenith': 41.11285376417578, 'azimuth': 182.5631874250523}
        >>> solar_irradiation(Z=1100.0, latitude=51.0486, longitude=-114.07,
        ... moment=pytz.timezone('America/Edmonton').localize(datetime(2018, 4, 15, 13, 43, 5)), surface_tilt=41.0,
        ... linke_turbidity=3, T=300, P=1E5,
        ... surface_azimuth=180.0, cache=cache)
        (1042.567770367, 918.237754854, 124.3300155131, 99.622865737, 24.7071497753)

        At night, there is no solar radiation and this function returns zeros:

        >>> solar_irradiation(Z=1100.0, latitude=51.0486, longitude=-114.07, linke_turbidity=3,
        ... moment=pytz.timezone('America/Edmonton').localize(datetime(2018, 4, 15, 2, 43, 5)), surface_tilt=41.0,
        ... surface_azimuth=180.0)
        (0.0, -0.0, 0.0, 0.0, 0.0)

        Notes
        -----
        The retrieval of `linke_turbidity` requires the pytables library (and
        Pandas); if it is not installed, specify a value of `linke_turbidity` to
        avoid the dependency.

        There is some redundancy of the calculated results, according to the
        following relations. The total irradiance is normally that desired for
        engineering calculations.

        poa_diffuse = poa_ground_diffuse + poa_sky_diffuse

        poa_global = poa_direct + poa_diffuse

        For a surface such as a pipe or vessel, an approach would be to split it
        into a number of rectangles and sum up the radiation absorbed by each.

        This calculation is fairly slow.

        References
        ----------
        .. [1] Will Holmgren, Calama-Consulting, Tony Lorenzo, Uwe Krien, bmu,
           DaCoEx, mayudong, et al. Pvlib/Pvlib-Python: 0.5.1. Zenodo, 2017.
           https://doi.org/10.5281/zenodo.1016425.
        '''
        # Atmospheric refraction at sunrise/sunset (0.5667 deg is an often used value)
        from fluids.optional.irradiance import get_absolute_airmass, get_relative_airmass, get_total_irradiance, ineichen

        moment_timetuple = moment.timetuple()
        moment_arg_dni = (moment_timetuple.tm_yday if
                          extraradiation_method == 'spencer' else moment)

        dni_extra = _get_extra_radiation_shim(moment_arg_dni, solar_constant=solar_constant,
                                   method=extraradiation_method,
                                   epoch_year=moment.year)

        if T is None or P is None:
            atmosphere = ATMOSPHERE_NRLMSISE00(Z=Z, latitude=latitude,
                                               longitude=longitude,
                                               day=moment_timetuple.tm_yday)
            if T is None:
                T = atmosphere.T
            if P is None:
                P = atmosphere.P

        if cache is not None and 'zenith' in cache:
            zenith = cache['zenith']
            apparent_zenith = cache['apparent_zenith']
            azimuth = cache['azimuth']
        else:
            apparent_zenith, zenith, _, _, azimuth, _ = solar_position(moment=moment,
                                                                       latitude=latitude,
                                                                       longitude=longitude,
                                                                       Z=Z, T=T, P=P,
                                                                       atmos_refract=atmos_refract)

        if linke_turbidity is None:
            try:
                import pvlib  # noqa: F401
            except:
                raise ImportError(PVLIB_MISSING_MSG)
            import pandas as pd
            from pvlib.clearsky import lookup_linke_turbidity
            linke_turbidity = float(lookup_linke_turbidity(
                pd.DatetimeIndex([moment]), latitude, longitude).values)

        if airmass_model in apparent_zenith_airmass_models:
            used_zenith = apparent_zenith
        elif airmass_model in true_zenith_airmass_models:
            used_zenith = zenith
        else:
            raise ValueError('Unrecognized airmass model')

        relative_airmass = get_relative_airmass(used_zenith, model=airmass_model)
        airmass_absolute = get_absolute_airmass(relative_airmass, pressure=P)

        ans = ineichen(apparent_zenith=apparent_zenith,
                       airmass_absolute=airmass_absolute,
                       linke_turbidity=linke_turbidity,
                       altitude=Z, dni_extra=solar_constant, perez_enhancement=True)
        ghi = ans['ghi']
        dni = ans['dni']
        dhi = ans['dhi']

    #    from pvlib.irradiance import get_total_irradiance
        ans = get_total_irradiance(surface_tilt=surface_tilt,
                          surface_azimuth=surface_azimuth,
                          solar_zenith=apparent_zenith, solar_azimuth=azimuth,
                          dni=dni, ghi=ghi, dhi=dhi, dni_extra=dni_extra,
                          airmass=airmass_absolute, albedo=albedo)
        poa_global = float(ans['poa_global'])
        poa_direct = float(ans['poa_direct'])
        poa_diffuse = float(ans['poa_diffuse'])
        poa_sky_diffuse = float(ans['poa_sky_diffuse'])
        poa_ground_diffuse = float(ans['poa_ground_diffuse'])
        return (poa_global, poa_direct, poa_diffuse, poa_sky_diffuse,
                poa_ground_diffuse)