CenterForTheBuiltEnvironment / pythermalcomfort
1 9
import warnings
2 9
from pythermalcomfort.psychrometrics import (
3
    units_converter,
4
    t_o,
5
    p_sat_torr,
6
    check_standard_compliance,
7
)
8 9
import math
9 9
from scipy import optimize
10 9
from numba import jit
11

12

13 9
def cooling_effect(tdb, tr, vr, rh, met, clo, wme=0, units="SI"):
14
    """
15
    Returns the value of the Cooling Effect (`CE`_) calculated in compliance with the
16
    ASHRAE 55 2017 Standard [1]_. The `CE`_ of the elevated air speed is the value that,
17
    when subtracted equally from both the average air temperature and the mean radiant
18
    temperature, yields the same `SET`_ under still air as in the first `SET`_ calculation
19
    under elevated air speed. The cooling effect is calculated only for air speed
20
    higher than 0.1 m/s.
21

22
    .. _CE: https://en.wikipedia.org/wiki/Thermal_comfort#Cooling_Effect
23

24
    Parameters
25
    ----------
26
    tdb : float
27
        dry bulb air temperature, default in [°C] in [°F] if `units` = 'IP'
28
    tr : float
29
        mean radiant temperature, default in [°C] in [°F] if `units` = 'IP'
30
    vr : float
31
        relative air velocity, default in [m/s] in [fps] if `units` = 'IP'
32

33
        Note: vr is the relative air velocity caused by body movement and not the air
34
        speed measured by the air velocity sensor.
35
        It can be calculate using the function
36
        :py:meth:`pythermalcomfort.psychrometrics.v_relative`.
37
    rh : float
38
        relative humidity, [%]
39
    met : float
40
        metabolic rate, [met]
41
    clo : float
42
        clothing insulation, [clo]
43
    wme : float
44
        external work, [met] default 0
45
    units: str default="SI"
46
        select the SI (International System of Units) or the IP (Imperial Units) system.
47

48
    Returns
49
    -------
50
    ce
51
        Cooling Effect, default in [°C] in [°F] if `units` = 'IP'
52

53
    Examples
54
    --------
55
    .. code-block:: python
56

57
        >>> from pythermalcomfort.models import cooling_effect
58
        >>> ce = cooling_effect(tdb=25, tr=25, vr=0.3, rh=50, met=1.2, clo=0.5)
59
        >>> print(ce)
60
        1.64
61

62
        >>> # for users who wants to use the IP system
63
        >>> ce = cooling_effect(tdb=77, tr=77, vr=1.64, rh=50, met=1, clo=0.6, units="IP")
64
        >>> print(ce)
65
        3.74
66

67
    Raises
68
    ------
69
    ValueError
70
        If the cooling effect could not be calculated
71
    """
72

73 9
    if units.lower() == "ip":
74 9
        tdb, tr, vr = units_converter(tdb=tdb, tr=tr, v=vr)
75

76 9
    if vr <= 0.1:
77 9
        return 0
78

79 9
    still_air_threshold = 0.1
80

81 9
    warnings.simplefilter("ignore")
82

83 9
    initial_set_tmp = set_tmp(tdb=tdb, tr=tr, v=vr, rh=rh, met=met, clo=clo, wme=wme)
84

85 9
    def function(x):
86 9
        return (
87
            set_tmp(
88
                tdb - x, tr - x, v=still_air_threshold, rh=rh, met=met, clo=clo, wme=wme
89
            )
90
            - initial_set_tmp
91
        )
92

93 9
    try:
94 9
        ce = optimize.brentq(function, 0.0, 15)
95 0
    except ValueError:
96 0
        ce = 0
97

98 9
    warnings.simplefilter("always")
99

100 9
    if ce == 0:
101 9
        warnings.warn(
102
            "The cooling effect could not be calculated, assuming ce = 0", UserWarning
103
        )
104

105 9
    if units.lower() == "ip":
106 9
        ce = ce / 1.8 * 3.28
107

108 9
    return round(ce, 2)
109

110

111 9
def pmv_ppd(tdb, tr, vr, rh, met, clo, wme=0, standard="ISO", units="SI"):
112
    """
113
    Returns Predicted Mean Vote (`PMV`_) and Predicted Percentage of Dissatisfied (
114
    `PPD`_) calculated in accordance to main thermal comfort Standards. The `PMV`_ is
115
    an index that
116
    predicts the mean value of the thermal sensation votes (self-reported perceptions)
117
    of a large group of people on a sensation scale expressed from –3 to +3
118
    corresponding to
119
    the categories \"cold,\" \"cool,\" \"slightly cool,\" \"neutral,\" \"slightly warm,
120
    \" \"warm,\" and \"hot.\"[1]_. The `PPD`_ is an index that establishes a quantitative
121
    prediction of the percentage of thermally dissatisfied people determined from
122
    `PMV`_ [1]_.
123

124
    Parameters
125
    ----------
126
    tdb : float
127
        dry bulb air temperature, default in [°C] in [°F] if `units` = 'IP'
128
    tr : float
129
        mean radiant temperature, default in [°C] in [°F] if `units` = 'IP'
130
    vr : float
131
        relative air velocity, default in [m/s] in [fps] if `units` = 'IP'
132

133
        Note: vr is the relative air velocity caused by body movement and not the air
134
        speed measured by the air velocity sensor.
135
        The relative air velocity can be calculate using the function
136
        :py:meth:`pythermalcomfort.psychrometrics.v_relative`.
137
    rh : float
138
        relative humidity, [%]
139
    met : float
140
        metabolic rate, [met]
141
    clo : float
142
        clothing insulation, [clo]
143

144
        Note: The ASHRAE 55 Standard suggests that the dynamic clothing insulation is
145
        used as input in the PMV model.
146
        The dynamic clothing insulation can be calculated using the function
147
        :py:meth:`pythermalcomfort.psychrometrics.clo_dynamic`.
148
    wme : float
149
        external work, [met] default 0
150
    standard: str (default="ISO")
151
        comfort standard used for calculation
152

153
        - If "ISO", then the ISO Equation is used
154
        - If "ASHRAE", then the ASHRAE Equation is used
155

156
        Note: While the PMV equation is the same for both the ISO and ASHRAE standards,
157
        the ASHRAE Standard Use of the PMV model is limited to air speeds below 0.20
158
        m/s (40 fpm).
159
        When air speeds exceed 0.20 m/s (40 fpm), the comfort zone boundaries are
160
        adjusted based on the SET model.
161
    units: str default="SI"
162
        select the SI (International System of Units) or the IP (Imperial Units) system.
163

164
    Returns
165
    -------
166
    PMV
167
        Predicted Mean Vote
168
    PPD
169
        Predicted Percentage of Dissatisfied occupants, [%]
170

171
    Notes
172
    -----
173
    You can use this function to calculate the `PMV`_ and `PPD`_ in accordance with
174
    either the ASHRAE 55 2017 Standard [1]_ or the ISO 7730 Standard [2]_.
175

176
    .. _PMV: https://en.wikipedia.org/wiki/Thermal_comfort#PMV/PPD_method
177
    .. _PPD: https://en.wikipedia.org/wiki/Thermal_comfort#PMV/PPD_method
178

179
    Examples
180
    --------
181
    .. code-block:: python
182

183
        >>> from pythermalcomfort.models import pmv_ppd
184
        >>> from pythermalcomfort.psychrometrics import v_relative
185
        >>> # calculate relative air velocity
186
        >>> v_r = v_relative(v=0.1, met=1.2)
187
        >>> # as you can see the relative air velocity is 0.16 m/s which is
188
        significantly higher than v
189
        >>> results = pmv_ppd(tdb=25, tr=25, vr=v_r, rh=50, met=1.2, clo=0.5, wme=0,
190
        standard="ISO")
191
        >>> print(results)
192
        {'pmv': -0.09, 'ppd': 5.2}
193

194
        >>> print(results['pmv'])
195
        -0.09
196

197
        >>> # for users who wants to use the IP system
198
        >>> results_ip = pmv_ppd(tdb=77, tr=77, vr=0.4, rh=50, met=1.2, clo=0.5,
199
        units="IP")
200
        >>> print(results_ip)
201
        {'pmv': 0.01, 'ppd': 5.0}
202

203
    Raises
204
    ------
205
    StopIteration
206
        Raised if the number of iterations exceeds the threshold
207
    ValueError
208
        The 'standard' function input parameter can only be 'ISO' or 'ASHRAE'
209
    """
210 9
    if units.lower() == "ip":
211 9
        tdb, tr, vr = units_converter(tdb=tdb, tr=tr, v=vr)
212

213 9
    standard = standard.lower()
214 9
    if standard not in ["iso", "ashrae"]:
215 9
        raise ValueError(
216
            "PMV calculations can only be performed in compliance with ISO or ASHRAE "
217
            "Standards"
218
        )
219

220 9
    check_standard_compliance(
221
        standard=standard, tdb=tdb, tr=tr, v=vr, rh=rh, met=met, clo=clo
222
    )
223

224
    # if the relative air velocity is higher than 0.1 then follow methodology ASHRAE
225
    # Appendix H, H3
226 9
    if standard == "ashrae" and vr >= 0.1:
227
        # calculate the cooling effect
228 9
        ce = cooling_effect(tdb=tdb, tr=tr, vr=vr, rh=rh, met=met, clo=clo, wme=wme)
229

230 9
        tdb = tdb - ce
231 9
        tr = tr - ce
232 9
        vr = 0.1
233

234 9
    _pmv = pmv_ppd_optimized(tdb, tr, vr, rh, met, clo, wme)
235

236 9
    _ppd = 100.0 - 95.0 * math.exp(-0.03353 * pow(_pmv, 4.0) - 0.2179 * pow(_pmv, 2.0))
237

238 9
    return {"pmv": round(_pmv, 2), "ppd": round(_ppd, 1)}
239

240

241 9
@jit(nopython=True)
242 3
def pmv_ppd_optimized(tdb, tr, vr, rh, met, clo, wme):
243

244 0
    pa = rh * 10 * math.exp(16.6536 - 4030.183 / (tdb + 235))
245

246 0
    icl = 0.155 * clo  # thermal insulation of the clothing in M2K/W
247 0
    m = met * 58.15  # metabolic rate in W/M2
248 0
    w = wme * 58.15  # external work in W/M2
249 0
    mw = m - w  # internal heat production in the human body
250
    # calculation of the clothing area factor
251 9
    if icl <= 0.078:
252 0
        fcl = 1 + (1.29 * icl)  # ratio of surface clothed body over nude body
253
    else:
254 0
        fcl = 1.05 + (0.645 * icl)
255

256
    # heat transfer coefficient by forced convection
257 0
    hcf = 12.1 * math.sqrt(vr)
258 0
    hc = hcf  # initialize variable
259 0
    taa = tdb + 273
260 0
    tra = tr + 273
261 0
    tcla = taa + (35.5 - tdb) / (3.5 * icl + 0.1)
262

263 0
    p1 = icl * fcl
264 0
    p2 = p1 * 3.96
265 0
    p3 = p1 * 100
266 0
    p4 = p1 * taa
267 0
    p5 = (308.7 - 0.028 * mw) + (p2 * (tra / 100.0) ** 4)
268 0
    xn = tcla / 100
269 0
    xf = tcla / 50
270 0
    eps = 0.00015
271

272 0
    n = 0
273 9
    while abs(xn - xf) > eps:
274 0
        xf = (xf + xn) / 2
275 0
        hcn = 2.38 * abs(100.0 * xf - taa) ** 0.25
276 9
        if hcf > hcn:
277 0
            hc = hcf
278
        else:
279 0
            hc = hcn
280 0
        xn = (p5 + p4 * hc - p2 * xf ** 4) / (100 + p3 * hc)
281 0
        n += 1
282 9
        if n > 150:
283 0
            raise StopIteration("Max iterations exceeded")
284

285 0
    tcl = 100 * xn - 273
286

287
    # heat loss diff. through skin
288 0
    hl1 = 3.05 * 0.001 * (5733 - (6.99 * mw) - pa)
289
    # heat loss by sweating
290 9
    if mw > 58.15:
291 0
        hl2 = 0.42 * (mw - 58.15)
292
    else:
293 0
        hl2 = 0
294
    # latent respiration heat loss
295 0
    hl3 = 1.7 * 0.00001 * m * (5867 - pa)
296
    # dry respiration heat loss
297 0
    hl4 = 0.0014 * m * (34 - tdb)
298
    # heat loss by radiation
299 0
    hl5 = 3.96 * fcl * (xn ** 4 - (tra / 100.0) ** 4)
300
    # heat loss by convection
301 0
    hl6 = fcl * hc * (tcl - tdb)
302

303 0
    ts = 0.303 * math.exp(-0.036 * m) + 0.028
304 0
    _pmv = ts * (mw - hl1 - hl2 - hl3 - hl4 - hl5 - hl6)
305

306 0
    return _pmv
307

308

309 9
def pmv(tdb, tr, vr, rh, met, clo, wme=0, standard="ISO", units="SI"):
310
    """
311
    Returns Predicted Mean Vote (`PMV`_) calculated in accordance to main thermal
312
    comfort Standards. The PMV is an index that predicts the mean value of the thermal
313
    sensation votes
314
    (self-reported perceptions) of a large group of people on a sensation scale
315
    expressed from –3 to +3 corresponding to the categories \"cold,\" \"cool,
316
    \" \"slightly cool,\"
317
    \"neutral,\" \"slightly warm,\" \"warm,\" and \"hot.\" [1]_
318

319
    Parameters
320
    ----------
321
    tdb : float
322
        dry bulb air temperature, default in [°C] in [°F] if `units` = 'IP'
323
    tr : float
324
        mean radiant temperature, default in [°C] in [°F] if `units` = 'IP'
325
    vr : float
326
        relative air velocity, default in [m/s] in [fps] if `units` = 'IP'
327

328
        Note: vr is the relative air velocity caused by body movement and not the air
329
        speed measured by the air velocity sensor.
330
        It can be calculate using the function
331
        :py:meth:`pythermalcomfort.psychrometrics.v_relative`.
332
    rh : float
333
        relative humidity, [%]
334
    met : float
335
        metabolic rate, [met]
336
    clo : float
337
        clothing insulation, [clo]
338

339
        Note: The ASHRAE 55 Standard suggests that the dynamic clothing insulation is
340
        used as input in the PMV model.
341
        The dynamic clothing insulation can be calculated using the function
342
        :py:meth:`pythermalcomfort.psychrometrics.clo_dynamic`.
343
    wme : float
344
        external work, [met] default 0
345
    standard: str (default="ISO")
346
        comfort standard used for calculation
347

348
        - If "ISO", then the ISO Equation is used
349
        - If "ASHRAE", then the ASHRAE Equation is used
350

351
        Note: While the PMV equation is the same for both the ISO and ASHRAE standards,
352
        the ASHRAE Standard Use of the PMV model is limited to air speeds below 0.20
353
        m/s (40 fpm).
354
        When air speeds exceed 0.20 m/s (40 fpm), the comfort zone boundaries are
355
        adjusted based on the SET model.
356
        See ASHRAE 55 2017 Appendix H for more information [1]_.
357
    units: str default="SI"
358
        select the SI (International System of Units) or the IP (Imperial Units) system.
359

360
    Returns
361
    -------
362
    PMV : float
363
        Predicted Mean Vote
364

365
    Notes
366
    -----
367
    You can use this function to calculate the `PMV`_ [1]_ [2]_.
368

369
    .. _PMV: https://en.wikipedia.org/wiki/Thermal_comfort#PMV/PPD_method
370

371
    Examples
372
    --------
373
    .. code-block:: python
374

375
        >>> from pythermalcomfort.models import pmv
376
        >>> from pythermalcomfort.psychrometrics import v_relative
377
        >>> # calculate relative air velocity
378
        >>> v_r = v_relative(v=0.1, met=1.2)
379
        >>> # as you can see the relative air velocity is 0.16 m/s which is
380
        significantly higher than v
381
        >>> results = pmv(tdb=25, tr=25, vr=v_r, rh=50, met=1.2, clo=0.5)
382
        >>> print(results)
383
        -0.09
384
    """
385

386 9
    return pmv_ppd(tdb, tr, vr, rh, met, clo, wme, standard=standard, units=units)[
387
        "pmv"
388
    ]
389

390

391 9
def set_tmp(
392
    tdb, tr, v, rh, met, clo, wme=0, body_surface_area=1.8258, patm=101325, units="SI",
393
):
394
    """
395
    Calculates the Standard Effective Temperature (SET). The SET is the temperature of
396
    an imaginary environment at 50% (rh), <0.1 m/s (20 fpm) average air speed (v),
397
    and tr = tdb ,
398
    in which the total heat loss from the skin of an imaginary occupant with an
399
    activity level of 1.0 met and a clothing level of 0.6 clo is the same as that
400
    from a person in the actual environment with actual clothing and activity level.
401

402
    Parameters
403
    ----------
404
    tdb : float
405
        dry bulb air temperature, default in [°C] in [°F] if `units` = 'IP'
406
    tr : float
407
        mean radiant temperature, default in [°C] in [°F] if `units` = 'IP'
408
    v : float
409
        air velocity, default in [m/s] in [fps] if `units` = 'IP'
410
    rh : float
411
        relative humidity, [%]
412
    met : float
413
        metabolic rate, [met]
414
    clo : float
415
        clothing insulation, [clo]
416
    wme : float
417
        external work, [met] default 0
418
    body_surface_area : float
419
        body surface area, default value 1.8258 [m2] in [ft2] if `units` = 'IP'
420
    patm : float
421
        atmospheric pressure, default value 101325 [Pa] in [atm] if `units` = 'IP'
422
    units: str default="SI"
423
        select the SI (International System of Units) or the IP (Imperial Units) system.
424

425
    Returns
426
    -------
427
    SET : float
428
        Standard effective temperature, [°C]
429

430
    Notes
431
    -----
432
    You can use this function to calculate the `SET`_ temperature in accordance with
433
    the ASHRAE 55 2017 Standard [1]_.
434

435
    .. _SET: https://en.wikipedia.org/wiki/Thermal_comfort#Standard_effective_temperature
436

437
    Examples
438
    --------
439
    .. code-block:: python
440

441
        >>> from pythermalcomfort.models import set_tmp
442
        >>> set_tmp(tdb=25, tr=25, v=0.1, rh=50, met=1.2, clo=.5)
443
        25.3
444

445
        >>> # for users who wants to use the IP system
446
        >>> set_tmp(tdb=77, tr=77, v=0.328, rh=50, met=1.2, clo=.5, units='IP')
447
        77.6
448

449
    """
450 9
    if units.lower() == "ip":
451 9
        if body_surface_area == 1.8258:
452 9
            body_surface_area = 19.65
453 9
        if patm == 101325:
454 9
            patm = 1
455 9
        tdb, tr, v, body_surface_area, patm = units_converter(
456
            tdb=tdb, tr=tr, v=v, area=body_surface_area, pressure=patm
457
        )
458

459 9
    check_standard_compliance(
460
        standard="ashrae", tdb=tdb, tr=tr, v=v, rh=rh, met=met, clo=clo
461
    )
462

463 9
    vapor_pressure = rh * p_sat_torr(tdb) / 100
464

465 9
    _set = set_optimized(
466
        tdb, tr, v, rh, met, clo, vapor_pressure, wme, body_surface_area, patm
467
    )
468

469 9
    if units.lower() == "ip":
470 9
        _set = units_converter(tmp=_set, from_units="si")[0]
471

472 9
    return round(_set, 1)
473

474

475 9
@jit(nopython=True)
476 3
def set_optimized(
477
    tdb, tr, v, rh, met, clo, vapor_pressure, wme, body_surface_area, patm,
478
):
479
    # Initial variables as defined in the ASHRAE 55-2017
480 0
    air_velocity = max(v, 0.1)
481 0
    k_clo = 0.25
482 0
    body_weight = 69.9
483 0
    met_factor = 58.2
484 0
    sbc = 0.000000056697  # Stefan-Boltzmann constant (W/m2K4)
485 0
    c_sw = 170  # driving coefficient for regulatory sweating
486 0
    c_dil = 120  # driving coefficient for vasodilation
487 0
    c_str = 0.5  # driving coefficient for vasoconstriction
488

489 0
    temp_skin_neutral = 33.7
490 0
    temp_core_neutral = 36.8
491 0
    temp_body_neutral = 36.49
492 0
    skin_blood_flow_neutral = 6.3
493

494 0
    temp_skin = temp_skin_neutral
495 0
    temp_core = temp_core_neutral
496 0
    skin_blood_flow = skin_blood_flow_neutral
497 0
    alfa = 0.1  # fractional skin mass
498 0
    e_sk = 0.1 * met  # total evaporative heat loss, W
499

500 0
    pressure_in_atmospheres = patm / 101325
501 0
    length_time_simulation = 60  # length time simulation
502 0
    r_clo = 0.155 * clo  # thermal resistance of clothing, °C M^2 /W
503

504 0
    f_a_cl = 1.0 + 0.15 * clo  # increase in body surface area due to clothing
505 0
    lr = 2.2 / pressure_in_atmospheres  # Lewis ratio
506 0
    rm = met * met_factor  # metabolic rate
507 0
    m = met * met_factor
508

509 9
    if clo <= 0:
510 0
        w_crit = 0.38 * pow(air_velocity, -0.29)  # evaporative efficiency
511 0
        i_cl = 1.0  # thermal resistance of clothing, clo
512
    else:
513 0
        w_crit = 0.59 * pow(air_velocity, -0.08)
514 0
        i_cl = 0.45
515

516
    # h_cc corrected convective heat transfer coefficient
517 0
    h_cc = 3.0 * pow(pressure_in_atmospheres, 0.53)
518
    # h_fc forced convective heat transfer coefficient, W/(m2 °C)
519 0
    h_fc = 8.600001 * pow((air_velocity * pressure_in_atmospheres), 0.53)
520 0
    h_cc = max(h_cc, h_fc)
521

522 0
    c_hr = 4.7  # linearized radiative heat transfer coefficient
523 0
    h_t = (
524
        c_hr + h_cc
525
    )  # sum of convective and radiant heat transfer coefficient W/(m2*K)
526 0
    r_a = 1.0 / (f_a_cl * h_t)  # resistance of air layer to dry heat
527 0
    t_op = (c_hr * tr + h_cc * tdb) / h_t  # operative temperature
528

529
    # initialize some variables
530 0
    dry = 0
531 0
    p_wet = 0
532 0
    _set = 0
533

534 0
    n_simulation = 0
535

536 9
    while n_simulation < length_time_simulation:
537

538 0
        n_simulation += 1
539

540 0
        iteration_limit = 150
541
        # t_cl temperature of the outer surface of clothing
542 0
        t_cl = (r_a * temp_skin + r_clo * t_op) / (r_a + r_clo)  # initial guess
543 0
        n_iterations = 0
544 0
        tc_converged = False
545

546 9
        while not tc_converged:
547

548 0
            c_hr = 4.0 * sbc * ((t_cl + tr) / 2.0 + 273.15) ** 3.0 * 0.72
549 0
            h_t = c_hr + h_cc
550 0
            r_a = 1.0 / (f_a_cl * h_t)
551 0
            t_op = (c_hr * tr + h_cc * tdb) / h_t
552 0
            t_cl_new = (r_a * temp_skin + r_clo * t_op) / (r_a + r_clo)
553 9
            if abs(t_cl_new - t_cl) <= 0.01:
554 0
                tc_converged = True
555 0
            t_cl = t_cl_new
556 0
            n_iterations += 1
557

558 9
            if n_iterations > iteration_limit:
559 0
                raise StopIteration("Max iterations exceeded")
560

561 0
        dry = (temp_skin - t_op) / (r_a + r_clo)  # total sensible heat loss, W
562
        # h_fcs rate of energy transport between core and skin, W
563 0
        h_fcs = (temp_core - temp_skin) * (5.28 + 1.163 * skin_blood_flow)
564 0
        q_res = (
565
            0.0023 * m * (44.0 - vapor_pressure)
566
        )  # latent heat loss due to respiration
567 0
        c_res = (
568
            0.0014 * m * (34.0 - tdb)
569
        )  # rate of convective heat loss from respiration, W/m2
570
        # todo maybe the iteration should stop when the following two terms are 0
571 0
        s_core = m - h_fcs - q_res - c_res - wme  # rate of energy storage in the core
572 0
        s_skin = h_fcs - dry - e_sk  # rate of energy storage in the skin
573 0
        TCSK = 0.97 * alfa * body_weight
574 0
        TCCR = 0.97 * (1 - alfa) * body_weight
575 0
        d_t_sk = (s_skin * body_surface_area) / (
576
            TCSK * 60.0
577
        )  # rate of change skin temperature °C per minute
578 0
        d_t_cr = (
579
            s_core * body_surface_area / (TCCR * 60.0)
580
        )  # rate of change core temperature °C per minute
581 0
        temp_skin = temp_skin + d_t_sk
582 0
        temp_core = temp_core + d_t_cr
583 0
        t_body = alfa * temp_skin + (1 - alfa) * temp_core  # mean body temperature, °C
584
        # sk_sig thermoregulatory control signal from the skin
585 0
        sk_sig = temp_skin - temp_skin_neutral
586 0
        warms = (sk_sig > 0) * sk_sig  # vasodialtion signal
587 0
        colds = ((-1.0 * sk_sig) > 0) * (-1.0 * sk_sig)  # vasoconstriction signal
588
        # c_reg_sig thermoregulatory control signal from the skin, °C
589 0
        c_reg_sig = temp_core - temp_core_neutral
590
        # c_warm vasodilation signal
591 0
        c_warm = (c_reg_sig > 0) * c_reg_sig
592
        # c_cold vasoconstriction signal
593 0
        c_cold = ((-1.0 * c_reg_sig) > 0) * (-1.0 * c_reg_sig)
594 0
        body_signal = t_body - temp_body_neutral
595 0
        WARMB = (body_signal > 0) * body_signal
596 0
        skin_blood_flow = (skin_blood_flow_neutral + c_dil * c_warm) / (
597
            1 + c_str * colds
598
        )
599 9
        if skin_blood_flow > 90.0:
600 0
            skin_blood_flow = 90.0
601 9
        if skin_blood_flow < 0.5:
602 0
            skin_blood_flow = 0.5
603 0
        regulatory_sweating = c_sw * WARMB * math.exp(warms / 10.7)
604 9
        if regulatory_sweating > 500.0:
605 0
            regulatory_sweating = 500.0
606 0
        e_rsw = 0.68 * regulatory_sweating  # heat lost by vaporization sweat
607 0
        r_ea = 1.0 / (lr * f_a_cl * h_cc)  # evaporative resistance air layer
608 0
        r_ecl = r_clo / (lr * i_cl)
609
        # e_max = maximum evaporative capacity
610 0
        e_max = (
611
            math.exp(18.6686 - 4030.183 / (temp_skin + 235.0)) - vapor_pressure
612
        ) / (r_ea + r_ecl)
613 0
        p_rsw = e_rsw / e_max  # ratio heat loss sweating to max heat loss sweating
614 0
        p_wet = 0.06 + 0.94 * p_rsw  # skin wetness
615 0
        e_diff = p_wet * e_max - e_rsw  # vapor diffusion through skin
616 9
        if p_wet > w_crit:
617 0
            p_wet = w_crit
618 0
            p_rsw = w_crit / 0.94
619 0
            e_rsw = p_rsw * e_max
620 0
            e_diff = 0.06 * (1.0 - p_rsw) * e_max
621 9
        if e_max < 0:
622 0
            e_diff = 0
623 0
            e_rsw = 0
624 0
            p_wet = w_crit
625 0
        e_sk = (
626
            e_rsw + e_diff
627
        )  # total evaporative heat loss sweating and vapor diffusion
628 0
        met_shivering = 19.4 * colds * c_cold  # met shivering W/m2
629 0
        m = rm + met_shivering
630 0
        alfa = 0.0417737 + 0.7451833 / (skin_blood_flow + 0.585417)
631

632 0
    hsk = dry + e_sk  # total heat loss from skin, W
633 0
    W = p_wet
634 0
    PSSK = math.exp(18.6686 - 4030.183 / (temp_skin + 235.0))
635 0
    CHRS = c_hr
636 9
    if met < 0.85:
637 0
        CHCS = 3.0
638
    else:
639 0
        CHCS = 5.66 * (met - 0.85) ** 0.39
640 9
    if CHCS < 3.0:
641 0
        CHCS = 3.0
642 0
    CTCS = CHCS + CHRS
643 0
    RCLOS = 1.52 / ((met - wme / met_factor) + 0.6944) - 0.1835
644 0
    RCLS = 0.155 * RCLOS
645 0
    FACLS = 1.0 + k_clo * RCLOS
646 0
    FCLS = 1.0 / (1.0 + 0.155 * FACLS * CTCS * RCLOS)
647 0
    IMS = 0.45
648 0
    ICLS = IMS * CHCS / CTCS * (1 - FCLS) / (CHCS / CTCS - FCLS * IMS)
649 0
    RAS = 1.0 / (FACLS * CTCS)
650 0
    REAS = 1.0 / (lr * FACLS * CHCS)
651 0
    RECLS = RCLS / (lr * ICLS)
652 0
    HD_S = 1.0 / (RAS + RCLS)
653 0
    HE_S = 1.0 / (REAS + RECLS)
654

655 0
    delta = 0.0001
656 0
    dx = 100.0
657 0
    set_old = round(temp_skin - hsk / HD_S, 2)
658 9
    while abs(dx) > 0.01:
659 0
        err_1 = (
660
            hsk
661
            - HD_S * (temp_skin - set_old)
662
            - W
663
            * HE_S
664
            * (PSSK - 0.5 * (math.exp(18.6686 - 4030.183 / (set_old + 235.0))))
665
        )
666 0
        err_2 = (
667
            hsk
668
            - HD_S * (temp_skin - (set_old + delta))
669
            - W
670
            * HE_S
671
            * (PSSK - 0.5 * (math.exp(18.6686 - 4030.183 / (set_old + delta + 235.0))))
672
        )
673 0
        _set = set_old - delta * err_1 / (err_2 - err_1)
674 0
        dx = _set - set_old
675 0
        set_old = _set
676

677 0
    return _set
678

679

680 9
def adaptive_ashrae(tdb, tr, t_running_mean, v, units="SI"):
681
    """
682
    Determines the adaptive thermal comfort based on ASHRAE 55. The adaptive model
683
    relates indoor design temperatures or acceptable temperature ranges to outdoor
684
    meteorological
685
    or climatological parameters.
686

687
    Parameters
688
    ----------
689
    tdb : float
690
        dry bulb air temperature, default in [°C] in [°F] if `units` = 'IP'
691
    tr : float
692
        mean radiant temperature, default in [°C] in [°F] if `units` = 'IP'
693
    t_running_mean: float
694
        running mean temperature, default in [°C] in [°C] in [°F] if `units` = 'IP'
695
    v : float
696
        air velocity, default in [m/s] in [fps] if `units` = 'IP'
697
    units: str default="SI"
698
        select the SI (International System of Units) or the IP (Imperial Units) system.
699

700
    Returns
701
    -------
702
    tmp_cmf : float
703
        Comfort temperature a that specific running mean temperature, default in [°C]
704
        or in [°F]
705
    tmp_cmf_80_low : float
706
        Lower acceptable comfort temperature for 80% occupants, default in [°C] or in [°F]
707
    tmp_cmf_80_up : float
708
        Upper acceptable comfort temperature for 80% occupants, default in [°C] or in [°F]
709
    tmp_cmf_90_low : float
710
        Lower acceptable comfort temperature for 90% occupants, default in [°C] or in [°F]
711
    tmp_cmf_90_up : float
712
        Upper acceptable comfort temperature for 90% occupants, default in [°C] or in [°F]
713
    acceptability_80 : bol
714
        Acceptability for 80% occupants
715
    acceptability_90 : bol
716
        Acceptability for 90% occupants
717

718
    Notes
719
    -----
720
    You can use this function to calculate if your conditions are within the `adaptive
721
    thermal comfort region`.
722
    Calculations with comply with the ASHRAE 55 2017 Standard [1]_.
723

724
    Examples
725
    --------
726
    .. code-block:: python
727

728
        >>> from pythermalcomfort.models import adaptive_ashrae
729
        >>> results = adaptive_ashrae(tdb=25, tr=25, t_running_mean=20, v=0.1)
730
        >>> print(results)
731
        {'tmp_cmf': 24.0, 'tmp_cmf_80_low': 20.5, 'tmp_cmf_80_up': 27.5,
732
        'tmp_cmf_90_low': 21.5, 'tmp_cmf_90_up': 26.5, 'acceptability_80': True,
733
        'acceptability_90': False}
734

735
        >>> print(results['acceptability_80'])
736
        True
737
        # The conditions you entered are considered to be comfortable for by 80% of the
738
        occupants
739

740
        >>> # for users who wants to use the IP system
741
        >>> results = adaptive_ashrae(tdb=77, tr=77, t_running_mean=68, v=0.3, units='ip')
742
        >>> print(results)
743
        {'tmp_cmf': 75.2, 'tmp_cmf_80_low': 68.9, 'tmp_cmf_80_up': 81.5,
744
        'tmp_cmf_90_low': 70.7, 'tmp_cmf_90_up': 79.7, 'acceptability_80': True,
745
        'acceptability_90': False}
746

747
        >>> results = adaptive_ashrae(tdb=25, tr=25, t_running_mean=9, v=0.1)
748
        ValueError: The running mean is outside the standards applicability limits
749
        # The adaptive thermal comfort model can only be used
750
        # if the running mean temperature is higher than 10°C
751

752
    Raises
753
    ------
754
    ValueError
755
        Raised if the input are outside the Standard's applicability limits
756

757
    """
758 9
    if units.lower() == "ip":
759 9
        tdb, tr, t_running_mean, vr = units_converter(
760
            tdb=tdb, tr=tr, tmp_running_mean=t_running_mean, v=v
761
        )
762

763 9
    check_standard_compliance(standard="ashrae", tdb=tdb, tr=tr, v=v)
764

765 9
    to = t_o(tdb, tr, v)
766

767
    # See if the running mean temperature is between 10 °C and 33.5 °C (the range where
768
    # the adaptive model is supposed to be used)
769 9
    if 10.0 <= t_running_mean <= 33.5:
770

771 9
        cooling_effect = 0
772
        # calculate cooling effect of elevated air speed when top > 25 degC.
773 9
        if v >= 0.6 and to >= 25:
774 9
            if v < 0.9:
775 0
                cooling_effect = 1.2
776 9
            elif v < 1.2:
777 0
                cooling_effect = 1.8
778
            else:
779 0
                cooling_effect = 2.2
780

781
        # Figure out the relation between comfort and outdoor temperature depending on
782
        # the level of conditioning.
783 9
        t_cmf = 0.31 * t_running_mean + 17.8
784 9
        tmp_cmf_80_low = t_cmf - 3.5
785 9
        tmp_cmf_90_low = t_cmf - 2.5
786 9
        tmp_cmf_80_up = t_cmf + 3.5 + cooling_effect
787 9
        tmp_cmf_90_up = t_cmf + 2.5 + cooling_effect
788

789 9
        def acceptability(t_cmf_lower, t_cmf_upper):
790
            # See if the conditions are comfortable.
791 9
            if t_cmf_lower < to < t_cmf_upper:
792 9
                return True
793
            else:
794 9
                return False
795

796 9
        acceptability_80 = acceptability(tmp_cmf_80_low, tmp_cmf_80_up)
797 9
        acceptability_90 = acceptability(tmp_cmf_90_low, tmp_cmf_90_up)
798

799 9
        if units.lower() == "ip":
800 9
            (
801
                t_cmf,
802
                tmp_cmf_80_low,
803
                tmp_cmf_80_up,
804
                tmp_cmf_90_low,
805
                tmp_cmf_90_up,
806
            ) = units_converter(
807
                from_units="si",
808
                tmp_cmf=t_cmf,
809
                tmp_cmf_80_low=tmp_cmf_80_low,
810
                tmp_cmf_80_up=tmp_cmf_80_up,
811
                tmp_cmf_90_low=tmp_cmf_90_low,
812
                tmp_cmf_90_up=tmp_cmf_90_up,
813
            )
814

815 9
        results = {
816
            "tmp_cmf": t_cmf,
817
            "tmp_cmf_80_low": tmp_cmf_80_low,
818
            "tmp_cmf_80_up": tmp_cmf_80_up,
819
            "tmp_cmf_90_low": tmp_cmf_90_low,
820
            "tmp_cmf_90_up": tmp_cmf_90_up,
821
            "acceptability_80": acceptability_80,
822
            "acceptability_90": acceptability_90,
823
        }
824

825
    else:
826 9
        raise ValueError(
827
            "The running mean is outside the standards applicability limits"
828
        )
829

830 9
    return results
831

832

833 9
def adaptive_en(tdb, tr, t_running_mean, v, units="SI"):
834
    """ Determines the adaptive thermal comfort based on EN 16798-1 2019 [3]_
835

836
    Parameters
837
    ----------
838
    tdb : float
839
        dry bulb air temperature, default in [°C] in [°F] if `units` = 'IP'
840
    tr : float
841
        mean radiant temperature, default in [°C] in [°F] if `units` = 'IP'
842
    t_running_mean: float
843
        running mean temperature, default in [°C] in [°C] in [°F] if `units` = 'IP'
844
    v : float
845
        air velocity, default in [m/s] in [fps] if `units` = 'IP'
846

847
        Note: Indoor operative temperature correction is applicable for buildings equipped
848
        with fans or personal systems providing building occupants with personal
849
        control over air speed at occupant level.
850
        For operative temperatures above 25°C the comfort zone upper limit can be
851
        increased by 1.2 °C (0.6 < v < 0.9 m/s), 1.8 °C (0.9 < v < 1.2 m/s), 2.2 °C ( v
852
        > 1.2 m/s)
853
    units: str default="SI"
854
        select the SI (International System of Units) or the IP (Imperial Units) system.
855

856
    Returns
857
    -------
858
    tmp_cmf : float
859
        Comfort temperature at that specific running mean temperature, default in [°C]
860
        or in [°F]
861
    acceptability_cat_i : bol
862
        If the indoor conditions comply with comfort category I
863
    acceptability_cat_ii : bol
864
        If the indoor conditions comply with comfort category II
865
    acceptability_cat_iii : bol
866
        If the indoor conditions comply with comfort category III
867
    tmp_cmf_cat_i_up : float
868
        Upper acceptable comfort temperature for category I, default in [°C] or in [°F]
869
    tmp_cmf_cat_ii_up : float
870
        Upper acceptable comfort temperature for category II, default in [°C] or in [°F]
871
    tmp_cmf_cat_iii_up : float
872
        Upper acceptable comfort temperature for category III, default in [°C] or in [°F]
873
    tmp_cmf_cat_i_low : float
874
        Lower acceptable comfort temperature for category I, default in [°C] or in [°F]
875
    tmp_cmf_cat_ii_low : float
876
        Lower acceptable comfort temperature for category II, default in [°C] or in [°F]
877
    tmp_cmf_cat_iii_low : float
878
        Lower acceptable comfort temperature for category III, default in [°C] or in [°F]
879

880
    Notes
881
    -----
882
    You can use this function to calculate if your conditions are within the EN
883
    adaptive thermal comfort region.
884
    Calculations with comply with the EN 16798-1 2019 [1]_.
885

886
    Examples
887
    --------
888
    .. code-block:: python
889

890
        >>> from pythermalcomfort.models import adaptive_en
891
        >>> results = adaptive_en(tdb=25, tr=25, t_running_mean=20, v=0.1)
892
        >>> print(results)
893
        {'tmp_cmf': 25.4, 'acceptability_cat_i': True, 'acceptability_cat_ii': True,
894
        'acceptability_cat_iii': True, ... }
895

896
        >>> print(results['acceptability_cat_i'])
897
        True
898
        # The conditions you entered are considered to comply with Category I
899

900
        >>> # for users who wants to use the IP system
901
        >>> results = adaptive_en(tdb=77, tr=77, t_running_mean=68, v=0.3, units='ip')
902
        >>> print(results)
903
        {'tmp_cmf': 77.7, 'acceptability_cat_i': True, 'acceptability_cat_ii': True,
904
        'acceptability_cat_iii': True, ... }
905

906
        >>> results = adaptive_en(tdb=25, tr=25, t_running_mean=9, v=0.1)
907
        ValueError: The running mean is outside the standards applicability limits
908
        # The adaptive thermal comfort model can only be used
909
        # if the running mean temperature is between 10 °C and 30 °C
910

911
    Raises
912
    ------
913
    ValueError
914
        Raised if the input are outside the Standard's applicability limits
915

916
    """
917

918 9
    if units.lower() == "ip":
919 0
        tdb, tr, t_running_mean, vr = units_converter(
920
            tdb=tdb, tr=tr, tmp_running_mean=t_running_mean, v=v
921
        )
922

923 9
    if (t_running_mean < 10) or (t_running_mean > 30):
924 0
        raise ValueError(
925
            "The running mean is outside the standards applicability limits"
926
        )
927

928 0
    to = t_o(tdb, tr, v)
929

930 0
    cooling_effect = 0
931
    # calculate cooling effect of elevated air speed when top > 25 degC.
932 9
    if v >= 0.6 and to >= 25:
933 9
        if v < 0.9:
934 0
            cooling_effect = 1.2
935 9
        elif v < 1.2:
936 0
            cooling_effect = 1.8
937
        else:
938 0
            cooling_effect = 2.2
939

940 0
    t_cmf = 0.33 * t_running_mean + 18.8
941

942 0
    t_cmf_i_lower = t_cmf - 3
943 0
    t_cmf_ii_lower = t_cmf - 4
944 0
    t_cmf_iii_lower = t_cmf - 5
945 0
    t_cmf_i_upper = t_cmf + 2 + cooling_effect
946 0
    t_cmf_ii_upper = t_cmf + 3 + cooling_effect
947 0
    t_cmf_iii_upper = t_cmf + 4 + cooling_effect
948

949 0
    def between(val, low, high):
950 0
        return low < val < high
951

952 9
    if between(to, t_cmf_i_lower, t_cmf_i_upper):
953 0
        acceptability_i, acceptability_ii, acceptability_iii = True, True, True
954 9
    elif between(to, t_cmf_ii_lower, t_cmf_ii_upper):
955 0
        acceptability_ii, acceptability_iii = True, True
956 0
        acceptability_i = False
957 9
    elif between(to, t_cmf_iii_lower, t_cmf_iii_upper):
958 0
        acceptability_iii = True
959 0
        acceptability_i, acceptability_ii = False, False
960
    else:
961 0
        acceptability_i, acceptability_ii, acceptability_iii = False, False, False
962

963 9
    if units.lower() == "ip":
964 0
        t_cmf, t_cmf_i_upper, t_cmf_ii_upper, t_cmf_iii_upper = units_converter(
965
            from_units="si",
966
            tmp_cmf=t_cmf,
967
            tmp_cmf_cat_i_up=t_cmf_i_upper,
968
            tmp_cmf_cat_ii_up=t_cmf_ii_upper,
969
            tmp_cmf_cat_iii_up=t_cmf_iii_upper,
970
        )
971 0
        t_cmf_i_lower, t_cmf_ii_lower, t_cmf_iii_lower = units_converter(
972
            from_units="si",
973
            tmp_cmf_cat_i_low=t_cmf_i_lower,
974
            tmp_cmf_cat_ii_low=t_cmf_ii_lower,
975
            tmp_cmf_cat_iii_low=t_cmf_iii_lower,
976
        )
977

978 0
    results = {
979
        "tmp_cmf": round(t_cmf, 1),
980
        "acceptability_cat_i": acceptability_i,
981
        "acceptability_cat_ii": acceptability_ii,
982
        "acceptability_cat_iii": acceptability_iii,
983
        "tmp_cmf_cat_i_up": round(t_cmf_i_upper, 1),
984
        "tmp_cmf_cat_ii_up": round(t_cmf_ii_upper, 1),
985
        "tmp_cmf_cat_iii_up": round(t_cmf_iii_upper, 1),
986
        "tmp_cmf_cat_i_low": round(t_cmf_i_lower, 1),
987
        "tmp_cmf_cat_ii_low": round(t_cmf_ii_lower, 1),
988
        "tmp_cmf_cat_iii_low": round(t_cmf_iii_lower, 1),
989
    }
990

991 0
    return results
992

993

994 9
def utci(tdb, tr, v, rh, units="SI"):
995
    """ Determines the Universal Thermal Climate Index (UTCI). The UTCI is the
996
    equivalent temperature for the environment derived from a reference environment.
997
    It is defined as the air temperature of the reference environment which produces
998
    the same strain index value in comparison with the reference individual's response
999
    to the real
1000
    environment. It is regarded as one of the most comprehensive indices for
1001
    calculating heat stress in outdoor spaces. The parameters that are taken into
1002
    account for calculating
1003
    UTCI involve dry-bulb temperature, mean radiation temperature, the pressure of
1004
    water vapor or relative humidity, and wind speed (at the elevation of 10 m) [7]_.
1005

1006
    Parameters
1007
    ----------
1008
    tdb : float
1009
        dry bulb air temperature, default in [°C] in [°F] if `units` = 'IP'
1010
    tr : float
1011
        mean radiant temperature, default in [°C] in [°F] if `units` = 'IP'
1012
    v : float
1013
        relative air velocity, default in [m/s] in [fps] if `units` = 'IP'
1014
    rh : float
1015
        relative humidity, [%]
1016
    units: str default="SI"
1017
        select the SI (International System of Units) or the IP (Imperial Units) system.
1018

1019
    Returns
1020
    -------
1021
    utci : float
1022
         Universal Thermal Climate Index, [°C] or in [°F]
1023

1024
    Notes
1025
    -----
1026
    You can use this function to calculate the Universal Thermal Climate Index (`UTCI`)
1027
    The applicability wind speed value must be between 0.5 and 17 m/s.
1028

1029
    .. _UTCI: http://www.utci.org/utcineu/utcineu.php
1030

1031
    Examples
1032
    --------
1033
    .. code-block:: python
1034

1035
        >>> from pythermalcomfort.models import utci
1036
        >>> utci(tdb=25, tr=25, v=1.0, rh=50)
1037
        24.6
1038

1039
        >>> # for users who wants to use the IP system
1040
        >>> utci(tdb=77, tr=77, v=3.28, rh=50, units='ip')
1041
        76.4
1042

1043
    Raises
1044
    ------
1045
    ValueError
1046
        Raised if the input are outside the Standard's applicability limits
1047

1048
    """
1049

1050 9
    if units.lower() == "ip":
1051 9
        tdb, tr, v = units_converter(tdb=tdb, tr=tr, v=v)
1052

1053 9
    check_standard_compliance(standard="utci", tdb=tdb, tr=tr, v=v)
1054

1055 9
    def exponential(t_db):
1056 9
        g = [
1057
            -2836.5744,
1058
            -6028.076559,
1059
            19.54263612,
1060
            -0.02737830188,
1061
            0.000016261698,
1062
            (7.0229056 * (10 ** (-10))),
1063
            (-1.8680009 * (10 ** (-13))),
1064
        ]
1065 9
        tk = t_db + 273.15  # air temp in K
1066 9
        es = 2.7150305 * math.log1p(tk)
1067 9
        for count, i in enumerate(g):
1068 9
            es = es + (i * (tk ** (count - 2)))
1069 9
        es = math.exp(es) * 0.01  # convert Pa to hPa
1070 9
        return es
1071

1072
    # Do a series of checks to be sure that the input values are within the bounds
1073
    # accepted by the model.
1074 9
    if (
1075
        (tdb < -50.0)
1076
        or (tdb > 50.0)
1077
        or (tr - tdb < -30.0)
1078
        or (tr - tdb > 70.0)
1079
        or (v < 0.5)
1080
        or (v > 17)
1081
    ):
1082 0
        raise ValueError(
1083
            "The value you entered are outside the equation applicability limits"
1084
        )
1085

1086
    # This is a python version of the UTCI_approx function
1087
    # Version a 0.002, October 2009
1088
    # tdb: air temperature, degrees Celsius
1089
    # ehPa: water vapour presure, hPa=hecto Pascal
1090
    # Tmrt: mean radiant temperature, degrees Celsius
1091
    # va10m: wind speed 10m above ground level in m/s
1092

1093 9
    eh_pa = exponential(tdb) * (rh / 100.0)
1094 9
    delta_t_tr = tr - tdb
1095 9
    pa = eh_pa / 10.0  # convert vapour pressure to kPa
1096

1097 9
    utci_approx = (
1098
        tdb
1099
        + (0.607562052)
1100
        + (-0.0227712343) * tdb
1101
        + (8.06470249 * (10 ** (-4))) * tdb * tdb
1102
        + (-1.54271372 * (10 ** (-4))) * tdb * tdb * tdb
1103
        + (-3.24651735 * (10 ** (-6))) * tdb * tdb * tdb * tdb
1104
        + (7.32602852 * (10 ** (-8))) * tdb * tdb * tdb * tdb * tdb
1105
        + (1.35959073 * (10 ** (-9))) * tdb * tdb * tdb * tdb * tdb * tdb
1106
        + (-2.25836520) * v
1107
        + (0.0880326035) * tdb * v
1108
        + (0.00216844454) * tdb * tdb * v
1109
        + (-1.53347087 * (10 ** (-5))) * tdb * tdb * tdb * v
1110
        + (-5.72983704 * (10 ** (-7))) * tdb * tdb * tdb * tdb * v
1111
        + (-2.55090145 * (10 ** (-9))) * tdb * tdb * tdb * tdb * tdb * v
1112
        + (-0.751269505) * v * v
1113
        + (-0.00408350271) * tdb * v * v
1114
        + (-5.21670675 * (10 ** (-5))) * tdb * tdb * v * v
1115
        + (1.94544667 * (10 ** (-6))) * tdb * tdb * tdb * v * v
1116
        + (1.14099531 * (10 ** (-8))) * tdb * tdb * tdb * tdb * v * v
1117
        + (0.158137256) * v * v * v
1118
        + (-6.57263143 * (10 ** (-5))) * tdb * v * v * v
1119
        + (2.22697524 * (10 ** (-7))) * tdb * tdb * v * v * v
1120
        + (-4.16117031 * (10 ** (-8))) * tdb * tdb * tdb * v * v * v
1121
        + (-0.0127762753) * v * v * v * v
1122
        + (9.66891875 * (10 ** (-6))) * tdb * v * v * v * v
1123
        + (2.52785852 * (10 ** (-9))) * tdb * tdb * v * v * v * v
1124
        + (4.56306672 * (10 ** (-4))) * v * v * v * v * v
1125
        + (-1.74202546 * (10 ** (-7))) * tdb * v * v * v * v * v
1126
        + (-5.91491269 * (10 ** (-6))) * v * v * v * v * v * v
1127
        + (0.398374029) * delta_t_tr
1128
        + (1.83945314 * (10 ** (-4))) * tdb * delta_t_tr
1129
        + (-1.73754510 * (10 ** (-4))) * tdb * tdb * delta_t_tr
1130
        + (-7.60781159 * (10 ** (-7))) * tdb * tdb * tdb * delta_t_tr
1131
        + (3.77830287 * (10 ** (-8))) * tdb * tdb * tdb * tdb * delta_t_tr
1132
        + (5.43079673 * (10 ** (-10))) * tdb * tdb * tdb * tdb * tdb * delta_t_tr
1133
        + (-0.0200518269) * v * delta_t_tr
1134
        + (8.92859837 * (10 ** (-4))) * tdb * v * delta_t_tr
1135
        + (3.45433048 * (10 ** (-6))) * tdb * tdb * v * delta_t_tr
1136
        + (-3.77925774 * (10 ** (-7))) * tdb * tdb * tdb * v * delta_t_tr
1137
        + (-1.69699377 * (10 ** (-9))) * tdb * tdb * tdb * tdb * v * delta_t_tr
1138
        + (1.69992415 * (10 ** (-4))) * v * v * delta_t_tr
1139
        + (-4.99204314 * (10 ** (-5))) * tdb * v * v * delta_t_tr
1140
        + (2.47417178 * (10 ** (-7))) * tdb * tdb * v * v * delta_t_tr
1141
        + (1.07596466 * (10 ** (-8))) * tdb * tdb * tdb * v * v * delta_t_tr
1142
        + (8.49242932 * (10 ** (-5))) * v * v * v * delta_t_tr
1143
        + (1.35191328 * (10 ** (-6))) * tdb * v * v * v * delta_t_tr
1144
        + (-6.21531254 * (10 ** (-9))) * tdb * tdb * v * v * v * delta_t_tr
1145
        + (-4.99410301 * (10 ** (-6))) * v * v * v * v * delta_t_tr
1146
        + (-1.89489258 * (10 ** (-8))) * tdb * v * v * v * v * delta_t_tr
1147
        + (8.15300114 * (10 ** (-8))) * v * v * v * v * v * delta_t_tr
1148
        + (7.55043090 * (10 ** (-4))) * delta_t_tr * delta_t_tr
1149
        + (-5.65095215 * (10 ** (-5))) * tdb * delta_t_tr * delta_t_tr
1150
        + (-4.52166564 * (10 ** (-7))) * tdb * tdb * delta_t_tr * delta_t_tr
1151
        + (2.46688878 * (10 ** (-8))) * tdb * tdb * tdb * delta_t_tr * delta_t_tr
1152
        + (2.42674348 * (10 ** (-10))) * tdb * tdb * tdb * tdb * delta_t_tr * delta_t_tr
1153
        + (1.54547250 * (10 ** (-4))) * v * delta_t_tr * delta_t_tr
1154
        + (5.24110970 * (10 ** (-6))) * tdb * v * delta_t_tr * delta_t_tr
1155
        + (-8.75874982 * (10 ** (-8))) * tdb * tdb * v * delta_t_tr * delta_t_tr
1156
        + (-1.50743064 * (10 ** (-9))) * tdb * tdb * tdb * v * delta_t_tr * delta_t_tr
1157
        + (-1.56236307 * (10 ** (-5))) * v * v * delta_t_tr * delta_t_tr
1158
        + (-1.33895614 * (10 ** (-7))) * tdb * v * v * delta_t_tr * delta_t_tr
1159
        + (2.49709824 * (10 ** (-9))) * tdb * tdb * v * v * delta_t_tr * delta_t_tr
1160
        + (6.51711721 * (10 ** (-7))) * v * v * v * delta_t_tr * delta_t_tr
1161
        + (1.94960053 * (10 ** (-9))) * tdb * v * v * v * delta_t_tr * delta_t_tr
1162
        + (-1.00361113 * (10 ** (-8))) * v * v * v * v * delta_t_tr * delta_t_tr
1163
        + (-1.21206673 * (10 ** (-5))) * delta_t_tr * delta_t_tr * delta_t_tr
1164
        + (-2.18203660 * (10 ** (-7))) * tdb * delta_t_tr * delta_t_tr * delta_t_tr
1165
        + (7.51269482 * (10 ** (-9))) * tdb * tdb * delta_t_tr * delta_t_tr * delta_t_tr
1166
        + (9.79063848 * (10 ** (-11)))
1167
        * tdb
1168
        * tdb
1169
        * tdb
1170
        * delta_t_tr
1171
        * delta_t_tr
1172
        * delta_t_tr
1173
        + (1.25006734 * (10 ** (-6))) * v * delta_t_tr * delta_t_tr * delta_t_tr
1174
        + (-1.81584736 * (10 ** (-9))) * tdb * v * delta_t_tr * delta_t_tr * delta_t_tr
1175
        + (-3.52197671 * (10 ** (-10)))
1176
        * tdb
1177
        * tdb
1178
        * v
1179
        * delta_t_tr
1180
        * delta_t_tr
1181
        * delta_t_tr
1182
        + (-3.36514630 * (10 ** (-8))) * v * v * delta_t_tr * delta_t_tr * delta_t_tr
1183
        + (1.35908359 * (10 ** (-10)))
1184
        * tdb
1185
        * v
1186
        * v
1187
        * delta_t_tr
1188
        * delta_t_tr
1189
        * delta_t_tr
1190
        + (4.17032620 * (10 ** (-10)))
1191
        * v
1192
        * v
1193
        * v
1194
        * delta_t_tr
1195
        * delta_t_tr
1196
        * delta_t_tr
1197
        + (-1.30369025 * (10 ** (-9)))
1198
        * delta_t_tr
1199
        * delta_t_tr
1200
        * delta_t_tr
1201
        * delta_t_tr
1202
        + (4.13908461 * (10 ** (-10)))
1203
        * tdb
1204
        * delta_t_tr
1205
        * delta_t_tr
1206
        * delta_t_tr
1207
        * delta_t_tr
1208
        + (9.22652254 * (10 ** (-12)))
1209
        * tdb
1210
        * tdb
1211
        * delta_t_tr
1212
        * delta_t_tr
1213
        * delta_t_tr
1214
        * delta_t_tr
1215
        + (-5.08220384 * (10 ** (-9)))
1216
        * v
1217
        * delta_t_tr
1218
        * delta_t_tr
1219
        * delta_t_tr
1220
        * delta_t_tr
1221
        + (-2.24730961 * (10 ** (-11)))
1222
        * tdb
1223
        * v
1224
        * delta_t_tr
1225
        * delta_t_tr
1226
        * delta_t_tr
1227
        * delta_t_tr
1228
        + (1.17139133 * (10 ** (-10)))
1229
        * v
1230
        * v
1231
        * delta_t_tr
1232
        * delta_t_tr
1233
        * delta_t_tr
1234
        * delta_t_tr
1235
        + (6.62154879 * (10 ** (-10)))
1236
        * delta_t_tr
1237
        * delta_t_tr
1238
        * delta_t_tr
1239
        * delta_t_tr
1240
        * delta_t_tr
1241
        + (4.03863260 * (10 ** (-13)))
1242
        * tdb
1243
        * delta_t_tr
1244
        * delta_t_tr
1245
        * delta_t_tr
1246
        * delta_t_tr
1247
        * delta_t_tr
1248
        + (1.95087203 * (10 ** (-12)))
1249
        * v
1250
        * delta_t_tr
1251
        * delta_t_tr
1252
        * delta_t_tr
1253
        * delta_t_tr
1254
        * delta_t_tr
1255
        + (-4.73602469 * (10 ** (-12)))
1256
        * delta_t_tr
1257
        * delta_t_tr
1258
        * delta_t_tr
1259
        * delta_t_tr
1260
        * delta_t_tr
1261
        * delta_t_tr
1262
        + (5.12733497) * pa
1263
        + (-0.312788561) * tdb * pa
1264
        + (-0.0196701861) * tdb * tdb * pa
1265
        + (9.99690870 * (10 ** (-4))) * tdb * tdb * tdb * pa
1266
        + (9.51738512 * (10 ** (-6))) * tdb * tdb * tdb * tdb * pa
1267
        + (-4.66426341 * (10 ** (-7))) * tdb * tdb * tdb * tdb * tdb * pa
1268
        + (0.548050612) * v * pa
1269
        + (-0.00330552823) * tdb * v * pa
1270
        + (-0.00164119440) * tdb * tdb * v * pa
1271
        + (-5.16670694 * (10 ** (-6))) * tdb * tdb * tdb * v * pa
1272
        + (9.52692432 * (10 ** (-7))) * tdb * tdb * tdb * tdb * v * pa
1273
        + (-0.0429223622) * v * v * pa
1274
        + (0.00500845667) * tdb * v * v * pa
1275
        + (1.00601257 * (10 ** (-6))) * tdb * tdb * v * v * pa
1276
        + (-1.81748644 * (10 ** (-6))) * tdb * tdb * tdb * v * v * pa
1277
        + (-1.25813502 * (10 ** (-3))) * v * v * v * pa
1278
        + (-1.79330391 * (10 ** (-4))) * tdb * v * v * v * pa
1279
        + (2.34994441 * (10 ** (-6))) * tdb * tdb * v * v * v * pa
1280
        + (1.29735808 * (10 ** (-4))) * v * v * v * v * pa
1281
        + (1.29064870 * (10 ** (-6))) * tdb * v * v * v * v * pa
1282
        + (-2.28558686 * (10 ** (-6))) * v * v * v * v * v * pa
1283
        + (-0.0369476348) * delta_t_tr * pa
1284
        + (0.00162325322) * tdb * delta_t_tr * pa
1285
        + (-3.14279680 * (10 ** (-5))) * tdb * tdb * delta_t_tr * pa
1286
        + (2.59835559 * (10 ** (-6))) * tdb * tdb * tdb * delta_t_tr * pa
1287
        + (-4.77136523 * (10 ** (-8))) * tdb * tdb * tdb * tdb * delta_t_tr * pa
1288
        + (8.64203390 * (10 ** (-3))) * v * delta_t_tr * pa
1289
        + (-6.87405181 * (10 ** (-4))) * tdb * v * delta_t_tr * pa
1290
        + (-9.13863872 * (10 ** (-6))) * tdb * tdb * v * delta_t_tr * pa
1291
        + (5.15916806 * (10 ** (-7))) * tdb * tdb * tdb * v * delta_t_tr * pa
1292
        + (-3.59217476 * (10 ** (-5))) * v * v * delta_t_tr * pa
1293
        + (3.28696511 * (10 ** (-5))) * tdb * v * v * delta_t_tr * pa
1294
        + (-7.10542454 * (10 ** (-7))) * tdb * tdb * v * v * delta_t_tr * pa
1295
        + (-1.24382300 * (10 ** (-5))) * v * v * v * delta_t_tr * pa
1296
        + (-7.38584400 * (10 ** (-9))) * tdb * v * v * v * delta_t_tr * pa
1297
        + (2.20609296 * (10 ** (-7))) * v * v * v * v * delta_t_tr * pa
1298
        + (-7.32469180 * (10 ** (-4))) * delta_t_tr * delta_t_tr * pa
1299
        + (-1.87381964 * (10 ** (-5))) * tdb * delta_t_tr * delta_t_tr * pa
1300
        + (4.80925239 * (10 ** (-6))) * tdb * tdb * delta_t_tr * delta_t_tr * pa
1301
        + (-8.75492040 * (10 ** (-8))) * tdb * tdb * tdb * delta_t_tr * delta_t_tr * pa
1302
        + (2.77862930 * (10 ** (-5))) * v * delta_t_tr * delta_t_tr * pa
1303
        + (-5.06004592 * (10 ** (-6))) * tdb * v * delta_t_tr * delta_t_tr * pa
1304
        + (1.14325367 * (10 ** (-7))) * tdb * tdb * v * delta_t_tr * delta_t_tr * pa
1305
        + (2.53016723 * (10 ** (-6))) * v * v * delta_t_tr * delta_t_tr * pa
1306
        + (-1.72857035 * (10 ** (-8))) * tdb * v * v * delta_t_tr * delta_t_tr * pa
1307
        + (-3.95079398 * (10 ** (-8))) * v * v * v * delta_t_tr * delta_t_tr * pa
1308
        + (-3.59413173 * (10 ** (-7))) * delta_t_tr * delta_t_tr * delta_t_tr * pa
1309
        + (7.04388046 * (10 ** (-7))) * tdb * delta_t_tr * delta_t_tr * delta_t_tr * pa
1310
        + (-1.89309167 * (10 ** (-8)))
1311
        * tdb
1312
        * tdb
1313
        * delta_t_tr
1314
        * delta_t_tr
1315
        * delta_t_tr
1316
        * pa
1317
        + (-4.79768731 * (10 ** (-7))) * v * delta_t_tr * delta_t_tr * delta_t_tr * pa
1318
        + (7.96079978 * (10 ** (-9)))
1319
        * tdb
1320
        * v
1321
        * delta_t_tr
1322
        * delta_t_tr
1323
        * delta_t_tr
1324
        * pa
1325
        + (1.62897058 * (10 ** (-9)))
1326
        * v
1327
        * v
1328
        * delta_t_tr
1329
        * delta_t_tr
1330
        * delta_t_tr
1331
        * pa
1332
        + (3.94367674 * (10 ** (-8)))
1333
        * delta_t_tr
1334
        * delta_t_tr
1335
        * delta_t_tr
1336
        * delta_t_tr
1337
        * pa
1338
        + (-1.18566247 * (10 ** (-9)))
1339
        * tdb
1340
        * delta_t_tr
1341
        * delta_t_tr
1342
        * delta_t_tr
1343
        * delta_t_tr
1344
        * pa
1345
        + (3.34678041 * (10 ** (-10)))
1346
        * v
1347
        * delta_t_tr
1348
        * delta_t_tr
1349
        * delta_t_tr
1350
        * delta_t_tr
1351
        * pa
1352
        + (-1.15606447 * (10 ** (-10)))
1353
        * delta_t_tr
1354
        * delta_t_tr
1355
        * delta_t_tr
1356
        * delta_t_tr
1357
        * delta_t_tr
1358
        * pa
1359
        + (-2.80626406) * pa * pa
1360
        + (0.548712484) * tdb * pa * pa
1361
        + (-0.00399428410) * tdb * tdb * pa * pa
1362
        + (-9.54009191 * (10 ** (-4))) * tdb * tdb * tdb * pa * pa
1363
        + (1.93090978 * (10 ** (-5))) * tdb * tdb * tdb * tdb * pa * pa
1364
        + (-0.308806365) * v * pa * pa
1365
        + (0.0116952364) * tdb * v * pa * pa
1366
        + (4.95271903 * (10 ** (-4))) * tdb * tdb * v * pa * pa
1367
        + (-1.90710882 * (10 ** (-5))) * tdb * tdb * tdb * v * pa * pa
1368
        + (0.00210787756) * v * v * pa * pa
1369
        + (-6.98445738 * (10 ** (-4))) * tdb * v * v * pa * pa
1370
        + (2.30109073 * (10 ** (-5))) * tdb * tdb * v * v * pa * pa
1371
        + (4.17856590 * (10 ** (-4))) * v * v * v * pa * pa
1372
        + (-1.27043871 * (10 ** (-5))) * tdb * v * v * v * pa * pa
1373
        + (-3.04620472 * (10 ** (-6))) * v * v * v * v * pa * pa
1374
        + (0.0514507424) * delta_t_tr * pa * pa
1375
        + (-0.00432510997) * tdb * delta_t_tr * pa * pa
1376
        + (8.99281156 * (10 ** (-5))) * tdb * tdb * delta_t_tr * pa * pa
1377
        + (-7.14663943 * (10 ** (-7))) * tdb * tdb * tdb * delta_t_tr * pa * pa
1378
        + (-2.66016305 * (10 ** (-4))) * v * delta_t_tr * pa * pa
1379
        + (2.63789586 * (10 ** (-4))) * tdb * v * delta_t_tr * pa * pa
1380
        + (-7.01199003 * (10 ** (-6))) * tdb * tdb * v * delta_t_tr * pa * pa
1381
        + (-1.06823306 * (10 ** (-4))) * v * v * delta_t_tr * pa * pa
1382
        + (3.61341136 * (10 ** (-6))) * tdb * v * v * delta_t_tr * pa * pa
1383
        + (2.29748967 * (10 ** (-7))) * v * v * v * delta_t_tr * pa * pa
1384
        + (3.04788893 * (10 ** (-4))) * delta_t_tr * delta_t_tr * pa * pa
1385
        + (-6.42070836 * (10 ** (-5))) * tdb * delta_t_tr * delta_t_tr * pa * pa
1386
        + (1.16257971 * (10 ** (-6))) * tdb * tdb * delta_t_tr * delta_t_tr * pa * pa
1387
        + (7.68023384 * (10 ** (-6))) * v * delta_t_tr * delta_t_tr * pa * pa
1388
        + (-5.47446896 * (10 ** (-7))) * tdb * v * delta_t_tr * delta_t_tr * pa * pa
1389
        + (-3.59937910 * (10 ** (-8))) * v * v * delta_t_tr * delta_t_tr * pa * pa
1390
        + (-4.36497725 * (10 ** (-6))) * delta_t_tr * delta_t_tr * delta_t_tr * pa * pa
1391
        + (1.68737969 * (10 ** (-7)))
1392
        * tdb
1393
        * delta_t_tr
1394
        * delta_t_tr
1395
        * delta_t_tr
1396
        * pa
1397
        * pa
1398
        + (2.67489271 * (10 ** (-8)))
1399
        * v
1400
        * delta_t_tr
1401
        * delta_t_tr
1402
        * delta_t_tr
1403
        * pa
1404
        * pa
1405
        + (3.23926897 * (10 ** (-9)))
1406
        * delta_t_tr
1407
        * delta_t_tr
1408
        * delta_t_tr
1409
        * delta_t_tr
1410
        * pa
1411
        * pa
1412
        + (-0.0353874123) * pa * pa * pa
1413
        + (-0.221201190) * tdb * pa * pa * pa
1414
        + (0.0155126038) * tdb * tdb * pa * pa * pa
1415
        + (-2.63917279 * (10 ** (-4))) * tdb * tdb * tdb * pa * pa * pa
1416
        + (0.0453433455) * v * pa * pa * pa
1417
        + (-0.00432943862) * tdb * v * pa * pa * pa
1418
        + (1.45389826 * (10 ** (-4))) * tdb * tdb * v * pa * pa * pa
1419
        + (2.17508610 * (10 ** (-4))) * v * v * pa * pa * pa
1420
        + (-6.66724702 * (10 ** (-5))) * tdb * v * v * pa * pa * pa
1421
        + (3.33217140 * (10 ** (-5))) * v * v * v * pa * pa * pa
1422
        + (-0.00226921615) * delta_t_tr * pa * pa * pa
1423
        + (3.80261982 * (10 ** (-4))) * tdb * delta_t_tr * pa * pa * pa
1424
        + (-5.45314314 * (10 ** (-9))) * tdb * tdb * delta_t_tr * pa * pa * pa
1425
        + (-7.96355448 * (10 ** (-4))) * v * delta_t_tr * pa * pa * pa
1426
        + (2.53458034 * (10 ** (-5))) * tdb * v * delta_t_tr * pa * pa * pa
1427
        + (-6.31223658 * (10 ** (-6))) * v * v * delta_t_tr * pa * pa * pa
1428
        + (3.02122035 * (10 ** (-4))) * delta_t_tr * delta_t_tr * pa * pa * pa
1429
        + (-4.77403547 * (10 ** (-6))) * tdb * delta_t_tr * delta_t_tr * pa * pa * pa
1430
        + (1.73825715 * (10 ** (-6))) * v * delta_t_tr * delta_t_tr * pa * pa * pa
1431
        + (-4.09087898 * (10 ** (-7)))
1432
        * delta_t_tr
1433
        * delta_t_tr
1434
        * delta_t_tr
1435
        * pa
1436
        * pa
1437
        * pa
1438
        + (0.614155345) * pa * pa * pa * pa
1439
        + (-0.0616755931) * tdb * pa * pa * pa * pa
1440
        + (0.00133374846) * tdb * tdb * pa * pa * pa * pa
1441
        + (0.00355375387) * v * pa * pa * pa * pa
1442
        + (-5.13027851 * (10 ** (-4))) * tdb * v * pa * pa * pa * pa
1443
        + (1.02449757 * (10 ** (-4))) * v * v * pa * pa * pa * pa
1444
        + (-0.00148526421) * delta_t_tr * pa * pa * pa * pa
1445
        + (-4.11469183 * (10 ** (-5))) * tdb * delta_t_tr * pa * pa * pa * pa
1446
        + (-6.80434415 * (10 ** (-6))) * v * delta_t_tr * pa * pa * pa * pa
1447
        + (-9.77675906 * (10 ** (-6))) * delta_t_tr * delta_t_tr * pa * pa * pa * pa
1448
        + (0.0882773108) * pa * pa * pa * pa * pa
1449
        + (-0.00301859306) * tdb * pa * pa * pa * pa * pa
1450
        + (0.00104452989) * v * pa * pa * pa * pa * pa
1451
        + (2.47090539 * (10 ** (-4))) * delta_t_tr * pa * pa * pa * pa * pa
1452
        + (0.00148348065) * pa * pa * pa * pa * pa * pa
1453
    )
1454

1455 9
    if units.lower() == "ip":
1456 9
        utci_approx = units_converter(tmp=utci_approx, from_units="si")[0]
1457

1458 9
    return round(utci_approx, 1)
1459

1460

1461 9
def clo_tout(tout, units="SI"):
1462
    """ Representative clothing insulation Icl as a function of outdoor air temperature
1463
    at 06:00 a.m [4]_.
1464

1465
    Parameters
1466
    ----------
1467
    tout : float
1468
        outdoor air temperature at 06:00 a.m., default in [°C] in [°F] if `units` = 'IP'
1469
    units: str default="SI"
1470
        select the SI (International System of Units) or the IP (Imperial Units) system.
1471

1472
    Returns
1473
    -------
1474
    clo : float
1475
         Representative clothing insulation Icl, [clo]
1476

1477
    Notes
1478
    -----
1479
    The ASHRAE 55 2017 states that it is acceptable to determine the clothing
1480
    insulation Icl using this equation in mechanically conditioned buildings [1]_.
1481

1482
    Examples
1483
    --------
1484
    .. code-block:: python
1485

1486
        >>> from pythermalcomfort.models import clo_tout
1487
        >>> clo_tout(tout=27)
1488
        0.46
1489

1490
    """
1491 9
    if units.lower() == "ip":
1492 9
        tout = units_converter(tmp=tout)[0]
1493

1494 9
    if tout < -5:
1495 0
        clo = 1
1496 9
    elif tout < 5:
1497 0
        clo = 0.818 - 0.0364 * tout
1498 9
    elif tout < 26:
1499 0
        clo = 10 ** (-0.1635 - 0.0066 * tout)
1500
    else:
1501 9
        clo = 0.46
1502

1503 9
    return round(clo, 2)
1504

1505

1506 9
def vertical_tmp_grad_ppd(tdb, tr, vr, rh, met, clo, vertical_tmp_grad, units="SI"):
1507
    """ Calculates the percentage of thermally dissatisfied people with a vertical
1508
    temperature gradient between feet and head [1]_.
1509
    This equation is only applicable for vr < 0.2 m/s (40 fps).
1510

1511
    Parameters
1512
    ----------
1513
    tdb : float
1514
        dry bulb air temperature, default in [°C] in [°F] if `units` = 'IP'
1515

1516
        Note: The air temperature is the average value over two heights: 0.6 m (24 in.)
1517
        and 1.1 m (43 in.) for seated occupants
1518
        and 1.1 m (43 in.) and 1.7 m (67 in.) for standing occupants.
1519
    tr : float
1520
        mean radiant temperature, default in [°C] in [°F] if `units` = 'IP'
1521
    vr : float
1522
        relative air velocity, default in [m/s] in [fps] if `units` = 'IP'
1523

1524
        Note: vr is the relative air velocity caused by body movement and not the air
1525
        speed measured by the air velocity sensor.
1526
        It can be calculate using the function
1527
        :py:meth:`pythermalcomfort.psychrometrics.v_relative`.
1528
    rh : float
1529
        relative humidity, [%]
1530
    met : float
1531
        metabolic rate, [met]
1532
    clo : float
1533
        clothing insulation, [clo]
1534
    vertical_tmp_grad : float
1535
        vertical temperature gradient between the feet and the head, default in [°C/m]
1536
        in [°F/ft] if `units` = 'IP'
1537
    units: str default="SI"
1538
        select the SI (International System of Units) or the IP (Imperial Units) system.
1539

1540
    Returns
1541
    -------
1542
    PPD_vg: float
1543
        Predicted Percentage of Dissatisfied occupants with vertical temperature
1544
        gradient, [%]
1545
    Acceptability: bol
1546
        The ASHRAE 55 2017 standard defines that the value of air speed at the ankle
1547
        level is acceptable if PPD_ad is lower or equal than 5 %
1548

1549
    Examples
1550
    --------
1551
    .. code-block:: python
1552

1553
        >>> from pythermalcomfort.models import vertical_tmp_grad_ppd
1554
        >>> results = vertical_tmp_grad_ppd(25, 25, 0.1, 50, 1.2, 0.5, 7)
1555
        >>> print(results)
1556
        {'PPD_vg': 12.6, 'Acceptability': False}
1557

1558
    """
1559 9
    if units.lower() == "ip":
1560 9
        tdb, tr, vr = units_converter(tdb=tdb, tr=tr, v=vr)
1561 9
        vertical_tmp_grad = vertical_tmp_grad / 1.8 * 3.28
1562

1563 9
    check_standard_compliance(
1564
        standard="ashrae", tdb=tdb, tr=tr, v_limited=vr, rh=rh, met=met, clo=clo
1565
    )
1566

1567 9
    tsv = pmv(tdb, tr, vr, rh, met, clo, standard="ashrae")
1568 9
    numerator = math.exp(0.13 * (tsv - 1.91) ** 2 + 0.15 * vertical_tmp_grad - 1.6)
1569 9
    ppd_val = round((numerator / (1 + numerator) - 0.345) * 100, 1)
1570 9
    acceptability = ppd_val <= 5
1571 9
    return {"PPD_vg": ppd_val, "Acceptability": acceptability}
1572

1573

1574 9
def ankle_draft(tdb, tr, vr, rh, met, clo, v_ankle, units="SI"):
1575
    """
1576
    Calculates the percentage of thermally dissatisfied people with the ankle draft (
1577
    0.1 m) above floor level [1]_.
1578
    This equation is only applicable for vr < 0.2 m/s (40 fps).
1579

1580
    Parameters
1581
    ----------
1582
    tdb : float
1583
        dry bulb air temperature, default in [°C] in [°F] if `units` = 'IP'
1584

1585
        Note: The air temperature is the average value over two heights: 0.6 m (24 in.)
1586
        and 1.1 m (43 in.) for seated occupants
1587
        and 1.1 m (43 in.) and 1.7 m (67 in.) for standing occupants.
1588
    tr : float
1589
        mean radiant temperature, default in [°C] in [°F] if `units` = 'IP'
1590
    vr : float
1591
        relative air velocity, default in [m/s] in [fps] if `units` = 'IP'
1592

1593
        Note: vr is the relative air velocity caused by body movement and not the air
1594
        speed measured by the air velocity sensor.
1595
        It can be calculate using the function
1596
        :py:meth:`pythermalcomfort.psychrometrics.v_relative`.
1597
    rh : float
1598
        relative humidity, [%]
1599
    met : float
1600
        metabolic rate, [met]
1601
    clo : float
1602
        clothing insulation, [clo]
1603
    v_ankle : float
1604
        air speed at the 0.1 m (4 in.) above the floor, default in [m/s] in [fps] if
1605
        `units` = 'IP'
1606
    units: str default="SI"
1607
        select the SI (International System of Units) or the IP (Imperial Units) system.
1608

1609
    Returns
1610
    -------
1611
    PPD_ad: float
1612
        Predicted Percentage of Dissatisfied occupants with ankle draft, [%]
1613
    Acceptability: bol
1614
        The ASHRAE 55 2017 standard defines that the value of air speed at the ankle
1615
        level is acceptable if PPD_ad is lower or equal than 20 %
1616

1617
    Examples
1618
    --------
1619
    .. code-block:: python
1620

1621
        >>> from pythermalcomfort.models import ankle_draft
1622
        >>> results = ankle_draft(25, 25, 0.2, 50, 1.2, 0.5, 0.3, units="SI")
1623
        >>> print(results)
1624
        {'PPD_ad': 18.6, 'Acceptability': True}
1625

1626
    """
1627 9
    if units.lower() == "ip":
1628 9
        tdb, tr, vr, v_ankle = units_converter(tdb=tdb, tr=tr, v=vr, vel=v_ankle)
1629

1630 9
    check_standard_compliance(
1631
        standard="ashrae", tdb=tdb, tr=tr, v_limited=vr, rh=rh, met=met, clo=clo
1632
    )
1633

1634 9
    tsv = pmv(tdb, tr, vr, rh, met, clo, standard="ashrae")
1635 9
    ppd_val = round(
1636
        math.exp(-2.58 + 3.05 * v_ankle - 1.06 * tsv)
1637
        / (1 + math.exp(-2.58 + 3.05 * v_ankle - 1.06 * tsv))
1638
        * 100,
1639
        1,
1640
    )
1641 9
    acceptability = ppd_val <= 20
1642 9
    return {"PPD_ad": ppd_val, "Acceptability": acceptability}
1643

1644

1645 9
def solar_gain(
1646
    sol_altitude,
1647
    sharp,
1648
    sol_radiation_dir,
1649
    sol_transmittance,
1650
    f_svv,
1651
    f_bes,
1652
    asw=0.7,
1653
    posture="seated",
1654
    floor_reflectance=0.6,
1655
):
1656
    """
1657
        Calculates the solar gain to the human body using the Effective Radiant Field (
1658
        ERF) [1]_. The ERF is a measure of the net energy flux to or from the human body.
1659
        ERF is expressed in W over human body surface area [w/m2]. In addition,
1660
        it calculates the delta mean radiant temperature. Which is the amount by which
1661
        the mean radiant
1662
        temperature of the space should be increased if no solar radiation is present.
1663

1664
        Parameters
1665
        ----------
1666
        sol_altitude : float
1667
            Solar altitude, degrees from horizontal [deg]. Ranges between 0 and 90.
1668
        sharp : float
1669
            Solar horizontal angle relative to the front of the person (SHARP) [deg].
1670
            Ranges between 0 and 180 and is symmetrical on either side. Zero (0) degrees
1671
            represents direct-beam radiation from the front, 90 degrees represents
1672
            direct-beam radiation from the side, and 180 degrees rep- resent direct-beam
1673
            radiation from the back. SHARP is the angle between the sun and the person
1674
            only. Orientation relative to compass or to room is not included in SHARP.
1675
        posture : str
1676
            Default 'seated' list of available options 'standing', 'supine' or 'seated'
1677
        sol_radiation_dir : float
1678
            Direct-beam solar radiation, [W/m2]. Ranges between 200 and 1000. See Table
1679
            C2-3 of ASHRAE 55 2017 [1]_.
1680
        sol_transmittance : float
1681
            Total solar transmittance, ranges from 0 to 1. The total solar
1682
            transmittance of window systems, including glazing unit, blinds, and other
1683
            façade treatments, shall be determined using one of the following methods:
1684
            i) Provided by manufacturer or from the National Fenestration Rating
1685
            Council approved Lawrence Berkeley National Lab International Glazing
1686
            Database.
1687
            ii) Glazing unit plus venetian blinds or other complex or unique shades
1688
            shall be calculated using National Fenestration Rating Council approved
1689
            software or Lawrence Berkeley National Lab Complex Glazing Database.
1690
        f_svv : float
1691
            Fraction of sky-vault view fraction exposed to body, ranges from 0 to 1.
1692
            It can be calculate using the function
1693
            :py:meth:`pythermalcomfort.psychrometrics.f_svv`.
1694
        f_bes : float
1695
            Fraction of the possible body surface exposed to sun, ranges from 0 to 1.
1696
            See Table C2-2 and equation C-7 ASHRAE 55 2017 [1]_.
1697
        asw: float
1698
            The average short-wave absorptivity of the occupant. It will range widely,
1699
            depending on the color of the occupant’s skin as well as the color and
1700
            amount of clothing covering the body.
1701
            A value of 0.7 shall be used unless more specific information about the
1702
            clothing or skin color of the occupants is available.
1703
            Note: Short-wave absorptivity typically ranges from 0.57 to 0.84, depending
1704
            on skin and clothing color. More information is available in Blum (1945).
1705
        floor_reflectance: float
1706
            Floor refectance. It is assumed to be constant and equal to 0.6.
1707

1708
        Notes
1709
        -----
1710
        More information on the calculation procedure can be found in Appendix C of [1]_.
1711

1712
        Returns
1713
        -------
1714
        erf: float
1715
            Solar gain to the human body using the Effective Radiant Field [W/m2]
1716
        delta_mrt: float
1717
            Delta mean radiant temperature. The amount by which the mean radiant
1718
            temperature of the space should be increased if no solar radiation is present.
1719

1720
        Examples
1721
        --------
1722
        .. code-block:: python
1723

1724
            >>> from pythermalcomfort.models import solar_gain
1725
            >>> results = solar_gain(sol_altitude=0, sharp=120,
1726
            sol_radiation_dir=800, sol_transmittance=0.5, f_svv=0.5, f_bes=0.5,
1727
            asw=0.7, posture='seated')
1728
            >>> print(results)
1729
            {'erf': 42.9, 'delta_mrt': 10.3}
1730

1731
        """
1732

1733 9
    posture = posture.lower()
1734 9
    if posture not in ["standing", "supine", "seated"]:
1735 0
        raise ValueError("Posture has to be either standing, supine or seated")
1736

1737 9
    def find_span(arr, x):
1738 9
        for i in range(0, len(arr)):
1739 9
            if arr[i + 1] >= x >= arr[i]:
1740 9
                return i
1741 0
        return -1
1742

1743 9
    deg_to_rad = 0.0174532925
1744 9
    hr = 6
1745 9
    i_diff = 0.2 * sol_radiation_dir
1746

1747 9
    fp_table = [
1748
        [0.35, 0.35, 0.314, 0.258, 0.206, 0.144, 0.082],
1749
        [0.342, 0.342, 0.31, 0.252, 0.2, 0.14, 0.082],
1750
        [0.33, 0.33, 0.3, 0.244, 0.19, 0.132, 0.082],
1751
        [0.31, 0.31, 0.275, 0.228, 0.175, 0.124, 0.082],
1752
        [0.283, 0.283, 0.251, 0.208, 0.16, 0.114, 0.082],
1753
        [0.252, 0.252, 0.228, 0.188, 0.15, 0.108, 0.082],
1754
        [0.23, 0.23, 0.214, 0.18, 0.148, 0.108, 0.082],
1755
        [0.242, 0.242, 0.222, 0.18, 0.153, 0.112, 0.082],
1756
        [0.274, 0.274, 0.245, 0.203, 0.165, 0.116, 0.082],
1757
        [0.304, 0.304, 0.27, 0.22, 0.174, 0.121, 0.082],
1758
        [0.328, 0.328, 0.29, 0.234, 0.183, 0.125, 0.082],
1759
        [0.344, 0.344, 0.304, 0.244, 0.19, 0.128, 0.082],
1760
        [0.347, 0.347, 0.308, 0.246, 0.191, 0.128, 0.082],
1761
    ]
1762 9
    if posture == "seated":
1763 9
        fp_table = [
1764
            [0.29, 0.324, 0.305, 0.303, 0.262, 0.224, 0.177],
1765
            [0.292, 0.328, 0.294, 0.288, 0.268, 0.227, 0.177],
1766
            [0.288, 0.332, 0.298, 0.29, 0.264, 0.222, 0.177],
1767
            [0.274, 0.326, 0.294, 0.289, 0.252, 0.214, 0.177],
1768
            [0.254, 0.308, 0.28, 0.276, 0.241, 0.202, 0.177],
1769
            [0.23, 0.282, 0.262, 0.26, 0.233, 0.193, 0.177],
1770
            [0.216, 0.26, 0.248, 0.244, 0.22, 0.186, 0.177],
1771
            [0.234, 0.258, 0.236, 0.227, 0.208, 0.18, 0.177],
1772
            [0.262, 0.26, 0.224, 0.208, 0.196, 0.176, 0.177],
1773
            [0.28, 0.26, 0.21, 0.192, 0.184, 0.17, 0.177],
1774
            [0.298, 0.256, 0.194, 0.174, 0.168, 0.168, 0.177],
1775
            [0.306, 0.25, 0.18, 0.156, 0.156, 0.166, 0.177],
1776
            [0.3, 0.24, 0.168, 0.152, 0.152, 0.164, 0.177],
1777
        ]
1778

1779 9
    if posture == "supine":
1780 0
        alt_temp = sol_altitude
1781 0
        sol_altitude = abs(90 - sharp)
1782 0
        sharp = alt_temp
1783

1784 9
    alt_range = [0, 15, 30, 45, 60, 75, 90]
1785 9
    az_range = [0, 15, 30, 45, 60, 75, 90, 105, 120, 135, 150, 165, 180]
1786 9
    alt_i = find_span(alt_range, sol_altitude)
1787 9
    az_i = find_span(az_range, sharp)
1788 9
    fp11 = fp_table[az_i][alt_i]
1789 9
    fp12 = fp_table[az_i][alt_i + 1]
1790 9
    fp21 = fp_table[az_i + 1][alt_i]
1791 9
    fp22 = fp_table[az_i + 1][alt_i + 1]
1792 9
    az1 = az_range[az_i]
1793 9
    az2 = az_range[az_i + 1]
1794 9
    alt1 = alt_range[alt_i]
1795 9
    alt2 = alt_range[alt_i + 1]
1796 9
    fp = fp11 * (az2 - sharp) * (alt2 - sol_altitude)
1797 9
    fp += fp21 * (sharp - az1) * (alt2 - sol_altitude)
1798 9
    fp += fp12 * (az2 - sharp) * (sol_altitude - alt1)
1799 9
    fp += fp22 * (sharp - az1) * (sol_altitude - alt1)
1800 9
    fp /= (az2 - az1) * (alt2 - alt1)
1801

1802 9
    f_eff = 0.725
1803 9
    if posture == "seated":
1804 9
        f_eff = 0.696
1805

1806 9
    sw_abs = asw
1807 9
    lw_abs = 0.95
1808

1809 9
    e_diff = f_eff * f_svv * 0.5 * sol_transmittance * i_diff
1810 9
    e_direct = f_eff * fp * sol_transmittance * f_bes * sol_radiation_dir
1811 9
    e_refl = (
1812
        f_eff
1813
        * f_svv
1814
        * 0.5
1815
        * sol_transmittance
1816
        * (sol_radiation_dir * math.sin(sol_altitude * deg_to_rad) + i_diff)
1817
        * floor_reflectance
1818
    )
1819

1820 9
    e_solar = e_diff + e_direct + e_refl
1821 9
    erf = e_solar * (sw_abs / lw_abs)
1822 9
    d_mrt = erf / (hr * f_eff)
1823

1824
    # print(fp, e_diff, e_direct, e_refl, e_solar, erf, d_mrt)
1825

1826 9
    return {"erf": round(erf, 1), "delta_mrt": round(d_mrt, 1)}
1827

1828

1829
# add the following models:
1830
# todo radiant_tmp_asymmetry
1831
# todo draft
1832
# todo floor_surface_tmp
1833
# todo effective_tmp
1834
# more info here: https://www.rdocumentation.org/packages/comf/versions/0.1.9
1835
# more info here: https://rdrr.io/cran/comf/man/

Read our documentation on viewing source code .

Loading