Spectral Modeling

def spectrl2(apparent_zenith, aoi, surface_tilt, ground_albedo, 
             surface_pressure, relative_humidity, precipitable_water, 
             ozone=0.31, nitrogen_dioxide=0.0002, aerosol_turbidity_500nm=0.1,
             dayofyear=1, wavelengths=None):
    """
    Calculate spectral irradiance using SPECTRL2 model.
    
    Parameters:
    - apparent_zenith: numeric, apparent solar zenith angle in degrees
    - aoi: numeric, angle of incidence on surface in degrees  
    - surface_tilt: numeric, surface tilt angle in degrees
    - ground_albedo: numeric, ground reflectance (0-1)
    - surface_pressure: numeric, surface pressure in millibars
    - relative_humidity: numeric, relative humidity (0-100)
    - precipitable_water: numeric, precipitable water in cm
    - ozone: numeric, ozone column thickness in atm-cm
    - nitrogen_dioxide: numeric, nitrogen dioxide column in atm-cm
    - aerosol_turbidity_500nm: numeric, aerosol optical depth at 500nm
    - dayofyear: int, day of year (1-366)
    - wavelengths: array-like, wavelengths in nm (280-4000)
    
    Returns:
    tuple: (wavelengths, total_irradiance, direct_irradiance, diffuse_irradiance, 
            global_horizontal_irradiance, direct_normal_irradiance)
    """

def get_reference_spectra(wavelengths=None, standard="ASTM G173-03"):
    """
    Get standard reference solar spectra.
    
    Parameters:
    - wavelengths: array-like, wavelengths in nm for interpolation
    - standard: str, reference standard ('ASTM G173-03', 'AM15G', 'AM15D', 'AM0')
    
    Returns:
    pandas.DataFrame with wavelengths and spectral irradiance columns
    """

def average_photon_energy(spectra):
    """
    Calculate average photon energy of spectrum.
    
    Parameters:
    - spectra: pandas.DataFrame with wavelengths (nm) and irradiance columns
    
    Returns:
    pandas.Series with average photon energies in eV
    """

def spectral_factor_sapm(airmass_absolute, module):
    """
    Calculate SAPM spectral correction factor.
    
    Parameters:
    - airmass_absolute: numeric, absolute airmass  
    - module: dict, module parameters containing spectral coefficients
      Required keys: 'a0', 'a1', 'a2', 'a3', 'a4'
    
    Returns:
    numeric, spectral correction factor (typically 0.8-1.2)
    """

def spectral_factor_firstsolar(precipitable_water, airmass_absolute, 
                               module_type, coefficients=None):
    """
    Calculate First Solar spectral correction factor.
    
    Parameters:
    - precipitable_water: numeric, precipitable water in cm
    - airmass_absolute: numeric, absolute airmass
    - module_type: str, module technology ('cdte', 'multijunction', 'cigs')
    - coefficients: dict, custom spectral coefficients
    
    Returns:
    numeric, spectral correction factor
    """

def spectral_factor_caballero(precipitable_water, airmass_absolute, aod500, 
                              alpha=1.14, module_type='csi'):
    """
    Calculate Caballero spectral correction factor.
    
    Parameters:
    - precipitable_water: numeric, precipitable water in cm
    - airmass_absolute: numeric, absolute airmass
    - aod500: numeric, aerosol optical depth at 500nm
    - alpha: numeric, Angstrom alpha parameter
    - module_type: str, module technology ('csi', 'asi', 'cdte', 'pc')
    
    Returns:
    numeric, spectral correction factor
    """

def spectral_factor_pvspec(airmass_absolute, clearsky_index, 
                           module_type='monosi', 
                           coefficients=None):
    """
    Calculate PVSPEC spectral correction factor.
    
    Parameters:
    - airmass_absolute: numeric, absolute airmass (1.0-10.0)
    - clearsky_index: numeric, clearsky index (0.0-1.0)
    - module_type: str, module type ('monosi', 'xsi', 'cdte', 'asi', 'cigs', 'perovskite')
    - coefficients: dict, custom model coefficients
    
    Returns:
    numeric, spectral correction factor
    """

def spectral_factor_jrc(airmass, clearsky_index, module_type=None, 
                        pwat=None, am_coeff=None):
    """
    Calculate JRC spectral correction factor.
    
    Parameters:
    - airmass: numeric, airmass (1.0-10.0)
    - clearsky_index: numeric, clearsky index (0.0-1.0)  
    - module_type: str, technology type ('csi', 'asi', 'cdte', 'cigs')
    - pwat: numeric, precipitable water in cm (optional)
    - am_coeff: dict, custom airmass coefficients
    
    Returns:
    numeric, spectral correction factor
    """

def calc_spectral_mismatch_field(sr, e_sun, e_ref=None):
    """
    Calculate spectral mismatch for field conditions.
    
    Parameters:
    - sr: pandas.Series, spectral response vs wavelength (nm)
    - e_sun: pandas.Series, solar spectral irradiance vs wavelength
    - e_ref: pandas.Series, reference spectral irradiance (default AM1.5G)
    
    Returns:
    numeric, spectral mismatch factor
    """

def get_example_spectral_response(wavelength=None):
    """
    Get example spectral response functions for common PV technologies.
    
    Parameters:
    - wavelength: array-like, wavelengths in nm for interpolation
    
    Returns:
    dict with technology names as keys and spectral response arrays as values
    Technologies: 'si_c', 'si_pc', 'cdte', 'cigs', 'asi', 'gaas', 'organic'
    """

def sr_to_qe(sr, wavelength=None, normalize=False):
    """
    Convert spectral response to quantum efficiency.
    
    Parameters:
    - sr: pandas.Series, spectral response (A/W) vs wavelength
    - wavelength: array-like, wavelengths in nm (default from sr index)
    - normalize: bool, normalize QE to peak value
    
    Returns:
    pandas.Series, quantum efficiency vs wavelength
    """

def qe_to_sr(qe, wavelength=None, normalize=False):
    """
    Convert quantum efficiency to spectral response.
    
    Parameters:
    - qe: pandas.Series, quantum efficiency vs wavelength  
    - wavelength: array-like, wavelengths in nm (default from qe index)
    - normalize: bool, normalize SR to peak value
    
    Returns:
    pandas.Series, spectral response (A/W) vs wavelength
    """

def spectral_factor_photovoltaic(airmass_absolute, precipitable_water, 
                                 aerosol_optical_depth, spectral_response,
                                 wavelengths=None):
    """
    Calculate spectral correction factor using custom spectral response.
    
    Parameters:
    - airmass_absolute: numeric, absolute airmass
    - precipitable_water: numeric, precipitable water in cm
    - aerosol_optical_depth: numeric, aerosol optical depth at 500nm
    - spectral_response: array-like, device spectral response
    - wavelengths: array-like, corresponding wavelengths in nm
    
    Returns:
    numeric, spectral correction factor
    """

def integrated_spectral_response(spectral_irradiance, spectral_response, 
                                 wavelengths=None):
    """
    Calculate integrated spectral response for given irradiance spectrum.
    
    Parameters:
    - spectral_irradiance: array-like, spectral irradiance (W/m²/nm)
    - spectral_response: array-like, device spectral response (A/W or QE)
    - wavelengths: array-like, wavelengths in nm
    
    Returns:
    numeric, integrated response
    """

def bandwidth_utilization_factor(spectral_response, reference_spectrum, 
                                 field_spectrum, wavelengths=None):
    """
    Calculate bandwidth utilization factor comparing field to reference conditions.
    
    Parameters:
    - spectral_response: array-like, device spectral response
    - reference_spectrum: array-like, reference spectral irradiance
    - field_spectrum: array-like, field spectral irradiance  
    - wavelengths: array-like, wavelengths in nm
    
    Returns:
    numeric, bandwidth utilization factor
    """

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

# Location and time
lat, lon = 37.8756, -122.2441  # Berkeley, CA
times = pd.date_range('2023-06-21 06:00', '2023-06-21 18:00', 
                     freq='H', tz='US/Pacific')

# Calculate solar position and atmospheric parameters
solar_pos = solarposition.get_solarposition(times, lat, lon)
airmass_rel = atmosphere.get_relative_airmass(solar_pos['zenith'])
airmass_abs = atmosphere.get_absolute_airmass(airmass_rel)

# Atmospheric conditions
precipitable_water = 2.5  # cm
aod500 = 0.15  # aerosol optical depth at 500nm

# Calculate spectral correction factors for different PV technologies
technologies = {
    'c-Si': {'type': 'csi', 'sapm_params': {'a0': 0.928, 'a1': 0.068, 'a2': -0.0077, 'a3': 0.0001, 'a4': -0.000002}},
    'CdTe': {'type': 'cdte', 'sapm_params': {'a0': 0.87, 'a1': 0.122, 'a2': -0.0047, 'a3': -0.000085, 'a4': 0.0000020}},
    'a-Si': {'type': 'asi', 'sapm_params': {'a0': 1.12, 'a1': -0.047, 'a2': -0.0085, 'a3': 0.000047, 'a4': 0.0000024}}
}

spectral_factors = {}

for tech_name, params in technologies.items():
    # SAPM spectral factor
    sapm_factor = spectrum.spectral_factor_sapm(airmass_abs, params['samp_params'])
    
    # First Solar factor (for CdTe)
    if tech_name == 'CdTe':
        fs_factor = spectrum.spectral_factor_firstsolar(
            precipitable_water, airmass_abs, 'cdte'
        )
    else:
        fs_factor = None
    
    # Caballero factor
    caballero_factor = spectrum.spectral_factor_caballero(
        precipitable_water, airmass_abs, aod500, module_type=params['type']
    )
    
    spectral_factors[tech_name] = {
        'sapm': samp_factor,
        'firstsolar': fs_factor,
        'caballero': caballero_factor
    }

# Plot spectral correction factors throughout the day
fig, axes = plt.subplots(2, 2, figsize=(15, 10))

# SAPM factors
for tech in technologies.keys():
    axes[0, 0].plot(times, spectral_factors[tech]['sapm'], 
                   label=tech, marker='o', linewidth=2)
axes[0, 0].set_title('SAPM Spectral Correction Factors')
axes[0, 0].set_ylabel('Spectral Factor')
axes[0, 0].legend()
axes[0, 0].grid(True)

# Caballero factors
for tech in technologies.keys():
    axes[0, 1].plot(times, spectral_factors[tech]['caballero'], 
                   label=tech, marker='s', linewidth=2)
axes[0, 1].set_title('Caballero Spectral Correction Factors')
axes[0, 1].set_ylabel('Spectral Factor')
axes[0, 1].legend()
axes[0, 1].grid(True)

# Airmass variation
axes[1, 0].plot(times, airmass_abs, 'ko-', linewidth=2)
axes[1, 0].set_title('Absolute Airmass')
axes[1, 0].set_ylabel('Airmass')
axes[1, 0].grid(True)

# Solar elevation
axes[1, 1].plot(times, 90 - solar_pos['zenith'], 'ro-', linewidth=2)
axes[1, 1].set_title('Solar Elevation Angle')
axes[1, 1].set_ylabel('Elevation (degrees)')
axes[1, 1].set_xlabel('Time')
axes[1, 1].grid(True)

plt.tight_layout()
plt.show()

# Print summary statistics
print("Daily Average Spectral Correction Factors:")
for tech in technologies.keys():
    sapm_avg = np.mean(spectral_factors[tech]['samp'])
    caballero_avg = np.mean(spectral_factors[tech]['caballero'])
    print(f"{tech:4s}: SAPM = {sapm_avg:.3f}, Caballero = {caballero_avg:.3f}")

import pvlib
from pvlib import spectrum, atmosphere, solarposition, clearsky
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

# Location and conditions
lat, lon = 40.0150, -105.2705  # Boulder, CO
elevation = 1655  # meters
time = pd.Timestamp('2023-06-21 12:00:00', tz='US/Mountain')

# Solar position
solar_pos = solarposition.get_solarposition(time, lat, lon)
zenith = solar_pos['zenith'].iloc[0]
azimuth = solar_pos['azimuth'].iloc[0]

# Atmospheric conditions
surface_pressure = atmosphere.alt2pres(elevation)
precipitable_water = 1.5  # cm
ozone = 0.31  # atm-cm
aod500 = 0.1  # aerosol optical depth
relative_humidity = 45  # percent

# Surface parameters
surface_tilt = 30  # degrees
surface_azimuth = 180  # degrees (south-facing)
ground_albedo = 0.2

# Calculate angle of incidence
from pvlib import irradiance
aoi = irradiance.aoi(surface_tilt, surface_azimuth, zenith, azimuth)

# Calculate spectral irradiance using SPECTRL2
wavelengths = np.arange(300, 1200, 5)  # 300-1200 nm in 5nm steps

spectral_result = spectrum.spectrl2(
    apparent_zenith=zenith,
    aoi=aoi,
    surface_tilt=surface_tilt,
    ground_albedo=ground_albedo,
    surface_pressure=surface_pressure/100,  # convert Pa to mbar
    relative_humidity=relative_humidity,
    precipitable_water=precipitable_water,
    ozone=ozone,
    aerosol_turbidity_500nm=aod500,
    dayofyear=172,  # June 21
    wavelengths=wavelengths
)

wavelengths_out = spectral_result[0]
poa_total = spectral_result[1]
poa_direct = spectral_result[2] 
poa_diffuse = spectral_result[3]
ghi_spectral = spectral_result[4]
dni_spectral = spectral_result[5]

# Get reference spectrum for comparison
reference_spectra = spectrum.get_reference_spectra()
am15g_spectrum = reference_spectra['wavelength'], reference_spectra['global_tilt_37']

# Plot spectral irradiance
fig, axes = plt.subplots(2, 2, figsize=(15, 10))

# POA spectral components
axes[0, 0].plot(wavelengths_out, poa_total, label='Total POA', linewidth=2)
axes[0, 0].plot(wavelengths_out, poa_direct, label='Direct POA', linewidth=2)
axes[0, 0].plot(wavelengths_out, poa_diffuse, label='Diffuse POA', linewidth=2)
axes[0, 0].set_xlabel('Wavelength (nm)')
axes[0, 0].set_ylabel('Spectral Irradiance (W/m²/nm)')
axes[0, 0].set_title('Plane-of-Array Spectral Irradiance Components')
axes[0, 0].legend()
axes[0, 0].grid(True)

# Horizontal spectral components  
axes[0, 1].plot(wavelengths_out, ghi_spectral, label='GHI', linewidth=2)
axes[0, 1].plot(wavelengths_out, dni_spectral, label='DNI', linewidth=2)
axes[0, 1].set_xlabel('Wavelength (nm)')
axes[0, 1].set_ylabel('Spectral Irradiance (W/m²/nm)')
axes[0, 1].set_title('Horizontal Spectral Irradiance')
axes[0, 1].legend()
axes[0, 1].grid(True)

# Comparison with AM1.5G reference
# Interpolate reference to same wavelengths for comparison
ref_interp = np.interp(wavelengths_out, am15g_spectrum[0], am15g_spectrum[1])
axes[1, 0].plot(wavelengths_out, poa_total, label='SPECTRL2 POA', linewidth=2)
axes[1, 0].plot(wavelengths_out, ref_interp, label='AM1.5G Reference', 
               linewidth=2, linestyle='--')
axes[1, 0].set_xlabel('Wavelength (nm)')
axes[1, 0].set_ylabel('Spectral Irradiance (W/m²/nm)')
axes[1, 0].set_title('Comparison with AM1.5G Reference')
axes[1, 0].legend()
axes[1, 0].grid(True)

# Spectral ratio
spectral_ratio = poa_total / ref_interp
axes[1, 1].plot(wavelengths_out, spectral_ratio, 'r-', linewidth=2)
axes[1, 1].set_xlabel('Wavelength (nm)')
axes[1, 1].set_ylabel('Ratio (Field/Reference)')
axes[1, 1].set_title('Spectral Ratio vs AM1.5G')
axes[1, 1].grid(True)
axes[1, 1].axhline(y=1.0, color='k', linestyle='--', alpha=0.5)

plt.tight_layout()
plt.show()

# Calculate integrated values
poa_integrated = np.trapz(poa_total, wavelengths_out)
ghi_integrated = np.trapz(ghi_spectral, wavelengths_out)
ref_integrated = np.trapz(ref_interp, wavelengths_out)

print(f"\nIntegrated Irradiance Values:")
print(f"POA Total: {poa_integrated:.1f} W/m²")
print(f"GHI: {ghi_integrated:.1f} W/m²")
print(f"AM1.5G Reference: {ref_integrated:.1f} W/m²")
print(f"POA/Reference Ratio: {poa_integrated/ref_integrated:.3f}")

import pvlib
from pvlib import spectrum
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

# Get example spectral response functions
example_sr = spectrum.get_example_spectral_response()

# Define wavelength range
wavelengths = np.arange(300, 1200, 2)  # 300-1200 nm in 2nm steps

# Get different spectral responses for analysis
technologies = ['si_c', 'si_pc', 'cdte', 'cigs', 'asi']
spectral_responses = {}

for tech in technologies:
    if tech in example_sr:
        # Interpolate to common wavelength grid
        sr_interp = np.interp(wavelengths, example_sr['wavelength'], 
                             example_sr[tech], left=0, right=0)
        spectral_responses[tech] = pd.Series(sr_interp, index=wavelengths)

# Get reference spectra
ref_spectra = spectrum.get_reference_spectra(wavelengths=wavelengths)
am15g = ref_spectra['global_tilt_37']
am15d = ref_spectra['direct_normal']

# Simulate different field conditions by modifying reference spectrum
# Condition 1: High airmass (red-shifted)
am_high = am15g * np.exp(-0.1 * (wavelengths - 600) / 300)  # Red shift

# Condition 2: High water vapor (IR absorption) 
water_absorption = 1 - 0.3 * np.exp(-((wavelengths - 940) / 40)**2)  # 940nm absorption
am_humid = am15g * water_absorption

# Condition 3: High aerosols (blue scattering)
blue_scattering = 1 - 0.2 * np.exp(-((wavelengths - 400) / 100)**2)  # Blue reduction
am_aerosol = am15g * blue_scattering

field_conditions = {
    'AM1.5G Reference': am15g,
    'High Airmass': am_high,  
    'High Humidity': am_humid,
    'High Aerosols': am_aerosol
}

# Calculate spectral mismatch factors for each technology and condition
mismatch_results = {}

for tech_name, sr in spectral_responses.items():
    mismatch_results[tech_name] = {}
    
    for condition_name, spectrum_field in field_conditions.items():
        if condition_name == 'AM1.5G Reference':
            mismatch_factor = 1.0  # Reference condition
        else:
            mismatch_factor = spectrum.calc_spectral_mismatch_field(
                sr, spectrum_field, am15g
            )
        
        mismatch_results[tech_name][condition_name] = mismatch_factor

# Create results dataframe
results_df = pd.DataFrame(mismatch_results).T
results_df = results_df.round(4)

print("Spectral Mismatch Factors:")
print(results_df)

# Plot spectral responses and field conditions
fig, axes = plt.subplots(2, 2, figsize=(15, 10))

# Spectral responses
for tech_name, sr in spectral_responses.items():
    axes[0, 0].plot(wavelengths, sr, label=tech_name, linewidth=2)
axes[0, 0].set_xlabel('Wavelength (nm)')
axes[0, 0].set_ylabel('Spectral Response')
axes[0, 0].set_title('PV Technology Spectral Response')
axes[0, 0].legend()
axes[0, 0].grid(True)

# Field spectra compared to reference
for condition_name, spectrum_field in field_conditions.items():
    if condition_name != 'AM1.5G Reference':
        ratio = spectrum_field / am15g
        axes[0, 1].plot(wavelengths, ratio, label=condition_name, linewidth=2)

axes[0, 1].set_xlabel('Wavelength (nm)')
axes[0, 1].set_ylabel('Ratio to AM1.5G')
axes[0, 1].set_title('Field Conditions vs AM1.5G Reference')
axes[0, 1].legend()
axes[0, 1].grid(True)
axes[0, 1].axhline(y=1.0, color='k', linestyle='--', alpha=0.5)

# Mismatch factors by technology
x_pos = np.arange(len(technologies))
width = 0.2
conditions = ['High Airmass', 'High Humidity', 'High Aerosols']
colors = ['red', 'blue', 'green']

for i, condition in enumerate(conditions):
    values = [mismatch_results[tech][condition] for tech in technologies]
    axes[1, 0].bar(x_pos + i*width, values, width, label=condition, 
                   color=colors[i], alpha=0.7)

axes[1, 0].set_xlabel('PV Technology')
axes[1, 0].set_ylabel('Spectral Mismatch Factor')
axes[1, 0].set_title('Spectral Mismatch by Technology and Condition')
axes[1, 0].set_xticks(x_pos + width)
axes[1, 0].set_xticklabels(technologies)
axes[1, 0].legend()
axes[1, 0].grid(True, axis='y')
axes[1, 0].axhline(y=1.0, color='k', linestyle='--', alpha=0.5)

# Technology ranking by spectral sensitivity
sensitivity = {}
for tech in technologies:
    # Calculate standard deviation of mismatch factors (excluding reference)
    values = [mismatch_results[tech][cond] for cond in conditions]
    sensitivity[tech] = np.std(values)

sorted_techs = sorted(sensitivity.items(), key=lambda x: x[1])
tech_names = [item[0] for item in sorted_techs]
sens_values = [item[1] for item in sorted_techs]

axes[1, 1].barh(tech_names, sens_values, color='orange', alpha=0.7)
axes[1, 1].set_xlabel('Spectral Sensitivity (Std Dev of Mismatch)')
axes[1, 1].set_title('Technology Ranking by Spectral Sensitivity')
axes[1, 1].grid(True, axis='x')

plt.tight_layout()
plt.show()

# Find best and worst performing technologies
print(f"\nMost spectrally stable technology: {sorted_techs[0][0]} (σ = {sorted_techs[0][1]:.4f})")
print(f"Most spectrally sensitive technology: {sorted_techs[-1][0]} (σ = {sorted_techs[-1][1]:.4f})")

# Calculate potential energy differences
print(f"\nPotential annual energy impact (assuming constant conditions):")
for tech in technologies:
    for condition in conditions:
        impact = (mismatch_results[tech][condition] - 1.0) * 100
        print(f"{tech} in {condition}: {impact:+.1f}%")

import pvlib
from pvlib import spectrum, atmosphere, clearsky, solarposition
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

# Location and time range
lat, lon = 35.0522, -106.5408  # Albuquerque, NM (high altitude, arid)
times = pd.date_range('2023-01-01', '2023-12-31', freq='MS')  # Monthly

# Calculate seasonal variations
seasonal_analysis = []

for time in times:
    # Solar position at solar noon
    solar_pos = solarposition.get_solarposition(time.replace(hour=12), lat, lon)
    zenith = solar_pos['zenith'].iloc[0]
    
    # Atmospheric parameters (seasonal variations)
    month = time.month
    
    # Seasonal precipitable water (higher in summer)
    pw_base = 1.0  # cm
    pw_seasonal = pw_base * (1 + 0.5 * np.sin(2 * np.pi * (month - 3) / 12))
    
    # Aerosol optical depth (higher in summer dust season)
    aod_base = 0.08
    aod_seasonal = aod_base * (1 + 0.3 * np.sin(2 * np.pi * (month - 6) / 12))
    
    # Calculate airmass
    airmass_rel = atmosphere.get_relative_airmass(zenith)
    airmass_abs = atmosphere.get_absolute_airmass(airmass_rel, pressure=85000)  # High altitude
    
    # Calculate clearsky index (simulate seasonal cloud variations)
    clearsky_model = clearsky.ineichen(zenith, airmass_abs, linke_turbidity=3.0)
    ghi_clear = clearsky_model['ghi']
    
    # Simulate measured GHI with seasonal cloud patterns
    cloud_factor = 0.9 - 0.2 * np.sin(2 * np.pi * (month - 1) / 12)  # More clouds in winter
    ghi_measured = ghi_clear * cloud_factor
    clearsky_idx = ghi_measured / ghi_clear if ghi_clear > 0 else 0
    
    # Calculate spectral correction factors for c-Si
    sapm_params = {'a0': 0.928, 'a1': 0.068, 'a2': -0.0077, 'a3': 0.0001, 'a4': -0.000002}
    
    spectral_sapm = spectrum.spectral_factor_sapm(airmass_abs, sapm_params)
    spectral_fs = spectrum.spectral_factor_firstsolar(pw_seasonal, airmass_abs, 'csi')
    spectral_caballero = spectrum.spectral_factor_caballero(
        pw_seasonal, airmass_abs, aod_seasonal, module_type='csi'
    )
    spectral_pvspec = spectrum.spectral_factor_pvspec(
        airmass_abs, clearsky_idx, module_type='monosi'
    )
    
    seasonal_analysis.append({
        'month': month,
        'zenith': zenith,
        'airmass': airmass_abs,
        'precipitable_water': pw_seasonal,
        'aod500': aod_seasonal,
        'clearsky_index': clearsky_idx,
        'spectral_sapm': spectral_sapm,
        'spectral_firstsolar': spectral_fs,
        'spectral_caballero': spectral_caballero,
        'spectral_pvspec': spectral_pvspec
    })

# Convert to DataFrame
seasonal_df = pd.DataFrame(seasonal_analysis)
seasonal_df.set_index('month', inplace=True)

# Plot seasonal variations
fig, axes = plt.subplots(3, 2, figsize=(15, 12))

# Atmospheric parameters
axes[0, 0].plot(seasonal_df.index, seasonal_df['airmass'], 'bo-', linewidth=2)
axes[0, 0].set_title('Seasonal Airmass Variation')
axes[0, 0].set_ylabel('Airmass')
axes[0, 0].grid(True)

axes[0, 1].plot(seasonal_df.index, seasonal_df['precipitable_water'], 'go-', linewidth=2, label='PW')
ax_twin = axes[0, 1].twinx()
ax_twin.plot(seasonal_df.index, seasonal_df['aod500'], 'ro-', linewidth=2, label='AOD')
axes[0, 1].set_title('Precipitable Water and AOD')
axes[0, 1].set_ylabel('Precipitable Water (cm)', color='g')
ax_twin.set_ylabel('AOD at 500nm', color='r')
axes[0, 1].grid(True)

# Clearsky conditions
axes[1, 0].plot(seasonal_df.index, seasonal_df['clearsky_index'], 'ko-', linewidth=2)
axes[1, 0].set_title('Seasonal Clearsky Index')
axes[1, 0].set_ylabel('Clearsky Index')
axes[1, 0].grid(True)

# Spectral correction factors
spectral_cols = ['spectral_sapm', 'spectral_firstsolar', 'spectral_caballero', 'spectral_pvspec']
colors = ['blue', 'red', 'green', 'orange']
labels = ['SAPM', 'First Solar', 'Caballero', 'PVSPEC']

for i, (col, color, label) in enumerate(zip(spectral_cols, colors, labels)):
    axes[1, 1].plot(seasonal_df.index, seasonal_df[col], 'o-', 
                   color=color, linewidth=2, label=label)

axes[1, 1].set_title('Spectral Correction Factors')
axes[1, 1].set_ylabel('Spectral Factor')
axes[1, 1].legend()
axes[1, 1].grid(True)

# Annual performance impact
reference_factor = 1.0
performance_impact = {}

for col in spectral_cols:
    # Calculate monthly energy impact
    monthly_impact = (seasonal_df[col] - reference_factor) * 100
    annual_avg_impact = monthly_impact.mean()
    performance_impact[col] = {
        'monthly': monthly_impact,
        'annual_avg': annual_avg_impact,
        'seasonal_range': monthly_impact.max() - monthly_impact.min()
    }

# Plot performance impacts
months = seasonal_df.index
for i, (col, color, label) in enumerate(zip(spectral_cols, colors, labels)):
    axes[2, 0].plot(months, performance_impact[col]['monthly'], 'o-', 
                   color=color, linewidth=2, label=label)

axes[2, 0].set_title('Monthly Performance Impact')
axes[2, 0].set_xlabel('Month')
axes[2, 0].set_ylabel('Performance Impact (%)')
axes[2, 0].legend()
axes[2, 0].grid(True)
axes[2, 0].axhline(y=0, color='k', linestyle='--', alpha=0.5)

# Summary statistics
model_names = [label.replace(' ', '_') for label in labels]
annual_avgs = [performance_impact[col]['annual_avg'] for col in spectral_cols]
seasonal_ranges = [performance_impact[col]['seasonal_range'] for col in spectral_cols]

x_pos = np.arange(len(model_names))
width = 0.35

bars1 = axes[2, 1].bar(x_pos - width/2, annual_avgs, width, label='Annual Avg', alpha=0.7)
bars2 = axes[2, 1].bar(x_pos + width/2, seasonal_ranges, width, label='Seasonal Range', alpha=0.7)

axes[2, 1].set_title('Annual Impact Summary')
axes[2, 1].set_xlabel('Spectral Model')
axes[2, 1].set_ylabel('Performance Impact (%)')
axes[2, 1].set_xticks(x_pos)
axes[2, 1].set_xticklabels(model_names, rotation=45)
axes[2, 1].legend()
axes[2, 1].grid(True, axis='y')

plt.tight_layout()
plt.show()

# Print summary statistics
print("\nAnnual Spectral Performance Analysis:")
print("="*50)
for i, col in enumerate(spectral_cols):
    label = labels[i]
    avg_impact = performance_impact[col]['annual_avg']
    range_impact = performance_impact[col]['seasonal_range'] 
    
    print(f"{label:12s}: Annual Avg = {avg_impact:+.2f}%, Seasonal Range = {range_impact:.2f}%")

print(f"\nLocation: {lat:.2f}°N, {lon:.2f}°W (High altitude, arid climate)")
print("Analysis shows seasonal variations in spectral performance due to:")
print("- Solar angle variations (airmass)")
print("- Seasonal atmospheric moisture")  
print("- Dust/aerosol loading")
print("- Cloud cover patterns")