def rd2xy(input,ra=None,dec=None,coordfile=None,colnames=None,
            precision=6,output=None,verbose=True):
    """ Primary interface to perform coordinate transformations from
        pixel to sky coordinates using STWCS and full distortion models
        read from the input image header.
    """
    single_coord = False
    if coordfile is not None:
        if colnames in blank_list:
            colnames = ['c1','c2']
        elif isinstance(colnames,type('a')):
            colnames = colnames.split(',')
        # convert input file coordinates to lists of decimal degrees values
        xlist,ylist = tweakutils.readcols(coordfile,cols=colnames)
    else:
        if isinstance(ra,np.ndarray):
            ralist = ra.tolist()
            declist = dec.tolist()
        elif not isinstance(ra, list):
            ralist = [ra]
            declist = [dec]
        else:
            ralist = ra
            declist = dec
        xlist  = [0]*len(ralist)
        ylist = [0]*len(ralist)
        if len(xlist) == 1:
            single_coord = True
        for i,(r,d) in enumerate(zip(ralist,declist)):
            # convert input value into decimal degrees value
            xval,yval = tweakutils.parse_skypos(r,d)
            xlist[i] = xval
            ylist[i] = yval

    # start by reading in WCS+distortion info for input image
    inwcs = wcsutil.HSTWCS(input)
    if inwcs.wcs.is_unity():
        print("####\nNo valid WCS found in {}.\n  Results may be invalid.\n####\n".format(input))

    # Now, convert pixel coordinates into sky coordinates
    try:
        outx,outy = inwcs.all_world2pix(xlist,ylist,1)
    except RuntimeError:
        outx,outy = inwcs.wcs_world2pix(xlist,ylist,1)

    # add formatting based on precision here...
    xstr = []
    ystr = []
    fmt = "%."+repr(precision)+"f"
    for x,y in zip(outx,outy):
        xstr.append(fmt%x)
        ystr.append(fmt%y)

    if verbose or (not verbose and util.is_blank(output)):
        print ('# Coordinate transformations for ',input)
        print('# X      Y         RA             Dec\n')
        for x,y,r,d in zip(xstr,ystr,xlist,ylist):
            print("%s  %s    %s  %s"%(x,y,r,d))

    # Create output file, if specified
    if output:
        f = open(output,mode='w')
        f.write("# Coordinates converted from %s\n"%input)
        for x,y in zip(xstr,ystr):
            f.write('%s    %s\n'%(x,y))
        f.close()
        print('Wrote out results to: ',output)

    if single_coord:
        outx = outx[0]
        outy = outy[0]

    return outx, outy