Skip to content

Commit d2b1ac5

Browse files
VolkrBVolkrB
committed
DeltaT with numba decorators, Unit tests
1 parent d3f46ff commit d2b1ac5

File tree

2 files changed

+158
-68
lines changed

2 files changed

+158
-68
lines changed

‎pvlib/spa.py

Lines changed: 137 additions & 68 deletions
Original file line numberDiff line numberDiff line change
@@ -25,25 +25,29 @@ def nocompile(*args, **kwargs):
2525

2626
if os.getenv('PVLIB_USE_NUMBA', '0') != '0':
2727
try:
28-
from numba import jit, __version__
28+
from numba import jit, __version__, vectorize
2929
except ImportError:
3030
warnings.warn('Could not import numba, falling back to numpy ' +
3131
'calculation')
32+
vcompile = np.vectorize
3233
jcompile = nocompile
3334
USE_NUMBA = False
3435
else:
3536
major, minor = __version__.split('.')[:2]
3637
if int(major + minor) >= 17:
3738
# need at least numba >= 0.17.0
39+
vcompile = vectorize
3840
jcompile = jit
3941
USE_NUMBA = True
4042
else:
4143
warnings.warn('Numba version must be >= 0.17.0, falling back to ' +
4244
'numpy')
45+
vcompile = np.vectorize
4346
jcompile = nocompile
4447
USE_NUMBA = False
4548
else:
46-
jcompile = nocompile
49+
vcompile = np.vectorize
50+
jcompile = nocompile
4751
USE_NUMBA = False
4852

4953

@@ -426,8 +430,7 @@ def nocompile(*args, **kwargs):
426430
])
427431

428432

429-
@jcompile('float64(int64, int64, int64, int64, int64, int64, int64)',
430-
nopython=True)
433+
@jcompile('float64(int64, int64, int64, int64, int64, int64, int64)', nopython=True)
431434
def julian_day_dt(year, month, day, hour, minute, second, microsecond):
432435
"""This is the original way to calculate the julian day from the NREL paper.
433436
However, it is much faster to convert to unix/epoch time and then convert
@@ -973,6 +976,130 @@ def solar_position_loop(unixtime, loc_args, out):
973976
out[5, i] = eot
974977

975978

979+
@jcompile(nopython=True)
980+
def delta_tcalc(year, month):
981+
982+
y = year + (month - 0.5)/12
983+
984+
if year < -500:
985+
delta_t = -20+32*((y-1820)/100)**2
986+
987+
elif -500 <= year < 500:
988+
delta_t = 10583.6-1014.41*(y/100)+33.78311*(y/100)**2-5.952053*(y/100)**3- \
989+
0.1798452*(y/100)**4+0.022174192*(y/100)**5+0.0090316521*(y/100)**6
990+
991+
elif 500 <= year < 1600:
992+
delta_t = 1574.2-556.01*((y-1000)/100)+71.23472*((y-1000)/100)**2+ \
993+
0.319781*((y-1000)/100)**3-0.8503463 *((y-1000)/100)**4- \
994+
0.005050998*((y-1000)/100)**5+0.0083572073*((y-1000)/100)**6
995+
996+
elif 1600 <= year < 1700:
997+
delta_t = 120-0.9808*(y-1600)-0.01532*(y-1600)**2+(y-1600)**3/7129
998+
999+
elif 1700 <= year < 1800:
1000+
delta_t = 8.83+0.1603*(y-1700)-0.0059285*(y-1700)**2+0.00013336*(y-1700)**3- \
1001+
(y-1700)**4/1174000
1002+
1003+
elif 1800 <= year < 1860:
1004+
delta_t = 13.72-0.332447*(y-1800)+0.0068612*(y-1800)**2+0.0041116*(y-1800)**3- \
1005+
0.00037436*(y-1800)**4+0.0000121272*(y-1800)**5-0.0000001699*(y-1800)**6+ \
1006+
0.000000000875*(y-1800)**7
1007+
1008+
elif 1860 <= year < 1900:
1009+
delta_t = 7.6+0.5737*(y-1860)-0.251754*(y-1860)**2+0.01680668*(y-1860)**3- \
1010+
0.0004473624*(y-1860)**4+(y-1860)**5/233174
1011+
1012+
elif 1900 <= year < 1920:
1013+
delta_t = -2.79+1.494119*(y-1900)-0.0598939*(y-1900)**2+0.0061966*(y-1900)**3- \
1014+
0.000197*(y-1900)**4
1015+
1016+
elif 1920 <= year < 1941:
1017+
delta_t = 21.20+0.84493*(y-1920)-0.076100*(y-1920)**2+0.0020936*(y-1920)**3
1018+
1019+
elif 1941 <= year < 1961:
1020+
delta_t = 29.07+0.407*(y-1950)-(y-1950)**2/233+(y-1950)**3/2547
1021+
1022+
elif 1961 <= year < 1986:
1023+
delta_t = 45.45+1.067*(y-1975)-(y-1975)**2/260-(y-1975)**3/718
1024+
1025+
elif 1986 <= year < 2005:
1026+
delta_t = 63.86+0.3345*(y-2000)-0.060374*(y-2000)**2+0.0017275*(y-2000)**3+ \
1027+
0.000651814*(y-2000)**4+0.00002373599*(y-2000)**5
1028+
1029+
elif 2005 <= year < 2050:
1030+
delta_t = 62.92+0.32217*(y-2000)+0.005589*(y-2000)**2
1031+
1032+
elif 2050 <= year < 2150:
1033+
delta_t = -20+32*((y-1820)/100)**2-0.5628*(2150-y)
1034+
1035+
elif year > 2150:
1036+
delta_t = -20+32*((y-1820)/100)**2
1037+
1038+
return delta_t
1039+
1040+
1041+
@vcompile
1042+
def delta_tcalc_array(year, month):
1043+
1044+
y = year + (month - 0.5)/12
1045+
1046+
if year < -500:
1047+
delta_t = -20+32*((y-1820)/100)**2
1048+
1049+
elif -500 <= year < 500:
1050+
delta_t = 10583.6-1014.41*(y/100)+33.78311*(y/100)**2-5.952053*(y/100)**3- \
1051+
0.1798452*(y/100)**4+0.022174192*(y/100)**5+0.0090316521*(y/100)**6
1052+
1053+
elif 500 <= year < 1600:
1054+
delta_t = 1574.2-556.01*((y-1000)/100)+71.23472*((y-1000)/100)**2+ \
1055+
0.319781*((y-1000)/100)**3-0.8503463 *((y-1000)/100)**4- \
1056+
0.005050998*((y-1000)/100)**5+0.0083572073*((y-1000)/100)**6
1057+
1058+
elif 1600 <= year < 1700:
1059+
delta_t = 120-0.9808*(y-1600)-0.01532*(y-1600)**2+(y-1600)**3/7129
1060+
1061+
elif 1700 <= year < 1800:
1062+
delta_t = 8.83+0.1603*(y-1700)-0.0059285*(y-1700)**2+0.00013336*(y-1700)**3- \
1063+
(y-1700)**4/1174000
1064+
1065+
elif 1800 <= year < 1860:
1066+
delta_t = 13.72-0.332447*(y-1800)+0.0068612*(y-1800)**2+0.0041116*(y-1800)**3- \
1067+
0.00037436*(y-1800)**4+0.0000121272*(y-1800)**5-0.0000001699*(y-1800)**6+ \
1068+
0.000000000875*(y-1800)**7
1069+
1070+
elif 1860 <= year < 1900:
1071+
delta_t = 7.6+0.5737*(y-1860)-0.251754*(y-1860)**2+0.01680668*(y-1860)**3- \
1072+
0.0004473624*(y-1860)**4+(y-1860)**5/233174
1073+
1074+
elif 1900 <= year < 1920:
1075+
delta_t = -2.79+1.494119*(y-1900)-0.0598939*(y-1900)**2+0.0061966*(y-1900)**3- \
1076+
0.000197*(y-1900)**4
1077+
1078+
elif 1920 <= year < 1941:
1079+
delta_t = 21.20+0.84493*(y-1920)-0.076100*(y-1920)**2+0.0020936*(y-1920)**3
1080+
1081+
elif 1941 <= year < 1961:
1082+
delta_t = 29.07+0.407*(y-1950)-(y-1950)**2/233+(y-1950)**3/2547
1083+
1084+
elif 1961 <= year < 1986:
1085+
delta_t = 45.45+1.067*(y-1975)-(y-1975)**2/260-(y-1975)**3/718
1086+
1087+
elif 1986 <= year < 2005:
1088+
delta_t = 63.86+0.3345*(y-2000)-0.060374*(y-2000)**2+0.0017275*(y-2000)**3+ \
1089+
0.000651814*(y-2000)**4+0.00002373599*(y-2000)**5
1090+
1091+
elif 2005 <= year < 2050:
1092+
delta_t = 62.92+0.32217*(y-2000)+0.005589*(y-2000)**2
1093+
1094+
elif 2050 <= year < 2150:
1095+
delta_t = -20+32*((y-1820)/100)**2-0.5628*(2150-y)
1096+
1097+
elif year > 2150:
1098+
delta_t = -20+32*((y-1820)/100)**2
1099+
1100+
return delta_t
1101+
1102+
9761103
def solar_position_numba(unixtime, lat, lon, elev, pressure, temp, delta_t,
9771104
atmos_refract, numthreads, sst=False, esd=False):
9781105
"""Calculate the solar position using the numba compiled functions
@@ -1301,70 +1428,12 @@ def earthsun_distance(unixtime, delta_t, numthreads):
13011428

13021429
def calculate_deltaT(year, month):
13031430

1304-
if year > 3000 or year < -500: pvl_logger.warning('Delta T'
1305-
' is unknown for years before -1999 and after 3000.'
1306-
' DeltaT values will be calculated, but the calculations'
1307-
' are not intended to be used for these years.')
1308-
1309-
#Taken from http://eclipse.gsfc.nasa.gov/SEcat5/deltatpoly.html
1310-
1311-
y = year + (month - 0.5)/12
1312-
1313-
1314-
if year < -500:
1315-
delta_t = -20+32*((y-1820)/100)**2
1316-
1317-
elif -500 <= year < 500:
1318-
delta_t = 10583.6-1014.41*(y/100)+33.78311*(y/100)**2-5.952053*(y/100)**3- \
1319-
0.1798452*(y/100)**4+0.022174192*(y/100)**5+0.0090316521*(y/100)**6
1431+
if np.isscalar(year) is True:
1432+
deltaT=delta_tcalc(year, month)
1433+
return deltaT
13201434

1321-
elif 500 <= year < 1600:
1322-
delta_t = 1574.2-556.01*((y-1000)/100)+71.23472*((y-1000)/100)**2+ \
1323-
0.319781*((y-1000)/100)**3-0.8503463 *((y-1000)/100)**4- \
1324-
0.005050998*((y-1000)/100)**5+0.0083572073*((y-1000)/100)**6
1325-
1326-
elif 1600 <= year < 1700:
1327-
delta_t = 120-0.9808*(y-1600)-0.01532*(y-1600)**2+(y-1600)**3/7129
1328-
1329-
elif 1700 <= year < 1800:
1330-
delta_t = 8.83+0.1603*(y-1700)-0.0059285*(y-1700)**2+0.00013336*(y-1700)**3- \
1331-
(y-1700)**4/1174000
1332-
1333-
elif 1800 <= year < 1860:
1334-
delta_t = 13.72-0.332447*(y-1800)+0.0068612*(y-1800)**2+0.0041116*(y-1800)**3- \
1335-
0.00037436*(y-1800)**4+0.0000121272*(y-1800)**5-0.0000001699*(y-1800)**6+ \
1336-
0.000000000875*(y-1800)**7
1337-
1338-
elif 1860 <= year < 1900:
1339-
delta_t = 7.6+0.5737*(y-1860)-0.251754*(y-1860)**2+0.01680668*(y-1860)**3- \
1340-
0.0004473624*(y-1860)**4+(y-1860)**5/233174
1341-
1342-
elif 1900 <= year < 1920:
1343-
delta_t = -2.79+1.494119*(y-1900)-0.0598939*(y-1900)**2+0.0061966*(y-1900)**3- \
1344-
0.000197*(y-1900)**4
1345-
1346-
elif 1920 <= year < 1941:
1347-
delta_t = 21.20+0.84493*(y-1920)-0.076100*(y-1920)**2+0.0020936*(y-1920)**3
1348-
1349-
elif 1941 <= year < 1961:
1350-
delta_t = 29.07+0.407*(y-1950)-(y-1950)**2/233+(y-1950)**3/2547
1351-
1352-
elif 1961 <= year < 1986:
1353-
delta_t = 45.45+1.067*(y-1975)-(y-1975)**2/260-(y-1975)**3/718
1354-
1355-
elif 1986 <= year < 2005:
1356-
delta_t = 63.86+0.3345*(y-2000)-0.060374*(y-2000)**2+0.0017275*(y-2000)**3+ \
1357-
0.000651814*(y-2000)**4+0.00002373599*(y-2000)**5
1358-
1359-
elif 2005 <= year < 2050:
1360-
delta_t = 62.92+0.32217*(y-2000)+0.005589*(y-2000)**2
1361-
1362-
elif 2050 <= year < 2150:
1363-
delta_t = -20+32*((y-1820)/100)**2-0.5628*(2150-y)
1364-
1365-
elif year > 2150:
1366-
delta_t = -20+32*((y-1820)/100)**2
1367-
1368-
return delta_t
1435+
else:
1436+
deltaT=delta_tcalc_array(year,month)
1437+
return deltaT
13691438

13701439

‎pvlib/test/test_spa.py

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -76,6 +76,13 @@
7676
theta0 = 90 - e0
7777
Gamma = 14.340241
7878
Phi = 194.340241
79+
yearArray = np.array([-499, 500, 1000, 1500, 1800, 1900, 1950, 1970, 1985, 1990, 2000, 2005])
80+
monthArray = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12])
81+
DTactual = 54.413442486
82+
DTactualArray = np.array([1.7184831e+04, 5.7088051e+03, 1.5730419e+03,
83+
1.9801820e+02, 1.3596506e+01, -2.1171894e+00,
84+
2.9289261e+01, 4.0824887e+01, 5.4724581e+01,
85+
5.7426651e+01, 6.4108015e+01, 6.5038015e+01])
7986

8087

8188
class SpaBase(object):
@@ -316,6 +323,20 @@ def test_earthsun_distance(self):
316323
result = self.spa.earthsun_distance(unixtimes, 64.0, 1)
317324
assert_almost_equal(R, result, 6)
318325

326+
def test_delta_tcalc_array(self):
327+
result = self.spa.delta_tcalc_array(yearArray, monthArray)
328+
assert_almost_equal(DTactualArray, result, 3)
329+
330+
def test_delta_tcalc(self):
331+
result = self.spa.delta_tcalc(1985 ,2)
332+
assert_almost_equal(DTactual, result, 8)
333+
334+
def test_calculate_deltaT(self):
335+
resultArray = self.spa.calculate_deltaT(yearArray, monthArray)
336+
assert_almost_equal(DTactualArray, resultArray, 3)
337+
resultScalar = self.spa.calculate_deltaT(1985 ,2)
338+
assert_almost_equal(DTactual, resultScalar, 8)
339+
319340
class NumpySpaTest(unittest.TestCase, SpaBase):
320341
"""Import spa without compiling to numba then run tests"""
321342
@classmethod

0 commit comments

Comments
 (0)