"""
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 math import acos, asin, atan, atan2, cos, degrees, radians, sin, tan

from fluids.constants import deg2rad, rad2deg
from fluids.numerics import sincos

__all__ = ['julian_day_dt', 'julian_day', 'julian_ephemeris_day', 'julian_century',
           'julian_ephemeris_century', 'julian_ephemeris_millennium', 'heliocentric_longitude',
           'heliocentric_latitude', 'heliocentric_radius_vector', 'geocentric_longitude',
           'geocentric_latitude', 'mean_elongation', 'mean_anomaly_sun', 'mean_anomaly_moon',
           'moon_argument_latitude', 'moon_ascending_longitude', 'longitude_nutation',
           'obliquity_nutation', 'mean_ecliptic_obliquity', 'true_ecliptic_obliquity',
           'aberration_correction', 'apparent_sun_longitude', 'mean_sidereal_time',
           'apparent_sidereal_time', 'geocentric_sun_right_ascension', 'geocentric_sun_declination',
           'local_hour_angle', 'equatorial_horizontal_parallax', 'uterm', 'xterm', 'yterm',
           'parallax_sun_right_ascension', 'topocentric_sun_right_ascension', 'topocentric_sun_declination',
           'topocentric_local_hour_angle', 'topocentric_elevation_angle_without_atmosphere',
           'atmospheric_refraction_correction', 'topocentric_elevation_angle', 'topocentric_zenith_angle',
           'topocentric_astronomers_azimuth', 'topocentric_azimuth_angle', 'sun_mean_longitude',
           'equation_of_time', 'calculate_deltat', 'longitude_obliquity_nutation',
           'transit_sunrise_sunset',
           ]
nan = float("nan")


HELIO_RADIUS_TABLE_LIST_0 = [[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]]

HELIO_RADIUS_TABLE_LIST_1 = [[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],
]
HELIO_RADIUS_TABLE_LIST_2 = [[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]]
HELIO_RADIUS_TABLE_LIST_3 = [[145.0, 4.273, 6283.076],
 [7.0, 3.92, 12566.15]]
HELIO_RADIUS_TABLE_LIST_4 = [[4.0, 2.56, 6283.08]]

NUTATION_YTERM_LIST_0 = [0.0, -2.0, 0.0, 0.0, 0.0, 0.0, -2.0, 0.0, 0.0, -2.0, -2.0, -2.0, 0.0, 2.0, 0.0, 2.0, 0.0, 0.0, -2.0, 0.0, 2.0, 0.0, 0.0, -2.0, 0.0, -2.0, 0.0, 0.0, 2.0, -2.0, 0.0, -2.0, 0.0, 0.0, 2.0, 2.0, 0.0, -2.0, 0.0, 2.0, 2.0, -2.0, -2.0, 2.0, 2.0, 0.0, -2.0, -2.0, 0.0, -2.0, -2.0, 0.0, -1.0, -2.0, 1.0, 0.0, 0.0, -1.0, 0.0, 0.0, 2.0, 0.0, 2.0]
NUTATION_YTERM_LIST_1 = [0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.0, 0.0, 2.0, 1.0, 0.0, -1.0, 0.0, 0.0, 0.0, 1.0, 1.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0, -1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, -1.0, 1.0, -1.0, -1.0, 0.0, -1.0]
NUTATION_YTERM_LIST_2 = [0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, -1.0, 0.0, 1.0, -1.0, -1.0, 1.0, 2.0, -2.0, 0.0, 2.0, 2.0, 1.0, 0.0, 0.0, -1.0, 0.0, -1.0, 0.0, 0.0, 1.0, 0.0, 2.0, -1.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 2.0, 1.0, -2.0, 0.0, 1.0, 0.0, 0.0, 2.0, 2.0, 0.0, 1.0, 1.0, 0.0, 0.0, 1.0, -2.0, 1.0, 1.0, 1.0, -1.0, 3.0, 0.0]
NUTATION_YTERM_LIST_3 = [0.0, 2.0, 2.0, 0.0, 0.0, 0.0, 2.0, 2.0, 2.0, 2.0, 0.0, 2.0, 2.0, 0.0, 0.0, 2.0, 0.0, 2.0, 0.0, 2.0, 2.0, 2.0, 0.0, 2.0, 2.0, 2.0, 2.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, -2.0, 2.0, 2.0, 2.0, 0.0, 2.0, 2.0, 0.0, 2.0, 2.0, 0.0, 0.0, 0.0, 2.0, 0.0, 2.0, 0.0, 2.0, -2.0, 0.0, 0.0, 0.0, 2.0, 2.0, 0.0, 0.0, 2.0, 2.0, 2.0, 2.0]
NUTATION_YTERM_LIST_4 = [1.0, 2.0, 2.0, 2.0, 0.0, 0.0, 2.0, 1.0, 2.0, 2.0, 0.0, 1.0, 2.0, 0.0, 1.0, 2.0, 1.0, 1.0, 0.0, 1.0, 2.0, 2.0, 0.0, 2.0, 0.0, 0.0, 1.0, 0.0, 1.0, 2.0, 1.0, 1.0, 1.0, 0.0, 1.0, 2.0, 2.0, 0.0, 2.0, 1.0, 0.0, 2.0, 1.0, 1.0, 1.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 2.0, 2.0, 2.0, 2.0]

NUTATION_ABCD_LIST = [[-171996.0, -174.2, 92025.0, 8.9],
 [-13187.0, -1.6, 5736.0, -3.1],
 [-2274.0, -0.2, 977.0, -0.5],
 [2062.0, 0.2, -895.0, 0.5],
 [1426.0, -3.4, 54.0, -0.1],
 [712.0, 0.1, -7.0, 0.0],
 [-517.0, 1.2, 224.0, -0.6],
 [-386.0, -0.4, 200.0, 0.0],
 [-301.0, 0.0, 129.0, -0.1],
 [217.0, -0.5, -95.0, 0.3],
 [-158.0, 0.0, 0.0, 0.0],
 [129.0, 0.1, -70.0, 0.0],
 [123.0, 0.0, -53.0, 0.0],
 [63.0, 0.0, 0.0, 0.0],
 [63.0, 0.1, -33.0, 0.0],
 [-59.0, 0.0, 26.0, 0.0],
 [-58.0, -0.1, 32.0, 0.0],
 [-51.0, 0.0, 27.0, 0.0],
 [48.0, 0.0, 0.0, 0.0],
 [46.0, 0.0, -24.0, 0.0],
 [-38.0, 0.0, 16.0, 0.0],
 [-31.0, 0.0, 13.0, 0.0],
 [29.0, 0.0, 0.0, 0.0],
 [29.0, 0.0, -12.0, 0.0],
 [26.0, 0.0, 0.0, 0.0],
 [-22.0, 0.0, 0.0, 0.0],
 [21.0, 0.0, -10.0, 0.0],
 [17.0, -0.1, 0.0, 0.0],
 [16.0, 0.0, -8.0, 0.0],
 [-16.0, 0.1, 7.0, 0.0],
 [-15.0, 0.0, 9.0, 0.0],
 [-13.0, 0.0, 7.0, 0.0],
 [-12.0, 0.0, 6.0, 0.0],
 [11.0, 0.0, 0.0, 0.0],
 [-10.0, 0.0, 5.0, 0.0],
 [-8.0, 0.0, 3.0, 0.0],
 [7.0, 0.0, -3.0, 0.0],
 [-7.0, 0.0, 0.0, 0.0],
 [-7.0, 0.0, 3.0, 0.0],
 [-7.0, 0.0, 3.0, 0.0],
 [6.0, 0.0, 0.0, 0.0],
 [6.0, 0.0, -3.0, 0.0],
 [6.0, 0.0, -3.0, 0.0],
 [-6.0, 0.0, 3.0, 0.0],
 [-6.0, 0.0, 3.0, 0.0],
 [5.0, 0.0, 0.0, 0.0],
 [-5.0, 0.0, 3.0, 0.0],
 [-5.0, 0.0, 3.0, 0.0],
 [-5.0, 0.0, 3.0, 0.0],
 [4.0, 0.0, 0.0, 0.0],
 [4.0, 0.0, 0.0, 0.0],
 [4.0, 0.0, 0.0, 0.0],
 [-4.0, 0.0, 0.0, 0.0],
 [-4.0, 0.0, 0.0, 0.0],
 [-4.0, 0.0, 0.0, 0.0],
 [3.0, 0.0, 0.0, 0.0],
 [-3.0, 0.0, 0.0, 0.0],
 [-3.0, 0.0, 0.0, 0.0],
 [-3.0, 0.0, 0.0, 0.0],
 [-3.0, 0.0, 0.0, 0.0],
 [-3.0, 0.0, 0.0, 0.0],
 [-3.0, 0.0, 0.0, 0.0],
 [-3.0, 0.0, 0.0, 0.0]]

HELIO_LAT_TABLE_LIST_0 = [[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]]

HELIO_LAT_TABLE_LIST_1 = [[9.0, 3.9, 5507.55],
 [6.0, 1.73, 5223.69]]

#HELIO_LONG_TABLE_LIST = HELIO_LONG_TABLE.tolist()
HELIO_LONG_TABLE_LIST_0 = [[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]]
HELIO_LONG_TABLE_LIST_1 = [[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],
 ]
HELIO_LONG_TABLE_LIST_2 = [[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]]

HELIO_LONG_TABLE_LIST_3 = [[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]
 ]
HELIO_LONG_TABLE_LIST_4 = [[114.0, 3.142, 0.0],
 [8.0, 4.13, 6283.08],
 [1.0, 3.84, 12566.15]]



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.
    """
    # Not used anywhere!
    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*1.1574074074074073e-05 + 2440587.5
#    jd = unixtime/86400.0 + 2440587.5
    return jd


def julian_ephemeris_day(julian_day, delta_t):
    jde = julian_day + delta_t*1.1574074074074073e-05
#    jde = julian_day + delta_t * 1.0 / 86400.0
    return jde


def julian_century(julian_day):
    jc = (julian_day - 2451545.0)*2.7378507871321012e-05# * 1.0 / 36525
    return jc


def julian_ephemeris_century(julian_ephemeris_day):
#    1/36525.0 =  2.7378507871321012e-05
    jce = (julian_ephemeris_day - 2451545.0)*2.7378507871321012e-05
    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
    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]))
    l5 = -0.9999987317275395
    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
    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)
               )
    HELIO_LAT_TABLE_LIST_1_ROW = HELIO_LAT_TABLE_LIST_1[0]
    b1 += (HELIO_LAT_TABLE_LIST_1_ROW[0]
           * cos(HELIO_LAT_TABLE_LIST_1_ROW[1]
                    + HELIO_LAT_TABLE_LIST_1_ROW[2] * jme))

    HELIO_LAT_TABLE_LIST_1_ROW = HELIO_LAT_TABLE_LIST_1[1]
    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
    # Would be possible to save a few multiplies of table1row[2]*jme, table1row[1]*jme as they are dups
    for row in range(40):
        table0row = HELIO_RADIUS_TABLE_LIST_0[row]
        r0 += (table0row[0]*cos(table0row[1] + table0row[2]*jme))
    for row in range(10):
        table1row = HELIO_RADIUS_TABLE_LIST_1[row]
        r1 += (table1row[0]*cos(table1row[1] + table1row[2]*jme))
    for row in range(6):
        table2row = HELIO_RADIUS_TABLE_LIST_2[row]
        r2 += (table2row[0]*cos(table2row[1] + table2row[2]*jme))

    table3row = HELIO_RADIUS_TABLE_LIST_3[0]
    r3 += (table3row[0]*cos(table3row[1] + table3row[2]*jme))
    table3row = HELIO_RADIUS_TABLE_LIST_3[1]
    r3 += (table3row[0]*cos(table3row[1] + table3row[2]*jme))

#    table4row = HELIO_RADIUS_TABLE_LIST_4[0]
#    r4 = (table4row[0]*cos(table4row[1] + table4row[2]*jme))
    r4 = (4.0*cos(2.56 + 6283.08*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):
    return (julian_ephemeris_century*(julian_ephemeris_century
            *(5.27776898149614e-6*julian_ephemeris_century - 0.0019142)
            + 445267.11148) + 297.85036)
#    x0 = (297.85036
#          + 445267.111480 * julian_ephemeris_century
#          - 0.0019142 * julian_ephemeris_century**2
#          + julian_ephemeris_century**3 / 189474.0)
#    return x0


def mean_anomaly_sun(julian_ephemeris_century):
    return (julian_ephemeris_century*(julian_ephemeris_century*(
            -3.33333333333333e-6*julian_ephemeris_century - 0.0001603)
            + 35999.05034) + 357.52772)
#    x1 = (357.52772
#          + 35999.050340 * julian_ephemeris_century
#          - 0.0001603 * julian_ephemeris_century**2
#          - julian_ephemeris_century**3 / 300000.0)
#    return x1


def mean_anomaly_moon(julian_ephemeris_century):
    return (julian_ephemeris_century*(julian_ephemeris_century*(
            1.77777777777778e-5*julian_ephemeris_century + 0.0086972)
        + 477198.867398) + 134.96298)
#    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):
    return julian_ephemeris_century*(julian_ephemeris_century*(
            3.05558101873071e-6*julian_ephemeris_century - 0.0036825)
        + 483202.017538) + 93.27191
#    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):
    return (julian_ephemeris_century*(julian_ephemeris_century*(
            2.22222222222222e-6*julian_ephemeris_century + 0.0020708)
            - 1934.136261) + 125.04452)
#    x4 = (125.04452
#          - 1934.136261 * julian_ephemeris_century
#          + 0.0020708 * julian_ephemeris_century**2
#          + julian_ephemeris_century**3 / 450000)
#    return x4


def longitude_obliquity_nutation(julian_ephemeris_century, x0, x1, x2, x3, x4):
    x0, x1, x2, x3, x4 = deg2rad*x0, deg2rad*x1, deg2rad*x2, deg2rad*x3, deg2rad*x4
    delta_psi_sum = 0.0
    delta_eps_sum = 0.0
    # If the sincos formulation is used, the speed up is ~8% with numba.
    for row in range(63):
        arg = (NUTATION_YTERM_LIST_0[row]*x0 +
               NUTATION_YTERM_LIST_1[row]*x1 +
               NUTATION_YTERM_LIST_2[row]*x2 +
               NUTATION_YTERM_LIST_3[row]*x3 +
               NUTATION_YTERM_LIST_4[row]*x4)
        arr = NUTATION_ABCD_LIST[row]
        sinarg, cosarg = sincos(arg)
#        sinarg = sin(arg)
#        cosarg = sqrt(1.0 - sinarg*sinarg)
        t0 = (arr[0] + julian_ephemeris_century*arr[1])
        delta_psi_sum += t0*sinarg
#        delta_psi_sum += t0*sin(arg)
        t0 = (arr[2] + julian_ephemeris_century*arr[3])
        delta_eps_sum += t0*cosarg
#        delta_eps_sum += t0*cos(arg)
    delta_psi = delta_psi_sum/36000000.0
    delta_eps = delta_eps_sum/36000000.0
    res = [0.0]*2
    res[0] = delta_psi
    res[1] = delta_eps
    return res

def longitude_nutation(julian_ephemeris_century, x0, x1, x2, x3, x4):
    x0, x1, x2, x3, x4 = deg2rad*x0, deg2rad*x1, deg2rad*x2, deg2rad*x3, deg2rad*x4
    delta_psi_sum = 0.0
    for row in range(63):
#       # None can be skipped but the multiplies can be with effort -2 to 2 with dict - just might be slower
        argsin = (NUTATION_YTERM_LIST_0[row]*x0 +
                  NUTATION_YTERM_LIST_1[row]*x1 +
                  NUTATION_YTERM_LIST_2[row]*x2 +
                  NUTATION_YTERM_LIST_3[row]*x3 +
                  NUTATION_YTERM_LIST_4[row]*x4)
        term = (NUTATION_ABCD_LIST[row][0] + NUTATION_ABCD_LIST[row][1]
                * 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
    x0, x1, x2, x3, x4 = deg2rad*x0, deg2rad*x1, deg2rad*x2, deg2rad*x3, deg2rad*x4
    for row in range(63):
        argcos = (NUTATION_YTERM_LIST_0[row]*x0 +
                  NUTATION_YTERM_LIST_1[row]*x1 +
                  NUTATION_YTERM_LIST_2[row]*x2 +
                  NUTATION_YTERM_LIST_3[row]*x3 +
                  NUTATION_YTERM_LIST_4[row]*x4)
        term = (NUTATION_ABCD_LIST[row][2]
               + NUTATION_ABCD_LIST[row][3]*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
    return mean_ecliptic_obliquity*0.0002777777777777778  + 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)
          + 0.000387933*julian_century2
          - julian_century2*julian_century/38710000.0)
    return v0 % 360.0


def apparent_sidereal_time(mean_sidereal_time, longitude_nutation,
                           true_ecliptic_obliquity):
    v = mean_sidereal_time + longitude_nutation*cos(deg2rad*true_ecliptic_obliquity)
    return v


def geocentric_sun_right_ascension(apparent_sun_longitude,
                                   true_ecliptic_obliquity,
                                   geocentric_latitude):
    num = (sin(deg2rad*apparent_sun_longitude)
           * cos(deg2rad*true_ecliptic_obliquity)
           - tan(deg2rad*geocentric_latitude)
           * sin(deg2rad*true_ecliptic_obliquity))
    alpha = degrees(atan2(num, cos(
        deg2rad*apparent_sun_longitude)))
    return alpha % 360


def geocentric_sun_declination(apparent_sun_longitude, true_ecliptic_obliquity,
                               geocentric_latitude):
    delta = degrees(asin(sin(deg2rad*geocentric_latitude) *
                                 cos(deg2rad*true_ecliptic_obliquity) +
                                 cos(deg2rad*geocentric_latitude) *
                                 sin(deg2rad*true_ecliptic_obliquity) *
                                 sin(deg2rad*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):
    return 8.794 / (3600.0 * earth_radius_vector)


def uterm(observer_latitude):
    u = atan(0.99664719*tan(deg2rad*observer_latitude))
    return u


def xterm(u, observer_latitude, observer_elevation):
    # 1/6378140.0 = const
    x = (cos(u) + observer_elevation*1.5678552054360676e-07*cos(deg2rad*observer_latitude))
    return x


def yterm(u, observer_latitude, observer_elevation):
    # 1/6378140.0 = const
    y = (0.99664719 * sin(u) + observer_elevation*1.5678552054360676e-07
         * sin(deg2rad*observer_latitude))
    return y


def parallax_sun_right_ascension(xterm, equatorial_horizontal_parallax,
                                 local_hour_angle, geocentric_sun_declination):
    x0 = sin(deg2rad*equatorial_horizontal_parallax)
    x1 = deg2rad*local_hour_angle
    num = -xterm*x0*sin(x1)
    denom = (cos(deg2rad*geocentric_sun_declination) - xterm*x0 * cos(x1))
    delta_alpha = degrees(atan2(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):
    x0 = sin(deg2rad*equatorial_horizontal_parallax)
    num = ((sin(deg2rad*geocentric_sun_declination) - yterm
            * x0)
           * cos(deg2rad*parallax_sun_right_ascension))
    denom = (cos(deg2rad*geocentric_sun_declination) - xterm
             * x0
             * cos(deg2rad*local_hour_angle))
    delta = degrees(atan2(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
                                                   ):

    r_observer_latitude = deg2rad*observer_latitude
    r_topocentric_sun_declination = deg2rad*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(deg2rad*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)
    delta_e = ((local_pressure / 1010.0) * (283.0 / (273.0 + local_temp))
               * 1.02 / (60.0 * tan(deg2rad*(
                   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):
    num = sin(deg2rad*topocentric_local_hour_angle)
    denom = (cos(deg2rad*topocentric_local_hour_angle)
             * sin(deg2rad*observer_latitude)
             - tan(deg2rad*topocentric_sun_declination)
             * cos(deg2rad*observer_latitude))
    gamma = degrees(atan2(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):
    term = cos(deg2rad*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):
    """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.

    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.
    """
    jd = julian_day(unixtime)
    jde = julian_ephemeris_day(jd, delta_t)
    jce = julian_ephemeris_century(jde)
    jme = julian_ephemeris_millennium(jce)
    R = heliocentric_radius_vector(jme)
    return R


def solar_position(unixtime, lat, lon, elev, pressure, temp, delta_t,
                   atmos_refract, sst=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 specifying 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.

    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)
    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, delta_epsilon = longitude_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: # numba: delete
        return v, alpha, delta # numba: delete

    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]


try:
    if IS_NUMBA:  # type: ignore # noqa: F821
        import threading

        import numba
        import numpy as np
        # This is 3x slower without nogil
        @numba.njit(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(len(unixtime)):
                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)
                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)
                delta_psi, delta_epsilon = longitude_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:
                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


except:
    pass


def transit_sunrise_sunset(dates, lat, lon, delta_t):
    """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.

    Returns
    -------
    tuple : (transit, sunrise, sunset) localized to UTC

    >>> transit_sunrise_sunset(1523836800, 51.0486, -114.07, 70.68302220312503)
    (1523907360.3863413, 1523882341.570479, 1523932345.7781625)
    """
    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.0
    ttdayp1 = ttday0 + 86400.0

    # index 0 is v, 1 is alpha, 2 is delta
    utday_res = solar_position(utday, 0, 0, 0, 0, 0, delta_t,
                               0, sst=True)
    v = utday_res[0]

    ttday0_res = solar_position(ttday0, 0, 0, 0, 0, 0, delta_t,
                                0, sst=True)
    ttdayn1_res = solar_position(ttdayn1, 0, 0, 0, 0, 0, delta_t,
                                 0, sst=True)
    ttdayp1_res = solar_position(ttdayp1, 0, 0, 0, 0, 0, delta_t,
                                 0, sst=True)
    m0 = (ttday0_res[1] - lon - v) / 360
    cos_arg = ((-0.014543315936696236 - sin(radians(lat)) # sin(radians(-0.8333)) = -0.0145...
               * sin(radians(ttday0_res[2]))) /
               (cos(radians(lat)) * cos(radians(ttday0_res[2]))))
    if abs(cos_arg) > 1:
        cos_arg = nan

    H0 = degrees(acos(cos_arg)) % 180

    m = [0.0]*3
    m[0] = m0 % 1
    m[1] = (m[0] - H0 / 360.0)
    m[2] = (m[0] + H0 / 360.0)

    # 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 = [0.0]*3
    for i in range(3):
        vs[i] = v + 360.985647*m[i]
    n = [0.0]*3
    for i in range(3):
        n[i] = m[i] + delta_t / 86400.0

    a = ttday0_res[1] - ttdayn1_res[1]

    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 = [0.0]*3
    delta_prime = [0.0]*3
    Hp = [0.0]*3
    for i in range(3):
        alpha_prime[i] = ttday0_res[1] + (n[i] * (a + b + c * n[i]))*0.5
        delta_prime[i] = ttday0_res[2] + (n[i] * (ap + bp + cp * n[i]))*0.5
        Hp[i] = (vs[i] + lon - alpha_prime[i]) % 360
        if Hp[i] >= 180.0:
            Hp[i] = Hp[i] - 360.0


    #alpha_prime = ttday0_res[1] + (n * (a + b + c * n)) / 2 # this is vect
    #delta_prime = ttday0_res[2] + (n * (ap + bp + cp * n)) / 2 # this is vect
    #Hp = (vs + lon - alpha_prime) % 360
    #Hp[Hp >= 180] = Hp[Hp >= 180] - 360
    x1 = sin(radians(lat))
    x2 = cos(radians(lat))

    h = [0.0]*3
    for i in range(3):
        h[i] = degrees(asin(x1*sin(radians(delta_prime[i])) + x2 * cos(radians(delta_prime[i])) * cos(radians(Hp[i]))))

    T = float((m[0] - Hp[0] / 360.0) * 86400.0)
    R = float((m[1] + (h[1] + 0.8333) / (360.0 * cos(radians(delta_prime[1])) *
                                   cos(radians(lat)) *
                                   sin(radians(Hp[1])))) * 86400.0)
    S = float((m[2] + (h[2] + 0.8333) / (360.0 * cos(radians(delta_prime[2])) *
                                   cos(radians(lat)) *
                                   sin(radians(Hp[2])))) * 86400.0)

    if add_a_day:
        S += 86400.0
    if sub_a_day:
        R -= 86400.0

    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 (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