"""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
Wind Models (requires Fortran compiler!)
----------------------------------------
.. autofunction:: hwm93
.. autofunction:: hwm14
"""
import os
from math import cos, exp, pi, radians, sin, sqrt
from fluids.constants import N_A, R, au
from fluids.numerics import brenth, quad
from fluids.numerics import numpy as np
try:
from datetime import datetime
except:
pass
__all__ = ['ATMOSPHERE_1976', 'ATMOSPHERE_NRLMSISE00', 'hwm93', 'hwm14',
'earthsun_distance', 'solar_position', 'solar_irradiation',
'sunrise_sunset']
no_gfortran_error = """This function uses f2py to encapsulate a fortran \
routine. However, f2py did not detect one on installation and could not compile \
this routine. """
try:
# Needed by hwm14
os.environ["HWMPATH"] = os.path.join(os.path.dirname(__file__), 'optional')
except:
pass
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
[docs]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.
[docs] @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))
[docs] @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)
[docs] @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
[docs] @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)
[docs] @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
[docs] @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 = brenth(H_for_P_ATMOSPHERE_1976_err, -610.0, 86000.0, 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]
[docs]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]
[docs]def hwm93(Z, latitude=0, longitude=0, day=0, seconds=0, f107=150.,
f107_avg=150., geomagnetic_disturbance_index=4):
r'''Horizontal Wind Model 1993, for calculating wind velocity in the
atmosphere as a function of time of year, longitude and latitude, solar
activity and earth's geomagnetic disturbance.
The model is described across the publications [1]_, [2]_, and [3]_.
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
[W/m^2/Hz]
f107_avg : float, optional
81-day sfu average; centered on specified day if possible, otherwise
use the previous days [W/m^2/Hz]
geomagnetic_disturbance_index : float, optional
Average daily `Ap` or also known as planetary magnetic index.
Returns
-------
v_north : float
Wind velocity, meridional (Northward) [m/s]
v_east : float
Wind velocity, zonal (Eastward) [m/s]
Examples
--------
>>> hwm93(5E5, 45, 50, 365) # doctest: +SKIP
(-73.00312042236328, 0.1485661268234253)
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/hwm93/.
F2PY auto-compilation support is not yet currently supported.
To compile this file, run the following command in a shell after navigating
to $FLUIDSPATH/fluids/optional/. This should generate the file hwm93.so
in that directory.
.. code-block:: bash
f2py -c hwm93.pyf hwm93.for --f77flags="-std=legacy"
If the module is not compiled, an import error will be raised.
References
----------
.. [1] Hedin, A. E., N. W. Spencer, and T. L. Killeen. "Empirical Global
Model of Upper Thermosphere Winds Based on Atmosphere and Dynamics
Explorer Satellite Data." Journal of Geophysical Research: Space Physics
93, no. A9 (September 1, 1988): 9959-78. doi:10.1029/JA093iA09p09959.
.. [2] Hedin, A. E., M. A. Biondi, R. G. Burnside, G. Hernandez, R. M.
Johnson, T. L. Killeen, C. Mazaudier, et al. "Revised Global Model of
Thermosphere Winds Using Satellite and Ground-Based Observations."
Journal of Geophysical Research: Space Physics 96, no. A5 (May 1, 1991):
7657-88. doi:10.1029/91JA00251.
.. [3] Hedin, A. E., E. L. Fleming, A. H. Manson, F. J. Schmidlin, S. K.
Avery, R. R. Clark, S. J. Franke, et al. "Empirical Wind Model for the
Upper, Middle and Lower Atmosphere." Journal of Atmospheric and
Terrestrial Physics 58, no. 13 (September 1996): 1421-47.
doi:10.1016/0021-9169(95)00122-0.
'''
try:
from fluids.optional.hwm93 import gws5
except: # pragma: no cover
raise ImportError(no_gfortran_error)
slt_hour = seconds/3600. + longitude/15.
ans = gws5(day, seconds, Z/1000., latitude, longitude, slt_hour, f107,
f107_avg, geomagnetic_disturbance_index)
return tuple(ans.tolist())
[docs]def hwm14(Z, latitude=0, longitude=0, day=0, seconds=0,
geomagnetic_disturbance_index=4):
r'''Horizontal Wind Model 2014, for calculating wind velocity in the
atmosphere as a function of time of year, longitude and latitude, and
earth's geomagnetic disturbance. The model is described in [1]_.
The model no longer accounts for solar flux.
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]
geomagnetic_disturbance_index : float, optional
Average daily `Ap` or also known as planetary magnetic index.
Returns
-------
v_north : float
Wind velocity, meridional (Northward) [m/s]
v_east : float
Wind velocity, zonal (Eastward) [m/s]
Examples
--------
>>> hwm14(5E5, 45, 50, 365) # doctest: +SKIP
(-38.64341354370117, 12.871272087097168)
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 http://onlinelibrary.wiley.com/store/10.1002/2014EA000089/asset/supinfo/ess224-sup-0002-supinfo.tgz?v=1&s=2a957ba70b7cf9dd0612d9430076297c3634ea75.
F2PY auto-compilation support is not yet currently supported.
To compile this file, run the following commands in a shell after navigating
to $FLUIDSPATH/fluids/optional/. This should generate the file hwm14.so
in that directory.
Generate a .pyf signature file:
.. code-block:: bash
f2py -m hwm14 -h hwm14.pyf hwm14.f90
Compile the interface:
.. code-block:: bash
f2py -c hwm14.pyf hwm14.f90
If the module is not compiled, an import error will be raised.
No patches were necessary to either the generated pyf or hwm14.f90 file,
as the authors of [1]_ have made it F2PY compatible.
Developed using 73 million data points taken by 44 instruments over 60
years.
References
----------
.. [1] Drob, Douglas P., John T. Emmert, John W. Meriwether, Jonathan J.
Makela, Eelco Doornbos, Mark Conde, Gonzalo Hernandez, et al. "An Update
to the Horizontal Wind Model (HWM): The Quiet Time Thermosphere." Earth
and Space Science 2, no. 7 (July 1, 2015): 2014EA000089.
doi:10.1002/2014EA000089.
'''
# Needed by hwm14
os.environ["HWMPATH"] = os.path.join(os.path.dirname(__file__), 'optional')
try:
try:
from fluids.optional import hwm14
except:
from optional import hwm14
except: # pragma: no cover
raise ImportError(no_gfortran_error)
ans = hwm14.hwm14(day, seconds, Z*1e-3, latitude, longitude, 0, 0,
0, np.array([np.nan, geomagnetic_disturbance_index]))
return tuple(ans.tolist())
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
[docs]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'
[docs]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
[docs]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
[docs]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)
[docs]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)