feat(user/zeropad.sh):Use find instead of ls; use bash functions
[BK-2020-03.git] / unitproc / python / bkt-humidity
1 #!/usr/bin/env python3
2 # Desc: Humidity conversion functions from EVA-2020-02-2
3 # Ref/Attrib: https://gitlab.com/baltakatei/ninfacyzga-01
4
5 def rel_to_abs(T,P,RH):
6 """Returns absolute humidity given relative humidity.
7
8 Inputs:
9 --------
10 T : float
11 Absolute temperature in units Kelvin (K).
12 P : float
13 Total pressure in units Pascals (Pa).
14 RH : float
15 Relative humidity in units percent (%).
16
17 Output:
18 --------
19 absolute_humidity : float
20 Absolute humidity in units [kg water vapor / kg dry air].
21
22 References:
23 --------
24 1. Sonntag, D. "Advancements in the field of hygrometry". 1994. https://doi.org/10.1127/metz/3/1994/51
25 2. Green, D. "Perry's Chemical Engineers' Handbook" (8th Edition). Page "12-4". McGraw-Hill Professional Publishing. 2007.
26
27 Version: 0.0.1
28 Author: Steven Baltakatei Sandoval
29 License: GPLv3+
30 """
31
32 import math;
33
34 # Check input types
35 T = float(T);
36 P = float(P);
37 RH = float(RH);
38
39 #debug
40 # print('DEBUG:Input Temperature (K) :' + str(T));
41 # print('DEBUG:Input Pressure (Pa) :' + str(P));
42 # print('DEBUG:Input Rel. Humidity (%) :' + str(RH));
43
44 # Set constants and initial conversions
45 epsilon = 0.62198 # (molar mass of water vapor) / (molar mass of dry air)
46 t = T - 273.15; # Celsius from Kelvin
47 P_hpa = P / 100; # hectoPascals (hPa) from Pascals (Pa)
48
49 # Calculate e_w(T), saturation vapor pressure of water in a pure phase, in Pascals
50 ln_e_w = -6096*T**-1 + 21.2409642 - 2.711193*10**-2*T + 1.673952*10**-5*T**2 + 2.433502*math.log(T); # Sonntag-1994 eq 7; e_w in Pascals
51 e_w = math.exp(ln_e_w);
52 e_w_hpa = e_w / 100; # also save e_w in hectoPascals (hPa)
53 # print('DEBUG:ln_e_w:' + str(ln_e_w)); # debug
54 # print('DEBUG:e_w:' + str(e_w)); # debug
55
56 # Calculate f_w(P,T), enhancement factor for water
57 f_w = 1 + (10**-4*e_w_hpa)/(273 + t)*(((38 + 173*math.exp(-t/43))*(1 - (e_w_hpa / P_hpa))) + ((6.39 + 4.28*math.exp(-t / 107))*((P_hpa / e_w_hpa) - 1))); # Sonntag-1994 eq 22.
58 # print('DEBUG:f_w:' + str(f_w)); # debug
59
60 # Calculate e_prime_w(P,T), saturation vapor pressure of water in air-water mixture, in Pascals
61 e_prime_w = f_w * e_w; # Sonntag-1994 eq 18
62 # print('DEBUG:e_prime_w:' + str(e_prime_w)); # debug
63
64 # Calculate e_prime, vapor pressure of water in air, in Pascals
65 e_prime = (RH / 100) * e_prime_w;
66 # print('DEBUG:e_prime:' + str(e_prime)); # debug
67
68 # Calculate r, the absolute humidity, in [kg water vapor / kg dry air]
69 r = (epsilon * e_prime) / (P - e_prime);
70 # print('DEBUG:r:' + str(r)); # debug
71
72 return float(r);
73
74 def rel_to_dpt(T,P,RH):
75 """Returns dew point temperature given relative humidity.
76
77 Inputs:
78 --------
79 T : float
80 Absolute temperature in units Kelvin (K).
81 P : float
82 Total pressure in units Pascals (Pa).
83 RH : float
84 Relative humidity in units percent (%).
85
86 Output:
87 --------
88 T_d : float
89 Dew point temperature in units Kelvin (K).
90
91 References:
92 --------
93 1. Sonntag, D. "Advancements in the field of hygrometry". 1994. https://doi.org/10.1127/metz/3/1994/51
94 2. Green, D. "Perry's Chemical Engineers' Handbook" (8th Edition). Page "12-4". McGraw-Hill Professional Publishing. 2007.
95
96 Version: 0.0.1
97 Author: Steven Baltakatei Sandoval
98 License: GPLv3+
99 """
100
101 import math;
102
103 # Check input types
104 T = float(T);
105 P = float(P);
106 RH = float(RH);
107
108 #debug
109 # print('DEBUG:Input Temperature (K) :' + str(T));
110 # print('DEBUG:Input Pressure (Pa) :' + str(P));
111 # print('DEBUG:Input Rel. Humidity (%) :' + str(RH));
112
113 # Set constants and initial conversions
114 epsilon = 0.62198 # (molar mass of water vapor) / (molar mass of dry air)
115 t = T - 273.15; # Celsius from Kelvin
116 P_hpa = P / 100; # hectoPascals (hPa) from Pascals (Pa)
117
118 # Calculate e_w(T), saturation vapor pressure of water in a pure phase, in Pascals
119 ln_e_w = -6096*T**-1 + 21.2409642 - 2.711193*10**-2*T + 1.673952*10**-5*T**2 + 2.433502*math.log(T); # Sonntag-1994 eq 7; e_w in Pascals
120 e_w = math.exp(ln_e_w);
121 e_w_hpa = e_w / 100; # also save e_w in hectoPascals (hPa)
122 # print('DEBUG:ln_e_w:' + str(ln_e_w)); # debug
123 # print('DEBUG:e_w:' + str(e_w)); # debug
124
125 # Calculate f_w(P,T), enhancement factor for water
126 f_w = 1 + (10**-4*e_w_hpa)/(273 + t)*(((38 + 173*math.exp(-t/43))*(1 - (e_w_hpa / P_hpa))) + ((6.39 + 4.28*math.exp(-t / 107))*((P_hpa / e_w_hpa) - 1))); # Sonntag-1994 eq 22.
127 # print('DEBUG:f_w:' + str(f_w)); # debug
128
129 # Calculate e_prime_w(P,T), saturation vapor pressure of water in air-water mixture, in Pascals
130 e_prime_w = f_w * e_w; # Sonntag-1994 eq 18
131 # print('DEBUG:e_prime_w:' + str(e_prime_w)); # debug
132
133 # Calculate e_prime, vapor pressure of water in air, in Pascals
134 e_prime = (RH / 100) * e_prime_w;
135 # print('DEBUG:e_prime:' + str(e_prime)); # debug
136
137 n = 0; repeat_flag = True;
138 while repeat_flag == True:
139 # print('DEBUG:n:' + str(n)); # debug
140
141 # Calculate f_w_td, the enhancement factor for water at dew point temperature.
142 if n == 0:
143 f = 1.0016 + 3.15*10**-6*P_hpa - (0.074 / P_hpa); # Sonntag-1994 eq 24
144 f_w_td = f; # initial approximation
145 elif n > 0:
146 t_d_prev = float(t_d); # save previous t_d value for later comparison
147 f_w_td = 1 + (10**-4*e_w_hpa)/(273 + t_d)*(((38 + 173*math.exp(-t_d/43))*(1 - (e_w_hpa / P_hpa))) + ((6.39 + 4.28*math.exp(-t_d / 107))*((P_hpa / e_w_hpa) - 1))); # Sonntag-1994 eq 22.
148 # print('DEBUG:f_w_td:' + str(f_w_td)); # debug
149
150 # Calculate e, the vapor pressure of water in the pure phase, in Pascals
151 e = (e_prime / f_w_td); # Sonntag-1994 eq 9 and 20
152 # print('DEBUG:e:' + str(e)); # debug
153
154 # Calculate y, an intermediate dew point calculation variable
155 y = math.log(e / 611.213);
156 # print('DEBUG:y:' + str(y)); # debug
157
158 # Calculate t_d, the dew point temperature in degrees Celsius
159 t_d = 13.715*y + 8.4262*10**-1*y**2 + 1.9048*10**-2*y**3 + 7.8158*10**-3*y**4;# Sonntag-1994 eq 10
160 # print('DEBUG:t_d:' + str(t_d)); # debug
161
162 if n == 0:
163 # First run
164 repeat_flag = True;
165 else:
166 # Test t_d accuracy
167 t_d_diff = math.fabs(t_d - t_d_prev);
168 # print('DEBUG:t_d :' + str(t_d)); # debug
169 # print('DEBUG:t_d_prev:' + str(t_d_prev)); # debug
170 # print('DEBUG:t_d_diff:' + str(t_d_diff)); # debug
171 if t_d_diff < 0.01:
172 repeat_flag = False;
173 else:
174 repeat_flag = True;
175
176 # Calculate T_d, the dew point temperature in Kelvin
177 T_d = 273.15 + t_d;
178 # print('DEBUG:T_d:' + str(T_d)); # debug
179
180 if n > 100:
181 return T_d; # good enough
182
183 # update loop counter
184 n += 1;
185 return T_d;