src/shading/solar_geometry.py

"""
ShadowModel: compute solar position, tracker tilt, and shadow projection
from a solar panel onto a vine canopy grid.

Coordinate system:
- Row-local: u = along-row direction (azimuth 315°), v = cross-row (perpendicular), z = up
- World: x = East, y = North, z = up
"""

from __future__ import annotations

import numpy as np
import pandas as pd
import pvlib

from config.settings import (
    CANOPY_HEIGHT,
    CANOPY_WIDTH,
    PANEL_HEIGHT,
    PANEL_WIDTH,
    ROW_AZIMUTH,
    ROW_SPACING,
    SITE_ALTITUDE,
    SITE_LATITUDE,
    SITE_LONGITUDE,
    TRACKER_GCR,
    TRACKER_MAX_ANGLE,
)


class ShadowModel:
    """3D shadow projection from a single-axis tracker onto vine canopy rows."""

    def __init__(
        self,
        panel_width: float = PANEL_WIDTH,
        panel_height: float = PANEL_HEIGHT,
        row_spacing: float = ROW_SPACING,
        canopy_height: float = CANOPY_HEIGHT,
        canopy_width: float = CANOPY_WIDTH,
        row_azimuth: float = ROW_AZIMUTH,
        lat: float = SITE_LATITUDE,
        lon: float = SITE_LONGITUDE,
        altitude: float = SITE_ALTITUDE,
        n_vertical: int = 3,
        n_horizontal: int = 10,
        n_rows: int = 5,
    ):
        self.panel_width = panel_width
        self.panel_height = panel_height
        self.row_spacing = row_spacing
        self.canopy_height = canopy_height
        self.canopy_width = canopy_width
        self.row_azimuth = row_azimuth
        self.lat = lat
        self.lon = lon
        self.altitude = altitude
        self.n_vertical = n_vertical
        self.n_horizontal = n_horizontal
        self.n_rows = n_rows

        # Row direction unit vectors in world coords (x=East, y=North)
        az_rad = np.radians(row_azimuth)
        self._row_u = np.array([np.sin(az_rad), np.cos(az_rad)])   # along row
        self._row_v = np.array([np.cos(az_rad), -np.sin(az_rad)])  # cross-row (perpendicular, right-hand)

        # Canopy grid in cross-row direction
        self._grid_v = np.linspace(
            -canopy_width / 2, canopy_width / 2, n_horizontal,
        )
        self._grid_z = np.linspace(
            canopy_height / n_vertical / 2,
            canopy_height - canopy_height / n_vertical / 2,
            n_vertical,
        )
        # LAI weight per vertical zone: top 50%, mid 35%, bottom 15%
        self.lai_weights = np.array([0.15, 0.35, 0.50])  # bottom to top

    def get_solar_position(self, times: pd.DatetimeIndex) -> pd.DataFrame:
        """Compute solar elevation and azimuth for site location."""
        loc = pvlib.location.Location(self.lat, self.lon, altitude=self.altitude)
        sp = loc.get_solarposition(times)
        return sp[["apparent_elevation", "azimuth"]].rename(
            columns={"apparent_elevation": "solar_elevation", "azimuth": "solar_azimuth"},
        )

    def compute_tracker_tilt(self, solar_azimuth: float, solar_elevation: float) -> dict:
        """
        Single-axis tracker angle using pvlib.tracking.singleaxis.
        Tracker axis is along row_azimuth, maximizing energy (normal tracking).
        Returns dict with tracker_theta, aoi, surface_tilt, surface_azimuth.
        """
        if solar_elevation <= 0:
            return {"tracker_theta": 0.0, "aoi": 90.0, "surface_tilt": 0.0, "surface_azimuth": 0.0}

        zenith = 90.0 - solar_elevation

        # Use positional args for zenith and azimuth so we work with both
        # pvlib < 0.13.1 (apparent_azimuth) and >= 0.13.1 (solar_azimuth)
        result = pvlib.tracking.singleaxis(
            zenith,
            solar_azimuth,
            axis_tilt=0,
            axis_azimuth=self.row_azimuth,
            max_angle=TRACKER_MAX_ANGLE,
            backtrack=True,
            gcr=TRACKER_GCR,
        )

        # pvlib can return arrays; extract scalar for float() and np.isnan()
        def _scalar(x):
            a = np.asarray(x)
            return a.flat[0] if a.size else np.nan

        theta = float(_scalar(result["tracker_theta"]))
        if np.isnan(theta):
            theta = 0.0
        aoi = _scalar(result["aoi"])
        surf_tilt = _scalar(result["surface_tilt"])
        surf_az = _scalar(result["surface_azimuth"])
        return {
            "tracker_theta": theta,
            "aoi": 90.0 if np.isnan(aoi) else float(aoi),
            "surface_tilt": 0.0 if np.isnan(surf_tilt) else float(surf_tilt),
            "surface_azimuth": 0.0 if np.isnan(surf_az) else float(surf_az),
        }

    def _panel_corners_local(self, tracker_tilt: float) -> np.ndarray:
        """
        4 panel corners in row-local frame (v=cross-row, u=along-row, z=up).
        Returns shape (4, 3) as [v, u, z].
        pvlib convention: positive theta = panel faces +v (NE for axis 315°).
        Negate theta so that positive theta lowers the +v edge (panel faces +v).
        """
        half_w = self.panel_width / 2
        tilt_rad = np.radians(-tracker_tilt)  # negate to match pvlib convention
        cos_t, sin_t = np.cos(tilt_rad), np.sin(tilt_rad)
        half_len = 1.0  # panel length along row (for visualization)

        # Corners: [v, u, z] — panel tilts in the v-z plane
        return np.array([
            [-half_w * cos_t, -half_len, self.panel_height - half_w * sin_t],
            [half_w * cos_t, -half_len, self.panel_height + half_w * sin_t],
            [half_w * cos_t, half_len, self.panel_height + half_w * sin_t],
            [-half_w * cos_t, half_len, self.panel_height - half_w * sin_t],
        ])

    def panel_corners_world(self, tracker_tilt: float, row_offset: float = 0.0) -> np.ndarray:
        """
        Panel corners in world coordinates for a row at cross-row offset.
        Returns shape (4, 3) as [x, y, z].
        """
        local = self._panel_corners_local(tracker_tilt)
        world = np.zeros_like(local)
        for i in range(4):
            # v -> cross-row, u -> along-row
            world[i, 0] = local[i, 0] * self._row_v[0] + local[i, 1] * self._row_u[0] + row_offset * self._row_v[0]
            world[i, 1] = local[i, 0] * self._row_v[1] + local[i, 1] * self._row_u[1] + row_offset * self._row_v[1]
            world[i, 2] = local[i, 2]
        return world

    def vine_box_world(self, row_offset: float = 0.0) -> np.ndarray:
        """
        8 vine canopy box corners in world coordinates.
        Returns shape (8, 3) as [x, y, z].
        """
        cw = self.canopy_width / 2
        ch = self.canopy_height
        half_len = 1.5  # length along row for visualization
        # Local corners [v, u, z]
        local = np.array([
            [-cw, -half_len, 0], [cw, -half_len, 0],
            [cw, half_len, 0], [-cw, half_len, 0],
            [-cw, -half_len, ch], [cw, -half_len, ch],
            [cw, half_len, ch], [-cw, half_len, ch],
        ])
        world = np.zeros_like(local)
        for i in range(8):
            world[i, 0] = local[i, 0] * self._row_v[0] + local[i, 1] * self._row_u[0] + row_offset * self._row_v[0]
            world[i, 1] = local[i, 0] * self._row_v[1] + local[i, 1] * self._row_u[1] + row_offset * self._row_v[1]
            world[i, 2] = local[i, 2]
        return world

    def project_shadow(
        self,
        solar_elevation: float,
        solar_azimuth: float,
        tracker_tilt: float | None = None,
    ) -> np.ndarray:
        """
        Combined shadow on a vine row from its own panel AND neighboring panels.
        Returns boolean mask of shape (n_vertical, n_horizontal).
        True = shaded, False = sunlit.
        """
        if solar_elevation <= 2.0:
            return np.ones((self.n_vertical, self.n_horizontal), dtype=bool)

        if tracker_tilt is None:
            tracker_tilt = self.compute_tracker_tilt(solar_azimuth, solar_elevation)["tracker_theta"]

        # Aggregate shadows from own panel + immediate neighbors (east & west)
        mask = np.zeros((self.n_vertical, self.n_horizontal), dtype=bool)
        for panel_offset_idx in (-1, 0, 1):
            source_offset = panel_offset_idx * self.row_spacing
            m = self.project_shadow_on_row(
                solar_elevation, solar_azimuth, tracker_tilt,
                source_row_offset=source_offset,
                target_row_offset=0.0,
            )
            mask |= m
        return mask

    def project_shadow_on_row(
        self,
        solar_elevation: float,
        solar_azimuth: float,
        tracker_tilt: float,
        source_row_offset: float,
        target_row_offset: float,
    ) -> np.ndarray:
        """
        Project shadow from panel at source_row_offset onto canopy at target_row_offset.
        Uses ray-segment intersection: for each canopy zone, trace a ray toward the sun
        and check if it hits the panel. This correctly handles shadow on the top face,
        east face, and west face of the canopy rectangle.
        Returns boolean mask (n_vertical, n_horizontal).
        """
        if solar_elevation <= 2.0:
            return np.zeros((self.n_vertical, self.n_horizontal), dtype=bool)

        elev_rad = np.radians(solar_elevation)
        sun_horiz = np.array([np.cos(elev_rad) * np.sin(np.radians(solar_azimuth)),
                              np.cos(elev_rad) * np.cos(np.radians(solar_azimuth))])
        sun_cross = float(np.dot(sun_horiz, self._row_v))
        sun_z = np.sin(elev_rad)

        if sun_z <= 0.01:
            return np.zeros((self.n_vertical, self.n_horizontal), dtype=bool)

        # Panel segment in cross-section (v-z plane), in target canopy frame
        row_shift = target_row_offset - source_row_offset
        corners_local = self._panel_corners_local(tracker_tilt)
        # Panel lower and upper edges (2 unique v-z points in cross-section)
        p1_v = corners_local[0, 0] - row_shift  # lower edge v
        p1_z = corners_local[0, 2]              # lower edge z
        p2_v = corners_local[1, 0] - row_shift  # upper edge v
        p2_z = corners_local[1, 2]              # upper edge z

        # Panel segment direction
        dv_p = p2_v - p1_v
        dz_p = p2_z - p1_z

        mask = np.zeros((self.n_vertical, self.n_horizontal), dtype=bool)

        # Ray-segment intersection for each canopy zone
        # Ray from (v_c, z_c) toward sun: (v_c + sun_cross*t, z_c + sun_z*t)
        # Panel segment: (p1_v + s*dv_p, p1_z + s*dz_p), s in [0,1]
        # Solve: v_c + sun_cross*t = p1_v + s*dv_p
        #        z_c + sun_z*t     = p1_z + s*dz_p
        det = sun_cross * dz_p - sun_z * dv_p
        if abs(det) < 1e-10:
            return mask  # sun ray parallel to panel

        for iz, z_c in enumerate(self._grid_z):
            for iv, v_c in enumerate(self._grid_v):
                dv = p1_v - v_c
                dz = p1_z - z_c
                t = (dv * dz_p - dz * dv_p) / det
                s = (dv * sun_z - dz * sun_cross) / det
                # t > 0: panel is above/toward the sun (not behind canopy)
                # 0 <= s <= 1: hit within panel width
                mask[iz, iv] = (t > 0) and (0.0 <= s <= 1.0)

        return mask

    def shadow_world_polygon(
        self,
        solar_elevation: float,
        solar_azimuth: float,
        tracker_tilt: float,
        row_offset: float = 0.0,
    ) -> np.ndarray | None:
        """
        Shadow polygon on canopy top in world coords for visualization.
        Returns (4, 3) array or None if sun below horizon.
        """
        if solar_elevation <= 2.0:
            return None

        elev_rad = np.radians(solar_elevation)
        az_rad = np.radians(solar_azimuth)
        sun_x = np.cos(elev_rad) * np.sin(az_rad)
        sun_y = np.cos(elev_rad) * np.cos(az_rad)
        sun_z = np.sin(elev_rad)

        if sun_z <= 0.01:
            return None

        corners = self.panel_corners_world(tracker_tilt, row_offset)
        ch = self.canopy_height
        shadow = np.zeros((4, 3))
        for i in range(4):
            dz = corners[i, 2] - ch
            shadow[i, 0] = corners[i, 0] - (sun_x / sun_z) * dz
            shadow[i, 1] = corners[i, 1] - (sun_y / sun_z) * dz
            shadow[i, 2] = ch + 0.01
        return shadow

    def compute_par_distribution(
        self,
        total_par: float,
        shadow_mask: np.ndarray,
        diffuse_fraction: float = 0.15,
        solar_elevation: float | None = None,
        solar_azimuth: float | None = None,
        tracker_tilt: float | None = None,
    ) -> np.ndarray:
        """PAR per zone accounting for overhead shadow, side light, and diffuse.

        When solar position is provided, zones near the row edge on the
        sun-facing side receive additional side light that enters between
        panel rows.  Vertical variation: top zones get more sky-view diffuse,
        bottom zones get less.

        Falls back to the simple sunlit/shaded model when solar position is
        not provided.
        """
        n_v, n_h = shadow_mask.shape
        par = np.full((n_v, n_h), total_par * diffuse_fraction, dtype=float)

        # --- Overhead direct beam: only where NOT shaded ---
        par[~shadow_mask] = total_par

        # --- Side light from row gaps (requires solar geometry) ---
        if solar_elevation is not None and solar_azimuth is not None and solar_elevation > 2:
            sun_cross, sun_z = self._sun_cross_component(solar_elevation, solar_azimuth)

            # Normalise horizontal positions to [-1, +1] across canopy width
            v_norm = np.linspace(-1.0, 1.0, n_h)

            # Side-light fraction: zones on the sun-facing edge of the canopy
            # receive direct light between panel rows.  The fraction decays
            # linearly from the edge toward the interior.
            #   sun_cross > 0 → sun from +v side → right edge lit
            #   sun_cross < 0 → sun from -v side → left edge lit
            if abs(sun_cross) > 0.01:
                # How far each horizontal position is from the sun-facing edge (0-1)
                if sun_cross > 0:
                    edge_proximity = np.clip((v_norm + 1) / 2, 0, 1)  # 1 at +v edge
                else:
                    edge_proximity = np.clip((1 - v_norm) / 2, 0, 1)  # 1 at -v edge

                # Side light penetration depends on sun elevation:
                # low sun → more horizontal light between rows
                # high sun → less side penetration
                horiz_factor = abs(sun_cross) / (abs(sun_cross) + sun_z)

                # Row-gap openness: fraction of sky visible from the side
                gap_fraction = 1.0 - self.panel_width / self.row_spacing  # ~0.62

                # Side PAR contribution for each horizontal position
                side_par = total_par * horiz_factor * gap_fraction * edge_proximity

                # Apply to shaded zones only (sunlit zones already have full PAR)
                for iz in range(n_v):
                    for ih in range(n_h):
                        if shadow_mask[iz, ih]:
                            par[iz, ih] += side_par[ih]

            # --- Vertical variation in diffuse sky view ---
            # Top zones see more open sky, bottom zones are self-shaded by
            # upper canopy.  Sky-view factor from 0.4 (bottom) to 1.0 (top).
            z_fracs = np.linspace(0.0, 1.0, n_v)
            sky_view = 0.4 + 0.6 * z_fracs  # bottom=0.4, top=1.0
            diffuse_sky = total_par * diffuse_fraction * sky_view
            for iz in range(n_v):
                for ih in range(n_h):
                    if shadow_mask[iz, ih]:
                        # Replace base diffuse with sky-view-weighted diffuse
                        par[iz, ih] = (par[iz, ih]
                                       - total_par * diffuse_fraction
                                       + diffuse_sky[iz])

        return np.clip(par, 0, total_par)

    def compute_face_par_zones(
        self,
        total_par: float,
        solar_elevation: float,
        solar_azimuth: float,
        tracker_tilt: float,
        diffuse_fraction: float = 0.15,
        include_panels: bool = True,
    ) -> dict:
        """Compute PAR reaching each canopy face at zone level.

        Returns dict with:
            east: array(3,) PAR at Bottom/Middle/Top zones on east face
            west: array(3,) PAR at Bottom/Middle/Top zones on west face
            top: array(3,) PAR at 3 positions across top face (W-edge, Centre, E-edge)
        """
        n_z = self.n_vertical  # 3
        ch = self.canopy_height
        cw2 = self.canopy_width / 2

        # Zone centres along canopy height
        zone_z = np.linspace(ch / n_z / 2, ch - ch / n_z / 2, n_z)

        # Diffuse PAR with sky-view factor per height zone
        z_fracs = np.linspace(0.0, 1.0, n_z)
        sky_view = 0.4 + 0.6 * z_fracs  # bottom=0.4, top=1.0
        diffuse_par = total_par * diffuse_fraction * sky_view

        east_par = diffuse_par.copy()
        west_par = diffuse_par.copy()

        if solar_elevation <= 2.0 or total_par <= 0:
            top_par = np.full(3, total_par * diffuse_fraction)
            return {"east": east_par, "west": west_par, "top": top_par}

        sun_cross, sun_z = self._sun_cross_component(solar_elevation, solar_azimuth)
        if sun_z <= 0.01:
            top_par = np.full(3, total_par * diffuse_fraction)
            return {"east": east_par, "west": west_par, "top": top_par}

        # Panel ray-tracing setup (reuses same logic as compute_face_shading)
        if include_panels:
            panels = []
            for off_idx in (-1, 0, 1):
                row_shift = -off_idx * self.row_spacing
                corners = self._panel_corners_local(tracker_tilt)
                panels.append((corners[0, 0] - row_shift, corners[0, 2],
                               corners[1, 0] - row_shift, corners[1, 2]))

            def _shaded(v_c, z_c):
                for p1_v, p1_z, p2_v, p2_z in panels:
                    dv_p, dz_p = p2_v - p1_v, p2_z - p1_z
                    det = sun_cross * dz_p - sun_z * dv_p
                    if abs(det) < 1e-10:
                        continue
                    dv, dz = p1_v - v_c, p1_z - z_c
                    t = (dv * dz_p - dz * dv_p) / det
                    s = (dv * sun_z - dz * sun_cross) / det
                    if t > 0 and 0.0 <= s <= 1.0:
                        return True
                return False
        else:
            def _shaded(v_c, z_c):
                return False

        # Side-light factor for vertical faces
        gap_fraction = 1.0 - self.panel_width / self.row_spacing
        horiz_factor = abs(sun_cross) / (abs(sun_cross) + sun_z)
        side_direct = total_par * horiz_factor * gap_fraction

        # Self-shading: face normal away from sun receives no direct beam
        east_self_shaded = (sun_cross <= 0)
        west_self_shaded = (sun_cross >= 0)

        # East face (v = +cw2, normal = +v)
        for iz in range(n_z):
            if east_self_shaded:
                east_par[iz] = diffuse_par[iz] * 0.3  # back face
            elif not _shaded(cw2, zone_z[iz]):
                east_par[iz] = side_direct + diffuse_par[iz]
            # else: panel-shaded, stays at diffuse_par[iz]

        # West face (v = -cw2, normal = -v)
        for iz in range(n_z):
            if west_self_shaded:
                west_par[iz] = diffuse_par[iz] * 0.3
            elif not _shaded(-cw2, zone_z[iz]):
                west_par[iz] = side_direct + diffuse_par[iz]

        # Top face: 3 positions across width (W-edge, Centre, E-edge)
        top_positions = np.linspace(-cw2, cw2, 3)
        top_par = np.empty(3)
        for ip, vp in enumerate(top_positions):
            # Sample a small region around each position
            offsets = np.linspace(vp - cw2 / 3, vp + cw2 / 3, 5)
            offsets = np.clip(offsets, -cw2, cw2)
            n_sunlit = sum(1 for v in offsets if not _shaded(v, ch))
            frac = n_sunlit / len(offsets)
            top_par[ip] = total_par * frac + total_par * diffuse_fraction * (1 - frac)

        east_par = np.clip(east_par, 0, None)
        west_par = np.clip(west_par, 0, None)

        return {"east": east_par, "west": west_par, "top": top_par}

    def sunlit_fraction(self, shadow_mask: np.ndarray) -> float:
        """Fraction of canopy zones that are sunlit (0-1)."""
        return float(1.0 - shadow_mask.mean())

    def fruiting_zone_shadow(
        self,
        shadow_mask: np.ndarray,
        fruiting_zone_idx: int | None = None,
    ) -> dict:
        """Report shading specifically on the fruiting zone (mid-canopy).

        Parameters
        ----------
        shadow_mask : ndarray of shape (n_vertical, n_horizontal)
            Boolean shadow mask (True = shaded).
        fruiting_zone_idx : int, optional
            Vertical zone index for the fruiting zone. Default from settings.

        Returns
        -------
        dict with fruiting_zone_shaded_pct, fruiting_zone_sunlit_pct,
        fruiting_zone_mask (boolean array of horizontal positions).
        """
        if fruiting_zone_idx is None:
            from config.settings import FRUITING_ZONE_INDEX
            fruiting_zone_idx = FRUITING_ZONE_INDEX

        fz_row = shadow_mask[fruiting_zone_idx, :]  # shape (n_horizontal,)
        fz_shaded_fraction = float(fz_row.mean())
        return {
            "fruiting_zone_shaded_pct": round(fz_shaded_fraction * 100, 1),
            "fruiting_zone_sunlit_pct": round((1 - fz_shaded_fraction) * 100, 1),
            "fruiting_zone_mask": fz_row,
        }

    def evaluate_candidate_offsets(
        self,
        solar_elevation: float,
        solar_azimuth: float,
        theta_astro: float,
        offsets: list[int | float] | None = None,
        total_par: float = 1500.0,
        fruiting_zone_idx: int | None = None,
    ) -> dict:
        """Evaluate shadow at astronomical angle + each candidate offset.

        For each offset, computes the shadow mask at tilt = theta_astro + offset,
        then derives PAR distribution, overall sunlit fraction, and fruiting zone
        metrics. Used by the TradeoffEngine to find the minimum effective dose.

        Parameters
        ----------
        solar_elevation, solar_azimuth : float
            Current sun position (degrees).
        theta_astro : float
            Astronomical (energy-maximizing) tracker tilt (degrees).
        offsets : list of int/float, optional
            Candidate offsets to evaluate. Default from settings.CANDIDATE_OFFSETS.
        total_par : float
            Incoming PAR (umol m-2 s-1).
        fruiting_zone_idx : int, optional
            Vertical zone index for the fruiting zone. Default from settings.

        Returns
        -------
        dict keyed by offset value, each containing:
            shadow_mask, par_distribution, sunlit_fraction, fruiting_zone,
            top_canopy_sunlit_pct.
        """
        if offsets is None:
            from config.settings import CANDIDATE_OFFSETS
            offsets = CANDIDATE_OFFSETS

        results = {}
        for offset in offsets:
            theta = theta_astro + offset
            mask = self.project_shadow(solar_elevation, solar_azimuth, theta)
            par_dist = self.compute_par_distribution(
                total_par, mask,
                solar_elevation=solar_elevation,
                solar_azimuth=solar_azimuth,
                tracker_tilt=theta,
            )
            fz = self.fruiting_zone_shadow(mask, fruiting_zone_idx)

            # Top canopy (zone 2) sunlit fraction — protect > 70% for photosynthesis
            top_zone_idx = self.n_vertical - 1
            top_sunlit_pct = round(float(1.0 - mask[top_zone_idx, :].mean()) * 100, 1)

            # Mean PAR across the fruiting zone
            fz_idx = fruiting_zone_idx if fruiting_zone_idx is not None else 1
            fz_mean_par = round(float(par_dist[fz_idx, :].mean()), 1)

            results[offset] = {
                "shadow_mask": mask,
                "par_distribution": par_dist,
                "sunlit_fraction": round(self.sunlit_fraction(mask), 3),
                "fruiting_zone": fz,
                "top_canopy_sunlit_pct": top_sunlit_pct,
                "fruiting_zone_mean_par": fz_mean_par,
            }

        return results

    def _sun_cross_component(self, solar_elevation: float, solar_azimuth: float) -> tuple[float, float]:
        """Return (sun_cross, sun_z) — sun direction projected onto the cross-row plane."""
        elev_rad = np.radians(solar_elevation)
        sun_horiz = np.array([np.cos(elev_rad) * np.sin(np.radians(solar_azimuth)),
                              np.cos(elev_rad) * np.cos(np.radians(solar_azimuth))])
        sun_cross = float(np.dot(sun_horiz, self._row_v))
        sun_z = np.sin(elev_rad)
        return sun_cross, sun_z

    def compute_face_shading(
        self,
        solar_elevation: float,
        solar_azimuth: float,
        tracker_tilt: float,
        n_samples: int = 20,
        include_panels: bool = True,
    ) -> dict:
        """
        Compute shading on top, east, and west faces of the canopy box.
        Includes self-shading: faces whose outward normal points away from the
        sun receive no direct light (the vine box is opaque).

        When include_panels=False, only self-shading is computed (reference case).

        Faces:
        - Top: z = canopy_height, width = canopy_width  (normal = +z)
        - East: v = +canopy_width/2, height = canopy_height  (normal = +v, NE-facing)
        - West: v = -canopy_width/2, height = canopy_height  (normal = -v, SW-facing)
        """
        cw2 = self.canopy_width / 2
        ch = self.canopy_height
        top_area = self.canopy_width  # 0.6 m
        east_area = ch               # 1.2 m
        west_area = ch               # 1.2 m
        total_area = top_area + east_area + west_area

        _zero = {
            "top_sunlit": 0.0, "east_sunlit": 0.0, "west_sunlit": 0.0,
            "top_area": top_area, "east_area": east_area, "west_area": west_area,
            "total_sunlit_area": 0.0, "total_area": total_area,
            "sunlit_fraction": 0.0,
        }

        if solar_elevation <= 2.0:
            return _zero

        sun_cross, sun_z = self._sun_cross_component(solar_elevation, solar_azimuth)

        if sun_z <= 0.01:
            return _zero

        # --- Self-shading: faces whose normal faces away from the sun ---
        # East face normal = +v → sunlit only when sun_cross > 0
        # West face normal = -v → sunlit only when sun_cross < 0
        # Top face normal = +z  → always sunlit during daytime
        east_self_shaded = (sun_cross <= 0)
        west_self_shaded = (sun_cross >= 0)

        # If a face is self-shaded, it's 0% sunlit regardless of panels
        if east_self_shaded:
            east_sunlit = 0.0
        if west_self_shaded:
            west_sunlit = 0.0

        if include_panels:
            # Panel segments from 3 rows (own + neighbors) in target frame
            panels = []
            for off_idx in (-1, 0, 1):
                row_shift = -off_idx * self.row_spacing
                corners = self._panel_corners_local(tracker_tilt)
                p1_v = corners[0, 0] - row_shift
                p1_z = corners[0, 2]
                p2_v = corners[1, 0] - row_shift
                p2_z = corners[1, 2]
                panels.append((p1_v, p1_z, p2_v, p2_z))

            def _is_panel_shaded(v_c, z_c):
                """Check if a point is shaded by any panel."""
                for p1_v, p1_z, p2_v, p2_z in panels:
                    dv_p = p2_v - p1_v
                    dz_p = p2_z - p1_z
                    det = sun_cross * dz_p - sun_z * dv_p
                    if abs(det) < 1e-10:
                        continue
                    dv = p1_v - v_c
                    dz = p1_z - z_c
                    t = (dv * dz_p - dz * dv_p) / det
                    s = (dv * sun_z - dz * sun_cross) / det
                    if t > 0 and 0.0 <= s <= 1.0:
                        return True
                return False
        else:
            def _is_panel_shaded(v_c, z_c):
                return False

        # Sample top face: z = ch, v from -cw2 to +cw2 (never self-shaded)
        top_v = np.linspace(-cw2, cw2, n_samples)
        top_sunlit_count = sum(1 for v in top_v if not _is_panel_shaded(v, ch))
        top_sunlit = top_sunlit_count / n_samples

        # East face: skip sampling if self-shaded
        if not east_self_shaded:
            east_z = np.linspace(0, ch, n_samples)
            east_sunlit_count = sum(1 for z in east_z if not _is_panel_shaded(cw2, z))
            east_sunlit = east_sunlit_count / n_samples

        # West face: skip sampling if self-shaded
        if not west_self_shaded:
            west_z = np.linspace(0, ch, n_samples)
            west_sunlit_count = sum(1 for z in west_z if not _is_panel_shaded(-cw2, z))
            west_sunlit = west_sunlit_count / n_samples

        total_sunlit_area = top_sunlit * top_area + east_sunlit * east_area + west_sunlit * west_area

        return {
            "top_sunlit": top_sunlit,
            "east_sunlit": east_sunlit,
            "west_sunlit": west_sunlit,
            "top_area": top_area,
            "east_area": east_area,
            "west_area": west_area,
            "total_sunlit_area": total_sunlit_area,
            "total_area": total_area,
            "sunlit_fraction": total_sunlit_area / total_area,
        }

    def compute_face_shadow_bounds(
        self,
        solar_elevation: float,
        solar_azimuth: float,
        tracker_tilt: float,
        n_samples: int = 30,
        include_panels: bool = True,
    ) -> dict:
        """
        Compute shadow boundaries on each face for 3D visualization.
        Includes self-shading: faces whose normal faces away from the sun
        are entirely shaded. When include_panels=False, only self-shading.
        Returns dict with shaded z-ranges on east/west faces and v-range on top face.
        """
        cw2 = self.canopy_width / 2
        ch = self.canopy_height

        _all_shaded = {"top_shaded": (-cw2, cw2), "east_shaded": (0, ch), "west_shaded": (0, ch)}

        if solar_elevation <= 2.0:
            return _all_shaded

        sun_cross, sun_z = self._sun_cross_component(solar_elevation, solar_azimuth)

        if sun_z <= 0.01:
            return _all_shaded

        # Self-shading: faces whose normal faces away from the sun
        east_self_shaded = (sun_cross <= 0)
        west_self_shaded = (sun_cross >= 0)

        if include_panels:
            panels = []
            for off_idx in (-1, 0, 1):
                row_shift = -off_idx * self.row_spacing
                corners = self._panel_corners_local(tracker_tilt)
                panels.append((corners[0, 0] - row_shift, corners[0, 2],
                               corners[1, 0] - row_shift, corners[1, 2]))

            def _is_panel_shaded(v_c, z_c):
                for p1_v, p1_z, p2_v, p2_z in panels:
                    dv_p, dz_p = p2_v - p1_v, p2_z - p1_z
                    det = sun_cross * dz_p - sun_z * dv_p
                    if abs(det) < 1e-10:
                        continue
                    dv, dz = p1_v - v_c, p1_z - z_c
                    t = (dv * dz_p - dz * dv_p) / det
                    s = (dv * sun_z - dz * sun_cross) / det
                    if t > 0 and 0.0 <= s <= 1.0:
                        return True
                return False
        else:
            def _is_panel_shaded(v_c, z_c):
                return False

        # Top face: find shaded v-range (never self-shaded)
        top_v = np.linspace(-cw2, cw2, n_samples)
        top_shaded = [v for v in top_v if _is_panel_shaded(v, ch)]
        top_bounds = (min(top_shaded), max(top_shaded)) if top_shaded else None

        # East face: entirely self-shaded when sun_cross <= 0
        if east_self_shaded:
            east_bounds = (0, ch)
        else:
            east_z = np.linspace(0, ch, n_samples)
            east_shaded = [z for z in east_z if _is_panel_shaded(cw2, z)]
            east_bounds = (min(east_shaded), max(east_shaded)) if east_shaded else None

        # West face: entirely self-shaded when sun_cross >= 0
        if west_self_shaded:
            west_bounds = (0, ch)
        else:
            west_z = np.linspace(0, ch, n_samples)
            west_shaded = [z for z in west_z if _is_panel_shaded(-cw2, z)]
            west_bounds = (min(west_shaded), max(west_shaded)) if west_shaded else None

        return {"top_shaded": top_bounds, "east_shaded": east_bounds, "west_shaded": west_bounds}