Index: trunk/psModules/src/detrend/pmDark.c
===================================================================
--- trunk/psModules/src/detrend/pmDark.c	(revision 16823)
+++ trunk/psModules/src/detrend/pmDark.c	(revision 16823)
@@ -0,0 +1,611 @@
+#include <stdio.h>
+#include <pslib.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)
+                           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;
+    }
+
+    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, 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);
+    for (int i = 0; i < numInputs; i++) {
+        values->data[i] = psVectorAlloc(numOrdinates, PS_TYPE_F32);
+    }
+    psVector *roMask = psVectorAlloc(numInputs, PS_TYPE_U8); // Mask for bad readouts
+    psVectorInit(roMask, 0);
+    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++) {
+            pmReadout *ro = inputs->data[j]; // Readout of interest
+            psVector *val = values->data[j]; // Value vector for readout
+
+            float value = NAN;          // Value of ordinate
+            if (!ordinateLookup(&value, ord->name, ro)) {
+                roMask->data.U8[j] = 0xff;
+                val->data.F32[i] = NAN;
+                numBadInputs++;
+                continue;
+            }
+            if (ord->scale) {
+                value = 2.0 * (value - ord->min) / (ord->max - ord->min) - 1.0;
+            }
+            val->data.F32[i] = value;
+        }
+    }
+    numInputs -= numBadInputs;
+
+    if (psTraceGetLevel("psModules.detrend") > 9) {
+        for (int i = 0; i < numInputs; i++) {
+            psVector *val = values->data[i];
+            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) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Insufficient inputs (%d) to fit polynomial terms (%d).",
+                numInputs, numTerms);
+        psFree(values);
+        psFree(roMask);
+        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.");
+        return false;
+    }
+    if (output->readouts->n != numTerms) {
+        output->readouts = psArrayRealloc(output->readouts, numTerms);
+    }
+    int outRow0, outCol0;               // 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);
+        }
+    }
+
+    // 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 (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];
+            }
+        }
+    }
+
+    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);
+
+    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;
+    }
+    PS_ASSERT_ARRAYS_SIZE_EQUAL(dark->readouts, ordinates, 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, ord->name, readout)) {
+            psError(PS_ERR_UNKNOWN, true, "Unable to find value for %s", ord->name);
+            psFree(values);
+            return false;
+        }
+        values->data.F32[i] = value;
+    }
+
+    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
+            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
+    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 *new = psMetadataLookupPtr(&mdok, cell->analysis, PM_DARK_ANALYSIS_ORDINATES);
+            if (!mdok) {
+                continue;
+            }
+            if (ordinates) {
+                if (new != ordinates) {
+                    psError(PS_ERR_UNKNOWN, true, "Dark ordinates differ across cells.");
+                    return false;
+                }
+            } else {
+                ordinates = new;
+            }
+        }
+    }
+
+    if (!ordinates) {
+        psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find dark ordinates.");
+        return false;
+    }
+
+    return pmDarkWrite(fits, hdu->header, ordinates);
+}
+
+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
+    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 *new = psMetadataLookupPtr(&mdok, cell->analysis, PM_DARK_ANALYSIS_ORDINATES);
+        if (!mdok) {
+            continue;
+        }
+        if (ordinates) {
+            if (new != ordinates) {
+                psError(PS_ERR_UNKNOWN, true, "Dark ordinates differ across cells.");
+                return false;
+            }
+        } else {
+            ordinates = new;
+        }
+    }
+
+    if (!ordinates) {
+        psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find dark ordinates.");
+        return false;
+    }
+
+    return pmDarkWrite(fits, hdu->header, ordinates);
+}
+
+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;
+    }
+
+    return pmDarkWrite(fits, hdu->header, ordinates);
+}
+
+bool pmDarkWrite(psFits *fits, const psMetadata *header, const psArray *ordinates)
+{
+    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;
+    }
+
+    // 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);
+        return false;
+    }
+    psFree(table);
+
+    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;
+    }
+
+    psArray *ordinates = pmDarkRead(fits);
+    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);
+        }
+    }
+    psFree(ordinates);
+
+    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;
+    }
+
+    psArray *ordinates = pmDarkRead(fits);
+    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);
+    }
+    psFree(ordinates);
+
+    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;
+    }
+
+    psArray *ordinates = pmDarkRead(fits);
+    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);
+    psFree(ordinates);
+
+    return true;
+}
+
+psArray *pmDarkRead(psFits *fits)
+{
+    PS_ASSERT_FITS_NON_NULL(fits, false);
+
+    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;
+    }
+
+    psArray *table = psFitsReadTable(fits); // FITS Table with ordinates
+    int numOrdinates = table->n;        // Number of ordinates
+    psArray *ordinates = psArrayAlloc(numOrdinates); // Parsed ordinates
+
+    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);
+
+    return ordinates;
+}
+