-
Notifications
You must be signed in to change notification settings - Fork 0
/
PET_FLindberg.py
357 lines (302 loc) · 10 KB
/
PET_FLindberg.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
# -*- coding: utf-8 -*-
"""
Created on Fri Jul 30 13:11:22 2021
@author: weedingb
"""
"""Adjusted code by Fredrik Lindberg, orginally from Dag Wastberg, Rambolls"""
import numpy as np
#from collections import namedtuple
#PET_person=namedtuple("PET_person","mbody age height activity sex clo")
def calculate_PET_index(Ta, Pa, Tmrt, va, pet):
mbody=pet.mbody
age=pet.age
height=pet.height
activity=pet.activity
sex=pet.sex
clo=pet.clo
pet_index=np.zeros_like(Tmrt)
for x in range(pet_index.shape[0]):
for y in range(pet_index.shape[1]):
pet_index[x][y]=_PET(Ta[x],Pa[x],Tmrt[x][y],va[x][y],mbody,age,height,activity,clo,sex)
def calculate_PET_index_vec(Ta, Pa, Tmrt, va,pet):
mbody=pet.mbody
age=pet.age
height=pet.height
activity=pet.activity
sex=pet.sex
clo=pet.clo
pet_index=_PET(Ta,Pa,Tmrt,va,mbody,age,height,activity,clo,sex)
def _PET(ta,RH,tmrt,v,mbody,age,ht,work,icl,sex):
"""
Args:
ta: air temperature
RH: relative humidity
tmrt: Mean Radiant temperature
v: wind at pedestrian heigh
mbody: body masss (kg)
age: person's age (years)
ht: height (meters)
work: activity level (W)
icl: clothing amount (0-5)
sex: 1=male 2=female
Returns:
"""
#my mod
if tmrt == np.nan:
print('Tmrt is nan')
return np.nan
# humidity conversion
vps = 6.107 * (10. ** (7.5 * ta / (238. + ta)))
vpa = RH * vps / 100 # water vapour presure, kPa
po = 1013.25 # Pressure
p = 1013.25 # Pressure
rob = 1.06
cb = 3.64 * 1000
food = 0
emsk = 0.99
emcl = 0.95
evap = 2.42e6
sigma = 5.67e-8
cair = 1.01 * 1000
eta = 0 # No idea what eta is
c_1 = 0.
c_2 = 0.
c_3 = 0.
c_4 = 0.
c_5 = 0.
c_6 = 0.
c_7 = 0.
c_8 = 0.
c_9 = 0.
c_10 = 0.
c_11 = 0.
# INBODY
metbf = 3.19 * mbody ** (3 / 4) * (1 + 0.004 * (30 - age) + 0.018 * ((ht * 100 / (mbody ** (1 / 3))) - 42.1))
metbm = 3.45 * mbody ** (3 / 4) * (1 + 0.004 * (30 - age) + 0.010 * ((ht * 100 / (mbody ** (1 / 3))) - 43.4))
if sex == 1:
met = metbm + work
else:
met = metbf + work
h = met * (1 - eta)
rtv = 1.44e-6 * met
# sensible respiration energy
tex = 0.47 * ta + 21.0
eres = cair * (ta - tex) * rtv
# latent respiration energy
vpex = 6.11 * 10 ** (7.45 * tex / (235 + tex))
erel = 0.623 * evap / p * (vpa - vpex) * rtv
# sum of the results
ere = eres + erel
# calcul constants
feff = 0.725
adu = 0.203 * mbody ** 0.425 * ht ** 0.725
facl = (-2.36 + 173.51 * icl - 100.76 * icl * icl + 19.28 * (icl ** 3)) / 100
if facl > 1:
facl = 1
rcl = (icl / 6.45) / facl
y = 1
# should these be else if statements?
if icl < 2:
y = (ht-0.2) / ht
if icl <= 0.6:
y = 0.5
if icl <= 0.3:
y = 0.1
fcl = 1 + 0.15 * icl
r2 = adu * (fcl - 1. + facl) / (2 * 3.14 * ht * y)
r1 = facl * adu / (2 * 3.14 * ht * y)
di = r2 - r1
acl = adu * facl + adu * (fcl - 1)
tcore = [0] * 8
wetsk = 0
hc = 2.67 + 6.5 * v ** 0.67
hc = hc * (p / po) ** 0.55
c_1 = h + ere
he = 0.633 * hc / (p * cair)
fec = 1 / (1 + 0.92 * hc * rcl)
htcl = 6.28 * ht * y * di / (rcl * np.log(r2 / r1) * acl)
aeff = adu * feff
c_2 = adu * rob * cb
c_5 = 0.0208 * c_2
c_6 = 0.76075 * c_2
rdsk = 0.79 * 10 ** 7
rdcl = 0
count2 = 0
j = 1
while count2 == 0 and j < 7:
tsk = 34
count1 = 0
tcl = (ta + tmrt + tsk) / 3
count3 = 1
enbal2 = 0
while count1 <= 3:
enbal = 0
while (enbal*enbal2) >= 0 and count3 < 200:
enbal2 = enbal
# 20
rclo2 = emcl * sigma * ((tcl + 273.2) ** 4 - (tmrt + 273.2) ** 4) * feff
tsk = 1 / htcl * (hc * (tcl - ta) + rclo2) + tcl
# radiation balance
rbare = aeff * (1 - facl) * emsk * sigma * ((tmrt + 273.2) ** 4 - (tsk + 273.2) ** 4)
rclo = feff * acl * emcl * sigma * ((tmrt + 273.2) ** 4 - (tcl + 273.2) ** 4)
rsum = rbare + rclo
# convection
cbare = hc * (ta - tsk) * adu * (1 - facl)
cclo = hc * (ta - tcl) * acl
csum = cbare + cclo
# core temperature
c_3 = 18 - 0.5 * tsk
c_4 = 5.28 * adu * c_3
c_7 = c_4 - c_6 - tsk * c_5
c_8 = -c_1 * c_3 - tsk * c_4 + tsk * c_6
c_9 = c_7 * c_7 - 4. * c_5 * c_8
c_10 = 5.28 * adu - c_6 - c_5 * tsk
c_11 = c_10 * c_10 - 4 * c_5 * (c_6 * tsk - c_1 - 5.28 * adu * tsk)
# tsk[tsk==36]=36.01
if tsk == 36:
tsk = 36.01
tcore[7] = c_1 / (5.28 * adu + c_2 * 6.3 / 3600) + tsk
tcore[3] = c_1 / (5.28 * adu + (c_2 * 6.3 / 3600) / (1 + 0.5 * (34 - tsk))) + tsk
if c_11 >= 0:
tcore[6] = (-c_10-c_11 ** 0.5) / (2 * c_5)
if c_11 >= 0:
tcore[1] = (-c_10+c_11 ** 0.5) / (2 * c_5)
if c_9 >= 0:
tcore[2] = (-c_7+abs(c_9) ** 0.5) / (2 * c_5)
if c_9 >= 0:
tcore[5] = (-c_7-abs(c_9) ** 0.5) / (2 * c_5)
tcore[4] = c_1 / (5.28 * adu + c_2 * 1 / 40) + tsk
# transpiration
tbody = 0.1 * tsk + 0.9 * tcore[j]
sw = 304.94 * (tbody - 36.6) * adu / 3600000
vpts = 6.11 * 10 ** (7.45 * tsk / (235. + tsk))
if tbody <= 36.6:
sw = 0
if sex == 2:
sw = 0.7 * sw
eswphy = -sw * evap
eswpot = he * (vpa - vpts) * adu * evap * fec
wetsk = eswphy / eswpot
if wetsk > 1:
wetsk = 1
eswdif = eswphy - eswpot
if eswdif <= 0:
esw = eswpot
else:
esw = eswphy
if esw > 0:
esw = 0
# diffusion
ed = evap / (rdsk + rdcl) * adu * (1 - wetsk) * (vpa - vpts)
# MAX VB
vb1 = 34 - tsk
vb2 = tcore[j] - 36.6
if vb2 < 0:
vb2 = 0
if vb1 < 0:
vb1 = 0
vb = (6.3 + 75 * vb2) / (1 + 0.5 * vb1)
# energy balance
enbal = h + ed + ere + esw + csum + rsum + food
# clothing's temperature
if count1 == 0:
xx = 1
if count1 == 1:
xx = 0.1
if count1 == 2:
xx = 0.01
if count1 == 3:
xx = 0.001
if enbal > 0:
tcl = tcl + xx
else:
tcl = tcl - xx
count3 = count3 + 1
count1 = count1 + 1
enbal2 = 0
if j == 2 or j == 5:
if c_9 >= 0:
if tcore[j] >= 36.6 and tsk <= 34.050:
if (j != 4 and vb >= 91) or (j == 4 and vb < 89):
pass
else:
if vb > 90:
vb = 90
count2 = 1
if j == 6 or j == 1:
if c_11 > 0:
if tcore[j] >= 36.6 and tsk > 33.850:
if (j != 4 and vb >= 91) or (j == 4 and vb < 89):
pass
else:
if vb > 90:
vb = 90
count2 = 1
if j == 3:
if tcore[j] < 36.6 and tsk <= 34.000:
if (j != 4 and vb >= 91) or (j == 4 and vb < 89):
pass
else:
if vb > 90:
vb = 90
count2 = 1
if j == 7:
if tcore[j] < 36.6 and tsk > 34.000:
if (j != 4 and vb >= 91) or (j == 4 and vb < 89):
pass
else:
if vb > 90:
vb = 90
count2 = 1
if j == 4:
if (j != 4 and vb >= 91) or (j == 4 and vb < 89):
pass
else:
if vb > 90:
vb = 90
count2 = 1
j = j + 1
# PET_cal
tx = ta
enbal2 = 0
count1 = 0
enbal = 0
hc = 2.67 + 6.5 * 0.1 ** 0.67
hc = hc * (p / po) ** 0.55
while count1 <= 3:
while (enbal * enbal2) >= 0:
enbal2 = enbal
# radiation balance
rbare = aeff * (1 - facl) * emsk * sigma * ((tx + 273.2) ** 4 - (tsk + 273.2) ** 4)
rclo = feff * acl * emcl * sigma * ((tx + 273.2) ** 4 - (tcl + 273.2) ** 4)
rsum = rbare + rclo
# convection
cbare = hc * (tx - tsk) * adu * (1 - facl)
cclo = hc * (tx - tcl) * acl
csum = cbare + cclo
# diffusion
ed = evap / (rdsk + rdcl) * adu * (1 - wetsk) * (12 - vpts)
# respiration
tex = 0.47 * tx + 21
eres = cair * (tx - tex) * rtv
vpex = 6.11 * 10 ** (7.45 * tex / (235 + tex))
erel = 0.623 * evap / p * (12 - vpex) * rtv
ere = eres + erel
# energy balance
enbal = h + ed + ere + esw + csum + rsum
# iteration concerning Tx
if count1 == 0:
xx = 1
if count1 == 1:
xx = 0.1
if count1 == 2:
xx = 0.01
if count1 == 3:
xx = 0.001
if enbal > 0:
tx = tx - xx
if enbal < 0:
tx = tx + xx
count1 = count1 + 1
enbal2 = 0
return tx