#!/usr/bin/env python

import os, sys, numpy, pyfits, math
import astLib.astWCS as astWCS

if __name__ == "__main__":

    catfn = sys.argv[1]
    fitsfn = sys.argv[2]
    size = float(sys.argv[3])
    output_dir = sys.argv[4]
    suffix = sys.argv[5]

    with open(catfn, "r") as cat:
        lines = cat.readlines()

    hdulist = pyfits.open(fitsfn)
    imghdu = hdulist[0]
    wcs = astWCS.WCS(imghdu.header, mode='pyfits')
    pixelscale = wcs.getPixelSizeDeg() * 3600.
    if (size < 0):
        npixels = int(numpy.fabs(size))
    else:
        npixels = int(size / pixelscale)


    for line in lines[2:]:
        items = line.split(",")
        
        try:
            id = int(items[0])
            objname = items[1].replace(" ", "_")
            ra = float(items[2])
            dec = float(items[3])
            vel = float(items[5])

            dust = items[17].startswith("x")
            interacting = items[18].startswith("x")
            debris = items[19].startswith("x")
            shells = items[25].startswith("x")

        except:
            continue

            
        if (interacting or debris or shells):
            print "% 5d %-30s %.5f %.5f %d" % (id, objname, ra, dec, vel), \
                " %5s"*4 % (interacting, debris, shells, (interacting| debris| shells)), \
                "###"
        else:
            continue

        info = "%s" % (objname)
        if (dust): info += " - dust"
        if (interacting): info += " - interacting"
        if (debris): info += " - debris"
        if (shells): info += " - shells"

        xy = wcs.wcs2pix(ra, dec)

        x1 = numpy.max([0, int(math.floor(xy[0] - npixels))])
        x2 = numpy.min([imghdu.data.shape[1], int(math.ceil(xy[0] + npixels))])
        y1 = numpy.max([0, int(math.floor(xy[1] - npixels))])
        y2 = numpy.min([imghdu.data.shape[0], int(math.ceil(xy[1] + npixels))])

        cutout = imghdu.data[y1:y2, x1:x2]

        out_hdu = pyfits.PrimaryHDU(data=cutout, header=imghdu.header)
        out_hdu.header['CRPIX1'] -= x1
        out_hdu.header['CRPIX2'] -= y1
        out_hdu.header['OBJECT'] = info

        output_fn = "%s/%04d_%s.%s.fits" % (output_dir, id, objname, suffix)
        if (os.path.isfile(output_fn)):
            os.remove(output_fn)
        print "Writing cutout to %s" % (output_fn)
        out_hdu.writeto(output_fn, clobber=True)