# 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