Failed unit test for nrel 2008 solar position algorithm

kurt-rhee commented 3 weeks ago

Describe the bug In the documentation for the nrel 2008 solar position algorithm, there are a set of "unit tests" in table A4.1

https://www.nrel.gov/docs/fy08osti/34302.pdf

To Reproduce Steps to reproduce the behavior:

run julian_day_dt() with the following parameters

year = 837,
month = 4,
day = 10,
hour = 7,
minute = 12,
second = 0,
microsecond = 0

result is 2026867.8

Expected behavior expected result in table A4.1 is 2026871.8

Screenshots

Versions:

`pvlib.10.5:
python: 3.11

Additional context Perhaps this has something to do with the gregorian / julian date cutoff? All of the other tests complete as expected.

kandersolar commented 3 weeks ago

Perhaps this has something to do with the gregorian / julian date cutoff?

It seems so. I reproduce the expected values exactly using this modified version of the function that follows the text description around equation 4 regarding the calendar difference:

import pvlib

def julian_day_dt(year, month, day, hour, minute, second, microsecond):
    """This is the original way to calculate the julian day from the NREL paper.
    However, it is much faster to convert to unix/epoch time and then convert
    to julian day. Note that the date must be UTC."""
    if month <= 2:
        year = year-1
        month = month+12
    a = int(year/100)
    b = 2 - a + int(a * 0.25)
    frac_of_day = (microsecond / 1e6 + (second + minute * 60 + hour * 3600)
                   ) * 1.0 / (3600*24)
    d = day + frac_of_day
    jd = int(365.25 * (year + 4716)) + int(30.6001 * (month + 1)) + d - 1524.5
    if jd > 2299160.0:
        jd += b

    return jd

tests = [
    ((2000, 1, 1, 12, 0, 0), 2451545.0),
    ((1999, 1, 1, 0, 0, 0), 2451179.5),
    ((1987, 1, 27, 0, 0, 0), 2446822.5),
    ((1987, 6, 19, 12, 0, 0), 2446966.0),
    ((1988, 1, 27, 0, 0, 0), 2447187.5),
    ((1988, 6, 19, 12, 0, 0), 2447332.0),
    ((1900, 1, 1, 0, 0, 0), 2415020.5),
    ((1600, 1, 1, 0, 0, 0), 2305447.5),

    ((1600, 12, 31, 0, 0, 0), 2305812.5),
    ((837, 4, 10, 7, 12, 0), 2026871.8),
    ((-123, 12, 31, 0, 0, 0), 1676496.5),
    ((-122, 1, 1, 0, 0, 0), 1676497.5),
    ((-1000, 7, 12, 12, 0, 0), 1356001.0),
    ((-1000, 2, 29, 0, 0, 0), 1355866.5),
    ((-1001, 8, 17, 21, 36, 0), 1355671.4),
    ((-4712, 1, 1, 12, 0, 0), 0.0),
]

print("pvlib\tmodified")
for inputs, expected in tests:
    ret1 = pvlib.spa.julian_day_dt(*inputs, microsecond=0)
    ret2 = julian_day_dt(*inputs, microsecond=0)
    print(ret1 - expected, '\t', ret2 - expected)

kurt-rhee commented 3 weeks ago

Nice catch @kandersolar! Quickest issue resolution of all time.

Do you want me to put in a PR for this or is it small enough that an outside commit isn't worth it.

kandersolar commented 3 weeks ago

It's worth fixing IMHO. We know our current implementation doesn't follow its reference, and we even have test cases to prove it. This function isn't used in pvlib's main SPA calculation flow, so I doubt this issue affects many (if any) real-life applications, but in service of pvlib's goal of providing "reference" implementations I think it deserves to be fixed. A PR would be welcome!

kurt-rhee commented 3 weeks ago

Alright! I won't be able to this month but when I'm back I'll put a PR in to get this fixed.

On Wed, Jun 5, 2024, 5:22 PM Kevin Anderson @.***> wrote:

It's worth fixing IMHO. We know our current implementation doesn't follow its reference, and we even have test cases to prove it. This function isn't used in pvlib's main SPA calculation flow, so I doubt this issue affects many (if any) real-life applications, but in service of pvlib's goal of providing "reference" implementations I think it deserves to be fixed. A PR would be welcome!

— Reply to this email directly, view it on GitHub https://github.com/pvlib/pvlib-python/issues/2077#issuecomment-2151161212, or unsubscribe https://github.com/notifications/unsubscribe-auth/AH4Y3NUD6TVS3I7SIZ5IBUTZF6TURAVCNFSM6AAAAABI3TN5Z2VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDCNJRGE3DCMRRGI . You are receiving this because you authored the thread.Message ID: @.***>

kurt-rhee commented 3 weeks ago

On second thought, @kylefmacdonald perhaps you would like a entry level introduction to contributing to pvlib?

pvlib / pvlib-python

Failed unit test for nrel 2008 solar position algorithm #2077