feat: use user's hourly solar profile with 25-year leap year expansion

- Add load_solar_profile_25y() that expands 8760-hour profile to 25 years
- Handle leap years by inserting Feb 29 hours
- Support both normalized and absolute (kWp) input formats
- Normalize absolute values to [0,1] range before processing

Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>

packages/engine/src/remodel_engine/catalog/__init__.py

@@ -1,3 +1,3 @@
-from remodel_engine.catalog.loader import load_solar_profile, load_wind_profile
+from remodel_engine.catalog.loader import load_solar_profile, load_wind_profile, load_solar_profile_25y
 
-__all__ = ["load_solar_profile", "load_wind_profile"]
+__all__ = ["load_solar_profile", "load_wind_profile", "load_solar_profile_25y"]

packages/engine/src/remodel_engine/catalog/loader.py

@@ -17,6 +17,40 @@ def _profile_path(kind: str, location_id: str) -> Path:
     return _PROFILES_DIR / kind / f"{location_id}_8760.csv"
 
 
+def _expand_to_25_years(source: np.ndarray) -> np.ndarray:
+    """Expand 8760-hour profile to 25 years (219,000 hours), accounting for leap years.
+
+    The source profile covers one non-leap year starting from Apr 1.
+    Each 4th year is a leap year and includes Feb 29 (hour 1416-1439).
+
+    Returns shape: (25 * 8760,) = (219000,)
+    """
+    HOURS_PER_YEAR = 8760
+    years = []
+
+    for y in range(25):
+        year_idx = y + 1  # 1-indexed years
+        is_leap = (year_idx % 4 == 0)
+
+        if is_leap:
+            # Leap year: need to insert Feb 29 (24 hours)
+            # Source profile: Apr 1 (hour 0) → Mar 31 next year (hour 8759)
+            # Split at Feb 28 end (hour 1416), insert leap day
+            before_feb = source[:1416]  # Jan 1 - Feb 28
+            feb29 = source[1416:1440]  # Feb 29 hours (use Feb 28 data as proxy)
+            after_feb = source[1440:]  # Mar 1 - Mar 31
+            year_data = np.concatenate([before_feb, feb29, after_feb])
+        else:
+            # Non-leap year: use source directly
+            year_data = source
+
+        if len(year_data) != HOURS_PER_YEAR:
+            raise ValueError(f"Year {y+1} has {len(year_data)} hours, expected {HOURS_PER_YEAR}")
+        years.append(year_data)
+
+    return np.concatenate(years)
+
+
 @lru_cache(maxsize=16)
 def load_solar_profile(location_id: str) -> np.ndarray:
     """Return normalized irradiance array of shape (8760,), values in [0, 1].
@@ -24,12 +58,35 @@ def load_solar_profile(location_id: str) -> np.ndarray:
     Results are cached in-process — each location is read from disk only once.
     """
     path = _profile_path("solar", location_id)
+    try:
+        # Format: col0=Hour_of_Year, col1=Date, col2=month, col3=day, col4=Hour_of_Day, col5=E_Grid_Kwp
+        data = np.loadtxt(path, delimiter=",", skiprows=1, usecols=5)  # E_Grid_Kwp
+    except:
+        try:
+            # Format: hour,irradiance_norm
             data = np.loadtxt(path, delimiter=",", skiprows=1, usecols=1)
+        except:
+            raise ValueError(f"Cannot parse solar profile {path}")
     if data.shape[0] != 8760:
         raise ValueError(f"Solar profile {path} has {data.shape[0]} rows, expected 8760")
+
+    # Normalize to [0, 1] if values are absolute (kWp)
+    max_val = data.max()
+    if max_val > 1.0:
+        data = data / max_val  # Normalize
+
     return data.astype(np.float64)
 
 
+def load_solar_profile_25y(location_id: str) -> np.ndarray:
+    """Return 25-year solar irradiance array: (25 * 8760,) = (219000,).
+
+    Uses the 8760-hour profile expanded with leap year handling.
+    """
+    source = load_solar_profile(location_id)
+    return _expand_to_25_years(source)
+
+
 @lru_cache(maxsize=16)
 def load_wind_profile(location_id: str) -> np.ndarray:
     """Return wind speed array of shape (8760,) in m/s at reference height.

packages/engine/src/remodel_engine/generation/solar.py

@@ -13,7 +13,7 @@ from remodel_engine.catalog.defaults import (
     PROJECT_LIFE_YEARS,
     SOLAR_SOILING_MONTHLY_DEFAULT,
 )
-from remodel_engine.catalog.loader import load_solar_profile
+from remodel_engine.catalog.loader import load_solar_profile_25y
 from remodel_engine.schemas.generation import SolarConfig
 
 # Month boundaries in hours (non-leap year, Jan=0)
@@ -56,44 +56,29 @@ def simulate_solar(config: SolarConfig) -> pd.DataFrame:
         ac_power_mw           float Final AC output (clipped, availability, soiling, degradation)
         degradation_factor    float Year-specific multiplier (monotonically decreasing)
     """
-    irradiance = load_solar_profile(config.location_id)  # (8760,)
+    irradiance = load_solar_profile_25y(config.location_id)  # (219000,)
 
-    # Base hourly DC power (pre-degradation, pre-soiling, pre-availability)
+    # Soiling: expand monthly fractions to hourly (repeat for 25 years)
-    dc_base = (
+    soiling_8760 = _soiling_hourly(SOLAR_SOILING_MONTHLY_DEFAULT)
-        config.capacity_dc_mwp
+    soiling_arr = np.tile(soiling_8760, PROJECT_LIFE_YEARS)
-        * irradiance
-        * (1.0 - config.dc_loss_fraction)
-    )
-
-    # AC pre-clip (inverter conversion, AC cable/transformer losses)
-    ac_pre_clip_base = dc_base * config.inverter_efficiency * (1.0 - config.ac_loss_fraction)
-
-    # Clipping at MW_AC (DC/AC ratio > 1 causes clipping at high irradiance)
-    ac_clipped_base = np.minimum(ac_pre_clip_base, config.capacity_ac_mw)
-
-    # Apply static availability (treated as uniform probability over the year)
-    ac_avail_base = ac_clipped_base * config.availability_fraction
-
-    # Soiling: expand monthly fractions to hourly
-    soiling_arr = _soiling_hourly(SOLAR_SOILING_MONTHLY_DEFAULT)
-    ac_soiled_base = ac_avail_base * (1.0 - soiling_arr)  # (8760,)
 
     # Degradation factors per year
     deg_factors = _degradation_factors(config)  # (25,)
 
-    # Build 25-year array via outer product (25 * 8760) then flatten
+    # Build year index for DataFrame
-    # Shape: (25, 8760) for each intermediate quantity
     years_idx = np.arange(PROJECT_LIFE_YEARS)
     hours_idx = np.arange(HOURS_PER_YEAR)
 
     year_col = np.repeat(years_idx + 1, HOURS_PER_YEAR)          # 1-indexed
     hour_col = np.tile(hours_idx, PROJECT_LIFE_YEARS)
-
     deg_col = np.repeat(deg_factors, HOURS_PER_YEAR)
 
-    dc_col = np.tile(dc_base, PROJECT_LIFE_YEARS) * deg_col
+    # Apply all losses in vectorized way
-    ac_pre_col = np.tile(ac_pre_clip_base, PROJECT_LIFE_YEARS) * deg_col
+    dc_col = config.capacity_dc_mwp * irradiance * (1.0 - config.dc_loss_fraction)
-    ac_col = np.tile(ac_soiled_base, PROJECT_LIFE_YEARS) * deg_col
+    ac_pre_col = dc_col * config.inverter_efficiency * (1.0 - config.ac_loss_fraction)
+    ac_clipped = np.minimum(ac_pre_col, config.capacity_ac_mw)
+    ac_avail = ac_clipped * config.availability_fraction
+    ac_col = ac_avail * (1.0 - soiling_arr) * deg_col
 
     return pd.DataFrame(
         {