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)
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))
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)
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]:
In [ ]: