Model/domain/sap10_calculator/worksheet/solar_gains.py

"""SAP 10.2 §6 + Appendix U §U3.2 — solar gains.

Two layers:

1. `surface_solar_flux_w_per_m2(orientation, pitch_deg, region, month)`
   implements the §U3.2 polynomial that converts the horizontal solar
   irradiance from Table U3 into a per-orientation per-pitch surface flux.
   Reads:
   - S_h,m from Appendix U Table U3 (already in `domain.sap10_calculator.climate.appendix_u`)
   - δ from Appendix U Table U3 footer (already in `appendix_u.solar_declination_deg`)
   - φ from Appendix U Table U4 (this module)
   - k1..k9 from Appendix U Table U5 (this module)
   - Formula:  S(orient, p, m) = S_h,m × R_h-inc
               R_h-inc = A cos²(φ-δ) + B cos(φ-δ) + C
               where A, B, C are cubics in sin(p/2) with coefficients k1-k9.

2. `window_solar_gain_w(area_m2, surface_flux, g_value, frame_factor,
   overshading)` implements §6.1 equation (5) `G = 0.9 × A × S × g⊥ × FF × Z`,
   plus an `_g_window` alternative for BFRC-rated windows (equation (6)).

Defaults for g⊥ (Table 6b), FF (Table 6c), Z (Table 6d) are deferred to
the cert→inputs mapper slice — this module takes them as caller inputs so
its physics is independent of cert-shape assumptions.

Reference: SAP 10.2 specification §6 + Appendix U §§U3.2-3.3 (pages
127-129); Table U5 columns map 8 cardinal cert orientations to 5
coefficient sets (N, NE/NW, E/W, SE/SW, S).
"""

from __future__ import annotations

from dataclasses import dataclass
from decimal import Decimal, ROUND_HALF_UP
from enum import Enum
from math import cos, radians, sin
from typing import Final

from datatypes.epc.domain.epc_property_data import EpcPropertyData, SapWindow
from domain.sap10_calculator.worksheet.conservatory import (
    CONSERVATORY_ROOF_PITCH_DEG,
    conservatory_geometry,
)
from domain.sap10_calculator.tables.pcdb.postcode_weather import PostcodeClimate
from domain.sap10_calculator.climate.appendix_u import (
    horizontal_solar_irradiance_w_per_m2,
    solar_declination_deg,
)
from domain.sap10_calculator.worksheet.internal_gains import OvershadingCategory


def _decimal_window_area_2dp(width: float, height: float) -> float:
    """W × H in Decimal arithmetic, HALF_UP-quantised at 2 d.p. Mirrors
    `_round_area_2dp(W × H)` but lands on the exact .005 spec boundary
    that float multiplication drops (e.g. 0.65 × 0.70 = 0.4550 exact /
    0.45499... in float — float rounding snaps to 0.45 vs spec 0.46).
    Matches `heat_transmission._decimal_round_half_up_product` so solar
    gains' per-window areas agree with the fabric cascade."""
    d = Decimal(str(width)) * Decimal(str(height))
    return float(d.quantize(Decimal("0.01"), rounding=ROUND_HALF_UP))


# Table 6d first column — winter solar access factor Z for heating gains.
# Distinct from the lighting Z_L (third column, §5) and cooling Z (second
# column, out of scope). SAP 10.2 spec p178.
_Z_SOLAR_BY_OVERSHADING: Final[dict[OvershadingCategory, float]] = {
    OvershadingCategory.HEAVY: 0.3,
    OvershadingCategory.MORE_THAN_AVERAGE: 0.54,
    OvershadingCategory.AVERAGE: 0.77,
    OvershadingCategory.VERY_LITTLE: 1.0,
}


def z_solar_for_overshading(overshading: OvershadingCategory) -> float:
    """SAP 10.2 Table 6d winter solar access factor Z for heating gains."""
    return _Z_SOLAR_BY_OVERSHADING[overshading]


class Orientation(Enum):
    """Cardinal compass directions per SAP10 orientation codes 1-8."""

    N = "N"
    NE = "NE"
    E = "E"
    SE = "SE"
    S = "S"
    SW = "SW"
    W = "W"
    NW = "NW"


# Appendix U Table U4 — representative latitude (°N) for each of 22 SAP
# climate regions (region 0 = UK average).
_LATITUDE_DEG: Final[tuple[float, ...]] = (
    53.5,   # 0 UK average
    51.6, 51.1, 50.9, 50.5, 51.5, 52.6, 53.5, 54.6, 55.2, 54.4,
    53.5, 52.1, 52.6, 55.9, 56.2, 57.3, 57.5, 57.7, 59.0, 60.1, 54.6,
)


# Appendix U Table U5 — k1..k9 polynomial constants per Table-U5 column.
# Columns: North, NE/NW, E/W, SE/SW, South.
_K_NORTH: Final[tuple[float, ...]] = (26.3, -38.5, 14.8, -16.5, 27.3, -11.9, -1.06, 0.0872, -0.191)
_K_NE_NW: Final[tuple[float, ...]] = (0.165, -3.68, 3.0, 6.38, -4.53, -0.405, -4.38, 4.89, -1.99)
_K_E_W: Final[tuple[float, ...]] = (1.44, -2.36, 1.07, -0.514, 1.89, -1.64, -0.542, -0.757, 0.604)
_K_SE_SW: Final[tuple[float, ...]] = (-2.95, 2.89, 1.17, 5.67, -3.54, -4.28, -2.72, -0.25, 3.07)
_K_SOUTH: Final[tuple[float, ...]] = (-0.66, -0.106, 2.93, 3.63, -0.374, -7.4, -2.71, -0.991, 4.59)


_ORIENTATION_TO_K: Final[dict[Orientation, tuple[float, ...]]] = {
    Orientation.N: _K_NORTH,
    Orientation.NE: _K_NE_NW,
    Orientation.E: _K_E_W,
    Orientation.SE: _K_SE_SW,
    Orientation.S: _K_SOUTH,
    Orientation.SW: _K_SE_SW,
    Orientation.W: _K_E_W,
    Orientation.NW: _K_NE_NW,
}


def _latitude_deg(region_or_climate: "int | PostcodeClimate") -> float:
    if isinstance(region_or_climate, PostcodeClimate):
        return region_or_climate.latitude_deg
    if not 0 <= region_or_climate < len(_LATITUDE_DEG):
        raise ValueError(
            f"region must be 0..{len(_LATITUDE_DEG) - 1}, got {region_or_climate}"
        )
    return _LATITUDE_DEG[region_or_climate]


def surface_solar_flux_w_per_m2(
    *,
    orientation: Orientation,
    pitch_deg: float,
    region: "int | PostcodeClimate",
    month: int,
) -> float:
    """Per-orientation per-pitch monthly solar flux on a surface (W/m²).

    SAP 10.2 Appendix U §U3.2 polynomial conversion from the horizontal
    irradiance in Table U3 to any orientation/tilt combination. Accepts
    either a SAP region index (0..21) or a `PostcodeClimate` record from
    PCDB Table 172 (demand cascade).
    """
    s_h = horizontal_solar_irradiance_w_per_m2(region, month)
    declination = solar_declination_deg(month)
    latitude = _latitude_deg(region)

    half_pitch = radians(pitch_deg) / 2.0
    s = sin(half_pitch)
    s2 = s * s
    s3 = s2 * s

    k1, k2, k3, k4, k5, k6, k7, k8, k9 = _ORIENTATION_TO_K[orientation]
    a = k1 * s3 + k2 * s2 + k3 * s
    b = k4 * s3 + k5 * s2 + k6 * s
    c = k7 * s3 + k8 * s2 + k9 * s + 1.0

    cos_phi_minus_delta = cos(radians(latitude - declination))
    r_h_inc = a * cos_phi_minus_delta * cos_phi_minus_delta + b * cos_phi_minus_delta + c
    return s_h * r_h_inc


def window_solar_gain_w(
    *,
    area_m2: float,
    surface_flux_w_per_m2: float,
    g_perpendicular: float,
    frame_factor: float,
    overshading_factor: float,
) -> float:
    """Solar gain through a window (W). SAP 10.2 §6.1 equation (5):

        G_solar = 0.9 × A_w × S × g⊥ × FF × Z

    where 0.9 corrects for the ratio of typical average transmittance to
    that at normal incidence. For BFRC-rated windows whose quoted g_window
    already bakes in 0.9 × g⊥ × FF, use `window_solar_gain_w_bfrc` instead.
    """
    return 0.9 * area_m2 * surface_flux_w_per_m2 * g_perpendicular * frame_factor * overshading_factor


# SAP 10.2 Table 6b — total solar energy transmittance g⊥ at normal incidence
# by SAP10 glazing_type code (p178). Default 0.76 (modal double-glazed UK
# stock) when the cert's glazing_type is missing or unrecognised.
#
# Codes 1-6 follow the SAP 10.2 Table 6b ordering. Codes 7-15 are RdSAP-21
# schema additions (per datatypes/epc/domain/epc_codes.csv); triple-glazed
# cohort certs lodge code 14 (triple 2022+) which previously fell through
# to the 0.76 DG default and over-counted solar gains.
_G_PERPENDICULAR_BY_GLAZING_TYPE: Final[dict[int, float]] = {
    1: 0.85,  # single glazed
    2: 0.76,  # double glazed (air or argon filled) — 2002-2022
    3: 0.76,  # double glazed (air or argon filled) — pre-2002
    4: 0.63,  # double glazed low-E soft-coat
    5: 0.76,  # window with secondary glazing
    6: 0.68,  # triple glazed (air or argon filled)
    # RdSAP 21 schema extensions — triple-glazed variants all share the
    # same Table 6b g⊥; secondary glazing mirrors the existing code 5 value.
    7: 0.76,   # double glazing, known data
    8: 0.68,   # triple glazing, known data
    9: 0.68,   # triple glazing, installed 2002-2022
    10: 0.68,  # triple glazing, installed pre-2002
    11: 0.76,  # secondary glazing, normal emissivity
    12: 0.76,  # secondary glazing, low emissivity
    13: 0.76,  # double glazing, installed 2022+
    14: 0.68,  # triple glazing, installed 2022+
    15: 0.85,  # single glazing, known data
}
_G_PERPENDICULAR_DEFAULT: Final[float] = 0.76


# SAP 10.2 Table 6c — frame factor (proportion of opening glazed) by frame
# material substring. Case-insensitive matching; cert lodges this as free
# text ("PVC", "Wood", "Metal").
_FRAME_FACTOR_BY_MATERIAL_SUBSTR: Final[tuple[tuple[str, float], ...]] = (
    ("metal", 0.8),
    ("aluminium", 0.8),
    ("wood", 0.7),
    ("timber", 0.7),
    ("pvc", 0.7),
    ("upvc", 0.7),
    ("composite", 0.7),
)
_FRAME_FACTOR_DEFAULT: Final[float] = 0.7


# SAP10 octant code → Orientation enum. Cert windows with a code outside 1..8
# (e.g. 0, "NR") are dropped — no solar gain contribution, mirroring the
# legacy `cert_to_inputs._window_inputs` shortcut.
ORIENTATION_BY_SAP10_CODE: Final[dict[int, Orientation]] = {
    1: Orientation.N,
    2: Orientation.NE,
    3: Orientation.E,
    4: Orientation.SE,
    5: Orientation.S,
    6: Orientation.SW,
    7: Orientation.W,
    8: Orientation.NW,
}


_ROOF_WINDOW_DEFAULT_PITCH_DEG: Final[float] = 45.0  # RdSAP10 Table 24
_ROOFLIGHT_PITCH_DEG: Final[float] = 0.0  # SAP10.2 §U3.2 p128
_HORIZONTAL_Z: Final[float] = 1.0  # Table 6d note 2 — roof windows + rooflights
_MONTHS: Final[range] = range(1, 13)


@dataclass(frozen=True)
class RoofWindowInput:
    """SAP 10.2 §6 "Roof window" — glazed opening in a pitched roof. Pitch
    < 70° is treated at its actual inclination; pitch ≥ 70° collapses to
    vertical per §U3.2. Z=1.0 regardless of overshading (Table 6d note 2).
    """

    area_m2: float
    orientation: Orientation
    g_perpendicular: float
    frame_factor: float
    pitch_deg: float = _ROOF_WINDOW_DEFAULT_PITCH_DEG


@dataclass(frozen=True)
class RooflightInput:
    """SAP 10.2 §6 "Rooflight" — horizontal glazed opening. Per §U3.2 p128
    rooflights are assumed horizontal (pitch = 0°); orientation is therefore
    immaterial (the §U3.2 polynomial degenerates to S_h,m at pitch 0).
    Z=1.0 regardless of overshading (Table 6d note 2).
    """

    area_m2: float
    g_perpendicular: float
    frame_factor: float


@dataclass(frozen=True)
class SolarGainsResult:
    """SAP 10.2 §6 line refs (74)..(83), each a 12-tuple of watts per month.

    Returned by `solar_gains_from_cert`. Downstream §7 utilisation factor
    + §9 heat-balance consume `total_solar_gains_monthly_w` directly; the
    per-orientation tuples are exposed for worksheet conformance + audit.
    Field names mirror the SAP 10.2 line refs.
    """

    north_monthly_w: tuple[float, ...]               # (74)
    northeast_monthly_w: tuple[float, ...]           # (75)
    east_monthly_w: tuple[float, ...]                # (76)
    southeast_monthly_w: tuple[float, ...]           # (77)
    south_monthly_w: tuple[float, ...]               # (78)
    southwest_monthly_w: tuple[float, ...]           # (79)
    west_monthly_w: tuple[float, ...]                # (80)
    northwest_monthly_w: tuple[float, ...]           # (81)
    roof_windows_monthly_w: tuple[float, ...]        # (82)
    rooflights_monthly_w: tuple[float, ...]          # (82a)
    total_solar_gains_monthly_w: tuple[float, ...]   # (83)


def _g_perpendicular(w: SapWindow) -> float:
    """Table 6b g⊥. Prefer the cert's measured value
    (`window_transmission_details.solar_transmittance`) when present —
    Elmhurst lodges it for every window; Table 6b is the cascade fallback.
    """
    if w.window_transmission_details is not None:
        return float(w.window_transmission_details.solar_transmittance)
    if isinstance(w.glazing_type, int) and w.glazing_type in _G_PERPENDICULAR_BY_GLAZING_TYPE:
        return _G_PERPENDICULAR_BY_GLAZING_TYPE[w.glazing_type]
    return _G_PERPENDICULAR_DEFAULT


def _frame_factor(w: SapWindow) -> float:
    """Table 6c frame factor. Prefer cert's `frame_factor`; else look up
    by `frame_material` substring."""
    if w.frame_factor is not None:
        return float(w.frame_factor)
    material = (w.frame_material or "").lower()
    for needle, ff in _FRAME_FACTOR_BY_MATERIAL_SUBSTR:
        if needle in material:
            return ff
    return _FRAME_FACTOR_DEFAULT


def _orientation(w: SapWindow) -> Orientation | None:
    """Map cert `orientation` code (1..8) to enum; None for unmapped."""
    if isinstance(w.orientation, int) and w.orientation in ORIENTATION_BY_SAP10_CODE:
        return ORIENTATION_BY_SAP10_CODE[w.orientation]
    return None


def _vertical_window_gain_monthly_w(
    *,
    w: SapWindow,
    orientation: Orientation,
    region: "int | PostcodeClimate",
    z_solar: float,
) -> tuple[float, ...]:
    """Compute the 12-tuple of monthly solar gain (W) for one vertical wall
    window. Pitch = 90° always; Table 6b/6c lookups derive g⊥ and FF."""
    # RdSAP 10 §15 "Rounding of data" (p.66): "All element areas (gross)
    # including window areas: 2 d.p." — matches heat_transmission's per-
    # window area rounding so solar gains and conduction agree on area.
    area = _decimal_window_area_2dp(float(w.window_width), float(w.window_height))
    g_perp = _g_perpendicular(w)
    ff = _frame_factor(w)
    return tuple(
        window_solar_gain_w(
            area_m2=area,
            surface_flux_w_per_m2=surface_solar_flux_w_per_m2(
                orientation=orientation, pitch_deg=90.0, region=region, month=m,
            ),
            g_perpendicular=g_perp,
            frame_factor=ff,
            overshading_factor=z_solar,
        )
        for m in _MONTHS
    )


def _sum_tuples(*tuples: tuple[float, ...]) -> tuple[float, ...]:
    """Element-wise sum of 12-tuples. Returns 12-tuple of zeros if empty."""
    if not tuples:
        return (0.0,) * 12
    return tuple(sum(t[i] for t in tuples) for i in range(12))


def solar_gains_from_cert(
    *,
    epc: EpcPropertyData,
    region: "int | PostcodeClimate",
    overshading: OvershadingCategory = OvershadingCategory.AVERAGE,
    roof_windows: tuple[RoofWindowInput, ...] = (),
    rooflights: tuple[RooflightInput, ...] = (),
) -> SolarGainsResult:
    """SAP 10.2 §6 orchestrator — chain (74)..(83) for the dwelling.

    Inputs:
      epc           cert (sap_windows feed (74)..(81); rooflights NOT lodged)
      region        SAP climate region (0 = UK avg for SAP rating pass)
      overshading   Table 6d bucket → Z for vertical wall windows
      roof_windows  RdSAP §6 "Roof window" — pitched, Z=1.0
      rooflights    RdSAP §6 "Rooflight" — horizontal, Z=1.0

    Coverage caveat: cert summaries do not lodge a distinct rooflight code,
    so `cert_to_inputs` passes both `roof_windows` and `rooflights` as empty
    tuples. Conformance against fixtures with roof glazing is exercised via
    explicit `_build_section_6_epc(fixture)` test wrappers.
    """
    z_vertical = z_solar_for_overshading(overshading)

    per_orientation: dict[Orientation, list[tuple[float, ...]]] = {
        o: [] for o in Orientation
    }
    for w in epc.sap_windows:
        orientation = _orientation(w)
        if orientation is None:
            continue
        per_orientation[orientation].append(
            _vertical_window_gain_monthly_w(
                w=w, orientation=orientation, region=region, z_solar=z_vertical,
            )
        )

    roof_windows_monthly = _sum_tuples(*(
        tuple(
            window_solar_gain_w(
                area_m2=rw.area_m2,
                surface_flux_w_per_m2=surface_solar_flux_w_per_m2(
                    orientation=rw.orientation, pitch_deg=rw.pitch_deg,
                    region=region, month=m,
                ),
                g_perpendicular=rw.g_perpendicular,
                frame_factor=rw.frame_factor,
                overshading_factor=_HORIZONTAL_Z,
            )
            for m in _MONTHS
        )
        for rw in roof_windows
    ))

    rooflights_monthly = _sum_tuples(*(
        tuple(
            window_solar_gain_w(
                area_m2=rl.area_m2,
                surface_flux_w_per_m2=surface_solar_flux_w_per_m2(
                    orientation=Orientation.N,  # immaterial at pitch 0
                    pitch_deg=_ROOFLIGHT_PITCH_DEG,
                    region=region, month=m,
                ),
                g_perpendicular=rl.g_perpendicular,
                frame_factor=rl.frame_factor,
                overshading_factor=_HORIZONTAL_Z,
            )
            for m in _MONTHS
        )
        for rl in rooflights
    ))

    per_orientation_summed = {
        o: _sum_tuples(*per_orientation[o]) for o in Orientation
    }

    # RdSAP 10 §6.1 (PDF p.49) + Table 25 note (p.51): "The orientation of
    # windows in a conservatory is not recorded, thus solar gains are
    # calculated using the default solar flux (East/West orientation, with
    # 20° pitch for roof windows)." Average overshading per §7 (Table 6d).
    # The glazed wall bills onto the (76) East line (vertical, Z=z_vertical);
    # the glazed roof onto the (82) roof-window line (20° pitch, Z=1.0).
    cons = conservatory_geometry(epc)
    if cons is not None:
        cons_wall_monthly = tuple(
            window_solar_gain_w(
                area_m2=cons.glazed_wall_area_m2,
                surface_flux_w_per_m2=surface_solar_flux_w_per_m2(
                    orientation=Orientation.E, pitch_deg=90.0,
                    region=region, month=m,
                ),
                g_perpendicular=cons.g_value,
                frame_factor=cons.frame_factor,
                overshading_factor=z_vertical,
            )
            for m in _MONTHS
        )
        cons_roof_monthly = tuple(
            window_solar_gain_w(
                area_m2=cons.glazed_roof_area_m2,
                surface_flux_w_per_m2=surface_solar_flux_w_per_m2(
                    orientation=Orientation.E,
                    pitch_deg=CONSERVATORY_ROOF_PITCH_DEG,
                    region=region, month=m,
                ),
                g_perpendicular=cons.g_value,
                frame_factor=cons.frame_factor,
                overshading_factor=_HORIZONTAL_Z,
            )
            for m in _MONTHS
        )
        per_orientation_summed[Orientation.E] = _sum_tuples(
            per_orientation_summed[Orientation.E], cons_wall_monthly,
        )
        roof_windows_monthly = _sum_tuples(
            roof_windows_monthly, cons_roof_monthly,
        )

    total = _sum_tuples(
        *per_orientation_summed.values(),
        roof_windows_monthly,
        rooflights_monthly,
    )

    return SolarGainsResult(
        north_monthly_w=per_orientation_summed[Orientation.N],
        northeast_monthly_w=per_orientation_summed[Orientation.NE],
        east_monthly_w=per_orientation_summed[Orientation.E],
        southeast_monthly_w=per_orientation_summed[Orientation.SE],
        south_monthly_w=per_orientation_summed[Orientation.S],
        southwest_monthly_w=per_orientation_summed[Orientation.SW],
        west_monthly_w=per_orientation_summed[Orientation.W],
        northwest_monthly_w=per_orientation_summed[Orientation.NW],
        roof_windows_monthly_w=roof_windows_monthly,
        rooflights_monthly_w=rooflights_monthly,
        total_solar_gains_monthly_w=total,
    )