Guest User

Untitled

a guest
Jul 14th, 2020
94
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()
RAW Paste Data