#include <stdio.h>
#include <pslib.h>

#include "ppMops.h"

int main(int argc, char *argv[])
{
    if (argc != 4) {
        fprintf(stderr, "Insufficient arguments.\n");
        fprintf(stderr, "Usage: %s DETECTIONS ZP OUTPUT\n", argv[0]);
        exit(PS_EXIT_CONFIG_ERROR);
    }

    ppMopsData *data = ppMopsDataAlloc(); // Configuration data
    data->detections = psStringCopy(argv[1]);
    data->zp = atof(argv[2]);
    data->output = psStringCopy(argv[3]);

    psArray *detections = NULL;         // Detections
    psMetadata *header = NULL;          // Header for detections
    {
        psFits *fits = psFitsOpen(data->detections, "r"); // FITS file
        header = psFitsReadHeader(NULL, fits);
        if (!header) {
            psErrorStackPrint(stderr, "Unable to read header");
            psFitsClose(fits);
            psFree(data);
            exit(PS_EXIT_DATA_ERROR);
        }
        if (!psFitsMoveExtName(fits, IN_EXTNAME)) {
            psErrorStackPrint(stderr, "Unable to move to extension %s", IN_EXTNAME);
            psFitsClose(fits);
            psFree(header);
            psFree(data);
            exit(PS_EXIT_DATA_ERROR);
        }
        detections = psFitsReadTable(fits);
        psFitsClose(fits);
        if (!detections) {
            psErrorStackPrint(stderr, "Unable to read detections");
            psFree(data);
            psFree(header);
            exit(PS_EXIT_DATA_ERROR);
        }
    }

    // Translate the columns
    int numRows = detections->n;        // Number of rows
    psArray *output = psArrayAlloc(numRows); // Output table
    double plateScale = 0.0;                 // Average plate scale
    for (int i = 0; i < numRows; i++) {
        psMetadata *inRow = detections->data[i]; // Input row
        psMetadata *outRow = output->data[i] = psMetadataAlloc(); // Output row

        double ra = psMetadataLookupF64(NULL, inRow, "RA_PSF");
        double dec = psMetadataLookupF64(NULL, inRow, "DEC_PSF");
        double mag = psMetadataLookupF64(NULL, inRow, "PSF_INST_MAG") + data->zp;
        double magErr = psMetadataLookupF64(NULL, inRow, "PSF_INST_MAG_SIG");
        double ext = psMetadataLookupF64(NULL, inRow, "EXT_NSIGMA");
        double xErr = psMetadataLookupF64(NULL, inRow, "X_PSF_SIG");
        double yErr = psMetadataLookupF64(NULL, inRow, "Y_PSF_SIG");
        double scale = psMetadataLookupF64(NULL, inRow, "PLTSCALE");
        double angle = psMetadataLookupF64(NULL, inRow, "POSANGLE");
        psU32 flags = psMetadataLookupU32(NULL, inRow, "FLAGS");

        plateScale += scale;

        // XXX Not at all sure I've got the angles around the right way here...
        double cosAngle = cos(angle), sinAngle = sin(angle);
        double raErr = scale * sqrt(PS_SQR(cosAngle) * PS_SQR(xErr) + PS_SQR(sinAngle) * PS_SQR(yErr));
        double decErr = scale * sqrt(PS_SQR(sinAngle) * PS_SQR(xErr) + PS_SQR(cosAngle) * PS_SQR(yErr));

        psMetadataAddF64(outRow, PS_LIST_TAIL, "RA_DEG", 0, "Right ascension (degrees)", ra);
        psMetadataAddF64(outRow, PS_LIST_TAIL, "DEC_DEG", 0, "Declination (degrees)", dec);
        psMetadataAddF64(outRow, PS_LIST_TAIL, "RA_SIG", 0, "Right ascension error (degrees)", raErr);
        psMetadataAddF64(outRow, PS_LIST_TAIL, "DEC_SIG", 0, "Declination error (degrees)", decErr);
        psMetadataAddF64(outRow, PS_LIST_TAIL, "MAG", 0, "Magnitude", mag);
        psMetadataAddF64(outRow, PS_LIST_TAIL, "MAG_SIG", 0, "Magnitude error", magErr);
        psMetadataAddU32(outRow, PS_LIST_TAIL, "FLAGS", 0, "Detection bit flags", flags);

        // The below need fixing
        psMetadataAddF64(outRow, PS_LIST_TAIL, "STARPSF", 0, "EXT_NSIGMA", ext);
        psMetadataAddF64(outRow, PS_LIST_TAIL, "ANG", 0, "Position angle of trail (degrees)", 0.0);
        psMetadataAddF64(outRow, PS_LIST_TAIL, "ANG_SIG", 0, "Position angle error (degrees)", 0.0);
        psMetadataAddF64(outRow, PS_LIST_TAIL, "LEN", 0, "Length of trail (arcsec)", 0.0);
        psMetadataAddF64(outRow, PS_LIST_TAIL, "LEN_SIG", 0, "Length error (arcsec)", 0.0);
    }
    plateScale /= numRows;
    psFree(detections);

    // Translate the header
    psMetadata *outHeader = psMetadataAlloc(); // Output header
    ppMopsVersionHeader(outHeader);
    {
        double ra = psMetadataLookupF64(NULL, header, "FPA.RA");
        double dec = psMetadataLookupF64(NULL, header, "FPA.DEC");
        const char *filter = psMetadataLookupStr(NULL, header, "FPA.FILTER");
        float exptime = psMetadataLookupF32(NULL, header, "EXPTIME");
        double angle = psMetadataLookupF64(NULL, header, "FPA.POSANGLE");
        double alt = psMetadataLookupF64(NULL, header, "FPA.ALT");
        double az = psMetadataLookupF64(NULL, header, "FPA.AZ");
        int imageid = psMetadataLookupS32(NULL, header, "IMAGEID");

        float psf = plateScale * 0.5 * (psMetadataLookupF32(NULL, header, "FWHM_MAJ") +
                                        psMetadataLookupF32(NULL, header, "FWHM_MIN"));

        // XXX This is wrong
        int fpaid = psMetadataLookupS32(NULL, header, "IMAGEID");


        psMetadataAddF64(outHeader, PS_LIST_TAIL, "RA", 0, "Right ascension of boresight", ra);
        psMetadataAddF64(outHeader, PS_LIST_TAIL, "DEC", 0, "Declination of boresight", dec);
        psMetadataAddStr(outHeader, PS_LIST_TAIL, "FILTER", 0, "Filter name", filter);
        psMetadataAddF64(outHeader, PS_LIST_TAIL, "EXPTIME", 0, "Exposure time (sec)", exptime);
        psMetadataAddF64(outHeader, PS_LIST_TAIL, "ROTANGLE", 0, "Rotator position angle", angle);
        psMetadataAddF64(outHeader, PS_LIST_TAIL, "TEL_ALT", 0, "Telescope altitude", alt);
        psMetadataAddF64(outHeader, PS_LIST_TAIL, "TEL_AZ", 0, "Telescope azimuth", az);
        psMetadataAddF64(outHeader, PS_LIST_TAIL, "DIFFIMID", 0, "Difference image identifier", imageid);
        psMetadataAddF64(outHeader, PS_LIST_TAIL, "FPA_ID", 0, "Exposure identifier", fpaid);
        psMetadataAddS32(outHeader, PS_LIST_TAIL, "OBSCODE", 0, "IAU Observatory code", OBSERVATORY_CODE);
        psMetadataAddF32(outHeader, PS_LIST_TAIL, "STARPSF", 0, "Stellar PSF (arcsec)", psf);

        // These are completely fake
        psMetadataAddF32(outHeader, PS_LIST_TAIL, "LIMITMAG", 0, "Limiting magnitude", 25.0);
        psMetadataAddF32(outHeader, PS_LIST_TAIL, "DE1", 0, "Detection efficiency (FAKE)", 0.0);
        psMetadataAddF32(outHeader, PS_LIST_TAIL, "DE2", 0, "Detection efficiency (FAKE)", 0.0);
        psMetadataAddF32(outHeader, PS_LIST_TAIL, "DE3", 0, "Detection efficiency (FAKE)", 0.0);
        psMetadataAddF32(outHeader, PS_LIST_TAIL, "DE4", 0, "Detection efficiency (FAKE)", 0.0);
        psMetadataAddF32(outHeader, PS_LIST_TAIL, "DE5", 0, "Detection efficiency (FAKE)", 0.0);
        psMetadataAddF32(outHeader, PS_LIST_TAIL, "DE6", 0, "Detection efficiency (FAKE)", 0.0);
        psMetadataAddF32(outHeader, PS_LIST_TAIL, "DE7", 0, "Detection efficiency (FAKE)", 0.0);
        psMetadataAddF32(outHeader, PS_LIST_TAIL, "DE8", 0, "Detection efficiency (FAKE)", 0.0);
        psMetadataAddF32(outHeader, PS_LIST_TAIL, "DE9", 0, "Detection efficiency (FAKE)", 0.0);
        psMetadataAddF32(outHeader, PS_LIST_TAIL, "DE10", 0, "Detection efficiency (FAKE)", 0.0);
    }
    psFree(header);

    // Write the new table
    {
        psFits *fits = psFitsOpen(data->output, "w"); // FITS file
        if (!fits) {
            psErrorStackPrint(stderr, "Unable to open %s for writing", data->output);
            psFree(outHeader);
            psFree(output);
            psFree(data);
            exit(PS_EXIT_SYS_ERROR);
        }
        if (!psFitsWriteTable(fits, outHeader, output, OUT_EXTNAME)) {
            psErrorStackPrint(stderr, "Unable to write table.");
            psFree(outHeader);
            psFree(output);
            psFree(data);
            exit(PS_EXIT_SYS_ERROR);
        }
        psFitsClose(fits);
    }

    psFree(outHeader);
    psFree(output);
    psFree(data);

    return PS_EXIT_SUCCESS;
}
