"""
US Standard Atmosphere 1976 thermophysical model, according to
:cite:`NASA1976USStandardAtmosphere`.
"""
from datetime import datetime
import numpy as np
import numpy.ma as ma
import xarray as xr
from scipy.integrate import cumtrapz
from scipy.interpolate import interp1d
from eradiate import __version__
from ..units import unit_registry as ureg
# ------------------------------------------------------------------------------
#
# Atmospheric vertical profile data set generator
#
# ------------------------------------------------------------------------------
[docs]@ureg.wraps(ret=None, args="m", strict=False)
def make_profile(levels=ureg.Quantity(np.linspace(0.0, 1e5, 51), "m")):
r"""
Makes an atmosphere vertical profile based on the US Standard Atmosphere
1976 thermophysical model.
Parameters
----------
levels : array-like, default: np.linspace(0, 1e5, 51) * ureg.m
Levels altitudes [m]. The values must be sorted by increasing order.
Valid range: 0 to 1e6 m.
Returns
-------
Dataset
Data set holding the values of the pressure, temperature,
total number density and number densities of the individual
gas species in each layer.
Notes
-----
The pressure, temperature and number densities given in each layer of
the altitude mesh are computed at the altitude of the layers centers.
In other words, the layer's middle is taken as the altitude
representative of the whole layer. For example, in a layer with lower
and upper altitudes of 1000 and 2000 m, the thermophysical variables
are computed at the altitude of 1500 m.
"""
if np.any(levels > 1e6) or np.any(levels < 0.0):
raise ValueError("Levels altitudes must be in [0, 1e6] m.")
z_layer = (levels[:-1] + levels[1:]) / 2
# create the US76 data set
ds = create(
ureg.Quantity(z_layer, "m"),
variables=["p", "t", "n", "n_tot"],
)
# derive atmospheric thermophysical properties profile data set
thermoprops_ds = (
xr.Dataset(
data_vars={
"p": ds.p,
"t": ds.t,
"n": ds.n_tot,
"mr": ds.n / ds.n_tot,
}
)
.rename_dims({"z": "z_layer"})
.drop_vars("z")
)
thermoprops_ds.data_vars["mr"].attrs.update(
dict(standard_name="mixing_ratio", long_name="mixing ratio", units="")
)
thermoprops_ds.coords["z_layer"] = (
"z_layer",
z_layer,
dict(
standard_name="layer_altitude",
long_name="layer altitude",
units="m",
),
)
thermoprops_ds.coords["z_level"] = (
"z_level",
levels,
dict(
standard_name="level_altitude",
long_name="level altitude",
units="m",
),
)
thermoprops_ds.attrs = dict(
convention="CF-1.8",
title="U.S. Standard Atmosphere 1976",
history=f"{datetime.now().strftime('%Y-%m-%d %H:%M:%S')} - "
f"data creation - eradiate.scenes.atmosphere.us76.make_profile",
source=f"eradiate, version {__version__}",
references="U.S. Standard Atmosphere, 1976, NASA-TM-X-74335, "
"NOAA-S/T-76-1562",
)
return thermoprops_ds
# ------------------------------------------------------------------------------
#
# Constants.
# As much as possible, constants names are chosen to be as close as possible to
# the notations used in :cite:`NASA1976USStandardAtmosphere`.
#
# ------------------------------------------------------------------------------
# Boltzmann constant
K = 1.380622e-23 # [J/K]
# Molar masses of the individual species
M = {
"N2": 0.0280134,
"O2": 0.0319988,
"Ar": 0.039948,
"CO2": 0.04400995,
"Ne": 0.020183,
"He": 0.0040026,
"Kr": 0.08380,
"Xe": 0.13130,
"CH4": 0.01604303,
"H2": 0.00201594,
"O": 0.01599939,
"H": 0.00100797,
} # [kg/mol]
# Sea level mean (mixture) molar mass
M0 = 0.028964425278793997 # [kg/mol]
# Avogadro number
NA = 6.022169e23 # [mol^-1]
# Universal gas constant
R = 8.31432 # [J/(mol*K)]
# Sea level volume fractions of the gas species present below 86 km
F = {
"N2": 0.78084,
"O2": 0.209476,
"Ar": 0.00934,
"CO2": 0.000314,
"Ne": 0.00001818,
"He": 0.00000524,
"Kr": 0.00000114,
"Xe": 0.000000087,
"CH4": 0.000002,
"H2": 0.0000005,
} # [dimensionless]
# Sea level gravity
G0 = 9.80665 # [m/s^2]
# Geopotential altitudes of the layers' boundaries (below 86 km)
H = [
0.0,
11e3,
20e3,
32e3,
47e3,
51e3,
71e3,
84852.05,
] # [m]
# Temperature gradients in the seven layers (below 86 km)
LK = [
-0.0065,
0.0,
0.0010,
0.0028,
0.0,
-0.0028,
-0.0020,
] # [K/m]
# Pressure at sea level
P0 = 101325.0 # [Pa]
# Effective Earth radius
R0 = 6.356766e6 # [m]
# Temperature at sea level
T0 = 288.15 # [K]
S = 110.4 # [K]
BETA = 1.458e6 # [kg/(m*s*K^1/2)]
GAMMA = 1.40 # [dimensionless]
SIGMA = 3.65e-10 # [m]
# Thermal diffusion constants of the individual species present above 86 km
ALPHA = {
"N2": 0.0, # [dimensionless]
"O": 0.0, # [dimensionless]
"O2": 0.0, # [dimensionless]
"Ar": 0.0, # [dimensionless]
"He": -0.4, # [dimensionless]
"H": -0.25, # [dimensionless]
}
A = {
"N2": None,
"O": 6.986e20, # [m^-1*s^-1]
"O2": 4.863e20, # [m^-1*s^-1]
"Ar": 4.487e20, # [m^-1*s^-1]
"He": 1.7e21, # [m^-1*s^-1]
"H": 3.305e21, # [m^-1*s^-1]
}
B = {
"N2": None,
"O": 0.75, # [dimensionless]
"O2": 0.75, # [dimensionless]
"Ar": 0.87, # [dimensionless]
"He": 0.691, # [dimensionless]
"H": 0.5, # [dimensionless]
}
# Eddy diffusion coefficients
K_7 = 1.2e2 # [m^2/s]
K_10 = 0.0 # [m^2/s]
# Vertical transport constants of the individual species present above 86 km
Q1 = {
"O": -5.809644e-4, # [km^-3]
"O2": 1.366212e-4, # [km^-3]
"Ar": 9.434079e-5, # [km^-3]
"He": -2.457369e-4, # [km^-3]
}
Q2 = {
"O": -3.416248e-3, # [km^-3], /!\ above 97 km, Q2 = 0.
"O2": 0.0, # [km^-3]
"Ar": 0.0, # [km^-3]
"He": 0.0, # [km^-3]
}
U1 = {
"O": 56.90311, # [km]
"O2": 86.0, # [km]
"Ar": 86.0, # [km]
"He": 86.0, # [km]
}
U2 = {"O": 97.0, "O2": None, "Ar": None, "He": None} # [km]
W1 = {
"O": 2.706240e-5, # [km^-3]
"O2": 8.333333e-5, # [km^-3]
"Ar": 8.333333e-5, # [km^-3]
"He": 6.666667e-4, # [km^-3]
}
W2 = {"O": 5.008765e-4, "O2": None, "Ar": None, "He": None} # [km^-3]
# Altitudes of the levels delimiting 5 layers above 86 km
Z7 = 86.0 # [km]
Z8 = 91.0 # [km]
Z9 = 110.0 # [km]
Z10 = 120.0 # [km]
Z11 = 500.0 # [km]
Z12 = 1000.0 # [km]
# Temperature at the different levels above 86 km
T7 = 186.8673 # [K]
T9 = 240.0 # [K]
T10 = 360.0 # [K]
T11 = 999.2356 # [K]
TINF = 1000.0 # [K]
LAMBDA = 0.01875 # [km^-1]
# Temperature gradients
LK7 = 0.0 # [K/km]
LK9 = 12.0 # [K/km]
# Molecular nitrogen at altitude = Z7
N2_7 = 1.129794e20 # [m^-3]
# Atomic oxygen at altitude = Z7
O_7 = 8.6e16 # [m^-3]
# Molecular oxygen at altitude = Z7
O2_7 = 3.030898e19 # [m^-3]
# Argon at altitude = Z7
AR_7 = 1.351400e18 # [m^-3]
# Helium at altitude = Z7 (assumes typo at page 13)
HE_7 = 7.5817e14 # [m^-3]
# Hydrogen at altitude = Z7
H_11 = 8.0e10 # [m^-3]
# Vertical flux
PHI = 7.2e11 # [m^-2 * s^-1]
# List of all gas species
SPECIES = ["N2", "O2", "Ar", "CO2", "Ne", "He", "Kr", "Xe", "CH4", "H2", "O", "H"]
# List of variables computed by the model
VARIABLES = [
"t",
"p",
"n",
"n_tot",
"rho",
"mv",
"hp",
"v",
"mfp",
"f",
"cs",
"mu",
"nu",
"kt",
]
# Variables standard names with respect to the Climate and Forecast (CF)
# convention
STANDARD_NAME = {
"t": "air_temperature",
"p": "air_pressure",
"n": "number_density",
"n_tot": "air_number_density",
"rho": "air_density",
"mv": "molar_volume",
"hp": "pressure_scale_height",
"v": "mean_air_particles_speed",
"mfp": "mean_free_path",
"f": "mean_collision_frequency",
"cs": "speed_of_sound_in_air",
"mu": "air_dynamic_viscosity",
"nu": "air_kinematic_viscosity",
"kt": "air_thermal_conductivity_coefficient",
"z": "altitude",
"h": "geopotential_height",
}
# Units of relevant quantities
UNITS = {
"t": "K",
"p": "Pa",
"n": "m^-3",
"n_tot": "m^-3",
"rho": "kg/m^3",
"mv": "m^3",
"hp": "m",
"v": "m/s",
"mfp": "m",
"f": "s^-1",
"cs": "m/s",
"mu": "kg/(m*s)",
"nu": "m^2/s",
"kt": "J/m",
"z": "m",
"h": "m",
"species": "",
}
# Variables dimensions
DIMS = {
"t": "z",
"p": "z",
"n": ("species", "z"),
"n_tot": "z",
"rho": "z",
"mv": "z",
"hp": "z",
"v": "z",
"mfp": "z",
"f": "z",
"cs": "z",
"mu": "z",
"nu": "z",
"kt": "z",
}
# ------------------------------------------------------------------------------
#
# Computational functions.
# The U.S. Standard Atmosphere 1976 model divides the atmosphere into two
# altitude regions:
# 1. the low-altitude region, from 0 to 86 kilometers
# 2. the high-altitude region, from 86 to 1000 kilometers.
# The majority of computational functions hereafter are specialised for one or
# the other altitude region and is valid only in that altitude region, not in
# the other.
#
# ------------------------------------------------------------------------------
[docs]@ureg.wraps(ret=None, args=("m", None), strict=False)
def create(z, variables=None):
"""
Creates a US Standard Atmosphere 1976 data set using specified altitudes
values.
Parameters
----------
z : array-like
1-D array with altitude values [m].
variables : list, optional
Names of the variables to compute.
Returns
-------
Dataset
Data set holding the values of the different atmospheric variables.
Warnings
--------
The returned U.S. Standard Atmosphere 1976 data set is not an
atmospheric vertical profile data set. See
:func:`eradiate.thermoprops.us76.make_profile`
if you are interested in generating an atmospheric vertical profile
based on the U.S. Standard Atmosphere 1976 model.
"""
if np.any(z < 0.0):
raise ValueError("altitude values must be greater than or equal to " "zero")
if np.any(z > 1000000.0):
raise ValueError("altitude values must be less then or equal to 1e6 m")
if variables is None:
variables = VARIABLES
else:
for var in variables:
if var not in VARIABLES:
raise ValueError(var, " is not a valid variable name")
# initialise data set
ds = init_data_set(ureg.Quantity(z, "m"))
# compute the model in the low-altitude region
compute_low_altitude(ds, ds.coords["z"] <= 86000.0, inplace=True)
# compute the model in the high-altitude region
compute_high_altitude(ds, ds.coords["z"] > 86000.0, inplace=True)
# replace all np.nan with 0. in number densities values
n = ds.n.values
n[np.isnan(n)] = 0.0
ds.n.values = n
# list names of variables to drop from the data set
names = []
for var in ds.data_vars:
if var not in variables:
names.append(var)
return ds.drop_vars(names)
[docs]def compute_low_altitude(data_set, mask=None, inplace=False):
r"""
Computes the US Standard Atmosphere 1976 in the low-altitude region.
Parameters
----------
data_set : Dataset
Data set to compute.
mask : DataArray, optional
Mask to select the region of the data set to compute.
By default, the mask selects the entire data set.
inplace : bool, default: False
If true, modifies ``data_set`` in place, else returns a copy of
``data_set``.
Default: False.
Returns
-------
Dataset or None
If ``inplace`` is True, returns nothing, else returns a copy of
``data_set``.
the values of the computed variables.
"""
if mask is None:
mask = xr.full_like(data_set.coords["z"], True, dtype=bool)
if inplace:
ds = data_set
else:
ds = data_set.copy(deep=True)
altitudes = ds.coords["z"][mask]
z = altitudes.values
# compute levels temperature and pressure values
tb, pb = compute_levels_temperature_and_pressure_low_altitude()
# compute geopotential height, temperature and pressure
h = to_geopotential_height(z)
t = compute_temperature_low_altitude(h, tb)
p = compute_pressure_low_altitude(h, pb, tb)
# compute the auxiliary atmospheric variables
n_tot = NA * p / (R * t)
rho = p * M0 / (R * t)
g = compute_gravity(z)
mu = BETA * np.power(t, 1.5) / (t + S)
# assign data set with computed values
ds["t"].loc[dict(z=altitudes)] = t
ds["p"].loc[dict(z=altitudes)] = p
ds["n_tot"].loc[dict(z=altitudes)] = n_tot
species = ["N2", "O2", "Ar", "CO2", "Ne", "He", "Kr", "Xe", "CH4", "H2"]
for i, s in enumerate(SPECIES):
if s in species:
ds["n"][i].loc[dict(z=altitudes)] = F[s] * n_tot
ds["rho"].loc[dict(z=altitudes)] = rho
ds["mv"].loc[dict(z=altitudes)] = NA / n_tot
ds["hp"].loc[dict(z=altitudes)] = R * t / (g * M0)
ds["v"].loc[dict(z=altitudes)] = np.sqrt(8.0 * R * t / (np.pi * M0))
ds["mfp"].loc[dict(z=altitudes)] = np.sqrt(2.0) / (
2.0 * np.pi * np.power(SIGMA, 2.0) * n_tot
)
ds["f"].loc[dict(z=altitudes)] = (
4.0
* NA
* np.power(SIGMA, 2.0)
* np.sqrt(np.pi * np.power(p, 2.0) / (R * M0 * t))
)
ds["cs"].loc[dict(z=altitudes)] = np.sqrt(GAMMA * R * t / M0)
ds["mu"].loc[dict(z=altitudes)] = mu
ds["nu"].loc[dict(z=altitudes)] = mu / rho
ds["kt"].loc[dict(z=altitudes)] = (
2.64638e-3 * np.power(t, 1.5) / (t + 245.4 * np.power(10.0, -12.0 / t))
)
if not inplace:
return ds
[docs]def compute_high_altitude(data_set, mask=None, inplace=False):
r"""
Computes the US Standard Atmosphere 1976 in the high-altitude region.
Parameters
----------
data_set : Dataset
Data set to compute.
mask : DataArray, optional
Mask to select the region of the data set to compute.
By default, the mask selects the entire data set.
inplace : bool, default: False
If true, modifies ``data_set`` in place, else returns a copy of
``data_set``.
Default: False.
Returns
-------
Dataset or None
If ``inplace`` is True, returns nothing, else returns a copy of
``data_set``.
"""
if mask is None:
mask = xr.full_like(data_set.coords["z"], True, dtype=bool)
if inplace:
ds = data_set
else:
ds = data_set.copy(deep=True)
altitudes = ds.coords["z"][mask]
if len(altitudes) == 0:
return ds
z = ureg.Quantity(altitudes.values, "m")
n = compute_number_densities_high_altitude(z)
species = ["N2", "O", "O2", "Ar", "He", "H"]
ni = np.array([n[s] for s in species])
n_tot = np.sum(ni, axis=0)
fi = ni / n_tot[np.newaxis, :]
mi = np.array([M[s] for s in species])
m = np.sum(fi * mi[:, np.newaxis], axis=0)
t = compute_temperature_high_altitude(z)
p = K * n_tot * t
rho = np.sum(ni * mi[:, np.newaxis], axis=0) / NA
g = compute_gravity(z)
# assign data set with computed values
ds["t"].loc[dict(z=altitudes)] = t
ds["p"].loc[dict(z=altitudes)] = p
ds["n_tot"].loc[dict(z=altitudes)] = n_tot
for i, s in enumerate(SPECIES):
if s in species:
ds["n"][i].loc[dict(z=altitudes)] = n[s]
ds["rho"].loc[dict(z=altitudes)] = rho
ds["mv"].loc[dict(z=altitudes)] = NA / n_tot
ds["hp"].loc[dict(z=altitudes)] = R * t / (g * m)
ds["v"].loc[dict(z=altitudes)] = np.sqrt(8.0 * R * t / (np.pi * m))
ds["mfp"].loc[dict(z=altitudes)] = np.sqrt(2.0) / (
2.0 * np.pi * np.power(SIGMA, 2.0) * n_tot
)
ds["f"].loc[dict(z=altitudes)] = (
4.0
* NA
* np.power(SIGMA, 2.0)
* np.sqrt(np.pi * np.power(p, 2.0) / (R * m * t))
)
if not inplace:
return ds
[docs]@ureg.wraps(ret=None, args="m", strict=False)
def init_data_set(z):
r"""
Initialises the data set.
Parameters
----------
z : array
Altitudes values [m]
Returns
-------
Dataset
Initialised data set.
"""
data_vars = {}
for var in VARIABLES:
if var != "n":
try:
data_vars[var] = (
DIMS[var],
np.full(z.shape, np.nan),
{"units": UNITS[var], "standard_name": STANDARD_NAME[var]},
)
except KeyError:
data_vars[var] = (
DIMS[var],
np.full(z.shape, np.nan),
{"units": UNITS[var], "standard_name": STANDARD_NAME[var]},
)
else:
data_vars[var] = (
DIMS[var],
np.full((len(SPECIES), len(z)), np.nan),
{"units": UNITS[var], "standard_name": STANDARD_NAME["n"]},
)
coords = {"z": ("z", z, {"units": UNITS["z"]}), "species": ("species", SPECIES)}
# TODO: set function name in history field dynamically
attrs = {
"convention": "CF-1.8",
"title": "U.S. Standard Atmosphere 1976",
"history": f"{datetime.now().strftime('%Y-%m-%d %H:%M:%S')} - data creation - "
f"eradiate.scenes.atmosphere.us76.create",
"source": f"eradiate, version {__version__}",
"references": "U.S. Standard Atmosphere, 1976, NASA-TM-X-74335, NOAA-S/T-76-1562",
}
return xr.Dataset(data_vars, coords, attrs)
[docs]def compute_levels_temperature_and_pressure_low_altitude():
r"""
Computes the temperature and the pressure values at the 8 levels
of the low-altitude model.
Returns
-------
array
Levels temperature values [K].
array
Levels pressure values [Pa].
"""
tb = [T0]
pb = [P0]
for i in range(1, 8):
t_next = tb[i - 1] + LK[i - 1] * (H[i] - H[i - 1])
tb.append(t_next)
if LK[i - 1] == 0:
p_next = compute_pressure_low_altitude_zero_gradient(
H[i], H[i - 1], pb[i - 1], tb[i - 1]
)
else:
p_next = compute_pressure_low_altitude_non_zero_gradient(
H[i], H[i - 1], pb[i - 1], tb[i - 1], LK[i - 1]
)
pb.append(p_next)
return tb, pb
[docs]@ureg.wraps(ret=None, args="km", strict=False)
def compute_number_densities_high_altitude(altitudes):
r"""
Computes the number density of the individual species in the
high-altitude region.
.. note::
A uniform altitude grid is generated and used for the computation of the
integral as well as for the computation of the number densities of the
individual species. This gridded data is then interpolated at the query
``altitudes`` using a linear interpolation scheme in logarithmic space.
Parameters
----------
altitudes : array-like
Altitude value(s) [km].
Returns
-------
ndarray
Number densities of the individual species and total number density at
the given altitudes [m^-3].
The number densities of the individual species are stored in a single
2-D array where the first dimension is the gas species and the second
dimension is the altitude. The species are in the following order:
===== ======
Row Species
===== ======
0 N2
1 O
2 O2
3 Ar
4 He
5 H
===== ======
"""
# altitude grid
grid = np.concatenate(
(
np.linspace(start=Z7, stop=150.0, num=640, endpoint=False),
np.geomspace(start=150.0, stop=Z12, num=100, endpoint=True),
)
) # [km]
# pre-computed variables
m = compute_mean_molar_mass_high_altitude(grid) # [kg/mol]
g = compute_gravity(ureg.Quantity(grid, "km")) # [m / s^2]
t = compute_temperature_high_altitude(grid) # [K]
dt_dz = compute_temperature_gradient_high_altitude(grid) # [K/m]
below_115 = grid < 115.0
k = eddy_diffusion_coefficient(grid[below_115]) # [m^2/s]
n_grid = {}
# molecular nitrogen
y = m * g / (R * t) # [m^-1]
n_grid["N2"] = (
N2_7 * (T7 / t) * np.exp(-cumtrapz(y, 1e3 * grid, initial=0.0))
) # the factor 1000 is to convert km to m
# atomic oxygen
d = thermal_diffusion_coefficient(
background=n_grid["N2"][below_115], temperature=t[below_115], a=A["O"], b=B["O"]
)
y = thermal_diffusion_term_atomic_oxygen(
grid, g, t, dt_dz, d, k
) + velocity_term_atomic_oxygen(grid)
n_grid["O"] = O_7 * (T7 / t) * np.exp(-cumtrapz(y, 1e3 * grid, initial=0.0))
# molecular oxygen
d = thermal_diffusion_coefficient(
background=n_grid["N2"][below_115],
temperature=t[below_115],
a=A["O2"],
b=B["O2"],
)
y = thermal_diffusion_term("O2", grid, g, t, dt_dz, m, d, k) + velocity_term(
"O2", grid
)
n_grid["O2"] = O2_7 * (T7 / t) * np.exp(-cumtrapz(y, 1e3 * grid, initial=0.0))
# argon
background = (
n_grid["N2"][below_115] + n_grid["O"][below_115] + n_grid["O2"][below_115]
)
d = thermal_diffusion_coefficient(
background=background, temperature=t[below_115], a=A["Ar"], b=B["Ar"]
)
y = thermal_diffusion_term("Ar", grid, g, t, dt_dz, m, d, k) + velocity_term(
"Ar", grid
)
n_grid["Ar"] = AR_7 * (T7 / t) * np.exp(-cumtrapz(y, 1e3 * grid, initial=0.0))
# helium
background = (
n_grid["N2"][below_115] + n_grid["O"][below_115] + n_grid["O2"][below_115]
)
d = thermal_diffusion_coefficient(
background=background, temperature=t[below_115], a=A["He"], b=B["He"]
)
y = thermal_diffusion_term("He", grid, g, t, dt_dz, m, d, k) + velocity_term(
"He", grid
)
n_grid["He"] = HE_7 * (T7 / t) * np.exp(-cumtrapz(y, 1e3 * grid, initial=0.0))
# hydrogen
# below 500 km
mask = (grid >= 150.0) & (grid <= 500.0)
background = (
n_grid["N2"][mask]
+ n_grid["O"][mask]
+ n_grid["O2"][mask]
+ n_grid["Ar"][mask]
+ n_grid["He"][mask]
)
d = thermal_diffusion_coefficient(background, t[mask], A["H"], B["H"])
alpha = ALPHA["H"]
_tau = tau_function(grid[mask], below_500=True)
y = (PHI / d) * np.power(t[mask] / T11, 1 + alpha) * np.exp(_tau)
integral_values = cumtrapz(y[::-1], 1e3 * grid[mask][::-1], initial=0.0)
integral_values = integral_values[::-1]
n_below_500 = (
(H_11 - integral_values) * np.power(T11 / t[mask], 1 + alpha) * np.exp(-_tau)
)
# above 500 km
_tau = tau_function(grid[grid > 500.0], below_500=False)
n_above_500 = H_11 * np.power(T11 / t[grid > 500.0], 1 + alpha) * np.exp(-_tau)
n_grid["H"] = np.concatenate((n_below_500, n_above_500))
n = {
s: log_interp1d(grid, n_grid[s])(altitudes)
for s in ["N2", "O", "O2", "Ar", "He"]
}
# Below 150 km, the number density of atomic hydrogen is zero.
n["H"] = np.concatenate(
(
np.zeros(len(altitudes[altitudes < 150.0])),
log_interp1d(grid[grid >= 150.0], n_grid["H"])(
altitudes[altitudes >= 150.0]
),
)
)
return n
[docs]@ureg.wraps(ret=None, args="km", strict=False)
def compute_mean_molar_mass_high_altitude(z):
r"""
Computes the mean molar mass in the high-altitude region.
Parameters
----------
z : array
Altitude [km].
Returns
-------
ndarray
Mean molar mass [kg/mol].
"""
return np.where(z <= 100.0, M0, M["N2"])
[docs]@ureg.wraps(ret=None, args="km", strict=False)
def compute_temperature_high_altitude(altitude):
r"""
Computes the temperature in the high-altitude region.
Parameters
----------
altitude : array
Altitude values [km].
Returns
-------
array:
Temperature values [K].
"""
r0 = R0 / 1e3 # km
a = -76.3232 # K
b = -19.9429 # km
tc = 263.1905 # K
def t(z):
r"""
Compute the temperature at a given altitude.
Parameters
----------
z : float
Altitude [km].
Returns
-------
float:
Temperature [K].
"""
if Z7 <= z <= Z8:
return T7
elif Z8 < z <= Z9:
return tc + a * np.sqrt(1.0 - np.power((z - Z8) / b, 2.0))
elif Z9 < z <= Z10:
return T9 + LK9 * (z - Z9)
elif Z10 < z <= Z12:
return TINF - (TINF - T10) * np.exp(
-LAMBDA * (z - Z10) * (r0 + Z10) / (r0 + z)
)
else:
raise ValueError("altitude value is out of range")
return np.vectorize(t)(altitude)
[docs]@ureg.wraps(ret=None, args="km", strict=False)
def compute_temperature_gradient_high_altitude(altitude):
r"""
Computes the temperature gradient in the high-altitude region.
Parameters
----------
altitude : array
Altitude values [km].
Returns
-------
ndarray
Temperature gradient values [K/m].
"""
a = -76.3232 # [dimensionless]
b = -19.9429 # km
def gradient(z):
r"""
Computes the temperature gradient at a given altitude.
Parameters
----------
z : float
Altitude [km].
Returns
-------
float:
Temperature gradient [K/km].
"""
if Z7 <= z <= Z8:
return LK7
elif Z8 < z <= Z9:
return -a / b * ((z - Z8) / b) / np.sqrt(1 - np.square((z - Z8) / b))
elif Z9 < z <= Z10:
return LK9
elif Z10 < z <= Z12:
zeta = (z - Z10) * (R0 + Z10) / (R0 + z)
return (
LAMBDA
* (TINF - T10)
* np.square((R0 + Z10) / (R0 + z))
* np.exp(-LAMBDA * zeta)
)
else:
raise ValueError(f"altitude z out of range, should be in " f"[{Z7}, {Z12}]")
return np.vectorize(gradient)(altitude) / 1e3 # converts K/km to K/m
[docs]@ureg.wraps(ret=None, args=("m^-3", "K", "m^-1*s^-1", "m^-1*s^-1"), strict=False)
def thermal_diffusion_coefficient(background, temperature, a, b):
r"""
Computes the thermal diffusion coefficient values in the high-altitude
region.
Parameters
----------
n : array
Background number density values [m^-3].
Parameters
----------
t : array
Temperature values [K].
a : float
Thermal diffusion constant a [m^-1*s^-1].
b : float
Thermal diffusion constant b [m^-1*s^-1].
Returns
-------
array
Values of the thermal diffusion coefficient [m^2/s].
"""
return (a / background) * np.power(temperature / 273.15, b)
[docs]@ureg.wraps(ret=None, args="km", strict=False)
def eddy_diffusion_coefficient(z):
r"""
Computes the values of the Eddy diffusion coefficient in the high-altitude
region.
Parameters
----------
z : array
Altitude values [km].
Returns
-------
array
Eddy diffusion coefficient values [m^2/s].
Notes
-----
Valid in the altitude region :math:`86 <= z <= 150` km.
"""
return np.where(
z < 95.0, K_7, K_7 * np.exp(1.0 - (400.0 / (400.0 - np.square(z - 95.0))))
)
[docs]@ureg.wraps(
ret=None,
args=("m/s^2", "K", "K/m", "kg/mol", "kg/mol", None, "m^2/s", "m^2/s"),
strict=False,
)
def f_below_115_km(g, t, dt_dz, m, mi, alpha, d, k):
r"""
Evaluates the function :math:`f` defined by equation (36) in
:cite:`NASA1976USStandardAtmosphere` in the altitude region :math:`86
<= z <= 115` km.
Parameters
----------
g : float or array-like
Values of gravity at the different altitudes [m / s^2].
t : float or array-like
Values of temperature at the different altitudes [K].
dt_dz : float or array-like
Values of temperature gradient at the different altitudes [K/m].
m : float
Molar mass [kg/mol].
mi : float
Species molar mass [kg/mol]
alpha : float
Alpha thermal diffusion constant [dimensionless].
d : float or array-like
Values of the thermal diffusion coefficient at the different altitudes
[m^2/s].
k : float or array-like
Values of the Eddy diffusion coefficient at the different altitudes
[m^2/s].
Returns
-------
array
Values of the function f at the different altitudes [m^-1].
"""
return (g / (R * t)) * (d / (d + k)) * (mi + (m * k) / d + (alpha * R * dt_dz) / g)
[docs]@ureg.wraps(ret=None, args=("m/s^2", "K", "K/m", "kg/mol", None), strict=False)
def f_above_115_km(g, t, dt_dz, mi, alpha):
r"""
Evaluates the function :math:`f` defined by equation (36) in
:cite:`NASA1976USStandardAtmosphere` in the altitude region :math:`115 <
z <= 1000` km.
Parameters
----------
g : float or array-like
Values of gravity at the different altitudes [m/s^2].
t : float or array-like
Values of temperature at the different altitudes [K].
dt_dz : float or array-like
Values of temperature gradient at the different altitudes [K/m].
mi : float
Species molar mass [kg/mol]
alpha : float
Alpha thermal diffusion constant [dimensionless].
Returns
-------
array
Values of the function f at the different altitudes [m^-1].
"""
return (g / (R * t)) * (mi + ((alpha * R) / g) * dt_dz)
[docs]@ureg.wraps(
ret=None,
args=(None, "km", "m/s^2", "K", "K/m", "kg/mol", "m^2/s", "m^2/s"),
strict=False,
)
def thermal_diffusion_term(species, grid, g, t, dt_dz, m, d, k):
r"""
Computes the thermal diffusion term of a given species in the
high-altitude region.
Parameters
----------
species : str
Species.
grid : array
Altitude grid [km].
g : array
Values of the gravity on the altitude grid [m/s^2].
t : array
Values of the temperature on the altitude grid [K].
dt_dz : array
Values of the temperature gradient on the altitude grid [K/m].
m : array
Values of the mean molar mass on the altitude grid [kg/mol].
d : array
Values of the molecular diffusion coefficient on the altitude grid, for altitudes < 115 km [m^2/s].
k : array
Values of the eddy diffusion coefficient on the altitude grid, for altitudes < 115 km [m^2/s].
Returns
-------
array
Values of the thermal diffusion term [km^-1].
"""
fo1 = f_below_115_km(
g[grid < 115.0],
t[grid < 115.0],
dt_dz[grid < 115.0],
m[grid < 115.0],
M[species],
ALPHA[species],
d,
k,
)
fo2 = f_above_115_km(
g[grid >= 115.0],
t[grid >= 115.0],
dt_dz[grid >= 115.0],
M[species],
ALPHA[species],
)
return np.concatenate((fo1, fo2))
[docs]@ureg.wraps(ret=None, args=("km", "m/s^2", "K", "K/m", "m^2/s", "m^2/s"), strict=False)
def thermal_diffusion_term_atomic_oxygen(grid, g, t, dt_dz, d, k):
r"""
Computes the thermal diffusion term of atomic oxygen in the high-altitude
region.
Parameters
----------
grid : array
Altitude grid [km].
g : array
Values of the gravity on the altitude grid [m/s^2].
t : array
Values of the temperature on the altitude grid [K].
dt_dz : array
Values of the temperature gradient on the altitude grid [K/m].
d : array
Values of thermal diffusion coefficient on the altitude grid [m^2/s].
k : array
Values of the Eddy diffusion coefficient on the altitude grid [m^2/s].
Returns
-------
array:
Values of the thermal diffusion term [km^-1].
"""
mask1, mask2 = grid < 115.0, grid >= 115.0
x1 = f_below_115_km(
g[mask1], t[mask1], dt_dz[mask1], M["N2"], M["O"], ALPHA["O"], d, k
)
x2 = f_above_115_km(
g[grid >= 115.0], t[grid >= 115.0], dt_dz[grid >= 115.0], M["O"], ALPHA["O"]
)
return np.concatenate((x1, x2))
[docs]@ureg.wraps(
ret=None, args=("m", "km^-3", "km^-3", "km", "km", "km^-3", "km^-3"), strict=False
)
def velocity_term_hump(z, q1, q2, u1, u2, w1, w2):
r"""
Computes the transport term given by equation (37) in
:cite:`NASA1976USStandardAtmosphere`.
Parameters
----------
z : float or array-like
Altitude [km].
q1 : float
Value of the Q constant [km^-3].
q2 : float
Value of the q constant [km^-3].
u1 : float
Value of the U constant [km].
u2 : floa
Value of the u constant [km].
w1 : floa
Value of the W constant [km^-3].
w2 : floa
Value of the w constant [km^-3].
Returns
-------
float or array-like:
Values of the transport term [km^-1].
Notes
-----
Valid in the altitude region: 86 km <= z <= 150 km
"""
# @formatter:off
return (
q1 * np.square(z - u1) * np.exp(-w1 * np.power(z - u1, 3.0))
+ q2 * np.square(u2 - z) * np.exp(-w2 * np.power(u2 - z, 3.0))
) / 1e3 # the factor 1e3 converts m^-1 to km^-1
# @formatter:on
[docs]@ureg.wraps(ret=None, args=("km", "km^-3", "km^-3", "km^-3"), strict=False)
def velocity_term_no_hump(z, q1, u1, w1):
r"""
Computes the transport term given by equation (37) in
:cite:`NASA1976USStandardAtmosphere` where the second term is zero.
Parameters
----------
z : float or array-like
Altitude [km].
q1 : float
Value of the Q constant [km^-3].
u1 : float
Value of the U constant [km].
w1 : floa
Value of the W constant [km^-3].
Returns
-------
float or array-like:
Values of the transport term [km^-1].
Notes
-----
Valid in the altitude region :math:`86 <= z <= 150` km.
"""
# @formatter:off
return (
q1 * np.square(z - u1) * np.exp(-w1 * np.power(z - u1, 3.0))
) / 1e3 # the factor 1e3 converts m^-1 to km^-1
# @formatter:on
[docs]@ureg.wraps(ret=None, args=(None, "km"), strict=False)
def velocity_term(species, grid):
r"""
Computes the velocity term of a given species in the high-altitude region.
Parameters
----------
species : str
Species.
grid : array
Altitude grid [km].
Returns
-------
array:
Values of the velocity terms [km^-1].
Notes
-----
Not valid for atomic oxygen. See :func:`velocity_term_atomic_oxygen`
"""
x1 = velocity_term_no_hump(
grid[grid <= 150.0], Q1[species], U1[species], W1[species]
)
# Above 150 km, the velocity term is neglected, as indicated at p. 14 in
# :cite:`NASA1976USStandardAtmosphere`
x2 = np.zeros(len(grid[grid > 150.0]))
return np.concatenate((x1, x2))
[docs]@ureg.wraps(ret=None, args="km", strict=False)
def velocity_term_atomic_oxygen(grid):
r"""
Computes the velocity term of atomic oxygen in the high-altitude region.
Parameters
----------
grid : array
Altitude grid [km].
Returns
-------
array:
Values of the velocity term [km^-1].
"""
mask1, mask2 = grid <= 150.0, grid > 150.0
x1 = np.where(
grid[mask1] <= 97.0,
velocity_term_hump(
grid[mask1], Q1["O"], Q2["O"], U1["O"], U2["O"], W1["O"], W2["O"]
),
velocity_term_no_hump(grid[mask1], Q1["O"], U1["O"], W1["O"]),
)
x2 = np.zeros(len(grid[mask2]))
return np.concatenate((x1, x2))
[docs]@ureg.wraps(ret=None, args=("km", None), strict=False)
def tau_function(z_grid, below_500=True):
r"""
Computes the integral given by equation (40) in
:cite:`NASA1976USStandardAtmosphere` at each point of an altitude grid.
Parameters
----------
z_grid : array-like
Altitude grid (values sorted by ascending order) to use for integration [km].
below_500 : bool
True if altitudes in z_grid are lower than 500 km, False otherwise.
Returns
-------
array:
Integral evaluations [dimensionless].
Notes
-----
Valid for altitudes between 150 km and 500 km.
"""
if below_500:
z_grid = z_grid[::-1]
y = (
M["H"]
* compute_gravity(ureg.Quantity(z_grid, "km"))
/ (R * compute_temperature_high_altitude(z_grid))
) # [m^-1]
integral_values = cumtrapz(y, 1e3 * z_grid, initial=0.0) # the factor 1e3
# converts z_grid to meters
if below_500:
return integral_values[::-1]
else:
return integral_values
[docs]def log_interp1d(x, y):
"""
Computes the linear interpolation of :math:`y(x)` in logarithmic space.
Parameters
----------
x : array
1-D array of real value.
y : array
N-D array of real values. The length of y along the interpolation axis
must be equal to the length of x.
Returns
-------
callable:
Function whose call method uses interpolation to find the value of new
points.
"""
logx = np.log10(x)
logy = np.log10(y)
lin_interp = interp1d(logx, logy, kind="linear")
def log_interp(z):
return np.power(10.0, lin_interp(np.log10(z)))
return log_interp
[docs]@ureg.wraps(ret=None, args=("m", "Pa", "K"), strict=False)
def compute_pressure_low_altitude(h, pb, tb):
r"""
Computes the pressure in the low-altitude region.
Parameters
----------
h : array
Geopotential height values [m].
p_levels : array-like
Levels pressure [Pa].
Returns
-------
array:
Pressure values [Pa].
"""
# we create a mask for each layer
masks = [ma.masked_inside(h, H[i - 1], H[i]).mask for i in range(1, 8)]
# for each layer, we evaluate the pressure based on whether the
# temperature gradient is zero or not
p = np.empty(len(h))
for i, mask in enumerate(masks):
if LK[i] == 0:
p[mask] = compute_pressure_low_altitude_zero_gradient(
h[mask], H[i], pb[i], tb[i]
)
else:
p[mask] = compute_pressure_low_altitude_non_zero_gradient(
h[mask], H[i], pb[i], tb[i], LK[i]
)
return p
[docs]@ureg.wraps(ret=None, args=("m", "m", "Pa", "K"), strict=False)
def compute_pressure_low_altitude_zero_gradient(h, hb, pb, tb):
r"""
Computes the pressure in the low-altitude region when the temperature
gradient is zero.
Parameters
----------
h : float or array-like
Geopotential height [m].
hb : float or array-like
Geopotential height at the bottom of the layer [m].
pb : float or array-like
Pressure at the bottom of the layer [Pa].
tb : float or array-like
Temperature at the bottom of the layer [K].
Returns
-------
float or array-like:
Pressure [Pa].
"""
return pb * np.exp(-G0 * M0 * (h - hb) / (R * tb))
[docs]@ureg.wraps(ret=None, args=("m", "m", "Pa", "K", "K/m"), strict=False)
def compute_pressure_low_altitude_non_zero_gradient(h, hb, pb, tb, lkb):
r"""
Computes the pressure in the low-altitude region when the temperature
gradient is non-zero.
Parameters
----------
h : float or array-like
Geopotential height [m].
hb : float or array-like
Geopotential height at the bottom of the layer [m].
pb : float or array-like
Pressure at the bottom of the layer [Pa].
tb : float or array-like
Temperature at the bottom of the layer [K].
Returns
-------
float or array-like:
Pressure [Pa].
"""
return pb * np.power(tb / (tb + lkb * (h - hb)), G0 * M0 / (R * lkb))
[docs]@ureg.wraps(ret=None, args=("m", "K"), strict=False)
def compute_temperature_low_altitude(h, tb):
r"""
Computes the temperature in the low-altitude region.
Parameters
----------
h : array
Geopotential height values [m].
tb : array-like
Levels temperature values [K].
Returns
-------
array:
Temperature [K].
"""
# we create a mask for each layer
masks = [ma.masked_inside(h, H[i - 1], H[i]).mask for i in range(1, 8)]
# for each layer, we evaluate the pressure based on whether the
# temperature gradient is zero or not
t = np.empty(len(h))
for i, mask in enumerate(masks):
if LK[i] == 0:
t[mask] = tb[i]
else:
t[mask] = tb[i] + LK[i] * (h[mask] - H[i])
return t
[docs]@ureg.wraps(ret=None, args="m", strict=False)
def to_altitude(h):
r"""
Converts geopotential height to (geometric) altitude.
Parameters
----------
h : float or array-like
Geopotential altitude [m].
Returns
-------
float or array-like:
Altitude [m]
"""
return R0 * h / (R0 - h)
[docs]@ureg.wraps(ret=None, args="m", strict=False)
def to_geopotential_height(z):
r"""
Converts altitude to geopotential height.
Parameters
----------
z : float or array-like
Altitude [m].
Returns
-------
float or array-like:
Geopotential height [m]
"""
return R0 * z / (R0 + z)
[docs]@ureg.wraps(ret=None, args="m", strict=False)
def compute_gravity(z):
r"""
Computes the gravity.
Parameters
----------
z : float or array-like
Altitude [m].
Returns
-------
float or array-like:
Gravity [m/s^2].
"""
return G0 * np.power((R0 / (R0 + z)), 2.0)
