In [72]:
%matplotlib inline
from __future__ import (absolute_import, division, print_function,
unicode_literals)
import numpy as np
from scipy.optimize import leastsq, fmin_l_bfgs_b
import matplotlib
import matplotlib.pyplot as plt
import datetime
import astropy.units as u
from astropy.coordinates import (EarthLocation, Latitude,
Longitude, get_sun)
from astropy.time import Time
from astropy.coordinates import AltAz
##################################################
# IERS table patch from @eteq ####################
from astropy.utils.data import download_file
from astropy.utils import iers
iers.IERS.iers_table = iers.IERS_A.open(download_file(iers.IERS_A_URL,
cache=True))
##################################################
# Some convenience functions
def get_altitude(target, time, location, kwargs={}):
'''
Compute altitude of a `target` at a particular time, given
an observer's `location` (EarthLocation).
'''
altaz = target.transform_to(AltAz(target.ra, target.dec,
location=location, obstime=time,
**kwargs))
return altaz.alt
def timebracket(t0, dt, N):
'''
Generate `N` times bracketing `t0` going
backwards and forwards `dt` hours
t0 : datetime.datetime
Initial time
dt : float
Number of hours before and after t0
N : int
Number of datetimes in range
'''
previous_dt = t0 + np.linspace(-dt, 0, N)*u.hour
next_dt = t0 + np.linspace(0, dt, N)*u.hour
return previous_dt, next_dt
def timelinspace(t0, dt, N):
'''
Generate `N` linearly spaced times between `t0`
and t0+dt where [dt] = day
t0 : datetime.datetime
Initial time
dt : float
Number of days before/after t0.
dt < 0 = before, dt > 0 = after
N : int
Number of times to generate in range
'''
# Keep array in time order
if dt > 0:
time_array = t0 + np.linspace(0, dt, N)*u.day
elif dt < 0:
time_array = t0 + np.linspace(dt, 0, N)*u.day
return time_array
def horiz_cross_rising(t, a, horizon=0, unit=u.degree,
return_bounds=False, return_a=False):
'''
Find time `t` when values in array `a` go from
negative to positive (exclude endpoints)
`return_nearest` will return nearest times to zero-crossing
'''
condition = (a[:-1] < horizon*unit) * (a[1:] > horizon*unit)
if sum(condition) < 1:
return []
if return_bounds:
nearest_index = np.argwhere(condition)[0][0]
if a[nearest_index] > horizon*unit:
lower_bound, upper_bound = t[nearest_index-1], t[nearest_index]
lower_a, upper_a = a[nearest_index-1].value, a[nearest_index].value
else:
lower_bound, upper_bound = t[nearest_index], t[nearest_index+1]
lower_a, upper_a = a[nearest_index].value, a[nearest_index+1].value
if return_a:
return (lower_bound, upper_bound), (lower_a, upper_a)
else:
return lower_bound, upper_bound
return t[condition][0]
def horiz_cross_setting(t, a, horizon=0, unit=u.degree, return_bounds=False):
'''
Find time `t` when values in array `a` go from
positive to negative (exclude endpoints)
'''
condition = (a[:-1] > horizon*unit) * (a[1:] < horizon*unit)
if sum(condition) < 1:
return []
if return_bounds:
nearest_index = np.argwhere(condition)[0][0]
if a[nearest_index] < horizon*unit:
lower_bound, upper_bound = t[nearest_index-1], t[nearest_index]
else:
lower_bound, upper_bound = t[nearest_index], t[nearest_index+1]
return lower_bound, upper_bound
return t[condition][0]
def get_rise(time, dt, location, N=1000, return_bounds=False, return_a=False):
times = timelinspace(time, dt, N)
sun = get_sun(times)
altitudes = get_altitude(sun, times, location)
return horiz_cross_rising(times, altitudes,
return_bounds=return_bounds,
return_a=return_a)
# A variety of ways to compute sunrise
def sunrise_doublegrid(time, location, which='nearest', N=300):
if which == 'nearest' or which == 'next':
# Calculate next sunrise
# Do a coarse search with error O(~5 minutes),
# return upper/lower bounds
rise_bounds = get_rise(time, 1, location, N=N,
return_bounds=True)
# Do a fine search with error O(~1 second),
# return the nearest time to sunrise
next_rise = get_rise(rise_bounds[0],
(rise_bounds[1] - rise_bounds[0]).value,
location, N=N)
if which == 'next':
return next_rise
if which == 'nearest' or which == 'previous':
# Calculate previous sunrise
# Do a coarse search with error O(~5 minutes),
# return upper/lower bounds
rise_bounds = get_rise(time, -1, location, N=N,
return_bounds=True)
# Do a fine search with error O(~1 second),
# return the nearest time to sunrise
previous_rise = get_rise(rise_bounds[0],
(rise_bounds[1] - rise_bounds[0]).value,
location, N=N)
if which == 'previous':
return previous_rise
if which == 'nearest':
if abs(previous_rise - time) < abs(next_rise - time):
return previous_rise
else:
return next_rise
raise ValueError, "which can be: 'next', 'previous', 'nearest'"
def two_point_interp(times, altitudes):
'''
Draw a straight line between bounding times, find x-intercept.
'''
slope = (altitudes[1] - altitudes[0])/(times[1].jd - times[0].jd)
return Time(times[1].jd - altitudes[1]/slope, format='jd')
def sunrise_dumbinterp(time, location, which='nearest', N=150):
if which == 'nearest' or which == 'next':
# Calculate next sunrise
# Do a coarse search with error O(~7.2 minutes),
# return upper/lower bounds, then interpolate between
# and find x-intercept
next_rise = two_point_interp(*get_rise(time, 1, location,
N=N, return_bounds=True,
return_a=True))
if False:
print(t_horizon, rise_bounds[0].jd,rise_bounds[1].jd)
x_range = np.linspace(rise_bounds[0].jd, rise_bounds[1].jd, 2)
y_range = slope*x_range + rise_alts[1] - slope*rise_bounds[1].jd
plt.plot(x_range, y_range)
plt.plot(x_range, rise_alts,'o')
plt.axvline(t_horizon)
plt.axhline(0)
plt.show()
if which == 'next':
return next_rise
if which == 'nearest' or which == 'previous':
# Calculate previous sunrise
# Do a coarse search with error O(~5 minutes),
# return upper/lower bounds
rise_bounds = get_rise(time, -1, location, N=N,
return_bounds=True)
# Do a fine search with error O(~1 second),
# return the nearest time to sunrise
previous_rise = get_rise(rise_bounds[0],
(rise_bounds[1] - rise_bounds[0]).value,
location, N=N)
if which == 'previous':
return previous_rise
if which == 'nearest':
if abs(previous_rise - time) < abs(next_rise - time):
return previous_rise
else:
return next_rise
raise ValueError, "which can be: 'next', 'previous', 'nearest'"
# Input observer location at Subaru
latitude = '19:49:42.600'
longitude = '-155:28:48.900'
elevation = 0*u.m
location = EarthLocation(lat=latitude, lon=longitude,
height=elevation)
time = Time('2015-05-29',
location=location, scale='utc')
print(sunrise_dumbinterp(time, location, which='next').iso)
In [73]:
%%timeit
sunrise_doublegrid(time, location)
In [74]:
%%timeit
sunrise_dumbinterp(time, location)
In [78]:
print('doublegrid:\t{}\ninterp:\t\t{}'.format(sunrise_doublegrid(time, location),
sunrise_dumbinterp(time, location)))
In [58]:
timeits = []
times = []
N_range = np.logspace(1, 2.77, 10) #np.arange(100, 1100, 100)
print(N_range)
for N in N_range:
times.append(sunrise_dumbinterp(time, location, N=N))
t = %timeit -o sunrise_dumbinterp(time, location, N=N)
timeits.append(t.best)
In [71]:
matplotlib.rcParams['font.size'] = 16
fig, ax = plt.subplots(1, 2, figsize=(16,8))
ax[0].plot_date(Time(times).plot_date, timeits, '-o')
ax[0].set(xlabel='Sunrise [UT]', ylabel='Runtime [s]')
[l.set(rotation=45, ha='right') for l in ax[0].get_xticklabels()]
for x, y, N_i in zip(Time(times).plot_date, timeits, N_range):
ax[0].annotate(int(N_i), xy=(x,y),
ha='right', va='bottom')
ax[1].semilogy(N_range,
np.abs(24*60*60*(Time(times) - times[-1]).value))
ax[1].set(xlabel='N grid points', ylabel='|error| [s]')
plt.show()
In [ ]: