# Untitled

a guest
Nov 30th, 2019
29
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1.     import numpy as np
2.     import matplotlib.pyplot as plt
3.
4.     class Thing(object):
5.         def __init__(self, name):
6.             self.name = name
7.
8.     names = ('Rosetta', 'comet_67P', 'Sun', 'Earth')
9.
10.     things = []
11.     for name in names:
12.         thing = Thing(name)
13.         things.append(thing)
14.         fname = 'Rosetta ' + name + ' BIG horizons_results.txt'
15.
16.         with open(fname, 'r') as infile:
18.
19.         iSOE      = [i for i, line in enumerate(lines) if "\$\$SOE" in line][0]
20.         iEOE      = [i for i, line in enumerate(lines) if "\$\$EOE" in line][0]
21.         lines     = [line.split(',') for line in lines[iSOE+1:iEOE]]
22.         thing.JD  = np.array([float(line[0])    for line in lines])
23.         state     = np.array([[float(x) for x in line[2:8]] for line in lines]).T.copy()
24.         thing.x   = state[0:3]
25.         thing.v   = state[3:6]
26.
27.     rosetta, comet_67P, sun, earth = things
28.
29.     i_2011 = [i for (i, line) in enumerate(lines) if '2011-Jan-01' in line[1]][0]
30.     i_in   = [i for (i, line) in enumerate(lines) if '2011-Jun-08' in line[1]][0]
31.     i_out  = [i for (i, line) in enumerate(lines) if '2014-Jan-20' in line[1]][0]
32.
33.     JD_2011, JD_in, JD_out = [rosetta.JD[i] for i in (i_2011, i_in, i_out)]
34.
35.     years   = 2011 + (rosetta.JD - JD_2011)/365.2564
36.
37.     bodies = (comet_67P, sun, earth)
38.     for body in bodies:
39.         body.distance = np.sqrt(((body.x - rosetta.x)**2).sum(axis=0))
40.         body.speed    = np.sqrt(((body.v - rosetta.v)**2).sum(axis=0))
41.
42.     AU = 1.495978707E+08 # km
43.
44.     print(sun.distance[[i_in, i_out]]/AU)
45.
46.     if True:
47.         plt.figure()
48.         for i, body in enumerate(bodies):
49.             plt.subplot(3, 1, i+1)
50.             plt.plot(years, body.distance/AU)
51.             plt.plot(years[[i_in, i_out]], body.distance[[i_in, i_out]]/AU, 'ok')
52.             plt.title(body.name, fontsize=14)
53.             if body is sun:
54.                 plt.text(2011.8, 4.2, '4.46 AU')
55.                 plt.text(2014.5, 4.2, '4.49 AU', fontsize=10)
56.         plt.suptitle('Rosetta spacecraft distances (AU) and hibernation period \n 2011-06-08 to 2014-01-20', fontsize=13)
57.         plt.show()
RAW Paste Data