Guest User

Untitled

a guest
Feb 25th, 2019
79
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.     GMe = 3.986E+14 # m^3/s^2
  4.     Re  = 6378137.  # meters
  5.  
  6.     halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]
  7.     degs, rads        = 180/pi, pi/180
  8.     Nhalf             = 200
  9.     theta             = np.linspace(0, twopi, 2*Nhalf+1)
  10.     peri, apo         = 1.1, 6
  11.     a, ep             = 0.5*(apo + peri), (apo-peri)/(apo+peri)
  12.     rcirc             = 1.1
  13.     rellipse          = a * (1-ep**2)/ (1-ep*np.cos(theta))
  14.     x0, y0            = [1.0   * f(theta) for f in (np.cos, np.sin)]
  15.     xc, yc            = [rcirc   * f(theta) for f in (np.cos, np.sin)]
  16.     xe, ye            = [rellipse* f(theta) for f in (np.cos, np.sin)]
  17.  
  18.     if True:
  19.         plt.plot(x0, y0, '-b', linewidth=2)
  20.         plt.plot(xc, yc, '-k')
  21.         plt.plot(xe, ye, '-k')
  22.         plt.plot(xe[::Nhalf], ye[::Nhalf], 'ok', markersize=8)
  23.         plt.xlim(-2, 7)
  24.         plt.show()
  25.  
  26.     rrel = np.arange(60) + 1.0
  27.     r    = Re * rrel
  28.  
  29.     vorb    = np.sqrt(GMe/r)
  30.     vesc    = np.sqrt(2) * vorb
  31.     aa      = 0.5*(r + Re*peri)
  32.     vapo    = np.sqrt(GMe*(2./r - 1./aa))
  33.     vexcrit = vapo + vesc
  34.  
  35.     if True:
  36.         plt.figure()
  37.         things = (vesc, vapo, vexcrit)
  38.         names  = ('vesc', 'vapo', 'vexcrit')
  39.         rplot = r * 0.001 # km
  40.         ilab  = 25
  41.         for thing, name in zip(things, names):
  42.             plt.plot(rplot, thing)
  43.             plt.text(rplot[ilab], thing[ilab]+100, name, fontsize=14)
  44.         plt.ylim(None, 6000)
  45.         plt.xlabel('r_apo (km)', fontsize=18)
  46.         plt.ylabel('velocity (m/s)', fontsize=18)
  47.         for v in (2000, 3000, 4000):
  48.             plt.plot(rplot, v*np.ones_like(rplot), '-k', linewidth=0.5)
  49.        
  50.         plt.show()
RAW Paste Data