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