Guest User

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:
  17.             lines = infile.readlines()
  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