#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 "pmDetrendThreads.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
                           bool *inRange, // is value within min : max range?
                           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);

    *inRange = true;

    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);
            *inRange = 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;
}

// this creates and saves: values, roMask, norm, orders, counts, sigma, and saves the on output->analysis
bool pmDarkCombinePrepare(pmCell *output, const psArray *inputs, psArray *ordinates, const char *normConcept)
{
    psArray *values = psArrayAlloc(inputs->n);
    psVector *roMask = psVectorAlloc(inputs->n, PS_TYPE_U8); // Mask for bad readouts
    psVector *norm = normConcept ? psVectorAlloc(inputs->n, PS_TYPE_F32) : NULL; // Normalizations for each
    psVector *orders = psVectorAlloc(ordinates->n, PS_TYPE_U8); // Orders for each concept

    psVectorInit(roMask, 0);

    bool inRange = false;
    int numBadInputs = 0;               // Number of bad inputs

    // build the 'norm' vector and the 'values' vectors, count the number of bad inputs
    for (int i = 0; i < inputs->n; i++) {
        values->data[i] = psVectorAlloc(ordinates->n, PS_TYPE_F32);
        if (!norm) continue;

        pmReadout *readout = inputs->data[i]; // Readout of interest
        float normValue;            // Normalisation value
        if (!ordinateLookup(&normValue, &inRange, normConcept, false, NAN, NAN, readout)) {
            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;
    }

    // build the 'orders' vector and set the array of 'values'
    for (int i = 0; i < ordinates->n; 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);
            psFree(norm);
            return false;
        }
        orders->data.U8[i] = ord->order;

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

            pmReadout *readout = inputs->data[j]; // Readout of interest
            float value = NAN;          // Value of ordinate
            if (!ordinateLookup(&value, &inRange, ord->name, ord->scale, ord->min, ord->max, readout)) {
                roMask->data.U8[j] = 0xff;
                val->data.F32[i] = NAN;
                numBadInputs++;
                continue;
            }
            if (!inRange) {
                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 < inputs->n; i++) {
            psVector *val = values->data[i];
            (void) val; // avoid unused variable message when tracing is compiled out
            for (int j = 0; j < ordinates->n; j++) {
                psTrace("psModules.detrend", 9, "Image %d, ordinate %d: %f\n", i, j, val->data.F32[j]);
            }
        }
    }

    int numTerms = 1;                   // Number of terms in polynomial
    for (int i = 0; i < orders->n; i++) {
        numTerms += orders->data.U8[i];
    }

    if (numTerms > inputs->n - numBadInputs) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Insufficient inputs (%ld) to fit polynomial terms (%d).",
                inputs->n - numBadInputs, numTerms);
        psFree(values);
        psFree(roMask);
        psFree(norm);
        return false;
    }

    // determine the output image size based on the input images
    int row0, col0, numCols, numRows;
    if (!pmReadoutStackSetOutputSize(&col0, &row0, &numCols, &numRows, inputs)) {
        psError(PS_ERR_UNKNOWN, false, "problem setting output readout size.");
        return false;
    }

    // the output is potentially a cube (depending on the dimensionality of the fit)
    if (output->readouts->n != numTerms) {
        output->readouts = psArrayRealloc(output->readouts, numTerms);
    }

    // generate the required output images based on the specified sizes
    for (int i = 0; i < numTerms; i++) {
        pmReadout *readout = output->readouts->data[i]; // Readout to update
        if (!readout) {
            readout = output->readouts->data[i] = pmReadoutAlloc(output);
        }

        pmReadoutStackDefineOutput(readout, col0, row0, numCols, numRows, false, false, 0);
        psTrace("psModules.imcombine", 7, "Output minimum: %d,%d\n", col0, row0);
    }

    // these calls allocate and save the requested images on the output analysis metadata
    psImage *counts = pmReadoutSetAnalysisImage(output->readouts->data[0], PM_READOUT_STACK_ANALYSIS_COUNT,
                                                numCols, numRows, PS_TYPE_U16, 0);
    if (!counts) {
        return false;
    }
    psImage *sigma = pmReadoutSetAnalysisImage(output->readouts->data[0], PM_READOUT_STACK_ANALYSIS_SIGMA,
                                               numCols, numRows, PS_TYPE_F32, NAN);
    if (!sigma) {
        return false;
    }

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

    psMetadataAddPtr(output->analysis, PS_LIST_TAIL, "DARK.VALUES",  PS_DATA_ARRAY  | PS_META_REPLACE,
                     "Dark values", values);
    psMetadataAddPtr(output->analysis, PS_LIST_TAIL, "DARK.RO.MASK", PS_DATA_VECTOR | PS_META_REPLACE,
                     "Dark Readout Mask", roMask);
    psMetadataAddPtr(output->analysis, PS_LIST_TAIL, "DARK.NORM",    PS_DATA_VECTOR | PS_META_REPLACE,
                     "Dark norm", norm);
    psMetadataAddPtr(output->analysis, PS_LIST_TAIL, "DARK.ORDERS",  PS_DATA_VECTOR | PS_META_REPLACE,
                     "Dark orders", orders);

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

    psFree(norm);
    psFree(roMask);
    psFree(orders);
    psFree(values);

    return true;
}

// do the combine work for this portion of the output (range is set by input data)
bool pmDarkCombine(pmCell *output, const psArray *inputs, int iter, float rej, psMaskType maskVal)
{
    PS_ASSERT_PTR_NON_NULL(output, false);
    PS_ASSERT_PTR_NON_NULL(output->readouts, false);
    PS_ASSERT_INT_NONNEGATIVE(output->readouts->n, false);
    PS_ASSERT_ARRAY_NON_NULL(inputs, false);
    PS_ASSERT_INT_NONNEGATIVE(iter, false);
    PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, false);

    pthread_t id = pthread_self();
    char name[64];
    sprintf (name, "%x", (unsigned int) id);
    psTimerStart (name);

    bool mdok = false;

    // retrieve the required parameter vectors
    psArray *values  = psMetadataLookupPtr(&mdok, output->analysis, "DARK.VALUES");
    psAssert(values, "values not supplied");
    psVector *roMask = psMetadataLookupPtr(&mdok, output->analysis, "DARK.RO.MASK");
    psAssert(roMask, "roMask not supplied");
    psVector *orders = psMetadataLookupPtr(&mdok, output->analysis, "DARK.ORDERS");
    psAssert(orders, "orders not supplied");

    psPolynomialMD *poly = psPolynomialMDAlloc(orders); // Polynomial for fitting

    // retrieve the norm vector, if supplied
    psVector *norm       = psMetadataLookupPtr(&mdok, output->analysis, "DARK.NORM");

    // retrieve the 'counts' and 'sigma' images
    psImage *counts = pmReadoutGetAnalysisImage(output->readouts->data[0], PM_READOUT_STACK_ANALYSIS_COUNT);
    if (!counts) {
        return false;
    }
    psImage *sigma = pmReadoutGetAnalysisImage(output->readouts->data[0], PM_READOUT_STACK_ANALYSIS_SIGMA);
    if (!sigma) {
        return false;
    }

    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.");
        return false;
    }

    pmReadout *outReadout = output->readouts->data[0];

    // Iterate over pixels, fitting polynomial
    psVector *pixels = psVectorAlloc(inputs->n, PS_TYPE_F32); // Stack of pixels
    psVector *mask   = psVectorAlloc(inputs->n, PS_TYPE_MASK); // Mask for stack
    for (int i = minInputRows; i < maxInputRows; i++) {
        int yOut = i - outReadout->row0; // 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 - outReadout->col0; // 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 < poly->coeff->n; 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(pixels);
    psFree(mask);

    return true;
}

bool pmDarkApplyScan_Threaded(psThreadJob *job)
{
    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);

    pmReadout *readout     = job->args->data[0]; // Readout to correct
    pmCell *dark           = job->args->data[1]; // Dark to apply
    const psVector *orders = job->args->data[2]; // Polynomial orders for each ordinate
    const psVector *values = job->args->data[3]; // Values for each ordinate

    psMaskType bad = PS_SCALAR_VALUE(job->args->data[4], U8); // Mask value to give bad pixels
    bool doNorm    = PS_SCALAR_VALUE(job->args->data[5], U8); // Normalise values?
    float norm     = PS_SCALAR_VALUE(job->args->data[6], F32); // Value by which to normalise
    int rowStart   = PS_SCALAR_VALUE(job->args->data[7], S32); // Starting row for scan
    int rowStop    = PS_SCALAR_VALUE(job->args->data[8], S32); // Stopping row for scan

    return pmDarkApplyScan(readout, dark, orders, values, bad, doNorm, norm, rowStart, rowStop);
}

bool pmDarkApplyScan(pmReadout *readout, const pmCell *dark, const psVector *orders, const psVector *values,
                     psMaskType bad, bool doNorm, float norm, int rowStart, int rowStop)
{
    int numCols = readout->image->numCols;
    int numTerms = dark->readouts->n;   // Number of polynomial terms

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

    for (int y = rowStart; y < rowStop; 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);

    return true;
}

bool pmDarkApply(pmReadout *readout, 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
    bool inRange = false;

    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, &inRange, 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, &inRange, 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;
    }

    // thread here by scan

    bool threaded = true;
    int scanRows = pmDetrendGetScanRows();
    if (scanRows == 0) {
        threaded = false;
        scanRows = readout->image->numRows;
    }

    for (int rowStart = 0; rowStart < readout->image->numRows; rowStart += scanRows) {
        int rowStop = PS_MIN(rowStart + scanRows, readout->image->numRows);

        if (threaded) {
            psThreadJob *job = psThreadJobAlloc("PSMODULES_DETREND_DARK");
            psArrayAdd(job->args, 1, readout);
            psArrayAdd(job->args, 1, dark);
            psArrayAdd(job->args, 1, orders);
            psArrayAdd(job->args, 1, values);
            PS_ARRAY_ADD_SCALAR(job->args, bad, PS_TYPE_MASK);
            PS_ARRAY_ADD_SCALAR(job->args, doNorm, PS_TYPE_U8);
            PS_ARRAY_ADD_SCALAR(job->args, norm, PS_TYPE_F32);
            PS_ARRAY_ADD_SCALAR(job->args, rowStart, PS_TYPE_S32);
            PS_ARRAY_ADD_SCALAR(job->args, rowStop, PS_TYPE_S32);

            if (!psThreadJobAddPending (job)) {
                psFree(job);
                psFree(orders);
                psFree(values);
                return false;
            }
            psFree(job);
        } else if (!pmDarkApplyScan(readout, dark, orders, values, bad, doNorm, norm, rowStart, rowStop)) {
            psError(PS_ERR_UNKNOWN, false, "Unable to apply dark.");
            psFree(orders);
            psFree(values);
            return false;
        }
    }

    if (threaded) {
        // wait here for the threaded jobs to finish
        if (!psThreadPoolWait(true)) {
            psError(PS_ERR_UNKNOWN, false, "Unable to apply dark.");
            psFree(orders);
            psFree(values);
            return false;
        }
    }

    psFree(orders);
    psFree(values);

    return true;
}

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

    if (!pmFPAWrite(fpa, fits, config, 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, pmConfig *config, bool blank, bool recurse)
{
    PS_ASSERT_PTR_NON_NULL(chip, false);
    PS_ASSERT_FITS_NON_NULL(fits, false);

    if (!pmChipWrite(chip, fits, config, 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, pmConfig *config, 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, config, 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, pmConfig *config)
{
    PS_ASSERT_PTR_NON_NULL(fpa, false);
    PS_ASSERT_FITS_NON_NULL(fits, false);

    if (!pmFPARead(fpa, fits, config)) {
        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, pmConfig *config)
{
    PS_ASSERT_PTR_NON_NULL(chip, false);
    PS_ASSERT_FITS_NON_NULL(fits, false);

    if (!pmChipRead(chip, fits, config)) {
        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, pmConfig *config)
{
    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, config)) {
        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;
}

