# "humid-air.py" Compute air properties.
# Copyright (C) 2016, 2023, 2026 Aubrey Jaffer

# Permission to copy this code, to modify it, to redistribute it,
# to distribute modified versions, and to use it for any purpose is
# granted, subject to the following restrictions and understandings.

# 1.  Any copy made of this dictionary must include this copyright notice
# in full.

# 2.  I have made no warranty or representation that this data is
# error-free, and I am under no obligation to provide any services, by
# way of maintenance, update, or otherwise.

# 3.  In conjunction with products arising from the use of this material,
# there shall be no use of my name in any advertising, promotional, or
# sales literature without prior written consent in each case.

# Appendix B: Thermal and Transport Properties of Humid Air
# Aubrey Jaffer.
# Skin-friction and forced convection from rough and smooth plates.
# Thermo, 3(4):711–775, 2023, doi:10.3390/thermo3040040.

import math

def humid_air_properties(T_C, T_f_C, P_Pa, RH_percent):
    """
    Python translation of the Scheme function `humid-air-properties`.

    Inputs:
        T_C      : plate temperature in °C
        T_f_C    : air temperature in °C (for viscosity, etc.)
        P_Pa     : total pressure in Pa
        RH_percent : relative humidity in %

    Returns:
        [rho, mu, k, c_p]
        rho : humid air density (kg/m³)
        mu  : dynamic viscosity of humid air (Pa·s)
        k   : thermal conductivity of humid air (W/(m·K))
        c_p : specific heat at constant pressure (J/(kg·K))
    """

    # Constants
    R   = 8.3144621          # J/(mol·K), universal gas constant
    M_A = 28.97e-3           # kg/mol, molar mass of dry air
    M_V = 18.01528e-3        # kg/mol, molar mass of water vapor

    r_m = M_A / M_V

    T_K   = T_C   + 273.15
    T_f_K = T_f_C + 273.15

    # RH*P_V: vapor pressure at T_f_K
    RH_P_V = (
        0.01 * RH_percent *
        math.exp(
            (-0.63536311e4 / T_f_K)
            + 0.3404926034e2
            + T_f_K * (
                -0.19509874e-1
                + T_f_K * (0.12811805e-4)
            )
        )
    )

    # Density of humid air at temperature T_K
    def rho(T_K_local):
        # (/ (+ (* (- P.Pa RH*P_V) M_A)
        #       (* RH*P_V M_V))
        #    R T.K)
        return ((P_Pa - RH_P_V) * M_A + RH_P_V * M_V) / (R * T_K_local)

    # absolute humidity * r_m --> vapor pressure per dry air pressure
    chi_m = RH_P_V / (P_Pa - RH_P_V)

    # Dynamic viscosity of dry air (mu_a6) as polynomial in T_K
    mu_a6 = (
        0.40401e-6
        + T_K * (
            0.074582e-6
            + T_K * (
                -5.7171e-11
                + T_K * (
                    2.9928e-14
                    + T_K * (-6.2524e-18)
                )
            )
        )
    )

    # Dynamic viscosity of water vapor (mu_v) using gamma = 647.27 / T_K
    gamma = 647.27 / T_K
    mu_v = (
        1e-6
        / math.sqrt(gamma)
        / (
            0.0181583
            + gamma * (
                0.0177624
                + gamma * (
                    0.0105287
                    + gamma * (-0.0036744)
                )
            )
        )
    )

    # Phi(M_r, mu_r) = sqrt(.125)*sq(1+sqrt(mu_r)/(M_r)**.25)/sqrt(1+M_r)
    def Phi(M_r, mu_r):
        return (
            (1 + math.sqrt(mu_r) / (M_r ** 0.25)) ** 2
            / math.sqrt(1 + M_r)
            / math.sqrt(8.0)
        )

    # mum(x,mua,muv) =
    #   mua/(1+Phi(r_m,mua/muv)*x*r_m) + muv/(1+Phi(1/r_m,muv/mua)/(x*r_m))
    def mum(chi_m_local, mua, muv):
        term1 = mua / (1 + Phi(r_m, mua / muv) * chi_m_local)
        term2 = muv / (1 + (Phi(1.0 / r_m, muv / mua) / chi_m_local))
        return term1 + term2

    # mu_m(P,T,RH) = mum(x(P,T,RH), mu_a6(T), mu_v(T))
    mu = mum(chi_m, mu_a6, mu_v)

    # Reduced density and temperature for thermal conductivity of dry air
    rho_r = P_Pa / T_K / 287.058 / 314.3
    T_r   = T_K / 132.5

    # Thermal conductivity of dry air, k_a
    k_a = (
        0.0259778
        * (
            0.239503 * T_r
            + 0.00649768 * math.sqrt(T_r)
            + 1.0
            - 1.92615 / T_r
            + 2.00383 / (T_r ** 2)
            - 1.07553 / (T_r ** 3)
            + 0.229414 / (T_r ** 4)
            + 0.402287 * rho_r
            + 0.356603 * (rho_r ** 2)
            - 0.163159 * (rho_r ** 3)
            + 0.138059 * (rho_r ** 4)
            - 0.0201725 * (rho_r ** 5)
        )
    )

    # Thermal conductivity of water vapor, k_v, as polynomial in T_C
    k_v = (
        1.74822e-2
        + 7.69127e-5 * T_C
        - 3.23464e-7 * (T_C ** 2)
        + 2.59524e-9 * (T_C ** 3)
        - 3.17650e-12 * (T_C ** 4)
    )

    # Ratio of viscosities
    r_u = mu_a6 / mu_v

    # Mixture thermal conductivity k
    k = (
        k_a / (1 + Phi(r_m, r_u) * chi_m)
        + k_v / (1 + (Phi(1.0 / r_m, 1.0 / r_u) / chi_m))
    )

    # Specific heats of vapor and air
    c_pV = 1869.0 + T_C * (-0.2578 + T_C * 1.941e-2)
    c_pA = (
        1034.0
        + T_K * (
            -0.2849
            + T_K * (
                0.7817e-3
                + T_K * (
                    -0.4971e-6
                    + T_K * 0.1077e-9
                )
            )
        )
    )

    # Vapor pressure fraction
    chi_p = RH_P_V / P_Pa

    # Mixture specific heat c_p
    numerator_cp = (
        c_pA * (1 - chi_p) * M_A
        + c_pV * chi_p * M_V
    )
    denominator_cp = (
        (1 - chi_p) * M_A
        + chi_p * M_V
    )
    c_p = numerator_cp / denominator_cp

    return [rho(T_K), mu, k, c_p]

# print(humid_air_properties(75, 25, 101325, 50))
## [1.0080661310059178, 2.0463460975363706e-05, 0.029599589056937527, 1019.4455914897032]
# print(humid_air_properties(75, 25, 101325, 25))
## [1.0110641138683436, 2.053730764585827e-05, 0.0296570855331149, 1014.7920286897102]