SHARE
TWEET

TimeDilation.py

krusader74 Mar 7th, 2016 71 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')
  127. ax = fig.add_subplot(111)
  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)
  136. ax = plt.gca().add_artist(ISS_legend)
  137.  
  138. plt.show()
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top