Model/domain/sap10_calculator/worksheet/mean_internal_temperature.py

"""SAP 10.2 mean internal temperature (Tables 9, 9b, 9c).

The dwelling has two temperature zones during heating periods:

- Living area at T_h1 = 21 °C
- Elsewhere at T_h2 (Table 9 formulas, depending on control type)

When the heating is off, both zones cool toward a 'steady-cool' temperature
T_sc; the time-weighted reduction over the off-period is u (Table 9b).
Mean temperature = T_h − (u1 + u2) summed over the day's off-periods, then
the two zones are blended by living-area fraction (Table 9c step 7).

Standard SAP heating schedule:
- 9 hours heating per day, two off-periods of 7 h and 8 h (control types 1, 2)
- Control type 3 zones the heating: 9 h and 8 h off in 'elsewhere'

Reference: SAP 10.2 specification (14-03-2025) Table 9 (page 183),
Table 9b (page 184), Table 9c (page 185).
"""

from __future__ import annotations

from dataclasses import dataclass
from typing import Final, Optional

from domain.sap10_calculator.worksheet.utilisation_factor import utilisation_factor


_T_H1_C: Final[float] = 21.0
_HLP_CLAMP_FOR_T_H2: Final[float] = 6.0
_MONTHS: Final[range] = range(12)
_TIME_CONSTANT_DIVISOR_KJ_TO_WH: Final[float] = 3.6


# SAP 10.2 PDF p.107 Table N5 — additional days at longer heating duration
# for variable heating duration packages. Each row gives (PSR, N24,9, N16,9):
# days per year operating at 24 / 16 hours respectively, instead of the
# standard 9 hours/day. "Use linear interpolation for intermediate values
# of plant size ratio, rounding the result to the nearest whole number of
# days." Clamped to the table's bounds per the same convention as PSR
# efficiency interpolation (PDF p.101 lines 6007-6008).
_TABLE_N5_VARIABLE_HEATING_DAYS: Final[tuple[tuple[float, int, int], ...]] = (
    (0.2, 218, 6),
    (0.3, 191, 22),
    (0.4, 168, 29),
    (0.5, 128, 56),
    (0.6, 94, 74),
    (0.7, 50, 95),
    (0.8, 26, 103),
    (0.9, 14, 92),
    (1.0, 8, 77),
    (1.1, 4, 55),
    (1.2, 3, 38),
)

# SAP 10.2 Table 1a — calendar days per month (non-leap), used for
# Equation N5 + the N24,9 / N16,9 day allocation algorithm (PDF p.107).
_DAYS_IN_MONTH: Final[tuple[int, ...]] = (31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31)
# SAP 10.2 PDF p.107 — month indices (0-based Jan..Dec) in cold-to-warm
# order: Jan, Dec, Feb, Mar, Nov, Apr, Oct, May. The remaining four
# months (Jun, Jul, Aug, Sep) never receive any allocation — for cohort
# PSRs the year totals never exceed the sum of these eight cold months.
_TABLE_N5_ALLOCATION_ORDER: Final[tuple[int, ...]] = (0, 11, 1, 2, 10, 3, 9, 4)


def extended_heating_days_from_psr_variable(*, psr: float) -> tuple[int, int]:
    """SAP 10.2 Appendix N3.5 + Table N5 (PDF p.107) — for heat-pump
    packages with `heating_duration_code = "V"` (Variable), linearly
    interpolate annual N24,9 and N16,9 totals between the bracketing
    Table N5 rows at the dwelling's plant size ratio, rounding the
    result to the nearest whole number of days.

    Clamps to the table's bounds (PSR ≤ 0.2 → first row; PSR ≥ 1.2 →
    last row) per the same convention as PSR efficiency interpolation
    in Appendix N (PDF p.101 lines 6007-6008).

    For the legacy fixed durations:
      "24" → (365, 0)
      "16" → (0, 365)
      "9"  → (0, 0)
    Those branches are caller responsibilities (Table N4) — this helper
    only covers the Variable case (the only duration lodged on modern
    PCDB Table 362 records per footnote 48).
    """
    if psr <= _TABLE_N5_VARIABLE_HEATING_DAYS[0][0]:
        return (_TABLE_N5_VARIABLE_HEATING_DAYS[0][1], _TABLE_N5_VARIABLE_HEATING_DAYS[0][2])
    if psr >= _TABLE_N5_VARIABLE_HEATING_DAYS[-1][0]:
        return (_TABLE_N5_VARIABLE_HEATING_DAYS[-1][1], _TABLE_N5_VARIABLE_HEATING_DAYS[-1][2])
    for low, high in zip(
        _TABLE_N5_VARIABLE_HEATING_DAYS,
        _TABLE_N5_VARIABLE_HEATING_DAYS[1:],
    ):
        low_psr, low_n24, low_n16 = low
        high_psr, high_n24, high_n16 = high
        if low_psr <= psr <= high_psr:
            span = high_psr - low_psr
            t = (psr - low_psr) / span if span > 0 else 0.0
            n24_f = low_n24 + (high_n24 - low_n24) * t
            n16_f = low_n16 + (high_n16 - low_n16) * t
            return (round(n24_f), round(n16_f))
    raise AssertionError("PSR bracket not found despite range check")


def allocate_extended_heating_days_to_months(
    *,
    n24_9_year: int,
    n16_9_year: int,
) -> tuple[tuple[int, int], ...]:
    """SAP 10.2 Appendix N3.5 (PDF p.107) — distribute the annual N24,9
    and N16,9 day counts to months following the spec's allocation
    order ("Jan, Dec, Feb, Mar, Nov, Apr, Oct, May (coldest to the
    warmest)"). N24,9 days are filled first across the order, then
    N16,9 days fill the days not yet claimed by N24,9.

    Returns a length-12 tuple of `(N24,9_m, N16,9_m)` Jan..Dec. The
    summer months (Jun, Jul, Aug, Sep) always return (0, 0) for the
    PSR/duration cases this codebase handles — they're outside the
    allocation order.
    """
    n24_remaining = n24_9_year
    n16_remaining = n16_9_year
    allocations: list[tuple[int, int]] = [(0, 0)] * 12

    # Sweep 1 — N24,9: fill each cold month up to its day count.
    for m_idx in _TABLE_N5_ALLOCATION_ORDER:
        if n24_remaining <= 0:
            break
        month_days = _DAYS_IN_MONTH[m_idx]
        take = min(n24_remaining, month_days)
        allocations[m_idx] = (take, 0)
        n24_remaining -= take

    # Sweep 2 — N16,9: fill remaining month space (month_days − N24,m).
    for m_idx in _TABLE_N5_ALLOCATION_ORDER:
        if n16_remaining <= 0:
            break
        month_days = _DAYS_IN_MONTH[m_idx]
        n24_m, _ = allocations[m_idx]
        space = month_days - n24_m
        if space <= 0:
            continue
        take = min(n16_remaining, space)
        allocations[m_idx] = (n24_m, take)
        n16_remaining -= take

    return tuple(allocations)


def extended_zone_mean_temperature_c(
    *,
    heating_temperature_c: float,
    t_bimodal_c: float,
    t_unimodal_c: float,
    n24_9_m: int,
    n16_9_m: int,
    days_in_month: int,
) -> float:
    """SAP 10.2 Appendix N3.5 Equation N5 (PDF p.107) — blend a zone's
    monthly mean temperature across three heating patterns:

        T = [N24,9 × Th + N16,9 × T_uni + (Nm − N16,9 − N24,9) × T_bi] / Nm

    The three patterns differ in their daily off-period structure:
      - 24-hour day: zero off periods → T = Th
      - Unimodal (16-hour day): one off period of 8h (0700-2300 heating
        period per Table N7 footnote b) → T = Th − u1(8h)
      - Bimodal (9-hour day): two off periods (7+8h for living-area and
        elsewhere control-type 1/2; 9+8h for elsewhere control-type 3
        per Table N7) → T = Th − u1 − u2

    The caller passes the pre-computed `t_bimodal_c` and `t_unimodal_c`
    (already reduced from Th by the relevant Table 9b u terms); this
    helper just does the day-weighted blend.

    When `n24_9_m = n16_9_m = 0` the formula collapses to `t_bimodal_c`,
    so non-HP certs and warm months flow through unchanged.
    """
    if days_in_month <= 0:
        return t_bimodal_c
    bimodal_days = days_in_month - n16_9_m - n24_9_m
    return (
        n24_9_m * heating_temperature_c
        + n16_9_m * t_unimodal_c
        + bimodal_days * t_bimodal_c
    ) / days_in_month


def elsewhere_heating_temperature_c(
    *,
    heat_loss_parameter: float,
    control_type: int,
) -> float:
    """Rest-of-dwelling temperature during heating periods (Table 9).

    Control type 1 (e.g. boiler system without sufficient controls):
        T_h2 = 21 − 0.5 × HLP
    Control type 2 or 3 (programmer + room thermostat or better):
        T_h2 = 21 − HLP + HLP² / 12

    HLP is clamped at 6.0 W/m²K for this computation per Table 9 note (e).
    """
    hlp = min(heat_loss_parameter, _HLP_CLAMP_FOR_T_H2)
    if control_type == 1:
        return _T_H1_C - 0.5 * hlp
    return _T_H1_C - hlp + (hlp * hlp) / 12.0


def off_period_temperature_reduction_c(
    *,
    off_period_hours: float,
    heating_temperature_c: float,
    external_temperature_c: float,
    responsiveness: float,
    total_gains_w: float,
    heat_transfer_coefficient_w_per_k: float,
    utilisation_factor: float,
    time_constant_h: float,
) -> float:
    """SAP 10.2 Table 9b — temperature reduction `u` during a heating-off
    period, in °C below the heating-period temperature.

        t_c   = 4 + 0.25 × τ
        T_sc  = (1 − R)(T_h − 2) + R × (T_e + η·G / H)
        if t_off ≤ t_c:  u = 0.5 × t_off² × (T_h − T_sc) / (24 × t_c)
        if t_off >  t_c:  u = (T_h − T_sc) × (t_off − 0.5·t_c) / 24
    """
    t_c = 4.0 + 0.25 * time_constant_h
    h_safe = heat_transfer_coefficient_w_per_k if heat_transfer_coefficient_w_per_k > 0 else 1.0
    t_sc = (1.0 - responsiveness) * (heating_temperature_c - 2.0) + responsiveness * (
        external_temperature_c + utilisation_factor * total_gains_w / h_safe
    )
    delta_t = heating_temperature_c - t_sc
    if off_period_hours <= t_c:
        return 0.5 * off_period_hours * off_period_hours * delta_t / (24.0 * t_c)
    return delta_t * (off_period_hours - 0.5 * t_c) / 24.0


# Standard SAP heating schedule (Table 9):
# Living area:               off (7, 8) regardless of control type.
# Elsewhere control type 1,2: off (7, 8).
# Elsewhere control type 3:   off (9, 8).
_LIVING_AREA_OFF_HOURS: Final[tuple[float, float]] = (7.0, 8.0)
_ELSEWHERE_OFF_HOURS_TYPE_12: Final[tuple[float, float]] = (7.0, 8.0)
_ELSEWHERE_OFF_HOURS_TYPE_3: Final[tuple[float, float]] = (9.0, 8.0)
# SAP 10.2 Appendix N3.5 Table N7 footnote (b): "heating 0700-2300" =
# 16 hours on, one 8-hour off period for both zones regardless of
# control type. Used by Equation N5's T_unimodal term.
_UNIMODAL_OFF_HOURS: Final[float] = 8.0


@dataclass(frozen=True)
class MeanInternalTemperatureResult:
    """SAP 10.2 §7 worksheet line refs (85)..(94).

    Returned by `mean_internal_temperature_monthly`. Downstream §8 space
    heating consumes `adjusted_mean_internal_temp_monthly` (93) and
    `utilisation_factor_whole_monthly` (94) directly; per-zone tuples are
    exposed for worksheet conformance + audit.
    """

    living_area_heating_temp_c: float                            # (85)
    utilisation_factor_living_monthly: tuple[float, ...]         # (86)
    mean_internal_temp_living_monthly: tuple[float, ...]         # (87)
    elsewhere_heating_temp_monthly: tuple[float, ...]            # (88)
    utilisation_factor_elsewhere_monthly: tuple[float, ...]      # (89)
    mean_internal_temp_elsewhere_monthly: tuple[float, ...]      # (90)
    living_area_fraction: float                                  # (91)
    mean_internal_temp_monthly: tuple[float, ...]                # (92)
    adjusted_mean_internal_temp_monthly: tuple[float, ...]       # (93)
    utilisation_factor_whole_monthly: tuple[float, ...]          # (94)


def _zone_mean_temp_with_per_zone_eta(
    *,
    heating_temperature_c: float,
    off_hours_first: float,
    off_hours_second: float,
    external_temp_c: float,
    responsiveness: float,
    total_gains_w: float,
    heat_transfer_coefficient_w_per_k: float,
    time_constant_h: float,
) -> tuple[float, float]:
    """SAP 10.2 Table 9c steps 2–4 (or 5–6) for one zone.

    Computes η using L = H(Th − Te) — the zone's own heating temperature —
    then folds it into Table 9b u1+u2 to yield zone mean. Returns (η, T_zone).
    Distinct from `_zone_mean_temperature_c` which takes a precomputed η
    that's shared across zones (the pre-fix shape).
    """
    loss_rate_w = max(0.0, heat_transfer_coefficient_w_per_k * (heating_temperature_c - external_temp_c))
    eta = utilisation_factor(
        total_gains_w=total_gains_w,
        heat_loss_rate_w=loss_rate_w,
        time_constant_h=time_constant_h,
    )
    common = dict(
        heating_temperature_c=heating_temperature_c,
        external_temperature_c=external_temp_c,
        responsiveness=responsiveness,
        total_gains_w=total_gains_w,
        heat_transfer_coefficient_w_per_k=heat_transfer_coefficient_w_per_k,
        utilisation_factor=eta,
        time_constant_h=time_constant_h,
    )
    u1 = off_period_temperature_reduction_c(off_period_hours=off_hours_first, **common)
    u2 = off_period_temperature_reduction_c(off_period_hours=off_hours_second, **common)
    return eta, heating_temperature_c - u1 - u2


def mean_internal_temperature_monthly(
    *,
    monthly_external_temp_c: tuple[float, ...],
    monthly_total_gains_w: tuple[float, ...],
    monthly_heat_transfer_coefficient_w_per_k: tuple[float, ...],
    thermal_mass_parameter_kj_per_m2_k: float,
    total_floor_area_m2: float,
    control_type: int,
    responsiveness: float,
    living_area_fraction: float,
    control_temperature_adjustment_c: float = 0.0,
    secondary_fraction: float = 0.0,
    secondary_responsiveness: float = 1.0,
    extended_heating_days_per_month: Optional[tuple[tuple[int, int], ...]] = None,
) -> MeanInternalTemperatureResult:
    """SAP 10.2 §7 orchestrator — chain Table 9c steps 1–9 for all 12 months.

    Per-zone η is the load-bearing detail vs the legacy single-η path:
      - (86) η_living = f(Ti = T_h1 = 21°C)
      - (89) η_elsewhere = f(Ti = T_h2 from Table 9)
      - (94) η_whole = f(Ti = (93) adjusted MIT) — feeds §8 step 9.

    Inputs:
      monthly_*           length-12 tuples covering Jan..Dec
      thermal_mass_parameter_kj_per_m2_k  TMP (35) for τ = TMP·TFA/(3.6·HLC)
      total_floor_area_m2                 (4) for τ derivation
      control_type        1/2/3 buckets feed Table 9 T_h2 + Table 9 off-hours
      responsiveness      Table 4d R for the Table 9b u-formula
      living_area_fraction (91) for (92) zone blending
      control_temperature_adjustment_c (93) Table 4e adj (defaults 0 — all 6
                          Elmhurst fixtures have (93) = (92), so this is a
                          known shortcut; cert-side mapping is a future slice).
      secondary_fraction  (203) fraction of main heat from second main system
                          when both heat the whole house (Table 9c case 1).
                          Defaults 0 (single-main); used to compute weighted R
                          per Table 9b: R_eff = (1-frac)·R_primary + frac·R_secondary.
                          Case 2 (different parts heated) deferred — no fixture.
      extended_heating_days_per_month  SAP 10.2 Appendix N3.5 — 12-tuple of
                          (N24,9_m, N16,9_m) for heat-pump packages with
                          PCDB data. When provided, the orchestrator
                          replaces Table 9c steps 3-4 with Equation N5
                          (blending Th, T_unimodal, T_bimodal). When
                          None (the default — every non-HP cert), the
                          standard SAP heating schedule applies: T_zone
                          = T_bimodal directly.
    """
    effective_responsiveness = (
        (1.0 - secondary_fraction) * responsiveness
        + secondary_fraction * secondary_responsiveness
    )
    elsewhere_off_hours = (
        _ELSEWHERE_OFF_HOURS_TYPE_3 if control_type == 3 else _ELSEWHERE_OFF_HOURS_TYPE_12
    )

    eta_living: list[float] = []
    t_1: list[float] = []
    t_h2: list[float] = []
    eta_elsewhere: list[float] = []
    t_2: list[float] = []
    t_internal: list[float] = []
    t_adj: list[float] = []
    eta_whole: list[float] = []

    for m in _MONTHS:
        ext = monthly_external_temp_c[m]
        gains = monthly_total_gains_w[m]
        h = monthly_heat_transfer_coefficient_w_per_k[m]
        hlp = h / total_floor_area_m2 if total_floor_area_m2 > 0 else 0.0
        tau = (
            thermal_mass_parameter_kj_per_m2_k * total_floor_area_m2
            / (_TIME_CONSTANT_DIVISOR_KJ_TO_WH * h)
            if h > 0 else float("inf")
        )

        # Living area — steps 1-4
        eta_l, t_l_bimodal = _zone_mean_temp_with_per_zone_eta(
            heating_temperature_c=_T_H1_C,
            off_hours_first=_LIVING_AREA_OFF_HOURS[0],
            off_hours_second=_LIVING_AREA_OFF_HOURS[1],
            external_temp_c=ext, responsiveness=effective_responsiveness,
            total_gains_w=gains, heat_transfer_coefficient_w_per_k=h,
            time_constant_h=tau,
        )
        eta_living.append(eta_l)

        # Elsewhere — steps 5-6
        t_h2_m = elsewhere_heating_temperature_c(
            heat_loss_parameter=hlp, control_type=control_type,
        )
        t_h2.append(t_h2_m)
        eta_e, t_e_bimodal = _zone_mean_temp_with_per_zone_eta(
            heating_temperature_c=t_h2_m,
            off_hours_first=elsewhere_off_hours[0],
            off_hours_second=elsewhere_off_hours[1],
            external_temp_c=ext, responsiveness=effective_responsiveness,
            total_gains_w=gains, heat_transfer_coefficient_w_per_k=h,
            time_constant_h=tau,
        )
        eta_elsewhere.append(eta_e)

        # SAP 10.2 Appendix N3.5 Equation N5 — when the caller provides
        # per-month (N24,9, N16,9) day allocations, blend Th / T_unimodal
        # / T_bimodal for each zone. T_unimodal applies one 8-hour off
        # period per Table N7 footnote (b) at the same η as the bimodal
        # path (η depends on Th + HLC only, not the heating schedule).
        if extended_heating_days_per_month is not None:
            n24_m, n16_m = extended_heating_days_per_month[m]
            days_m = _DAYS_IN_MONTH[m]
            u_uni_living = off_period_temperature_reduction_c(
                off_period_hours=_UNIMODAL_OFF_HOURS,
                heating_temperature_c=_T_H1_C,
                external_temperature_c=ext,
                responsiveness=effective_responsiveness,
                total_gains_w=gains,
                heat_transfer_coefficient_w_per_k=h,
                utilisation_factor=eta_l,
                time_constant_h=tau,
            )
            t_l_unimodal = _T_H1_C - u_uni_living
            u_uni_elsewhere = off_period_temperature_reduction_c(
                off_period_hours=_UNIMODAL_OFF_HOURS,
                heating_temperature_c=t_h2_m,
                external_temperature_c=ext,
                responsiveness=effective_responsiveness,
                total_gains_w=gains,
                heat_transfer_coefficient_w_per_k=h,
                utilisation_factor=eta_e,
                time_constant_h=tau,
            )
            t_e_unimodal = t_h2_m - u_uni_elsewhere
            t_l = extended_zone_mean_temperature_c(
                heating_temperature_c=_T_H1_C,
                t_bimodal_c=t_l_bimodal,
                t_unimodal_c=t_l_unimodal,
                n24_9_m=n24_m,
                n16_9_m=n16_m,
                days_in_month=days_m,
            )
            t_e = extended_zone_mean_temperature_c(
                heating_temperature_c=t_h2_m,
                t_bimodal_c=t_e_bimodal,
                t_unimodal_c=t_e_unimodal,
                n24_9_m=n24_m,
                n16_9_m=n16_m,
                days_in_month=days_m,
            )
        else:
            t_l = t_l_bimodal
            t_e = t_e_bimodal
        t_1.append(t_l)
        t_2.append(t_e)

        # Blend (step 7) + Table 4e adj (step 8)
        t_int = living_area_fraction * t_l + (1.0 - living_area_fraction) * t_e
        t_internal.append(t_int)
        t_adjusted = t_int + control_temperature_adjustment_c
        t_adj.append(t_adjusted)

        # Step 9 — η_whole at the adjusted MIT for §8 space heating
        loss_whole_w = max(0.0, h * (t_adjusted - ext))
        eta_whole.append(
            utilisation_factor(
                total_gains_w=gains,
                heat_loss_rate_w=loss_whole_w,
                time_constant_h=tau,
            )
        )

    return MeanInternalTemperatureResult(
        living_area_heating_temp_c=_T_H1_C,
        utilisation_factor_living_monthly=tuple(eta_living),
        mean_internal_temp_living_monthly=tuple(t_1),
        elsewhere_heating_temp_monthly=tuple(t_h2),
        utilisation_factor_elsewhere_monthly=tuple(eta_elsewhere),
        mean_internal_temp_elsewhere_monthly=tuple(t_2),
        living_area_fraction=living_area_fraction,
        mean_internal_temp_monthly=tuple(t_internal),
        adjusted_mean_internal_temp_monthly=tuple(t_adj),
        utilisation_factor_whole_monthly=tuple(eta_whole),
    )