Untitled

a guest
Jul 14th, 2020
266
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1.     import numpy as np
2.     import matplotlib.pyplot as plt
3.
4.     halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]
5.     to_degs, to_rads = 180/pi, pi/180
6.
7.     e = 0.055
8.     a = 1.0
9.     mu = 1.0
10.     T = twopi # period for a=1 and mu=1
11.
12.     # from https://space.stackexchange.com/a/43237/12102
13.
14.     theta = np.linspace(-pi, pi, 501)
15.     E = 2. * np.arctan(np.sqrt((1.-e)/(1.+e)) * np.tan(theta/2))
16.     t = a * np.sqrt(a/mu) * (E - e * np.sin(E))
17.
18.     libration = theta - twopi * t/T
19.
20.     if True:
21.         plt.figure()
22.         fs = 14
23.         plt.subplot(2, 2, 1)
24.         plt.plot(t/T, theta*to_degs)
25.         plt.xlabel('time (fraction of period)', fontsize=fs)
26.         plt.ylabel('theta (deg)', fontsize=fs)
27.         plt.xlim(-0.5, 0.5)
28.         plt.ylim(-180, 180)
29.         plt.subplot(2, 2, 2)
30.         plt.plot(theta*to_degs, t/T)
31.         plt.xlabel('theta (deg)', fontsize=fs)
32.         plt.ylabel('time (fraction of period)', fontsize=fs)
33.         plt.xlim(-180, 180)
34.         plt.ylim(-0.5, 0.5)
35.         plt.subplot(2, 1, 2)
36.         plt.plot(t/T, libration*to_degs)
37.         plt.xlabel('time (fraction of period)', fontsize=fs)
38.         plt.ylabel('apparent libration (deg)', fontsize=fs)
39.         plt.xlim(-0.5, 0.5)
40.         plt.text(0.15, 0, 'max: ' + str(round((libration*to_degs).max(), 2))  + u'\xb0', fontsize=fs)
41.         # https://stackoverflow.com/questions/3215168/how-to-get-character-in-a-string-in-python
42.         plt.show()