Bifacial PV Modeling

def get_irradiance_poa(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth,
                       gcr, height, pitch, ghi, dhi, dni, albedo,
                       model='isotropic', npoints=100):
    """
    Calculate plane-of-array irradiance for infinite sheds bifacial configuration.
    
    Parameters:
    - surface_tilt: numeric, surface tilt angle in degrees
    - surface_azimuth: numeric, surface azimuth in degrees  
    - solar_zenith: numeric, solar zenith angle in degrees
    - solar_azimuth: numeric, solar azimuth in degrees
    - gcr: numeric, ground coverage ratio (0-1)
    - height: numeric, module height above ground in meters
    - pitch: numeric, row spacing in meters
    - ghi: numeric, global horizontal irradiance in W/m²
    - dhi: numeric, diffuse horizontal irradiance in W/m²
    - dni: numeric, direct normal irradiance in W/m²
    - albedo: numeric, ground reflectance (0-1)
    - model: str, sky diffuse model ('isotropic', 'perez')
    - npoints: int, number of discretization points
    
    Returns:
    dict with keys:
    - poa_front: front-side plane-of-array irradiance
    - poa_rear: rear-side plane-of-array irradiance  
    - poa_total: total bifacial plane-of-array irradiance
    """

def get_irradiance(surface_tilt, surface_azimuth, solar_zenith, solar_azimuth,
                   gcr, height, pitch, ghi, dhi, dni, albedo, 
                   iam_front=1.0, iam_rear=1.0, bifaciality=0.8):
    """
    Calculate irradiance components for infinite sheds bifacial system.
    
    Parameters:
    - surface_tilt: numeric, surface tilt angle in degrees
    - surface_azimuth: numeric, surface azimuth in degrees
    - solar_zenith: numeric, solar zenith angle in degrees  
    - solar_azimuth: numeric, solar azimuth in degrees
    - gcr: numeric, ground coverage ratio
    - height: numeric, module height above ground in meters
    - pitch: numeric, row spacing in meters
    - ghi: numeric, global horizontal irradiance in W/m²
    - dhi: numeric, diffuse horizontal irradiance in W/m²
    - dni: numeric, direct normal irradiance in W/m²
    - albedo: numeric, ground reflectance (0-1)
    - iam_front: numeric, front-side incidence angle modifier
    - iam_rear: numeric, rear-side incidence angle modifier  
    - bifaciality: numeric, rear/front power ratio under same irradiance
    
    Returns:
    dict with detailed irradiance breakdown including:
    - front_direct, front_diffuse, front_reflected
    - rear_direct, rear_diffuse, rear_reflected  
    - total_absorbed_front, total_absorbed_rear
    """

def pvfactors_timeseries(solar_azimuth, solar_zenith, surface_azimuth, surface_tilt,
                         axis_azimuth, timestamps, dni, dhi, gcr, pvrow_height, 
                         pvrow_width, albedo, n_pvrows=3, index_observed_pvrow=1,
                         rho_front_pvrow=0.03, rho_back_pvrow=0.05, 
                         horizon_band_angle=15.0, run_parallel_calculations=True,
                         n_workers_for_parallel_calcs=2):
    """
    Run PVFactors timeseries simulation for detailed bifacial modeling.
    
    Parameters:
    - solar_azimuth: array-like, solar azimuth angles in degrees
    - solar_zenith: array-like, solar zenith angles in degrees
    - surface_azimuth: numeric, PV surface azimuth in degrees
    - surface_tilt: numeric, PV surface tilt in degrees
    - axis_azimuth: numeric, axis of rotation azimuth in degrees
    - timestamps: pandas.DatetimeIndex, simulation timestamps
    - dni: array-like, direct normal irradiance in W/m²
    - dhi: array-like, diffuse horizontal irradiance in W/m²
    - gcr: numeric, ground coverage ratio
    - pvrow_height: numeric, PV row height in meters
    - pvrow_width: numeric, PV row width in meters  
    - albedo: numeric, ground albedo (0-1)
    - n_pvrows: int, number of PV rows to model
    - index_observed_pvrow: int, index of PV row to report (0-based)
    - rho_front_pvrow: numeric, front surface reflectance
    - rho_back_pvrow: numeric, rear surface reflectance
    - horizon_band_angle: numeric, discretization angle for horizon band
    - run_parallel_calculations: bool, enable parallel processing
    - n_workers_for_parallel_calcs: int, number of parallel workers
    
    Returns:
    pandas.DataFrame with columns:
    - total_inc_front, total_inc_rear: incident irradiance on front/rear
    - total_abs_front, total_abs_rear: absorbed irradiance on front/rear
    """

def vf_ground_sky_2d(rotation, gcr, x, pitch, height, max_rows=10):
    """
    Calculate 2D view factors between ground and sky.
    
    Parameters:
    - rotation: numeric, surface rotation angle in degrees
    - gcr: numeric, ground coverage ratio
    - x: numeric, position along ground (normalized by pitch)
    - pitch: numeric, row spacing in meters
    - height: numeric, PV module height in meters
    - max_rows: int, maximum number of rows to consider
    
    Returns:
    numeric, view factor from ground point to sky
    """

def vf_ground_sky_2d_integ(surface_tilt, gcr, height, pitch, max_rows=10, 
                           x0=0, x1=1):
    """
    Calculate integrated 2D view factors between ground and sky.
    
    Parameters:
    - surface_tilt: numeric, surface tilt angle in degrees
    - gcr: numeric, ground coverage ratio
    - height: numeric, PV module height in meters
    - pitch: numeric, row spacing in meters
    - max_rows: int, maximum number of rows to consider
    - x0: numeric, integration start position (normalized)
    - x1: numeric, integration end position (normalized)
    
    Returns:
    numeric, integrated view factor
    """

def vf_row_sky_2d(surface_tilt, gcr, x):
    """
    Calculate 2D view factors between PV row and sky.
    
    Parameters:
    - surface_tilt: numeric, surface tilt angle in degrees
    - gcr: numeric, ground coverage ratio  
    - x: numeric, position along PV row (normalized)
    
    Returns:
    numeric, view factor from PV row point to sky
    """

def vf_row_ground_2d(surface_tilt, gcr, x):
    """
    Calculate 2D view factors between PV row and ground.
    
    Parameters:
    - surface_tilt: numeric, surface tilt angle in degrees
    - gcr: numeric, ground coverage ratio
    - x: numeric, position along PV row (normalized)
    
    Returns:  
    numeric, view factor from PV row point to ground
    """

def vf_row_sky_2d_integ(surface_tilt, gcr, x0=0, x1=1):
    """
    Calculate integrated 2D view factors between PV row and sky.
    
    Parameters:
    - surface_tilt: numeric, surface tilt angle in degrees
    - gcr: numeric, ground coverage ratio
    - x0: numeric, integration start position (normalized)
    - x1: numeric, integration end position (normalized)
    
    Returns:
    numeric, integrated view factor from PV row to sky
    """

def vf_row_ground_2d_integ(surface_tilt, gcr, x0=0, x1=1):
    """
    Calculate integrated 2D view factors between PV row and ground.
    
    Parameters:
    - surface_tilt: numeric, surface tilt angle in degrees
    - gcr: numeric, ground coverage ratio
    - x0: numeric, integration start position (normalized)  
    - x1: numeric, integration end position (normalized)
    
    Returns:
    numeric, integrated view factor from PV row to ground
    """

def power_mismatch_deline(poa_front, poa_rear, bifaciality=0.8, 
                          n_rear_strings=None, mismatch_factor=0.98):
    """
    Calculate power mismatch loss for bifacial modules using Deline model.
    
    Parameters:
    - poa_front: array-like, front-side plane-of-array irradiance
    - poa_rear: array-like, rear-side plane-of-array irradiance  
    - bifaciality: numeric, rear/front power ratio under same irradiance
    - n_rear_strings: int, number of rear-illuminated strings (optional)
    - mismatch_factor: numeric, electrical mismatch factor
    
    Returns:
    dict with keys:
    - p_mp_front: front-side power at maximum power point
    - p_mp_rear: rear-side power at maximum power point
    - p_mp_total: total bifacial power output
    - mismatch_loss: power loss due to mismatch
    """

def effective_irradiance_bifacial(poa_front, poa_rear, bifaciality=0.8):
    """
    Calculate effective irradiance for bifacial modules.
    
    Parameters:
    - poa_front: numeric, front-side plane-of-array irradiance in W/m²
    - poa_rear: numeric, rear-side plane-of-array irradiance in W/m²
    - bifaciality: numeric, rear/front power ratio
    
    Returns:
    numeric, effective irradiance in W/m²
    """

def bifacial_gain(poa_front, poa_rear, bifaciality=0.8):
    """
    Calculate bifacial gain compared to monofacial equivalent.
    
    Parameters:
    - poa_front: numeric, front-side irradiance in W/m²
    - poa_rear: numeric, rear-side irradiance in W/m²  
    - bifaciality: numeric, rear/front efficiency ratio
    
    Returns:
    numeric, bifacial gain factor (>1.0 indicates benefit)
    """

import pvlib
from pvlib import bifacial, solarposition, irradiance, atmosphere
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

# System parameters
lat, lon = 40.0150, -105.2705  # Boulder, CO
surface_tilt = 30  # degrees
surface_azimuth = 180  # south-facing
gcr = 0.4  # ground coverage ratio
height = 1.5  # module height above ground in meters
pitch = height / (gcr * np.cos(np.radians(surface_tilt)))  # calculated pitch
albedo = 0.25  # ground reflectance
bifaciality = 0.8  # rear/front efficiency ratio

# Time series for one day
times = pd.date_range('2023-06-21 05:00', '2023-06-21 19:00', 
                     freq='H', tz='US/Mountain')

# Calculate solar position
solar_pos = solarposition.get_solarposition(times, lat, lon)

# Get clear sky irradiance
clear_sky = pvlib.clearsky.ineichen(times, lat, lon, altitude=1655)

# Calculate front-side POA irradiance (traditional method)
poa_front_traditional = irradiance.get_total_irradiance(
    surface_tilt, surface_azimuth,
    solar_pos['zenith'], solar_pos['azimuth'],
    clear_sky['dni'], clear_sky['ghi'], clear_sky['dhi'],
    albedo=albedo
)

# Calculate bifacial irradiance using infinite sheds model
bifacial_results = []

for i, time in enumerate(times):
    # Get irradiance for this timestep
    result = bifacial.infinite_sheds.get_irradiance_poa(
        surface_tilt=surface_tilt,
        surface_azimuth=surface_azimuth,
        solar_zenith=solar_pos['zenith'].iloc[i],
        solar_azimuth=solar_pos['azimuth'].iloc[i],
        gcr=gcr,
        height=height,
        pitch=pitch,
        ghi=clear_sky['ghi'].iloc[i],
        dhi=clear_sky['dhi'].iloc[i], 
        dni=clear_sky['dni'].iloc[i],
        albedo=albedo,
        model='isotropic'
    )
    
    bifacial_results.append({
        'time': time,
        'poa_front': result['poa_front'],
        'poa_rear': result['poa_rear'],
        'poa_total': result['poa_total']
    })

# Convert to DataFrame
bifacial_df = pd.DataFrame(bifacial_results)
bifacial_df.set_index('time', inplace=True)

# Calculate bifacial gains
bifacial_df['effective_irradiance'] = bifacial.effective_irradiance_bifacial(
    bifacial_df['poa_front'], bifacial_df['poa_rear'], bifaciality
)

bifacial_df['bifacial_gain'] = bifacial.bifacial_gain(
    bifacial_df['poa_front'], bifacial_df['poa_rear'], bifaciality
)

# Calculate daily totals
daily_totals = {
    'Front POA': bifacial_df['poa_front'].sum() / 1000,  # kWh/m²
    'Rear POA': bifacial_df['poa_rear'].sum() / 1000,
    'Effective': bifacial_df['effective_irradiance'].sum() / 1000,
    'Traditional': poa_front_traditional['poa_global'].sum() / 1000
}

print("Daily Irradiation Summary (kWh/m²):")
for key, value in daily_totals.items():
    print(f"{key:12s}: {value:.2f}")

daily_gain = (daily_totals['Effective'] - daily_totals['Traditional']) / daily_totals['Traditional'] * 100
print(f"\nDaily Bifacial Gain: {daily_gain:.1f}%")

# Plot results
fig, axes = plt.subplots(3, 1, figsize=(12, 10))

# Irradiance components
axes[0].plot(bifacial_df.index, bifacial_df['poa_front'], 
            label='Front POA', linewidth=2, color='blue')
axes[0].plot(bifacial_df.index, bifacial_df['poa_rear'], 
            label='Rear POA', linewidth=2, color='red')
axes[0].plot(bifacial_df.index, poa_front_traditional['poa_global'], 
            label='Traditional POA', linewidth=2, linestyle='--', color='gray')
axes[0].set_ylabel('Irradiance (W/m²)')
axes[0].set_title('Bifacial vs Traditional Irradiance')
axes[0].legend()
axes[0].grid(True)

# Effective irradiance and gain
axes[1].plot(bifacial_df.index, bifacial_df['effective_irradiance'], 
            label='Effective Irradiance', linewidth=2, color='green')
axes[1].set_ylabel('Effective Irradiance (W/m²)')
axes[1].set_title('Bifacial Effective Irradiance')
axes[1].legend()
axes[1].grid(True)

# Bifacial gain factor
axes[2].plot(bifacial_df.index, bifacial_df['bifacial_gain'], 
            'o-', linewidth=2, color='purple')
axes[2].set_ylabel('Bifacial Gain Factor')
axes[2].set_xlabel('Time')
axes[2].set_title('Bifacial Gain Throughout Day')
axes[2].grid(True)
axes[2].axhline(y=1.0, color='k', linestyle='--', alpha=0.5)

plt.tight_layout()
plt.show()

# Analyze rear irradiance sources
print(f"\nRear Irradiance Analysis:")
print(f"Peak rear irradiance: {bifacial_df['poa_rear'].max():.1f} W/m²")
print(f"Average rear/front ratio: {(bifacial_df['poa_rear'] / bifacial_df['poa_front']).mean():.3f}")
print(f"Peak bifacial gain: {bifacial_df['bifacial_gain'].max():.3f}")

import pvlib
from pvlib import bifacial, solarposition, iotools
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

# Load weather data
lat, lon = 37.8756, -122.2441  # Berkeley, CA
api_key = 'DEMO_KEY'  # Replace with actual API key
email = 'user@example.com'

# Get one week of weather data
try:
    weather_data, metadata = iotools.get_psm3(
        lat, lon, api_key, email, 
        names=2023, 
        attributes=['ghi', 'dni', 'dhi'],
        leap_day=False
    )
    # Select one week in summer
    start_date = '2023-06-15'
    end_date = '2023-06-22'
    weather_week = weather_data[start_date:end_date]
    
except:
    # Fallback to generated clear sky data
    print("Using clear sky data (replace with actual weather data)")
    times = pd.date_range('2023-06-15', '2023-06-22', freq='H', tz='US/Pacific')
    clear_sky = pvlib.clearsky.ineichen(times, lat, lon)
    weather_week = clear_sky.rename(columns={'ghi': 'ghi', 'dni': 'dni', 'dhi': 'dhi'})

# System configuration
system_params = {
    'surface_tilt': 25,         # degrees
    'surface_azimuth': 180,     # south-facing
    'axis_azimuth': 180,        # tracking axis orientation
    'gcr': 0.35,               # ground coverage ratio
    'pvrow_height': 2.0,        # meters
    'pvrow_width': 4.0,         # meters
    'albedo': 0.20,            # ground reflectance
    'n_pvrows': 5,             # number of rows to model
    'index_observed_pvrow': 2,  # middle row (0-indexed)
    'bifaciality': 0.85        # rear/front efficiency ratio
}

# Calculate solar position
solar_pos = solarposition.get_solarposition(weather_week.index, lat, lon)

# Run PVFactors simulation
print("Running PVFactors simulation...")
try:
    pvfactors_results = bifacial.pvfactors.pvfactors_timeseries(
        solar_azimuth=solar_pos['azimuth'],
        solar_zenith=solar_pos['zenith'],
        surface_azimuth=system_params['surface_azimuth'],
        surface_tilt=system_params['surface_tilt'],
        axis_azimuth=system_params['axis_azimuth'],
        timestamps=weather_week.index,
        dni=weather_week['dni'],
        dhi=weather_week['dhi'],
        gcr=system_params['gcr'],
        pvrow_height=system_params['pvrow_height'],
        pvrow_width=system_params['pvrow_width'],
        albedo=system_params['albedo'],
        n_pvrows=system_params['n_pvrows'],
        index_observed_pvrow=system_params['index_observed_pvrow'],
        run_parallel_calculations=True
    )
    
    pvfactors_available = True
    
except Exception as e:
    print(f"PVFactors not available: {e}")
    print("Using simplified infinite sheds model...")
    pvfactors_available = False

# Compare with infinite sheds model for validation
infinite_sheds_results = []

for i, (timestamp, weather_row) in enumerate(weather_week.iterrows()):
    if i % 24 == 0:  # Progress indicator
        print(f"Processing day {i//24 + 1}/7...")
    
    result = bifacial.infinite_sheds.get_irradiance(
        surface_tilt=system_params['surface_tilt'],
        surface_azimuth=system_params['surface_azimuth'],
        solar_zenith=solar_pos['zenith'].iloc[i],
        solar_azimuth=solar_pos['azimuth'].iloc[i],
        gcr=system_params['gcr'],
        height=system_params['pvrow_height'],
        pitch=system_params['pvrow_width'] / system_params['gcr'],
        ghi=weather_row['ghi'],
        dhi=weather_row['dhi'],
        dni=weather_row['dni'],
        albedo=system_params['albedo'],
        bifaciality=system_params['bifaciality']
    )
    
    infinite_sheds_results.append({
        'timestamp': timestamp,
        'front_poa': result['total_absorbed_front'],
        'rear_poa': result['total_absorbed_rear']
    })

infinite_df = pd.DataFrame(infinite_sheds_results)
infinite_df.set_index('timestamp', inplace=True)

# Calculate performance metrics
if pvfactors_available:
    # PVFactors results
    pvf_front = pvfactors_results['total_abs_front']
    pvf_rear = pvfactors_results['total_abs_rear']
    
    # Effective irradiance
    pvf_effective = pvf_front + system_params['bifaciality'] * pvf_rear
    
    # Compare models
    comparison_df = pd.DataFrame({
        'pvf_front': pvf_front,
        'pvf_rear': pvf_rear,
        'pvf_effective': pvf_effective,
        'inf_front': infinite_df['front_poa'],
        'inf_rear': infinite_df['rear_poa'],
        'inf_effective': infinite_df['front_poa'] + system_params['bifaciality'] * infinite_df['rear_poa']
    })
    
    # Statistical comparison
    print("\nModel Comparison Statistics:")
    print("="*40)
    
    for component in ['front', 'rear', 'effective']:
        pvf_col = f'pvf_{component}'
        inf_col = f'inf_{component}'
        
        mae = np.mean(np.abs(comparison_df[pvf_col] - comparison_df[inf_col]))
        rmse = np.sqrt(np.mean((comparison_df[pvf_col] - comparison_df[inf_col])**2))
        mbe = np.mean(comparison_df[pvf_col] - comparison_df[inf_col])
        
        print(f"{component.title()} Irradiance:")
        print(f"  MAE:  {mae:.2f} W/m²")
        print(f"  RMSE: {rmse:.2f} W/m²") 
        print(f"  MBE:  {mbe:.2f} W/m²")

# Plot detailed analysis
if pvfactors_available:
    fig, axes = plt.subplots(4, 1, figsize=(15, 12))
    
    # Daily profiles comparison
    sample_day = comparison_df.index.date[len(comparison_df)//2]  # Middle day
    day_data = comparison_df[comparison_df.index.date == sample_day]
    
    axes[0].plot(day_data.index.hour, day_data['pvf_front'], 
                'b-', label='PVFactors Front', linewidth=2)
    axes[0].plot(day_data.index.hour, day_data['inf_front'], 
                'b--', label='Infinite Sheds Front', linewidth=2)
    axes[0].set_ylabel('Front Irradiance (W/m²)')
    axes[0].set_title(f'Front Irradiance Comparison - {sample_day}')
    axes[0].legend()
    axes[0].grid(True)
    
    axes[1].plot(day_data.index.hour, day_data['pvf_rear'], 
                'r-', label='PVFactors Rear', linewidth=2)
    axes[1].plot(day_data.index.hour, day_data['inf_rear'], 
                'r--', label='Infinite Sheds Rear', linewidth=2)
    axes[1].set_ylabel('Rear Irradiance (W/m²)')
    axes[1].set_title('Rear Irradiance Comparison')
    axes[1].legend()
    axes[1].grid(True)
    
    # Scatter plots
    axes[2].scatter(comparison_df['inf_front'], comparison_df['pvf_front'], 
                   alpha=0.6, s=10, label='Front')
    axes[2].scatter(comparison_df['inf_rear'], comparison_df['pvf_rear'], 
                   alpha=0.6, s=10, label='Rear')
    
    max_val = max(comparison_df[['pvf_front', 'pvf_rear', 'inf_front', 'inf_rear']].max())
    axes[2].plot([0, max_val], [0, max_val], 'k--', alpha=0.5)
    axes[2].set_xlabel('Infinite Sheds (W/m²)')
    axes[2].set_ylabel('PVFactors (W/m²)')
    axes[2].set_title('Model Correlation')
    axes[2].legend()
    axes[2].grid(True)
    
    # Weekly energy totals
    daily_totals = comparison_df.resample('D').sum() / 1000  # kWh/m²/day
    
    x_days = range(len(daily_totals))
    width = 0.35
    
    axes[3].bar([x - width/2 for x in x_days], daily_totals['pvf_effective'], 
               width, label='PVFactors', alpha=0.7)
    axes[3].bar([x + width/2 for x in x_days], daily_totals['inf_effective'], 
               width, label='Infinite Sheds', alpha=0.7)
    
    axes[3].set_xlabel('Day')
    axes[3].set_ylabel('Daily Energy (kWh/m²/day)')
    axes[3].set_title('Daily Effective Energy Comparison')
    axes[3].legend()
    axes[3].grid(True, axis='y')
    
    plt.tight_layout()
    plt.show()
    
    # Weekly summary
    weekly_totals = {
        'PVFactors Front': comparison_df['pvf_front'].sum() / 1000,
        'PVFactors Rear': comparison_df['pvf_rear'].sum() / 1000,
        'PVFactors Effective': comparison_df['pvf_effective'].sum() / 1000,
        'Infinite Sheds Front': comparison_df['inf_front'].sum() / 1000,
        'Infinite Sheds Rear': comparison_df['inf_rear'].sum() / 1000,
        'Infinite Sheds Effective': comparison_df['inf_effective'].sum() / 1000
    }
    
    print(f"\nWeekly Energy Summary (kWh/m²):")
    print("="*40)
    for key, value in weekly_totals.items():
        print(f"{key:25s}: {value:.2f}")

else:
    # Plot infinite sheds results only
    fig, axes = plt.subplots(3, 1, figsize=(12, 10))
    
    # Calculate effective irradiance
    infinite_df['effective'] = (infinite_df['front_poa'] + 
                               system_params['bifaciality'] * infinite_df['rear_poa'])
    
    # Time series
    axes[0].plot(infinite_df.index, infinite_df['front_poa'], 
                label='Front POA', linewidth=1.5)
    axes[0].plot(infinite_df.index, infinite_df['rear_poa'], 
                label='Rear POA', linewidth=1.5)
    axes[0].set_ylabel('Irradiance (W/m²)')
    axes[0].set_title('Bifacial Irradiance - Infinite Sheds Model')
    axes[0].legend()
    axes[0].grid(True)
    
    # Daily totals
    daily_totals = infinite_df.resample('D').sum() / 1000
    
    axes[1].bar(range(len(daily_totals)), daily_totals['front_poa'], 
               alpha=0.7, label='Front')
    axes[1].bar(range(len(daily_totals)), daily_totals['rear_poa'], 
               bottom=daily_totals['front_poa'], alpha=0.7, label='Rear')
    axes[1].set_xlabel('Day')
    axes[1].set_ylabel('Daily Energy (kWh/m²/day)')
    axes[1].set_title('Daily Energy Breakdown')
    axes[1].legend()
    axes[1].grid(True, axis='y')
    
    # Bifacial gain analysis
    bifacial_gain = infinite_df['effective'] / infinite_df['front_poa']
    axes[2].plot(infinite_df.index, bifacial_gain, 'g-', linewidth=1.5)
    axes[2].set_ylabel('Bifacial Gain Factor')
    axes[2].set_xlabel('Time')
    axes[2].set_title('Bifacial Gain Throughout Week')
    axes[2].grid(True)
    axes[2].axhline(y=1.0, color='k', linestyle='--', alpha=0.5)
    
    plt.tight_layout()
    plt.show()
    
    print(f"\nWeekly Summary (Infinite Sheds Model):")
    print("="*40)
    print(f"Front Energy:     {daily_totals['front_poa'].sum():.2f} kWh/m²")
    print(f"Rear Energy:      {daily_totals['rear_poa'].sum():.2f} kWh/m²") 
    print(f"Effective Energy: {daily_totals['effective'].sum():.2f} kWh/m²")
    print(f"Average Bifacial Gain: {bifacial_gain.mean():.3f}")
    print(f"Peak Bifacial Gain:    {bifacial_gain.max():.3f}")

import pvlib
from pvlib import bifacial
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize_scalar

# Parameter ranges for analysis
surface_tilts = np.arange(0, 61, 5)  # 0 to 60 degrees
gcrs = np.arange(0.1, 0.81, 0.05)   # 0.1 to 0.8
heights = np.arange(0.5, 4.1, 0.25) # 0.5 to 4.0 meters

# Fixed parameters
pitch_base = 5.0  # meters

def analyze_view_factors(surface_tilt, gcr, height):
    """
    Analyze view factors for given system configuration.
    """
    pitch = height / (gcr * np.cos(np.radians(surface_tilt)))
    
    # Ground view factors at different positions
    x_positions = np.linspace(0, 1, 21)  # 0 to 1 (normalized by pitch)
    vf_ground_sky = []
    
    for x in x_positions:
        vf = bifacial.utils.vf_ground_sky_2d(
            rotation=surface_tilt, gcr=gcr, x=x, 
            pitch=pitch, height=height, max_rows=10
        )
        vf_ground_sky.append(vf)
    
    # Row view factors
    vf_row_sky_integrated = bifacial.utils.vf_row_sky_2d_integ(surface_tilt, gcr)
    vf_row_ground_integrated = bifacial.utils.vf_row_ground_2d_integ(surface_tilt, gcr)
    
    # Ground integrated view factor
    vf_ground_sky_integrated = bifacial.utils.vf_ground_sky_2d_integ(
        surface_tilt, gcr, height, pitch, max_rows=10
    )
    
    return {
        'vf_ground_sky_positions': vf_ground_sky,
        'x_positions': x_positions,
        'vf_row_sky': vf_row_sky_integrated,
        'vf_row_ground': vf_row_ground_integrated,
        'vf_ground_sky_avg': vf_ground_sky_integrated,
        'pitch': pitch
    }

# Analyze view factor sensitivity
print("Analyzing view factor sensitivity...")

# Effect of tilt angle
tilt_analysis = []
for tilt in surface_tilts:
    result = analyze_view_factors(tilt, gcr=0.4, height=2.0)
    tilt_analysis.append({
        'tilt': tilt,
        'vf_row_sky': result['vf_row_sky'],
        'vf_row_ground': result['vf_row_ground'],
        'vf_ground_sky': result['vf_ground_sky_avg']
    })

# Effect of GCR
gcr_analysis = []
for gcr in gcrs:
    result = analyze_view_factors(surface_tilt=30, gcr=gcr, height=2.0)
    gcr_analysis.append({
        'gcr': gcr,
        'vf_row_sky': result['vf_row_sky'],
        'vf_row_ground': result['vf_row_ground'], 
        'vf_ground_sky': result['vf_ground_sky_avg']
    })

# Effect of height
height_analysis = []
for height in heights:
    result = analyze_view_factors(surface_tilt=30, gcr=0.4, height=height)
    height_analysis.append({
        'height': height,
        'vf_row_sky': result['vf_row_sky'],
        'vf_row_ground': result['vf_row_ground'],
        'vf_ground_sky': result['vf_ground_sky_avg']
    })

# Plot view factor analysis
fig, axes = plt.subplots(2, 3, figsize=(18, 12))

# Tilt angle effects
tilts = [item['tilt'] for item in tilt_analysis]
axes[0, 0].plot(tilts, [item['vf_row_sky'] for item in tilt_analysis], 
               'b-o', label='Row→Sky', linewidth=2)
axes[0, 0].plot(tilts, [item['vf_row_ground'] for item in tilt_analysis], 
               'r-s', label='Row→Ground', linewidth=2)
axes[0, 0].plot(tilts, [item['vf_ground_sky'] for item in tilt_analysis], 
               'g-^', label='Ground→Sky', linewidth=2)
axes[0, 0].set_xlabel('Surface Tilt (degrees)')
axes[0, 0].set_ylabel('View Factor')
axes[0, 0].set_title('View Factors vs Surface Tilt')
axes[0, 0].legend()
axes[0, 0].grid(True)

# GCR effects
gcr_values = [item['gcr'] for item in gcr_analysis]
axes[0, 1].plot(gcr_values, [item['vf_row_sky'] for item in gcr_analysis], 
               'b-o', label='Row→Sky', linewidth=2)
axes[0, 1].plot(gcr_values, [item['vf_row_ground'] for item in gcr_analysis], 
               'r-s', label='Row→Ground', linewidth=2)
axes[0, 1].plot(gcr_values, [item['vf_ground_sky'] for item in gcr_analysis], 
               'g-^', label='Ground→Sky', linewidth=2)
axes[0, 1].set_xlabel('Ground Coverage Ratio')
axes[0, 1].set_ylabel('View Factor')
axes[0, 1].set_title('View Factors vs GCR')
axes[0, 1].legend()
axes[0, 1].grid(True)

# Height effects
height_values = [item['height'] for item in height_analysis]
axes[0, 2].plot(height_values, [item['vf_row_sky'] for item in height_analysis], 
               'b-o', label='Row→Sky', linewidth=2)
axes[0, 2].plot(height_values, [item['vf_row_ground'] for item in height_analysis], 
               'r-s', label='Row→Ground', linewidth=2)
axes[0, 2].plot(height_values, [item['vf_ground_sky'] for item in height_analysis], 
               'g-^', label='Ground→Sky', linewidth=2)
axes[0, 2].set_xlabel('Module Height (m)')
axes[0, 2].set_ylabel('View Factor')
axes[0, 2].set_title('View Factors vs Height')
axes[0, 2].legend()
axes[0, 2].grid(True)

# Detailed position analysis for specific configuration
detail_config = analyze_view_factors(surface_tilt=30, gcr=0.4, height=2.0)

# Ground view factors along pitch
axes[1, 0].plot(detail_config['x_positions'], detail_config['vf_ground_sky_positions'], 
               'ko-', linewidth=2, markersize=6)
axes[1, 0].set_xlabel('Position (normalized by pitch)')
axes[1, 0].set_ylabel('Ground→Sky View Factor')
axes[1, 0].set_title('Ground View Factor Variation')
axes[1, 0].grid(True)

# 2D heatmap of rear irradiance benefit (simplified model)
tilt_grid, gcr_grid = np.meshgrid(np.arange(10, 51, 5), np.arange(0.2, 0.71, 0.05))
benefit_matrix = np.zeros_like(tilt_grid)

for i, gcr_val in enumerate(np.arange(0.2, 0.71, 0.05)):
    for j, tilt_val in enumerate(np.arange(10, 51, 5)):
        result = analyze_view_factors(tilt_val, gcr_val, height=2.0)
        # Simplified rear irradiance benefit metric
        benefit = result['vf_row_ground'] * result['vf_ground_sky_avg']
        benefit_matrix[i, j] = benefit

contour = axes[1, 1].contourf(tilt_grid, gcr_grid, benefit_matrix, levels=20, cmap='viridis')
axes[1, 1].set_xlabel('Surface Tilt (degrees)')
axes[1, 1].set_ylabel('Ground Coverage Ratio')
axes[1, 1].set_title('Rear Irradiance Benefit (Relative)')
plt.colorbar(contour, ax=axes[1, 1])

# Optimization example - find optimal GCR for maximum rear benefit
def rear_benefit_objective(gcr, tilt=30, height=2.0):
    """
    Objective function for rear irradiance benefit (to maximize).
    Returns negative value for minimization.
    """
    result = analyze_view_factors(tilt, gcr, height)
    # Simple benefit metric: product of relevant view factors
    benefit = result['vf_row_ground'] * result['vf_ground_sky_avg']
    return -benefit  # Negative for minimization

# Find optimal GCR for different tilt angles
optimal_gcrs = []
for tilt in np.arange(10, 51, 10):
    opt_result = minimize_scalar(
        rear_benefit_objective, 
        bounds=(0.1, 0.8), 
        method='bounded',
        args=(tilt, 2.0)
    )
    optimal_gcrs.append({
        'tilt': tilt,
        'optimal_gcr': opt_result.x,
        'benefit': -opt_result.fun
    })

opt_tilts = [item['tilt'] for item in optimal_gcrs]
opt_gcr_values = [item['optimal_gcr'] for item in optimal_gcrs]

axes[1, 2].plot(opt_tilts, opt_gcr_values, 'ro-', linewidth=2, markersize=8)
axes[1, 2].set_xlabel('Surface Tilt (degrees)')
axes[1, 2].set_ylabel('Optimal GCR')
axes[1, 2].set_title('Optimal GCR vs Tilt for Max Rear Benefit')
axes[1, 2].grid(True)

plt.tight_layout()
plt.show()

# Print optimization results
print("\nView Factor Optimization Results:")
print("="*50)
print("Optimal GCR for maximum rear irradiance benefit:")
for item in optimal_gcrs:
    print(f"Tilt {item['tilt']:2.0f}°: Optimal GCR = {item['optimal_gcr']:.3f}, "
          f"Benefit = {item['benefit']:.4f}")

# Configuration recommendations
print(f"\nConfiguration Recommendations:")
print("="*30)
print("High rear irradiance configurations:")
print("- Lower GCR (0.3-0.5) increases ground view of sky")
print("- Moderate tilt (20-35°) balances front and rear performance")  
print("- Higher mounting (2-3m) reduces inter-row shading")
print("- Light-colored ground surface increases reflected irradiance")

# Practical example with recommendations
print(f"\nPractical Design Example:")
print("="*25)
recommended_config = analyze_view_factors(surface_tilt=25, gcr=0.35, height=2.5)
print(f"Recommended: 25° tilt, GCR=0.35, height=2.5m")
print(f"Row→Ground VF: {recommended_config['vf_row_ground']:.3f}")
print(f"Ground→Sky VF: {recommended_config['vf_ground_sky_avg']:.3f}")
print(f"Row→Sky VF:    {recommended_config['vf_row_sky']:.3f}")
print(f"Pitch:         {recommended_config['pitch']:.1f} m")