import time, math
#
# Find the azimuth and elevation of an object
#
def rad_to_deg(rad):
    "Converts radians to degrees"
    return (rad*180.0)/math.pi

def deg_to_rad(deg):
    "Converts degrees to radians"
    return (deg/180.0)*math.pi

def eq_to_hor(ha_angle,dec_deg,lat_deg):
#
# returns azimuth and altitude
#
# from starlink slalib  sla_DE2H
#
    ha = deg_to_rad(15.*ha_angle)
    dec= deg_to_rad(dec_deg)
    phi= deg_to_rad(lat_deg)
    sh = math.sin(ha)
    ch = math.cos(ha)
    sd = math.sin(dec)
    cd = math.cos(dec)
    sp = math.sin(phi)
    cp = math.cos(phi)
    x  = -1.*ch*cd*sp + sd*cp
    y  = -1.*sh*cd
    z  = ch*cd*cp + sd*sp
    r  = math.sqrt(x*x+y*y)
    if r == 0.0:
       a = 0.0
    else:
       a = math.atan2(y,x)
    if a < 0.0: a = a + 2.*math.pi
    az_rad = a
    el_rad = math.atan2(z,r)
    az = rad_to_deg(az_rad)
    el = rad_to_deg(el_rad)
#    print "---------------->",ha_angle, dec_deg, lat_deg, az,el
    return az,el


def find_lst(tel_long_deg,year,month,mday,hour,min,sec):
#
# returns local sidereal time give the observer's longitude
#
    (sid_hour,sid_min,sid_sec) = find_sid(year,month,mday,hour,min,sec)
    sid_dhour = float(sid_hour) + float(sid_min)/60.0 + float(sid_sec)/3600.
    lst_hours = sid_dhour - tel_long_deg/15.
    if lst_hours > 24.0:
        lst_hours = lst_hours - 24.0
    if lst_hours < 0.0:
        lst_hours = lst_hours + 24.0
 
    return lst_hours
 
def find_sid(year,month,mday,hour,min,sec):
#
# returns greenwich sidereal time
#
# from Duffet-Smith's book "Practical Astronomy with your Calculator"
#
    ut_hours = hour + (min + sec/60.)/60.
    mday = mday + ut_hours/24.
 
    if  month <= 2 :
        year = year - 1
        month = month + 12
 
    a = int(year/100.)
    b = 2. - a + int(a/4.)
 
    if year < 0 :
        c = int((365.2500*year)-0.7500)
    else:
        c = int(365.2500*year)
 
    d = int(30.600100*(month+1))
    jd = b + c + d + int(mday) + 1720994.500
 
    s = jd - 2451545.000
    t = s/36525.000
    t0 = 6.697374558 + (2400.051336*t) + (0.000025862*(t*t));
    t0 = (t0 - int(t0/24.)*24)
    if t0 < 0.0: t0 = t0 + 24.
    ut = 1.002737909*ut_hours
    tmp = int((ut + t0)/24.)
    gst = ut + t0 - tmp*24.
    gst_hour = int(gst)
    gst_min  = int((gst - gst_hour)*60.)
    gst_sec  = int((gst - gst_hour - gst_min/60.)*3600.)
    return (gst_hour, gst_min, gst_sec)

##########################################################################

DURHAM_LONG_DEG =   1.583
DURHAM_LAT_DEG  =  54.767

year  = 2007
month = 11
mday = 29
hour = 17
min = 54
sec  = 0
print " "
print "Positions at ", year, month, mday, hour, min, sec

lst_hr = find_lst(DURHAM_LONG_DEG,year,month,mday,hour,min,sec)
print "Local Sidereal Time is ", lst_hr, "hours"

print " "
print "Vega"
ra_hours  = 18. + 36./60. + 56.3/3600.
dec_deg   = 38. + 47./60. + 01./3600.
print "RA is  ", ra_hours," hours"
print "Dec is ",dec_deg, "deg"
hour_angle = lst_hr - ra_hours
(az,el) = eq_to_hor(hour_angle,dec_deg,DURHAM_LAT_DEG)
print "Azimuth is   ", az," deg"
print "Elevation is ", el, "deg"

print " "
print "Mirfak"
ra_hours  = 03. + 24./60. + 19.4/3600.
dec_deg   = 49. + 51./60. + 40./3600.
print "RA is  ", ra_hours," hours"
print "Dec is ",dec_deg, "deg"
hour_angle = lst_hr - ra_hours
(az,el) = eq_to_hor(hour_angle,dec_deg,DURHAM_LAT_DEG)
print "Azimuth is   ", az," deg"
print "Elevation is ", el, "deg"

print " "
print "Altair"
ra_hours  = 19. + 50./60. + 47.0/3600.
dec_deg   =  8. + 52./60. + 6./3600.
print "RA is  ", ra_hours," hours"
print "Dec is ",dec_deg, "deg"
hour_angle = lst_hr - ra_hours
(az,el) = eq_to_hor(hour_angle,dec_deg,DURHAM_LAT_DEG)
print "Azimuth is   ", az," deg"
print "Elevation is ", el, "deg"

print " "
print "Dubhe"
ra_hours  = 11. + 03./60. + 43.7/3600.
dec_deg   = 61. + 45./60. + 3.7/3600.
print "RA is  ", ra_hours," hours"
print "Dec is ",dec_deg, "deg"
hour_angle = lst_hr - ra_hours
(az,el) = eq_to_hor(hour_angle,dec_deg,DURHAM_LAT_DEG)
print "Azimuth is   ", az," deg"
print "Elevation is ", el, "deg"