From f5bd6349e1172782a42473c467a4d5422af34622 Mon Sep 17 00:00:00 2001 From: "Adam R. Jensen" <39184289+AdamRJensen@users.noreply.github.com> Date: Tue, 29 Sep 2026 10:29:40 +0200 Subject: [PATCH 1/7] Add pvlib.spa.transit_sunrise_sunset to documentation --- docs/sphinx/source/reference/solarposition.rst | 1 + docs/sphinx/source/whatsnew/v0.16.2.rst | 2 ++ pvlib/spa.py | 7 ++++++- 3 files changed, 9 insertions(+), 1 deletion(-) diff --git a/docs/sphinx/source/reference/solarposition.rst b/docs/sphinx/source/reference/solarposition.rst index dbe93d5392..42d02c5de1 100644 --- a/docs/sphinx/source/reference/solarposition.rst +++ b/docs/sphinx/source/reference/solarposition.rst @@ -41,6 +41,7 @@ Functions for calculating sunrise, sunset and transit times. solarposition.sun_rise_set_transit_ephem solarposition.sun_rise_set_transit_spa solarposition.sun_rise_set_transit_geometric + spa.transit_sunrise_sunset The spa module contains the implementation of the built-in NREL SPA diff --git a/docs/sphinx/source/whatsnew/v0.16.2.rst b/docs/sphinx/source/whatsnew/v0.16.2.rst index f704525a12..833da40157 100644 --- a/docs/sphinx/source/whatsnew/v0.16.2.rst +++ b/docs/sphinx/source/whatsnew/v0.16.2.rst @@ -22,6 +22,8 @@ Enhancements Documentation ~~~~~~~~~~~~~ +* Add :py:func:`pvlib.spa.transit_sunrise_sunset` to the API reference. + (:issue:`2255`) Testing diff --git a/pvlib/spa.py b/pvlib/spa.py index fe072b2cc9..18c7843f45 100644 --- a/pvlib/spa.py +++ b/pvlib/spa.py @@ -1108,7 +1108,7 @@ def solar_position(unixtime, lat, lon, elev, pressure, temp, delta_t, def transit_sunrise_sunset(dates, lat, lon, delta_t, numthreads): """ Calculate the sun transit, sunrise, and sunset - for a set of dates at a given location. + for a set of dates at a given location using the NREL SPA [1]_. Parameters ---------- @@ -1129,6 +1129,11 @@ def transit_sunrise_sunset(dates, lat, lon, delta_t, numthreads): ------- tuple : (transit, sunrise, sunset) localized to UTC + References + ---------- + .. [1] I. Reda and A. Andreas, Solar position algorithm for solar + radiation applications. Solar Energy, vol. 76, no. 5, pp. 577-589, + 2004. :doi:`10.1016/j.solener.2003.12.003` """ if ((dates % 86400) != 0.0).any(): From 48f95c56a89bb15cb038f9899b102c3ea4a42156 Mon Sep 17 00:00:00 2001 From: "Adam R. Jensen" <39184289+AdamRJensen@users.noreply.github.com> Date: Tue, 29 Sep 2026 10:31:27 +0200 Subject: [PATCH 2/7] Update v0.16.2.rst --- docs/sphinx/source/whatsnew/v0.16.2.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/sphinx/source/whatsnew/v0.16.2.rst b/docs/sphinx/source/whatsnew/v0.16.2.rst index 833da40157..6cc4056a0a 100644 --- a/docs/sphinx/source/whatsnew/v0.16.2.rst +++ b/docs/sphinx/source/whatsnew/v0.16.2.rst @@ -23,7 +23,7 @@ Enhancements Documentation ~~~~~~~~~~~~~ * Add :py:func:`pvlib.spa.transit_sunrise_sunset` to the API reference. - (:issue:`2255`) + (:issue:`2255`, :pull:`2873`) Testing From 969a7703e8522c5d347452f3c42ee1b7026e6acf Mon Sep 17 00:00:00 2001 From: Kevin Anderson Date: Mon, 26 Oct 2026 18:20:17 +0100 Subject: [PATCH 3/7] Apply batched suggestions from code review Co-authored-by: Rajiv Daxini <143435106+RDaxini@users.noreply.github.com> --- pvlib/spa.py | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/pvlib/spa.py b/pvlib/spa.py index 18c7843f45..deb6a07e38 100644 --- a/pvlib/spa.py +++ b/pvlib/spa.py @@ -1108,7 +1108,7 @@ def solar_position(unixtime, lat, lon, elev, pressure, temp, delta_t, def transit_sunrise_sunset(dates, lat, lon, delta_t, numthreads): """ Calculate the sun transit, sunrise, and sunset - for a set of dates at a given location using the NREL SPA [1]_. + for a set of dates at a given location using the NREL SPA [1]_, [2]_. Parameters ---------- @@ -1134,7 +1134,16 @@ def transit_sunrise_sunset(dates, lat, lon, delta_t, numthreads): .. [1] I. Reda and A. Andreas, Solar position algorithm for solar radiation applications. Solar Energy, vol. 76, no. 5, pp. 577-589, 2004. :doi:`10.1016/j.solener.2003.12.003` - """ + References + ---------- + .. [1] I. Reda and A. Andreas, Solar position algorithm for solar + radiation applications. Solar Energy, vol. 76, no. 5, pp. 577-589, + 2004. :doi:`10.1016/j.solener.2003.12.003` + .. [2] I. Reda, and A. Andreas. "Corrigendum to “Solar position + algorithm for solar radiation applications." [Solar Energy 76 + (2004) 577–589] Solar Energy, vol. 81, no. 6, pp. 838, 2007. + :doi:`10.1016/j.solener.2007.01.003` + """ if ((dates % 86400) != 0.0).any(): raise ValueError('Input dates must be at 00:00 UTC') From 2a420ba6280b4085a784f4d0ab7693c691a177cd Mon Sep 17 00:00:00 2001 From: Kevin Anderson Date: Mon, 26 Oct 2026 18:27:03 +0100 Subject: [PATCH 4/7] Apply batched suggestions from code review Co-authored-by: Kevin Anderson --- pvlib/spa.py | 2753 +++++++++++++++++++++++++------------------------- 1 file changed, 1375 insertions(+), 1378 deletions(-) diff --git a/pvlib/spa.py b/pvlib/spa.py index 6e36a084cd..6536b7747c 100644 --- a/pvlib/spa.py +++ b/pvlib/spa.py @@ -1,1140 +1,1137 @@ -""" -Calculate the sun's position using the NREL Solar Position Algorithm -(SPA) either using numpy arrays or compiling the code to machine language -with numba. -""" - -# Contributors: -# Created by Tony Lorenzo (@alorenzo175), Univ. of Arizona, 2015 - -import os -import threading -import warnings - -import numpy as np - - -# this block is a way to use an environment variable to switch between -# compiling the functions with numba or just use numpy -def nocompile(*args, **kwargs): - return lambda func: func - - -if os.getenv('PVLIB_USE_NUMBA', '0') != '0': - try: - from numba import jit - except ImportError: - warnings.warn('Could not import numba, falling back to numpy ' + - 'calculation') - jcompile = nocompile - USE_NUMBA = False - else: - jcompile = jit - USE_NUMBA = True -else: - jcompile = nocompile - USE_NUMBA = False - - -# heliocentric longitude coefficients -L0 = np.array([ - [175347046.0, 0.0, 0.0], - [3341656.0, 4.6692568, 6283.07585], - [34894.0, 4.6261, 12566.1517], - [3497.0, 2.7441, 5753.3849], - [3418.0, 2.8289, 3.5231], - [3136.0, 3.6277, 77713.7715], - [2676.0, 4.4181, 7860.4194], - [2343.0, 6.1352, 3930.2097], - [1324.0, 0.7425, 11506.7698], - [1273.0, 2.0371, 529.691], - [1199.0, 1.1096, 1577.3435], - [990.0, 5.233, 5884.927], - [902.0, 2.045, 26.298], - [857.0, 3.508, 398.149], - [780.0, 1.179, 5223.694], - [753.0, 2.533, 5507.553], - [505.0, 4.583, 18849.228], - [492.0, 4.205, 775.523], - [357.0, 2.92, 0.067], - [317.0, 5.849, 11790.629], - [284.0, 1.899, 796.298], - [271.0, 0.315, 10977.079], - [243.0, 0.345, 5486.778], - [206.0, 4.806, 2544.314], - [205.0, 1.869, 5573.143], - [202.0, 2.458, 6069.777], - [156.0, 0.833, 213.299], - [132.0, 3.411, 2942.463], - [126.0, 1.083, 20.775], - [115.0, 0.645, 0.98], - [103.0, 0.636, 4694.003], - [102.0, 0.976, 15720.839], - [102.0, 4.267, 7.114], - [99.0, 6.21, 2146.17], - [98.0, 0.68, 155.42], - [86.0, 5.98, 161000.69], - [85.0, 1.3, 6275.96], - [85.0, 3.67, 71430.7], - [80.0, 1.81, 17260.15], - [79.0, 3.04, 12036.46], - [75.0, 1.76, 5088.63], - [74.0, 3.5, 3154.69], - [74.0, 4.68, 801.82], - [70.0, 0.83, 9437.76], - [62.0, 3.98, 8827.39], - [61.0, 1.82, 7084.9], - [57.0, 2.78, 6286.6], - [56.0, 4.39, 14143.5], - [56.0, 3.47, 6279.55], - [52.0, 0.19, 12139.55], - [52.0, 1.33, 1748.02], - [51.0, 0.28, 5856.48], - [49.0, 0.49, 1194.45], - [41.0, 5.37, 8429.24], - [41.0, 2.4, 19651.05], - [39.0, 6.17, 10447.39], - [37.0, 6.04, 10213.29], - [37.0, 2.57, 1059.38], - [36.0, 1.71, 2352.87], - [36.0, 1.78, 6812.77], - [33.0, 0.59, 17789.85], - [30.0, 0.44, 83996.85], - [30.0, 2.74, 1349.87], - [25.0, 3.16, 4690.48] -]) -L1 = np.array([ - [628331966747.0, 0.0, 0.0], - [206059.0, 2.678235, 6283.07585], - [4303.0, 2.6351, 12566.1517], - [425.0, 1.59, 3.523], - [119.0, 5.796, 26.298], - [109.0, 2.966, 1577.344], - [93.0, 2.59, 18849.23], - [72.0, 1.14, 529.69], - [68.0, 1.87, 398.15], - [67.0, 4.41, 5507.55], - [59.0, 2.89, 5223.69], - [56.0, 2.17, 155.42], - [45.0, 0.4, 796.3], - [36.0, 0.47, 775.52], - [29.0, 2.65, 7.11], - [21.0, 5.34, 0.98], - [19.0, 1.85, 5486.78], - [19.0, 4.97, 213.3], - [17.0, 2.99, 6275.96], - [16.0, 0.03, 2544.31], - [16.0, 1.43, 2146.17], - [15.0, 1.21, 10977.08], - [12.0, 2.83, 1748.02], - [12.0, 3.26, 5088.63], - [12.0, 5.27, 1194.45], - [12.0, 2.08, 4694.0], - [11.0, 0.77, 553.57], - [10.0, 1.3, 6286.6], - [10.0, 4.24, 1349.87], - [9.0, 2.7, 242.73], - [9.0, 5.64, 951.72], - [8.0, 5.3, 2352.87], - [6.0, 2.65, 9437.76], - [6.0, 4.67, 4690.48] -]) -L2 = np.array([ - [52919.0, 0.0, 0.0], - [8720.0, 1.0721, 6283.0758], - [309.0, 0.867, 12566.152], - [27.0, 0.05, 3.52], - [16.0, 5.19, 26.3], - [16.0, 3.68, 155.42], - [10.0, 0.76, 18849.23], - [9.0, 2.06, 77713.77], - [7.0, 0.83, 775.52], - [5.0, 4.66, 1577.34], - [4.0, 1.03, 7.11], - [4.0, 3.44, 5573.14], - [3.0, 5.14, 796.3], - [3.0, 6.05, 5507.55], - [3.0, 1.19, 242.73], - [3.0, 6.12, 529.69], - [3.0, 0.31, 398.15], - [3.0, 2.28, 553.57], - [2.0, 4.38, 5223.69], - [2.0, 3.75, 0.98] -]) -L3 = np.array([ - [289.0, 5.844, 6283.076], - [35.0, 0.0, 0.0], - [17.0, 5.49, 12566.15], - [3.0, 5.2, 155.42], - [1.0, 4.72, 3.52], - [1.0, 5.3, 18849.23], - [1.0, 5.97, 242.73] -]) -L4 = np.array([ - [114.0, 3.142, 0.0], - [8.0, 4.13, 6283.08], - [1.0, 3.84, 12566.15] -]) -L5 = np.array([ - [1.0, 3.14, 0.0] -]) - - -# heliocentric latitude coefficients -B0 = np.array([ - [280.0, 3.199, 84334.662], - [102.0, 5.422, 5507.553], - [80.0, 3.88, 5223.69], - [44.0, 3.7, 2352.87], - [32.0, 4.0, 1577.34] -]) -B1 = np.array([ - [9.0, 3.9, 5507.55], - [6.0, 1.73, 5223.69] -]) - - -# heliocentric radius coefficients -R0 = np.array([ - [100013989.0, 0.0, 0.0], - [1670700.0, 3.0984635, 6283.07585], - [13956.0, 3.05525, 12566.1517], - [3084.0, 5.1985, 77713.7715], - [1628.0, 1.1739, 5753.3849], - [1576.0, 2.8469, 7860.4194], - [925.0, 5.453, 11506.77], - [542.0, 4.564, 3930.21], - [472.0, 3.661, 5884.927], - [346.0, 0.964, 5507.553], - [329.0, 5.9, 5223.694], - [307.0, 0.299, 5573.143], - [243.0, 4.273, 11790.629], - [212.0, 5.847, 1577.344], - [186.0, 5.022, 10977.079], - [175.0, 3.012, 18849.228], - [110.0, 5.055, 5486.778], - [98.0, 0.89, 6069.78], - [86.0, 5.69, 15720.84], - [86.0, 1.27, 161000.69], - [65.0, 0.27, 17260.15], - [63.0, 0.92, 529.69], - [57.0, 2.01, 83996.85], - [56.0, 5.24, 71430.7], - [49.0, 3.25, 2544.31], - [47.0, 2.58, 775.52], - [45.0, 5.54, 9437.76], - [43.0, 6.01, 6275.96], - [39.0, 5.36, 4694.0], - [38.0, 2.39, 8827.39], - [37.0, 0.83, 19651.05], - [37.0, 4.9, 12139.55], - [36.0, 1.67, 12036.46], - [35.0, 1.84, 2942.46], - [33.0, 0.24, 7084.9], - [32.0, 0.18, 5088.63], - [32.0, 1.78, 398.15], - [28.0, 1.21, 6286.6], - [28.0, 1.9, 6279.55], - [26.0, 4.59, 10447.39] -]) -R1 = np.array([ - [103019.0, 1.10749, 6283.07585], - [1721.0, 1.0644, 12566.1517], - [702.0, 3.142, 0.0], - [32.0, 1.02, 18849.23], - [31.0, 2.84, 5507.55], - [25.0, 1.32, 5223.69], - [18.0, 1.42, 1577.34], - [10.0, 5.91, 10977.08], - [9.0, 1.42, 6275.96], - [9.0, 0.27, 5486.78] -]) -R2 = np.array([ - [4359.0, 5.7846, 6283.0758], - [124.0, 5.579, 12566.152], - [12.0, 3.14, 0.0], - [9.0, 3.63, 77713.77], - [6.0, 1.87, 5573.14], - [3.0, 5.47, 18849.23] -]) -R3 = np.array([ - [145.0, 4.273, 6283.076], - [7.0, 3.92, 12566.15] -]) -R4 = np.array([ - [4.0, 2.56, 6283.08] -]) - - -# longitude and obliquity nutation coefficients -NUTATION_ABCD_ARRAY = np.array([ - [-171996, -174.2, 92025, 8.9], - [-13187, -1.6, 5736, -3.1], - [-2274, -0.2, 977, -0.5], - [2062, 0.2, -895, 0.5], - [1426, -3.4, 54, -0.1], - [712, 0.1, -7, 0], - [-517, 1.2, 224, -0.6], - [-386, -0.4, 200, 0], - [-301, 0, 129, -0.1], - [217, -0.5, -95, 0.3], - [-158, 0, 0, 0], - [129, 0.1, -70, 0], - [123, 0, -53, 0], - [63, 0, 0, 0], - [63, 0.1, -33, 0], - [-59, 0, 26, 0], - [-58, -0.1, 32, 0], - [-51, 0, 27, 0], - [48, 0, 0, 0], - [46, 0, -24, 0], - [-38, 0, 16, 0], - [-31, 0, 13, 0], - [29, 0, 0, 0], - [29, 0, -12, 0], - [26, 0, 0, 0], - [-22, 0, 0, 0], - [21, 0, -10, 0], - [17, -0.1, 0, 0], - [16, 0, -8, 0], - [-16, 0.1, 7, 0], - [-15, 0, 9, 0], - [-13, 0, 7, 0], - [-12, 0, 6, 0], - [11, 0, 0, 0], - [-10, 0, 5, 0], - [-8, 0, 3, 0], - [7, 0, -3, 0], - [-7, 0, 0, 0], - [-7, 0, 3, 0], - [-7, 0, 3, 0], - [6, 0, 0, 0], - [6, 0, -3, 0], - [6, 0, -3, 0], - [-6, 0, 3, 0], - [-6, 0, 3, 0], - [5, 0, 0, 0], - [-5, 0, 3, 0], - [-5, 0, 3, 0], - [-5, 0, 3, 0], - [4, 0, 0, 0], - [4, 0, 0, 0], - [4, 0, 0, 0], - [-4, 0, 0, 0], - [-4, 0, 0, 0], - [-4, 0, 0, 0], - [3, 0, 0, 0], - [-3, 0, 0, 0], - [-3, 0, 0, 0], - [-3, 0, 0, 0], - [-3, 0, 0, 0], - [-3, 0, 0, 0], - [-3, 0, 0, 0], - [-3, 0, 0, 0], -]) - -NUTATION_YTERM_ARRAY = np.array([ - [0, 0, 0, 0, 1], - [-2, 0, 0, 2, 2], - [0, 0, 0, 2, 2], - [0, 0, 0, 0, 2], - [0, 1, 0, 0, 0], - [0, 0, 1, 0, 0], - [-2, 1, 0, 2, 2], - [0, 0, 0, 2, 1], - [0, 0, 1, 2, 2], - [-2, -1, 0, 2, 2], - [-2, 0, 1, 0, 0], - [-2, 0, 0, 2, 1], - [0, 0, -1, 2, 2], - [2, 0, 0, 0, 0], - [0, 0, 1, 0, 1], - [2, 0, -1, 2, 2], - [0, 0, -1, 0, 1], - [0, 0, 1, 2, 1], - [-2, 0, 2, 0, 0], - [0, 0, -2, 2, 1], - [2, 0, 0, 2, 2], - [0, 0, 2, 2, 2], - [0, 0, 2, 0, 0], - [-2, 0, 1, 2, 2], - [0, 0, 0, 2, 0], - [-2, 0, 0, 2, 0], - [0, 0, -1, 2, 1], - [0, 2, 0, 0, 0], - [2, 0, -1, 0, 1], - [-2, 2, 0, 2, 2], - [0, 1, 0, 0, 1], - [-2, 0, 1, 0, 1], - [0, -1, 0, 0, 1], - [0, 0, 2, -2, 0], - [2, 0, -1, 2, 1], - [2, 0, 1, 2, 2], - [0, 1, 0, 2, 2], - [-2, 1, 1, 0, 0], - [0, -1, 0, 2, 2], - [2, 0, 0, 2, 1], - [2, 0, 1, 0, 0], - [-2, 0, 2, 2, 2], - [-2, 0, 1, 2, 1], - [2, 0, -2, 0, 1], - [2, 0, 0, 0, 1], - [0, -1, 1, 0, 0], - [-2, -1, 0, 2, 1], - [-2, 0, 0, 0, 1], - [0, 0, 2, 2, 1], - [-2, 0, 2, 0, 1], - [-2, 1, 0, 2, 1], - [0, 0, 1, -2, 0], - [-1, 0, 1, 0, 0], - [-2, 1, 0, 0, 0], - [1, 0, 0, 0, 0], - [0, 0, 1, 2, 0], - [0, 0, -2, 2, 2], - [-1, -1, 1, 0, 0], - [0, 1, 1, 0, 0], - [0, -1, 1, 2, 2], - [2, -1, -1, 2, 2], - [0, 0, 3, 2, 2], - [2, -1, 0, 2, 2], -]) - - -@jcompile('float64(int64, int64, int64, int64, int64, int64, int64)', - nopython=True) -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 - - -@jcompile('float64(float64)', nopython=True) -def julian_day(unixtime): - jd = unixtime * 1.0 / 86400 + 2440587.5 - return jd - - -@jcompile('float64(float64, float64)', nopython=True) -def julian_ephemeris_day(julian_day, delta_t): - jde = julian_day + delta_t * 1.0 / 86400 - return jde - - -@jcompile('float64(float64)', nopython=True) -def julian_century(julian_day): - jc = (julian_day - 2451545) * 1.0 / 36525 - return jc - - -@jcompile('float64(float64)', nopython=True) -def julian_ephemeris_century(julian_ephemeris_day): - jce = (julian_ephemeris_day - 2451545) * 1.0 / 36525 - return jce - - -@jcompile('float64(float64)', nopython=True) -def julian_ephemeris_millennium(julian_ephemeris_century): - jme = julian_ephemeris_century * 1.0 / 10 - return jme - - -# omit type signature here; specifying read-only arrays requires use of the -# numba.types API, meaning numba must be available to import. -# https://github.com/numba/numba/issues/4511 -@jcompile(nopython=True) -def sum_mult_cos_add_mult(arr, x): - # shared calculation used for heliocentric longitude, latitude, and radius - s = 0. - for row in range(arr.shape[0]): - s += arr[row, 0] * np.cos(arr[row, 1] + arr[row, 2] * x) - return s - -@jcompile('float64(float64)', nopython=True) -def heliocentric_longitude(jme): - l0 = sum_mult_cos_add_mult(L0, jme) - l1 = sum_mult_cos_add_mult(L1, jme) - l2 = sum_mult_cos_add_mult(L2, jme) - l3 = sum_mult_cos_add_mult(L3, jme) - l4 = sum_mult_cos_add_mult(L4, jme) - l5 = sum_mult_cos_add_mult(L5, jme) - - l_rad = (l0 + l1 * jme + l2 * jme**2 + l3 * jme**3 + l4 * jme**4 + - l5 * jme**5)/10**8 - l = np.rad2deg(l_rad) - return l % 360 - -@jcompile('float64(float64)', nopython=True) -def heliocentric_latitude(jme): - b0 = sum_mult_cos_add_mult(B0, jme) - b1 = sum_mult_cos_add_mult(B1, jme) - - b_rad = (b0 + b1 * jme)/10**8 - b = np.rad2deg(b_rad) - return b - - -@jcompile('float64(float64)', nopython=True) -def heliocentric_radius_vector(jme): - r0 = sum_mult_cos_add_mult(R0, jme) - r1 = sum_mult_cos_add_mult(R1, jme) - r2 = sum_mult_cos_add_mult(R2, jme) - r3 = sum_mult_cos_add_mult(R3, jme) - r4 = sum_mult_cos_add_mult(R4, jme) - - r = (r0 + r1 * jme + r2 * jme**2 + r3 * jme**3 + r4 * jme**4)/10**8 - return r - - -@jcompile('float64(float64)', nopython=True) -def geocentric_longitude(heliocentric_longitude): - theta = heliocentric_longitude + 180.0 - return theta % 360 - - -@jcompile('float64(float64)', nopython=True) -def geocentric_latitude(heliocentric_latitude): - beta = -1.0*heliocentric_latitude - return beta - - -@jcompile('float64(float64)', nopython=True) -def mean_elongation(julian_ephemeris_century): - x0 = (297.85036 - + 445267.111480 * julian_ephemeris_century - - 0.0019142 * julian_ephemeris_century**2 - + julian_ephemeris_century**3 / 189474) - return x0 - - -@jcompile('float64(float64)', nopython=True) -def mean_anomaly_sun(julian_ephemeris_century): - x1 = (357.52772 - + 35999.050340 * julian_ephemeris_century - - 0.0001603 * julian_ephemeris_century**2 - - julian_ephemeris_century**3 / 300000) - return x1 - - -@jcompile('float64(float64)', nopython=True) -def mean_anomaly_moon(julian_ephemeris_century): - x2 = (134.96298 - + 477198.867398 * julian_ephemeris_century - + 0.0086972 * julian_ephemeris_century**2 - + julian_ephemeris_century**3 / 56250) - return x2 - - -@jcompile('float64(float64)', nopython=True) -def moon_argument_latitude(julian_ephemeris_century): - x3 = (93.27191 - + 483202.017538 * julian_ephemeris_century - - 0.0036825 * julian_ephemeris_century**2 - + julian_ephemeris_century**3 / 327270) - return x3 - - -@jcompile('float64(float64)', nopython=True) -def moon_ascending_longitude(julian_ephemeris_century): - x4 = (125.04452 - - 1934.136261 * julian_ephemeris_century - + 0.0020708 * julian_ephemeris_century**2 - + julian_ephemeris_century**3 / 450000) - return x4 - - -@jcompile( - 'void(float64, float64, float64, float64, float64, float64, float64[:])', - nopython=True) -def longitude_obliquity_nutation(julian_ephemeris_century, x0, x1, x2, x3, x4, - out): - delta_psi_sum = 0.0 - delta_eps_sum = 0.0 - for row in range(NUTATION_YTERM_ARRAY.shape[0]): - a = NUTATION_ABCD_ARRAY[row, 0] - b = NUTATION_ABCD_ARRAY[row, 1] - c = NUTATION_ABCD_ARRAY[row, 2] - d = NUTATION_ABCD_ARRAY[row, 3] - arg = np.radians( - NUTATION_YTERM_ARRAY[row, 0]*x0 + - NUTATION_YTERM_ARRAY[row, 1]*x1 + - NUTATION_YTERM_ARRAY[row, 2]*x2 + - NUTATION_YTERM_ARRAY[row, 3]*x3 + - NUTATION_YTERM_ARRAY[row, 4]*x4 - ) - delta_psi_sum += (a + b * julian_ephemeris_century) * np.sin(arg) - delta_eps_sum += (c + d * julian_ephemeris_century) * np.cos(arg) - delta_psi = delta_psi_sum*1.0/36000000 - delta_eps = delta_eps_sum*1.0/36000000 - # seems like we ought to be able to return a tuple here instead - # of resorting to `out`, but returning a UniTuple from this - # function caused calculations elsewhere to give the wrong result. - # very difficult to investigate since it did not occur when using - # object mode. issue was observed on numba 0.56.4 - out[0] = delta_psi - out[1] = delta_eps - - -@jcompile('float64(float64)', nopython=True) -def mean_ecliptic_obliquity(julian_ephemeris_millennium): - U = 1.0*julian_ephemeris_millennium/10 - e0 = (84381.448 - 4680.93 * U - 1.55 * U**2 - + 1999.25 * U**3 - 51.38 * U**4 - 249.67 * U**5 - - 39.05 * U**6 + 7.12 * U**7 + 27.87 * U**8 - + 5.79 * U**9 + 2.45 * U**10) - return e0 - - -@jcompile('float64(float64, float64)', nopython=True) -def true_ecliptic_obliquity(mean_ecliptic_obliquity, obliquity_nutation): - e0 = mean_ecliptic_obliquity - deleps = obliquity_nutation - e = e0*1.0/3600 + deleps - return e - - -@jcompile('float64(float64)', nopython=True) -def aberration_correction(earth_radius_vector): - deltau = -20.4898 / (3600 * earth_radius_vector) - return deltau - - -@jcompile('float64(float64, float64, float64)', nopython=True) -def apparent_sun_longitude(geocentric_longitude, longitude_nutation, - aberration_correction): - lamd = geocentric_longitude + longitude_nutation + aberration_correction - return lamd - - -@jcompile('float64(float64, float64)', nopython=True) -def mean_sidereal_time(julian_day, julian_century): - v0 = (280.46061837 + 360.98564736629 * (julian_day - 2451545) - + 0.000387933 * julian_century**2 - julian_century**3 / 38710000) - return v0 % 360.0 - - -@jcompile('float64(float64, float64, float64)', nopython=True) -def apparent_sidereal_time(mean_sidereal_time, longitude_nutation, - true_ecliptic_obliquity): - v = mean_sidereal_time + longitude_nutation * np.cos( - np.radians(true_ecliptic_obliquity)) - return v - - -@jcompile('float64(float64, float64, float64)', nopython=True) -def geocentric_sun_right_ascension(apparent_sun_longitude, - true_ecliptic_obliquity, - geocentric_latitude): - true_ecliptic_obliquity_rad = np.radians(true_ecliptic_obliquity) - apparent_sun_longitude_rad = np.radians(apparent_sun_longitude) - - num = (np.sin(apparent_sun_longitude_rad) - * np.cos(true_ecliptic_obliquity_rad) - - np.tan(np.radians(geocentric_latitude)) - * np.sin(true_ecliptic_obliquity_rad)) - alpha = np.degrees(np.arctan2(num, np.cos(apparent_sun_longitude_rad))) - return alpha % 360 - - -@jcompile('float64(float64, float64, float64)', nopython=True) -def geocentric_sun_declination(apparent_sun_longitude, true_ecliptic_obliquity, - geocentric_latitude): - geocentric_latitude_rad = np.radians(geocentric_latitude) - true_ecliptic_obliquity_rad = np.radians(true_ecliptic_obliquity) - - delta = np.degrees(np.arcsin(np.sin(geocentric_latitude_rad) * - np.cos(true_ecliptic_obliquity_rad) + - np.cos(geocentric_latitude_rad) * - np.sin(true_ecliptic_obliquity_rad) * - np.sin(np.radians(apparent_sun_longitude)))) - return delta - - -@jcompile('float64(float64, float64, float64)', nopython=True) -def local_hour_angle(apparent_sidereal_time, observer_longitude, - sun_right_ascension): - """Measured westward from south""" - H = apparent_sidereal_time + observer_longitude - sun_right_ascension - return H % 360 - - -@jcompile('float64(float64)', nopython=True) -def equatorial_horizontal_parallax(earth_radius_vector): - xi = 8.794 / (3600 * earth_radius_vector) - return xi - - -@jcompile('float64(float64)', nopython=True) -def uterm(observer_latitude): - u = np.arctan(0.99664719 * np.tan(np.radians(observer_latitude))) - return u - - -@jcompile('float64(float64, float64, float64)', nopython=True) -def xterm(u, observer_latitude, observer_elevation): - x = (np.cos(u) + observer_elevation / 6378140 - * np.cos(np.radians(observer_latitude))) - return x - - -@jcompile('float64(float64, float64, float64)', nopython=True) -def yterm(u, observer_latitude, observer_elevation): - y = (0.99664719 * np.sin(u) + observer_elevation / 6378140 - * np.sin(np.radians(observer_latitude))) - return y - - -@jcompile('float64(float64, float64,float64, float64)', nopython=True) -def parallax_sun_right_ascension(xterm, equatorial_horizontal_parallax, - local_hour_angle, geocentric_sun_declination): - equatorial_horizontal_parallax_rad = \ - np.radians(equatorial_horizontal_parallax) - local_hour_angle_rad = np.radians(local_hour_angle) - - num = (-xterm * np.sin(equatorial_horizontal_parallax_rad) - * np.sin(local_hour_angle_rad)) - denom = (np.cos(np.radians(geocentric_sun_declination)) - - xterm * np.sin(equatorial_horizontal_parallax_rad) - * np.cos(local_hour_angle_rad)) - delta_alpha = np.degrees(np.arctan2(num, denom)) - return delta_alpha - - -@jcompile('float64(float64, float64)', nopython=True) -def topocentric_sun_right_ascension(geocentric_sun_right_ascension, - parallax_sun_right_ascension): - alpha_prime = geocentric_sun_right_ascension + parallax_sun_right_ascension - return alpha_prime - - -@jcompile('float64(float64, float64, float64, float64, float64, float64)', - nopython=True) -def topocentric_sun_declination(geocentric_sun_declination, xterm, yterm, - equatorial_horizontal_parallax, - parallax_sun_right_ascension, - local_hour_angle): - geocentric_sun_declination_rad = np.radians(geocentric_sun_declination) - equatorial_horizontal_parallax_rad = \ - np.radians(equatorial_horizontal_parallax) - - num = ((np.sin(geocentric_sun_declination_rad) - yterm - * np.sin(equatorial_horizontal_parallax_rad)) - * np.cos(np.radians(parallax_sun_right_ascension))) - denom = (np.cos(geocentric_sun_declination_rad) - xterm - * np.sin(equatorial_horizontal_parallax_rad) - * np.cos(np.radians(local_hour_angle))) - delta = np.degrees(np.arctan2(num, denom)) - return delta - - -@jcompile('float64(float64, float64)', nopython=True) -def topocentric_local_hour_angle(local_hour_angle, - parallax_sun_right_ascension): - H_prime = local_hour_angle - parallax_sun_right_ascension - return H_prime - - -@jcompile('float64(float64, float64, float64)', nopython=True) -def topocentric_elevation_angle_without_atmosphere(observer_latitude, - topocentric_sun_declination, - topocentric_local_hour_angle - ): - observer_latitude_rad = np.radians(observer_latitude) - topocentric_sun_declination_rad = np.radians(topocentric_sun_declination) - e0 = np.degrees(np.arcsin( - np.sin(observer_latitude_rad) - * np.sin(topocentric_sun_declination_rad) - + np.cos(observer_latitude_rad) - * np.cos(topocentric_sun_declination_rad) - * np.cos(np.radians(topocentric_local_hour_angle)))) - return e0 - - -@jcompile('float64(float64, float64, float64, float64)', nopython=True) -def atmospheric_refraction_correction(local_pressure, local_temp, - topocentric_elevation_angle_wo_atmosphere, - atmos_refract): - # switch sets delta_e when the sun is below the horizon - switch = topocentric_elevation_angle_wo_atmosphere >= -1.0 * ( - 0.26667 + atmos_refract) - delta_e = ((local_pressure / 1010.0) * (283.0 / (273 + local_temp)) - * 1.02 / (60 * np.tan(np.radians( - topocentric_elevation_angle_wo_atmosphere - + 10.3 / (topocentric_elevation_angle_wo_atmosphere - + 5.11))))) * switch - return delta_e - - -@jcompile('float64(float64, float64)', nopython=True) -def topocentric_elevation_angle(topocentric_elevation_angle_without_atmosphere, - atmospheric_refraction_correction): - e = (topocentric_elevation_angle_without_atmosphere - + atmospheric_refraction_correction) - return e - - -@jcompile('float64(float64)', nopython=True) -def topocentric_zenith_angle(topocentric_elevation_angle): - theta = 90 - topocentric_elevation_angle - return theta - - -@jcompile('float64(float64, float64, float64)', nopython=True) -def topocentric_astronomers_azimuth(topocentric_local_hour_angle, - topocentric_sun_declination, - observer_latitude): - topocentric_local_hour_angle_rad = np.radians(topocentric_local_hour_angle) - observer_latitude_rad = np.radians(observer_latitude) - num = np.sin(topocentric_local_hour_angle_rad) - denom = (np.cos(topocentric_local_hour_angle_rad) - * np.sin(observer_latitude_rad) - - np.tan(np.radians(topocentric_sun_declination)) - * np.cos(observer_latitude_rad)) - gamma = np.degrees(np.arctan2(num, denom)) - return gamma % 360 - - -@jcompile('float64(float64)', nopython=True) -def topocentric_azimuth_angle(topocentric_astronomers_azimuth): - phi = topocentric_astronomers_azimuth + 180 - return phi % 360 - - -@jcompile('float64(float64)', nopython=True) -def sun_mean_longitude(julian_ephemeris_millennium): - M = (280.4664567 + 360007.6982779 * julian_ephemeris_millennium - + 0.03032028 * julian_ephemeris_millennium**2 - + julian_ephemeris_millennium**3 / 49931 - - julian_ephemeris_millennium**4 / 15300 - - julian_ephemeris_millennium**5 / 2000000) - return M - - -@jcompile('float64(float64, float64, float64, float64)', nopython=True) -def equation_of_time(sun_mean_longitude, geocentric_sun_right_ascension, - longitude_nutation, true_ecliptic_obliquity): - E = (sun_mean_longitude - 0.0057183 - geocentric_sun_right_ascension + - longitude_nutation * np.cos(np.radians(true_ecliptic_obliquity))) - # limit between 0 and 360 - E = E % 360 - # convert to minutes - E *= 4 - greater = E > 20 - less = E < -20 - other = (E <= 20) & (E >= -20) - E = greater * (E - 1440) + less * (E + 1440) + other * E - return E - - -@jcompile('void(float64[:], float64[:], float64[:], float64[:,:])', - nopython=True, nogil=True) -def solar_position_loop(unixtime, delta_t, loc_args, out): - """Loop through the time array and calculate the solar position""" - lat = loc_args[0] - lon = loc_args[1] - elev = loc_args[2] - pressure = loc_args[3] - temp = loc_args[4] - atmos_refract = loc_args[5] - sst = loc_args[6] - esd = loc_args[7] - - for i in range(unixtime.shape[0]): - utime = unixtime[i] - dT = delta_t[i] - jd = julian_day(utime) - jde = julian_ephemeris_day(jd, dT) - jc = julian_century(jd) - jce = julian_ephemeris_century(jde) - jme = julian_ephemeris_millennium(jce) - R = heliocentric_radius_vector(jme) - if esd: - out[0, i] = R - continue - L = heliocentric_longitude(jme) - B = heliocentric_latitude(jme) - Theta = geocentric_longitude(L) - beta = geocentric_latitude(B) - x0 = mean_elongation(jce) - x1 = mean_anomaly_sun(jce) - x2 = mean_anomaly_moon(jce) - x3 = moon_argument_latitude(jce) - x4 = moon_ascending_longitude(jce) - l_o_nutation = np.empty((2,)) - longitude_obliquity_nutation(jce, x0, x1, x2, x3, x4, l_o_nutation) - delta_psi = l_o_nutation[0] - delta_epsilon = l_o_nutation[1] - epsilon0 = mean_ecliptic_obliquity(jme) - epsilon = true_ecliptic_obliquity(epsilon0, delta_epsilon) - delta_tau = aberration_correction(R) - lamd = apparent_sun_longitude(Theta, delta_psi, delta_tau) - v0 = mean_sidereal_time(jd, jc) - v = apparent_sidereal_time(v0, delta_psi, epsilon) - alpha = geocentric_sun_right_ascension(lamd, epsilon, beta) - delta = geocentric_sun_declination(lamd, epsilon, beta) - if sst: - out[0, i] = v - out[1, i] = alpha - out[2, i] = delta - continue - m = sun_mean_longitude(jme) - eot = equation_of_time(m, alpha, delta_psi, epsilon) - H = local_hour_angle(v, lon, alpha) - xi = equatorial_horizontal_parallax(R) - u = uterm(lat) - x = xterm(u, lat, elev) - y = yterm(u, lat, elev) - delta_alpha = parallax_sun_right_ascension(x, xi, H, delta) - delta_prime = topocentric_sun_declination(delta, x, y, xi, delta_alpha, - H) - H_prime = topocentric_local_hour_angle(H, delta_alpha) - e0 = topocentric_elevation_angle_without_atmosphere(lat, delta_prime, - H_prime) - delta_e = atmospheric_refraction_correction(pressure, temp, e0, - atmos_refract) - e = topocentric_elevation_angle(e0, delta_e) - theta = topocentric_zenith_angle(e) - theta0 = topocentric_zenith_angle(e0) - gamma = topocentric_astronomers_azimuth(H_prime, delta_prime, lat) - phi = topocentric_azimuth_angle(gamma) - out[0, i] = theta - out[1, i] = theta0 - out[2, i] = e - out[3, i] = e0 - out[4, i] = phi - out[5, i] = eot - - -def solar_position_numba(unixtime, lat, lon, elev, pressure, temp, delta_t, - atmos_refract, numthreads, sst=False, esd=False): - """Calculate the solar position using the numba compiled functions - and multiple threads. Very slow if functions are not numba compiled. - """ - # these args are the same for each thread - loc_args = np.array([lat, lon, elev, pressure, temp, - atmos_refract, sst, esd], dtype=np.float64) - - # turn delta_t into an array if it isn't already - delta_t = np.full_like(unixtime, delta_t, dtype=np.float64) - - # construct dims x ulength array to put the results in - ulength = unixtime.shape[0] - if sst: - dims = 3 - elif esd: - dims = 1 - else: - dims = 6 - result = np.empty((dims, ulength), dtype=np.float64) - - if unixtime.dtype != np.float64: - unixtime = unixtime.astype(np.float64) - - if ulength < numthreads: - numthreads = ulength - - if numthreads <= 1: - solar_position_loop(unixtime, delta_t, loc_args, result) - return result - - # split the input and output arrays into numthreads chunks - split0 = np.array_split(unixtime, numthreads) - split1 = np.array_split(delta_t, numthreads) - split2 = np.array_split(result, numthreads, axis=1) - chunks = [ - [a0, a1, loc_args, a2] - for a0, a1, a2 in zip(split0, split1, split2) - ] - # Spawn one thread per chunk - threads = [threading.Thread(target=solar_position_loop, args=chunk) - for chunk in chunks] - for thread in threads: - thread.start() - for thread in threads: - thread.join() - return result - - -def solar_position_numpy(unixtime, lat, lon, elev, pressure, temp, delta_t, - atmos_refract, numthreads, sst=False, esd=False): - """Calculate the solar position assuming unixtime is a numpy array. Note - this function will not work if the solar position functions were - compiled with numba. - """ - - jd = julian_day(unixtime) - jde = julian_ephemeris_day(jd, delta_t) - jc = julian_century(jd) - jce = julian_ephemeris_century(jde) - jme = julian_ephemeris_millennium(jce) - R = heliocentric_radius_vector(jme) - if esd: - return (R, ) - L = heliocentric_longitude(jme) - B = heliocentric_latitude(jme) - Theta = geocentric_longitude(L) - beta = geocentric_latitude(B) - x0 = mean_elongation(jce) - x1 = mean_anomaly_sun(jce) - x2 = mean_anomaly_moon(jce) - x3 = moon_argument_latitude(jce) - x4 = moon_ascending_longitude(jce) - l_o_nutation = np.empty((2, len(x0))) - longitude_obliquity_nutation(jce, x0, x1, x2, x3, x4, l_o_nutation) - delta_psi = l_o_nutation[0] - delta_epsilon = l_o_nutation[1] - epsilon0 = mean_ecliptic_obliquity(jme) - epsilon = true_ecliptic_obliquity(epsilon0, delta_epsilon) - delta_tau = aberration_correction(R) - lamd = apparent_sun_longitude(Theta, delta_psi, delta_tau) - v0 = mean_sidereal_time(jd, jc) - v = apparent_sidereal_time(v0, delta_psi, epsilon) - alpha = geocentric_sun_right_ascension(lamd, epsilon, beta) - delta = geocentric_sun_declination(lamd, epsilon, beta) - if sst: - return v, alpha, delta - m = sun_mean_longitude(jme) - eot = equation_of_time(m, alpha, delta_psi, epsilon) - H = local_hour_angle(v, lon, alpha) - xi = equatorial_horizontal_parallax(R) - u = uterm(lat) - x = xterm(u, lat, elev) - y = yterm(u, lat, elev) - delta_alpha = parallax_sun_right_ascension(x, xi, H, delta) - delta_prime = topocentric_sun_declination(delta, x, y, xi, delta_alpha, H) - H_prime = topocentric_local_hour_angle(H, delta_alpha) - e0 = topocentric_elevation_angle_without_atmosphere(lat, delta_prime, - H_prime) - delta_e = atmospheric_refraction_correction(pressure, temp, e0, - atmos_refract) - e = topocentric_elevation_angle(e0, delta_e) - theta = topocentric_zenith_angle(e) - theta0 = topocentric_zenith_angle(e0) - gamma = topocentric_astronomers_azimuth(H_prime, delta_prime, lat) - phi = topocentric_azimuth_angle(gamma) - return theta, theta0, e, e0, phi, eot - - -def solar_position(unixtime, lat, lon, elev, pressure, temp, delta_t, - atmos_refract, numthreads=8, sst=False, esd=False): - - """ - Calculate the solar position using the - NREL SPA described in [1]_, [2]_. - - If numba is installed, the functions can be compiled - and the code runs quickly. If not, the functions - still evaluate but use numpy instead. - - Parameters - ---------- - unixtime : numpy array - Array of unix/epoch timestamps to calculate solar position for. - Unixtime is the number of seconds since Jan. 1, 1970 00:00:00 UTC. - A pandas.DatetimeIndex is easily converted using .view(np.int64)/10**9 - lat : float - Latitude to calculate solar position for - lon : float - Longitude to calculate solar position for - elev : float - Elevation of location in meters - pressure : int or float - avg. yearly pressure at location in millibars; - used for atmospheric correction - temp : int or float - avg. yearly temperature at location in - degrees C; used for atmospheric correction - delta_t : float or array - Difference between terrestrial time and UT1. - atmos_refract : float - The approximate atmospheric refraction (in degrees) - at sunrise and sunset. - numthreads: int, optional, default 8 - Number of threads to use for computation if numba>=0.17 - is installed. - sst : bool, default False - If True, return only data needed for sunrise, sunset, and transit - calculations. - esd : bool, default False - If True, return only Earth-Sun distance in AU - - Returns - ------- - Numpy Array with elements: - apparent zenith, - zenith, - apparent_elevation, - elevation, - azimuth, - equation_of_time - - References - ---------- - [1] I. Reda and A. Andreas, Solar position algorithm for solar radiation - applications. Solar Energy, vol. 76, no. 5, pp. 577-589, 2004. - - [2] I. Reda and A. Andreas, Corrigendum to Solar position algorithm for - solar radiation applications. Solar Energy, vol. 81, no. 6, p. 838, 2007. - """ - if USE_NUMBA: - do_calc = solar_position_numba - else: - do_calc = solar_position_numpy - - result = do_calc(unixtime, lat, lon, elev, pressure, - temp, delta_t, atmos_refract, numthreads, - sst, esd) - - if not isinstance(result, np.ndarray): - try: - result = np.array(result) - except Exception: - pass - - return result - - -def transit_sunrise_sunset(dates, lat, lon, delta_t, numthreads): - """ - Calculate the sun transit, sunrise, and sunset - for a set of dates at a given location using the NREL SPA [1]_, [2]_. - - Parameters - ---------- - dates : array - Numpy array of ints/floats corresponding to the Unix time - for the dates of interest, must be midnight UTC (00:00+00:00) - on the day of interest. - lat : float - Latitude of location to perform calculation for - lon : float - Longitude of location - delta_t : float or array - Difference between terrestrial time and UT. USNO has tables. - numthreads : int - Number to threads to use for calculation (if using numba) - - Returns - ------- - tuple : (transit, sunrise, sunset) localized to UTC - - References - ---------- - .. [1] I. Reda and A. Andreas, Solar position algorithm for solar - radiation applications. Solar Energy, vol. 76, no. 5, pp. 577-589, - 2004. :doi:`10.1016/j.solener.2003.12.003` +""" +Calculate the sun's position using the NREL Solar Position Algorithm +(SPA) either using numpy arrays or compiling the code to machine language +with numba. +""" + +# Contributors: +# Created by Tony Lorenzo (@alorenzo175), Univ. of Arizona, 2015 + +import os +import threading +import warnings + +import numpy as np + + +# this block is a way to use an environment variable to switch between +# compiling the functions with numba or just use numpy +def nocompile(*args, **kwargs): + return lambda func: func + + +if os.getenv('PVLIB_USE_NUMBA', '0') != '0': + try: + from numba import jit + except ImportError: + warnings.warn('Could not import numba, falling back to numpy ' + + 'calculation') + jcompile = nocompile + USE_NUMBA = False + else: + jcompile = jit + USE_NUMBA = True +else: + jcompile = nocompile + USE_NUMBA = False + + +# heliocentric longitude coefficients +L0 = np.array([ + [175347046.0, 0.0, 0.0], + [3341656.0, 4.6692568, 6283.07585], + [34894.0, 4.6261, 12566.1517], + [3497.0, 2.7441, 5753.3849], + [3418.0, 2.8289, 3.5231], + [3136.0, 3.6277, 77713.7715], + [2676.0, 4.4181, 7860.4194], + [2343.0, 6.1352, 3930.2097], + [1324.0, 0.7425, 11506.7698], + [1273.0, 2.0371, 529.691], + [1199.0, 1.1096, 1577.3435], + [990.0, 5.233, 5884.927], + [902.0, 2.045, 26.298], + [857.0, 3.508, 398.149], + [780.0, 1.179, 5223.694], + [753.0, 2.533, 5507.553], + [505.0, 4.583, 18849.228], + [492.0, 4.205, 775.523], + [357.0, 2.92, 0.067], + [317.0, 5.849, 11790.629], + [284.0, 1.899, 796.298], + [271.0, 0.315, 10977.079], + [243.0, 0.345, 5486.778], + [206.0, 4.806, 2544.314], + [205.0, 1.869, 5573.143], + [202.0, 2.458, 6069.777], + [156.0, 0.833, 213.299], + [132.0, 3.411, 2942.463], + [126.0, 1.083, 20.775], + [115.0, 0.645, 0.98], + [103.0, 0.636, 4694.003], + [102.0, 0.976, 15720.839], + [102.0, 4.267, 7.114], + [99.0, 6.21, 2146.17], + [98.0, 0.68, 155.42], + [86.0, 5.98, 161000.69], + [85.0, 1.3, 6275.96], + [85.0, 3.67, 71430.7], + [80.0, 1.81, 17260.15], + [79.0, 3.04, 12036.46], + [75.0, 1.76, 5088.63], + [74.0, 3.5, 3154.69], + [74.0, 4.68, 801.82], + [70.0, 0.83, 9437.76], + [62.0, 3.98, 8827.39], + [61.0, 1.82, 7084.9], + [57.0, 2.78, 6286.6], + [56.0, 4.39, 14143.5], + [56.0, 3.47, 6279.55], + [52.0, 0.19, 12139.55], + [52.0, 1.33, 1748.02], + [51.0, 0.28, 5856.48], + [49.0, 0.49, 1194.45], + [41.0, 5.37, 8429.24], + [41.0, 2.4, 19651.05], + [39.0, 6.17, 10447.39], + [37.0, 6.04, 10213.29], + [37.0, 2.57, 1059.38], + [36.0, 1.71, 2352.87], + [36.0, 1.78, 6812.77], + [33.0, 0.59, 17789.85], + [30.0, 0.44, 83996.85], + [30.0, 2.74, 1349.87], + [25.0, 3.16, 4690.48] +]) +L1 = np.array([ + [628331966747.0, 0.0, 0.0], + [206059.0, 2.678235, 6283.07585], + [4303.0, 2.6351, 12566.1517], + [425.0, 1.59, 3.523], + [119.0, 5.796, 26.298], + [109.0, 2.966, 1577.344], + [93.0, 2.59, 18849.23], + [72.0, 1.14, 529.69], + [68.0, 1.87, 398.15], + [67.0, 4.41, 5507.55], + [59.0, 2.89, 5223.69], + [56.0, 2.17, 155.42], + [45.0, 0.4, 796.3], + [36.0, 0.47, 775.52], + [29.0, 2.65, 7.11], + [21.0, 5.34, 0.98], + [19.0, 1.85, 5486.78], + [19.0, 4.97, 213.3], + [17.0, 2.99, 6275.96], + [16.0, 0.03, 2544.31], + [16.0, 1.43, 2146.17], + [15.0, 1.21, 10977.08], + [12.0, 2.83, 1748.02], + [12.0, 3.26, 5088.63], + [12.0, 5.27, 1194.45], + [12.0, 2.08, 4694.0], + [11.0, 0.77, 553.57], + [10.0, 1.3, 6286.6], + [10.0, 4.24, 1349.87], + [9.0, 2.7, 242.73], + [9.0, 5.64, 951.72], + [8.0, 5.3, 2352.87], + [6.0, 2.65, 9437.76], + [6.0, 4.67, 4690.48] +]) +L2 = np.array([ + [52919.0, 0.0, 0.0], + [8720.0, 1.0721, 6283.0758], + [309.0, 0.867, 12566.152], + [27.0, 0.05, 3.52], + [16.0, 5.19, 26.3], + [16.0, 3.68, 155.42], + [10.0, 0.76, 18849.23], + [9.0, 2.06, 77713.77], + [7.0, 0.83, 775.52], + [5.0, 4.66, 1577.34], + [4.0, 1.03, 7.11], + [4.0, 3.44, 5573.14], + [3.0, 5.14, 796.3], + [3.0, 6.05, 5507.55], + [3.0, 1.19, 242.73], + [3.0, 6.12, 529.69], + [3.0, 0.31, 398.15], + [3.0, 2.28, 553.57], + [2.0, 4.38, 5223.69], + [2.0, 3.75, 0.98] +]) +L3 = np.array([ + [289.0, 5.844, 6283.076], + [35.0, 0.0, 0.0], + [17.0, 5.49, 12566.15], + [3.0, 5.2, 155.42], + [1.0, 4.72, 3.52], + [1.0, 5.3, 18849.23], + [1.0, 5.97, 242.73] +]) +L4 = np.array([ + [114.0, 3.142, 0.0], + [8.0, 4.13, 6283.08], + [1.0, 3.84, 12566.15] +]) +L5 = np.array([ + [1.0, 3.14, 0.0] +]) + + +# heliocentric latitude coefficients +B0 = np.array([ + [280.0, 3.199, 84334.662], + [102.0, 5.422, 5507.553], + [80.0, 3.88, 5223.69], + [44.0, 3.7, 2352.87], + [32.0, 4.0, 1577.34] +]) +B1 = np.array([ + [9.0, 3.9, 5507.55], + [6.0, 1.73, 5223.69] +]) + + +# heliocentric radius coefficients +R0 = np.array([ + [100013989.0, 0.0, 0.0], + [1670700.0, 3.0984635, 6283.07585], + [13956.0, 3.05525, 12566.1517], + [3084.0, 5.1985, 77713.7715], + [1628.0, 1.1739, 5753.3849], + [1576.0, 2.8469, 7860.4194], + [925.0, 5.453, 11506.77], + [542.0, 4.564, 3930.21], + [472.0, 3.661, 5884.927], + [346.0, 0.964, 5507.553], + [329.0, 5.9, 5223.694], + [307.0, 0.299, 5573.143], + [243.0, 4.273, 11790.629], + [212.0, 5.847, 1577.344], + [186.0, 5.022, 10977.079], + [175.0, 3.012, 18849.228], + [110.0, 5.055, 5486.778], + [98.0, 0.89, 6069.78], + [86.0, 5.69, 15720.84], + [86.0, 1.27, 161000.69], + [65.0, 0.27, 17260.15], + [63.0, 0.92, 529.69], + [57.0, 2.01, 83996.85], + [56.0, 5.24, 71430.7], + [49.0, 3.25, 2544.31], + [47.0, 2.58, 775.52], + [45.0, 5.54, 9437.76], + [43.0, 6.01, 6275.96], + [39.0, 5.36, 4694.0], + [38.0, 2.39, 8827.39], + [37.0, 0.83, 19651.05], + [37.0, 4.9, 12139.55], + [36.0, 1.67, 12036.46], + [35.0, 1.84, 2942.46], + [33.0, 0.24, 7084.9], + [32.0, 0.18, 5088.63], + [32.0, 1.78, 398.15], + [28.0, 1.21, 6286.6], + [28.0, 1.9, 6279.55], + [26.0, 4.59, 10447.39] +]) +R1 = np.array([ + [103019.0, 1.10749, 6283.07585], + [1721.0, 1.0644, 12566.1517], + [702.0, 3.142, 0.0], + [32.0, 1.02, 18849.23], + [31.0, 2.84, 5507.55], + [25.0, 1.32, 5223.69], + [18.0, 1.42, 1577.34], + [10.0, 5.91, 10977.08], + [9.0, 1.42, 6275.96], + [9.0, 0.27, 5486.78] +]) +R2 = np.array([ + [4359.0, 5.7846, 6283.0758], + [124.0, 5.579, 12566.152], + [12.0, 3.14, 0.0], + [9.0, 3.63, 77713.77], + [6.0, 1.87, 5573.14], + [3.0, 5.47, 18849.23] +]) +R3 = np.array([ + [145.0, 4.273, 6283.076], + [7.0, 3.92, 12566.15] +]) +R4 = np.array([ + [4.0, 2.56, 6283.08] +]) + + +# longitude and obliquity nutation coefficients +NUTATION_ABCD_ARRAY = np.array([ + [-171996, -174.2, 92025, 8.9], + [-13187, -1.6, 5736, -3.1], + [-2274, -0.2, 977, -0.5], + [2062, 0.2, -895, 0.5], + [1426, -3.4, 54, -0.1], + [712, 0.1, -7, 0], + [-517, 1.2, 224, -0.6], + [-386, -0.4, 200, 0], + [-301, 0, 129, -0.1], + [217, -0.5, -95, 0.3], + [-158, 0, 0, 0], + [129, 0.1, -70, 0], + [123, 0, -53, 0], + [63, 0, 0, 0], + [63, 0.1, -33, 0], + [-59, 0, 26, 0], + [-58, -0.1, 32, 0], + [-51, 0, 27, 0], + [48, 0, 0, 0], + [46, 0, -24, 0], + [-38, 0, 16, 0], + [-31, 0, 13, 0], + [29, 0, 0, 0], + [29, 0, -12, 0], + [26, 0, 0, 0], + [-22, 0, 0, 0], + [21, 0, -10, 0], + [17, -0.1, 0, 0], + [16, 0, -8, 0], + [-16, 0.1, 7, 0], + [-15, 0, 9, 0], + [-13, 0, 7, 0], + [-12, 0, 6, 0], + [11, 0, 0, 0], + [-10, 0, 5, 0], + [-8, 0, 3, 0], + [7, 0, -3, 0], + [-7, 0, 0, 0], + [-7, 0, 3, 0], + [-7, 0, 3, 0], + [6, 0, 0, 0], + [6, 0, -3, 0], + [6, 0, -3, 0], + [-6, 0, 3, 0], + [-6, 0, 3, 0], + [5, 0, 0, 0], + [-5, 0, 3, 0], + [-5, 0, 3, 0], + [-5, 0, 3, 0], + [4, 0, 0, 0], + [4, 0, 0, 0], + [4, 0, 0, 0], + [-4, 0, 0, 0], + [-4, 0, 0, 0], + [-4, 0, 0, 0], + [3, 0, 0, 0], + [-3, 0, 0, 0], + [-3, 0, 0, 0], + [-3, 0, 0, 0], + [-3, 0, 0, 0], + [-3, 0, 0, 0], + [-3, 0, 0, 0], + [-3, 0, 0, 0], +]) + +NUTATION_YTERM_ARRAY = np.array([ + [0, 0, 0, 0, 1], + [-2, 0, 0, 2, 2], + [0, 0, 0, 2, 2], + [0, 0, 0, 0, 2], + [0, 1, 0, 0, 0], + [0, 0, 1, 0, 0], + [-2, 1, 0, 2, 2], + [0, 0, 0, 2, 1], + [0, 0, 1, 2, 2], + [-2, -1, 0, 2, 2], + [-2, 0, 1, 0, 0], + [-2, 0, 0, 2, 1], + [0, 0, -1, 2, 2], + [2, 0, 0, 0, 0], + [0, 0, 1, 0, 1], + [2, 0, -1, 2, 2], + [0, 0, -1, 0, 1], + [0, 0, 1, 2, 1], + [-2, 0, 2, 0, 0], + [0, 0, -2, 2, 1], + [2, 0, 0, 2, 2], + [0, 0, 2, 2, 2], + [0, 0, 2, 0, 0], + [-2, 0, 1, 2, 2], + [0, 0, 0, 2, 0], + [-2, 0, 0, 2, 0], + [0, 0, -1, 2, 1], + [0, 2, 0, 0, 0], + [2, 0, -1, 0, 1], + [-2, 2, 0, 2, 2], + [0, 1, 0, 0, 1], + [-2, 0, 1, 0, 1], + [0, -1, 0, 0, 1], + [0, 0, 2, -2, 0], + [2, 0, -1, 2, 1], + [2, 0, 1, 2, 2], + [0, 1, 0, 2, 2], + [-2, 1, 1, 0, 0], + [0, -1, 0, 2, 2], + [2, 0, 0, 2, 1], + [2, 0, 1, 0, 0], + [-2, 0, 2, 2, 2], + [-2, 0, 1, 2, 1], + [2, 0, -2, 0, 1], + [2, 0, 0, 0, 1], + [0, -1, 1, 0, 0], + [-2, -1, 0, 2, 1], + [-2, 0, 0, 0, 1], + [0, 0, 2, 2, 1], + [-2, 0, 2, 0, 1], + [-2, 1, 0, 2, 1], + [0, 0, 1, -2, 0], + [-1, 0, 1, 0, 0], + [-2, 1, 0, 0, 0], + [1, 0, 0, 0, 0], + [0, 0, 1, 2, 0], + [0, 0, -2, 2, 2], + [-1, -1, 1, 0, 0], + [0, 1, 1, 0, 0], + [0, -1, 1, 2, 2], + [2, -1, -1, 2, 2], + [0, 0, 3, 2, 2], + [2, -1, 0, 2, 2], +]) + + +@jcompile('float64(int64, int64, int64, int64, int64, int64, int64)', + nopython=True) +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 + + +@jcompile('float64(float64)', nopython=True) +def julian_day(unixtime): + jd = unixtime * 1.0 / 86400 + 2440587.5 + return jd + + +@jcompile('float64(float64, float64)', nopython=True) +def julian_ephemeris_day(julian_day, delta_t): + jde = julian_day + delta_t * 1.0 / 86400 + return jde + + +@jcompile('float64(float64)', nopython=True) +def julian_century(julian_day): + jc = (julian_day - 2451545) * 1.0 / 36525 + return jc + + +@jcompile('float64(float64)', nopython=True) +def julian_ephemeris_century(julian_ephemeris_day): + jce = (julian_ephemeris_day - 2451545) * 1.0 / 36525 + return jce + + +@jcompile('float64(float64)', nopython=True) +def julian_ephemeris_millennium(julian_ephemeris_century): + jme = julian_ephemeris_century * 1.0 / 10 + return jme + + +# omit type signature here; specifying read-only arrays requires use of the +# numba.types API, meaning numba must be available to import. +# https://github.com/numba/numba/issues/4511 +@jcompile(nopython=True) +def sum_mult_cos_add_mult(arr, x): + # shared calculation used for heliocentric longitude, latitude, and radius + s = 0. + for row in range(arr.shape[0]): + s += arr[row, 0] * np.cos(arr[row, 1] + arr[row, 2] * x) + return s + +@jcompile('float64(float64)', nopython=True) +def heliocentric_longitude(jme): + l0 = sum_mult_cos_add_mult(L0, jme) + l1 = sum_mult_cos_add_mult(L1, jme) + l2 = sum_mult_cos_add_mult(L2, jme) + l3 = sum_mult_cos_add_mult(L3, jme) + l4 = sum_mult_cos_add_mult(L4, jme) + l5 = sum_mult_cos_add_mult(L5, jme) + + l_rad = (l0 + l1 * jme + l2 * jme**2 + l3 * jme**3 + l4 * jme**4 + + l5 * jme**5)/10**8 + l = np.rad2deg(l_rad) + return l % 360 + +@jcompile('float64(float64)', nopython=True) +def heliocentric_latitude(jme): + b0 = sum_mult_cos_add_mult(B0, jme) + b1 = sum_mult_cos_add_mult(B1, jme) + + b_rad = (b0 + b1 * jme)/10**8 + b = np.rad2deg(b_rad) + return b + + +@jcompile('float64(float64)', nopython=True) +def heliocentric_radius_vector(jme): + r0 = sum_mult_cos_add_mult(R0, jme) + r1 = sum_mult_cos_add_mult(R1, jme) + r2 = sum_mult_cos_add_mult(R2, jme) + r3 = sum_mult_cos_add_mult(R3, jme) + r4 = sum_mult_cos_add_mult(R4, jme) + + r = (r0 + r1 * jme + r2 * jme**2 + r3 * jme**3 + r4 * jme**4)/10**8 + return r + + +@jcompile('float64(float64)', nopython=True) +def geocentric_longitude(heliocentric_longitude): + theta = heliocentric_longitude + 180.0 + return theta % 360 + + +@jcompile('float64(float64)', nopython=True) +def geocentric_latitude(heliocentric_latitude): + beta = -1.0*heliocentric_latitude + return beta + + +@jcompile('float64(float64)', nopython=True) +def mean_elongation(julian_ephemeris_century): + x0 = (297.85036 + + 445267.111480 * julian_ephemeris_century + - 0.0019142 * julian_ephemeris_century**2 + + julian_ephemeris_century**3 / 189474) + return x0 + + +@jcompile('float64(float64)', nopython=True) +def mean_anomaly_sun(julian_ephemeris_century): + x1 = (357.52772 + + 35999.050340 * julian_ephemeris_century + - 0.0001603 * julian_ephemeris_century**2 + - julian_ephemeris_century**3 / 300000) + return x1 + + +@jcompile('float64(float64)', nopython=True) +def mean_anomaly_moon(julian_ephemeris_century): + x2 = (134.96298 + + 477198.867398 * julian_ephemeris_century + + 0.0086972 * julian_ephemeris_century**2 + + julian_ephemeris_century**3 / 56250) + return x2 + + +@jcompile('float64(float64)', nopython=True) +def moon_argument_latitude(julian_ephemeris_century): + x3 = (93.27191 + + 483202.017538 * julian_ephemeris_century + - 0.0036825 * julian_ephemeris_century**2 + + julian_ephemeris_century**3 / 327270) + return x3 + + +@jcompile('float64(float64)', nopython=True) +def moon_ascending_longitude(julian_ephemeris_century): + x4 = (125.04452 + - 1934.136261 * julian_ephemeris_century + + 0.0020708 * julian_ephemeris_century**2 + + julian_ephemeris_century**3 / 450000) + return x4 + + +@jcompile( + 'void(float64, float64, float64, float64, float64, float64, float64[:])', + nopython=True) +def longitude_obliquity_nutation(julian_ephemeris_century, x0, x1, x2, x3, x4, + out): + delta_psi_sum = 0.0 + delta_eps_sum = 0.0 + for row in range(NUTATION_YTERM_ARRAY.shape[0]): + a = NUTATION_ABCD_ARRAY[row, 0] + b = NUTATION_ABCD_ARRAY[row, 1] + c = NUTATION_ABCD_ARRAY[row, 2] + d = NUTATION_ABCD_ARRAY[row, 3] + arg = np.radians( + NUTATION_YTERM_ARRAY[row, 0]*x0 + + NUTATION_YTERM_ARRAY[row, 1]*x1 + + NUTATION_YTERM_ARRAY[row, 2]*x2 + + NUTATION_YTERM_ARRAY[row, 3]*x3 + + NUTATION_YTERM_ARRAY[row, 4]*x4 + ) + delta_psi_sum += (a + b * julian_ephemeris_century) * np.sin(arg) + delta_eps_sum += (c + d * julian_ephemeris_century) * np.cos(arg) + delta_psi = delta_psi_sum*1.0/36000000 + delta_eps = delta_eps_sum*1.0/36000000 + # seems like we ought to be able to return a tuple here instead + # of resorting to `out`, but returning a UniTuple from this + # function caused calculations elsewhere to give the wrong result. + # very difficult to investigate since it did not occur when using + # object mode. issue was observed on numba 0.56.4 + out[0] = delta_psi + out[1] = delta_eps + + +@jcompile('float64(float64)', nopython=True) +def mean_ecliptic_obliquity(julian_ephemeris_millennium): + U = 1.0*julian_ephemeris_millennium/10 + e0 = (84381.448 - 4680.93 * U - 1.55 * U**2 + + 1999.25 * U**3 - 51.38 * U**4 - 249.67 * U**5 + - 39.05 * U**6 + 7.12 * U**7 + 27.87 * U**8 + + 5.79 * U**9 + 2.45 * U**10) + return e0 + + +@jcompile('float64(float64, float64)', nopython=True) +def true_ecliptic_obliquity(mean_ecliptic_obliquity, obliquity_nutation): + e0 = mean_ecliptic_obliquity + deleps = obliquity_nutation + e = e0*1.0/3600 + deleps + return e + + +@jcompile('float64(float64)', nopython=True) +def aberration_correction(earth_radius_vector): + deltau = -20.4898 / (3600 * earth_radius_vector) + return deltau + + +@jcompile('float64(float64, float64, float64)', nopython=True) +def apparent_sun_longitude(geocentric_longitude, longitude_nutation, + aberration_correction): + lamd = geocentric_longitude + longitude_nutation + aberration_correction + return lamd + + +@jcompile('float64(float64, float64)', nopython=True) +def mean_sidereal_time(julian_day, julian_century): + v0 = (280.46061837 + 360.98564736629 * (julian_day - 2451545) + + 0.000387933 * julian_century**2 - julian_century**3 / 38710000) + return v0 % 360.0 + + +@jcompile('float64(float64, float64, float64)', nopython=True) +def apparent_sidereal_time(mean_sidereal_time, longitude_nutation, + true_ecliptic_obliquity): + v = mean_sidereal_time + longitude_nutation * np.cos( + np.radians(true_ecliptic_obliquity)) + return v + + +@jcompile('float64(float64, float64, float64)', nopython=True) +def geocentric_sun_right_ascension(apparent_sun_longitude, + true_ecliptic_obliquity, + geocentric_latitude): + true_ecliptic_obliquity_rad = np.radians(true_ecliptic_obliquity) + apparent_sun_longitude_rad = np.radians(apparent_sun_longitude) + + num = (np.sin(apparent_sun_longitude_rad) + * np.cos(true_ecliptic_obliquity_rad) + - np.tan(np.radians(geocentric_latitude)) + * np.sin(true_ecliptic_obliquity_rad)) + alpha = np.degrees(np.arctan2(num, np.cos(apparent_sun_longitude_rad))) + return alpha % 360 + + +@jcompile('float64(float64, float64, float64)', nopython=True) +def geocentric_sun_declination(apparent_sun_longitude, true_ecliptic_obliquity, + geocentric_latitude): + geocentric_latitude_rad = np.radians(geocentric_latitude) + true_ecliptic_obliquity_rad = np.radians(true_ecliptic_obliquity) + + delta = np.degrees(np.arcsin(np.sin(geocentric_latitude_rad) * + np.cos(true_ecliptic_obliquity_rad) + + np.cos(geocentric_latitude_rad) * + np.sin(true_ecliptic_obliquity_rad) * + np.sin(np.radians(apparent_sun_longitude)))) + return delta + + +@jcompile('float64(float64, float64, float64)', nopython=True) +def local_hour_angle(apparent_sidereal_time, observer_longitude, + sun_right_ascension): + """Measured westward from south""" + H = apparent_sidereal_time + observer_longitude - sun_right_ascension + return H % 360 + + +@jcompile('float64(float64)', nopython=True) +def equatorial_horizontal_parallax(earth_radius_vector): + xi = 8.794 / (3600 * earth_radius_vector) + return xi + + +@jcompile('float64(float64)', nopython=True) +def uterm(observer_latitude): + u = np.arctan(0.99664719 * np.tan(np.radians(observer_latitude))) + return u + + +@jcompile('float64(float64, float64, float64)', nopython=True) +def xterm(u, observer_latitude, observer_elevation): + x = (np.cos(u) + observer_elevation / 6378140 + * np.cos(np.radians(observer_latitude))) + return x + + +@jcompile('float64(float64, float64, float64)', nopython=True) +def yterm(u, observer_latitude, observer_elevation): + y = (0.99664719 * np.sin(u) + observer_elevation / 6378140 + * np.sin(np.radians(observer_latitude))) + return y + + +@jcompile('float64(float64, float64,float64, float64)', nopython=True) +def parallax_sun_right_ascension(xterm, equatorial_horizontal_parallax, + local_hour_angle, geocentric_sun_declination): + equatorial_horizontal_parallax_rad = \ + np.radians(equatorial_horizontal_parallax) + local_hour_angle_rad = np.radians(local_hour_angle) + + num = (-xterm * np.sin(equatorial_horizontal_parallax_rad) + * np.sin(local_hour_angle_rad)) + denom = (np.cos(np.radians(geocentric_sun_declination)) + - xterm * np.sin(equatorial_horizontal_parallax_rad) + * np.cos(local_hour_angle_rad)) + delta_alpha = np.degrees(np.arctan2(num, denom)) + return delta_alpha + + +@jcompile('float64(float64, float64)', nopython=True) +def topocentric_sun_right_ascension(geocentric_sun_right_ascension, + parallax_sun_right_ascension): + alpha_prime = geocentric_sun_right_ascension + parallax_sun_right_ascension + return alpha_prime + + +@jcompile('float64(float64, float64, float64, float64, float64, float64)', + nopython=True) +def topocentric_sun_declination(geocentric_sun_declination, xterm, yterm, + equatorial_horizontal_parallax, + parallax_sun_right_ascension, + local_hour_angle): + geocentric_sun_declination_rad = np.radians(geocentric_sun_declination) + equatorial_horizontal_parallax_rad = \ + np.radians(equatorial_horizontal_parallax) + + num = ((np.sin(geocentric_sun_declination_rad) - yterm + * np.sin(equatorial_horizontal_parallax_rad)) + * np.cos(np.radians(parallax_sun_right_ascension))) + denom = (np.cos(geocentric_sun_declination_rad) - xterm + * np.sin(equatorial_horizontal_parallax_rad) + * np.cos(np.radians(local_hour_angle))) + delta = np.degrees(np.arctan2(num, denom)) + return delta + + +@jcompile('float64(float64, float64)', nopython=True) +def topocentric_local_hour_angle(local_hour_angle, + parallax_sun_right_ascension): + H_prime = local_hour_angle - parallax_sun_right_ascension + return H_prime + + +@jcompile('float64(float64, float64, float64)', nopython=True) +def topocentric_elevation_angle_without_atmosphere(observer_latitude, + topocentric_sun_declination, + topocentric_local_hour_angle + ): + observer_latitude_rad = np.radians(observer_latitude) + topocentric_sun_declination_rad = np.radians(topocentric_sun_declination) + e0 = np.degrees(np.arcsin( + np.sin(observer_latitude_rad) + * np.sin(topocentric_sun_declination_rad) + + np.cos(observer_latitude_rad) + * np.cos(topocentric_sun_declination_rad) + * np.cos(np.radians(topocentric_local_hour_angle)))) + return e0 + + +@jcompile('float64(float64, float64, float64, float64)', nopython=True) +def atmospheric_refraction_correction(local_pressure, local_temp, + topocentric_elevation_angle_wo_atmosphere, + atmos_refract): + # switch sets delta_e when the sun is below the horizon + switch = topocentric_elevation_angle_wo_atmosphere >= -1.0 * ( + 0.26667 + atmos_refract) + delta_e = ((local_pressure / 1010.0) * (283.0 / (273 + local_temp)) + * 1.02 / (60 * np.tan(np.radians( + topocentric_elevation_angle_wo_atmosphere + + 10.3 / (topocentric_elevation_angle_wo_atmosphere + + 5.11))))) * switch + return delta_e + + +@jcompile('float64(float64, float64)', nopython=True) +def topocentric_elevation_angle(topocentric_elevation_angle_without_atmosphere, + atmospheric_refraction_correction): + e = (topocentric_elevation_angle_without_atmosphere + + atmospheric_refraction_correction) + return e + + +@jcompile('float64(float64)', nopython=True) +def topocentric_zenith_angle(topocentric_elevation_angle): + theta = 90 - topocentric_elevation_angle + return theta + + +@jcompile('float64(float64, float64, float64)', nopython=True) +def topocentric_astronomers_azimuth(topocentric_local_hour_angle, + topocentric_sun_declination, + observer_latitude): + topocentric_local_hour_angle_rad = np.radians(topocentric_local_hour_angle) + observer_latitude_rad = np.radians(observer_latitude) + num = np.sin(topocentric_local_hour_angle_rad) + denom = (np.cos(topocentric_local_hour_angle_rad) + * np.sin(observer_latitude_rad) + - np.tan(np.radians(topocentric_sun_declination)) + * np.cos(observer_latitude_rad)) + gamma = np.degrees(np.arctan2(num, denom)) + return gamma % 360 + + +@jcompile('float64(float64)', nopython=True) +def topocentric_azimuth_angle(topocentric_astronomers_azimuth): + phi = topocentric_astronomers_azimuth + 180 + return phi % 360 + + +@jcompile('float64(float64)', nopython=True) +def sun_mean_longitude(julian_ephemeris_millennium): + M = (280.4664567 + 360007.6982779 * julian_ephemeris_millennium + + 0.03032028 * julian_ephemeris_millennium**2 + + julian_ephemeris_millennium**3 / 49931 + - julian_ephemeris_millennium**4 / 15300 + - julian_ephemeris_millennium**5 / 2000000) + return M + + +@jcompile('float64(float64, float64, float64, float64)', nopython=True) +def equation_of_time(sun_mean_longitude, geocentric_sun_right_ascension, + longitude_nutation, true_ecliptic_obliquity): + E = (sun_mean_longitude - 0.0057183 - geocentric_sun_right_ascension + + longitude_nutation * np.cos(np.radians(true_ecliptic_obliquity))) + # limit between 0 and 360 + E = E % 360 + # convert to minutes + E *= 4 + greater = E > 20 + less = E < -20 + other = (E <= 20) & (E >= -20) + E = greater * (E - 1440) + less * (E + 1440) + other * E + return E + + +@jcompile('void(float64[:], float64[:], float64[:], float64[:,:])', + nopython=True, nogil=True) +def solar_position_loop(unixtime, delta_t, loc_args, out): + """Loop through the time array and calculate the solar position""" + lat = loc_args[0] + lon = loc_args[1] + elev = loc_args[2] + pressure = loc_args[3] + temp = loc_args[4] + atmos_refract = loc_args[5] + sst = loc_args[6] + esd = loc_args[7] + + for i in range(unixtime.shape[0]): + utime = unixtime[i] + dT = delta_t[i] + jd = julian_day(utime) + jde = julian_ephemeris_day(jd, dT) + jc = julian_century(jd) + jce = julian_ephemeris_century(jde) + jme = julian_ephemeris_millennium(jce) + R = heliocentric_radius_vector(jme) + if esd: + out[0, i] = R + continue + L = heliocentric_longitude(jme) + B = heliocentric_latitude(jme) + Theta = geocentric_longitude(L) + beta = geocentric_latitude(B) + x0 = mean_elongation(jce) + x1 = mean_anomaly_sun(jce) + x2 = mean_anomaly_moon(jce) + x3 = moon_argument_latitude(jce) + x4 = moon_ascending_longitude(jce) + l_o_nutation = np.empty((2,)) + longitude_obliquity_nutation(jce, x0, x1, x2, x3, x4, l_o_nutation) + delta_psi = l_o_nutation[0] + delta_epsilon = l_o_nutation[1] + epsilon0 = mean_ecliptic_obliquity(jme) + epsilon = true_ecliptic_obliquity(epsilon0, delta_epsilon) + delta_tau = aberration_correction(R) + lamd = apparent_sun_longitude(Theta, delta_psi, delta_tau) + v0 = mean_sidereal_time(jd, jc) + v = apparent_sidereal_time(v0, delta_psi, epsilon) + alpha = geocentric_sun_right_ascension(lamd, epsilon, beta) + delta = geocentric_sun_declination(lamd, epsilon, beta) + if sst: + out[0, i] = v + out[1, i] = alpha + out[2, i] = delta + continue + m = sun_mean_longitude(jme) + eot = equation_of_time(m, alpha, delta_psi, epsilon) + H = local_hour_angle(v, lon, alpha) + xi = equatorial_horizontal_parallax(R) + u = uterm(lat) + x = xterm(u, lat, elev) + y = yterm(u, lat, elev) + delta_alpha = parallax_sun_right_ascension(x, xi, H, delta) + delta_prime = topocentric_sun_declination(delta, x, y, xi, delta_alpha, + H) + H_prime = topocentric_local_hour_angle(H, delta_alpha) + e0 = topocentric_elevation_angle_without_atmosphere(lat, delta_prime, + H_prime) + delta_e = atmospheric_refraction_correction(pressure, temp, e0, + atmos_refract) + e = topocentric_elevation_angle(e0, delta_e) + theta = topocentric_zenith_angle(e) + theta0 = topocentric_zenith_angle(e0) + gamma = topocentric_astronomers_azimuth(H_prime, delta_prime, lat) + phi = topocentric_azimuth_angle(gamma) + out[0, i] = theta + out[1, i] = theta0 + out[2, i] = e + out[3, i] = e0 + out[4, i] = phi + out[5, i] = eot + + +def solar_position_numba(unixtime, lat, lon, elev, pressure, temp, delta_t, + atmos_refract, numthreads, sst=False, esd=False): + """Calculate the solar position using the numba compiled functions + and multiple threads. Very slow if functions are not numba compiled. + """ + # these args are the same for each thread + loc_args = np.array([lat, lon, elev, pressure, temp, + atmos_refract, sst, esd], dtype=np.float64) + + # turn delta_t into an array if it isn't already + delta_t = np.full_like(unixtime, delta_t, dtype=np.float64) + + # construct dims x ulength array to put the results in + ulength = unixtime.shape[0] + if sst: + dims = 3 + elif esd: + dims = 1 + else: + dims = 6 + result = np.empty((dims, ulength), dtype=np.float64) + + if unixtime.dtype != np.float64: + unixtime = unixtime.astype(np.float64) + + if ulength < numthreads: + numthreads = ulength + + if numthreads <= 1: + solar_position_loop(unixtime, delta_t, loc_args, result) + return result + + # split the input and output arrays into numthreads chunks + split0 = np.array_split(unixtime, numthreads) + split1 = np.array_split(delta_t, numthreads) + split2 = np.array_split(result, numthreads, axis=1) + chunks = [ + [a0, a1, loc_args, a2] + for a0, a1, a2 in zip(split0, split1, split2) + ] + # Spawn one thread per chunk + threads = [threading.Thread(target=solar_position_loop, args=chunk) + for chunk in chunks] + for thread in threads: + thread.start() + for thread in threads: + thread.join() + return result + + +def solar_position_numpy(unixtime, lat, lon, elev, pressure, temp, delta_t, + atmos_refract, numthreads, sst=False, esd=False): + """Calculate the solar position assuming unixtime is a numpy array. Note + this function will not work if the solar position functions were + compiled with numba. + """ + + jd = julian_day(unixtime) + jde = julian_ephemeris_day(jd, delta_t) + jc = julian_century(jd) + jce = julian_ephemeris_century(jde) + jme = julian_ephemeris_millennium(jce) + R = heliocentric_radius_vector(jme) + if esd: + return (R, ) + L = heliocentric_longitude(jme) + B = heliocentric_latitude(jme) + Theta = geocentric_longitude(L) + beta = geocentric_latitude(B) + x0 = mean_elongation(jce) + x1 = mean_anomaly_sun(jce) + x2 = mean_anomaly_moon(jce) + x3 = moon_argument_latitude(jce) + x4 = moon_ascending_longitude(jce) + l_o_nutation = np.empty((2, len(x0))) + longitude_obliquity_nutation(jce, x0, x1, x2, x3, x4, l_o_nutation) + delta_psi = l_o_nutation[0] + delta_epsilon = l_o_nutation[1] + epsilon0 = mean_ecliptic_obliquity(jme) + epsilon = true_ecliptic_obliquity(epsilon0, delta_epsilon) + delta_tau = aberration_correction(R) + lamd = apparent_sun_longitude(Theta, delta_psi, delta_tau) + v0 = mean_sidereal_time(jd, jc) + v = apparent_sidereal_time(v0, delta_psi, epsilon) + alpha = geocentric_sun_right_ascension(lamd, epsilon, beta) + delta = geocentric_sun_declination(lamd, epsilon, beta) + if sst: + return v, alpha, delta + m = sun_mean_longitude(jme) + eot = equation_of_time(m, alpha, delta_psi, epsilon) + H = local_hour_angle(v, lon, alpha) + xi = equatorial_horizontal_parallax(R) + u = uterm(lat) + x = xterm(u, lat, elev) + y = yterm(u, lat, elev) + delta_alpha = parallax_sun_right_ascension(x, xi, H, delta) + delta_prime = topocentric_sun_declination(delta, x, y, xi, delta_alpha, H) + H_prime = topocentric_local_hour_angle(H, delta_alpha) + e0 = topocentric_elevation_angle_without_atmosphere(lat, delta_prime, + H_prime) + delta_e = atmospheric_refraction_correction(pressure, temp, e0, + atmos_refract) + e = topocentric_elevation_angle(e0, delta_e) + theta = topocentric_zenith_angle(e) + theta0 = topocentric_zenith_angle(e0) + gamma = topocentric_astronomers_azimuth(H_prime, delta_prime, lat) + phi = topocentric_azimuth_angle(gamma) + return theta, theta0, e, e0, phi, eot + + +def solar_position(unixtime, lat, lon, elev, pressure, temp, delta_t, + atmos_refract, numthreads=8, sst=False, esd=False): + + """ + Calculate the solar position using the + NREL SPA described in [1]_, [2]_. + + If numba is installed, the functions can be compiled + and the code runs quickly. If not, the functions + still evaluate but use numpy instead. + + Parameters + ---------- + unixtime : numpy array + Array of unix/epoch timestamps to calculate solar position for. + Unixtime is the number of seconds since Jan. 1, 1970 00:00:00 UTC. + A pandas.DatetimeIndex is easily converted using .view(np.int64)/10**9 + lat : float + Latitude to calculate solar position for + lon : float + Longitude to calculate solar position for + elev : float + Elevation of location in meters + pressure : int or float + avg. yearly pressure at location in millibars; + used for atmospheric correction + temp : int or float + avg. yearly temperature at location in + degrees C; used for atmospheric correction + delta_t : float or array + Difference between terrestrial time and UT1. + atmos_refract : float + The approximate atmospheric refraction (in degrees) + at sunrise and sunset. + numthreads: int, optional, default 8 + Number of threads to use for computation if numba>=0.17 + is installed. + sst : bool, default False + If True, return only data needed for sunrise, sunset, and transit + calculations. + esd : bool, default False + If True, return only Earth-Sun distance in AU + + Returns + ------- + Numpy Array with elements: + apparent zenith, + zenith, + apparent_elevation, + elevation, + azimuth, + equation_of_time + + References + ---------- + [1] I. Reda and A. Andreas, Solar position algorithm for solar radiation + applications. Solar Energy, vol. 76, no. 5, pp. 577-589, 2004. + + [2] I. Reda and A. Andreas, Corrigendum to Solar position algorithm for + solar radiation applications. Solar Energy, vol. 81, no. 6, p. 838, 2007. + """ + if USE_NUMBA: + do_calc = solar_position_numba + else: + do_calc = solar_position_numpy + + result = do_calc(unixtime, lat, lon, elev, pressure, + temp, delta_t, atmos_refract, numthreads, + sst, esd) + + if not isinstance(result, np.ndarray): + try: + result = np.array(result) + except Exception: + pass + + return result + + +def transit_sunrise_sunset(dates, lat, lon, delta_t, numthreads): + """ + Calculate the sun transit, sunrise, and sunset + for a set of dates at a given location using the NREL SPA. + + The NREL SPA is documented in [1]_, [2]_. + + Parameters + ---------- + dates : array + Numpy array of ints/floats corresponding to the Unix time + for the dates of interest, must be midnight UTC (00:00+00:00) + on the day of interest. + lat : float + Latitude of location to perform calculation for + lon : float + Longitude of location + delta_t : float or array + Difference between terrestrial time and UT. USNO has tables. + numthreads : int + Number to threads to use for calculation (if using numba) + + Returns + ------- + tuple : (transit, sunrise, sunset) localized to UTC + References ---------- .. [1] I. Reda and A. Andreas, Solar position algorithm for solar @@ -1145,244 +1142,244 @@ def transit_sunrise_sunset(dates, lat, lon, delta_t, numthreads): (2004) 577–589] Solar Energy, vol. 81, no. 6, pp. 838, 2007. :doi:`10.1016/j.solener.2007.01.003` """ - - if ((dates % 86400) != 0.0).any(): - raise ValueError('Input dates must be at 00:00 UTC') - - utday = (dates // 86400) * 86400 - ttday0 = utday - delta_t - ttdayn1 = ttday0 - 86400 - ttdayp1 = ttday0 + 86400 - - # index 0 is v, 1 is alpha, 2 is delta - utday_res = solar_position(utday, 0, 0, 0, 0, 0, delta_t, - 0, numthreads, sst=True) - v = utday_res[0] - - ttday0_res = solar_position(ttday0, 0, 0, 0, 0, 0, delta_t, - 0, numthreads, sst=True) - ttdayn1_res = solar_position(ttdayn1, 0, 0, 0, 0, 0, delta_t, - 0, numthreads, sst=True) - ttdayp1_res = solar_position(ttdayp1, 0, 0, 0, 0, 0, delta_t, - 0, numthreads, sst=True) - m0 = (ttday0_res[1] - lon - v) / 360 - cos_arg = ((np.sin(np.radians(-0.8333)) - np.sin(np.radians(lat)) - * np.sin(np.radians(ttday0_res[2]))) / - (np.cos(np.radians(lat)) * np.cos(np.radians(ttday0_res[2])))) - cos_arg[abs(cos_arg) > 1] = np.nan - H0 = np.degrees(np.arccos(cos_arg)) % 180 - - m = np.empty((3, len(utday))) - m[0] = m0 % 1 - m[1] = (m[0] - H0 / 360) - m[2] = (m[0] + H0 / 360) - - # need to account for fractions of day that may be the next or previous - # day in UTC - add_a_day = m[2] >= 1 - sub_a_day = m[1] < 0 - m[1] = m[1] % 1 - m[2] = m[2] % 1 - vs = v + 360.985647 * m - n = m + delta_t / 86400 - - a = ttday0_res[1] - ttdayn1_res[1] - a[abs(a) > 2] = a[abs(a) > 2] % 1 - ap = ttday0_res[2] - ttdayn1_res[2] - ap[abs(ap) > 2] = ap[abs(ap) > 2] % 1 - b = ttdayp1_res[1] - ttday0_res[1] - b[abs(b) > 2] = b[abs(b) > 2] % 1 - bp = ttdayp1_res[2] - ttday0_res[2] - bp[abs(bp) > 2] = bp[abs(bp) > 2] % 1 - c = b - a - cp = bp - ap - - alpha_prime = ttday0_res[1] + (n * (a + b + c * n)) / 2 - delta_prime = ttday0_res[2] + (n * (ap + bp + cp * n)) / 2 - Hp = (vs + lon - alpha_prime) % 360 - Hp[Hp >= 180] = Hp[Hp >= 180] - 360 - - h = np.degrees(np.arcsin(np.sin(np.radians(lat)) * - np.sin(np.radians(delta_prime)) + - np.cos(np.radians(lat)) * - np.cos(np.radians(delta_prime)) - * np.cos(np.radians(Hp)))) - - T = (m[0] - Hp[0] / 360) * 86400 - R = (m[1] + (h[1] + 0.8333) / (360 * np.cos(np.radians(delta_prime[1])) * - np.cos(np.radians(lat)) * - np.sin(np.radians(Hp[1])))) * 86400 - S = (m[2] + (h[2] + 0.8333) / (360 * np.cos(np.radians(delta_prime[2])) * - np.cos(np.radians(lat)) * - np.sin(np.radians(Hp[2])))) * 86400 - - S[add_a_day] += 86400 - R[sub_a_day] -= 86400 - - transit = T + utday - sunrise = R + utday - sunset = S + utday - - return transit, sunrise, sunset - - -def earthsun_distance(unixtime, delta_t, numthreads): - """ - Calculates the distance from the earth to the sun using the - NREL SPA algorithm described in [1]. - - Parameters - ---------- - unixtime : numpy array - Array of unix/epoch timestamps to calculate solar position for. - Unixtime is the number of seconds since Jan. 1, 1970 00:00:00 UTC. - A pandas.DatetimeIndex is easily converted using .view(np.int64)/10**9 - delta_t : float or array - Difference between terrestrial time and UT. USNO has tables. - numthreads : int - Number to threads to use for calculation (if using numba) - - Returns - ------- - R : array - Earth-Sun distance in AU. - - References - ---------- - [1] Reda, I., Andreas, A., 2003. Solar position algorithm for solar - radiation applications. Technical report: NREL/TP-560-34302. Golden, - USA, http://www.nlr.gov. - """ - - R = solar_position(unixtime, 0, 0, 0, 0, 0, delta_t, - 0, numthreads, esd=True)[0] - - return R - - -def calculate_deltat(year, month): - """Calculate the difference between Terrestrial Dynamical Time (TD) - and Universal Time (UT). - - Equations taken from http://eclipse.gsfc.nasa.gov/SEcat5/deltatpoly.html - """ - - plw = 'Deltat is unknown for years before -1999 and after 3000. ' \ - 'Delta values will be calculated, but the calculations ' \ - 'are not intended to be used for these years.' - - try: - if np.any((year > 3000) | (year < -1999)): - warnings.warn(plw) - except ValueError: - if (year > 3000) | (year < -1999): - warnings.warn(plw) - except TypeError: - return 0 - - y = year + (month - 0.5)/12 - - deltat = np.where(year < -500, - - -20+32*((y-1820)/100)**2, 0) - - deltat = np.where((-500 <= year) & (year < 500), - - 10583.6-1014.41*(y/100) - + 33.78311*(y/100)**2 - - 5.952053*(y/100)**3 - - 0.1798452*(y/100)**4 - + 0.022174192*(y/100)**5 - + 0.0090316521*(y/100)**6, deltat) - - deltat = np.where((500 <= year) & (year < 1600), - - 1574.2-556.01*((y-1000)/100) - + 71.23472*((y-1000)/100)**2 - + 0.319781*((y-1000)/100)**3 - - 0.8503463*((y-1000)/100)**4 - - 0.005050998*((y-1000)/100)**5 - + 0.0083572073*((y-1000)/100)**6, deltat) - - deltat = np.where((1600 <= year) & (year < 1700), - - 120-0.9808*(y-1600) - - 0.01532*(y-1600)**2 - + (y-1600)**3/7129, deltat) - - deltat = np.where((1700 <= year) & (year < 1800), - - 8.83+0.1603*(y-1700) - - 0.0059285*(y-1700)**2 - + 0.00013336*(y-1700)**3 - - (y-1700)**4/1174000, deltat) - - deltat = np.where((1800 <= year) & (year < 1860), - - 13.72-0.332447*(y-1800) - + 0.0068612*(y-1800)**2 - + 0.0041116*(y-1800)**3 - - 0.00037436*(y-1800)**4 - + 0.0000121272*(y-1800)**5 - - 0.0000001699*(y-1800)**6 - + 0.000000000875*(y-1800)**7, deltat) - - deltat = np.where((1860 <= year) & (year < 1900), - - 7.62+0.5737*(y-1860) - - 0.251754*(y-1860)**2 - + 0.01680668*(y-1860)**3 - - 0.0004473624*(y-1860)**4 - + (y-1860)**5/233174, deltat) - - deltat = np.where((1900 <= year) & (year < 1920), - - -2.79+1.494119*(y-1900) - - 0.0598939*(y-1900)**2 - + 0.0061966*(y-1900)**3 - - 0.000197*(y-1900)**4, deltat) - - deltat = np.where((1920 <= year) & (year < 1941), - - 21.20+0.84493*(y-1920) - - 0.076100*(y-1920)**2 - + 0.0020936*(y-1920)**3, deltat) - - deltat = np.where((1941 <= year) & (year < 1961), - - 29.07+0.407*(y-1950) - - (y-1950)**2/233 - + (y-1950)**3/2547, deltat) - - deltat = np.where((1961 <= year) & (year < 1986), - - 45.45+1.067*(y-1975) - - (y-1975)**2/260 - - (y-1975)**3/718, deltat) - - deltat = np.where((1986 <= year) & (year < 2005), - - 63.86+0.3345*(y-2000) - - 0.060374*(y-2000)**2 - + 0.0017275*(y-2000)**3 - + 0.000651814*(y-2000)**4 - + 0.00002373599*(y-2000)**5, deltat) - - deltat = np.where((2005 <= year) & (year < 2050), - - 62.92+0.32217*(y-2000) - + 0.005589*(y-2000)**2, deltat) - - deltat = np.where((2050 <= year) & (year < 2150), - - -20+32*((y-1820)/100)**2 - - 0.5628*(2150-y), deltat) - - deltat = np.where(year >= 2150, - - -20+32*((y-1820)/100)**2, deltat) - - deltat = deltat.item() if np.isscalar(year) & np.isscalar(month)\ - else deltat - - return deltat + + if ((dates % 86400) != 0.0).any(): + raise ValueError('Input dates must be at 00:00 UTC') + + utday = (dates // 86400) * 86400 + ttday0 = utday - delta_t + ttdayn1 = ttday0 - 86400 + ttdayp1 = ttday0 + 86400 + + # index 0 is v, 1 is alpha, 2 is delta + utday_res = solar_position(utday, 0, 0, 0, 0, 0, delta_t, + 0, numthreads, sst=True) + v = utday_res[0] + + ttday0_res = solar_position(ttday0, 0, 0, 0, 0, 0, delta_t, + 0, numthreads, sst=True) + ttdayn1_res = solar_position(ttdayn1, 0, 0, 0, 0, 0, delta_t, + 0, numthreads, sst=True) + ttdayp1_res = solar_position(ttdayp1, 0, 0, 0, 0, 0, delta_t, + 0, numthreads, sst=True) + m0 = (ttday0_res[1] - lon - v) / 360 + cos_arg = ((np.sin(np.radians(-0.8333)) - np.sin(np.radians(lat)) + * np.sin(np.radians(ttday0_res[2]))) / + (np.cos(np.radians(lat)) * np.cos(np.radians(ttday0_res[2])))) + cos_arg[abs(cos_arg) > 1] = np.nan + H0 = np.degrees(np.arccos(cos_arg)) % 180 + + m = np.empty((3, len(utday))) + m[0] = m0 % 1 + m[1] = (m[0] - H0 / 360) + m[2] = (m[0] + H0 / 360) + + # need to account for fractions of day that may be the next or previous + # day in UTC + add_a_day = m[2] >= 1 + sub_a_day = m[1] < 0 + m[1] = m[1] % 1 + m[2] = m[2] % 1 + vs = v + 360.985647 * m + n = m + delta_t / 86400 + + a = ttday0_res[1] - ttdayn1_res[1] + a[abs(a) > 2] = a[abs(a) > 2] % 1 + ap = ttday0_res[2] - ttdayn1_res[2] + ap[abs(ap) > 2] = ap[abs(ap) > 2] % 1 + b = ttdayp1_res[1] - ttday0_res[1] + b[abs(b) > 2] = b[abs(b) > 2] % 1 + bp = ttdayp1_res[2] - ttday0_res[2] + bp[abs(bp) > 2] = bp[abs(bp) > 2] % 1 + c = b - a + cp = bp - ap + + alpha_prime = ttday0_res[1] + (n * (a + b + c * n)) / 2 + delta_prime = ttday0_res[2] + (n * (ap + bp + cp * n)) / 2 + Hp = (vs + lon - alpha_prime) % 360 + Hp[Hp >= 180] = Hp[Hp >= 180] - 360 + + h = np.degrees(np.arcsin(np.sin(np.radians(lat)) * + np.sin(np.radians(delta_prime)) + + np.cos(np.radians(lat)) * + np.cos(np.radians(delta_prime)) + * np.cos(np.radians(Hp)))) + + T = (m[0] - Hp[0] / 360) * 86400 + R = (m[1] + (h[1] + 0.8333) / (360 * np.cos(np.radians(delta_prime[1])) * + np.cos(np.radians(lat)) * + np.sin(np.radians(Hp[1])))) * 86400 + S = (m[2] + (h[2] + 0.8333) / (360 * np.cos(np.radians(delta_prime[2])) * + np.cos(np.radians(lat)) * + np.sin(np.radians(Hp[2])))) * 86400 + + S[add_a_day] += 86400 + R[sub_a_day] -= 86400 + + transit = T + utday + sunrise = R + utday + sunset = S + utday + + return transit, sunrise, sunset + + +def earthsun_distance(unixtime, delta_t, numthreads): + """ + Calculates the distance from the earth to the sun using the + NREL SPA algorithm described in [1]. + + Parameters + ---------- + unixtime : numpy array + Array of unix/epoch timestamps to calculate solar position for. + Unixtime is the number of seconds since Jan. 1, 1970 00:00:00 UTC. + A pandas.DatetimeIndex is easily converted using .view(np.int64)/10**9 + delta_t : float or array + Difference between terrestrial time and UT. USNO has tables. + numthreads : int + Number to threads to use for calculation (if using numba) + + Returns + ------- + R : array + Earth-Sun distance in AU. + + References + ---------- + [1] Reda, I., Andreas, A., 2003. Solar position algorithm for solar + radiation applications. Technical report: NREL/TP-560-34302. Golden, + USA, http://www.nlr.gov. + """ + + R = solar_position(unixtime, 0, 0, 0, 0, 0, delta_t, + 0, numthreads, esd=True)[0] + + return R + + +def calculate_deltat(year, month): + """Calculate the difference between Terrestrial Dynamical Time (TD) + and Universal Time (UT). + + Equations taken from http://eclipse.gsfc.nasa.gov/SEcat5/deltatpoly.html + """ + + plw = 'Deltat is unknown for years before -1999 and after 3000. ' \ + 'Delta values will be calculated, but the calculations ' \ + 'are not intended to be used for these years.' + + try: + if np.any((year > 3000) | (year < -1999)): + warnings.warn(plw) + except ValueError: + if (year > 3000) | (year < -1999): + warnings.warn(plw) + except TypeError: + return 0 + + y = year + (month - 0.5)/12 + + deltat = np.where(year < -500, + + -20+32*((y-1820)/100)**2, 0) + + deltat = np.where((-500 <= year) & (year < 500), + + 10583.6-1014.41*(y/100) + + 33.78311*(y/100)**2 + - 5.952053*(y/100)**3 + - 0.1798452*(y/100)**4 + + 0.022174192*(y/100)**5 + + 0.0090316521*(y/100)**6, deltat) + + deltat = np.where((500 <= year) & (year < 1600), + + 1574.2-556.01*((y-1000)/100) + + 71.23472*((y-1000)/100)**2 + + 0.319781*((y-1000)/100)**3 + - 0.8503463*((y-1000)/100)**4 + - 0.005050998*((y-1000)/100)**5 + + 0.0083572073*((y-1000)/100)**6, deltat) + + deltat = np.where((1600 <= year) & (year < 1700), + + 120-0.9808*(y-1600) + - 0.01532*(y-1600)**2 + + (y-1600)**3/7129, deltat) + + deltat = np.where((1700 <= year) & (year < 1800), + + 8.83+0.1603*(y-1700) + - 0.0059285*(y-1700)**2 + + 0.00013336*(y-1700)**3 + - (y-1700)**4/1174000, deltat) + + deltat = np.where((1800 <= year) & (year < 1860), + + 13.72-0.332447*(y-1800) + + 0.0068612*(y-1800)**2 + + 0.0041116*(y-1800)**3 + - 0.00037436*(y-1800)**4 + + 0.0000121272*(y-1800)**5 + - 0.0000001699*(y-1800)**6 + + 0.000000000875*(y-1800)**7, deltat) + + deltat = np.where((1860 <= year) & (year < 1900), + + 7.62+0.5737*(y-1860) + - 0.251754*(y-1860)**2 + + 0.01680668*(y-1860)**3 + - 0.0004473624*(y-1860)**4 + + (y-1860)**5/233174, deltat) + + deltat = np.where((1900 <= year) & (year < 1920), + + -2.79+1.494119*(y-1900) + - 0.0598939*(y-1900)**2 + + 0.0061966*(y-1900)**3 + - 0.000197*(y-1900)**4, deltat) + + deltat = np.where((1920 <= year) & (year < 1941), + + 21.20+0.84493*(y-1920) + - 0.076100*(y-1920)**2 + + 0.0020936*(y-1920)**3, deltat) + + deltat = np.where((1941 <= year) & (year < 1961), + + 29.07+0.407*(y-1950) + - (y-1950)**2/233 + + (y-1950)**3/2547, deltat) + + deltat = np.where((1961 <= year) & (year < 1986), + + 45.45+1.067*(y-1975) + - (y-1975)**2/260 + - (y-1975)**3/718, deltat) + + deltat = np.where((1986 <= year) & (year < 2005), + + 63.86+0.3345*(y-2000) + - 0.060374*(y-2000)**2 + + 0.0017275*(y-2000)**3 + + 0.000651814*(y-2000)**4 + + 0.00002373599*(y-2000)**5, deltat) + + deltat = np.where((2005 <= year) & (year < 2050), + + 62.92+0.32217*(y-2000) + + 0.005589*(y-2000)**2, deltat) + + deltat = np.where((2050 <= year) & (year < 2150), + + -20+32*((y-1820)/100)**2 + - 0.5628*(2150-y), deltat) + + deltat = np.where(year >= 2150, + + -20+32*((y-1820)/100)**2, deltat) + + deltat = deltat.item() if np.isscalar(year) & np.isscalar(month)\ + else deltat + + return deltat From fc56783a4d8e846e6a04ad143865de021e367235 Mon Sep 17 00:00:00 2001 From: Kevin Anderson Date: Wed, 28 Oct 2026 10:01:39 -0400 Subject: [PATCH 5/7] Revert "Apply batched suggestions from code review" This reverts commit 2a420ba6280b4085a784f4d0ab7693c691a177cd. --- pvlib/spa.py | 2753 +++++++++++++++++++++++++------------------------- 1 file changed, 1378 insertions(+), 1375 deletions(-) diff --git a/pvlib/spa.py b/pvlib/spa.py index 6536b7747c..6e36a084cd 100644 --- a/pvlib/spa.py +++ b/pvlib/spa.py @@ -1,1137 +1,1140 @@ -""" -Calculate the sun's position using the NREL Solar Position Algorithm -(SPA) either using numpy arrays or compiling the code to machine language -with numba. -""" - -# Contributors: -# Created by Tony Lorenzo (@alorenzo175), Univ. of Arizona, 2015 - -import os -import threading -import warnings - -import numpy as np - - -# this block is a way to use an environment variable to switch between -# compiling the functions with numba or just use numpy -def nocompile(*args, **kwargs): - return lambda func: func - - -if os.getenv('PVLIB_USE_NUMBA', '0') != '0': - try: - from numba import jit - except ImportError: - warnings.warn('Could not import numba, falling back to numpy ' + - 'calculation') - jcompile = nocompile - USE_NUMBA = False - else: - jcompile = jit - USE_NUMBA = True -else: - jcompile = nocompile - USE_NUMBA = False - - -# heliocentric longitude coefficients -L0 = np.array([ - [175347046.0, 0.0, 0.0], - [3341656.0, 4.6692568, 6283.07585], - [34894.0, 4.6261, 12566.1517], - [3497.0, 2.7441, 5753.3849], - [3418.0, 2.8289, 3.5231], - [3136.0, 3.6277, 77713.7715], - [2676.0, 4.4181, 7860.4194], - [2343.0, 6.1352, 3930.2097], - [1324.0, 0.7425, 11506.7698], - [1273.0, 2.0371, 529.691], - [1199.0, 1.1096, 1577.3435], - [990.0, 5.233, 5884.927], - [902.0, 2.045, 26.298], - [857.0, 3.508, 398.149], - [780.0, 1.179, 5223.694], - [753.0, 2.533, 5507.553], - [505.0, 4.583, 18849.228], - [492.0, 4.205, 775.523], - [357.0, 2.92, 0.067], - [317.0, 5.849, 11790.629], - [284.0, 1.899, 796.298], - [271.0, 0.315, 10977.079], - [243.0, 0.345, 5486.778], - [206.0, 4.806, 2544.314], - [205.0, 1.869, 5573.143], - [202.0, 2.458, 6069.777], - [156.0, 0.833, 213.299], - [132.0, 3.411, 2942.463], - [126.0, 1.083, 20.775], - [115.0, 0.645, 0.98], - [103.0, 0.636, 4694.003], - [102.0, 0.976, 15720.839], - [102.0, 4.267, 7.114], - [99.0, 6.21, 2146.17], - [98.0, 0.68, 155.42], - [86.0, 5.98, 161000.69], - [85.0, 1.3, 6275.96], - [85.0, 3.67, 71430.7], - [80.0, 1.81, 17260.15], - [79.0, 3.04, 12036.46], - [75.0, 1.76, 5088.63], - [74.0, 3.5, 3154.69], - [74.0, 4.68, 801.82], - [70.0, 0.83, 9437.76], - [62.0, 3.98, 8827.39], - [61.0, 1.82, 7084.9], - [57.0, 2.78, 6286.6], - [56.0, 4.39, 14143.5], - [56.0, 3.47, 6279.55], - [52.0, 0.19, 12139.55], - [52.0, 1.33, 1748.02], - [51.0, 0.28, 5856.48], - [49.0, 0.49, 1194.45], - [41.0, 5.37, 8429.24], - [41.0, 2.4, 19651.05], - [39.0, 6.17, 10447.39], - [37.0, 6.04, 10213.29], - [37.0, 2.57, 1059.38], - [36.0, 1.71, 2352.87], - [36.0, 1.78, 6812.77], - [33.0, 0.59, 17789.85], - [30.0, 0.44, 83996.85], - [30.0, 2.74, 1349.87], - [25.0, 3.16, 4690.48] -]) -L1 = np.array([ - [628331966747.0, 0.0, 0.0], - [206059.0, 2.678235, 6283.07585], - [4303.0, 2.6351, 12566.1517], - [425.0, 1.59, 3.523], - [119.0, 5.796, 26.298], - [109.0, 2.966, 1577.344], - [93.0, 2.59, 18849.23], - [72.0, 1.14, 529.69], - [68.0, 1.87, 398.15], - [67.0, 4.41, 5507.55], - [59.0, 2.89, 5223.69], - [56.0, 2.17, 155.42], - [45.0, 0.4, 796.3], - [36.0, 0.47, 775.52], - [29.0, 2.65, 7.11], - [21.0, 5.34, 0.98], - [19.0, 1.85, 5486.78], - [19.0, 4.97, 213.3], - [17.0, 2.99, 6275.96], - [16.0, 0.03, 2544.31], - [16.0, 1.43, 2146.17], - [15.0, 1.21, 10977.08], - [12.0, 2.83, 1748.02], - [12.0, 3.26, 5088.63], - [12.0, 5.27, 1194.45], - [12.0, 2.08, 4694.0], - [11.0, 0.77, 553.57], - [10.0, 1.3, 6286.6], - [10.0, 4.24, 1349.87], - [9.0, 2.7, 242.73], - [9.0, 5.64, 951.72], - [8.0, 5.3, 2352.87], - [6.0, 2.65, 9437.76], - [6.0, 4.67, 4690.48] -]) -L2 = np.array([ - [52919.0, 0.0, 0.0], - [8720.0, 1.0721, 6283.0758], - [309.0, 0.867, 12566.152], - [27.0, 0.05, 3.52], - [16.0, 5.19, 26.3], - [16.0, 3.68, 155.42], - [10.0, 0.76, 18849.23], - [9.0, 2.06, 77713.77], - [7.0, 0.83, 775.52], - [5.0, 4.66, 1577.34], - [4.0, 1.03, 7.11], - [4.0, 3.44, 5573.14], - [3.0, 5.14, 796.3], - [3.0, 6.05, 5507.55], - [3.0, 1.19, 242.73], - [3.0, 6.12, 529.69], - [3.0, 0.31, 398.15], - [3.0, 2.28, 553.57], - [2.0, 4.38, 5223.69], - [2.0, 3.75, 0.98] -]) -L3 = np.array([ - [289.0, 5.844, 6283.076], - [35.0, 0.0, 0.0], - [17.0, 5.49, 12566.15], - [3.0, 5.2, 155.42], - [1.0, 4.72, 3.52], - [1.0, 5.3, 18849.23], - [1.0, 5.97, 242.73] -]) -L4 = np.array([ - [114.0, 3.142, 0.0], - [8.0, 4.13, 6283.08], - [1.0, 3.84, 12566.15] -]) -L5 = np.array([ - [1.0, 3.14, 0.0] -]) - - -# heliocentric latitude coefficients -B0 = np.array([ - [280.0, 3.199, 84334.662], - [102.0, 5.422, 5507.553], - [80.0, 3.88, 5223.69], - [44.0, 3.7, 2352.87], - [32.0, 4.0, 1577.34] -]) -B1 = np.array([ - [9.0, 3.9, 5507.55], - [6.0, 1.73, 5223.69] -]) - - -# heliocentric radius coefficients -R0 = np.array([ - [100013989.0, 0.0, 0.0], - [1670700.0, 3.0984635, 6283.07585], - [13956.0, 3.05525, 12566.1517], - [3084.0, 5.1985, 77713.7715], - [1628.0, 1.1739, 5753.3849], - [1576.0, 2.8469, 7860.4194], - [925.0, 5.453, 11506.77], - [542.0, 4.564, 3930.21], - [472.0, 3.661, 5884.927], - [346.0, 0.964, 5507.553], - [329.0, 5.9, 5223.694], - [307.0, 0.299, 5573.143], - [243.0, 4.273, 11790.629], - [212.0, 5.847, 1577.344], - [186.0, 5.022, 10977.079], - [175.0, 3.012, 18849.228], - [110.0, 5.055, 5486.778], - [98.0, 0.89, 6069.78], - [86.0, 5.69, 15720.84], - [86.0, 1.27, 161000.69], - [65.0, 0.27, 17260.15], - [63.0, 0.92, 529.69], - [57.0, 2.01, 83996.85], - [56.0, 5.24, 71430.7], - [49.0, 3.25, 2544.31], - [47.0, 2.58, 775.52], - [45.0, 5.54, 9437.76], - [43.0, 6.01, 6275.96], - [39.0, 5.36, 4694.0], - [38.0, 2.39, 8827.39], - [37.0, 0.83, 19651.05], - [37.0, 4.9, 12139.55], - [36.0, 1.67, 12036.46], - [35.0, 1.84, 2942.46], - [33.0, 0.24, 7084.9], - [32.0, 0.18, 5088.63], - [32.0, 1.78, 398.15], - [28.0, 1.21, 6286.6], - [28.0, 1.9, 6279.55], - [26.0, 4.59, 10447.39] -]) -R1 = np.array([ - [103019.0, 1.10749, 6283.07585], - [1721.0, 1.0644, 12566.1517], - [702.0, 3.142, 0.0], - [32.0, 1.02, 18849.23], - [31.0, 2.84, 5507.55], - [25.0, 1.32, 5223.69], - [18.0, 1.42, 1577.34], - [10.0, 5.91, 10977.08], - [9.0, 1.42, 6275.96], - [9.0, 0.27, 5486.78] -]) -R2 = np.array([ - [4359.0, 5.7846, 6283.0758], - [124.0, 5.579, 12566.152], - [12.0, 3.14, 0.0], - [9.0, 3.63, 77713.77], - [6.0, 1.87, 5573.14], - [3.0, 5.47, 18849.23] -]) -R3 = np.array([ - [145.0, 4.273, 6283.076], - [7.0, 3.92, 12566.15] -]) -R4 = np.array([ - [4.0, 2.56, 6283.08] -]) - - -# longitude and obliquity nutation coefficients -NUTATION_ABCD_ARRAY = np.array([ - [-171996, -174.2, 92025, 8.9], - [-13187, -1.6, 5736, -3.1], - [-2274, -0.2, 977, -0.5], - [2062, 0.2, -895, 0.5], - [1426, -3.4, 54, -0.1], - [712, 0.1, -7, 0], - [-517, 1.2, 224, -0.6], - [-386, -0.4, 200, 0], - [-301, 0, 129, -0.1], - [217, -0.5, -95, 0.3], - [-158, 0, 0, 0], - [129, 0.1, -70, 0], - [123, 0, -53, 0], - [63, 0, 0, 0], - [63, 0.1, -33, 0], - [-59, 0, 26, 0], - [-58, -0.1, 32, 0], - [-51, 0, 27, 0], - [48, 0, 0, 0], - [46, 0, -24, 0], - [-38, 0, 16, 0], - [-31, 0, 13, 0], - [29, 0, 0, 0], - [29, 0, -12, 0], - [26, 0, 0, 0], - [-22, 0, 0, 0], - [21, 0, -10, 0], - [17, -0.1, 0, 0], - [16, 0, -8, 0], - [-16, 0.1, 7, 0], - [-15, 0, 9, 0], - [-13, 0, 7, 0], - [-12, 0, 6, 0], - [11, 0, 0, 0], - [-10, 0, 5, 0], - [-8, 0, 3, 0], - [7, 0, -3, 0], - [-7, 0, 0, 0], - [-7, 0, 3, 0], - [-7, 0, 3, 0], - [6, 0, 0, 0], - [6, 0, -3, 0], - [6, 0, -3, 0], - [-6, 0, 3, 0], - [-6, 0, 3, 0], - [5, 0, 0, 0], - [-5, 0, 3, 0], - [-5, 0, 3, 0], - [-5, 0, 3, 0], - [4, 0, 0, 0], - [4, 0, 0, 0], - [4, 0, 0, 0], - [-4, 0, 0, 0], - [-4, 0, 0, 0], - [-4, 0, 0, 0], - [3, 0, 0, 0], - [-3, 0, 0, 0], - [-3, 0, 0, 0], - [-3, 0, 0, 0], - [-3, 0, 0, 0], - [-3, 0, 0, 0], - [-3, 0, 0, 0], - [-3, 0, 0, 0], -]) - -NUTATION_YTERM_ARRAY = np.array([ - [0, 0, 0, 0, 1], - [-2, 0, 0, 2, 2], - [0, 0, 0, 2, 2], - [0, 0, 0, 0, 2], - [0, 1, 0, 0, 0], - [0, 0, 1, 0, 0], - [-2, 1, 0, 2, 2], - [0, 0, 0, 2, 1], - [0, 0, 1, 2, 2], - [-2, -1, 0, 2, 2], - [-2, 0, 1, 0, 0], - [-2, 0, 0, 2, 1], - [0, 0, -1, 2, 2], - [2, 0, 0, 0, 0], - [0, 0, 1, 0, 1], - [2, 0, -1, 2, 2], - [0, 0, -1, 0, 1], - [0, 0, 1, 2, 1], - [-2, 0, 2, 0, 0], - [0, 0, -2, 2, 1], - [2, 0, 0, 2, 2], - [0, 0, 2, 2, 2], - [0, 0, 2, 0, 0], - [-2, 0, 1, 2, 2], - [0, 0, 0, 2, 0], - [-2, 0, 0, 2, 0], - [0, 0, -1, 2, 1], - [0, 2, 0, 0, 0], - [2, 0, -1, 0, 1], - [-2, 2, 0, 2, 2], - [0, 1, 0, 0, 1], - [-2, 0, 1, 0, 1], - [0, -1, 0, 0, 1], - [0, 0, 2, -2, 0], - [2, 0, -1, 2, 1], - [2, 0, 1, 2, 2], - [0, 1, 0, 2, 2], - [-2, 1, 1, 0, 0], - [0, -1, 0, 2, 2], - [2, 0, 0, 2, 1], - [2, 0, 1, 0, 0], - [-2, 0, 2, 2, 2], - [-2, 0, 1, 2, 1], - [2, 0, -2, 0, 1], - [2, 0, 0, 0, 1], - [0, -1, 1, 0, 0], - [-2, -1, 0, 2, 1], - [-2, 0, 0, 0, 1], - [0, 0, 2, 2, 1], - [-2, 0, 2, 0, 1], - [-2, 1, 0, 2, 1], - [0, 0, 1, -2, 0], - [-1, 0, 1, 0, 0], - [-2, 1, 0, 0, 0], - [1, 0, 0, 0, 0], - [0, 0, 1, 2, 0], - [0, 0, -2, 2, 2], - [-1, -1, 1, 0, 0], - [0, 1, 1, 0, 0], - [0, -1, 1, 2, 2], - [2, -1, -1, 2, 2], - [0, 0, 3, 2, 2], - [2, -1, 0, 2, 2], -]) - - -@jcompile('float64(int64, int64, int64, int64, int64, int64, int64)', - nopython=True) -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 - - -@jcompile('float64(float64)', nopython=True) -def julian_day(unixtime): - jd = unixtime * 1.0 / 86400 + 2440587.5 - return jd - - -@jcompile('float64(float64, float64)', nopython=True) -def julian_ephemeris_day(julian_day, delta_t): - jde = julian_day + delta_t * 1.0 / 86400 - return jde - - -@jcompile('float64(float64)', nopython=True) -def julian_century(julian_day): - jc = (julian_day - 2451545) * 1.0 / 36525 - return jc - - -@jcompile('float64(float64)', nopython=True) -def julian_ephemeris_century(julian_ephemeris_day): - jce = (julian_ephemeris_day - 2451545) * 1.0 / 36525 - return jce - - -@jcompile('float64(float64)', nopython=True) -def julian_ephemeris_millennium(julian_ephemeris_century): - jme = julian_ephemeris_century * 1.0 / 10 - return jme - - -# omit type signature here; specifying read-only arrays requires use of the -# numba.types API, meaning numba must be available to import. -# https://github.com/numba/numba/issues/4511 -@jcompile(nopython=True) -def sum_mult_cos_add_mult(arr, x): - # shared calculation used for heliocentric longitude, latitude, and radius - s = 0. - for row in range(arr.shape[0]): - s += arr[row, 0] * np.cos(arr[row, 1] + arr[row, 2] * x) - return s - -@jcompile('float64(float64)', nopython=True) -def heliocentric_longitude(jme): - l0 = sum_mult_cos_add_mult(L0, jme) - l1 = sum_mult_cos_add_mult(L1, jme) - l2 = sum_mult_cos_add_mult(L2, jme) - l3 = sum_mult_cos_add_mult(L3, jme) - l4 = sum_mult_cos_add_mult(L4, jme) - l5 = sum_mult_cos_add_mult(L5, jme) - - l_rad = (l0 + l1 * jme + l2 * jme**2 + l3 * jme**3 + l4 * jme**4 + - l5 * jme**5)/10**8 - l = np.rad2deg(l_rad) - return l % 360 - -@jcompile('float64(float64)', nopython=True) -def heliocentric_latitude(jme): - b0 = sum_mult_cos_add_mult(B0, jme) - b1 = sum_mult_cos_add_mult(B1, jme) - - b_rad = (b0 + b1 * jme)/10**8 - b = np.rad2deg(b_rad) - return b - - -@jcompile('float64(float64)', nopython=True) -def heliocentric_radius_vector(jme): - r0 = sum_mult_cos_add_mult(R0, jme) - r1 = sum_mult_cos_add_mult(R1, jme) - r2 = sum_mult_cos_add_mult(R2, jme) - r3 = sum_mult_cos_add_mult(R3, jme) - r4 = sum_mult_cos_add_mult(R4, jme) - - r = (r0 + r1 * jme + r2 * jme**2 + r3 * jme**3 + r4 * jme**4)/10**8 - return r - - -@jcompile('float64(float64)', nopython=True) -def geocentric_longitude(heliocentric_longitude): - theta = heliocentric_longitude + 180.0 - return theta % 360 - - -@jcompile('float64(float64)', nopython=True) -def geocentric_latitude(heliocentric_latitude): - beta = -1.0*heliocentric_latitude - return beta - - -@jcompile('float64(float64)', nopython=True) -def mean_elongation(julian_ephemeris_century): - x0 = (297.85036 - + 445267.111480 * julian_ephemeris_century - - 0.0019142 * julian_ephemeris_century**2 - + julian_ephemeris_century**3 / 189474) - return x0 - - -@jcompile('float64(float64)', nopython=True) -def mean_anomaly_sun(julian_ephemeris_century): - x1 = (357.52772 - + 35999.050340 * julian_ephemeris_century - - 0.0001603 * julian_ephemeris_century**2 - - julian_ephemeris_century**3 / 300000) - return x1 - - -@jcompile('float64(float64)', nopython=True) -def mean_anomaly_moon(julian_ephemeris_century): - x2 = (134.96298 - + 477198.867398 * julian_ephemeris_century - + 0.0086972 * julian_ephemeris_century**2 - + julian_ephemeris_century**3 / 56250) - return x2 - - -@jcompile('float64(float64)', nopython=True) -def moon_argument_latitude(julian_ephemeris_century): - x3 = (93.27191 - + 483202.017538 * julian_ephemeris_century - - 0.0036825 * julian_ephemeris_century**2 - + julian_ephemeris_century**3 / 327270) - return x3 - - -@jcompile('float64(float64)', nopython=True) -def moon_ascending_longitude(julian_ephemeris_century): - x4 = (125.04452 - - 1934.136261 * julian_ephemeris_century - + 0.0020708 * julian_ephemeris_century**2 - + julian_ephemeris_century**3 / 450000) - return x4 - - -@jcompile( - 'void(float64, float64, float64, float64, float64, float64, float64[:])', - nopython=True) -def longitude_obliquity_nutation(julian_ephemeris_century, x0, x1, x2, x3, x4, - out): - delta_psi_sum = 0.0 - delta_eps_sum = 0.0 - for row in range(NUTATION_YTERM_ARRAY.shape[0]): - a = NUTATION_ABCD_ARRAY[row, 0] - b = NUTATION_ABCD_ARRAY[row, 1] - c = NUTATION_ABCD_ARRAY[row, 2] - d = NUTATION_ABCD_ARRAY[row, 3] - arg = np.radians( - NUTATION_YTERM_ARRAY[row, 0]*x0 + - NUTATION_YTERM_ARRAY[row, 1]*x1 + - NUTATION_YTERM_ARRAY[row, 2]*x2 + - NUTATION_YTERM_ARRAY[row, 3]*x3 + - NUTATION_YTERM_ARRAY[row, 4]*x4 - ) - delta_psi_sum += (a + b * julian_ephemeris_century) * np.sin(arg) - delta_eps_sum += (c + d * julian_ephemeris_century) * np.cos(arg) - delta_psi = delta_psi_sum*1.0/36000000 - delta_eps = delta_eps_sum*1.0/36000000 - # seems like we ought to be able to return a tuple here instead - # of resorting to `out`, but returning a UniTuple from this - # function caused calculations elsewhere to give the wrong result. - # very difficult to investigate since it did not occur when using - # object mode. issue was observed on numba 0.56.4 - out[0] = delta_psi - out[1] = delta_eps - - -@jcompile('float64(float64)', nopython=True) -def mean_ecliptic_obliquity(julian_ephemeris_millennium): - U = 1.0*julian_ephemeris_millennium/10 - e0 = (84381.448 - 4680.93 * U - 1.55 * U**2 - + 1999.25 * U**3 - 51.38 * U**4 - 249.67 * U**5 - - 39.05 * U**6 + 7.12 * U**7 + 27.87 * U**8 - + 5.79 * U**9 + 2.45 * U**10) - return e0 - - -@jcompile('float64(float64, float64)', nopython=True) -def true_ecliptic_obliquity(mean_ecliptic_obliquity, obliquity_nutation): - e0 = mean_ecliptic_obliquity - deleps = obliquity_nutation - e = e0*1.0/3600 + deleps - return e - - -@jcompile('float64(float64)', nopython=True) -def aberration_correction(earth_radius_vector): - deltau = -20.4898 / (3600 * earth_radius_vector) - return deltau - - -@jcompile('float64(float64, float64, float64)', nopython=True) -def apparent_sun_longitude(geocentric_longitude, longitude_nutation, - aberration_correction): - lamd = geocentric_longitude + longitude_nutation + aberration_correction - return lamd - - -@jcompile('float64(float64, float64)', nopython=True) -def mean_sidereal_time(julian_day, julian_century): - v0 = (280.46061837 + 360.98564736629 * (julian_day - 2451545) - + 0.000387933 * julian_century**2 - julian_century**3 / 38710000) - return v0 % 360.0 - - -@jcompile('float64(float64, float64, float64)', nopython=True) -def apparent_sidereal_time(mean_sidereal_time, longitude_nutation, - true_ecliptic_obliquity): - v = mean_sidereal_time + longitude_nutation * np.cos( - np.radians(true_ecliptic_obliquity)) - return v - - -@jcompile('float64(float64, float64, float64)', nopython=True) -def geocentric_sun_right_ascension(apparent_sun_longitude, - true_ecliptic_obliquity, - geocentric_latitude): - true_ecliptic_obliquity_rad = np.radians(true_ecliptic_obliquity) - apparent_sun_longitude_rad = np.radians(apparent_sun_longitude) - - num = (np.sin(apparent_sun_longitude_rad) - * np.cos(true_ecliptic_obliquity_rad) - - np.tan(np.radians(geocentric_latitude)) - * np.sin(true_ecliptic_obliquity_rad)) - alpha = np.degrees(np.arctan2(num, np.cos(apparent_sun_longitude_rad))) - return alpha % 360 - - -@jcompile('float64(float64, float64, float64)', nopython=True) -def geocentric_sun_declination(apparent_sun_longitude, true_ecliptic_obliquity, - geocentric_latitude): - geocentric_latitude_rad = np.radians(geocentric_latitude) - true_ecliptic_obliquity_rad = np.radians(true_ecliptic_obliquity) - - delta = np.degrees(np.arcsin(np.sin(geocentric_latitude_rad) * - np.cos(true_ecliptic_obliquity_rad) + - np.cos(geocentric_latitude_rad) * - np.sin(true_ecliptic_obliquity_rad) * - np.sin(np.radians(apparent_sun_longitude)))) - return delta - - -@jcompile('float64(float64, float64, float64)', nopython=True) -def local_hour_angle(apparent_sidereal_time, observer_longitude, - sun_right_ascension): - """Measured westward from south""" - H = apparent_sidereal_time + observer_longitude - sun_right_ascension - return H % 360 - - -@jcompile('float64(float64)', nopython=True) -def equatorial_horizontal_parallax(earth_radius_vector): - xi = 8.794 / (3600 * earth_radius_vector) - return xi - - -@jcompile('float64(float64)', nopython=True) -def uterm(observer_latitude): - u = np.arctan(0.99664719 * np.tan(np.radians(observer_latitude))) - return u - - -@jcompile('float64(float64, float64, float64)', nopython=True) -def xterm(u, observer_latitude, observer_elevation): - x = (np.cos(u) + observer_elevation / 6378140 - * np.cos(np.radians(observer_latitude))) - return x - - -@jcompile('float64(float64, float64, float64)', nopython=True) -def yterm(u, observer_latitude, observer_elevation): - y = (0.99664719 * np.sin(u) + observer_elevation / 6378140 - * np.sin(np.radians(observer_latitude))) - return y - - -@jcompile('float64(float64, float64,float64, float64)', nopython=True) -def parallax_sun_right_ascension(xterm, equatorial_horizontal_parallax, - local_hour_angle, geocentric_sun_declination): - equatorial_horizontal_parallax_rad = \ - np.radians(equatorial_horizontal_parallax) - local_hour_angle_rad = np.radians(local_hour_angle) - - num = (-xterm * np.sin(equatorial_horizontal_parallax_rad) - * np.sin(local_hour_angle_rad)) - denom = (np.cos(np.radians(geocentric_sun_declination)) - - xterm * np.sin(equatorial_horizontal_parallax_rad) - * np.cos(local_hour_angle_rad)) - delta_alpha = np.degrees(np.arctan2(num, denom)) - return delta_alpha - - -@jcompile('float64(float64, float64)', nopython=True) -def topocentric_sun_right_ascension(geocentric_sun_right_ascension, - parallax_sun_right_ascension): - alpha_prime = geocentric_sun_right_ascension + parallax_sun_right_ascension - return alpha_prime - - -@jcompile('float64(float64, float64, float64, float64, float64, float64)', - nopython=True) -def topocentric_sun_declination(geocentric_sun_declination, xterm, yterm, - equatorial_horizontal_parallax, - parallax_sun_right_ascension, - local_hour_angle): - geocentric_sun_declination_rad = np.radians(geocentric_sun_declination) - equatorial_horizontal_parallax_rad = \ - np.radians(equatorial_horizontal_parallax) - - num = ((np.sin(geocentric_sun_declination_rad) - yterm - * np.sin(equatorial_horizontal_parallax_rad)) - * np.cos(np.radians(parallax_sun_right_ascension))) - denom = (np.cos(geocentric_sun_declination_rad) - xterm - * np.sin(equatorial_horizontal_parallax_rad) - * np.cos(np.radians(local_hour_angle))) - delta = np.degrees(np.arctan2(num, denom)) - return delta - - -@jcompile('float64(float64, float64)', nopython=True) -def topocentric_local_hour_angle(local_hour_angle, - parallax_sun_right_ascension): - H_prime = local_hour_angle - parallax_sun_right_ascension - return H_prime - - -@jcompile('float64(float64, float64, float64)', nopython=True) -def topocentric_elevation_angle_without_atmosphere(observer_latitude, - topocentric_sun_declination, - topocentric_local_hour_angle - ): - observer_latitude_rad = np.radians(observer_latitude) - topocentric_sun_declination_rad = np.radians(topocentric_sun_declination) - e0 = np.degrees(np.arcsin( - np.sin(observer_latitude_rad) - * np.sin(topocentric_sun_declination_rad) - + np.cos(observer_latitude_rad) - * np.cos(topocentric_sun_declination_rad) - * np.cos(np.radians(topocentric_local_hour_angle)))) - return e0 - - -@jcompile('float64(float64, float64, float64, float64)', nopython=True) -def atmospheric_refraction_correction(local_pressure, local_temp, - topocentric_elevation_angle_wo_atmosphere, - atmos_refract): - # switch sets delta_e when the sun is below the horizon - switch = topocentric_elevation_angle_wo_atmosphere >= -1.0 * ( - 0.26667 + atmos_refract) - delta_e = ((local_pressure / 1010.0) * (283.0 / (273 + local_temp)) - * 1.02 / (60 * np.tan(np.radians( - topocentric_elevation_angle_wo_atmosphere - + 10.3 / (topocentric_elevation_angle_wo_atmosphere - + 5.11))))) * switch - return delta_e - - -@jcompile('float64(float64, float64)', nopython=True) -def topocentric_elevation_angle(topocentric_elevation_angle_without_atmosphere, - atmospheric_refraction_correction): - e = (topocentric_elevation_angle_without_atmosphere - + atmospheric_refraction_correction) - return e - - -@jcompile('float64(float64)', nopython=True) -def topocentric_zenith_angle(topocentric_elevation_angle): - theta = 90 - topocentric_elevation_angle - return theta - - -@jcompile('float64(float64, float64, float64)', nopython=True) -def topocentric_astronomers_azimuth(topocentric_local_hour_angle, - topocentric_sun_declination, - observer_latitude): - topocentric_local_hour_angle_rad = np.radians(topocentric_local_hour_angle) - observer_latitude_rad = np.radians(observer_latitude) - num = np.sin(topocentric_local_hour_angle_rad) - denom = (np.cos(topocentric_local_hour_angle_rad) - * np.sin(observer_latitude_rad) - - np.tan(np.radians(topocentric_sun_declination)) - * np.cos(observer_latitude_rad)) - gamma = np.degrees(np.arctan2(num, denom)) - return gamma % 360 - - -@jcompile('float64(float64)', nopython=True) -def topocentric_azimuth_angle(topocentric_astronomers_azimuth): - phi = topocentric_astronomers_azimuth + 180 - return phi % 360 - - -@jcompile('float64(float64)', nopython=True) -def sun_mean_longitude(julian_ephemeris_millennium): - M = (280.4664567 + 360007.6982779 * julian_ephemeris_millennium - + 0.03032028 * julian_ephemeris_millennium**2 - + julian_ephemeris_millennium**3 / 49931 - - julian_ephemeris_millennium**4 / 15300 - - julian_ephemeris_millennium**5 / 2000000) - return M - - -@jcompile('float64(float64, float64, float64, float64)', nopython=True) -def equation_of_time(sun_mean_longitude, geocentric_sun_right_ascension, - longitude_nutation, true_ecliptic_obliquity): - E = (sun_mean_longitude - 0.0057183 - geocentric_sun_right_ascension + - longitude_nutation * np.cos(np.radians(true_ecliptic_obliquity))) - # limit between 0 and 360 - E = E % 360 - # convert to minutes - E *= 4 - greater = E > 20 - less = E < -20 - other = (E <= 20) & (E >= -20) - E = greater * (E - 1440) + less * (E + 1440) + other * E - return E - - -@jcompile('void(float64[:], float64[:], float64[:], float64[:,:])', - nopython=True, nogil=True) -def solar_position_loop(unixtime, delta_t, loc_args, out): - """Loop through the time array and calculate the solar position""" - lat = loc_args[0] - lon = loc_args[1] - elev = loc_args[2] - pressure = loc_args[3] - temp = loc_args[4] - atmos_refract = loc_args[5] - sst = loc_args[6] - esd = loc_args[7] - - for i in range(unixtime.shape[0]): - utime = unixtime[i] - dT = delta_t[i] - jd = julian_day(utime) - jde = julian_ephemeris_day(jd, dT) - jc = julian_century(jd) - jce = julian_ephemeris_century(jde) - jme = julian_ephemeris_millennium(jce) - R = heliocentric_radius_vector(jme) - if esd: - out[0, i] = R - continue - L = heliocentric_longitude(jme) - B = heliocentric_latitude(jme) - Theta = geocentric_longitude(L) - beta = geocentric_latitude(B) - x0 = mean_elongation(jce) - x1 = mean_anomaly_sun(jce) - x2 = mean_anomaly_moon(jce) - x3 = moon_argument_latitude(jce) - x4 = moon_ascending_longitude(jce) - l_o_nutation = np.empty((2,)) - longitude_obliquity_nutation(jce, x0, x1, x2, x3, x4, l_o_nutation) - delta_psi = l_o_nutation[0] - delta_epsilon = l_o_nutation[1] - epsilon0 = mean_ecliptic_obliquity(jme) - epsilon = true_ecliptic_obliquity(epsilon0, delta_epsilon) - delta_tau = aberration_correction(R) - lamd = apparent_sun_longitude(Theta, delta_psi, delta_tau) - v0 = mean_sidereal_time(jd, jc) - v = apparent_sidereal_time(v0, delta_psi, epsilon) - alpha = geocentric_sun_right_ascension(lamd, epsilon, beta) - delta = geocentric_sun_declination(lamd, epsilon, beta) - if sst: - out[0, i] = v - out[1, i] = alpha - out[2, i] = delta - continue - m = sun_mean_longitude(jme) - eot = equation_of_time(m, alpha, delta_psi, epsilon) - H = local_hour_angle(v, lon, alpha) - xi = equatorial_horizontal_parallax(R) - u = uterm(lat) - x = xterm(u, lat, elev) - y = yterm(u, lat, elev) - delta_alpha = parallax_sun_right_ascension(x, xi, H, delta) - delta_prime = topocentric_sun_declination(delta, x, y, xi, delta_alpha, - H) - H_prime = topocentric_local_hour_angle(H, delta_alpha) - e0 = topocentric_elevation_angle_without_atmosphere(lat, delta_prime, - H_prime) - delta_e = atmospheric_refraction_correction(pressure, temp, e0, - atmos_refract) - e = topocentric_elevation_angle(e0, delta_e) - theta = topocentric_zenith_angle(e) - theta0 = topocentric_zenith_angle(e0) - gamma = topocentric_astronomers_azimuth(H_prime, delta_prime, lat) - phi = topocentric_azimuth_angle(gamma) - out[0, i] = theta - out[1, i] = theta0 - out[2, i] = e - out[3, i] = e0 - out[4, i] = phi - out[5, i] = eot - - -def solar_position_numba(unixtime, lat, lon, elev, pressure, temp, delta_t, - atmos_refract, numthreads, sst=False, esd=False): - """Calculate the solar position using the numba compiled functions - and multiple threads. Very slow if functions are not numba compiled. - """ - # these args are the same for each thread - loc_args = np.array([lat, lon, elev, pressure, temp, - atmos_refract, sst, esd], dtype=np.float64) - - # turn delta_t into an array if it isn't already - delta_t = np.full_like(unixtime, delta_t, dtype=np.float64) - - # construct dims x ulength array to put the results in - ulength = unixtime.shape[0] - if sst: - dims = 3 - elif esd: - dims = 1 - else: - dims = 6 - result = np.empty((dims, ulength), dtype=np.float64) - - if unixtime.dtype != np.float64: - unixtime = unixtime.astype(np.float64) - - if ulength < numthreads: - numthreads = ulength - - if numthreads <= 1: - solar_position_loop(unixtime, delta_t, loc_args, result) - return result - - # split the input and output arrays into numthreads chunks - split0 = np.array_split(unixtime, numthreads) - split1 = np.array_split(delta_t, numthreads) - split2 = np.array_split(result, numthreads, axis=1) - chunks = [ - [a0, a1, loc_args, a2] - for a0, a1, a2 in zip(split0, split1, split2) - ] - # Spawn one thread per chunk - threads = [threading.Thread(target=solar_position_loop, args=chunk) - for chunk in chunks] - for thread in threads: - thread.start() - for thread in threads: - thread.join() - return result - - -def solar_position_numpy(unixtime, lat, lon, elev, pressure, temp, delta_t, - atmos_refract, numthreads, sst=False, esd=False): - """Calculate the solar position assuming unixtime is a numpy array. Note - this function will not work if the solar position functions were - compiled with numba. - """ - - jd = julian_day(unixtime) - jde = julian_ephemeris_day(jd, delta_t) - jc = julian_century(jd) - jce = julian_ephemeris_century(jde) - jme = julian_ephemeris_millennium(jce) - R = heliocentric_radius_vector(jme) - if esd: - return (R, ) - L = heliocentric_longitude(jme) - B = heliocentric_latitude(jme) - Theta = geocentric_longitude(L) - beta = geocentric_latitude(B) - x0 = mean_elongation(jce) - x1 = mean_anomaly_sun(jce) - x2 = mean_anomaly_moon(jce) - x3 = moon_argument_latitude(jce) - x4 = moon_ascending_longitude(jce) - l_o_nutation = np.empty((2, len(x0))) - longitude_obliquity_nutation(jce, x0, x1, x2, x3, x4, l_o_nutation) - delta_psi = l_o_nutation[0] - delta_epsilon = l_o_nutation[1] - epsilon0 = mean_ecliptic_obliquity(jme) - epsilon = true_ecliptic_obliquity(epsilon0, delta_epsilon) - delta_tau = aberration_correction(R) - lamd = apparent_sun_longitude(Theta, delta_psi, delta_tau) - v0 = mean_sidereal_time(jd, jc) - v = apparent_sidereal_time(v0, delta_psi, epsilon) - alpha = geocentric_sun_right_ascension(lamd, epsilon, beta) - delta = geocentric_sun_declination(lamd, epsilon, beta) - if sst: - return v, alpha, delta - m = sun_mean_longitude(jme) - eot = equation_of_time(m, alpha, delta_psi, epsilon) - H = local_hour_angle(v, lon, alpha) - xi = equatorial_horizontal_parallax(R) - u = uterm(lat) - x = xterm(u, lat, elev) - y = yterm(u, lat, elev) - delta_alpha = parallax_sun_right_ascension(x, xi, H, delta) - delta_prime = topocentric_sun_declination(delta, x, y, xi, delta_alpha, H) - H_prime = topocentric_local_hour_angle(H, delta_alpha) - e0 = topocentric_elevation_angle_without_atmosphere(lat, delta_prime, - H_prime) - delta_e = atmospheric_refraction_correction(pressure, temp, e0, - atmos_refract) - e = topocentric_elevation_angle(e0, delta_e) - theta = topocentric_zenith_angle(e) - theta0 = topocentric_zenith_angle(e0) - gamma = topocentric_astronomers_azimuth(H_prime, delta_prime, lat) - phi = topocentric_azimuth_angle(gamma) - return theta, theta0, e, e0, phi, eot - - -def solar_position(unixtime, lat, lon, elev, pressure, temp, delta_t, - atmos_refract, numthreads=8, sst=False, esd=False): - - """ - Calculate the solar position using the - NREL SPA described in [1]_, [2]_. - - If numba is installed, the functions can be compiled - and the code runs quickly. If not, the functions - still evaluate but use numpy instead. - - Parameters - ---------- - unixtime : numpy array - Array of unix/epoch timestamps to calculate solar position for. - Unixtime is the number of seconds since Jan. 1, 1970 00:00:00 UTC. - A pandas.DatetimeIndex is easily converted using .view(np.int64)/10**9 - lat : float - Latitude to calculate solar position for - lon : float - Longitude to calculate solar position for - elev : float - Elevation of location in meters - pressure : int or float - avg. yearly pressure at location in millibars; - used for atmospheric correction - temp : int or float - avg. yearly temperature at location in - degrees C; used for atmospheric correction - delta_t : float or array - Difference between terrestrial time and UT1. - atmos_refract : float - The approximate atmospheric refraction (in degrees) - at sunrise and sunset. - numthreads: int, optional, default 8 - Number of threads to use for computation if numba>=0.17 - is installed. - sst : bool, default False - If True, return only data needed for sunrise, sunset, and transit - calculations. - esd : bool, default False - If True, return only Earth-Sun distance in AU - - Returns - ------- - Numpy Array with elements: - apparent zenith, - zenith, - apparent_elevation, - elevation, - azimuth, - equation_of_time - - References - ---------- - [1] I. Reda and A. Andreas, Solar position algorithm for solar radiation - applications. Solar Energy, vol. 76, no. 5, pp. 577-589, 2004. - - [2] I. Reda and A. Andreas, Corrigendum to Solar position algorithm for - solar radiation applications. Solar Energy, vol. 81, no. 6, p. 838, 2007. - """ - if USE_NUMBA: - do_calc = solar_position_numba - else: - do_calc = solar_position_numpy - - result = do_calc(unixtime, lat, lon, elev, pressure, - temp, delta_t, atmos_refract, numthreads, - sst, esd) - - if not isinstance(result, np.ndarray): - try: - result = np.array(result) - except Exception: - pass - - return result - - -def transit_sunrise_sunset(dates, lat, lon, delta_t, numthreads): - """ - Calculate the sun transit, sunrise, and sunset - for a set of dates at a given location using the NREL SPA. - - The NREL SPA is documented in [1]_, [2]_. - - Parameters - ---------- - dates : array - Numpy array of ints/floats corresponding to the Unix time - for the dates of interest, must be midnight UTC (00:00+00:00) - on the day of interest. - lat : float - Latitude of location to perform calculation for - lon : float - Longitude of location - delta_t : float or array - Difference between terrestrial time and UT. USNO has tables. - numthreads : int - Number to threads to use for calculation (if using numba) - - Returns - ------- - tuple : (transit, sunrise, sunset) localized to UTC - +""" +Calculate the sun's position using the NREL Solar Position Algorithm +(SPA) either using numpy arrays or compiling the code to machine language +with numba. +""" + +# Contributors: +# Created by Tony Lorenzo (@alorenzo175), Univ. of Arizona, 2015 + +import os +import threading +import warnings + +import numpy as np + + +# this block is a way to use an environment variable to switch between +# compiling the functions with numba or just use numpy +def nocompile(*args, **kwargs): + return lambda func: func + + +if os.getenv('PVLIB_USE_NUMBA', '0') != '0': + try: + from numba import jit + except ImportError: + warnings.warn('Could not import numba, falling back to numpy ' + + 'calculation') + jcompile = nocompile + USE_NUMBA = False + else: + jcompile = jit + USE_NUMBA = True +else: + jcompile = nocompile + USE_NUMBA = False + + +# heliocentric longitude coefficients +L0 = np.array([ + [175347046.0, 0.0, 0.0], + [3341656.0, 4.6692568, 6283.07585], + [34894.0, 4.6261, 12566.1517], + [3497.0, 2.7441, 5753.3849], + [3418.0, 2.8289, 3.5231], + [3136.0, 3.6277, 77713.7715], + [2676.0, 4.4181, 7860.4194], + [2343.0, 6.1352, 3930.2097], + [1324.0, 0.7425, 11506.7698], + [1273.0, 2.0371, 529.691], + [1199.0, 1.1096, 1577.3435], + [990.0, 5.233, 5884.927], + [902.0, 2.045, 26.298], + [857.0, 3.508, 398.149], + [780.0, 1.179, 5223.694], + [753.0, 2.533, 5507.553], + [505.0, 4.583, 18849.228], + [492.0, 4.205, 775.523], + [357.0, 2.92, 0.067], + [317.0, 5.849, 11790.629], + [284.0, 1.899, 796.298], + [271.0, 0.315, 10977.079], + [243.0, 0.345, 5486.778], + [206.0, 4.806, 2544.314], + [205.0, 1.869, 5573.143], + [202.0, 2.458, 6069.777], + [156.0, 0.833, 213.299], + [132.0, 3.411, 2942.463], + [126.0, 1.083, 20.775], + [115.0, 0.645, 0.98], + [103.0, 0.636, 4694.003], + [102.0, 0.976, 15720.839], + [102.0, 4.267, 7.114], + [99.0, 6.21, 2146.17], + [98.0, 0.68, 155.42], + [86.0, 5.98, 161000.69], + [85.0, 1.3, 6275.96], + [85.0, 3.67, 71430.7], + [80.0, 1.81, 17260.15], + [79.0, 3.04, 12036.46], + [75.0, 1.76, 5088.63], + [74.0, 3.5, 3154.69], + [74.0, 4.68, 801.82], + [70.0, 0.83, 9437.76], + [62.0, 3.98, 8827.39], + [61.0, 1.82, 7084.9], + [57.0, 2.78, 6286.6], + [56.0, 4.39, 14143.5], + [56.0, 3.47, 6279.55], + [52.0, 0.19, 12139.55], + [52.0, 1.33, 1748.02], + [51.0, 0.28, 5856.48], + [49.0, 0.49, 1194.45], + [41.0, 5.37, 8429.24], + [41.0, 2.4, 19651.05], + [39.0, 6.17, 10447.39], + [37.0, 6.04, 10213.29], + [37.0, 2.57, 1059.38], + [36.0, 1.71, 2352.87], + [36.0, 1.78, 6812.77], + [33.0, 0.59, 17789.85], + [30.0, 0.44, 83996.85], + [30.0, 2.74, 1349.87], + [25.0, 3.16, 4690.48] +]) +L1 = np.array([ + [628331966747.0, 0.0, 0.0], + [206059.0, 2.678235, 6283.07585], + [4303.0, 2.6351, 12566.1517], + [425.0, 1.59, 3.523], + [119.0, 5.796, 26.298], + [109.0, 2.966, 1577.344], + [93.0, 2.59, 18849.23], + [72.0, 1.14, 529.69], + [68.0, 1.87, 398.15], + [67.0, 4.41, 5507.55], + [59.0, 2.89, 5223.69], + [56.0, 2.17, 155.42], + [45.0, 0.4, 796.3], + [36.0, 0.47, 775.52], + [29.0, 2.65, 7.11], + [21.0, 5.34, 0.98], + [19.0, 1.85, 5486.78], + [19.0, 4.97, 213.3], + [17.0, 2.99, 6275.96], + [16.0, 0.03, 2544.31], + [16.0, 1.43, 2146.17], + [15.0, 1.21, 10977.08], + [12.0, 2.83, 1748.02], + [12.0, 3.26, 5088.63], + [12.0, 5.27, 1194.45], + [12.0, 2.08, 4694.0], + [11.0, 0.77, 553.57], + [10.0, 1.3, 6286.6], + [10.0, 4.24, 1349.87], + [9.0, 2.7, 242.73], + [9.0, 5.64, 951.72], + [8.0, 5.3, 2352.87], + [6.0, 2.65, 9437.76], + [6.0, 4.67, 4690.48] +]) +L2 = np.array([ + [52919.0, 0.0, 0.0], + [8720.0, 1.0721, 6283.0758], + [309.0, 0.867, 12566.152], + [27.0, 0.05, 3.52], + [16.0, 5.19, 26.3], + [16.0, 3.68, 155.42], + [10.0, 0.76, 18849.23], + [9.0, 2.06, 77713.77], + [7.0, 0.83, 775.52], + [5.0, 4.66, 1577.34], + [4.0, 1.03, 7.11], + [4.0, 3.44, 5573.14], + [3.0, 5.14, 796.3], + [3.0, 6.05, 5507.55], + [3.0, 1.19, 242.73], + [3.0, 6.12, 529.69], + [3.0, 0.31, 398.15], + [3.0, 2.28, 553.57], + [2.0, 4.38, 5223.69], + [2.0, 3.75, 0.98] +]) +L3 = np.array([ + [289.0, 5.844, 6283.076], + [35.0, 0.0, 0.0], + [17.0, 5.49, 12566.15], + [3.0, 5.2, 155.42], + [1.0, 4.72, 3.52], + [1.0, 5.3, 18849.23], + [1.0, 5.97, 242.73] +]) +L4 = np.array([ + [114.0, 3.142, 0.0], + [8.0, 4.13, 6283.08], + [1.0, 3.84, 12566.15] +]) +L5 = np.array([ + [1.0, 3.14, 0.0] +]) + + +# heliocentric latitude coefficients +B0 = np.array([ + [280.0, 3.199, 84334.662], + [102.0, 5.422, 5507.553], + [80.0, 3.88, 5223.69], + [44.0, 3.7, 2352.87], + [32.0, 4.0, 1577.34] +]) +B1 = np.array([ + [9.0, 3.9, 5507.55], + [6.0, 1.73, 5223.69] +]) + + +# heliocentric radius coefficients +R0 = np.array([ + [100013989.0, 0.0, 0.0], + [1670700.0, 3.0984635, 6283.07585], + [13956.0, 3.05525, 12566.1517], + [3084.0, 5.1985, 77713.7715], + [1628.0, 1.1739, 5753.3849], + [1576.0, 2.8469, 7860.4194], + [925.0, 5.453, 11506.77], + [542.0, 4.564, 3930.21], + [472.0, 3.661, 5884.927], + [346.0, 0.964, 5507.553], + [329.0, 5.9, 5223.694], + [307.0, 0.299, 5573.143], + [243.0, 4.273, 11790.629], + [212.0, 5.847, 1577.344], + [186.0, 5.022, 10977.079], + [175.0, 3.012, 18849.228], + [110.0, 5.055, 5486.778], + [98.0, 0.89, 6069.78], + [86.0, 5.69, 15720.84], + [86.0, 1.27, 161000.69], + [65.0, 0.27, 17260.15], + [63.0, 0.92, 529.69], + [57.0, 2.01, 83996.85], + [56.0, 5.24, 71430.7], + [49.0, 3.25, 2544.31], + [47.0, 2.58, 775.52], + [45.0, 5.54, 9437.76], + [43.0, 6.01, 6275.96], + [39.0, 5.36, 4694.0], + [38.0, 2.39, 8827.39], + [37.0, 0.83, 19651.05], + [37.0, 4.9, 12139.55], + [36.0, 1.67, 12036.46], + [35.0, 1.84, 2942.46], + [33.0, 0.24, 7084.9], + [32.0, 0.18, 5088.63], + [32.0, 1.78, 398.15], + [28.0, 1.21, 6286.6], + [28.0, 1.9, 6279.55], + [26.0, 4.59, 10447.39] +]) +R1 = np.array([ + [103019.0, 1.10749, 6283.07585], + [1721.0, 1.0644, 12566.1517], + [702.0, 3.142, 0.0], + [32.0, 1.02, 18849.23], + [31.0, 2.84, 5507.55], + [25.0, 1.32, 5223.69], + [18.0, 1.42, 1577.34], + [10.0, 5.91, 10977.08], + [9.0, 1.42, 6275.96], + [9.0, 0.27, 5486.78] +]) +R2 = np.array([ + [4359.0, 5.7846, 6283.0758], + [124.0, 5.579, 12566.152], + [12.0, 3.14, 0.0], + [9.0, 3.63, 77713.77], + [6.0, 1.87, 5573.14], + [3.0, 5.47, 18849.23] +]) +R3 = np.array([ + [145.0, 4.273, 6283.076], + [7.0, 3.92, 12566.15] +]) +R4 = np.array([ + [4.0, 2.56, 6283.08] +]) + + +# longitude and obliquity nutation coefficients +NUTATION_ABCD_ARRAY = np.array([ + [-171996, -174.2, 92025, 8.9], + [-13187, -1.6, 5736, -3.1], + [-2274, -0.2, 977, -0.5], + [2062, 0.2, -895, 0.5], + [1426, -3.4, 54, -0.1], + [712, 0.1, -7, 0], + [-517, 1.2, 224, -0.6], + [-386, -0.4, 200, 0], + [-301, 0, 129, -0.1], + [217, -0.5, -95, 0.3], + [-158, 0, 0, 0], + [129, 0.1, -70, 0], + [123, 0, -53, 0], + [63, 0, 0, 0], + [63, 0.1, -33, 0], + [-59, 0, 26, 0], + [-58, -0.1, 32, 0], + [-51, 0, 27, 0], + [48, 0, 0, 0], + [46, 0, -24, 0], + [-38, 0, 16, 0], + [-31, 0, 13, 0], + [29, 0, 0, 0], + [29, 0, -12, 0], + [26, 0, 0, 0], + [-22, 0, 0, 0], + [21, 0, -10, 0], + [17, -0.1, 0, 0], + [16, 0, -8, 0], + [-16, 0.1, 7, 0], + [-15, 0, 9, 0], + [-13, 0, 7, 0], + [-12, 0, 6, 0], + [11, 0, 0, 0], + [-10, 0, 5, 0], + [-8, 0, 3, 0], + [7, 0, -3, 0], + [-7, 0, 0, 0], + [-7, 0, 3, 0], + [-7, 0, 3, 0], + [6, 0, 0, 0], + [6, 0, -3, 0], + [6, 0, -3, 0], + [-6, 0, 3, 0], + [-6, 0, 3, 0], + [5, 0, 0, 0], + [-5, 0, 3, 0], + [-5, 0, 3, 0], + [-5, 0, 3, 0], + [4, 0, 0, 0], + [4, 0, 0, 0], + [4, 0, 0, 0], + [-4, 0, 0, 0], + [-4, 0, 0, 0], + [-4, 0, 0, 0], + [3, 0, 0, 0], + [-3, 0, 0, 0], + [-3, 0, 0, 0], + [-3, 0, 0, 0], + [-3, 0, 0, 0], + [-3, 0, 0, 0], + [-3, 0, 0, 0], + [-3, 0, 0, 0], +]) + +NUTATION_YTERM_ARRAY = np.array([ + [0, 0, 0, 0, 1], + [-2, 0, 0, 2, 2], + [0, 0, 0, 2, 2], + [0, 0, 0, 0, 2], + [0, 1, 0, 0, 0], + [0, 0, 1, 0, 0], + [-2, 1, 0, 2, 2], + [0, 0, 0, 2, 1], + [0, 0, 1, 2, 2], + [-2, -1, 0, 2, 2], + [-2, 0, 1, 0, 0], + [-2, 0, 0, 2, 1], + [0, 0, -1, 2, 2], + [2, 0, 0, 0, 0], + [0, 0, 1, 0, 1], + [2, 0, -1, 2, 2], + [0, 0, -1, 0, 1], + [0, 0, 1, 2, 1], + [-2, 0, 2, 0, 0], + [0, 0, -2, 2, 1], + [2, 0, 0, 2, 2], + [0, 0, 2, 2, 2], + [0, 0, 2, 0, 0], + [-2, 0, 1, 2, 2], + [0, 0, 0, 2, 0], + [-2, 0, 0, 2, 0], + [0, 0, -1, 2, 1], + [0, 2, 0, 0, 0], + [2, 0, -1, 0, 1], + [-2, 2, 0, 2, 2], + [0, 1, 0, 0, 1], + [-2, 0, 1, 0, 1], + [0, -1, 0, 0, 1], + [0, 0, 2, -2, 0], + [2, 0, -1, 2, 1], + [2, 0, 1, 2, 2], + [0, 1, 0, 2, 2], + [-2, 1, 1, 0, 0], + [0, -1, 0, 2, 2], + [2, 0, 0, 2, 1], + [2, 0, 1, 0, 0], + [-2, 0, 2, 2, 2], + [-2, 0, 1, 2, 1], + [2, 0, -2, 0, 1], + [2, 0, 0, 0, 1], + [0, -1, 1, 0, 0], + [-2, -1, 0, 2, 1], + [-2, 0, 0, 0, 1], + [0, 0, 2, 2, 1], + [-2, 0, 2, 0, 1], + [-2, 1, 0, 2, 1], + [0, 0, 1, -2, 0], + [-1, 0, 1, 0, 0], + [-2, 1, 0, 0, 0], + [1, 0, 0, 0, 0], + [0, 0, 1, 2, 0], + [0, 0, -2, 2, 2], + [-1, -1, 1, 0, 0], + [0, 1, 1, 0, 0], + [0, -1, 1, 2, 2], + [2, -1, -1, 2, 2], + [0, 0, 3, 2, 2], + [2, -1, 0, 2, 2], +]) + + +@jcompile('float64(int64, int64, int64, int64, int64, int64, int64)', + nopython=True) +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 + + +@jcompile('float64(float64)', nopython=True) +def julian_day(unixtime): + jd = unixtime * 1.0 / 86400 + 2440587.5 + return jd + + +@jcompile('float64(float64, float64)', nopython=True) +def julian_ephemeris_day(julian_day, delta_t): + jde = julian_day + delta_t * 1.0 / 86400 + return jde + + +@jcompile('float64(float64)', nopython=True) +def julian_century(julian_day): + jc = (julian_day - 2451545) * 1.0 / 36525 + return jc + + +@jcompile('float64(float64)', nopython=True) +def julian_ephemeris_century(julian_ephemeris_day): + jce = (julian_ephemeris_day - 2451545) * 1.0 / 36525 + return jce + + +@jcompile('float64(float64)', nopython=True) +def julian_ephemeris_millennium(julian_ephemeris_century): + jme = julian_ephemeris_century * 1.0 / 10 + return jme + + +# omit type signature here; specifying read-only arrays requires use of the +# numba.types API, meaning numba must be available to import. +# https://github.com/numba/numba/issues/4511 +@jcompile(nopython=True) +def sum_mult_cos_add_mult(arr, x): + # shared calculation used for heliocentric longitude, latitude, and radius + s = 0. + for row in range(arr.shape[0]): + s += arr[row, 0] * np.cos(arr[row, 1] + arr[row, 2] * x) + return s + +@jcompile('float64(float64)', nopython=True) +def heliocentric_longitude(jme): + l0 = sum_mult_cos_add_mult(L0, jme) + l1 = sum_mult_cos_add_mult(L1, jme) + l2 = sum_mult_cos_add_mult(L2, jme) + l3 = sum_mult_cos_add_mult(L3, jme) + l4 = sum_mult_cos_add_mult(L4, jme) + l5 = sum_mult_cos_add_mult(L5, jme) + + l_rad = (l0 + l1 * jme + l2 * jme**2 + l3 * jme**3 + l4 * jme**4 + + l5 * jme**5)/10**8 + l = np.rad2deg(l_rad) + return l % 360 + +@jcompile('float64(float64)', nopython=True) +def heliocentric_latitude(jme): + b0 = sum_mult_cos_add_mult(B0, jme) + b1 = sum_mult_cos_add_mult(B1, jme) + + b_rad = (b0 + b1 * jme)/10**8 + b = np.rad2deg(b_rad) + return b + + +@jcompile('float64(float64)', nopython=True) +def heliocentric_radius_vector(jme): + r0 = sum_mult_cos_add_mult(R0, jme) + r1 = sum_mult_cos_add_mult(R1, jme) + r2 = sum_mult_cos_add_mult(R2, jme) + r3 = sum_mult_cos_add_mult(R3, jme) + r4 = sum_mult_cos_add_mult(R4, jme) + + r = (r0 + r1 * jme + r2 * jme**2 + r3 * jme**3 + r4 * jme**4)/10**8 + return r + + +@jcompile('float64(float64)', nopython=True) +def geocentric_longitude(heliocentric_longitude): + theta = heliocentric_longitude + 180.0 + return theta % 360 + + +@jcompile('float64(float64)', nopython=True) +def geocentric_latitude(heliocentric_latitude): + beta = -1.0*heliocentric_latitude + return beta + + +@jcompile('float64(float64)', nopython=True) +def mean_elongation(julian_ephemeris_century): + x0 = (297.85036 + + 445267.111480 * julian_ephemeris_century + - 0.0019142 * julian_ephemeris_century**2 + + julian_ephemeris_century**3 / 189474) + return x0 + + +@jcompile('float64(float64)', nopython=True) +def mean_anomaly_sun(julian_ephemeris_century): + x1 = (357.52772 + + 35999.050340 * julian_ephemeris_century + - 0.0001603 * julian_ephemeris_century**2 + - julian_ephemeris_century**3 / 300000) + return x1 + + +@jcompile('float64(float64)', nopython=True) +def mean_anomaly_moon(julian_ephemeris_century): + x2 = (134.96298 + + 477198.867398 * julian_ephemeris_century + + 0.0086972 * julian_ephemeris_century**2 + + julian_ephemeris_century**3 / 56250) + return x2 + + +@jcompile('float64(float64)', nopython=True) +def moon_argument_latitude(julian_ephemeris_century): + x3 = (93.27191 + + 483202.017538 * julian_ephemeris_century + - 0.0036825 * julian_ephemeris_century**2 + + julian_ephemeris_century**3 / 327270) + return x3 + + +@jcompile('float64(float64)', nopython=True) +def moon_ascending_longitude(julian_ephemeris_century): + x4 = (125.04452 + - 1934.136261 * julian_ephemeris_century + + 0.0020708 * julian_ephemeris_century**2 + + julian_ephemeris_century**3 / 450000) + return x4 + + +@jcompile( + 'void(float64, float64, float64, float64, float64, float64, float64[:])', + nopython=True) +def longitude_obliquity_nutation(julian_ephemeris_century, x0, x1, x2, x3, x4, + out): + delta_psi_sum = 0.0 + delta_eps_sum = 0.0 + for row in range(NUTATION_YTERM_ARRAY.shape[0]): + a = NUTATION_ABCD_ARRAY[row, 0] + b = NUTATION_ABCD_ARRAY[row, 1] + c = NUTATION_ABCD_ARRAY[row, 2] + d = NUTATION_ABCD_ARRAY[row, 3] + arg = np.radians( + NUTATION_YTERM_ARRAY[row, 0]*x0 + + NUTATION_YTERM_ARRAY[row, 1]*x1 + + NUTATION_YTERM_ARRAY[row, 2]*x2 + + NUTATION_YTERM_ARRAY[row, 3]*x3 + + NUTATION_YTERM_ARRAY[row, 4]*x4 + ) + delta_psi_sum += (a + b * julian_ephemeris_century) * np.sin(arg) + delta_eps_sum += (c + d * julian_ephemeris_century) * np.cos(arg) + delta_psi = delta_psi_sum*1.0/36000000 + delta_eps = delta_eps_sum*1.0/36000000 + # seems like we ought to be able to return a tuple here instead + # of resorting to `out`, but returning a UniTuple from this + # function caused calculations elsewhere to give the wrong result. + # very difficult to investigate since it did not occur when using + # object mode. issue was observed on numba 0.56.4 + out[0] = delta_psi + out[1] = delta_eps + + +@jcompile('float64(float64)', nopython=True) +def mean_ecliptic_obliquity(julian_ephemeris_millennium): + U = 1.0*julian_ephemeris_millennium/10 + e0 = (84381.448 - 4680.93 * U - 1.55 * U**2 + + 1999.25 * U**3 - 51.38 * U**4 - 249.67 * U**5 + - 39.05 * U**6 + 7.12 * U**7 + 27.87 * U**8 + + 5.79 * U**9 + 2.45 * U**10) + return e0 + + +@jcompile('float64(float64, float64)', nopython=True) +def true_ecliptic_obliquity(mean_ecliptic_obliquity, obliquity_nutation): + e0 = mean_ecliptic_obliquity + deleps = obliquity_nutation + e = e0*1.0/3600 + deleps + return e + + +@jcompile('float64(float64)', nopython=True) +def aberration_correction(earth_radius_vector): + deltau = -20.4898 / (3600 * earth_radius_vector) + return deltau + + +@jcompile('float64(float64, float64, float64)', nopython=True) +def apparent_sun_longitude(geocentric_longitude, longitude_nutation, + aberration_correction): + lamd = geocentric_longitude + longitude_nutation + aberration_correction + return lamd + + +@jcompile('float64(float64, float64)', nopython=True) +def mean_sidereal_time(julian_day, julian_century): + v0 = (280.46061837 + 360.98564736629 * (julian_day - 2451545) + + 0.000387933 * julian_century**2 - julian_century**3 / 38710000) + return v0 % 360.0 + + +@jcompile('float64(float64, float64, float64)', nopython=True) +def apparent_sidereal_time(mean_sidereal_time, longitude_nutation, + true_ecliptic_obliquity): + v = mean_sidereal_time + longitude_nutation * np.cos( + np.radians(true_ecliptic_obliquity)) + return v + + +@jcompile('float64(float64, float64, float64)', nopython=True) +def geocentric_sun_right_ascension(apparent_sun_longitude, + true_ecliptic_obliquity, + geocentric_latitude): + true_ecliptic_obliquity_rad = np.radians(true_ecliptic_obliquity) + apparent_sun_longitude_rad = np.radians(apparent_sun_longitude) + + num = (np.sin(apparent_sun_longitude_rad) + * np.cos(true_ecliptic_obliquity_rad) + - np.tan(np.radians(geocentric_latitude)) + * np.sin(true_ecliptic_obliquity_rad)) + alpha = np.degrees(np.arctan2(num, np.cos(apparent_sun_longitude_rad))) + return alpha % 360 + + +@jcompile('float64(float64, float64, float64)', nopython=True) +def geocentric_sun_declination(apparent_sun_longitude, true_ecliptic_obliquity, + geocentric_latitude): + geocentric_latitude_rad = np.radians(geocentric_latitude) + true_ecliptic_obliquity_rad = np.radians(true_ecliptic_obliquity) + + delta = np.degrees(np.arcsin(np.sin(geocentric_latitude_rad) * + np.cos(true_ecliptic_obliquity_rad) + + np.cos(geocentric_latitude_rad) * + np.sin(true_ecliptic_obliquity_rad) * + np.sin(np.radians(apparent_sun_longitude)))) + return delta + + +@jcompile('float64(float64, float64, float64)', nopython=True) +def local_hour_angle(apparent_sidereal_time, observer_longitude, + sun_right_ascension): + """Measured westward from south""" + H = apparent_sidereal_time + observer_longitude - sun_right_ascension + return H % 360 + + +@jcompile('float64(float64)', nopython=True) +def equatorial_horizontal_parallax(earth_radius_vector): + xi = 8.794 / (3600 * earth_radius_vector) + return xi + + +@jcompile('float64(float64)', nopython=True) +def uterm(observer_latitude): + u = np.arctan(0.99664719 * np.tan(np.radians(observer_latitude))) + return u + + +@jcompile('float64(float64, float64, float64)', nopython=True) +def xterm(u, observer_latitude, observer_elevation): + x = (np.cos(u) + observer_elevation / 6378140 + * np.cos(np.radians(observer_latitude))) + return x + + +@jcompile('float64(float64, float64, float64)', nopython=True) +def yterm(u, observer_latitude, observer_elevation): + y = (0.99664719 * np.sin(u) + observer_elevation / 6378140 + * np.sin(np.radians(observer_latitude))) + return y + + +@jcompile('float64(float64, float64,float64, float64)', nopython=True) +def parallax_sun_right_ascension(xterm, equatorial_horizontal_parallax, + local_hour_angle, geocentric_sun_declination): + equatorial_horizontal_parallax_rad = \ + np.radians(equatorial_horizontal_parallax) + local_hour_angle_rad = np.radians(local_hour_angle) + + num = (-xterm * np.sin(equatorial_horizontal_parallax_rad) + * np.sin(local_hour_angle_rad)) + denom = (np.cos(np.radians(geocentric_sun_declination)) + - xterm * np.sin(equatorial_horizontal_parallax_rad) + * np.cos(local_hour_angle_rad)) + delta_alpha = np.degrees(np.arctan2(num, denom)) + return delta_alpha + + +@jcompile('float64(float64, float64)', nopython=True) +def topocentric_sun_right_ascension(geocentric_sun_right_ascension, + parallax_sun_right_ascension): + alpha_prime = geocentric_sun_right_ascension + parallax_sun_right_ascension + return alpha_prime + + +@jcompile('float64(float64, float64, float64, float64, float64, float64)', + nopython=True) +def topocentric_sun_declination(geocentric_sun_declination, xterm, yterm, + equatorial_horizontal_parallax, + parallax_sun_right_ascension, + local_hour_angle): + geocentric_sun_declination_rad = np.radians(geocentric_sun_declination) + equatorial_horizontal_parallax_rad = \ + np.radians(equatorial_horizontal_parallax) + + num = ((np.sin(geocentric_sun_declination_rad) - yterm + * np.sin(equatorial_horizontal_parallax_rad)) + * np.cos(np.radians(parallax_sun_right_ascension))) + denom = (np.cos(geocentric_sun_declination_rad) - xterm + * np.sin(equatorial_horizontal_parallax_rad) + * np.cos(np.radians(local_hour_angle))) + delta = np.degrees(np.arctan2(num, denom)) + return delta + + +@jcompile('float64(float64, float64)', nopython=True) +def topocentric_local_hour_angle(local_hour_angle, + parallax_sun_right_ascension): + H_prime = local_hour_angle - parallax_sun_right_ascension + return H_prime + + +@jcompile('float64(float64, float64, float64)', nopython=True) +def topocentric_elevation_angle_without_atmosphere(observer_latitude, + topocentric_sun_declination, + topocentric_local_hour_angle + ): + observer_latitude_rad = np.radians(observer_latitude) + topocentric_sun_declination_rad = np.radians(topocentric_sun_declination) + e0 = np.degrees(np.arcsin( + np.sin(observer_latitude_rad) + * np.sin(topocentric_sun_declination_rad) + + np.cos(observer_latitude_rad) + * np.cos(topocentric_sun_declination_rad) + * np.cos(np.radians(topocentric_local_hour_angle)))) + return e0 + + +@jcompile('float64(float64, float64, float64, float64)', nopython=True) +def atmospheric_refraction_correction(local_pressure, local_temp, + topocentric_elevation_angle_wo_atmosphere, + atmos_refract): + # switch sets delta_e when the sun is below the horizon + switch = topocentric_elevation_angle_wo_atmosphere >= -1.0 * ( + 0.26667 + atmos_refract) + delta_e = ((local_pressure / 1010.0) * (283.0 / (273 + local_temp)) + * 1.02 / (60 * np.tan(np.radians( + topocentric_elevation_angle_wo_atmosphere + + 10.3 / (topocentric_elevation_angle_wo_atmosphere + + 5.11))))) * switch + return delta_e + + +@jcompile('float64(float64, float64)', nopython=True) +def topocentric_elevation_angle(topocentric_elevation_angle_without_atmosphere, + atmospheric_refraction_correction): + e = (topocentric_elevation_angle_without_atmosphere + + atmospheric_refraction_correction) + return e + + +@jcompile('float64(float64)', nopython=True) +def topocentric_zenith_angle(topocentric_elevation_angle): + theta = 90 - topocentric_elevation_angle + return theta + + +@jcompile('float64(float64, float64, float64)', nopython=True) +def topocentric_astronomers_azimuth(topocentric_local_hour_angle, + topocentric_sun_declination, + observer_latitude): + topocentric_local_hour_angle_rad = np.radians(topocentric_local_hour_angle) + observer_latitude_rad = np.radians(observer_latitude) + num = np.sin(topocentric_local_hour_angle_rad) + denom = (np.cos(topocentric_local_hour_angle_rad) + * np.sin(observer_latitude_rad) + - np.tan(np.radians(topocentric_sun_declination)) + * np.cos(observer_latitude_rad)) + gamma = np.degrees(np.arctan2(num, denom)) + return gamma % 360 + + +@jcompile('float64(float64)', nopython=True) +def topocentric_azimuth_angle(topocentric_astronomers_azimuth): + phi = topocentric_astronomers_azimuth + 180 + return phi % 360 + + +@jcompile('float64(float64)', nopython=True) +def sun_mean_longitude(julian_ephemeris_millennium): + M = (280.4664567 + 360007.6982779 * julian_ephemeris_millennium + + 0.03032028 * julian_ephemeris_millennium**2 + + julian_ephemeris_millennium**3 / 49931 + - julian_ephemeris_millennium**4 / 15300 + - julian_ephemeris_millennium**5 / 2000000) + return M + + +@jcompile('float64(float64, float64, float64, float64)', nopython=True) +def equation_of_time(sun_mean_longitude, geocentric_sun_right_ascension, + longitude_nutation, true_ecliptic_obliquity): + E = (sun_mean_longitude - 0.0057183 - geocentric_sun_right_ascension + + longitude_nutation * np.cos(np.radians(true_ecliptic_obliquity))) + # limit between 0 and 360 + E = E % 360 + # convert to minutes + E *= 4 + greater = E > 20 + less = E < -20 + other = (E <= 20) & (E >= -20) + E = greater * (E - 1440) + less * (E + 1440) + other * E + return E + + +@jcompile('void(float64[:], float64[:], float64[:], float64[:,:])', + nopython=True, nogil=True) +def solar_position_loop(unixtime, delta_t, loc_args, out): + """Loop through the time array and calculate the solar position""" + lat = loc_args[0] + lon = loc_args[1] + elev = loc_args[2] + pressure = loc_args[3] + temp = loc_args[4] + atmos_refract = loc_args[5] + sst = loc_args[6] + esd = loc_args[7] + + for i in range(unixtime.shape[0]): + utime = unixtime[i] + dT = delta_t[i] + jd = julian_day(utime) + jde = julian_ephemeris_day(jd, dT) + jc = julian_century(jd) + jce = julian_ephemeris_century(jde) + jme = julian_ephemeris_millennium(jce) + R = heliocentric_radius_vector(jme) + if esd: + out[0, i] = R + continue + L = heliocentric_longitude(jme) + B = heliocentric_latitude(jme) + Theta = geocentric_longitude(L) + beta = geocentric_latitude(B) + x0 = mean_elongation(jce) + x1 = mean_anomaly_sun(jce) + x2 = mean_anomaly_moon(jce) + x3 = moon_argument_latitude(jce) + x4 = moon_ascending_longitude(jce) + l_o_nutation = np.empty((2,)) + longitude_obliquity_nutation(jce, x0, x1, x2, x3, x4, l_o_nutation) + delta_psi = l_o_nutation[0] + delta_epsilon = l_o_nutation[1] + epsilon0 = mean_ecliptic_obliquity(jme) + epsilon = true_ecliptic_obliquity(epsilon0, delta_epsilon) + delta_tau = aberration_correction(R) + lamd = apparent_sun_longitude(Theta, delta_psi, delta_tau) + v0 = mean_sidereal_time(jd, jc) + v = apparent_sidereal_time(v0, delta_psi, epsilon) + alpha = geocentric_sun_right_ascension(lamd, epsilon, beta) + delta = geocentric_sun_declination(lamd, epsilon, beta) + if sst: + out[0, i] = v + out[1, i] = alpha + out[2, i] = delta + continue + m = sun_mean_longitude(jme) + eot = equation_of_time(m, alpha, delta_psi, epsilon) + H = local_hour_angle(v, lon, alpha) + xi = equatorial_horizontal_parallax(R) + u = uterm(lat) + x = xterm(u, lat, elev) + y = yterm(u, lat, elev) + delta_alpha = parallax_sun_right_ascension(x, xi, H, delta) + delta_prime = topocentric_sun_declination(delta, x, y, xi, delta_alpha, + H) + H_prime = topocentric_local_hour_angle(H, delta_alpha) + e0 = topocentric_elevation_angle_without_atmosphere(lat, delta_prime, + H_prime) + delta_e = atmospheric_refraction_correction(pressure, temp, e0, + atmos_refract) + e = topocentric_elevation_angle(e0, delta_e) + theta = topocentric_zenith_angle(e) + theta0 = topocentric_zenith_angle(e0) + gamma = topocentric_astronomers_azimuth(H_prime, delta_prime, lat) + phi = topocentric_azimuth_angle(gamma) + out[0, i] = theta + out[1, i] = theta0 + out[2, i] = e + out[3, i] = e0 + out[4, i] = phi + out[5, i] = eot + + +def solar_position_numba(unixtime, lat, lon, elev, pressure, temp, delta_t, + atmos_refract, numthreads, sst=False, esd=False): + """Calculate the solar position using the numba compiled functions + and multiple threads. Very slow if functions are not numba compiled. + """ + # these args are the same for each thread + loc_args = np.array([lat, lon, elev, pressure, temp, + atmos_refract, sst, esd], dtype=np.float64) + + # turn delta_t into an array if it isn't already + delta_t = np.full_like(unixtime, delta_t, dtype=np.float64) + + # construct dims x ulength array to put the results in + ulength = unixtime.shape[0] + if sst: + dims = 3 + elif esd: + dims = 1 + else: + dims = 6 + result = np.empty((dims, ulength), dtype=np.float64) + + if unixtime.dtype != np.float64: + unixtime = unixtime.astype(np.float64) + + if ulength < numthreads: + numthreads = ulength + + if numthreads <= 1: + solar_position_loop(unixtime, delta_t, loc_args, result) + return result + + # split the input and output arrays into numthreads chunks + split0 = np.array_split(unixtime, numthreads) + split1 = np.array_split(delta_t, numthreads) + split2 = np.array_split(result, numthreads, axis=1) + chunks = [ + [a0, a1, loc_args, a2] + for a0, a1, a2 in zip(split0, split1, split2) + ] + # Spawn one thread per chunk + threads = [threading.Thread(target=solar_position_loop, args=chunk) + for chunk in chunks] + for thread in threads: + thread.start() + for thread in threads: + thread.join() + return result + + +def solar_position_numpy(unixtime, lat, lon, elev, pressure, temp, delta_t, + atmos_refract, numthreads, sst=False, esd=False): + """Calculate the solar position assuming unixtime is a numpy array. Note + this function will not work if the solar position functions were + compiled with numba. + """ + + jd = julian_day(unixtime) + jde = julian_ephemeris_day(jd, delta_t) + jc = julian_century(jd) + jce = julian_ephemeris_century(jde) + jme = julian_ephemeris_millennium(jce) + R = heliocentric_radius_vector(jme) + if esd: + return (R, ) + L = heliocentric_longitude(jme) + B = heliocentric_latitude(jme) + Theta = geocentric_longitude(L) + beta = geocentric_latitude(B) + x0 = mean_elongation(jce) + x1 = mean_anomaly_sun(jce) + x2 = mean_anomaly_moon(jce) + x3 = moon_argument_latitude(jce) + x4 = moon_ascending_longitude(jce) + l_o_nutation = np.empty((2, len(x0))) + longitude_obliquity_nutation(jce, x0, x1, x2, x3, x4, l_o_nutation) + delta_psi = l_o_nutation[0] + delta_epsilon = l_o_nutation[1] + epsilon0 = mean_ecliptic_obliquity(jme) + epsilon = true_ecliptic_obliquity(epsilon0, delta_epsilon) + delta_tau = aberration_correction(R) + lamd = apparent_sun_longitude(Theta, delta_psi, delta_tau) + v0 = mean_sidereal_time(jd, jc) + v = apparent_sidereal_time(v0, delta_psi, epsilon) + alpha = geocentric_sun_right_ascension(lamd, epsilon, beta) + delta = geocentric_sun_declination(lamd, epsilon, beta) + if sst: + return v, alpha, delta + m = sun_mean_longitude(jme) + eot = equation_of_time(m, alpha, delta_psi, epsilon) + H = local_hour_angle(v, lon, alpha) + xi = equatorial_horizontal_parallax(R) + u = uterm(lat) + x = xterm(u, lat, elev) + y = yterm(u, lat, elev) + delta_alpha = parallax_sun_right_ascension(x, xi, H, delta) + delta_prime = topocentric_sun_declination(delta, x, y, xi, delta_alpha, H) + H_prime = topocentric_local_hour_angle(H, delta_alpha) + e0 = topocentric_elevation_angle_without_atmosphere(lat, delta_prime, + H_prime) + delta_e = atmospheric_refraction_correction(pressure, temp, e0, + atmos_refract) + e = topocentric_elevation_angle(e0, delta_e) + theta = topocentric_zenith_angle(e) + theta0 = topocentric_zenith_angle(e0) + gamma = topocentric_astronomers_azimuth(H_prime, delta_prime, lat) + phi = topocentric_azimuth_angle(gamma) + return theta, theta0, e, e0, phi, eot + + +def solar_position(unixtime, lat, lon, elev, pressure, temp, delta_t, + atmos_refract, numthreads=8, sst=False, esd=False): + + """ + Calculate the solar position using the + NREL SPA described in [1]_, [2]_. + + If numba is installed, the functions can be compiled + and the code runs quickly. If not, the functions + still evaluate but use numpy instead. + + Parameters + ---------- + unixtime : numpy array + Array of unix/epoch timestamps to calculate solar position for. + Unixtime is the number of seconds since Jan. 1, 1970 00:00:00 UTC. + A pandas.DatetimeIndex is easily converted using .view(np.int64)/10**9 + lat : float + Latitude to calculate solar position for + lon : float + Longitude to calculate solar position for + elev : float + Elevation of location in meters + pressure : int or float + avg. yearly pressure at location in millibars; + used for atmospheric correction + temp : int or float + avg. yearly temperature at location in + degrees C; used for atmospheric correction + delta_t : float or array + Difference between terrestrial time and UT1. + atmos_refract : float + The approximate atmospheric refraction (in degrees) + at sunrise and sunset. + numthreads: int, optional, default 8 + Number of threads to use for computation if numba>=0.17 + is installed. + sst : bool, default False + If True, return only data needed for sunrise, sunset, and transit + calculations. + esd : bool, default False + If True, return only Earth-Sun distance in AU + + Returns + ------- + Numpy Array with elements: + apparent zenith, + zenith, + apparent_elevation, + elevation, + azimuth, + equation_of_time + + References + ---------- + [1] I. Reda and A. Andreas, Solar position algorithm for solar radiation + applications. Solar Energy, vol. 76, no. 5, pp. 577-589, 2004. + + [2] I. Reda and A. Andreas, Corrigendum to Solar position algorithm for + solar radiation applications. Solar Energy, vol. 81, no. 6, p. 838, 2007. + """ + if USE_NUMBA: + do_calc = solar_position_numba + else: + do_calc = solar_position_numpy + + result = do_calc(unixtime, lat, lon, elev, pressure, + temp, delta_t, atmos_refract, numthreads, + sst, esd) + + if not isinstance(result, np.ndarray): + try: + result = np.array(result) + except Exception: + pass + + return result + + +def transit_sunrise_sunset(dates, lat, lon, delta_t, numthreads): + """ + Calculate the sun transit, sunrise, and sunset + for a set of dates at a given location using the NREL SPA [1]_, [2]_. + + Parameters + ---------- + dates : array + Numpy array of ints/floats corresponding to the Unix time + for the dates of interest, must be midnight UTC (00:00+00:00) + on the day of interest. + lat : float + Latitude of location to perform calculation for + lon : float + Longitude of location + delta_t : float or array + Difference between terrestrial time and UT. USNO has tables. + numthreads : int + Number to threads to use for calculation (if using numba) + + Returns + ------- + tuple : (transit, sunrise, sunset) localized to UTC + + References + ---------- + .. [1] I. Reda and A. Andreas, Solar position algorithm for solar + radiation applications. Solar Energy, vol. 76, no. 5, pp. 577-589, + 2004. :doi:`10.1016/j.solener.2003.12.003` References ---------- .. [1] I. Reda and A. Andreas, Solar position algorithm for solar @@ -1142,244 +1145,244 @@ def transit_sunrise_sunset(dates, lat, lon, delta_t, numthreads): (2004) 577–589] Solar Energy, vol. 81, no. 6, pp. 838, 2007. :doi:`10.1016/j.solener.2007.01.003` """ - - if ((dates % 86400) != 0.0).any(): - raise ValueError('Input dates must be at 00:00 UTC') - - utday = (dates // 86400) * 86400 - ttday0 = utday - delta_t - ttdayn1 = ttday0 - 86400 - ttdayp1 = ttday0 + 86400 - - # index 0 is v, 1 is alpha, 2 is delta - utday_res = solar_position(utday, 0, 0, 0, 0, 0, delta_t, - 0, numthreads, sst=True) - v = utday_res[0] - - ttday0_res = solar_position(ttday0, 0, 0, 0, 0, 0, delta_t, - 0, numthreads, sst=True) - ttdayn1_res = solar_position(ttdayn1, 0, 0, 0, 0, 0, delta_t, - 0, numthreads, sst=True) - ttdayp1_res = solar_position(ttdayp1, 0, 0, 0, 0, 0, delta_t, - 0, numthreads, sst=True) - m0 = (ttday0_res[1] - lon - v) / 360 - cos_arg = ((np.sin(np.radians(-0.8333)) - np.sin(np.radians(lat)) - * np.sin(np.radians(ttday0_res[2]))) / - (np.cos(np.radians(lat)) * np.cos(np.radians(ttday0_res[2])))) - cos_arg[abs(cos_arg) > 1] = np.nan - H0 = np.degrees(np.arccos(cos_arg)) % 180 - - m = np.empty((3, len(utday))) - m[0] = m0 % 1 - m[1] = (m[0] - H0 / 360) - m[2] = (m[0] + H0 / 360) - - # need to account for fractions of day that may be the next or previous - # day in UTC - add_a_day = m[2] >= 1 - sub_a_day = m[1] < 0 - m[1] = m[1] % 1 - m[2] = m[2] % 1 - vs = v + 360.985647 * m - n = m + delta_t / 86400 - - a = ttday0_res[1] - ttdayn1_res[1] - a[abs(a) > 2] = a[abs(a) > 2] % 1 - ap = ttday0_res[2] - ttdayn1_res[2] - ap[abs(ap) > 2] = ap[abs(ap) > 2] % 1 - b = ttdayp1_res[1] - ttday0_res[1] - b[abs(b) > 2] = b[abs(b) > 2] % 1 - bp = ttdayp1_res[2] - ttday0_res[2] - bp[abs(bp) > 2] = bp[abs(bp) > 2] % 1 - c = b - a - cp = bp - ap - - alpha_prime = ttday0_res[1] + (n * (a + b + c * n)) / 2 - delta_prime = ttday0_res[2] + (n * (ap + bp + cp * n)) / 2 - Hp = (vs + lon - alpha_prime) % 360 - Hp[Hp >= 180] = Hp[Hp >= 180] - 360 - - h = np.degrees(np.arcsin(np.sin(np.radians(lat)) * - np.sin(np.radians(delta_prime)) + - np.cos(np.radians(lat)) * - np.cos(np.radians(delta_prime)) - * np.cos(np.radians(Hp)))) - - T = (m[0] - Hp[0] / 360) * 86400 - R = (m[1] + (h[1] + 0.8333) / (360 * np.cos(np.radians(delta_prime[1])) * - np.cos(np.radians(lat)) * - np.sin(np.radians(Hp[1])))) * 86400 - S = (m[2] + (h[2] + 0.8333) / (360 * np.cos(np.radians(delta_prime[2])) * - np.cos(np.radians(lat)) * - np.sin(np.radians(Hp[2])))) * 86400 - - S[add_a_day] += 86400 - R[sub_a_day] -= 86400 - - transit = T + utday - sunrise = R + utday - sunset = S + utday - - return transit, sunrise, sunset - - -def earthsun_distance(unixtime, delta_t, numthreads): - """ - Calculates the distance from the earth to the sun using the - NREL SPA algorithm described in [1]. - - Parameters - ---------- - unixtime : numpy array - Array of unix/epoch timestamps to calculate solar position for. - Unixtime is the number of seconds since Jan. 1, 1970 00:00:00 UTC. - A pandas.DatetimeIndex is easily converted using .view(np.int64)/10**9 - delta_t : float or array - Difference between terrestrial time and UT. USNO has tables. - numthreads : int - Number to threads to use for calculation (if using numba) - - Returns - ------- - R : array - Earth-Sun distance in AU. - - References - ---------- - [1] Reda, I., Andreas, A., 2003. Solar position algorithm for solar - radiation applications. Technical report: NREL/TP-560-34302. Golden, - USA, http://www.nlr.gov. - """ - - R = solar_position(unixtime, 0, 0, 0, 0, 0, delta_t, - 0, numthreads, esd=True)[0] - - return R - - -def calculate_deltat(year, month): - """Calculate the difference between Terrestrial Dynamical Time (TD) - and Universal Time (UT). - - Equations taken from http://eclipse.gsfc.nasa.gov/SEcat5/deltatpoly.html - """ - - plw = 'Deltat is unknown for years before -1999 and after 3000. ' \ - 'Delta values will be calculated, but the calculations ' \ - 'are not intended to be used for these years.' - - try: - if np.any((year > 3000) | (year < -1999)): - warnings.warn(plw) - except ValueError: - if (year > 3000) | (year < -1999): - warnings.warn(plw) - except TypeError: - return 0 - - y = year + (month - 0.5)/12 - - deltat = np.where(year < -500, - - -20+32*((y-1820)/100)**2, 0) - - deltat = np.where((-500 <= year) & (year < 500), - - 10583.6-1014.41*(y/100) - + 33.78311*(y/100)**2 - - 5.952053*(y/100)**3 - - 0.1798452*(y/100)**4 - + 0.022174192*(y/100)**5 - + 0.0090316521*(y/100)**6, deltat) - - deltat = np.where((500 <= year) & (year < 1600), - - 1574.2-556.01*((y-1000)/100) - + 71.23472*((y-1000)/100)**2 - + 0.319781*((y-1000)/100)**3 - - 0.8503463*((y-1000)/100)**4 - - 0.005050998*((y-1000)/100)**5 - + 0.0083572073*((y-1000)/100)**6, deltat) - - deltat = np.where((1600 <= year) & (year < 1700), - - 120-0.9808*(y-1600) - - 0.01532*(y-1600)**2 - + (y-1600)**3/7129, deltat) - - deltat = np.where((1700 <= year) & (year < 1800), - - 8.83+0.1603*(y-1700) - - 0.0059285*(y-1700)**2 - + 0.00013336*(y-1700)**3 - - (y-1700)**4/1174000, deltat) - - deltat = np.where((1800 <= year) & (year < 1860), - - 13.72-0.332447*(y-1800) - + 0.0068612*(y-1800)**2 - + 0.0041116*(y-1800)**3 - - 0.00037436*(y-1800)**4 - + 0.0000121272*(y-1800)**5 - - 0.0000001699*(y-1800)**6 - + 0.000000000875*(y-1800)**7, deltat) - - deltat = np.where((1860 <= year) & (year < 1900), - - 7.62+0.5737*(y-1860) - - 0.251754*(y-1860)**2 - + 0.01680668*(y-1860)**3 - - 0.0004473624*(y-1860)**4 - + (y-1860)**5/233174, deltat) - - deltat = np.where((1900 <= year) & (year < 1920), - - -2.79+1.494119*(y-1900) - - 0.0598939*(y-1900)**2 - + 0.0061966*(y-1900)**3 - - 0.000197*(y-1900)**4, deltat) - - deltat = np.where((1920 <= year) & (year < 1941), - - 21.20+0.84493*(y-1920) - - 0.076100*(y-1920)**2 - + 0.0020936*(y-1920)**3, deltat) - - deltat = np.where((1941 <= year) & (year < 1961), - - 29.07+0.407*(y-1950) - - (y-1950)**2/233 - + (y-1950)**3/2547, deltat) - - deltat = np.where((1961 <= year) & (year < 1986), - - 45.45+1.067*(y-1975) - - (y-1975)**2/260 - - (y-1975)**3/718, deltat) - - deltat = np.where((1986 <= year) & (year < 2005), - - 63.86+0.3345*(y-2000) - - 0.060374*(y-2000)**2 - + 0.0017275*(y-2000)**3 - + 0.000651814*(y-2000)**4 - + 0.00002373599*(y-2000)**5, deltat) - - deltat = np.where((2005 <= year) & (year < 2050), - - 62.92+0.32217*(y-2000) - + 0.005589*(y-2000)**2, deltat) - - deltat = np.where((2050 <= year) & (year < 2150), - - -20+32*((y-1820)/100)**2 - - 0.5628*(2150-y), deltat) - - deltat = np.where(year >= 2150, - - -20+32*((y-1820)/100)**2, deltat) - - deltat = deltat.item() if np.isscalar(year) & np.isscalar(month)\ - else deltat - - return deltat + + if ((dates % 86400) != 0.0).any(): + raise ValueError('Input dates must be at 00:00 UTC') + + utday = (dates // 86400) * 86400 + ttday0 = utday - delta_t + ttdayn1 = ttday0 - 86400 + ttdayp1 = ttday0 + 86400 + + # index 0 is v, 1 is alpha, 2 is delta + utday_res = solar_position(utday, 0, 0, 0, 0, 0, delta_t, + 0, numthreads, sst=True) + v = utday_res[0] + + ttday0_res = solar_position(ttday0, 0, 0, 0, 0, 0, delta_t, + 0, numthreads, sst=True) + ttdayn1_res = solar_position(ttdayn1, 0, 0, 0, 0, 0, delta_t, + 0, numthreads, sst=True) + ttdayp1_res = solar_position(ttdayp1, 0, 0, 0, 0, 0, delta_t, + 0, numthreads, sst=True) + m0 = (ttday0_res[1] - lon - v) / 360 + cos_arg = ((np.sin(np.radians(-0.8333)) - np.sin(np.radians(lat)) + * np.sin(np.radians(ttday0_res[2]))) / + (np.cos(np.radians(lat)) * np.cos(np.radians(ttday0_res[2])))) + cos_arg[abs(cos_arg) > 1] = np.nan + H0 = np.degrees(np.arccos(cos_arg)) % 180 + + m = np.empty((3, len(utday))) + m[0] = m0 % 1 + m[1] = (m[0] - H0 / 360) + m[2] = (m[0] + H0 / 360) + + # need to account for fractions of day that may be the next or previous + # day in UTC + add_a_day = m[2] >= 1 + sub_a_day = m[1] < 0 + m[1] = m[1] % 1 + m[2] = m[2] % 1 + vs = v + 360.985647 * m + n = m + delta_t / 86400 + + a = ttday0_res[1] - ttdayn1_res[1] + a[abs(a) > 2] = a[abs(a) > 2] % 1 + ap = ttday0_res[2] - ttdayn1_res[2] + ap[abs(ap) > 2] = ap[abs(ap) > 2] % 1 + b = ttdayp1_res[1] - ttday0_res[1] + b[abs(b) > 2] = b[abs(b) > 2] % 1 + bp = ttdayp1_res[2] - ttday0_res[2] + bp[abs(bp) > 2] = bp[abs(bp) > 2] % 1 + c = b - a + cp = bp - ap + + alpha_prime = ttday0_res[1] + (n * (a + b + c * n)) / 2 + delta_prime = ttday0_res[2] + (n * (ap + bp + cp * n)) / 2 + Hp = (vs + lon - alpha_prime) % 360 + Hp[Hp >= 180] = Hp[Hp >= 180] - 360 + + h = np.degrees(np.arcsin(np.sin(np.radians(lat)) * + np.sin(np.radians(delta_prime)) + + np.cos(np.radians(lat)) * + np.cos(np.radians(delta_prime)) + * np.cos(np.radians(Hp)))) + + T = (m[0] - Hp[0] / 360) * 86400 + R = (m[1] + (h[1] + 0.8333) / (360 * np.cos(np.radians(delta_prime[1])) * + np.cos(np.radians(lat)) * + np.sin(np.radians(Hp[1])))) * 86400 + S = (m[2] + (h[2] + 0.8333) / (360 * np.cos(np.radians(delta_prime[2])) * + np.cos(np.radians(lat)) * + np.sin(np.radians(Hp[2])))) * 86400 + + S[add_a_day] += 86400 + R[sub_a_day] -= 86400 + + transit = T + utday + sunrise = R + utday + sunset = S + utday + + return transit, sunrise, sunset + + +def earthsun_distance(unixtime, delta_t, numthreads): + """ + Calculates the distance from the earth to the sun using the + NREL SPA algorithm described in [1]. + + Parameters + ---------- + unixtime : numpy array + Array of unix/epoch timestamps to calculate solar position for. + Unixtime is the number of seconds since Jan. 1, 1970 00:00:00 UTC. + A pandas.DatetimeIndex is easily converted using .view(np.int64)/10**9 + delta_t : float or array + Difference between terrestrial time and UT. USNO has tables. + numthreads : int + Number to threads to use for calculation (if using numba) + + Returns + ------- + R : array + Earth-Sun distance in AU. + + References + ---------- + [1] Reda, I., Andreas, A., 2003. Solar position algorithm for solar + radiation applications. Technical report: NREL/TP-560-34302. Golden, + USA, http://www.nlr.gov. + """ + + R = solar_position(unixtime, 0, 0, 0, 0, 0, delta_t, + 0, numthreads, esd=True)[0] + + return R + + +def calculate_deltat(year, month): + """Calculate the difference between Terrestrial Dynamical Time (TD) + and Universal Time (UT). + + Equations taken from http://eclipse.gsfc.nasa.gov/SEcat5/deltatpoly.html + """ + + plw = 'Deltat is unknown for years before -1999 and after 3000. ' \ + 'Delta values will be calculated, but the calculations ' \ + 'are not intended to be used for these years.' + + try: + if np.any((year > 3000) | (year < -1999)): + warnings.warn(plw) + except ValueError: + if (year > 3000) | (year < -1999): + warnings.warn(plw) + except TypeError: + return 0 + + y = year + (month - 0.5)/12 + + deltat = np.where(year < -500, + + -20+32*((y-1820)/100)**2, 0) + + deltat = np.where((-500 <= year) & (year < 500), + + 10583.6-1014.41*(y/100) + + 33.78311*(y/100)**2 + - 5.952053*(y/100)**3 + - 0.1798452*(y/100)**4 + + 0.022174192*(y/100)**5 + + 0.0090316521*(y/100)**6, deltat) + + deltat = np.where((500 <= year) & (year < 1600), + + 1574.2-556.01*((y-1000)/100) + + 71.23472*((y-1000)/100)**2 + + 0.319781*((y-1000)/100)**3 + - 0.8503463*((y-1000)/100)**4 + - 0.005050998*((y-1000)/100)**5 + + 0.0083572073*((y-1000)/100)**6, deltat) + + deltat = np.where((1600 <= year) & (year < 1700), + + 120-0.9808*(y-1600) + - 0.01532*(y-1600)**2 + + (y-1600)**3/7129, deltat) + + deltat = np.where((1700 <= year) & (year < 1800), + + 8.83+0.1603*(y-1700) + - 0.0059285*(y-1700)**2 + + 0.00013336*(y-1700)**3 + - (y-1700)**4/1174000, deltat) + + deltat = np.where((1800 <= year) & (year < 1860), + + 13.72-0.332447*(y-1800) + + 0.0068612*(y-1800)**2 + + 0.0041116*(y-1800)**3 + - 0.00037436*(y-1800)**4 + + 0.0000121272*(y-1800)**5 + - 0.0000001699*(y-1800)**6 + + 0.000000000875*(y-1800)**7, deltat) + + deltat = np.where((1860 <= year) & (year < 1900), + + 7.62+0.5737*(y-1860) + - 0.251754*(y-1860)**2 + + 0.01680668*(y-1860)**3 + - 0.0004473624*(y-1860)**4 + + (y-1860)**5/233174, deltat) + + deltat = np.where((1900 <= year) & (year < 1920), + + -2.79+1.494119*(y-1900) + - 0.0598939*(y-1900)**2 + + 0.0061966*(y-1900)**3 + - 0.000197*(y-1900)**4, deltat) + + deltat = np.where((1920 <= year) & (year < 1941), + + 21.20+0.84493*(y-1920) + - 0.076100*(y-1920)**2 + + 0.0020936*(y-1920)**3, deltat) + + deltat = np.where((1941 <= year) & (year < 1961), + + 29.07+0.407*(y-1950) + - (y-1950)**2/233 + + (y-1950)**3/2547, deltat) + + deltat = np.where((1961 <= year) & (year < 1986), + + 45.45+1.067*(y-1975) + - (y-1975)**2/260 + - (y-1975)**3/718, deltat) + + deltat = np.where((1986 <= year) & (year < 2005), + + 63.86+0.3345*(y-2000) + - 0.060374*(y-2000)**2 + + 0.0017275*(y-2000)**3 + + 0.000651814*(y-2000)**4 + + 0.00002373599*(y-2000)**5, deltat) + + deltat = np.where((2005 <= year) & (year < 2050), + + 62.92+0.32217*(y-2000) + + 0.005589*(y-2000)**2, deltat) + + deltat = np.where((2050 <= year) & (year < 2150), + + -20+32*((y-1820)/100)**2 + - 0.5628*(2150-y), deltat) + + deltat = np.where(year >= 2150, + + -20+32*((y-1820)/100)**2, deltat) + + deltat = deltat.item() if np.isscalar(year) & np.isscalar(month)\ + else deltat + + return deltat From 11303bb53080c0df6f55f4bd93cfcbb6d0690cbe Mon Sep 17 00:00:00 2001 From: Kevin Anderson Date: Wed, 28 Oct 2026 10:07:02 -0400 Subject: [PATCH 6/7] re-apply suggestion --- pvlib/spa.py | 17 ++++++----------- 1 file changed, 6 insertions(+), 11 deletions(-) diff --git a/pvlib/spa.py b/pvlib/spa.py index 6e36a084cd..07582b9d22 100644 --- a/pvlib/spa.py +++ b/pvlib/spa.py @@ -1109,7 +1109,7 @@ def solar_position(unixtime, lat, lon, elev, pressure, temp, delta_t, def transit_sunrise_sunset(dates, lat, lon, delta_t, numthreads): """ Calculate the sun transit, sunrise, and sunset - for a set of dates at a given location using the NREL SPA [1]_, [2]_. + for a set of dates at a given location using the NREL SPA [1]_, [2]_. Parameters ---------- @@ -1135,16 +1135,11 @@ def transit_sunrise_sunset(dates, lat, lon, delta_t, numthreads): .. [1] I. Reda and A. Andreas, Solar position algorithm for solar radiation applications. Solar Energy, vol. 76, no. 5, pp. 577-589, 2004. :doi:`10.1016/j.solener.2003.12.003` - References - ---------- - .. [1] I. Reda and A. Andreas, Solar position algorithm for solar - radiation applications. Solar Energy, vol. 76, no. 5, pp. 577-589, - 2004. :doi:`10.1016/j.solener.2003.12.003` - .. [2] I. Reda, and A. Andreas. "Corrigendum to “Solar position - algorithm for solar radiation applications." [Solar Energy 76 - (2004) 577–589] Solar Energy, vol. 81, no. 6, pp. 838, 2007. - :doi:`10.1016/j.solener.2007.01.003` - """ + .. [2] I. Reda, and A. Andreas. "Corrigendum to “Solar position + algorithm for solar radiation applications." [Solar Energy 76 + (2004) 577–589] Solar Energy, vol. 81, no. 6, pp. 838, 2007. + :doi:`10.1016/j.solener.2007.01.003` + """ if ((dates % 86400) != 0.0).any(): raise ValueError('Input dates must be at 00:00 UTC') From 831c0821088f27cfc446796bfc4e1498ec0030cd Mon Sep 17 00:00:00 2001 From: Kevin Anderson Date: Wed, 28 Oct 2026 10:09:32 -0400 Subject: [PATCH 7/7] lint and rst fix --- pvlib/spa.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/pvlib/spa.py b/pvlib/spa.py index 07582b9d22..9f1a3a345c 100644 --- a/pvlib/spa.py +++ b/pvlib/spa.py @@ -1109,7 +1109,9 @@ def solar_position(unixtime, lat, lon, elev, pressure, temp, delta_t, def transit_sunrise_sunset(dates, lat, lon, delta_t, numthreads): """ Calculate the sun transit, sunrise, and sunset - for a set of dates at a given location using the NREL SPA [1]_, [2]_. + for a set of dates at a given location using the NREL SPA. + + The NREL SPA is documented in [1]_, [2]_. Parameters ---------- @@ -1135,7 +1137,7 @@ def transit_sunrise_sunset(dates, lat, lon, delta_t, numthreads): .. [1] I. Reda and A. Andreas, Solar position algorithm for solar radiation applications. Solar Energy, vol. 76, no. 5, pp. 577-589, 2004. :doi:`10.1016/j.solener.2003.12.003` - .. [2] I. Reda, and A. Andreas. "Corrigendum to “Solar position + .. [2] I. Reda, and A. Andreas. "Corrigendum to “Solar position algorithm for solar radiation applications." [Solar Energy 76 (2004) 577–589] Solar Energy, vol. 81, no. 6, pp. 838, 2007. :doi:`10.1016/j.solener.2007.01.003`