pvlib
Failed unit test for nrel 2008 solar position algorithm
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.
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)
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.
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!
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: @.***>
On second thought, @kylefmacdonald perhaps you would like a entry level introduction to contributing to pvlib?