Index: /branches/eam_rel9_p0/psModules/src/astrom/Makefile.am
===================================================================
--- /branches/eam_rel9_p0/psModules/src/astrom/Makefile.am	(revision 6061)
+++ /branches/eam_rel9_p0/psModules/src/astrom/Makefile.am	(revision 6062)
@@ -4,5 +4,6 @@
 libpsmoduleastrom_la_LDFLAGS  = -release $(PACKAGE_VERSION)
 libpsmoduleastrom_la_SOURCES  = \
-	pmAstrometry.c \
+	pmFPA.c \
+	pmFPAAstrometry.c \
 	pmAstrometryObjects.c \
 	pmFPAConstruct.c \
@@ -16,10 +17,9 @@
 	pmConceptsStandard.c \
 	pmChipMosaic.c
-###	pmFPAConceptsGet.c
-###	pmFPAConceptsSet.c
 
 psmoduleincludedir = $(includedir)
 psmoduleinclude_HEADERS = \
-	pmAstrometry.h \
+	pmFPA.h \
+	pmFPAAstrometry.h \
 	pmAstrometryObjects.h \
 	pmFPAConstruct.h \
@@ -32,5 +32,3 @@
 	pmConceptsWrite.h \
 	pmConceptsStandard.h \
-	pmChipMosaic.h
-###	pmFPAConceptsGet.h
-###	pmFPAConceptsSet.h
+	pmChipMosaic.h 
Index: /branches/eam_rel9_p0/psModules/src/astrom/pmAstrometry.c
===================================================================
--- /branches/eam_rel9_p0/psModules/src/astrom/pmAstrometry.c	(revision 6061)
+++ /branches/eam_rel9_p0/psModules/src/astrom/pmAstrometry.c	(revision 6062)
@@ -13,6 +13,6 @@
 * XXX: Should we implement non-linear cell->chip transforms?
 *
-*  @version $Revision: 1.11.2.1.2.1 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2006-01-14 03:10:19 $
+*  @version $Revision: 1.11.2.1.2.2 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2006-01-20 02:36:41 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -63,5 +63,10 @@
         psFree(readout->weight);
         psFree(readout->analysis);
+        #if 0
+
         psFree(readout->parent);
+        #endif
+
+        readout->parent = NULL;
     }
 }
@@ -80,4 +85,6 @@
         // in order to avoid memory reference counter problems.
         //
+        #if 0
+
         for (psS32 i = 0 ; i < cell->readouts->n ; i++) {
             pmReadout *tmpReadout = (pmReadout *) cell->readouts->data[i];
@@ -87,7 +94,12 @@
             }
         }
+        psFree(cell->parent);
+        #endif
+
+        cell->parent = NULL;
+
         psFree(cell->readouts);
-        psFree(cell->parent);
         psFree(cell->hdu);
+
     }
 }
@@ -104,4 +116,6 @@
         // in order to avoid memory reference counter problems.
         //
+        #if 0
+
         for (psS32 i = 0 ; i < chip->cells->n ; i++) {
             pmCell *tmpCell = (pmCell *) chip->cells->data[i];
@@ -111,6 +125,9 @@
             }
         }
+        psFree(chip->parent);
+        #endif
+
+        chip->parent = NULL;
         psFree(chip->cells);
-        psFree(chip->parent);
         psFree(chip->hdu);
     }
@@ -131,4 +148,6 @@
         // in order to avoid memory reference counter problems.
         //
+        #if 0
+
         for (psS32 i = 0 ; i < fpa->chips->n ; i++) {
             pmChip *tmpChip = (pmChip *) fpa->chips->data[i];
@@ -138,4 +157,5 @@
             }
         }
+        #endif
         psFree(fpa->chips);
         psFree(fpa->hdu);
@@ -199,5 +219,5 @@
     tmpCell->analysis = psMetadataAlloc();
     tmpCell->readouts = psArrayAlloc(0);
-    tmpCell->parent = psMemIncrRefCounter(chip);
+    tmpCell->parent = chip;
     if (chip != NULL) {
         chip->cells = psArrayAdd(chip->cells, 1, (psPtr) tmpCell);
@@ -232,5 +252,5 @@
     tmpChip->analysis = psMetadataAlloc();
     tmpChip->cells = psArrayAlloc(0);
-    tmpChip->parent = psMemIncrRefCounter(fpa);
+    tmpChip->parent = fpa;
     if (fpa != NULL) {
         fpa->chips = psArrayAdd(fpa->chips, 1, (psPtr) tmpChip);
@@ -338,488 +358,4 @@
 }
 
-
-
-/*****************************************************************************/
-/* FUNCTION IMPLEMENTATION - PUBLIC                                          */
-/*****************************************************************************/
-
-pmCell* pmCellInFPA(
-    const psPlane* fpaCoord,
-    const pmFPA* FPA)
-{
-    PS_ASSERT_PTR_NON_NULL(fpaCoord, NULL);
-    PS_ASSERT_PTR_NON_NULL(FPA, NULL);
-
-    pmChip* tmpChip = NULL;
-    psPlane chipCoord;
-    pmCell* outCell = NULL;
-
-    // Determine which chip contains the fpaCoords.
-    tmpChip = pmChipInFPA(fpaCoord, FPA);
-    if (tmpChip == NULL) {
-        return(NULL);
-    }
-
-    // Convert to those chip coordinates.
-    psPlane *rc = pmCoordFPAToChip(&chipCoord, fpaCoord, tmpChip);
-    if (rc == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine Chip coords.\n");
-        return(NULL);
-    }
-
-    // Determine which cell contains those chip coordinates.
-    outCell = pmCellInChip(&chipCoord, tmpChip);
-    if (outCell == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine the cell.\n");
-        return(NULL);
-    }
-
-    return (outCell);
-}
-
-pmChip* pmChipInFPA(
-    const psPlane* fpaCoord,
-    const pmFPA* FPA)
-{
-    PS_ASSERT_PTR_NON_NULL(fpaCoord, NULL);
-    PS_ASSERT_PTR_NON_NULL(FPA, NULL);
-    PS_ASSERT_PTR_NON_NULL(FPA->chips, NULL);
-
-    psArray* chips = FPA->chips;
-    psS32 nChips = chips->n;
-    psPlane chipCoord;
-    pmCell *tmpCell = NULL;
-
-    //
-    // Loop through every chip in this FPA.  Convert the original FPA
-    // coordinates to chip coordinates for that chip.  Then, determine if any
-    // cells in that chip contain those chip coordinates.
-    // XXX: Depending on the number of chips, and their topology, there may be
-    // a much more efficient way of doing this.
-    //
-    for (psS32 i = 0; i < nChips; i++) {
-        pmChip* tmpChip = chips->data[i];
-        PS_ASSERT_PTR_NON_NULL(tmpChip, NULL);
-        PS_ASSERT_PTR_NON_NULL(tmpChip->fromFPA, NULL);
-
-        psPlaneTransformApply(&chipCoord, tmpChip->fromFPA, fpaCoord);
-
-        tmpCell = pmCellInChip(&chipCoord, tmpChip);
-        if (tmpCell != NULL) {
-            return(tmpChip);
-        }
-    }
-
-    psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine the chip.\n");
-    return (NULL);
-}
-
-
-pmCell* pmCellInChip(
-    const psPlane* chipCoord,
-    const pmChip* chip)
-{
-    PS_ASSERT_PTR_NON_NULL(chipCoord, NULL);
-    PS_ASSERT_PTR_NON_NULL(chip, NULL);
-
-    psPlane cellCoord;
-    psArray* cells;
-
-    cells = chip->cells;
-    if (cells == NULL) {
-        return NULL;
-    }
-
-    //
-    // We loop over each cell in the chip.  We transform the chipCoord into
-    // a cellCoord for that cell and determine if that cellCoord is valid.
-    // If so, then we return that cell.
-    // XXX: Depending on the number of cells, and their topology, there may be
-    // a much more efficient way of doing this.
-    //
-    for (psS32 i = 0; i < cells->n; i++) {
-        pmCell* tmpCell = (pmCell* ) cells->data[i];
-        PS_ASSERT_PTR_NON_NULL(tmpCell, NULL);
-
-        psPlaneTransform *chipToCell = NULL;
-        if (true ==  p_psIsProjectionLinear(tmpCell->toChip)) {
-            chipToCell = p_psPlaneTransformLinearInvert(tmpCell->toChip);
-        } else {
-            psLogMsg(__func__, PS_LOG_WARN, "WARNING: non-linear cell->chip transforms are not yet implemented.\n");
-            //chipToCell = psPlaneTransformInvert(NULL, tmpCell->toChip, NULL, -1);
-            chipToCell = NULL;
-        }
-        if (chipToCell == NULL) {
-            psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not invert the Cell->toChip transform.\n");
-            return(NULL);
-        }
-        psArray* readouts = tmpCell->readouts;
-
-        if (readouts != NULL) {
-            for (psS32 j = 0; j < readouts->n; j++) {
-                pmReadout* tmpReadout = readouts->data[j];
-                PS_ASSERT_READOUT_NON_NULL(tmpReadout, NULL);
-
-                psPlaneTransformApply(&cellCoord,
-                                      chipToCell,
-                                      chipCoord);
-
-                if (checkValidImageCoords(cellCoord.x,
-                                          cellCoord.y,
-                                          tmpReadout->image)) {
-                    psFree(chipToCell);
-                    return (tmpCell);
-                }
-            }
-        }
-        psFree(chipToCell);
-    }
-
-    //psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine the cell.\n");
-    return (NULL);
-}
-
-
-psPlane* pmCoordCellToFPA(
-    psPlane* fpaCoord,
-    const psPlane* cellCoord,
-    const pmCell* cell)
-{
-    PS_ASSERT_PTR_NON_NULL(cellCoord, NULL);
-    PS_ASSERT_PTR_NON_NULL(cell, NULL);
-
-    psPlane *rc = psPlaneTransformApply(fpaCoord, cell->toFPA, cellCoord);
-    if (rc == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform cell coords to FPA coords.\n");
-    }
-    return(rc);
-}
-
-
-psPlane* pmCoordChipToFPA(
-    psPlane* outCoord,
-    const psPlane* inCoord,
-    const pmChip* chip)
-{
-    PS_ASSERT_PTR_NON_NULL(inCoord, NULL);
-    PS_ASSERT_PTR_NON_NULL(chip, NULL);
-
-    psPlane *rc = psPlaneTransformApply(outCoord, chip->toFPA, inCoord);
-    if (rc == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform chip coords to FPA coords.\n");
-    }
-    return(rc);
-}
-
-
-psPlane* pmCoordFPAToChip(
-    psPlane* chipCoord,
-    const psPlane* fpaCoord,
-    const pmChip* chip)
-{
-    PS_ASSERT_PTR_NON_NULL(fpaCoord, NULL);
-    PS_ASSERT_PTR_NON_NULL(chip, NULL);
-    PS_ASSERT_PTR_NON_NULL(chip->fromFPA, NULL);
-
-    psPlane *rc = psPlaneTransformApply(chipCoord, chip->fromFPA, fpaCoord);
-    if (rc == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform FPA coords to Chip coords.\n");
-    }
-    return(rc);
-}
-
-psPlane* pmCoordCellToChip(
-    psPlane* outCoord,
-    const psPlane* inCoord,
-    const pmCell* cell)
-{
-    PS_ASSERT_PTR_NON_NULL(inCoord, NULL);
-    PS_ASSERT_PTR_NON_NULL(cell, NULL);
-
-    psPlane *rc = psPlaneTransformApply(outCoord, cell->toChip, inCoord);
-    if (rc == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform Cell coords to Chip coords.\n");
-    }
-    return(rc);
-}
-
-psPlane* pmCoordChipToCell(
-    psPlane* cellCoord,
-    const psPlane* chipCoord,
-    const pmCell* cell)
-{
-    PS_ASSERT_PTR_NON_NULL(chipCoord, NULL);
-    PS_ASSERT_PTR_NON_NULL(cell, NULL);
-    PS_ASSERT_PTR_NON_NULL(cell->parent, NULL);
-
-    pmCell *tmpCell = pmCellInChip(chipCoord, cell->parent);
-    if (tmpCell == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine the proper cell.\n");
-        return(NULL);
-    }
-
-    psPlaneTransform *tmpChipToCell = NULL;
-    PS_ASSERT_PTR_NON_NULL(tmpCell->toChip, NULL);
-    if (true ==  p_psIsProjectionLinear(tmpCell->toChip)) {
-        tmpChipToCell = p_psPlaneTransformLinearInvert(tmpCell->toChip);
-    } else {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: non-linear cell->chip transforms are not yet implemented.\n");
-        // XXX: tmpChipToCell = psPlaneTransformInvert(NULL, tmpCell->toChip, NULL, -1);
-        tmpChipToCell = NULL;
-    }
-    if (tmpChipToCell == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not invert the Cell->toChip transform.\n");
-        return(NULL);
-    }
-
-    psPlane *rc = psPlaneTransformApply(cellCoord, tmpChipToCell, chipCoord);
-    if (rc == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform Chip coords to Cell coords.\n");
-    }
-    psFree(tmpChipToCell);
-    return(rc);
-}
-
-psPlane* pmCoordFPAToTP(
-    psPlane* outCoord,
-    const psPlane* inCoord,
-    double color,
-    double magnitude,
-    const pmFPA* fpa)
-{
-    PS_ASSERT_PTR_NON_NULL(inCoord, NULL);
-    PS_ASSERT_PTR_NON_NULL(fpa, NULL);
-
-    psPlane *rc = psPlaneDistortApply(outCoord, fpa->toTangentPlane, inCoord, color, magnitude);
-    if (rc == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform FPA coords to tangent plane coords.\n");
-    }
-    return(rc);
-}
-
-psPlane* pmCoordTPToFPA(
-    psPlane* fpaCoord,
-    const psPlane* tpCoord,
-    double color,
-    double magnitude,
-    const pmFPA* fpa)
-{
-    PS_ASSERT_PTR_NON_NULL(tpCoord, NULL);
-    PS_ASSERT_PTR_NON_NULL(fpa, NULL);
-    PS_ASSERT_PTR_NON_NULL(fpa->fromTangentPlane, NULL);
-
-    psPlane *rc = psPlaneDistortApply(fpaCoord, fpa->fromTangentPlane, tpCoord, color, magnitude);
-    if (rc == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform tangent plane coords to FPA coords.\n");
-    }
-    return(rc);
-}
-
-
-/*****************************************************************************
-XXXDeproject(outSphere, coord, projection): This private routine is a wrapper
-for p_psDeproject().  The reason: p_psDeproject() and p_psProject() combined
-do not seem to produce the original coordinates when they even though they
-should.  XXXDeproject() simply negates the ->r and ->d members of the output
-psSphere if the input ->y is larger than 0.0.  I don't know why it works.
- 
-I'm guessing the p_psProject() and p_psDeproject() functions have bugs.
- 
-XXX: It appears that p_psProject() and p_psDeproject() have been fixed.
-Remove this.
- *****************************************************************************/
-psSphere* XXXDeproject(
-    psSphere *outSphere,
-    const psPlane* coord,
-    const psProjection* projection)
-{
-    psSphere *rc = p_psDeproject(outSphere, coord, projection);
-
-    if (coord->y >= 0.0) {
-        rc->d = -rc->d;
-        rc->r = -rc->r;
-    }
-
-    return(rc);
-}
-
-/*****************************************************************************
-  *****************************************************************************/
-psSphere* pmCoordTPToSky(
-    psSphere* outSphere,
-    const psPlane* tpCoord,
-    const psProjection *projection)
-{
-    PS_ASSERT_PTR_NON_NULL(tpCoord, NULL);
-    PS_ASSERT_PTR_NON_NULL(projection, NULL);
-
-    //    psSphere *rc = XXXDeproject(outSphere, tpCoord, projection);
-    psSphere *rc = p_psDeproject(outSphere, tpCoord, projection);
-    if (rc == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform tangent plane coords to sky coords.\n");
-    }
-    return(rc);
-}
-
-/*****************************************************************************
- *****************************************************************************/
-psPlane* pmCoordSkyToTP(
-    psPlane* tpCoord,
-    const psSphere* in,
-    const psProjection *projection)
-{
-    PS_ASSERT_PTR_NON_NULL(in, NULL);
-    PS_ASSERT_PTR_NON_NULL(projection, NULL);
-
-    psPlane *rc = p_psProject(tpCoord, in, projection);
-    if (rc == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform sky to tangent plane coords.\n");
-    }
-    return(rc);
-}
-
-/*****************************************************************************
- *****************************************************************************/
-psSphere* pmCoordCellToSky(
-    psSphere* skyCoord,
-    const psPlane* cellCoord,
-    double color,
-    double magnitude,
-    const pmCell* cell)
-{
-    PS_ASSERT_PTR_NON_NULL(cellCoord, NULL);
-    PS_ASSERT_PTR_NON_NULL(cell, NULL);
-    PS_ASSERT_PTR_NON_NULL(cell->toFPA, NULL);
-    PS_ASSERT_PTR_NON_NULL(cell->parent, NULL);
-    PS_ASSERT_PTR_NON_NULL(cell->parent->parent, NULL);
-    PS_ASSERT_PTR_NON_NULL(cell->parent->parent->toTangentPlane, NULL);
-    PS_ASSERT_PTR_NON_NULL(cell->parent->parent->projection, NULL);
-    psPlane fpaCoord;
-    psPlane tpCoord;
-    psPlane *rc;
-    pmFPA* parFPA = (cell->parent)->parent;
-
-    // Convert the input cell coordinates to FPA coordinates.
-    rc = psPlaneTransformApply(&fpaCoord, cell->toFPA, cellCoord);
-    if (rc == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could transform cell coords to FPA coords.\n");
-        return(NULL);
-    }
-
-    // Convert the FPA coordinates to tangent plane Coordinates.
-    rc = psPlaneDistortApply(&tpCoord, parFPA->toTangentPlane, &fpaCoord, color, magnitude);
-    if (rc == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could transform FPA coords to tangent plane coords.\n");
-        return(NULL);
-    }
-
-    // Convert the tangent plane Coordinates to sky coordinates.
-    psSphere *rc2 = pmCoordTPToSky(skyCoord, &tpCoord, parFPA->projection);
-    if (rc2 == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform cell coords to sky coords.\n");
-    }
-
-    return(rc2);
-}
-
-/*****************************************************************************
- *****************************************************************************/
-psPlane* pmCoordSkyToCell(
-    psPlane* cellCoord,
-    const psSphere* skyCoord,
-    float color,
-    float magnitude,
-    const pmCell* cell)
-{
-    PS_ASSERT_PTR_NON_NULL(skyCoord, NULL);
-    PS_ASSERT_PTR_NON_NULL(cell, NULL);
-    PS_ASSERT_PTR_NON_NULL(cell->parent, NULL);
-    PS_ASSERT_PTR_NON_NULL(cell->parent->parent, NULL);
-    pmChip *parChip = cell->parent;
-    pmFPA *parFPA = parChip->parent;
-    psPlane tpCoord;
-    psPlane fpaCoord;
-    psPlane chipCoord;
-    psPlane *rc;
-
-    // Convert the skyCoords to tangent plane coords.
-    rc = pmCoordSkyToTP(&tpCoord, skyCoord, parFPA->projection);
-    if (rc == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine tangent plane coords.\n");
-        return(NULL);
-    }
-
-    // Convert the tangent plane coords to FPA coords.
-    rc = pmCoordTPToFPA(&fpaCoord, &tpCoord, color, magnitude, parFPA);
-    if (rc == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine FPA coords.\n");
-        return(NULL);
-    }
-
-    // Convert the FPA coords to chip coords.
-    rc = pmCoordFPAToChip(&chipCoord, &fpaCoord, parChip);
-    if (rc == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine chip coords.\n");
-        return(NULL);
-    }
-
-    // Convert the chip coords to cell coords.
-    rc = pmCoordChipToCell(cellCoord, &chipCoord, cell);
-    if (rc == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine cell coords.\n");
-        return(NULL);
-    }
-
-    return (cellCoord);
-}
-
-/*****************************************************************************
- *****************************************************************************/
-psSphere* pmCoordCellToSkyQuick(
-    psSphere* outSphere,
-    const psPlane* cellCoord,
-    const pmCell* cell)
-{
-    PS_ASSERT_PTR_NON_NULL(cellCoord, NULL);
-    PS_ASSERT_PTR_NON_NULL(cell, NULL);
-    PS_ASSERT_PTR_NON_NULL(cell->toSky, NULL);
-    psPlane outPlane;
-    psPlane *rc;
-    rc = psPlaneTransformApply(&outPlane, cell->toSky, cellCoord);
-    if (rc == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could transform cell coords to sky coords.\n");
-        return(NULL);
-    }
-
-    psSphere *out = outSphere;
-    if (out == NULL) {
-        out = psSphereAlloc();
-    }
-    out->r = outPlane.y;
-    out->d = outPlane.x;
-
-    return(out);
-}
-
-/*****************************************************************************
- *****************************************************************************/
-psPlane* pmCoordSkyToCellQuick(
-    psPlane* cellCoord,
-    const psSphere* skyCoord,
-    const pmCell* cell)
-{
-    PS_ASSERT_PTR_NON_NULL(skyCoord, NULL);
-    PS_ASSERT_PTR_NON_NULL(cell, NULL);
-    PS_ASSERT_PTR_NON_NULL(cell->toSky, NULL);
-    psPlane skyPlane;
-    skyPlane.y = skyCoord->r;
-    skyPlane.x = skyCoord->d;
-
-    psPlane *rc = psPlaneTransformApply(cellCoord, cell->toSky, &skyPlane);
-    if (rc == NULL) {
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform sky to cell coords.\n");
-    }
-    return(cellCoord);
-}
 
 
@@ -925,2 +461,74 @@
     return(numChips);
 }
+
+
+bool pmCellSetWeights(pmCell *cell // Cell for which to set weights
+                     )
+{
+    float gain = psMetadataLookupF32(NULL, cell->concepts, "CELL.GAIN"); // Cell gain
+    float readnoise = psMetadataLookupF32(NULL, cell->concepts, "CELL.READNOISE"); // Cell read noise
+    psRegion *trimsec = psMetadataLookupPtr(NULL, cell->concepts, "CELL.TRIMSEC"); // Trim section
+
+    p_pmHDU *hdu = cell->hdu;           // The data unit, containing the weight and mask originals
+    if (!hdu) {
+        pmChip *chip = cell->parent;    // The parent chip
+        if (chip->hdu) {
+            hdu = chip->hdu;
+        } else {
+            pmFPA *fpa = chip->parent;  // The parent FPA
+            hdu = fpa->hdu;
+            if (!hdu) {
+                psError(PS_ERR_UNKNOWN, true, "Unable to find the HDU in the hierarchy!\n");
+                return false;
+            }
+        }
+    }
+
+    psArray *pixels = hdu->images;      // Array of images
+    psArray *weights = hdu->weights;    // Array of weight images
+    psArray *masks = hdu->masks;        // Array of mask images
+    // Generate the weights and masks if required
+    if (! weights) {
+        weights = psArrayAlloc(pixels->n);
+        for (int i = 0; i < pixels->n; i++) {
+            psImage *image = pixels->data[i];
+            weights->data[i] = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32);
+            psImageInit(weights->data[i], 0.0);
+        }
+        hdu->weights = weights;
+    }
+    if (! masks) {
+        masks = psArrayAlloc(pixels->n);
+        for (int i = 0; i < pixels->n; i++) {
+            psImage *image = pixels->data[i];
+            masks->data[i] = psImageAlloc(image->numCols, image->numRows, PS_TYPE_U8);
+            psImageInit(masks->data[i], 0);
+        }
+        hdu->masks = masks;
+    }
+
+    // Set the pixels
+    psArray *readouts = cell->readouts; // Array of readouts
+    for (int i = 0; i < readouts->n; i++) {
+        pmReadout *readout = readouts->data[i]; // The readout of interest
+
+        if (! readout->weight) {
+            readout->weight = weights->data[i];
+        }
+        if (! readout->mask) {
+            readout->mask = masks->data[i];
+        }
+
+        // Mask is already set to 0
+
+        // Set weight image to the variance = g*f + rn^2
+        psImage *image = psImageSubset(readout->image, *trimsec); // The pixels
+        psImage *weight = psImageSubset(readout->weight, *trimsec); // The weight map
+        psBinaryOp(weight, image, "*", psScalarAlloc(gain, PS_TYPE_F32));
+        psBinaryOp(weight, weight, "+", psScalarAlloc(readnoise*readnoise, PS_TYPE_F32));
+    }
+
+    return true;
+}
+
+
Index: /branches/eam_rel9_p0/psModules/src/astrom/pmAstrometryObjects.c
===================================================================
--- /branches/eam_rel9_p0/psModules/src/astrom/pmAstrometryObjects.c	(revision 6061)
+++ /branches/eam_rel9_p0/psModules/src/astrom/pmAstrometryObjects.c	(revision 6062)
@@ -8,6 +8,6 @@
 *  @author EAM, IfA
 *
-*  @version $Revision: 1.1.10.1 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2006-01-15 18:19:38 $
+*  @version $Revision: 1.1.10.2 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2006-01-20 02:36:41 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -20,5 +20,5 @@
 #include <math.h>
 #include "pslib.h"
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 #include "pmAstrometryObjects.h"
 
Index: /branches/eam_rel9_p0/psModules/src/astrom/pmAstrometryObjects.h
===================================================================
--- /branches/eam_rel9_p0/psModules/src/astrom/pmAstrometryObjects.h	(revision 6061)
+++ /branches/eam_rel9_p0/psModules/src/astrom/pmAstrometryObjects.h	(revision 6062)
@@ -8,6 +8,6 @@
 *  @author EAM, IfA
 *
-*  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2005-12-05 20:49:30 $
+*  @version $Revision: 1.1.10.1 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2006-01-20 02:36:41 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -21,8 +21,8 @@
 # include <unistd.h>   // for unlink
 # include <pslib.h>
-# include <pmAstrometry.h>
-
-/*
- * 
+# include <pmFPA.h>
+
+/*
+ *
  * This structure specifies the coordinate of the detection in each of the
  * four necessary coordinate frames: pix defines the position in the psReadout
@@ -35,7 +35,7 @@
  * coordinates, while the reference detections will be projected to the other
  * frames from the sky coordinates.
- * 
+ *
  * XXX: There are more members here than in the SDRS.
- * 
+ *
  */
 typedef struct
@@ -55,9 +55,9 @@
 
 /*
- * 
+ *
  * The pmAstromMatch structure defines the cross-correlation between two
  * arrays. A single such data item specifies that item number pmAstromMatch.idx1
  * in the first list corresponds to pmAstromMatch.idx2 in the second list.
- * 
+ *
  */
 typedef struct
@@ -70,7 +70,7 @@
 
 /*
- * 
+ *
  * XXX: Not in SDRS.
- * 
+ *
  */
 typedef struct
@@ -88,5 +88,5 @@
 
 /*
- * 
+ *
  * If the two sets of coordinates are expected to agree very well (ie, the
  * current best-guess astrometric solution is quite close to the radius. The
@@ -97,5 +97,5 @@
  * ASTROM.MATCH.RADIUS). The output consists an array of pmAstromMatch values,
  * defined below.
- * 
+ *
  */
 psArray *pmAstromRadiusMatch(
@@ -108,11 +108,11 @@
 
 /*
- * 
+ *
  * This function accepts an array of pmAstromObj objects and rotates them by
  * the given angle about the given center coordinate pCenter,qCenter in the Focal
  * Plane Array coordinates.
- * 
+ *
  * XXX: This differs from the SDRS
- * 
+ *
  */
 /* SDRS
@@ -132,5 +132,5 @@
 
 /*
- * 
+ *
  * If the two sets of coordinates are not known to agree well, but the
  * relative scale and approximate relative rotation is known, then a much faster
@@ -147,7 +147,7 @@
  * allowing the procedure to scan over a range of rotations. We define the
  * following function to apply this matching algorithm:
- * 
+ *
  * XXX: In the SDRS, this function is a pointer.
- * 
+ *
  */
 pmAstromStats pmAstromGridMatch(
@@ -159,5 +159,5 @@
 
 /*
- * 
+ *
  * The result of a pmAstromGridMatch may be used to modify the astrometry
  * transformation information for a pmFPA image hierarchy structure. The result
@@ -167,7 +167,7 @@
  * the linear terms of the pmFPA.toTangentPlane transformation. These two
  * adjustments are made using the function:
- * 
+ *
  * XXX: This function name is different in the SDRS.
- * 
+ *
  */
 psPlaneTransform *pmAstromGridApply(
@@ -178,11 +178,11 @@
 
 /*
- * 
+ *
  * This function is identical to pmAstromGridMatch, but is valid for only a
  * single relative rotation. The input config information need not contain any of
  * the GRID.*.ANGLE entries (they will be ignored).
- * 
+ *
  * XXX: This function name is different in the SDRS.
- * 
+ *
  */
 /* in pmAstromGrid.c */
@@ -195,5 +195,5 @@
 
 /*
- * 
+ *
  * This function accepts the raw and reference source lists and the list of
  * matched entries. It uses the matched list to determine a polynomial
@@ -209,7 +209,7 @@
  * modifications to pmFPA.toTangentPlane incorporate the rotation component of
  * the linear terms and the higher-order terms of the polynomial fits.
- * 
+ *
  * XXX: No prototype code.
- * 
+ *
  */
 bool pmAstromFitFPA(
@@ -235,7 +235,7 @@
  *ASTROM.ORDER).  The result of this fit is a set of modifications of the
  *components of the pmChip.toFPA transformation.
- * 
+ *
  * XXX: No prototype code.
- * 
+ *
  */
 bool pmAstromFitChip(
@@ -249,5 +249,5 @@
 
 /*
- * 
+ *
  * The following function determines the position residual, in the tangent
  * plane, as a function of position in the focal plane, for a collection of raw
@@ -255,7 +255,7 @@
  * the bin size over which the gradient is measured (keyword: ASTROM.GRAD.BOX).
  * The function returns an array of pmAstromGradient structures, defined below.
- * 
+ *
  * XXX: No prototype code.
- * 
+ *
  */
 psArray pmAstromMeasureGradients(
@@ -269,9 +269,9 @@
 
 /*
- * 
+ *
  * The following data structure carries the information about the residual
  * gradient of source positions in the tangent plane (pmAstromObj.TP) as a
  * function of position in the focal plane (pmAstromObj.FP).
- * 
+ *
  */
 typedef struct
@@ -285,5 +285,5 @@
 
 /*
- * 
+ *
  * The gradient set measured above can be fitted with a pair of 2D
  * polynomials. The resulting fits can then be related back to the implied
@@ -292,7 +292,7 @@
  * supplied pmFPA structure. The configuration variable supplies the polynomial
  * order (keyword: ASTROM.DISTORT.ORDER).
- * 
+ *
  * XXX: No prototype code.
- * 
+ *
  */
 psArray pmAstromFitDistortion(
@@ -310,15 +310,15 @@
  ******************************************************************************/
 /*
- * 
- * 
- * 
- * 
- */
-
-
-/*
- * 
+ *
+ *
+ *
+ *
+ */
+
+
+/*
+ *
  * Allocates a pmAstromObj struct.
- * 
+ *
  */
 pmAstromObj *pmAstromObjAlloc (void);
@@ -327,7 +327,7 @@
 
 /*
- * 
+ *
  * Copies a pmAstromObj struct.
- * 
+ *
  */
 pmAstromObj *pmAstromObjCopy(
@@ -338,7 +338,7 @@
 
 /*
- * 
- * 
- * 
+ *
+ *
+ *
  */
 pmAstromMatch *pmAstromMatchAlloc(
@@ -351,7 +351,7 @@
 
 /*
- * 
- * 
- * 
+ *
+ *
+ *
  */
 psPlaneTransform *pmAstromMatchFit(
@@ -366,7 +366,7 @@
 
 /*
- * 
- * 
- * 
+ *
+ *
+ *
  */
 int pmAstromObjSortByFPX(
@@ -378,7 +378,7 @@
 
 /*
- * 
- * 
- * 
+ *
+ *
+ *
  */
 int pmAstromObjSortByMag(
Index: /branches/eam_rel9_p0/psModules/src/astrom/pmChipMosaic.c
===================================================================
--- /branches/eam_rel9_p0/psModules/src/astrom/pmChipMosaic.c	(revision 6061)
+++ /branches/eam_rel9_p0/psModules/src/astrom/pmChipMosaic.c	(revision 6062)
@@ -3,5 +3,5 @@
 
 #include "pslib.h"
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 #include "pmChipMosaic.h"
 
@@ -34,9 +34,12 @@
     // Get the maximum extent of the mosaic image
     int xMin = INT_MAX;
-    int xMax = 0;
+    int xMax = - INT_MAX;
     int yMin = INT_MAX;
-    int yMax = 0;
+    int yMax = - INT_MAX;
     for (int i = 0; i < source->n; i++) {
         psImage *image = source->data[i]; // The image of interest
+        if (!image) {
+            continue;
+        }
 
         assert(image->type.type == PS_TYPE_F32); // Only implemented for F32 images so far.
@@ -74,4 +77,7 @@
     for (int i = 0; i < source->n; i++) {
         psImage *image = source->data[i]; // The image of interest
+        if (!image) {
+            continue;
+        }
         if (xBinSource->data.S32[i] == xBinTarget && yBinSource->data.S32[i] == yBinTarget &&
                 xFlip->data.U8[i] == 0 && yFlip->data.U8[i] == 0) {
@@ -107,11 +113,13 @@
 
 // Mosaic a chip together into a single image
-static psImage *mosaicChip(pmChip *chip,// Chip to mosaic
-                           int xBinChip, int yBinChip // Binning of mosaic image in x and y
-                          )
+static bool mosaicChip(pmChip *chip,// Chip to mosaic
+                       int xBinChip, int yBinChip // Binning of mosaic image in x and y
+                      )
 {
 
     psArray *cells = chip->cells;       // The array of cells
     psArray *images = psArrayAlloc(cells->n); // Array of images that will be mosaicked
+    psArray *weights = psArrayAlloc(cells->n); // Array of weight images to be mosaicked
+    psArray *masks = psArrayAlloc(cells->n); // Array of mask images to be mosaicked
     psVector *x0 = psVectorAlloc(cells->n, PS_TYPE_S32); // Origin x coordinates
     psVector *y0 = psVectorAlloc(cells->n, PS_TYPE_S32); // Origin y coordinates
@@ -126,4 +134,5 @@
         x0->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.X0");
         y0->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.Y0");
+        psTrace(__func__, 9, "Cell %d: x0=%d y0=%d\n", i, x0->data.S32[i], y0->data.S32[i]);
         xBin->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN");
         yBin->data.S32[i] = psMetadataLookupS32(NULL, cell->concepts, "CELL.XBIN");
@@ -159,9 +168,14 @@
         }
         psImage *image = ((pmReadout*)readouts->data[0])->image; // The image to put into the mosaic
+        psImage *weight = ((pmReadout*)readouts->data[0])->weight;
+        psImage *mask = ((pmReadout*)readouts->data[0])->weight;
         images->data[i] = psImageSubset(image, *trimsec); // Trimmed image
-    }
-
+        weights->data[i] = weight ? psImageSubset(weight, *trimsec) : NULL;
+        masks->data[i] = mask ? psImageSubset(mask, *trimsec) : NULL;
+    }
     // Mosaic the images together and we're done
-    psImage *mosaic = p_pmImageMosaic(images, xFlip, yFlip, xBin, yBin, xBinChip, yBinChip, x0, y0);
+    psImage *image = p_pmImageMosaic(images, xFlip, yFlip, xBin, yBin, xBinChip, yBinChip, x0, y0);
+    psImage *weight = p_pmImageMosaic(weights, xFlip, yFlip, xBin, yBin, xBinChip, yBinChip, x0, y0);
+    psImage *mask = p_pmImageMosaic(masks, xFlip, yFlip, xBin, yBin, xBinChip, yBinChip, x0, y0);
 
     // Clean up
@@ -173,6 +187,31 @@
     psFree(yFlip);
     psFree(images);
-
-    return mosaic;
+    psFree(weights);
+    psFree(masks);
+
+    // Chop off all the component cells, and put in a new one
+
+    #if 1
+    // XXX For the sake of getting something working, I am not going to bother with sorting out where
+    // the double free is coming from.  I'm going to drop the pointers on the array and create a memory
+    // leak.  We can clean this up later, when we're not under as much pressure.
+    chip->cells = psArrayAlloc(0);
+    #else
+
+    psArrayElementsFree(chip->cells);
+    chip->cells = psArrayRealloc(chip->cells, 0);
+    chip->cells->n = 0;
+    #endif
+
+    pmCell *cell = pmCellAlloc(chip, NULL, __func__); // New cell
+    cell->exists = true;
+    cell->process = true;
+    pmReadout *readout = pmReadoutAlloc(cell); // New readout
+    readout->image = image;
+    readout->weight = weight;
+    readout->mask = mask;
+    //psFree(readout);
+
+    return true;
 }
 
@@ -193,13 +232,6 @@
 
         numChips++;
-        psImage *mosaic = mosaicChip(chip, xBinChip, yBinChip); // Mosaic of cells within the chip
-
-        // Chop off all the component cells, and put in a new one
-        psArrayElementsFree(chip->cells);
-        chip->cells = psArrayRealloc(chip->cells, 1);
-        pmCell *cell = pmCellAlloc(chip, NULL, __func__); // New cell
-        pmReadout *readout = pmReadoutAlloc(cell); // New readout
-        readout->image = mosaic;
-        psFree(readout);
+        mosaicChip(chip, xBinChip, yBinChip); // Mosaic of cells within the chip
+
     }
 
Index: /branches/eam_rel9_p0/psModules/src/astrom/pmChipMosaic.h
===================================================================
--- /branches/eam_rel9_p0/psModules/src/astrom/pmChipMosaic.h	(revision 6061)
+++ /branches/eam_rel9_p0/psModules/src/astrom/pmChipMosaic.h	(revision 6062)
@@ -3,5 +3,5 @@
 
 #include "pslib.h"
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 
 // Mosaic multiple images, with flips, binning and offsets
Index: /branches/eam_rel9_p0/psModules/src/astrom/pmConcepts.c
===================================================================
--- /branches/eam_rel9_p0/psModules/src/astrom/pmConcepts.c	(revision 6061)
+++ /branches/eam_rel9_p0/psModules/src/astrom/pmConcepts.c	(revision 6062)
@@ -99,10 +99,10 @@
 
 // Set all registered concepts to blank value for the specified level
-static bool conceptsBlank(psMetadata *specs,  // One of the concepts specifications
+static bool conceptsBlank(psMetadata **specs,  // One of the concepts specifications
                           psMetadata *target // Place to install the concepts
                          )
 {
     pmConceptsInit();
-    psMetadataIterator *specsIter = psMetadataIteratorAlloc(specs, PS_LIST_HEAD, NULL); // Iterator on specs
+    psMetadataIterator *specsIter = psMetadataIteratorAlloc(*specs, PS_LIST_HEAD, NULL); // Iterator on specs
     psMetadataItem *specItem = NULL;    // Item from the specs metadata
     while ((specItem = psMetadataGetAndIncrement(specsIter))) {
@@ -117,5 +117,5 @@
 
 // Read all registered concepts for the specified level
-static bool conceptsRead(psMetadata *specs, // One of the concepts specifications
+static bool conceptsRead(psMetadata **specs, // One of the concepts specifications
                          pmFPA *fpa,    // The FPA
                          pmChip *chip,  // The chip
@@ -126,5 +126,5 @@
 {
     pmConceptsInit();
-    psMetadataIterator *specsIter = psMetadataIteratorAlloc(specs, PS_LIST_HEAD, NULL); // Iterator on specs
+    psMetadataIterator *specsIter = psMetadataIteratorAlloc(*specs, PS_LIST_HEAD, NULL); // Iterator on specs
     psMetadataItem *specItem = NULL;    // Item from the specs metadata
     while ((specItem = psMetadataGetAndIncrement(specsIter))) {
@@ -147,5 +147,5 @@
 
 // Write all registered concepts for the specified level
-static bool conceptsWrite(psMetadata *specs, // One of the concepts specifications
+static bool conceptsWrite(psMetadata **specs, // One of the concepts specifications
                           pmFPA *fpa,   // The FPA
                           pmChip *chip, // The chip
@@ -160,5 +160,5 @@
     while ((item = psMetadataGetAndIncrement(iter))) {
         const char *name = item->name;  // Name of the concept
-        psMetadataItem *specItem = psMetadataLookup(specs, name); // Specification for the concept
+        psMetadataItem *specItem = psMetadataLookup(*specs, name); // Specification for the concept
         if (specItem) {
             pmConceptSpec *spec = specItem->data.V; // The specification
@@ -182,5 +182,6 @@
                        )
 {
-    return conceptsBlank(conceptsFPA, fpa->concepts);
+    psTrace("psModule.concepts", 5, "Blanking FPA concepts: %x %x\n", conceptsFPA, fpa->concepts);
+    return conceptsBlank(&conceptsFPA, fpa->concepts);
 }
 
@@ -190,5 +191,6 @@
                       )
 {
-    return conceptsRead(conceptsFPA, fpa, NULL, NULL, db, fpa->concepts);
+    psTrace("psModule.concepts", 5, "Reading FPA concepts: %x %x\n", conceptsFPA, fpa->concepts);
+    return conceptsRead(&conceptsFPA, fpa, NULL, NULL, db, fpa->concepts);
 }
 
@@ -198,5 +200,6 @@
                        )
 {
-    return conceptsWrite(conceptsFPA, fpa, NULL, NULL, db, fpa->concepts);
+    psTrace("psModule.concepts", 5, "Writing FPA concepts: %x %x\n", conceptsFPA, fpa->concepts);
+    return conceptsWrite(&conceptsFPA, fpa, NULL, NULL, db, fpa->concepts);
 }
 
@@ -205,5 +208,6 @@
                         )
 {
-    return conceptsBlank(conceptsChip, chip->concepts);
+    psTrace("psModule.concepts", 5, "Blanking chip concepts: %x %x\n", conceptsChip, chip->concepts);
+    return conceptsBlank(&conceptsChip, chip->concepts);
 }
 
@@ -213,6 +217,7 @@
                        )
 {
+    psTrace("psModule.concepts", 5, "Reading chip concepts: %x %x\n", conceptsChip, chip->concepts);
     pmFPA *fpa = chip->parent;          // FPA to which the chip belongs
-    return conceptsRead(conceptsChip, fpa, chip, NULL, db, chip->concepts);
+    return conceptsRead(&conceptsChip, fpa, chip, NULL, db, chip->concepts);
 }
 
@@ -222,6 +227,7 @@
                         )
 {
+    psTrace("psModule.concepts", 5, "Writing chip concepts: %x %x\n", conceptsChip, chip->concepts);
     pmFPA *fpa = chip->parent;          // FPA to which the chip belongs
-    return conceptsWrite(conceptsChip, fpa, chip, NULL, db, chip->concepts);
+    return conceptsWrite(&conceptsChip, fpa, chip, NULL, db, chip->concepts);
 }
 
@@ -230,5 +236,6 @@
                         )
 {
-    return conceptsBlank(conceptsCell, cell->concepts);
+    psTrace("psModule.concepts", 5, "Blanking cell concepts: %x %x\n", conceptsCell, cell->concepts);
+    return conceptsBlank(&conceptsCell, cell->concepts);
 }
 
@@ -238,7 +245,8 @@
                        )
 {
+    psTrace("psModule.concepts", 5, "Writing cell concepts: %x %x\n", conceptsCell, cell->concepts);
     pmChip *chip = cell->parent;        // Chip to which the cell belongs
     pmFPA *fpa = chip->parent;          // FPA to which the chip belongs
-    return conceptsRead(conceptsCell, fpa, chip, cell, db, cell->concepts);
+    return conceptsRead(&conceptsCell, fpa, chip, cell, db, cell->concepts);
 }
 
@@ -248,7 +256,8 @@
                         )
 {
+    psTrace("psModule.concepts", 5, "Writing cell concepts: %x %x\n", conceptsCell, cell->concepts);
     pmChip *chip = cell->parent;        // Chip to which the cell belongs
     pmFPA *fpa = chip->parent;          // FPA to which the chip belongs
-    return conceptsWrite(conceptsCell, fpa, chip, cell, db, cell->concepts);
+    return conceptsWrite(&conceptsCell, fpa, chip, cell, db, cell->concepts);
 }
 
@@ -297,5 +306,5 @@
         conceptsChip = psMetadataAlloc();
         init = true;
-        // XXX Insert here the native Chip concepts
+        // There are no standard concepts at the chip level to be installed
     }
     if (! conceptsCell) {
@@ -397,4 +406,14 @@
         }
 
+        // CELL.X0
+        pmConceptRegister(psMetadataItemAllocS32("CELL.X0", "Position of (0,0) on the chip", 0),
+                          (pmConceptReadFunc)pmConceptRead_CELL_X0,
+                          (pmConceptWriteFunc)pmConceptWrite_CELL_X0, PM_CONCEPT_LEVEL_CELL);
+
+        // CELL.Y0
+        pmConceptRegister(psMetadataItemAllocS32("CELL.Y0", "Position of (0,0) on the chip", 0),
+                          (pmConceptReadFunc)pmConceptRead_CELL_Y0,
+                          (pmConceptWriteFunc)pmConceptWrite_CELL_Y0, PM_CONCEPT_LEVEL_CELL);
+
     }
 
Index: /branches/eam_rel9_p0/psModules/src/astrom/pmConcepts.h
===================================================================
--- /branches/eam_rel9_p0/psModules/src/astrom/pmConcepts.h	(revision 6061)
+++ /branches/eam_rel9_p0/psModules/src/astrom/pmConcepts.h	(revision 6062)
@@ -4,5 +4,5 @@
 #include "pslib.h"
 
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 
 // Function to call to read a concept
Index: /branches/eam_rel9_p0/psModules/src/astrom/pmConceptsRead.c
===================================================================
--- /branches/eam_rel9_p0/psModules/src/astrom/pmConceptsRead.c	(revision 6061)
+++ /branches/eam_rel9_p0/psModules/src/astrom/pmConceptsRead.c	(revision 6062)
@@ -3,5 +3,5 @@
 #include "pslib.h"
 
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 #include "pmConceptsRead.h"
 #include "psAdditionals.h"
Index: /branches/eam_rel9_p0/psModules/src/astrom/pmConceptsRead.h
===================================================================
--- /branches/eam_rel9_p0/psModules/src/astrom/pmConceptsRead.h	(revision 6061)
+++ /branches/eam_rel9_p0/psModules/src/astrom/pmConceptsRead.h	(revision 6062)
@@ -2,5 +2,5 @@
 #define PM_CONCEPTS_READ_H
 
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 
 psMetadataItem *pmConceptReadFromCamera(pmCell *cell, // The cell
Index: /branches/eam_rel9_p0/psModules/src/astrom/pmConceptsStandard.c
===================================================================
--- /branches/eam_rel9_p0/psModules/src/astrom/pmConceptsStandard.c	(revision 6061)
+++ /branches/eam_rel9_p0/psModules/src/astrom/pmConceptsStandard.c	(revision 6062)
@@ -5,5 +5,5 @@
 #include "pmConceptsRead.h"
 #include "pmConceptsWrite.h"
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 #include "pmConceptsStandard.h"
 #include "psAdditionals.h"
@@ -987,6 +987,7 @@
 {
     psMetadataItem *x0item = psMetadataLookup(cell->concepts, "CELL.X0");
-    x0item->data.S32 -= fortranCorr(fpa, "CELL.X0");
-    return pmConceptWriteItem(fpa, chip, cell, db, x0item);
+    psMetadataItem *newItem = psMetadataItemAllocS32("CELL.X0", "Position of (0,0) on the chip",
+                              x0item->data.S32 - fortranCorr(fpa, "CELL.X0"));
+    return pmConceptWriteItem(fpa, chip, cell, db, newItem);
 }
 
@@ -994,6 +995,7 @@
 {
     psMetadataItem *y0item = psMetadataLookup(cell->concepts, "CELL.Y0");
-    y0item->data.S32 -= fortranCorr(fpa, "CELL.Y0");
-    return pmConceptWriteItem(fpa, chip, cell, db, y0item);
-}
-
+    psMetadataItem *newItem = psMetadataItemAllocS32("CELL.Y0", "Position of (0,0) on the chip",
+                              y0item->data.S32 - fortranCorr(fpa, "CELL.Y0"));
+    return pmConceptWriteItem(fpa, chip, cell, db, newItem);
+}
+
Index: /branches/eam_rel9_p0/psModules/src/astrom/pmConceptsStandard.h
===================================================================
--- /branches/eam_rel9_p0/psModules/src/astrom/pmConceptsStandard.h	(revision 6061)
+++ /branches/eam_rel9_p0/psModules/src/astrom/pmConceptsStandard.h	(revision 6062)
@@ -3,5 +3,5 @@
 
 #include "pslib.h"
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 
 psMetadataItem *pmConceptRead_FPA_RA(pmFPA *fpa, pmChip *chip, pmCell *cell, psDB *db);
Index: /branches/eam_rel9_p0/psModules/src/astrom/pmConceptsWrite.c
===================================================================
--- /branches/eam_rel9_p0/psModules/src/astrom/pmConceptsWrite.c	(revision 6061)
+++ /branches/eam_rel9_p0/psModules/src/astrom/pmConceptsWrite.c	(revision 6062)
@@ -3,5 +3,5 @@
 #include "pslib.h"
 
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 #include "pmConceptsRead.h"
 #include "psAdditionals.h"
Index: /branches/eam_rel9_p0/psModules/src/astrom/pmConceptsWrite.h
===================================================================
--- /branches/eam_rel9_p0/psModules/src/astrom/pmConceptsWrite.h	(revision 6061)
+++ /branches/eam_rel9_p0/psModules/src/astrom/pmConceptsWrite.h	(revision 6062)
@@ -3,5 +3,5 @@
 
 #include "pslib.h"
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 
 // Well, not really "write", but check to make sure it's there and matches
Index: /branches/eam_rel9_p0/psModules/src/astrom/pmFPA.c
===================================================================
--- /branches/eam_rel9_p0/psModules/src/astrom/pmFPA.c	(revision 6062)
+++ /branches/eam_rel9_p0/psModules/src/astrom/pmFPA.c	(revision 6062)
@@ -0,0 +1,509 @@
+/** @file  pmFPA.c
+*
+*  @brief This file defines the basic types for the FPA hierarchy
+*
+*  @ingroup AstroImage
+*
+*  @author GLG, MHPCC
+*
+* XXX: We should review the extent of the warning messages on these functions
+* when the transformations are not successful.
+*
+* XXX: Should we implement non-linear cell->chip transforms?
+*
+*  @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2006-01-20 02:36:41 $
+*
+*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
+*/
+
+/******************************************************************************/
+/*  INCLUDE FILES                                                             */
+/******************************************************************************/
+#include <string.h>
+#include <math.h>
+#include "pslib.h"
+
+#include "pmFPA.h"
+#include "pmConcepts.h"
+
+/******************************************************************************
+ *****************************************************************************/
+#define PS_FREE_HIERARCHY 1
+static void readoutFree(pmReadout *readout)
+{
+    if (readout != NULL) {
+        psFree(readout->image);
+        psFree(readout->mask);
+        psFree(readout->weight);
+        psFree(readout->analysis);
+        #if 0
+
+        psFree(readout->parent);
+        #endif
+
+        readout->parent = NULL;
+    }
+}
+
+static void cellFree(pmCell *cell)
+{
+    if (cell != NULL) {
+        psFree(cell->toChip);
+        psFree(cell->toFPA);
+        psFree(cell->toSky);
+        psFree(cell->concepts);
+        psFree(cell->analysis);
+        psFree(cell->camera);
+        //
+        // Set the parent to NULL in all cell->readouts before psFree(cell->readouts)
+        // in order to avoid memory reference counter problems.
+        //
+        #if 0
+
+        for (psS32 i = 0 ; i < cell->readouts->n ; i++) {
+            pmReadout *tmpReadout = (pmReadout *) cell->readouts->data[i];
+            tmpReadout->parent = NULL;
+            if (PS_FREE_HIERARCHY == 1) {
+                psFree(tmpReadout);
+            }
+        }
+        psFree(cell->parent);
+        #endif
+
+        cell->parent = NULL;
+
+        psFree(cell->readouts);
+        psFree(cell->hdu);
+
+    }
+}
+
+static void chipFree(pmChip* chip)
+{
+    if (chip != NULL) {
+        psFree(chip->toFPA);
+        psFree(chip->fromFPA);
+        psFree(chip->concepts);
+        psFree(chip->analysis);
+        //
+        // Set the parent to NULL in all chip->cells before psFree(chip->cells)
+        // in order to avoid memory reference counter problems.
+        //
+        #if 0
+
+        for (psS32 i = 0 ; i < chip->cells->n ; i++) {
+            pmCell *tmpCell = (pmCell *) chip->cells->data[i];
+            tmpCell->parent = NULL;
+            if (PS_FREE_HIERARCHY == 1) {
+                psFree(tmpCell);
+            }
+        }
+        psFree(chip->parent);
+        #endif
+
+        chip->parent = NULL;
+        psFree(chip->cells);
+        psFree(chip->hdu);
+    }
+}
+
+
+static void FPAFree(pmFPA *fpa)
+{
+    if (fpa != NULL) {
+        psFree(fpa->fromTangentPlane);
+        psFree(fpa->toTangentPlane);
+        psFree(fpa->projection);
+        psFree(fpa->concepts);
+        psFree(fpa->analysis);
+        psFree(fpa->camera);
+        //
+        // Set the parent to NULL in all fpa->chips before psFree(fpa->chips)
+        // in order to avoid memory reference counter problems.
+        //
+        #if 0
+
+        for (psS32 i = 0 ; i < fpa->chips->n ; i++) {
+            pmChip *tmpChip = (pmChip *) fpa->chips->data[i];
+            tmpChip->parent = NULL;
+            if (PS_FREE_HIERARCHY == 1) {
+                psFree(tmpChip);
+            }
+        }
+        #endif
+        psFree(fpa->chips);
+        psFree(fpa->hdu);
+        psFree(fpa->phu);
+    }
+}
+
+void p_pmHDUFree(p_pmHDU *hdu)
+{
+    if (hdu) {
+        psFree(hdu->extname);
+        psFree(hdu->header);
+        psFree(hdu->images);
+        psFree(hdu->masks);
+        psFree(hdu->weights);
+    }
+}
+
+// XXX: Verify these default values for row0, col0, rowBins, colBins
+// PAP: These values may disappear in the future in favour of values in parent->concepts?
+pmReadout *pmReadoutAlloc(pmCell *cell)
+{
+    pmReadout *tmpReadout = (pmReadout *) psAlloc(sizeof(pmReadout));
+
+    tmpReadout->col0 = 0;
+    tmpReadout->row0 = 0;
+    tmpReadout->colBins = 0;
+    tmpReadout->rowBins = 0;
+    tmpReadout->image = NULL;
+    tmpReadout->mask = NULL;
+    tmpReadout->weight = NULL;
+    tmpReadout->analysis = psMetadataAlloc();
+    tmpReadout->parent = cell;
+    if (cell != NULL) {
+        cell->readouts = psArrayAdd(cell->readouts, 1, (psPtr) tmpReadout);
+    }
+    psMemSetDeallocator(tmpReadout, (psFreeFunc) readoutFree);
+    return(tmpReadout);
+}
+
+// XXX: Verify these default values for row0, col0.
+// PAP: These values may disappear in the future in favour of values in the "concepts"?
+pmCell *pmCellAlloc(
+    pmChip *chip,
+    psMetadata *cameraData,
+    const char *name)
+{
+    pmCell *tmpCell = (pmCell *) psAlloc(sizeof(pmCell));
+
+    tmpCell->col0 = 0;
+    tmpCell->row0 = 0;
+    tmpCell->toChip = NULL;
+    tmpCell->toFPA = NULL;
+    tmpCell->toSky = NULL;
+    tmpCell->concepts = psMetadataAlloc();
+    psBool rc = psMetadataAddStr(tmpCell->concepts, PS_LIST_HEAD, "CELL.NAME", 0, NULL, name);
+    if (rc == false) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not add CELL.NAME to metadata.\n");
+    }
+    tmpCell->camera = psMemIncrRefCounter(cameraData);
+    tmpCell->analysis = psMetadataAlloc();
+    tmpCell->readouts = psArrayAlloc(0);
+    tmpCell->parent = chip;
+    if (chip != NULL) {
+        chip->cells = psArrayAdd(chip->cells, 1, (psPtr) tmpCell);
+    }
+    tmpCell->process = true;            // All cells are processed by default
+    tmpCell->exists = false;            // Not yet read in
+    tmpCell->hdu = NULL;
+
+    pmConceptsBlankCell(tmpCell);
+
+    psMemSetDeallocator(tmpCell, (psFreeFunc) cellFree);
+    return(tmpCell);
+}
+
+// XXX: Verify these default values for row0, col0.
+// PAP: row0, col0 may disappear in the future in favour of storing values in the "concepts".
+pmChip *pmChipAlloc(
+    pmFPA *fpa,
+    const char *name)
+{
+    pmChip *tmpChip = (pmChip *) psAlloc(sizeof(pmChip));
+
+    tmpChip->col0 = 0;
+    tmpChip->row0 = 0;
+    tmpChip->toFPA = NULL;
+    tmpChip->fromFPA = NULL;
+    tmpChip->concepts = psMetadataAlloc();
+    psBool rc = psMetadataAddStr(tmpChip->concepts, PS_LIST_HEAD, "CHIP.NAME", 0, NULL, name);
+    if (rc == false) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not add CHIP.NAME to metadata.\n");
+    }
+    tmpChip->analysis = psMetadataAlloc();
+    tmpChip->cells = psArrayAlloc(0);
+    tmpChip->parent = fpa;
+    if (fpa != NULL) {
+        fpa->chips = psArrayAdd(fpa->chips, 1, (psPtr) tmpChip);
+    }
+    tmpChip->process = true;            // Work on all chips, by default
+    tmpChip->exists = false;            // Not read in yet
+    tmpChip->hdu = NULL;
+
+    pmConceptsBlankChip(tmpChip);
+
+    psMemSetDeallocator(tmpChip, (psFreeFunc) chipFree);
+    return(tmpChip);
+}
+
+pmFPA *pmFPAAlloc(const psMetadata *camera)
+{
+    pmFPA *tmpFPA = (pmFPA *) psAlloc(sizeof(pmFPA));
+
+    tmpFPA->fromTangentPlane = NULL;
+    tmpFPA->toTangentPlane = NULL;
+    tmpFPA->projection = NULL;
+    tmpFPA->concepts = psMetadataAlloc();
+    tmpFPA->analysis = NULL;
+    tmpFPA->camera = psMemIncrRefCounter((psPtr)camera);
+    tmpFPA->chips = psArrayAlloc(0);
+    tmpFPA->hdu = NULL;
+    tmpFPA->phu = NULL;
+
+    pmConceptsBlankFPA(tmpFPA);
+
+    psMemSetDeallocator(tmpFPA, (psFreeFunc) FPAFree);
+    return(tmpFPA);
+}
+
+p_pmHDU *p_pmHDUAlloc(const char *extname)
+{
+    p_pmHDU *hdu = psAlloc(sizeof(p_pmHDU));
+    psMemSetDeallocator(hdu, (psFreeFunc)p_pmHDUFree);
+
+    hdu->extname = psStringCopy(extname);
+    hdu->header = NULL;
+    hdu->images = NULL;
+    hdu->masks = NULL;
+    hdu->weights = NULL;
+
+    return hdu;
+}
+
+static psBool cellCheckParents(pmCell *cell)
+{
+    if (cell == NULL) {
+        return(true);
+    }
+    psBool flag = true;
+
+    for (psS32 i = 0 ; i < cell->readouts->n ; i++) {
+        pmReadout *tmpReadout = (pmReadout *) cell->readouts->data[i];
+        PS_ASSERT_PTR_NON_NULL(tmpReadout, false);
+        if (tmpReadout->parent != cell) {
+            tmpReadout->parent = cell;
+            flag = false;
+        }
+    }
+    return(flag);
+}
+
+static psBool chipCheckParents(pmChip *chip)
+{
+    if (chip == NULL) {
+        return(true);
+    }
+    psBool flag = true;
+
+    for (psS32 i = 0 ; i < chip->cells->n ; i++) {
+        pmCell *tmpCell = (pmCell *) chip->cells->data[i];
+        PS_ASSERT_PTR_NON_NULL(tmpCell, false);
+        if (tmpCell->parent != chip) {
+            tmpCell->parent = chip;
+            flag = false;
+        }
+
+        flag &= cellCheckParents(tmpCell);
+    }
+    return(flag);
+}
+
+psBool pmFPACheckParents(pmFPA *fpa)
+{
+    if (fpa == NULL) {
+        return(true);
+    }
+    psBool flag = true;
+
+    for (psS32 i = 0 ; i < fpa->chips->n ; i++) {
+        pmChip *tmpChip = (pmChip *) fpa->chips->data[i];
+        PS_ASSERT_PTR_NON_NULL(tmpChip, false);
+        if (tmpChip->parent != fpa) {
+            tmpChip->parent = fpa;
+            flag = false;
+        }
+
+        flag &= chipCheckParents(tmpChip);
+    }
+    return(flag);
+}
+
+
+
+/*****************************************************************************
+ *****************************************************************************/
+
+// Set cells within a chip to be processed or not
+static bool setCellsProcess(const pmChip *chip, // Chip of interest
+                            bool process  // Process this chip?
+                           )
+{
+    PS_ASSERT_PTR_NON_NULL(chip, false);
+
+    psArray *cells = chip->cells;       // Component cells
+    if (! cells) {
+        return false;
+    }
+    for (int i = 0; i < cells->n; i++) {
+        pmCell *tmpCell = cells->data[i]; // Cell of interest
+        if (tmpCell) {
+            tmpCell->process = process;
+        }
+    }
+
+    return true;
+}
+
+/*****************************************************************************
+ *****************************************************************************/
+bool pmFPASelectChip(
+    pmFPA *fpa,
+    int chipNum)
+{
+    PS_ASSERT_PTR_NON_NULL(fpa, false);
+
+    psArray *chips = fpa->chips;        // Component chips
+    if ((chips == NULL) || (chipNum >= chips->n)) {
+        return(false);
+    }
+
+    for (int i = 0 ; i < chips->n ; i++) {
+        pmChip *tmpChip = (pmChip *) chips->data[i];
+        if (tmpChip == NULL) {
+            continue;
+        }
+        if (i == chipNum) {
+            tmpChip->process = true;
+            setCellsProcess(tmpChip, true);
+        } else {
+            tmpChip->process = false;
+            setCellsProcess(tmpChip, false);
+        }
+
+    }
+
+    return true;
+}
+
+
+/*****************************************************************************
+XXX: The SDRS is ambiguous on a few things:
+    Whether or not the other chips should be set process=true. [PAP: No]
+    Should we return the number of chip process=true before or after they're set, [PAP: After]
+ *****************************************************************************/
+/**
+ *
+ * pmFPAExcludeChip shall set process to false only for the specified chip
+ * number (chipNum). In the event that the specified chip number does not exist
+ * within the fpa, the function shall generate a warning, and perform no action.
+ * The function shall return the number of chips within the fpa that have process
+ * set to true.
+ *
+ */
+int pmFPAExcludeChip(
+    pmFPA *fpa,
+    int chipNum)
+{
+    PS_ASSERT_PTR_NON_NULL(fpa, false);
+
+    psArray *chips = fpa->chips;        // Component chips
+    if (chips == NULL) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: fpa->chips == NULL\n");
+        return(0);
+    }
+    if ((chipNum >= chips->n) || (NULL == (pmChip *) chips->data[chipNum])) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: the specified chip (%d) does not exist.\n", chipNum);
+        return(0);
+    }
+
+    int numChips = 0;                   // Number of chips to be processed
+    for (int i = 0 ; i < chips->n ; i++) {
+        pmChip *tmpChip = (pmChip *) chips->data[i]; // Chip of interest
+        if (tmpChip != NULL) {
+            if (i == chipNum) {
+                tmpChip->process = false;
+                setCellsProcess(tmpChip, false); // Wipe out the cell as well
+            } else if (tmpChip->process) {
+                numChips++;
+            }
+        }
+    }
+
+    return(numChips);
+}
+
+
+bool pmCellSetWeights(pmCell *cell // Cell for which to set weights
+                     )
+{
+    float gain = psMetadataLookupF32(NULL, cell->concepts, "CELL.GAIN"); // Cell gain
+    float readnoise = psMetadataLookupF32(NULL, cell->concepts, "CELL.READNOISE"); // Cell read noise
+    psRegion *trimsec = psMetadataLookupPtr(NULL, cell->concepts, "CELL.TRIMSEC"); // Trim section
+
+    p_pmHDU *hdu = cell->hdu;           // The data unit, containing the weight and mask originals
+    if (!hdu) {
+        pmChip *chip = cell->parent;    // The parent chip
+        if (chip->hdu) {
+            hdu = chip->hdu;
+        } else {
+            pmFPA *fpa = chip->parent;  // The parent FPA
+            hdu = fpa->hdu;
+            if (!hdu) {
+                psError(PS_ERR_UNKNOWN, true, "Unable to find the HDU in the hierarchy!\n");
+                return false;
+            }
+        }
+    }
+
+    psArray *pixels = hdu->images;      // Array of images
+    psArray *weights = hdu->weights;    // Array of weight images
+    psArray *masks = hdu->masks;        // Array of mask images
+    // Generate the weights and masks if required
+    if (! weights) {
+        weights = psArrayAlloc(pixels->n);
+        for (int i = 0; i < pixels->n; i++) {
+            psImage *image = pixels->data[i];
+            weights->data[i] = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32);
+            psImageInit(weights->data[i], 0.0);
+        }
+        hdu->weights = weights;
+    }
+    if (! masks) {
+        masks = psArrayAlloc(pixels->n);
+        for (int i = 0; i < pixels->n; i++) {
+            psImage *image = pixels->data[i];
+            masks->data[i] = psImageAlloc(image->numCols, image->numRows, PS_TYPE_U8);
+            psImageInit(masks->data[i], 0);
+        }
+        hdu->masks = masks;
+    }
+
+    // Set the pixels
+    psArray *readouts = cell->readouts; // Array of readouts
+    for (int i = 0; i < readouts->n; i++) {
+        pmReadout *readout = readouts->data[i]; // The readout of interest
+
+        if (! readout->weight) {
+            readout->weight = weights->data[i];
+        }
+        if (! readout->mask) {
+            readout->mask = masks->data[i];
+        }
+
+        // Mask is already set to 0
+
+        // Set weight image to the variance = g*f + rn^2
+        psImage *image = psImageSubset(readout->image, *trimsec); // The pixels
+        psImage *weight = psImageSubset(readout->weight, *trimsec); // The weight map
+        psBinaryOp(weight, image, "*", psScalarAlloc(gain, PS_TYPE_F32));
+        psBinaryOp(weight, weight, "+", psScalarAlloc(readnoise*readnoise, PS_TYPE_F32));
+    }
+
+    return true;
+}
+
+
Index: /branches/eam_rel9_p0/psModules/src/astrom/pmFPA.h
===================================================================
--- /branches/eam_rel9_p0/psModules/src/astrom/pmFPA.h	(revision 6062)
+++ /branches/eam_rel9_p0/psModules/src/astrom/pmFPA.h	(revision 6062)
@@ -0,0 +1,284 @@
+/** @file  pmFPA.h
+*
+*  @brief This file defines the basic types the focal plane hierarchy.
+*
+*  @ingroup AstroImage
+*
+*  @author GLG, MHPCC
+*
+*  @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2006-01-20 02:36:41 $
+*
+*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
+*/
+
+#ifndef PM_FPA_H
+#define PM_FPA_H
+#if HAVE_CONFIG_H
+#include <config.h>
+#endif
+#include "pslib.h"
+#include "psDB.h"
+
+/// @addtogroup AstroImage
+/// @{
+
+// XXX: Is this correct?  Must determine what p_pmHDU is.
+// XXX: Create the p_pmHDU alloc/free.
+typedef struct
+{
+    const char *extname;                // Extension name, if it corresponds to this level
+    psMetadata *header;                 // The FITS header, if it corresponds to this level
+    psArray *images;                    // The pixel data, if it corresponds to this level
+    psArray *masks;                     // The mask data, if it corresponds to this level
+    psArray *weights;                   // The weight data, if it corresponds to this level
+}
+p_pmHDU;
+
+/** Focal plane data structure
+ *
+ *  A focal plane consists of one or more chips (according to the number of
+ *  pieces of contiguous silicon). It contains metadata containers for the
+ *  concepts and analysis, a link to the parent, and pointers to the FITS header,
+ *  if that corresponds to this level (the FPA may be the PHU, but will not ever
+ *  contain pixels). For astrometry, it contains a transformation from the focal
+ *  plane to the tangent plane and the fixed pattern residuals. It is expected
+ *  that the transformation will consist of two 4D polynomials (i.e. a function
+ *  of two coordinates in position, the magnitude of the object, and the color of
+ *  the object) in order to correct for optical distortions and the effects of
+ *  the atmosphere; hence we think that it is prudent to include a reverse
+ *  transformation which will be derived from numerically inverting the forward
+ *  transformation.
+ *
+ */
+typedef struct
+{
+    // Astrometric transformations
+    psPlaneDistort* fromTangentPlane;   ///< Transformation from tangent plane to focal plane
+    psPlaneDistort* toTangentPlane;     ///< Transformation from focal plane to tangent plane
+    psProjection *projection;           ///< Projection from tangent plane to sky
+    // Information
+    psMetadata *concepts;               ///< Cache for PS concepts
+    psMetadata *analysis;               ///< FPA-level analysis metadata
+    const psMetadata *camera;           ///< Camera configuration
+    psArray *chips;                     ///< The chips
+    p_pmHDU *hdu;                       ///< FITS data
+    psMetadata *phu;                    ///< Primary Header
+}
+pmFPA;
+
+/** Chip data structure
+ *
+ *  A chip consists of one or more cells (according to the number of amplifiers
+ *  on the device). The chip contains metadata containers for the concepts and
+ *  analysis, a link to the parent, and pointers to the pointers to the various
+ *  FITS data, if that corresponds to this level. For astrometry, in addition to
+ *  the rough positioning information, it contains a coordinate transform from
+ *  the chip to the focal plane. It is expected that this transform will consist
+ *  of two second-order 2D polynomials; hence we think that it is prudent to
+ *  include a reverse transformation which will be derived from numerically
+ *  inverting the forward transformation. A boolean indicates whether the chip is
+ *  of interest, allowing it to be excluded from analysis.
+ *
+ */
+typedef struct
+{
+    // Offset specifying position on focal plane
+    int col0;                           ///< Offset from the left of FPA.
+    int row0;                           ///< Offset from the bottom of FPA.
+    // Astrometric transformations
+    psPlaneTransform* toFPA;            ///< Transformation from chip to FPA coordinates
+    psPlaneTransform* fromFPA;          ///< Transformation from FPA to chip coordinates
+    // Information
+    psMetadata *concepts;               ///< Cache for PS concepts
+    psMetadata *analysis;               ///< Chip-level analysis metadata
+    psArray *cells;                     ///< The cells (referred to by name)
+    pmFPA *parent;                      ///< Parent FPA
+    bool process;                       ///< Do we bother about reading and working with this chip?
+    bool exists;                        ///< Does the chip exist --- has it been read in?
+    p_pmHDU *hdu;                       ///< FITS data
+}
+pmChip;
+
+/** Cell data structure
+ *
+ *  A cell consists of one or more readouts.  It also contains a pointer to the
+ *  cell's metadata, and its parent chip.  On the astrometry side, it also
+ *  contains coordinate transforms from the cell to chip, from the cell to
+ *  focal-plane, as well as a "quick and dirty" tranform from the cell to
+ *  sky coordinates.
+ *
+ */
+typedef struct
+{
+    // Offset specifying position on chip
+    int col0;                           ///< Offset from the left of chip.
+    int row0;                           ///< Offset from the bottom of chip.
+    // Astrometric transformations
+    psPlaneTransform* toChip;           ///< Transformations from cell to chip coordinates
+    psPlaneTransform* toFPA;            ///< Transformations from cell to FPA coordinates
+    psPlaneTransform* toSky;            ///< Transformations from cell to sky coordinates
+    // Information
+    psMetadata *concepts;               ///< Cache for PS concepts
+    psMetadata *camera;                 ///< Camera Info
+    psMetadata *analysis;               ///< Cell-level analysis metadata
+    psArray *readouts;                  ///< The readouts (referred to by number)
+    pmChip *parent;                     ///< Parent chip
+    bool process;                       ///< Do we bother about reading and working with this cell?
+    bool exists;                        ///< Does the cell exist --- has it been read in?
+    p_pmHDU *hdu;                       ///< FITS data
+}
+pmCell;
+
+/** Readout data structure.
+ *
+ *  A readout is the result of a single read of a cell (or a portion thereof).
+ *  It contains the offset from the lower-left corner of the chip, in the case
+ *  that the CCD was windowed, as well as the binning factors and parity (if the
+ *  binning value is negative, then the parity is reversed). It also contains the
+ *  pixel data, metadata containers for the concepts and analysis, and a link to
+ *  the parent.
+ *
+ */
+typedef struct
+{
+    // Position on the cell
+    // XXX: These may be removed in the future; use parent->concepts instead?
+    int col0;                           ///< Offset from the left of chip.
+    int row0;                           ///< Offset from the bottom of chip.
+    int colBins;                        ///< Amount of binning in x-dimension
+    int rowBins;                        ///< Amount of binning in y-dimension
+    // Information
+    psImage *image;                     ///< Imaging area of readout
+    psImage *mask;                      ///< Mask of input image
+    psImage *weight;                    ///< Weight of input image
+    psMetadata *analysis;               ///< Readout-level analysis metadata
+    pmCell *parent;                     ///< Parent cell
+}
+pmReadout;
+
+
+/** Allocates a pmReadout
+ *
+ *  The constructor shall make an empty pmReadout. If the parent cell is not
+ *  NULL, the parent link is made and the readout shall be placed in the
+ *  parents array of readouts. The metadata containers shall be allocated. All
+ *  other pointers in the structure shall be initialized to NULL.
+ *
+ *  @return pmReadout*    newly allocated pmReadout with all internal pointers set to NULL
+ */
+pmReadout *pmReadoutAlloc(
+    pmCell *cell                        ///< Parent cell
+);
+
+/** Allocates a pmCell
+ *
+ *  The constructor shall make an empty pmCell. If the parent chip is not NULL,
+ *  the parent link is made and the cell shall be placed in the parents array of
+ *  cells. The readouts array shall be allocated with a zero size, and the
+ *  metadata containers constructed. All other pointers in the structure shall be
+ *  initialized to NULL.
+ *
+ *  @return pmCell*    newly allocated pmCell
+ */
+pmCell *pmCellAlloc(
+    pmChip *chip,       ///< Parent chip
+    psMetadata *cameraData, ///< Camera data
+    const char *name    ///< Name of cell
+);
+
+/** Allocates a pmChip
+ *
+ *  The constructor shall make an empty pmChip. If the parent fpa is not NULL,
+ *  the parent link is made and the chip shall be placed in the parent's array
+ *  of chips. The cells array shall be allocated with a zero size, and the
+ *  metadata containers constructed. All other pointers in the structure shall be
+ *  initialized to NULL.
+ *
+ *  @return pmChip*    newly allocated pmChip
+ */
+pmChip *pmChipAlloc(
+    pmFPA *fpa,                         ///< FPA to which the chip belongs
+    const char *name                    ///< Name of chip
+);
+
+/** Allocates a pmFPA
+ *
+ *  The constructor shall make an empty pmFPA. The chips array shall be
+ *  allocated with a zero size, the camera and db pointers set to the values
+ *  provided, and the concepts metadata constructed. All other pointers in the
+ *  structure shall be initialized to NULL.
+ *
+ */
+pmFPA *pmFPAAlloc(
+    const psMetadata *camera            ///< Camera configuration
+);
+
+/** Allocates a p_pmHDU
+ *
+ * XXX: More detailed description
+ *
+ * @return p_pmHDU*    newly allocated p_pmHDU
+ */
+p_pmHDU *p_pmHDUAlloc(const char *extname // Extension name
+                     );
+
+
+/** Verify parent links.
+ *
+ *  This function checks the validity of the parent links in the FPA hierarchy.
+ *  If a parent link is not set (or not set correctly), it is corrected, and the
+ *  function shall return false. If all the parent pointers were correct, the
+ *  function shall return true.
+ *
+ */
+bool pmFPACheckParents(
+    pmFPA *fpa
+);
+
+
+
+/** FUNC DESC
+ *
+ *
+ *
+ *
+ */
+
+
+
+/**
+ *
+ * pmFPASelectChip shall set valid to true for the specified chip number
+ * (chipNum), and all other chips shall have valid set to false. In the event
+ * that the specified chip number does not exist within the fpa, the function
+ * shall return false.
+ *
+ */
+bool pmFPASelectChip(
+    pmFPA *fpa,
+    int chipNum
+);
+
+/**
+ *
+ * pmFPAExcludeChip shall set valid to false only for the specified chip
+ * number (chipNum). In the event that the specified chip number does not exist
+ * within the fpa, the function shall generate a warning, and perform no action.
+ * The function shall return the number of chips within the fpa that have valid
+ * set to true.
+ *
+ */
+int pmFPAExcludeChip(
+    pmFPA *fpa,
+    int chipNum
+);
+
+
+// Set the weights and masks within a cell, based on the gain and RN
+bool pmCellSetWeights(pmCell *cell // Cell for which to set weights
+                     );
+
+
+
+#endif // #ifndef PM_FPA_H
Index: /branches/eam_rel9_p0/psModules/src/astrom/pmFPAAstrometry.c
===================================================================
--- /branches/eam_rel9_p0/psModules/src/astrom/pmFPAAstrometry.c	(revision 6062)
+++ /branches/eam_rel9_p0/psModules/src/astrom/pmFPAAstrometry.c	(revision 6062)
@@ -0,0 +1,512 @@
+#include <stdio.h>
+#include "pslib.h"
+#include "pmFPAAstrometry.h"
+#include "pmFPA.h"
+
+
+/*****************************************************************************
+checkValidImageCoords(): this is a private function which simply determines if
+the supplied x,y coordinates are in the range for the supplied psImage.
+ 
+XXX: What about col0 and row0
+XXX: This should return a psBool.
+XXX: Macro this for speed.
+ *****************************************************************************/
+static psS32 checkValidImageCoords(
+    double x,
+    double y,
+    psImage* tmpImage)
+{
+    PS_ASSERT_IMAGE_NON_NULL(tmpImage, 0);
+
+    // The FLT_EPSILON is because -0.0 was failing this.
+    if (((x+FLT_EPSILON) < 0.0) || (x > (double)tmpImage->numCols) ||
+            ((y+FLT_EPSILON) < 0.0) || (y > (double)tmpImage->numRows)) {
+        return (0);
+    }
+
+    return (1);
+}
+
+/*****************************************************************************/
+/* FUNCTION IMPLEMENTATION - PUBLIC                                          */
+/*****************************************************************************/
+
+pmCell* pmCellInFPA(
+    const psPlane* fpaCoord,
+    const pmFPA* FPA)
+{
+    PS_ASSERT_PTR_NON_NULL(fpaCoord, NULL);
+    PS_ASSERT_PTR_NON_NULL(FPA, NULL);
+
+    pmChip* tmpChip = NULL;
+    psPlane chipCoord;
+    pmCell* outCell = NULL;
+
+    // Determine which chip contains the fpaCoords.
+    tmpChip = pmChipInFPA(fpaCoord, FPA);
+    if (tmpChip == NULL) {
+        return(NULL);
+    }
+
+    // Convert to those chip coordinates.
+    psPlane *rc = pmCoordFPAToChip(&chipCoord, fpaCoord, tmpChip);
+    if (rc == NULL) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine Chip coords.\n");
+        return(NULL);
+    }
+
+    // Determine which cell contains those chip coordinates.
+    outCell = pmCellInChip(&chipCoord, tmpChip);
+    if (outCell == NULL) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine the cell.\n");
+        return(NULL);
+    }
+
+    return (outCell);
+}
+
+pmChip* pmChipInFPA(
+    const psPlane* fpaCoord,
+    const pmFPA* FPA)
+{
+    PS_ASSERT_PTR_NON_NULL(fpaCoord, NULL);
+    PS_ASSERT_PTR_NON_NULL(FPA, NULL);
+    PS_ASSERT_PTR_NON_NULL(FPA->chips, NULL);
+
+    psArray* chips = FPA->chips;
+    psS32 nChips = chips->n;
+    psPlane chipCoord;
+    pmCell *tmpCell = NULL;
+
+    //
+    // Loop through every chip in this FPA.  Convert the original FPA
+    // coordinates to chip coordinates for that chip.  Then, determine if any
+    // cells in that chip contain those chip coordinates.
+    // XXX: Depending on the number of chips, and their topology, there may be
+    // a much more efficient way of doing this.
+    //
+    for (psS32 i = 0; i < nChips; i++) {
+        pmChip* tmpChip = chips->data[i];
+        PS_ASSERT_PTR_NON_NULL(tmpChip, NULL);
+        PS_ASSERT_PTR_NON_NULL(tmpChip->fromFPA, NULL);
+
+        psPlaneTransformApply(&chipCoord, tmpChip->fromFPA, fpaCoord);
+
+        tmpCell = pmCellInChip(&chipCoord, tmpChip);
+        if (tmpCell != NULL) {
+            return(tmpChip);
+        }
+    }
+
+    psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine the chip.\n");
+    return (NULL);
+}
+
+
+pmCell* pmCellInChip(
+    const psPlane* chipCoord,
+    const pmChip* chip)
+{
+    PS_ASSERT_PTR_NON_NULL(chipCoord, NULL);
+    PS_ASSERT_PTR_NON_NULL(chip, NULL);
+
+    psPlane cellCoord;
+    psArray* cells;
+
+    cells = chip->cells;
+    if (cells == NULL) {
+        return NULL;
+    }
+
+    //
+    // We loop over each cell in the chip.  We transform the chipCoord into
+    // a cellCoord for that cell and determine if that cellCoord is valid.
+    // If so, then we return that cell.
+    // XXX: Depending on the number of cells, and their topology, there may be
+    // a much more efficient way of doing this.
+    //
+    for (psS32 i = 0; i < cells->n; i++) {
+        pmCell* tmpCell = (pmCell* ) cells->data[i];
+        PS_ASSERT_PTR_NON_NULL(tmpCell, NULL);
+
+        psPlaneTransform *chipToCell = NULL;
+        if (true ==  p_psIsProjectionLinear(tmpCell->toChip)) {
+            chipToCell = p_psPlaneTransformLinearInvert(tmpCell->toChip);
+        } else {
+            psLogMsg(__func__, PS_LOG_WARN, "WARNING: non-linear cell->chip transforms are not yet implemented.\n");
+            //chipToCell = psPlaneTransformInvert(NULL, tmpCell->toChip, NULL, -1);
+            chipToCell = NULL;
+        }
+        if (chipToCell == NULL) {
+            psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not invert the Cell->toChip transform.\n");
+            return(NULL);
+        }
+        psArray* readouts = tmpCell->readouts;
+
+        if (readouts != NULL) {
+            for (psS32 j = 0; j < readouts->n; j++) {
+                pmReadout* tmpReadout = readouts->data[j];
+                PS_ASSERT_READOUT_NON_NULL(tmpReadout, NULL);
+
+                psPlaneTransformApply(&cellCoord,
+                                      chipToCell,
+                                      chipCoord);
+
+                if (checkValidImageCoords(cellCoord.x,
+                                          cellCoord.y,
+                                          tmpReadout->image)) {
+                    psFree(chipToCell);
+                    return (tmpCell);
+                }
+            }
+        }
+        psFree(chipToCell);
+    }
+
+    //psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine the cell.\n");
+    return (NULL);
+}
+
+
+psPlane* pmCoordCellToFPA(
+    psPlane* fpaCoord,
+    const psPlane* cellCoord,
+    const pmCell* cell)
+{
+    PS_ASSERT_PTR_NON_NULL(cellCoord, NULL);
+    PS_ASSERT_PTR_NON_NULL(cell, NULL);
+
+    psPlane *rc = psPlaneTransformApply(fpaCoord, cell->toFPA, cellCoord);
+    if (rc == NULL) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform cell coords to FPA coords.\n");
+    }
+    return(rc);
+}
+
+
+psPlane* pmCoordChipToFPA(
+    psPlane* outCoord,
+    const psPlane* inCoord,
+    const pmChip* chip)
+{
+    PS_ASSERT_PTR_NON_NULL(inCoord, NULL);
+    PS_ASSERT_PTR_NON_NULL(chip, NULL);
+
+    psPlane *rc = psPlaneTransformApply(outCoord, chip->toFPA, inCoord);
+    if (rc == NULL) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform chip coords to FPA coords.\n");
+    }
+    return(rc);
+}
+
+
+psPlane* pmCoordFPAToChip(
+    psPlane* chipCoord,
+    const psPlane* fpaCoord,
+    const pmChip* chip)
+{
+    PS_ASSERT_PTR_NON_NULL(fpaCoord, NULL);
+    PS_ASSERT_PTR_NON_NULL(chip, NULL);
+    PS_ASSERT_PTR_NON_NULL(chip->fromFPA, NULL);
+
+    psPlane *rc = psPlaneTransformApply(chipCoord, chip->fromFPA, fpaCoord);
+    if (rc == NULL) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform FPA coords to Chip coords.\n");
+    }
+    return(rc);
+}
+
+psPlane* pmCoordCellToChip(
+    psPlane* outCoord,
+    const psPlane* inCoord,
+    const pmCell* cell)
+{
+    PS_ASSERT_PTR_NON_NULL(inCoord, NULL);
+    PS_ASSERT_PTR_NON_NULL(cell, NULL);
+
+    psPlane *rc = psPlaneTransformApply(outCoord, cell->toChip, inCoord);
+    if (rc == NULL) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform Cell coords to Chip coords.\n");
+    }
+    return(rc);
+}
+
+psPlane* pmCoordChipToCell(
+    psPlane* cellCoord,
+    const psPlane* chipCoord,
+    const pmCell* cell)
+{
+    PS_ASSERT_PTR_NON_NULL(chipCoord, NULL);
+    PS_ASSERT_PTR_NON_NULL(cell, NULL);
+    PS_ASSERT_PTR_NON_NULL(cell->parent, NULL);
+
+    pmCell *tmpCell = pmCellInChip(chipCoord, cell->parent);
+    if (tmpCell == NULL) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine the proper cell.\n");
+        return(NULL);
+    }
+
+    psPlaneTransform *tmpChipToCell = NULL;
+    PS_ASSERT_PTR_NON_NULL(tmpCell->toChip, NULL);
+    if (true ==  p_psIsProjectionLinear(tmpCell->toChip)) {
+        tmpChipToCell = p_psPlaneTransformLinearInvert(tmpCell->toChip);
+    } else {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: non-linear cell->chip transforms are not yet implemented.\n");
+        // XXX: tmpChipToCell = psPlaneTransformInvert(NULL, tmpCell->toChip, NULL, -1);
+        tmpChipToCell = NULL;
+    }
+    if (tmpChipToCell == NULL) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not invert the Cell->toChip transform.\n");
+        return(NULL);
+    }
+
+    psPlane *rc = psPlaneTransformApply(cellCoord, tmpChipToCell, chipCoord);
+    if (rc == NULL) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform Chip coords to Cell coords.\n");
+    }
+    psFree(tmpChipToCell);
+    return(rc);
+}
+
+psPlane* pmCoordFPAToTP(
+    psPlane* outCoord,
+    const psPlane* inCoord,
+    double color,
+    double magnitude,
+    const pmFPA* fpa)
+{
+    PS_ASSERT_PTR_NON_NULL(inCoord, NULL);
+    PS_ASSERT_PTR_NON_NULL(fpa, NULL);
+
+    psPlane *rc = psPlaneDistortApply(outCoord, fpa->toTangentPlane, inCoord, color, magnitude);
+    if (rc == NULL) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform FPA coords to tangent plane coords.\n");
+    }
+    return(rc);
+}
+
+psPlane* pmCoordTPToFPA(
+    psPlane* fpaCoord,
+    const psPlane* tpCoord,
+    double color,
+    double magnitude,
+    const pmFPA* fpa)
+{
+    PS_ASSERT_PTR_NON_NULL(tpCoord, NULL);
+    PS_ASSERT_PTR_NON_NULL(fpa, NULL);
+    PS_ASSERT_PTR_NON_NULL(fpa->fromTangentPlane, NULL);
+
+    psPlane *rc = psPlaneDistortApply(fpaCoord, fpa->fromTangentPlane, tpCoord, color, magnitude);
+    if (rc == NULL) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform tangent plane coords to FPA coords.\n");
+    }
+    return(rc);
+}
+
+
+/*****************************************************************************
+XXXDeproject(outSphere, coord, projection): This private routine is a wrapper
+for p_psDeproject().  The reason: p_psDeproject() and p_psProject() combined
+do not seem to produce the original coordinates when they even though they
+should.  XXXDeproject() simply negates the ->r and ->d members of the output
+psSphere if the input ->y is larger than 0.0.  I don't know why it works.
+ 
+I'm guessing the p_psProject() and p_psDeproject() functions have bugs.
+ 
+XXX: It appears that p_psProject() and p_psDeproject() have been fixed.
+Remove this.
+ *****************************************************************************/
+psSphere* XXXDeproject(
+    psSphere *outSphere,
+    const psPlane* coord,
+    const psProjection* projection)
+{
+    psSphere *rc = p_psDeproject(outSphere, coord, projection);
+
+    if (coord->y >= 0.0) {
+        rc->d = -rc->d;
+        rc->r = -rc->r;
+    }
+
+    return(rc);
+}
+
+/*****************************************************************************
+  *****************************************************************************/
+psSphere* pmCoordTPToSky(
+    psSphere* outSphere,
+    const psPlane* tpCoord,
+    const psProjection *projection)
+{
+    PS_ASSERT_PTR_NON_NULL(tpCoord, NULL);
+    PS_ASSERT_PTR_NON_NULL(projection, NULL);
+
+    //    psSphere *rc = XXXDeproject(outSphere, tpCoord, projection);
+    psSphere *rc = p_psDeproject(outSphere, tpCoord, projection);
+    if (rc == NULL) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform tangent plane coords to sky coords.\n");
+    }
+    return(rc);
+}
+
+/*****************************************************************************
+ *****************************************************************************/
+psPlane* pmCoordSkyToTP(
+    psPlane* tpCoord,
+    const psSphere* in,
+    const psProjection *projection)
+{
+    PS_ASSERT_PTR_NON_NULL(in, NULL);
+    PS_ASSERT_PTR_NON_NULL(projection, NULL);
+
+    psPlane *rc = p_psProject(tpCoord, in, projection);
+    if (rc == NULL) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform sky to tangent plane coords.\n");
+    }
+    return(rc);
+}
+
+/*****************************************************************************
+ *****************************************************************************/
+psSphere* pmCoordCellToSky(
+    psSphere* skyCoord,
+    const psPlane* cellCoord,
+    double color,
+    double magnitude,
+    const pmCell* cell)
+{
+    PS_ASSERT_PTR_NON_NULL(cellCoord, NULL);
+    PS_ASSERT_PTR_NON_NULL(cell, NULL);
+    PS_ASSERT_PTR_NON_NULL(cell->toFPA, NULL);
+    PS_ASSERT_PTR_NON_NULL(cell->parent, NULL);
+    PS_ASSERT_PTR_NON_NULL(cell->parent->parent, NULL);
+    PS_ASSERT_PTR_NON_NULL(cell->parent->parent->toTangentPlane, NULL);
+    PS_ASSERT_PTR_NON_NULL(cell->parent->parent->projection, NULL);
+    psPlane fpaCoord;
+    psPlane tpCoord;
+    psPlane *rc;
+    pmFPA* parFPA = (cell->parent)->parent;
+
+    // Convert the input cell coordinates to FPA coordinates.
+    rc = psPlaneTransformApply(&fpaCoord, cell->toFPA, cellCoord);
+    if (rc == NULL) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could transform cell coords to FPA coords.\n");
+        return(NULL);
+    }
+
+    // Convert the FPA coordinates to tangent plane Coordinates.
+    rc = psPlaneDistortApply(&tpCoord, parFPA->toTangentPlane, &fpaCoord, color, magnitude);
+    if (rc == NULL) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could transform FPA coords to tangent plane coords.\n");
+        return(NULL);
+    }
+
+    // Convert the tangent plane Coordinates to sky coordinates.
+    psSphere *rc2 = pmCoordTPToSky(skyCoord, &tpCoord, parFPA->projection);
+    if (rc2 == NULL) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform cell coords to sky coords.\n");
+    }
+
+    return(rc2);
+}
+
+/*****************************************************************************
+ *****************************************************************************/
+psPlane* pmCoordSkyToCell(
+    psPlane* cellCoord,
+    const psSphere* skyCoord,
+    float color,
+    float magnitude,
+    const pmCell* cell)
+{
+    PS_ASSERT_PTR_NON_NULL(skyCoord, NULL);
+    PS_ASSERT_PTR_NON_NULL(cell, NULL);
+    PS_ASSERT_PTR_NON_NULL(cell->parent, NULL);
+    PS_ASSERT_PTR_NON_NULL(cell->parent->parent, NULL);
+    pmChip *parChip = cell->parent;
+    pmFPA *parFPA = parChip->parent;
+    psPlane tpCoord;
+    psPlane fpaCoord;
+    psPlane chipCoord;
+    psPlane *rc;
+
+    // Convert the skyCoords to tangent plane coords.
+    rc = pmCoordSkyToTP(&tpCoord, skyCoord, parFPA->projection);
+    if (rc == NULL) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine tangent plane coords.\n");
+        return(NULL);
+    }
+
+    // Convert the tangent plane coords to FPA coords.
+    rc = pmCoordTPToFPA(&fpaCoord, &tpCoord, color, magnitude, parFPA);
+    if (rc == NULL) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine FPA coords.\n");
+        return(NULL);
+    }
+
+    // Convert the FPA coords to chip coords.
+    rc = pmCoordFPAToChip(&chipCoord, &fpaCoord, parChip);
+    if (rc == NULL) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine chip coords.\n");
+        return(NULL);
+    }
+
+    // Convert the chip coords to cell coords.
+    rc = pmCoordChipToCell(cellCoord, &chipCoord, cell);
+    if (rc == NULL) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not determine cell coords.\n");
+        return(NULL);
+    }
+
+    return (cellCoord);
+}
+
+/*****************************************************************************
+ *****************************************************************************/
+psSphere* pmCoordCellToSkyQuick(
+    psSphere* outSphere,
+    const psPlane* cellCoord,
+    const pmCell* cell)
+{
+    PS_ASSERT_PTR_NON_NULL(cellCoord, NULL);
+    PS_ASSERT_PTR_NON_NULL(cell, NULL);
+    PS_ASSERT_PTR_NON_NULL(cell->toSky, NULL);
+    psPlane outPlane;
+    psPlane *rc;
+    rc = psPlaneTransformApply(&outPlane, cell->toSky, cellCoord);
+    if (rc == NULL) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could transform cell coords to sky coords.\n");
+        return(NULL);
+    }
+
+    psSphere *out = outSphere;
+    if (out == NULL) {
+        out = psSphereAlloc();
+    }
+    out->r = outPlane.y;
+    out->d = outPlane.x;
+
+    return(out);
+}
+
+/*****************************************************************************
+ *****************************************************************************/
+psPlane* pmCoordSkyToCellQuick(
+    psPlane* cellCoord,
+    const psSphere* skyCoord,
+    const pmCell* cell)
+{
+    PS_ASSERT_PTR_NON_NULL(skyCoord, NULL);
+    PS_ASSERT_PTR_NON_NULL(cell, NULL);
+    PS_ASSERT_PTR_NON_NULL(cell->toSky, NULL);
+    psPlane skyPlane;
+    skyPlane.y = skyCoord->r;
+    skyPlane.x = skyCoord->d;
+
+    psPlane *rc = psPlaneTransformApply(cellCoord, cell->toSky, &skyPlane);
+    if (rc == NULL) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not transform sky to cell coords.\n");
+    }
+    return(cellCoord);
+}
Index: /branches/eam_rel9_p0/psModules/src/astrom/pmFPAAstrometry.h
===================================================================
--- /branches/eam_rel9_p0/psModules/src/astrom/pmFPAAstrometry.h	(revision 6062)
+++ /branches/eam_rel9_p0/psModules/src/astrom/pmFPAAstrometry.h	(revision 6062)
@@ -0,0 +1,188 @@
+#ifndef PM_FPA_ASTROMETRY_H
+#define PM_FPA_ASTROMETRY_H
+
+#include "pslib.h"
+#include "pmFPA.h"
+
+/** Find cooresponding cell for given FPA coordinate
+ *
+ *  @return pmCell*    the cell cooresponding to the coord in FPA
+ */
+pmCell* pmCellInFPA(
+    const psPlane* coord,              ///< the coordinate in FPA plane
+    const pmFPA* FPA                   ///< the FPA to search for the cell
+);
+
+
+/** Find cooresponding chip for given FPA coordinate
+ *
+ *  @return pmChip*    the chip cooresponding to coord
+ */
+pmChip* pmChipInFPA(
+    const psPlane* coord,              ///< the coordinate in FPA plane
+    const pmFPA* FPA                   ///< the FPA to search for the cell
+);
+
+
+/** Find cooresponding cell for given Chip coordinate
+ *
+ *  @return pmCell*    the cell cooresponding to coord
+ */
+pmCell* pmCellInChip(
+    const psPlane* coord,              ///< the coordinate in Chip plane
+    const pmChip* chip                 ///< the chip to search for the cell
+);
+
+
+/** Translate a cell coordinate into a chip coordinate
+ *
+ *  @return psPlane*    the resulting chip coordinate
+ */
+psPlane* pmCoordCellToChip(
+    psPlane* out,                      ///< a plane struct to recycle. If NULL, a new struct is created
+    const psPlane* in,                 ///< the coordinate within Cell
+    const pmCell* cell                 ///< the Cell in interest
+);
+
+
+/** Translate a chip coordinate into a FPA coordinate
+ *
+ *  @return psPlane*    the resulting FPA coordinate
+ */
+psPlane* pmCoordChipToFPA(
+    psPlane* out,                      ///< a plane struct to recycle. If NULL, a new struct is created
+    const psPlane* in,                 ///< the coordinate within Chip
+    const pmChip* chip                 ///< the chip in interest
+);
+
+
+/** Translate a FPA coordinate into a Tangent Plane coordinate
+ *
+ *  @return psPlane*    the resulting Tangent Plane coordinate
+ */
+psPlane* pmCoordFPAToTP(
+    psPlane* out,                      ///< a plane struct to recycle. If NULL, a new struct is created
+    const psPlane* in,                 ///< the coordinate within FPA
+    double color,                      ///< Color of source
+    double magnitude,                  ///< Magnitude of source
+    const pmFPA* fpa                   ///< the FPA in interest
+);
+
+
+/** Translate a Tangent Plane coordinate into a Sky coordinate
+ *
+ *  @return psSphere*    the resulting Sky coordinate
+ */
+psSphere* pmCoordTPToSky(
+    psSphere* out,                     ///< a sphere struct to recycle. If NULL, a new struct is created
+    const psPlane* in,                ///< the coordinate within Tangent Plane
+    const psProjection *projection
+);
+
+/** Translate a cell coordinate into a FPA coordinate
+ *
+ *  @return psPlane*    the resulting FPA coordinate
+ */
+psPlane* pmCoordCellToFPA(
+    psPlane* out,                      ///< a plane struct to recycle. If NULL, a new struct is created
+    const psPlane* in,                 ///< the coordinate within cell
+    const pmCell* cell                 ///< the cell in interest
+);
+
+
+/** Translate a cell coordinate into a Sky coordinate
+ *
+ *  @return psSphere*    the resulting Sky coordinate
+ */
+psSphere* pmCoordCellToSky(
+    psSphere* out,                     ///< a sphere struct to recycle. If NULL, a new struct is created
+    const psPlane* in,                 ///< the coordinate within cell
+    double color,                      ///< Color of source
+    double magnitude,                  ///< Magnitude of source
+    const pmCell* cell                 ///< the cell in interest
+);
+
+
+/** Translate a cell coordinate into a Sky coordinate using a 'quick and
+ *  dirty' method
+ *
+ *  @return psSphere*    the resulting Sky coordinate
+ */
+psSphere* pmCoordCellToSkyQuick(
+    psSphere* out,                     ///< a sphere struct to recycle. If NULL, a new struct is created
+    const psPlane* in,                 ///< the coordinate within cell
+    const pmCell* cell                 ///< the cell in interest
+);
+
+
+/** Translate a Sky coordinate into a Tangent Plane coordinate
+ *
+ *  @return psPlane*    the resulting Tangent Plane coordinate
+ */
+psPlane* pmCoordSkyToTP(
+    psPlane* out,                      ///< a plane struct to recycle. If NULL, a new struct is created
+    const psSphere* in,                ///< the sky coordinate
+    const psProjection *projection
+);
+
+/** Translate a Tangent Plane coordinate into a FPA coordinate
+ *
+ *  @return psPlane*    the resulting FPA coordinate
+ */
+psPlane* pmCoordTPToFPA(
+    psPlane* out,                      ///< a plane struct to recycle. If NULL, a new struct is created
+    const psPlane* in,                 ///< the coordinate within tangent plane
+    double color,                      ///< Color of source
+    double magnitude,                  ///< Magnitude of source
+    const pmFPA* fpa                   ///< the FPA of interest
+);
+
+
+/** Translate a FPA coordinate into a chip coordinate
+ *
+ *  @return psPlane*    the resulting chip coordinate
+ */
+psPlane* pmCoordFPAToChip(
+    psPlane* out,                      ///< a plane struct to recycle. If NULL, a new struct is created
+    const psPlane* in,                 ///< the FPA coordinate
+    const pmChip* chip                 ///< the chip of interest
+);
+
+
+/** Translate a chip coordinate into a cell coordinate
+ *
+ *  @return psPlane*    the resulting cell coordinate
+ */
+psPlane* pmCoordChipToCell(
+    psPlane* out,                      ///< a plane struct to recycle. If NULL, a new struct is created
+    const psPlane* in,                 ///< the Chip coordinate
+    const pmCell* cell                 ///< the cell of interest
+);
+
+
+/** Translate a sky coordinate into a cell coordinate
+ *
+ *  @return psPlane*    the resulting cell coordinate
+ */
+psPlane* pmCoordSkyToCell(
+    psPlane* out,                      ///< a plane struct to recycle. If NULL, a new struct is created
+    const psSphere* in,                ///< the Sky coordinate
+    float color,                       ///< Color of source
+    float magnitude,                   ///< Magnitude of source
+    const pmCell* cell                 ///< the cell of interest
+);
+
+
+/** Translate a sky coordinate into a cell coordinate using a 'quick and
+ *  dirty' method
+ *
+ *  @return psPlane*    the resulting cell coordinate
+ */
+psPlane* pmCoordSkyToCellQuick(
+    psPlane* out,                      ///< a plane struct to recycle. If NULL, a new struct is created
+    const psSphere* in,                ///< the Sky coordinate
+    const pmCell* cell                 ///< the cell of interest
+);
+
+
+#endif
Index: /branches/eam_rel9_p0/psModules/src/astrom/pmFPAConstruct.c
===================================================================
--- /branches/eam_rel9_p0/psModules/src/astrom/pmFPAConstruct.c	(revision 6061)
+++ /branches/eam_rel9_p0/psModules/src/astrom/pmFPAConstruct.c	(revision 6062)
@@ -3,8 +3,7 @@
 #include "pslib.h"
 
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 #include "pmFPAConstruct.h"
 #include "psAdditionals.h"
-#include "pmFPAConceptsGet.h"
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
Index: /branches/eam_rel9_p0/psModules/src/astrom/pmFPAConstruct.h
===================================================================
--- /branches/eam_rel9_p0/psModules/src/astrom/pmFPAConstruct.h	(revision 6061)
+++ /branches/eam_rel9_p0/psModules/src/astrom/pmFPAConstruct.h	(revision 6062)
@@ -3,5 +3,5 @@
 
 #include "pslib.h"
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 
 // Read the contents of a FITS file (format specified by the camera configuration) into memory
Index: /branches/eam_rel9_p0/psModules/src/astrom/pmFPARead.c
===================================================================
--- /branches/eam_rel9_p0/psModules/src/astrom/pmFPARead.c	(revision 6061)
+++ /branches/eam_rel9_p0/psModules/src/astrom/pmFPARead.c	(revision 6062)
@@ -3,5 +3,5 @@
 #include "pslib.h"
 
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 #include "pmFPARead.h"
 #include "pmConcepts.h"
@@ -94,5 +94,7 @@
     hdu->header = header;
 
-    // XXX: Insert mask and weight reading here
+    // Mask and weight not set yet
+    hdu->weights = NULL;
+    hdu->masks = NULL;
 
     return true;
Index: /branches/eam_rel9_p0/psModules/src/astrom/pmFPARead.h
===================================================================
--- /branches/eam_rel9_p0/psModules/src/astrom/pmFPARead.h	(revision 6061)
+++ /branches/eam_rel9_p0/psModules/src/astrom/pmFPARead.h	(revision 6062)
@@ -3,5 +3,5 @@
 
 #include "pslib.h"
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 
 bool pmFPARead(pmFPA *fpa,              // FPA to read into
Index: /branches/eam_rel9_p0/psModules/src/astrom/pmFPAWrite.c
===================================================================
--- /branches/eam_rel9_p0/psModules/src/astrom/pmFPAWrite.c	(revision 6061)
+++ /branches/eam_rel9_p0/psModules/src/astrom/pmFPAWrite.c	(revision 6062)
@@ -3,5 +3,5 @@
 #include "pslib.h"
 
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 #include "pmFPARead.h"
 #include "pmConcepts.h"
Index: /branches/eam_rel9_p0/psModules/src/astrom/pmFPAWrite.h
===================================================================
--- /branches/eam_rel9_p0/psModules/src/astrom/pmFPAWrite.h	(revision 6061)
+++ /branches/eam_rel9_p0/psModules/src/astrom/pmFPAWrite.h	(revision 6062)
@@ -3,5 +3,5 @@
 
 #include "pslib.h"
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 
 bool pmFPAWrite(psFits *fits,           // FITS file to which to write
Index: /branches/eam_rel9_p0/psModules/src/astrom/pmReadout.c
===================================================================
--- /branches/eam_rel9_p0/psModules/src/astrom/pmReadout.c	(revision 6061)
+++ /branches/eam_rel9_p0/psModules/src/astrom/pmReadout.c	(revision 6062)
@@ -2,5 +2,5 @@
 #include "pslib.h"
 
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 
 // Get the bias images for a readout, using the CELL.BIASSEC
Index: /branches/eam_rel9_p0/psModules/src/astrom/pmReadout.h
===================================================================
--- /branches/eam_rel9_p0/psModules/src/astrom/pmReadout.h	(revision 6061)
+++ /branches/eam_rel9_p0/psModules/src/astrom/pmReadout.h	(revision 6062)
@@ -3,5 +3,5 @@
 
 #include "pslib.h"
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 
 // Get the bias sections for a specific readout
Index: /branches/eam_rel9_p0/psModules/src/detrend/pmFlatField.c
===================================================================
--- /branches/eam_rel9_p0/psModules/src/detrend/pmFlatField.c	(revision 6061)
+++ /branches/eam_rel9_p0/psModules/src/detrend/pmFlatField.c	(revision 6062)
@@ -24,6 +24,6 @@
  *  @author Ross Harman, MHPCC
  *
- *  @version $Revision: 1.4.8.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-17 03:18:39 $
+ *  @version $Revision: 1.4.8.1.2.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-20 02:38:28 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -44,8 +44,6 @@
 
 
-bool pmFlatField(pmReadout *in, pmReadout *mask, const pmReadout *flat)
+bool pmFlatField(pmReadout *in, const pmReadout *flat)
 {
-    // XXX: Not sure if this is correct.  Must consult with IfA.
-    PS_ASSERT_PTR_NON_NULL(mask, false);
     int i = 0;
     int j = 0;
@@ -80,5 +78,5 @@
         return false;
     }
-    inMask = mask->image;
+    inMask = in->mask;
 
     // Check input image and its mask are not larger than flat image
@@ -90,5 +88,5 @@
         return false;
     }
-    if (inMask->numRows > flatImage->numRows || inMask->numCols > flatImage->numCols) {
+    if (inMask && (inMask->numRows > flatImage->numRows || inMask->numCols > flatImage->numCols)) {
         psError( PS_ERR_BAD_PARAMETER_SIZE, true,
                  PS_ERRORTEXT_pmFlatField_SIZE_MASK_IMAGE,
@@ -112,5 +110,5 @@
                  totOffRow, totOffCol, inImage->numRows, inImage->numCols);
         return false;
-    } else if(totOffRow>=inMask->numRows || totOffCol>=inMask->numCols) {
+    } else if(inMask && (totOffRow>=inMask->numRows || totOffCol>=inMask->numCols)) {
         psError( PS_ERR_BAD_PARAMETER_SIZE, true,
                  PS_ERRORTEXT_pmFlatField_OFFSET_MASK_IMAGE,
@@ -122,5 +120,4 @@
     inType = inImage->type.type;
     flatType = flatImage->type.type;
-    maskType = inMask->type.type;
     if(PS_IS_PSELEMTYPE_COMPLEX(inType)) {
         psError( PS_ERR_BAD_PARAMETER_TYPE, true,
@@ -133,5 +130,5 @@
                  flatType);
         return false;
-    } else if(maskType != PS_TYPE_MASK) {
+    } else if(inMask && inMask->type.type != PS_TYPE_MASK) {
         psError( PS_ERR_BAD_PARAMETER_TYPE, true,
                  PS_ERRORTEXT_pmFlatField_TYPE_MASK_IMAGE,
@@ -153,5 +150,7 @@
             if(flatImage->data.TYPE[j][i] <= 0.0) {                                                          \
                 /* Negative or zero flat pixels shall be masked in input image as  PM_MASK_FLAT */           \
-                inMask->data.PS_TYPE_MASK_DATA[j][i] |= PM_MASK_FLAT;                                        \
+                if (inMask) {                                                                                \
+                    inMask->data.PS_TYPE_MASK_DATA[j][i] |= PM_MASK_FLAT;                                    \
+                }                                                                                            \
                 flatImage->data.TYPE[j][i] = 0.0;                                                            \
             }                                                                                                \
@@ -160,5 +159,5 @@
     for(j = totOffRow; j < inImage->numRows; j++) {                                                          \
         for(i = totOffCol; i < inImage->numCols; i++) {                                                      \
-            if(!inMask->data.PS_TYPE_MASK_DATA[j][i]) {                                                      \
+            if(inMask && !inMask->data.PS_TYPE_MASK_DATA[j][i]) {                                            \
                 /* Module shall divide the input image by the flat-fielded image */                          \
                 inImage->data.TYPE[j][i] /= flatImage->data.TYPE[j][i];                                      \
Index: /branches/eam_rel9_p0/psModules/src/detrend/pmFlatField.h
===================================================================
--- /branches/eam_rel9_p0/psModules/src/detrend/pmFlatField.h	(revision 6061)
+++ /branches/eam_rel9_p0/psModules/src/detrend/pmFlatField.h	(revision 6062)
@@ -24,6 +24,6 @@
  *  @author Ross Harman, MHPCC
  *
- *  @version $Revision: 1.2.8.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-17 03:18:39 $
+ *  @version $Revision: 1.2.8.1.2.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-20 02:38:28 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -31,5 +31,5 @@
 
 #include "pslib.h"
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 
 
@@ -43,5 +43,4 @@
 bool pmFlatField(
     pmReadout *in,          ///< Readout with input image
-    pmReadout *mask,        ///< Input image mask
     const pmReadout *flat   ///< Readout with flat image
 );
Index: /branches/eam_rel9_p0/psModules/src/detrend/pmMaskBadPixels.h
===================================================================
--- /branches/eam_rel9_p0/psModules/src/detrend/pmMaskBadPixels.h	(revision 6061)
+++ /branches/eam_rel9_p0/psModules/src/detrend/pmMaskBadPixels.h	(revision 6062)
@@ -24,6 +24,6 @@
  *  @author Ross Harman, MHPCC
  *
- *  @version $Revision: 1.2.8.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-17 03:18:39 $
+ *  @version $Revision: 1.2.8.1.2.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-20 02:38:28 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -31,5 +31,5 @@
 
 #include "pslib.h"
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 
 /** Mask values */
Index: /branches/eam_rel9_p0/psModules/src/detrend/pmNonLinear.h
===================================================================
--- /branches/eam_rel9_p0/psModules/src/detrend/pmNonLinear.h	(revision 6061)
+++ /branches/eam_rel9_p0/psModules/src/detrend/pmNonLinear.h	(revision 6062)
@@ -10,6 +10,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.2.8.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-17 03:18:39 $
+ *  @version $Revision: 1.2.8.1.2.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-20 02:38:28 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -21,5 +21,5 @@
 
 #include "pslib.h"
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 
 pmReadout *pmNonLinearityPolynomial(pmReadout *in,
Index: /branches/eam_rel9_p0/psModules/src/imcombine/pmReadoutCombine.h
===================================================================
--- /branches/eam_rel9_p0/psModules/src/imcombine/pmReadoutCombine.h	(revision 6061)
+++ /branches/eam_rel9_p0/psModules/src/imcombine/pmReadoutCombine.h	(revision 6062)
@@ -5,6 +5,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.2.10.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-31 04:34:52 $
+ *  @version $Revision: 1.2.10.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-20 02:38:28 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -23,5 +23,5 @@
 #include "pslib.h"
 #include "psConstants.h"
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 
 typedef struct
Index: /branches/eam_rel9_p0/psModules/src/imsubtract/pmSubtractBias.c
===================================================================
--- /branches/eam_rel9_p0/psModules/src/imsubtract/pmSubtractBias.c	(revision 6061)
+++ /branches/eam_rel9_p0/psModules/src/imsubtract/pmSubtractBias.c	(revision 6062)
@@ -11,6 +11,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.6.8.1.2.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-31 04:35:58 $
+ *  @version $Revision: 1.6.8.1.2.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-20 02:38:28 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -22,8 +22,14 @@
 #endif
 
+#include <assert.h>
 #include "pmSubtractBias.h"
 
 #define PM_SUBTRACT_BIAS_POLYNOMIAL_ORDER 2
 #define PM_SUBTRACT_BIAS_SPLINE_ORDER 3
+
+
+#define MAX(a,b) ((a) > (b) ? (a) : (b))
+#define MIN(a,b) ((a) < (b) ? (a) : (b))
+
 
 // XXX: put these in psConstants.h
@@ -56,4 +62,31 @@
 }\
 
+
+void overscanOptionsFree(pmOverscanOptions *options)
+{
+    psFree(options->stat);
+    psFree(options->poly);
+    psFree(options->spline);
+}
+
+pmOverscanOptions *pmOverscanOptionsAlloc(bool single, pmFit fitType, unsigned int order, psStats *stat)
+{
+    pmOverscanOptions *opts = psAlloc(sizeof(pmOverscanOptions));
+    psMemSetDeallocator(opts, (psFreeFunc)overscanOptionsFree);
+
+    // Inputs
+    opts->single = single;
+    opts->fitType = fitType;
+    opts->order = order;
+    opts->stat = psMemIncrRefCounter(stat);
+
+    // Outputs
+    opts->poly = NULL;
+    opts->spline = NULL;
+
+    return opts;
+}
+
+
 /******************************************************************************
 psSubtractFrame(): this routine will take as input a readout for the input
@@ -61,40 +94,60 @@
 place from the input image.
 *****************************************************************************/
-static pmReadout *SubtractFrame(pmReadout *in,
-                                const pmReadout *bias)
-{
-    psS32 i;
-    psS32 j;
-
-    if (bias == NULL) {
-        psLogMsg(__func__, PS_LOG_NOTE,
-                 "WARNING: pmSubtractBias.c: SubtractFrame(): bias frame is NULL.  Returning original image.\n");
-        return(in);
-    }
-
-
-    if ((in->image->numRows + in->row0 - bias->row0) > bias->image->numRows) {
-        psError(PS_ERR_UNKNOWN,true, "bias image does not have enough rows.  Returning in image\n");
-        return(in);
-    }
-    if ((in->image->numCols + in->col0 - bias->col0) > bias->image->numCols) {
-        psError(PS_ERR_UNKNOWN,true, "bias image does not have enough columns.  Returning in image\n");
-        return(in);
-    }
-
-    for (i=0;i<in->image->numRows;i++) {
-        for (j=0;j<in->image->numCols;j++) {
-            in->image->data.F32[i][j]-=
-                bias->image->data.F32[i+in->row0-bias->row0][j+in->col0-bias->col0];
-            if ((in->mask != NULL) && (bias->mask != NULL)) {
-                (in->mask->data.U8[i][j])|=
-                    bias->mask->data.U8[i+in->row0-bias->row0][j+in->col0-bias->col0];
+static bool SubtractFrame(pmReadout *in,// Input readout
+                          const pmReadout *sub, // Readout to be subtracted from input
+                          float scale   // Scale to apply before subtracting
+                         )
+{
+    assert(in);
+    assert(sub);
+
+    // Get the trim sections
+    psRegion *inTrimsec = psMetadataLookupPtr(NULL, in->parent->concepts, "CELL.TRIMSEC");
+    psRegion *subTrimsec = psMetadataLookupPtr(NULL, sub->parent->concepts, "CELL.TRIMSEC");
+    psImage *inImage = psImageSubset(in->image, *inTrimsec); // The input image
+    psImage *subImage = psImageSubset(sub->image, *subTrimsec); // The image to be subtracted
+    psImage *inMask = in->mask ? psImageSubset(in->mask, *inTrimsec) : NULL; // The input mask
+    psImage *subMask = sub->mask ? psImageSubset(sub->mask, *subTrimsec) : NULL; // The input mask
+
+    // Offsets of the cells
+    int x0in = psMetadataLookupS32(NULL, in->parent->concepts, "CELL.X0");
+    int y0in = psMetadataLookupS32(NULL, in->parent->concepts, "CELL.Y0");
+    int x0sub = psMetadataLookupS32(NULL, sub->parent->concepts, "CELL.X0");
+    int y0sub = psMetadataLookupS32(NULL, sub->parent->concepts, "CELL.Y0");
+
+    if ((inImage->numCols + x0in - x0sub) > subImage->numCols) {
+        psError(PS_ERR_UNKNOWN, true, "Image does not have enough columns for subtraction.\n");
+        return false;
+    }
+    if ((inImage->numRows + y0in - y0sub) > subImage->numRows) {
+        psError(PS_ERR_UNKNOWN, true, "Image does not have enough rows for subtraction.\n");
+        return false;
+    }
+
+    if (scale == 1.0) {
+        for (int i = 0; i < inImage->numRows; i++) {
+            for (int j = 0; j < inImage->numCols; j++) {
+                inImage->data.F32[i][j] -= subImage->data.F32[i+y0in-y0sub][j+x0in-x0sub];
+                if (inMask && subMask) {
+                    inMask->data.U8[i][j] |= subMask->data.U8[i+y0in-y0sub][j+x0in-x0sub];
+                }
             }
         }
-    }
-
-    return(in);
-}
-
+    } else {
+        for (int i = 0; i < inImage->numRows; i++) {
+            for (int j = 0; j < inImage->numCols; j++) {
+                inImage->data.F32[i][j] -= subImage->data.F32[i+y0in-y0sub][j+x0in-x0sub] * scale;
+                if (inMask && subMask) {
+                    inMask->data.U8[i][j] |= subMask->data.U8[i+y0in-y0sub][j+x0in-x0sub];
+                }
+            }
+        }
+    }
+
+    return true;
+}
+
+
+#if 0
 /******************************************************************************
 ImageSubtractScalar(): subtract a scalar from the input image.
@@ -114,4 +167,5 @@
     return(image);
 }
+#endif
 
 /******************************************************************************
@@ -180,4 +234,5 @@
 
 
+#if 0
 /******************************************************************************
 ScaleOverscanVector(): this routine takes as input an arbitrary vector,
@@ -283,15 +338,87 @@
 }
 
+#endif
+
+// Produce an overscan vector from an array of pixels
+static psVector *overscanVector(pmOverscanOptions *overscanOpts, // Overscan options
+                                const psArray *pixels, // Array of vectors containing the pixel values
+                                psStats *myStats // Statistic to use in reducing the overscan
+                               )
+{
+    // Reduce the overscans
+    psVector *reduced = psVectorAlloc(pixels->n, PS_TYPE_F32); // Overscan for each row
+    psVector *ordinate = psVectorAlloc(pixels->n, PS_TYPE_F32); // Ordinate
+    psVector *mask = psVectorAlloc(pixels->n, PS_TYPE_U8); // Mask for fitting
+    for (int i = 0; i < pixels->n; i++) {
+        psVector *values = pixels->data[i]; // Vector with overscan values
+        if (values->n > 0) {
+            mask->data.U8[i] = 0;
+            ordinate->data.F32[i] = 2.0*(float)i/(float)pixels->n - 1.0; // Scale to [-1,1]
+            psVectorStats(myStats, values, NULL, NULL, 0);
+            double reducedVal = NAN; // Result of statistics
+            if (! p_psGetStatValue(myStats, &reducedVal)) {
+                psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result "
+                        "of statistics on row %d.\n", i);
+                return NULL;
+            }
+            reduced->data.F32[i] = reducedVal;
+            //            printf("%f ", reducedVal);
+        } else if (overscanOpts->fitType == PM_FIT_NONE) {
+            psError(PS_ERR_UNKNOWN, true, "The overscan is not supplied for all points on the "
+                    "image, and no fit is requested.\n");
+            return NULL;
+        } else {
+            // We'll fit this one out
+            mask->data.U8[i] = 1;
+        }
+
+        //        printf("\n");
+    }
+
+    // Fit the overscan, if required
+    switch (overscanOpts->fitType) {
+    case PM_FIT_NONE:
+        // No fitting --- that's easy.
+        break;
+    case PM_FIT_POLY_ORD:
+        overscanOpts->poly = psPolynomial1DAlloc(overscanOpts->order, PS_POLYNOMIAL_ORD);
+        overscanOpts->poly = psVectorFitPolynomial1D(overscanOpts->poly, mask, 1, reduced, NULL,
+                             ordinate);
+        psFree(reduced);
+        reduced = psPolynomial1DEvalVector(overscanOpts->poly, ordinate);
+        break;
+    case PM_FIT_POLY_CHEBY:
+        overscanOpts->poly = psPolynomial1DAlloc(overscanOpts->order, PS_POLYNOMIAL_CHEB);
+        overscanOpts->poly = psVectorFitPolynomial1D(overscanOpts->poly, mask, 1, reduced, NULL,
+                             ordinate);
+        psFree(reduced);
+        reduced = psPolynomial1DEvalVector(overscanOpts->poly, ordinate);
+        break;
+    case PM_FIT_SPLINE:
+        // XXX I don't think psSpline1D is up to scratch yet --- it has no mask, and requires an
+        // input spline
+        overscanOpts->spline = psVectorFitSpline1D(reduced, ordinate);
+        psFree(reduced);
+        reduced = psSpline1DEvalVector(overscanOpts->spline, ordinate);
+        break;
+    default:
+        psError(PS_ERR_UNKNOWN, true, "Unknown value for the fitting type: %d\n", overscanOpts->fitType);
+        return NULL;
+        break;
+    }
+
+    psFree(ordinate);
+    psFree(mask);
+
+    return reduced;
+}
+
+
+
 /******************************************************************************
 XXX: The SDRS does not specify type support.  F32 is implemented here.
  *****************************************************************************/
-pmReadout *pmSubtractBias(pmReadout *in,
-                          void *fitSpec,
-                          const psList *overscans,
-                          pmOverscanAxis overScanAxis,
-                          psStats *stat,
-                          psS32 nBinOrig,
-                          pmFit fit,
-                          const pmReadout *bias)
+pmReadout *pmSubtractBias(pmReadout *in, pmOverscanOptions *overscanOpts,
+                          const pmReadout *bias, const pmReadout *dark)
 {
     psTrace(".psModule.pmSubtracBias.pmSubtractBias", 4,
@@ -301,390 +428,166 @@
     PS_ASSERT_READOUT_TYPE(in, PS_TYPE_F32, NULL);
 
-    //
-    // If the overscans != NULL, then check the type of each image.
-    //
-    if (overscans != NULL) {
-        psListElem *tmpOverscan = (psListElem *) overscans->head;
-        while (NULL != tmpOverscan) {
-            psImage *myOverscanImage = (psImage *) tmpOverscan->data;
-            PS_ASSERT_IMAGE_TYPE(myOverscanImage, PS_TYPE_F32, NULL);
-            tmpOverscan = tmpOverscan->next;
-        }
-    }
-
-    if ((overscans == NULL) && (overScanAxis != PM_OVERSCAN_NONE)) {
-        psError(PS_ERR_UNKNOWN,true, "(overscans == NULL) && (overScanAxis != PM_OVERSCAN_NONE).  Returning in image\n");
-        return(in);
-    }
-
-    // Check for an unallowable pmFit.
-    if ((fit != PM_OVERSCAN_NONE) &&
-            (fit != PM_OVERSCAN_ROWS) &&
-            (fit != PM_OVERSCAN_COLUMNS) &&
-            (fit != PM_OVERSCAN_ALL)) {
-        psError(PS_ERR_UNKNOWN, true, "fit is unallowable (%d).  Returning in image.\n", fit);
-        return(in);
-    }
-    // Check for an unallowable pmOverscanAxis.
-    if ((overScanAxis != PM_OVERSCAN_NONE) &&
-            (overScanAxis != PM_OVERSCAN_ROWS) &&
-            (overScanAxis != PM_OVERSCAN_COLUMNS) &&
-            (overScanAxis != PM_OVERSCAN_ALL)) {
-        psError(PS_ERR_UNKNOWN, true, "overScanAxis is unallowable (%d).  Returning in image.\n", overScanAxis);
-        return(in);
-    }
-    psS32 i;
-    psS32 j;
-    psS32 numBins = 0;
-    static psVector *overscanVector = NULL;
-    psVector *tmpRow = NULL;
-    psVector *tmpCol = NULL;
-    psVector *myBin = NULL;
-    psVector *binVec = NULL;
-    psListElem *tmpOverscan = NULL;
-    double statValue;
-    psImage *myOverscanImage = NULL;
-    psPolynomial1D *myPoly = NULL;
-    psSpline1D *mySpline = NULL;
-    psS32 nBin;
-
-    //
-    //  Create a static stats data structure and determine the highest
-    //  priority stats option.
-    //
-    static psStats *myStats = NULL;
-    if (myStats == NULL) {
-        myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
-        p_psMemSetPersistent(myStats, true);
-    }
-    if (stat != NULL) {
-        myStats->options = GenNewStatOptions(stat);
-    }
-
-
-    if (overScanAxis == PM_OVERSCAN_NONE) {
-        if (fit != PM_FIT_NONE) {
-            psLogMsg(__func__, PS_LOG_WARN,
-                     "WARNING: pmSubtractBias.(): overScanAxis equals NONE, and fit does not equal NONE.  Proceeding to full fram subtraction.\n");
-        }
-
-        if (overscans != NULL) {
-            psLogMsg(__func__, PS_LOG_WARN,
-                     "WARNING: pmSubtractBias.(): overScanAxis equals NONE and overscans does not equal NULL.  Proceeding to full fram subtraction.\n");
-        }
-        return(SubtractFrame(in, bias));
-    }
-
-    if ((overScanAxis == PM_OVERSCAN_ALL) && (fit != PM_FIT_NONE)) {
-        psLogMsg(__func__, PS_LOG_WARN,
-                 "WARNING: pmSubtractBias.(): overScanAxis equals ALL, and fit does not equal NONE.  Proceeding with the rest of the module.\n");
-    }
-
-
-    //
-    // We subtract each overscan region from the image data.
-    // If we get here we know that overscans != NULL.
-    //
-
-    if (overScanAxis == PM_OVERSCAN_ALL) {
-        tmpOverscan = (psListElem *) overscans->head;
-        while (NULL != tmpOverscan) {
-            myOverscanImage = (psImage *) tmpOverscan->data;
-
-            PS_ASSERT_IMAGE_TYPE(myOverscanImage, PS_TYPE_F32, NULL);
-            psStats *rc = psImageStats(myStats, myOverscanImage, NULL, (psMaskType)0xffffffff);
-            if (rc == NULL) {
-                psError(PS_ERR_UNKNOWN, false, "psImageStats(): could not perform requested statistical operation.  Returning in image.\n");
+    pmCell *cell = in->parent;      // The parent cell
+    psRegion *trimsec = psMetadataLookupPtr(NULL, cell->concepts, "CELL.TRIMSEC"); // The trim region
+    psImage *image = psImageSubset(in->image, *trimsec); // The image corresponding to the trim region
+
+    // Overscan processing
+    if (overscanOpts) {
+        // Check for an unallowable pmFit.
+        if (overscanOpts->fitType != PM_FIT_NONE && overscanOpts->fitType != PM_FIT_POLY_ORD &&
+                overscanOpts->fitType != PM_FIT_POLY_CHEBY && overscanOpts->fitType != PM_FIT_SPLINE) {
+            psError(PS_ERR_UNKNOWN, true, "Invalid fit type (%d).  Returning original image.\n", overscanOpts->fitType);
+            return(in);
+        }
+
+        // Get the list of overscans
+        psList *overscanRegions = psMetadataLookupPtr(NULL, cell->concepts, "CELL.BIASSEC");
+        psList *overscans = psListAlloc(NULL); // List of the overscan images
+        psListIterator *iter = psListIteratorAlloc(overscanRegions, PS_LIST_HEAD, false); // Iterator
+        psRegion *biassec = NULL;       // A BIASSEC region from the list
+        while ((biassec = psListGetAndIncrement(iter))) {
+            psImage *overscan = psImageSubset(in->image, *biassec);
+            psListAdd(overscans, PS_LIST_TAIL, overscan);
+            psFree(overscan);
+        }
+        psFree(iter);
+
+        psStats *myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); // A new psStats, to avoid clobbering original
+        myStats->options = GenNewStatOptions(overscanOpts->stat);
+
+        // Reduce all overscan pixels to a single value
+        if (overscanOpts->single) {
+            psVector *pixels = psVectorAlloc(0, PS_TYPE_F32);
+            pixels->n = 0;
+            psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator
+            psImage *overscan = NULL;   // Overscan image from iterator
+            while ((overscan = psListGetAndIncrement(iter))) {
+                int index = pixels->n;  // Index
+                pixels = psVectorRealloc(pixels, pixels->n + overscan->numRows * overscan->numCols);
+                // XXX Reimplement with memcpy
+                for (int i = 0; i < overscan->numRows; i++) {
+                    for (int j = 0; j < overscan->numCols; j++) {
+                        pixels->data.F32[index++] = overscan->data.F32[i][j];
+                    }
+                }
+
+            }
+            psFree(iter);
+
+            (void)psVectorStats(myStats, pixels, NULL, NULL, 0);
+            double reduced = NAN;     // Result of statistics
+            if (! p_psGetStatValue(myStats, &reduced)) {
+                psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning input image.\n");
                 return(in);
             }
-            if (false == p_psGetStatValue(myStats, &statValue)) {
-                psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning in image.\n");
-                return(in);
+            (void)psBinaryOp(image, image, "-", psScalarAlloc((float)reduced, PS_TYPE_F32));
+        } else {
+
+            // We do the regular overscan subtraction
+            bool readRows = psMetadataLookupBool(NULL, cell->concepts, "CELL.READDIR"); // Read direction
+
+            if (readRows) {
+                // The read direction is rows
+                psArray *pixels = psArrayAlloc(image->numRows); // Array of vectors containing pixels
+                for (int i = 0; i < pixels->n; i++) {
+                    psVector *values = psVectorAlloc(0, PS_TYPE_F32);
+                    values->n = 0;
+                    pixels->data[i] = values;
+                }
+
+                // Pull the pixels out into the vectors
+                psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator
+                psImage *overscan = NULL; // Overscan image from iterator
+                while ((overscan = psListGetAndIncrement(iter))) {
+                    int diff = image->row0 - overscan->row0; // Offset between the two regions
+                    for (int i = MAX(0,diff); i < MIN(image->numRows, overscan->numRows + diff); i++) {
+                        // i is row on overscan
+                        // XXX Reimplement with memcpy
+                        psVector *values = pixels->data[i];
+                        int index = values->n; // Index in the vector
+                        values = psVectorRealloc(values, values->n + overscan->numCols);
+                        for (int j = 0; j < overscan->numCols; j++) {
+                            values->data.F32[index++] = overscan->data.F32[i][j];
+                        }
+                        values->n += overscan->numCols;
+                        pixels->data[i] = values; // Update the pointer in case it's moved
+                    }
+                }
+
+                // Reduce the overscans
+                psVector *reduced = overscanVector(overscanOpts, pixels, myStats);
+                if (! reduced) {
+                    return in;
+                }
+
+                // Subtract row by row
+                for (int i = 0; i < image->numRows; i++) {
+                    for (int j = 0; j < image->numCols; j++) {
+                        image->data.F32[i][j] -= reduced->data.F32[i];
+                    }
+                }
+                psFree(reduced);
+
+            } else {
+                // The read direction is columns
+                psArray *pixels = psArrayAlloc(image->numCols); // Array of vectors containing pixels
+                for (int i = 0; i < pixels->n; i++) {
+                    psVector *values = psVectorAlloc(0, PS_TYPE_F32);
+                    values->n = 0;
+                    pixels->data[i] = values;
+                }
+
+                // Pull the pixels out into the vectors
+                psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator
+                psImage *overscan = NULL; // Overscan image from iterator
+                while ((overscan = psListGetAndIncrement(iter))) {
+                    int diff = image->col0 - overscan->col0; // Offset between the two regions
+                    for (int i = MAX(0,diff); i < MIN(image->numCols, overscan->numCols + diff); i++) {
+                        // i is column on overscan
+                        // XXX Reimplement with memcpy
+                        psVector *values = pixels->data[i];
+                        int index = values->n; // Index in the vector
+                        values = psVectorRealloc(values, values->n + overscan->numRows);
+                        for (int j = 0; j < overscan->numRows; j++) {
+                            values->data.F32[index++] = overscan->data.F32[i][j];
+                        }
+                        values->n += overscan->numRows;
+                        pixels->data[i] = values; // Update the pointer in case it's moved
+                    }
+                }
+
+                // Reduce the overscans
+                psVector *reduced = overscanVector(overscanOpts, pixels, myStats);
+                if (! reduced) {
+                    return in;
+                }
+
+                // Subtract column by column
+                for (int i = 0; i < image->numCols; i++) {
+                    for (int j = 0; j < image->numRows; j++) {
+                        image->data.F32[j][i] -= reduced->data.F32[i];
+                    }
+                }
+                psFree(reduced);
             }
-            ImageSubtractScalar(in->image, statValue);
-
-            tmpOverscan = tmpOverscan->next;
-        }
-        return(in);
-    }
-
-    // This check is redundant with above code.
-    if (!((overScanAxis == PM_OVERSCAN_ROWS) || (overScanAxis == PM_OVERSCAN_COLUMNS))) {
-        psError(PS_ERR_UNKNOWN, true, "overScanAxis is unallowable (%d).\nReturning in image.\n", overScanAxis);
-        return(in);
-    }
-
-    tmpOverscan = (psListElem *) overscans->head;
-    while (NULL != tmpOverscan) {
-        //        PS_IMAGE_PRINT_F32_HIDEF(in->image);
-        myOverscanImage = (psImage *) tmpOverscan->data;
-
-        if (overScanAxis == PM_OVERSCAN_ROWS) {
-            if (myOverscanImage->numCols != (in->image)->numCols) {
-                psLogMsg(__func__, PS_LOG_WARN,
-                         "WARNING: pmSubtractBias.(): overscan image has %d columns, input image has %d columns\n",
-                         myOverscanImage->numCols, in->image->numCols);
-            }
-
-            // We create a row vector and subtract this vector from image.
-            // XXX: Is there a better way to extract a psVector from a psImage without
-            // having to copy every element in that vector?
-            overscanVector = psVectorAlloc(myOverscanImage->numCols, PS_TYPE_F32);
-            for (i=0;i<overscanVector->n;i++) {
-                overscanVector->data.F32[i] = 0.0;
-            }
-            tmpRow = psVectorAlloc(myOverscanImage->numRows, PS_TYPE_F32);
-
-            // For each column of the input image, loop through every row,
-            // collect the pixel in that row, then performed the specified
-            // statistical op on those pixels.  Store this in overscanVector.
-            for (i=0;i<myOverscanImage->numCols;i++) {
-                for (j=0;j<myOverscanImage->numRows;j++) {
-                    tmpRow->data.F32[j] = myOverscanImage->data.F32[j][i];
-                }
-                psStats *rc = psVectorStats(myStats, tmpRow, NULL, NULL, 0);
-                if (rc == NULL) {
-                    psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation.  Returning in image.\n");
-                    return(in);
-                }
-                if (false ==  p_psGetStatValue(rc, &statValue)) {
-                    psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning in image.\n");
-                    return(in);
-                }
-                overscanVector->data.F32[i] = statValue;
-            }
-            psFree(tmpRow);
-
-            // Scale the overscan vector to the size of the input image.
-            if (overscanVector->n != in->image->numCols) {
-                if ((fit == PM_FIT_POLYNOMIAL) || (fit == PM_FIT_SPLINE)) {
-                    psVector *newVec = ScaleOverscanVector(overscanVector,
-                                                           in->image->numCols,
-                                                           fitSpec, fit);
-                    if (newVec == NULL) {
-                        psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector(): could not scale the overscan vector.  Returning in image.\n");
-                        return(in);
-                    }
-                    psFree(overscanVector);
-                    overscanVector = newVec;
-                } else {
-                    psError(PS_ERR_UNKNOWN, true, "Don't know how to scale the overscan vector.  Set fit to PM_FIT_SPLINE or PM_FIT_POLYNOMIAL.  Returning in image.\n");
-                    psFree(overscanVector);
-                    return(in);
-                }
-            }
-        }
-
-        if (overScanAxis == PM_OVERSCAN_COLUMNS) {
-            if (myOverscanImage->numRows != (in->image)->numRows) {
-                psLogMsg(__func__, PS_LOG_WARN,
-                         "WARNING: pmSubtractBias.(): overscan image has %d rows, input image has %d rows\n",
-                         myOverscanImage->numRows, in->image->numRows);
-            }
-
-            // We create a column vector and subtract this vector from image.
-            overscanVector = psVectorAlloc(myOverscanImage->numRows, PS_TYPE_F32);
-            for (i=0;i<overscanVector->n;i++) {
-                overscanVector->data.F32[i] = 0.0;
-            }
-            tmpCol = psVectorAlloc(myOverscanImage->numCols, PS_TYPE_F32);
-
-            // For each row of the input image, loop through every column,
-            // collect the pixel in that row, then performed the specified
-            // statistical op on those pixels.  Store this in overscanVector.
-            for (i=0;i<myOverscanImage->numRows;i++) {
-                for (j=0;j<myOverscanImage->numCols;j++) {
-                    tmpCol->data.F32[j] = myOverscanImage->data.F32[i][j];
-                }
-                psStats *rc = psVectorStats(myStats, tmpCol, NULL, NULL, 0);
-                if (rc == NULL) {
-                    psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation.  Returning in image.\n");
-                    return(in);
-                }
-                if (false ==  p_psGetStatValue(rc, &statValue)) {
-                    psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning in image.\n");
-                    return(in);
-                }
-                overscanVector->data.F32[i] = statValue;
-            }
-            psFree(tmpCol);
-
-            // Scale the overscan vector to the size of the input image.
-            if (overscanVector->n != in->image->numRows) {
-                if ((fit == PM_FIT_POLYNOMIAL) || (fit == PM_FIT_SPLINE)) {
-                    psVector *newVec = ScaleOverscanVector(overscanVector,
-                                                           in->image->numRows,
-                                                           fitSpec, fit);
-                    if (newVec == NULL) {
-                        psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector(): could not scale the overscan vector.  Returning in image.\n");
-                        return(in);
-                    }
-                    psFree(overscanVector);
-                    overscanVector = newVec;
-                } else {
-                    psError(PS_ERR_UNKNOWN, true, "Don't know how to scale the overscan vector.  Set fit to PM_FIT_SPLINE or PM_FIT_POLYNOMIAL.  Returning in image.\n");
-                    psFree(overscanVector);
-                    return(in);
-                }
-            }
-        }
-
-        //
-        // Re-bin the overscan vector (change its length).
-        //
-        // Only if nBinOrig > 1.
-        if ((nBinOrig > 1) && (nBinOrig < overscanVector->n)) {
-            numBins = 1+((overscanVector->n)/nBinOrig);
-            myBin = psVectorAlloc(numBins, PS_TYPE_F32);
-            binVec = psVectorAlloc(nBinOrig, PS_TYPE_F32);
-
-            for (i=0;i<numBins;i++) {
-                for(j=0;j<nBinOrig;j++) {
-                    if (overscanVector->n > ((i*nBinOrig)+j)) {
-                        binVec->data.F32[j] = overscanVector->data.F32[(i*nBinOrig)+j];
-                    } else {
-                        // XXX: we get here if nBinOrig does not evenly divide
-                        // the overscanVector vector.  This is the last bin.  Should
-                        // we change the binVec->n to acknowledge that?
-                        binVec->n = j;
-                    }
-                }
-                psStats *rc = psVectorStats(myStats, binVec, NULL, NULL, 0);
-                if (rc == NULL) {
-                    psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation.  Returning in image.\n");
-                    return(in);
-                }
-                if (false ==  p_psGetStatValue(rc, &statValue)) {
-                    psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation.  Returning in image.\n");
-                    return(in);
-                }
-                myBin->data.F32[i] = statValue;
-            }
-
-            // Change the effective size of overscanVector.
-            overscanVector->n = numBins;
-            for (i=0;i<numBins;i++) {
-                overscanVector->data.F32[i] = myBin->data.F32[i];
-            }
-            psFree(binVec);
-            psFree(myBin);
-            nBin = nBinOrig;
-        } else {
-            nBin = 1;
-        }
-
-        // At this point the number of data points in overscanVector should be
-        // equal to the number of rows/columns (whatever is appropriate) in the
-        // image divided by numBins.
-        //
-
-
-        //
-        // This doesn't seem right.  The only way to do a spline fit is if,
-        // by SDRS requirements, fitSpec is not-NULL>  But in order for it
-        // to be non-NULL, someone must have called psSpline1DAlloc() with
-        // the min, max, and number of splines.
-        //
-        if (!((fitSpec == NULL) || (fit == PM_FIT_NONE))) {
-            //
-            // Fit a polynomial or spline to the overscan vector.
-            //
-            if (fit == PM_FIT_POLYNOMIAL) {
-                myPoly = (psPolynomial1D *) fitSpec;
-                myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, overscanVector, NULL, NULL);
-                if (myPoly == NULL) {
-                    psError(PS_ERR_UNKNOWN, false, "(3) Could not fit a polynomial to overscan vector.  Returning in image.\n");
-                    psFree(overscanVector);
-                    return(in);
-                }
-            } else if (fit == PM_FIT_SPLINE) {
-                // XXX: This makes no sense
-                // XXX: must free mySpline?
-                mySpline = (psSpline1D *) fitSpec;
-                mySpline = psVectorFitSpline1D(NULL, overscanVector);
-                if (mySpline == NULL) {
-                    psError(PS_ERR_UNKNOWN, false, "Could not fit a spline to overscan vector.  Returning in image.\n");
-                    psFree(overscanVector);
-                    return(in);
-                }
-            }
-
-            //
-            // Subtract fitted overscan vector row-wise from the image.
-            //
-            if (overScanAxis == PM_OVERSCAN_ROWS) {
-                for (i=0;i<(in->image)->numCols;i++) {
-                    psF32 tmpF32 = 0.0;
-                    if (fit == PM_FIT_POLYNOMIAL) {
-                        tmpF32 = psPolynomial1DEval(myPoly, ((psF32) i) / ((psF32) nBin));
-                    } else if (fit == PM_FIT_SPLINE) {
-                        tmpF32 = psSpline1DEval(mySpline, ((psF32) i) / ((psF32) nBin));
-                    }
-                    for (j=0;j<(in->image)->numRows;j++) {
-                        (in->image)->data.F32[j][i]-= tmpF32;
-                    }
-                }
-            }
-
-            //
-            // Subtract fitted overscan vector column-wise from the image.
-            //
-            if (overScanAxis == PM_OVERSCAN_COLUMNS) {
-                for (i=0;i<(in->image)->numRows;i++) {
-                    psF32 tmpF32 = 0.0;
-                    if (fit == PM_FIT_POLYNOMIAL) {
-                        tmpF32 = psPolynomial1DEval(myPoly, ((psF32) i) / ((psF32) nBin));
-                    } else if (fit == PM_FIT_SPLINE) {
-                        tmpF32 = psSpline1DEval(mySpline, ((psF32) i) / ((psF32) nBin));
-                    }
-
-                    for (j=0;j<(in->image)->numCols;j++) {
-                        (in->image)->data.F32[i][j]-= tmpF32;
-                    }
-                }
-            }
-        } else {
-            //
-            // If we get here, then no polynomials were fit to the overscan
-            // vector.  We simply subtract it, taking into account binning,
-            // from the image.
-            //
-
-            //
-            // Subtract overscan vector row-wise from the image.
-            //
-            if (overScanAxis == PM_OVERSCAN_ROWS) {
-                for (i=0;i<(in->image)->numCols;i++) {
-                    for (j=0;j<(in->image)->numRows;j++) {
-                        (in->image)->data.F32[j][i]-= overscanVector->data.F32[i/nBin];
-                    }
-                }
-            }
-
-            //
-            // Subtract overscan vector column-wise from the image.
-            //
-            if (overScanAxis == PM_OVERSCAN_COLUMNS) {
-                for (i=0;i<(in->image)->numRows;i++) {
-                    for (j=0;j<(in->image)->numCols;j++) {
-                        (in->image)->data.F32[i][j]-= overscanVector->data.F32[i/nBin];
-                    }
-                }
-            }
-        }
-
-        psFree(overscanVector);
-
-        tmpOverscan = tmpOverscan->next;
-    }
-
-    psTrace(".psModule.pmSubtracBias.pmSubtractBias", 4,
-            "---- pmSubtractBias() exit ----\n");
-
-    if (bias != NULL) {
-        return(SubtractFrame(in, bias));
-    }
-    return(in);
-}
-
-
+        }
+        psFree(myStats);
+        psFree(overscans);
+    } // End of overscan subtraction
+
+    // Bias frame subtraction
+    if (bias) {
+        SubtractFrame(in, bias, 1.0);
+    }
+
+    if (dark) {
+        // Get the scaling
+        float inTime = psMetadataLookupF32(NULL, in->parent->concepts, "CELL.DARKTIME");
+        float darkTime = psMetadataLookupF32(NULL, dark->parent->concepts, "CELL.DARKTIME");
+        SubtractFrame(in, dark, inTime/darkTime);
+    }
+
+    return in;
+}
+
+
Index: /branches/eam_rel9_p0/psModules/src/imsubtract/pmSubtractBias.h
===================================================================
--- /branches/eam_rel9_p0/psModules/src/imsubtract/pmSubtractBias.h	(revision 6061)
+++ /branches/eam_rel9_p0/psModules/src/imsubtract/pmSubtractBias.h	(revision 6062)
@@ -11,6 +11,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.4.8.1.2.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-12-31 04:35:36 $
+ *  @version $Revision: 1.4.8.1.2.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-20 02:38:28 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -29,5 +29,5 @@
 #include "pslib.h"
 
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 
 typedef enum {
@@ -40,9 +40,29 @@
 
 typedef enum {
-    PM_FIT_NONE,                              ///< No fit
-    PM_FIT_POLYNOMIAL,                        ///< Fit polynomial
-    PM_FIT_SPLINE                             ///< Fit cubic splines
+    PM_FIT_NONE,                        ///< No fit
+    PM_FIT_POLY_ORD,                    ///< Fit ordinary polynomial
+    PM_FIT_POLY_CHEBY,                  ///< Fit Chebyshev polynomial
+    PM_FIT_SPLINE                       ///< Fit cubic splines
 } pmFit;
 
+typedef struct
+{
+    // Inputs
+    bool single;                // Reduce all overscan regions to a single value?
+    pmFit fitType;              // Type of fit to overscan
+    unsigned int order;         // Order of polynomial, or number of spline pieces
+    psStats *stat;              // Statistic to use when reducing the minor direction
+    // Outputs
+    psPolynomial1D *poly;       // Result of polynomial fit
+    psSpline1D *spline;         // Result of spline fit
+}
+pmOverscanOptions;
+
+pmOverscanOptions *pmOverscanOptionsAlloc(bool single, pmFit fitType, unsigned int order, psStats *stat);
+
+pmReadout *pmSubtractBias(pmReadout *in, pmOverscanOptions *overscanOpts,
+                          const pmReadout *bias, const pmReadout *dark);
+
+#if 0
 pmReadout *pmSubtractBias(pmReadout *in,                ///< The input pmReadout image
                           void *fitSpec,                ///< A polynomial or spline, defining the fit type.
@@ -53,4 +73,5 @@
                           pmFit fit,                    ///< PM_FIT_SPLINE, PM_FIT_POLYNOMIAL, or PM_FIT_NONE
                           const pmReadout *bias);       ///< A possibly NULL bias pmReadout which is to be subtracted
+#endif
 
 #endif
Index: /branches/eam_rel9_p0/psModules/src/imsubtract/pmSubtractSky.h
===================================================================
--- /branches/eam_rel9_p0/psModules/src/imsubtract/pmSubtractSky.h	(revision 6061)
+++ /branches/eam_rel9_p0/psModules/src/imsubtract/pmSubtractSky.h	(revision 6062)
@@ -6,6 +6,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-09-28 20:43:52 $
+ *  @version $Revision: 1.1.18.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-20 02:38:28 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -23,5 +23,5 @@
 #include<math.h>
 #include "pslib.h"
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 
 // XXX: this is pmFit in pmSubtractBias.c, named psFit here.
Index: /branches/eam_rel9_p0/psModules/src/objects/pmObjects.h
===================================================================
--- /branches/eam_rel9_p0/psModules/src/objects/pmObjects.h	(revision 6061)
+++ /branches/eam_rel9_p0/psModules/src/objects/pmObjects.h	(revision 6062)
@@ -4,5 +4,5 @@
  * images is one of the critical tasks of the IPP or any astronomical software
  * system. This file will define structures and functions related to the task
- * of source detection and measurement. The elements defined in this section 
+ * of source detection and measurement. The elements defined in this section
  * are generally low-level components which can be connected together to
  * construct a complete object measurement suite.
@@ -10,6 +10,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.4.4.3 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-01-15 18:21:19 $
+ *  @version $Revision: 1.4.4.4 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-20 02:38:28 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -27,5 +27,5 @@
 #include <math.h>
 #include "pslib.h"
-#include "pmAstrometry.h"
+#include "pmFPA.h"
 /**
  * In the object analysis process, we will use specific mask values to mark the
@@ -45,5 +45,5 @@
 
 /** pmPeakType
- * 
+ *
  *  A peak pixel may have several features which may be determined when the
  *  peak is found or measured. These are specified by the pmPeakType enum.
@@ -52,5 +52,5 @@
  *  edge. The PM_PEAK_FLAT represents a peak pixel which has more than a specific
  *  number of neighbors at the same value, within some tolarence:
- * 
+ *
  */
 typedef enum {
@@ -63,9 +63,9 @@
 
 /** pmPeak data structure
- *  
+ *
  *  A source has the capacity for several types of measurements. The
  *  simplest measurement of a source is the location and flux of the peak pixel
  *  associated with the source:
- *  
+ *
  */
 typedef struct
@@ -80,9 +80,9 @@
 
 /** pmMoments data structure
- *  
+ *
  * One of the simplest measurements which can be made quickly for an object
  * are the object moments. We specify a structure to carry the moment information
  * for a specific source:
- *  
+ *
  */
 typedef struct
@@ -103,9 +103,9 @@
 
 /** pmPSFClump data structure
- * 
+ *
  * A collection of object moment measurements can be used to determine
  * approximate object classes. The key to this analysis is the location and
  * statistics (in the second-moment plane,
- *  
+ *
  */
 typedef struct
@@ -130,5 +130,5 @@
 
 /** pmModel data structure
- * 
+ *
  * Every source may have two types of models: a PSF model and a EXT (extended-source)
  * model. The PSF model represents the best fit of the image PSF to the specific
@@ -136,5 +136,5 @@
  * object by the PSF, not by the fit. The EXT model represents the best fit of
  * the given model to the object, with all shape parameters floating in the fit.
- *  
+ *
  */
 typedef struct
@@ -152,12 +152,12 @@
 
 /** pmSourceType enumeration
- * 
+ *
  * A given source may be identified as most-likely to be one of several source
  * types. The pmSource entry pmSourceType defines the current best-guess for this
  * source.
- * 
+ *
  * XXX: The values given below are currently illustrative and will require
  * some modification as the source classification code is developed. (TBD)
- * 
+ *
  */
 typedef enum {
@@ -186,9 +186,9 @@
 
 /** pmSource data structure
- *  
+ *
  *  This source has the capacity for several types of measurements. The
  *  simplest measurement of a source is the location and flux of the peak pixel
  *  associated with the source:
- *  
+ *
  */
 typedef struct
@@ -224,5 +224,5 @@
 
 /** pmMomentsAlloc()
- * 
+ *
  */
 pmMoments *pmMomentsAlloc();
@@ -230,5 +230,5 @@
 
 /** pmModelAlloc()
- * 
+ *
  */
 pmModel *pmModelAlloc(pmModelType type);
@@ -236,5 +236,5 @@
 
 /** pmSourceAlloc()
- * 
+ *
  */
 pmSource  *pmSourceAlloc();
@@ -242,5 +242,5 @@
 
 /** pmFindVectorPeaks()
- * 
+ *
  * Find all local peaks in the given vector above the given threshold. A peak
  * is defined as any element with a value greater than its two neighbors and with
@@ -256,5 +256,5 @@
  * a vector containing the coordinates (element number) of the detected peaks
  * (type psU32).
- * 
+ *
  */
 psVector *pmFindVectorPeaks(
@@ -265,5 +265,5 @@
 
 /** pmFindImagePeaks()
- * 
+ *
  * Find all local peaks in the given image above the given threshold. This
  * function should find all row peaks using pmFindVectorPeaks, then test each row
@@ -276,5 +276,5 @@
  * peaks at the (+x,+y) corner. When selecting the peaks, their type must also be
  * set. The result of this function is an array of pmPeak entries.
- * 
+ *
  */
 psArray *pmFindImagePeaks(
@@ -285,8 +285,8 @@
 
 /** pmCullPeaks()
- * 
+ *
  * Eliminate peaks from the psList that have a peak value above the given
  * maximum, or fall outside the valid region.
- * 
+ *
  */
 psList *pmCullPeaks(
@@ -298,5 +298,5 @@
 
 /** pmPeaksSubset()
- * 
+ *
  * Create a new peaks array, removing certain types of peaks from the input
  * array of peaks based on the given criteria. Peaks should be eliminated if they
@@ -304,5 +304,5 @@
  * the valid region.  The result of the function is a new array with a reduced
  * number of peaks.
- * 
+ *
  */
 psArray *pmPeaksSubset(
@@ -314,5 +314,5 @@
 
 /** pmSourceDefinePixels()
- * 
+ *
  * Define psImage subarrays for the source located at coordinates x,y on the
  * image set defined by readout. The pixels defined by this operation consist of
@@ -326,7 +326,7 @@
  * example). This function should be used to define a region of interest around a
  * source, including both source and sky pixels.
- * 
+ *
  * XXX: must code this.
- * 
+ *
  */
 // XXX: Uncommenting the pmReadout causes compile errors.
@@ -341,5 +341,5 @@
 
 /** pmSourceLocalSky()
- * 
+ *
  * Measure the local sky in the vicinity of the given source. The Radius
  * defines the square aperture in which the moments will be measured. This
@@ -351,5 +351,5 @@
  * This function allocates the pmMoments structure. The resulting sky is used to
  * set the value of the pmMoments.sky element of the provided pmSource structure.
- * 
+ *
  */
 bool pmSourceLocalSky(
@@ -361,5 +361,5 @@
 
 /** pmSourceMoments()
- * 
+ *
  * Measure source moments for the given source, using the value of
  * source.moments.sky provided as the local background value and the peak
@@ -368,5 +368,5 @@
  * are measured within the given circular radius of the source.peak coordinates.
  * The return value indicates the success (TRUE) of the operation.
- * 
+ *
  */
 bool pmSourceMoments(
@@ -377,5 +377,5 @@
 
 /** pmSourcePSFClump()
- * 
+ *
  * We use the source moments to make an initial, approximate source
  * classification, and as part of the information needed to build a PSF model for
@@ -386,5 +386,5 @@
  * similar. The function returns a pmPSFClump structure, representing the
  * centroid and size of the clump in the sigma_x, sigma_y second-moment plane.
- * 
+ *
  * The goal is to identify and characterize the stellar clump within the
  * sigma_x, sigma_y second-moment plane.  To do this, an image is constructed to
@@ -397,10 +397,10 @@
  *  * used to calculate the median and standard deviation of the sigma_x, sigma_y
  * values. These resulting values are returned via the pmPSFClump structure.
- * 
+ *
  * The return value indicates the success (TRUE) of the operation.
- * 
+ *
  * XXX: Limit the S/N of the candidate sources (part of Metadata)? (TBD).
  * XXX: Save the clump parameters on the Metadata (TBD)
- * 
+ *
  */
 pmPSFClump pmSourcePSFClump(
@@ -411,5 +411,5 @@
 
 /** pmSourceRoughClass()
- * 
+ *
  * Based on the specified data values, make a guess at the source
  * classification. The sources are provides as a psArray of pmSource entries.
@@ -418,5 +418,5 @@
  * can be extracted from the metadata using the given keywords. Except as noted,
  * the data type for these parameters are psF32.
- * 
+ *
  */
 bool pmSourceRoughClass(
@@ -428,5 +428,5 @@
 
 /** pmSourceModelGuess()
- * 
+ *
  * Convert available data to an initial guess for the given model. This
  * function allocates a pmModel entry for the pmSource structure based on the
@@ -434,5 +434,5 @@
  * are specified for each model below. The guess values are placed in the model
  * parameters. The function returns TRUE on success or FALSE on failure.
- * 
+ *
  */
 pmModel *pmSourceModelGuess(
@@ -443,7 +443,7 @@
 
 /** pmContourType
- * 
+ *
  * Only one type is defined at present.
- * 
+ *
  */
 typedef enum {
@@ -455,5 +455,5 @@
 
 /** pmSourceContour()
- * 
+ *
  * Find points in a contour for the given source at the given level. If type
  * is PM_CONTOUR_CRUDE, the contour is found by starting at the source peak,
@@ -465,5 +465,5 @@
  * This function may be used as part of the model guess inputs.  Other contour
  * types may be specified in the future for more refined contours (TBD)
- * 
+ *
  */
 psArray *pmSourceContour(
@@ -476,5 +476,5 @@
 
 /** pmSourceFitModel()
- * 
+ *
  * Fit the requested model to the specified source. The starting guess for the
  * model is given by the input source.model parameter values. The pixels of
@@ -482,5 +482,5 @@
  * function calls psMinimizeLMChi2() on the image data. The function returns TRUE
  * on success or FALSE on failure.
- * 
+ *
  */
 bool pmSourceFitModel(
@@ -492,12 +492,12 @@
 
 /** pmModelFitStatus()
- * 
+ *
  * This function wraps the call to the model-specific function returned by
  * pmModelFitStatusFunc_GetFunction.  The model-specific function examines the
  * model parameters, parameter errors, Chisq, S/N, and other parameters available
  * from model to decide if the particular fit was successful or not.
- * 
+ *
  * XXX: Must code this.
- * 
+ *
  */
 bool pmModelFitStatus(
@@ -507,5 +507,5 @@
 
 /** pmSourceAddModel()
- * 
+ *
  * Add the given source model flux to/from the provided image. The boolean
  * option center selects if the source is re-centered to the image center or if
@@ -514,5 +514,5 @@
  * is at most the pixel range specified by the source.pixels image. The success
  * status is returned.
- * 
+ *
  */
 bool pmSourceAddModel(
@@ -526,5 +526,5 @@
 
 /** pmSourceSubModel()
- * 
+ *
  * Subtract the given source model flux to/from the provided image. The
  * boolean option center selects if the source is re-centered to the image center
@@ -533,5 +533,5 @@
  * image is at most the pixel range specified by the source.pixels image. The
  * success status is returned.
- * 
+ *
  */
 bool pmSourceSubModel(
@@ -545,5 +545,5 @@
 
 /**
- * 
+ *
  * The function returns both the magnitude of the fit, defined as -2.5log(flux),
  * where the flux is integrated under the model, theoretically from a radius of 0
@@ -552,7 +552,7 @@
  * not excluded by the aperture mask. The model flux is calculated by calling the
  * model-specific function provided by pmModelFlux_GetFunction.
- * 
+ *
  * XXX: must code this.
- * 
+ *
  */
 bool pmSourcePhotometry(
@@ -566,9 +566,9 @@
 
 /**
- * 
+ *
  * This function converts the source classification into the closest available
  * approximation to the Dophot classification scheme:
  * XXX EAM : fix this to use current source classification scheme
- * 
+ *
  * PM_SOURCE_DEFECT: 8
  * PM_SOURCE_SATURATED: 8
@@ -585,5 +585,5 @@
  * PM_SOURCE_POOR_FIT_GAL: 2
  * PM_SOURCE_OTHER: ?
- * 
+ *
  */
 int pmSourceDophotType(
@@ -593,11 +593,11 @@
 
 /** pmSourceSextractType()
- * 
+ *
  * This function converts the source classification into the closest available
  * approximation to the Sextractor classification scheme. the correspondence is
  * not yet defined (TBD) .
- * 
+ *
  * XXX: Must code this.
- * 
+ *
  */
 int pmSourceSextractType(
@@ -606,5 +606,5 @@
 
 /** pmSourceFitModel_v5()
- * 
+ *
  * Fit the requested model to the specified source. The starting guess for the
  * model is given by the input source.model parameter values. The pixels of
@@ -612,5 +612,5 @@
  * function calls psMinimizeLMChi2() on the image data. The function returns TRUE
  * on success or FALSE on failure.
- * 
+ *
  */
 bool pmSourceFitModel_v5(
@@ -622,5 +622,5 @@
 
 /** pmSourceFitModel_v7()
- * 
+ *
  * Fit the requested model to the specified source. The starting guess for the
  * model is given by the input source.model parameter values. The pixels of
@@ -628,5 +628,5 @@
  * function calls psMinimizeLMChi2() on the image data. The function returns TRUE
  * on success or FALSE on failure.
- * 
+ *
  */
 bool pmSourceFitModel_v7(
@@ -638,7 +638,7 @@
 
 /** pmSourcePhotometry()
- * 
+ *
  * XXX: Need descriptions
- * 
+ *
  */
 bool pmSourcePhotometry(