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

#include "psPolynomialMD.h"

#include "pmHDU.h"
#include "pmFPA.h"
#include "pmHDUUtils.h"
#include "pmFPARead.h"
#include "pmFPAWrite.h"
#include "pmReadoutStack.h"

#include "pmDark.h"

#define PM_DARK_FITS_EXTNAME "PS_DARK"  // FITS extension name for ordinates table
#define PM_DARK_FITS_NAME    "NAME"     // Column name for concept name in ordinates table
#define PM_DARK_FITS_ORDER   "ORDER"    // Column name for polynomial order in ordinates table
#define PM_DARK_FITS_SCALE   "SCALE"    // Column name for scaling option in ordinates table
#define PM_DARK_FITS_MIN     "MIN"      // Column name for minimum value in ordinates table
#define PM_DARK_FITS_MAX     "MAX"      // Column name for maximum value in ordinates table



// Look up the value of an ordinate in a readout
static bool ordinateLookup(float *value, // Value of ordinate, to return
                           const char *name, // Name of ordinate (concept name)
                           bool scale,  // Scale the value?
                           float min, float max, // Minimum and maximum values for scaling
                           const pmReadout *ro // Readout of interest
                           )
{
    assert(value);
    assert(name);
    assert(ro);

    pmCell *cell = ro->parent; // Parent cell
    if (!cell) {
        return false;
    }
    psMetadataItem *item = psMetadataLookup(cell->concepts, name);
    if (!item) {
        pmChip *chip = cell->parent; // Parent chip
        if (!chip) {
            return false;
        }
        item = psMetadataLookup(chip->concepts, name);
        if (!item) {
            pmFPA *fpa = chip->parent; // Parent FPA
            if (!fpa) {
                return false;
            }
            item = psMetadataLookup(fpa->concepts, name);
            if (!item) {
                psWarning("Unable to find concept %s in readout", name);
                return false;
            }
        }
    }

    *value = psMetadataItemParseF32(item); // Value of interest
    if (!isfinite(*value)) {
        psWarning("Non-finite value (%f) of concept %s in readout", *value, name);
        return false;
    }
    if (scale) {
        if (*value < min || *value > max) {
            psWarning("Value of concept %s (%f) outside range (%f:%f)", name, *value, min, max);
            return false;
        }
        *value = 2.0 * (*value - min) / (max - min) - 1.0;
    }

    return true;
}

static void darkOrdinateFree(pmDarkOrdinate *ord)
{
    psFree(ord->name);
    return;
}

pmDarkOrdinate *pmDarkOrdinateAlloc(const char *name, int order)
{
    pmDarkOrdinate *ord = psAlloc(sizeof(pmDarkOrdinate)); // Ordinate data, to return
    psMemSetDeallocator(ord, (psFreeFunc)darkOrdinateFree);

    ord->name = psStringCopy(name);
    ord->order = order;
    ord->scale = false;
    ord->min = NAN;
    ord->max = NAN;

    return ord;
}


bool pmDarkCombine(pmCell *output, const psArray *inputs, psArray *ordinates, const char *normConcept,
                   int iter, float rej, psMaskType maskVal)
{
    PS_ASSERT_PTR_NON_NULL(output, false);
    PS_ASSERT_ARRAY_NON_NULL(inputs, false);
    PS_ASSERT_ARRAY_NON_NULL(ordinates, false);
    PS_ASSERT_INT_NONNEGATIVE(iter, false);
    PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, false);

    // Extract fitting orders
    int numOrdinates = ordinates->n;    // Number of ordinates
    int numInputs = inputs->n;          // Number of inputs
    int numBadInputs = 0;               // Number of bad inputs
    psArray *values = psArrayAlloc(numInputs);
    psVector *roMask = psVectorAlloc(numInputs, PS_TYPE_U8); // Mask for bad readouts
    psVectorInit(roMask, 0);
    psVector *norm = normConcept ? psVectorAlloc(numInputs, PS_TYPE_F32) : NULL; // Normalisations for each
    for (int i = 0; i < numInputs; i++) {
        values->data[i] = psVectorAlloc(numOrdinates, PS_TYPE_F32);
        if (norm) {
            pmReadout *ro = inputs->data[i]; // Readout of interest
            float normValue;            // Normalisation value
            if (!ordinateLookup(&normValue, normConcept, false, NAN, NAN, ro)) {
                psWarning("Unable to find value of %s for readout %d", normConcept, i);
                roMask->data.U8[i] = 0xff;
                norm->data.F32[i] = NAN;
                numBadInputs++;
                continue;
            }
            if (normValue == 0.0) {
                psWarning("Normalisation value (%s) for readout %d is zero", normConcept, i);
                roMask->data.U8[i] = 0xff;
                norm->data.F32[i] = NAN;
                numBadInputs++;
                continue;
            }
            norm->data.F32[i] = 1.0 / normValue;
        }
    }
    psVector *orders = psVectorAlloc(numOrdinates, PS_TYPE_U8); // Orders for each concept
    for (int i = 0; i < numOrdinates; i++) {
        pmDarkOrdinate *ord = ordinates->data[i]; // Ordinate information
        if (ord->order <= 0) {
            psError(PS_ERR_UNKNOWN, true, "Bad order for concept %s (%d) --- ignored", ord->name, ord->order);
            psFree(values);
            psFree(roMask);
            psFree(orders);
            return false;
        }
        orders->data.U8[i] = ord->order;

        // Mask the readout and move on
        #define MASK_READOUT_VALUE { \
            roMask->data.U8[j] = 0xff; \
            val->data.F32[i] = NAN; \
            numBadInputs++; \
            continue; \
        }

        for (int j = 0; j < numInputs; j++) {
            psVector *val = values->data[j]; // Value vector for readout
            if (roMask->data.U8[j]) {
                val->data.F32[i] = NAN;
                continue;
            }

            pmReadout *ro = inputs->data[j]; // Readout of interest
            float value = NAN;          // Value of ordinate
            if (!ordinateLookup(&value, ord->name, ord->scale, ord->min, ord->max, ro)) {
                roMask->data.U8[j] = 0xff;
                val->data.F32[i] = NAN;
                numBadInputs++;
                continue;
            }
            val->data.F32[i] = value;
        }
    }

    if (psTraceGetLevel("psModules.detrend") > 9) {
        for (int i = 0; i < numInputs; i++) {
            psVector *val = values->data[i];
            (void) val; // avoid unused variable message when tracing is compiled out
            for (int j = 0; j < numOrdinates; j++) {
                psTrace("psModules.detrend", 9, "Image %d, ordinate %d: %f\n", i, j, val->data.F32[j]);
            }
        }
    }

    psPolynomialMD *poly = psPolynomialMDAlloc(orders); // Polynomial for fitting
    psFree(orders);
    int numTerms = poly->coeff->n;      // Number of terms in polynomial
    if (numTerms > numInputs - numBadInputs) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Insufficient inputs (%d) to fit polynomial terms (%d).",
                numInputs - numBadInputs, numTerms);
        psFree(values);
        psFree(roMask);
        psFree(norm);
        return false;
    }

    // Set up output
    int minInputCols, maxInputCols, minInputRows, maxInputRows; // Smallest and largest values to combine
    int xSize, ySize;                   // Size of the output image
    if (!pmReadoutStackValidate(&minInputCols, &maxInputCols, &minInputRows, &maxInputRows, &xSize, &ySize,
                                inputs)) {
        psError(PS_ERR_UNKNOWN, false, "No valid input readouts.");
        psFree(values);
        psFree(roMask);
        psFree(norm);
        return false;
    }
    if (output->readouts->n != numTerms) {
        output->readouts = psArrayRealloc(output->readouts, numTerms);
    }
    int outRow0 = 0, outCol0 = 0;       // Output row0 and col0
    for (int i = 0; i < numTerms; i++) {
        pmReadout *ro = output->readouts->data[i]; // Readout to update
        if (!ro) {
            ro = output->readouts->data[i] = pmReadoutAlloc(output);
        }

        pmReadoutUpdateSize(ro, minInputCols, minInputRows, xSize, ySize, false, false, 0);
        if (i == 0) {
            outRow0 = ro->row0;
            outCol0 = ro->col0;
        } else {
            assert(ro->row0 == outRow0);
            assert(ro->col0 == outCol0);
        }
    }

    psImage *counts = pmReadoutAnalysisImage(output->readouts->data[0], PM_READOUT_STACK_ANALYSIS_COUNT,
                                             xSize, ySize, PS_TYPE_U16, 0);
    if (!counts) {
        return false;
    }
    psImage *sigma = pmReadoutAnalysisImage(output->readouts->data[0], PM_READOUT_STACK_ANALYSIS_SIGMA,
                                            xSize, ySize, PS_TYPE_F32, NAN);
    if (!sigma) {
        psFree(counts);
        return false;
    }

    // Iterate over pixels, fitting polynomial
    psVector *pixels = psVectorAlloc(numInputs, PS_TYPE_F32); // Stack of pixels
    psVector *mask   = psVectorAlloc(numInputs, PS_TYPE_MASK); // Mask for stack
    for (int i = minInputRows; i < maxInputRows; i++) {
        int yOut = i - outRow0; // y position on output readout
#ifdef SHOW_BUSY
        if (psTraceGetLevel("psModules.detrend") > 9) {
            printf("Processing row %d\r", i);
            fflush(stdout);
        }
#endif

        for (int j = minInputCols; j < maxInputCols; j++) {
            int xOut = j - outCol0; // x position on output readout

            psVectorInit(mask, 0);
            for (int r = 0; r < inputs->n; r++) {
                if (roMask->data.U8[r]) {
                    mask->data.PS_TYPE_MASK_DATA[r] = 0xff;
                    continue;
                }
                pmReadout *readout = inputs->data[r]; // Input readout
                int yIn = i - readout->row0; // y position on input readout
                int xIn = j - readout->col0; // x position on input readout

                pixels->data.F32[r] = readout->image->data.F32[yIn][xIn];
                if (norm) {
                    pixels->data.F32[r] *= norm->data.F32[r];
                }
                if (readout->mask) {
                    mask->data.PS_TYPE_MASK_DATA[r] = readout->mask->data.PS_TYPE_MASK_DATA[yIn][xIn];
                }
            }

            if (!psPolynomialMDClipFit(poly, pixels, NULL, mask, maskVal, values, iter, rej)) {
                psErrorClear();         // Nothing we can do about it
                psVectorInit(poly->coeff, NAN);
            }
            for (int k = 0; k < numTerms; k++) {
                pmReadout *ro = output->readouts->data[k]; // Readout of interest
                ro->image->data.F32[yOut][xOut] = poly->coeff->data.F64[k];
            }
            counts->data.U16[yOut][xOut] = poly->numFit;
            sigma->data.F32[yOut][xOut] = poly->stdevFit;
        }
    }

    psFree(norm);
    psFree(roMask);
    psFree(poly);
    psFree(pixels);
    psFree(mask);
    psFree(values);

    psMetadataAddPtr(output->analysis, PS_LIST_TAIL, PM_DARK_ANALYSIS_ORDINATES,
                     PS_DATA_ARRAY | PS_META_REPLACE, "Dark ordinates", ordinates);
    psMetadataAddStr(output->analysis, PS_LIST_TAIL, PM_DARK_ANALYSIS_NORM,
                     PS_META_REPLACE, "Dark normalisation", normConcept);

    for (int i = 0; i < numTerms; i++) {
        pmReadout *ro = output->readouts->data[i]; // Readout to update
        ro->data_exists = true;
    }
    output->data_exists = true;
    output->parent->data_exists = true;

    return true;
}



bool pmDarkApply(pmReadout *readout, const pmCell *dark, psMaskType bad)
{
    PS_ASSERT_PTR_NON_NULL(readout, false);
    PS_ASSERT_PTR_NON_NULL(dark, false);
    PS_ASSERT_IMAGE_NON_NULL(readout->image, false);
    PS_ASSERT_IMAGE_TYPE(readout->image, PS_TYPE_F32, false);
    int numCols = readout->image->numCols, numRows = readout->image->numRows; // Size of image
    if (readout->mask) {
        PS_ASSERT_IMAGE_NON_NULL(readout->mask, false);
        PS_ASSERT_IMAGES_SIZE_EQUAL(readout->mask, readout->image, false);
        PS_ASSERT_IMAGE_TYPE(readout->mask, PS_TYPE_MASK, false);
    }
    int numTerms = dark->readouts->n;   // Number of polynomial terms
    for (int i = 0; i < numTerms; i++) {
        pmReadout *ro = dark->readouts->data[i]; // Dark readout
        PS_ASSERT_PTR_NON_NULL(ro, false);
        PS_ASSERT_IMAGE_NON_NULL(ro->image, false);
        PS_ASSERT_IMAGE_SIZE(ro->image, numCols, numRows, false);
        PS_ASSERT_IMAGE_TYPE(ro->image, PS_TYPE_F32, false);
    }
    psArray *ordinates = psMetadataLookupPtr(NULL, dark->analysis, PM_DARK_ANALYSIS_ORDINATES); // Ordinates
    if (!ordinates) {
        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find dark ordinates.");
        return false;
    }
    bool mdok;                          // Status of MD lookup
    psString normConcept = psMetadataLookupStr(&mdok, dark->analysis, PM_DARK_ANALYSIS_NORM); // Normalisation

    int numOrdinates = ordinates->n;    // Number of ordinates
    psVector *values = psVectorAlloc(numOrdinates, PS_TYPE_F32); // Values of ordinates
    for (int i = 0; i < numOrdinates; i++) {
        pmDarkOrdinate *ord = ordinates->data[i]; // Ordinate of interest
        float value = NAN;              // Value for ordinate
        if (!ordinateLookup(&value, ord->name, ord->scale, ord->min, ord->max, readout)) {
            psError(PS_ERR_UNKNOWN, true, "Unable to find value for %s", ord->name);
            psFree(values);
            return false;
        }
        values->data.F32[i] = value;
    }
    float norm = NAN;                   // Normalisation value
    bool doNorm = false;                // Do normalisation?
    if (normConcept && strlen(normConcept) > 0) {
        if (!ordinateLookup(&norm, normConcept, false, NAN, NAN, readout)) {
            psError(PS_ERR_UNKNOWN, true, "Unable to find value for %s", normConcept);
            psFree(values);
            return false;
        }
        doNorm = true;
    }

    psVector *orders = psVectorAlloc(numOrdinates, PS_TYPE_U8); // Order for each polynomial
    for (int i = 0; i < numOrdinates; i++) {
        pmDarkOrdinate *ord = ordinates->data[i]; // Ordinate information
        if (ord->order <= 0) {
            psError(PS_ERR_UNKNOWN, true, "Bad order for concept %s (%d) --- ignored", ord->name, ord->order);
            psFree(values);
            psFree(orders);
            return false;
        }
        orders->data.U8[i] = ord->order;
    }

    psPolynomialMD *poly = psPolynomialMDAlloc(orders); // Polynomial to apply
    psFree(orders);

    for (int y = 0; y < numRows; y++) {
        for (int x = 0; x < numCols; x++) {
            for (int i = 0; i < numTerms; i++) {
                pmReadout *ro = dark->readouts->data[i]; // Dark readout
                poly->coeff->data.F64[i] = ro->image->data.F32[y][x];
            }
            float value = psPolynomialMDEval(poly, values); // Value of dark current
            if (doNorm) {
                value *= norm;
            }
            readout->image->data.F32[y][x] -= value;
            if (readout->mask && !isfinite(value)) {
                readout->mask->data.PS_TYPE_MASK_DATA[y][x] = bad;
            }
        }
    }

    psFree(poly);
    psFree(values);

    return true;
}

bool pmFPAWriteDark(pmFPA *fpa, psFits *fits, psDB *db, bool blank, bool recurse)
{
    PS_ASSERT_PTR_NON_NULL(fpa, false);
    PS_ASSERT_FITS_NON_NULL(fits, false);

    if (!pmFPAWrite(fpa, fits, db, blank, recurse)) {
        psError(PS_ERR_IO, false, "Unable to write FPA dark images");
        return false;
    }

    pmHDU *hdu = fpa->hdu;
    if (blank || !hdu) {
        // No more to do
        return true;
    }

    psArray *ordinates = NULL;      // Dark ordinates, to write
    const char *normConcept = NULL;     // Normalisation concept
    psArray *chips = fpa->chips;    // Component chips
    for (int i = 0; i < chips->n; i++) {
        pmChip *chip = chips->data[i]; // Chip of interest
        psArray *cells = chip->cells; // Component cells
        for (int j = 0; j < cells->n; j++) {
            pmCell *cell = cells->data[j]; // Cell of interest
            bool mdok;              // Status of MD lookup
            psArray *newOrd = psMetadataLookupPtr(&mdok, cell->analysis, PM_DARK_ANALYSIS_ORDINATES); // Ords
            if (!mdok) {
                continue;
            }
            psString newNorm = psMetadataLookupPtr(&mdok, cell->analysis, PM_DARK_ANALYSIS_NORM); // Norm
            if (ordinates) {
                if (newOrd != ordinates) {
                    psError(PS_ERR_UNKNOWN, true, "Dark ordinates differ across cells.");
                    return false;
                }
                if ((normConcept && (!newNorm || strcmp(normConcept, newNorm) != 0)) ||
                    (!normConcept && newNorm)) {
                    psError(PS_ERR_UNKNOWN, true, "Dark normalisations differ across cells.");
                    return false;
                }
            } else {
                ordinates = newOrd;
                normConcept = newNorm;
            }
        }
    }

    if (!ordinates) {
        psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find dark ordinates.");
        return false;
    }

    return pmDarkWrite(fits, hdu->header, ordinates, normConcept);
}

bool pmChipWriteDark(pmChip *chip, psFits *fits, psDB *db, bool blank, bool recurse)
{
    PS_ASSERT_PTR_NON_NULL(chip, false);
    PS_ASSERT_FITS_NON_NULL(fits, false);

    if (!pmChipWrite(chip, fits, db, blank, recurse)) {
        psError(PS_ERR_IO, false, "Unable to write chip dark images");
        return false;
    }

    pmHDU *hdu = pmHDUFromChip(chip);
    if (blank || !hdu) {
        // No more to do
        return true;
    }

    psArray *ordinates = NULL;          // Dark ordinates, to write
    const char *normConcept = NULL;     // Normalisation concept
    psArray *cells = chip->cells;       // Component cells
    for (int j = 0; j < cells->n; j++) {
        pmCell *cell = cells->data[j]; // Cell of interest
        bool mdok;                      // Status of MD lookup
        psArray *newOrd = psMetadataLookupPtr(&mdok, cell->analysis, PM_DARK_ANALYSIS_ORDINATES); // Ordinates
        if (!mdok) {
            continue;
        }
        psString newNorm = psMetadataLookupPtr(&mdok, cell->analysis, PM_DARK_ANALYSIS_NORM); // Normalisation
        if (ordinates) {
            if (newOrd != ordinates) {
                psError(PS_ERR_UNKNOWN, true, "Dark ordinates differ across cells.");
                return false;
            }
            if ((normConcept && (!newNorm || strcmp(normConcept, newNorm) != 0)) ||
                (!normConcept && newNorm)) {
                psError(PS_ERR_UNKNOWN, true, "Dark normalisations differ across cells.");
                return false;
            }
        } else {
            ordinates = newOrd;
            normConcept = newNorm;
        }
    }

    if (!ordinates) {
        psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find dark ordinates.");
        return false;
    }

    return pmDarkWrite(fits, hdu->header, ordinates, normConcept);
}

bool pmCellWriteDark(pmCell *cell, psFits *fits, psDB *db, bool blank)
{
    PS_ASSERT_PTR_NON_NULL(cell, false);
    PS_ASSERT_FITS_NON_NULL(fits, false);

    // Allow the usual pmFPAWrite functions to handle the heavy lifting for the images
    if (!pmCellWrite(cell, fits, db, blank)) {
        psError(PS_ERR_UNKNOWN, false, "Unable to write dark cell.");
        return false;
    }

    pmHDU *hdu = pmHDUFromCell(cell);
    if (blank || !hdu) {
        // No more to do
        return true;
    }

    psArray *ordinates = psMetadataLookupPtr(NULL, cell->analysis, PM_DARK_ANALYSIS_ORDINATES);
    if (!ordinates) {
        psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find dark ordinates.");
        return false;
    }
    bool mdok;                          // Status of MD lookup
    psString normConcept = psMetadataLookupPtr(&mdok, cell->analysis, PM_DARK_ANALYSIS_NORM); // Normalisation

    return pmDarkWrite(fits, hdu->header, ordinates, normConcept);
}

bool pmDarkWrite(psFits *fits, psMetadata *header, const psArray *ordinates, const char *normConcept)
{
    PS_ASSERT_FITS_NON_NULL(fits, false);
    PS_ASSERT_FITS_WRITABLE(fits, false);
    PS_ASSERT_ARRAY_NON_NULL(ordinates, false);

    // Only write table once per FITS file
    bool write = false;             // Write table?
    if (psFitsGetSize(fits) <= 1) {
        write = true;
    } else {
        psMetadata *headers = psFitsReadHeaderSet(NULL, fits); // FITS headers
        if (!psMetadataLookup(headers, PM_DARK_FITS_EXTNAME)) {
            write = true;
        }
        psFree(headers);
    }
    if (!write) {
        return true;
    }

    if (!psMemIncrRefCounter((psMetadata*)header)) {
        header = psMetadataAlloc();
    }
    psMetadataAddStr(header, PS_LIST_TAIL, PM_DARK_HEADER_NORM, PS_META_REPLACE,
                     "Dark normalisation concept", normConcept);

    if (ordinates->n > 0) {
        // Format ordinates into FITS table
        int numOrdinates = ordinates->n;// Number of ordinates
        psArray *table = psArrayAlloc(numOrdinates); // FITS table, constructed from ordinates
        for (int i = 0; i < ordinates->n; i++) {
            pmDarkOrdinate *ord = ordinates->data[i]; // Ordinate of interest
            psMetadata *row = psMetadataAlloc(); // FITS table row
            psMetadataAddStr(row, PS_LIST_TAIL, PM_DARK_FITS_NAME, 0, "Concept name", ord->name);
            psMetadataAddS32(row, PS_LIST_TAIL, PM_DARK_FITS_ORDER, 0, "Polynomial order", ord->order);
            psMetadataAddBool(row, PS_LIST_TAIL, PM_DARK_FITS_SCALE, 0, "Scale values?", ord->scale);
            psMetadataAddF32(row, PS_LIST_TAIL, PM_DARK_FITS_MIN, 0, "Minimum value", ord->min);
            psMetadataAddF32(row, PS_LIST_TAIL, PM_DARK_FITS_MAX, 0, "Maximum value", ord->max);
            table->data[i] = row;
        }

        if (!psFitsWriteTable(fits, header, table, PM_DARK_FITS_EXTNAME)) {
            psError(PS_ERR_IO, false, "Unable to write dark ordinates.");
            psFree(table);
            psFree(header);
            return false;
        }
        psFree(table);
    } else {
        // No ordinates to write to a table, so write to a blank header.
        if (!psFitsWriteBlank(fits, header, PM_DARK_FITS_EXTNAME)) {
            psError(PS_ERR_IO, false, "Unable to write dark header.");
            psFree(header);
            return false;
        }
    }
    psFree(header);

    return true;
}


bool pmFPAReadDark(pmFPA *fpa, psFits *fits, psDB *db)
{
    PS_ASSERT_PTR_NON_NULL(fpa, false);
    PS_ASSERT_FITS_NON_NULL(fits, false);

    if (!pmFPARead(fpa, fits, db)) {
        psError(PS_ERR_UNKNOWN, false, "Unable to read dark FPA.");
        return false;
    }

    psString normConcept = NULL;        // Normalisation concept
    psArray *ordinates = pmDarkRead(&normConcept, fits); // Dark ordinates
    if (!ordinates) {
        psError(PS_ERR_IO, false, "Unable to read dark ordinates.");
        return false;
    }

    psArray *chips = fpa->chips;        // Component chips
    for (int i = 0; i < chips->n; i++) {
        pmChip *chip = chips->data[i];  // Chip of interest
        psArray *cells = chip->cells;   // Component cells
        for (int j = 0; j < cells->n; j++) {
            pmCell *cell = cells->data[j]; // Cell of interest
            psMetadataAddPtr(cell->analysis, PS_LIST_TAIL, PM_DARK_ANALYSIS_ORDINATES,
                             PS_DATA_ARRAY | PS_META_REPLACE, "Dark ordinates", ordinates);
            psMetadataAddStr(cell->analysis, PS_LIST_TAIL, PM_DARK_ANALYSIS_NORM, PS_META_REPLACE,
                             "Dark normalisation", normConcept);
        }
    }
    psFree(ordinates);
    psFree(normConcept);

    return true;
}

bool pmChipReadDark(pmChip *chip, psFits *fits, psDB *db)
{
    PS_ASSERT_PTR_NON_NULL(chip, false);
    PS_ASSERT_FITS_NON_NULL(fits, false);

    if (!pmChipRead(chip, fits, db)) {
        psError(PS_ERR_UNKNOWN, false, "Unable to read dark chip.");
        return false;
    }

    psString normConcept = NULL;        // Normalisation concept
    psArray *ordinates = pmDarkRead(&normConcept, fits); // Dark ordinates
    if (!ordinates) {
        psError(PS_ERR_IO, false, "Unable to read dark ordinates.");
        return false;
    }

    psArray *cells = chip->cells;       // Component cells
    for (int i = 0; i < cells->n; i++) {
        pmCell *cell = cells->data[i];  // Cell of interest
        psMetadataAddPtr(cell->analysis, PS_LIST_TAIL, PM_DARK_ANALYSIS_ORDINATES,
                         PS_DATA_ARRAY | PS_META_REPLACE, "Dark ordinates", ordinates);
        psMetadataAddStr(cell->analysis, PS_LIST_TAIL, PM_DARK_ANALYSIS_NORM, PS_META_REPLACE,
                         "Dark normalisation", normConcept);
    }
    psFree(ordinates);
    psFree(normConcept);

    return true;
}


bool pmCellReadDark(pmCell *cell, psFits *fits, psDB *db)
{
    PS_ASSERT_PTR_NON_NULL(cell, false);
    PS_ASSERT_FITS_NON_NULL(fits, false);

    // Allow the usual pmFPARead functions to handle the heavy lifting for the images
    if (!pmCellRead(cell, fits, db)) {
        psError(PS_ERR_UNKNOWN, false, "Unable to read dark cell.");
        return false;
    }

    psString normConcept = NULL;        // Normalisation concept
    psArray *ordinates = pmDarkRead(&normConcept, fits); // Dark ordinates
    if (!ordinates) {
        psError(PS_ERR_IO, false, "Unable to read dark ordinates.");
        return false;
    }
    psMetadataAddPtr(cell->analysis, PS_LIST_TAIL, PM_DARK_ANALYSIS_ORDINATES,
                     PS_DATA_ARRAY | PS_META_REPLACE, "Dark ordinates", ordinates);
    psMetadataAddStr(cell->analysis, PS_LIST_TAIL, PM_DARK_ANALYSIS_NORM, PS_META_REPLACE,
                     "Dark normalisation", normConcept);
    psFree(ordinates);
    psFree(normConcept);

    return true;
}

psArray *pmDarkRead(psString *normConcept, psFits *fits)
{
    PS_ASSERT_PTR_NON_NULL(normConcept, NULL);
    PS_ASSERT_PTR_NULL(*normConcept, NULL);
    PS_ASSERT_FITS_NON_NULL(fits, NULL);

    if (!psFitsMoveExtName(fits, PM_DARK_FITS_EXTNAME)) {
        psError(PS_ERR_IO, false, "Unable to find extension with dark ordinates table (%s).",
                PM_DARK_FITS_EXTNAME);
        return false;
    }

    psMetadata *header = psFitsReadHeader(NULL, fits); // Header
    bool mdok;                          // Status of MD lookup
    *normConcept = psMemIncrRefCounter(psMetadataLookupStr(&mdok, header, PM_DARK_HEADER_NORM));

    psArray *ordinates = NULL;          // Dark ordinates to return

    psFitsType type = psFitsGetExtType(fits); // Type of FITS extension
    switch (type) {
      case PS_FITS_TYPE_IMAGE: {
          // Check that it's of zero size; otherwise we might have some conflict
          int numCols = psMetadataLookupS32(&mdok, header, "NAXIS1");
          int numRows = psMetadataLookupS32(&mdok, header, "NAXIS2");
          if (numCols != 0 || numRows != 0) {
              psError(PS_ERR_UNKNOWN, true, "Extension %s is not a DARK table.", PM_DARK_FITS_EXTNAME);
              psFree(header);
              return NULL;
          }
          // No ordinates fit --- only a constant term
          ordinates = psArrayAlloc(0);
          break;
      }
      case PS_FITS_TYPE_BINARY_TABLE:
      case PS_FITS_TYPE_ASCII_TABLE: {
          psArray *table = psFitsReadTable(fits); // FITS Table with ordinates
          int numOrdinates = table->n;        // Number of ordinates
          ordinates = psArrayAlloc(numOrdinates);

          for (int i = 0; i < numOrdinates; i++) {
              psMetadata *row = table->data[i]; // Row of interest
              const char *name = psMetadataLookupStr(NULL, row, PM_DARK_FITS_NAME); // Concept name
              int order = psMetadataLookupS32(NULL, row, PM_DARK_FITS_ORDER); // Polynomial order
              if (!name || order <= 0) {
                  psError(PS_ERR_UNKNOWN, false, "Bad value reading dark ordinates table.");
                  psFree(table);
                  psFree(ordinates);
                  return false;
              }
              pmDarkOrdinate *ord = pmDarkOrdinateAlloc(name, order); // Ordinate data
              ord->scale = psMetadataLookupBool(NULL, row, PM_DARK_FITS_SCALE);
              ord->min = psMetadataLookupF32(NULL, row, PM_DARK_FITS_MIN);
              ord->max = psMetadataLookupF32(NULL, row, PM_DARK_FITS_MAX);

              ordinates->data[i] = ord;
          }
          psFree(table);
          break;
      }
      default:
        psError(PS_ERR_UNKNOWN, true, "Unrecognised FITS extension type.");
        return NULL;
    }
    psFree(header);

    return ordinates;
}

