In [235]:
import ephem
import matplotlib.pyplot as plt
from datetime import datetime, timedelta
from ephem import newton
import matplotlib.dates as mdates
%matplotlib inline

In [242]:
from matplotlib.dates import HourLocator

In [122]:
sun = ephem.Sun()
crab = ephem.readdb('Crab Nebula,f|R,5:34:31.9,22:0:52,8.4,2000,360')
fort_sumner = ephem.Observer()
fort_sumner.lat = '34.4731'
fort_sumner.lon = '-104.2422'
fort_sumner.elevation = 4000
fort_sumner.pressure = 0

In [252]:
utc_offset = timedelta(hours=6)
launch_time = datetime(2016, 9, 12)
launch_time_utc = launch_time + utc_offset
print(launch_time_utc)


2016-09-12 06:00:00

from http://www.wolframalpha.com/input/?i=2016%2F9%2F15+sun+set+in+fort+sumner%2C+nm Sun set time is 7:04 pm Sun rise time is 6:42 am


In [253]:
fort_sumner.date = '2016/09/12 18:00'
sun.compute(fort_sumner)
print(sun.rise_time.datetime() - utc_offset)
print(sun.set_time.datetime() - utc_offset)
fort_sumner.date = sun.rise_time
sun.compute(fort_sumner)
print('%f' % np.rad2deg(sun.alt))


2016-09-12 09:09:42.714069
2016-09-12 16:35:51.711481
29.999629

In [254]:
times = [launch_time_utc + timedelta(0, 3*d) for d in range(1,10*6000)]
crab_elevation = []
sun_elevation = []
for d in times:
    fort_sumner.date = d
    crab.compute(fort_sumner)
    sun.compute(fort_sumner)
    sun_elevation.append(sun.alt)
    crab.compute(fort_sumner)
    crab_elevation.append(crab.alt)

In [255]:
times_local = [t - utc_offset for t in times]

fig = plt.figure(figsize=(15, 10))
ax1 = fig.add_subplot(111)
plt.axvspan(start_obs, end_obs, facecolor='gray', alpha=0.2, label='At float')
plt.axvline(launch_time, label='Launch')
ax2 = ax1.twiny()
ax2.set_xlim(launch_time + utc_offset, launch_time+timedelta(2) + utc_offset)
ax2.set_xlabel('UTC')
ax2.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H:%S'))

ax1.plot(times_local, np.rad2deg(sun_elevation), '-', label='Sun')
ax1.plot(times_local, np.rad2deg(crab_elevation), '-', label='Crab')
ax1.set_xlim(launch_time, launch_time+timedelta(2))
plt.grid(True)
plt.axhspan(30, 70, facecolor='gray', alpha=0.2, label='Observable')

start_obs = launch_time + timedelta(hours=6) + timedelta(hours=3)
end_obs = launch_time + timedelta(1) + timedelta(hours=15)
ax1.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H:%S'))
fig.autofmt_xdate()
ax1.set_xlabel('Local Time (MT)')
plt.ylim(0)
plt.legend()
plt.show()



In [256]:
def f(x):
    fort_sumner.date = x
    sun.compute(fort_sumner)
    return sun.alt - ephem.degrees('30')

In [264]:
fort_sumner.horizon = '30:00'
guess_list = [launch_date - timedelta(2), launch_date - timedelta(1), launch_date + timedelta(1)]
for guess in guess_list:
    fort_sumner.date = guess
    print 'Sun rising', fort_sumner.previous_rising(ephem.Sun()), 'Sun setting', fort_sumner.previous_setting(ephem.Sun())
    print 'Crab rising', fort_sumner.previous_rising(crab), 'Crab setting', fort_sumner.previous_setting(crab)


Sun rising 2016/9/12 15:08:22 Sun setting 2016/9/12 22:37:12
Crab rising 2016/9/12 08:33:36 Crab setting 2016/9/12 17:35:40
Sun rising 2016/9/13 15:09:13 Sun setting 2016/9/13 22:35:38
Crab rising 2016/9/13 08:29:40 Crab setting 2016/9/13 17:31:44
Sun rising 2016/9/15 15:10:57 Sun setting 2016/9/15 22:32:28
Crab rising 2016/9/15 08:21:49 Crab setting 2016/9/15 17:23:52

In [266]:
solar_obs_time = datetime(2016,9,12,16,37)-datetime(2016,9,12,10,0)+datetime(2016,9,12,16,35)-datetime(2016,9,12,9,9)

In [271]:
solar_obs_time.total_seconds()/(60.0*60.0)


Out[271]:
14.05

In [ ]: