# lucht.py - CO2, O2 en waterdamp in de woonmodule op een winterdag, met vraaggestuurde ventilatie
import math

R, T = 8.314, 293.15                 # J/(mol K), binnen 20 °C
V = 6.0 * 4.0 * 2.6                  # woonmodule (m³)
n_lucht = 101325 * V / (R * T)       # mol lucht in de module
ps = lambda t: 611.2 * math.exp(17.62 * t / (243.12 + t))   # verzadigingsdruk water (Pa)


def activiteit(uur):                 # twee bewoners: slapen van 23 tot 7 uur, overdag actief
    return 0.72 if (uur >= 23 or uur < 7) else 1.14            # gemiddeld precies 1


co2, o2, pw = 420e-6, 0.2095, 0.40 * ps(20.0)                # molfracties en dampdruk (Pa)
co2_b, o2_b, pw_b = 420e-6, 0.2095, 0.90 * ps(-3.0)          # buitenlucht: -3 °C, 90 %
Q, dt = 100.0, 1 / 60                                        # m³/h en een stap van 1 minuut (h)
log = []
for stap in range(48 * 60):                                  # twee dagen
    uur = (stap * dt) % 24
    a = activiteit(uur)
    # ventilatie in drie trappen met hysterese, gestuurd door de CO2-sensor
    if co2 > 900e-6:
        Q = 150.0
    elif co2 < 750e-6 and Q == 150.0:
        Q = 100.0
    if (uur >= 23 or uur < 7) and co2 < 850e-6:
        Q = 70.0                                             # 's nachts stiller, zolang het kan
    elif Q == 70.0 and not (uur >= 23 or uur < 7):
        Q = 100.0
    ververs = Q / V                                          # per uur ververste fractie
    co2 += dt * (2 * a * 1.0 / 44.009e-3 / 24 / n_lucht - ververs * (co2 - co2_b))
    o2 += dt * (-2 * a * 0.84 / 31.998e-3 / 24 / n_lucht - ververs * (o2 - o2_b))
    pw += dt * (2 * a * 2000 / 18.015 / 24 * R * T / V - ververs * (pw - pw_b))
    log.append((stap * dt, co2 * 1e6, o2 * 100, 100 * pw / ps(20.0), Q))

dag2 = log[24 * 60:]
print(f"CO2 {min(r[1] for r in dag2):.0f} tot {max(r[1] for r in dag2):.0f} ppm, "
      f"O2 hoogstens {20.95 - min(r[2] for r in dag2):.3f} procentpunt onder buiten")
print(f"RV {min(r[3] for r in dag2):.0f} tot {max(r[3] for r in dag2):.0f} %, "
      f"gemiddeld debiet {sum(r[4] for r in dag2) / len(dag2):.0f} m³/h")
print(f"minuten boven 1000 ppm: {sum(1 for r in dag2 if r[1] > 1000)}")