SHARE
TWEET

Untitled

a guest May 20th, 2019 226 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:
  25.         x_L1 = answer
  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.  
  60. EMB_pos = np.load('EMB_pos.npy')
  61. SOHOH_pos = np.load('SOHOH_pos.npy')
  62. SUN_pos = np.load('SUN_pos.npy')
  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
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
Top