Guest User

Untitled

a guest
Aug 30th, 2018
153
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. def deriv(X, t):
  2.     x, v = X.reshape(2, -1)
  3.     acc = -GM * x * ((x**2).sum())**-1.5
  4.     return np.hstack((v, acc))
  5.  
  6. import numpy as np
  7. import matplotlib.pyplot as plt
  8. from scipy.integrate import odeint as ODEint
  9.  
  10. # https://physics.stackexchange.com/questions/348854/parker-solar-probe-passing-extremely-close-to-the-sun-what-relativistic-effects
  11.  
  12. halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]
  13.  
  14. c    = 2.9979E+08
  15. a    = 57.8E+06 * 1000.
  16. GM   = 1.327E+20
  17. peri = 6.6E+06 * 1000.
  18. apo  = 2*a - peri
  19. vapo, vperi = [np.sqrt(GM*(2./r - 1./a)) for r in (apo, peri)]
  20.  
  21. T    = twopi * np.sqrt(a**3/GM)
  22.  
  23. X0   = np.hstack([apo, 0, 0, vapo])
  24.  
  25. time = np.linspace(0, T, 1000001)
  26. days = time/(24.*3600.)
  27.  
  28. answer, info = ODEint(deriv, X0, time, full_output=True, rtol=1E-10)
  29.  
  30. x, v = answer.T.reshape(2, 2, -1)
  31. r    = np.sqrt((x**2).sum(axis=0))
  32.  
  33. dfof = -(GM/c**2) * (2./r - 0.5/a)
  34. dt   = time[1] - time[0]
  35.  
  36. delta_t = dt * dfof
  37. DT      = delta_t.cumsum()
  38.  
  39. if True:
  40.     plt.figure()
  41.     plt.subplot(2, 1, 1)
  42.     plt.plot(x[0]/1000., x[1]/1000.)
  43.     plt.plot([0], [0], 'ok')
  44.     plt.xlim(-0.1E+08, 1.1E+08)
  45.     plt.ylim(-0.6E+08, 0.6E+08)
  46.     plt.title('heliocentric orbit (km)', fontsize=14)
  47.     plt.subplot(4, 1, 3)
  48.     plt.plot(days, dfof)
  49.     plt.title('delta f/f vs time (days)', fontsize=14)
  50.     plt.subplot(4, 1, 4)
  51.     plt.plot(days, DT)
  52.     plt.title('cumulative time (sec) vs time (days)', fontsize=14)
  53.     plt.show()
RAW Paste Data