In [158]:

    

import ephem
import matplotlib.pyplot as plt
from datetime import datetime, timedelta
from ephem import newton
%matplotlib inline


    

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 [141]:

    

utc_offset = timedelta(hours=6)
launch_time = datetime(2016, 9, 15)
launch_time_utc = datetime(2016, 9, 15) + utc_offset
print(launch_time_utc)


    

    




2016-09-15 06:00:00

In [142]:

    

fort_sumner.date = '2016/09/15 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-15 06:44:28.174546
2016-09-14 19:00:11.678369
0.000073

In [145]:

    

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 [169]:

    

fig = plt.figure(figsize=(15, 10))
times_local = [t - utc_offset for t in times]
ax = plt.plot(times_local, np.rad2deg(sun_elevation), '-', label='Sun')
ax = plt.plot(times_local, np.rad2deg(crab_elevation), '-', label='Crab')
plt.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)
plt.axvspan(start_obs, end_obs, facecolor='gray', alpha=0.2, label='At float')
plt.axvline(launch_time, label='Launch')
fig.autofmt_xdate()
plt.xlabel(times[0])
plt.ylim(0)
plt.title('Fort Sumner')
plt.legend()
plt.show()


    

In [166]:

    

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


    

In [167]:

    

x = fort_sumner.date
ephem.newton(f, x, x + 0.01)


    

    
Out[167]:





42629.43602009599

In [168]:

    

print 'At the time and date', fort_sumner.date
print 'the solar altitude is', sun.alt


    

    




At the time and date 2016/9/17 22:27:52
the solar altitude is 30:00:00.0

In [ ]: