Model/domain/sap10_calculator/worksheet/water_heating.py

"""SAP 10.2 §4 — water heating energy requirements.

Worksheet line refs (xlsx rows 207-304, sheet `NonRegionalWeather`):
  (42)    assumed occupancy N from Appendix J Table 1b
  (42a)m  monthly hot water usage for mixer showers
  (42b)m  monthly hot water usage for baths
  (42c)m  monthly hot water usage for other uses
  (43)    annual average hot water usage (litres/day)
  (44)m   daily hot water usage by month
  (45)m   energy content of monthly hot water demand
  (46)m   distribution loss = 0.15 × (45)m
  (47)–(56)  storage volume + water storage / HIU loss
  (57)    dedicated solar storage adjustment
  (58)–(61)  primary loss + combi loss
  (62)m   total monthly water heat requirement
  (63a)–(63d)  WWHRS / PV-diverter / Solar / FGHRS reductions
  (64)m   output from water heater
  (65)m   heat gains from water heating

Reference: SAP 10.2 specification §4 (pages 22-31) + Appendix J (pages
84-90); canonical xlsx worked example at the repo root.
"""

from __future__ import annotations

from dataclasses import dataclass
from math import exp
from typing import Final, Literal, Optional

from datatypes.epc.domain.epc_property_data import EpcPropertyData


_OCCUPANCY_TFA_FLOOR_M2: Final[float] = 13.9


@dataclass(frozen=True)
class WaterHeatingResult:
    """SAP 10.2 §4 worksheet outputs broken down per line ref so callers
    can audit the cascade against the canonical xlsx or an Elmhurst
    worksheet. Annual totals are days-weighted where appropriate.

    Field-to-line-ref mapping:
      (42)   occupancy
      (43)   annual_avg_hot_water_l_per_day
      (44)m  daily_hot_water_l_per_day_monthly
      (45)m  energy_content_monthly_kwh
      (46)m  distribution_loss_monthly_kwh
      (61)m  combi_loss_monthly_kwh
      (62)m  total_demand_monthly_kwh
      (64)m  output_monthly_kwh
      (65)m  heat_gains_monthly_kwh
    Annual sum of (64)m is exposed as `output_kwh_per_yr` for the
    calculator's `hot_water_kwh_per_yr` slot.
    """

    occupancy: float
    annual_avg_hot_water_l_per_day: float
    daily_hot_water_l_per_day_monthly: tuple[float, ...]
    energy_content_monthly_kwh: tuple[float, ...]
    distribution_loss_monthly_kwh: tuple[float, ...]
    solar_storage_monthly_kwh: tuple[float, ...]  # (57)m — Tables 2/2a/2b
    primary_loss_monthly_kwh: tuple[float, ...]   # (59)m — Table 3
    combi_loss_monthly_kwh: tuple[float, ...]
    total_demand_monthly_kwh: tuple[float, ...]
    output_monthly_kwh: tuple[float, ...]
    heat_gains_monthly_kwh: tuple[float, ...]
    electric_shower_monthly_kwh: tuple[float, ...]  # (64a)m — App J step 8
    output_kwh_per_yr: float
    electric_shower_kwh_per_yr: float  # Σ (64a)m — feeds §10a (247a)

# Table J2 — monthly factors for hot water use (also used by Appendix J
# equation J11 for "other uses"). Symmetric about the year midpoint.
_TABLE_J2_MONTHLY_FACTORS: Final[tuple[float, ...]] = (
    1.10, 1.06, 1.02, 0.98, 0.94, 0.90, 0.90, 0.94, 0.98, 1.02, 1.06, 1.10,
)

# Appendix J J11 footnote: -5% reduction in V_d,other,ave for dwellings
# designed to ≤125 L/person/day total water use.
_LOW_WATER_USE_REDUCTION: Final[float] = 0.05

# Table J1 — cold-water inlet temperatures (°C). Two columns per the spec:
# from a header tank (cooler) vs from the mains. The cert lodges which
# applies; both populations need spec-exact monthly arrays.
TABLE_J1_TCOLD_FROM_HEADER_TANK_C: Final[tuple[float, ...]] = (
    11.1, 11.3, 12.3, 14.5, 16.2, 18.8, 21.3, 19.3, 18.7, 16.2, 13.2, 11.2,
)
TABLE_J1_TCOLD_FROM_MAINS_C: Final[tuple[float, ...]] = (
    8.0, 8.2, 9.3, 12.7, 14.6, 16.7, 18.4, 17.6, 16.6, 14.3, 11.1, 8.5,
)

# Table J5 — behavioural variation factor for showers AND baths. Used by
# (42a)m showers (Appendix J step 1d) and (42b)m baths (step 2b) alike.
TABLE_J5_BEHAVIOURAL_FACTOR: Final[tuple[float, ...]] = (
    1.035, 1.021, 1.007, 0.993, 0.979, 0.965,
    0.965, 0.979, 0.993, 1.007, 1.021, 1.035,
)

# Appendix J equation J7 constants.
_BATH_VOLUME_L: Final[float] = 73.0
# Appendix J equation J8 / J3 mixing temperatures.
_HOT_DELIVERY_TEMPERATURE_C: Final[float] = 52.0
_WARM_BATH_TEMPERATURE_C: Final[float] = 42.0
_WARM_SHOWER_TEMPERATURE_C: Final[float] = 41.0
# Appendix J step 1d: assumed shower duration in minutes per event.
_SHOWER_DURATION_MIN: Final[float] = 6.0

# Days per month (non-leap year) — used by Appendix J J4 / J9 / J12 to
# turn monthly daily-rate arrays into days-weighted annual averages.
_DAYS_IN_MONTH: Final[tuple[int, ...]] = (
    31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31,
)
_DAYS_IN_YEAR: Final[int] = sum(_DAYS_IN_MONTH)


def assumed_occupancy(total_floor_area_m2: float) -> float:
    """SAP 10.2 §4 line (42) / Appendix J Table 1b.

    Piecewise occupancy by total floor area:
        TFA ≤ 13.9 m² : N = 1
        TFA  > 13.9 m² : N = 1 + 1.76 × (1 − exp(−0.000349 × (TFA − 13.9)²))
                              + 0.0013 × (TFA − 13.9)
    """
    if total_floor_area_m2 <= _OCCUPANCY_TFA_FLOOR_M2:
        return 1.0
    x = total_floor_area_m2 - _OCCUPANCY_TFA_FLOOR_M2
    return 1.0 + 1.76 * (1.0 - exp(-0.000349 * x * x)) + 0.0013 * x


def hot_water_mixer_showers_monthly_l_per_day(
    *,
    n_occupants: float,
    has_bath: bool,
    mixer_shower_flow_rates_l_per_min: tuple[float, ...],
    cold_water_temps_c: tuple[float, ...],
    n_electric_showers: int = 0,
) -> tuple[float, ...]:
    """SAP 10.2 §4 line (42a)m via Appendix J equations J1, J2, J3.

    Mixer showers draw from the main hot water system (cylinder or combi)
    and mix with cold to a 41 °C delivery. The per-day shower count is

        N_shower = 0.45 × N + 0.65    (bath present)
                 = 0.58 × N + 0.83    (no bath)
                 = 0                  (no shower outlets)

    Each outlet's warm-water draw for a month is the flow rate × 6 min ×
    Table J5 fbeh. The hot fraction is (41 − Tcold[m])/(52 − Tcold[m]).
    Per-outlet daily warm water is scaled by N_shower / N_outlets, then
    summed across the MIXER outlets to give (42a)m.

    N_outlets is the dwelling's TOTAL shower-outlet count: SAP 10.2
    Appendix J step 1a is explicit that the count INCLUDES "any
    instantaneous electric showers" (which then bill their own energy via
    (64a)m). So `mixer_shower_flow_rates_l_per_min` lists only the mixer
    outlets (iterated for the warm-water draw), but `n_electric_showers`
    must be added to the denominator — otherwise a dwelling with both a
    mixer and an electric shower assigns the FULL N_shower to the mixer
    system AND bills the electric shower on top, double-counting shower
    demand (over-counts main HW → under-rates the dwelling).
    """
    n_mixer_outlets = len(mixer_shower_flow_rates_l_per_min)
    if n_mixer_outlets == 0:
        return tuple(0.0 for _ in range(12))
    # Appendix J step 1a: Noutlets includes instantaneous electric showers.
    n_outlets = n_mixer_outlets + n_electric_showers
    if has_bath:
        n_shower = 0.45 * n_occupants + 0.65
    else:
        n_shower = 0.58 * n_occupants + 0.83
    monthly: list[float] = []
    for fbeh, tcold in zip(TABLE_J5_BEHAVIOURAL_FACTOR, cold_water_temps_c):
        f_hot = (_WARM_SHOWER_TEMPERATURE_C - tcold) / (
            _HOT_DELIVERY_TEMPERATURE_C - tcold
        )
        v_hot_total = 0.0
        for flow in mixer_shower_flow_rates_l_per_min:
            v_warm_per_outlet = flow * _SHOWER_DURATION_MIN * fbeh
            v_d_warm = v_warm_per_outlet * n_shower / n_outlets
            v_hot_total += v_d_warm * f_hot
        monthly.append(v_hot_total)
    return tuple(monthly)


def hot_water_baths_monthly_l_per_day(
    *,
    n_occupants: float,
    has_bath: bool,
    has_shower: bool,
    cold_water_temps_c: tuple[float, ...],
    low_water_use: bool,
) -> tuple[float, ...]:
    """SAP 10.2 §4 line (42b)m via Appendix J equations J6, J7, J8.

    Per-day hot water for bath fills across the 12 months. Per J6:

        N_bath = 0                   if no bath (but a shower is present)
               = 0.13 × N + 0.19     if shower is also present
               = 0.35 × N + 0.50     if no shower present, or no bath
                                     and no shower

    Each bath fills 73 L of warm water at 42 °C; the hot fraction depends
    on monthly cold-water temperature (Table J1, either header tank or
    mains depending on the cert). `low_water_use` knocks 5% off the
    warm-water term per the J7 footnote.

    `cold_water_temps_c` must be a 12-tuple of monthly Tcold values —
    pass `TABLE_J1_TCOLD_FROM_MAINS_C` for the common case.
    """
    if not has_bath and has_shower:
        return tuple(0.0 for _ in range(12))
    if has_bath and has_shower:
        n_bath = 0.13 * n_occupants + 0.19
    else:
        n_bath = 0.35 * n_occupants + 0.50
    lwu_factor = 1.0 - _LOW_WATER_USE_REDUCTION if low_water_use else 1.0
    monthly: list[float] = []
    for fbeh, tcold in zip(TABLE_J5_BEHAVIOURAL_FACTOR, cold_water_temps_c):
        v_warm = n_bath * _BATH_VOLUME_L * fbeh * lwu_factor
        f_hot = (_WARM_BATH_TEMPERATURE_C - tcold) / (
            _HOT_DELIVERY_TEMPERATURE_C - tcold
        )
        monthly.append(v_warm * f_hot)
    return tuple(monthly)


def total_hot_water_monthly_l_per_day(
    *,
    showers: tuple[float, ...],
    baths: tuple[float, ...],
    other_uses: tuple[float, ...],
) -> tuple[float, ...]:
    """SAP 10.2 §4 line (44)m via Appendix J equation J13:
        V_d,m = V_d,shower[m] + V_d,bath[m] + V_d,other[m]

    A pure element-wise sum of the three monthly demand streams. All
    three inputs must be 12-tuples — caller is responsible for ensuring
    they were computed against the same Tcold table.
    """
    return tuple(s + b + o for s, b, o in zip(showers, baths, other_uses))


def energy_content_of_hot_water_monthly_kwh(
    *,
    monthly_hot_water_l_per_day: tuple[float, ...],
    cold_water_temps_c: tuple[float, ...],
) -> tuple[float, ...]:
    """SAP 10.2 §4 line (45)m via Appendix J equation J14:
        (45)m = 4.18 × V_d,m × n_m × (52 − Tcold[m]) / 3600    [kWh/month]

    Sensible heat to raise the monthly hot water volume from Tcold[m] to
    the 52 °C delivery temperature. 4.18 J/(g·K) is the specific heat of
    water; dividing by 3600 converts J/g to Wh/g (= kWh/kg, since 1 L of
    water ≈ 1 kg).
    """
    return tuple(
        4.18 * vd * n * (_HOT_DELIVERY_TEMPERATURE_C - tcold) / 3600.0
        for vd, n, tcold in zip(
            monthly_hot_water_l_per_day, _DAYS_IN_MONTH, cold_water_temps_c
        )
    )


def distribution_loss_monthly_kwh(
    *,
    monthly_energy_content_kwh: tuple[float, ...],
    is_instantaneous_at_point_of_use: bool,
) -> tuple[float, ...]:
    """SAP 10.2 §4 line (46)m.

    Distribution loss is a flat 15% of (45)m unless the water heating is
    instantaneous at the point of use (Table 4a hot water codes 907 and
    909) — those have no distribution pipework so the loss is zero. Heat
    networks still incur the 15% loss whether or not a hot water cylinder
    is present (spec §4 step 7).
    """
    if is_instantaneous_at_point_of_use:
        return tuple(0.0 for _ in range(12))
    return tuple(0.15 * e for e in monthly_energy_content_kwh)


def water_efficiency_monthly_via_equation_d1(
    *,
    winter_efficiency_pct: float,
    summer_efficiency_pct: float,
    space_heating_monthly_useful_kwh: tuple[float, ...],
    water_heating_output_monthly_kwh: tuple[float, ...],
) -> tuple[float, ...]:
    """SAP 10.2 Appendix D §D2.1 (2) Equation D1 — monthly water-heating
    efficiency cascade for combi boilers and CPSUs that provide both
    space and water heating:

        η_water,monthly = (Q_space + Q_water) / (Q_space/η_winter + Q_water/η_summer)

    where Q_space (kWh/month) = (98c)m × (204) and Q_water (kWh/month)
    = (64)m. η_winter is the raw PCDB winter seasonal efficiency
    (Appendix D §D2.1 (2) note: "η_winter does not include any
    efficiency adjustment due to design flow temperature or controls").

    Two early-out rules per spec:
    - If summer_efficiency ≥ winter_efficiency (or the boiler is water-
      heating-only): η_water,monthly = η_summer for every month.
    - If both Q_space[m] and Q_water[m] = 0 in any month: η_water,
      monthly[m] = η_summer.
    """
    if summer_efficiency_pct >= winter_efficiency_pct:
        return (summer_efficiency_pct / 100.0,) * 12
    eff_winter = winter_efficiency_pct / 100.0
    eff_summer = summer_efficiency_pct / 100.0
    monthly: list[float] = []
    for q_space, q_water in zip(
        space_heating_monthly_useful_kwh, water_heating_output_monthly_kwh
    ):
        if q_space == 0.0 and q_water == 0.0:
            monthly.append(eff_summer)
            continue
        numerator = q_space + q_water
        denominator = q_space / eff_winter + q_water / eff_summer
        monthly.append(numerator / denominator)
    return tuple(monthly)


def combi_loss_monthly_kwh_table_3b_row_1_instantaneous(
    *,
    rejected_energy_proportion_r1: float,
    loss_factor_f1_kwh_per_day: float,
    energy_content_monthly_kwh: tuple[float, ...],
    daily_hot_water_monthly_l_per_day: tuple[float, ...],
) -> tuple[float, ...]:
    """SAP 10.2 Appendix J Table 3b row 1 (Instantaneous combi with non-
    storage FGHRS or without FGHRS, profile M only):

        (61)m = (45)m × r1 × fu + [F1 × n_m]

    where r1 = rejected energy proportion (PCDB Table 105 field 51),
    F1 = loss factor in kWh/day (PCDB field 52), and fu = V_d,m / 100
    when daily hot-water usage V_d,m < 100 L/day, else fu = 1.0.

    Applies only to combi boilers EN 13203-2 / OPS 26 tested with one
    profile (Separate DHW tests = 1). Other Table 3b rows (storage
    combis, storage-FGHRS variants) and Table 3c (two-profile tests)
    are deferred until a fixture exercises them. Untested combis fall
    back to the existing Table 3a path.
    """
    return tuple(
        e * rejected_energy_proportion_r1 * (v / 100.0 if v < 100.0 else 1.0)
        + loss_factor_f1_kwh_per_day * n
        for e, v, n in zip(
            energy_content_monthly_kwh,
            daily_hot_water_monthly_l_per_day,
            _DAYS_IN_MONTH,
        )
    )


_DVF_M_AND_L_LOWER_V_L_PER_DAY: Final[float] = 100.2
_DVF_M_AND_L_UPPER_V_L_PER_DAY: Final[float] = 199.8
_DVF_M_AND_L_UPPER_CLAMP: Final[float] = -99.6
_DVF_M_AND_S_LOWER_V_L_PER_DAY: Final[float] = 36.0
_DVF_M_AND_S_UPPER_V_L_PER_DAY: Final[float] = 100.2
_DVF_M_AND_S_LOWER_CLAMP: Final[float] = 64.2
_DVF_LINEAR_INTERCEPT: Final[float] = 100.2


def _table_3c_dvf(
    daily_hot_water_l_per_day: float,
    profile_pair: Literal["M+L", "M+S"],
) -> float:
    """SAP 10.2 Appendix J Table 3c — Daily Volume Factor.

    Piecewise function of V_d,m gated on which two EN 13203-2 / OPS 26
    profiles the combi was tested with (encoded in PCDF Table 105 field
    48 = separate_dhw_tests: 2 → schedules 2 and 3 = M+L; 3 → schedules
    2 and 1 = M+S).
    """
    v = daily_hot_water_l_per_day
    if profile_pair == "M+L":
        if v < _DVF_M_AND_L_LOWER_V_L_PER_DAY:
            return 0.0
        if v > _DVF_M_AND_L_UPPER_V_L_PER_DAY:
            return _DVF_M_AND_L_UPPER_CLAMP
        return _DVF_LINEAR_INTERCEPT - v
    if v < _DVF_M_AND_S_LOWER_V_L_PER_DAY:
        return _DVF_M_AND_S_LOWER_CLAMP
    if v > _DVF_M_AND_S_UPPER_V_L_PER_DAY:
        return 0.0
    return _DVF_LINEAR_INTERCEPT - v


def combi_loss_monthly_kwh_table_3c_two_profile_instantaneous(
    *,
    rejected_energy_proportion_r1: float,
    loss_factor_f2_kwh_per_day: float,
    rejected_factor_f3_per_litre: float,
    profile_pair: Literal["M+L", "M+S"],
    energy_content_monthly_kwh: tuple[float, ...],
    daily_hot_water_monthly_l_per_day: tuple[float, ...],
) -> tuple[float, ...]:
    """SAP 10.2 Appendix J Table 3c row 1 (Instantaneous combi with non-
    storage FGHRS or without FGHRS), two-profile tests:

        (61)m = (45)m × [r1 + DVF × F3] × fu + [F2 × n_m]

    where r1, F2 (kWh/day), F3 (litres⁻¹, can be negative) are PCDF
    Table 105 fields 51 / 56 / 57, DVF is the per-month Daily Volume
    Factor from `_table_3c_dvf`, and fu = V_d,m / 100 when V_d,m < 100
    else 1.0.

    Applies to PCDB combis with separate_dhw_tests ∈ {2, 3}: =2 means
    schedules 2 and 3 (profiles M and L); =3 means schedules 2 and 1
    (profiles M and S). The profile pair selects the DVF piecewise
    branch — see `_table_3c_dvf`. r2 (PCDF field 55) is lodged but the
    spec explicitly excludes it from SAP assessments ("only r1").

    Storage-FGHRS / storage-combi variants (Table 3c rows 2-5) are
    deferred until a fixture exercises them — mirrors the row-1-only
    coverage of Table 3b (see `pcdb_combi_loss_override`).
    """
    return tuple(
        e
        * (rejected_energy_proportion_r1 + _table_3c_dvf(v, profile_pair) * rejected_factor_f3_per_litre)
        * (v / 100.0 if v < 100.0 else 1.0)
        + loss_factor_f2_kwh_per_day * n
        for e, v, n in zip(
            energy_content_monthly_kwh,
            daily_hot_water_monthly_l_per_day,
            _DAYS_IN_MONTH,
        )
    )


def combi_loss_monthly_kwh_table_3a_keep_hot_time_clock() -> tuple[float, ...]:
    """SAP 10.2 §4 line (61)m — Table 3a row "Instantaneous, with keep-hot
    facility controlled by time clock": 600 × n_m / 365 kWh/month.

    A flat 600 kWh/year combi loss, prorated by month length. Unlike the
    "without keep-hot" rows there is no `fu` adjustment for low-volume
    draw — the loss is constant. Suitable for any standard non-PCDB-
    tested combi that the cert lodges as having a time-clock keep-hot
    facility (Appendix D section D1.16).
    """
    return tuple(600.0 * n / _DAYS_IN_YEAR for n in _DAYS_IN_MONTH)


def combi_loss_monthly_kwh_table_3a_row_1_no_keep_hot(
    *,
    daily_hot_water_monthly_l_per_day: tuple[float, ...],
) -> tuple[float, ...]:
    """SAP 10.2 §4 line (61)m — Table 3a row 1 "Instantaneous, without
    keep-hot facility": 600 × fu × n_m / 365 kWh/month, where fu = V_d,m
    / 100 when V_d,m < 100 L/day, else fu = 1.0 (SAP 10.2 spec p.160).

    Differs from the keep-hot time-clock row by the fu volume-scaling
    factor — for low-volume dwellings (V_d < 100 L/day on average ≈ N <
    2.5 occupants with no electric showers) the loss is proportionally
    less than 600 kWh/yr. For V_d ≥ 100 every month, fu collapses to 1.0
    and this row coincides with `..._keep_hot_time_clock()` (600 kWh/yr
    flat).

    Origin: BRE STP09-B04 §5.3 derived the 600 kWh/yr keep-hot baseline
    from observed cycling losses; the no-keep-hot variant scales by fu
    because instantaneous combis only cycle when actually drawing hot
    water, and low-draw dwellings stand idle.
    """
    return tuple(
        600.0 * min(1.0, v / 100.0) * n / _DAYS_IN_YEAR
        for v, n in zip(daily_hot_water_monthly_l_per_day, _DAYS_IN_MONTH)
    )


def combi_loss_monthly_kwh_table_3a_row_4_keep_hot_no_time_clock() -> tuple[float, ...]:
    """SAP 10.2 §4 line (61)m — Table 3a row "Instantaneous, with keep-hot
    facility not controlled by time clock": 900 × n_m / 365 kWh/month
    (SAP 10.2 spec p.160).

    A flat 900 kWh/year — 50% larger than the time-clocked variant
    because the keep-hot heater cycles around the clock rather than only
    during scheduled windows. No fu adjustment per spec: the keep-hot
    facility maintains store temperature regardless of draw.
    """
    return tuple(900.0 * n / _DAYS_IN_YEAR for n in _DAYS_IN_MONTH)


# SAP 10.2 Table 2 (PDF p.158) hot water storage loss factor L kWh/litre/day.
# Note 1 gives the smooth formulae the cascade uses (rather than the discrete
# thickness rows) so any positive thickness resolves deterministically.
_CYLINDER_INSULATION_FACTORY = "factory_insulated"
_CYLINDER_INSULATION_LOOSE_JACKET = "loose_jacket"


def cylinder_storage_loss_factor_table_2(
    *,
    insulation_type: Literal["factory_insulated", "loose_jacket"],
    thickness_mm: float,
) -> float:
    """SAP 10.2 Table 2 (PDF p.158) — hot water storage loss factor L
    in kWh/litre/day. Note 1 supplies the smooth formula:
        Cylinder, factory insulated: L = 0.005 + 0.55 / (t + 4.0)
        Cylinder, loose jacket:      L = 0.005 + 1.76 / (t + 12.8)
    where t is the insulation thickness in mm. Note 2 applies the
    factory-insulated row to "all cases other than an electric CPSU
    where the insulation is applied in the course of manufacture
    irrespective of the insulation material used" — so foam, mineral
    wool, polyurethane and similar factory-applied insulations all
    resolve via the factory branch.
    """
    if insulation_type == _CYLINDER_INSULATION_FACTORY:
        return 0.005 + 0.55 / (thickness_mm + 4.0)
    return 0.005 + 1.76 / (thickness_mm + 12.8)


def cylinder_volume_factor_table_2a(volume_l: float) -> float:
    """SAP 10.2 Table 2a (PDF p.158) — volume factor VF using Note 2's
    closed form `VF = (120 / Vc)^(1/3)`. The closed form matches the
    tabulated rows to 4 d.p. (V=160 → VF=0.9086 in the worksheet vs the
    table's 0.908 — Elmhurst computes via formula).
    """
    return (120.0 / volume_l) ** (1.0 / 3.0)


def cylinder_temperature_factor_table_2b(
    *,
    has_cylinder_thermostat: bool,
    separately_timed_dhw: bool,
) -> float:
    """SAP 10.2 Table 2b (PDF p.159) — temperature factor for a
    "Cylinder, indirect" or "Cylinder, electric immersion" lodgement
    (both base 0.60 in the "loss from Table 2" column). Multipliers per
    Notes a) / b):
        × 1.3 if cylinder thermostat is absent
        × 0.9 if domestic hot water is separately timed
    """
    factor = 0.60
    if not has_cylinder_thermostat:
        factor *= 1.3
    if separately_timed_dhw:
        factor *= 0.9
    return factor


# SAP 10.2 Table 3 (PDF p.159) — primary circuit loss for boilers and
# heat pumps connected to a hot water cylinder via insulated or
# uninsulated primary pipework. The spec lists the zero-loss
# configurations explicitly (combi boilers, integral-vessel heat pumps,
# CPSUs, thermal stores within 1.5 m insulated pipe, etc.); callers
# must gate this helper on those exemptions.
PIPEWORK_INSULATED_UNINSULATED: Final[float] = 0.0
PIPEWORK_INSULATED_FIRST_METRE: Final[float] = 0.1
PIPEWORK_INSULATED_ALL_ACCESSIBLE: Final[float] = 0.3
PIPEWORK_INSULATED_FULLY: Final[float] = 1.0

# Per Table 3 hours-per-day table: 5 winter / 3 summer if cylinder
# thermostat present and water heating not separately timed; 3 / 3 if
# cylinder thermostat present AND separately timed; 11 / 3 if no
# cylinder thermostat. "Use summer value for June, July, August and
# September and winter value for other months."
_SUMMER_MONTH_INDICES: Final[tuple[int, ...]] = (5, 6, 7, 8)  # Jun..Sep


def primary_circuit_hours_per_day_table_3(
    *,
    has_cylinder_thermostat: bool,
    separately_timed_dhw: bool,
    heat_network: bool = False,
) -> tuple[float, float]:
    """SAP 10.2 Table 3 (PDF p.159-160) — hours of primary circulation
    per day, returned as `(winter_hours, summer_hours)`:
        no thermostat                       → (11, 3)
        thermostat, not separately timed    → ( 5, 3)
        thermostat, separately timed        → ( 3, 3)

    SAP 10.2 Table 3 (PDF p.160) verbatim: "For heat networks apply the
    formula above with p = 1.0 and h = 3 for all months." So a heat-
    network main uses h=3 winter and summer regardless of the cylinder-
    thermostat / separate-timing lodgement.
    """
    if heat_network:
        return (3.0, 3.0)
    if not has_cylinder_thermostat:
        return (11.0, 3.0)
    if separately_timed_dhw:
        return (3.0, 3.0)
    return (5.0, 3.0)


def primary_loss_monthly_kwh(
    *,
    pipework_insulation_fraction: float,
    has_cylinder_thermostat: bool,
    separately_timed_dhw: bool,
    heat_network: bool = False,
) -> tuple[float, ...]:
    """SAP 10.2 §4 line (59)m via Table 3 (PDF p.159):
        (59)m = n_m × 14 × [{0.0091 × p + 0.0245 × (1 − p)} × h + 0.0263]
    where p is the fraction of primary pipework insulated and h is the
    hours of primary circulation per day (winter / summer split per
    `primary_circuit_hours_per_day_table_3`).

    `heat_network=True` selects the SAP 10.2 Table 3 (PDF p.160) heat-
    network row — h=3 for all months — regardless of the cylinder-
    thermostat / separate-timing lodgement.

    Returns 12 monthly values in calendar order Jan..Dec. Callers must
    gate this helper on the spec's zero-loss configurations
    (combi boilers, integral-vessel HPs, CPSUs, thermal stores ≤ 1.5 m
    insulated pipe, etc.) — the formula assumes the configuration
    incurs the loss.
    """
    p = pipework_insulation_fraction
    pipework_term = 0.0091 * p + 0.0245 * (1.0 - p)
    winter_h, summer_h = primary_circuit_hours_per_day_table_3(
        has_cylinder_thermostat=has_cylinder_thermostat,
        separately_timed_dhw=separately_timed_dhw,
        heat_network=heat_network,
    )
    return tuple(
        n * 14.0 * (
            pipework_term * (summer_h if m in _SUMMER_MONTH_INDICES else winter_h)
            + 0.0263
        )
        for m, n in enumerate(_DAYS_IN_MONTH)
    )


def cylinder_storage_loss_monthly_kwh(
    *,
    volume_l: float,
    insulation_type: Literal["factory_insulated", "loose_jacket"],
    thickness_mm: float,
    has_cylinder_thermostat: bool,
    separately_timed_dhw: bool,
    declared_loss_kwh_per_day: Optional[float] = None,
) -> tuple[float, ...]:
    """SAP 10.2 §4 line (56)m water storage loss per spec (PDF p.136).

    Two branches, selected by whether the manufacturer's declared loss
    factor is lodged:

    a) declared loss known (`declared_loss_kwh_per_day` set):
        (50) = (48) declared loss (kWh/day) × (49) Table-2b temperature factor
       → `volume_l` / `insulation_type` / `thickness_mm` are unused.
    b) declared loss not known (the default):
        (54) = (47) V × (51) L × (52) VF × (53) TF

        (55) = (50) or (54)
        (56)m = (55) × n_m           (n_m = days in month)

    Returns 12 monthly values in calendar order Jan..Dec. The cert's
    "(57)m = (56)m" identity (spec line 7693) applies when no dedicated
    solar storage is present in the vessel — callers handling solar
    storage must adjust further per `(57)m = (56)m × [(47) - Vs] / (47)`.
    """
    TF = cylinder_temperature_factor_table_2b(
        has_cylinder_thermostat=has_cylinder_thermostat,
        separately_timed_dhw=separately_timed_dhw,
    )
    if declared_loss_kwh_per_day is not None:
        # SAP 10.2 §4 (PDF p.136) branch a) — the lodged manufacturer's
        # declared loss (kWh/day) replaces the Table 2 V×L×VF computation;
        # the Table-2b temperature factor still applies (line (49)→(50)).
        combined_55 = declared_loss_kwh_per_day * TF
    else:
        L = cylinder_storage_loss_factor_table_2(
            insulation_type=insulation_type, thickness_mm=thickness_mm,
        )
        VF = cylinder_volume_factor_table_2a(volume_l)
        combined_55 = volume_l * L * VF * TF
    return tuple(combined_55 * n for n in _DAYS_IN_MONTH)


def total_water_heating_demand_monthly_kwh(
    *,
    energy_content_monthly_kwh: tuple[float, ...],
    distribution_loss_monthly_kwh: tuple[float, ...],
    solar_storage_monthly_kwh: tuple[float, ...],
    primary_loss_monthly_kwh: tuple[float, ...],
    combi_loss_monthly_kwh: tuple[float, ...],
) -> tuple[float, ...]:
    """SAP 10.2 §4 line (62)m via the spec's formula:
        (62)m = 0.85 × (45)m + (46)m + (57)m + (59)m + (61)m

    Note (56)m water-storage loss does NOT appear here — it's accounted
    for through the storage system's efficiency in the (63a-d) inputs and
    the (64)m output line. (46)m distribution loss appears in both (62)m
    and (65)m heat gains (with weight 0.8); that's intentional per spec.

    All five monthly arrays must be 12-tuples in calendar order Jan..Dec.
    """
    return tuple(
        0.85 * e + d + s + p + c
        for e, d, s, p, c in zip(
            energy_content_monthly_kwh,
            distribution_loss_monthly_kwh,
            solar_storage_monthly_kwh,
            primary_loss_monthly_kwh,
            combi_loss_monthly_kwh,
        )
    )


def output_from_water_heater_monthly_kwh(
    *,
    total_demand_monthly_kwh: tuple[float, ...],
    wwhrs_monthly_kwh: tuple[float, ...],
    pv_diverter_monthly_kwh: tuple[float, ...],
    solar_monthly_kwh: tuple[float, ...],
    fghrs_monthly_kwh: tuple[float, ...],
) -> tuple[float, ...]:
    """SAP 10.2 §4 line (64)m:
        (64)m = max(0, (62)m + (63a)m + (63b)m + (63c)m + (63d)m)

    Output from the water heater after subtracting renewable / heat-
    recovery contributions. The four reduction inputs are entered as
    negative quantities (heat displaced FROM the boiler/cylinder), so
    the formula uses + not -. Spec note "if (64)m < 0 then set to 0"
    floors the result month-by-month — a renewable-heavy system can't
    show negative delivered heat for the warmest months.
    """
    return tuple(
        max(0.0, t + w + pv + s + f)
        for t, w, pv, s, f in zip(
            total_demand_monthly_kwh,
            wwhrs_monthly_kwh,
            pv_diverter_monthly_kwh,
            solar_monthly_kwh,
            fghrs_monthly_kwh,
        )
    )


def heat_gains_from_water_heating_monthly_kwh(
    *,
    energy_content_monthly_kwh: tuple[float, ...],
    distribution_loss_monthly_kwh: tuple[float, ...],
    solar_storage_monthly_kwh: tuple[float, ...],
    primary_loss_monthly_kwh: tuple[float, ...],
    combi_loss_monthly_kwh: tuple[float, ...],
    electric_shower_monthly_kwh: tuple[float, ...],
) -> tuple[float, ...]:
    """SAP 10.2 §4 line (65)m heat gains released into the heated space:
        (65)m = 0.25 × [0.85 × (45)m + (61)m + (64a)m]
              + 0.80 × [(46)m + (57)m + (59)m]

    First bracket: delivered-heat losses (hot water at the tap + combi
    cycling losses + electric-shower waste heat) at 25 % recovery into
    the dwelling. Second bracket: pipe-side losses (distribution +
    solar storage + primary circuit) at 80 % recovery — these run
    through the heated envelope so most of the loss heats it.

    Per spec footnote at xlsx row 302, include (57)m only if the hot
    water store is in the dwelling. Callers should pass zero for
    out-of-dwelling stores (e.g. communal heat networks).
    """
    return tuple(
        0.25 * (0.85 * e + c + es) + 0.80 * (d + s + p)
        for e, d, s, p, c, es in zip(
            energy_content_monthly_kwh,
            distribution_loss_monthly_kwh,
            solar_storage_monthly_kwh,
            primary_loss_monthly_kwh,
            combi_loss_monthly_kwh,
            electric_shower_monthly_kwh,
        )
    )


def _days_weighted_average(monthly: tuple[float, ...]) -> float:
    """Σ value[m] × n_m / 365 — used by Appendix J equations J4 and J9."""
    return sum(v * d for v, d in zip(monthly, _DAYS_IN_MONTH)) / _DAYS_IN_YEAR


def annual_average_hot_water_other_uses_l_per_day(
    *, n_occupants: float, low_water_use: bool
) -> float:
    """SAP 10.2 §4 — V_d,other,ave per Appendix J step 3a: 9.8 × N + 14,
    less 5% for the low-water-use target. The annual average is computed
    BEFORE Table J2 monthly modulation, so it's the unmodulated baseline
    value (the days-weighted mean of (42c)m would drift slightly off
    because Table J2 doesn't days-average to exactly 1)."""
    annual = 9.8 * n_occupants + 14.0
    if low_water_use:
        annual *= 1.0 - _LOW_WATER_USE_REDUCTION
    return annual


def annual_average_hot_water_l_per_day(
    *,
    showers_monthly: tuple[float, ...],
    baths_monthly: tuple[float, ...],
    other_uses_annual_avg: float,
) -> float:
    """SAP 10.2 §4 line (43) via Appendix J equation J12:
        V_d,ave = V_d,shower,ave + V_d,bath,ave + V_d,other,ave

    Per the spec text after J12, (43) is the sum of the three component
    annual averages — NOT the days-weighted average of (44)m. The
    distinction only matters for "other uses": its monthly array (42c)m
    is the unmodulated annual baseline times Table J2 factors, and the
    days-weighted average of those factors is 0.9997 (not exactly 1.0),
    so taking the days-weighted mean of (42c)m would drift slightly low.
    Showers and baths only have annual averages via J4 / J9.
    """
    return (
        _days_weighted_average(showers_monthly)
        + _days_weighted_average(baths_monthly)
        + other_uses_annual_avg
    )


def hot_water_other_uses_monthly_l_per_day(
    *, n_occupants: float, low_water_use: bool
) -> tuple[float, ...]:
    """SAP 10.2 §4 line (42c)m via Appendix J equation J11.

    Annual-average daily hot water use for "other purposes" (i.e. not
    showers, not baths — sinks, dishwashers, etc.):

        V_d,other,ave = 9.8 × N + 14

    reduced by 5% if `low_water_use` is True (dwelling designed for ≤125
    L/person/day total water use). The monthly array applies Table J2's
    factor sequence so each entry is daily L for that month.
    """
    annual_average = 9.8 * n_occupants + 14.0
    if low_water_use:
        annual_average *= 1.0 - _LOW_WATER_USE_REDUCTION
    return tuple(annual_average * f for f in _TABLE_J2_MONTHLY_FACTORS)


_ELECTRIC_SHOWER_DEFAULT_KW: Final[float] = 9.3
_ELECTRIC_SHOWER_DURATION_HOURS: Final[float] = 0.1


def electric_shower_monthly_kwh(
    *,
    n_occupants: float,
    has_bath: bool,
    n_outlets: int,
    n_electric_showers: int,
    rated_power_kw: float = _ELECTRIC_SHOWER_DEFAULT_KW,
) -> tuple[float, ...]:
    """SAP 10.2 §4 line (64a)m via Appendix J (p.82) step 8.

    Per outlet for each month:
        N_ES = N_shower / N_outlets    (eq J16)
        E_ES,j,m = N_ES × f_beh × P_ES,j × 0.1 × n_m  (eq J17)
    Summed across electric-shower outlets j (eq J18).

    N_shower from step 1c (same branch as `hot_water_mixer_showers_monthly_
    l_per_day`); N_outlets is the cert-lodged total of mixer + electric
    outlets. P_ES,j defaults to Table J4 row "Instantaneous electric
    shower" = 9.3 kW for assessments of existing dwellings.

    Returns 12-tuple of zeros when there are no electric showers."""
    if n_electric_showers <= 0 or n_outlets <= 0:
        return tuple(0.0 for _ in range(12))
    if has_bath:
        n_shower = 0.45 * n_occupants + 0.65
    else:
        n_shower = 0.58 * n_occupants + 0.83
    n_es_per_outlet = n_shower / n_outlets
    return tuple(
        n_electric_showers
        * n_es_per_outlet
        * fbeh
        * rated_power_kw
        * _ELECTRIC_SHOWER_DURATION_HOURS
        * n_m
        for fbeh, n_m in zip(TABLE_J5_BEHAVIOURAL_FACTOR, _DAYS_IN_MONTH)
    )


def water_heating_from_cert(
    *,
    epc: EpcPropertyData,
    mixer_shower_flow_rates_l_per_min: tuple[float, ...],
    has_bath: bool,
    cold_water_temps_c: tuple[float, ...],
    low_water_use: bool,
    combi_loss_monthly_kwh_override: Optional[tuple[float, ...]] = None,
    solar_storage_monthly_kwh_override: Optional[tuple[float, ...]] = None,
    primary_loss_monthly_kwh_override: Optional[tuple[float, ...]] = None,
    solar_water_heating_monthly_kwh_override: Optional[tuple[float, ...]] = None,
    electric_shower_monthly_kwh_override: Optional[tuple[float, ...]] = None,
    has_electric_shower: bool = False,
    electric_shower_count: int = 0,
    is_instantaneous_at_point_of_use: bool = False,
) -> WaterHeatingResult:
    """SAP 10.2 §4 orchestrator — chain every line ref from (42) through
    (65) for a combi-gas dwelling with optional PCDB-backed combi loss.

    Inputs the cert / site notes contribute:
      - TFA            → occupancy (line 42)
      - bath presence  → bath formula branch (J6)
      - shower flow rates per mixer outlet → (42a)m
      - cold water source (mains / header tank) → Tcold table
      - low-water-use target flag → J7/J11 -5% reduction

    `combi_loss_monthly_kwh_override` lets callers inject a (61)m array
    derived from PCDB Table 3b/3c (tested boilers). When omitted the
    cascade defaults to Table 3a row "Instantaneous, with keep-hot
    facility controlled by time clock" — the modal lodging for non-PCDB
    combis.

    `electric_shower_monthly_kwh_override` lets callers inject (64a)m for
    dwellings with an instantaneous electric shower. The orchestrator
    doesn't yet derive (64a)m from the cert (cert shower-type code →
    Appendix J electric-shower equation), so callers with a worksheet to
    hand can pass the value through. Affects (65)m heat gains via the
    25%-recovery first bracket.

    All remaining (47)–(60), (63a-d) branches default to zero — suits
    the combi-no-storage-no-solar-no-renewables population. Cylinder +
    solar + WWHRS / PV-diverter / FGHRS paths land in future slices.
    """
    if epc.total_floor_area_m2 is None:
        raise ValueError("EpcPropertyData.total_floor_area_m2 is required for §4")
    n = assumed_occupancy(epc.total_floor_area_m2)
    showers = hot_water_mixer_showers_monthly_l_per_day(
        n_occupants=n,
        has_bath=has_bath,
        mixer_shower_flow_rates_l_per_min=mixer_shower_flow_rates_l_per_min,
        cold_water_temps_c=cold_water_temps_c,
        # SAP 10.2 Appendix J step 1a — Noutlets includes electric showers.
        n_electric_showers=electric_shower_count if has_electric_shower else 0,
    )
    has_shower = (
        len(mixer_shower_flow_rates_l_per_min) > 0
        or electric_shower_monthly_kwh_override is not None
        or has_electric_shower
    )
    baths = hot_water_baths_monthly_l_per_day(
        n_occupants=n,
        has_bath=has_bath,
        has_shower=has_shower,
        cold_water_temps_c=cold_water_temps_c,
        low_water_use=low_water_use,
    )
    other = hot_water_other_uses_monthly_l_per_day(
        n_occupants=n, low_water_use=low_water_use,
    )
    daily_total = total_hot_water_monthly_l_per_day(
        showers=showers, baths=baths, other_uses=other,
    )
    other_avg = annual_average_hot_water_other_uses_l_per_day(
        n_occupants=n, low_water_use=low_water_use,
    )
    annual_avg = annual_average_hot_water_l_per_day(
        showers_monthly=showers,
        baths_monthly=baths,
        other_uses_annual_avg=other_avg,
    )
    energy_content = energy_content_of_hot_water_monthly_kwh(
        monthly_hot_water_l_per_day=daily_total,
        cold_water_temps_c=cold_water_temps_c,
    )
    distribution = distribution_loss_monthly_kwh(
        monthly_energy_content_kwh=energy_content,
        is_instantaneous_at_point_of_use=is_instantaneous_at_point_of_use,
    )
    combi = (
        combi_loss_monthly_kwh_override
        if combi_loss_monthly_kwh_override is not None
        else combi_loss_monthly_kwh_table_3a_row_1_no_keep_hot(
            daily_hot_water_monthly_l_per_day=daily_total,
        )
    )
    zero12 = (0.0,) * 12
    solar_storage = (
        solar_storage_monthly_kwh_override
        if solar_storage_monthly_kwh_override is not None
        else zero12
    )
    primary_loss = (
        primary_loss_monthly_kwh_override
        if primary_loss_monthly_kwh_override is not None
        else zero12
    )
    total_demand = total_water_heating_demand_monthly_kwh(
        energy_content_monthly_kwh=energy_content,
        distribution_loss_monthly_kwh=distribution,
        solar_storage_monthly_kwh=solar_storage,
        primary_loss_monthly_kwh=primary_loss,
        combi_loss_monthly_kwh=combi,
    )
    solar_hw = (
        solar_water_heating_monthly_kwh_override
        if solar_water_heating_monthly_kwh_override is not None
        else zero12
    )
    output = output_from_water_heater_monthly_kwh(
        total_demand_monthly_kwh=total_demand,
        wwhrs_monthly_kwh=zero12,
        pv_diverter_monthly_kwh=zero12,
        solar_monthly_kwh=solar_hw,
        fghrs_monthly_kwh=zero12,
    )
    if electric_shower_monthly_kwh_override is not None:
        electric_shower = electric_shower_monthly_kwh_override
    elif electric_shower_count > 0:
        # Appendix J step 8 — N_outlets counts mixer + electric outlets
        # together (eq J16).
        n_outlets_total = len(mixer_shower_flow_rates_l_per_min) + electric_shower_count
        electric_shower = electric_shower_monthly_kwh(
            n_occupants=n,
            has_bath=has_bath,
            n_outlets=n_outlets_total,
            n_electric_showers=electric_shower_count,
        )
    else:
        electric_shower = zero12
    gains = heat_gains_from_water_heating_monthly_kwh(
        energy_content_monthly_kwh=energy_content,
        distribution_loss_monthly_kwh=distribution,
        solar_storage_monthly_kwh=solar_storage,
        primary_loss_monthly_kwh=primary_loss,
        combi_loss_monthly_kwh=combi,
        electric_shower_monthly_kwh=electric_shower,
    )
    return WaterHeatingResult(
        occupancy=n,
        annual_avg_hot_water_l_per_day=annual_avg,
        daily_hot_water_l_per_day_monthly=daily_total,
        energy_content_monthly_kwh=energy_content,
        distribution_loss_monthly_kwh=distribution,
        solar_storage_monthly_kwh=solar_storage,
        primary_loss_monthly_kwh=primary_loss,
        combi_loss_monthly_kwh=combi,
        total_demand_monthly_kwh=total_demand,
        output_monthly_kwh=output,
        heat_gains_monthly_kwh=gains,
        electric_shower_monthly_kwh=electric_shower,
        output_kwh_per_yr=sum(output),
        electric_shower_kwh_per_yr=sum(electric_shower),
    )