Skip to content

Commit

Permalink
move calculate_deltat from pvlib.spa to pvlib.solarposition
Browse files Browse the repository at this point in the history
  • Loading branch information
@kandersolar
kandersolar committed Jun 12, 2023
1 parent ef11aac commit f02266f
Show file tree
Hide file tree
Showing 3 changed files with 152 additions and 136 deletions.
2 changes: 1 addition & 1 deletion docs/sphinx/source/reference/solarposition.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@ Additional functions for quantities closely related to solar position.
solarposition.calc_time
solarposition.pyephem_earthsun_distance
solarposition.nrel_earthsun_distance
spa.calculate_deltat
solarposition.calculate_deltat


Functions for calculating sunrise, sunset and transit times.
Expand Down
157 changes: 151 additions & 6 deletions pvlib/solarposition.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -306,7 +306,7 @@ def spa_python(time, latitude, longitude,
avg. yearly air temperature in degrees C.
delta_t : float, optional, default 67.0
Difference between terrestrial time and UT1.
If delta_t is None, uses spa.calculate_deltat
If delta_t is None, uses :py:func:`calculate_deltat`
using time.year and time.month from pandas.DatetimeIndex.
For most simulations the default delta_t is sufficient.
*Note: delta_t = None will break code using nrel_numba,
Expand Down Expand Up @@ -370,7 +370,7 @@ def spa_python(time, latitude, longitude,

spa = _spa_python_import(how)

delta_t = delta_t or spa.calculate_deltat(time.year, time.month)
delta_t = delta_t or calculate_deltat(time.year, time.month)

app_zenith, zenith, app_elevation, elevation, azimuth, eot = \
spa.solar_position(unixtime, lat, lon, elev, pressure, temperature,
Expand Down Expand Up @@ -412,7 +412,7 @@ def sun_rise_set_transit_spa(times, latitude, longitude, how='numpy',
to machine code and run them multithreaded.
delta_t : float, optional, default 67.0
Difference between terrestrial time and UT1.
If delta_t is None, uses spa.calculate_deltat
If delta_t is None, uses :py:func:`calculate_deltat`
using times.year and times.month from pandas.DatetimeIndex.
For most simulations the default delta_t is sufficient.
*Note: delta_t = None will break code using nrel_numba,
Expand Down Expand Up @@ -449,7 +449,7 @@ def sun_rise_set_transit_spa(times, latitude, longitude, how='numpy',

spa = _spa_python_import(how)

delta_t = delta_t or spa.calculate_deltat(times.year, times.month)
delta_t = delta_t or calculate_deltat(times.year, times.month)

transit, sunrise, sunset = spa.transit_sunrise_sunset(
unixtime, lat, lon, delta_t, numthreads)
Expand Down Expand Up @@ -974,7 +974,7 @@ def nrel_earthsun_distance(time, how='numpy', delta_t=67.0, numthreads=4):
delta_t : float, optional, default 67.0
Difference between terrestrial time and UT1.
If delta_t is None, uses spa.calculate_deltat
If delta_t is None, uses :py:func:`calculate_deltat`
using time.year and time.month from pandas.DatetimeIndex.
For most simulations the default delta_t is sufficient.
*Note: delta_t = None will break code using nrel_numba,
Expand Down Expand Up @@ -1005,7 +1005,7 @@ def nrel_earthsun_distance(time, how='numpy', delta_t=67.0, numthreads=4):

spa = _spa_python_import(how)

delta_t = delta_t or spa.calculate_deltat(time.year, time.month)
delta_t = delta_t or calculate_deltat(time.year, time.month)

dist = spa.earthsun_distance(unixtime, delta_t, numthreads)

Expand Down Expand Up @@ -1476,3 +1476,148 @@ def sun_rise_set_transit_geometric(times, latitude, longitude, declination,
sunset = _local_times_from_hours_since_midnight(times, sunset_hour)
transit = _local_times_from_hours_since_midnight(times, transit_hour)
return sunrise, sunset, transit


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.
Parameters
----------
year : numeric
Calendar year for which to calculate the time offset
month : numeric
Calendar month (1-12) for which to calculate the time offset
Returns
-------
deltat : numeric
References
----------
.. [1] `NASA GSFC: Polynomial Expressions for Delta T (ΔT)
<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
129 changes: 0 additions & 129 deletions pvlib/spa.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1246,132 +1246,3 @@ def earthsun_distance(unixtime, delta_t, numthreads):
return R


def calculate_deltat(year, month):
"""Calculate the difference between Terrestrial Dynamical Time (TD)
and Universal Time (UT).
Note: This function is not yet compatible for calculations using
Numba.
Equations taken from http://eclipse.gsfc.nasa.gov/SEcat5/deltatpoly.html
"""

plw = 'Deltat is unknown for years before -1999 and after 3000. ' \
'Delta values will be calculated, but the calculations ' \
'are not intended to be used for these years.'

try:
if np.any((year > 3000) | (year < -1999)):
warnings.warn(plw)
except ValueError:
if (year > 3000) | (year < -1999):
warnings.warn(plw)
except TypeError:
return 0

y = year + (month - 0.5)/12

deltat = np.where(year < -500,

-20+32*((y-1820)/100)**2, 0)

deltat = np.where((-500 <= year) & (year < 500),

10583.6-1014.41*(y/100)
+ 33.78311*(y/100)**2
- 5.952053*(y/100)**3
- 0.1798452*(y/100)**4
+ 0.022174192*(y/100)**5
+ 0.0090316521*(y/100)**6, deltat)

deltat = np.where((500 <= year) & (year < 1600),

1574.2-556.01*((y-1000)/100)
+ 71.23472*((y-1000)/100)**2
+ 0.319781*((y-1000)/100)**3
- 0.8503463*((y-1000)/100)**4
- 0.005050998*((y-1000)/100)**5
+ 0.0083572073*((y-1000)/100)**6, deltat)

deltat = np.where((1600 <= year) & (year < 1700),

120-0.9808*(y-1600)
- 0.01532*(y-1600)**2
+ (y-1600)**3/7129, deltat)

deltat = np.where((1700 <= year) & (year < 1800),

8.83+0.1603*(y-1700)
- 0.0059285*(y-1700)**2
+ 0.00013336*(y-1700)**3
- (y-1700)**4/1174000, deltat)

deltat = np.where((1800 <= year) & (year < 1860),

13.72-0.332447*(y-1800)
+ 0.0068612*(y-1800)**2
+ 0.0041116*(y-1800)**3
- 0.00037436*(y-1800)**4
+ 0.0000121272*(y-1800)**5
- 0.0000001699*(y-1800)**6
+ 0.000000000875*(y-1800)**7, deltat)

deltat = np.where((1860 <= year) & (year < 1900),

7.62+0.5737*(y-1860)
- 0.251754*(y-1860)**2
+ 0.01680668*(y-1860)**3
- 0.0004473624*(y-1860)**4
+ (y-1860)**5/233174, deltat)

deltat = np.where((1900 <= year) & (year < 1920),

-2.79+1.494119*(y-1900)
- 0.0598939*(y-1900)**2
+ 0.0061966*(y-1900)**3
- 0.000197*(y-1900)**4, deltat)

deltat = np.where((1920 <= year) & (year < 1941),

21.20+0.84493*(y-1920)
- 0.076100*(y-1920)**2
+ 0.0020936*(y-1920)**3, deltat)

deltat = np.where((1941 <= year) & (year < 1961),

29.07+0.407*(y-1950)
- (y-1950)**2/233
+ (y-1950)**3/2547, deltat)

deltat = np.where((1961 <= year) & (year < 1986),

45.45+1.067*(y-1975)
- (y-1975)**2/260
- (y-1975)**3/718, deltat)

deltat = np.where((1986 <= year) & (year < 2005),

63.86+0.3345*(y-2000)
- 0.060374*(y-2000)**2
+ 0.0017275*(y-2000)**3
+ 0.000651814*(y-2000)**4
+ 0.00002373599*(y-2000)**5, deltat)

deltat = np.where((2005 <= year) & (year < 2050),

62.92+0.32217*(y-2000)
+ 0.005589*(y-2000)**2, deltat)

deltat = np.where((2050 <= year) & (year < 2150),

-20+32*((y-1820)/100)**2
- 0.5628*(2150-y), deltat)

deltat = np.where(year >= 2150,

-20+32*((y-1820)/100)**2, deltat)

deltat = deltat.item() if np.isscalar(year) & np.isscalar(month)\
else deltat

return deltat

0 comments on commit f02266f

Please sign in to comment.