Since the last post I've been playing more with PyEphem. One of the many interesting features of this astronomical computations module is the ability to read external orbital elements.

Here is an example on how to create an ISON comet object.

We can get the ISON comet orbital data from this link.

In [2]:

import ephem
import numpy as np
import matplotlib.pyplot as plt

ISON = "C/2012 S1 (ISON),h,11/28.7744/2013,62.3994,295.6528,345.5644,1.000002,0.012444,2000,7.5,3.2"
comet = ephem.readdb(ISON)
print(comet)

<ephem.HyperbolicBody 'C/2012 S1 (ISON)' at 0x7f90bdd4f710>

Now let's create a date vector from today up to the end of year and check ISON distance from the Earth and the Sun.

In [3]:

from pandas import DataFrame, date_range

days = date_range(start='2013/10/28', end='2013/12/31', freq='D')

earth_d, sun_d = [], [] 
for day in days:
    comet.compute(day)
    earth_d.append(comet.earth_distance)
    sun_d.append(comet.sun_distance)
    
df = DataFrame(np.c_[earth_d, sun_d], index=days, columns=['Earth', 'Sun'])

from scipy import constants
conversion = constants.au * 1e-3
df['Sun km'] = df['Sun'] * conversion
df['Earth km'] = df['Earth'] * conversion

df.head(5)

Out[3]:

	Earth	Sun	Sun km	Earth km
2013-10-28	1.343978	1.091559	1.632950e+08	2.010563e+08
2013-10-29	1.315059	1.068280	1.598124e+08	1.967301e+08
2013-10-30	1.286399	1.044747	1.562919e+08	1.924426e+08
2013-10-31	1.258032	1.020949	1.527318e+08	1.881990e+08
2013-11-01	1.229998	0.996874	1.491302e+08	1.840051e+08

In [4]:

kw = dict(markeredgecolor='k', marker='o', linestyle='none')
arrowprops = dict(arrowstyle="->", connectionstyle="arc3,rad=.2", color='k')

fig, ax = plt.subplots(figsize=(8, 4))
sx, sy, sm = df['Sun'].argmin(), df['Sun'].min(), df['Sun km'].min()
ex, ey, em = df['Earth'].argmin(), df['Earth'].min(), df['Earth km'].min()

ax.plot(ex, ey, markerfacecolor='#339933', **kw)
ax.plot(sx, sy, markerfacecolor='#ff3333', **kw)

ax.annotate('Closest to Earth is\n%s km' % em,
            xy=(ex, ey), xycoords='data',
            xytext=(-20, 30), textcoords='offset points', arrowprops=arrowprops)

ax.annotate('Closest to the Sun is\n%s km' % sm,
            xy=(sx, sy), xycoords='data',
            xytext=(-180, 30), textcoords='offset points', arrowprops=arrowprops)

ax = df[['Earth', 'Sun']].plot(ax=ax, zorder=0, color=['#339933', '#ff3333'], linewidth=2)
_ = ax.set_ylabel('Distance [AU]')

Finally, if ISON survives its "sungrazing orbit," we can try to see it when it will be closest to the Earth.

In [5]:

obs = ephem.Observer()
obs.lon = np.deg2rad(-46.665)
obs.lat = np.deg2rad(-23.473)
obs.date = ex
comet.compute(obs)
print("Date: %s\nRight ascension: %s\nDeclination: %s\nMagnitude: %s" %
      (obs.date, comet.ra, comet.dec, comet.mag))

Date: 2013/12/27 00:00:00
Right ascension: 16:22:39.01
Declination: 52:59:54.6
Magnitude: 5.69

In [6]:

from IPython.display import YouTubeVideo
YouTubeVideo("A1yH_DuC88M", autoplay=0, theme="light", color="red")

Out[6]:

In [7]:

HTML(html)

Out[7]:

This post was written as an IPython notebook. It is available for download or as a static html.

python4oceanographers by Filipe Fernandes is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.
Based on a work at https://ocefpaf.github.io/.

python4oceanographers

Turning ripples into waves

Catching ISON fever!

Comments