Index: trunk/psphot/src/psphotStackImageLoop.c
===================================================================
--- trunk/psphot/src/psphotStackImageLoop.c	(revision 27848)
+++ trunk/psphot/src/psphotStackImageLoop.c	(revision 28013)
@@ -1,9 +1,14 @@
 # include "psphotStandAlone.h"
-
-# define ESCAPE(MESSAGE) { \
-  psError(PSPHOT_ERR_DATA, false, MESSAGE); \
-  psFree (view); \
-  return false; \
-}
+#define WCS_NONLIN_TOL 0.001            // Non-linear tolerance for header WCS
+
+# define ESCAPE(MESSAGE) {				\
+	psError(PSPHOT_ERR_DATA, false, MESSAGE);	\
+	psFree (view);					\
+	return false;					\
+    }
+
+bool GetAstrometryFPA (pmConfig *config, pmFPAview *view);
+bool SetAstrometryFPA (pmConfig *config, pmFPAview *view, bool bilevelAstrometry);
+bool GetAstrometryChip (pmConfig *config, pmFPAview *view, bool bilevelAstrometry);
 
 bool psphotStackImageLoop (pmConfig *config) {
@@ -14,6 +19,9 @@
     pmReadout *readout;
 
-    pmFPAfile *input = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT");
-    if (!status) {
+    pmFPAfile *inputRaw = psMetadataLookupPtr (&status, config->files, "PSPHOT.STACK.INPUT.RAW");
+    pmFPAfile *inputCnv = psMetadataLookupPtr (&status, config->files, "PSPHOT.STACK.INPUT.CNV");
+    pmFPAfile *input = inputRaw ? inputRaw : inputCnv;
+
+    if (!input) {
         psError(PSPHOT_ERR_PROG, false, "Can't find input data!");
         return false;
@@ -24,4 +32,6 @@
     // XXX for now, just load the full set of images up front
     if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for fpa in psphot.");
+
+    bool bilevelAstrometry = GetAstrometryFPA(config, view);
 
     // for psphot, we force data to be read at the chip level
@@ -41,6 +51,5 @@
                 if (! readout->data_exists) { continue; }
 
-# if (0)		
-		// uncomment to generate matched psfs
+		// PSF matching
 		if (!psphotStackMatchPSFs (config, view)) {
                     psError(psErrorCodeLast(), false, "failure in psphotStackMatchPSFs for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout);
@@ -48,5 +57,4 @@
                     return false;
 		}
-# endif
 
 		// XXX for now, we assume there is only a single chip in the PHU:
@@ -69,7 +77,13 @@
 	    }
 	}
+
+	GetAstrometryChip(config, view, bilevelAstrometry);
+
 	// save output which is saved at the chip level
 	if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed output for Chip in psphot.");
     }
+
+    SetAstrometryFPA(config, view, bilevelAstrometry);
+    
     // save output which is saved at the fpa level
     if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed ouput for FPA in psphot.");
@@ -87,2 +101,151 @@
 */
 
+# define UPDATE_HEADER 0
+
+bool GetAstrometryFPA (pmConfig *config, pmFPAview *view) {
+
+    bool status = false;
+
+    bool foundAstrometry = false;
+    bool bilevelAstrometry = false;
+
+    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
+    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
+
+    // loop over the available readouts
+    for (int i = 0; i < num; i++) {
+
+	// find the currently selected readout
+	pmFPAfile *output = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.OUTPUT.IMAGE", i); // File of interest
+	psAssert (output, "missing file?");
+
+	pmFPAfile *inputRaw = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", i); // File of interest
+	pmFPAfile *inputCnv = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.CNV", i); // File of interest
+	pmFPAfile *input = inputRaw ? inputRaw : inputCnv;
+	psAssert (input, "missing input file");
+
+	// find the FPA phu
+	pmHDU *phu = pmFPAviewThisPHU(view, input->fpa);
+	if (!phu) {
+	    psWarning("Unable to read bilevel mosaic astrometry for input FPA entry %d", i);
+	    continue;
+	}
+
+	char *ctype = psMetadataLookupStr(NULL, phu->header, "CTYPE1");
+	if (!ctype) {
+	    psWarning("Error in WCS keywords for input FPA entry %d", i);
+	    continue;
+	}
+
+	if (!foundAstrometry) {
+	    bilevelAstrometry = !strcmp (&ctype[4], "-DIS");
+	} else {
+	    if (bilevelAstrometry != !strcmp (&ctype[4], "-DIS")) {
+		psAbort("astrometry type mis-match");
+	    }
+	}
+
+	if (bilevelAstrometry) {
+	    // update the output structures
+	    if (!pmAstromReadBilevelMosaic(output->fpa, phu->header)) {
+		psWarning("Unable to read bilevel mosaic astrometry for input FPA.");
+	    }
+	}
+    }
+    return bilevelAstrometry;
+}
+
+bool GetAstrometryChip (pmConfig *config, pmFPAview *view, bool bilevelAstrometry) {
+
+    bool status = false;
+
+    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
+    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
+
+    // loop over the available readouts
+    for (int i = 0; i < num; i++) {
+
+	// find the currently selected readout
+	pmFPAfile *output = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.OUTPUT.IMAGE", i); // File of interest
+	psAssert (output, "missing file?");
+
+	pmFPAfile *inputRaw = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.RAW", i); // File of interest
+	pmFPAfile *inputCnv = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.INPUT.CNV", i); // File of interest
+	pmFPAfile *input = inputRaw ? inputRaw : inputCnv;
+	psAssert (input, "missing input file");
+
+	// Need to update the output for astrometry.  Read WCS data from the corresponding
+        // header and save in the output fpa
+        pmHDU *inHDU = pmFPAviewThisHDU (view, input->fpa);
+	pmChip *outChip = pmFPAviewThisChip(view, output->fpa); ///< Chip in the output
+
+# if (UPDATE_HEADER)
+	pmHDU *outHDU = pmFPAviewThisHDU (view, output->fpa);
+	if (!outHDU) {
+	    pmFPAAddSourceFromView(output->fpa, "name", view, output->format);
+	    outHDU = pmFPAviewThisHDU (view, output->fpa);
+	    psAssert (outHDU, "failed to make HDU");
+	}
+# endif
+
+        if (bilevelAstrometry) {
+            if (!pmAstromReadBilevelChip (outChip, inHDU->header)) {
+                psWarning("Unable to read bilevel chip astrometry for input FPA.");
+		continue;
+	    }
+# if (UPDATE_HEADER)
+	    if (!pmAstromWriteBilevelChip(outHDU->header, outChip, WCS_NONLIN_TOL)) {
+		psWarning("Unable to generate WCS header.");
+		continue;
+	    }
+# endif
+        } else {
+            // we use a default FPA pixel scale of 1.0
+            if (!pmAstromReadWCS (output->fpa, outChip, inHDU->header, 1.0)) {
+                psWarning("Unable to read WCS astrometry for input FPA.");
+		continue;
+            }
+# if (UPDATE_HEADER)
+	    if (UPDATE_HEADER && !pmAstromWriteWCS(outHDU->header, output->fpa, outChip, WCS_NONLIN_TOL)) {
+		psWarning("Unable to generate WCS header.");
+		continue;
+	    }
+# endif
+	}
+    }
+
+    return true;
+}
+
+bool SetAstrometryFPA (pmConfig *config, pmFPAview *view, bool bilevelAstrometry) {
+
+    bool status = false;
+
+    if (!bilevelAstrometry) return true;
+
+    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
+    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
+
+    // loop over the available readouts
+    for (int i = 0; i < num; i++) {
+
+	// find the currently selected readout
+	pmFPAfile *output = pmFPAfileSelectSingle(config->files, "PSPHOT.STACK.OUTPUT.IMAGE", i); // File of interest
+	psAssert (output, "missing file?");
+
+# if (UPDATE_HEADER)
+	pmHDU *PHU = pmFPAviewThisPHU(view, output->fpa);
+	if (!PHU) {
+	    pmFPAAddSourceFromView(output->fpa, "name", view, output->format);
+	    PHU = pmFPAviewThisPHU (view, output->fpa);
+	    psAssert (PHU, "failed to make PHU");
+	}
+
+	if (!pmAstromWriteBilevelMosaic(PHU->header, output->fpa, WCS_NONLIN_TOL)) {
+	    psWarning("Unable to generate WCS header.");
+	}
+# endif
+    }
+
+    return true;
+}