import numpy as np
import matplotlib.pyplot as plt
##############################################
# Code is functioning; ignores annual motion #
# but was a useful exercise. Not sure why it #
# doesn't match with Max's graphs exactly w/ #
# phase, but it's not really important. #
# - 17 Jun 2016 #
##############################################
# This code assumes that the Earth is moving parallel to
# the source in its annual motion. Year implemented in time.
# Constants (in SI units)
c = 299792458 # Speed of light
rE = 6371000 # Radius of Earth
rOrb = 149600000000 # Radius of Earth's orbit
Td = 86400 # Seconds in a day
# Position of the observatory (e.g. Baton Rouge)
lat0 = 30.4583
# Position of the source (e.g. Crab Nebula)
alp0 = 5.5755389 # right ascension is in hours
dlt0 = 22.0145 # Declination is in degrees
# Convert to radians
alp = alp0 * np.pi / 12
dlt = dlt0 * np.pi / 180
lat = lat0 * np.pi / 180
secs = np.arange(0,Td)
nHat = np.array([[np.cos(alp) * np.cos(dlt)],
[np.sin(alp) * np.cos(dlt)],
[np.sin(dlt)]])
def r_orbit(t):
'''Gives rectangular coordinates of planet in orbit;
we will assume roughly constant over a day. '''
return rOrb * np.array([[0], [0], [0]])
def r_revol(t):
'''Gives rectangular coordinates of observatory relative center of Earth'''
lngi = 2 * np.pi * t / Td # Angle over time
return rE * np.array([[np.cos(lat) * np.cos(lngi)],
[np.sin(lat) * np.cos(lngi)],
[np.sin(lngi)]])
def r(t):
return r_orbit(t) + r_revol(t)
def romerDelay(t):
return r(t).T.dot(nHat) / c
# I want to get rid of this loop, but for now it works; must deal with later.
romer = np.zeros(Td)
for s in secs:
romer[s] = romerDelay(s)
plt.title("Romer Delay over 1 day")
plt.xlabel("time (s)")
plt.ylabel("Delay (s)")
plt.plot(romer)
plt.show()