Index: /trunk/psModules/src/objects/pmPSF_IO.c
===================================================================
--- /trunk/psModules/src/objects/pmPSF_IO.c	(revision 15227)
+++ /trunk/psModules/src/objects/pmPSF_IO.c	(revision 15228)
@@ -6,6 +6,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-09-28 00:38:42 $
+ *  @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-10-05 22:47:04 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -126,12 +126,17 @@
 bool pmPSFmodelWriteFPA (pmFPA *fpa, const pmFPAview *view, pmFPAfile *file, const pmConfig *config)
 {
+    pmFPAview *thisView = pmFPAviewAlloc (view->nRows);
+    *thisView = *view;
 
     for (int i = 0; i < fpa->chips->n; i++) {
         pmChip *chip = fpa->chips->data[i];
-        if (!pmPSFmodelWriteChip (chip, view, file, config)) {
+        thisView->chip = i;
+        if (!pmPSFmodelWriteChip (chip, thisView, file, config)) {
             psError(PS_ERR_IO, false, "Failed to write PSF for %dth chip", i);
-            return false;
-        }
-    }
+            psFree(thisView);
+            return false;
+        }
+    }
+    psFree(thisView);
     return true;
 }
@@ -141,6 +146,6 @@
 {
     if (!pmPSFmodelWrite (chip->analysis, view, file, config)) {
-	psError(PS_ERR_IO, false, "Failed to write PSF for chip");
-	return false;
+        psError(PS_ERR_IO, false, "Failed to write PSF for chip");
+        return false;
     }
     return true;
@@ -151,5 +156,5 @@
 //   - a PHU blank header
 // - image header        : FITS Image NAXIS = 0
-// if (trendMode == MAP) 
+// if (trendMode == MAP)
 //   - psf resid (+header) : FITS Image
 // else
@@ -186,40 +191,40 @@
     // define the EXTNAME values used for image header, table data, and residual image segments
     {
-	// lookup the EXTNAME values used for table data and image header segments
-	char *rule = NULL;
-
-	// Menu of EXTNAME rules
-	psMetadata *menu = psMetadataLookupMetadata(&status, file->camera, "EXTNAME.RULES");
-	if (!menu) {
-	    psError(PS_ERR_UNKNOWN, true, "missing EXTNAME.RULES in camera.config");
-	    return false;
-	}
-
-	// EXTNAME for image header
-	rule = psMetadataLookupStr(&status, menu, "PSF.HEAD");
-	if (!rule) {
-	    psError(PS_ERR_UNKNOWN, false, "missing entry for PSF.HEAD in EXTNAME.RULES in camera.config");
-	    return false;
-	}
-	headName = pmFPAfileNameFromRule (rule, file, view);
-
-	// EXTNAME for table data
-	rule = psMetadataLookupStr(&status, menu, "PSF.TABLE");
-	if (!rule) {
-	    psError(PS_ERR_UNKNOWN, false, "missing entry for PSF.TABLE in EXTNAME.RULES in camera.config");
-	    psFree (headName);
-	    return false;
-	}
-	tableName = pmFPAfileNameFromRule (rule, file, view);
-
-	// EXTNAME for resid data
-	rule = psMetadataLookupStr(&status, menu, "PSF.RESID");
-	if (!rule) {
-	    psError(PS_ERR_UNKNOWN, false, "missing entry for PSF.RESID in EXTNAME.RULES in camera.config");
-	    psFree (headName);
-	    psFree (tableName);
-	    return false;
-	}
-	residName = pmFPAfileNameFromRule (rule, file, view);
+        // lookup the EXTNAME values used for table data and image header segments
+        char *rule = NULL;
+
+        // Menu of EXTNAME rules
+        psMetadata *menu = psMetadataLookupMetadata(&status, file->camera, "EXTNAME.RULES");
+        if (!menu) {
+            psError(PS_ERR_UNKNOWN, true, "missing EXTNAME.RULES in camera.config");
+            return false;
+        }
+
+        // EXTNAME for image header
+        rule = psMetadataLookupStr(&status, menu, "PSF.HEAD");
+        if (!rule) {
+            psError(PS_ERR_UNKNOWN, false, "missing entry for PSF.HEAD in EXTNAME.RULES in camera.config");
+            return false;
+        }
+        headName = pmFPAfileNameFromRule (rule, file, view);
+
+        // EXTNAME for table data
+        rule = psMetadataLookupStr(&status, menu, "PSF.TABLE");
+        if (!rule) {
+            psError(PS_ERR_UNKNOWN, false, "missing entry for PSF.TABLE in EXTNAME.RULES in camera.config");
+            psFree (headName);
+            return false;
+        }
+        tableName = pmFPAfileNameFromRule (rule, file, view);
+
+        // EXTNAME for resid data
+        rule = psMetadataLookupStr(&status, menu, "PSF.RESID");
+        if (!rule) {
+            psError(PS_ERR_UNKNOWN, false, "missing entry for PSF.RESID in EXTNAME.RULES in camera.config");
+            psFree (headName);
+            psFree (tableName);
+            return false;
+        }
+        residName = pmFPAfileNameFromRule (rule, file, view);
     }
 
@@ -227,18 +232,18 @@
     // if this header block is new, write it to disk
     if (hdu->header != file->header) {
-	// add EXTNAME, EXTHEAD, EXTTYPE to header
-	psMetadataAddStr (hdu->header, PS_LIST_TAIL, "EXTTABLE", PS_META_REPLACE, "name of table extension", tableName);
-	psMetadataAddStr (hdu->header, PS_LIST_TAIL, "EXTRESID", PS_META_REPLACE, "name of resid extension", residName);
-	psMetadataAddStr (hdu->header, PS_LIST_TAIL, "EXTTYPE", PS_META_REPLACE, "extension type", "IMAGE");
-	if (!file->wrote_phu) {
-	    // this hdu->header acts as the PHU: set EXTEND to be true
-	    psMetadataAddBool (hdu->header, PS_LIST_TAIL, "EXTEND", PS_META_REPLACE, "this file has extensions", true);
-	    file->wrote_phu = true;
-	}
-
-	psFitsWriteBlank (file->fits, hdu->header, headName);
-	psTrace ("pmFPAfile", 5, "wrote ext head %s (type: %d)\n", file->filename, file->type);
-	file->header = hdu->header;
-	psFree (headName);
+        // add EXTNAME, EXTHEAD, EXTTYPE to header
+        psMetadataAddStr (hdu->header, PS_LIST_TAIL, "EXTTABLE", PS_META_REPLACE, "name of table extension", tableName);
+        psMetadataAddStr (hdu->header, PS_LIST_TAIL, "EXTRESID", PS_META_REPLACE, "name of resid extension", residName);
+        psMetadataAddStr (hdu->header, PS_LIST_TAIL, "EXTTYPE", PS_META_REPLACE, "extension type", "IMAGE");
+        if (!file->wrote_phu) {
+            // this hdu->header acts as the PHU: set EXTEND to be true
+            psMetadataAddBool (hdu->header, PS_LIST_TAIL, "EXTEND", PS_META_REPLACE, "this file has extensions", true);
+            file->wrote_phu = true;
+        }
+
+        psFitsWriteBlank (file->fits, hdu->header, headName);
+        psTrace ("pmFPAfile", 5, "wrote ext head %s (type: %d)\n", file->filename, file->type);
+        file->header = hdu->header;
+        psFree (headName);
     }
 
@@ -246,32 +251,32 @@
     pmPSF *psf = psMetadataLookupPtr (&status, analysis, "PSPHOT.PSF");
     if (!psf) {
-	psError(PS_ERR_UNKNOWN, true, "missing PSF for this analysis metadata");
-	psFree (tableName);
-	psFree (residName);
-	return false;
+        psError(PS_ERR_UNKNOWN, true, "missing PSF for this analysis metadata");
+        psFree (tableName);
+        psFree (residName);
+        return false;
     }
 
     // write the PSF model parameters in a FITS table
-    { 
-	// we need to write a header for the table,
-	psMetadata *header = psMetadataAlloc();
-
-	char *modelName = pmModelClassGetName (psf->type);
-	psMetadataAddStr (header, PS_LIST_TAIL, "PSF_NAME", 0, "PSF model name", modelName);
-
-	psMetadataAddBool (header, PS_LIST_TAIL, "ERR_LMM",  0, "Use Poisson errors in fits?", psf->poissonErrorsPhotLMM);
-	psMetadataAddBool (header, PS_LIST_TAIL, "ERR_LIN",  0, "Use Poisson errors in fits?", psf->poissonErrorsPhotLin);
-	psMetadataAddBool (header, PS_LIST_TAIL, "ERR_PAR",  0, "Use Poisson errors in fits?", psf->poissonErrorsParams);
-
-	int nPar = pmModelClassParameterCount (psf->type)    ;
-	psMetadataAdd (header, PS_LIST_TAIL, "PSF_NPAR", PS_DATA_S32, "PSF model parameter count", nPar);
-
-	psMetadataAddS32 (header, PS_LIST_TAIL, "IMAXIS1", 0, "Image X Size", psf->fieldNx);
-	psMetadataAddS32 (header, PS_LIST_TAIL, "IMAXIS2", 0, "Image Y Size", psf->fieldNy);
-	psMetadataAddS32 (header, PS_LIST_TAIL, "IMREF1",  0, "Image X Ref",  psf->fieldXo);
-	psMetadataAddS32 (header, PS_LIST_TAIL, "IMREF2",  0, "Image Y Ref",  psf->fieldYo);
-
-	// extract PSF Clump info
-	pmPSFClump psfClump;
+    {
+        // we need to write a header for the table,
+        psMetadata *header = psMetadataAlloc();
+
+        char *modelName = pmModelClassGetName (psf->type);
+        psMetadataAddStr (header, PS_LIST_TAIL, "PSF_NAME", 0, "PSF model name", modelName);
+
+        psMetadataAddBool (header, PS_LIST_TAIL, "ERR_LMM",  0, "Use Poisson errors in fits?", psf->poissonErrorsPhotLMM);
+        psMetadataAddBool (header, PS_LIST_TAIL, "ERR_LIN",  0, "Use Poisson errors in fits?", psf->poissonErrorsPhotLin);
+        psMetadataAddBool (header, PS_LIST_TAIL, "ERR_PAR",  0, "Use Poisson errors in fits?", psf->poissonErrorsParams);
+
+        int nPar = pmModelClassParameterCount (psf->type)    ;
+        psMetadataAdd (header, PS_LIST_TAIL, "PSF_NPAR", PS_DATA_S32, "PSF model parameter count", nPar);
+
+        psMetadataAddS32 (header, PS_LIST_TAIL, "IMAXIS1", 0, "Image X Size", psf->fieldNx);
+        psMetadataAddS32 (header, PS_LIST_TAIL, "IMAXIS2", 0, "Image Y Size", psf->fieldNy);
+        psMetadataAddS32 (header, PS_LIST_TAIL, "IMREF1",  0, "Image X Ref",  psf->fieldXo);
+        psMetadataAddS32 (header, PS_LIST_TAIL, "IMREF2",  0, "Image Y Ref",  psf->fieldYo);
+
+        // extract PSF Clump info
+        pmPSFClump psfClump;
 
         psfClump.X  = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.X");   assert (status);
@@ -280,139 +285,139 @@
         psfClump.dY = psMetadataLookupF32 (&status, recipe, "PSF.CLUMP.DY");  assert (status);
 
-	psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_CLX", 0, "psf clump center", psfClump.X);
-	psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_CLY", 0, "psf clump center", psfClump.Y);
-	psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_CLDX", 0, "psf clump size", psfClump.dX);
-	psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_CLDY", 0, "psf clump size", psfClump.dY);
-
-	// save the dimensions of each parameter
-	for (int i = 0; i < nPar; i++) {
-	    char name[9];
-	    int nX, nY;
-
-	    pmTrend2D *trend = psf->params->data[i];
-	    if (trend == NULL) continue;
-
-	    if (trend->mode == PM_TREND_MAP) {
-	      nX = trend->map->map->numCols;
-	      nY = trend->map->map->numRows;
-	    } else {
-	      nX = trend->poly->nX;
-	      nY = trend->poly->nY;
-	    }
-	    snprintf (name, 9, "PAR%02d_NX", i);
-	    psMetadataAddS32 (header, PS_LIST_TAIL, name, 0, "", nX);
-	    snprintf (name, 9, "PAR%02d_NY", i);
-	    psMetadataAddS32 (header, PS_LIST_TAIL, name, 0, "", nY);
-	    snprintf (name, 9, "PAR%02d_MD", i);
-	    char *modeName = pmTrend2DModeToString (trend->mode);
-	    psMetadataAddStr (header, PS_LIST_TAIL, name, 0, "", modeName);
-	    psFree (modeName);
-	}
-
-	// other required information describing the PSF
-	psMetadataAddF32 (header, PS_LIST_TAIL, "AP_RESID", 0, "aperture residual", psf->ApResid);
-	psMetadataAddF32 (header, PS_LIST_TAIL, "AP_ERROR", 0, "aperture residual scatter", psf->dApResid);
-	psMetadataAddF32 (header, PS_LIST_TAIL, "CHISQ",    0, "chi-square for fit", psf->chisq);
-	psMetadataAddS32 (header, PS_LIST_TAIL, "NSTARS",   0, "number of stars used to measure PSF", psf->nPSFstars);
-
-	// XXX can we drop this now?
-	psMetadataAddF32 (header, PS_LIST_TAIL, "SKY_BIAS", PS_DATA_F32, "sky bias level", psf->skyBias);
-
-	// build a FITS table of the PSF parameters
-	psArray *psfTable = psArrayAllocEmpty (100);
-	for (int i = 0; i < nPar; i++) {
-	    pmTrend2D *trend = psf->params->data[i];
-	    if (trend == NULL) continue; // skip unset parameters (eg, XPOS)
-
-	    if (trend->mode == PM_TREND_MAP) {
-		// write the image components into a table: this is needed because they may each be a different size
-		psImageMap *map = trend->map;
-		for (int ix = 0; ix < map->map->numCols; ix++) {
-		    for (int iy = 0; iy < map->map->numRows; iy++) {
-			psMetadata *row = psMetadataAlloc ();
-			psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TERM", 0, "", i);
-			psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER",    0, "", ix);
-			psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER",    0, "", iy);
-			psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE",      0, "", map->map->data.F32[iy][ix]);
-			psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR",      0, "", map->error->data.F32[iy][ix]);
-			psMetadataAddU8  (row, PS_LIST_TAIL, "MASK",       0, "", 0); // no cells are masked
-
-			psArrayAdd (psfTable, 100, row);
-			psFree (row);
-		    }
-		}
-	    } else {
-		// write the polynomial components into a table
-		psPolynomial2D *poly = trend->poly;
-		for (int ix = 0; ix <= poly->nX; ix++) {
-		    for (int iy = 0; iy <= poly->nY; iy++) {
-			psMetadata *row = psMetadataAlloc ();
-			psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TERM", 0, "", i);
-			psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER",    0, "", ix);
-			psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER",    0, "", iy);
-			psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE",      0, "", poly->coeff[ix][iy]);
-			psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR",      0, "", poly->coeffErr[ix][iy]);
-			psMetadataAddU8  (row, PS_LIST_TAIL, "MASK",       0, "", poly->mask[ix][iy]);
-
-			psArrayAdd (psfTable, 100, row);
-			psFree (row);
-		    }
-		}
-	    }
-	}
-
-	// write an empty FITS segment if we have no PSF information
-	if (psfTable->n == 0) {
-	    // XXX this is probably an error (if we have a PSF, how do we have no data?)
-	    psFitsWriteBlank (file->fits, header, tableName);
-	} else {
-	    psTrace ("pmFPAfile", 5, "writing psf data %s\n", tableName);
-	    if (!psFitsWriteTable (file->fits, header, psfTable, tableName)) {
-		psError(PS_ERR_IO, false, "writing psf table data %s\n", tableName);
-		psFree (tableName);
-		psFree (residName);
-		psFree (psfTable);
-		psFree (header);
-		return false;
-	    }
-	}
-	psFree (tableName);
-	psFree (psfTable);
-	psFree (header);
+        psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_CLX", 0, "psf clump center", psfClump.X);
+        psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_CLY", 0, "psf clump center", psfClump.Y);
+        psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_CLDX", 0, "psf clump size", psfClump.dX);
+        psMetadataAddF32 (header, PS_LIST_TAIL, "PSF_CLDY", 0, "psf clump size", psfClump.dY);
+
+        // save the dimensions of each parameter
+        for (int i = 0; i < nPar; i++) {
+            char name[9];
+            int nX, nY;
+
+            pmTrend2D *trend = psf->params->data[i];
+            if (trend == NULL) continue;
+
+            if (trend->mode == PM_TREND_MAP) {
+              nX = trend->map->map->numCols;
+              nY = trend->map->map->numRows;
+            } else {
+              nX = trend->poly->nX;
+              nY = trend->poly->nY;
+            }
+            snprintf (name, 9, "PAR%02d_NX", i);
+            psMetadataAddS32 (header, PS_LIST_TAIL, name, 0, "", nX);
+            snprintf (name, 9, "PAR%02d_NY", i);
+            psMetadataAddS32 (header, PS_LIST_TAIL, name, 0, "", nY);
+            snprintf (name, 9, "PAR%02d_MD", i);
+            char *modeName = pmTrend2DModeToString (trend->mode);
+            psMetadataAddStr (header, PS_LIST_TAIL, name, 0, "", modeName);
+            psFree (modeName);
+        }
+
+        // other required information describing the PSF
+        psMetadataAddF32 (header, PS_LIST_TAIL, "AP_RESID", 0, "aperture residual", psf->ApResid);
+        psMetadataAddF32 (header, PS_LIST_TAIL, "AP_ERROR", 0, "aperture residual scatter", psf->dApResid);
+        psMetadataAddF32 (header, PS_LIST_TAIL, "CHISQ",    0, "chi-square for fit", psf->chisq);
+        psMetadataAddS32 (header, PS_LIST_TAIL, "NSTARS",   0, "number of stars used to measure PSF", psf->nPSFstars);
+
+        // XXX can we drop this now?
+        psMetadataAddF32 (header, PS_LIST_TAIL, "SKY_BIAS", PS_DATA_F32, "sky bias level", psf->skyBias);
+
+        // build a FITS table of the PSF parameters
+        psArray *psfTable = psArrayAllocEmpty (100);
+        for (int i = 0; i < nPar; i++) {
+            pmTrend2D *trend = psf->params->data[i];
+            if (trend == NULL) continue; // skip unset parameters (eg, XPOS)
+
+            if (trend->mode == PM_TREND_MAP) {
+                // write the image components into a table: this is needed because they may each be a different size
+                psImageMap *map = trend->map;
+                for (int ix = 0; ix < map->map->numCols; ix++) {
+                    for (int iy = 0; iy < map->map->numRows; iy++) {
+                        psMetadata *row = psMetadataAlloc ();
+                        psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TERM", 0, "", i);
+                        psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER",    0, "", ix);
+                        psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER",    0, "", iy);
+                        psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE",      0, "", map->map->data.F32[iy][ix]);
+                        psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR",      0, "", map->error->data.F32[iy][ix]);
+                        psMetadataAddU8  (row, PS_LIST_TAIL, "MASK",       0, "", 0); // no cells are masked
+
+                        psArrayAdd (psfTable, 100, row);
+                        psFree (row);
+                    }
+                }
+            } else {
+                // write the polynomial components into a table
+                psPolynomial2D *poly = trend->poly;
+                for (int ix = 0; ix <= poly->nX; ix++) {
+                    for (int iy = 0; iy <= poly->nY; iy++) {
+                        psMetadata *row = psMetadataAlloc ();
+                        psMetadataAddS32 (row, PS_LIST_TAIL, "MODEL_TERM", 0, "", i);
+                        psMetadataAddS32 (row, PS_LIST_TAIL, "X_POWER",    0, "", ix);
+                        psMetadataAddS32 (row, PS_LIST_TAIL, "Y_POWER",    0, "", iy);
+                        psMetadataAddF32 (row, PS_LIST_TAIL, "VALUE",      0, "", poly->coeff[ix][iy]);
+                        psMetadataAddF32 (row, PS_LIST_TAIL, "ERROR",      0, "", poly->coeffErr[ix][iy]);
+                        psMetadataAddU8  (row, PS_LIST_TAIL, "MASK",       0, "", poly->mask[ix][iy]);
+
+                        psArrayAdd (psfTable, 100, row);
+                        psFree (row);
+                    }
+                }
+            }
+        }
+
+        // write an empty FITS segment if we have no PSF information
+        if (psfTable->n == 0) {
+            // XXX this is probably an error (if we have a PSF, how do we have no data?)
+            psFitsWriteBlank (file->fits, header, tableName);
+        } else {
+            psTrace ("pmFPAfile", 5, "writing psf data %s\n", tableName);
+            if (!psFitsWriteTable (file->fits, header, psfTable, tableName)) {
+                psError(PS_ERR_IO, false, "writing psf table data %s\n", tableName);
+                psFree (tableName);
+                psFree (residName);
+                psFree (psfTable);
+                psFree (header);
+                return false;
+            }
+        }
+        psFree (tableName);
+        psFree (psfTable);
+        psFree (header);
     }
 
     // write the residual images (3D)
     {
-	psMetadata *header = psMetadataAlloc ();
-	if (psf->residuals == NULL) {
-	    // set some header keywords to make it clear there are no residuals?
-	    psFitsWriteBlank (file->fits, header, residName);
-	    psFree (residName);
-	    psFree (header);
-	    return true;
-	}
-
-	psMetadataAddS32 (header, PS_LIST_TAIL, "XBIN",    0, "", psf->residuals->xBin);
-	psMetadataAddS32 (header, PS_LIST_TAIL, "YBIN",    0, "", psf->residuals->yBin);
-	psMetadataAddS32 (header, PS_LIST_TAIL, "XCENTER", 0, "", psf->residuals->xCenter);
-	psMetadataAddS32 (header, PS_LIST_TAIL, "YCENTER", 0, "", psf->residuals->yCenter);
-
-	// write the residuals as three planes of the image
-	// this call creates an extension with NAXIS3 = 3
-	if (psf->residuals->Rx) {
-	    // this call creates an extension with NAXIS3 = 3
-	    psArray *images = psArrayAllocEmpty (3);
-	    psArrayAdd (images, 1, psf->residuals->Ro);
-	    psArrayAdd (images, 1, psf->residuals->Rx);
-	    psArrayAdd (images, 1, psf->residuals->Ry);
-
-	    psFitsWriteImageCube (file->fits, header, images, residName);
-	    psFree (images);
-	} else {
-	    // this call creates an extension with NAXIS3 = 1
-	    psFitsWriteImage(file->fits, header, psf->residuals->Ro, 0, residName);
-	}
-	psFree (residName);
-	psFree (header);
+        psMetadata *header = psMetadataAlloc ();
+        if (psf->residuals == NULL) {
+            // set some header keywords to make it clear there are no residuals?
+            psFitsWriteBlank (file->fits, header, residName);
+            psFree (residName);
+            psFree (header);
+            return true;
+        }
+
+        psMetadataAddS32 (header, PS_LIST_TAIL, "XBIN",    0, "", psf->residuals->xBin);
+        psMetadataAddS32 (header, PS_LIST_TAIL, "YBIN",    0, "", psf->residuals->yBin);
+        psMetadataAddS32 (header, PS_LIST_TAIL, "XCENTER", 0, "", psf->residuals->xCenter);
+        psMetadataAddS32 (header, PS_LIST_TAIL, "YCENTER", 0, "", psf->residuals->yCenter);
+
+        // write the residuals as three planes of the image
+        // this call creates an extension with NAXIS3 = 3
+        if (psf->residuals->Rx) {
+            // this call creates an extension with NAXIS3 = 3
+            psArray *images = psArrayAllocEmpty (3);
+            psArrayAdd (images, 1, psf->residuals->Ro);
+            psArrayAdd (images, 1, psf->residuals->Rx);
+            psArrayAdd (images, 1, psf->residuals->Ry);
+
+            psFitsWriteImageCube (file->fits, header, images, residName);
+            psFree (images);
+        } else {
+            // this call creates an extension with NAXIS3 = 1
+            psFitsWriteImage(file->fits, header, psf->residuals->Ro, 0, residName);
+        }
+        psFree (residName);
+        psFree (header);
     }
 
@@ -465,14 +470,14 @@
     psMetadata *outhead = psMetadataAlloc();
     if (phu) {
-	psMetadataCopy (outhead, phu->header);
+        psMetadataCopy (outhead, phu->header);
     } else {
-	pmConfigConformHeader (outhead, file->format);
-
-	psMetadata *fileData = psMetadataLookupMetadata(NULL, file->format, "FILE"); // File information
-	const char *fpaNameHdr = psMetadataLookupStr(NULL, fileData, "FPA.NAME");
-	if (fpaNameHdr && strlen(fpaNameHdr) > 0) {
-	    const char *fpaName = psMetadataLookupStr(NULL, file->fpa->concepts, "FPA.NAME");
-	    psMetadataAddStr(outhead, PS_LIST_TAIL, fpaNameHdr, PS_META_REPLACE, "FPA name", fpaName);
-	}
+        pmConfigConformHeader (outhead, file->format);
+
+        psMetadata *fileData = psMetadataLookupMetadata(NULL, file->format, "FILE"); // File information
+        const char *fpaNameHdr = psMetadataLookupStr(NULL, fileData, "FPA.NAME");
+        if (fpaNameHdr && strlen(fpaNameHdr) > 0) {
+            const char *fpaName = psMetadataLookupStr(NULL, file->fpa->concepts, "FPA.NAME");
+            psMetadataAddStr(outhead, PS_LIST_TAIL, fpaNameHdr, PS_META_REPLACE, "FPA name", fpaName);
+        }
     }
 
@@ -524,6 +529,6 @@
 {
     if (!pmPSFmodelRead (chip->analysis, view, file, config)) {
-	psError(PS_ERR_IO, false, "Failed to write PSF for chip");
-	return false;
+        psError(PS_ERR_IO, false, "Failed to write PSF for chip");
+        return false;
     }
     return true;
@@ -628,27 +633,27 @@
         snprintf (name, 9, "PAR%02d_NX", i);
         binning->nXruff = psMetadataLookupS32 (&status, header, name);
-	if (!status) continue;		// not all parameters are defined
+        if (!status) continue;          // not all parameters are defined
 
         snprintf (name, 9, "PAR%02d_NY", i);
         binning->nYruff = psMetadataLookupS32 (&status, header, name);
-	if (!status) {
-	    psError(PS_ERR_UNKNOWN, true, "inconsistent PSF header: NX defined for PAR %d, but not NY", i);
-	    return false;
-	}
+        if (!status) {
+            psError(PS_ERR_UNKNOWN, true, "inconsistent PSF header: NX defined for PAR %d, but not NY", i);
+            return false;
+        }
 
         snprintf (name, 9, "PAR%02d_MD", i);
         char *modeName = psMetadataLookupStr (&status, header, name);
-	if (!status) {
-	    psError(PS_ERR_UNKNOWN, true, "inconsistent PSF header: NX & NY defined for PAR %d, but not MD", i);
-	    return false;
-	}
-	pmTrend2DMode psfTrendMode = pmTrend2DModeFromString (modeName);
-	if (psfTrendMode == PM_TREND_NONE) {
-	    psfTrendMode = PM_TREND_POLY_ORD;
-	}
-
-	psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER);
-	psImageBinningSetSkipByOffset (binning, options->psfFieldXo, options->psfFieldYo);
-	psf->params->data[i] = pmTrend2DNoImageAlloc (psfTrendMode, binning, NULL);
+        if (!status) {
+            psError(PS_ERR_UNKNOWN, true, "inconsistent PSF header: NX & NY defined for PAR %d, but not MD", i);
+            return false;
+        }
+        pmTrend2DMode psfTrendMode = pmTrend2DModeFromString (modeName);
+        if (psfTrendMode == PM_TREND_NONE) {
+            psfTrendMode = PM_TREND_POLY_ORD;
+        }
+
+        psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER);
+        psImageBinningSetSkipByOffset (binning, options->psfFieldXo, options->psfFieldYo);
+        psf->params->data[i] = pmTrend2DNoImageAlloc (psfTrendMode, binning, NULL);
     }
     psFree (binning);
@@ -674,31 +679,31 @@
 
         pmTrend2D *trend = psf->params->data[iPar];
-	if (trend == NULL) {
-	    psError(PS_ERR_UNKNOWN, true, "parameter %d not available", iPar);
-	    return false;
-	}
-
-	if (trend->mode == PM_TREND_MAP) {
-	    psImageMap *map = trend->map;
-	    assert (map);
-	    assert (map->map);
-	    assert (map->error);
-	    assert (xPow >= 0);
-	    assert (yPow >= 0);
-	    assert (xPow < map->map->numCols);
-	    assert (yPow < map->map->numRows);
-	    map->map->data.F32[yPow][xPow]    = psMetadataLookupF32 (&status, row, "VALUE");
-	    map->error->data.F32[yPow][xPow]  = psMetadataLookupF32 (&status, row, "ERROR");
-	} else {
-	    psPolynomial2D *poly = trend->poly;
-	    assert (poly);
-	    assert (xPow >= 0);
-	    assert (yPow >= 0);
-	    assert (xPow <= poly->nX);
-	    assert (yPow <= poly->nY);
-	    poly->coeff[xPow][yPow]    = psMetadataLookupF32 (&status, row, "VALUE");
-	    poly->coeffErr[xPow][yPow] = psMetadataLookupF32 (&status, row, "ERROR");
-	    poly->mask[xPow][yPow]     = psMetadataLookupU8  (&status, row, "MASK");
-	}
+        if (trend == NULL) {
+            psError(PS_ERR_UNKNOWN, true, "parameter %d not available", iPar);
+            return false;
+        }
+
+        if (trend->mode == PM_TREND_MAP) {
+            psImageMap *map = trend->map;
+            assert (map);
+            assert (map->map);
+            assert (map->error);
+            assert (xPow >= 0);
+            assert (yPow >= 0);
+            assert (xPow < map->map->numCols);
+            assert (yPow < map->map->numRows);
+            map->map->data.F32[yPow][xPow]    = psMetadataLookupF32 (&status, row, "VALUE");
+            map->error->data.F32[yPow][xPow]  = psMetadataLookupF32 (&status, row, "ERROR");
+        } else {
+            psPolynomial2D *poly = trend->poly;
+            assert (poly);
+            assert (xPow >= 0);
+            assert (yPow >= 0);
+            assert (xPow <= poly->nX);
+            assert (yPow <= poly->nY);
+            poly->coeff[xPow][yPow]    = psMetadataLookupF32 (&status, row, "VALUE");
+            poly->coeffErr[xPow][yPow] = psMetadataLookupF32 (&status, row, "ERROR");
+            poly->mask[xPow][yPow]     = psMetadataLookupU8  (&status, row, "MASK");
+        }
     }
     psFree (header);