﻿

# Untitled

a guest
Feb 25th, 2019
87
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