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

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

#include "ppMops.h"

/*
  ppMopsRead possibly modifies the args->version if the user did not
  set it explicitely.
 */
psArray *ppMopsRead(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];

        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;
        }

        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;
        }
	int skyChipPsfVersion = ppMopsGetSkyChipPsfVersion(fits);
	if (args->version == 0) {
	  psTrace("ppMops.read", 1, "Changing args->version to %d\n", skyChipPsfVersion);
	  args->version = (unsigned short) skyChipPsfVersion;
	}
	if (skyChipPsfVersion == 0) {
	  // Try to read with the user specified version?
	  skyChipPsfVersion = args->version;
	}
	/* Display a warning message if there are version
	   inconsistencies between the file and the flag (note that
	   those inconsistencies might be wanted) */
	if (skyChipPsfVersion != args->version) {
	  if (skyChipPsfVersion > args->version) {
	    psWarning("The FITS data will be downgraded from PS1_DV%d to PS1_DV%d\n",
		      skyChipPsfVersion, args->version);
	  } else { // Necessarily: skyChipPsfVersion > args->version
	    psWarning("The FITS data will be upgraded from PS1_DV%d to PS1_DV%d (new values set to default 0, NaN...)\n",
		      skyChipPsfVersion, args->version);	    
	  }
	}

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

        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 = (float) 0.5 * (psMetadataLookupF32(NULL, header, "FWHM_MAJ") +
				     psMetadataLookupF32(NULL, header, "FWHM_MIN"));

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

        psFree(header);

        psMetadata *table = psFitsReadTableAllColumns(fits); // Table of interest
        if (!table) {
            psError(PS_ERR_IO, false, "Unable to read table %d", i);
            return NULL;
        }
        det->table = table;
        psFitsClose(fits);
        if (args->version == 0) {
          if (skyChipPsfVersion < 2) {
            // XXX: TODO: Do we need to add dummy vectors for the missing columns?
           }
        }

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

        det->raErr = psVectorAlloc(size, PS_TYPE_F64);
        det->decErr = psVectorAlloc(size, PS_TYPE_F64);
        det->mask = psVectorAlloc(size, PS_TYPE_U8);

        // convert ra and dec to radians for use in the purge duplicates function
        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 NULL;
        }

        // Add our new vectors to the table so that duplicates and masked items may be purged
        psMetadataAddVector(table, PS_LIST_HEAD, "DEC_ERR", 0, NULL, det->decErr);
        psMetadataAddVector(table, PS_LIST_HEAD, "RA_ERR", 0, NULL, det->raErr);

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

        psVector *mag    = psMetadataLookupVector(NULL, table, "PSF_INST_MAG");
        psVector *magErr = psMetadataLookupVector(NULL, table, "PSF_INST_MAG_SIG");
        psVector *xErrV  = psMetadataLookupVector(NULL, table, "X_PSF_SIG");
        psVector *yErrV  = psMetadataLookupVector(NULL, table, "Y_PSF_SIG");
        psVector *scaleV = psMetadataLookupVector(NULL, table, "PLTSCALE");
        psVector *angleV = psMetadataLookupVector(NULL, table, "POSANGLE");
        psVector *flagsV = psMetadataLookupVector(NULL, table, "FLAGS");

        double plateScale = 0.0;        // Plate scale
        long numGood = 0;               // Number of good rows
        for (long row = 0; row < size; row++) {

            psU32 flags = flagsV->data.U32[row]; // psFitsTableGetU32(NULL, table, row, "FLAGS");
            if (flags & SOURCE_MASK) {
                psTrace("ppMops.read", 10, "Discarding row %ld from input %d because of flags: %ud", row, i, flags);
                det->mask->data.U8[row] = 0xFF;
                continue;
            }

            // Calculate error in RA, Dec

            double xErr = xErrV->data.F32[row];
            double yErr = yErrV->data.F32[row];
            double scale = scaleV->data.F32[row];
            double angle = angleV->data.F32[row];

            if (!isfinite(det->x->data.F32[row]) || !isfinite(det->y->data.F32[row]) ||
                !isfinite(det->ra->data.F64[row]) || !isfinite(det->dec->data.F64[row]) ||
                !isfinite(mag->data.F32[row]) || !isfinite(magErr->data.F32[row]) ||
                !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",
                        row, i,
                        det->x->data.F32[row], det->y->data.F32[row],
                        det->ra->data.F64[row], det->dec->data.F64[row],
                        mag->data.F32[row], magErr->data.F32[row],
                        xErr, yErr, scale, angle);
                det->mask->data.U8[row] = 0xFF;
                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[row] = errScale * sqrt(cosAngle2 * xErr2 + sinAngle2 * yErr2);
            det->decErr->data.F64[row] = errScale * sqrt(sinAngle2 * xErr2 + cosAngle2 * yErr2);

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

        // Are we using numGood for anything outside of this function?
        det->numGood = numGood;

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

        psTrace("ppMops.read", 2, "Read %ld good rows from %s\n", numGood, (const char*)name);
    }

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

    return detections;
}
