Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- class Lpoint(object):
- def __init__(self, name, X0, mus, opts, tol):
- self.name = name
- self.X0 = X0
- self.mus = mus
- self.mu0 = self.mus[0]
- self.opts = opts
- self.tol = tol
- self.Xs = []
- self.res = []
- self.solve_mus()
- self.X = np.vstack(self.Xs)
- self.oks = [res['success'] for res in self.res]
- self.ok = all(self.oks)
- def solve_mus(self):
- for mu in self.mus:
- results = minimize(self.H, self.X0, method='Nelder-Mead',
- options=self.opts, tol=self.tol,
- args=(mu))
- X = results['x']
- self.X0 = X
- self.Xs.append(X)
- self.res.append(results)
- def H(self, X, mu):
- x, y = X
- r1 = np.sqrt((x + mu)**2 + (y)**2)
- r2 = np.sqrt((x + mu - 1)**2 + (y)**2)
- # zero velocity forces # m1, m2 = 1-mu, mu
- xddot = -( (1-mu) * (x+mu) / r1**3 + mu * (x+mu-1) / r2**3 ) + x
- yddot = -( (1-mu) * y / r1**3 + mu * y / r2**3 ) + y
- return np.sqrt(xddot**2 + yddot**2)
- import numpy as np
- import matplotlib.pyplot as plt
- from scipy.optimize import minimize
- mus = [0.01, 0.02, 0.05, 0.1, 0.2, 0.5]
- X0s = ((0.7, 0), (1.2, 0), (-1, 0), (0.5, 1), (0.5, -1)) # arbitrary starting point
- X0s = [np.array(X0) for X0 in X0s]
- names = ('L1', 'L2', 'L3', 'L4', 'L5')
- opts = {'xtol': 1E-12, 'disp': False, 'maxiter':1000, 'maxfev':2000}
- Lpts = []
- for (X0, name) in zip(X0s, names):
- Lpt = Lpoint(name, X0, mus, opts, tol=1E-06)
- Lpts.append(Lpt)
- L1, L2, L3, L4, L5 = Lpts
- for L in Lpts:
- L.solve_mus()
- all_oks = sum([L.oks for L in Lpts], [])
- print ('All okay: ', all(all_oks), str(sum(all_oks)) + '/' + str(sum(all_oks)))
- if True:
- plt.figure()
- for i, mu in enumerate(mus):
- plt.subplot(3, 2, i+1)
- x1, x2 = -mu, 1-mu
- m1, m2 = 1-mu, mu
- for L in Lpts:
- x, y = L.X[i]
- plt.plot([x], [y], '.k')
- plt.plot([x1], [0], 'or', markersize = max(16 * m1, 6))
- plt.plot([x2], [0], 'ob', markersize = max(16 * m2, 6))
- plt.xlim(-1.5, 1.5)
- plt.ylim(-1.0, 1.0)
- plt.text(-1.4, 0.8, 'mu='+str(mu), fontsize=16)
- hw = 0.1
- plt.plot([-hw, hw], [0, 0], '-k', linewidth=0.5)
- plt.plot([0, 0], [-hw, hw], '-k', linewidth=0.5)
- plt.show()
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement