<jupyter_start><jupyter_text>An Import Loop and a Fiery Reentry<NAME>PyBay, San Francisco2018 August 18<jupyter_code>from IPython.display import HTML

style = """
<style>
#header, div.unselected { display: none !important }
#notebook, .container { padding: 0 !important }
#site { height: 100% !important }
</style>
"""
style = ""

HTML(style)<jupyter_output><empty_output><jupyter_text>1. The design of Skyfield2. The art of the search3. Tiangong-1 before dawn<jupyter_code>from pytz import timezone
eastern = timezone('US/Eastern')
pacific = timezone('US/Pacific')

EDT = -4
PDT = -7<jupyter_output><empty_output><jupyter_text>The design of Skyfield Pure Python No object mutation ```>>> my_list = [...]>>> my_list.sort()>>> print(my_list)```*vs*```>>> print(sorted(my_list))```<jupyter_code>import ephem
boston = ephem.Observer()
boston.lat = '42.37'
boston.lon = '-71.03'
mars = ephem.Mars()
altitudes = []
for hour in range(6):
    boston.date = '2018/8/18 %d' % hour  # awkward
    mars.compute(boston)                 # awkward
    altitudes.append(mars.alt)           # awkward
altitudes
from skyfield.api import load, Topos

ts = load.timescale()
t6 = ts.utc(2018, 8, 19, 6 - PDT)

planets = load('de421.bsp')
earth, sun = planets['earth'], planets['sun']
sf = earth + Topos('37.7749 N', '122.4194 W')
alt, az, distance = sf.at(t6).observe(sun).apparent().altaz()
print(alt)<jupyter_output>-06deg 12' 11.4"<jupyter_text>Separate positions from coordinates*Symptom: “PyEphem’s results don’t match…”* *two axes: which position? which coordinates?* ```astrometric.radec(J2000)astrometric.radec('date')astrometric.ecliptic_latlon()apparent = astrometric.apparent()apparent.radec(J2000)apparent.radec('date')apparent.ecliptic_latlon()``` Explicit units NumPy Arrays Method chaining```apparent = astrometric.apparent()geocentric = boston.at(t)``` ```topos = geocentric.subpoint()``` ```Traceback (most recent call last):  ...  File "/home/brandon/src/skyfield/skyfield/toposlib.py",  line 7, in     from .vectorlib import VectorFunctionImportError: cannot import name 'VectorFunction'``` *flaw in Python's design* http://rhodesmill.org/skyfield/  The art of the search1. Create a number2. Draw a picture3. Launch a search<jupyter_code>def f(hour):
    t = ts.utc(2018, 8, 19, array(hour) - PDT)
    alt, az, distance = sf.at(t).observe(sun).apparent().altaz()
    return alt.degrees

f(6)
plot([6, 7], [f(6), f(7)], 'bo--')
plot([6, 7], [0, 0], 'r-')
h = [6, 6.5, 7]
plot(h, f(h), 'bo--')
plot(h, zeros(len(h)), 'r-')
h = [6, 6.5, 6.75, 7]
plot(h, f(h), 'bo--')
plot(h, zeros(len(h)), 'r-')
h = [6, 6.5, 6.625, 6.75, 7]
plot(h, f(h), 'bo--')
plot(h, zeros(len(h)), 'r-')
h = [6, 6.5, 6.5625, 6.625, 6.75, 7]
plot(h, f(h), 'bo--')
plot(h, zeros(len(h)), 'r-')
interval_start = 7 - 6
interval_end = 1/3600
print(log(interval_start / interval_end) / log(2))
h0 = 6
h1 = 7
while h1 - h0 > 1/3600:
    h = (h1 + h0) / 2
    print(h, h1 - h0)
    if f(h) > 0:
        h1 = h
    else:
        h0 = h
h = [6, 6.5, 6.5625, 6.625, 6.75, 7]
plot(h, f(h), 'bo--')
plot(h, zeros(len(h)), 'r-')
# What if we guess the answer is at the intercept
# of our current two guesses with the f() == 0 line?
# "Secant Method"

h0 = 6
h1 = 7
while abs(h1 - h0) > 1/3600:
    h = (h0 * f(h1) - h1 * f(h0)) / (f(h1) - f(h0))
    h0, h1 = h1, h
    print(h0, h1, h1 - h0)<jupyter_output>7 6.54694463859 -0.453055361411
6.54694463859 6.5528044219 0.00585978330672
6.5528044219 6.55285986302 5.54411210123e-05<jupyter_text>*finesse*<jupyter_code>h0 = 13
h1 = 14
while abs(h1 - h0) > 1/3600:
    h = (h0 * f(h1) - h1 * f(h0)) / (f(h1) - f(h0))
    h0, h1 = h1, h
    print(h0, h1, h1 - h0)
def draw_arrow(x0, y0, x1, y1):
    annotate(
        '', (x1, y1), (x0, y0),
        xycoords='data',
        arrowprops={
            'arrowstyle':'->, head_width=0.5, head_length=1',
                    'linewidth': 1,
                    'color':'k'},
    )
fig, ax = plt.subplots()
h0 = 13
h1 = 14
axhline(0, color='r')
while abs(h1 - h0) > 1/3600:
    h = (h0 * f(h1) - h1 * f(h0)) / (f(h1) - f(h0))
    plot(h0, f(h0), 'bo')
    draw_arrow(h1, f(h1), h, 0)
    h0, h1 = h1, h
# The secant method can only "see"
# tiny fragments of the whole curve.

t = arange(*ax.get_xlim(), 0.1)
ax.plot(t, f(t))
fig<jupyter_output><empty_output><jupyter_text>secant method isdifficult to control*but*efficient does safety have to be expensive?<jupyter_code>%timeit f(6)
%timeit f([6])
%timeit f([6, 7, 8, 9])
h0 = 6
h1 = 7
h = linspace(h0, h1, 5)
h
f(h)
f(h) > 0
argmax(f(h) > 0)
h0, h1 = 6, 7
while h1 - h0 > 1/3600:
    h = linspace(h0, h1, 5)
    print(h)
    i = argmax(f(h) > 0)
    h0, h1 = h[i-1], h[i]
def secant_search(h0, h1, epsilon, f):
    while abs(h1 - h0) > epsilon:
        f0, f1 = f(h0), f(h1)
        h = (h0 * f1 - h1 * f0) / (f1 - f0)
        h0, h1 = h1, h
    return h
def phalanx_search(h0, h1, epsilon, length, f):
    while h1 - h0 > epsilon:
        h = linspace(h0, h1, length)
        i = argmax(f(h) > 0)
        h0, h1 = h[i-1], h[i]
    return h0
%timeit secant_search(6, 7, 1/3600, f)
%timeit phalanx_search(6, 7, 1/3600, 4, f)
%timeit phalanx_search(6, 7, 1/3600, 8, f)
%timeit phalanx_search(6, 7, 1/3600, 16, f)
def secant_search(h0, h1, epsilon, f):
    while abs(h1 - h0) > epsilon:
        f0, f1 = f(h0), f(h1)
        h = (h0 * f1 - h1 * f0) / (f1 - f0)
        h0, h1 = h1, h
    return h
def secant_search2(h0, h1, epsilon, f):
    while abs(h1 - h0) > epsilon:
        f0, f1 = f([h0, h1])
        h = (h0 * f1 - h1 * f0) / (f1 - f0)
        h0, h1 = h1, h
    return h
from scipy.optimize import brentq
print(brentq(f, 6, 7, xtol=1/3600))
print(secant_search(6, 7, 1/3600, f))
print(secant_search2(6, 7, 1/3600, f))
print(phalanx_search(6, 7, 1/3600, 3, f))
print(phalanx_search(6, 7, 1/3600, 5, f))
print(phalanx_search(6, 7, 1/3600, 9, f))
%timeit brentq(f, 0, 12, xtol=1/3600)
%timeit secant_search(0, 12, 1/3600, f)
%timeit secant_search2(0, 12, 1/3600, f)
%timeit phalanx_search(0, 12, 1/3600, 4, f)
%timeit phalanx_search(0, 12, 1/3600, 8, f)
%timeit phalanx_search(0, 12, 1/3600, 16, f)
from time import time
for i in range(3, 30):
    t0 = time()
    phalanx_search(6, 7, 1/3600, i, f)
    dt = time() - t0
    plot(i, dt, 'o')
ylim(0);<jupyter_output><empty_output><jupyter_text>Have we found sunrise?<jupyter_code>h0 // 1, h0 * 60 // 1 % 60, h0 * 3600 // 1 % 60<jupyter_output><empty_output><jupyter_text>Sunday, August 19, 2018            Pacific Daylight Time                    Sun                      Begin civil twilight6:01 a.m.            Sunrise6:29 a.m.            Sun transit1:13 p.m.            Sunset7:57 p.m.            End civil twilight8:25 p.m.Moon               Moonset12:50 a.m.               Moonrise3:15 p.m.               Moon transit8:24 p.m.         Moonset1:30 a.m. on following day “For computational purposes, sunrise or sunset is defined to occur when the geometric zenith distance of center of the Sun is 90.8333 degrees.”http://aa.usno.navy.mil/faq/docs/RST_defs.php<jupyter_code>def phalanx_search(t0, t1, epsilon, f):
    jd0, jd1 = t0.tt, t1.tt
    while jd1 - jd0 > epsilon:       # <NAME>
        jd = linspace(jd0, jd1, 16)
        t = ts.tt(jd=jd)
        i = argmax(f(t) > 0)         # search still looks for zero
        jd0, jd1 = jd[i-1], jd[i]
    return ts.tt(jd=jd0)

def sunrise(topos):
    def sunrise_altitude(t):
        alt, az, distance = topos.at(t).observe(sun).apparent().altaz()
        return alt.degrees + 0.8333  # it's f() that I've adjusted
    return sunrise_altitude
t = phalanx_search(
    ts.utc(2018, 8, 19, 6 - PDT),
    ts.utc(2018, 8, 19, 7 - PDT),
    1e-8,
    sunrise(sf),
)

t.astimezone(pacific).isoformat(' ')<jupyter_output><empty_output><jupyter_text>Sun                      Sunrise6:29 a.m.<jupyter_code># What if we want the sun's maximum altitude?

for h in 10, 11, 12, 13, 14:
    t = ts.utc(2018, 8, 19, h - PDT)
    alt, az, distance = sf.at(t).observe(sun).apparent().altaz()
    print(h, alt)
t = ts.utc(2018, 8, 19, [h - PDT for h in range(10, 15)])
alt, az, distance = sf.at(t).observe(sun).apparent().altaz()
d = alt.degrees
print(d)
print(d)
i = argmax(d)
print(i)
print(d[i-1], d[i], d[i+1])
def sunalt(t):
    alt, az, distance = sf.at(t).observe(sun).apparent().altaz()
    return alt.degrees
def maximize(t0, t1, epsilon, f):
    jd0, jd1 = t0.tt, t1.tt
    while jd1 - jd0 > epsilon:
        jd = linspace(jd0, jd1, 16)
        t = ts.tt(jd=jd)
        i = argmax(f(t))
        i = clip(i, 1, 14)
        jd0, jd1 = jd[i-1], jd[i+1]
    return ts.tt(jd=jd0)
t = maximize(ts.utc(2018, 8, 19, 10 - PDT),
             ts.utc(2018, 8, 19, 14 - PDT), 1e-8, sunalt)
print(t.astimezone(pacific))

day = ts.utc(2018, 8, 19, arange(-PDT, 24-PDT, 0.1))
plot(day.utc_datetime(), sunalt(day))
plot(t.utc_datetime(), sunalt(t), 'ro')<jupyter_output>2018-08-19 13:13:07.075000-07:00<jupyter_text>Tiangong-1 before dawn<jupyter_code>satellites = load.tle('stations-2018-03-31.tle')
sats = set(sat for sat in satellites.values())
tg = satellites['TIANGONG 1']
tg
# When was sunrise?

bluffton = Topos('40.8953 N', '83.8888 W')
EST = -4
t = phalanx_search(
    ts.utc(2018, 8, 19, 4 - EDT),
    ts.utc(2018, 8, 19, 10 - EDT),
    1e-8,
    sunrise(earth + bluffton),
)
t.astimezone(eastern).isoformat(' ')
# How high will Tiangong-1 pass overhead?
# I knew better than to do this:

t0, t1 = ts.utc(2018, 4, 1, [3-EDT, 7-EDT])
f = lambda t: (tg - bluffton).at(t).altaz()[0].degrees
tm = maximize(t0, t1, 1/24/3600, f)
print(tm.astimezone(eastern))
print(f(tm))
t = ts.utc(2018, 4, 1, arange(0-PDT, 12-PDT, 1/60))
alt, az, distance = (tg - bluffton).at(t).altaz()
plot(t.astimezone(eastern), alt.degrees)
plot(tm.astimezone(eastern), f(tm), 'ro')
print(alt)
above_horizon = alt.degrees > 0
print(above_horizon[300:320])
boundaries, = diff(above_horizon).nonzero()
print(boundaries)

passes = boundaries.reshape(len(boundaries) // 2, 2)
print(passes)
for i, j in passes:
    tm = maximize(t[i], t[j], 1/24/3600, f)
    dt = tm.astimezone(eastern).replace(microsecond=0)
    print(i, j, '-', dt, ' max:', f(tm))
def plot_sky(pass_indices):
    i, j = pass_indices
    
    # Set up the polar plot.
    ax = plt.subplot(111, projection='polar')
    ax.set_rlim([0, 90])
    ax.set_theta_zero_location('N')
    ax.set_theta_direction(-1)
    
    # Draw line and labels.
    θ = az.radians
    r = 90 - alt.degrees
    ax.plot(θ[i:j], r[i:j], 'ro--')
    for k in range(i, j):
        if r[k] < 90:
            text = t[k].astimezone(eastern).strftime('%H:%M')
            ax.text(θ[k], r[k], text, ha='right', va='bottom')
plot_sky(passes[0])
plot_sky(passes[1])
ranking = []
for sat in sats:
    position = sat.at(t)
    v = position.velocity.km_per_s
    speed = sqrt((v * v).sum(axis=0))
    if any(isnan(speed)):
        continue
    ranking.append((speed.max(), sat.name))
sorted(ranking, reverse=True)[:12]
from skyfield.constants import ERAD
ranking = []
for sat in sats:
    position = sat.at(t)
    p = position.position.km
    elevation = sqrt((p * p).sum(axis=0)) - ERAD / 1e3
    if isnan(elevation.min()):
        continue
    ranking.append((elevation.min(), sat.name))
sorted(ranking)[:12]
speeds = []
for sat in sats:
    position = sat.at(t)
    p = position.position.km
    v = position.velocity.km_per_s
    speed = sqrt((v * v).sum(axis=0))
    elevation = sqrt((p * p).sum(axis=0)) - ERAD / 1e3
    if any(isnan(speed)):
        continue
    degrees_per_second = (speed / elevation) / tau * 360
    speeds.append((degrees_per_second.max(), sat.name))
sorted(speeds, reverse=True)[:12]
# <NAME> - @brandon_rhodes - PyBay 2018

plot_sky(passes[1])
from IPython.display import HTML
HTML("""
<style>
h1, h2, h3, .rendered_html p {
  text-align: center;
}
.anchor-link {
  display: none;  /* ruins centering of titles */
}
li, pre {
  text-align: left;
}
.text_cell_render.rendered_html {
  font-size: larger;
  text-align: center;
}
</style>
<script src="slide_timer.js"></script>
<script>
setup_timer();
</script>
""")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17745683930942896, 'fraction_numerical': 0.07818725099601594, 'mean_word_length': 3.1692041522491348, 'pattern_counts': {'":': 0, '<': 50, '<?xml version=': 0, '>': 80, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '456415', 'license': 'MIT', 'n_tokens_mistral': 5735, 'n_tokens_neox': 5156, 'n_words': 1242, 'path': 'Talks/2018-08-pybay.ipynb', 'repo_name': 'Grant-Steinfeld/astronomy-notebooks', 'size': '12070'}