Guest User

Untitled

a guest
Oct 31st, 2020
174
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.     a_list = [0.723, 1.0, 1.524, 5.204]
  5.     T_list = [a**1.5 for a in a_list]
  6.  
  7.     times = 2 * np.pi * np.linspace(0, 10, 2001)
  8.  
  9.     phases = [times/T for T in T_list]
  10.  
  11.     positions = np.array([[a * f(phase) for f in (np.cos, np.sin, np.zeros_like)]
  12.                           for (a, phase) in zip(a_list, phases)])
  13.  
  14.     venus, earth, mars, jupiter = positions
  15.  
  16.  
  17.     T_moon = 0.5
  18.     a_moon = 1.3 / 150  # 1.3 million km divided by 1 AU
  19.  
  20.     phase_moon = times / T_moon
  21.  
  22.     moon_p = np.array([a_moon * f(phase_moon) for
  23.                        f in (np.cos, np.sin, np.zeros_like)])
  24.     moon_r = np.array([a_moon * f(-phase_moon) for
  25.                        f in (np.cos, np.sin, np.zeros_like)])
  26.  
  27.     incs = [f * np.pi / 180 for f in (10, 19.5, 20.5, 10, 10)]
  28.     for inc, thing in zip(incs, (venus, moon_p, moon_r, mars, jupiter)):
  29.         sinc, cinc = np.sin(inc), np.cos(inc)
  30.         thing1 = thing.copy()
  31.         thing1[0] = cinc * thing[0] - sinc * thing[2]
  32.         thing1[2] = sinc * thing[0] + cinc * thing[2]
  33.         thing[[0, 2]] = thing1[[0, 2]]
  34.  
  35.     moon_p +=  earth
  36.     moon_r +=  earth
  37.  
  38.     if True:
  39.         xe, ye, ze = earth
  40.         to_degs = 180 / np.pi
  41.         for i, body in enumerate((venus, moon_p, moon_r, mars, jupiter)):
  42.             d = body - earth
  43.             lon = np.arctan2(d[1], d[0])
  44.             lat = np.arctan2(d[2], np.sqrt(d[0]**2 + d[1]**2))
  45.             bk1 = (lon[1:] > 2) * (lon[:-1] < -2)
  46.             bk2 = (lon[:-1] > 2) * (lon[1:] < -2)
  47.             bk = bk1 + bk2
  48.             lon, lat = lon[:-1], lat[:-1]
  49.             lon[bk] = np.nan
  50.             plt.plot(to_degs * lon, to_degs * lat, linewidth=0.75)
  51.             plt.xlim(-180, 180)
  52.             plt.ylim(-90, 90)
  53.             plt.gca().set_aspect('equal')
  54.         plt.show()
  55.  
RAW Paste Data