Failed unit test for nrel 2008 solar position algorithm

[kurt-rhee]

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:

1. run julian_day_dt() with the following parameters

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

2. 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]

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]

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]

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]

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 <#2077 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AH4Y3NUD6TVS3I7SIZ5IBUTZF6TURAVCNFSM6AAAAABI3TN5Z2VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDCNJRGE3DCMRRGI> . You are receiving this because you authored the thread.Message ID: ***@***.***>

[kurt-rhee]

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