sunposition.py

# The MIT License (MIT)
#
# Copyright (c) 2021 Samuel Bear Powell
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

import os
import numpy as np
from datetime import datetime

_JIT = not os.environ.get('NUMBA_DISABLE_JIT',False)

if __name__ == '__main__':
    from argparse import ArgumentParser
    import sys
    parser = ArgumentParser(prog='sunposition',description='Compute sun position parameters given the time and location')
    parser.add_argument('--test',dest='test',action='store_true',help='Run tests')
    parser.add_argument('--version',action='version',version='%(prog)s 1.0')
    parser.add_argument('--citation',dest='cite',action='store_true',help='Print citation information')
    parser.add_argument('-t,--time',dest='t',type=str,default='now',help='"now" or date and time (UTC) in "YYYY-MM-DD hh:mm:ss.ssssss" format or a (UTC) POSIX timestamp')
    parser.add_argument('-lat,--latitude',dest='lat',type=float,default=51.48,help='latitude, in decimal degrees, positive for north')
    parser.add_argument('-lon,--longitude',dest='lon',type=float,default=0.0,help='longitude, in decimal degrees, positive for east')
    parser.add_argument('-e,--elevation',dest='elev',type=float,default=0,help='elevation, in meters')
    parser.add_argument('-T,--temperature',dest='temp',type=float,default=14.6,help='temperature, in degrees celcius')
    parser.add_argument('-p,--pressure',dest='p',type=float,default=1013.0,help='atmospheric pressure, in millibar')
    parser.add_argument('-a,--atmos_refract',dest='a',type=float,default=0.5667,help='atmospheric refraction at sunrise and sunset, in degrees')
    parser.add_argument('-dt',type=float,default=0.0,help='difference between earth\'s rotation time (TT) and universal time (UT1)')
    parser.add_argument('-r,--radians',dest='rad',action='store_true',help='Output in radians instead of degrees')
    parser.add_argument('--csv',dest='csv',action='store_true',help='Comma separated values (time,dt,lat,lon,elev,temp,pressure,az,zen,RA,dec,H)')
    parser.add_argument('--jit',dest='jit',action='store_true',help='Enable Numba acceleration (jit compilation time may overwhelm speed-up)')
    args = parser.parse_args()
    if args.cite:
        print("Implementation: Samuel Bear Powell, 2016")
        print("Algorithm:")
        print("Ibrahim Reda, Afshin Andreas, \"Solar position algorithm for solar radiation applications\", Solar Energy, Volume 76, Issue 5, 2004, Pages 577-589, ISSN 0038-092X, doi:10.1016/j.solener.2003.12.003")
        sys.exit(0)

    _JIT = args.jit

    if args.t == "now":
        args.t = datetime.utcnow()
    elif ":" in args.t and "-" in args.t:
        try:
            args.t = datetime.strptime(args.t,'%Y-%m-%d %H:%M:%S.%f') #with microseconds
        except:
            try:
                args.t = datetime.strptime(args.t,'%Y-%m-%d %H:%M:%S.') #without microseconds
            except:
                args.t = datetime.strptime(args.t,'%Y-%m-%d %H:%M:%S')
    else:
        args.t = datetime.utcfromtimestamp(int(args.t))

if _JIT:
    try:
        import scipy  # seemingly pointless but i'd advise keeping
        import numba
        from numba import jit
        _JIT = (numba.config.DISABLE_JIT == 0)

        @jit(nopython=True)
        def _polyval(p, x):
            y = 0.0
            for i,v in enumerate(p):
                y = y*x + v
            return y
    except:
        _JIT = False

if not _JIT:
    def jit(*args, **kwargs):
        #decorator that does nothing
        if len(args) == 1 and len(kwargs) == 0 and callable(args[0]):
            # called as @decorator
            return args[0]
        else:
            # called as @decorator(*args, **kwargs)
            return jit
    _polyval = np.polyval

def _calendar_time(dt):
    try:
        x = dt.year, dt.month, dt.day, dt.hour, dt.minute, dt.second, dt.microsecond
        return x
    except AttributeError:
        try:
            return _calendar_time(datetime.utcfromtimestamp(dt)) #will raise OSError if dt is not acceptable
        except:
            raise TypeError('dt must be datetime object or POSIX timestamp')

@jit(nopython=True)
def _julian_day(dt):
    """Calculate the Julian Day from a (year, month, day, hour, minute, second, microsecond) tuple"""
    # year and month numbers
    yr, mo, dy, hr, mn, sc, us = dt
    if mo <= 2:  # From paper: "if M = 1 or 2, then Y = Y - 1 and M = M + 12"
        mo += 12
        yr -= 1
    # day of the month with decimal time
    dy = dy + hr/24.0 + mn/(24.0*60.0) + sc/(24.0*60.0*60.0) + us/(24.0*60.0*60.0*1e6)
    # b is equal to 0 for the julian calendar and is equal to (2- A +
    # INT(A/4)), A = INT(Y/100), for the gregorian calendar
    a = int(yr / 100)
    b = 2 - a + int(a / 4)
    jd = int(365.25 * (yr + 4716)) + int(30.6001 * (mo + 1)) + dy + b - 1524.5
    return jd

@jit(nopython=True)
def _julian_ephemeris_day(jd, deltat):
    """Calculate the Julian Ephemeris Day from the Julian Day and delta-time = (terrestrial time - universal time) in seconds"""
    return jd + deltat / 86400.0

@jit(nopython=True)
def _julian_century(jd):
    """Caluclate the Julian Century from Julian Day or Julian Ephemeris Day"""
    return (jd - 2451545.0) / 36525.0

@jit(nopython=True)
def _julian_millennium(jc):
    """Calculate the Julian Millennium from Julian Ephemeris Century"""
    return jc / 10.0

@jit(nopython=True)
def _cos_sum(x, coeffs):
    y = np.zeros(len(coeffs))
    for i, abc in enumerate(coeffs):
        for a,b,c in abc:
            y[i] += a*np.cos(b + c*x)
    return y

# Earth Heliocentric Longitude coefficients (L0, L1, L2, L3, L4, and L5 in paper)
_EHL = (
    #L5:
    np.array([(1.0, 3.14, 0.0)]),
    #L4:
    np.array([(114.0, 3.142, 0.0), (8.0, 4.13, 6283.08), (1.0, 3.84, 12566.15)]),
    #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)]),
    #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)]),
    #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), (11.0, 0.77, 553.57),
    (10.0, 1.3, 3286.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)]),
    #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.4458, 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), (71.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)])
)

@jit(nopython=True)
def _heliocentric_longitude(jme):
    """Compute the Earth Heliocentric Longitude (L) in degrees given the Julian Ephemeris Millennium"""
    #L5, ..., L0
    Li = _cos_sum(jme, _EHL)
    L = _polyval(Li, jme) / 1e8
    L = np.rad2deg(L) % 360
    return L

#Earth Heliocentric Latitude coefficients (B0 and B1 in paper)
_EHB = (
    #B1:
    np.array([(9.0, 3.9, 5507.55), (6.0, 1.73, 5223.69)]),
    #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)])
)

@jit(nopython=True)
def _heliocentric_latitude(jme):
    """Compute the Earth Heliocentric Latitude (B) in degrees given the Julian Ephemeris Millennium"""
    Bi = _cos_sum(jme, _EHB)
    B = _polyval(Bi, jme) / 1e8
    B = np.rad2deg(B) % 360
    return B

#Earth Heliocentric Radius coefficients (R0, R1, R2, R3, R4)
_EHR = (
    #R4:
    np.array([(4.0, 2.56, 6283.08)]),
    #R3:
    np.array([(145.0, 4.273, 6283.076), (7.0, 3.92, 12566.15)]),
    #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)]),
    #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)]),
    #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), (85.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), (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)])
)

@jit(nopython=True)
def _heliocentric_radius(jme):
    """Compute the Earth Heliocentric Radius (R) in astronimical units given the Julian Ephemeris Millennium"""

    Ri = _cos_sum(jme, _EHR)
    R = _polyval(Ri, jme) / 1e8
    return R

@jit(nopython=True)
def _heliocentric_position(jme):
    """Compute the Earth Heliocentric Longitude, Latitude, and Radius given the Julian Ephemeris Millennium
        Returns (L, B, R) where L = longitude in degrees, B = latitude in degrees, and R = radius in astronimical units
    """
    return _heliocentric_longitude(jme), _heliocentric_latitude(jme), _heliocentric_radius(jme)

@jit(nopython=True)
def _geocentric_position(helio_pos):
    """Compute the geocentric latitude (Theta) and longitude (beta) (in degrees) of the sun given the earth's heliocentric position (L, B, R)"""
    L,B,R = helio_pos
    th = L + 180
    b = -B
    return (th, b)

@jit(nopython=True)
def _ecliptic_obliquity(jme, delta_epsilon):
    """Calculate the true obliquity of the ecliptic (epsilon, in degrees) given the Julian Ephemeris Millennium and the obliquity"""
    u = jme/10
    e0 = _polyval((2.45, 5.79, 27.87, 7.12, -39.05, -249.67, -51.38, 1999.25, -1.55, -4680.93, 84381.448), u)
    e = e0/3600.0 + delta_epsilon
    return e

#Nutation Longitude and Obliquity coefficients (Y)
_NLO_Y = np.array([(0.0,   0.0,   0.0,   0.0,   1.0), (-2.0,  0.0,   0.0,   2.0,   2.0), (0.0,   0.0,   0.0,   2.0,   2.0),
        (0.0,   0.0,   0.0,   0.0,   2.0), (0.0,   1.0,   0.0,   0.0,   0.0), (0.0,   0.0,   1.0,   0.0,   0.0),
        (-2.0,  1.0,   0.0,   2.0,   2.0), (0.0,   0.0,   0.0,   2.0,   1.0), (0.0,   0.0,   1.0,   2.0,   2.0),
        (-2.0,  -1.0,  0.0,   2.0,   2.0), (-2.0,  0.0,   1.0,   0.0,   0.0), (-2.0,  0.0,   0.0,   2.0,   1.0),
        (0.0,   0.0,   -1.0,  2.0,   2.0), (2.0,   0.0,   0.0,   0.0,   0.0), (0.0,   0.0,   1.0,   0.0,   1.0),
        (2.0,   0.0,   -1.0,  2.0,   2.0), (0.0,   0.0,   -1.0,  0.0,   1.0), (0.0,   0.0,   1.0,   2.0,   1.0),
        (-2.0,  0.0,   2.0,   0.0,   0.0), (0.0,   0.0,   -2.0,  2.0,   1.0), (2.0,   0.0,   0.0,   2.0,   2.0),
        (0.0,   0.0,   2.0,   2.0,   2.0), (0.0,   0.0,   2.0,   0.0,   0.0), (-2.0,  0.0,   1.0,   2.0,   2.0),
        (0.0,   0.0,   0.0,   2.0,   0.0), (-2.0,  0.0,   0.0,   2.0,   0.0), (0.0,   0.0,   -1.0,  2.0,   1.0),
        (0.0,   2.0,   0.0,   0.0,   0.0), (2.0,   0.0,   -1.0,  0.0,   1.0), (-2.0,  2.0,   0.0,   2.0,   2.0),
        (0.0,   1.0,   0.0,   0.0,   1.0), (-2.0,  0.0,   1.0,   0.0,   1.0), (0.0,   -1.0,  0.0,   0.0,   1.0),
        (0.0,   0.0,   2.0,   -2.0,  0.0), (2.0,   0.0,   -1.0,  2.0,   1.0), (2.0,   0.0,   1.0,   2.0,   2.0),
        (0.0,   1.0,   0.0,   2.0,   2.0), (-2.0,  1.0,   1.0,   0.0,   0.0), (0.0,   -1.0,  0.0,   2.0,   2.0),
        (2.0,   0.0,   0.0,   2.0,   1.0), (2.0,   0.0,   1.0,   0.0,   0.0), (-2.0,  0.0,   2.0,   2.0,   2.0),
        (-2.0,  0.0,   1.0,   2.0,   1.0), (2.0,   0.0,   -2.0,  0.0,   1.0), (2.0,   0.0,   0.0,   0.0,   1.0),
        (0.0,   -1.0,  1.0,   0.0,   0.0), (-2.0,  -1.0,  0.0,   2.0,   1.0), (-2.0,  0.0,   0.0,   0.0,   1.0),
        (0.0,   0.0,   2.0,   2.0,   1.0), (-2.0,  0.0,   2.0,   0.0,   1.0), (-2.0,  1.0,   0.0,   2.0,   1.0),
        (0.0,   0.0,   1.0,   -2.0,  0.0), (-1.0,  0.0,   1.0,   0.0,   0.0), (-2.0,  1.0,   0.0,   0.0,   0.0),
        (1.0,   0.0,   0.0,   0.0,   0.0), (0.0,   0.0,   1.0,   2.0,   0.0), (0.0,   0.0,   -2.0,  2.0,   2.0),
        (-1.0,  -1.0,  1.0,   0.0,   0.0), (0.0,   1.0,   1.0,   0.0,   0.0), (0.0,   -1.0,  1.0,   2.0,   2.0),
        (2.0,   -1.0,  -1.0,  2.0,   2.0), (0.0,   0.0,   3.0,   2.0,   2.0), (2.0,   -1.0,  0.0,   2.0,   2.0)])

#Nutation Longitude and Obliquity coefficients (a,b)
_NLO_AB = np.array([(-171996.0, -174.2), (-13187.0, -1.6), (-2274.0, -0.2), (2062.0, 0.2), (1426.0, -3.4), (712.0, 0.1),
        (-517.0, 1.2), (-386.0, -0.4), (-301.0, 0.0), (217.0, -0.5), (-158.0, 0.0), (129.0, 0.1),
        (123.0, 0.0), (63.0,  0.0), (63.0,  0.1), (-59.0, 0.0), (-58.0, -0.1), (-51.0, 0.0),
        (48.0,  0.0), (46.0,  0.0), (-38.0, 0.0), (-31.0, 0.0), (29.0,  0.0), (29.0,  0.0),
        (26.0,  0.0), (-22.0, 0.0), (21.0,  0.0), (17.0,  -0.1), (16.0,  0.0), (-16.0, 0.1),
        (-15.0, 0.0), (-13.0, 0.0), (-12.0, 0.0), (11.0,  0.0), (-10.0, 0.0), (-8.0,  0.0),
        (7.0,   0.0), (-7.0,  0.0), (-7.0,  0.0), (-7.0,  0.0), (6.0,   0.0), (6.0,   0.0),
        (6.0,   0.0), (-6.0,  0.0), (-6.0,  0.0), (5.0,   0.0), (-5.0,  0.0), (-5.0,  0.0),
        (-5.0,  0.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 Longitude and Obliquity coefficients (c,d)
_NLO_CD = np.array([(92025.0,   8.9), (5736.0,    -3.1), (977.0, -0.5), (-895.0,    0.5),
        (54.0,  -0.1), (-7.0,  0.0), (224.0, -0.6), (200.0, 0.0),
        (129.0, -0.1), (-95.0, 0.3), (0.0,   0.0), (-70.0, 0.0),
        (-53.0, 0.0), (0.0,   0.0), (-33.0, 0.0), (26.0,  0.0),
        (32.0,  0.0), (27.0,  0.0), (0.0,   0.0), (-24.0, 0.0),
        (16.0,  0.0), (13.0,  0.0), (0.0,   0.0), (-12.0, 0.0),
        (0.0,   0.0), (0.0,   0.0), (-10.0, 0.0), (0.0,   0.0),
        (-8.0,  0.0), (7.0,   0.0), (9.0,   0.0), (7.0,   0.0),
        (6.0,   0.0), (0.0,   0.0), (5.0,   0.0), (3.0,   0.0),
        (-3.0,  0.0), (0.0,   0.0), (3.0,   0.0), (3.0,   0.0),
        (0.0,   0.0), (-3.0,  0.0), (-3.0,  0.0), (3.0,   0.0),
        (3.0,   0.0), (0.0,   0.0), (3.0,   0.0), (3.0,   0.0),
        (3.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), (0.0, 0.0), (0.0, 0.0), (0.0, 0.0),
        (0.0, 0.0), (0.0, 0.0), (0.0, 0.0)])

@jit(nopython=True)
def _nutation_obliquity(jce):
    """compute the nutation in longitude (delta_psi) and the true obliquity (epsilon) given the Julian Ephemeris Century"""
    #mean elongation of the moon from the sun, in radians:
    #x0 = 297.85036 + 445267.111480*jce - 0.0019142*(jce**2) + (jce**3)/189474
    x0 = np.deg2rad(_polyval((1./189474, -0.0019142, 445267.111480, 297.85036),jce))
    #mean anomaly of the sun (Earth), in radians:
    x1 = np.deg2rad(_polyval((-1/3e5, -0.0001603, 35999.050340, 357.52772), jce))
    #mean anomaly of the moon, in radians:
    x2 = np.deg2rad(_polyval((1./56250, 0.0086972, 477198.867398, 134.96298), jce))
    #moon's argument of latitude, in radians:
    x3 = np.deg2rad(_polyval((1./327270, -0.0036825, 483202.017538, 93.27191), jce))
    #Longitude of the ascending node of the moon's mean orbit on the ecliptic
    # measured from the mean equinox of the date, in radians
    x4 = np.deg2rad(_polyval((1./45e4, 0.0020708, -1934.136261, 125.04452), jce))

    x = np.array([x0, x1, x2, x3, x4])

    a,b = _NLO_AB.T
    c,d = _NLO_CD.T
    dp = np.sum((a + b*jce)*np.sin(np.dot(_NLO_Y, x)))/36e6
    de = np.sum((c + d*jce)*np.cos(np.dot(_NLO_Y, x)))/36e6

    e = _ecliptic_obliquity(_julian_millennium(jce), de)

    return dp, e

@jit(nopython=True)
def _abberation_correction(R):
    """Calculate the abberation correction (delta_tau, in degrees) given the Earth Heliocentric Radius (in AU)"""
    return -20.4898/(3600*R)

@jit(nopython=True)
def _sun_longitude(helio_pos, delta_psi):
    """Calculate the apparent sun longitude (lambda, in degrees) and geocentric latitude (beta, in degrees) given the earth heliocentric position and delta_psi"""
    L,B,R = helio_pos
    theta = L + 180 #geocentric longitude
    beta = -B #geocentric latitude
    ll = theta + delta_psi + _abberation_correction(R)
    return ll, beta

@jit(nopython=True)
def _greenwich_sidereal_time(jd, delta_psi, epsilon):
    """Calculate the apparent Greenwich sidereal time (v, in degrees) given the Julian Day"""
    jc = _julian_century(jd)
    #mean sidereal time at greenwich, in degrees:
    v0 = (280.46061837 + 360.98564736629*(jd - 2451545) + 0.000387933*(jc**2) - (jc**3)/38710000) % 360
    v = v0 + delta_psi*np.cos(np.deg2rad(epsilon))
    return v

@jit(nopython=True)
def _sun_ra_decl(llambda, epsilon, beta):
    """Calculate the sun's geocentric right ascension (alpha, in degrees) and declination (delta, in degrees)"""
    l = np.deg2rad(llambda)
    e = np.deg2rad(epsilon)
    b = np.deg2rad(beta)
    alpha = np.arctan2(np.sin(l)*np.cos(e) - np.tan(b)*np.sin(e), np.cos(l)) #x1 / x2
    alpha = np.rad2deg(alpha) % 360
    delta = np.arcsin(np.sin(b)*np.cos(e) + np.cos(b)*np.sin(e)*np.sin(l))
    delta = np.rad2deg(delta)
    return alpha, delta

@jit(nopython=True)
def _sun_topo_ra_decl_hour(latitude, longitude, elevation, jd, delta_t = 0):
    """Calculate the sun's topocentric right ascension (alpha'), declination (delta'), and hour angle (H')"""

    jde = _julian_ephemeris_day(jd, delta_t)
    jce = _julian_century(jde)
    jme = _julian_millennium(jce)

    helio_pos = _heliocentric_position(jme)
    R = helio_pos[-1]
    phi, E = np.deg2rad(latitude), elevation
    #equatorial horizontal parallax of the sun, in radians
    xi = np.deg2rad(8.794/(3600*R)) #
    #rho = distance from center of earth in units of the equatorial radius
    #phi-prime = geocentric latitude
    #NB: These equations look like their based on WGS-84, but are rounded slightly
    # The WGS-84 reference ellipsoid has major axis a = 6378137 m, and flattening factor 1/f = 298.257223563
    # minor axis b = a*(1-f) = 6356752.3142 = 0.996647189335*a
    u = np.arctan(0.99664719*np.tan(phi)) #
    x = np.cos(u) + E*np.cos(phi)/6378140 #rho sin(phi-prime)
    y = 0.99664719*np.sin(u) + E*np.sin(phi)/6378140 #rho cos(phi-prime)

    delta_psi, epsilon = _nutation_obliquity(jce) #

    llambda, beta = _sun_longitude(helio_pos, delta_psi) #

    alpha, delta = _sun_ra_decl(llambda, epsilon, beta) #

    v = _greenwich_sidereal_time(jd, delta_psi, epsilon) #

    H = v + longitude - alpha #
    Hr, dr = np.deg2rad(H), np.deg2rad(delta)

    dar = np.arctan2(-x*np.sin(xi)*np.sin(Hr), np.cos(dr)-x*np.sin(xi)*np.cos(Hr))
    delta_alpha = np.rad2deg(dar) #

    alpha_prime = alpha + delta_alpha #
    delta_prime = np.rad2deg(np.arctan2((np.sin(dr) - y*np.sin(xi))*np.cos(dar), np.cos(dr) - y*np.sin(xi)*np.cos(Hr))) #
    H_prime = H - delta_alpha #

    return alpha_prime, delta_prime, H_prime

@jit(nopython=True)
def _sun_topo_azimuth_zenith(latitude, delta_prime, H_prime, temperature=14.6, pressure=1013, atmos_refract=0.5667):
    """Compute the sun's topocentric azimuth and zenith angles
    azimuth is measured eastward from north, zenith from vertical
    temperature = average temperature in C (default is 14.6 = global average in 2013)
    pressure = average pressure in mBar (default 1013 = global average)
    """
    SUN_RADIUS = 0.26667
    phi = np.deg2rad(latitude)
    dr, Hr = np.deg2rad(delta_prime), np.deg2rad(H_prime)
    P, T = pressure, temperature
    e0 = np.rad2deg(np.arcsin(np.sin(phi)*np.sin(dr) + np.cos(phi)*np.cos(dr)*np.cos(Hr)))
    delta_e = 0.0
    if e0 >= -1*(SUN_RADIUS + atmos_refract):
        tmp = np.deg2rad(e0 + 10.3/(e0+5.11))
        delta_e = (P/1010.0)*(283.0/(273+T))*(1.02/(60*np.tan(tmp)))

    e = e0 + delta_e
    zenith = 90 - e

    gamma = np.rad2deg(np.arctan2(np.sin(Hr), np.cos(Hr)*np.sin(phi) - np.tan(dr)*np.cos(phi))) % 360
    Phi = (gamma + 180) % 360  # azimuth from north
    return Phi, 90 - zenith, delta_e  # return elevation from horizon

@jit(nopython=True)
def _norm_lat_lon(lat,lon):
    if lat < -90 or lat > 90:
        #convert to cartesian and back
        x = np.cos(np.deg2rad(lon))*np.cos(np.deg2rad(lat))
        y = np.sin(np.deg2rad(lon))*np.cos(np.deg2rad(lat))
        z = np.sin(np.deg2rad(lat))
        r = np.sqrt(x**2 + y**2 + z**2)
        lon = np.rad2deg(np.arctan2(y,x)) % 360
        lat = np.rad2deg(np.arcsin(z/r))
    elif lon < 0 or lon > 360:
        lon = lon % 360
    return lat,lon

@jit(nopython=True)
def _topo_pos(jd,lat,lon,elev,dt,radians):
    """compute RA,dec,H, all in degrees"""
    lat,lon = _norm_lat_lon(lat,lon)
    RA, dec, H = _sun_topo_ra_decl_hour(lat, lon, elev, jd, dt)
    if radians:
        return np.deg2rad(RA), np.deg2rad(dec), np.deg2rad(H)
    else:
        return RA, dec, H

_topo_pos_v = np.vectorize(_topo_pos)

@jit(nopython=True)
def _pos(jd,lat,lon,elev,temp,press,atmos_refract,dt,radians):
    """Compute azimuth,zenith,RA,dec,H"""
    lat,lon = _norm_lat_lon(lat,lon)
    RA, dec, H = _sun_topo_ra_decl_hour(lat, lon, elev, jd, dt)
    azimuth, zenith, delta_e = _sun_topo_azimuth_zenith(lat, dec, H, temp, press, atmos_refract)
    if radians:
        return np.deg2rad(azimuth), np.deg2rad(zenith), np.deg2rad(RA), np.deg2rad(dec), np.deg2rad(H)
    else:
        return azimuth,zenith,RA,dec,H

_pos_v = np.vectorize(_pos)

@np.vectorize
def julian_day(dt):
    """Convert UTC datetimes or UTC timestamps to Julian days

    Parameters
    ----------
    dt : array_like
        UTC datetime objects or UTC timestamps (as per datetime.utcfromtimestamp)

    Returns
    -------
    jd : ndarray
        datetimes converted to fractional Julian days
    """
    t = _calendar_time(dt)
    return _julian_day(t)

@jit
def arcdist(p0,p1,radians=False):
    """Angular distance between azimuth,zenith pairs

    Parameters
    ----------
    p0 : array_like, shape (..., 2)
    p1 : array_like, shape (..., 2)
        p[...,0] = azimuth angles, p[...,1] = zenith angles
    radians : boolean (default False)
        If False, angles are in degrees, otherwise in radians

    Returns
    -------
    ad :  array_like, shape is broadcast(p0,p1).shape
        Arcdistances between corresponding pairs in p0,p1
        In degrees by default, in radians if radians=True
    """
    #formula comes from translating points into cartesian coordinates
    #taking the dot product to get the cosine between the two vectors
    #then arccos to return to angle, and simplify everything assuming real inputs
    p0,p1 = np.asarray(p0), np.asarray(p1)
    if not radians:
        p0,p1 = np.deg2rad(p0), np.deg2rad(p1)
    a0,z0 = p0[...,0], p0[...,1]
    a1,z1 = p1[...,0], p1[...,1]
    d = np.arccos(np.cos(z0)*np.cos(z1)+np.cos(a0-a1)*np.sin(z0)*np.sin(z1))
    if radians:
        return d
    else:
        return np.rad2deg(d)

def topocentric_sunpos(dt, latitude, longitude, elevation, delta_t=0, radians=False):
    """Compute the topocentric coordinates of the sun as viewed at the given time and location.

    Parameters
    ----------
    dt : array_like of datetime or float
        UTC datetime objects or UTC timestamps (as per datetime.utcfromtimestamp) representing the times of observations
    latitude, longitude : array_like of float
        decimal degrees, positive for north of the equator and east of Greenwich
    elevation : array_like of float
        meters, relative to the WGS-84 ellipsoid
    delta_t : array_like of float, optional
        seconds, default is 0, difference between the earth's rotation time (TT) and universal time (UT)
    radians : bool, optional
        return results in radians if True, degrees if False (default)

    Returns
    -------
    right_ascension : ndarray, topocentric
    declination : ndarray, topocentric
    hour_angle : ndarray, topocentric
    """

    jd = julian_day(dt)
    return _topo_pos_v(jd, latitude, longitude, elevation, delta_t, radians)

def sunpos(dt, latitude, longitude, elevation, temperature=None, pressure=None, atmos_refract=None, delta_t=0, radians=False):
    """Compute the observed and topocentric coordinates of the sun as viewed at the given time and location.

    Parameters
    ----------
    dt : array_like of datetime or float
        UTC datetime objects or UTC timestamps (as per datetime.utcfromtimestamp) representing the times of observations
    latitude, longitude : array_like of float
        decimal degrees, positive for north of the equator and east of Greenwich
    elevation : array_like of float
        meters, relative to the WGS-84 ellipsoid
    temperature : None or array_like of float, optional
        celcius, default is 14.6 (global average in 2013)
    pressure : None or array_like of float, optional
        millibar, default is 1013 (global average in ??)
    atmos_refract : None or array_like of float, optional
        Atmospheric refraction at sunrise and sunset, in degrees. Default is 0.5667
    delta_t : array_like of float, optional
        seconds, default is 0, difference between the earth's rotation time (TT) and universal time (UT)
    radians : bool, optional
        return results in radians if True, degrees if False (default)

    Returns
    -------
    azimuth_angle : ndarray, measured eastward from north
    zenith_angle : ndarray, measured down from vertical
    right_ascension : ndarray, topocentric
    declination : ndarray, topocentric
    hour_angle : ndarray, topocentric
    """

    if temperature is None:
        temperature = 14.6
    if pressure is None:
        pressure = 1013
    if atmos_refract is None:
        atmos_refract = 0.5667

    jd = julian_day(dt)
    return _pos_v(jd, latitude, longitude, elevation, temperature, pressure, atmos_refract, delta_t, radians)

def observed_sunpos(dt, latitude, longitude, elevation, temperature=None, pressure=None, atmos_refract=None, delta_t=0, radians=False):
    """Compute the observed coordinates of the sun as viewed at the given time and location.

    Parameters
    ----------
    dt : array_like of datetime or float
        UTC datetime objects or UTC timestamps (as per datetime.utcfromtimestamp) representing the times of observations
    latitude, longitude : array_like of float
        decimal degrees, positive for north of the equator and east of Greenwich
    elevation : array_like of float
        meters, relative to the WGS-84 ellipsoid
    temperature : None or array_like of float, optional
        celcius, default is 14.6 (global average in 2013)
    pressure : None or array_like of float, optional
        millibar, default is 1013 (global average in ??)
    atmos_refract : None or array_like of float, optional
        Atmospheric refraction at sunrise and sunset, in degrees. Default is 0.5667
    delta_t : array_like of float, optional
        seconds, default is 0, difference between the earth's rotation time (TT) and universal time (UT)
    radians : bool, optional
        return results in radians if True, degrees if False (default)

    Returns
    -------
    azimuth_angle : ndarray, measured eastward from north
    zenith_angle : ndarray, measured down from vertical
    """
    return sunpos(dt, latitude, longitude, elevation, temperature, pressure, atmos_refract, delta_t, radians)[:2]

def test():
    test_file = 'test_1.txt'
    # Parse and compare results from https://midcdmz.nrel.gov/solpos/spa.html
    param_names = ['syear','smonth','sday','eyear','emonth','eday','otype','step','stepunit','hr','min','sec','latitude','longitude','timezone','elev','press','temp','dut1','deltat','azmrot','slope','refract']
    param_dtype = np.dtype([(name, float) for name in param_names])
    params = np.loadtxt(test_file, param_dtype, delimiter=',', skiprows=2, max_rows=1)

    row_type = np.dtype([
            ('Date_M/D/YYYY', 'S10'),
            ('Time_H:MM:SS', 'S8'),
            ('Topo_zen', float),
            ('Topo_az', float),
            ('Julian_day', float),
            ('Julian_century', float),
            ('Julian_ephemeris_day', float),
            ('Julian_ephemeris_century', float),
            ('Julian_ephemeris_millennium', float),
            ('Earth_heliocentric_longitude', float),
            ('Earth_heliocentric_latitude', float),
            ('Earth_radius_vector', float),
            ('Geocentric_longitude', float),
            ('Geocentric_latitude', float),
            ('Mean_elongation', float),
            ('Mean_anomaly_sun', float),
            ('Mean_anomaly_moon', float),
            ('Argument_latitude_moon', float),
            ('Ascending_longitude_moon', float),
            ('Nutation_longitude', float),
            ('Nutation_obliquity', float),
            ('Ecliptic_mean_obliquity', float),
            ('Ecliptic_true_obliquity', float),
            ('Aberration_correction', float),
            ('Apparent_sun_longitude', float),
            ('Greenwich_mean_sidereal_time', float),
            ('Greenwich_sidereal_time', float),
            ('Geocentric_sun_right_ascension', float),
            ('Geocentric_sun_declination', float),
            ('Observer_hour_angle', float),
            ('Sun_equatorial_horizontal_parallax', float),
            ('Sun_right_ascension_parallax', float),
            ('Topo_sun_declination', float),
            ('Topo_sun_right_ascension', float),
            ('Topo_local_hour_angle', float),
            ('Topo_elevation_angle_uncorrected', float),
            ('Atmospheric_refraction_correction', float),
            ('Topo_elevation_angle_corrected', float),
            ('Equation_of_time', float),
            ('Sunrise_hour_angle', float),
            ('Sunset_hour_angle', float),
            ('Sun_transit_altitude', float)])

    true_data = np.loadtxt(test_file, row_type, delimiter=',', skiprows=4)

    def to_datetime(date_time_pair):
        s = str(b' '.join(date_time_pair),'UTF-8')
        return datetime.strptime(s, '%m/%d/%Y %H:%M:%S')

    def angle_diff(a1, a2, period=2*np.pi):
        """(a1 - a2 + d) % (2*d) - d; d = period/2"""
        d = period/2
        return ((a1 - a2 + d) % (period)) - d

    dts = [to_datetime(dt_pair) for dt_pair in true_data[['Date_M/D/YYYY','Time_H:MM:SS']]]
    lat,lon,elev,temp,press,deltat = params['latitude'],params['longitude'],params['elev'],params['temp'],params['press'],params['deltat']
    all_errs = []
    for dt,truth in zip(dts,true_data):
        t = _calendar_time(dt)
        jd = _julian_day(t) #Julian_day
        jde = _julian_ephemeris_day(jd, deltat) #Julian_ephemeris_day
        jce = _julian_century(jde) #Julian_ephemeris_century
        jme = _julian_millennium(jce) #Julian_ephemeris_millenium
        L,B,R = _heliocentric_position(jme) #Earth_heliocentric_longitude, Earth_heliocentric_latitude, Earth_radius_vector
        delta_psi, epsilon = _nutation_obliquity(jce) #Nutation_longitude, Ecliptic_true_obliquity
        theta,beta = _geocentric_position((L,B,R)) #Geocentric_longitude, Geocentric_latitude
        delta_tau = _abberation_correction(R) #Aberration_correction
        llambda, beta = _sun_longitude((L,B,R), delta_psi) #Apparent_sun_longitude, Geocentric_latitude (identical to previous)
        v = _greenwich_sidereal_time(jd, delta_psi, epsilon) #Greenwich_sidereal_time
        alpha, delta = _sun_ra_decl(llambda, epsilon, beta) #Geocentric_sun_right_ascension, Geocentric_sun_declination
        alpha_p, delta_p, H_p = _sun_topo_ra_decl_hour(lat,lon,elev,jd,deltat) #Topo_sun_right_ascension, Topo_sun_declination, Topo_local_hour_angle
        az, zen, delta_e = _sun_topo_azimuth_zenith(lat,delta_p,H_p,temp,press) #Topo_az, Topo_zen, Atmospheric_refraction_correction

        jd_err = jd - truth['Julian_day']
        jde_err = jde - truth['Julian_ephemeris_day']
        jce_err = jce - truth['Julian_ephemeris_century']
        jme_err = jme - truth['Julian_ephemeris_millennium']
        L_err = L - truth['Earth_heliocentric_longitude']
        B_err = B - truth['Earth_heliocentric_latitude']
        R_err = R - truth['Earth_radius_vector']
        delta_psi_err = delta_psi - truth['Nutation_longitude']
        epsilon_err = epsilon - truth['Ecliptic_true_obliquity']
        theta_err = theta - truth['Geocentric_longitude']
        beta_err = beta - truth['Geocentric_latitude']
        delta_tau_err = delta_tau - truth['Aberration_correction']
        lambda_err = llambda - truth['Apparent_sun_longitude']
        v_err = v - truth['Greenwich_sidereal_time']
        alpha_err = alpha - truth['Geocentric_sun_right_ascension']
        delta_err = delta - truth['Geocentric_sun_declination']
        alpha_prime_err = alpha_p - truth['Topo_sun_right_ascension']
        delta_prime_err = delta_p - truth['Topo_sun_declination']
        H_prime_err = angle_diff(H_p, truth['Topo_local_hour_angle'], 360)
        az_err = angle_diff(az, truth['Topo_az'], 360)
        delta_e_err = delta_e - truth['Atmospheric_refraction_correction']
        zen_err = zen - truth['Topo_zen']
        all_errs.append([jd_err,jde_err,jce_err,jme_err,L_err,B_err,R_err,delta_psi_err,
                    epsilon_err,theta_err,beta_err,delta_tau_err,lambda_err,
                    v_err,alpha_err,delta_err,alpha_prime_err,delta_prime_err,
                    H_prime_err,az_err,delta_e_err, zen_err])
    rms_err = np.sqrt(np.mean(np.array(all_errs)**2,0))
    err_names = ['Julian day', 'Julian ephemeris day', 'Julian ephemeris century', 'Julian ephemeris millennium', 'Earth heliocentric longitude', 'Earth heliocentric latitude', 'Earth radius vector', 'Nutation longitude', 'Ecliptic true obliquity', 'Geocentric longitude', 'Geocentric latitude', 'Aberration correction', 'Apparent sun longitude', 'Greenwich sidereal time', 'Geocentric sun right ascension', 'Geocentric sun declination', 'Topo sun right ascension', 'Topo sun declination', 'Topo local hour angle', 'Topo az', 'Atmospheric_refraction_correction','Topo zen']
    print('RMS Errors')
    for n, e in zip(err_names, rms_err):
        print('{}: {}'.format(n, e))

def main(args):
    az, zen, ra, dec, h = sunpos(args.t, args.lat, args.lon, args.elev, args.temp, args.p, args.dt, args.rad)
    if args.csv:
        #machine readable
        print('{t}, {dt}, {lat}, {lon}, {elev}, {temp}, {p}, {az}, {zen}, {ra}, {dec}, {h}'.format(t=args.t, dt=args.dt, lat=args.lat, lon=args.lon, elev=args.elev,temp=args.temp, p=args.p,az=az, zen=zen, ra=ra, dec=dec, h=h))
    else:
        dr='deg'
        if args.rad:
            dr='rad'
        print("Computing sun position at T = {t} + {dt} s".format(t=args.t, dt=args.dt))
        print("Lat, Lon, Elev = {lat} deg, {lon} deg, {elev} m".format(lat=args.lat, lon=args.lon, elev=args.elev))
        print("T, P = {temp} C, {press} mbar".format(temp=args.temp, press=args.p))
        print("Results:")
        print("Azimuth, zenith = {az} {dr}, {zen} {dr}".format(az=az,zen=zen,dr=dr))
        print("RA, dec, H = {ra} {dr}, {dec} {dr}, {h} {dr}".format(ra=ra, dec=dec, h=h, dr=dr))

if __name__ == '__main__':
    if args.test:
        test()
    else:
        main(args)