#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

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

#include "ppMops.h"

psArray *ppMopsRead(const ppMopsArguments *args)
{
    psTrace("ppMops.read", 1, "Reading input detections\n");

    psArray *inNames = args->input;          // Input names
    long num = inNames->n;                   // Number of inputs
    psArray *detections = psArrayAlloc(num); // Array of detections, to return
    for (int i = 0; i < num; i++) {
        psFits *fits = psFitsOpen(inNames->data[i], "r"); // FITS file
        if (!fits) {
            psError(PS_ERR_IO, false, "Unable to open input %d", i);
            return false;
        }
        psMetadata *header = psFitsReadHeader(NULL, fits); // Primary header
        if (!header) {
            psError(PS_ERR_IO, false, "Unable to read header %d", i);
            return false;
        }

        psS64 diffSkyfileId = psMetadataLookupS64(NULL, header, "IMAGEID"); // Identifier for image
        if (diffSkyfileId == 0) {
            psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find identifier for image %d", i);
            return false;
        }

        if (!psFitsMoveExtName(fits, "SkyChip.psf")) {
            psError(PS_ERR_IO, false, "Unable to move to HDU with detections");
            return false;
        }

        long size = psFitsTableSize(fits); // Size of table
        if (size <= 0) {
            psErrorStackPrint(stderr, "Unable to determine size of table %d", i);
            psErrorClear();
            psWarning("Ignoring input %d", i);
            psFree(header);
            psFitsClose(fits);
            continue;
        }
        ppMopsDetections *det = ppMopsDetectionsAlloc(size);

        psTrace("ppMops.read", 3, "Reading %ld rows from %s\n", size, (const char*)inNames->data[i]);

        det->raBoresight = psMemIncrRefCounter(psMetadataLookupStr(NULL, header, "FPA.RA"));
        det->decBoresight = psMemIncrRefCounter(psMetadataLookupStr(NULL, header, "FPA.DEC"));
        det->filter = psMemIncrRefCounter(psMetadataLookupStr(NULL, header, "FPA.FILTER"));
        det->airmass = psMetadataLookupF32(NULL, header, "AIRMASS");
        det->exptime = psMetadataLookupF32(NULL, header, "EXPTIME");
        det->posangle = psMetadataLookupF64(NULL, header, "FPA.POSANGLE");
        det->alt = psMetadataLookupF64(NULL, header, "FPA.ALT");
        det->az = psMetadataLookupF64(NULL, header, "FPA.AZ");
        det->mjd = psMetadataLookupF64(NULL, header, "MJD-OBS") + det->exptime / 2.0 / 3600 / 24;

        det->seeing = 0.5 * (psMetadataLookupF32(NULL, header, "FWHM_MAJ") +
                             psMetadataLookupF32(NULL, header, "FWHM_MIN"));

        int naxis1 = psMetadataLookupS32(NULL, header, "IMNAXIS1"); // Number of columns
        int naxis2 = psMetadataLookupS32(NULL, header, "IMNAXIS2"); // Number of rows

        psFree(header);

        psArray *table = psFitsReadTable(fits); // Table of interest
        if (!table) {
            psError(PS_ERR_IO, false, "Unable to read table %d", i);
            return false;
        }
        psFitsClose(fits);

        double plateScale = 0.0;        // Plate scale
        long numGood = 0;               // Number of good rows
        for (long j = 0; j < size; j++) {
            psMetadata *row = table->data[j]; // Row of interest

            psU32 flags = psMetadataLookupU32(NULL, row, "FLAGS");
            if (flags & SOURCE_MASK) {
                psTrace("ppMops.read", 10, "Discarding row %ld from input %d because of flags: %ud", j, i, flags);
                continue;
            }

            det->x->data.F32[numGood] = psMetadataLookupF32(NULL, row, "X_PSF");
            det->y->data.F32[numGood] = psMetadataLookupF32(NULL, row, "Y_PSF");
            det->ra->data.F64[numGood] = DEG_TO_RAD(psMetadataLookupF64(NULL, row, "RA_PSF"));
            det->dec->data.F64[numGood] = DEG_TO_RAD(psMetadataLookupF64(NULL, row, "DEC_PSF"));
            det->mag->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_INST_MAG");
            det->magErr->data.F32[numGood] = psMetadataLookupF32(NULL, row, "PSF_INST_MAG_SIG");
            det->extended->data.F32[numGood] = psMetadataLookupF32(NULL, row, "EXT_NSIGMA");
            det->angle->data.F32[numGood] = 0.0;
            det->angleErr->data.F32[numGood] = 0.0;
            det->length->data.F32[numGood] = 0.0;
            det->lengthErr->data.F32[numGood] = 0.0;
            det->flags->data.U32[numGood] = psMetadataLookupU32(NULL, row, "FLAGS");
            det->diffSkyfileId->data.S64[numGood] = diffSkyfileId;
            det->naxis1->data.S32[numGood] = naxis1;
            det->naxis2->data.S32[numGood] = naxis2;

            // Calculate error in RA, Dec
            double xErr = psMetadataLookupF64(NULL, row, "X_PSF_SIG");
            double yErr = psMetadataLookupF64(NULL, row, "Y_PSF_SIG");
            double scale = psMetadataLookupF64(NULL, row, "PLTSCALE");
            double angle = psMetadataLookupF64(NULL, row, "POSANGLE");

            if (!isfinite(det->x->data.F32[numGood]) || !isfinite(det->y->data.F32[numGood]) ||
                !isfinite(det->ra->data.F64[numGood]) || !isfinite(det->dec->data.F64[numGood]) ||
                !isfinite(det->mag->data.F32[numGood]) || !isfinite(det->magErr->data.F32[numGood]) ||
                !isfinite(xErr) || !isfinite(yErr) || !isfinite(scale) || !isfinite(angle)) {
                psTrace("ppMops.read", 10,
                        "Discarding row %ld from input %d because of non-finite values: "
                        "%f %f %lf %lf %f %f %f %f %f %f",
                        j, i,
                        det->x->data.F32[numGood], det->y->data.F32[numGood],
                        det->ra->data.F64[numGood], det->dec->data.F64[numGood],
                        det->mag->data.F32[numGood], det->magErr->data.F32[numGood],
                        xErr, yErr, scale, angle);
                continue;
            }

            // 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;
            det->raErr->data.F64[numGood] = errScale * sqrt(cosAngle2 * xErr2 + sinAngle2 * yErr2);
            det->decErr->data.F64[numGood] = errScale * sqrt(sinAngle2 * xErr2 + cosAngle2 * yErr2);

            det->mask->data.U8[numGood] = 0;
            plateScale += scale;
            numGood++;
        }
        det->seeing *= plateScale / numGood;

        det->x->n = numGood;
        det->y->n = numGood;
        det->ra->n = numGood;
        det->dec->n = numGood;
        det->raErr->n = numGood;
        det->decErr->n = numGood;
        det->mag->n = numGood;
        det->magErr->n = numGood;
        det->extended->n = numGood;
        det->angle->n = numGood;
        det->angleErr->n = numGood;
        det->length->n = numGood;
        det->lengthErr->n = numGood;
        det->flags->n = numGood;
        det->diffSkyfileId->n = numGood;
        det->naxis1->n = numGood;
        det->naxis2->n = numGood;
        det->mask->n = numGood;

        det->num = numGood;

        if (isfinite(args->zp) && numGood > 0) {
            psBinaryOp(det->mag, det->mag, "+", psScalarAlloc(args->zp, PS_TYPE_F32));
        }

        psTrace("ppMops.read", 2, "Read %ld good rows from %s\n", numGood, (const char*)inNames->data[i]);

        psFree(table);
        detections->data[i] = det;
    }

    psTrace("ppMops.read", 1, "Done reading input detections\n");

    return detections;
}
