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()