import types
from taskutil import get_global_namespace

def readantenna(antab=None, pl=None):
    """
    Read a text antenna table antab with columns X, Y, Z, and Diameter (all in m).
    """
    if not isinstance(pl, types.ModuleType):
        myf = get_global_namespace()
        pl = myf['pl']
    
    f = open(antab)
    stnx = []
    stny = []
    stnz = []
    stnd = []
    nant = 0
    line = '    '
    while (len(line) > 0):
        try: 
            line = f.readline()
            if (line.find('#') != 0): # If no # at start...
                ###ignoring line that has less than 4 elements
                if(len(line.split()) >3):
                    splitline = line.split()
                    stnx.append(float(splitline[0]))
                    stny.append(float(splitline[1]))
                    stnz.append(float(splitline[2]))
                    stnd.append(float(splitline[3]))
                    nant += 1					
        except:
            break

    f.close()
    return (pl.array(stnx), pl.array(stny), pl.array(stnz), pl.array(stnd), nant)

def mul_scal_qa(scal, q, qa=None):
    """
    Returns the product of a dimensionless number scal and quantity q.
    """
    if not qa:
        myf = get_global_namespace()
        qa = myf['qa']
    
    return qa.quantity(scal * q['value'], q['unit'])

def calc_pointings(direction, spacing, imsize, cell, relmargin=0.5):
    """
    If direction is a list, simply returns direction and the number of
    pointings in it.

    Otherwise, returns a hexagonally packed list of pointings separated by
    spacing and fitting inside an image specified by direction, imsize and
    cell, and the number of pointings.  The hexagonal packing starts with a
    horizontal row centered on direction, and the other rows alternate
    being horizontally offset by a half spacing.  All of the pointings will
    be within a rectangle relmargin * spacing smaller than the image on all
    sides, unless none can fit.  In that case direction is returned
    (i.e. single pointing).
    """
    myf = get_global_namespace()
    pl = myf['pl']
    qa = myf['qa']
    
    if type(direction) == list:
        return len(direction), direction

    epoch, centx, centy = direction_splitter(direction, qa)

    spacing  = qa.quantity(spacing)
    yspacing = mul_scal_qa(0.866025404, spacing, qa)

    if type(cell) == list:
        cellx, celly = map(qa.quantity, cell)
    else:
        cellx = qa.quantity(cell)
        celly = cellx

    ysize = mul_scal_qa(imsize[1], celly)
    nrows = 1+ int(pl.floor(qa.convert(qa.div(ysize, yspacing), '')['value']
                            - 2.309401077 * relmargin))

    xsize = mul_scal_qa(imsize[0], cellx, qa)
    availcols = 1 + qa.convert(qa.div(xsize, spacing),
                               '')['value'] - 2.0 * relmargin
    ncols = int(pl.floor(availcols))

    # By making the even rows shifted spacing/2 ahead, and possibly shorter,
    # the top and bottom rows (nrows odd), are guaranteed to be short.
    if availcols - ncols >= 0.5:                            # O O O
        evencols = ncols                                    #  O O O
        ncolstomin = 0.5 * (ncols - 0.5)
    else:
        evencols = ncols - 1                                #  O O 
        ncolstomin = 0.5 * (ncols - 1)                      # O O O
    pointings = []

    # Start from the top because in the Southern hemisphere it sets first.
    y = qa.add(centy, mul_scal_qa(0.5 * (nrows - 1), yspacing, qa))
    for row in xrange(0, nrows):         # xrange stops early.
        xspacing = mul_scal_qa(1.0 / pl.cos(qa.convert(y, 'rad')['value']),
                               spacing, qa)
        ystr = qa.formxxx(y, format='dms')

        if row % 2:                             # Odd
            xmin = qa.sub(centx, mul_scal_qa(ncolstomin, xspacing, qa))
            stopcolp1 = ncols
        else:                                   # Even (including 0)
            xmin = qa.sub(centx, mul_scal_qa(ncolstomin - 0.5, xspacing, qa))
            stopcolp1 = evencols
        for col in xrange(0, stopcolp1):        # xrange stops early.
            x = qa.formxxx(qa.add(xmin, mul_scal_qa(col, xspacing, qa)),
                           format='hms')
            pointings.append("%s%s %s" % (epoch, x, ystr))
        y = qa.sub(y, yspacing)

    if len(pointings) < 1:
        print "No pointings would fit within relmargin =", relmargin
        print "Using single pointing at", direction
        pointings = [direction]
    else:
        print "Using generated pointings:"
        print pointings

    return len(pointings), pointings

def plot_pointings(pointings, img, beamsize, dolab = True):
    pass

def average_direction(directions):
    """
    Returns the average of directions as a string.
    """
    myf = get_global_namespace()
    pl = myf['pl']
    qa = myf['qa']

    epoch0, x, y = direction_splitter(directions[0], qa)
    i = 1
    avgx = 0.0
    avgy = 0.0
    for drn in directions:
        epoch, x, y = direction_splitter(drn, qa)
        x = x['value']
        y = y['value']
        if epoch != epoch0:                     # Paranoia
            print "Warning: precession not handled by average_direction()"
        x = wrapang(x, avgx, 360.0, pl)
        avgx += (x - avgx) / i
        avgy += (y - avgy) / i
        i += 1
    avgx = qa.toangle('%fdeg' % avgx)
    avgy = qa.toangle('%fdeg' % avgy)
    avgx = qa.formxxx(avgx, format='hms')
    avgy = qa.formxxx(avgy, format='dms')
    return "%s%s %s" % (epoch0, avgx, avgy)

def direction_splitter(direction, qa=None):
    """
    Given a direction, return its epoch, x, and y parts.  Epoch will be ''
    if absent, or '%s ' % epoch if present.  x and y will be angle qa's in
    degrees.
    """
    dirl = direction.split()
    if len(dirl) == 3:
        epoch = dirl[0] + ' '
    else:
        epoch = ''
    if not qa:
        myf = get_global_namespace()
        qa = myf['qa']
    x, y = map(qa.toangle, dirl[-2:])
    return epoch, qa.convert(x, 'deg'), qa.convert(y, 'deg')

def wrapang(ang, target, period = 360.0, pl=None):
    """
    Returns ang wrapped so that it is within +-period/2 of target.
    """
    if not isinstance(pl, types.ModuleType):
        myf = get_global_namespace()
        pl = myf['pl']
    
    dang       = ang - target
    period     = pl.absolute(period)
    halfperiod = 0.5 * period
    if pl.absolute(dang) > halfperiod:
        nwraps = pl.floor(0.5 + float(dang) / period)
        ang -= nwraps * period
    return ang