PET2017_refactored.py

"""PET calculation routine (2017)
Refactored to follow PEP 8, add type hints, and encapsulate logic.

Original algorithm: Peter Hoeppe (VDI 3787‑2) + Djordje Spasic (Ladybug).
Refactor: OpenAI ChatGPT (o3), 2025‑07‑29.
"""
from __future__ import annotations

import math
from dataclasses import dataclass
from typing import List, Tuple

import numpy as np
from scipy.optimize import fsolve

# ---------------------------------------------------------------------------
# Module‑level constants (SI units unless noted)
# ---------------------------------------------------------------------------
PO_HPA: float = 1013.25        # Standard air pressure [hPa]
BLOOD_SPECIFIC_HEAT: float = 3.64e3   # Blood specific heat [J kg⁻¹ K⁻¹]
AIR_SPECIFIC_HEAT: float = 1.01e3     # Air specific heat  [J kg⁻¹ K⁻¹]
SKIN_EMISSIVITY: float = 0.99
CLOTHING_EMISSIVITY: float = 0.95
LATENT_HEAT_EVAP: float = 2.42e6      # Latent heat of evaporation [J kg⁻¹]
STEFAN_BOLTZMANN: float = 5.67e-8     # Stefan–Boltzmann constant [W m⁻² K⁻⁴]
MECH_EFFICIENCY: float = 0.0          # Fraction of metabolism converted to work

# Set‑point temperatures (°C)
TC_SET: float = 36.6
TSK_SET: float = 34.0
TBODY_SET: float = 0.1 * TSK_SET + 0.9 * TC_SET

# ---------------------------------------------------------------------------
# Data classes
# ---------------------------------------------------------------------------
@dataclass(slots=True)
class Person:
    """Anthropometric and demographic data."""
    age: float      # years
    sex: int        # 1 = male, 2 = female
    height: float   # m
    mass: float     # kg
    posture: int = 1  # 1 = standing, 2 = seated, 3 = lying


@dataclass(slots=True)
class Environment:
    """Micro‑climatic variables."""
    air_temp: float        # °C
    mrt: float             # mean radiant temperature, °C
    rel_humidity: float    # %
    wind_speed: float      # m s⁻¹
    pressure: float = PO_HPA  # hPa


@dataclass(slots=True)
class ActivityClothing:
    """Metabolic activity and apparel insulation."""
    metabolic_rate: float  # W
    icl: float             # clothing insulation [clo]


# ---------------------------------------------------------------------------
# Helper (physiological) functions
# ---------------------------------------------------------------------------

def _vasoconstriction(t_core: float, t_skin: float) -> Tuple[float, float]:
    """Return (skin_blood_flow [L m⁻² h⁻¹], alpha weighting for T_body)."""
    sig_skin = max(TSK_SET - t_skin, 0.0)
    sig_core = max(t_core - TC_SET, 0.0)
    flow = (6.3 + 75.0 * sig_core) / (1.0 + 0.5 * sig_skin)
    flow = min(flow, 90.0)
    alpha = 0.1  # fixed for steady‑state version
    return flow, alpha


def _sweating(t_body: float, t_skin: float) -> float:
    """Sweat rate [g m⁻² h⁻¹] using a linearised Gagge relation."""
    sig_body = max(t_body - TBODY_SET, 0.0)
    q_sw = 0.30494 * sig_body  # g m⁻² h⁻¹
    return min(q_sw, 500.0)


# ---------------------------------------------------------------------------
# Core calculator class
# ---------------------------------------------------------------------------
class PETCalculator:
    """Three‑node thermo‑physiological model with PET evaluation."""

    def __init__(self, person: Person, env: Environment, act_clo: ActivityClothing):
        self.person = person
        self.env = env
        self.act_clo = act_clo

    # --------------------------------------------------------------------
    # Public API
    # --------------------------------------------------------------------

    def solve_node_temperatures(
        self,
        initial_guess: List[float] | Tuple[float, float, float] = (38.0, 40.0, 40.0),
    ) -> np.ndarray:
        """Return equilibrium temperatures [T_core, T_skin, T_clo] (°C)."""
        sol, *_ = fsolve(self._energy_balance_vector, x0=initial_guess, args=(True,), full_output=True)
        return sol

    def compute_pet(
        self,
        stable_nodes: np.ndarray,
        t_min: float = -40.0,
        t_max: float = 60.0,
        tol: float = 1e-6,
    ) -> float:
        """Compute Physiological Equivalent Temperature (PET)."""

        def _scalar_balance(t_air: float) -> float:
            std_env = Environment(
                air_temp=t_air,
                mrt=t_air,
                rel_humidity=50.0,
                wind_speed=0.1,
                pressure=self.env.pressure,
            )
            std_clo = ActivityClothing(metabolic_rate=self.act_clo.metabolic_rate, icl=0.9)
            return PETCalculator(self.person, std_env, std_clo)._energy_balance_scalar(stable_nodes, t_air)

        lo, hi = t_min, t_max
        pet = (lo + hi) / 2.0
        while hi - lo > tol:
            if _scalar_balance(lo) * _scalar_balance(pet) < 0:
                hi = pet
            else:
                lo = pet
            pet = (lo + hi) / 2.0
        return pet

    def system_energy_balance_vector(self, nodes: np.ndarray) -> np.ndarray:
        """Return [core, skin, clothing] energy balances (W m⁻²)."""
        return np.asarray(self._energy_balance_vector(nodes, True), dtype=float)

    # --------------------------------------------------------------------
    # Internal machinery
    # --------------------------------------------------------------------

    def _energy_balance_vector(self, nodes: List[float], mode_actual: bool) -> List[float]:
        return self._energy_balance_generic(nodes, mode_actual, vector_out=True)

    def _energy_balance_scalar(self, nodes: np.ndarray, air_temp_ref: float) -> float:
        return self._energy_balance_generic(nodes, mode_actual=False, vector_out=False, air_override=air_temp_ref)

    # pylint: disable=too‑many‑locals,too‑many‑statements
    def _energy_balance_generic(
        self,
        nodes: List[float] | np.ndarray,
        mode_actual: bool,
        *,
        vector_out: bool,
        air_override: float | None = None,
    ) -> List[float] | float:
        """Vector or scalar energy balance based on the VDI three‑node model."""
        t_core, t_skin, t_clo = nodes
        env, person, act = self.env, self.person, self.act_clo

        # Allow overriding air temperature / MRT (for PET dichotomy)
        ta = env.air_temp if air_override is None else air_override
        tmrt = env.mrt if air_override is None else air_override
        hr, v, p = env.rel_humidity, env.wind_speed, env.pressure

        # Body surface area & posture coefficient
        adu = 0.203 * person.mass ** 0.425 * person.height ** 0.725
        feff = 0.725 if person.posture in (1, 3) else 0.696

        # Clothing area factors
        f_cl = 1.0 + 0.31 * act.icl
        facl = (173.51 * act.icl - 2.36 - 100.76 * act.icl**2 + 19.28 * act.icl**3) / 100.0
        a_clo = adu * facl + adu * (f_cl - 1.0)
        a_eff_r = adu * feff

        # Vapour pressure of water in air [hPa]
        vpa = hr / 100.0 * 6.105 * math.exp(17.27 * ta / (237.7 + ta)) if mode_actual else 12.0

        # Convective heat‑transfer coefficient
        hc_map = {
            1: 2.67 + 6.5 * v ** 0.67,  # standing
            2: 2.26 + 7.42 * v ** 0.67, # sitting
            3: 8.6 * v ** 0.513,        # lying
        }
        hc = hc_map[person.posture]
        if person.posture == 3:
            hc *= (p / PO_HPA) ** 0.55  # pressure correction when lying

        # Basal metabolism (WHO/ICBN)
        metab_f = 3.19 * person.mass ** 0.75 * (
            1 + 0.004 * (30 - person.age) + 0.018 * (person.height * 100 / person.mass ** (1/3) - 42.1)
        )
        metab_m = 3.45 * person.mass ** 0.75 * (
            1 + 0.004 * (30 - person.age) + 0.01  * (person.height * 100 / person.mass ** (1/3) - 43.4)
        )

        if mode_actual:
            metab = (act.metabolic_rate + metab_m) / adu
            fec   = (act.metabolic_rate + metab_f) / adu
        else:
            metab = (80 + metab_m) / adu
            fec   = (80 + metab_f) / adu

        h_met = (metab if person.sex == 1 else fec) * (1.0 - MECH_EFFICIENCY)

        # Respiratory heat losses
        t_exp = 0.47 * ta + 21.0
        d_vent = metab * 1.44e-6
        e_res_sens = AIR_SPECIFIC_HEAT * (ta - t_exp) * d_vent
        vp_exp = 6.11 * 10 ** (7.45 * t_exp / (235.0 + t_exp))
        e_res_lat = 0.623 * LATENT_HEAT_EVAP / p * (vpa - vp_exp) * d_vent
        e_res = e_res_sens + e_res_lat

        # Clothing & geometry
        r_cl = act.icl / 6.45  # m² K W⁻¹
        if facl > 1.0:
            facl = 1.0
        if act.icl >= 2.0:
            y = 1.0
        elif 0.6 < act.icl < 2.0:
            y = (person.height - 0.2) / person.height
        elif 0.3 < act.icl <= 0.6:
            y = 0.5
        elif 0.0 < act.icl <= 0.3:
            y = 0.1
        else:
            y = 0.0

        r2 = adu * (f_cl - 1.0 + facl) / (6.28 * person.height * y)
        r1 = facl * adu / (6.28 * person.height * y)
        di = r2 - r1

        q_flow, alpha = _vasoconstriction(t_core, t_skin)
        t_body = alpha * t_skin + (1.0 - alpha) * t_core

        h_tcl = 6.28 * person.height * y * di / (r_cl * math.log(r2 / r1) * a_clo)

        q_sw = _sweating(t_body, t_skin)
        e_sw = LATENT_HEAT_EVAP / 1000.0 * q_sw / 3600.0

        pv_sk = 6.105 * math.exp((17.27 * (t_skin + 273.15) - 4717.03) / (237.7 + t_skin))
        l_w = 1.67
        h_e_diff = hc * l_w
        f_ecl = 1.0 / (1.0 + 0.92 * hc * r_cl)
        e_max = h_e_diff * f_ecl * (pv_sk - vpa)
        w = min(e_sw / e_max if e_max else 0.0, 1.0)
        if e_sw < 0:
            e_sw = 0.0
        i_m = 0.38
        r_ecl = (1.0 / (f_cl * hc) + r_cl) / (l_w * i_m)
        e_diff = (1.0 - w) * (pv_sk - vpa) / r_ecl
        evap = -(e_diff + e_sw)

        # Radiation
        r_bare = a_eff_r * (1.0 - facl) * SKIN_EMISSIVITY * STEFAN_BOLTZMANN * (
            (tmrt + 273.15) ** 4 - (t_skin + 273.15) ** 4
        ) / adu
        r_clo = feff * a_clo * CLOTHING_EMISSIVITY * STEFAN_BOLTZMANN * (
            (tmrt + 273.15) ** 4 - (t_clo + 273.15) ** 4
        ) / adu
        r_sum = r_bare + r_clo

        # Convection
        c_bare = hc * (ta - t_skin) * adu * (1.0 - facl) / adu
        c_clo = hc * (ta - t_clo) * a_clo / adu
        c_sum = c_bare + c_clo

        # Energy balance equations
        core_bal = h_met + e_res - (q_flow / 3600.0 * BLOOD_SPECIFIC_HEAT + 5.28) * (t_core - t_skin)
        skin_bal = (
            r_bare + c_bare + evap + (q_flow / 3600.0 * BLOOD_SPECIFIC_HEAT + 5.28) * (t_core - t_skin) - h_tcl * (t_skin - t_clo)
        )
        clo_bal = c_clo + r_clo + h_tcl * (t_skin - t_clo)
        scalar_bal = h_met + e_res + r_sum + c_sum + evap

        if vector_out:
            return [core_bal, skin_bal, clo_bal]
        return scalar_bal


# ---------------------------------------------------------------------------
# Demonstration / self‑test
# ---------------------------------------------------------------------------
if __name__ == "__main__":
    person = Person(age=35, sex=1, height=1.80, mass=75.0, posture=1)
    env = Environment(air_temp=30.0, mrt=20.0, rel_humidity=50.0, wind_speed=1.0)
    act = ActivityClothing(metabolic_rate=80.0, icl=0.5)

    calc = PETCalculator(person, env, act)
    stable = calc.solve_node_temperatures()
    print("Nodes temperature [T_core, T_skin, T_clo]", stable)
    print("Thermal Balance", calc.system_energy_balance_vector(stable)[0])
    print("PET:", round(calc.compute_pet(stable), 2))