TimeDilation.py

import math
import matplotlib.pyplot as plt
from numpy import *
from scipy.optimize import brentq


########################
# Time Dilation on ISS #
########################

# ---------
# Constants
# ---------

G       = 6.67408E-11           # universal gravitational constant in m^3 kg^-1 s^-2
c       = 299792458.0           # speed of light in m s^-1
M       = 5.9742E24             # Mass of the earth in kg
Rp      = 6356752.0             # Radius of earth at poles in m
Req     = 6378137.0             # radius of earth at equator in m
h       = 405000.0              # average height of ISS above earth's surface in m
Rs      = 2*G*M/c**2            # Schwarzschild radius for earth--radius for earth to be blackhole
period  = 86162.4               # sidereal period of the earth in s
v_eq    = 2*math.pi*Req/period  # velocity of earth at equator in m s^-1
N       = 342                   # number of days Scott Kelly spent on ISS


# ------------------------------------------------------------------------------------------
# Mathematics of Satellite Motion
# A satellite's linear speed, acceleration and period
# can be computed from its distance R from earth
# See http://www.physicsclassroom.com/class/circles/Lesson-4/Mathematics-of-Satellite-Motion
# ------------------------------------------------------------------------------------------

# Orbital Speed Equation
def v(R):
    return math.sqrt( G * M / R)

# The Acceleration Equation
def a(R):
    return G * M / R**2

# Orbital Period Equation
def T(R):
    return math.sqrt(R**3 * 4 * math.pi**2 / (G * M))

v_ISS = v(Req + h) # velocity of ISS above equator, approx 7660 m s^-1


# ------------------
# General Relativity
# ------------------

# Using the Schwarzschild metric, compute $\frac{d\tau}{dt}$:
def dtau_dt(R, v):
    beta = v/c # relative velocity, compared to light
    return math.sqrt((1 - Rs/R) - (1 - Rs/R)**-1 * beta**2)

# 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.
def GR(R):
    return 10.0**6 * period * (dtau_dt(R,v(R)) / dtau_dt(Req,v_eq) - 1)

GR_ISS = GR(Req + h)
print("* Daily time dilation (gain or loss if negative) in microseconds of ISS according to GR (gravity + velocity) = %f" % (GR_ISS))


# ------------------
# Special Relativity
# ------------------

# The GR calculation above already comprises
# 1. Special Relativistic (SR) effects due to velocity only
# 2. additional effects due to gravity only

# In this section, we compute the SR effects separately.
# Then we subtract them from the GR effects, to find the effects of gravity alone...

# alpha is the reciprocal of gamma---the Lorentz contraction factor
def alpha(v):
    beta = v/c # relative velocity, compared to light
    return math.sqrt(1 - beta**2)

# Daily time dilation (gain or loss if negative) in microseconds due to orbital velocity only (SR):
def SR(R):
    return 10.0**6 * period * (alpha(v(R)) - 1)

print("* Daily time dilation (gain or loss if negative) in microseconds due to ISS's velocity only (SR) = %f" % (SR(Req + h)))

# Effect of gravity only (GR - SR)
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)))

# ------------------
# Scott Kelly
# ------------------

# 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):
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)))


# ------------------
# Zero Time Dilation
# ------------------

# At what orbital radius do the effects of gravity and velocity cancel each other out, i.e., R s.t. GR(R) = 0?

# Helper function to compute time dilation as multiple of earth's equatorial radius
def mult_to_GR(mult):
    return GR(mult*Req)

# Use scipy to find the zero of the GR function...
one_root = brentq(mult_to_GR, 1, 2)

# print result
print("\nRoot of GR time dilation function at %f radii in interval [%f,%f] " % (one_root, 1, 2))

# check solution
print("Daily time dilation (gain or loss if negative) in microseconds at R = %f * Req is %f" % (one_root, GR(one_root*Req)))


# ------------------
# Plot
# ------------------

factor = linspace(1,7,400)

fig = plt.figure()
fig.suptitle('Time Dilation', fontsize=14, fontweight='bold')
ax = fig.add_subplot(111)
ax.set_xlabel('Multiples of earth equatorial radius')
ax.set_ylabel('Daily time dilation (gain or loss if negative) in microseconds')

line1, = plt.plot(factor,[mult_to_GR(f) for f in factor], color='g')
line2, = plt.plot(factor,[0 for f in factor], color='k')
line3, = plt.plot(factor,[GR_ISS for f in factor], color='r', label="ISS")

ISS_legend = plt.legend(handles=[line3], loc=4)
ax = plt.gca().add_artist(ISS_legend)

plt.show()