test_irradiance.py

import datetime
from collections import OrderedDict
import warnings

import numpy as np
from numpy import array, nan
import pandas as pd

import pytest
from numpy.testing import (
    assert_almost_equal, assert_allclose, assert_equal, assert_raises
)
from pvlib import irradiance, albedo

from .conftest import (
    assert_frame_equal,
    assert_series_equal,
    requires_ephem,
    requires_numba,
)

from pvlib._deprecation import pvlibDeprecationWarning

# fixtures create realistic test input data
# test input data generated at Location(32.2, -111, 'US/Arizona', 700)
# test input data is hard coded to avoid dependencies on other parts of pvlib


@pytest.fixture
def times():
    # must include night values
    return pd.date_range(start='20140624', freq='6h', periods=4,
                         tz='US/Arizona')


@pytest.fixture
def irrad_data(times):
    return pd.DataFrame(np.array(
        [[0.,    0.,    0.],
         [79.73860422,  316.1949056,   40.46149818],
         [1042.48031487,  939.95469881,  118.45831879],
         [257.20751138,  646.22886049,   62.03376265]]),
        columns=['ghi', 'dni', 'dhi'], index=times)


@pytest.fixture
def ephem_data(times):
    return pd.DataFrame(np.array(
        [[124.0390863, 124.0390863, -34.0390863, -34.0390863,
          352.69550699,  -2.36677158],
         [82.85457044,  82.97705621,   7.14542956,   7.02294379,
          66.71410338,  -2.42072165],
         [10.56413562,  10.56725766,  79.43586438,  79.43274234,
          144.76567754,  -2.47457321],
         [72.41687122,  72.46903556,  17.58312878,  17.53096444,
          287.04104128,  -2.52831909]]),
        columns=['apparent_zenith', 'zenith', 'apparent_elevation',
                 'elevation', 'azimuth', 'equation_of_time'],
        index=times)


@pytest.fixture
def dni_et():
    return np.array(
        [1321.1655834833093, 1321.1655834833093, 1321.1655834833093,
         1321.1655834833093])


@pytest.fixture
def relative_airmass(times):
    return pd.Series([np.nan, 7.58831596, 1.01688136, 3.27930443], times)


# setup for et rad test. put it here for readability
timestamp = pd.Timestamp('20161026')
dt_index = pd.DatetimeIndex([timestamp])
doy = timestamp.dayofyear
dt_date = timestamp.date()
dt_datetime = datetime.datetime.combine(dt_date, datetime.time(0))
dt_np64 = np.datetime64(dt_datetime)
value = 1383.636203


@pytest.mark.parametrize('testval, expected', [
    (doy, value),
    (np.float64(doy), value),
    (dt_date, value),
    (dt_datetime, value),
    (dt_np64, value),
    (np.array([doy]), np.array([value])),
    (pd.Series([doy]), np.array([value])),
    (dt_index, pd.Series([value], index=dt_index)),
    (timestamp, value)
])
@pytest.mark.parametrize('method', [
    'asce', 'spencer', 'nrel', pytest.param('pyephem', marks=requires_ephem)])
def test_get_extra_radiation(testval, expected, method):
    out = irradiance.get_extra_radiation(testval, method=method)
    assert_allclose(out, expected, atol=10)


def test_get_extra_radiation_epoch_year():
    out = irradiance.get_extra_radiation(doy, method='nrel', epoch_year=2012)
    assert_allclose(out, 1382.4926804890767, atol=0.1)


@requires_numba
def test_get_extra_radiation_nrel_numba(times):
    with warnings.catch_warnings():
        # don't warn on method reload
        warnings.simplefilter("ignore")
        result = irradiance.get_extra_radiation(
            times, method='nrel', how='numba', numthreads=4)
        # and reset to no-numba state
        irradiance.get_extra_radiation(times, method='nrel')
    assert_allclose(result,
                    [1322.332316, 1322.296282, 1322.261205, 1322.227091])


def test_get_extra_radiation_invalid():
    with pytest.raises(ValueError):
        irradiance.get_extra_radiation(300, method='invalid')


def test_get_ground_diffuse_simple_float():
    result = irradiance.get_ground_diffuse(40, 900)
    assert_allclose(result, 26.32000014911496)


def test_get_ground_diffuse_simple_series(irrad_data):
    ground_irrad = irradiance.get_ground_diffuse(40, irrad_data['ghi'])
    assert ground_irrad.name == 'diffuse_ground'


def test_get_ground_diffuse_albedo_0(irrad_data):
    ground_irrad = irradiance.get_ground_diffuse(
        40, irrad_data['ghi'], albedo=0)
    assert (0 == ground_irrad).all()


def test_get_ground_diffuse_albedo_series(times):
    albedo = pd.Series(0.2, index=times)
    ground_irrad = irradiance.get_ground_diffuse(
        45, pd.Series(1000, index=times), albedo)
    expected = albedo * 0.5 * (1 - np.sqrt(2) / 2.) * 1000
    expected.name = 'diffuse_ground'
    assert_series_equal(ground_irrad, expected)


def test_grounddiffuse_albedo_invalid_surface(irrad_data):
    with pytest.raises(KeyError):
        irradiance.get_ground_diffuse(
            40, irrad_data['ghi'], surface_type='invalid')


def test_get_ground_diffuse_albedo_surface(irrad_data):
    result = irradiance.get_ground_diffuse(40, irrad_data['ghi'],
                                           surface_type='sand')
    assert_allclose(result, [0, 3.731058, 48.778813, 12.035025], atol=1e-4)


def test_isotropic_float():
    result = irradiance.isotropic(40, 100)
    assert_allclose(result, 88.30222215594891)


def test_isotropic_series(irrad_data):
    result = irradiance.isotropic(40, irrad_data['dhi'])
    assert_allclose(result, [0, 35.728402, 104.601328, 54.777191], atol=1e-4)


def test_klucher_series_float():
    # klucher inputs
    surface_tilt, surface_azimuth = 40.0, 180.0
    dhi, ghi = 100.0, 900.0
    solar_zenith, solar_azimuth = 20.0, 180.0
    # expect same result for floats and pd.Series
    expected = irradiance.klucher(
        surface_tilt, surface_azimuth,
        pd.Series(dhi), pd.Series(ghi),
        pd.Series(solar_zenith), pd.Series(solar_azimuth)
    )  # 94.99429931664851
    result = irradiance.klucher(
        surface_tilt, surface_azimuth, dhi, ghi, solar_zenith, solar_azimuth
    )
    assert_allclose(result, expected[0])


def test_klucher_series(irrad_data, ephem_data):
    result = irradiance.klucher(40, 180, irrad_data['dhi'], irrad_data['ghi'],
                                ephem_data['apparent_zenith'],
                                ephem_data['azimuth'])
    # pvlib matlab 1.4 does not contain the max(cos_tt, 0) correction
    # so, these values are different
    assert_allclose(result, [0., 36.789794, 109.209347, 56.965916], atol=1e-4)
    # expect same result for np.array and pd.Series
    expected = irradiance.klucher(
        40, 180, irrad_data['dhi'].values, irrad_data['ghi'].values,
        ephem_data['apparent_zenith'].values, ephem_data['azimuth'].values
    )
    assert_allclose(result, expected, atol=1e-4)


def test_haydavies(irrad_data, ephem_data, dni_et):
    result = irradiance.haydavies(
        40, 180, irrad_data['dhi'], irrad_data['dni'], dni_et,
        ephem_data['apparent_zenith'], ephem_data['azimuth'])
    # values from matlab 1.4 code
    assert_allclose(result, [0, 27.1775, 102.9949, 33.1909], atol=1e-4)


def test_haydavies_components(irrad_data, ephem_data, dni_et):
    keys = ['poa_sky_diffuse', 'poa_isotropic', 'poa_circumsolar',
            'poa_horizon']
    expected = pd.DataFrame(np.array(
        [[0, 27.1775, 102.9949, 33.1909],
         [0, 27.1775, 30.1818, 27.9837],
         [0, 0, 72.8130, 5.2071],
         [0, 0, 0, 0]]).T,
        columns=keys,
        index=irrad_data.index
    )
    # pandas
    result = irradiance.haydavies(
        40, 180, irrad_data['dhi'], irrad_data['dni'], dni_et,
        ephem_data['apparent_zenith'], ephem_data['azimuth'],
        return_components=True)
    assert_frame_equal(result, expected, check_less_precise=4)
    # numpy
    result = irradiance.haydavies(
        40, 180, irrad_data['dhi'].values, irrad_data['dni'].values, dni_et,
        ephem_data['apparent_zenith'].values, ephem_data['azimuth'].values,
        return_components=True)
    for key in keys:
        assert_allclose(result[key], expected[key], atol=1e-4)
    assert isinstance(result, dict)
    # scalar
    result = irradiance.haydavies(
        40, 180, irrad_data['dhi'].values[-1], irrad_data['dni'].values[-1],
        dni_et[-1], ephem_data['apparent_zenith'].values[-1],
        ephem_data['azimuth'].values[-1], return_components=True)
    for key in keys:
        assert_allclose(result[key], expected[key].iloc[-1], atol=1e-4)
    assert isinstance(result, dict)


def test_reindl(irrad_data, ephem_data, dni_et):
    result = irradiance.reindl(
        40, 180, irrad_data['dhi'], irrad_data['dni'], irrad_data['ghi'],
        dni_et, ephem_data['apparent_zenith'], ephem_data['azimuth'])
    # values from matlab 1.4 code
    assert_allclose(result, [0., 27.9412, 104.1317, 34.1663], atol=1e-4)


def test_king(irrad_data, ephem_data):
    result = irradiance.king(40, irrad_data['dhi'], irrad_data['ghi'],
                             ephem_data['apparent_zenith'])
    assert_allclose(result, [0, 44.629352, 115.182626, 79.719855], atol=1e-4)


def test_perez(irrad_data, ephem_data, dni_et, relative_airmass):
    dni = irrad_data['dni'].copy()
    dni.iloc[2] = np.nan
    out = irradiance.perez(40, 180, irrad_data['dhi'], dni,
                           dni_et, ephem_data['apparent_zenith'],
                           ephem_data['azimuth'], relative_airmass)
    expected = pd.Series(np.array(
        [0.,   31.46046871,  np.nan,   45.45539877]),
        index=irrad_data.index)
    assert_series_equal(out, expected, check_less_precise=2)


def test_perez_driesse(irrad_data, ephem_data, dni_et, relative_airmass):
    dni = irrad_data['dni'].copy()
    dni.iloc[2] = np.nan
    out = irradiance.perez_driesse(40, 180, irrad_data['dhi'], dni,
                                   dni_et, ephem_data['apparent_zenith'],
                                   ephem_data['azimuth'], relative_airmass)
    expected = pd.Series(np.array(
        [0.,   29.991,  np.nan,   47.397]),
        index=irrad_data.index)
    assert_series_equal(out, expected, check_less_precise=2)


def test_perez_driesse_airmass(irrad_data, ephem_data, dni_et):
    dni = irrad_data['dni'].copy()
    dni.iloc[2] = np.nan
    out = irradiance.perez_driesse(40, 180, irrad_data['dhi'], dni,
                                   dni_et, ephem_data['apparent_zenith'],
                                   ephem_data['azimuth'], airmass=None)
    expected = pd.Series(np.array(
        [0.,   29.991,  np.nan,   47.397]),
        index=irrad_data.index)
    assert_series_equal(out, expected, check_less_precise=2)


def test_perez_components(irrad_data, ephem_data, dni_et, relative_airmass):
    dni = irrad_data['dni'].copy()
    dni.iloc[2] = np.nan
    out = irradiance.perez(40, 180, irrad_data['dhi'], dni,
                           dni_et, ephem_data['apparent_zenith'],
                           ephem_data['azimuth'], relative_airmass,
                           return_components=True)
    expected = pd.DataFrame(np.array(
        [[0.,   31.46046871,  np.nan,   45.45539877],
         [0.,  26.84138589,          np.nan,  31.72696071],
         [0.,  0.,         np.nan,  4.47966439],
         [0.,  4.62212181,         np.nan,  9.25316454]]).T,
        columns=['poa_sky_diffuse', 'poa_isotropic', 'poa_circumsolar',
                 'poa_horizon'],
        index=irrad_data.index
    )
    expected_for_sum = expected['poa_sky_diffuse'].copy()
    expected_for_sum.iloc[2] = 0
    sum_components = out.iloc[:, 1:].sum(axis=1)
    sum_components.name = 'poa_sky_diffuse'

    assert_frame_equal(out, expected, check_less_precise=2)
    assert_series_equal(sum_components, expected_for_sum, check_less_precise=2)


def test_perez_driesse_components(irrad_data, ephem_data, dni_et,
                                  relative_airmass):
    dni = irrad_data['dni'].copy()
    dni.iloc[2] = np.nan
    out = irradiance.perez_driesse(40, 180, irrad_data['dhi'], dni,
                                   dni_et, ephem_data['apparent_zenith'],
                                   ephem_data['azimuth'], relative_airmass,
                                   return_components=True)

    expected = pd.DataFrame(np.array(
        [[0., 29.991, np.nan, 47.397],
         [0., 25.806, np.nan, 33.181],
         [0.,  0.000, np.nan,  4.197],
         [0.,  4.184, np.nan, 10.018]]).T,
        columns=['poa_sky_diffuse', 'poa_isotropic', 'poa_circumsolar',
                 'poa_horizon'],
        index=irrad_data.index
    )
    expected_for_sum = expected['poa_sky_diffuse'].copy()
    expected_for_sum.iloc[2] = 0
    sum_components = out.iloc[:, 1:].sum(axis=1)
    sum_components.name = 'poa_sky_diffuse'

    assert_frame_equal(out, expected, check_less_precise=2)
    assert_series_equal(sum_components, expected_for_sum, check_less_precise=2)


def test_perez_negative_horizon():
    times = pd.date_range(start='20190101 11:30:00', freq='1h',
                          periods=5, tz='US/Central')

    # Avoid test dependencies on functionality not being tested by hard-coding
    # the inputs. This data corresponds to Goodwin Creek in the afternoon on
    # 1/1/2019.
    # dni_e is slightly rounded from irradiance.get_extra_radiation
    # airmass from atmosphere.get_relative_airmas
    inputs = pd.DataFrame(np.array(
        [[158,         19,          1,          0,          0],
         [249,        165,        136,         93,         50],
         [57.746951,  57.564205,  60.813841,  66.989435,  75.353368],
         [171.003315, 187.346924, 202.974357, 216.725599, 228.317233],
         [1414,       1414,       1414,       1414,       1414],
         [1.869315,   1.859981,   2.044429,   2.544943,   3.900136]]).T,
        columns=['dni', 'dhi', 'solar_zenith',
                 'solar_azimuth', 'dni_extra', 'airmass'],
        index=times
    )

    out = irradiance.perez(34, 180, inputs['dhi'], inputs['dni'],
                           inputs['dni_extra'], inputs['solar_zenith'],
                           inputs['solar_azimuth'], inputs['airmass'],
                           model='allsitescomposite1990',
                           return_components=True)

    # sky_diffuse can be less than isotropic under certain conditions as
    # horizon goes negative
    expected = pd.DataFrame(np.array(
        [[281.410185, 152.20879, 123.867898, 82.836412, 43.517015],
         [166.785419, 142.24475, 119.173875, 83.525150, 45.725931],
         [113.548755,  16.09757,   9.956174,  3.142467,  0],
         [1.076010,  -6.13353,  -5.262151, -3.831230, -2.208923]]).T,
        columns=['poa_sky_diffuse', 'poa_isotropic', 'poa_circumsolar',
                 'poa_horizon'],
        index=times
    )

    expected_for_sum = expected['poa_sky_diffuse'].copy()
    sum_components = out.iloc[:, 1:].sum(axis=1)
    sum_components.name = 'poa_sky_diffuse'

    assert_frame_equal(out, expected, check_less_precise=2)
    assert_series_equal(sum_components, expected_for_sum, check_less_precise=2)


def test_perez_zero_dhi_and_dni_scalar():
    # Divides zero by zero.
    args = (20, 180, 0.0, 0.0, 1366.1, 89.96, 256.28, 37.32)

    out = irradiance.perez(*args)
    poa_sky_diffuse_expected = 0.0
    assert_equal(out, poa_sky_diffuse_expected)

    out = irradiance.perez(*args, return_components=True)
    expected = {
        "poa_sky_diffuse": poa_sky_diffuse_expected,
        "poa_isotropic": 0.0,
        "poa_circumsolar": 0.0,
        "poa_horizon": 0.0,
    }
    assert len(out) == len(expected)
    for key in expected.keys():
        assert_equal(out[key], expected[key], err_msg=key)

    assert_equal(
        out["poa_sky_diffuse"],
        out["poa_isotropic"] + out["poa_circumsolar"] + out["poa_horizon"],
    )


def test_perez_array_dhi_and_dni_combos():
    # Divides zero and non-zero by zero and various NaN division combos.
    args = (
        20, 180,
        np.array([0.0, 10.0, np.nan, np.nan, 0.0, 100.0, np.nan]),
        np.array([0.0, 0.0, 0.0, 100.0, np.nan, np.nan, np.nan]),
        1366.1, 89.96, 256.28, 37.32
    )

    out = irradiance.perez(*args)
    poa_sky_diffuse_expected = np.array(
        [0.0, 9.424924186619206, np.nan, np.nan, np.nan, np.nan, np.nan]
    )
    assert_allclose(out, poa_sky_diffuse_expected)

    out = irradiance.perez(*args, return_components=True)
    expected = {
        "poa_sky_diffuse": poa_sky_diffuse_expected,
        "poa_isotropic": np.array(
            [0.0, 9.162258932459126, np.nan, np.nan, np.nan, np.nan, np.nan]
        ),
        "poa_circumsolar": np.array(
            [0.0, 0.5187450944545264, np.nan, np.nan, np.nan, np.nan, np.nan]
        ),
        "poa_horizon": np.array(
            [0.0, -0.2560798402944465, np.nan, np.nan, np.nan, np.nan, np.nan]
        ),
    }
    assert len(out) == len(expected)
    for key in expected.keys():
        assert_allclose(out[key], expected[key], err_msg=key)

    assert_almost_equal(
        out["poa_sky_diffuse"],
        out["poa_isotropic"] + out["poa_circumsolar"] + out["poa_horizon"],
    )


def test_perez_array_dhi_and_dni_combos_nan_airmass():
    # Divides zero and non-zero by zero and various NaN division combos, when
    # airmass is NaN (e.g., for sun below horizon).
    args = (
        20, 180,
        np.array([0.0, 10.0, np.nan, np.nan, 0.0, 100.0, np.nan]),
        np.array([0.0, 0.0, 0.0, 100.0, np.nan, np.nan, np.nan]),
        1366.1, 91, 256.28, np.nan
    )

    out = irradiance.perez(*args)
    poa_sky_diffuse_expected = np.array(
        [0.0, 0.0, np.nan, np.nan, np.nan, np.nan, np.nan]
    )
    assert_allclose(out, poa_sky_diffuse_expected)

    out = irradiance.perez(*args, return_components=True)
    expected = {
        "poa_sky_diffuse": poa_sky_diffuse_expected,
        "poa_isotropic": np.array(
            [0.0, 0.0, np.nan, np.nan, np.nan, np.nan, np.nan]
        ),
        "poa_circumsolar": np.array(
            [0.0, 0.0, np.nan, np.nan, np.nan, np.nan, np.nan]
        ),
        "poa_horizon": np.array(
            [0.0, 0.0, np.nan, np.nan, np.nan, np.nan, np.nan]
        ),
    }
    assert len(out) == len(expected)
    for key in expected.keys():
        assert_allclose(out[key], expected[key], err_msg=key)

    assert_almost_equal(
        out["poa_sky_diffuse"],
        out["poa_isotropic"] + out["poa_circumsolar"] + out["poa_horizon"],
    )


def test_perez_zero_dhi_nonzero_dni_scalar():
    # Divides nonzero by zero.
    args = (20, 180, 0.0, 100.0, 1366.1, 89.96, 256.28, 37.32)

    with assert_raises(FloatingPointError):
        irradiance.perez(*args)

    with assert_raises(FloatingPointError):
        irradiance.perez(*args, return_components=True)


def test_perez_zero_dhi_nonzero_dni_array():
    # Divides nonzero by zero.
    args = (
        20, 180, np.array([0.0, 10.0, np.nan]), 100.0, 1366.1, 89.96,
        256.28, 37.32
    )

    with assert_raises(FloatingPointError):
        irradiance.perez(*args)

    with assert_raises(FloatingPointError):
        irradiance.perez(*args, return_components=True)


def test_perez_arrays(irrad_data, ephem_data, dni_et, relative_airmass):
    dni = irrad_data['dni'].copy()
    dni.iloc[2] = np.nan
    out = irradiance.perez(40, 180, irrad_data['dhi'].values, dni.values,
                           dni_et, ephem_data['apparent_zenith'].values,
                           ephem_data['azimuth'].values,
                           relative_airmass.values)
    expected = np.array(
        [0.,   31.46046871,  np.nan,   45.45539877])
    assert_allclose(out, expected, atol=1e-2)
    assert isinstance(out, np.ndarray)


def test_perez_driesse_arrays(irrad_data, ephem_data, dni_et,
                              relative_airmass):
    dni = irrad_data['dni'].copy()
    dni.iloc[2] = np.nan
    out = irradiance.perez_driesse(40, 180, irrad_data['dhi'].values,
                                   dni.values, dni_et,
                                   ephem_data['apparent_zenith'].values,
                                   ephem_data['azimuth'].values,
                                   relative_airmass.values)
    expected = np.array(
        [0.,   29.990,  np.nan,   47.396])
    assert_allclose(out, expected, atol=1e-2)
    assert isinstance(out, np.ndarray)


def test_perez_scalar():
    # copied values from fixtures
    out = irradiance.perez(40, 180, 118.45831879, 939.95469881,
                           1321.1655834833093, 10.56413562, 144.76567754,
                           1.01688136)
    # this will fail. out is ndarry with ndim == 0. fix in future version.
    # assert np.isscalar(out)
    assert_allclose(out, 109.084332)


def test_perez_driesse_scalar():
    # copied values from fixtures
    out = irradiance.perez_driesse(40, 180, 118.458, 939.954,
                                   1321.165, 10.564, 144.765, 1.016)
    # this will fail. out is ndarry with ndim == 0. fix in future version.
    # assert np.isscalar(out)
    assert_allclose(out, 110.341, atol=1e-2)


@pytest.mark.parametrize('model', ['isotropic', 'klucher', 'haydavies',
                                   'reindl', 'king', 'perez', 'perez-driesse'])
def test_sky_diffuse_zenith_close_to_90(model):
    # GH 432
    sky_diffuse = irradiance.get_sky_diffuse(
        30, 180, 89.999, 230,
        dni=10, ghi=51, dhi=50, dni_extra=1360, airmass=12, model=model)
    assert sky_diffuse < 100


def test_get_sky_diffuse_model_invalid():
    with pytest.raises(ValueError):
        irradiance.get_sky_diffuse(
            30, 180, 0, 180, 1000, 1100, 100, dni_extra=1360, airmass=1,
            model='invalid')


def test_get_sky_diffuse_missing_dni_extra():
    msg = 'dni_extra is required'
    with pytest.raises(ValueError, match=msg):
        irradiance.get_sky_diffuse(
            30, 180, 0, 180, 1000, 1100, 100, airmass=1,
            model='haydavies')


def test_get_sky_diffuse_missing_airmass(irrad_data, ephem_data, dni_et):
    # test assumes location is Tucson, AZ
    # calculated airmass should be the equivalent to fixture airmass
    dni = irrad_data['dni'].copy()
    dni.iloc[2] = np.nan
    out = irradiance.get_sky_diffuse(
        40, 180, ephem_data['apparent_zenith'], ephem_data['azimuth'], dni,
        irrad_data['ghi'], irrad_data['dhi'], dni_et,  model='perez')
    expected = pd.Series(np.array(
        [0., 31.46046871, np.nan, 45.45539877]),
        index=irrad_data.index)
    assert_series_equal(out, expected, check_less_precise=2)


def test_campbell_norman():
    expected = pd.DataFrame(np.array(
        [[863.859736967, 653.123094076, 220.65905025]]),
        columns=['ghi', 'dni', 'dhi'],
        index=[0])
    out = irradiance.campbell_norman(
        pd.Series([10]), pd.Series([0.5]), pd.Series([109764.21013135818]),
        dni_extra=1400)
    assert_frame_equal(out, expected)


def test_get_total_irradiance(irrad_data, ephem_data, dni_et,
                              relative_airmass):
    models = ['isotropic', 'klucher',
              'haydavies', 'reindl', 'king', 'perez', 'perez-driesse']

    for model in models:
        total = irradiance.get_total_irradiance(
            32, 180,
            ephem_data['apparent_zenith'], ephem_data['azimuth'],
            dni=irrad_data['dni'], ghi=irrad_data['ghi'],
            dhi=irrad_data['dhi'],
            dni_extra=dni_et, airmass=relative_airmass,
            model=model,
            surface_type='urban')

        assert total.columns.tolist() == ['poa_global', 'poa_direct',
                                          'poa_diffuse', 'poa_sky_diffuse',
                                          'poa_ground_diffuse']


@pytest.mark.parametrize('model', ['isotropic', 'klucher',
                                   'haydavies', 'reindl', 'king',
                                   'perez', 'perez-driesse'])
def test_get_total_irradiance_albedo(
        irrad_data, ephem_data, dni_et, relative_airmass, model):
    albedo = pd.Series(0.2, index=ephem_data.index)
    total = irradiance.get_total_irradiance(
        32, 180,
        ephem_data['apparent_zenith'], ephem_data['azimuth'],
        dni=irrad_data['dni'], ghi=irrad_data['ghi'],
        dhi=irrad_data['dhi'],
        dni_extra=dni_et, airmass=relative_airmass,
        model=model,
        albedo=albedo)

    assert total.columns.tolist() == ['poa_global', 'poa_direct',
                                      'poa_diffuse', 'poa_sky_diffuse',
                                      'poa_ground_diffuse']


@pytest.mark.parametrize('model', ['isotropic', 'klucher',
                                   'haydavies', 'reindl', 'king',
                                   'perez', 'perez-driesse'])
def test_get_total_irradiance_scalars(model):
    total = irradiance.get_total_irradiance(
        32, 180,
        10, 180,
        dni=1000, ghi=1100,
        dhi=100,
        dni_extra=1400, airmass=1,
        model=model,
        surface_type='urban')

    assert list(total.keys()) == ['poa_global', 'poa_direct',
                                  'poa_diffuse', 'poa_sky_diffuse',
                                  'poa_ground_diffuse']
    # test that none of the values are nan
    assert np.isnan(np.array(list(total.values()))).sum() == 0


def test_get_total_irradiance_missing_dni_extra():
    msg = 'dni_extra is required'
    with pytest.raises(ValueError, match=msg):
        irradiance.get_total_irradiance(
            32, 180,
            10, 180,
            dni=1000, ghi=1100,
            dhi=100,
            model='haydavies')


def test_get_total_irradiance_missing_airmass():
    total = irradiance.get_total_irradiance(
        32, 180,
        10, 180,
        dni=1000, ghi=1100,
        dhi=100,
        dni_extra=1400,
        model='perez')
    assert list(total.keys()) == ['poa_global', 'poa_direct',
                                  'poa_diffuse', 'poa_sky_diffuse',
                                  'poa_ground_diffuse']


def test_poa_components(irrad_data, ephem_data, dni_et, relative_airmass):
    aoi = irradiance.aoi(40, 180, ephem_data['apparent_zenith'],
                         ephem_data['azimuth'])
    gr_sand = irradiance.get_ground_diffuse(40, irrad_data['ghi'],
                                            surface_type='sand')
    diff_perez = irradiance.perez(
        40, 180, irrad_data['dhi'], irrad_data['dni'], dni_et,
        ephem_data['apparent_zenith'], ephem_data['azimuth'], relative_airmass)
    out = irradiance.poa_components(
        aoi, irrad_data['dni'], diff_perez, gr_sand)
    expected = pd.DataFrame(np.array(
        [[0.,  -0.,   0.,   0.,
            0.],
         [35.19456561,   0.,  35.19456561,  31.4635077,
            3.73105791],
         [956.18253696, 798.31939281, 157.86314414, 109.08433162,
          48.77881252],
         [90.99624896,  33.50143401,  57.49481495,  45.45978964,
          12.03502531]]),
        columns=['poa_global', 'poa_direct', 'poa_diffuse', 'poa_sky_diffuse',
                 'poa_ground_diffuse'],
        index=irrad_data.index)
    assert_frame_equal(out, expected)


@pytest.mark.parametrize('pressure,expected', [
    (93193,  [[830.46567,   0.79742,   0.93505],
              [676.18340,   0.63782,   3.02102]]),
    (None,   [[868.72425,   0.79742,   1.01664],
              [680.73800,   0.63782,   3.28463]]),
    (101325, [[868.72425,   0.79742,   1.01664],
              [680.73800,   0.63782,   3.28463]])
])
def test_disc_value(pressure, expected):
    # see GH 449 for pressure=None vs. 101325.
    columns = ['dni', 'kt', 'airmass']
    times = pd.DatetimeIndex(['2014-06-24T1200', '2014-06-24T1800'],
                             tz='America/Phoenix')
    ghi = pd.Series([1038.62, 254.53], index=times)
    zenith = pd.Series([10.567, 72.469], index=times)
    out = irradiance.disc(ghi, zenith, times, pressure=pressure)
    expected_values = np.array(expected)
    expected = pd.DataFrame(expected_values, columns=columns, index=times)
    # check the pandas dataframe. check_less_precise is weird
    assert_frame_equal(out, expected, check_less_precise=True)
    # use np.assert_allclose to check values more clearly
    assert_allclose(out.values, expected_values, atol=1e-5)


def test_disc_overirradiance():
    columns = ['dni', 'kt', 'airmass']
    ghi = np.array([3000])
    solar_zenith = np.full_like(ghi, 0)
    times = pd.date_range(start='2016-07-19 12:00:00', freq='1s',
                          periods=len(ghi), tz='America/Phoenix')
    out = irradiance.disc(ghi=ghi, solar_zenith=solar_zenith,
                          datetime_or_doy=times)
    expected = pd.DataFrame(np.array(
        [[8.72544336e+02, 1.00000000e+00, 9.99493933e-01]]),
        columns=columns, index=times)
    assert_frame_equal(out, expected)


def test_disc_min_cos_zenith_max_zenith():
    # map out behavior under difficult conditions with various
    # limiting kwargs settings
    columns = ['dni', 'kt', 'airmass']
    times = pd.DatetimeIndex(['2016-07-19 06:11:00'], tz='America/Phoenix')
    out = irradiance.disc(ghi=1.0, solar_zenith=89.99, datetime_or_doy=times)
    expected = pd.DataFrame(np.array(
        [[0.00000000e+00, 1.16046346e-02, 12.0]]),
        columns=columns, index=times)
    assert_frame_equal(out, expected)

    # max_zenith and/or max_airmass keep these results reasonable
    out = irradiance.disc(ghi=1.0, solar_zenith=89.99, datetime_or_doy=times,
                          min_cos_zenith=0)
    expected = pd.DataFrame(np.array(
        [[0.00000000e+00, 1.0, 12.0]]),
        columns=columns, index=times)
    assert_frame_equal(out, expected)

    # still get reasonable values because of max_airmass=12 limit
    out = irradiance.disc(ghi=1.0, solar_zenith=89.99, datetime_or_doy=times,
                          max_zenith=100)
    expected = pd.DataFrame(np.array(
        [[0., 1.16046346e-02, 12.0]]),
        columns=columns, index=times)
    assert_frame_equal(out, expected)

    # still get reasonable values because of max_airmass=12 limit
    out = irradiance.disc(ghi=1.0, solar_zenith=89.99, datetime_or_doy=times,
                          min_cos_zenith=0, max_zenith=100)
    expected = pd.DataFrame(np.array(
        [[277.50185968, 1.0, 12.0]]),
        columns=columns, index=times)
    assert_frame_equal(out, expected)

    # max_zenith keeps this result reasonable
    out = irradiance.disc(ghi=1.0, solar_zenith=89.99, datetime_or_doy=times,
                          min_cos_zenith=0, max_airmass=100)
    expected = pd.DataFrame(np.array(
        [[0.00000000e+00, 1.0, 36.39544757]]),
        columns=columns, index=times)
    assert_frame_equal(out, expected)

    # allow zenith to be close to 90 and airmass to be infinite
    # and we get crazy values
    out = irradiance.disc(ghi=1.0, solar_zenith=89.99, datetime_or_doy=times,
                          max_zenith=100, max_airmass=100)
    expected = pd.DataFrame(np.array(
        [[6.68577449e+03, 1.16046346e-02, 3.63954476e+01]]),
        columns=columns, index=times)
    assert_frame_equal(out, expected)

    # allow min cos zenith to be 0, zenith to be close to 90,
    # and airmass to be very big and we get even higher DNI values
    out = irradiance.disc(ghi=1.0, solar_zenith=89.99, datetime_or_doy=times,
                          min_cos_zenith=0, max_zenith=100, max_airmass=100)
    expected = pd.DataFrame(np.array(
        [[7.21238390e+03, 1., 3.63954476e+01]]),
        columns=columns, index=times)
    assert_frame_equal(out, expected)


def test_dirint_value():
    times = pd.DatetimeIndex(['2014-06-24T12-0700', '2014-06-24T18-0700'])
    ghi = pd.Series([1038.62, 254.53], index=times)
    zenith = pd.Series([10.567, 72.469], index=times)
    pressure = 93193.
    dirint_data = irradiance.dirint(ghi, zenith, times, pressure=pressure)
    assert_almost_equal(dirint_data.values,
                        np.array([868.8,  699.7]), 1)


def test_dirint_nans():
    times = pd.date_range(start='2014-06-24T12-0700', periods=5, freq='6h')
    ghi = pd.Series([np.nan, 1038.62, 1038.62, 1038.62, 1038.62], index=times)
    zenith = pd.Series([10.567, np.nan, 10.567, 10.567, 10.567], index=times)
    pressure = pd.Series([93193., 93193., np.nan, 93193., 93193.], index=times)
    temp_dew = pd.Series([10, 10, 10, np.nan, 10], index=times)
    dirint_data = irradiance.dirint(ghi, zenith, times, pressure=pressure,
                                    temp_dew=temp_dew)
    assert_almost_equal(dirint_data.values,
                        np.array([np.nan, np.nan, np.nan, np.nan, 893.1]), 1)


def test_dirint_tdew():
    times = pd.DatetimeIndex(['2014-06-24T12-0700', '2014-06-24T18-0700'])
    ghi = pd.Series([1038.62, 254.53], index=times)
    zenith = pd.Series([10.567, 72.469], index=times)
    pressure = 93193.
    dirint_data = irradiance.dirint(ghi, zenith, times, pressure=pressure,
                                    temp_dew=10)
    assert_almost_equal(dirint_data.values,
                        np.array([882.1,  672.6]), 1)


def test_dirint_no_delta_kt():
    times = pd.DatetimeIndex(['2014-06-24T12-0700', '2014-06-24T18-0700'])
    ghi = pd.Series([1038.62, 254.53], index=times)
    zenith = pd.Series([10.567, 72.469], index=times)
    pressure = 93193.
    dirint_data = irradiance.dirint(ghi, zenith, times, pressure=pressure,
                                    use_delta_kt_prime=False)
    assert_almost_equal(dirint_data.values,
                        np.array([861.9,  670.4]), 1)


def test_dirint_coeffs():
    coeffs = irradiance._get_dirint_coeffs()
    assert coeffs[0, 0, 0, 0] == 0.385230
    assert coeffs[0, 1, 2, 1] == 0.229970
    assert coeffs[3, 2, 6, 3] == 1.032260


def test_dirint_min_cos_zenith_max_zenith():
    # map out behavior under difficult conditions with various
    # limiting kwargs settings
    # times don't have any physical relevance
    times = pd.DatetimeIndex(['2014-06-24T12-0700', '2014-06-24T18-0700'])
    ghi = pd.Series([0, 1], index=times)
    solar_zenith = pd.Series([90, 89.99], index=times)

    out = irradiance.dirint(ghi, solar_zenith, times)
    expected = pd.Series([0.0, 0.0], index=times, name='dni')
    assert_series_equal(out, expected)

    out = irradiance.dirint(ghi, solar_zenith, times, min_cos_zenith=0)
    expected = pd.Series([0.0, 0.0], index=times, name='dni')
    assert_series_equal(out, expected)

    out = irradiance.dirint(ghi, solar_zenith, times, max_zenith=90)
    expected = pd.Series([0.0, 0.0], index=times, name='dni')
    assert_series_equal(out, expected, check_less_precise=True)

    out = irradiance.dirint(ghi, solar_zenith, times, min_cos_zenith=0,
                            max_zenith=90)
    expected = pd.Series([0.0, 144.264507], index=times, name='dni')
    assert_series_equal(out, expected, check_less_precise=True)

    out = irradiance.dirint(ghi, solar_zenith, times, min_cos_zenith=0,
                            max_zenith=100)
    expected = pd.Series([0.0, 144.264507], index=times, name='dni')
    assert_series_equal(out, expected, check_less_precise=True)


def test_ghi_from_poa_driesse(mocker):
    # inputs copied from test_gti_dirint
    times = pd.DatetimeIndex(
        ['2014-06-24T06-0700', '2014-06-24T09-0700', '2014-06-24T12-0700'])
    poa_global = np.array([20, 300, 1000])
    zenith = np.array([80, 45, 20])
    azimuth = np.array([90, 135, 180])
    surface_tilt = 30
    surface_azimuth = 180

    # test core function
    output = irradiance.ghi_from_poa_driesse_2024(
        surface_tilt, surface_azimuth, zenith, azimuth,
        poa_global, dni_extra=1366.1)

    expected = [22.089, 304.088, 931.143]
    assert_allclose(expected, output, atol=0.001)

    # test series output
    poa_global = pd.Series([20, 300, 1000], index=times)

    output = irradiance.ghi_from_poa_driesse_2024(
        surface_tilt, surface_azimuth, zenith, azimuth,
        poa_global, dni_extra=1366.1)

    assert isinstance(output, pd.Series)

    # test full_output option and special cases
    poa_global = np.array([0, 1500, np.nan])

    ghi, conv, niter = irradiance.ghi_from_poa_driesse_2024(
        surface_tilt, surface_azimuth, zenith, azimuth,
        poa_global, dni_extra=1366.1, full_output=True)

    expected = [0, np.nan, np.nan]
    assert_allclose(expected, ghi, atol=0.001)

    expected = [True, False, False]
    assert_allclose(expected, conv)

    expected = [0, -1, 0]
    assert_allclose(expected, niter)

    # test xtol argument
    poa_global = pd.Series([20, 300, 1000], index=times)
    # test exception
    xtol = -3.14159  # negative value raises exception in scipy.optimize.bisect
    with pytest.raises(ValueError, match=rf"xtol too small \({xtol} <= 0\)"):
        output = irradiance.ghi_from_poa_driesse_2024(
            surface_tilt, surface_azimuth, zenith, azimuth,
            poa_global, dni_extra=1366.1, xtol=xtol)
    # test propagation
    xtol = 3.141592
    bisect_spy = mocker.spy(irradiance, "bisect")
    output = irradiance.ghi_from_poa_driesse_2024(
        surface_tilt, surface_azimuth, zenith, azimuth,
        poa_global, dni_extra=1366.1, xtol=xtol)
    assert bisect_spy.call_args[1]["xtol"] == xtol


def test_ghi_from_poa_driesse_2023_deprecated():
    with pytest.warns(pvlibDeprecationWarning,
                      match="ghi_from_poa_driesse_2024"):
        irradiance.ghi_from_poa_driesse_2023(
            30, 180, 20, 180, 500, dni_extra=1366.1)


def test_gti_dirint():
    times = pd.DatetimeIndex(
        ['2014-06-24T06-0700', '2014-06-24T09-0700', '2014-06-24T12-0700'])
    poa_global = np.array([20, 300, 1000])
    aoi = np.array([100, 70, 10])
    zenith = np.array([80, 45, 20])
    azimuth = np.array([90, 135, 180])
    surface_tilt = 30
    surface_azimuth = 180

    # test defaults
    output = irradiance.gti_dirint(
        poa_global, aoi, zenith, azimuth, times, surface_tilt, surface_azimuth)

    expected_col_order = ['ghi', 'dni', 'dhi']
    expected = pd.DataFrame(array(
        [[21.05796198,    0.,   21.05796198],
         [291.40037163,   63.41290679,  246.56067523],
         [931.04078010,  695.94965324,  277.06172442]]),
        columns=expected_col_order, index=times)

    assert_frame_equal(output, expected)

    # test ignore calculate_gt_90
    output = irradiance.gti_dirint(
        poa_global, aoi, zenith, azimuth, times, surface_tilt, surface_azimuth,
        calculate_gt_90=False)

    expected_no_90 = expected.copy()
    expected_no_90.iloc[0, :] = np.nan

    assert_frame_equal(output, expected_no_90)

    # test pressure input
    pressure = 93193.
    output = irradiance.gti_dirint(
        poa_global, aoi, zenith, azimuth, times, surface_tilt, surface_azimuth,
        pressure=pressure)

    expected = pd.DataFrame(array(
        [[21.05796198,    0.,   21.05796198],
         [293.21310935,   63.27500913,  248.47092131],
         [932.46756378,  648.05001357,  323.49974813]]),
        columns=expected_col_order, index=times)

    assert_frame_equal(output, expected)

    # test albedo input
    albedo = 0.05
    output = irradiance.gti_dirint(
        poa_global, aoi, zenith, azimuth, times, surface_tilt, surface_azimuth,
        albedo=albedo)

    expected = pd.DataFrame(array(
        [[21.3592591,    0.,   21.3592591],
         [294.4985420,   66.25848451,  247.64671830],
         [941.7943404,  727.50552952,  258.16276278]]),
        columns=expected_col_order, index=times)

    assert_frame_equal(output, expected)

    # test with albedo as a Series
    albedo = pd.Series(0.05, index=times)
    output = irradiance.gti_dirint(
        poa_global, aoi, zenith, azimuth, times, surface_tilt, surface_azimuth,
        albedo=albedo)

    assert_frame_equal(output, expected)

    # test temp_dew input
    temp_dew = np.array([70, 80, 20])
    output = irradiance.gti_dirint(
        poa_global, aoi, zenith, azimuth, times, surface_tilt, surface_azimuth,
        temp_dew=temp_dew)

    expected = pd.DataFrame(array(
        [[21.05796198,    0.,           21.05796198],
         [295.06070190,   38.20346345,  268.0467738],
         [931.79627208,  689.81549269,  283.5817439]]),
        columns=expected_col_order, index=times)

    assert_frame_equal(output, expected)


def test_gti_dirint_data_error():
    times = pd.DatetimeIndex(
        ['2014-06-24T06-0700', '2014-06-24T09-0700', '2014-06-24T12-0700'])
    poa_global = np.array([20, 300, 1000])
    zenith = np.array([80, 45, 20])
    azimuth = np.array([90, 135, 180])
    surface_tilt = 30
    surface_azimuth = 180

    aoi = np.array([100, 170, 110])
    with pytest.raises(
            ValueError, match="There are no times with AOI < 90"):
        irradiance.gti_dirint(
            poa_global, aoi, zenith, azimuth, times, surface_tilt,
            surface_azimuth)


def test_erbs():
    index = pd.DatetimeIndex(['20190101']*3 + ['20190620'])
    ghi = pd.Series([0, 50, 1000, 1000], index=index)
    zenith = pd.Series([120, 85, 10, 10], index=index)
    expected = pd.DataFrame(np.array(
        [[0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
         [9.67192672e+01, 4.15703604e+01, 4.05723511e-01],
         [7.94205651e+02, 2.17860117e+02, 7.18132729e-01],
         [8.42001578e+02, 1.70790318e+02, 7.68214312e-01]]),
        columns=['dni', 'dhi', 'kt'], index=index)

    out = irradiance.erbs(ghi, zenith, index)

    assert_frame_equal(np.round(out, 0), np.round(expected, 0))


def test_erbs_driesse():
    index = pd.DatetimeIndex(['20190101']*3 + ['20190620'])
    ghi = pd.Series([0, 50, 1000, 1000], index=index)
    zenith = pd.Series([120, 85, 10, 10], index=index)
    # expected values are the same as for erbs original test
    expected = pd.DataFrame(np.array(
        [[0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
         [9.67192672e+01, 4.15703604e+01, 4.05723511e-01],
         [7.94205651e+02, 2.17860117e+02, 7.18132729e-01],
         [8.42001578e+02, 1.70790318e+02, 7.68214312e-01]]),
        columns=['dni', 'dhi', 'kt'], index=index)

    out = irradiance.erbs_driesse(ghi, zenith, index)

    assert_frame_equal(np.round(out, 0), np.round(expected, 0))

    # test with the new optional dni_extra argument
    dni_extra = irradiance.get_extra_radiation(index)

    out = irradiance.erbs_driesse(ghi, zenith, dni_extra=dni_extra)

    assert_frame_equal(np.round(out, 0), np.round(expected, 0))

    # test for required inputs
    with pytest.raises(ValueError):
        irradiance.erbs_driesse(ghi, zenith)


def test_boland():
    index = pd.DatetimeIndex(['20190101']*3 + ['20190620'])
    ghi = pd.Series([0, 50, 1000, 1000], index=index)
    zenith = pd.Series([120, 85, 10, 10], index=index)
    expected = pd.DataFrame(np.array(
        [[0.0,        0.0,        0.0],
         [81.9448546, 42.8580353, 0.405723511],
         [723.764990, 287.230626, 0.718132729],
         [805.020419, 207.209650, 0.768214312]]),
        columns=['dni', 'dhi', 'kt'], index=index)

    out = irradiance.boland(ghi, zenith, index)

    assert np.allclose(out, expected)


def test_orgill_hollands():
    index = pd.DatetimeIndex(['20190101']*3 + ['20190620'])
    ghi = pd.Series([0, 50, 1000, 1000], index=index)
    zenith = pd.Series([120, 85, 10, 10], index=index)
    expected = pd.DataFrame(np.array(
        [[0.0,        0.0,        0.0],
         [108.731366, 40.5234370, 0.405723511],
         [776.155771, 235.635779, 0.718132729],
         [835.696102, 177.000000, 0.768214312]]),
        columns=['dni', 'dhi', 'kt'], index=index)

    out = irradiance.orgill_hollands(ghi, zenith, index)

    assert np.allclose(out, expected)


def test_erbs_min_cos_zenith_max_zenith():
    # map out behavior under difficult conditions with various
    # limiting kwargs settings
    columns = ['dni', 'dhi', 'kt']
    times = pd.DatetimeIndex(['2016-07-19 06:11:00'], tz='America/Phoenix')

    # max_zenith keeps these results reasonable
    out = irradiance.erbs(ghi=1.0, zenith=89.99999,
                          datetime_or_doy=times, min_cos_zenith=0)
    expected = pd.DataFrame(np.array(
        [[0., 1., 1.]]),
        columns=columns, index=times)
    assert_frame_equal(out, expected)

    # 4-5 9s will produce bad behavior without max_zenith limit
    out = irradiance.erbs(ghi=1.0, zenith=89.99999,
                          datetime_or_doy=times, max_zenith=100)
    expected = pd.DataFrame(np.array(
        [[6.00115286e+03, 9.98952601e-01, 1.16377640e-02]]),
        columns=columns, index=times)
    assert_frame_equal(out, expected)

    # 1-2 9s will produce bad behavior without either limit
    out = irradiance.erbs(ghi=1.0, zenith=89.99, datetime_or_doy=times,
                          min_cos_zenith=0, max_zenith=100)
    expected = pd.DataFrame(np.array(
        [[4.78419761e+03, 1.65000000e-01, 1.00000000e+00]]),
        columns=columns, index=times)
    assert_frame_equal(out, expected)

    # check default behavior under hardest condition
    out = irradiance.erbs(ghi=1.0, zenith=90, datetime_or_doy=times)
    expected = pd.DataFrame(np.array(
        [[0., 1., 0.01163776]]),
        columns=columns, index=times)
    assert_frame_equal(out, expected)


def test_erbs_all_scalar():
    ghi = 1000
    zenith = 10
    doy = 180

    expected = OrderedDict()
    expected['dni'] = 8.42358014e+02
    expected['dhi'] = 1.70439297e+02
    expected['kt'] = 7.68919470e-01

    out = irradiance.erbs(ghi, zenith, doy)

    for k, v in out.items():
        assert_allclose(v, expected[k], 1e-2)


def test_dirindex(times):
    ghi = pd.Series([0, 0, 1038.62, 254.53], index=times)
    ghi_clearsky = pd.Series(
        np.array([0., 79.73860422, 1042.48031487, 257.20751138]),
        index=times
    )
    dni_clear = pd.Series(
        np.array([0., 316.1949056, 939.95469881, 646.22886049]),
        index=times
    )
    zenith = pd.Series(
        np.array([124.0390863, 82.85457044, 10.56413562, 72.41687122]),
        index=times
    )
    pressure = 93193.
    tdew = 10.
    out = irradiance.dirindex(ghi, ghi_clearsky, dni_clear,
                              zenith, times, pressure=pressure,
                              temp_dew=tdew)
    dirint_close_values = irradiance.dirint(ghi, zenith, times,
                                            pressure=pressure,
                                            use_delta_kt_prime=True,
                                            temp_dew=tdew).values
    expected_out = np.array([np.nan, 0., 748.31562800, 630.73752100])

    tolerance = 1e-8
    assert np.allclose(out, expected_out, rtol=tolerance, atol=0,
                       equal_nan=True)
    tol_dirint = 0.2
    assert np.allclose(
        out.values,
        dirint_close_values,
        rtol=tol_dirint,
        atol=0,
        equal_nan=True)


def test_dirindex_min_cos_zenith_max_zenith():
    # map out behavior under difficult conditions with various
    # limiting kwargs settings
    # times don't have any physical relevance
    times = pd.DatetimeIndex(['2014-06-24T12-0700', '2014-06-24T18-0700'])
    ghi = pd.Series([0, 1], index=times)
    ghi_clearsky = pd.Series([0, 1], index=times)
    dni_clear = pd.Series([0, 5], index=times)
    solar_zenith = pd.Series([90, 89.99], index=times)

    out = irradiance.dirindex(ghi, ghi_clearsky, dni_clear, solar_zenith,
                              times)
    expected = pd.Series([nan, nan], index=times)
    assert_series_equal(out, expected)

    out = irradiance.dirindex(ghi, ghi_clearsky, dni_clear, solar_zenith,
                              times, min_cos_zenith=0)
    expected = pd.Series([nan, nan], index=times)
    assert_series_equal(out, expected)

    out = irradiance.dirindex(ghi, ghi_clearsky, dni_clear, solar_zenith,
                              times, max_zenith=90)
    expected = pd.Series([nan, nan], index=times)
    assert_series_equal(out, expected)

    out = irradiance.dirindex(ghi, ghi_clearsky, dni_clear, solar_zenith,
                              times, min_cos_zenith=0, max_zenith=100)
    expected = pd.Series([nan, 5.], index=times)
    assert_series_equal(out, expected)


def test_dirint_array_inputs():
    """np.array and pd.Series inputs work correctly. GH #2751"""
    times = pd.date_range('2023-06-21 10:00', periods=2, freq='h', tz='UTC')

    # np.array input → should return pd.Series
    result = irradiance.dirint(
        ghi=np.array([500.0, 400.0]),
        solar_zenith=np.array([45.0, 50.0]),
        times=times,
        use_delta_kt_prime=False
    )
    assert isinstance(result, pd.Series)
    assert result.iloc[0] > 0

    # pd.Series input → should return pd.Series
    times2 = pd.date_range('2023-06-21 10:00', periods=3, freq='h', tz='UTC')
    result2 = irradiance.dirint(
        ghi=pd.Series([400, 500, 300], index=times2),
        solar_zenith=pd.Series([50, 40, 60], index=times2),
        times=times2
    )
    assert isinstance(result2, pd.Series)
    assert (result2 >= 0).all()


def test_dni():
    ghi = pd.Series([90, 100, 100, 100, 100])
    dhi = pd.Series([100, 90, 50, 50, 50])
    zenith = pd.Series([80, 100, 85, 70, 85])
    dni_clear = pd.Series([50, 50, 200, 50, 300])

    dni = irradiance.dni(ghi, dhi, zenith,
                         dni_clear=dni_clear, clearsky_tolerance=2)
    assert_series_equal(dni,
                        pd.Series([float('nan'), float('nan'), 400,
                                   146.190220008, 573.685662283]))

    dni = irradiance.dni(ghi, dhi, zenith)
    assert_series_equal(dni,
                        pd.Series([float('nan'), float('nan'), 573.685662283,
                                   146.190220008, 573.685662283]))


@pytest.mark.parametrize(
    'surface_tilt,surface_azimuth,solar_zenith,' +
    'solar_azimuth,aoi_expected,aoi_proj_expected',
    [(0, 0, 0, 0, 0, 1),
     (30, 180, 30, 180, 0, 1),
     (30, 180, 150, 0, 180, -1),
     (90, 0, 30, 60, 75.5224878, 0.25),
     (90, 0, 30, 170, 119.4987042, -0.4924038)])
def test_aoi_and_aoi_projection(surface_tilt, surface_azimuth, solar_zenith,
                                solar_azimuth, aoi_expected,
                                aoi_proj_expected):
    aoi = irradiance.aoi(surface_tilt, surface_azimuth, solar_zenith,
                         solar_azimuth)
    assert_allclose(aoi, aoi_expected, atol=1e-5)

    aoi_projection = irradiance.aoi_projection(
        surface_tilt, surface_azimuth, solar_zenith, solar_azimuth)
    assert_allclose(aoi_projection, aoi_proj_expected, atol=1e-6)


def test_aoi_projection_precision():
    # GH 1185 -- test that aoi_projection does not exceed 1.0, and when
    # given identical inputs, the returned projection is very close to 1.0

    # scalars
    zenith = 89.26778228223463
    azimuth = 60.932028605997004
    projection = irradiance.aoi_projection(zenith, azimuth, zenith, azimuth)
    assert projection <= 1
    assert np.isclose(projection, 1)

    # arrays
    zeniths = np.array([zenith])
    azimuths = np.array([azimuth])
    projections = irradiance.aoi_projection(zeniths, azimuths,
                                            zeniths, azimuths)
    assert all(projections <= 1)
    assert all(np.isclose(projections, 1))
    assert projections.dtype == np.dtype('float64')


@pytest.fixture
def airmass_kt():
    # disc algorithm stopped at am=12. test am > 12 for out of range behavior
    return np.array([1, 5, 12, 20])


def test_kt_kt_prime_factor(airmass_kt):
    out = irradiance._kt_kt_prime_factor(airmass_kt)
    expected = np.array([0.999971,  0.723088,  0.548811,  0.471068])
    assert_allclose(out, expected, atol=1e-5)


def test_clearsky_index():
    ghi = np.array([-1., 0., 1., 500., 1000., np.nan])
    ghi_measured, ghi_modeled = np.meshgrid(ghi, ghi)
    # default max_clearsky_index
    with np.errstate(invalid='ignore', divide='ignore'):
        out = irradiance.clearsky_index(ghi_measured, ghi_modeled)
    expected = np.array(
        [[1., 0., 0., 0., 0., np.nan],
         [0., 0., 0., 0., 0., np.nan],
         [0., 0., 1., 2., 2., np.nan],
         [0., 0., 0.002, 1., 2., np.nan],
         [0., 0., 0.001, 0.5, 1., np.nan],
         [np.nan, np.nan, np.nan, np.nan, np.nan, np.nan]])
    assert_allclose(out, expected, atol=0.001)
    # specify max_clearsky_index
    with np.errstate(invalid='ignore', divide='ignore'):
        out = irradiance.clearsky_index(ghi_measured, ghi_modeled,
                                        max_clearsky_index=1.5)
    expected = np.array(
        [[1., 0., 0., 0., 0., np.nan],
         [0., 0., 0., 0., 0., np.nan],
         [0., 0., 1., 1.5, 1.5, np.nan],
         [0., 0., 0.002, 1., 1.5, np.nan],
         [0., 0., 0.001, 0.5, 1., np.nan],
         [np.nan, np.nan, np.nan, np.nan, np.nan, np.nan]])
    assert_allclose(out, expected, atol=0.001)
    # scalars
    out = irradiance.clearsky_index(10, 1000)
    expected = 0.01
    assert_allclose(out, expected, atol=0.001)
    # series
    times = pd.date_range(start='20180601', periods=2, freq='12h')
    ghi_measured = pd.Series([100,  500], index=times)
    ghi_modeled = pd.Series([500, 1000], index=times)
    out = irradiance.clearsky_index(ghi_measured, ghi_modeled)
    expected = pd.Series([0.2, 0.5], index=times)
    assert_series_equal(out, expected)


def test_clearness_index():
    ghi = np.array([-1, 0, 1, 1000])
    solar_zenith = np.array([180, 90, 89.999, 0])
    ghi, solar_zenith = np.meshgrid(ghi, solar_zenith)
    # default min_cos_zenith
    out = irradiance.clearness_index(ghi, solar_zenith, 1370)
    # np.set_printoptions(precision=3, floatmode='maxprec', suppress=True)
    expected = np.array(
        [[0., 0., 0.011, 2.],
         [0., 0., 0.011, 2.],
         [0., 0., 0.011, 2.],
         [0., 0., 0.001, 0.73]])
    assert_allclose(out, expected, atol=0.001)
    # specify min_cos_zenith
    with np.errstate(invalid='ignore', divide='ignore'):
        out = irradiance.clearness_index(ghi, solar_zenith, 1400,
                                         min_cos_zenith=0)
    expected = np.array(
        [[0.,   nan, 2., 2.],
         [0., 0., 2., 2.],
         [0., 0., 2., 2.],
         [0., 0., 0.001, 0.714]])
    assert_allclose(out, expected, atol=0.001)
    # specify max_clearness_index
    out = irradiance.clearness_index(ghi, solar_zenith, 1370,
                                     max_clearness_index=0.82)
    expected = np.array(
        [[0.,  0.,  0.011,  0.82],
         [0.,  0.,  0.011,  0.82],
         [0.,  0.,  0.011,  0.82],
         [0.,  0.,  0.001,  0.73]])
    assert_allclose(out, expected, atol=0.001)
    # specify min_cos_zenith and max_clearness_index
    with np.errstate(invalid='ignore', divide='ignore'):
        out = irradiance.clearness_index(ghi, solar_zenith, 1400,
                                         min_cos_zenith=0,
                                         max_clearness_index=0.82)
    expected = np.array(
        [[0.,    nan,  0.82,  0.82],
         [0.,  0.,  0.82,  0.82],
         [0.,  0.,  0.82,  0.82],
         [0.,  0.,  0.001,  0.714]])
    assert_allclose(out, expected, atol=0.001)
    # scalars
    out = irradiance.clearness_index(1000, 10, 1400)
    expected = 0.725
    assert_allclose(out, expected, atol=0.001)
    # series
    times = pd.date_range(start='20180601', periods=2, freq='12h')
    ghi = pd.Series([0, 1000], index=times)
    solar_zenith = pd.Series([90, 0], index=times)
    extra_radiation = pd.Series([1360, 1400], index=times)
    out = irradiance.clearness_index(ghi, solar_zenith, extra_radiation)
    expected = pd.Series([0, 0.714285714286], index=times)
    assert_series_equal(out, expected)


def test_clearness_index_zenith_independent(airmass_kt):
    clearness_index = np.array([-1, 0, .1, 1])
    clearness_index, airmass_kt = np.meshgrid(clearness_index, airmass_kt)
    out = irradiance.clearness_index_zenith_independent(clearness_index,
                                                        airmass_kt)
    expected = np.array(
        [[0., 0., 0.1, 1.],
         [0., 0., 0.138, 1.383],
         [0., 0., 0.182, 1.822],
         [0., 0., 0.212, 2.]])
    assert_allclose(out, expected, atol=0.001)
    # test max_clearness_index
    out = irradiance.clearness_index_zenith_independent(
        clearness_index, airmass_kt, max_clearness_index=0.82)
    expected = np.array(
        [[0.,  0.,  0.1,  0.82],
         [0.,  0.,  0.138,  0.82],
         [0.,  0.,  0.182,  0.82],
         [0.,  0.,  0.212,  0.82]])
    assert_allclose(out, expected, atol=0.001)
    # scalars
    out = irradiance.clearness_index_zenith_independent(.4, 2)
    expected = 0.443
    assert_allclose(out, expected, atol=0.001)
    # series
    times = pd.date_range(start='20180601', periods=2, freq='12h')
    clearness_index = pd.Series([0, .5], index=times)
    airmass = pd.Series([np.nan, 2], index=times)
    out = irradiance.clearness_index_zenith_independent(clearness_index,
                                                        airmass)
    expected = pd.Series([np.nan, 0.553744437562], index=times)
    assert_series_equal(out, expected)


def test_complete_irradiance():
    # Generate dataframe to test on
    times = pd.date_range('2010-07-05 7:00:00-0700', periods=2, freq='h')
    i = pd.DataFrame({'ghi': [372.103976116, 497.087579068],
                      'dhi': [356.543700, 465.44400],
                      'dni': [49.63565561689957, 62.10624908037814]},
                     index=times)
    # Define the solar position and clearsky dataframe
    solar_position = pd.DataFrame({'apparent_zenith': [71.7303262449161,
                                                       59.369],
                                   'zenith': [71.7764, 59.395]},
                                  index=pd.DatetimeIndex([
                                      '2010-07-05 07:00:00-0700',
                                      '2010-07-05 08:00:00-0700']))
    clearsky = pd.DataFrame({'dni': [625.5254880160008, 778.7766443075865],
                             'ghi': [246.3508023804681, 469.461381740857],
                             'dhi': [50.25488725346631, 72.66909939636372]},
                            index=pd.DatetimeIndex([
                                '2010-07-05 07:00:00-0700',
                                '2010-07-05 08:00:00-0700']))
    # Test scenario where DNI is generated via component sum equation
    complete_df = irradiance.complete_irradiance(
        solar_position.apparent_zenith,
        ghi=i.ghi,
        dhi=i.dhi,
        dni=None,
        dni_clear=clearsky.dni)
    # Assert that the ghi, dhi, and dni series match the original dataframe
    # values
    assert_frame_equal(complete_df, i)
    # Test scenario where GHI is generated via component sum equation
    complete_df = irradiance.complete_irradiance(
        solar_position.apparent_zenith,
        ghi=None,
        dhi=i.dhi,
        dni=i.dni,
        dni_clear=clearsky.dni)
    # Assert that the ghi, dhi, and dni series match the original dataframe
    # values
    assert_frame_equal(complete_df, i)
    # Test scenario where DHI is generated via component sum equation
    complete_df = irradiance.complete_irradiance(
        solar_position.apparent_zenith,
        ghi=i.ghi,
        dhi=None,
        dni=i.dni,
        dni_clear=clearsky.dni)
    # Assert that the ghi, dhi, and dni series match the original dataframe
    # values
    assert_frame_equal(complete_df, i)
    # Test scenario where all parameters are passed (throw error)
    with pytest.raises(ValueError):
        irradiance.complete_irradiance(solar_position.apparent_zenith,
                                       ghi=i.ghi,
                                       dhi=i.dhi,
                                       dni=i.dni,
                                       dni_clear=clearsky.dni)
    # Test scenario where only one parameter is passed (throw error)
    with pytest.raises(ValueError):
        irradiance.complete_irradiance(solar_position.apparent_zenith,
                                       ghi=None,
                                       dhi=None,
                                       dni=i.dni,
                                       dni_clear=clearsky.dni)


def test_louche():

    index = pd.DatetimeIndex(['20190101']*3 + ['20190620']*1)
    ghi = pd.Series([0, 50, 1000, 1000], index=index)
    zenith = pd.Series([91, 85, 10, 10], index=index)
    expected = pd.DataFrame(np.array(
        [[0, 0, 0],
         [130.089669, 38.661938, 0.405724],
         [828.498650, 184.088106, 0.718133],
         [887.407348, 126.074364, 0.768214]]),
        columns=['dni', 'dhi', 'kt'], index=index)

    out = irradiance.louche(ghi, zenith, index)

    assert_frame_equal(out, expected)


def test_SURFACE_ALBEDOS_deprecated():
    with pytest.warns(pvlibDeprecationWarning, match='SURFACE_ALBEDOS has been'
                      ' moved to the albedo module as of v0.11.0. Please use'
                      ' pvlib.albedo.SURFACE_ALBEDOS.'):
        irradiance.SURFACE_ALBEDOS


@pytest.mark.filterwarnings("ignore:SURFACE_ALBEDOS")
def test_SURFACE_ALBEDO_equals():
    assert irradiance.SURFACE_ALBEDOS == albedo.SURFACE_ALBEDOS


def test_diffuse_par_spitters():
    solar_zenith, global_diffuse_fraction = np.meshgrid(
        [90, 85, 75, 60, 40, 30, 10, 0], [0.01, 0.1, 0.3, 0.6, 0.8, 0.99]
    )
    solar_zenith = solar_zenith.ravel()
    global_diffuse_fraction = global_diffuse_fraction.ravel()
    result = irradiance.diffuse_par_spitters(
        solar_zenith, global_diffuse_fraction
    )
    expected = np.array([
        0.01300, 0.01290, 0.01226, 0.01118, 0.01125, 0.01189, 0.01293, 0.01300,
        0.12970, 0.12874, 0.12239, 0.11174, 0.11236, 0.11868, 0.12905, 0.12970,
        0.38190, 0.37931, 0.36201, 0.33273, 0.33446, 0.35188, 0.38014, 0.38190,
        0.71520, 0.71178, 0.68859, 0.64787, 0.65033, 0.67472, 0.71288, 0.71520,
        0.88640, 0.88401, 0.86755, 0.83745, 0.83931, 0.85746, 0.88478, 0.88640,
        0.99591, 0.99576, 0.99472, 0.99270, 0.99283, 0.99406, 0.99581, 0.99591,
    ])  # fmt: skip
    assert_allclose(result, expected, atol=1e-5)