Index: /trunk/psphot/Makefile
===================================================================
--- /trunk/psphot/Makefile	(revision 4250)
+++ /trunk/psphot/Makefile	(revision 4251)
@@ -23,4 +23,5 @@
 $(SRC)/psPolynomials.$(ARCH).o \
 $(SRC)/psUtils.$(ARCH).o \
+$(SRC)/load_args.$(ARCH).o \
 $(SRC)/setup.$(ARCH).o \
 $(SRC)/image_stats.$(ARCH).o \
@@ -37,4 +38,5 @@
 $(SRC)/fit_galaxies.$(ARCH).o \
 $(SRC)/subtract_galaxies.$(ARCH).o \
+$(SRC)/LocalSky.$(ARCH).o \
 $(SRC)/find_defects.$(ARCH).o
 
Index: /trunk/psphot/src/LocalSky.c
===================================================================
--- /trunk/psphot/src/LocalSky.c	(revision 4250)
+++ /trunk/psphot/src/LocalSky.c	(revision 4251)
@@ -7,11 +7,17 @@
 			  psF32 Radius)
 {
+    psRegion *srcRegion;
+
     // Grab a subimage of the original image of size (2 * outerRadius).
     srcRegion = psRegionSquare (x, y, Radius);
-    srcRegion = psRegionForImage (srcRegion, imdata->image);
+    srcRegion = psRegionForImage (srcRegion, imdata->image, srcRegion);
 
-    psImage *subImage = psImageSubset(imdata->image, srcRegion);
-    psImage *subImageMask = psImageSubset(imdata->mask, srcRegion);
-    psImage *subImageNoise = psImageSubset(imdata->noise, srcRegion);
+//  psImage *subImage = psImageSubset(imdata->image, srcRegion);
+//  psImage *subImageMask = psImageSubset(imdata->mask, srcRegion);
+//  psImage *subImageNoise = psImageSubset(imdata->noise, srcRegion);
+
+    psImage *subImage = psImageSubset(imdata->image, srcRegion->x0, srcRegion->y0, srcRegion->x1, srcRegion->y1);
+    psImage *subImageMask = psImageSubset(imdata->mask, srcRegion->x0, srcRegion->y0, srcRegion->x1, srcRegion->y1);
+    psImage *subImageNoise = psImageSubset(imdata->noise, srcRegion->x0, srcRegion->y0, srcRegion->x1, srcRegion->y1);
 
     mySource->pixels = subImage;
@@ -30,25 +36,21 @@
     psImage *image = source->pixels;
     psImage *mask  = source->mask;
-    psImage *peak  = source->peak;
+    psPeak *peak  = source->peak;
+    psRegion *srcRegion;
 
     srcRegion = psRegionSquare (peak->x, peak->y, Radius);
-    srcRegion = psRegionForImage (srcRegion, imdata->mask);
-    psImageMaskRegion (imdata->mask, region, OR, 0x80);
+    srcRegion = psRegionForImage (srcRegion, mask, srcRegion);
+    psImageMaskRegion (mask, srcRegion, OR, 0x80);
 
     psStats *myStats = psStatsAlloc(statsOptions);
     myStats = psImageStats(myStats, image, mask, 0xff);
-    psImageMaskRegion (imdata->mask, region, AND, 0x7f);
+    psImageMaskRegion (mask, srcRegion, AND, 0x7f);
 
-    mySource->moments = pmMomentsAlloc();
+    source->moments = pmMomentsAlloc();
     psF64 tmpF64;
     p_psGetStatValue(myStats, &tmpF64);
-    mySource->moments->Sky = (psF32) tmpF64;
+    source->moments->Sky = (psF32) tmpF64;
 
     return (true);
 }
 
-// define a square region centered on the given coordinate
-psRegion *psRegionSquare (psF32 x, psF32 y, psF32 radius) {
-    region = psRegionAlloc (x - Radius, y - Radius, x + Radius + 1, y + Radius + 1);
-    return (region);
-}
Index: /trunk/psphot/src/apply_psf_model.c
===================================================================
--- /trunk/psphot/src/apply_psf_model.c	(revision 4250)
+++ /trunk/psphot/src/apply_psf_model.c	(revision 4251)
@@ -5,7 +5,9 @@
 // run this function to a specific flux limit?
 
-bool apply_psf_model (psImage *image, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) 
+bool apply_psf_model (psImageData *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) 
 { 
     bool  status;
+    float x;
+    float y;
     int   Nfit = 0;
     int   Niter = 0;
@@ -42,4 +44,5 @@
 	x = model->params->data.F32[2];
 	y = model->params->data.F32[3];
+	// XXX I need to check if the model center has moved too much relative to the peak
 
 	// set the fit radius based on the object flux limit and the model
@@ -51,7 +54,4 @@
 	}
 	
-	// mask off saturated pixels (move this to setup?)
-	pmSourceMaskSaturated (source, SATURATE);
-
 	// check if we need to redefine the pixels
 	if (model->radius > OUTER) {
Index: /trunk/psphot/src/choose_psf_model.c
===================================================================
--- /trunk/psphot/src/choose_psf_model.c	(revision 4250)
+++ /trunk/psphot/src/choose_psf_model.c	(revision 4251)
@@ -3,5 +3,5 @@
 // test PSF models and select best option
 
-pmPSF *choose_psf_model (psImage *image, psMetadata *config, psArray *sources) 
+pmPSF *choose_psf_model (psMetadata *config, psArray *sources) 
 { 
     bool        status;
Index: /trunk/psphot/src/find_defects.c
===================================================================
--- /trunk/psphot/src/find_defects.c	(revision 4250)
+++ /trunk/psphot/src/find_defects.c	(revision 4251)
@@ -1,6 +1,4 @@
 # include "psphot.h"
 # define NSIGMA 3.0
-
-void dump_image (psImage *image, char *filename);
 
 bool find_defects (psImage *zapmask, psArray *sources, psMetadata *config, pmPSF *psf) 
@@ -34,5 +32,5 @@
   // generate the image flux
   pmSourceAddModel (flux, NULL, model, true);
-  dump_image (flux, "psf.fits");
+  DumpImage (flux, "psf.fits");
 
   fprintf (stderr, "sx: %f, sy: %f\n", model->params->data.F32[4], model->params->data.F32[5]);
Index: /trunk/psphot/src/find_defects_new.c
===================================================================
--- /trunk/psphot/src/find_defects_new.c	(revision 4251)
+++ /trunk/psphot/src/find_defects_new.c	(revision 4251)
@@ -0,0 +1,12 @@
+# include "psphot.h"
+
+bool find_defects (psImage *zapmask, psArray *sources, psMetadata *config, pmPSF *psf) 
+{ 
+
+  psTimerStart ("psphot");
+
+  // for each source, search for all row peaks
+  // measure 
+
+  return (true);
+}
Index: /trunk/psphot/src/fit_galaxies.c
===================================================================
--- /trunk/psphot/src/fit_galaxies.c	(revision 4250)
+++ /trunk/psphot/src/fit_galaxies.c	(revision 4251)
@@ -8,9 +8,4 @@
     int   Nfit = 0;
     int   Niter = 0;
-
-    float XBORDER  = psMetadataLookupF32 (&status, config, "XBORDER");
-    float YBORDER  = psMetadataLookupF32 (&status, config, "YBORDER");
-    psRegion *keep = psRegionAlloc (XBORDER, -XBORDER, YBORDER, -YBORDER);
-    keep           = psRegionForImage (keep, image, keep);
 
     float RADIUS   = psMetadataLookupF32 (&status, config, "FIT_RADIUS");
@@ -25,6 +20,11 @@
 
 	// need a better model guess and a better radius choice
-	pmSourceSetPixelCircle (source, image, RADIUS);
-	pmSourceMaskRegion (source, keep);
+	// when radius is not fixed, we will need to check if new radius fits on image
+# if (0)
+	if (model->radius > OUTER) {
+	  // allocate image, noise, mask arrays for each peak (square of radius OUTER)
+	  pmSourceDefinePixels (source, imdata, x, y, OUTER);
+	}
+# endif
 	psModel  *model  = pmSourceModelGuess (source, modelType); 
 
Index: /trunk/psphot/src/load_args.c
===================================================================
--- /trunk/psphot/src/load_args.c	(revision 4250)
+++ /trunk/psphot/src/load_args.c	(revision 4251)
@@ -2,4 +2,11 @@
 
 psMetadata *load_args (int *argc, char **argv) {
+
+  int N, Nfail;
+
+  // basic pslib options
+  psLogSetFormat ("M");
+  psLogArguments (argc, argv);
+  psTraceArguments (argc, argv);
 
   // identify options in args list - these go on config 
@@ -19,5 +26,5 @@
   }
 
-  if (argc != 4) usage ();
+  if (*argc != 4) usage ();
 
   psMetadata *config = psMetadataParseConfig (NULL, &Nfail, argv[3], FALSE);
Index: /trunk/psphot/src/psUtils.c
===================================================================
--- /trunk/psphot/src/psUtils.c	(revision 4250)
+++ /trunk/psphot/src/psUtils.c	(revision 4251)
@@ -1,21 +1,3 @@
 # include "psphot.h"
-
-// set actual region based on image parameters
-// XXX EAM : this needs to be changes to use psRegion rather than psRegion*
-psRegion *psRegionForImage (psRegion *out, psImage *image, psRegion *in) {
-    
-    if (out == NULL) {
-	out = psRegionAlloc(in->x0, in->x1, in->y0, in->y1);
-    } else {
-	*out = *in;
-    }
-    if (out->x1 <= 0) {
-	out->x1 = image->numCols + out->x1;
-    }
-    if (out->y1 <= 0) {
-	out->y1 = image->numRows + out->y1;
-    }
-    return (out);
-}
 
 static psHash *timers = NULL;
@@ -336,2 +318,132 @@
   return (false);
 }	    
+
+// define a square region centered on the given coordinate
+psRegion *psRegionSquare (psF32 x, psF32 y, psF32 radius) {
+    psRegion *region;
+    region = psRegionAlloc (x - radius, x + radius + 1,    
+			    y - radius, y + radius + 1);
+    return (region);
+}
+
+// set actual region based on image parameters:
+// compensate for negative upper limits
+// XXX this is inconsistent: the coordindates should always be in the parent
+//     frame, which means the negative values should subtract from Nx,Ny of
+//     the parent, not the child.  but, we don't carry the dimensions of the 
+//     parent in the psImage container.  for now, us the child Nx,Ny
+// force range to be on this subimage 
+// XXX EAM : this needs to be changes to use psRegion rather than psRegion*
+psRegion *psRegionForImage (psRegion *out, psImage *image, psRegion *in) {
+    
+    // x0,y0, x1,y1 are in *parent* units
+
+    if (out == NULL) {
+	out = psRegionAlloc(in->x0, in->x1, in->y0, in->y1);
+    } else {
+	*out = *in;
+    }
+    // XXX these are probably wrong (see above)
+    if (out->x1 <= 0) {
+	out->x1 = image->col0 + image->numCols + out->x1;
+    }
+    if (out->y1 <= 0) {
+	out->y1 = image->row0 + image->numRows + out->y1;
+    }
+
+    // force the lower-limits to be on the child
+    out->x0 = PS_MAX(image->col0, out->x0);
+    out->y0 = PS_MAX(image->row0, out->y0);
+
+    // force the upper-limits to be on the child
+    out->x1 = PS_MIN(image->col0 + image->numCols, out->x1);
+    out->y1 = PS_MIN(image->row0 + image->numRows, out->y1);
+    return (out);
+}
+
+// mask the area contained by the region
+// the region is defined wrt the parent image
+void psImageMaskRegion (psImage *image, psRegion *region, bool logical_and, int maskValue) {
+
+    for (int iy = 0; iy < image->numRows; iy++) {
+	for (int ix = 0; ix < image->numCols; ix++) {
+	    if (ix + image->col0 <  region->x0) continue;
+	    if (ix + image->col0 >= region->x1) continue;
+	    if (iy + image->row0 <  region->y0) continue;
+	    if (iy + image->row0 >= region->y1) continue;
+	    if (logical_and) {
+		image->data.U8[iy][ix] &= maskValue;
+	    } else {
+		image->data.U8[iy][ix] |= maskValue;
+	    }
+	}
+    }
+}
+
+// mask the area not contained by the region
+// the region is defined wrt the parent image
+void psImageKeepRegion (psImage *image, psRegion *region, bool logical_and, int maskValue) {
+
+    for (int iy = 0; iy < image->numRows; iy++) {
+	for (int ix = 0; ix < image->numCols; ix++) {
+	    if (ix + image->col0 <  region->x0) goto maskit;
+	    if (ix + image->col0 >= region->x1) goto maskit;
+	    if (iy + image->row0 <  region->y0) goto maskit;
+	    if (iy + image->row0 >= region->y1) goto maskit;
+	    continue;
+	maskit:
+	    if (logical_and) {
+		image->data.U8[iy][ix] &= maskValue;
+	    } else {
+		image->data.U8[iy][ix] |= maskValue;
+	    }
+	}
+    }
+}
+
+// mask the area contained by the region
+// the region is defined wrt the parent image
+void psImageMaskCircle (psImage *image, double x, double y, double radius, bool logical_and, int maskValue) {
+
+    double dx, dy, r2, R2;
+
+    R2 = PS_SQR(radius);
+
+    for (int iy = 0; iy < image->numRows; iy++) {
+	for (int ix = 0; ix < image->numCols; ix++) {
+	    dx = ix + image->col0 - x;
+	    dy = iy + image->row0 - y;
+	    r2 = PS_SQR(dx) + PS_SQR(dy);
+	    if (r2 > R2) continue;
+	    if (logical_and) {
+		image->data.U8[iy][ix] &= maskValue;
+	    } else {
+		image->data.U8[iy][ix] |= maskValue;
+	    }
+	}
+    }
+}
+
+// mask the area contained by the region
+// the region is defined wrt the parent image
+void psImageKeepCircle (psImage *image, double x, double y, double radius, bool logical_and, int maskValue) {
+
+    double dx, dy, r2, R2;
+
+    R2 = PS_SQR(radius);
+
+    for (int iy = 0; iy < image->numRows; iy++) {
+	for (int ix = 0; ix < image->numCols; ix++) {
+	    dx = ix + image->col0 - x;
+	    dy = iy + image->row0 - y;
+	    r2 = PS_SQR(dx) + PS_SQR(dy);
+	    if (r2 < R2) continue;
+	    if (logical_and) {
+		image->data.U8[iy][ix] &= maskValue;
+	    } else {
+		image->data.U8[iy][ix] |= maskValue;
+	    }
+	}
+    }
+}
+
Index: /trunk/psphot/src/psphot-utils.c
===================================================================
--- /trunk/psphot/src/psphot-utils.c	(revision 4250)
+++ /trunk/psphot/src/psphot-utils.c	(revision 4251)
@@ -135,7 +135,8 @@
     psFitsWriteImage (fits, NULL, image, 0, NULL);
     psFree (fits);
-    return;
+    return true;
 }
 
+# if (0)
 void pmSourceMaskSaturated (psSource *source, float saturate) {
 
@@ -150,86 +151,5 @@
     }
 }
-
-// mask the area contained by the region
-// the region is defined wrt the parent image
-void psImageMaskRegion (psImage *image, psRegion *region, bool logical_and, int maskValue) {
-
-    for (int iy = 0; iy < image->numRows; iy++) {
-	for (int ix = 0; ix < image->numCols; ix++) {
-	    if (ix + image->col0 <  region->x0) continue;
-	    if (ix + image->col0 >= region->x1) continue;
-	    if (iy + image->row0 <  region->y0) continue;
-	    if (iy + image->row0 >= region->y1) continue;
-	    if (logical_and) {
-		image->data.U8[iy][ix] &= maskValue;
-	    } else {
-		image->data.U8[iy][ix] |= maskValue;
-	    }
-	}
-    }
-}
-
-// mask the area not contained by the region
-// the region is defined wrt the parent image
-void psImageKeepRegion (psImage *image, psRegion *region, bool logical_and, int maskValue) {
-
-    for (int iy = 0; iy < image->numRows; iy++) {
-	for (int ix = 0; ix < image->numCols; ix++) {
-	    if (ix + image->col0 <  region->x0) goto maskit;
-	    if (ix + image->col0 >= region->x1) goto maskit;
-	    if (iy + image->row0 <  region->y0) goto maskit;
-	    if (iy + image->row0 >= region->y1) goto maskit;
-	    continue;
-	maskit:
-	    if (logical_and) {
-		image->data.U8[iy][ix] &= maskValue;
-	    } else {
-		image->data.U8[iy][ix] |= maskValue;
-	    }
-	}
-    }
-}
-
-// mask the area contained by the region
-// the region is defined wrt the parent image
-void psImageMaskCircle (psImage *image, double x, double y, double radius, bool logical_and, int maskValue) {
-
-    double R2 = PS_SQR(radius);
-
-    for (int iy = 0; iy < image->numRows; iy++) {
-	for (int ix = 0; ix < image->numCols; ix++) {
-	    dx = ix + image->col0 - x;
-	    dy = iy + image->row0 - y;
-	    r2 = PS_SQR(dx) + PS_SQR(dy);
-	    if (r2 > R2) continue;
-	    if (logical_and) {
-		image->data.U8[iy][ix] &= maskValue;
-	    } else {
-		image->data.U8[iy][ix] |= maskValue;
-	    }
-	}
-    }
-}
-
-// mask the area contained by the region
-// the region is defined wrt the parent image
-void psImageKeepCircle (psImage *image, double x, double y, double radius, bool logical_and, int maskValue) {
-
-    double R2 = PS_SQR(radius);
-
-    for (int iy = 0; iy < image->numRows; iy++) {
-	for (int ix = 0; ix < image->numCols; ix++) {
-	    dx = ix + image->col0 - x;
-	    dy = iy + image->row0 - y;
-	    r2 = PS_SQR(dx) + PS_SQR(dy);
-	    if (r2 < R2) continue;
-	    if (logical_and) {
-		image->data.U8[iy][ix] &= maskValue;
-	    } else {
-		image->data.U8[iy][ix] |= maskValue;
-	    }
-	}
-    }
-}
+# endif
 
 // create a subset of sources for the PS_SOURCE_PSFSTAR entries
Index: /trunk/psphot/src/psphot.c
===================================================================
--- /trunk/psphot/src/psphot.c	(revision 4250)
+++ /trunk/psphot/src/psphot.c	(revision 4251)
@@ -10,6 +10,4 @@
     psStats     *sky     = NULL;
 
-    psLogArguments (&argc, argv);
-    psTraceArguments (&argc, argv);
     config = load_args (&argc, argv);
 
@@ -31,13 +29,8 @@
 
     // use stellar objects SN > PSF_SN_LIM
-    psf = choose_psf_model (imdata, config, sources);
+    psf = choose_psf_model (config, sources);
 
-    if (0) { 
-      psImage *zap = psImageAlloc (imdata->image->numCols, imdata->image->numRows, PS_TYPE_F32);
-      psImageInit (zap, 0);
-      find_defects (zap, sources, config, psf);
-      DumpImage (zap, argv[2]);
-      exit (0);
-    }
+    // identify defects based on 1-pixel wide x and y direction 2nd moments
+    find_defects (zap, sources, config, psf);
 
     // test PSF on all except SATURATE and DEFECT (artifacts)
Index: /trunk/psphot/src/psphot.h
===================================================================
--- /trunk/psphot/src/psphot.h	(revision 4250)
+++ /trunk/psphot/src/psphot.h	(revision 4251)
@@ -41,4 +41,5 @@
 int get_argument (int argc, char **argv, char *arg);
 int remove_argument (int N, int *argc, char **argv);
+psMetadata *load_args (int *argc, char **argv);
 
 psArray *pmPeaksSubset(psArray *peaks, psF32 maxValue, const psRegion *valid);
@@ -66,5 +67,5 @@
 pmPSF *pmPSFAlloc (psModelType type);
 pmPSF_Test *pmPSF_TestAlloc (psArray *stars, char *modelName);
-pmPSF_Test *pmPSF_TestModel (psArray *stars, char *modelName);
+pmPSF_Test *pmPSF_TestModel (psArray *sources, char *modelName, float radius);
 
 bool pmPSFFromModels (pmPSF *psf, psArray *models, psVector *mask);
@@ -74,10 +75,10 @@
 
 // top-level psphot functions
-bool setup (psImage **image, psMetadata **header, psMetadata **config, char **argv);
-psStats *image_stats (psImage *image, psMetadata *config);
-psArray *find_peaks (psImage *image, psMetadata *config, psStats *sky);
-psArray *source_moments (psImage *image, psMetadata *config, psArray *allpeaks);
-pmPSF *choose_psf_model (psImage *image, psMetadata *config, psArray *sources);
-bool apply_psf_model (psImage *image, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
+psImageData *setup (psMetadata *config);
+psStats *image_stats (psImageData *imdata, psMetadata *config);
+psArray *find_peaks (psImageData *imdata, psMetadata *config, psStats *sky);
+psArray *source_moments (psImageData *imdata, psMetadata *config, psArray *allpeaks);
+pmPSF *choose_psf_model (psMetadata *config, psArray *sources);
+bool apply_psf_model (psImageData *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
 bool test_psf_scatter (psArray *sources);
 bool mark_psf_sources (psArray *sources, psMetadata *config);
@@ -90,2 +91,15 @@
 bool mark_psf_source (psSource *source, float shapeNsigma);
 bool find_defects (psImage *zapmask, psArray *sources, psMetadata *config, pmPSF *psf);
+
+bool pmSourceDefinePixels(psSource *mySource, const psImageData *imdata, psF32 x, psF32 y, psF32 Radius);
+bool pmSourceLocalSky_EAM (psSource *source, psStatsOptions statsOptions, psF32 Radius);
+psRegion *psRegionSquare (psF32 x, psF32 y, psF32 radius);
+
+// image mask functions
+void psImageMaskRegion (psImage *image, psRegion *region, bool logical_and, int maskValue);
+void psImageKeepRegion (psImage *image, psRegion *region, bool logical_and, int maskValue);
+void psImageMaskCircle (psImage *image, double x, double y, double radius, bool logical_and, int maskValue);
+void psImageKeepCircle (psImage *image, double x, double y, double radius, bool logical_and, int maskValue);
+
+# define AND true
+# define OR false
Index: /trunk/psphot/src/pspsf.c
===================================================================
--- /trunk/psphot/src/pspsf.c	(revision 4250)
+++ /trunk/psphot/src/pspsf.c	(revision 4251)
@@ -46,6 +46,10 @@
 }
 
-pmPSF_Test *pmPSF_TestModel (psArray *sources, char *modelName, float RADIUS) {
+pmPSF_Test *pmPSF_TestModel (psArray *sources, char *modelName, float RADIUS) 
+{
   
+    float x;
+    float y;
+    bool status;
     int Nflt = 0;
     int Npsf = 0;
@@ -102,8 +106,10 @@
 	}
 	psModel  *modelPSF = psModelFromPSF (modelFLT, test->psf); // set shape for this model
-
-	psImageKeepCircle (source->mask, RADIUS, OR, 0x80);
+	x = source->peak->x;
+	y = source->peak->y;
+
+	psImageKeepCircle (source->mask, x, y, RADIUS, OR, 0x80);
 	status = pmSourceFitModel (source, modelPSF, true);
-	psImageKeepCircle (source->mask, RADIUS, AND, 0x7f);
+	psImageKeepCircle (source->mask, x, y, RADIUS, AND, 0x7f);
 
 	// skip poor fits
Index: /trunk/psphot/src/setup.c
===================================================================
--- /trunk/psphot/src/setup.c	(revision 4250)
+++ /trunk/psphot/src/setup.c	(revision 4251)
@@ -13,6 +13,6 @@
     psImage *noise = NULL;
     psImage *mask = NULL;
-
-    psRegion  region = {0,0,0,0};	// a region representing the entire array
+    psRegion region = {0,0,0,0};	// a region representing the entire array
+    bool status;
 
     psTimerStart ("psphot");
@@ -24,5 +24,5 @@
 
     // read header and image data from INPUT 
-    char *input = psMetadataLookupSTR (&status, config, "INPUT");
+    char *input = psMetadataLookupPtr (&status, config, "INPUT");
     psFits *file = psFitsAlloc (input);
     header = psFitsReadHeader (header, file);
@@ -31,5 +31,4 @@
 
     // grab these values from the approrpiate location (image header if necessary)
-    bool  status   = false;
     float RDNOISE  = pmConfigLookupF32 (&status, config, header, "RDNOISE");
     float GAIN     = pmConfigLookupF32 (&status, config, header, "GAIN");
@@ -38,5 +37,5 @@
 
     // load the noise image if it is supplied, otherwise build from input
-    char *noiseName = psMetadataLookupSTR (&status, config, "NOISE");
+    char *noiseName = psMetadataLookupPtr (&status, config, "NOISE");
     if (status == true) {
       file = psFitsAlloc (noiseName);
@@ -44,14 +43,26 @@
       psFree (file);
     } else {
-      psScalar *value = psScalarAlloc (1.0 / GAIN, PS_TYPE_F64);
-      noise = BinaryOp (NULL, image, '/', value);
+      psScalar *value;
+
+      noise = psImageAlloc (image->numCols, image->numRows, PS_TYPE_F32);
+      for (int iy = 0; iy < image->numRows; iy++) {
+	for (int ix = 0; ix < image->numCols; ix++) {
+	  noise->data.F32[iy][ix] = image->data.F32[iy][ix] / GAIN + PS_SQR(RDNOISE/GAIN);
+	}
+      }
+
+      // XXX EAM it looks like psBinaryOp is broken for (image op scalar)
+      # if (0)
+      value = psScalarAlloc (1.0 / GAIN, PS_TYPE_F64);
+      noise = (psImage *) psBinaryOp (NULL, image, "/", value);
       psFree (value);
 
-      psScalar *value = psScalarAlloc (PS_SQR(RDNOISE/GAIN), PS_TYPE_F64);
-      noise = BinaryOp (noise, noise, '+', value);
+      value = psScalarAlloc (PS_SQR(RDNOISE/GAIN), PS_TYPE_F64);
+      noise = (psImage *) psBinaryOp (noise, noise, "+", value);
       psFree (value);
+      # endif
     }
 
-    char *maskName = psMetadataLookupSTR (&status, config, "MASK");
+    char *maskName = psMetadataLookupPtr (&status, config, "MASK");
     if (status == true) {
       file = psFitsAlloc (maskName);
@@ -67,5 +78,5 @@
     psRegion *keep = psRegionAlloc (XBORDER, -XBORDER, YBORDER, -YBORDER);
     keep           = psRegionForImage (keep, image, keep);
-    psImageKeepRegion (mask, keep, 0x01);
+    psImageKeepRegion (mask, keep, OR, 0x01);
 
     // mask the saturated pixels
@@ -87,4 +98,4 @@
     imdata->mask = mask;
 
-    return (true);
+    return (imdata);
 }
Index: /trunk/psphot/src/source_moments.c
===================================================================
--- /trunk/psphot/src/source_moments.c	(revision 4250)
+++ /trunk/psphot/src/source_moments.c	(revision 4251)
@@ -1,9 +1,8 @@
 # include "psphot.h"
 
-psArray *source_moments (psImageData *imdata, psMetadata *config, psArray *allpeaks) 
+psArray *source_moments (psImageData *imdata, psMetadata *config, psArray *peaks) 
 {
     bool     status  = false;
     psArray *sources = NULL;
-    psArray *peaks   = NULL;
 
     psTimerStart ("psphot");
Index: /trunk/psphot/src/subtract_psf_source.c
===================================================================
--- /trunk/psphot/src/subtract_psf_source.c	(revision 4250)
+++ /trunk/psphot/src/subtract_psf_source.c	(revision 4251)
@@ -25,4 +25,7 @@
   pmSourceSubModel (pixels, source->mask, model, false);
   model->params->data.F32[0] = sky;
+
+  // XXX adjust (enhance) the noise matrix a la dophot?
+
   return (true);
 }