"""
The ``shading`` module contains functions that model module shading and the
associated effects on PV module output
"""
import numpy as np
import pandas as pd
from pvlib.tools import sind, cosd
from pvlib._deprecation import renamed_kwarg_warning
[docs]
def ground_angle(surface_tilt, gcr, slant_height):
    """
    Angle from horizontal of the line from a point on the row slant length
    to the bottom of the facing row.
    The angles are clockwise from horizontal, rather than the usual
    counterclockwise direction.
    Parameters
    ----------
    surface_tilt : numeric
        Surface tilt angle in degrees from horizontal, e.g., surface facing up
        = 0, surface facing horizon = 90. [degree]
    gcr : float
        ground coverage ratio, ratio of row slant length to row spacing.
        [unitless]
    slant_height : numeric
        The distance up the module's slant height to evaluate the ground
        angle, as a fraction [0-1] of the module slant height [unitless].
    Returns
    -------
    psi : numeric
        Angle [degree].
    """
    #  : \\            \
    #  :  \\            \
    #  :   \\            \
    #  :    \\            \  facing row
    #  :     \\.___________\
    #  :       \  ^*-.  psi \
    #  :        \  x   *-.   \
    #  :         \  v      *-.\
    #  :          \<-----P---->\
    x1 = gcr * slant_height * sind(surface_tilt)
    x2 = gcr * slant_height * cosd(surface_tilt) + 1
    psi = np.arctan2(x1, x2)  # do this before rad2deg because it handles 0 / 0
    return np.rad2deg(psi) 
[docs]
def masking_angle(surface_tilt, gcr, slant_height):
    """
    The elevation angle below which diffuse irradiance is blocked.
    The ``height`` parameter determines how far up the module's surface to
    evaluate the masking angle.  The lower the point, the steeper the masking
    angle [1]_.  SAM uses a "worst-case" approach where the masking angle
    is calculated for the bottom of the array (i.e. ``slant_height=0``) [2]_.
    Parameters
    ----------
    surface_tilt : numeric
        Panel tilt from horizontal [degrees].
    gcr : float
        The ground coverage ratio of the array [unitless].
    slant_height : numeric
        The distance up the module's slant height to evaluate the masking
        angle, as a fraction [0-1] of the module slant height [unitless].
    Returns
    -------
    mask_angle : numeric
        Angle from horizontal where diffuse light is blocked by the
        preceding row [degrees].
    See Also
    --------
    masking_angle_passias
    sky_diffuse_passias
    References
    ----------
    .. [1] D. Passias and B. Källbäck, "Shading effects in rows of solar cell
       panels", Solar Cells, Volume 11, Pages 281-291.  1984.
       :doi:`10.1016/0379-6787(84)90017-6`
    .. [2] Gilman, P. et al., (2018). "SAM Photovoltaic Model Technical
       Reference Update", NREL Technical Report NREL/TP-6A20-67399.
       Available at https://www.nrel.gov/docs/fy18osti/67399.pdf
    """
    # The original equation (8 in [1]) requires pitch and collector width,
    # but it's easy to non-dimensionalize it to make it a function of GCR
    # by factoring out B from the argument to arctan.
    numerator = gcr * (1 - slant_height) * sind(surface_tilt)
    denominator = 1 - gcr * (1 - slant_height) * cosd(surface_tilt)
    phi = np.arctan(numerator / denominator)
    return np.degrees(phi) 
[docs]
def masking_angle_passias(surface_tilt, gcr):
    r"""
    The average masking angle over the slant height of a row.
    The masking angle is the angle from horizontal where the sky dome is
    blocked by the row in front. The masking angle is larger near the lower
    edge of a row than near the upper edge. This function calculates the
    average masking angle as described in [1]_.
    Parameters
    ----------
    surface_tilt : numeric
        Panel tilt from horizontal [degrees].
    gcr : float
        The ground coverage ratio of the array [unitless].
    Returns
    ----------
    mask_angle : numeric
        Average angle from horizontal where diffuse light is blocked by the
        preceding row [degrees].
    See Also
    --------
    masking_angle
    sky_diffuse_passias
    Notes
    -----
    The pvlib-python authors believe that Eqn. 9 in [1]_ is incorrect.
    Here we use an independent equation.  First, Eqn. 8 is non-dimensionalized
    (recasting in terms of GCR):
    .. math::
        \psi(z') = \arctan \left [
            \frac{(1 - z') \sin \beta}
                 {\mathrm{GCR}^{-1} + (z' - 1) \cos \beta}
        \right ]
    Where :math:`GCR = B/C` and :math:`z' = z/B`. The average masking angle
    :math:`\overline{\psi} = \int_0^1 \psi(z') \mathrm{d}z'` is then
    evaluated symbolically using Maxima (using :math:`X = 1/\mathrm{GCR}`):
    .. code-block:: none
        load(scifac)    /* for the gcfac function */
        assume(X>0, cos(beta)>0, cos(beta)-X<0);   /* X is 1/GCR */
        gcfac(integrate(atan((1-z)*sin(beta)/(X+(z-1)*cos(beta))), z, 0, 1))
    This yields the equation implemented by this function:
    .. math::
        \overline{\psi} = \
            &-\frac{X}{2} \sin\beta \log | 2 X \cos\beta - (X^2 + 1)| \\
            &+ (X \cos\beta - 1) \arctan \frac{X \cos\beta - 1}{X \sin\beta} \\
            &+ (1 - X \cos\beta) \arctan \frac{\cos\beta}{\sin\beta} \\
            &+ X \log X \sin\beta
    The pvlib-python authors have validated this equation against numerical
    integration of :math:`\overline{\psi} = \int_0^1 \psi(z') \mathrm{d}z'`.
    References
    ----------
    .. [1] D. Passias and B. Källbäck, "Shading effects in rows of solar cell
       panels", Solar Cells, Volume 11, Pages 281-291.  1984.
       :doi:`10.1016/0379-6787(84)90017-6`
    """
    # wrap it in an array so that division by zero is handled well
    beta = np.radians(np.array(surface_tilt))
    sin_b = np.sin(beta)
    cos_b = np.cos(beta)
    X = 1/gcr
    with np.errstate(divide='ignore', invalid='ignore'):  # ignore beta=0
        term1 = -X * sin_b * np.log(np.abs(2 * X * cos_b - (X**2 + 1))) / 2
        term2 = (X * cos_b - 1) * np.arctan((X * cos_b - 1) / (X * sin_b))
        term3 = (1 - X * cos_b) * np.arctan(cos_b / sin_b)
        term4 = X * np.log(X) * sin_b
    psi_avg = term1 + term2 + term3 + term4
    # when beta=0, divide by zero makes psi_avg NaN.  replace with 0:
    psi_avg = np.where(np.isfinite(psi_avg), psi_avg, 0)
    if isinstance(surface_tilt, pd.Series):
        psi_avg = pd.Series(psi_avg, index=surface_tilt.index)
    return np.degrees(psi_avg) 
[docs]
def sky_diffuse_passias(masking_angle):
    r"""
    The diffuse irradiance loss caused by row-to-row sky diffuse shading.
    Even when the sun is high in the sky, a row's view of the sky dome will
    be partially blocked by the row in front. This causes a reduction in the
    diffuse irradiance incident on the module. The reduction depends on the
    masking angle, the elevation angle from a point on the shaded module to
    the top of the shading row. In [1]_ the masking angle is calculated as
    the average across the module height. SAM assumes the "worst-case" loss
    where the masking angle is calculated for the bottom of the array [2]_.
    This function, as in [1]_, makes the assumption that sky diffuse
    irradiance is isotropic.
    Parameters
    ----------
    masking_angle : numeric
        The elevation angle below which diffuse irradiance is blocked
        [degrees].
    Returns
    -------
    derate : numeric
        The fraction [0-1] of blocked sky diffuse irradiance.
    See Also
    --------
    masking_angle
    masking_angle_passias
    References
    ----------
    .. [1] D. Passias and B. Källbäck, "Shading effects in rows of solar cell
       panels", Solar Cells, Volume 11, Pages 281-291.  1984.
       :doi:`10.1016/0379-6787(84)90017-6`
    .. [2] Gilman, P. et al., (2018). "SAM Photovoltaic Model Technical
       Reference Update", NREL Technical Report NREL/TP-6A20-67399.
       Available at https://www.nrel.gov/docs/fy18osti/67399.pdf
    """
    return 1 - cosd(masking_angle/2)**2 
[docs]
@renamed_kwarg_warning(
    since="0.13.1",
    old_param_name="axis_tilt",
    new_param_name="axis_slope",
)
def projected_solar_zenith_angle(solar_zenith, solar_azimuth,
                                 axis_slope, axis_azimuth):
    r"""
    Calculate projected solar zenith angle in degrees.
    This solar zenith angle is projected onto the plane whose normal vector is
    defined by ``axis_slope`` and ``axis_azimuth``. The normal vector is in the
    direction of ``axis_azimuth`` (clockwise from north) and tilted from
    horizontal by ``axis_slope``. See Figure 5 in [1]_:
    .. figure:: ../../_images/Anderson_Mikofski_2020_Fig5.jpg
       :alt: Wire diagram of coordinates systems to obtain the projected angle.
       :align: center
       :scale: 50 %
       Fig. 5, [1]_: Solar coordinates projection onto tracker rotation plane.
    Parameters
    ----------
    solar_zenith : numeric
        Sun's apparent zenith in degrees.
    solar_azimuth : numeric
        Sun's azimuth in degrees.
    axis_slope : numeric
        Axis tilt angle in degrees. From horizontal plane to array plane.
        .. versionchanged:: 0.13.1
            Renamed from ``axis_tilt`` to ``axis_slope``
    axis_azimuth : numeric
        Axis azimuth angle in degrees.
        North = 0°; East = 90°; South = 180°; West = 270°
    Returns
    -------
    Projected_solar_zenith : numeric
        In degrees.
    Notes
    -----
    This projection has a variety of applications in PV. For example:
    - Projecting the sun's position onto the plane perpendicular to
      the axis of a single-axis tracker (i.e. the plane
      whose normal vector coincides with the tracker torque tube)
      yields the tracker rotation angle that maximizes direct irradiance
      capture. This tracking strategy is called *true-tracking*. Learn more
      about tracking in
      :ref:`sphx_glr_gallery_solar-tracking_plot_single_axis_tracking.py`.
    - Self-shading in large PV arrays is often modeled by assuming
      a simplified 2-D array geometry where the sun's position is
      projected onto the plane perpendicular to the PV rows.
      The projected zenith angle is then used for calculations
      regarding row-to-row shading.
    Examples
    --------
    Calculate the ideal true-tracking angle for a horizontal north-south
    single-axis tracker:
    >>> rotation = projected_solar_zenith_angle(solar_zenith, solar_azimuth,
    >>>                                         axis_slope=0, axis_azimuth=180)
    Calculate the projected zenith angle in a south-facing fixed tilt array
    (note: the ``axis_azimuth`` of a fixed-tilt row points along the length
    of the row):
    >>> psza = projected_solar_zenith_angle(solar_zenith, solar_azimuth,
    >>>                                     axis_slope=0, axis_azimuth=90)
    References
    ----------
    .. [1] K. Anderson and M. Mikofski, 'Slope-Aware Backtracking for
       Single-Axis Trackers', National Renewable Energy Lab. (NREL), Golden,
       CO (United States);
       NREL/TP-5K00-76626, Jul. 2020. :doi:`10.2172/1660126`.
    See Also
    --------
    pvlib.solarposition.get_solarposition
    """
    # Assume the tracker reference frame is right-handed. Positive y-axis is
    # oriented along tracking axis; from north, the y-axis is rotated clockwise
    # by the axis azimuth and tilted from horizontal by the axis tilt. The
    # positive x-axis is 90 deg clockwise from the y-axis and parallel to
    # horizontal (e.g., if the y-axis is south, the x-axis is west); the
    # positive z-axis is normal to the x and y axes, pointed upward.
    # Since elevation = 90 - zenith, sin(90-x) = cos(x) & cos(90-x) = sin(x):
    # Notation from [1], modified to use zenith instead of elevation
    # cos(elevation) = sin(zenith) and sin(elevation) = cos(zenith)
    # Avoid recalculating these values
    sind_solar_zenith = sind(solar_zenith)
    cosd_axis_azimuth = cosd(axis_azimuth)
    sind_axis_azimuth = sind(axis_azimuth)
    sind_axis_slope = sind(axis_slope)
    # Sun's x, y, z coords
    sx = sind_solar_zenith * sind(solar_azimuth)
    sy = sind_solar_zenith * cosd(solar_azimuth)
    sz = cosd(solar_zenith)
    # Eq. (4); sx', sz' values from sun coordinates projected onto surface
    sx_prime = sx * cosd_axis_azimuth - sy * sind_axis_azimuth
    sz_prime = (
        sx * sind_axis_azimuth * sind_axis_slope
        + sy * sind_axis_slope * cosd_axis_azimuth
        + sz * cosd(axis_slope)
    )
    # Eq. (5); angle between sun's beam and surface
    theta_T = np.degrees(np.arctan2(sx_prime, sz_prime))
    return theta_T 
[docs]
@renamed_kwarg_warning(
    since="0.13.1",
    old_param_name="axis_tilt",
    new_param_name="axis_slope",
)
@renamed_kwarg_warning(
    since="0.13.1",
    old_param_name="cross_axis_tilt",
    new_param_name="cross_axis_slope",
)
def shaded_fraction1d(
    solar_zenith,
    solar_azimuth,
    axis_azimuth,
    shaded_row_rotation,
    *,
    collector_width,
    pitch,
    axis_slope=0,
    surface_to_axis_offset=0,
    cross_axis_slope=0,
    shading_row_rotation=None,
):
    r"""
    Shaded fraction in the vertical dimension of tilted rows, or perpendicular
    to the axis of horizontal rows.
    If ``shading_row_rotation`` isn't provided, it is assumed that
    both the shaded row and the shading row (the one blocking the
    direct beam) have the same rotation and azimuth values.
    .. warning::
        The function assumes that the roles of the shaded and shading rows
        remain the same during the day. In the case where the shading and
        shaded rows change throughout the day, e.g. a N-S single-axis tracker,
        the inputs must be switched depending on the sign of the projected
        solar zenith angle. See the Examples section below.
    .. versionadded:: 0.11.0
    Parameters
    ----------
    solar_zenith : numeric
        Solar zenith angle, in degrees.
    solar_azimuth : numeric
        Solar azimuth angle, in degrees.
    axis_azimuth : numeric
        Axis azimuth of the rotation axis of a tracker, in degrees.
        Fixed-tilt arrays can be considered as a particular case of a tracker.
        North=0º, South=180º, East=90º, West=270º.
    shaded_row_rotation : numeric
        Right-handed rotation of the row receiving the shade, with respect
        to ``axis_azimuth``. In degrees :math:`^{\circ}`.
    collector_width : numeric
        Vertical length of a tilted row. The returned ``shaded_fraction``
        is the ratio of the shadow over this value.
    pitch : numeric
        Axis-to-axis horizontal spacing of the row.
    axis_slope : numeric, default 0
        Tilt of the rows axis from horizontal. In degrees :math:`^{\circ}`.
        .. versionchanged:: 0.13.1
            Renamed from ``axis_tilt`` to ``axis_slope``
    surface_to_axis_offset : numeric, default 0
        Distance between the rotating axis and the collector surface.
        May be used to account for a torque tube offset.
    cross_axis_slope : numeric, default 0
        Angle of the plane containing the rows' axes relative to
        horizontal. Right-handed rotation with respect to the rows axes.
        See :term:`cross_axis_slope`. In degrees :math:`^{\circ}`.
        .. versionchanged:: 0.13.1
            Renamed from ``cross_axis_tilt`` to ``cross_axis_slope``
    shading_row_rotation : numeric, optional
        Right-handed rotation of the row casting the shadow, with respect
        to the row axis. In degrees :math:`^{\circ}`.
    Returns
    -------
    shaded_fraction : numeric
        The fraction of the collector width shaded by an adjacent row. A
        value of 1 is completely shaded and 0 is no shade.
    Notes
    -----
    All length parameters must have the same units.
    Parameters are defined as follow:
    .. figure:: ../../_images/Anderson_Jensen_2024_Fig3.png
        :alt: Diagram showing the two rows and the parameters of the model.
        Figure 3 of [1]_. See correspondence between this nomenclature and the
        function parameters in the table below.
    +------------------+----------------------------+---------------------+
    | Symbol           | Parameter                  | Units               |
    +==================+============================+=====================+
    | :math:`\theta_1` | ``shading_row_rotation``   |                     |
    +------------------+----------------------------+                     |
    | :math:`\theta_2` | ``shaded_row_rotation``    | Degrees             |
    +------------------+----------------------------+ :math:`^{\circ}`    |
    | :math:`\beta_c`  | ``cross_axis_slope``       |                     |
    +------------------+----------------------------+---------------------+
    | :math:`p`        | ``pitch``                  | Any consistent      |
    +------------------+----------------------------+ length unit across  |
    | :math:`\ell`     | ``collector_width``        | all these           |
    +------------------+----------------------------+ parameters, e.g.    |
    | :math:`z_0`      | ``surface_to_axis_offset`` | :math:`m`.          |
    +------------------+----------------------------+---------------------+
    | :math:`f_s`      | Return value               | Dimensionless       |
    +------------------+----------------------------+---------------------+
    Examples
    --------
    **Fixed-tilt south-facing array on flat terrain**
    Tilted row with a pitch of 3 m, a collector width of
    2 m, and row rotations of 30°. In the morning.
    >>> shaded_fraction1d(solar_zenith=80, solar_azimuth=135,
    ...     axis_azimuth=90, shaded_row_rotation=30, shading_row_rotation=30,
    ...     collector_width=2, pitch=3, axis_slope=0,
    ...     surface_to_axis_offset=0.05, cross_axis_slope=0)
    0.47755694708090535
    **Fixed-tilt north-facing array on sloped terrain**
    Tilted row with a pitch of 4 m, a collector width of
    2.5 m, and row rotations of 50° for the shaded
    row and 30° for the shading row. The rows are on a
    10° slope, where their axis is on the most inclined
    direction (zero cross-axis slope). Shaded in the morning.
    >>> shaded_fraction1d(solar_zenith=80, solar_azimuth=75.5,
    ...     axis_azimuth=270, shaded_row_rotation=50, shading_row_rotation=30,
    ...     collector_width=2.5, pitch=4, axis_slope=10,
    ...     surface_to_axis_offset=0.05, cross_axis_slope=0)
    0.793244836197256
    **N-S single-axis tracker on sloped terrain**
    Horizontal trackers with a pitch of 3 m, a collector width of
    1.4 m, and tracker rotations of 30° pointing east,
    in the morning. Terrain slope is 7° west-east (east-most
    tracker is higher than the west-most tracker).
    >>> shaded_fraction1d(solar_zenith=80, solar_azimuth=90, axis_azimuth=180,
    ...     shaded_row_rotation=-30, collector_width=1.4, pitch=3,
    ...     axis_slope=0, surface_to_axis_offset=0.10, cross_axis_slope=7)
    0.8242176864434579
    Note the previous example only is valid for the shaded fraction of the
    west-most tracker in the morning, and assuming it is the
    shaded tracker during all the day is incorrect.
    During the afternoon, it is the one casting the shadow onto the
    east-most tracker.
    To calculate the shaded fraction for the east-most
    tracker, you must input the corresponding ``shaded_row_rotation``
    in the afternoon.
    >>> shaded_fraction1d(solar_zenith=80, solar_azimuth=270, axis_azimuth=180,
    ...     shaded_row_rotation=30, collector_width=1.4, pitch=3, axis_slope=0,
    ...     surface_to_axis_offset=0.10, cross_axis_slope=7)
    0.018002567182254348
    You must switch the input/output depending on the
    sign of the projected solar zenith angle. See
    :py:func:`~pvlib.shading.projected_solar_zenith_angle` and the example
    :ref:`sphx_glr_gallery_shading_plot_shaded_fraction1d_ns_hsat_example.py`
    See also
    --------
    pvlib.shading.projected_solar_zenith_angle
    References
    ----------
    .. [1] Kevin S. Anderson, Adam R. Jensen; Shaded fraction and backtracking
        in single-axis trackers on rolling terrain. J. Renewable Sustainable
        Energy 1 March 2024; 16 (2): 023504. :doi:`10.1063/5.0202220`
    """
    # For nomenclature you may refer to [1].
    # rotation of row casting the shadow defaults to shaded row's one
    if shading_row_rotation is None:
        shading_row_rotation = shaded_row_rotation
    # projected solar zenith angle
    projected_solar_zenith = projected_solar_zenith_angle(
        solar_zenith,
        solar_azimuth,
        axis_slope,
        axis_azimuth,
    )
    # calculate repeated elements
    thetas_1_S_diff = shading_row_rotation - projected_solar_zenith
    thetas_2_S_diff = shaded_row_rotation - projected_solar_zenith
    thetaS_rotation_diff = projected_solar_zenith - cross_axis_slope
    cos_theta_2_S_diff_abs = np.abs(cosd(thetas_2_S_diff))
    # Eq. (12) of [1]
    t_asterisk = (
        0.5
        + np.abs(cosd(thetas_1_S_diff)) / cos_theta_2_S_diff_abs / 2
        + (
            np.sign(projected_solar_zenith)
            * surface_to_axis_offset
            / collector_width
            / cos_theta_2_S_diff_abs
            * (sind(thetas_2_S_diff) - sind(thetas_1_S_diff))
        )
        - (
            pitch
            / collector_width
            * cosd(thetaS_rotation_diff)
            / cos_theta_2_S_diff_abs
            / cosd(cross_axis_slope)
        )
    )
    return np.clip(t_asterisk, 0, 1) 
[docs]
def direct_martinez(
    poa_global,
    poa_direct,
    shaded_fraction,
    shaded_blocks,
    total_blocks,
):
    r"""
    A shading loss power factor for non-monolithic silicon
    modules and arrays with an arbitrary number of bypass diodes.
    This experimental model reduces the direct and circumsolar
    irradiance reaching the module's cells based on the number of *blocks*
    affected by the shadow.
    More on blocks in the *Notes* section and in [1]_.
    .. versionadded:: 0.11.0
    Parameters
    ----------
    poa_global : numeric
        Plane of array global irradiance. [W/m²].
    poa_direct : numeric
        Plane of array direct and circumsolar irradiance. [W/m²].
    shaded_fraction : numeric
        Fraction of module surface area that is shaded. [Unitless].
    shaded_blocks : numeric
        Number of blocks affected by the shadow. [Unitless].
        If a floating point number is provided, it will be rounded up.
    total_blocks : int
        Number of total blocks. Unitless.
    Returns
    -------
    shading_losses : numeric
        Fraction of DC power lost due to shading. [Unitless]
    Notes
    -----
    The implemented equations are (6) and (8) from [1]_:
    .. math::
        (1 - F_{ES}) = (1 - F_{GS}) \left(1 - \frac{N_{SB}}{N_{TB} + 1}\right)
        \quad \text{(6)}
        \left(1 - \frac{P_{S}}{P_{NS}}\right) = \left(1 -
        \frac{\left[(B + D^{CIR})(1 - F_{ES}) + D^{ISO} + R\right]}{G}\right)
        \quad \text{(8)}
    In (6), :math:`(1 - F_{ES})` is the correction factor to be multiplied by
    the direct and circumsolar irradiance, :math:`F_{GS}` is the shaded
    fraction of the collector, :math:`N_{SB}` is the number of shaded blocks
    and :math:`N_{TB}` is the number of total blocks.
    In (8), :math:`\frac{P_{S}}{P_{NS}}` is the fraction of DC power lost due
    to shading, :math:`P_{S}` is the power output of the shaded module,
    :math:`P_{NS}` is the power output of the non-shaded module,
    :math:`B + D^{CIR}` is the beam and circumsolar irradiance,
    :math:`D^{ISO} + R` is the sum of diffuse and albedo irradiances and
    :math:`G` is the global irradiance.
    **Blocks terminology:**
    A *block* is defined in [1]_ as a group of solar cells protected by a
    bypass diode. Also, a *block* is shaded when at least one of its
    cells is partially shaded.
    The total number of blocks and their layout depend on the module(s) used.
    Many manufacturers don't specify this information explicitly.
    However, these values can be inferred from:
    - the number of bypass diodes
    - where and how many junction boxes are present on the back of the module
    - whether or not the module is comprised of *half-cut cells*
    The latter two are heavily correlated.
    For example:
    1. A module with 1 bypass diode behaves as 1 block.
    2. A module with 3 bypass diodes and 1 junction box is likely to have 3
       blocks.
    3. A half-cut module with 3 junction boxes (split junction boxes) is
       likely to have 3x2 blocks. The number of blocks along the longest
       side of the module is 2 and along the shortest side is 3.
    4. A module without bypass diodes doesn't constitute a block, but may be
       part of one.
    Examples
    --------
    Minimal example. For a complete example, see
    :ref:`sphx_glr_gallery_shading_plot_martinez_shade_loss.py`.
    >>> import numpy as np
    >>> from pvlib import shading
    >>> total_blocks = 3  # blocks along the vertical of the module
    >>> POA_direct_and_circumsolar, POA_diffuse = 600, 80  # W/m²
    >>> POA_global = POA_direct_and_circumsolar + POA_diffuse
    >>> P_out_unshaded = 3000  # W
    >>> # calculation of the shaded fraction for the collector
    >>> shaded_fraction = shading.shaded_fraction1d(
    >>>     solar_zenith=80, solar_azimuth=180,
    >>>     axis_azimuth=90, shaded_row_rotation=25,
    >>>     collector_width=0.5, pitch=1, surface_to_axis_offset=0,
    >>>     cross_axis_slope=5.711, shading_row_rotation=50)
    >>> # calculation of the number of shaded blocks
    >>> shaded_blocks = np.ceil(total_blocks*shaded_fraction)
    >>> # apply the Martinez power losses to the calculated shading
    >>> loss_fraction = shading.direct_martinez(
    >>>     POA_global, POA_direct_and_circumsolar,
    >>>     shaded_fraction, shaded_blocks, total_blocks)
    >>> P_out_corrected = P_out_unshaded * (1 - loss_fraction)
    See Also
    --------
    shaded_fraction1d : to calculate 1-dimensional shaded fraction
    References
    ----------
    .. [1] F. Martínez-Moreno, J. Muñoz, and E. Lorenzo, 'Experimental model
        to estimate shading losses on PV arrays', Solar Energy Materials and
        Solar Cells, vol. 94, no. 12, pp. 2298-2303, Dec. 2010,
        :doi:`10.1016/j.solmat.2010.07.029`.
    """  # Contributed by Echedey Luis, 2024
    beam_factor = (  # Eq. (6) of [1]
        (1 - shaded_fraction)
        * (1 - np.ceil(shaded_blocks) / (1 + total_blocks))
    )
    return (  # Eq. (8) of [1]
        1
        - (
            poa_direct * beam_factor
            + (poa_global - poa_direct)  # diffuse and albedo
        )
        / poa_global
    )