# 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.
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