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

#include <stdio.h>
#include <string.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++) {
        const char *name = inNames->data[i]; // File name
        psFits *fits = psFitsOpen(name, "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;
        }

        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 = detections->data[i] = ppMopsDetectionsAlloc();
        det->component = psStringNCopy(name, strrchr(name, '.') - name); // Strip off extension
        det->header = header;
        det->num = size;

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

        det->naxis1 = psMetadataLookupS32(NULL, header, "IMNAXIS1");
        det->naxis2 = psMetadataLookupS32(NULL, header, "IMNAXIS2");

        det->psfHeader = psFitsReadHeader(NULL, fits);
        if (!det->psfHeader) {
            psError(psErrorCodeLast(), false, "Unable to read header %d", i);
            return false;
        }
        psMetadata *table = det->table = psFitsReadTableAllColumns(fits); // Table of interest
        if (!table) {
            psError(psErrorCodeLast(), false, "Unable to read table %d", i);
            return false;
        }

        psVector *ra = psMetadataLookupVector(NULL, table, "RA_PSF");
        psVector *dec = psMetadataLookupVector(NULL, table, "DEC_PSF");

        det->ra = (psVector*)psBinaryOp(NULL, ra, "*", psScalarAlloc(DEG_TO_RAD(1.0), PS_TYPE_F64));
        det->dec = (psVector*)psBinaryOp(NULL, dec, "*", psScalarAlloc(DEG_TO_RAD(1.0), PS_TYPE_F64));
        det->x = psMemIncrRefCounter(psMetadataLookupVector(NULL, table, "X_PSF"));
        det->y = psMemIncrRefCounter(psMetadataLookupVector(NULL, table, "Y_PSF"));
        if (!det->ra || !det->dec || !det->x || !det->y) {
            psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find all of RA, Dec, X and Y columns");
            return false;
        }

        psTrace("ppMops.read", 2, "Read %ld rows from %s\n", det->num, det->component);

        if (!psFitsMoveExtName(fits, "SkyChip.deteff")) {
            psWarning("No detection efficiencies included in %s", det->component);
            psErrorStackPrint(stderr, "No detection efficiencies included in %s", det->component);
            psErrorClear();
            continue;
        }

        det->deteffHeader = psFitsReadHeader(NULL, fits);
        if (!det->deteffHeader) {
            psWarning("Unable to read detection efficiency header in %s", det->component);
            psErrorClear();
            continue;
        }
        det->deteffTable = psFitsReadTableAllColumns(fits);
        if (!det->deteffTable) {
            psFree(det->deteffHeader);
            det->deteffHeader = NULL;
            psWarning("Unable to read detection efficiency table in %s", det->component);
            psErrorClear();
            continue;
        }
        psTrace("ppMops.read", 2, "Read detection efficiency from %s\n", det->component);
        psFitsClose(fits);
    }

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

    return detections;
}
