#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]);

    if (!isfinite(data->zp)) {
        psErrorStackPrint(stderr, "Zero point is unknown\n");
        exit(PS_EXIT_CONFIG_ERROR);
    }

    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 numIn = detections->n;        // Number of rows in input
    int numOut = 0;                   // Number of rows in output
    psArray *output = psArrayAllocEmpty(numIn); // Output table
    double plateScale = 0.0;                 // Average plate scale
    for (int i = 0; i < numIn; i++) {
        psMetadata *inRow = detections->data[i]; // Input 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");
        float ext = psMetadataLookupF32(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");

        if (!isfinite(mag) || !isfinite(magErr) || (flags & SOURCE_MASK) ||
            ((flags & PM_SOURCE_MODE_DEFECT) && !(flags & PM_SOURCE_MODE_MOMENTS_FAILURE))
            ) {
            // DEFECT can be due to bad moments
            continue;
        }

        psMetadata *outRow = output->data[numOut] = psMetadataAlloc(); // Output row

        numOut++;
        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 cosAngle2 = PS_SQR(cosAngle), sinAngle2 = PS_SQR(sinAngle);
        double xErr2 = PS_SQR(xErr), yErr2 = PS_SQR(yErr);
        double errScale = scale / 3600.0;
        double raErr = errScale * sqrt(cosAngle2 * xErr2 + sinAngle2 * yErr2);
        double decErr = errScale * sqrt(sinAngle2 * xErr2 + cosAngle2 * yErr2);

        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);
    }
    output->n = numOut;
    plateScale /= numOut;
    psFree(detections);

    // Translate the header
    psMetadata *outHeader = psMetadataAlloc(); // Output header
    ppMopsVersionHeader(outHeader);
    {
        const char *ra = psMetadataLookupStr(NULL, header, "FPA.RA");
        const char *dec = psMetadataLookupStr(NULL, header, "FPA.DEC");
        const char *filter = psMetadataLookupStr(NULL, header, "FPA.FILTER");
        float airmass = psMetadataLookupF32(NULL, header, "FPA.AIRMASS");
        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");
        double mjd = psMetadataLookupF64(NULL, header, "MJD-OBS") + exptime / 2.0 / 3600 / 24;

        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");


        psMetadataAddStr(outHeader, PS_LIST_TAIL, "RA", 0, "Right ascension of boresight", ra);
        psMetadataAddStr(outHeader, PS_LIST_TAIL, "DEC", 0, "Declination of boresight", dec);
        psMetadataAddF32(outHeader, PS_LIST_TAIL, "AIRMASS", 0, "Airmass of exposure", airmass);
        psMetadataAddF64(outHeader, PS_LIST_TAIL, "MJD-OBS", 0, "MJD of exposure midpoint", mjd);
        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);
        psMetadataAddS32(outHeader, PS_LIST_TAIL, "DIFFIMID", 0, "Difference image identifier", imageid);
        psMetadataAddS32(outHeader, PS_LIST_TAIL, "FPA_ID", 0, "Exposure identifier", fpaid);
        psMetadataAddStr(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 (FAKE)", 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 (numOut == 0) {
            // Write dummy table
            psMetadata *outRow = psMetadataAlloc(); // Dummy output row
            psMetadataAddF64(outRow, PS_LIST_TAIL, "RA_DEG", 0, "Right ascension (degrees)", NAN);
            psMetadataAddF64(outRow, PS_LIST_TAIL, "DEC_DEG", 0, "Declination (degrees)", NAN);
            psMetadataAddF64(outRow, PS_LIST_TAIL, "RA_SIG", 0, "Right ascension error (degrees)", NAN);
            psMetadataAddF64(outRow, PS_LIST_TAIL, "DEC_SIG", 0, "Declination error (degrees)", NAN);
            psMetadataAddF64(outRow, PS_LIST_TAIL, "MAG", 0, "Magnitude", NAN);
            psMetadataAddF64(outRow, PS_LIST_TAIL, "MAG_SIG", 0, "Magnitude error", NAN);
            psMetadataAddU32(outRow, PS_LIST_TAIL, "FLAGS", 0, "Detection bit flags", 0);

            // The below need fixing
            psMetadataAddF64(outRow, PS_LIST_TAIL, "STARPSF", 0, "EXT_NSIGMA", NAN);
            psMetadataAddF64(outRow, PS_LIST_TAIL, "ANG", 0, "Position angle of trail (degrees)", NAN);
            psMetadataAddF64(outRow, PS_LIST_TAIL, "ANG_SIG", 0, "Position angle error (degrees)", NAN);
            psMetadataAddF64(outRow, PS_LIST_TAIL, "LEN", 0, "Length of trail (arcsec)", NAN);
            psMetadataAddF64(outRow, PS_LIST_TAIL, "LEN_SIG", 0, "Length error (arcsec)", NAN);
            if (!psFitsWriteTableEmpty(fits, outHeader, outRow, OUT_EXTNAME)) {
                psErrorStackPrint(stderr, "Unable to write table.");
                psFree(outHeader);
                psFree(output);
                psFree(outRow);
                psFree(data);
                exit(PS_EXIT_SYS_ERROR);
            }
            psFree(outRow);
        } else 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;
}
