""" Calculate the solar position using the NREL SPA algorithm 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, __version__ except ImportError: warnings.warn('Could not import numba, falling back to numpy ' + 'calculation') jcompile = nocompile USE_NUMBA = False else: major, minor = __version__.split('.')[:2] if int(major + minor) >= 17: # need at least numba >= 0.17.0 jcompile = jit USE_NUMBA = True else: warnings.warn('Numba version must be >= 0.17.0, falling back to ' + 'numpy') jcompile = nocompile USE_NUMBA = False else: jcompile = nocompile USE_NUMBA = False TABLE_1_DICT = { '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]]), '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]]), '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]]) } resize_mapping = { 'L1': (64, 3), 'L2': (64, 3), 'L3': (64, 3), 'L4': (64, 3), 'L5': (64, 3), 'B1': (5, 3), 'R1': (40, 3), 'R2': (40, 3), 'R3': (40, 3), 'R4': (40, 3)} # make arrays uniform size for efficient broadcasting in numba, fill with 0s # np.resize does not work because it fills with repeated copies for key, dims in resize_mapping.items(): new_rows = dims[0] - TABLE_1_DICT[key].shape[0] TABLE_1_DICT[key] = np.append(TABLE_1_DICT[key], np.zeros((new_rows, 3)), axis=0) HELIO_LONG_TABLE = np.array([TABLE_1_DICT['L0'], TABLE_1_DICT['L1'], TABLE_1_DICT['L2'], TABLE_1_DICT['L3'], TABLE_1_DICT['L4'], TABLE_1_DICT['L5']]) HELIO_LAT_TABLE = np.array([TABLE_1_DICT['B0'], TABLE_1_DICT['B1']]) HELIO_RADIUS_TABLE = np.array([TABLE_1_DICT['R0'], TABLE_1_DICT['R1'], TABLE_1_DICT['R2'], TABLE_1_DICT['R3'], TABLE_1_DICT['R4']]) 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 + b - 1524.5) 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 @jcompile('float64(float64)', nopython=True) def heliocentric_longitude(jme): l0 = 0.0 l1 = 0.0 l2 = 0.0 l3 = 0.0 l4 = 0.0 l5 = 0.0 for row in range(HELIO_LONG_TABLE.shape[1]): l0 += (HELIO_LONG_TABLE[0, row, 0] * np.cos(HELIO_LONG_TABLE[0, row, 1] + HELIO_LONG_TABLE[0, row, 2] * jme) ) l1 += (HELIO_LONG_TABLE[1, row, 0] * np.cos(HELIO_LONG_TABLE[1, row, 1] + HELIO_LONG_TABLE[1, row, 2] * jme) ) l2 += (HELIO_LONG_TABLE[2, row, 0] * np.cos(HELIO_LONG_TABLE[2, row, 1] + HELIO_LONG_TABLE[2, row, 2] * jme) ) l3 += (HELIO_LONG_TABLE[3, row, 0] * np.cos(HELIO_LONG_TABLE[3, row, 1] + HELIO_LONG_TABLE[3, row, 2] * jme) ) l4 += (HELIO_LONG_TABLE[4, row, 0] * np.cos(HELIO_LONG_TABLE[4, row, 1] + HELIO_LONG_TABLE[4, row, 2] * jme) ) l5 += (HELIO_LONG_TABLE[5, row, 0] * np.cos(HELIO_LONG_TABLE[5, row, 1] + HELIO_LONG_TABLE[5, row, 2] * 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 = 0.0 b1 = 0.0 for row in range(HELIO_LAT_TABLE.shape[1]): b0 += (HELIO_LAT_TABLE[0, row, 0] * np.cos(HELIO_LAT_TABLE[0, row, 1] + HELIO_LAT_TABLE[0, row, 2] * jme) ) b1 += (HELIO_LAT_TABLE[1, row, 0] * np.cos(HELIO_LAT_TABLE[1, row, 1] + HELIO_LAT_TABLE[1, row, 2] * 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 = 0.0 r1 = 0.0 r2 = 0.0 r3 = 0.0 r4 = 0.0 for row in range(HELIO_RADIUS_TABLE.shape[1]): r0 += (HELIO_RADIUS_TABLE[0, row, 0] * np.cos(HELIO_RADIUS_TABLE[0, row, 1] + HELIO_RADIUS_TABLE[0, row, 2] * jme) ) r1 += (HELIO_RADIUS_TABLE[1, row, 0] * np.cos(HELIO_RADIUS_TABLE[1, row, 1] + HELIO_RADIUS_TABLE[1, row, 2] * jme) ) r2 += (HELIO_RADIUS_TABLE[2, row, 0] * np.cos(HELIO_RADIUS_TABLE[2, row, 1] + HELIO_RADIUS_TABLE[2, row, 2] * jme) ) r3 += (HELIO_RADIUS_TABLE[3, row, 0] * np.cos(HELIO_RADIUS_TABLE[3, row, 1] + HELIO_RADIUS_TABLE[3, row, 2] * jme) ) r4 += (HELIO_RADIUS_TABLE[4, row, 0] * np.cos(HELIO_RADIUS_TABLE[4, row, 1] + HELIO_RADIUS_TABLE[4, row, 2] * 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('float64(float64, float64, float64, float64, float64, float64)', nopython=True) def longitude_nutation(julian_ephemeris_century, x0, x1, x2, x3, x4): delta_psi_sum = 0 for row in range(NUTATION_YTERM_ARRAY.shape[0]): a = NUTATION_ABCD_ARRAY[row, 0] b = NUTATION_ABCD_ARRAY[row, 1] argsin = (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) term = (a + b * julian_ephemeris_century) * np.sin(np.radians(argsin)) delta_psi_sum += term delta_psi = delta_psi_sum*1.0/36000000 return delta_psi @jcompile('float64(float64, float64, float64, float64, float64, float64)', nopython=True) def obliquity_nutation(julian_ephemeris_century, x0, x1, x2, x3, x4): delta_eps_sum = 0.0 for row in range(NUTATION_YTERM_ARRAY.shape[0]): c = NUTATION_ABCD_ARRAY[row, 2] d = NUTATION_ABCD_ARRAY[row, 3] argcos = (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) term = (c + d * julian_ephemeris_century) * np.cos(np.radians(argcos)) delta_eps_sum += term delta_eps = delta_eps_sum*1.0/36000000 return 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): num = (np.sin(np.radians(apparent_sun_longitude)) * np.cos(np.radians(true_ecliptic_obliquity)) - np.tan(np.radians(geocentric_latitude)) * np.sin(np.radians(true_ecliptic_obliquity))) alpha = np.degrees(np.arctan2(num, np.cos( np.radians(apparent_sun_longitude)))) return alpha % 360 @jcompile('float64(float64, float64, float64)', nopython=True) def geocentric_sun_declination(apparent_sun_longitude, true_ecliptic_obliquity, geocentric_latitude): delta = np.degrees(np.arcsin(np.sin(np.radians(geocentric_latitude)) * np.cos(np.radians(true_ecliptic_obliquity)) + np.cos(np.radians(geocentric_latitude)) * np.sin(np.radians(true_ecliptic_obliquity)) * 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): num = (-xterm * np.sin(np.radians(equatorial_horizontal_parallax)) * np.sin(np.radians(local_hour_angle))) denom = (np.cos(np.radians(geocentric_sun_declination)) - xterm * np.sin(np.radians(equatorial_horizontal_parallax)) * np.cos(np.radians(local_hour_angle))) 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): num = ((np.sin(np.radians(geocentric_sun_declination)) - yterm * np.sin(np.radians(equatorial_horizontal_parallax))) * np.cos(np.radians(parallax_sun_right_ascension))) denom = (np.cos(np.radians(geocentric_sun_declination)) - xterm * np.sin(np.radians(equatorial_horizontal_parallax)) * 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 ): e0 = np.degrees(np.arcsin( np.sin(np.radians(observer_latitude)) * np.sin(np.radians(topocentric_sun_declination)) + np.cos(np.radians(observer_latitude)) * np.cos(np.radians(topocentric_sun_declination)) * 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): num = np.sin(np.radians(topocentric_local_hour_angle)) denom = (np.cos(np.radians(topocentric_local_hour_angle)) * np.sin(np.radians(observer_latitude)) - np.tan(np.radians(topocentric_sun_declination)) * np.cos(np.radians(observer_latitude))) 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[:,:])', nopython=True, nogil=True) def solar_position_loop(unixtime, 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] delta_t = loc_args[5] atmos_refract = loc_args[6] sst = loc_args[7] esd = loc_args[8] for i in range(unixtime.shape[0]): utime = unixtime[i] jd = julian_day(utime) 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: 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) delta_psi = longitude_nutation(jce, x0, x1, x2, x3, x4) delta_epsilon = obliquity_nutation(jce, x0, x1, x2, x3, x4) 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, delta_t, atmos_refract, sst, esd]) # 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: warnings.warn('The number of threads is more than the length of ' 'the time array. Only using %s threads.'.format(ulength)) numthreads = ulength if numthreads <= 1: solar_position_loop(unixtime, loc_args, result) return result # split the input and output arrays into numthreads chunks split0 = np.array_split(unixtime, numthreads) split2 = np.array_split(result, numthreads, axis=1) chunks = [[a0, loc_args, split2[i]] for i, a0 in enumerate(split0)] # 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) delta_psi = longitude_nutation(jce, x0, x1, x2, x3, x4) delta_epsilon = obliquity_nutation(jce, x0, x1, x2, x3, x4) 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 algorithm described in [1]. 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 .astype(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, optional If delta_t is None, uses spa.calculate_deltat using time.year and time.month from pandas.DatetimeIndex. For most simulations specifing delta_t is sufficient. Difference between terrestrial time and UT1. *Note: delta_t = None will break code using nrel_numba, this will be fixed in a future version. By default, use USNO historical data and predictions atmos_refrac : float, optional 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, elevation, apparent_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. 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 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 """ 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 .astype(np.int64)/10**9 delta_t : float 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.nrel.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). Note: This function is not yet compatible for calculations using Numba. 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