Comparisons of Equation 9 in Passias & Kallbach, "Shading effects in rows of solar cell panels"

passias_masking_angle_comparison.py

import numpy as np
import matplotlib.pyplot as plt

def passias_masking_angle_original(slant_height, pitch, surface_tilt,
                                   reverse_minus):
    """
    implementation of Equation 9.  `correct_minus` is a boolean that
    controls whether an unclear minus sign is included or not.
    """
    B = slant_height
    C = pitch
    beta = np.radians(surface_tilt)
    sin_beta = np.sin(beta)
    cos_beta = np.cos(beta)
    f = B * sin_beta
    g = B * cos_beta
    arctan_f_C_g = np.arctan(f / (C - g))
    g_C_g = g / (C - g)
    sign = -1 if reverse_minus else 1
    psi_avg = sin_beta * C/B * (
        np.tan(beta) * arctan_f_C_g +
        np.log1p(g_C_g) -
        - sign * 0.5 * np.log1p(g_C_g**2) -
        arctan_f_C_g / (sin_beta * cos_beta * (1 + g_C_g))
    )
    return np.degrees(psi_avg)

def passias_masking_angle(gcr, surface_tilt):
    """ As implemented in the pvlib PR """
    beta = np.radians(surface_tilt)
    sin_beta = np.sin(beta)
    cos_beta = np.cos(beta)
    X = 1/gcr

    term1 = -X * sin_beta * np.log(np.abs(2 * X * cos_beta - (X**2 + 1))) / 2
    term2 = (X * cos_beta - 1) * np.arctan((X * cos_beta - 1) / (X * sin_beta))
    term3 = (1 - X * cos_beta) * np.arctan(cos_beta / sin_beta)
    term4 = X * np.log(X) * sin_beta

    psi_avg = term1 + term2 + term3 + term4

    return np.degrees(psi_avg)

def passias_masking_angle_numerical(slant_height, pitch, surface_tilt):
    """
    naive numerical implementation for comparison to the two analytic methods
    """
    B = slant_height
    C = pitch
    beta = np.radians(surface_tilt)
    psi = lambda z: np.arctan((B-z) * np.sin(beta) / (C - B*np.cos(beta) + z*np.cos(beta)))
    dz = B/1000
    z = np.arange(0, B, dz)
    psi_avg = np.sum([psi(zp) for zp in z], axis=0)*dz / B
    return np.degrees(psi_avg)


surface_tilt = np.arange(0.1, 90, 0.1)

fig, axes = plt.subplots(1, 4, figsize=(12, 4))
for k in [1, 1.5, 2, 2.5, 3, 4, 5, 7, 10]:
    gcr = 1/k

    # pvlib PR version
    psi = passias_masking_angle(gcr, surface_tilt)
    axes[0].plot(surface_tilt, psi, label='k={}'.format(k))

    # paper implementation, version 1
    pitch = 1.0
    slant_height = gcr * pitch
    psi = passias_masking_angle_original(slant_height, pitch, surface_tilt,
                                         True)
    axes[1].plot(surface_tilt, psi, label='k={}'.format(k))

    # paper implementation, version 2
    psi = passias_masking_angle_original(slant_height, pitch, surface_tilt,
                                         False)
    axes[2].plot(surface_tilt, psi, label='k={}'.format(k))

    # numerical integration
    psi = passias_masking_angle_numerical(slant_height, pitch, surface_tilt)
    axes[3].plot(surface_tilt, psi, label='k={}'.format(k))

axes[0].set_title('pvlib PR')
axes[1].set_title('passias v1')
axes[2].set_title('passias v2')
axes[3].set_title('numerical')

plt.show()