Skip to content
GitLab
Commit 1fcc53bb authored by Scott Griffiths's avatar Scott Griffiths
Browse files

Added astropy.units to humidity.py

parent b3a86220
Hide whitespace changes
Inline Side-by-side
humidity.py
View file @ 1fcc53bb
......@@ -12,11 +12,13 @@ Adapted from MHumidity, written by Markus Gaug <markus.gaug@uab.cat>, 04/2013
import numpy as np
from math import exp, sqrt, log10
from astropy import units as u
MOLAR_MASS_WATER_VAPOR = 0.018015 # molar mass of water vapor [kg/mol]
GAS_CONSTANT = 8.31451 # gas constant [J/mol/K]
MOLAR_MASS_WATER_VAPOR = 0.018015*u.kg/u.mol
GAS_CONSTANT = 8.31451*u.J/u.mol/u.K
@u.quantity_input(p=u.hPa, T=u.K)
def Compressibility(p, T, Xw):
"""
Calculates the compressibility of moist air, according to:
......@@ -35,21 +37,25 @@ def Compressibility(p, T, Xw):
Returns:
compressibility of moist air (dimensionless constant, 0 < Z < 1)
"""
Tc = T - 273.15 # temperature in deg. C
pT = 100*p/T # ratio of pressure to temperature, with pressure in Pascals
a0 = 1.58123E-6 # K Pa^-1
a1 = -2.9331E-8 # Pa^-1
a2 = 1.1043E-10 # K^-1 Pa^-1
b0 = 5.707E-6 # K Pa^-1
b1 = -2.051E-8 # Pa^-1
c0 = 1.9898E-4 # K Pa^-1
c1 = -2.376E-6 # Pa^-1
d0 = 1.83E-11 # K^2 Pa^-2
d1 = -7.65E-9 # K^2 Pa^-2
Z = 1 - pT*(a0 + a1*Tc + a2*Tc*Tc + b0*Xw + b1*Xw*Tc + c0*Xw*Xw + c1*Xw*Xw*Tc) + pT*pT*(d0 + d1*Xw*Xw)
Tc = T.to(u.deg_C, equivalencies=u.temperature())
pT = p.to(u.Pa)/T # ratio of pressure to temperature
a0 = 1.58123E-6*u.K*u.Pa**-1
a1 = -2.9331E-8*u.Pa**-1
a2 = 1.1043E-10*u.K**-1*u.Pa**-1
b0 = 5.707E-6*u.K*u.Pa**-1
b1 = -2.051E-8*u.Pa**-1
c0 = 1.9898E-4*u.K*u.Pa**-1
c1 = -2.376E-6*u.Pa**-1
d0 = 1.83E-11*u.K**2*u.Pa**-2
d1 = -7.65E-9*u.K**2*u.Pa**-2
Z = 1 - \
pT*(a0 + a1*Tc + a2*Tc*Tc + b0*Xw + b1*Xw*Tc + c0*Xw*Xw + c1*Xw*Xw*Tc) \
+ pT*pT*(d0 + d1*Xw*Xw)
assert Z.unit is u.dimensionless
return Z
@u.quantity_input(p=u.Pa, T=u.K)
def EnhancementFactor(p, T):
"""
Calculates the enhancement factor of water vapor in air.
......@@ -65,19 +71,22 @@ def EnhancementFactor(p, T):
Returns:
enhancement factor (dimensionless constant)
"""
Tc = T - 273.15 # temparture in deg. C
f = 1.00062 + 3.14E-8*p + 5.6E-7*pow(Tc, 2)
Tc = T.to(u.deg_C, equivalencies=u.temperature())
f = 1.00062 + (3.14E-8*u.Pa**-1)*p + (5.6E-7*u.K**-2)*pow(Tc, 2)
assert f.unit is u.dimensionless
return f
# TODO: This function needs a simple unit test to check which function is called
@u.quantity_input(T=u.K)
def SaturationVaporPressure(T):
if T > 273.15:
if T > 273.15*u.K:
return SaturationVaporPressureDavis(T)
else:
return SaturationVaporPressureGoffGratch(T)
@u.quantity_input(T=u.K)
def SaturationVaporPressureDavis(T):
"""
Calculates the vapor pressure at saturation, according to:
......@@ -94,11 +103,17 @@ def SaturationVaporPressureDavis(T):
Returns:
saturation vapor pressure in hPa
"""
E = exp(1.2378847e-5*T*T - 1.9121316e-2*T + 33.93711047 - 6343.1645/T)
return E/100 # divide by 100 for hPa
c1 = 1.2378847e-5*u.K**-2
c2 = 1.9121316e-2*u.K**-1
c3 = 33.93711047*u.dimensionless
c4 = 6343.1645*u.K
E = exp(1.2378847e-5*T*T - 1.9121316e-2*T + 33.93711047 - 6343.1645/T)*u.Pa
return E.to(u.hPa)
# TODO: Goff-Gratch doesn't seem to be as good as Buck, or maybe Wexler
@u.quantity_input(T=u.K)
def SaturationVaporPressureGoffGratch(T):
"""
Calculates the vapor pressure at saturation, according to:
......@@ -112,7 +127,7 @@ def SaturationVaporPressureGoffGratch(T):
Returns:
saturation vapor pressure in hPa
"""
theta = 373.16/T # ratio of steam point (100 deg C) to temperature
theta = (373.16*u.K)/T # ratio of steam point (100 deg C) to temperature
c = [-7.90298*(theta - 1),
5.02808*log10(theta),
-1.3816e-7*(pow(10, 11.344*(1 - 1/theta)) - 1),
......@@ -120,7 +135,7 @@ def SaturationVaporPressureGoffGratch(T):
log10(1013.246)]
log10_ew = np.sum(c)
ew = pow(10, log10_ew)
return ew
return ew*u.hPa
# TODO: Check this function and write unit test
......@@ -140,16 +155,17 @@ def SaturationVaporPressureOverWater(T):
Returns:
saturation vapor pressure in hPa
"""
omega = T - 2.38555575678E-01/(T - 6.50175348448E+02)
omega = T - 2.38555575678e-1/(T - 6.50175348448e2)
omega2 = omega*omega
A = omega2 + 1.16705214528E+03*omega - 7.24213167032E+05
B = -1.70738469401E+01*omega2 + 1.20208247025E+04*omega - 3.23255503223E+06
C = 1.49151086135E+01*omega2 - 4.82326573616E+03*omega + 4.05113405421E+05
A = omega2 + 1.16705214528e3*omega - 7.24213167032e5
B = -1.70738469401e1*omega2 + 1.20208247025e4*omega - 3.23255503223e6
C = 1.49151086135e1*omega2 - 4.82326573616e3*omega + 4.05113405421e5
X = -B + sqrt(B*B - 4*A*C)
return 1e4*pow(2*C/X, 4)
# TODO: Check this function and write unit test
@u.quantity_input(T=u.K)
def SaturationVaporPressureOverIce(T):
"""
Calculates the vapor pressure at saturation over ice, according to IAPWS:
......@@ -166,11 +182,12 @@ def SaturationVaporPressureOverIce(T):
Returns:
saturation vapor pressure in hPa
"""
theta = T/273.16
theta = T/(273.16*u.K)
Y = -13.928169*(1 - pow(theta, -1.5)) + 34.7078238*(1 - pow(theta, -1.25))
return 6.11657*exp(Y)
return 6.11657*exp(Y)*u.hPa
@u.quantity_input(p=u.hPa, T=u.K, RH=u.percent)
def MolarFractionWaterVapor(p, T, RH):
"""
Calculates the molar fraction of water vapor in moist air.
......@@ -187,12 +204,13 @@ def MolarFractionWaterVapor(p, T, RH):
Returns:
Xw : molar fraction of water vapor in moist air (dimensionless)
"""
f = EnhancementFactor(100*p, T)
f = EnhancementFactor(p.to(u.Pa), T)
psv = SaturationVaporPressure(T)
Xw = f * RH/100 * psv/p # Xw = f * h * E(T)/p
Xw = f * RH.to(u.dimensionless_unscaled) * psv/p # Xw = f * h * E(T)/p
return Xw
@u.quantity_input(p=u.hPa, T=u.K)
def DensityMoistAir(p, T, Z, Xw, C):
"""
Density equation of moist air, according to:
......@@ -208,15 +226,16 @@ def DensityMoistAir(p, T, Z, Xw, C):
Returns:
density of moist air (kg m^-3)
"""
p *= 100 # convert pressure to Pa
p = p.to(u.Pa)
R = GAS_CONSTANT
Mw = MOLAR_MASS_WATER_VAPOR
Ma = 1e-3*(28.9635 + 12.011e-6*(C - 400)) # molar mass of dry air [kg/mol]
Ma = (28.9635 + 12.011e-6*(C - 400))*1e-3*u.kg/u.mol
rho = p*Ma/(Z*R*T) * (1 - Xw*(1 - Mw/Ma)) # Tomasi eq. 12
return rho
# TODO: This function has no unit test
@u.quantity_input(T=u.K, RH=u.percent)
def PartialPressureWaterVapor(T, RH):
"""
Calculates the partial pressure of water vapor in the air.
......@@ -228,5 +247,5 @@ def PartialPressureWaterVapor(T, RH):
ew : water vapor partial pressure in hPa
"""
# water vapor pressure: e = h * E(T)
ew = (RH/100)*SaturationVaporPressureDavis(T)
ew = RH.to(u.dimensionless_unscaled)*SaturationVaporPressureDavis(T)
return ew
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment