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