Guest User

Untitled

a guest
Nov 20th, 2018
114
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. class Thing(object):
  2.     def __init__(self, name):
  3.         self.name = name
  4.  
  5. import numpy as np
  6. import matplotlib.pyplot as plt
  7. from mpl_toolkits.mplot3d import Axes3D
  8.  
  9. names  = ('Earth', 'Huygens', 'Saturn', 'Sun', 'Titan')
  10. fnames = ['Huygens to Titan ' + n + ' horizons_results.txt' for n in names]
  11.  
  12. halfpi, pi, twopi = [f*np.pi for f in [0.5, 1.0, 2.0]]
  13. degs, rads        = 180./pi, pi/180.
  14. AU                = 149597870.700  # kilometers
  15.  
  16. JDs, posns, vels, linez = [], [], [], []
  17. for fname in fnames:
  18.     with open(fname, 'r') as infile:
  19.  
  20.         lines = infile.read().splitlines()
  21.  
  22.     iSOE = [i for i, line in enumerate(lines) if "$$SOE" in line][0]
  23.     iEOE = [i for i, line in enumerate(lines) if "$$EOE" in line][0]
  24.  
  25.     print iSOE, iEOE, lines[iSOE], lines[iEOE]
  26.  
  27.     lines = [line.split(',') for line in lines[iSOE+1:iEOE]]
  28.     JD  = np.array([float(line[0]) for line in lines])
  29.     pos = np.array([[float(item) for item in line[2:5]] for line in lines])
  30.     vel = np.array([[float(item) for item in line[5:8]] for line in lines])
  31.  
  32.     pos, vel = [thing.T for thing in pos, vel]
  33.  
  34.     JDs.append(JD)
  35.     posns.append(pos)
  36.     vels.append(vel)
  37.     linez.append(lines)
  38.  
  39. things = []
  40. for JD, pos, vel, name in zip(JDs, posns, vels, names):
  41.     thing = Thing(name)
  42.     thing.JD  = JD
  43.     thing.pos = pos
  44.     thing.vel = vel
  45.     things.append(thing)
  46.  
  47. Earth, Huygens, Saturn, Sun, Titan = things
  48.  
  49. minutes = np.arange(0-240, 0+30)
  50.  
  51. d = Huygens.pos - Titan.pos
  52. r = np.sqrt((d**2).sum(axis=0))
  53. x, y, z = d
  54. rxy     = np.sqrt(x**2 + y**2)
  55. ph = np.arctan2(y, x)
  56. th = np.arctan2(z, rxy)
  57.  
  58. v = Huygens.vel - Titan.vel
  59. s = np.sqrt((v**2).sum(axis=0))
  60. vn  = v/s
  61. rn  = d/r
  62. vdotr = (vn*-rn).sum(axis=0)
  63.    
  64. if True:
  65.     rr = r[3430-240:3430+30]
  66.     plt.figure()
  67.     plt.plot(minutes, ph[3430-240:3430+30])
  68.     plt.plot(minutes, th[3430-240:3430+30])
  69.     plt.title('Huygens approach theta and phi', fontsize=18)
  70.     plt.show()
  71.  
  72. if True:
  73.     plt.figure()
  74.     plt.subplot(2, 1, 1)
  75.     plt.title('Huygens vs Titan v dot r vs time (min)', fontsize=20)
  76.     plt.plot(minutes, vdotr[3430-240:3430+30])
  77.     plt.subplot(2, 1, 2)
  78.     plt.title('Huygens vs Titan nadir angle', fontsize=20)
  79.     plt.plot(minutes, degs*np.arccos(vdotr[3430-240:3430+30]))
  80.     plt.show()
  81.  
  82. if True:
  83.     rr = r[3430-240:3430+30]
  84.     plt.figure()
  85.     plt.subplot(3, 1, 1)
  86.     plt.plot(minutes[1:], rr[1:]-rr[:-1])
  87.     plt.ylim(-360, 10)
  88.     plt.title('Huygens vs Titan dr/dt (km/min)', fontsize=18)
  89.     plt.subplot(3, 1, 2)
  90.     plt.plot(minutes[1:], -(rr[1:]-rr[:-1]))
  91.     plt.yscale('log')
  92.     plt.title('Huygens vs Titan dr/dt (km/min)', fontsize=18)
  93.     plt.subplot(3, 1, 3)
  94.     plt.plot(minutes, degs*np.arccos(vdotr[3430-240:3430+30]))
  95.     plt.title('Huygens vs Titan nadir angle', fontsize=18)
  96.     plt.show()
  97.  
  98. if True:
  99.     rr = r[3430-240:3430+30]
  100.     plt.figure()
  101.     plt.subplot(2, 1, 1)
  102.     plt.plot(minutes, rr-rr.min())
  103.     plt.subplot(2, 1, 2)
  104.     plt.plot(minutes, rr-rr.min())
  105.     plt.ylim(-5, 200)    
  106.     plt.suptitle('Huygens altitude (km) vs time (min)', fontsize=18)
  107.     plt.show()
RAW Paste Data