# -*- coding: utf-8 -*- """ irradiance.py from pvlib ======================== Stripped down, vendorized version from: https://github.com/pvlib/pvlib-python/ Calculate the solar position using the NREL SPA algorithm either using numpy arrays or compiling the code to machine language with numba. The rational for not including this library as a strict dependency is to avoid including a dependency on pandas, keeping load time low, and PyPy compatibility Created by Tony Lorenzo (@alorenzo175), Univ. of Arizona, 2015 For a full list of contributors to this file, see the `pvlib` repository. The copyright notice (BSD-3 clause) is as follows: BSD 3-Clause License Copyright (c) 2013-2018, Sandia National Laboratories and pvlib python Development Team All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. Neither the name of the {organization} nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. """ from __future__ import division import os import time from datetime import datetime import math from math import degrees, sin, cos, tan, radians, asin, atan2 import numpy as np ndarray = np.ndarray 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]]) } # Resizing is just adding zeros... not sure why #print(TABLE_1_DICT['L1'].shape) #print(TABLE_1_DICT['L1'].tolist()) TABLE_1_DICT['L1'].resize((64, 3), refcheck=False) #print(TABLE_1_DICT['L1'].shape) #print(TABLE_1_DICT['L1'].tolist()) TABLE_1_DICT['L2'].resize((64, 3), refcheck=False) TABLE_1_DICT['L3'].resize((64, 3), refcheck=False) TABLE_1_DICT['L4'].resize((64, 3), refcheck=False) TABLE_1_DICT['L5'].resize((64, 3), refcheck=False) TABLE_1_DICT['B1'].resize((5, 3), refcheck=False) TABLE_1_DICT['R1'].resize((40, 3), refcheck=False) TABLE_1_DICT['R2'].resize((40, 3), refcheck=False) TABLE_1_DICT['R3'].resize((40, 3), refcheck=False) TABLE_1_DICT['R4'].resize((40, 3), refcheck=False) 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], ]) HELIO_LONG_TABLE_LIST = HELIO_LONG_TABLE.tolist() HELIO_RADIUS_TABLE_LIST = HELIO_RADIUS_TABLE.tolist() NUTATION_YTERM_LIST = NUTATION_YTERM_ARRAY.tolist() NUTATION_ABCD_LIST = NUTATION_ABCD_ARRAY.tolist() HELIO_LAT_TABLE_LIST = HELIO_LAT_TABLE.tolist() 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 + (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 def julian_day(unixtime): jd = unixtime/86400.0 + 2440587.5 return jd def julian_ephemeris_day(julian_day, delta_t): jde = julian_day + delta_t * 1.0 / 86400 return jde def julian_century(julian_day): jc = (julian_day - 2451545) * 1.0 / 36525 return jc def julian_ephemeris_century(julian_ephemeris_day): jce = (julian_ephemeris_day - 2451545)/36525.0 return jce def julian_ephemeris_millennium(julian_ephemeris_century): jme = julian_ephemeris_century*0.1 return jme def heliocentric_longitude(jme): # Might be able to replace this with a pade approximation? # Looping over rows is probably still faster than (a, b, c) # Maximum optimization l0 = 0.0 l1 = 0.0 l2 = 0.0 l3 = 0.0 l4 = 0.0 l5 = 0.0 if isinstance(jme, ndarray): cos = np.cos rad2deg = np.rad2deg else: cos = math.cos rad2deg = math.degrees HELIO_LONG_TABLE_LIST_0 = HELIO_LONG_TABLE_LIST[0] HELIO_LONG_TABLE_LIST_1 = HELIO_LONG_TABLE_LIST[1] HELIO_LONG_TABLE_LIST_2 = HELIO_LONG_TABLE_LIST[2] HELIO_LONG_TABLE_LIST_3 = HELIO_LONG_TABLE_LIST[3] HELIO_LONG_TABLE_LIST_4 = HELIO_LONG_TABLE_LIST[4] for row in range(64): HELIO_LONG_TABLE_LIST_0_ROW = HELIO_LONG_TABLE_LIST_0[row] l0 += (HELIO_LONG_TABLE_LIST_0_ROW[0] * cos(HELIO_LONG_TABLE_LIST_0_ROW[1] + HELIO_LONG_TABLE_LIST_0_ROW[2] * jme) ) for row in range(34): HELIO_LONG_TABLE_LIST_1_ROW = HELIO_LONG_TABLE_LIST_1[row] l1 += (HELIO_LONG_TABLE_LIST_1_ROW[0] * cos(HELIO_LONG_TABLE_LIST_1_ROW[1] + HELIO_LONG_TABLE_LIST_1_ROW[2] * jme) ) for row in range(20): HELIO_LONG_TABLE_LIST_2_ROW = HELIO_LONG_TABLE_LIST_2[row] l2 += (HELIO_LONG_TABLE_LIST_2_ROW[0] * cos(HELIO_LONG_TABLE_LIST_2_ROW[1] + HELIO_LONG_TABLE_LIST_2_ROW[2] * jme) ) for row in range(7): HELIO_LONG_TABLE_LIST_3_ROW = HELIO_LONG_TABLE_LIST_3[row] l3 += (HELIO_LONG_TABLE_LIST_3_ROW[0] * cos(HELIO_LONG_TABLE_LIST_3_ROW[1] + HELIO_LONG_TABLE_LIST_3_ROW[2] * jme) ) for row in range(3): HELIO_LONG_TABLE_LIST_4_ROW = HELIO_LONG_TABLE_LIST_4[row] l4 += (HELIO_LONG_TABLE_LIST_4_ROW[0] * cos(HELIO_LONG_TABLE_LIST_4_ROW[1] + HELIO_LONG_TABLE_LIST_4_ROW[2] * jme) ) l5 = (HELIO_LONG_TABLE_LIST[5][0][0] * cos(HELIO_LONG_TABLE_LIST[5][0][1])) l_rad = (jme*(jme*(jme*(jme*(jme*l5 + l4) + l3) + l2) + l1) + l0)*1E-8 l = rad2deg(l_rad) return l % 360 def heliocentric_latitude(jme): b0 = 0.0 b1 = 0.0 if isinstance(jme, ndarray): cos = np.cos rad2deg = np.rad2deg else: cos = math.cos rad2deg = math.degrees HELIO_LAT_TABLE_LIST_0 = HELIO_LAT_TABLE_LIST[0] HELIO_LAT_TABLE_LIST_1 = HELIO_LAT_TABLE_LIST[1] for row in range(5): HELIO_LAT_TABLE_LIST_0_ROW = HELIO_LAT_TABLE_LIST[0][row] b0 += (HELIO_LAT_TABLE_LIST_0_ROW[0] * cos(HELIO_LAT_TABLE_LIST_0_ROW[1] + HELIO_LAT_TABLE_LIST_0_ROW[2] * jme) ) for row in range(2): HELIO_LAT_TABLE_LIST_1_ROW = HELIO_LAT_TABLE_LIST[1][row] b1 += (HELIO_LAT_TABLE_LIST_1_ROW[0] * cos(HELIO_LAT_TABLE_LIST_1_ROW[1] + HELIO_LAT_TABLE_LIST_1_ROW[2] * jme) ) b_rad = (b0 + b1 * jme)*1E-8 b = rad2deg(b_rad) return b def heliocentric_radius_vector(jme): # no optimizations can be thought of r0 = 0.0 r1 = 0.0 r2 = 0.0 r3 = 0.0 r4 = 0.0 if isinstance(jme, ndarray): cos = np.cos else: cos = math.cos # Would be possible to save a few multiplies of table1row[2]*jme, table1row[1]*jme as they are dups table0, table1, table2, table3, table4 = HELIO_RADIUS_TABLE_LIST[0:5] for row in range(40): table0row = table0[row] r0 += (table0row[0]*cos(table0row[1] + table0row[2]*jme)) for row in range(10): table1row = table1[row] r1 += (table1row[0]*cos(table1row[1] + table1row[2]*jme)) for row in range(6): table2row = table2[row] r2 += (table2row[0]*cos(table2row[1] + table2row[2]*jme)) for row in range(2): table3row = table3[row] r3 += (table3row[0]*cos(table3row[1] + table3row[2]*jme)) table4row = table4[0] r4 = (table4row[0]*cos(table4row[1] + table4row[2]*jme)) return (jme*(jme*(jme*(jme*r4 + r3) + r2) + r1) + r0)*1E-8 def geocentric_longitude(heliocentric_longitude): theta = heliocentric_longitude + 180.0 return theta % 360 def geocentric_latitude(heliocentric_latitude): beta = -heliocentric_latitude return beta 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 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 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 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 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 def longitude_nutation(julian_ephemeris_century, x0, x1, x2, x3, x4): if isinstance(julian_ephemeris_century, ndarray): sin = np.sin radians = np.radians else: sin = math.sin radians = math.radians x0, x1, x2, x3, x4 = radians(x0), radians(x1), radians(x2), radians(x3), radians(x4) delta_psi_sum = 0.0 for row in range(63): a = NUTATION_ABCD_LIST[row][0] b = NUTATION_ABCD_LIST[row][1] # # None can be skipped but the multiplies can be with effort -2 to 2 with dict - just might be slower NUTATION_YTERM_LIST_ROW = NUTATION_YTERM_LIST[row] argsin = (NUTATION_YTERM_LIST_ROW[0]*x0 + NUTATION_YTERM_LIST_ROW[1]*x1 + NUTATION_YTERM_LIST_ROW[2]*x2 + NUTATION_YTERM_LIST_ROW[3]*x3 + NUTATION_YTERM_LIST_ROW[4]*x4) term = (a + b * julian_ephemeris_century)*sin(argsin) delta_psi_sum += term delta_psi = delta_psi_sum/36000000.0 return delta_psi def obliquity_nutation(julian_ephemeris_century, x0, x1, x2, x3, x4): delta_eps_sum = 0.0 if isinstance(julian_ephemeris_century, ndarray): cos = np.cos radians = np.radians else: cos = math.cos radians = math.radians x0, x1, x2, x3, x4 = radians(x0), radians(x1), radians(x2), radians(x3), radians(x4) for row in range(63): c = NUTATION_ABCD_LIST[row][2] d = NUTATION_ABCD_LIST[row][3] NUTATION_YTERM_LIST_ROW = NUTATION_YTERM_LIST[row] argcos = (NUTATION_YTERM_LIST_ROW[0]*x0 + NUTATION_YTERM_LIST_ROW[1]*x1 + NUTATION_YTERM_LIST_ROW[2]*x2 + NUTATION_YTERM_LIST_ROW[3]*x3 + NUTATION_YTERM_LIST_ROW[4]*x4) term = (c + d * julian_ephemeris_century)*cos(argcos) delta_eps_sum += term delta_eps = delta_eps_sum/36000000.0 return delta_eps def mean_ecliptic_obliquity(julian_ephemeris_millennium): U = 0.1*julian_ephemeris_millennium e0 = (U*(U*(U*(U*(U*(U*(U*(U*(U*(2.45*U + 5.79) + 27.87) + 7.12) - 39.05) - 249.67) - 51.38) + 1999.25) - 1.55) - 4680.93) + 84381.448) return e0 def true_ecliptic_obliquity(mean_ecliptic_obliquity, obliquity_nutation): e0 = mean_ecliptic_obliquity deleps = obliquity_nutation e = e0/3600.0 + deleps return e def aberration_correction(earth_radius_vector): # -20.4898 / (3600) deltau = -0.005691611111111111/earth_radius_vector return deltau def apparent_sun_longitude(geocentric_longitude, longitude_nutation, aberration_correction): lamd = geocentric_longitude + longitude_nutation + aberration_correction return lamd def mean_sidereal_time(julian_day, julian_century): julian_century2 = julian_century*julian_century v0 = (280.46061837 + 360.98564736629*(julian_day - 2451545) + 0.000387933*julian_century2 - julian_century2*julian_century/38710000) return v0 % 360.0 def apparent_sidereal_time(mean_sidereal_time, longitude_nutation, true_ecliptic_obliquity): if isinstance(true_ecliptic_obliquity, ndarray): cos = np.cos radians = np.radians else: cos = math.cos radians = math.radians v = mean_sidereal_time + longitude_nutation*cos(radians(true_ecliptic_obliquity)) return v def geocentric_sun_right_ascension(apparent_sun_longitude, true_ecliptic_obliquity, geocentric_latitude): if isinstance(apparent_sun_longitude, ndarray): sin = np.sin cos = np.cos tan = np.tan radians = np.radians degrees = np.degrees arctan2 = np.arctan2 else: sin = math.sin cos = math.cos tan = math.tan radians = math.radians degrees = math.degrees arctan2 = math.atan2 num = (sin(radians(apparent_sun_longitude)) * cos(radians(true_ecliptic_obliquity)) - tan(radians(geocentric_latitude)) * sin(radians(true_ecliptic_obliquity))) alpha = degrees(arctan2(num, cos( radians(apparent_sun_longitude)))) return alpha % 360 def geocentric_sun_declination(apparent_sun_longitude, true_ecliptic_obliquity, geocentric_latitude): if isinstance(apparent_sun_longitude, ndarray): sin = np.sin cos = np.cos radians = np.radians degrees = np.degrees arcsin = np.arcsin else: sin = math.sin cos = math.cos radians = math.radians degrees = math.degrees arcsin = math.asin delta = degrees(arcsin(sin(radians(geocentric_latitude)) * cos(radians(true_ecliptic_obliquity)) + cos(radians(geocentric_latitude)) * sin(radians(true_ecliptic_obliquity)) * sin(radians(apparent_sun_longitude)))) return delta 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 def equatorial_horizontal_parallax(earth_radius_vector): xi = 8.794 / (3600 * earth_radius_vector) return xi def uterm(observer_latitude): if isinstance(observer_latitude, ndarray): tan = np.tan arctan = np.arctan radians = np.radians else: tan = math.tan radians = math.radians arctan = math.atan u = arctan(0.99664719*tan(radians(observer_latitude))) return u def xterm(u, observer_latitude, observer_elevation): if isinstance(observer_latitude, ndarray) or isinstance(u, ndarray): cos = np.cos radians = np.radians else: cos = math.cos radians = math.radians x = (cos(u) + observer_elevation/6378140.0*cos(radians(observer_latitude))) return x def yterm(u, observer_latitude, observer_elevation): if isinstance(observer_latitude, ndarray) or isinstance(u, ndarray): sin = np.sin radians = np.radians else: sin = math.sin radians = math.radians y = (0.99664719 * sin(u) + observer_elevation/6378140.0 * sin(radians(observer_latitude))) return y def parallax_sun_right_ascension(xterm, equatorial_horizontal_parallax, local_hour_angle, geocentric_sun_declination): if (isinstance(xterm, ndarray) or isinstance(equatorial_horizontal_parallax, ndarray) or isinstance(local_hour_angle, ndarray) or isinstance(geocentric_sun_declination, ndarray)): sin = np.sin cos = np.cos radians = np.radians degrees = np.degrees arctan2 = np.arctan2 else: sin = math.sin cos = math.cos radians = math.radians degrees = math.degrees arctan2 = math.atan2 num = (-xterm * sin(radians(equatorial_horizontal_parallax)) * sin(radians(local_hour_angle))) denom = (cos(radians(geocentric_sun_declination)) - xterm * sin(radians(equatorial_horizontal_parallax)) * cos(radians(local_hour_angle))) delta_alpha = degrees(arctan2(num, denom)) return delta_alpha 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 def topocentric_sun_declination(geocentric_sun_declination, xterm, yterm, equatorial_horizontal_parallax, parallax_sun_right_ascension, local_hour_angle): if isinstance(geocentric_sun_declination, ndarray): sin = np.sin cos = np.cos radians = np.radians degrees = np.degrees arctan2 = np.arctan2 else: sin = math.sin cos = math.cos radians = math.radians degrees = math.degrees arctan2 = math.atan2 num = ((sin(radians(geocentric_sun_declination)) - yterm * sin(radians(equatorial_horizontal_parallax))) * cos(radians(parallax_sun_right_ascension))) denom = (cos(radians(geocentric_sun_declination)) - xterm * sin(radians(equatorial_horizontal_parallax)) * cos(radians(local_hour_angle))) delta = degrees(arctan2(num, denom)) return delta def topocentric_local_hour_angle(local_hour_angle, parallax_sun_right_ascension): H_prime = local_hour_angle - parallax_sun_right_ascension return H_prime def topocentric_elevation_angle_without_atmosphere(observer_latitude, topocentric_sun_declination, topocentric_local_hour_angle ): if (isinstance(observer_latitude, ndarray) or isinstance(topocentric_sun_declination, ndarray) or isinstance(topocentric_local_hour_angle, ndarray)): 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)))) else: observer_latitude = float(observer_latitude) topocentric_sun_declination = float(topocentric_sun_declination) topocentric_local_hour_angle = float(topocentric_local_hour_angle) r_observer_latitude = radians(observer_latitude) r_topocentric_sun_declination = radians(topocentric_sun_declination) e0 = degrees(asin( sin(r_observer_latitude) * sin(r_topocentric_sun_declination) + cos(r_observer_latitude) * cos(r_topocentric_sun_declination) * cos(radians(topocentric_local_hour_angle)))) return e0 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) if isinstance(topocentric_elevation_angle_wo_atmosphere, ndarray): tan, radians = np.tan, np.radians else: tan, radians = math.tan, math.radians delta_e = ((local_pressure / 1010.0) * (283.0 / (273.0 + local_temp)) * 1.02 / (60 * tan(radians( topocentric_elevation_angle_wo_atmosphere + 10.3 / (topocentric_elevation_angle_wo_atmosphere + 5.11))))) * switch return delta_e def topocentric_elevation_angle(topocentric_elevation_angle_without_atmosphere, atmospheric_refraction_correction): e = (topocentric_elevation_angle_without_atmosphere + atmospheric_refraction_correction) return e def topocentric_zenith_angle(topocentric_elevation_angle): theta = 90.0 - topocentric_elevation_angle return theta def topocentric_astronomers_azimuth(topocentric_local_hour_angle, topocentric_sun_declination, observer_latitude): if (isinstance(topocentric_local_hour_angle, ndarray) or isinstance(topocentric_sun_declination, ndarray) or isinstance(observer_latitude, ndarray)): sin = np.sin cos = np.cos tan = np.tan radians = np.radians degrees = np.degrees arctan2 = np.arctan2 else: sin = math.sin cos = math.cos tan = math.tan radians = math.radians degrees = math.degrees arctan2 = math.atan2 num = sin(radians(topocentric_local_hour_angle)) denom = (cos(radians(topocentric_local_hour_angle)) * sin(radians(observer_latitude)) - tan(radians(topocentric_sun_declination)) * cos(radians(observer_latitude))) gamma = degrees(arctan2(num, denom)) return gamma % 360.0 def topocentric_azimuth_angle(topocentric_astronomers_azimuth): phi = topocentric_astronomers_azimuth + 180.0 return phi % 360.0 def sun_mean_longitude(julian_ephemeris_millennium): M = julian_ephemeris_millennium*(julian_ephemeris_millennium*( julian_ephemeris_millennium*(julian_ephemeris_millennium*( -5.0e-7*julian_ephemeris_millennium - 6.5359477124183e-5) + 2.00276381406341e-5) + 0.03032028) + 360007.6982779) + 280.4664567 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): if isinstance(true_ecliptic_obliquity, ndarray): term = np.cos(np.radians(true_ecliptic_obliquity)) else: term = cos(radians(true_ecliptic_obliquity)) E = (sun_mean_longitude - 0.0057183 - geocentric_sun_right_ascension + longitude_nutation * term) # limit between 0 and 360 E = E % 360 # convert to minutes E *= 4.0 greater = E > 20.0 less = E < -20.0 other = (E <= 20.0) & (E >= -20.0) E = greater * (E - 1440.0) + less * (E + 1440.0) + other * E return E def earthsun_distance(unixtime, delta_t, numthreads=1): """ 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 solar_position_numpy(unixtime, lat, lon, elev, pressure, temp, delta_t, atmos_refract, numthreads=None, 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 None 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 ------- list 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. """ 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) # allows ndarrays here # Clean of mostly numpy though here if esd: # This is all that is used in earth sun distance return [R] # numpy dirty L = heliocentric_longitude(jme) B = heliocentric_latitude(jme) # nbumpy dirty 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, array=False): result = solar_position_numpy(unixtime, lat, lon, elev, pressure, temp, delta_t, atmos_refract, numthreads, sst, esd) if not array: return result if not isinstance(result, 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 """ isnumpy = isinstance(dates, ndarray) if isnumpy: condition = ((dates % 86400) != 0.0).any() else: condition = (dates % 86400) != 0.0 if condition: 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])))) if isnumpy: cos_arg[abs(cos_arg) > 1] = np.nan else: if abs(cos_arg) > 1: cos_arg = np.nan H0 = np.degrees(np.arccos(cos_arg)) % 180 if isnumpy: m = np.empty((3, len(utday))) else: m = np.empty((3, 1)) 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] if isnumpy: 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 else: if abs(a) > 2: a = a %1 ap = ttday0_res[2] - ttdayn1_res[2] if (abs(ap) > 2): ap = ap % 1 b = ttdayp1_res[1] - ttday0_res[1] if (abs(b) > 2): b = b % 1 bp = ttdayp1_res[2] - ttday0_res[2] if abs(bp) > 2: bp = bp % 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 calculate_deltat(year, month): y = year + (month - 0.5)/12 if isinstance(y, ndarray): return np.vectorize(calculate_deltat)(year, month) if (2005 <= year) & (year < 2050): t1 = (y-2000.0) deltat = (62.92+0.32217*t1 + 0.005589*t1*t1) elif (1986 <= year) & (year < 2005): t1 = y - 2000.0 deltat = (63.86+0.3345*t1 - 0.060374*t1**2 + 0.0017275*t1**3 + 0.000651814*t1**4 + 0.00002373599*t1**5) elif (2050 <= year) & (year < 2150): deltat = (-20+32*((y-1820)/100)**2 - 0.5628*(2150-y)) elif year < -500.0: deltat = -20.0 + 32*(0.01*(y-1820.0))**2 elif (-500 <= year) & (year < 500): t1 = y/100 deltat = (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) elif (500 <= year) & (year < 1600): t1 = (y-1000)/100 deltat = (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) elif (1600 <= year) & (year < 1700): t1 = (y-1600.0) deltat = (120-0.9808*(y-1600) - 0.01532*(y-1600)**2 + (y-1600)**3/7129) elif (1700 <= year) & (year < 1800): t1 = (y - 1700.0) deltat = (8.83+0.1603*(y-1700) - 0.0059285*(y-1700)**2 + 0.00013336*(y-1700)**3 - (y-1700)**4/1174000) elif (1800 <= year) & (year < 1860): t1 = y - 1800.0 deltat = (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) elif (1860 <= year) & (year < 1900): t1 = y-1860.0 deltat = (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) elif (1900 <= year) & (year < 1920): t1 = y - 1900.0 deltat = (-2.79+1.494119*(y-1900) - 0.0598939*(y-1900)**2 + 0.0061966*(y-1900)**3 - 0.000197*(y-1900)**4) elif (1920 <= year) & (year < 1941): t1 = y - 1920.0 deltat = (21.20+0.84493*(y-1920) - 0.076100*(y-1920)**2 + 0.0020936*(y-1920)**3) elif (1941 <= year) & (year < 1961): t1 = y - 1950.0 deltat = (29.07+0.407*(y-1950) - (y-1950)**2/233 + (y-1950)**3/2547) elif (1961 <= year) & (year < 1986): t1 = y-1975 deltat = (45.45+1.067*(y-1975) - (y-1975)**2/260 - (y-1975)**3/718) elif year >= 2150: deltat = -20+32*((y-1820)/100)**2 return deltat