Advertisement
Guest User

Untitled

a guest
Jun 19th, 2019
58
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 2.39 KB | None | 0 0
  1. class Lpoint(object):
  2. def __init__(self, name, X0, mus, opts, tol):
  3. self.name = name
  4. self.X0 = X0
  5. self.mus = mus
  6. self.mu0 = self.mus[0]
  7. self.opts = opts
  8. self.tol = tol
  9. self.Xs = []
  10. self.res = []
  11.  
  12. self.solve_mus()
  13. self.X = np.vstack(self.Xs)
  14. self.oks = [res['success'] for res in self.res]
  15. self.ok = all(self.oks)
  16.  
  17. def solve_mus(self):
  18. for mu in self.mus:
  19. results = minimize(self.H, self.X0, method='Nelder-Mead',
  20. options=self.opts, tol=self.tol,
  21. args=(mu))
  22. X = results['x']
  23. self.X0 = X
  24. self.Xs.append(X)
  25. self.res.append(results)
  26.  
  27. def H(self, X, mu):
  28. x, y = X
  29. r1 = np.sqrt((x + mu)**2 + (y)**2)
  30. r2 = np.sqrt((x + mu - 1)**2 + (y)**2)
  31. # zero velocity forces # m1, m2 = 1-mu, mu
  32. xddot = -( (1-mu) * (x+mu) / r1**3 + mu * (x+mu-1) / r2**3 ) + x
  33. yddot = -( (1-mu) * y / r1**3 + mu * y / r2**3 ) + y
  34. return np.sqrt(xddot**2 + yddot**2)
  35.  
  36. import numpy as np
  37. import matplotlib.pyplot as plt
  38. from scipy.optimize import minimize
  39.  
  40. mus = [0.01, 0.02, 0.05, 0.1, 0.2, 0.5]
  41.  
  42. X0s = ((0.7, 0), (1.2, 0), (-1, 0), (0.5, 1), (0.5, -1)) # arbitrary starting point
  43. X0s = [np.array(X0) for X0 in X0s]
  44. names = ('L1', 'L2', 'L3', 'L4', 'L5')
  45.  
  46. opts = {'xtol': 1E-12, 'disp': False, 'maxiter':1000, 'maxfev':2000}
  47.  
  48. Lpts = []
  49. for (X0, name) in zip(X0s, names):
  50. Lpt = Lpoint(name, X0, mus, opts, tol=1E-06)
  51. Lpts.append(Lpt)
  52.  
  53. L1, L2, L3, L4, L5 = Lpts
  54.  
  55. for L in Lpts:
  56. L.solve_mus()
  57.  
  58. all_oks = sum([L.oks for L in Lpts], [])
  59.  
  60. print ('All okay: ', all(all_oks), str(sum(all_oks)) + '/' + str(sum(all_oks)))
  61.  
  62. if True:
  63. plt.figure()
  64. for i, mu in enumerate(mus):
  65. plt.subplot(3, 2, i+1)
  66. x1, x2 = -mu, 1-mu
  67. m1, m2 = 1-mu, mu
  68. for L in Lpts:
  69. x, y = L.X[i]
  70. plt.plot([x], [y], '.k')
  71. plt.plot([x1], [0], 'or', markersize = max(16 * m1, 6))
  72. plt.plot([x2], [0], 'ob', markersize = max(16 * m2, 6))
  73. plt.xlim(-1.5, 1.5)
  74. plt.ylim(-1.0, 1.0)
  75. plt.text(-1.4, 0.8, 'mu='+str(mu), fontsize=16)
  76. hw = 0.1
  77. plt.plot([-hw, hw], [0, 0], '-k', linewidth=0.5)
  78. plt.plot([0, 0], [-hw, hw], '-k', linewidth=0.5)
  79. plt.show()
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement