# seizoen.py - de warmtepomp van het station over een rekenjaar, uur per uur
import numpy as np
from kringproces import cyclus                    # het model van hierboven (Peng-Robinson, propaan)

UA, WINST = 50.0, 533.3                          # warmteverlies (W/K) en interne winst (W) van het station
VD = 0.00162                                     # slagvolume maal toerental van de compressor bij 90 Hz (m³/s)


def stooklijn(Tb):                               # aanvoertemperatuur van het water (°C)
    return min(40.0, max(25.0, 20 + 15 * (20 - Tb) / 28))


def ontdooien(Tb):                               # extra energie om ijs van de verdamper te smelten (fractie)
    if Tb <= -7 or Tb >= 7:
        return 0.0
    return 0.10 * (Tb + 7) / 9 if Tb <= 2 else 0.10 * (7 - Tb) / 5


def warmtepomp(Tb):
    c = cyclus(Tb - 6.0, stooklijn(Tb) + 5.0)     # verdampen 6 K onder de lucht, condenseren 5 K boven het water
    lam = 0.95 - 0.025 * c["pc"] / c["pe"]        # volumetrisch rendement van de compressor
    return VD * c["rho1"] * lam * c["qc"], c["cop"] * 0.92 * (1 - ontdooien(Tb))   # Qmax (W), COP


T = np.loadtxt(open("rekenjaar.csv"))            # 8760 uurtemperaturen van het rekenjaar (°C)
rooster = np.arange(-20.0, 21.0)                 # tabel om de 1 °C, daarna interpoleren
tab = np.array([warmtepomp(t) for t in rooster])
Qmax, COP = np.interp(T, rooster, tab[:, 0]), np.interp(T, rooster, tab[:, 1])

Q = np.clip(UA * (20 - T) - WINST, 0, None)      # warmtevraag per uur (W)
P = Q / COP + np.where(Q > 0, 20 + 30 * np.minimum(Q / Qmax, 1), 0)   # compressor + pomp + ventilator (W)
tekort = np.clip(Q - Qmax, 0, None)              # wat de noodverwarming zou moeten bijleggen

print(f"warmte {Q.sum() / 1e3:.0f} kWh, elektriciteit {P.sum() / 1e3:.0f} kWh, SCOP {Q.sum() / P.sum():.2f}")
print(f"waarvan hulpenergie (pomp en ventilator): {(P - Q / COP).sum() / 1e3:.0f} kWh")
print(f"hoogste warmtevraag {Q.max():.0f} W, hoogste opgenomen vermogen {P.max():.0f} W")
print(f"uren met warmtevraag {int((Q > 0).sum())}, uren onder de kleinste stand (pendelen): {int(((Q > 0) & (Q < Qmax / 3)).sum())}")
print(f"tekort voor de noodverwarming: {tekort.sum() / 1e3:.1f} kWh")