﻿

# Untitled

a guest
May 20th, 2019
346
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. def derivxx(x, mu):
2.     r1 = np.sqrt((x      + mu)**2)
3.     r2 = np.sqrt((x - 1. + mu)**2)
4.     term_1 = x
5.     term_2 = -(1.-mu) * (x + mu) / r1**3
6.     term_3 =     -mu  * (x - 1. + mu) / r2**3
7.
8.     xddot  = term_1 + term_2 + term_3
9.
10.     return xddot
11.
12. def get_x_L1():
13.
14.     mu_EM  = GMm / (GMm + GMe)
15.     mu_SE  = GMe / (GMe + GMs)
16.     mu_SEb = (GMe + GMm) / (GMe + GMm + GMs)
17.
18.     mu = mu_SEb
19.
20.     xtol = 1E-10
21.     answer, blob = spo.brentq(derivxx, 0.9, 0.9999, args=(mu),
22.                               full_output = True, xtol=xtol)
23.     print "converged = ", blob.converged
24.     if blob.converged:
26.     else:
27.         x_L1 = None
28.
29.     return x_L1
30.
31. def ROTIT(X, th):
32.
33.     sth, cth, zth, oth = ( np.sin(th), np.cos(th),
34.                            np.zeros_like(th), np.ones_like(th) )
35.
36.     xR = (X * np.vstack((cth, -sth,  zth))).sum(axis=0)
37.     yR = (X * np.vstack((sth,  cth,  zth))).sum(axis=0)
38.     zR = (X * np.vstack((zth,  zth,  oth))).sum(axis=0)
39.
40.     return np.vstack([xR, yR, zR])
41.
42. def cosine_ramp(a, b, n):
43.     x    = np.linspace(0, np.pi, n)
44.     ramp = 0.5*(1. - np.cos(x))
45.     y    = a + (b-a)*ramp
46.
47.     return y
48.
49. import numpy as np
50. import matplotlib.pyplot as plt
51. from mpl_toolkits.mplot3d import Axes3D
52.
53. halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]
54. GMs, GMe, GMm     = (1.32712440018E+20, 3.98600418E+14, 4.9048695E+12)  # m^3 s^-2
55.
56. import json
57. import scipy.optimize as spo
58. from scipy.optimize import minimize as mini
59.
63.
64. x_L1 = get_x_L1()  # about 0.99  # I'm going to say the sun doesn't move here
65.
66. EMB_from_SUN  = EMB_pos - SUN_pos
67.
68. L1_from_SUN   = x_L1 * EMB_from_SUN
69.
70. SOHO_from_SUN = SOHOH_pos - SUN_pos
71.
72. SOHO_from_EMB = SOHO_from_SUN - EMB_from_SUN
73. SOHO_from_L1  = SOHO_from_SUN - L1_from_SUN
74.
75. EMB_from_L1   = EMB_from_SUN - L1_from_SUN
76.
77. theta = np.arctan2(EMB_from_SUN[1], EMB_from_SUN[0])  # OK here is the angle!
78.
79. SOHO_from_L1_rot  = ROTIT(SOHO_from_L1, -theta)
80. EMB_from_L1_rot   = ROTIT(EMB_from_L1,  -theta)
81.
82. fig = plt.figure(figsize=[10, 8])
83.
84. ax = fig.add_subplot(1, 1, 1, projection='3d')
85. x, y, z = SOHO_from_L1_rot
86. ax.plot(x, y, z, '-c', linewidth=0.7)
87.
88. nn1, nn2 = 20, 20  # first and last twenty point
89. ax.plot(x[:nn1], y[:nn1], z[:nn1], '.k')
90. ax.plot(x[-nn2:], y[-nn2:], z[-nn2:], '.k')
91. ax.plot([0], [0], [0], 'or')
92. x, y, z = EMB_from_L1_rot
93. ax.plot(x, y, z, marker='o', fillstyle='full', markeredgecolor='blue',
94.         color='blue', markeredgewidth=0.0)
95.
96. if True:
97.     # Build a set of views
98.     ns, nm = 5, 31
99.     azis   = ns*[0] + cosine_ramp(0, -90, nm).tolist() + ns*[-90] + nm*[-90] + ns*[-90]
100.     altis  = ns*[0] + nm*[0] + ns*[0] + cosine_ramp(0, 90, nm).tolist() + ns*[90]
101.     views  = zip(altis, azis)
102.
103.     for i, (alt, az) in enumerate(views[::10]):
104.         ax.view_init(alt, az)
105.         plt.savefig('SOHO_3D_' + str(10000+i)[1:])
106.
107. plt.show()
RAW Paste Data