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

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

#include "ppMops.h"
#include "ppTranslateVersion.h"

static bool addOutputColumn(psMetadata *table, const psArray *detections, long total, char *outColName, char *inColName, bool convertTo32);
static bool addSkyfileIDColumn(psMetadata *table, const psArray *detections, long total, char *colName);

bool ppMopsWrite(const psArray *detections, const ppMopsArguments *args)
{

    psFits *fits = psFitsOpen(args->output, "w"); // FITS file
    if (!fits) {
        psError(PS_ERR_IO, false, "Unable to open output file.");
        return false;
    }

    psMetadata *header = psMetadataAlloc(); // Header to write
    psString source = ppTranslateSource(), version = ppTranslateVersion();
    psMetadataAddStr(header, PS_LIST_TAIL, "SWSOURCE", 0, "Software source", source);
    psMetadataAddStr(header, PS_LIST_TAIL, "SWVERSN", 0, "Software version", version);
    ppTranslateVersionHeader(header);
    psFree(source);
    psFree(version);

    psMetadataAddStr(header, PS_LIST_TAIL, "EXP_NAME", 0, "Exposure name", args->exp_name);
    psMetadataAddS64(header, PS_LIST_TAIL, "EXP_ID", 0, "Exposure identifier", args->exp_id);
    psMetadataAddS64(header, PS_LIST_TAIL, "CHIP_ID", 0, "Chip stage identifier", args->chip_id);
    psMetadataAddS64(header, PS_LIST_TAIL, "CAM_ID", 0, "Cam stage identifier", args->cam_id);
    psMetadataAddS64(header, PS_LIST_TAIL, "FAKE_ID", 0, "Fake stage identifier", args->fake_id);
    psMetadataAddS64(header, PS_LIST_TAIL, "WARP_ID", 0, "Warp stage identifier", args->warp_id);
    psMetadataAddS64(header, PS_LIST_TAIL, "DIFF_ID", 0, "Diff stage identifier", args->diff_id);
    psMetadataAddBool(header, PS_LIST_TAIL, "DIFF_POS", 0, "Positive subtraction?", args->positive);

    // Get these header words from the first input non null input
    ppMopsDetections *det = NULL;
    for (int d = 0; d < detections->n && det == NULL; d++) {
        det = detections->data[d];
    }
    if (det != NULL) {
        psMetadataAddF64(header, PS_LIST_TAIL, "MJD-OBS", 0, "MJD of exposure midpoint", det->mjd);
        psMetadataAddStr(header, PS_LIST_TAIL, "RA", 0, "Right Ascension of boresight", det->raBoresight);
        psMetadataAddStr(header, PS_LIST_TAIL, "DEC", 0, "Declination of boresight", det->decBoresight);
	psMetadataAddF64(header, PS_LIST_TAIL, "TEL_ALT", 0, "Telescope altitude", det->alt);
	psMetadataAddF64(header, PS_LIST_TAIL, "TEL_AZ", 0, "Telescope azimuth", det->az);
	psMetadataAddF64(header, PS_LIST_TAIL, "EXPTIME", 0, "Exposure time (sec)", det->exptime);
	psMetadataAddF64(header, PS_LIST_TAIL, "ROTANGLE", 0, "Rotator position angle", det->posangle);
	psMetadataAddStr(header, PS_LIST_TAIL, "FILTER", 0, "Filter name", det->filter);
	psMetadataAddF32(header, PS_LIST_TAIL, "AIRMASS", 0, "Airmass of exposure", det->airmass);
        psMetadataAddF32(header, PS_LIST_TAIL, "SEEING", 0, "Mean seeing", det->seeing);
    } else {
	psWarning("no inputs with surviving detections. output header will be incomplete");
    }
    psMetadataAddStr(header, PS_LIST_TAIL, "OBSCODE", 0, "IAU Observatory code", OBSERVATORY_CODE);
    psMetadataAddF32(header, PS_LIST_TAIL, "MAGZP", 0, "Magnitude zero point", args->zp);
    psMetadataAddF32(header, PS_LIST_TAIL, "MAGZPERR", 0, "Error in magnitude zero point", args->zpErr);
    psMetadataAddF32(header, PS_LIST_TAIL, "ASTRORMS", 0, "RMS of astrometric fit", args->rmsAstrom);

    //field in header that tells about the CMF version
    char cmfVersion[8];
    sprintf(cmfVersion, PS1_DV_FORMAT, args->version);
    psMetadataAddStr(header, PS_LIST_TAIL, "CMFVERSION", 0, "CMF version", cmfVersion);

    // Find the total number of detections

    long total = 0;
    for (long i=0; i<detections->n; i++) {
        ppMopsDetections *det = detections->data[i];
	if (!det) {
	    continue;
	}
        total += det->num;
    }

    psTrace("ppMops.write", 1, "Writing %ld rows to %s", total, args->output);

    if (total == 0) {
        // Write dummy table
        psMetadata *row = psMetadataAlloc(); // Output row
        psMetadataAddF64(row, PS_LIST_TAIL, "RA", 0, "Right ascension (degrees)", NAN);
        psMetadataAddF64(row, PS_LIST_TAIL, "RA_ERR", 0, "Right ascension error (degrees)", NAN);
        psMetadataAddF64(row, PS_LIST_TAIL, "DEC", 0, "Declination (degrees)", NAN);
        psMetadataAddF64(row, PS_LIST_TAIL, "DEC_ERR", 0, "Declination error (degrees)", NAN);
        psMetadataAddF32(row, PS_LIST_TAIL, "MAG", 0, "Magnitude", NAN);
        psMetadataAddF32(row, PS_LIST_TAIL, "MAG_ERR", 0, "Magnitude error", NAN);
        psMetadataAddF32(row, PS_LIST_TAIL, "PSF_CHI2", 0, "chi^2 of PSF fit", NAN);
        psMetadataAddS32(row, PS_LIST_TAIL, "PSF_DOF", 0, "Degrees of freedom of PSF fit", 0);
        psMetadataAddF32(row, PS_LIST_TAIL, "CR_SIGNIFICANCE", 0, "Significance of CR", NAN);
        psMetadataAddF32(row, PS_LIST_TAIL, "EXT_SIGNIFICANCE", 0, "Significance of extendedness", NAN);
        psMetadataAddF32(row, PS_LIST_TAIL, "PSF_MAJOR", 0, "PSF major axis (pixels)", NAN);
        psMetadataAddF32(row, PS_LIST_TAIL, "PSF_MINOR", 0, "PSF minor axis (pixels)", NAN);
        psMetadataAddF32(row, PS_LIST_TAIL, "PSF_THETA", 0, "PSF position angle (deg on chip)", NAN);
        psMetadataAddF32(row, PS_LIST_TAIL, "PSF_QUALITY", 0, "PSF quality factor", NAN);
        psMetadataAddS32(row, PS_LIST_TAIL, "PSF_NPIX", 0, "Number of pixels in PSF", 0);
        psMetadataAddF32(row, PS_LIST_TAIL, "MOMENTS_XX", 0, "xx moment", NAN);
        psMetadataAddF32(row, PS_LIST_TAIL, "MOMENTS_XY", 0, "xy moment", NAN);
        psMetadataAddF32(row, PS_LIST_TAIL, "MOMENTS_YY", 0, "yy moment", NAN);
        psMetadataAddU32(row, PS_LIST_TAIL, "FLAGS", 0, "Detection bit flags", 0);
        psMetadataAddS64(row, PS_LIST_TAIL, "DIFF_SKYFILE_ID", 0, "Identifier for diff skyfile", 0);

        psMetadataAddS32(row, PS_LIST_TAIL, "N_POS", 0, "Number of positive pixels", 0);
        psMetadataAddF32(row, PS_LIST_TAIL, "F_POS", 0, "Fraction of positive pixels", NAN);
        psMetadataAddF32(row, PS_LIST_TAIL, "RATIO_BAD", 0, "Ratio of positive pixels to negative", NAN);
        psMetadataAddF32(row, PS_LIST_TAIL, "RATIO_MASK", 0, "Ratio of positive pixels to masked", NAN);
        psMetadataAddF32(row, PS_LIST_TAIL, "RATIO_ALL", 0, "Ratio of positive pixels to all", NAN);

	if (args->version == 2) {
	      // Write data of version 2 (see ICD)
	      psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET",         PS_DATA_U32, "IPP detection identifier index",
			     NAN);
	      psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX",    PS_DATA_F32, "PSF fit instrumental magnitude",
			     NAN);
	      psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG",PS_DATA_F32, "Sigma of PSF instrumental magnitude",
			     NAN);
	      psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG",           PS_DATA_F32, "magnitude in standard aperture",
			     NAN);
	      psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RAW",       PS_DATA_F32, "magnitude in real aperture",
			     NAN);
	      psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",
			     NAN);
	      psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX",          PS_DATA_F32, "instrumental flux in standard aperture",
			     NAN);
	      psMetadataAdd (row, PS_LIST_TAIL, "AP_FLUX_SIG",      PS_DATA_F32, "aperture flux error",
			     NAN);
	      psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",
			     NAN);
	      psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",
			     NAN);
	      psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG",  PS_DATA_F32, "measured scatter of zero point calibration",
			     NAN);
	      psMetadataAdd (row, PS_LIST_TAIL, "SKY",              PS_DATA_F32, "Sky level",
			     NAN);
	      psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA",        PS_DATA_F32, "Sigma of sky level",
			     NAN);
	      psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF_PERFECT",   PS_DATA_F32, "PSF coverage/quality factor (poor)",
			     NAN);
	      psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_R1",       PS_DATA_F32, "first radial moment",
			     NAN);
	      psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_RH",       PS_DATA_F32, "half radial moment",
			     NAN);
	      psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX",        PS_DATA_F32, "Kron Flux (in 2.5 R1)",
			     NAN);
	      psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_ERR",    PS_DATA_F32, "Kron Flux Error",
			     NAN);
	      psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_INNER",  PS_DATA_F32, "Kron Flux (in 1.0 R1)", 
			     NAN);
	      psMetadataAdd (row, PS_LIST_TAIL, "KRON_FLUX_OUTER",  PS_DATA_F32, "Kron Flux (in 4.0 R1)",
			     NAN);
	      psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_P",         PS_DATA_F32, "distance to positive match source",
			     NAN);
	      psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_P",        PS_DATA_F32, "signal-to-noise of pos match src",
			     NAN);
	      psMetadataAdd (row, PS_LIST_TAIL, "DIFF_R_M",         PS_DATA_F32, "distance to negative match source",
			     NAN);
	      psMetadataAdd (row, PS_LIST_TAIL, "DIFF_SN_M",        PS_DATA_F32, "signal-to-noise of neg match src",
			     NAN);
	      psMetadataAdd (row, PS_LIST_TAIL, "FLAGS2",           PS_DATA_U32, "psphot analysis flags (group 2)",
			     0);
	      psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES",         PS_DATA_U16, "Number of frames overlapping source center", 
			     0);
	      psMetadataAdd (row, PS_LIST_TAIL, "PADDING",          PS_DATA_S16, "padding", 
			     0);
	}

        if (!psFitsWriteTableEmpty(fits, header, row, OUT_EXTNAME)) {
            psErrorStackPrint(stderr, "Unable to write empty table.");
            psFree(header);
            psFree(row);
            return false;
        }
        psFree(row);
    } else {

#define addColumn(_outName, _inName, _convertTo32) \
        if (!addOutputColumn(table, detections, total, _outName, _inName, _convertTo32)) { \
            psError(PS_ERR_UNKNOWN, false, "Failed to add column %s", _outName); \
            return false; \
        }

        // Allocate the output table
        psMetadata *table = psMetadataAlloc();
        addColumn("RA", "RA_PSF", 0);
        addColumn("RA_ERR", NULL, 0);      // calculated from various parameters including X_PSF_SIG Y_PSF_SIG and POSANG
        addColumn("DEC", "DEC_PSF", 0);
        addColumn("DEC_ERR", NULL, 0);     // calculated from various parameters including X_PSF_SIG Y_PSF_SIG and POSANG
        addColumn("MAG", "PSF_INST_MAG", 0);
        addColumn("MAG_ERR", "PSF_INST_MAG_SIG", 0);
        addColumn("PSF_CHI2", "PSF_CHISQ", 0);
        addColumn("PSF_DOF", "PSF_NDOF", 1);
        addColumn("CR_SIGNIFICANCE", "CR_NSIGMA", 0);
        addColumn("EXT_SIGNIFICANCE", "EXT_NSIGMA", 0);
        addColumn("PSF_MAJOR", NULL, 0);
        addColumn("PSF_MINOR", NULL, 0);     
        addColumn("PSF_THETA", NULL, 0);    
        addColumn("PSF_QUALITY", "PSF_QF", 0); 
        addColumn("PSF_NPIX", NULL, 1);        
        addColumn("MOMENTS_XX", NULL, 0);
        addColumn("MOMENTS_XY", NULL, 0);
        addColumn("MOMENTS_YY", NULL, 0);
        addColumn("N_POS", "DIFF_NPOS", 1);
        addColumn("F_POS", "DIFF_FRATIO", 0);
        addColumn("RATIO_BAD", "DIFF_NRATIO_BAD", 0);
        addColumn("RATIO_MASK", "DIFF_NRATIO_MASK", 0);
        addColumn("RATIO_ALL", "DIFF_NRATIO_ALL", 0);
        addColumn("FLAGS", "FLAGS", 1);
        addSkyfileIDColumn(table, detections, total, "DIFF_SKYFILE_ID");
        if (args->version == 2) {
            addColumn("IPP_IDET", NULL, 1);
            addColumn("PSF_INST_FLUX", NULL, 0);
            addColumn("PSF_INST_FLUX_SIG", NULL, 0);
            addColumn("AP_MAG", NULL, 0);
            addColumn("AP_MAG_RAW", NULL, 0);
            addColumn("AP_MAG_RADIUS", NULL, 0);
            addColumn("AP_FLUX", NULL, 0);
            addColumn("AP_FLUX_SIG", NULL, 0);
            addColumn("PEAK_FLUX_AS_MAG", NULL, 0);
            addColumn("CAL_PSF_MAG", NULL, 0);
            addColumn("CAL_PSF_MAG_SIG", NULL, 0);
            addColumn("SKY", NULL, 0);
            addColumn("SKY_SIGMA", NULL, 0);
            addColumn("PSF_QF_PERFECT", NULL, 0);
            addColumn("MOMENTS_R1", NULL, 0);
            addColumn("MOMENTS_RH", NULL, 0);
            addColumn("KRON_FLUX", NULL, 0);
            addColumn("KRON_FLUX_ERR", NULL, 0);
            addColumn("KRON_FLUX_INNER", NULL, 0);
            addColumn("KRON_FLUX_OUTER", NULL, 0);
            addColumn("DIFF_R_P", NULL, 0);
            addColumn("DIFF_SN_P", NULL, 0);
            addColumn("DIFF_R_M", NULL, 0);
            addColumn("DIFF_SN_M", NULL, 0);
            addColumn("FLAGS2", NULL, 1);
            addColumn("IPP_IDET", NULL, 0);
            addColumn("N_FRAMES", NULL, 0);
            addColumn("PADDING", NULL, 0);
        }
        if (!psFitsWriteTableAllColumns(fits, header, table, OUT_EXTNAME)) {
            psError(psErrorCodeLast(), false, "Unable to write table");
            return false;
        }
        psFree(table);
    }

    psFree(header);
    psFitsClose(fits);

    psTrace("ppMops.write", 1, "Done writing %ld rows to %s", det->num, args->output);

    return true;
}

static bool addOutputColumn(psMetadata *table, const psArray *detections, long outputSize, char *outColumnName, char *inColumnName, bool convertTo32)
{
    if (inColumnName == NULL) {
        inColumnName = outColumnName;
    }

    psVector *out = NULL;
    if (convertTo32) {
        // psFitsReadTableAllColumns reads columns of cfitsio type LONG and ULONG into a 64 bit integers
        // We want to write 32 bits to the output.
        int next = 0;
        for (long i=0; i<detections->n; i++) {
            ppMopsDetections *det = detections->data[i];
            if (!det || det->num == 0) {
                // no detections survived for this input
                continue;
            }
            psVector *in = psMetadataLookupVector(NULL, det->table, inColumnName);
            if (!in) {
                psError(PS_ERR_PROGRAMMING, true, "failed to find input column: %s", inColumnName);
                return false;
            }
            if (in->type.type != PS_TYPE_S64 && in->type.type != PS_TYPE_U64) {
                psError(PS_ERR_PROGRAMMING, true, "input column to convert is not S64 or U64: %s %d",
                    inColumnName, in->type.type);
                return false;
            }
            if (out == NULL) {
                // First time through set up the output vector and the copy parameters
                if (in->type.type == PS_TYPE_S64) {
                    out = psVectorAlloc(outputSize, PS_TYPE_S32);
                } else {
                    out = psVectorAlloc(outputSize, PS_TYPE_U32);
                }
            }
            for (long d=0; d < det->num; d++) {
                if (in->type.type == PS_TYPE_S64) {
                    out->data.S32[next++] = in->data.S64[d];
                } else {
                    out->data.U32[next++] = in->data.U64[d];
                }
            }
        }
    } else {
        void *next = NULL;
        int elementSize = 0;    // size of elements in vector... We are making assumptions here about the organization of primitives in memory so we can use memcopy
        for (long i=0; i<detections->n; i++) {
            ppMopsDetections *det = detections->data[i];
            if (!det || det->num == 0) {
                // no detections survived for this input
                continue;
            }
            psVector *in = psMetadataLookupVector(NULL, det->table, inColumnName);
            if (!in) {
                psError(PS_ERR_PROGRAMMING, true, "failed to find input column: %s", inColumnName);
                return false;
            }
            if (out == NULL) {
                // First time through set up the output vector and the copy parameters
                out = psVectorAlloc(outputSize, in->type.type);
                next = (void *) out->data.U8;
                switch (in->type.type) {
                    case PS_TYPE_S8:
                    elementSize = sizeof(psS8);
                    break;
                case PS_TYPE_U8:
                    elementSize = sizeof(psU8);
                    break;
                case PS_TYPE_S16:
                    elementSize = sizeof(psS16);
                    break;
                case PS_TYPE_U16:
                    elementSize = sizeof(psU16);
                    break;
                case PS_TYPE_S32:
                    elementSize = sizeof(psS32);
                    break;
                case PS_TYPE_U32:
                    elementSize = sizeof(psU32);
                    break;
                case PS_TYPE_S64:
                    elementSize = sizeof(psS64);
                    break;
                case PS_TYPE_U64:
                    elementSize = sizeof(psU64);
                    break;
                case PS_TYPE_F32:
                    elementSize = sizeof(psF32);
                    break;
                case PS_TYPE_F64:
                    elementSize = sizeof(psF64);
                    break;
                default:
                    psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Unknown vector type %d", in->type.type);
                    return false;

                }
            }
            // We are doing nasty things here so we can use memcpy. 
            // It would be safer to do a proper loop over the elements.
            long toCopy = det->num * elementSize;
            memcpy(next, in->data.U8, toCopy);
            next += toCopy;
        }
    }

    // Finally add the new column to the output table
    psMetadataAddVector(table, PS_LIST_TAIL, outColumnName, 0, NULL, out);
    psFree(out);    // drop reference

    return true;
}
static bool addSkyfileIDColumn(psMetadata *table, const psArray *detections, long total, char *colName)
{
    psVector *out = psVectorAlloc(total, PS_TYPE_S64);
    long next = 0;
    for (long i = 0; i<detections->n; i++) {
        ppMopsDetections *det = detections->data[i];
	if (!det) {
	    continue;
	}
        psS64 diffSkyfileId = det->diffSkyfileId;
        for (long j = 0; j < det->num; j++) {
            out->data.S64[next++] = diffSkyfileId;
        }
    }
    psMetadataAddVector(table, PS_LIST_TAIL, colName, 0, NULL, out);
    psFree(out);
    return true;
}
