| """ |
| 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 |
|
|
| |
| az_rad = np.radians(row_azimuth) |
| self._row_u = np.array([np.sin(az_rad), np.cos(az_rad)]) |
| self._row_v = np.array([np.cos(az_rad), -np.sin(az_rad)]) |
|
|
| |
| 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, |
| ) |
| |
| self.lai_weights = np.array([0.15, 0.35, 0.50]) |
|
|
| 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 |
|
|
| |
| |
| 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, |
| ) |
|
|
| |
| 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) |
| cos_t, sin_t = np.cos(tilt_rad), np.sin(tilt_rad) |
| half_len = 1.0 |
|
|
| |
| 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): |
| |
| 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 |
| |
| 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"] |
|
|
| |
| 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) |
|
|
| |
| row_shift = target_row_offset - source_row_offset |
| corners_local = self._panel_corners_local(tracker_tilt) |
| |
| p1_v = corners_local[0, 0] - row_shift |
| p1_z = corners_local[0, 2] |
| p2_v = corners_local[1, 0] - row_shift |
| p2_z = corners_local[1, 2] |
|
|
| |
| dv_p = p2_v - p1_v |
| dz_p = p2_z - p1_z |
|
|
| mask = np.zeros((self.n_vertical, self.n_horizontal), dtype=bool) |
|
|
| |
| |
| |
| |
| |
| det = sun_cross * dz_p - sun_z * dv_p |
| if abs(det) < 1e-10: |
| return mask |
|
|
| 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 |
| |
| |
| 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) |
|
|
| |
| par[~shadow_mask] = total_par |
|
|
| |
| 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) |
|
|
| |
| v_norm = np.linspace(-1.0, 1.0, n_h) |
|
|
| |
| |
| |
| |
| |
| if abs(sun_cross) > 0.01: |
| |
| if sun_cross > 0: |
| edge_proximity = np.clip((v_norm + 1) / 2, 0, 1) |
| else: |
| edge_proximity = np.clip((1 - v_norm) / 2, 0, 1) |
|
|
| |
| |
| |
| horiz_factor = abs(sun_cross) / (abs(sun_cross) + sun_z) |
|
|
| |
| gap_fraction = 1.0 - self.panel_width / self.row_spacing |
|
|
| |
| side_par = total_par * horiz_factor * gap_fraction * edge_proximity |
|
|
| |
| for iz in range(n_v): |
| for ih in range(n_h): |
| if shadow_mask[iz, ih]: |
| par[iz, ih] += side_par[ih] |
|
|
| |
| |
| |
| z_fracs = np.linspace(0.0, 1.0, n_v) |
| sky_view = 0.4 + 0.6 * z_fracs |
| 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]: |
| |
| 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 |
| ch = self.canopy_height |
| cw2 = self.canopy_width / 2 |
|
|
| |
| zone_z = np.linspace(ch / n_z / 2, ch - ch / n_z / 2, n_z) |
|
|
| |
| z_fracs = np.linspace(0.0, 1.0, n_z) |
| sky_view = 0.4 + 0.6 * z_fracs |
| 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} |
|
|
| |
| 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 |
|
|
| |
| 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 |
|
|
| |
| east_self_shaded = (sun_cross <= 0) |
| west_self_shaded = (sun_cross >= 0) |
|
|
| |
| for iz in range(n_z): |
| if east_self_shaded: |
| east_par[iz] = diffuse_par[iz] * 0.3 |
| elif not _shaded(cw2, zone_z[iz]): |
| east_par[iz] = side_direct + diffuse_par[iz] |
| |
|
|
| |
| 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_positions = np.linspace(-cw2, cw2, 3) |
| top_par = np.empty(3) |
| for ip, vp in enumerate(top_positions): |
| |
| 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, :] |
| 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_zone_idx = self.n_vertical - 1 |
| top_sunlit_pct = round(float(1.0 - mask[top_zone_idx, :].mean()) * 100, 1) |
|
|
| |
| 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 |
| east_area = ch |
| west_area = ch |
| 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 |
|
|
| |
| |
| |
| |
| east_self_shaded = (sun_cross <= 0) |
| west_self_shaded = (sun_cross >= 0) |
|
|
| |
| if east_self_shaded: |
| east_sunlit = 0.0 |
| if west_self_shaded: |
| west_sunlit = 0.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) |
| 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 |
|
|
| |
| 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 |
|
|
| |
| 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 |
|
|
| |
| 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 |
|
|
| |
| 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_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 |
|
|
| |
| 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 |
|
|
| |
| 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} |
|
|
