from slalib import *
import numpy as np
import string as st
import time as tm
# http://www.python.org/doc/current/lib/socket-example.html
import socket
from control20m import *
ant = Control20m()
ant.take_control()
ant.sync_unix_time()

USE_SKT = False
#USE_SKT = True
if USE_SKT :
    #HOST = '192.33.116.227'
    #HOST = '10.16.96.163'
    HOST = 'byu20meter.gb.nrao.edu'
    PORT = 5010
    socket.setdefaulttimeout(10)
    _skt = socket.socket(socket.AF_INET, socket.SOCK_STREAM)
    _skt.connect((HOST, PORT))
#_skt.close()

print '->->-> Need to look at settling time. <-<-<-'

def take_antenna_control ( ) :
    """
    take_antenna_control()
    Put the 20-meter antenna under Python program control.
    """
    ant.take_control()
    
def release_antenna ( ) :
    """
    release_antenna()
    Release the 20-meter antenna from Python program control.
    """
    ant.release_control()

def hadec_to_azel ( ha_hrs, dec_deg ) :
    """
    hadec_to_azel ( ha_hrs, dec_deg )
    Convert Hour Angle - Declination coordinates to Azimuth and Elevation
    without making any pointing or refraction corrections.  HA is in decimal
    hours and Dec, Az, and El are in decimal degrees.
    The returned value is a tuple (az, el)
    """
    ha_rad = np.pi * ha_hrs / 12.0
    dec_rad = np.pi * dec_deg / 180.0
    vec = sla_dcs2c(ha_rad, dec_rad)
    # use the latitude of the 140-ft, which is nearly due west of the 20-m
    desig, name, of_long, of_lat, elev = sla_obs(0, 'GBVA140')
    #print of_lat * 180 / np.pi
    rmat = sla_deuler('yxx', np.pi / 2.0 - of_lat, 0.0, 0.0)
    new_vec = sla_dmxv(rmat, vec)
    az, el = sla_dcc2s(new_vec)
    return (180.0 * (1.0 + az / np.pi), 180.0 * el / np.pi)

#print hadec_to_azel(0.0, 0.0)
#print hadec_to_azel(-1.0, 50.0)

def str2frac ( dt_str ) :
    """
    str2frac(dt_str)
    Convert a colon separated hexagesimal string to decimal.
    'xx:yy:zz.zz' -> xx.xxxx, xx can be signed
    """
    parsed = st.split(dt_str, ':')
    sign = 1.0
    if parsed[0][0] == '-' :
        sign = -1.0
    return sign * (abs(st.atof(parsed[0])) + (st.atof(parsed[1]) +
                                 st.atof(parsed[2]) / 60.0) / 60.0)

#print str2frac('11:22:33')
#print str2frac('-11:22:33')

def j2000_to_epoch ( ra_str, dec_str, mjd ) :
    """
    j2000_to_epoch('HH:MM:SS.ss', '-DD:MM:SS.s', mjd)
    Convert hexagesimal J2000 R.A.and Dec. string values to R.A.and Dec
    for the given Modified Julian Date using J2000 precession and nutation.
    The returned value is a tuple (ra, dec), where ra is in decimal hours,
    and dec is in decimal degrees.
    """
    rmat = sla_prenut(2000, mjd)
    ra = np.pi * str2frac(ra_str) / 12.0
    dec = np.pi * str2frac(dec_str) / 180.0
    vec = sla_dcs2c(ra, dec)
    new_vec = sla_dmxv(rmat, vec)
    pra, pdec = sla_dcc2s(new_vec)
    if pra < 0.0 :
        pra += 2.0 * np.pi
    return (12.0 * pra / np.pi, 180.0 * pdec / np.pi)

#print 'Precessed:', j2000_to_epoch('1:00:00', '1:00:00', 54630)

def get_current_mjd ( ) :
    """
    get_current_mjd()
    Reads the system clock and converts to fractional Modified Julian Date
    with fractional second precision to the accuracy of the system clock.
    """
    frac_sec1 = 1.0; frac_sec2 = 0.0
    while (frac_sec1 > frac_sec2) :
        ftime1 = tm.time()
        gmt = tm.gmtime()
        ftime2 = tm.time()
        frac_sec1 = ftime1 - int(ftime1)
        frac_sec2 = ftime2 - int(ftime2)
    #print gmt
    yr = gmt[0]
    mjd = 54465 + gmt[7] + (gmt[3] + (gmt[4] + gmt[5] / 60.0) / \
                             60.0) / 24.0 + (frac_sec1 + frac_sec2) / 172800.0
    for i in range(2009, 2052) :
        if i > yr :
            break
        if (i - 1) % 4 == 0 :
            mjd += 366
        else :
            mjd += 365
    return mjd

#print get_current_mjd()

_RH = 0.8; _ABS_TEMP = 290.0; _PRESS = 1000
print 'RH:', _RH * 100, '%,', '_ABS_TEMP:', _ABS_TEMP, 'K,', \
      'Pressure:', _PRESS, 'mBar'
_REFA, _REFB = sla_refco(881.0, _ABS_TEMP, _PRESS, _RH,
                         1000.0, 0.67086, 0.0065, 1e-10)

def get_refracted_el ( el ) :
    """
    get_refracted_el(el)
    Returns the apparent elevation of an object after refraction
    correction based on the relative humidity, absolute temperature,
    barometric pressure and other site parameters hardwired into this code.
    The elevation is in decimal degrees.
    """
    zu = np.pi * (90.0 - el) / 180.0
    zr = sla_refz(zu, _REFA, _REFB)
    return 90.0 - 180.0 * zr / np.pi

#print get_refracted_el(10.0)

def get_corrected_azel ( az, el, with_refraction=True ) :
    """
    get_corrected_azel(az, el, with_refraction=True)
    Returns azimuth and elevation in decimal degrees corrected for
    telescope pointing errors and, if the with_refraction parameter is True,
    for atmospheric refraction.  The telescope pointing error coefficients
    and equations are hardwired into code in the file, control20m.py.  The
    input parameters are true or apparent azimuth and elevation, and the
    output values are azimuth and elevation to be used as telescope position
    commands.
    """
    if with_refraction :
        rel = get_refracted_el(el)
    else :
        rel = el
    #encoder_az = az
    encoder_az = az + ant.get_az_pointing_offset(az, rel)
    if encoder_az >= 360.0 :
        encoder_az -= 360.0
    if encoder_az < 0.0 :
        encoder_az += 360.0
    #encoder_el = rel
    encoder_el = rel + ant.get_el_pointing_offset(az, rel)
    return (encoder_az, encoder_el)

#print get_corrected_azel(0.0, 10.0, True)

def get_object_posn ( object ) :
    """
    get_object_posn(object)
    Returns the RA and Dec for the current epoch corrected for precession
    and nutation.  The input parasmeter may be one of two source name
    strings, 'CasA', 'CygA', or a J2000 hexagesimal string value pair in
    the format 'HH:MM:SS.ss,-DD:MM:SS.s'  The returned value is a tuple
    (ra, dec) in decimal hours and decimal degrees.
    """
    if object == 'CasA' :
        object = '23:23:26.7,+58:49:03'
    elif object == 'CygA' :
        object = '19:59:28.3,+40:44:02'
    elif object == 'W3OH' :
        object = '02:27:04.1, 61:52:22.0'
    elif object == 'W49N' :
        object = '19:10:13.2, 09:06:12.0'
    pos = st.split(object, ',')
    if len(pos) != 2 :
        print 'get_objet_posn(): Object position error', object
        return (False, False)
    return j2000_to_epoch(pos[0], pos[1], get_current_mjd())

#print get_object_posn('11:22:33,-5:33:77')

def map_azel_to_azel ( map_az_deg, map_el_deg, src_az_deg, src_el_deg ) :
    """
    map_azel_to_azel ( map_az_deg, map_el_deg, src_az_deg, src_el_deg )
    This routine is used in mapping a telescope beam pattern, where the beam
    coordinates are specified in a near-rectangular grid (actually spherical
    coordinates, azimuth and elevation, centered on az=0.0, el=0.0).  The
    last two parameters are the apparent azimuth and elevation of the radio
    source, without telecope pointing corrections, being used to map the
    beam.  The first two parameters specify where this radio source is to be
    placed in the beam map.  All four parameters are in decimal degrees.
    The returned value is a tuple (az, el) of coordinates where the primary
    telescope axis is to be pointed.  Positive map_az drives the telescope
    to smaller azimuth, and positive map_el drives the telescope to lower
    elevation.  All arguments and return valuies are in degrees.  The
    returned (az, el) are corrected for atmospheric refraction but not
    telescope pointing errors.
    """
    # convert degrees to radians
    map_az_rad = np.pi * map_az_deg / 180.0
    map_el_rad = np.pi * map_el_deg / 180.0
    src_az_rad = np.pi * src_az_deg / 180.0
    src_el_rad = np.pi * src_el_deg / 180.0
    # convert map spherical coordinates into a 3-D cartesian unit vector
    vec = sla_dcs2c(map_az_rad, map_el_rad)
    #print vec
    # the offset telescope position is determined iteratively, starting
    # with a guess
    elg = src_el_rad - map_el_rad
    azg = src_az_rad - map_az_rad / np.cos(elg)
    # position tolerance of 2 arcsec in radians
    err_tol = 2.0 * np.pi / (180.0 * 3600.0)
    num_iter = 10  # maximum number of iterations
    for i in range(num_iter) :
        # determine the rotation matrix for moving the map origin to the
        # guessed azimuth and elevation
        rmat = sla_deuler('yzx', elg, -azg, 0.0)
        # do that coordinate rotation
        new_vec = sla_dmxv(rmat, vec)
        #print new_vec
        # covert the rotated map vector to az, el coordinates
        az, el = sla_dcc2s(new_vec)
        #print i, '%10.5f %10.5f %10.5f %10.5f' % \
        #      (azg * 180 / np.pi, elg * 180 / np.pi,
        #       az * 180 / np.pi, el * 180 / np.pi)
        # compute the differences between map position az, el and the source
        # az, el
        az_diff = src_az_rad - az
        if az_diff > np.pi :
            az_diff -= 2.0 * np.pi
        el_diff = src_el_rad - el
        #print '-->', az_diff, el_diff
        # if the differences are within tolerance, quit
        if abs(az_diff) < err_tol and abs(el_diff) < err_tol :
            break
        # if not, apply the differences for a new guess
        azg += az_diff
        elg += el_diff
        # stop, if taking too many iterations
        if i >= num_iter - 1 :
            print 'map_azel_to_azel(%4.2f, %4.2f, %4.2f, %4.2f) did not converge' % (map_az_deg, map_el_deg, src_az_deg, src_el_deg)
            break
    if azg < 0.0 :
        azg += 2.0 * np.pi
    if azg > 2.0 * np.pi :
        azg -= 2.0 * np.pi
    return {'tel_az':180.0 * azg / np.pi, 'tel_el':180.0 * elg / np.pi}

#for az in np.arange(0.0, 361, 5.0) :
#    for el in np.arange(0.0, 79.0, 5.0) :
#        map_azel_to_azel(10.0, -10.0, az, el)

#srcaz = 0.0; srcel = 5.0
#mapaz = 5.0; mapel = 5.0
#t = map_azel_to_azel(mapaz, mapel, srcaz, srcel)
#print 't:', t
#vec1 = sla_dcs2c(srcaz * np.pi / 180.0, srcel * np.pi / 180.0)
#print 'vec1:', vec1
#vec2 = sla_dcs2c(t['tel_az'] * np.pi / 180.0, t['tel_el'] * np.pi / 180.0)
#print 'vec2:', vec2
#dv = np.sqrt((vec1[0] - vec2[0])**2 + (vec1[1] - vec2[1])**2 + (vec1[2] - vec2[2])**2)
#print 'dv, acsin(dv):', dv, np.arcsin(dv / 2.0) * 360.0 / np.pi
#vec3 = sla_dcs2c(mapel * np.pi / 180.0, mapel * np.pi / 180.0)
#print 'vec3:', vec3
#dvm = np.sqrt((1.0 - vec3[0])**2 + (vec3[1])**2 + (vec3[2])**2)
#print 'dvm, acsin(dvm):', dvm, np.arcsin(dvm / 2.0) * 360.0 / np.pi

DUT1 = -0.43  # UT1 - UTC in seconds  UT1 = UTC + DUT1

def get_last_str ( mjd, utc_string ) :
    """
    get_last_str(mjd, utc_string)
    Converts the specified integer Modified Julian Date and hexagesimal
    UTC, 'HH:MM:SS.ss', to Local Apparent Sidereal Time at the 20-meter
    telescope using the DUT1 = UT1 - UTC value hardwired into this code.
    The current value of DUT1 can be found in IERS Bulletin A.
    http://maia.usno.navy,mil/ser7/ser7.dat
    The returned value is in fraction of a sidereal day.
    """
    utc_frac = str2frac(utc_string) / 24.0
    mjd_frac = mjd + utc_frac
    #desig, name, gbt_long, lat, elev = sla_obs(0, 'GBT')
    desig, name, of_long, lat, elev = sla_obs(0, 'GBVA140')
    lon_frac =  of_long / (2.0 * np.pi)
    gast_rad = sla_gmst(mjd_frac + DUT1 / 86400.0) + sla_eqeqx(mjd_frac)
    diff140to20m = 2.26  # seconds
    last = (gast_rad - of_long) / (2.0 * np.pi) + diff140to20m / 86400.0
    if last >= 1.0 :
        last -= 1.0
    if last < 0.0 :
        last += 1.0
    return last

#print (get_last_str(54630.0, '0:0:0') - 0.50511288448786651)  * 86400.0
#print (get_last_str(54630.0, '12:12:12') - 0.014977253263429011)  * 86400.0

def get_last ( mjd_frac ) :
    """
    get_last(mjd_frac)
    Converts the specified fractional Modified Julian Date (fractional part
    is UTC fraction of a day), to Local Apparent Sidereal Time at the 20-meter
    telescope using the DUT1 = UT1 - UTC value hardwired into this code.
    The current value of DUT1 can be found in IERS Bulletin A.
    http://maia.usno.navy,mil/ser7/ser7.dat
    The returned value is in fraction of a sidereal day.  To get the current
    LAST use get_last(get_current_mjd()).
    """
    #desig, name, gbt_long, lat, elev = sla_obs(0, 'GBT')
    desig, name, of_long, lat, elev = sla_obs(0, 'GBVA140')
    lon_frac =  of_long / (2.0 * np.pi)
    gast_rad = sla_gmst(mjd_frac + DUT1 / 86400.0) + sla_eqeqx(mjd_frac)
    diff140to20m = 2.26  # seconds
    return (gast_rad - of_long) / (2.0 * np.pi) + diff140to20m / 86400.0

#print (get_last(54630.5084722222222) - 0.014977253263429011)  * 86400.0
#print get_last(get_current_mjd())

def make_track ( object, start_mjd, stop_mjd, cross_el_offset, el_offset ) :
    """
    make_track(object, start_mjd, stop_mjd, cross_el_offset, el_offset)
    Creates a series of time, position and velocity values to be followed by
    the 20-meter telescope drive system for the purpose of mapping the
    telescope beam pattern.  The first input parameter is the radio source
    used to map the beam and can be one of three name strings, 'CasA', 'CygA',
    or 'GOES 12' or a J2000 coordinate pair as a hexagesimal string of the
    following format: 'HH:MM:SS.s,-DD:MM:SS.s'.  The second and third
    parameters are start and stop dates and times in fractional Modified
    Julian Date, where the fractional part is UTC day fraction.  The last two
    parameters specify the coordinates of the beam map where the radio source
    is to be located and tracked.  See the __doc__ string for
    map_azel_to_azel(). The offsets are in decimal degrees.  The returned
    value is a dictionary of equal-length numpy arrays of the form
    {'mjd':[], 'az':[], 'az_vel':[], 'el':[], 'el_vel':[]} in fractional MJD,
    decimal degrees, and degrees per second.  The (az, el) values are
    corrected for refraction and telescope pointing errors so may be used
    directly by the telescope drive system.  The last array values are not
    the last point in the track.  This is reached by continuing from the
    last position at the specified rates until the stop_mjd.
    """
    if object == 'GOES 12' :
        src_az = 180.0 - 7.75
        el = 45.20
        src_el = get_refracted_el(el)
        tel_azel = map_azel_to_azel(cross_el_offset, el_offset, src_az, src_el)
        tel_az, tel_el = get_corrected_azel(tel_azel['tel_az'],
                                            tel_azel['tel_el'], False)
        return {'mjd':np.array([start_mjd]),
                'az':np.array([tel_az]), 'az_vel':np.array([0.0]),
                'el':np.array([tel_el]), 'el_vel':np.array([0.0])}
        
    ra, dec = get_object_posn(object)
    if ra == False :
        return
    start_last = get_last(start_mjd)
    start_ha = 24.0 * start_last - ra
    if start_ha < -12.0 :
        start_ha += 24.0
    if start_ha > 12.0 :
        start_ha -= 24.0
    stop_last = get_last(stop_mjd)
    stop_ha = 24.0 * stop_last - ra
    if stop_ha < -12.0 :
        stop_ha += 24.0
    if stop_ha > 12.0 :
        stop_ha -= 24.0
    if stop_ha < start_ha :
        stop_ha += 24.0
    duration = (stop_ha - start_ha) * 3600.0
    min_steps = 1
    while duration / min_steps > 12.0 :
        min_steps *= 2
    
    if stop_last < start_last :
        stop_last += 1.0
    num_iter = 3
    max_num_pts = min_steps * (2**num_iter)
    pt_mjd = np.zeros(max_num_pts + 1, dtype=np.float64)
    pt_az = np.zeros(max_num_pts + 1, dtype=np.float64)
    pt_az_vel = np.zeros(max_num_pts, dtype=np.float64)
    pt_el = np.zeros(max_num_pts + 1, dtype=np.float64)
    pt_el_vel = np.zeros(max_num_pts, dtype=np.float64)
    posn_tol = 2.0 / 3600.0  # degrees
    for iter in range(num_iter) :
        num_pts = min_steps * (2**(iter + 1))
        ha_step = (stop_ha - start_ha) / num_pts
        for pt in range(num_pts + 1) :
            src_az, el = hadec_to_azel(start_ha + ha_step * pt, dec)
            src_el = get_refracted_el(el)
            tel_azel = map_azel_to_azel(cross_el_offset, el_offset,
                                        src_az, src_el)
            pt_mjd[pt] = start_mjd + pt * ha_step * 0.997269 / 24.0
            tel_az, tel_el = get_corrected_azel(tel_azel['tel_az'],
                                                tel_azel['tel_el'], False)
            pt_az[pt] = tel_az
            pt_el[pt] = tel_el
        tol_flag = True
        for pt in range(1, num_pts + 1, 2) :
            az_err = pt_az[pt] - (pt_az[pt - 1] + pt_az[pt + 1]) / 2.0
            az_err *= np.cos(pt_el[pt] * np.pi / 180.0)
            el_err = pt_el[pt] - (pt_el[pt - 1] + pt_el[pt + 1]) / 2.0
            err = np.sqrt(az_err**2 + el_err**2)
            if err > posn_tol :
                tol_flag = False
        if tol_flag :
            break
    if not tol_flag :
        print 'WARNING: make_track() did not converge!'
    num_pts_out = num_pts / 2
    step_sec = 2.0 * ha_step * 0.997269 * 3600.0  # seconds
    for pt in range(num_pts_out) :
        pt_mjd[pt] = pt_mjd[pt * 2]
        pt_az[pt] = pt_az[pt * 2]
        pt_el[pt] = pt_el[pt * 2]
        pt_az_vel[pt] = (pt_az[pt * 2 + 2] - pt_az[pt * 2]) / step_sec
        pt_el_vel[pt] = (pt_el[pt * 2 + 2] - pt_el[pt * 2]) / step_sec
    return {'mjd':pt_mjd[:num_pts_out],
            'az':pt_az[:num_pts_out], 'az_vel':pt_az_vel[:num_pts_out],
            'el':pt_el[:num_pts_out], 'el_vel':pt_el_vel[:num_pts_out]}

#print make_track('19:22:33,65:33:77', 54642.7803896, 54642.7803896 + 30.0/86400, 5.0, 5.0)
#print make_track('GOES 12', 54642.7803896, 54642.7803896 + 30.0/86400, 5.0, 5.0)
#begin_mjd = get_current_mjd()
#for i in range(100) :
#    start_mjd = get_current_mjd() + 0.5 / 86400.0
#    stop_mjd = start_mjd + 10.0 / 86400.0
#    trk = make_track('08:50:00.0,+11:00:007.0', start_mjd, stop_mjd, 5.0, 5.0)
#end_mjd = get_current_mjd()
#print '100 make_track()s:', (end_mjd - begin_mjd) * 86400.0

def wait_for_slew ( cmd_az, cmd_el, tol_deg=0.02 ) :
    """
    wait_for_slew(cmd_az, cmd_el, tol_deg=0.02)
    Returns only when the telescope encoder positions match the parameter
    values within the specified tolerance.  All parameters are in degrees.
    Positions are checked once per second.
    """
    old_tel_stat = False
    old_mjd = 0.0
    tm.sleep(1.0)
    while True :
        new_mjd = get_current_mjd()
        #print '1111'
        tel_stat = ant.get_time_tagged_status()
        #print tel_stat['az'], cmd_az, tel_stat['el'], cmd_el
        if (abs(tel_stat['az'] - cmd_az) < tol_deg and
            abs(tel_stat['el'] - cmd_el) < tol_deg) :
            break
        if old_tel_stat :
            delta_t = (new_mjd - old_mjd) * 86400.0
            az_vel = (tel_stat['az'] - old_tel_stat['az']) / delta_t
            el_vel = (tel_stat['el'] - old_tel_stat['el']) / delta_t
            #print 'az_vel: %5.3f, el_vel: %5.3f deg/sec' % (az_vel, el_vel)
        old_mjd = new_mjd
        old_tel_stat = tel_stat
        log_ant_posn()
        tm.sleep(1.0)

def execute_track ( trk, stop_mjd, cross_el_offset, el_offset, dflag=False ) :
    """
    execute_track(trk, stop_mjd)
    Commands the 20-meter antenna to follow the track specified in the
    first parameter until the stop time specified in the second parameter.
    See the make_track() __doc__ for a description of the first parameter
    format.  The stop_mjd is in fractional Modified Julian Date (MJD +
    UTC fraction of a day).  The return value is Boolean True unless an
    error occurred.  This function is normally called by the track()
    function, which has called the make_track() function to get the first
    parameter.  The fifth parameter determines whether a start data
    message is sent to the data computer.
    """
    for i in range(len(trk['mjd'])) :
        while (trk['mjd'][i] - get_current_mjd()) * 86400.0 > 3.0 :
            tm.sleep(0.5)
            log_ant_posn()

        utc_sec = (trk['mjd'][i] - int(trk['mjd'][i])) * 86400.0
        #time_str = mjd_to_time_str(trk['mjd'][i])
        status = ant.set_azel(trk['az'][i], trk['el'][i], tim=utc_sec,
                              azvel=trk['az_vel'][i], elvel=trk['el_vel'][i])
        if not status :
            return status
        #print get_current_mjd(), '+++++++++'
        if i == 0 and dflag :
            tm.sleep(1.0)
            msg1 = 'Az=%5.3f,El=%5.3f,' % (trk['az'][i], trk['el'][i])
            msg2 = 'Xeloff=%5.3f,Eloff=%5.3f' % (cross_el_offset, el_offset)
            _skt.send('Arrived at ' + msg1 + msg2 + '\r\n')
    ctr = 0
    while stop_mjd > get_current_mjd() :
        tm.sleep(0.5)
        ctr += 1
        if ctr % 2 :
            log_ant_posn()
    if dflag :
        data = _skt.recv(1024)
        print 'Received:', repr(data)
    return True

def track ( object, duration, cross_el_offset=0.0, el_offset=0.0,
            dflag=False ) :
    """
    track(object, duration, cross_el_offset=0.0, el_offset=0.0, dflag=False)
    This is the primary user function for sending the 20-meter telescope to
    a radio source or J2000 R.A./Dec position on the sky and tracking it for
    the number of seconds given in the second parameter.  The first parameter
    can be one of three name strings, 'CasA', 'CygA', or 'GOES 12', or a
    J2000 coordinate pair as a hexagesimal string of the following format:
    'HH:MM:SS.s,-DD:MM:SS.s'.  The tracking duration begins after the
    telescope has reached the commanded position.  The pointing may be offset
    in near-rectangular coordinates from the commanded position by the last
    two parameters, which are in degrees.  Positive cross_el_offset drives
    the telescope to smaller azimuth, and positive el_offset drives the
    telescope to lower elevation.  The track is adjusted to keep the
    commanded R.A./Dec in a fixed point of the antenna beam pattern for the
    purpose of beam mapping.  The fifth parameter determines whether a
    start data message is sent to the data computer.
    """
    tel_stat = ant.get_time_tagged_status()
    cur_az = tel_stat['az']
    cur_el = tel_stat['el']
    #mjd = get_current_mjd()  # includes day fraction
    #ra, dec = get_object_posn(object)  # precessed & nutated
    #last = get_last(mjd)  # in fraction of a day
    #obj_ha = 24.0 * last - ra
    #if obj_ha < -12.0 :
    #    obj_ha += 24.0
    #if obj_ha > 12.0 :
    #    obj_ha -= 24.0
    #az, el = hadec_to_azel(obj_ha, dec)
    #az_dist = abs(az - cur_az)

    start_mjd = get_current_mjd() + 0.5 / 86400.0
    log_ant_command(object, start_mjd, 0.0, cross_el_offset, el_offset)
    stop_mjd = start_mjd + 5.0 / 86400.0
    trk = make_track(object, start_mjd, stop_mjd,
                     cross_el_offset, el_offset)
    az_dist = abs(trk['az'][0] - cur_az)
    if az_dist > 180 :
        az_dist = abs(360 - az_dist)
    el_dist = abs(trk['el'][0] - cur_el)
    slew_time = max(az_dist / 2.0, el_dist / 2.0)  # seconds
    print 'slew time:', slew_time, 'sec to', trk['az'][0], trk['el'][0]

    start_mjd = get_current_mjd() + slew_time / 86400.0
    stop_mjd = start_mjd + 5.0 / 86400.0
    trk = make_track(object, start_mjd, stop_mjd,
                     cross_el_offset, el_offset)

    status = ant.set_azel(trk['az'][0], trk['el'][0])
    if not status :
        return status

    wait_for_slew(trk['az'][0], trk['el'][0])

    start_mjd = get_current_mjd() + 0.5 / 86400.0
    stop_mjd = start_mjd + duration / 86400.0
    trk = make_track(object, start_mjd, stop_mjd,
                     cross_el_offset, el_offset)
    #print trk
    log_ant_command(object, start_mjd, stop_mjd, cross_el_offset, el_offset)
    status = execute_track(trk, stop_mjd, cross_el_offset, el_offset, dflag)
    ant.stop()
    return status

_st = tm.gmtime()
_fn = 'log%4d_%02d_%02d_%02d:%02d:%02d' % (_st[0], _st[1], _st[2],
                                           _st[3], _st[4], _st[5])
_log_fd = file(_fn, 'wb')
#print 'No Log File'

def log_ant_posn ( ) :
    """
    log_ant_posn() sends the current fractional MJD, antenna controller time,
    and antenna position to a text log file, which was opened when the
    'array_pointing' module was imported into python.  The file name contains
    the UTC date and time when it was created, e.g., log2008_07_01_15:56:00.
    """
    mjd = get_current_mjd()
    tel_stat = ant.get_time_tagged_status()
    _log_fd.write('%02d:%02d:%02d %15.9f %8.4f %8.4f\n' % \
                      (tel_stat['UTC_hr'], tel_stat['UTC_min'],
                       tel_stat['UTC_sec'], mjd, tel_stat['az'],
                       tel_stat['el']))
    _log_fd.flush()

def log_ant_command ( object, start_mjd, stop_mjd,
                      cross_el_offset, el_offset ) :
    """
    log_ant_command(object, start_mjd, stop_mjd, cross_el_offset, el_offset)
    Sends the given parameters to a text log file, which was opened when the
    'array_pointing' module was imported into python.  The file name contains
    the UTC date and time when it was created, e.g., log2008_07_01_15:56:00.
    To distinguish these command parameters from the actual antenna times
    and positions in the same file, these commands are preceded by the '#'
    character
    """
    _log_fd.write('# %s %15.9f %15.9f %6.4f %6.4f\n' %
                  (object, start_mjd, stop_mjd, cross_el_offset, el_offset))
    _log_fd.flush()

def mjd_to_time_str ( mjd ) :
    """
    mjd_to_time_str(mjd) converts fractional MJD to 'HH:MM:SS.ss' format.
    """
    utc_sec = (mjd - int(mjd)) * 86400.0
    print 'mjd sec:', utc_sec
    hrs = int(utc_sec / 3600.0)
    utc_sec -= 3600.0 * hrs
    mint = int(utc_sec / 60.0)
    utc_sec -= 60.0 * mint
    return '%02d:%02d:%05.2f' % (hrs, mint, utc_sec)

#print mjd_to_time_str(get_current_mjd())

def box ( object, increment, box_num, duration, xoff=0.0, eloff=0.0,
          dflag=False ) :
    """
    box(object, increment, box_num, duration, dflag=False)
    Executes a set of pointings in a square box around the given source
    position.  The first input parameter can be one of three name strings,
    'CasA', 'CygA', or 'GOES 12' or a J2000 coordinate pair as hexagesimal
    strings in the following format: 'HH:MM:SS.s,-DD:MM:SS.s'.  The second
    parameter is the spacing between pointings in decimal degrees.  The third
    parameter is the size of the box in number of increments.  box_num=0 will
    do one pointing on the object.  box_num=1 will execute the offset sequence,
    (-1,1), (0,1), (1,1), (1,0), (1,-1), (0,-1), (-1,-1), (-1,0).
    box_num=2 will execute the sequence (-2,2), -1,2), (0,2), etc.  The
    third parameter is the dwell time on each point in seconds.  The fifth
    parameter determines whether a start data message is sent to the data
    computer.
    """
    if box_num == 0 :
        track(object, duration, xoff, eloff, dflag)
        #print '0,0'
        return

    # top
    for x in range(-box_num, box_num+1) :
        track(object, duration, x * increment + xoff,
              box_num * increment + eloff, dflag)
        #print x, box_num
    # + side
    for y in range(box_num-1, -box_num, -1) :
        track(object, duration, box_num * increment + xoff,
              y * increment + eloff, dflag)
        #print box_num, y
    # bottom
    for x in range(box_num, -box_num-1, -1) :
        track(object, duration, x * increment + xoff,
              -box_num * increment + eloff, dflag)
        #print x, -box_num
    # - side
    for y in range(-box_num+1, box_num) :
        track(object, duration, -box_num * increment + xoff,
              y * increment + eloff, dflag)
        #print -box_num, y

def box_set ( object, increment, num_boxes, duration, xoff=0.0, eloff=0.0,
              repeat_on=True, dflag=False ) :
    """
    box_set(object, increment, num_boxes, duration, dflag=False)
    Executes a set of concentric boxes around the object specified in the
    first parameter.  See the box.__doc__ for the definition of a box and
    the second and fourth parameters.  If num_boxes=0, only one pointing
    on the object will be done.  If num_boxes=1, the on-source pointing and
    the smallest box will be done.  If the last parameter is True (default),
    an on-source pointing will be done after each box.  The fifth
    parameter determines whether a start data message is sent to the data
    computer.
    """
    box(object, increment, 0, duration, xoff, eloff, dflag)

    for box_num in range(1, num_boxes + 1) :
        box(object, increment, box_num, duration, xoff, eloff, dflag)
        if repeat_on :
            box(object, increment, 0, duration, xoff, eloff, dflag)

def raster ( object, increment, max_x, max_y, duration, xoff=0.0, eloff=0.0,
             repeat_on=True, dflag=False) :
    """
    raster(object, increment, max_x, max_y, duration, repeat_on=True)
    Executes a symmetric raster scan around the given source position.
    The first input parameter can be one of three name strings, 'CasA', 'CygA',
    or 'GOES 12', or a J2000 coordinate pair as hexagesimal strings in the
    following format: 'HH:MM:SS.s,-DD:MM:SS.s'.  The second parameter is
    the spacing between pointings in decimal degrees. The third and fourth
    parameters are the number of increments from the center position to the
    sides of the raster, e.g. max_x=2, max_y=1 would do 3 rows of 5
    pointings.  The fifth parameter is the dwell time on each pointing.
    If the last parameter is True (default), an on-source pointing will be
    done after each row.  The center row is done first, then one increment
    above and below, two increments above and below, etc.  The seventh
    parameter determines whether a start data message is sent to the data
    computer.
    """
    for y in range(max_y + 1) :
        for x in range(-max_x, max_x + 1) :
            track(object, duration,
                  x * increment + xoff, y * increment + eloff, dflag)
        if repeat_on :
            track(object, duration, xoff, eloff, dflag)
        if y != 0 :
            for x in range(-max_x, max_x + 1) :
                track(object, duration,
                      x * increment + xoff, -y * increment + eloff, dflag)
            if repeat_on :
                track(object, duration, xoff, eloff, dflag)

def column ( object, increment, max_y, duration, xoff=0.0, eloff=0.0,
             dflag=False) :
    """
    column(object, increment, max_y, duration, dflag=False)
    Executes a symmetric raster scan around the given source position.
    The first input parameter can be one of three name strings, 'CasA', 'CygA',
    or 'GOES 12', or a J2000 coordinate pair as hexagesimal strings in the
    following format: 'HH:MM:SS.s,-DD:MM:SS.s'.  The second parameter is
    the spacing between pointings in decimal degrees. The third and fourth
    parameters are the number of increments from the center position to the
    sides of the raster, e.g. max_x=2, max_y=1 would do 3 rows of 5
    pointings.  The fifth parameter is the dwell time on each pointing.
    If the last parameter is True (default), an on-source pointing will be
    done after each row.  The center row is done first, then one increment
    above and below, two increments above and below, etc.  The seventh
    parameter determines whether a start data message is sent to the data
    computer.
    """
    for y in range(-max_y, max_y + 1) :
        track(object, duration, xoff, y * increment + eloff, dflag)
    
def five_point ( object, increment, duration, xoff=0.0, eloff=0.0,
                 dflag=False ) :
    """
    five_point(object, increment, duration, dflag=False)
    Executes a symmetric five-point scan around the given source position.
    The first input parameter can be one of three name strings, 'CasA', 'CygA',
    or 'GOES 12', or a J2000 coordinate pair as hexagesimal strings in the
    following format: 'HH:MM:SS.s,-DD:MM:SS.s'.  The second parameter is
    the spacing between pointings in decimal degrees.  The on-source scan is
    done first, then -az, -el, +az, +el.  The fourth
    parameter determines whether a start data message is sent to the data
    computer.
    """
    track(object, duration, xoff, eloff, dflag)
    track(object, duration, increment + xoff, eloff, dflag)
    track(object, duration, xoff, increment + eloff, dflag)
    track(object, duration, -increment + xoff, eloff, dflag)
    track(object, duration, xoff, -increment + eloff, dflag)

def five_pointx ( object, increment, duration, xoff=0.0, eloff=0.0,
                  dflag=False ) :
    """
    five_point(object, increment, duration, dflag=False)
    Executes a symmetric five-point scan around the given source position.
    The first input parameter can be one of three name strings, 'CasA', 'CygA',
    or 'GOES 12', or a J2000 coordinate pair as hexagesimal strings in the
    following format: 'HH:MM:SS.s,-DD:MM:SS.s'.  The second parameter is
    the spacing between pointings in decimal degrees.  The on-source scan is
    done first, then -az, -el, +az, +el.  The fourth
    parameter determines whether a start data message is sent to the data
    computer.
    """
    track(object, duration, 0.0, 0.0, dflag)
    track(object, duration, increment + xoff, 0.0, dflag)
    track(object, duration, 0.0, increment + eloff, dflag)
    track(object, duration, -increment + xoff, 0.0, dflag)
    track(object, duration, 0.0, -increment + eloff, dflag)

def test_one_cmd ( az, el, az_vel=0.0, el_vel=0.0 ) :
    ant.set_computer_timeout(20.0)
    mjd = get_current_mjd() + 3.0 / 86400.0
    time_str = mjd_to_time_str(mjd)
    #print time_str
    ant.set_azel(az, el, tim=time_str, azvel=az_vel, elvel=el_vel)
    while (mjd + 13.0 / 86400.0) > get_current_mjd() :
        log_ant_posn()
        tm.sleep(0.5)
    ant.stop()

def test_two_cmds ( az, el, az_vel=0.0, el_vel=0.0 ) :
    ant.set_computer_timeout(30.0)
    mjd = get_current_mjd() + 3.0 / 86400.0
    utc_sec = (mjd - int(mjd)) * 86400.0
    #time_str = mjd_to_time_str(mjd)
    ant.set_azel(az, el, tim=utc_sec, azvel=az_vel, elvel=el_vel)

    mjd += 10.0 / 86400.0
    #while mjd > get_current_mjd() :
    #    tm.sleep(1.0)
    #time_str = mjd_to_time_str(mjd)
    utc_sec = (mjd - int(mjd)) * 86400.0
    ant.set_azel(az + 1.0, el + 1.0, tim=utc_sec, azvel=az_vel, elvel=el_vel)
    while (mjd + 13.0 / 86400.0) > get_current_mjd() :
        log_ant_posn()
        tm.sleep(0.5)
    ant.stop()

def test_three_cmds ( az, el, az_vel=0.0, el_vel=0.0 ) :
    ant.set_computer_timeout(30.0)
    mjd = get_current_mjd() + 3.0 / 86400.0
    time_str = mjd_to_time_str(mjd)
    ant.set_azel(az, el, tim=time_str, azvel=az_vel, elvel=el_vel)

    mjd += 10.0 / 86400.0
    #while mjd > get_current_mjd() :
    #    tm.sleep(1.0)
    time_str = mjd_to_time_str(mjd)
    ant.set_azel(az + 1.0, el + 1.0, tim=time_str, azvel=az_vel, elvel=el_vel)

    mjd += 10.0 / 86400.0
    #while mjd > get_current_mjd() :
    #    tm.sleep(1.0)
    time_str = mjd_to_time_str(mjd)
    ant.set_azel(az + 2.0, el + 2.0, tim=time_str, azvel=az_vel, elvel=el_vel)

    while (mjd + 13.0 / 86400.0) > get_current_mjd() :
        log_ant_posn()
        tm.sleep(0.5)
    ant.stop()

def test_slew ( az, el ) :
    ant.set_computer_timeout(30.0)
    ant.set_azel(az, el)
    wait_for_slew(az, el, tol_deg=0.02)
    ant.stop()

def test_execute_track ( az, el ) :
    num_pts = 3
    tmjd = np.zeros(num_pts, dtype=np.float64)
    taz = np.zeros(num_pts, dtype=np.float64)
    taz_vel = np.zeros(num_pts, dtype=np.float64)
    tel = np.zeros(num_pts, dtype=np.float64)
    tel_vel = np.zeros(num_pts, dtype=np.float64)
    delta_t = 5.0
    delta_p = 0.1
    mjd = get_current_mjd() + 3.0 / 86400.0
    for i in range(num_pts) :
        tmjd[i] = mjd + i * delta_t
        taz[i] = az + i * delta_p
        taz_vel[i] = delta_p / delta_t
        tel[i] = el + i * delta_p
        tel_vel[i] = delta_p / delta_t
    stop_mjd = tmjd[-1] + delta_t
    trk = {'mjd':tmjd, 'az':taz, 'az_vel':taz_vel, 'el':tel, 'el_vel':tel_vel}
    execute_track(trk, stop_mjd)