import astropy.units as u 
import astropy.coordinates as c
import astropy.time as t
import numpy as np
import math

boston = c.EarthLocation.from_geodetic(-71.0833*u.deg,42.3333*u.deg)
midnight = t.Time( '1988-03-20', scale='utc')
h0 = -0.5667 *u.deg
venus = c.SkyCoord(ra=41.73129*u.deg, dec=18.44092*u.deg)

sidereal_time = c.Angle( (11,50,58.10), unit=u.hour)
print sidereal_time - midnight.sidereal_time('apparent','greenwich')

m0 = (venus.ra - boston.longitude - sidereal_time)/(360*u.deg/u.sday)
m0+=1*u.sday
print m0 
print (m0 - 0.81965*u.sday).to(u.s)

def fix_day(day) : 
    if day > 1 * u.sday : 
        day -= 1 *u.sday
    elif day < 0 * u.sday : 
        day += 1 * u.sday
    return day


cos_H0 = (np.sin(h0) - np.sin(boston.latitude)*np.sin(venus.dec))/(np.cos(boston.latitude)*np.cos(venus.dec))
print cos_H0 # -0.3178735
H0 = np.arccos(cos_H0)
m1 = m0 - H0/(360*u.deg/u.sday)
m2 = m0 + H0/(360*u.deg/u.sday)
m = u.Quantity([fix_day(m0),fix_day(m1),fix_day(m2)])
ref_m = [0.81965, 0.51817, 0.12113] * u.sday
print H0.to(u.deg) # 108.5344
print m
print (m - ref_m)
print (m - ref_m).to(u.s)

#m_solar =  m * 1.0027379093604878 * u.day/u.sday
m_solar =  m * 360.985647 * u.deg/u.sday 
#m_solar = m_solar * u.day/u.sday
theta0 = sidereal_time + m_solar
theta0.wrap_at(360*u.deg, inplace=True)
print theta0.deg
ref_theta0 = [113.62397, 4.79401, 221.46827] * u.deg
print (theta0 - ref_theta0).deg
print (theta0 - ref_theta0).to(u.arcsec)

cr_m_solar = ref_m * 360.985647 * u.deg/u.sday
cr_theta0 = sidereal_time + cr_m_solar
cr_theta0.wrap_at(360*u.deg, inplace=True)
print (cr_theta0-ref_theta0).deg
print (cr_theta0-ref_theta0).to(u.arcsec)

test_theta0 = sidereal_time + 360.985647 * 0.12113 * u.deg
print (test_theta0 - 221.46827 * u.deg).to(u.arcsec)

venus_interp = c.SkyCoord(ra=[42.59324,42.27648,41.85927], dec=[0,18.64229,18.48835], unit=u.deg)
print venus_interp
H = theta0 + boston.longitude - venus_interp.ra
print H.deg
ref_H = [-0.05257, -108.56577, 108.52570] * u.deg
print (H - ref_H).deg
print (H - ref_H).to(u.arcsec)


cr_H = c.Angle(ref_theta0 + boston.longitude - venus_interp.ra)
print cr_H.deg
print (cr_H - ref_H).deg
print (cr_H - ref_H).to(u.arcsec)

sin_h = np.sin(boston.latitude)*np.sin(venus_interp.dec) + (np.cos(boston.latitude)*np.cos(venus_interp.dec)*np.cos(H))
h = np.arcsin(sin_h)
print h.to(u.deg)
ref_h = [90, -0.44393, -0.52711] * u.deg
print (h - ref_h).to(u.deg)
print (h - ref_h).to(u.arcsec)

cr_sin_h = np.sin(boston.latitude)*np.sin(venus_interp.dec) + (np.cos(boston.latitude)*np.cos(venus_interp.dec)*np.cos(ref_H))
cr_h = np.arcsin(cr_sin_h)
print cr_h.to(u.deg)
print (cr_h - ref_h).to(u.deg)
print (cr_h - ref_h).to(u.arcsec)

dm0 = - (H[0] / (360*(u.deg/u.sday)))
dm1 = (h[1]-h0)/(360*u.deg/u.sday * np.cos(venus_interp[1].dec) * np.cos(boston.latitude) * np.sin(H[1]))
dm2 = (h[2]-h0)/(360*u.deg/u.sday * np.cos(venus_interp[2].dec) * np.cos(boston.latitude) * np.sin(H[2]))
delta_m = u.Quantity([dm0, dm1, dm2])
ref_delta_m = [0.00015, -0.00051, 0.00017] *u.sday
print delta_m.to(u.sday)
print (delta_m - ref_delta_m).to(u.min)
print (delta_m - ref_delta_m).to(u.s)

dm0 = - (ref_H[0] / (360*(u.deg/u.sday)))
dm1 = (ref_h[1]-h0)/(360*u.deg/u.sday * np.cos(venus_interp[1].dec) * np.cos(boston.latitude) * np.sin(ref_H[1]))
dm2 = (ref_h[2]-h0)/(360*u.deg/u.sday * np.cos(venus_interp[2].dec) * np.cos(boston.latitude) * np.sin(ref_H[2]))
cr_delta_m = u.Quantity([dm0, dm1, dm2])
print cr_delta_m.to(u.min)
print (cr_delta_m - ref_delta_m).to(u.min)
print (cr_delta_m - ref_delta_m).to(u.s)