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

// Configuration for Phase 4
typedef struct {
    // Input parameters
    int nInputs;			// Number of input files
    const char *camera;			// Name of the camera configuration file
    const psList *inNames;		// Input file names
    unsigned long skyId;		// Sky tangent plane identification number
    // Output parameters
    const char *combined;		// Combined image filename
    const char *subtracted;		// Subtracted image filename
    const char *variables;		// Output variable candidates filename
    const char *recipeFile;		// Recipe filename
} phase4Config;


// A tangent plane (for combined image and static sky)
// Assume we're dealing with a rectangle only for now.
typedef struct {
    long int id;			// Identification number --- key in MDDB
    double ra0, dec0;			// Tangent point --- from MDDB
    double sizeRA, sizeDec;		// Size in RA and Dec --- from MDDB
    double scale;			// Pixel scale --- from MDDB
    char *name;				// File name --- from MDDB or image server
    psImage *pixels;			// The pixels that comprise the tangent plane
    psImage *mask;			// Mask image
    psSomeCollectionOfObjects *objects;	// Objects in the image
} psTangentPlane;

void main (int argc, char **argv)
{
    // Bury the command line parsing into some function that just works
    phase4Config *config = parseCommandLine(argc, argv); // Parse the command line

    // Read the recipe metadata
    psMetadata *recipe = psMetadataParseConfig(NULL, config->recipeFile, false);

    // Don't particularly care about the camera and site details --- Phase 4 doesn't care where
    // the data comes from, but it does need the following for each image:
    // * Gain and read noise for noise model --- from the camera configuration
    // * Background value for noise model (assume constant over the region of interest) --- from image MD
    // * Photometric calibration for relative and absolute scalings --- from image MD
    // * Astrometric calibration to map to the sky --- from file specified in image MD
    // We need to be told which camera; assume all data is from the same camera (easy to generalise
    // if this is not true).

    // Read the camera configuration
    psMetadata *camera = psMetadataParseConfig(NULL, config->camera, false);
    // What are the keywords for gain, readnoise, photometric calibration, background?
    // This basically treats the "camera" configuration as a translation table.
    const char *gainHeader = psMetadataLookup(camera, "GAIN");
    const char *readnoiseHeader = psMetadataLookup(camera, "READNOISE");
    const char *photcalHeader = psMetadataLookup(camera, "PHOTCAL");
    const char *backgroundHeader = psMetadataLookup(camera, "BACKGRND");
    const char *astromHeader = psMetadataLookup(camera, "ASTROM");

    // Initialise the various arrays/vectors we need
    psArray *headers = psArrayAlloc(config->nInputs); // The image headers
    psArray *inputs = psArrayAlloc(config->nInputs); // The images
    psArray *astrometry = psArrayAlloc(config->nInputs); // The astrometry data for each image
    psVector *gain = psVectorAlloc(config->nInputs, PS_TYPE_F32); // Gains
    psVector *readnoise = psVectorAlloc(config->nInputs, PS_TYPE_F32); // Read noises
    psVector *photcal = psVectorAlloc(config->nInputs, PS_TYPE_F32); // Photometric calibrations
    psVector *background = psVectorAlloc(config->nInputs, PS_TYPE_F32); // Background flux values

    // This function queries the MDDB to get sizeRA,sizeDec,ra0,dec0,scale,name for the tangent plane of
    // interest
    psTangentPlane *staticSky = psTangentPlaneFromDB(config->sky);

    // Read in each of the inputs
    for (psListSetIterator(config->inNames, PS_LIST_HEAD), // Set the iterator and get the first one
	     char *inName = psListGetNext(input),
	     int i = 0;
	 inName != NULL;		// Keep going until nothing left in the list
	 inName = psListGetNext(input),	// Get the next one on the list
	     i++) {

	// Read the header
	psMetadata *header = psFITSReadHeader(inName);
	// Put it in the array of headers for future reference
	headers->data[i] = header;

	// Check the validity of the given header based on expectactions for this camera
	status = pmCameraValidateHeaders(headers, camera);
	// Need to act if status is bad --- generate a warning and ignore this image

	// Get the gain, readnoise, photometric calibration, background
	// I think this is not strictly correct --- doesn't psMetadataLookup return a void* ?
	gain->data.F32[i] = (psF32) psMetadataLookup(header, gainHeader);
	readnoise->data.F32[i] = (psF32) psMetadataLookup(header, readnoiseHeader);
	photcal->data.F32[i] = (psF32) psMetadataLookup(header, photcalHeader);
	background->data.F32[i] = (psF32) psMetadataLookup(header, backgroundHeader);

	// Get the name of the astrometry file (XML) and generate the FPA structure from it (no pixels)
	const char *astromFile = psMetadataLookup(header, astromHeader);
	psFPA *fpa = psFPAFromAstrometry(astromFile);
	astrometry->data[i] = fpa;

	// Iterate over the FPA structure
	for (int chipnum = 0, int nChips = fpa->chips->n; chipnum < nChips; chipnum++) {
	    psChip *chip = fpa->chips->data[chipnum]; // The chip of interest
	    for (int cellnum = 0, int nCells = chip->cells->n; cellnum < nCells; cellnum++) {
		psCell *cell = chip->cells->data[cellnum]; // The cell of interest

		// QUESTION: Do we want to bother about the cells that were used for fast guiding?
		// Assume that we care --- may as well.
		
		// QUESTION: Are multiple readouts stacked in phase 4 or earlier?
		// Assume that they should be stacked here.

		// QUESTION: How do we know the names of the chip, cell, readout in the FITS file?
		// Assume that each of psFPA, psChip, psCell, psReadout contains a "const char *source"
		// which contains some sort of URI, e.g.:
		// FPA: fpa12345_12.fits --- This particular FPA contains only a single OTA
		// Chip: fpa12345_12.fits --- The OTA is entirely contained in the FITS file
		// Cell: fpa12345_12.fits|cell34 --- The cell has extension name "cell34"
		// Readout: fpa12345_12.fits|cell34|56 --- The readout is the 56th slice in the 3rd axis
		// Assume further that, given the "source", there exists a module that loads the pixels.


		// If there are multiple readouts, then stack them
		if (cell->readouts->n > 1) {
		    psList *readouts = psListAlloc(NULL); // List of readouts for stacking

		    // Read the images in
		    for (int readoutNum = 0; readoutNum < cell->readouts->n; readoutNum++) {
			(void)psReadoutRead(cell->readouts->data[readoutNum]); // Read the readout into memory
			psListAdd(readouts, PS_LIST_TAIL, cell->readouts->data[readoutNum]); // Put on list
		    }

		    // QUESTION: Do we want to bother about combining these "properly", or just sum them?
		    // Assume that we just sum them, which saves measuring the zeros and scales.  CR rejection
		    // can be done when we combine multiple images (consecutive or from multiple telescopes).

		    // Perhaps we want a pmReadoutStack module, which would convert an array of readouts from
		    // fast guiding to a single pseudo-readout.

		    // Just do a simple mean
		    psStats *combineStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); // Statistics for combination
		    pmCombineParams *combine = pmCombineParamsAlloc(combineStats, 0, 0.0, 0.0, 0); // How to
												   // combine
		    psImage *stack = pmReadoutCombine(NULL, readouts, combine, 0.0, 1.0, false, 1.0, 0.0);
		    // And now multiply by the number to get the sum
		    (void)psImageBinaryOp(stack, stack, "*", psScalarAlloc(cell->readouts->n, PS_TYPE_U16));

		    // Probably also have to multiply by some exposure time ratio to scale with the rest of
		    // the focal plane --- neglected here.

		    // Clean up the temporary steps
		    psFree(readouts);
		    psFree(combine);
		    psFree(combineStats);

		    // Stuff the stack back into the cell --- probably should play with the metadata, but I'm
		    // going to ignore that for now.

		    // Clean out the old structure
		    for (int readoutNum = 0; readoutNum < cell->readouts->n; readoutNum++) {
			psFree(cell->readouts->data[readoutNum]);
		    }
		    // Resize and stuff the stack in.
		    psArrayRealloc(cell->readouts, 1);
		    cell->readouts->data[0] = stack;
		} else {
		    // Otherwise, just read the cell in
		    cell->readouts->data[0] = psReadoutRead(cell->readouts->data[0]);
		} // Reading in readouts
	    } // Cells
	} // Chips
    } // FPAs

    // Now all the data is in memory.  Combine the multiple images
    psTangentPlane *combined = pmFPACombine(NULL, staticSky, inputs, gain, readnoise, photcal, background);

    // Write the combined image out
    psTangentPlaneWriteFITS(combined, config->combined);

    // Run the object detection/measurement --- this is not well thought out yet
    (void)psDetectAndMeasureObjects(combined, otherParameters);

    // Load the static sky into the structure already defined
    (void)psTangentPlaneRead(staticSky);

    // Get the seeing ratio --- involves matching stars and measuring the mean ratio of FWHM
    float seeingRatio = pmSeeingRatio(combined, staticSky, otherParameters);

    // Do the subtraction --- involves measuring the convolution kernel, applying the kernel and subtracting
    psTangentPlane *subtracted = NULL;
    if (seeingRatio < 1.0) {
	subtracted = pmImageSubtraction(NULL, combined, staticSky, otherParameters);
    } else {
	subtracted = pmImageSubtraction(NULL, staticSky, combined, otherParameters);
    }

    // Write out the subtracted image
    psTangentPlaneWriteFITS(subtraction, config->subtracted);
    
    // Find and measure objects on the subtracted image
    (void)psDetectAndMeasureObjects(subtracted, otherParameters);

    // Filter the list of candidate variable sources
    psList *variables = pmFilterVariables(subtracted, combined, staticSky, otherParameters);

    // Output the list of variable sources
    (void)pmOutputVariables(variables, config->variables);

    // All done.
}


