Incidence Angle Modifier Models

def ashrae(aoi, b=0.05):
    """
    ASHRAE incidence angle modifier model.
    
    Parameters:
    - aoi: numeric, angle of incidence in degrees
    - b: float, parameter to adjust IAM vs AOI (default 0.05)
    
    Returns:
    - iam: numeric, incident angle modifier (0-1)
    """

def physical(aoi, n=1.526, K=4.0, L=0.002, *, n_ar=None):
    """
    Physical IAM model based on Fresnel reflections and absorption.
    
    Parameters:
    - aoi: numeric, angle of incidence in degrees
    - n: numeric, effective refractive index (default 1.526 for glass)
    - K: numeric, glazing extinction coefficient in 1/m (default 4.0)
    - L: numeric, glazing thickness in meters (default 0.002)
    - n_ar: numeric, refractive index of AR coating (optional)
    
    Returns:
    - iam: numeric, incident angle modifier
    """

def martin_ruiz(aoi, a_r=0.16):
    """
    Martin and Ruiz incidence angle model.
    
    Parameters:
    - aoi: numeric, angle of incidence in degrees
    - a_r: numeric, angular losses coefficient (0.08-0.25 typical)
    
    Returns:
    - iam: numeric, incident angle modifier
    """

def martin_ruiz_diffuse(surface_tilt, a_r=0.16, c1=0.4244, c2=None):
    """
    Martin-Ruiz diffuse IAM factors for sky and ground irradiance.
    
    Parameters:
    - surface_tilt: numeric, surface tilt angle in degrees
    - a_r: numeric, angular losses coefficient
    - c1: float, first fitting parameter (default 0.4244)
    - c2: float, second fitting parameter (calculated if None)
    
    Returns:
    - iam_sky: numeric, IAM for sky diffuse irradiance
    - iam_ground: numeric, IAM for ground-reflected diffuse irradiance
    """

def sapm(aoi, module, upper=None):
    """
    SAPM incidence angle modifier model.
    
    Parameters:
    - aoi: numeric, angle of incidence in degrees
    - module: dict, module parameters containing B0-B5 coefficients
    - upper: float, upper limit on results (optional)
    
    Returns:
    - iam: numeric, SAPM angle of incidence loss coefficient (F2)
    """

def interp(aoi, theta_ref, iam_ref, method='linear', normalize=True):
    """
    IAM by interpolating reference measurements.
    
    Parameters:
    - aoi: numeric, angle of incidence in degrees
    - theta_ref: numeric, reference angles where IAM is known
    - iam_ref: numeric, IAM values at reference angles
    - method: str, interpolation method ('linear', 'quadratic', 'cubic')
    - normalize: bool, normalize to IAM=1 at normal incidence
    
    Returns:
    - iam: numeric, interpolated incident angle modifier
    """

def schlick(aoi):
    """
    Schlick approximation to Fresnel equations for IAM.
    
    Parameters:
    - aoi: numeric, angle of incidence in degrees
    
    Returns:
    - iam: numeric, incident angle modifier
    """

def schlick_diffuse(surface_tilt):
    """
    Analytical integration of Schlick model for diffuse IAM.
    
    Parameters:
    - surface_tilt: numeric, surface tilt angle in degrees
    
    Returns:
    - iam_sky: numeric, IAM for sky diffuse irradiance
    - iam_ground: numeric, IAM for ground-reflected diffuse irradiance
    """

def marion_diffuse(model, surface_tilt, **kwargs):
    """
    Diffuse IAM using Marion's integration method.
    
    Parameters:
    - model: str, IAM model name ('ashrae', 'physical', 'martin_ruiz', 'sapm', 'schlick')
    - surface_tilt: numeric, surface tilt angle in degrees
    - **kwargs: additional parameters for the IAM model
    
    Returns:
    - iam: dict with keys 'sky', 'horizon', 'ground' containing IAM values
    """

def marion_integrate(function, surface_tilt, region, num=None):
    """
    Integrate IAM function over solid angle region.
    
    Parameters:
    - function: callable, IAM function to integrate
    - surface_tilt: numeric, surface tilt angle in degrees
    - region: str, integration region ('sky', 'horizon', 'ground')
    - num: int, number of integration increments
    
    Returns:
    - iam: numeric, diffuse correction factor for specified region
    """

def convert(source_name, source_params, target_name, weight=None, 
           fix_n=True, xtol=None):
    """
    Convert parameters from one IAM model to another.
    
    Parameters:
    - source_name: str, source model ('ashrae', 'martin_ruiz', 'physical')
    - source_params: dict, parameters for source model
    - target_name: str, target model ('ashrae', 'martin_ruiz', 'physical')
    - weight: function, weighting function for residuals
    - fix_n: bool, fix refractive index when converting to physical model
    - xtol: float, optimization tolerance
    
    Returns:
    - dict, parameters for target model
    """

def fit(measured_aoi, measured_iam, model_name, weight=None, xtol=None):
    """
    Fit IAM model parameters to measured data.
    
    Parameters:
    - measured_aoi: array-like, measured angles of incidence in degrees
    - measured_iam: array-like, measured IAM values
    - model_name: str, model to fit ('ashrae', 'martin_ruiz', 'physical')
    - weight: function, weighting function for residuals
    - xtol: float, optimization tolerance
    
    Returns:
    - dict, fitted model parameters
    """

import pvlib
from pvlib import iam
import numpy as np
import matplotlib.pyplot as plt

# Angle of incidence range
aoi = np.linspace(0, 90, 91)

# Calculate IAM using different models
iam_ashrae = iam.ashrae(aoi, b=0.05)
iam_physical = iam.physical(aoi, n=1.526, K=4.0, L=0.002)
iam_martin_ruiz = iam.martin_ruiz(aoi, a_r=0.16)

# Plot comparison
plt.figure(figsize=(10, 6))
plt.plot(aoi, iam_ashrae, 'b-', label='ASHRAE (b=0.05)', linewidth=2)
plt.plot(aoi, iam_physical, 'r--', label='Physical (n=1.526, K=4, L=2mm)', linewidth=2)
plt.plot(aoi, iam_martin_ruiz, 'g:', label='Martin-Ruiz (a_r=0.16)', linewidth=2)
plt.xlabel('Angle of Incidence (degrees)')
plt.ylabel('Incidence Angle Modifier')
plt.title('IAM Model Comparison')
plt.legend()
plt.grid(True)
plt.xlim(0, 90)
plt.ylim(0, 1)
plt.show()

print("IAM Values at Key Angles:")
print("AOI (°)  ASHRAE  Physical  Martin-Ruiz")
for angle in [0, 30, 45, 60, 75, 90]:
    idx = angle
    print(f"{angle:5.0f}   {iam_ashrae[idx]:.3f}    {iam_physical[idx]:.3f}      {iam_martin_ruiz[idx]:.3f}")

import pvlib
from pvlib import iam
import numpy as np
import matplotlib.pyplot as plt

# Angle of incidence range
aoi = np.linspace(0, 90, 91)

# Compare with and without AR coating
iam_no_ar = iam.physical(aoi, n=1.526, K=4.0, L=0.002)
iam_with_ar = iam.physical(aoi, n=1.526, K=4.0, L=0.002, n_ar=1.29)

# Different glass types
iam_low_iron = iam.physical(aoi, n=1.526, K=2.0, L=0.002)  # Low-iron glass
iam_thick_glass = iam.physical(aoi, n=1.526, K=4.0, L=0.004)  # Thicker glass

plt.figure(figsize=(12, 8))

plt.subplot(2, 2, 1)
plt.plot(aoi, iam_no_ar, 'b-', label='No AR coating', linewidth=2)
plt.plot(aoi, iam_with_ar, 'r--', label='With AR coating (n=1.29)', linewidth=2)
plt.xlabel('Angle of Incidence (degrees)')
plt.ylabel('IAM')
plt.title('Effect of Anti-Reflective Coating')
plt.legend()
plt.grid(True)

plt.subplot(2, 2, 2)
plt.plot(aoi, iam_no_ar, 'b-', label='Standard glass (K=4)', linewidth=2)
plt.plot(aoi, iam_low_iron, 'g--', label='Low-iron glass (K=2)', linewidth=2)
plt.xlabel('Angle of Incidence (degrees)')
plt.ylabel('IAM')
plt.title('Effect of Glass Type')
plt.legend()
plt.grid(True)

plt.subplot(2, 2, 3)
plt.plot(aoi, iam_no_ar, 'b-', label='2mm glass', linewidth=2)
plt.plot(aoi, iam_thick_glass, 'm--', label='4mm glass', linewidth=2)
plt.xlabel('Angle of Incidence (degrees)')
plt.ylabel('IAM')
plt.title('Effect of Glass Thickness')
plt.legend()
plt.grid(True)

plt.subplot(2, 2, 4)
# Show transmission improvement with AR coating
improvement = iam_with_ar - iam_no_ar
plt.plot(aoi, improvement * 100, 'k-', linewidth=2)
plt.xlabel('Angle of Incidence (degrees)')
plt.ylabel('IAM Improvement (%)')
plt.title('AR Coating Improvement')
plt.grid(True)

plt.tight_layout()
plt.show()

print("AR Coating Benefits:")
print("AOI (°)  No AR   With AR  Improvement (%)")
for angle in [0, 30, 45, 60, 75]:
    no_ar = iam_no_ar[angle]
    with_ar = iam_with_ar[angle]
    improvement = (with_ar - no_ar) * 100
    print(f"{angle:5.0f}   {no_ar:.3f}   {with_ar:.3f}    {improvement:10.2f}")

import pvlib
from pvlib import iam, pvsystem
import numpy as np
import matplotlib.pyplot as plt

# Get module parameters from SAM database
modules = pvsystem.retrieve_sam('SandiaMod')

# Select a few different modules for comparison
module_names = [
    'Canadian_Solar_CS5P_220M___2009_',
    'SunPower_SPR_315E_WHT_D__2013_',
    'First_Solar_FS_380__2010_'
]

aoi = np.linspace(0, 90, 91)

plt.figure(figsize=(12, 8))

# Plot IAM curves for different modules
for i, module_name in enumerate(module_names):
    if module_name in modules:
        module = modules[module_name]
        iam_sapm = iam.sapm(aoi, module)
        
        # Clean up module name for display
        display_name = module_name.replace('_', ' ').replace('   ', ' - ')
        plt.plot(aoi, iam_sapm, linewidth=2, label=display_name)
        
        print(f"\nModule: {display_name}")
        print("SAPM IAM Coefficients:")
        for coeff in ['B0', 'B1', 'B2', 'B3', 'B4', 'B5']:
            if coeff in module:
                print(f"  {coeff}: {module[coeff]:.6f}")

plt.xlabel('Angle of Incidence (degrees)')
plt.ylabel('SAPM IAM (F2)')
plt.title('SAPM IAM Curves for Different Module Technologies')
plt.legend()
plt.grid(True)
plt.xlim(0, 90)
plt.ylim(0, 1.1)
plt.show()

# Show how SAPM can exceed 1.0 at some angles
print("\nSAMP IAM Values (can exceed 1.0):")
print("AOI (°)  " + "  ".join([name[:15] for name in module_names if name in modules]))
for angle in [0, 15, 30, 45, 60, 75, 90]:
    row = f"{angle:5.0f}   "
    for name in module_names:
        if name in modules:
            module = modules[name]
            iam_val = iam.sapm(angle, module)
            row += f"{iam_val:.3f}        "
    print(row)

import pvlib
from pvlib import iam
import numpy as np
import matplotlib.pyplot as plt

# Surface tilt angles
surface_tilts = np.arange(0, 91, 5)

# Calculate diffuse IAM using different approaches
results = {}

# Martin-Ruiz diffuse (analytical)
for tilt in surface_tilts:
    iam_sky, iam_ground = iam.martin_ruiz_diffuse(tilt, a_r=0.16)
    if 'martin_ruiz_sky' not in results:
        results['martin_ruiz_sky'] = []
        results['martin_ruiz_ground'] = []
    results['martin_ruiz_sky'].append(iam_sky)
    results['martin_ruiz_ground'].append(iam_ground)

# Schlick diffuse (analytical)  
for tilt in surface_tilts:
    iam_sky, iam_ground = iam.schlick_diffuse(tilt)
    if 'schlick_sky' not in results:
        results['schlick_sky'] = []
        results['schlick_ground'] = []
    results['schlick_sky'].append(iam_sky)
    results['schlick_ground'].append(iam_ground)

# Marion integration (numerical) - sample a few points
sample_tilts = [0, 20, 45, 70, 90]
marion_results = {}
for tilt in sample_tilts:
    # Physical model
    marion_physical = iam.marion_diffuse('physical', tilt)
    # ASHRAE model
    marion_ashrae = iam.marion_diffuse('ashrae', tilt, b=0.05)
    
    if tilt not in marion_results:
        marion_results[tilt] = {}
    marion_results[tilt]['physical'] = marion_physical
    marion_results[tilt]['ashrae'] = marion_ashrae

# Plot results
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))

# Sky diffuse
ax1.plot(surface_tilts, results['martin_ruiz_sky'], 'b-', 
         label='Martin-Ruiz', linewidth=2)
ax1.plot(surface_tilts, results['schlick_sky'], 'r--', 
         label='Schlick', linewidth=2)

# Add Marion integration points
for tilt in sample_tilts:
    if tilt in marion_results:
        ax1.plot(tilt, marion_results[tilt]['physical']['sky'], 
                'go', markersize=8, label='Marion Physical' if tilt == sample_tilts[0] else "")
        ax1.plot(tilt, marion_results[tilt]['ashrae']['sky'], 
                'mo', markersize=8, label='Marion ASHRAE' if tilt == sample_tilts[0] else "")

ax1.set_xlabel('Surface Tilt (degrees)')
ax1.set_ylabel('Sky Diffuse IAM')
ax1.set_title('Sky Diffuse IAM vs Surface Tilt')
ax1.legend()
ax1.grid(True)

# Ground diffuse
ax2.plot(surface_tilts, results['martin_ruiz_ground'], 'b-', 
         label='Martin-Ruiz', linewidth=2)
ax2.plot(surface_tilts, results['schlick_ground'], 'r--', 
         label='Schlick', linewidth=2)

# Add Marion integration points
for tilt in sample_tilts:
    if tilt in marion_results:
        ax2.plot(tilt, marion_results[tilt]['physical']['ground'], 
                'go', markersize=8, label='Marion Physical' if tilt == sample_tilts[0] else "")
        ax2.plot(tilt, marion_results[tilt]['ashrae']['ground'], 
                'mo', markersize=8, label='Marion ASHRAE' if tilt == sample_tilts[0] else "")

ax2.set_xlabel('Surface Tilt (degrees)')
ax2.set_ylabel('Ground Diffuse IAM')
ax2.set_title('Ground Diffuse IAM vs Surface Tilt')
ax2.legend()
ax2.grid(True)

plt.tight_layout()
plt.show()

print("Diffuse IAM Comparison at Key Tilt Angles:")
print("Tilt (°)  Martin-Ruiz Sky  Schlick Sky  Martin-Ruiz Ground  Schlick Ground")
for i, tilt in enumerate([0, 20, 45, 70, 90]):
    idx = tilt // 5  # Index in results arrays
    mr_sky = results['martin_ruiz_sky'][idx]
    schlick_sky = results['schlick_sky'][idx]
    mr_ground = results['martin_ruiz_ground'][idx]
    schlick_ground = results['schlick_ground'][idx]
    print(f"{tilt:6.0f}    {mr_sky:11.3f}    {schlick_sky:9.3f}    {mr_ground:13.3f}    {schlick_ground:11.3f}")

import pvlib
from pvlib import iam
import numpy as np
import matplotlib.pyplot as plt

# Generate synthetic measured data (simulating real measurements)
np.random.seed(42)
aoi_measured = np.array([0, 10, 20, 30, 40, 50, 60, 70, 80])

# "True" model (physical with known parameters)
true_params = {'n': 1.526, 'K': 4.0, 'L': 0.002}
iam_true = iam.physical(aoi_measured, **true_params)

# Add measurement noise
iam_measured = iam_true + np.random.normal(0, 0.01, len(iam_true))
iam_measured = np.clip(iam_measured, 0, 1)  # Keep physical bounds

print(f"Synthetic Measured Data:")
print("AOI (°)  True IAM  Measured IAM")
for i, angle in enumerate(aoi_measured):
    print(f"{angle:5.0f}   {iam_true[i]:.3f}     {iam_measured[i]:.3f}")

# Fit different models to the measured data
fitted_ashrae = iam.fit(aoi_measured, iam_measured, 'ashrae')
fitted_martin_ruiz = iam.fit(aoi_measured, iam_measured, 'martin_ruiz')
fitted_physical = iam.fit(aoi_measured, iam_measured, 'physical')

print(f"\nFitted Model Parameters:")
print(f"ASHRAE: {fitted_ashrae}")
print(f"Martin-Ruiz: {fitted_martin_ruiz}")
print(f"Physical: {fitted_physical}")
print(f"True Physical: {true_params}")

# Convert between models
ashrae_to_physical = iam.convert('ashrae', fitted_ashrae, 'physical')
martin_ruiz_to_ashrae = iam.convert('martin_ruiz', fitted_martin_ruiz, 'ashrae')

print(f"\nModel Conversions:")
print(f"ASHRAE → Physical: {ashrae_to_physical}")
print(f"Martin-Ruiz → ASHRAE: {martin_ruiz_to_ashrae}")

# Evaluate fitted models over full AOI range
aoi_full = np.linspace(0, 90, 91)
iam_fitted_ashrae = iam.ashrae(aoi_full, **fitted_ashrae)
iam_fitted_martin_ruiz = iam.martin_ruiz(aoi_full, **fitted_martin_ruiz)
iam_fitted_physical = iam.physical(aoi_full, **fitted_physical)
iam_true_full = iam.physical(aoi_full, **true_params)

# Plot results
plt.figure(figsize=(12, 8))

plt.subplot(2, 2, 1)
plt.plot(aoi_full, iam_true_full, 'k-', label='True Physical', linewidth=3)
plt.plot(aoi_measured, iam_measured, 'ro', label='Measured Data', markersize=8)
plt.plot(aoi_full, iam_fitted_ashrae, 'b--', label='Fitted ASHRAE', linewidth=2)
plt.plot(aoi_full, iam_fitted_martin_ruiz, 'g:', label='Fitted Martin-Ruiz', linewidth=2)
plt.plot(aoi_full, iam_fitted_physical, 'r-.', label='Fitted Physical', linewidth=2)
plt.xlabel('Angle of Incidence (degrees)')
plt.ylabel('IAM')
plt.title('Model Fitting Comparison')
plt.legend()
plt.grid(True)

plt.subplot(2, 2, 2)
# Show residuals
residuals_ashrae = iam_true_full - iam_fitted_ashrae
residuals_martin_ruiz = iam_true_full - iam_fitted_martin_ruiz
residuals_physical = iam_true_full - iam_fitted_physical

plt.plot(aoi_full, residuals_ashrae * 100, 'b--', label='ASHRAE', linewidth=2)
plt.plot(aoi_full, residuals_martin_ruiz * 100, 'g:', label='Martin-Ruiz', linewidth=2)
plt.plot(aoi_full, residuals_physical * 100, 'r-.', label='Physical', linewidth=2)
plt.xlabel('Angle of Incidence (degrees)')
plt.ylabel('Residual (%)')
plt.title('Fitting Residuals (True - Fitted)')
plt.legend()
plt.grid(True)

plt.subplot(2, 2, 3)
# Compare converted models
iam_converted_physical = iam.physical(aoi_full, **ashrae_to_physical)
iam_converted_ashrae = iam.ashrae(aoi_full, **martin_ruiz_to_ashrae)

plt.plot(aoi_full, iam_fitted_ashrae, 'b-', label='Original ASHRAE', linewidth=2)
plt.plot(aoi_full, iam_converted_physical, 'b--', label='ASHRAE→Physical', linewidth=2)
plt.plot(aoi_full, iam_fitted_martin_ruiz, 'g-', label='Original Martin-Ruiz', linewidth=2)
plt.plot(aoi_full, iam_converted_ashrae, 'g--', label='Martin-Ruiz→ASHRAE', linewidth=2)
plt.xlabel('Angle of Incidence (degrees)')
plt.ylabel('IAM')
plt.title('Model Conversions')
plt.legend()
plt.grid(True)

plt.subplot(2, 2, 4)
# Error metrics
models = ['ASHRAE', 'Martin-Ruiz', 'Physical']
rmse_values = []
mae_values = []

for residuals in [residuals_ashrae, residuals_martin_ruiz, residuals_physical]:
    rmse = np.sqrt(np.mean(residuals**2)) * 100
    mae = np.mean(np.abs(residuals)) * 100
    rmse_values.append(rmse)
    mae_values.append(mae)

x = np.arange(len(models))
width = 0.35

plt.bar(x - width/2, rmse_values, width, label='RMSE', alpha=0.8)
plt.bar(x + width/2, mae_values, width, label='MAE', alpha=0.8)
plt.xlabel('Model')
plt.ylabel('Error (%)')
plt.title('Fitting Error Metrics')
plt.xticks(x, models)
plt.legend()
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print(f"\nFitting Error Summary:")
print("Model         RMSE (%)  MAE (%)")
for i, model in enumerate(models):
    print(f"{model:12s}  {rmse_values[i]:6.3f}   {mae_values[i]:5.3f}")

import pvlib
from pvlib import iam
import numpy as np
import matplotlib.pyplot as plt

# Reference data points (typically from manufacturer measurements)
theta_ref = np.array([0, 10, 20, 30, 40, 50, 60, 70, 80, 90])
iam_ref = np.array([1.000, 1.000, 0.999, 0.996, 0.992, 0.985, 0.974, 0.955, 0.921, 0.000])

# AOI range for interpolation
aoi = np.linspace(0, 90, 181)

# Different interpolation methods
iam_linear = iam.interp(aoi, theta_ref, iam_ref, method='linear')
iam_quadratic = iam.interp(aoi, theta_ref, iam_ref, method='quadratic')
iam_cubic = iam.interp(aoi, theta_ref, iam_ref, method='cubic')

# Compare with parametric models fitted to same reference data
fitted_ashrae = iam.fit(theta_ref, iam_ref, 'ashrae')
fitted_martin_ruiz = iam.fit(theta_ref, iam_ref, 'martin_ruiz')

iam_ashrae_fit = iam.ashrae(aoi, **fitted_ashrae)
iam_martin_ruiz_fit = iam.martin_ruiz(aoi, **fitted_martin_ruiz)

# Plot comparison
plt.figure(figsize=(12, 8))

plt.subplot(2, 1, 1)
plt.plot(theta_ref, iam_ref, 'ko', markersize=8, label='Reference Data')
plt.plot(aoi, iam_linear, 'b-', label='Linear Interpolation', linewidth=2)
plt.plot(aoi, iam_quadratic, 'r--', label='Quadratic Interpolation', linewidth=2)
plt.plot(aoi, iam_cubic, 'g:', label='Cubic Interpolation', linewidth=2)
plt.xlabel('Angle of Incidence (degrees)')
plt.ylabel('IAM')
plt.title('Interpolation Methods Comparison')
plt.legend()
plt.grid(True)
plt.xlim(0, 90)
plt.ylim(0.8, 1.05)

plt.subplot(2, 1, 2)
plt.plot(theta_ref, iam_ref, 'ko', markersize=8, label='Reference Data')
plt.plot(aoi, iam_cubic, 'g-', label='Cubic Interpolation', linewidth=2)
plt.plot(aoi, iam_ashrae_fit, 'm--', label=f'ASHRAE Fit (b={fitted_ashrae["b"]:.3f})', linewidth=2)
plt.plot(aoi, iam_martin_ruiz_fit, 'c:', label=f'Martin-Ruiz Fit (a_r={fitted_martin_ruiz["a_r"]:.3f})', linewidth=2)
plt.xlabel('Angle of Incidence (degrees)')
plt.ylabel('IAM')
plt.title('Interpolation vs Parametric Model Fits')
plt.legend()
plt.grid(True)
plt.xlim(0, 90)
plt.ylim(0.8, 1.05)

plt.tight_layout()
plt.show()

print("Interpolation vs Parametric Models at Key Angles:")
print("AOI (°)  Linear  Quadratic  Cubic   ASHRAE  Martin-Ruiz")
for angle in [0, 15, 30, 45, 60, 75, 90]:
    lin = np.interp(angle, aoi, iam_linear)
    quad = np.interp(angle, aoi, iam_quadratic)
    cub = np.interp(angle, aoi, iam_cubic)
    ash = np.interp(angle, aoi, iam_ashrae_fit)
    mr = np.interp(angle, aoi, iam_martin_ruiz_fit)
    print(f"{angle:5.0f}   {lin:.3f}    {quad:.3f}     {cub:.3f}   {ash:.3f}     {mr:.3f}")

# Calculate RMS differences from cubic interpolation (reference)
rms_linear = np.sqrt(np.mean((iam_linear - iam_cubic)**2)) * 100
rms_quadratic = np.sqrt(np.mean((iam_quadratic - iam_cubic)**2)) * 100
rms_ashrae = np.sqrt(np.mean((iam_ashrae_fit - iam_cubic)**2)) * 100
rms_martin_ruiz = np.sqrt(np.mean((iam_martin_ruiz_fit - iam_cubic)**2)) * 100

print(f"\nRMS Differences from Cubic Interpolation:")
print(f"Linear:      {rms_linear:.4f}%")
print(f"Quadratic:   {rms_quadratic:.4f}%")
print(f"ASHRAE:      {rms_ashrae:.4f}%")
print(f"Martin-Ruiz: {rms_martin_ruiz:.4f}%")