• API
• FAQ
• Tools
• Archive
SHARE
TWEET

# TimeDilation.py

krusader74 Mar 7th, 2016 81 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. import math
2. import matplotlib.pyplot as plt
3. from numpy import *
4. from scipy.optimize import brentq
5.
6.
7. ########################
8. # Time Dilation on ISS #
9. ########################
10.
11. # ---------
12. # Constants
13. # ---------
14.
15. G       = 6.67408E-11           # universal gravitational constant in m^3 kg^-1 s^-2
16. c       = 299792458.0           # speed of light in m s^-1
17. M       = 5.9742E24             # Mass of the earth in kg
18. Rp      = 6356752.0             # Radius of earth at poles in m
19. Req     = 6378137.0             # radius of earth at equator in m
20. h       = 405000.0              # average height of ISS above earth's surface in m
21. Rs      = 2*G*M/c**2            # Schwarzschild radius for earth--radius for earth to be blackhole
22. period  = 86162.4               # sidereal period of the earth in s
23. v_eq    = 2*math.pi*Req/period  # velocity of earth at equator in m s^-1
24. N       = 342                   # number of days Scott Kelly spent on ISS
25.
26.
27. # ------------------------------------------------------------------------------------------
28. # Mathematics of Satellite Motion
29. # A satellite's linear speed, acceleration and period
30. # can be computed from its distance R from earth
31. # See http://www.physicsclassroom.com/class/circles/Lesson-4/Mathematics-of-Satellite-Motion
32. # ------------------------------------------------------------------------------------------
33.
34. # Orbital Speed Equation
35. def v(R):
36.     return math.sqrt( G * M / R)
37.
38. # The Acceleration Equation
39. def a(R):
40.     return G * M / R**2
41.
42. # Orbital Period Equation
43. def T(R):
44.     return math.sqrt(R**3 * 4 * math.pi**2 / (G * M))
45.
46. v_ISS = v(Req + h) # velocity of ISS above equator, approx 7660 m s^-1
47.
48.
49. # ------------------
50. # General Relativity
51. # ------------------
52.
53. # Using the Schwarzschild metric, compute $\frac{d\tau}{dt}$:
54. def dtau_dt(R, v):
55.     beta = v/c # relative velocity, compared to light
56.     return math.sqrt((1 - Rs/R) - (1 - Rs/R)**-1 * beta**2)
57.
58. # Daily time dilation (gain or loss if negative) in microseconds as a function of orbit radius R. At R = 1.49741*Req approx. there is no time dilation.
59. def GR(R):
60.     return 10.0**6 * period * (dtau_dt(R,v(R)) / dtau_dt(Req,v_eq) - 1)
61.
62. GR_ISS = GR(Req + h)
63. print("* Daily time dilation (gain or loss if negative) in microseconds of ISS according to GR (gravity + velocity) = %f" % (GR_ISS))
64.
65.
66. # ------------------
67. # Special Relativity
68. # ------------------
69.
70. # The GR calculation above already comprises
71. # 1. Special Relativistic (SR) effects due to velocity only
72. # 2. additional effects due to gravity only
73.
74. # In this section, we compute the SR effects separately.
75. # Then we subtract them from the GR effects, to find the effects of gravity alone...
76.
77. # alpha is the reciprocal of gamma---the Lorentz contraction factor
78. def alpha(v):
79.     beta = v/c # relative velocity, compared to light
80.     return math.sqrt(1 - beta**2)
81.
82. # Daily time dilation (gain or loss if negative) in microseconds due to orbital velocity only (SR):
83. def SR(R):
84.     return 10.0**6 * period * (alpha(v(R)) - 1)
85.
86. print("* Daily time dilation (gain or loss if negative) in microseconds due to ISS's velocity only (SR) = %f" % (SR(Req + h)))
87.
88. # Effect of gravity only (GR - SR)
89. print("* Daily time dilation (gain or loss if negative) in microseconds of ISS due to gravity only (GR - SR) = %f" % (GR_ISS - SR(Req + h)))
90.
91. # ------------------
92. # Scott Kelly
93. # ------------------
94.
95. # Number of microseconds older (+) / younger (-) Scott Kelly is, compared to his twin, Mark, after N days on the ISS as a result of GR (gravity + velocity):
96. print("\nNumber of microseconds older (+) / younger (-) Scott Kelly is, compared to his twin, Mark, after %d days on the ISS as a result of GR (gravity + velocity) = %f" % (N, N * (GR_ISS)))
97.
98.
99. # ------------------
100. # Zero Time Dilation
101. # ------------------
102.
103. # At what orbital radius do the effects of gravity and velocity cancel each other out, i.e., R s.t. GR(R) = 0?
104.
105. # Helper function to compute time dilation as multiple of earth's equatorial radius
106. def mult_to_GR(mult):
107.     return GR(mult*Req)
108.
109. # Use scipy to find the zero of the GR function...
110. one_root = brentq(mult_to_GR, 1, 2)
111.
112. # print result
113. print("\nRoot of GR time dilation function at %f radii in interval [%f,%f] " % (one_root, 1, 2))
114.
115. # check solution
116. print("Daily time dilation (gain or loss if negative) in microseconds at R = %f * Req is %f" % (one_root, GR(one_root*Req)))
117.
118.
119. # ------------------
120. # Plot
121. # ------------------
122.
123. factor = linspace(1,7,400)
124.
125. fig = plt.figure()
126. fig.suptitle('Time Dilation', fontsize=14, fontweight='bold')
128. ax.set_xlabel('Multiples of earth equatorial radius')
129. ax.set_ylabel('Daily time dilation (gain or loss if negative) in microseconds')
130.
131. line1, = plt.plot(factor,[mult_to_GR(f) for f in factor], color='g')
132. line2, = plt.plot(factor,[0 for f in factor], color='k')
133. line3, = plt.plot(factor,[GR_ISS for f in factor], color='r', label="ISS")
134.
135. ISS_legend = plt.legend(handles=[line3], loc=4)