Index: /tags/eam_tags/psphot_dev_01/psphot/.cvsignore
===================================================================
--- /tags/eam_tags/psphot_dev_01/psphot/.cvsignore	(revision 22334)
+++ /tags/eam_tags/psphot_dev_01/psphot/.cvsignore	(revision 22334)
@@ -0,0 +1,1 @@
+bin
Index: /tags/eam_tags/psphot_dev_01/psphot/Makefile
===================================================================
--- /tags/eam_tags/psphot_dev_01/psphot/Makefile	(revision 22334)
+++ /tags/eam_tags/psphot_dev_01/psphot/Makefile	(revision 22334)
@@ -0,0 +1,98 @@
+default: psphot
+help:
+	@echo "USAGE: make psphot"
+
+# CC      =       gcc -g -Wall -O2
+CC      =       gcc -g -Wall 
+SRC     =       src
+BIN     =       bin
+
+DESTBIN =       $(HOME)/src/bin/$(ARCH)
+
+LPSLIB  :=      $(shell pslib-config --libs)
+IPSLIB  :=      $(shell pslib-config --cflags)
+
+INCS	= 	$(IPSLIB)
+LIBS	= 	-lpsmodule $(LPSLIB)
+CFLAGS	=	$(INCS) -std=c99
+LFLAGS	=	$(LIBS) 
+
+PSPHOT = \
+$(SRC)/psphot.$(ARCH).o            \
+$(SRC)/psphotArguments.$(ARCH).o   \
+$(SRC)/psphotSetup.$(ARCH).o	   \
+$(SRC)/psphotImageStats.$(ARCH).o  \
+$(SRC)/psphotSourceStats.$(ARCH).o \
+$(SRC)/psphotChoosePSF.$(ARCH).o   \
+$(SRC)/psphotApplyPSF.$(ARCH).o	   \
+$(SRC)/psphotFitGalaxies.$(ARCH).o \
+$(SRC)/psphotOutput.$(ARCH).o      \
+$(SRC)/psphotMarkPSF.$(ARCH).o     \
+$(SRC)/psphotSubtractPSF.$(ARCH).o \
+$(SRC)/pspsf.$(ARCH).o             \
+$(SRC)/psLibUtils.$(ARCH).o	   \
+$(SRC)/psModulesUtils.$(ARCH).o	   \
+$(SRC)/pmSourceUtils.$(ARCH).o	   \
+$(SRC)/pmPeaksSigmaLimit.$(ARCH).o \
+$(SRC)/psLine.$(ARCH).o		   \
+$(SRC)/psPolynomials.$(ARCH).o	   \
+$(SRC)/psMinimize.$(ARCH).o	   \
+$(SRC)/psImageData.$(ARCH).o	   \
+$(SRC)/psphotSortBySN.$(ARCH).o
+
+FITSOURCE = \
+$(SRC)/fitsource.$(ARCH).o \
+$(SRC)/onesource.$(ARCH).o \
+$(SRC)/fs_args.$(ARCH).o \
+$(SRC)/psphot-utils.$(ARCH).o \
+$(SRC)/psPolynomials.$(ARCH).o \
+$(SRC)/psUtils.$(ARCH).o \
+$(SRC)/setup.$(ARCH).o \
+$(SRC)/LocalSky.$(ARCH).o
+
+psphot: $(BIN)/psphot.$(ARCH)
+$(BIN)/psphot.$(ARCH) : $(PSPHOT)
+$(PSPHOT): $(SRC)/psphot.h
+
+fitsource: $(BIN)/fitsource.$(ARCH)
+$(BIN)/fitsource.$(ARCH) : $(FITSOURCE)
+$(FITSOURCE): $(SRC)/psphot.h
+
+INSTALL = psphot
+
+# dependancy rules for binary code #########################
+.PRECIOUS: %.$(ARCH).o
+.PRECIOUS: $(BIN)/%.$(ARCH)
+
+%.$(ARCH).o : %.c
+	$(CC) $(CFLAGS) -c $< -o $@
+
+$(BIN)/%.$(ARCH) : $(SRC)/%.$(ARCH).o
+	@if [ ! -d $(BIN) ]; then mkdir -p $(BIN); fi
+	$(CC) $^ -o $@ $(LFLAGS)
+	@echo "done with $@"
+
+$(DESTBIN)/%: $(BIN)/%.$(ARCH)
+	@if [ ! -d $(DESTBIN) ]; then mkdir -p $(DESTBIN); fi
+	rm -f $(DESTBIN)/$*
+	cp $(BIN)/$*.$(ARCH) $(DESTBIN)/$*
+
+$(INSTALL): % : $(BIN)/%.$(ARCH)
+
+%.clean :
+	rm -f $(BIN)/$*.$(ARCH)
+	rm -f $(SRC)/*.$(ARCH).o
+
+%.install:
+	make $(DESTBIN)/$*
+
+# utilities #################################################
+
+install:
+	for i in $(INSTALL); do make $$i.install; done
+
+clean:	
+	rm -f $(BIN)/*.$(ARCH)
+	rm -f `find . -name "*.o"`
+	rm -f `find . -name "*~"`
+	rm -f `find . -name "#*"`
Index: /tags/eam_tags/psphot_dev_01/psphot/doc/notes.txt
===================================================================
--- /tags/eam_tags/psphot_dev_01/psphot/doc/notes.txt	(revision 22334)
+++ /tags/eam_tags/psphot_dev_01/psphot/doc/notes.txt	(revision 22334)
@@ -0,0 +1,176 @@
+
+Notes on psphot
+
+2005.09.05
+
+  I have a working version of PSPhot which handles PSF and galaxy
+  models.  I am tagging the module with an alpha version number. this
+  version compiles with psLib-0.5.0, at least my version of it.  My
+  current goal is to make it work with the current psLib head (nearly
+  0.7.0) with a mimimum of _EAM versions of functions.
+
+  the tag for the working version is dev_01
+
+2005.06.04
+
+  progress on psphot has moved along well.  At this point, the process
+  loads the image, finds the peaks, determines a best PSF model,
+  identifies objects which are consistent with that model, and
+  attempts to fit a galaxy model to the objects which are not
+  consistent with the model. 
+
+  a variety of minor issues remain, as well as some major issues. 
+
+  - use a proper noise image to keep the fits honest after other
+    objects have been subtracted.  
+
+  - define the masks in terms of a single mask image.  this image
+    could be provided initially by the user.   
+
+  - add noise enhancement (couple with mask marking subtraction)
+
+  - 
+
+2005.04.12
+
+- psPeak, psSource, etc: should be pmPeak, pmSource, etc
+- pmCullPeaks: should be pmCullImagePeaks
+
+- does pmSourceMoments need an image argument?
+
+- psLibInit / p_psTimeInit
+
+  we should not have to know about astronomy time concepts (ie, have a
+  timeConfig file) in order to use the time functions for basic work.
+  This means having a PS_TIME_UNIX which just works with the UNIX
+  clock and avoids all initialization issues.
+
+- psImageStats is much too slow for basic median
+  - this function does a complete sort on all pixel values.  we need
+  to use the histogram method (spec a different type of median)
+
+- psImageStats, ROBUST_MEDIAN fails
+- there are errors in psVectorStats for ROBUST MEAN, STDEV (& median?)
+  * this is ill-defined.  revise the ADD: too sensitive to the
+  binning, smoothing scales 
+
+- pmFindImagePeaks is not finding any peaks:
+  typo in interior row section: see bug 359
+  after this fix: PASS
+
+* p_psGetRowVectorFromImage : PASS
+* psGetRowVectorFromImage : PASS 
+- both have about the same time (~15 us +/- 5us) for 2048 pix on alala
+
+* pmFindVectorPeaks : PASS (returned row data in data.U32)
+
+- # include <string.h> : needed in psMetadata.c
+
+* MyListAddPeak / pmPeakAlloc are inconsistent wrt col,row vs x,y
+
+- why does psMetadataItemAllocV allocate a string of length
+  MAX_STRING_LENGTH for the comment, then use strncpy?
+
+  should be:
+
+    // set metadata item comment
+    if (comment == NULL) {
+        // Per SDRS, null isn't allowed, must use "" instead
+        metadataItem->comment = psStringCopy ("");
+    } else {
+        metadataItem->comment = psStringCopy (comment);
+    }
+
+psMetadataAdd -> psMetadataAddItem is all confused:
+ - we end up with two psMetadataItemAlloc calls, when only one should
+   be used
+ - psMetadataItemAlloc was doing the wrong thing if the incoming data
+   was NULL
+ - I think the logic in this block is wrong as well.  
+
+ * I fixed this in psMetadata.c
+
+2005.04.14 : 
+
+ - psImageStats is still much too slow
+
+ - pmObjects : changed to using psArray to carry peaks / sources,
+   rather than psList
+
+ - added function pmPeaksSubset to replace pmCullPeaks
+
+ - pmSourceLocalSky: added min/max fncs to force subimage to stay in
+   image
+
+ - Sx, Sy can go negative: have forced limit to 0,0
+
+ - pmSourceRoughClass: implemented sigmaX, sigmaY search for stellar
+   clump (uses pmComparePeakDescend)
+
+   - the search ignores the 0,0 pixel
+
+   - pmSourceRoughClass uses all detected peaks.  allow an exclusion
+     for peak > max, peak < min?
+
+   - need to push clump stats on the metadata
+
+   - the source classes are assigned so they are mutually exclusive 
+
+   - change PS_SOURCE_OTHER to STAR? 
+
+ - pmSourceSetPixelCircle:
+
+   - name should be pmSourceSetPixelsCircle (geometry could be an
+     option...)
+
+   - I cleaned the defined region to match the convensions of
+     pmSourceLocalSky
+
+   - truncate and force subimage to lie on image
+
+   - CheckRadius2 seems like an inefficient coding
+
+   - valid pixels have mask value of 1 (inverse of other funcs)
+
+ - pmSourceFitModel:
+
+   - num iterations is much too large
+
+   - 
+
+mrq:
+
+ - allocate alpha, beta arrays
+ - mrq2dcof (pars, alpha, beta)
+ - make guess pars', alpha', beta'
+ - mrq2dcof (pars', alpha', beta')
+ - if (chisq < old chisq), keep new guess
+
+
+2005.04.19:
+
+- psMinimizeLMChi2 has some serious problems.  I re-wrote it starting
+  from the mrqmin example in ohana.   after a variety little bugs, it
+  seems to work quite well.  the program 'modeltest.c' runs two
+  different gaussians through psMinimizeLMChi2, and gets the right
+  answer quite quickly.  The first pass at the implementation had a
+  variety of problems.  to avoid any worry about errors in the
+  psMatrix code, I implemented psGaussJordan, basing the code on the
+  elixir gaussj.c code.  this worked fine, and let me find the errors
+  elsewhere in the code. the LUD version now works just the same as
+  the GaussJ version.
+
+
+2005.04.21
+
+- psMinimizeLMChi2 does an excellent job now in most cases.  I am able
+  to run pmSourceFitModel on the PSF stars quite well.  I added in the
+  (nearly) correct errors (actually, just using sqrt(N), limiting to
+  N>=1).  I see one problem object, which is quite faint, and does not
+  converge well: I get very large values for x,y, Sx, Sy, etc.  I need
+  to think about constraining these fits from running way off scale.
+
+  fit:   0.123879 sec for 277 stars (gauss)
+  fit:   0.155940 sec for 277 stars (pgauss) (!)
+  fit:   0.153292 sec for 277 stars
+  (0.5 msec per star)
Index: /tags/eam_tags/psphot_dev_01/psphot/doc/psphot.txt
===================================================================
--- /tags/eam_tags/psphot_dev_01/psphot/doc/psphot.txt	(revision 22334)
+++ /tags/eam_tags/psphot_dev_01/psphot/doc/psphot.txt	(revision 22334)
@@ -0,0 +1,5 @@
+
+psphot is a stand-alone program which implements the Pan-STARRS IPP
+object detection and analysis functions for using with individual
+images.  
+
Index: /tags/eam_tags/psphot_dev_01/psphot/src/fitsource.c
===================================================================
--- /tags/eam_tags/psphot_dev_01/psphot/src/fitsource.c	(revision 22334)
+++ /tags/eam_tags/psphot_dev_01/psphot/src/fitsource.c	(revision 22334)
@@ -0,0 +1,21 @@
+# include "psphot.h"
+
+int main (int argc, char **argv) {
+
+    psMetadata  *config  = NULL;		// user-provided configuration information
+    psImageData *imdata  = NULL;
+    float x, y;
+
+    config = fs_args (&argc, argv);
+
+    // load input data (image and config)
+    // create or load mask and noise images
+    imdata = setup (config);
+
+    x = atof (argv[4]);
+    y = atof (argv[5]);
+
+    onesource (imdata, config, argv[3], x, y);
+
+    exit (0);
+}
Index: /tags/eam_tags/psphot_dev_01/psphot/src/fs_args.c
===================================================================
--- /tags/eam_tags/psphot_dev_01/psphot/src/fs_args.c	(revision 22334)
+++ /tags/eam_tags/psphot_dev_01/psphot/src/fs_args.c	(revision 22334)
@@ -0,0 +1,49 @@
+# include "psphot.h"
+
+psMetadata *fs_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 
+
+  char *mask = NULL;
+  if ((N = get_argument (*argc, argv, "--mask"))) {
+    remove_argument (N, argc, argv);
+    mask = psStringCopy (argv[N]);
+    remove_argument (N, argc, argv);
+  }
+
+  char *noise = NULL;
+  if ((N = get_argument (*argc, argv, "--noise"))) {
+    remove_argument (N, argc, argv);
+    noise = psStringCopy (argv[N]);
+    remove_argument (N, argc, argv);
+  }
+
+  if (*argc != 6) fs_usage ();
+
+  psMetadata *config = psMetadataParseConfig (NULL, &Nfail, argv[2], FALSE);
+
+  psMetadataAddStr (config, PS_LIST_HEAD, "INPUT",  "", argv[1]);
+  if (mask != NULL) {
+    psMetadataAddStr (config, PS_LIST_HEAD, "MASK", "", mask);
+  }
+  if (noise != NULL) {
+    psMetadataAddStr (config, PS_LIST_HEAD, "NOISE", "", noise);
+  }
+  return (config);
+}
+
+int fs_usage () {
+
+    fprintf (stderr, "USAGE: psphot (image.fits) (config) (model) (x) (y)\n");
+    fprintf (stderr, "options: \n");
+    fprintf (stderr, "  --mask  (filename)\n");
+    fprintf (stderr, "  --noise (filename)\n");
+    exit (2);
+}
Index: /tags/eam_tags/psphot_dev_01/psphot/src/onesource.c
===================================================================
--- /tags/eam_tags/psphot_dev_01/psphot/src/onesource.c	(revision 22334)
+++ /tags/eam_tags/psphot_dev_01/psphot/src/onesource.c	(revision 22334)
@@ -0,0 +1,43 @@
+# include "psphot.h"
+
+bool onesource (psImageData *imdata, psMetadata *config, char *modelName, float x, float y) {
+
+    bool status;
+    int modelType;
+    float sky, obsMag, fitMag;
+
+    modelType   = psModelSetType (modelName);
+    if (modelType < 0) psAbort ("fitsource", "unknown model %s", modelName);
+
+    float INNER  = psMetadataLookupF32 (&status, config, "INNER_RADIUS");
+    float OUTER  = psMetadataLookupF32 (&status, config, "OUTER_RADIUS");
+    float MRAD   = psMetadataLookupF32 (&status, config, "GAL_MOMENTS_RADIUS");
+    float RADIUS = psMetadataLookupF32 (&status, config, "FIT_RADIUS");
+
+    psSource *source = pmSourceAlloc();
+    source->peak = pmPeakAlloc (x, y, 0, 0);
+    pmSourceDefinePixels (source, imdata, x, y, OUTER);
+
+    status = pmSourceLocalSky_EAM (source, PS_STAT_SAMPLE_MEDIAN, INNER);
+    status = pmSourceMoments_EAM (source, MRAD);
+    source->peak->counts = source->moments->Peak;
+
+    psModel *model = pmSourceModelGuess (source, modelType); 
+    fprintf (stderr, "guess slope: %f\n", model->params->data.F32[7]);
+
+    psImageKeepCircle (source->mask, x, y, RADIUS, OR, 0x80);
+    status = pmSourceFitModel_EAM (source, model, false);
+
+    pmSourcePhotometry (&fitMag, &obsMag, model, source->pixels, source->mask);
+    fprintf (stderr, "ap: %f, fit: %f, apmifit: %f\n", obsMag, fitMag, obsMag - fitMag);
+
+    // subtract object, leave local sky
+    sky = model->params->data.F32[0];
+    model->params->data.F32[0] = 0;
+    pmSourceSubModel (source->pixels, source->mask, model, false);
+    model->params->data.F32[0] = sky;
+    
+    DumpImage (source->pixels, "resid.fits");
+    DumpImage (source->mask, "mask.fits");
+    return true;
+}
Index: /tags/eam_tags/psphot_dev_01/psphot/src/pmPeaksSigmaLimit.c
===================================================================
--- /tags/eam_tags/psphot_dev_01/psphot/src/pmPeaksSigmaLimit.c	(revision 22334)
+++ /tags/eam_tags/psphot_dev_01/psphot/src/pmPeaksSigmaLimit.c	(revision 22334)
@@ -0,0 +1,43 @@
+# include "psphot.h"
+
+psArray *pmPeaksSigmaLimit (psImageData *imdata, psMetadata *config, psStats *sky) 
+{
+    bool  status = false;
+    float NSIGMA;
+    float SIGMA;
+    float threshold;
+
+    // smooth the image 
+
+    psTimerStart ("psphot");
+
+    SIGMA  = psMetadataLookupF32 (&status, config, "PEAKS_SMOOTH_SIGMA");
+    NSIGMA = psMetadataLookupF32 (&status, config, "PEAKS_SMOOTH_NSIGMA");
+
+    psImage *smooth = psImageCopy (NULL, imdata->image, PS_TYPE_F32);
+    psImageSmooth (smooth, SIGMA, NSIGMA);
+    psLogMsg ("psphot", 4, "smooth: %f sec\n", psTimerMark ("psphot"));
+
+    psTimerStart ("psphot");
+
+    // set peak threshold
+    NSIGMA    = psMetadataLookupF32 (&status, config, "PEAKS_NSIGMA_LIMIT");
+    threshold = NSIGMA*sky->sampleStdev + sky->sampleMean;
+    psLogMsg ("psphot", 3, "threshold: %f DN\n", threshold);
+
+    // find the peaks in the smoothed image 
+    psArray *peaks = pmFindImagePeaks (smooth, threshold);
+    if (peaks == NULL) psAbort ("find peaks", "no peaks found");
+    psFree (smooth);
+
+    // make this optional
+    // DumpPeaks (peaks, "peaks.dat");
+
+    psLogMsg ("psphot", 3, "%d peaks: %f sec\n", peaks->n, psTimerMark ("psphot"));
+
+    return (peaks);
+}
+
+// In this function, we smooth the image, then search for the peaks 
+// Should we also subtract a super-binned image? (as an option?)
+// We need to gracefully handle no source detections
Index: /tags/eam_tags/psphot_dev_01/psphot/src/pmSourceUtils.c
===================================================================
--- /tags/eam_tags/psphot_dev_01/psphot/src/pmSourceUtils.c	(revision 22334)
+++ /tags/eam_tags/psphot_dev_01/psphot/src/pmSourceUtils.c	(revision 22334)
@@ -0,0 +1,51 @@
+# include "psphot.h"
+
+bool pmSourceDefinePixels(psSource *mySource, 
+			  const psImageData *imdata,
+			  psF32 x, 
+			  psF32 y,
+			  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);
+
+    //  use these when psImageSubset is updated
+    //  psImage *subImage = psImageSubset(imdata->image, srcRegion);
+    //  psImage *subImageMask = psImageSubset(imdata->mask, srcRegion);
+    //  psImage *subImageNoise = psImageSubset(imdata->noise, srcRegion);
+
+    mySource->pixels = psImageSubset(imdata->image, srcRegion->x0, srcRegion->y0, srcRegion->x1, srcRegion->y1);
+    mySource->noise  = psImageSubset(imdata->noise, srcRegion->x0, srcRegion->y0, srcRegion->x1, srcRegion->y1);
+    mySource->mask   = psImageSubset(imdata->mask,  srcRegion->x0, srcRegion->y0, srcRegion->x1, srcRegion->y1);
+
+    psFree (srcRegion);
+    return(mySource);
+}
+
+bool pmSourcePhotometry (float *fitMag, float *obsMag, psModel *model, psImage *image, psImage *mask) {
+
+    float obsSum = 0;
+    float fitSum = 0;
+    float sky = model->params->data.F32[0];
+
+    psModelFlux modelFluxFunc = psModelFlux_GetFunction (model->type);
+    fitSum = modelFluxFunc (model->params);
+
+    for (int ix = 0; ix < image->numCols; ix++) {
+	for (int iy = 0; iy < image->numRows; iy++) {
+	    if (mask->data.U8[iy][ix]) continue;
+	    obsSum += image->data.F32[iy][ix] - sky;
+	    // fitSum += psModelEval (model, image, ix, iy) - sky;
+	}
+    }
+    if (obsSum <= 0) return false;
+    if (fitSum <= 0) return false;
+
+    *fitMag = -2.5*log10(fitSum);
+    *obsMag = -2.5*log10(obsSum);
+    return (true);
+} 
+
Index: /tags/eam_tags/psphot_dev_01/psphot/src/psImageData.c
===================================================================
--- /tags/eam_tags/psphot_dev_01/psphot/src/psImageData.c	(revision 22334)
+++ /tags/eam_tags/psphot_dev_01/psphot/src/psImageData.c	(revision 22334)
@@ -0,0 +1,24 @@
+# include "psphot.h"
+
+static void psImageDataFree (psImageData *imdata) {
+
+  if (imdata == NULL) return;
+
+  psFree (imdata->image);
+  psFree (imdata->header);
+  psFree (imdata->noise);
+  psFree (imdata->mask);
+  return;
+}
+
+psImageData *psImageDataAlloc (psImage *image, psImage *noise, psImage *mask, psMetadata *header) {
+  
+    psImageData *imdata = psAlloc(sizeof(psImageData));
+    imdata->image = image;
+    imdata->header = header;
+    imdata->noise = noise;
+    imdata->mask = mask;
+
+    p_psMemSetDeallocator(imdata, (psFreeFcn) psImageDataFree);
+    return (imdata);
+}
Index: /tags/eam_tags/psphot_dev_01/psphot/src/psLibUtils.c
===================================================================
--- /tags/eam_tags/psphot_dev_01/psphot/src/psLibUtils.c	(revision 22334)
+++ /tags/eam_tags/psphot_dev_01/psphot/src/psLibUtils.c	(revision 22334)
@@ -0,0 +1,456 @@
+# include "psphot.h"
+
+// XXX EAM : these utility functions should be added back into PSLib
+
+static psHash *timers = NULL;
+
+bool psTimerClear (char *name) {
+
+  bool status;
+
+  if (name == NULL) return false;
+
+  status = psHashRemove (timers, name);
+  return (status);
+}
+
+void psTimerFree () {
+
+  psFree (timers);
+  p_psTimeFinalize();
+  return;
+}
+
+// start/restart a named timer
+bool psTimerStart (char *name) {
+
+    psTime *start;
+
+    if (timers == NULL) timers = psHashAlloc (16);
+
+    start = psTimeGetTime (PS_TIME_UTC);
+    psHashAdd (timers, name, start);
+    psFree (start); 
+    return (TRUE);
+}
+
+// get current elapsed time on named timer
+psF64 psTimerMark (char *name) {
+
+    psTime *start;
+    psTime *mark;
+    psF64   delta;
+
+    if (timers == NULL) return (0);
+
+    start = psHashLookup (timers, name);
+    if (start == NULL) return (0);
+
+    mark = psTimeGetTime (PS_TIME_UTC);
+    delta = psTimeDelta (mark, start);
+    psFree (mark);
+    // psFree (start); -- XXX is psHashLookup not psMemCopying?
+
+    return (delta);
+}
+
+// we have log levels 1 (Error), 2 (Warning), 3 (Info), 4 (Details), 5 (Minutiae)
+// 2 = default, -v = 3, -vv = 4, -vvv = 5 
+psS32 psLogArguments (int *argc, char **argv) {
+  
+    int N, level;
+
+    // default log level is 2
+    level = 2;
+
+    // set in order, so that -vvv overrides -vv overrides -v
+    if ((N = psArgumentGet (*argc, argv, "-v"))) {
+	psArgumentRemove (N, argc, argv);
+	level = 3;
+    }
+    if ((N = psArgumentGet (*argc, argv, "-vv"))) {
+	psArgumentRemove (N, argc, argv);
+	level = 4;
+    }
+    if ((N = psArgumentGet (*argc, argv, "-vvv"))) {
+	psArgumentRemove (N, argc, argv);
+	level = 5;
+    }
+
+    if ((N = psArgumentGet (*argc, argv, "-logfmt"))) {
+	if (*argc < N + 2) {
+	    psAbort ("psLogArguments", "USAGE: -logfmt (format)");
+	}
+	psArgumentRemove (N, argc, argv);
+	psLogSetFormat (argv[N]); // XXX EAM : this function should return an error if the log format is invalid
+	psArgumentRemove (N, argc, argv);
+    }
+
+    // set the level, return the level
+    psLogSetLevel (level);
+    return (level);
+}
+
+// set trace levels by facility
+psS32 psTraceArguments (int *argc, char **argv) {
+  
+    int N;
+
+    // argument format is: -trace (facil) (level) 
+    while ((N = psArgumentGet (*argc, argv, "-trace"))) {
+	if (*argc < N + 3) {
+	    psAbort ("psTraceArguments", "USAGE: -trace (facility) (level)");
+	}
+	psArgumentRemove (N, argc, argv);
+	psTraceSetLevel (argv[N], atoi(argv[N+1]));
+	psArgumentRemove (N, argc, argv);
+	psArgumentRemove (N, argc, argv);
+    }
+    if ((N = psArgumentGet (*argc, argv, "-trace-levels"))) {
+	psTracePrintLevels ();
+	exit (2);
+    }
+    return (TRUE);
+}
+
+// find the location of the specified argument
+int psArgumentGet (int argc, char **argv, char *arg) {
+
+    int i;
+
+    for (i = 0; i < argc; i++) {
+	if (!strcmp(argv[i], arg))
+	    return (i);
+    }
+  
+    return ((int)NULL);
+}
+
+// remove the specified argument (by location)
+int psArgumentRemove (int N, int *argc, char **argv) {
+
+    int i;
+
+    if ((N != (int)NULL) && (N != 0)) {
+	(*argc)--;
+	for (i = N; i < *argc; i++) {
+	    argv[i] = argv[i+1];
+	}
+    }
+
+    return (N);
+}
+
+// alternate implementation of this function from pmObjects.c
+psVector *psGetRowVectorFromImage(psImage *image,
+				  psU32 row)
+{
+    PS_IMAGE_CHECK_NULL(image, NULL);
+    PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, NULL);
+
+    psVector *tmpVector = psVectorAlloc(image->numCols, PS_TYPE_F32);
+    memcpy (tmpVector->data.F32, image->data.F32[row], image->numCols*sizeof(psF32));
+    return(tmpVector);
+}
+
+// extract config informatin from config data or from image header, if specified
+psF32 pmConfigLookupF32 (bool *status, psMetadata *config, psMetadata *header, char *name) {
+
+    char *source;
+    char *keyword;
+    psF32 value;
+    psMetadataItem *item;
+
+    // look for the entry in the config collection
+    item = psMetadataLookup (config, name);
+    if (item == NULL) {
+	psLogMsg ("pmConfigLookupF32", 2, "no key %s in config data, trying as header keyword", name);
+	value = psMetadataLookupF32 (status, header, name);
+	if (*status == false) {
+	    psAbort ("pmConfigLookupF32", "no key %s in header", name);
+	}
+	*status = true;
+	return (value);
+    }	
+
+    // I'm either expecting a string, with the name "HD:keyword"...
+    if (item->type == PS_META_STR) {
+	source = item->data.V;
+	if (!strncasecmp (source, "HD:", 3)) {
+	    keyword = &source[3]; 
+	    value = psMetadataLookupF32 (status, header, keyword);
+	    if (*status == false) {
+		psAbort ("pmConfigLookupF32", "no key %s in config", name);
+	    }
+	    *status = true;
+	    // psFree (item);
+	    return (value);
+	}	
+    }
+
+    //  ... or a value (F32?)
+    if (item->type == PS_META_F32) {
+	value = item->data.F32;
+	// psFree (item);
+	return (value);
+    }
+
+    *status = false;
+    psError(PS_ERR_UNKNOWN, true, "invalid item");
+    return (0);
+}
+
+void psImageSmooth (psImage *image, float sigma, float Nsigma) {
+
+    int Nx, Ny, Npixel, Nrange;
+    float factor, g, s;
+    psVector *temp;
+
+    // relevant terms
+    Nrange = sigma*Nsigma + 0.5;
+    Npixel = 2*Nrange + 1;
+    factor = -0.5/(sigma*sigma);
+
+    Nx = image->numCols;
+    Ny = image->numRows;
+
+    // generate gaussian
+    psVector *gaussnorm = psVectorAlloc (Npixel, PS_TYPE_F32);
+    for (int i = -Nrange; i < Nrange + 1; i++) {
+	gaussnorm->data.F32[i+Nrange] = exp (factor*i*i);
+    }
+    psF32 *gauss = &gaussnorm->data.F32[Nrange];
+
+    // smooth in X direction
+    temp = psVectorAlloc (Nx, PS_TYPE_F32);
+    for (int j = 0; j < Ny; j++) {
+	psF32 *vi = image->data.F32[j];
+	psF32 *vo = temp->data.F32;
+	for (int i = 0; i < Nx; i++) {
+	    g = s = 0;
+	    for (int n = -Nrange; n < Nrange + 1; n++) {
+		if (i+n < 0) continue;
+		if (i+n >= Nx) continue;
+		s += gauss[n]*vi[i+n];
+		g += gauss[n];
+	    }
+	    vo[i] = s / g;
+	}
+	memcpy (image->data.F32[j], temp->data.F32, Nx*sizeof(psF32));
+    }
+    psFree (temp);
+
+    // smooth in Y direction
+    temp = psVectorAlloc (image->numRows, PS_TYPE_F32);
+    for (int i = 0; i < Nx; i++) {
+	psF32  *vo = temp->data.F32;
+	psF32 **vi = image->data.F32;
+	for (int j = 0; j < Ny; j++) {
+	    g = s = 0;
+	    for (int n = -Nrange; n < Nrange + 1; n++) {
+		if (j+n < 0) continue;
+		if (j+n >= Ny) continue;
+		s += gauss[n]*vi[j+n][i];
+		g += gauss[n];
+	    }
+	    vo[j] = s / g;
+	}
+	// replace temp in image
+	for (int j = 0; j < Ny; j++) {
+	    vi[j][i] = vo[j];
+	}
+    }
+    psFree (temp);
+    psFree (gaussnorm);
+}
+    
+bool psImageInit (psImage *image,...) {
+
+  va_list argp;
+  psU8  vU8;
+  psF32 vF32;
+  psF64 vF64;
+
+  if (image == NULL) return (false);
+
+  va_start (argp, image);
+
+  switch (image->type.type) {
+    case PS_TYPE_U8:
+      vU8 = va_arg (argp, psU32);
+
+      for (int iy = 0; iy < image->numRows; iy++) {
+	for (int ix = 0; ix < image->numCols; ix++) {
+	  image->data.U8[iy][ix] = vU8;
+	}
+      }
+      break;
+
+    case PS_TYPE_F32:
+      vF32 = va_arg (argp, psF64);
+      
+      for (int iy = 0; iy < image->numRows; iy++) {
+	for (int ix = 0; ix < image->numCols; ix++) {
+	  image->data.F32[iy][ix] = vF32;
+	}
+      }
+      return (true);
+
+    case PS_TYPE_F64:
+      vF64 = va_arg (argp, psF64);
+
+      for (int iy = 0; iy < image->numRows; iy++) {
+	for (int ix = 0; ix < image->numCols; ix++) {
+	  image->data.F64[iy][ix] = vF64;
+	}
+      }
+      return (true);
+
+    default:
+      psAbort ("psphot.psUtils", "datatype %d not defined in psImageInit\n", image->type);
+      return (false);
+  }
+  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;
+	    }
+	}
+    }
+}
+
+psVector *psVectorCreate (double lower, double upper, double delta, psElemType type) {
+
+  int nBin = (upper - lower) / delta;
+
+  psVector *out = psVectorAlloc (nBin, type);
+  
+  for (int i = 0; i < nBin; i++) {
+    out->data.F64[i] = lower + i * delta;
+  }
+
+  return (out);
+}
Index: /tags/eam_tags/psphot_dev_01/psphot/src/psLine.c
===================================================================
--- /tags/eam_tags/psphot_dev_01/psphot/src/psLine.c	(revision 22334)
+++ /tags/eam_tags/psphot_dev_01/psphot/src/psLine.c	(revision 22334)
@@ -0,0 +1,46 @@
+# include "psphot.h"
+
+static void psLineFree (psLine *line) {
+
+  if (line == NULL) return;
+
+  psFree (line->line);
+  return;
+}
+
+// allocate a psLine structrue
+psLine *psLineAlloc (int Nline) {
+
+  psLine *line;
+  line = psAlloc (sizeof(psLine));
+  p_psMemSetDeallocator(line, (psFreeFcn) psLineFree);
+
+  line->Nline = 0;
+  line->NLINE = Nline;
+  line->line = psAlloc (Nline);
+  return (line);
+}
+
+bool psLineInit (psLine *line) {
+  if (line == NULL) return (false);
+  line->Nline = 0;
+  return (true);
+}
+
+bool psLineAdd (psLine *line, char *format, ...) {
+
+    int Nchar;
+    va_list ap;
+
+    if (line == NULL) return (false);
+
+    int nMax = line->NLINE - line->Nline;
+
+    va_start (ap, format);
+    Nchar = vsnprintf (&line->line[line->Nline], nMax, format, ap);
+    line->Nline += PS_MIN (nMax - 1, Nchar);
+
+    if (Nchar >= nMax) return (false);
+    return (true);
+}
+
Index: /tags/eam_tags/psphot_dev_01/psphot/src/psMinimize.c
===================================================================
--- /tags/eam_tags/psphot_dev_01/psphot/src/psMinimize.c	(revision 22334)
+++ /tags/eam_tags/psphot_dev_01/psphot/src/psMinimize.c	(revision 22334)
@@ -0,0 +1,250 @@
+# include "psphot.h"
+
+// XXX EAM : this file contains my alternate implementations of psMinimizeLMChi2
+
+// XXX EAM this implementation of MinLM includes limits on params & dparams
+psBool psMinimizeLMChi2_EAM(psMinimization *min,
+			    psImage *covar,
+			    psVector *params,
+			    const psVector *paramMask,
+			    const psArray *x,
+			    const psVector *y,
+			    const psVector *yErr,
+			    psMinimizeLMChi2Func func)
+{
+    PS_PTR_CHECK_NULL(min, NULL);
+    PS_VECTOR_CHECK_NULL(params, NULL);
+    PS_VECTOR_CHECK_EMPTY(params, NULL);
+    PS_PTR_CHECK_NULL(x, NULL);
+    PS_VECTOR_CHECK_NULL(y, NULL);
+    PS_VECTOR_CHECK_EMPTY(y, NULL);
+    PS_VECTOR_CHECK_SIZE_EQUAL(x, y, NULL);
+    PS_PTR_CHECK_NULL(func, NULL);
+
+    // this function has test and current values for several things
+    // the current best value is in lower case
+    // the next guess value is in upper case
+
+    // allocate internal arrays (current vs Guess)
+    psImage *alpha   = psImageAlloc (params->n, params->n, PS_TYPE_F64);
+    psImage *Alpha   = psImageAlloc (params->n, params->n, PS_TYPE_F64);
+    psVector *beta   = psVectorAlloc (params->n, PS_TYPE_F64);
+    psVector *Beta   = psVectorAlloc (params->n, PS_TYPE_F64);
+    psVector *Params = psVectorAlloc (params->n, PS_TYPE_F32);
+    psVector *dy     = NULL;
+    psF64 Chisq = 0.0;
+    psF64 lambda = 0.001;
+
+    psVector *beta_lim = NULL;
+    psVector *param_min = NULL;
+    psVector *param_max = NULL;
+
+    // if we are provided a covar image, we expect to find these three vectors in first three rows
+    if (covar != NULL) {
+	beta_lim  = psVectorAlloc (params->n, PS_TYPE_F32);
+	param_min = psVectorAlloc (params->n, PS_TYPE_F32);
+	param_max = psVectorAlloc (params->n, PS_TYPE_F32);
+	for (int i = 0; i < params->n; i++) {
+	    beta_lim->data.F32[i] = covar->data.F64[0][i];
+	    param_min->data.F32[i] = covar->data.F64[1][i];
+	    param_max->data.F32[i] = covar->data.F64[2][i];
+	}
+	psImageRecycle (covar, params->n, params->n, PS_TYPE_F64);
+    }
+	
+
+    // why is this needed here??? the initial guess on params is provided by the user
+    Params = psVectorCopy (Params, params, PS_TYPE_F32);
+
+    // the user provides the error or NULL.  we need to convert
+    // to appropriate weights
+    dy = psVectorAlloc (y->n, PS_TYPE_F32);
+    if (yErr != NULL) {
+        for (int i = 0; i < dy->n; i++) {
+	    if (yErr->data.F32[i] == 0.0) {
+		dy->data.F32[i] = 1.0;
+		// mask this?  bad pixel, obviously...
+	    } else {
+		dy->data.F32[i] = 1.0 / PS_SQR (yErr->data.F32[i]);
+	    }
+        }
+    } else {
+        for (int i = 0; i < dy->n; i++) {
+            dy->data.F32[i] = 1.0;
+        }
+    }
+
+    // calculate initial alpha and beta, set chisq (min->value)
+    min->value = p_psMinLM_SetABX (alpha, beta, params, paramMask, x, y, dy, func);
+    if (isnan(min->value)) {
+	min->iter = min->maxIter;
+	return (false);
+    }
+    # ifndef PS_NO_TRACE
+    // dump some useful info if trace is defined
+    if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") >= 5) {
+        p_psImagePrint  (psTraceGetDestination(), alpha, "alpha guess");
+        p_psVectorPrint (psTraceGetDestination(), beta, "beta guess");
+        p_psVectorPrint (psTraceGetDestination(), params, "params guess");
+    }
+    if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") == 4) {
+	p_psVectorPrintRow (psTraceGetDestination(), Params, "params guess");
+    }
+    # endif /* PS_NO_TRACE */
+
+    // iterate until the tolerance is reached, or give up
+    while ((min->lastDelta > min->tol) && (min->iter < min->maxIter)) {
+
+        // set a new guess for Alpha, Beta, Params
+        p_psMinLM_GuessABP_EAM (Alpha, Beta, Params, alpha, beta, params, paramMask, beta_lim, param_min, param_max, lambda);
+
+        // measure linear model prediction
+        psF64 dLinear = p_psMinLM_dLinear (Beta, beta, lambda);
+
+        # ifndef PS_NO_TRACE
+        // dump some useful info if trace is defined
+        if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") >= 5) {
+            p_psImagePrint  (psTraceGetDestination(), Alpha, "alpha guess");
+            p_psVectorPrint (psTraceGetDestination(), Beta, "beta guess");
+            p_psVectorPrint (psTraceGetDestination(), Params, "params guess");
+        }
+        if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") == 4) {
+            p_psVectorPrintRow (psTraceGetDestination(), Params, "params guess");
+        }
+        # endif /* PS_NO_TRACE */
+
+        // calculate Chisq for new guess, update Alpha & Beta
+        Chisq = p_psMinLM_SetABX (Alpha, Beta, Params, paramMask, x, y, dy, func);
+
+	// XXX EAM alternate convergence criterion:
+	// compare the delta (min->value - Chisq) with the
+	// expected delta from the linear model (dLinear)
+        // accept new guess (if improvement), or increase lambda
+	psF64 rho = (min->value - Chisq) / dLinear;
+
+        psTrace (".psLib.dataManip.psMinimizeLMChi2", 4, "last chisq: %f, new chisq %f, delta: %f, rho: %f\n", min->value, Chisq, min->lastDelta, rho);
+        # ifndef PS_NO_TRACE
+        // dump some useful info if trace is defined
+        if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") >= 5) {
+            p_psImagePrint  (psTraceGetDestination(), Alpha, "alpha guess");
+            p_psVectorPrint (psTraceGetDestination(), Beta, "beta guess");
+            p_psVectorPrint (psTraceGetDestination(), Params, "params guess");
+        }
+        # endif /* PS_NO_TRACE */
+
+        /* if (Chisq < min->value) {  */
+	if (rho > 0.0) {  
+            min->lastDelta = (min->value - Chisq) / (dy->n - params->n);
+            min->value = Chisq;
+            alpha  = psImageCopy (alpha, Alpha, PS_TYPE_F64);
+            beta   = psVectorCopy (beta, Beta, PS_TYPE_F64);
+            params = psVectorCopy (params, Params, PS_TYPE_F32);
+	    lambda *= 0.1;
+        } else {
+            lambda *= 10.0;
+        }
+        min->iter ++;
+    }
+    psTrace (".psLib.dataManip.psMinimizeLMChi2", 3, "chisq: %f, last delta: %f, Niter: %d\n", min->value, min->lastDelta, min->iter);
+
+    // construct & return the covariance matrix (if requested)
+    if (covar != NULL) {
+      p_psMinLM_GuessABP_EAM (covar, Beta, Params, alpha, beta, params, paramMask, beta_lim, param_min, param_max, 0.0);
+    }      
+
+    // free the internal temporary data
+    psFree (alpha);
+    psFree (Alpha);
+    psFree (beta);
+    psFree (Beta);
+    psFree (Params);
+    psFree (dy);
+
+    if (min->iter == min->maxIter) {
+      return (false);
+    } 
+    return (true);
+}
+
+// XXX EAM : can we use static copies of LUv, LUm, A?
+psBool p_psMinLM_GuessABP_EAM (psImage  *Alpha,
+			       psVector *Beta,
+			       psVector *Params,
+			       const psImage  *alpha,
+			       const psVector *beta,
+			       const psVector *params,
+			       const psVector *paramMask,
+			       const psVector *beta_lim,
+			       const psVector *params_min,
+			       const psVector *params_max,
+			       psF64 lambda)
+{
+
+    # define USE_LU_DECOMP 1
+    # if (USE_LU_DECOMP)
+    psVector *LUv = NULL;
+    psImage  *LUm = NULL;
+    psImage  *A   = NULL;
+    psF32    det;
+
+    // LU decomposition version
+    psTrace (".psLib.dataManip.psMinLM_GuessABP", 5, "using LUD version\n");
+
+    // set new guess values (creates matrix A)
+    A = psImageCopy (NULL, alpha, PS_TYPE_F64);
+    for (int j = 0; j < params->n; j++) {
+	if ((paramMask != NULL) && (paramMask->data.U8[j])) continue;
+        A->data.F64[j][j] = alpha->data.F64[j][j] * (1.0 + lambda);
+    }
+
+    // solve A*beta = Beta (Alpha = 1/A)
+    // these operations do not modify the input values (creates LUm, LUv)
+    LUm   = psMatrixLUD (NULL, &LUv, A);
+    Beta  = psMatrixLUSolve (Beta, LUm, beta, LUv);
+    Alpha = psMatrixInvert (Alpha, A, &det);
+
+    # else
+    // gauss-jordan version
+    psTrace (".psLib.dataManip.psMinLM_GuessABP", 5, "using Gauss-J version");
+
+    // set new guess values (creates matrix A)
+    Beta = psVectorCopy (Beta, beta, PS_TYPE_F64);
+    Alpha = psImageCopy (Alpha, alpha, PS_TYPE_F64);
+    for (int j = 0; j < params->n; j++) {
+	if ((paramMask != NULL) && (paramMask->data.U8[j])) continue;
+        Alpha->data.F64[j][j] = alpha->data.F64[j][j] * (1.0 + lambda);
+    }
+
+    psGaussJordan (Alpha, Beta);
+    # endif
+
+    // apply Beta to get new Params values
+    for (int j = 0; j < params->n; j++) {
+	if ((paramMask != NULL) && (paramMask->data.U8[j])) continue;
+        // Params->data.F32[j] = params->data.F32[j] - Beta->data.F64[j];
+	// continue;
+	// compare Beta to beta limits
+	if (beta_lim != NULL) {
+	    if (fabs(Beta->data.F64[j]) > fabs(beta_lim->data.F32[j])) {
+		Beta->data.F64[j] = (Beta->data.F64[j] > 0) ? fabs(beta_lim->data.F32[j]) : -fabs(beta_lim->data.F32[j]);
+	    }
+	}
+        Params->data.F32[j] = params->data.F32[j] - Beta->data.F64[j];
+	// compare new params to param limits
+	if (params_max != NULL) {
+	    Params->data.F32[j] = PS_MIN (Params->data.F32[j], params_max->data.F32[j]);
+	}
+	if (params_min != NULL) {
+	    Params->data.F32[j] = PS_MAX (Params->data.F32[j], params_min->data.F32[j]);
+	}
+    }
+
+    # if (USE_LU_DECOMP)
+    psFree (A);
+    psFree (LUm);
+    psFree (LUv);
+    # endif
+
+    return true;
+}
+
Index: /tags/eam_tags/psphot_dev_01/psphot/src/psModulesUtils.c
===================================================================
--- /tags/eam_tags/psphot_dev_01/psphot/src/psModulesUtils.c	(revision 22334)
+++ /tags/eam_tags/psphot_dev_01/psphot/src/psModulesUtils.c	(revision 22334)
@@ -0,0 +1,320 @@
+# include "psphot.h"
+
+// XXX EAM : these are my alternate implementations of psModule functions
+
+// this sets and clears bit 0x80
+bool pmSourceLocalSky_EAM (psSource *source,
+			   psStatsOptions statsOptions,
+			   psF32 Radius)
+{
+
+    psImage *image = source->pixels;
+    psImage *mask  = source->mask;
+    psPeak *peak  = source->peak;
+    psRegion *srcRegion;
+
+    srcRegion = psRegionSquare (peak->x, peak->y, Radius);
+    srcRegion = psRegionForImage (srcRegion, mask, srcRegion);
+    psImageMaskRegion (mask, srcRegion, OR, 0x80);
+
+    psStats *myStats = psStatsAlloc(statsOptions);
+    myStats = psImageStats(myStats, image, mask, 0xff);
+    psImageMaskRegion (mask, srcRegion, AND, 0x7f);
+
+    psF64 tmpF64;
+    p_psGetStatValue(myStats, &tmpF64);
+
+    psFree (myStats);
+    psFree (srcRegion);
+
+    if (isnan(tmpF64)) return (false);
+    source->moments = pmMomentsAlloc();
+    source->moments->Sky = (psF32) tmpF64;
+    return (true);
+}
+
+# define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 15
+# define PM_SOURCE_FIT_MODEL_TOLERANCE 0.1
+bool pmSourceFitModel_EAM(psSource *source,
+			  psModel *model,
+			  const bool PSF)
+{
+    PS_PTR_CHECK_NULL(source, false);
+    PS_PTR_CHECK_NULL(source->moments, false);
+    PS_PTR_CHECK_NULL(source->peak, false);
+    PS_PTR_CHECK_NULL(source->pixels, false);
+    PS_PTR_CHECK_NULL(source->mask, false);
+    PS_PTR_CHECK_NULL(source->noise, false);
+    psBool fitStatus = true;
+    psBool onPic     = true;
+    psBool rc        = true;
+    psF32  Ro, ymodel;
+
+
+    // XXX EAM : is it necessary for the mask & noise to exist?  the
+    //           tests below could be conditions (!NULL)
+
+    psModelFunc modelFunc = psModelFunc_GetFunction (model->type);
+    psModelLimits modelLimits = psModelLimits_GetFunction (model->type);
+
+    psVector *params = model->params;
+    psVector *dparams = model->dparams;
+    psVector *paramMask = NULL;
+
+    psVector *beta_lim = NULL;
+    psVector *params_min = NULL;
+    psVector *params_max = NULL;
+
+    int nParams = PSF ? params->n - 4 : params->n;
+    psF32 So = params->data.F32[0];
+
+    // find the number of valid pixels
+    // XXX EAM : this loop and the loop below could just be one pass
+    //           using the psArrayAdd and psVectorExtend functions
+    psS32 count = 0;
+    for (psS32 i = 0; i < source->pixels->numRows; i++) {
+        for (psS32 j = 0; j < source->pixels->numCols; j++) {
+            if (source->mask->data.U8[i][j] == 0) {
+                count++;
+            }
+        }
+    }
+    if (count <  nParams + 1) {
+      psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n");
+      return(false);
+    }
+
+    // construct the coordinate and value entries
+    psArray *x = psArrayAlloc(count);
+    psVector *y = psVectorAlloc(count, PS_TYPE_F32);
+    psVector *yErr = psVectorAlloc(count, PS_TYPE_F32);
+    psS32 tmpCnt = 0;
+    for (psS32 i = 0; i < source->pixels->numRows; i++) {
+        for (psS32 j = 0; j < source->pixels->numCols; j++) {
+            if (source->mask->data.U8[i][j] == 0) {
+                psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
+                // XXX: Convert i/j to image space:
+                // XXX EAM: coord order is (x,y) == (col,row)
+                coord->data.F32[0] = (psF32) (j + source->pixels->col0);
+                coord->data.F32[1] = (psF32) (i + source->pixels->row0);
+                x->data[tmpCnt] = (psPtr *) coord;
+                y->data.F32[tmpCnt] = source->pixels->data.F32[i][j];
+
+		ymodel = modelFunc (NULL, model->params, coord);
+		
+		// this test enhances the noise based on deviation from the model flux
+		Ro = 1.0 + fabs (y->data.F32[tmpCnt] - ymodel) / sqrt(PS_SQR(ymodel - So) + PS_SQR(So));
+                yErr->data.F32[tmpCnt] = sqrt (source->noise->data.F32[i][j] * Ro);
+		// XXX EAM : note the wasted effort: we carry dY^2, take sqrt(), then 
+		//           the minimization function calculates sq()
+                tmpCnt++;
+            }
+        }
+    }
+
+    psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS,
+                            PM_SOURCE_FIT_MODEL_TOLERANCE);
+
+    // PSF model only fits first 4 parameters, FLT model fits all
+    if (PSF) {
+      paramMask = psVectorAlloc (params->n, PS_TYPE_U8);
+      for (int i = 0; i < 4; i++) {
+	paramMask->data.U8[i] = 0;
+      }
+      for (int i = 4; i < paramMask->n; i++) {
+	paramMask->data.U8[i] = 1;
+      }
+    }  
+
+    psImage *covar = psImageAlloc (params->n, 3, PS_TYPE_F64);
+    modelLimits (&beta_lim, &params_min, &params_max);
+    for (int i = 0; i < params->n; i++) {
+	covar->data.F64[0][i] = beta_lim->data.F32[i];
+	covar->data.F64[1][i] = params_min->data.F32[i];
+	covar->data.F64[2][i] = params_max->data.F32[i];
+    }
+
+    psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n");
+    fitStatus = psMinimizeLMChi2_EAM(myMin, covar, params, paramMask, x, y, yErr, modelFunc);
+    for (int i = 0; i < dparams->n; i++) {
+	if ((paramMask != NULL) && paramMask->data.U8[i]) continue;
+	dparams->data.F32[i] = sqrt(covar->data.F64[i][i]);
+    }
+ 
+    // XXX EAM: we need to do something (give an error?) if rc is false
+    // XXX EAM: psMinimizeLMChi2 does not check convergence
+
+    // XXX models can go insane: reject these
+    onPic &= (params->data.F32[2] >= source->pixels->col0);
+    onPic &= (params->data.F32[2] <  source->pixels->col0 + source->pixels->numCols);
+    onPic &= (params->data.F32[3] >= source->pixels->row0);
+    onPic &= (params->data.F32[3] <  source->pixels->row0 + source->pixels->numRows);
+
+    // XXX EAM: save the resulting chisq, nDOF, nIter
+    model->chisq = myMin->value;
+    model->nIter = myMin->iter;
+    model->nDOF  = y->n - nParams;
+
+    // XXX EAM get the Gauss-Newton distance for fixed model parameters
+    if (paramMask != NULL) {
+	psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);
+	psMinimizeGaussNewtonDelta (delta, params, NULL, x, y, yErr, modelFunc);
+	for (int i = 0; i < dparams->n; i++) {
+	    if (!paramMask->data.U8[i]) continue;
+	    dparams->data.F32[i] = delta->data.F64[i];
+	}
+    }
+
+    psFree(paramMask);
+    psFree(x);
+    psFree(y);
+    psFree(myMin);
+
+    rc = (onPic && fitStatus);
+    return(rc);
+}
+
+/******************************************************************************
+pmSourceMoments(source, radius): this function takes a subImage defined in the
+psSource data structure, along with the peak location, and determines the
+various moments associated with that peak.
+ 
+Requires the following to have been created:
+    psSource
+    psSource->peak
+    psSource->pixels
+ 
+XXX: The peak calculations are done in image coords, not subImage coords.
+ 
+XXX: mask values?
+*****************************************************************************/
+# define VALID_RADIUS(X,Y) (((R2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0)
+
+bool pmSourceMoments_EAM(psSource *source,
+			 psF32 radius)
+{
+    PS_PTR_CHECK_NULL(source, NULL);
+    PS_PTR_CHECK_NULL(source->peak, NULL);
+    PS_PTR_CHECK_NULL(source->pixels, NULL);
+    PS_PTR_CHECK_NULL(source->noise, NULL);
+    PS_FLOAT_COMPARE(0.0, radius, NULL);
+
+    //
+    // XXX: Verify the setting for sky if source->moments == NULL.
+    //
+    psF32 sky = 0.0;
+    if (source->moments == NULL) {
+        source->moments = pmMomentsAlloc();
+    } else {
+        sky = source->moments->Sky;
+    }
+
+    //
+    // Sum = SUM (z - sky)
+    // X1  = SUM (x - xc)*(z - sky)
+    // X2  = SUM (x - xc)^2 * (z - sky)
+    // XY  = SUM (x - xc)*(y - yc)*(z - sky)
+    //
+    psF32 peakPixel = -PS_MAX_F32;
+    psS32 numPixels = 0;
+    psF32 Sum = 0.0;
+    psF32 X1 = 0.0;
+    psF32 Y1 = 0.0;
+    psF32 X2 = 0.0;
+    psF32 Y2 = 0.0;
+    psF32 XY = 0.0;
+    psF32 x  = 0;
+    psF32 y  = 0;
+    psF32 R2 = PS_SQR(radius);
+
+    psF32 xPeak = source->peak->x;
+    psF32 yPeak = source->peak->y;
+
+    // XXX why do I get different results for these two methods of finding Sx?
+    // XXX Sx, Sy would be better measured if we clip pixels close to sky
+    // XXX Sx, Sy can still be imaginary, so we probably need to keep Sx^2?
+    // We loop through all pixels in this subimage (source->pixels), and for each
+    // pixel that is not masked, AND within the radius of the peak pixel, we
+    // proceed with the moments calculation.  need to do two loops for a
+    // numerically stable result.  first loop: get the sums.
+    // XXX EAM : mask == 0 is valid
+
+    for (psS32 row = 0; row < source->pixels->numRows ; row++) {
+        for (psS32 col = 0; col < source->pixels->numCols ; col++) {
+            if ((source->mask != NULL) && (source->mask->data.U8[row][col])) continue;
+
+	    psF32 xDiff = col + source->pixels->col0 - xPeak;
+	    psF32 yDiff = row + source->pixels->row0 - yPeak;
+
+	    if (!VALID_RADIUS(xDiff, yDiff)) continue;
+
+	    psF32 pDiff = source->pixels->data.F32[row][col] - sky;
+
+	    // XXX EAM : check for valid S/N in pixel
+	    if (pDiff / sqrt(source->noise->data.F32[row][col]) < 1) continue;
+	    
+	    Sum += pDiff;
+	    X1  += xDiff * pDiff;
+	    Y1  += yDiff * pDiff;
+	    XY  += xDiff * yDiff * pDiff;
+	    
+	    X2  += PS_SQR(xDiff) * pDiff;
+	    Y2  += PS_SQR(yDiff) * pDiff;
+	    
+	    peakPixel = PS_MAX (source->pixels->data.F32[row][col], peakPixel);
+	    numPixels++;
+        }
+    }
+    // XXX EAM - the limit is a bit arbitrary.  make it user defined?
+    if ((numPixels < 3) || (Sum <= 0)) {
+      psTrace (".psModules.pmSourceMoments", 5, "no valid pixels for source\n");
+      return (false);
+    }
+
+    psTrace (".psModules.pmSourceMoments", 5, 
+	     "sky: %f  Sum: %f  X1: %f  Y1: %f  X2: %f  Y2: %f  XY: %f  Npix: %d\n", 
+	     sky, Sum, X1, Y1, X2, Y2, XY, numPixels);
+
+    //
+    // first moment X  = X1/Sum + xc
+    // second moment X = sqrt (X2/Sum - (X1/Sum)^2)
+    // Sxy             = XY / Sum
+    //
+    x = X1/Sum;
+    y = Y1/Sum;
+    if ((fabs(x) > radius) || (fabs(y) > radius)) {
+      psTrace (".psModules.pmSourceMoments", 5, 
+	       "large centroid swing; invalid peak %d, %d\n", 
+	       source->peak->x, source->peak->y);
+      return (false);
+    }
+
+    source->moments->x = x + xPeak;
+    source->moments->y = y + yPeak;
+
+    source->moments->Sxy = XY/Sum - x*y;
+    source->moments->Sum = Sum;
+    source->moments->Peak = peakPixel;
+    source->moments->nPixels = numPixels;
+
+    // XXX EAM : these values can be negative, so we need to limit the range
+    source->moments->Sx = sqrt(PS_MAX(X2/Sum - PS_SQR(x), 0));
+    source->moments->Sy = sqrt(PS_MAX(Y2/Sum - PS_SQR(y), 0));
+
+    psTrace (".psModules.pmSourceMoments", 4, 
+	     "sky: %f  Sum: %f  x: %f  y: %f  Sx: %f  Sy: %f  Sxy: %f\n", 
+	     sky, Sum, source->moments->x, source->moments->y, 
+	     source->moments->Sx, source->moments->Sy, source->moments->Sxy);
+
+    return(true);
+}
+
+bool pmModelFitStatus (psModel *model) {
+
+    bool status;
+
+    psModelFitStatusFunc statusFunc = psModelFitStatusFunc_GetFunction (model->type);
+    status = statusFunc (model);
+
+    return (status);
+}
Index: /tags/eam_tags/psphot_dev_01/psphot/src/psPolynomials.c
===================================================================
--- /tags/eam_tags/psphot_dev_01/psphot/src/psPolynomials.c	(revision 22334)
+++ /tags/eam_tags/psphot_dev_01/psphot/src/psPolynomials.c	(revision 22334)
@@ -0,0 +1,580 @@
+# include "psphot.h"
+
+// write out the terms of the given 1D polynomial
+void psPolynomial1DDump (psPolynomial1D *poly) {
+
+    for (int i = 0; i < poly->n + 1; i++) {
+	fprintf (stderr, "x^%d : %g +/- %g\n", i, poly->coeff[i], poly->coeffErr[i]);
+    }
+}    
+
+psF32 Polynomial1DEval_EAM(psF32 x, const psPolynomial1D* myPoly)
+{
+    psS32 loop_x = 0;
+    psF32 polySum = 0.0;
+    psF32 xSum = 1.0;
+
+    for (loop_x = 0; loop_x < myPoly->n + 1; loop_x++) {
+        if (myPoly->mask[loop_x] == 0) {
+            polySum += xSum * myPoly->coeff[loop_x];
+        }
+        xSum *= x;
+    }
+
+    return(polySum);
+}
+
+psVector *Polynomial1DEvalVector_EAM(const psPolynomial1D *myPoly,
+				       const psVector *x)
+{
+    PS_POLY_CHECK_NULL(myPoly, NULL);
+    PS_VECTOR_CHECK_NULL(x, NULL);
+    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
+
+    psVector *tmp;
+
+    tmp = psVectorAlloc(x->n, PS_TYPE_F64);
+    for (psS32 i=0;i<x->n;i++) {
+        tmp->data.F64[i] = Polynomial1DEval_EAM(x->data.F64[i], myPoly);
+    }
+
+    return(tmp);
+}
+
+static void psPolynomial1DFree (psPolynomial1D *poly) {
+
+  if (poly == NULL) return;
+  psFree (poly->coeff);
+  psFree (poly->coeffErr);
+  psFree (poly->mask);
+  return;
+}
+
+// XXX EAM : use Nterm = Norder + 1 definition  
+// XXX EAM : should we provide both order and nterms in struct?
+psPolynomial1D* Polynomial1DAlloc(psS32 nOrder,
+				  psPolynomialType type)
+{
+    PS_INT_CHECK_NON_NEGATIVE(nOrder, NULL);
+
+    psS32 i = 0;
+    psS32 nTerm = nOrder + 1;
+    psPolynomial1D* newPoly = NULL;
+
+    newPoly = (psPolynomial1D* ) psAlloc(sizeof(psPolynomial1D));
+    newPoly->type = type;
+    newPoly->n = nOrder;
+    newPoly->coeff = (psF32 *)psAlloc(nTerm * sizeof(psF32));
+    newPoly->coeffErr = (psF32 *)psAlloc(nTerm * sizeof(psF32));
+    newPoly->mask = (psU8 *)psAlloc(nTerm * sizeof(psU8));
+    for (i = 0; i < nTerm; i++) {
+        newPoly->coeff[i] = 0.0;
+        newPoly->coeffErr[i] = 0.0;
+        newPoly->mask[i] = 0;
+    }
+    p_psMemSetDeallocator(newPoly, (psFreeFcn) psPolynomial1DFree);
+    return(newPoly);
+}
+
+// XXX EAM : my alternate BuildSums1D
+static psVector *BuildSums1D(psVector* sums, 
+			     psF64 x,
+			     psS32 nTerm)
+{
+    psS32 nSum = 0;
+    psF64 xSum = 0.0;
+
+    nSum = 2*nTerm;
+    if (sums == NULL) {
+        sums = psVectorAlloc(nSum, PS_TYPE_F64);
+    }
+    if (nSum > sums->n) {
+        sums = psVectorRealloc(sums, nSum);
+    }
+
+    xSum = 1.0;
+    for (int i = 0; i < nSum; i++) {
+        sums->data.F64[i] = xSum;
+        xSum *= x;
+    }
+    return (sums);
+}
+
+// XXX EAM : test version of 1d fitting
+psPolynomial1D* VectorFitPolynomial1DOrd_EAM(psPolynomial1D* myPoly,
+					     psVector *mask,
+					     const psVector *x,
+					     const psVector *y,
+					     const psVector *yErr)
+{
+    // I think this is 1 dimension down
+    psImage*  	  A = NULL;
+    psVector* 	  B = NULL;
+    psVector* xSums = NULL;
+    psS32 nTerm;
+    psF64 wt;
+
+    psTrace(".psLib.dataManip.VectorFitPolynomial1DOrd", 4,
+            "---- VectorFitPolynomial1DOrd() begin ----\n");
+
+    // dump minutiae
+# ifndef PS_NO_TRACE
+    if (psTraceGetLevel (".psLib.dataManip.VectorFitPolynomial1DOrd") >= 5) {
+      FILE *f = psTraceGetDestination ();
+      fprintf (f, "VectorFitPolynomial1D()\n");
+      for (int i = 0; i < x->n; i++) {
+	fprintf (f, "(x, y, yErr) is (%f, %f, %f)\n", x->data.F64[i], y->data.F64[i], yErr->data.F64[i]);
+      }
+    }
+# endif
+      
+    nTerm = myPoly->n + 1;
+    A     = psImageAlloc(nTerm, nTerm, PS_TYPE_F64);
+    B     = psVectorAlloc(nTerm, PS_TYPE_F64);
+
+    // Initialize data structures (why is this not a function!)
+    for (int i = 0; i < nTerm; i++) {
+        B->data.F64[i] = 0.0;
+        for (int j = 0; j < nTerm; j++) {
+            A->data.F64[i][j] = 0.0;
+        }
+    }
+
+    // xSums look like: 1, x, x^2, ... x^(2n+1)
+
+    // Build the B and A data structs.
+    for (int k = 0; k < x->n; k++) {
+        if ((mask != NULL) && mask->data.U8[k]) continue;
+	xSums = BuildSums1D(xSums, x->data.F64[k], nTerm);
+      
+	if (yErr == NULL) {
+	    wt = 1.0;
+	} else {
+	    // this should probably by yErr^2 !!
+	    // this should filter yErr == 0 values
+	    wt = 1.0 / PS_SQR(yErr->data.F64[k]);
+	}
+	for (int i = 0; i < nTerm; i++) {
+	    B->data.F64[i] += y->data.F64[k] * xSums->data.F64[i] * wt;
+	}
+
+	// we could skip half of the array and assign at the end
+	// we must handle masked orders
+	for (int i = 0; i < nTerm; i++) {
+	    for (int j = 0; j < nTerm; j++) {
+		A->data.F64[i][j] += xSums->data.F64[i + j] * wt;
+	    }
+	}
+    }
+ 
+    // GaussJordan version
+    if (0) {
+	// does the solution in place
+	psGaussJordan (A, B);
+    
+	// the first nTerm entries in B correspond directly to the desired 
+	// polynomial coefficients.  this is only true for the 1D case
+	for (int k = 0; k < nTerm; k++) {
+	    myPoly->coeff[k] = B->data.F64[k];
+	}
+    } 
+    else 
+    // LUD version of the fit
+    {
+	psImage *ALUD = NULL;
+	psVector* outPerm = NULL;
+	psVector* coeffs = NULL;
+    
+	ALUD = psImageAlloc(nTerm, nTerm, PS_TYPE_F64);
+	ALUD = psMatrixLUD(ALUD, &outPerm, A);
+	coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm);
+	for (int k = 0; k < nTerm; k++) {
+	    myPoly->coeff[k] = coeffs->data.F64[k];
+	}
+	psFree (ALUD);
+	psFree (outPerm);
+	psFree (coeffs);
+    }
+
+    psFree(A);
+    psFree(B);
+    psFree(xSums);
+
+    psTrace(".psLib.dataManip.VectorFitPolynomial1DOrd", 4,
+            "---- VectorFitPolynomial1DOrd() begin ----\n");
+    return (myPoly);
+}
+
+// ********************** 2D polynomial functions ******************
+
+// XXX EAM : this version uses myPoly->nX as Norder, not Nterms
+psVector *Polynomial2DEvalVector(const psPolynomial2D *myPoly,
+				 const psVector *x,
+				 const psVector *y)
+
+{
+    PS_POLY_CHECK_NULL(myPoly, NULL);
+    PS_VECTOR_CHECK_NULL(x, NULL);
+    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NULL);
+    PS_VECTOR_CHECK_NULL(y, NULL);
+    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F32, NULL);
+
+    psVector *tmp;
+    psS32 vecLen=x->n;
+
+    // Determine the length of the output vector to by the minimum of the x,y vectors
+    if (y->n < vecLen) {
+        vecLen = y->n;
+    }
+
+    // Create output vector to return
+    tmp = psVectorAlloc(vecLen, PS_TYPE_F32);
+
+    // Evaluate the polynomial at the specified points
+    for (psS32 i=0; i<vecLen; i++) {
+        tmp->data.F32[i] = Polynomial2DEval(myPoly,x->data.F32[i],y->data.F32[i]);
+    }
+
+    // Return output vector
+    return(tmp);
+}
+
+// XXX EAM : this version uses the F64 vectors
+psVector *Polynomial2DEvalVectorD(const psPolynomial2D *myPoly,
+				  const psVector *x,
+				  const psVector *y)
+
+{
+    PS_POLY_CHECK_NULL(myPoly, NULL);
+    PS_VECTOR_CHECK_NULL(x, NULL);
+    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
+    PS_VECTOR_CHECK_NULL(y, NULL);
+    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F64, NULL);
+
+    psVector *tmp;
+    psS32 vecLen=x->n;
+
+    // Determine the length of the output vector to by the minimum of the x,y vectors
+    if (y->n < vecLen) {
+        vecLen = y->n;
+    }
+
+    // Create output vector to return
+    tmp = psVectorAlloc(vecLen, PS_TYPE_F64);
+
+    // Evaluate the polynomial at the specified points
+    for (psS32 i=0; i<vecLen; i++) {
+        tmp->data.F64[i] = Polynomial2DEval(myPoly,x->data.F64[i],y->data.F64[i]);
+    }
+
+    // Return output vector
+    return(tmp);
+}
+
+static void psPolynomial2DFree(psPolynomial2D *poly) {
+
+  if (poly == NULL) return;
+  for (int i = 0; i < poly->nX + 1; i++) {
+    psFree (poly->coeff[i]);
+    psFree (poly->coeffErr[i]);
+    psFree (poly->mask[i]);
+  }
+  psFree (poly->coeff);
+  psFree (poly->coeffErr);
+  psFree (poly->mask);
+  return;
+}
+
+// XXX EAM : use Nterm = Norder + 1 definition  
+// the user requests a polynomial of order Norder
+psPolynomial2D* Polynomial2DAlloc(psS32 nXorder, psS32 nYorder,
+                                    psPolynomialType type)
+{
+    PS_INT_CHECK_NON_NEGATIVE(nXorder, NULL);
+    PS_INT_CHECK_NON_NEGATIVE(nYorder, NULL);
+
+    psS32 x = 0;
+    psS32 y = 0;
+    psS32 nXterm = nXorder + 1;
+    psS32 nYterm = nYorder + 1;
+    psPolynomial2D* newPoly = NULL;
+
+    newPoly = (psPolynomial2D* ) psAlloc(sizeof(psPolynomial2D));
+    newPoly->type = type;
+    newPoly->nX = nXorder;
+    newPoly->nY = nYorder;
+
+    newPoly->coeff = (psF32 **)psAlloc(nXterm * sizeof(psF32 *));
+    newPoly->coeffErr = (psF32 **)psAlloc(nXterm * sizeof(psF32 *));
+    newPoly->mask = (psU8 **)psAlloc(nXterm * sizeof(psU8 *));
+    for (x = 0; x < nXterm; x++) {
+        newPoly->coeff[x] = (psF32 *)psAlloc(nYterm * sizeof(psF32));
+        newPoly->coeffErr[x] = (psF32 *)psAlloc(nYterm * sizeof(psF32));
+        newPoly->mask[x] = (psU8 *)psAlloc(nYterm * sizeof(psU8));
+    }
+    for (x = 0; x < nXterm; x++) {
+        for (y = 0; y < nYterm; y++) {
+            newPoly->coeff[x][y] = 0.0;
+            newPoly->coeffErr[x][y] = 0.0;
+            newPoly->mask[x][y] = 0;
+        }
+    }
+    p_psMemSetDeallocator(newPoly, (psFreeFcn) psPolynomial2DFree);
+    return(newPoly);
+}
+
+// XXX EAM : BuildSums2D in analogy with BuildSums1D
+static psImage *BuildSums2D(psImage* sums,
+			    psF64 x,      psF64 y,
+			    psS32 nXterm, psS32 nYterm)
+{
+    psS32 nXsum = 0;
+    psS32 nYsum = 0;
+    psF64 xSum = 1.0;
+    psF64 ySum = 1.0;
+
+    nXsum = 2*nXterm;
+    nYsum = 2*nYterm;
+    if (sums == NULL) {
+        sums = psImageAlloc(nXsum, nYsum, PS_TYPE_F64);
+    }
+    if ((nXsum != sums->numCols) || (nYsum != sums->numRows)) {
+	psFree (sums);
+        sums = psImageAlloc(nXsum, nYsum, PS_TYPE_F64);
+    }
+
+    ySum = 1.0;
+    for (int j = 0; j < nYsum; j++) {
+	xSum = ySum;
+	for (int i = 0; i < nXsum; i++) {
+	    sums->data.F64[i][j] = xSum;
+	    xSum *= x;
+	}
+	ySum *= y;
+    }
+    return (sums);
+}
+
+// XXX EAM : test version of 2d fitting
+psPolynomial2D* VectorFitPolynomial2DOrd_EAM(psPolynomial2D* myPoly,
+					     psVector* mask,
+					     const psVector* x,
+					     const psVector* y,
+					     const psVector* z,
+					     const psVector* zErr)
+{
+    // I think this is 1 dimension down
+    psImage*  	  A = NULL;
+    psVector* 	  B = NULL;
+    psImage*   Sums = NULL;
+    psF64 wt;
+    psS32 nTerm, nXterm, nYterm;
+
+    nXterm = myPoly->nX + 1;
+    nYterm = myPoly->nY + 1;
+    nTerm = nXterm * nYterm;
+
+    A = psImageAlloc(nTerm, nTerm, PS_TYPE_F64);
+    B = psVectorAlloc(nTerm, PS_TYPE_F64);
+
+    // Initialize data structures (why is this not a function!)
+    for (int i = 0; i < nTerm; i++) {
+        B->data.F64[i] = 0.0;
+        for (int j = 0; j < nTerm; j++) {
+            A->data.F64[i][j] = 0.0;
+        }
+    }
+
+    // Sums look like: 1, x, x^2, ... x^(2n+1), y, xy, x^2y, ... x^(2n+1)
+
+    // Build the B and A data structs.
+    for (int k  = 0; k < x->n; k++) {
+      if ((mask != NULL) && mask->data.U8[k]) continue;
+	Sums = BuildSums2D(Sums, x->data.F64[k], y->data.F64[k], nXterm, nYterm);
+      
+	if (zErr == NULL) {
+	    wt = 1.0;
+	} else {
+	    // this should probably by zErr^2 !!
+	    // this should filter zErr == 0 values
+	    wt = 1.0 / zErr->data.F64[k];
+	}
+
+	// we could skip half of the array and assign at the end
+	// we must handle masked orders
+	for (int n = 0; n < nXterm; n++) {
+	    for (int m = 0; m < nYterm; m++) {
+		B->data.F64[n+m*nXterm] += z->data.F64[k] * Sums->data.F64[n][m] * wt;
+	    }
+	}
+
+	for (int i = 0; i < nXterm; i++) {
+	    for (int j = 0; j < nYterm; j++) {
+		for (int n = 0; n < nXterm; n++) {
+		    for (int m = 0; m < nYterm; m++) {
+			A->data.F64[i+j*nXterm][n+m*nXterm] += Sums->data.F64[i+n][j+m] * wt;
+		    }
+		}
+	    }
+	}
+    }
+ 
+    // does the solution in place
+    psGaussJordan (A, B);
+    
+    // XXX: How do we know if these routines were successful?
+    // ALUD = psMatrixLUD(ALUD, &outPerm, A);
+    // coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm);
+
+    for (int n = 0; n < nXterm; n++) {
+	for (int m = 0; m < nYterm; m++) {
+	    myPoly->coeff[n][m] = B->data.F64[n+m*nXterm];
+	}
+    }
+
+    psFree(A);
+    psFree(B);
+    psFree(Sums);
+
+    psTrace(".psLib.dataManip.VectorFitPolynomial2DOrd", 4,
+            "---- VectorFitPolynomial2DOrd() begin ----\n");
+    return (myPoly);
+}
+
+// write out the terms of the given 2D polynomial
+void psPolynomial2DDump (psPolynomial2D *poly) {
+
+    for (int i = 0; i < poly->nX + 1; i++) {
+	for (int j = 0; j < poly->nY + 1; j++) {
+	    fprintf (stderr, "x^%d y^%d : %g +/- %g\n", i, j, poly->coeff[i][j], poly->coeffErr[i][j]);
+	}
+    }
+}    
+
+psPolynomial2D* RobustFit2D_nomask(psPolynomial2D* poly,
+			    const psVector* x,
+			    const psVector* y,
+			    const psVector* z,
+			    const psVector* dz)
+{
+    psVector *X;
+    psVector *Y;
+    psVector *Z;
+    psVector *dZ;
+
+    psVector *zFit   = NULL;
+    psVector *zResid = NULL;
+    psStats  *stats  = NULL;
+
+    X  = psVectorCopy (NULL, x, PS_TYPE_F64);
+    Y  = psVectorCopy (NULL, y, PS_TYPE_F64);
+    Z  = psVectorCopy (NULL, z, PS_TYPE_F64);
+    dZ = psVectorCopy (NULL, dz, PS_TYPE_F64);
+
+    for (int N = 0; N < 3; N++) {
+	// XXX EAM : this would be better defined with an element mask
+	poly   = VectorFitPolynomial2DOrd_EAM (poly, NULL, X, Y, Z, dZ);
+	zFit   = Polynomial2DEvalVectorD (poly, x, y);
+	zResid = (psVector *) psBinaryOp (NULL, (void *) z, "-", (void *) zFit);
+
+	stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
+	stats  = psVectorStats (stats, zResid, NULL, NULL, 0);
+	psTrace (".psphot.RobustFit", 4, "residual stats for robust fit:  %g +/- %g (%d pts)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues);
+
+	// re-create X, Y, Z, dZ if pts are valid
+	int n = 0;
+	for (int i = 0; i < zResid->n; i++) {
+	    if (fabs(zResid->data.F64[i] - stats->clippedMean) > 3*stats->clippedStdev) {
+	      continue;
+	    } 
+	    X->data.F64[n]  =  x->data.F64[i];
+	    Y->data.F64[n]  =  y->data.F64[i];
+	    Z->data.F64[n]  =  z->data.F64[i];
+	    dZ->data.F64[n] = dz->data.F64[i];
+	    n++;
+	}
+	X->n = n;
+	Y->n = n;
+	Z->n = n;
+	dZ->n = n;
+    }
+    return (poly);
+}
+
+// XXX EAM : be careful here with F32 vs F64 vectors 
+psPolynomial2D* RobustFit2D(psPolynomial2D* poly,
+			    psVector* mask,
+			    const psVector* x,
+			    const psVector* y,
+			    const psVector* z,
+			    const psVector* dz)
+{
+    PS_VECTOR_CHECK_NULL(mask, NULL);
+    PS_VECTOR_CHECK_NULL(x, NULL);
+    PS_VECTOR_CHECK_NULL(y, NULL);
+    PS_VECTOR_CHECK_NULL(z, NULL);
+    PS_VECTOR_CHECK_NULL(dz, NULL);
+
+    psVector *zFit   = NULL;
+    psVector *zResid = psVectorAlloc (x->n, PS_TYPE_F64);
+    psStats  *stats  = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
+
+    for (int N = 0; N < 3; N++) {
+	poly   = VectorFitPolynomial2DOrd_EAM (poly, mask, x, y, z, dz);
+	zFit   = Polynomial2DEvalVectorD (poly, x, y);
+	zResid = (psVector *) psBinaryOp (zResid, (void *) z, "-", (void *) zFit);
+
+	stats  = psVectorStats (stats, zResid, NULL, mask, 1);
+	psTrace (".psphot.RobustFit", 4, "residual stats for robust fit:  %g +/- %g\n", 
+		 stats->sampleMean, stats->sampleStdev);
+
+	// set mask if pts are not valid
+	// we are masking out any point which is out of range 
+	// recovery is not allowed with this scheme
+	for (int i = 0; i < zResid->n; i++) {
+	  if (mask->data.U8[i]) continue;
+	  if (fabs(zResid->data.F64[i] - stats->sampleMean) > 3*stats->sampleStdev) {
+	    mask->data.U8[i] = 1;
+	    continue;
+	  }       
+	}
+	psFree (zFit);
+    }
+    psFree (zResid);
+    psFree (stats);
+    return (poly);
+}
+
+// XXX EAM : VectorFitPolynomial2DOrd and Polynomial2DEvalVector require different types (F32 vs F64)
+
+# if (0) // moved to psLib
+// XXX EAM : this version uses myPoly->nX as Norder, not Nterms
+psF32 Polynomial2DEval(const psPolynomial2D* myPoly,
+		       psF32 x,
+		       psF32 y)
+{
+    PS_POLY_CHECK_NULL(myPoly, NAN);
+
+    psS32 loop_x = 0;
+    psS32 loop_y = 0;
+    psF32 polySum = 0.0;
+    psF32 xSum = 1.0;
+    psF32 ySum = 1.0;
+
+    // XXX EAM : nX is order, not nTerms
+    for (loop_x = 0; loop_x < myPoly->nX + 1; loop_x++) {
+        ySum = xSum;
+	// XXX EAM : nX is order, not nTerms
+        for (loop_y = 0; loop_y < myPoly->nY + 1; loop_y++) {
+            if (myPoly->mask[loop_x][loop_y] == 0) {
+                polySum += ySum * myPoly->coeff[loop_x][loop_y];
+            }
+            ySum *= y;
+        }
+        xSum *= x;
+    }
+
+    return(polySum);
+}
+# endif
+
Index: /tags/eam_tags/psphot_dev_01/psphot/src/psphot.c
===================================================================
--- /tags/eam_tags/psphot_dev_01/psphot/src/psphot.c	(revision 22334)
+++ /tags/eam_tags/psphot_dev_01/psphot/src/psphot.c	(revision 22334)
@@ -0,0 +1,51 @@
+# include "psphot.h"
+
+int main (int argc, char **argv) {
+
+    psMetadata  *config  = NULL;
+    psImageData *imdata  = NULL;
+    psArray     *sources = NULL;
+    psArray     *peaks   = NULL;
+    pmPSF       *psf     = NULL;
+    psStats     *sky     = NULL;
+    pmPSFClump   psfClump;
+
+    config = psphotArguments (&argc, argv);
+
+    // load input data (image and config)
+    // create or load mask and noise images
+    // we have memory leaks here -- may be from psMetadata
+    imdata = psphotSetup (config);
+
+    // measure image stats for initial guess 
+    sky = psphotImageStats (imdata, config);
+
+    // find the peaks in the image
+    peaks = pmPeaksSigmaLimit (imdata, config, sky);
+
+    // construct sources and measure basic stats
+    sources = psphotSourceStats (imdata, config, peaks);
+
+    // classify sources based on moments, brightness
+    // XXX - this has 1 leak
+    psfClump = pmSourcePSFClump (sources, config);
+
+    // group into STAR, COSMIC, GALAXY, SATURATED
+    pmSourceRoughClass (sources, config, psfClump);
+
+    // source analysis is done in S/N order (brightest first)
+    sources = psArraySort (sources, psphotSortBySN);
+
+    // use bright stellar objects to measure PSF
+    psf = psphotChoosePSF (config, sources, sky);
+
+    // test PSF on all except SATURATE and DEFECT (artifacts)
+    psphotApplyPSF (imdata, config, sources, psf, sky);
+
+    // fit extended objects with galaxy models
+    psphotFitGalaxies (imdata, config, sources, sky);
+
+    // write out data in appropriate format
+    psphotOutput (imdata, config, sources);
+    exit (0);
+}
Index: /tags/eam_tags/psphot_dev_01/psphot/src/psphot.h
===================================================================
--- /tags/eam_tags/psphot_dev_01/psphot/src/psphot.h	(revision 22334)
+++ /tags/eam_tags/psphot_dev_01/psphot/src/psphot.h	(revision 22334)
@@ -0,0 +1,158 @@
+# include <pslib.h>
+# include <pmObjects.h>
+# include <string.h>
+# include <strings.h>
+# include <unistd.h>
+# include <stdlib.h>
+# include <math.h>
+
+# define  M_PI 3.14159265358979323846264338328      /* pi */
+
+typedef struct {
+  psImage *image;
+  psImage *mask;
+  psImage *noise;
+  psMetadata *header;
+} psImageData;
+
+// data to test a given PSF model type
+typedef struct {
+  psModelType modelType;
+  pmPSF      *psf;
+  psArray    *sources;      // pointers to the original sources
+  psArray    *modelFLT;     // model fits, floating parameters 
+  psArray    *modelPSF;     // model fits, PSF parameters
+  psVector   *mask;
+  psVector   *metric;
+  psVector   *fitMag;
+  float       ApResid;
+  float       dApResid;
+  float       skyBias;
+} pmPSF_Test;
+
+// structure to carry a dynamic string
+typedef struct {
+  int NLINE;
+  int Nline;
+  char *line;
+} psLine;
+
+enum {
+    PSPHOT_MASK_INVALID   = 0x01,
+    PSPHOT_MASK_SATURATED = 0x02,
+    PSPHOT_MASK_KEEP      = 0x08,
+} MaskValues;
+
+// used by mask operations
+# define AND true
+# define OR false
+
+# define psMemCopy(A)(psMemIncrRefCounter((A)))
+
+// top-level psphot functions
+psMetadata  *psphotArguments (int *argc, char **argv);
+psImageData *psphotSetup (psMetadata *config);
+psStats     *psphotImageStats (psImageData *imdata, psMetadata *config);
+psArray     *pmPeaksSigmaLimit (psImageData *imdata, psMetadata *config, psStats *sky);
+psArray     *psphotSourceStats (psImageData *imdata, psMetadata *config, psArray *allpeaks);
+pmPSF       *psphotChoosePSF (psMetadata *config, psArray *sources, psStats *sky);
+bool         psphotApplyPSF (psImageData *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
+bool         psphotFitGalaxies (psImageData *imdata, psMetadata *config, psArray *sources, psStats *skyStats);
+void         psphotOutput (psImageData *imdata, psMetadata *config, psArray *sources);
+
+bool         psphotMarkPSF (psSource *source, float shapeNsigma, float minSN, float maxChi, float SATURATE);
+bool         psphotSubtractPSF (psSource *source);
+int 	     psphotSortBySN (const void **a, const void **b);
+
+// psImageData functions
+psImageData *psImageDataAlloc (psImage *image, psImage *noise, psImage *mask, psMetadata *header);
+
+// psf utilities
+pmPSF       *pmPSFAlloc (psModelType type);
+pmPSF_Test  *pmPSF_TestAlloc (psArray *stars, char *modelName);
+pmPSF_Test  *pmPSF_TestModel (psArray *sources, char *modelName, float radius);
+bool	     pmPSFFromModels (pmPSF *psf, psArray *models, psVector *mask);
+psModel	    *psModelFromPSF (psModel *model, pmPSF *psf);
+bool	     pmSourcePhotometry (float *fitMag, float *obsMag, psModel *model, psImage *image, psImage *mask);
+bool	     pmPSFMetricModel (pmPSF_Test *test, float RADIUS);
+
+// output functions
+bool 	     pmSourcesWriteText (psImageData *imdata, char *filename, psArray *sources);
+bool 	     pmSourcesWriteOBJ (psImageData *imdata, char *filename, psArray *sources);
+bool 	     pmSourcesWriteCMP (psImageData *imdata, char *filename, psArray *sources);
+bool 	     pmSourcesWriteCMF (psImageData *imdata, char *filename, psArray *sources);
+bool 	     pmSourcesWriteSX (psImageData *imdata, char *filename, psArray *sources);
+int  	     pmSourcesDophotType (psSource *source);
+bool 	     psPeaksWriteText (psArray *sources, char *filename);
+bool 	     psMomentsWriteText (psArray *sources, char *filename);
+bool 	     psModelWritePSFs (psArray *sources, char *filename);
+bool 	     psModelWriteFLTs (psArray *sources, char *filename);
+bool 	     psModelWriteNULLs (psArray *sources, char *filename);
+
+// psModule extra utilities
+bool 	     pmSourceDefinePixels(psSource *mySource, const psImageData *imdata, psF32 x, psF32 y, psF32 Radius);
+bool 	     pmSourceLocalSky_EAM (psSource *source, psStatsOptions statsOptions, psF32 Radius);
+bool 	     pmSourceFitModel_EAM(psSource *source, psModel *model, const bool PSF);
+bool 	     pmSourceMoments_EAM(psSource *source, psF32 radius);
+bool 	     pmModelFitStatus (psModel *model);
+int	     pmSourceDophotType (psSource *source);
+
+// minimize 
+psBool 	     p_psMinLM_GuessABP_EAM (psImage  *Alpha, psVector *Beta, psVector *Params, const psImage  *alpha, const psVector *beta, const psVector *params, const psVector *paramMask, const psVector *beta_lim, const psVector *params_min, const psVector *params_max, psF64 lambda);
+psBool 	     psMinimizeLMChi2_EAM(psMinimization *min, psImage *covar, psVector *params, const psVector *paramMask, const psArray *x, const psVector *y, const psVector *yErr, psMinimizeLMChi2Func func);
+psF64        p_psMinLM_dLinear (const psVector *Beta, const psVector *beta, psF64 lambda);
+
+// psLib extra utilities
+bool 	     psTimerStart (char *name);
+void 	     psTimerFree ();
+bool 	     psTimerClear (char *name);
+psF64 	     psTimerMark (char *name);
+psS32 	     psLogArguments (int *argc, char **argv);
+psS32 	     psTraceArguments (int *argc, char **argv);
+int	     psArgumentGet (int argc, char **argv, char *arg);
+int	     psArgumentRemove (int N, int *argc, char **argv);
+psF32        pmConfigLookupF32 (bool *status, psMetadata *config, psMetadata *header, char *name);
+psVector    *psVectorCreate (double lower, double upper, double delta, psElemType type);
+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);
+psVector    *psGetRowVectorFromImage(psImage *image, psU32 row);
+
+// basic image functions
+bool         psImageInit (psImage *image,...);
+void	     psImageSmooth (psImage *image, float sigma, float Nsigma);
+psRegion    *psRegionForImage (psRegion *out, psImage *image, psRegion *in);
+psRegion    *psRegionSquare (psF32 x, psF32 y, psF32 radius);
+int          psphotSaveImage (psMetadata *header, psImage *image, char *filename);
+
+// psLine functions
+psLine      *psLineAlloc (int Nline);
+bool	     psLineInit (psLine *line);
+bool	     psLineAdd (psLine *line, char *format, ...);
+
+// polynomial functions
+psF32	        Polynomial2DEval(const psPolynomial2D* myPoly, psF32 x, psF32 y);
+psImage	       *psBuildSums2D(psImage* sums,psF64 x,psF64 y,psS32 nXterm, psS32 nYterm);
+psPolynomial2D *VectorFitPolynomial2DOrd_EAM(psPolynomial2D* myPoly, psVector* mask, const psVector* x, const psVector* y, const psVector* z, const psVector* zErr);
+psPolynomial2D *Polynomial2DAlloc(psS32 nXorder, psS32 nYorder, psPolynomialType type);
+void	        psPolynomial2DDump (psPolynomial2D *poly);
+psPolynomial2D *RobustFit2D_nomask(psPolynomial2D* poly, const psVector* x, const psVector* y, const psVector* z, const psVector* zErr);
+psPolynomial2D *RobustFit2D(psPolynomial2D* poly, psVector* mask, const psVector* x, const psVector* y, const psVector* z, const psVector* zErr);
+psVector       *Polynomial2DEvalVector(const psPolynomial2D *myPoly, const psVector *x,const psVector *y);
+psVector       *psBuildSums1D(psVector* sums, psF64 x,psS32 nTerm);
+void	        psPolynomial1DDump (psPolynomial1D *poly);
+psF32		Polynomial1DEval_EAM(psF32 x, const psPolynomial1D* myPoly);
+psVector       *Polynomial1DEvalVector_EAM(const psPolynomial1D *myPoly, const psVector *x);
+psPolynomial1D *Polynomial1DAlloc(psS32 nOrder, psPolynomialType type);
+psPolynomial1D *VectorFitPolynomial1DOrd_EAM(psPolynomial1D* myPoly, psVector* mask, const psVector* x, const psVector* y, const psVector* yErr);
+
+// fitsource
+bool	        onesource (psImageData *imdata, psMetadata *config, char *modelName, float x, float y);
+psMetadata     *fs_args (int *argc, char **argv);
+int             fs_usage ();
+bool		DumpImage (psImage *image, char *filename);
+
+// deprecated
+// void pmSourceMaskRegion (psSource *source, psRegion *region);
+//psVector    *p_psGetRowVectorFromImage(psImage *image, psU32 row);
Index: /tags/eam_tags/psphot_dev_01/psphot/src/psphotApplyPSF.c
===================================================================
--- /tags/eam_tags/psphot_dev_01/psphot/src/psphotApplyPSF.c	(revision 22334)
+++ /tags/eam_tags/psphot_dev_01/psphot/src/psphotApplyPSF.c	(revision 22334)
@@ -0,0 +1,93 @@
+# include "psphot.h"
+
+// fit psf model to all objects 
+// PSFSTAR objects will be refitted
+// run this function to a specific flux limit?
+
+bool psphotApplyPSF (psImageData *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) 
+{ 
+    bool  status;
+    float x;
+    float y;
+    int   Nfit = 0;
+    int   Nsub = 0;
+    int   Niter = 0;
+
+    psTimerStart ("psphot");
+
+    // we may set this differently here from the value used to mark likely saturated stars
+    float SATURATION   = psMetadataLookupF32 (&status, config, "SATURATION");
+    float OUTER_RADIUS = psMetadataLookupF32 (&status, config, "SKY_OUTER_RADIUS");
+
+    float PSF_MIN_SN   	   = psMetadataLookupF32 (&status, config, "PSF_MIN_SN");
+    float PSF_MAX_CHI  	   = psMetadataLookupF32 (&status, config, "PSF_MAX_CHI");
+    float PSF_FIT_NSIGMA   = psMetadataLookupF32 (&status, config, "PSF_FIT_NSIGMA");
+    float PSF_FIT_PADDING  = psMetadataLookupF32 (&status, config, "PSF_FIT_PADDING");
+    float PSF_SHAPE_NSIGMA = psMetadataLookupF32 (&status, config, "PSF_SHAPE_NSIGMA");
+
+    // set the object surface-brightness limit for fitted pixels
+    float FLUX_LIMIT  = PSF_FIT_NSIGMA * sky->sampleStdev;
+    psLogMsg ("psphot.apply_psf_model", 3, "fitting pixels with at least %f object counts\n", FLUX_LIMIT);
+
+    // this function specifies the radius at this the model hits the given flux
+    psModelRadius modelRadius = psModelRadius_GetFunction (psf->type);
+
+    for (int i = 0; i < sources->n; i++) {
+
+	psSource *source = sources->data[i];
+
+	// skip non-astronomical objects (very likely defects)
+	// XXX EAM : should we try these anyway?
+	if (source->type == PS_SOURCE_DEFECT) continue; 
+	if (source->type == PS_SOURCE_SATURATED) continue;
+
+	// use the source moments, etc to guess basic model parameters
+	psModel *modelFLT = pmSourceModelGuess (source, psf->type); 
+
+	// set PSF parameters for this model
+	psModel *model = psModelFromPSF (modelFLT, psf);
+	x = model->params->data.F32[2];
+	y = model->params->data.F32[3];
+	psFree (modelFLT);
+
+	// set the fit radius based on the object flux limit and the model
+	// XXX EAM : FLUX_LIMIT should be set based on local sky model (not global median)
+	model->radius = modelRadius (model->params, FLUX_LIMIT) + PSF_FIT_PADDING;
+	if (isnan(model->radius)) {
+	  psAbort ("apply_psf_model", "error in radius");
+	}
+	
+	// check if we need to redefine the pixels
+	// XXX EAM : a better test would examine the source pixels
+	if (model->radius > OUTER_RADIUS) {
+	  // (re)-allocate image, noise, mask arrays for each peak (square of radius OUTER)
+	  pmSourceDefinePixels (source, imdata, x, y, model->radius);
+	}
+
+	// fit PSF model (set/unset the pixel mask)
+	psImageKeepCircle (source->mask, x, y, model->radius, OR, PSPHOT_MASK_KEEP);
+	status = pmSourceFitModel (source, model, true);
+	psImageKeepCircle (source->mask, x, y, model->radius, AND, ~PSPHOT_MASK_KEEP);
+	if (!status || (model->params->data.F32[1] < 0)) {
+	  psLogMsg ("psphot", 3, "PSF fit failed for %f, %f (%d iterations)\n", x, y, model->nIter);
+	  source->type = PS_SOURCE_FAIL_FIT_PSF;  // better choice?
+	  psFree (model);
+	  continue;
+	}
+
+	// model succeeded : tentatively keep it
+	source->modelPSF = model;
+	Niter += model[0].nIter;
+	Nfit ++;
+
+	// is it a good model?
+	psphotMarkPSF (source, PSF_SHAPE_NSIGMA, PSF_MIN_SN, PSF_MAX_CHI, SATURATION);
+	if (psphotSubtractPSF (source)) {
+	  Nsub ++;
+	}
+    }
+
+    psLogMsg ("psphot", 3, "fit PSF models: %f sec for %d objects (%d total iterations)\n", psTimerMark ("psphot"), Nfit, Niter);
+    psLogMsg ("psphot", 3, "subtracted %d PSF objects\n", Nsub);
+    return (true);
+}
Index: /tags/eam_tags/psphot_dev_01/psphot/src/psphotArguments.c
===================================================================
--- /tags/eam_tags/psphot_dev_01/psphot/src/psphotArguments.c	(revision 22334)
+++ /tags/eam_tags/psphot_dev_01/psphot/src/psphotArguments.c	(revision 22334)
@@ -0,0 +1,73 @@
+# include "psphot.h"
+
+static void usage (void) ;
+
+psMetadata *psphotArguments (int *argc, char **argv) {
+
+  int N, Nfail;
+  int mode = PS_META_STR | PS_META_REPLACE;
+
+  // basic pslib options
+  fprintf (stderr, "starting...\n");
+  psLogSetFormat ("M");
+  psLogArguments (argc, argv);
+  psTraceArguments (argc, argv);
+
+  // optional mask image - add to config
+  char *mask = NULL;
+  if ((N = psArgumentGet (*argc, argv, "-mask"))) {
+    psArgumentRemove (N, argc, argv);
+    mask = psStringCopy (argv[N]);
+    psArgumentRemove (N, argc, argv);
+  }
+
+  // optional noise image - add to config
+  char *noise = NULL;
+  if ((N = psArgumentGet (*argc, argv, "-noise"))) {
+    psArgumentRemove (N, argc, argv);
+    noise = psStringCopy (argv[N]);
+    psArgumentRemove (N, argc, argv);
+  }
+
+  // optional output residual image - add to config
+  char *resid = NULL;
+  if ((N = psArgumentGet (*argc, argv, "-resid"))) {
+    psArgumentRemove (N, argc, argv);
+    resid = psStringCopy (argv[N]);
+    psArgumentRemove (N, argc, argv);
+  }
+
+  if (*argc != 4) usage ();
+
+  // load config information
+  psMetadata *config = psMetadataAlloc ();
+  psMetadataAdd (config, PS_LIST_HEAD, "PSF_MODEL", PS_META_MULTI, "folder for psf model entries", NULL);
+  config = psMetadataParseConfig (config, &Nfail, argv[3], FALSE);
+  fprintf (stderr, "loaded config...\n");
+
+  // identify input image & optional noise & mask images
+  // command-line entries override config-file entries
+  psMetadataAdd (config, PS_LIST_HEAD, "IMAGE",       mode, "", argv[1]);
+  psMetadataAdd (config, PS_LIST_HEAD, "OUTPUT_FILE", mode, "", argv[2]);
+
+  if (mask != NULL) {
+    psMetadataAdd (config, PS_LIST_HEAD, "MASK_IMAGE", mode, "", mask);
+  }
+  if (noise != NULL) {
+    psMetadataAdd (config, PS_LIST_HEAD, "NOISE_IMAGE", mode, "", noise);
+  }
+  if (resid != NULL) {
+    psMetadataAdd (config, PS_LIST_HEAD, "RESID_IMAGE", mode, "", resid);
+  }
+  return (config);
+}
+
+static void usage (void) {
+
+    fprintf (stderr, "USAGE: psphot (image) (output) (config)\n");
+    fprintf (stderr, "options: \n");
+    fprintf (stderr, "  -mask  (filename)\n");
+    fprintf (stderr, "  -noise (filename)\n");
+    fprintf (stderr, "  -resid (filename)\n");
+    exit (2);
+}
Index: /tags/eam_tags/psphot_dev_01/psphot/src/psphotChoosePSF.c
===================================================================
--- /tags/eam_tags/psphot_dev_01/psphot/src/psphotChoosePSF.c	(revision 22334)
+++ /tags/eam_tags/psphot_dev_01/psphot/src/psphotChoosePSF.c	(revision 22334)
@@ -0,0 +1,94 @@
+# include "psphot.h"
+
+// test PSF models and select best option
+
+pmPSF *psphotChoosePSF (psMetadata *config, psArray *sources, psStats *skystats) 
+{ 
+    bool       	    status;
+    char       	   *modelName;
+    pmPSF_Test 	   *test  = NULL;
+    psArray    	   *stars = NULL;
+    psMetadataItem *item  = NULL;
+
+    // array to store candidate PSF stars
+    stars = psArrayAlloc (sources->n);
+    stars->n = 0;
+
+    // select the candidate PSF stars (pointers to original sources)
+    for (int i = 0; i < sources->n; i++) {
+	psSource *source = sources->data[i];
+	if (source->type != PS_SOURCE_PSFSTAR) continue;
+	psArrayAdd (stars, 200, source);
+    }
+    psTrace (".psphot.pspsf", 3, "selected candidate %d PSF objects\n", stars->n);
+
+    // get the fixed PSF fit radius
+    // XXX EAM : check that PSF_FIT_RADIUS < SKY_OUTER_RADIUS
+    float RADIUS = psMetadataLookupF32 (&status, config, "PSF_FIT_RADIUS");
+
+    // get the list pointers for the PSF_MODEL entries
+    psMetadataItem *mdi = psMetadataLookup (config, "PSF_MODEL");
+    if (mdi == NULL) psAbort ("psphotChoosePSF", "missing PSF_MODEL selection");
+
+    psList *list = (psList *) mdi->data.list;
+    psListIterator *iter = psListIteratorAlloc (list, PS_LIST_HEAD, FALSE); 
+
+    // set up an array to store the test results
+    psArray *tests = psArrayAlloc (list->size);
+    // XXX EAM : new value : psArray *tests = psArrayAlloc (list->n);
+
+    // test each model option listed in config
+    for (int i = 0; i < tests->n; i++) { 
+
+	item = psListGetAndIncrement (iter);
+	modelName = item->data.V;
+
+	tests->data[i] = pmPSF_TestModel (stars, modelName, RADIUS);
+	psFree (modelName);
+	psFree (item);
+    }
+    psFree (iter);
+    // psFree (list); XXX EAM - is list freed with iter?
+    psFree (stars);
+
+    // select the best of the test models
+    // here we are using the clippedStdev on the metric as the indicator
+    test = tests->data[0];
+    int bestN = 0;
+    float bestM = test->dApResid;
+    for (int i = 1; i < tests->n; i++) {
+	test = tests->data[i];
+	float M = test->dApResid;
+	if (M < bestM) {
+	    bestM = M;
+	    bestN = i;
+	}
+    }
+
+    // keep only the selected test:
+    test = tests->data[bestN];
+    modelName = psModelGetType (test->modelType);
+    psLogMsg ("psphot.pspsf", 3, "selected psf model %s, ApResid: %f +/- %f\n", modelName, test->ApResid, test->dApResid);
+
+    // XXX I am having trouble tracking down the double frees vs leaks
+    # if (0)
+    // we are keeping the modelPSF fits from the PSF test, 
+    // but we replace these later on (apply_psf_model)
+    // set source->model based on best psf model fit
+    for (int i = 0; i < test->sources->n; i++) {
+	psSource *source = test->sources->data[i];
+	// drop masked sources from PSFSTAR list
+	if (test->mask->data.U8[i]) {
+	  source->type = PS_SOURCE_OTHER;
+	  source->modelPSF = NULL;
+	} else {
+	  // source->modelPSF = psMemCopy(test->modelPSF->data[i]);
+	  source->modelPSF = test->modelPSF->data[i];
+	}
+    }
+    # endif
+
+    pmPSF *psf = psMemCopy(test->psf);
+    psFree (tests);
+    return (psf);
+}
Index: /tags/eam_tags/psphot_dev_01/psphot/src/psphotFitGalaxies.c
===================================================================
--- /tags/eam_tags/psphot_dev_01/psphot/src/psphotFitGalaxies.c	(revision 22334)
+++ /tags/eam_tags/psphot_dev_01/psphot/src/psphotFitGalaxies.c	(revision 22334)
@@ -0,0 +1,114 @@
+# include "psphot.h"
+
+bool psphotFitGalaxies (psImageData *imdata, psMetadata *config, psArray *sources, psStats *skyStats) 
+{ 
+    bool  status, goodfit;
+    float x;
+    float y;
+    float sky;
+    int   Nfit = 0;
+    int   Nfail = 0;
+    int   Niter = 0;
+
+    float OUTER_RADIUS	   = psMetadataLookupF32 (&status, config, "SKY_OUTER_RADIUS");
+
+    float GAL_MIN_SN  	   = psMetadataLookupF32 (&status, config, "GAL_MIN_SN");
+    float GAL_FIT_NSIGMA   = psMetadataLookupF32 (&status, config, "GAL_FIT_NSIGMA");
+    float GAL_FIT_PADDING  = psMetadataLookupF32 (&status, config, "GAL_FIT_PADDING");
+    float GAL_MOMENTS_RAD  = psMetadataLookupF32 (&status, config, "GAL_MOMENTS_RADIUS");
+
+    float FLUX_LIMIT       = GAL_FIT_NSIGMA * skyStats->sampleStdev;
+
+    char         *modelName   = psMetadataLookupPtr (&status, config, "GAL_MODEL");
+    psModelType   modelType   = psModelSetType (modelName);
+    psModelRadius modelRadius = psModelRadius_GetFunction (modelType);
+
+    psTimerStart ("psphot");
+
+    for (int i = 0; i < sources->n; i++) {
+	psSource *source = sources->data[i];
+
+	// sources which should not be fitted
+	// skip all valid stars
+	if (source->type == PS_SOURCE_PSFSTAR) continue;
+	if (source->type == PS_SOURCE_SATSTAR) continue;
+	if (source->type == PS_SOURCE_GOODSTAR) continue;
+
+	// skip all likely defects
+	if (source->type == PS_SOURCE_DEFECT) continue;
+	if (source->type == PS_SOURCE_SATURATED) continue;
+
+	// skip poorly fitted stars
+	if (source->type == PS_SOURCE_FAINTSTAR) continue;
+	if (source->type == PS_SOURCE_POOR_FIT_PSF) continue;
+
+	// XXX EAM when do we pick these up again?
+	if (source->moments->SN < GAL_MIN_SN) {
+	  source->type = PS_SOURCE_FAINT_GALAXY;  // better choice?
+	  continue;
+	}
+
+	// recalculate the source moments using the larger galaxy moments radius
+	status = pmSourceMoments_EAM (source, GAL_MOMENTS_RAD);
+	if (!status) {
+	  source->type = PS_SOURCE_DROP_GALAXY;  // better choice?
+	  continue;
+	}
+
+	// use the source moments, etc to guess basic model parameters
+	psModel  *model  = pmSourceModelGuess (source, modelType); 
+	x = model->params->data.F32[2];
+	y = model->params->data.F32[3];
+
+	// set the fit radius based on the object flux limit and the model
+	// XXX EAM : FLUX_LIMIT should be set based on local sky model (not global median)
+	model->radius = modelRadius (model->params, FLUX_LIMIT) + GAL_FIT_PADDING;
+	if (isnan(model->radius)) psAbort ("fit_galaxies", "error in radius");
+
+	// check if we need to redefine the pixels
+	// XXX EAM : a better test would examine the source pixels
+	if (model->radius > OUTER_RADIUS) {
+	  // (re)-allocate image, noise, mask arrays for each peak (square of radius OUTER)
+	  pmSourceDefinePixels (source, imdata, x, y, model->radius);
+	}
+
+	// fit FLT (not PSF) model (set/unset the pixel mask)
+	psImageKeepCircle (source->mask, x, y, model->radius, OR, PSPHOT_MASK_KEEP);
+	status = pmSourceFitModel_EAM (source, model, false);
+	psImageKeepCircle (source->mask, x, y, model->radius, AND, ~PSPHOT_MASK_KEEP);
+	if (!status) {
+	  // if the fit fails, we need to change the classification
+	  psLogMsg ("psphot", 3, "GAL fit failed for %f, %f (%d iterations, %f radius)\n", x, y, model->nIter, model->radius);
+	  source->type = PS_SOURCE_FAIL_FIT_GAL;  // better choice?
+	  source->modelFLT = model;
+	  Nfail ++;
+	  continue;
+	}
+
+	// check if model fit is acceptable
+	goodfit = pmModelFitStatus (model);
+	if (!goodfit) {
+	  // if the fit fails, we need to change the classification
+	  psLogMsg ("psphot", 3, "GAL fit poor for %f, %f (%d iterations, %f radius)\n", x, y, model->nIter, model->radius);
+	  source->type = PS_SOURCE_POOR_FIT_GAL;  // better choice?
+	  source->modelFLT = model;
+	  Nfail ++;
+	  continue;
+	}
+
+	// accept the model
+	source->type = PS_SOURCE_GALAXY;
+	source->modelFLT = model;
+	Niter += model[0].nIter;
+	Nfit++;
+
+	// subtract object, leave local sky
+	sky = model->params->data.F32[0];
+	model->params->data.F32[0] = 0;
+	pmSourceSubModel (source->pixels, source->mask, model, false);
+	model->params->data.F32[0] = sky;
+
+    }
+    psLogMsg ("psphot", 3, "fit galaxies: %f sec for %d galaxies (%d failures, %d total iterations)\n", psTimerMark ("psphot"), Nfit, Nfail, Niter);
+    return (true);
+}
Index: /tags/eam_tags/psphot_dev_01/psphot/src/psphotImageStats.c
===================================================================
--- /tags/eam_tags/psphot_dev_01/psphot/src/psphotImageStats.c	(revision 22334)
+++ /tags/eam_tags/psphot_dev_01/psphot/src/psphotImageStats.c	(revision 22334)
@@ -0,0 +1,56 @@
+# include "psphot.h"
+
+psStats *psphotImageStats (psImageData *imdata, psMetadata *config) 
+{ 
+    psStats *stats = NULL;
+    psStats *sky   = NULL;
+
+    psImage *image = imdata->image;
+
+    // get image stats on a subset of the image (random 1e5 pts)
+    psTimerStart ("psphot");
+
+    // this section fixes a bug in the psImageStats definition:
+    // it should limit the number of elements used for the initial 
+    // pass on the stats
+    { 
+	bool     status  = false;
+	int      Npixels = image->numRows*image->numCols;
+	int      Nsubset = PS_MIN (Npixels, psMetadataLookupF32 (&status, config, "IMSTATS_NPIX"));
+	psVector *subset = psVectorAlloc (Nsubset, PS_TYPE_F32);
+	psRandom *rnd    = psRandomAlloc (PS_RANDOM_TAUS, 0);
+
+	// choose Nsubset points between 0 and Nx*Ny, convert to coords
+	Npixels = image->numRows*image->numCols;
+	for (int i = 0; i < Nsubset; i++) {
+	  double frnd = psRandomUniform (rnd);
+	  int pixel = Npixels * frnd;
+	  int ix = pixel / image->numCols;
+	  int iy = pixel % image->numCols;
+	  subset->data.F32[i] = image->data.F32[iy][ix];
+	}
+
+	// this should use ROBUST not SAMPLE median
+	// robust median is broken in pslib (0.5)
+	stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN); 
+	stats = psVectorStats (stats, subset, NULL, NULL, 0); // use a mask?
+	psFree (subset);
+	psFree (rnd);
+    }
+
+    // we store these values in mean,stdev 
+    bool status = false;
+    float NOISE = psMetadataLookupF32 (&status, config, "RDNOISE");
+    float GAIN  = psMetadataLookupF32 (&status, config, "GAIN");
+
+    // convert instrumental background to poisson stats
+    sky = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
+    sky->sampleMean   = stats->sampleMedian;
+    sky->sampleStdev  = sqrt(sky->sampleMean/GAIN + PS_SQR(NOISE));
+
+    psLogMsg ("psphot", 3, "stats: %f sec\n", psTimerMark ("psphot"));
+    psLogMsg ("psphot", 3, "background: %f +/- %f\n", sky->sampleMean, sky->sampleStdev);
+
+    psFree (stats);
+    return (sky);
+}
Index: /tags/eam_tags/psphot_dev_01/psphot/src/psphotMarkPSF.c
===================================================================
--- /tags/eam_tags/psphot_dev_01/psphot/src/psphotMarkPSF.c	(revision 22334)
+++ /tags/eam_tags/psphot_dev_01/psphot/src/psphotMarkPSF.c	(revision 22334)
@@ -0,0 +1,96 @@
+# include "psphot.h"
+
+// identify objects consistent with PSF shape/magnitude distribution
+// we expect dparams[4],dparams[5] to have a scatter of:
+// sigma_x / (S/N) * sqrt(2)
+// 1 / (params[4],params[5])*(S/N)
+
+// sigma_x : 1 / SX
+// dsx : 1 / (SX * SN)
+// dsx_o = hypot (1/(SX*SN), MIN_DSX)
+
+// any objects which is consistent with the PSF should have 
+// abs(dparams[5]) < N * dsxLine(mag) & abs(dparams[6]) < N * dsyLine(mag)
+// this includes a minimum buffer (DS) for the brighter objects
+
+// saturated stars should fall outside (but are already IDed)
+// galaxies should be larger, cosmic rays smaller, but need to test?
+// we also reject objects with S/N too low or ChiSquare to high
+
+// any object which meets the criterion is marked as PS_SOURCE_BRIGHTSTAR 
+
+// floor for DS value 
+// XXX EAM : add to configuration?
+# define MIN_DS 0.01
+
+bool psphotMarkPSF (psSource *source, float shapeNsigma, float minSN, float maxChi, float SATURATE)
+{ 
+    int keep;
+    float dSX, dSY, SX, SY, SN;
+    float nSx, nSy, Chi;
+
+    if (source->modelPSF == NULL) return (false);
+
+    // if object has fitted peak above saturation, label as SATSTAR
+    // remember: fit does not use saturated pixels (masked)
+    if (source->modelPSF->params->data.F32[1] >= SATURATE) {
+	if (source->type == PS_SOURCE_PSFSTAR) {
+	    psLogMsg ("psphot", 3, "PSFSTAR marked SATSTAR\n");
+	}
+	source->type = PS_SOURCE_SATSTAR;
+	return (true);
+    } 
+    if (source->type == PS_SOURCE_SATSTAR) {
+	psLogMsg ("psphot", 4, "SATSTAR marked GOODSTAR (fitted peak below saturation)\n");
+	source->type = PS_SOURCE_GOODSTAR;
+    }
+
+    SN  = source->modelPSF->params->data.F32[1]/source->modelPSF->dparams->data.F32[1];
+    SX  = source->modelPSF->params->data.F32[4];
+    SY  = source->modelPSF->params->data.F32[5];
+    dSX = source->modelPSF->dparams->data.F32[4];
+    dSY = source->modelPSF->dparams->data.F32[5];
+    Chi = source->modelPSF->chisq / source->modelPSF->nDOF;
+
+    // swing of sigma_x,y in sigmas
+    nSx = dSX / hypot (MIN_DS, 1 / (SX * SN));
+    nSy = dSY / hypot (MIN_DS, 1 / (SY * SN));
+
+    // assign PS_SOURCE_GOODSTAR to bright objects within PSF region of dparams[]
+    keep = TRUE;
+    keep &= (fabs(nSx) < shapeNsigma);
+    keep &= (fabs(nSy) < shapeNsigma);
+    keep &= (SN > minSN);
+    keep &= (Chi < maxChi);
+    if (keep) {
+	if (source->type == PS_SOURCE_PSFSTAR) return (true);
+	source->type = PS_SOURCE_GOODSTAR;
+	return (true);
+    }
+    
+    if (source->type == PS_SOURCE_PSFSTAR) {
+	psLogMsg ("psphot", 3, "PSFSTAR demoted based on fit quality\n");
+    }
+
+    // object appears to be small, suspected defect
+    if ((nSx <= -shapeNsigma) || (nSy <= -shapeNsigma)) {
+	source->type = PS_SOURCE_DEFECT;
+	return (false);
+    }
+
+    // object appears to be large, suspected galaxy
+    if ((nSx >= shapeNsigma) || (nSy >= shapeNsigma)) {
+	source->type = PS_SOURCE_GALAXY;
+	return (false);
+    }
+
+    // object appears to be extremely faint: what is this?
+    if (SN < minSN) {
+	source->type = PS_SOURCE_FAINTSTAR;
+	return (false);
+    }
+
+    // these are pooly fitted, probable stars near other stars?
+    source->type = PS_SOURCE_POOR_FIT_PSF;
+    return (false);
+}	
Index: /tags/eam_tags/psphot_dev_01/psphot/src/psphotOutput.c
===================================================================
--- /tags/eam_tags/psphot_dev_01/psphot/src/psphotOutput.c	(revision 22334)
+++ /tags/eam_tags/psphot_dev_01/psphot/src/psphotOutput.c	(revision 22334)
@@ -0,0 +1,643 @@
+# include "psphot.h"
+
+// output functions: we have several fixed modes (sx, obj, cmp)
+void psphotOutput (psImageData *imdata, psMetadata *config, psArray *sources) {
+
+    bool status;
+
+    char *outputMode = psMetadataLookupPtr (&status, config, "OUTPUT_MODE");
+    char *outputFile = psMetadataLookupPtr (&status, config, "OUTPUT_FILE");
+    char *residImage = psMetadataLookupPtr (&status, config, "RESID_IMAGE");
+
+    if (residImage != NULL) psphotSaveImage (imdata->header, imdata->image, residImage);
+
+    if (outputFile == NULL) {
+	psLogMsg ("output", 3, "no data output file selected");
+	return;
+    }
+    if (outputMode == NULL) {
+	psLogMsg ("output", 3, "no data output mode selected");
+	return;
+    }
+    if (!strcasecmp (outputMode, "TEXT")) {
+	pmSourcesWriteText (imdata, outputFile, sources);
+	return;
+    }
+
+    if (!strcasecmp (outputMode, "OBJ")) {
+	pmSourcesWriteOBJ (imdata, outputFile, sources);
+	return;
+    }
+  
+    if (!strcasecmp (outputMode, "SX")) {
+	pmSourcesWriteSX (imdata, outputFile, sources);
+	return;
+    }
+  
+    if (!strcasecmp (outputMode, "CMP")) {
+	pmSourcesWriteCMP (imdata, outputFile, sources);
+	return;
+    }
+  
+    if (!strcasecmp (outputMode, "CMF")) {
+	pmSourcesWriteCMF (imdata, outputFile, sources);
+	return;
+    }
+    psAbort ("psphot", "unknown output mode %s", outputMode);
+}
+
+bool pmSourcesWriteText (psImageData *imdata, char *filename, psArray *sources) {
+
+    char *name = (char *) psAlloc (strlen(filename) + 10);
+
+    sprintf (name, "%s.psf.dat", filename);
+    psModelWritePSFs (sources, name);
+
+    sprintf (name, "%s.flt.dat", filename);
+    psModelWriteFLTs (sources, name);
+
+    sprintf (name, "%s.nul.dat", filename);
+    psModelWriteNULLs (sources, name);
+
+    sprintf (name, "%s.mnt.dat", filename);
+    psMomentsWriteText (sources, name);
+
+    psFree (name);
+    return true;
+}
+
+// dophot-style output list with fixed line width
+bool pmSourcesWriteOBJ (psImageData *imdata, char *filename, psArray *sources) {
+
+    int i, type, status;
+    psModel *model;
+    psF32 *PAR, *dPAR;
+    float sky, dmag, apMag, fitMag;
+
+    psLine *line = psLineAlloc (104);  // 104 is dophot-defined line length
+
+    FILE *f = fopen (filename, "w");
+    if (f == NULL) {
+	psLogMsg ("WriteSourceOBJ", 3, "can't open output file for output %s\n", filename);
+	return false;
+    }
+
+    // write sources with models first 
+    for (i = 0; i < sources->n; i++) {
+	psSource *source = (psSource *) sources->data[i];
+
+	// use the correct model
+	switch (source->type) {
+
+	    // use PSF model with stars
+	  case PS_SOURCE_PSFSTAR:
+	  case PS_SOURCE_SATSTAR:
+	  case PS_SOURCE_GOODSTAR:
+	    model = source->modelPSF;
+	    break;
+
+	    // use FLT model with galaxies
+	  case PS_SOURCE_GALAXY:
+	    model = source->modelFLT;
+	    break;
+	    
+	    // skip defects and poorly fitted stars or galaxies
+	  case PS_SOURCE_DEFECT:
+	  case PS_SOURCE_SATURATED:
+	  case PS_SOURCE_FAINTSTAR:
+	  case PS_SOURCE_POOR_FIT_PSF:
+	  case PS_SOURCE_POOR_FIT_GAL:
+	  case PS_SOURCE_FAIL_FIT_GAL:
+	  default:
+	    model = NULL;
+	}
+	if (model == NULL) continue;
+
+	PAR = model->params->data.F32;
+	dPAR = model->dparams->data.F32;
+
+	// save local sky for later
+	sky = model->params->data.F32[0];
+
+	// replace source flux
+	// XXX EAM : need to add 'sky' boolean option here 
+	model->params->data.F32[0] = 0;
+	pmSourceAddModel (source->pixels, source->mask, model, false);
+	model->params->data.F32[0] = sky;
+
+	// measure object photometry
+	status = pmSourcePhotometry (&fitMag, &apMag, model, source->pixels, source->mask);
+
+	// subtract object, leave local sky
+	// XXX EAM : need to add 'sky' boolean option here 
+	model->params->data.F32[0] = 0;
+	pmSourceSubModel (source->pixels, source->mask, model, false);
+	model->params->data.F32[0] = sky;
+
+	if (status == FALSE) continue;
+
+	dmag = dPAR[1] / PAR[1];
+	type = pmSourceDophotType (source);
+	psLineInit (line);
+	psLineAdd (line, "%3d", source->type);
+	psLineAdd (line, "%8.2f", PAR[2]);
+	psLineAdd (line, "%8.2f", PAR[3]);
+	psLineAdd (line, "%8.3f", fitMag);
+	psLineAdd (line, "%6.3f", dmag);
+	psLineAdd (line, "%9.2f", PAR[0]);
+	psLineAdd (line, "%9.3f", PAR[4]);
+	psLineAdd (line, "%9.3f", PAR[5]);
+	psLineAdd (line, "%7.2f", PAR[6]);
+	psLineAdd (line, "%8.3f", 32.0);
+	psLineAdd (line, "%8.3f", apMag);
+	psLineAdd (line, "%8.2f\n", apMag - fitMag);
+	fwrite (line->line, 1, line->Nline, f);
+    }
+    fclose (f);
+    return true;
+}
+
+// elixir/sextractor-style output list with fixed line width
+bool pmSourcesWriteSX (psImageData *imdata, char *filename, psArray *sources) {
+
+    int i, status, flags;
+    psModel *model;
+    psF32 *PAR, *dPAR;
+    float sky, dmag, apMag, fitMag;
+
+    psLine *line = psLineAlloc (110);  // 104 is dophot-defined line length
+
+    FILE *f = fopen (filename, "w");
+    if (f == NULL) {
+	psLogMsg ("WriteSourceOBJ", 3, "can't open output file for output %s\n", filename);
+	return false;
+    }
+
+    // write sources with models first 
+    for (i = 0; i < sources->n; i++) {
+	psSource *source = (psSource *) sources->data[i];
+	flags = 0;
+
+	// use the correct model
+	switch (source->type) {
+
+	    // use PSF model with stars
+	  case PS_SOURCE_SATSTAR:
+	    flags = 4;
+	    model = source->modelPSF;
+	    break;
+
+	  case PS_SOURCE_PSFSTAR:
+	  case PS_SOURCE_GOODSTAR:
+	    model = source->modelPSF;
+	    break;
+
+	    // use FLT model with galaxies
+	  case PS_SOURCE_GALAXY:
+	    model = source->modelFLT;
+	    break;
+	    
+	    // skip defects and poorly fitted stars or galaxies
+	  case PS_SOURCE_DEFECT:
+	  case PS_SOURCE_SATURATED:
+	  case PS_SOURCE_FAINTSTAR:
+	  case PS_SOURCE_POOR_FIT_PSF:
+	  case PS_SOURCE_POOR_FIT_GAL:
+	  case PS_SOURCE_FAIL_FIT_GAL:
+	  default:
+	    model = NULL;
+	}
+	if (model == NULL) continue;
+
+	PAR = model->params->data.F32;
+	dPAR = model->dparams->data.F32;
+
+	// save local sky for later
+	sky = model->params->data.F32[0];
+
+	// replace source flux
+	// XXX EAM : need to add 'sky' boolean option here 
+	model->params->data.F32[0] = 0;
+	pmSourceAddModel (source->pixels, source->mask, model, false);
+	model->params->data.F32[0] = sky;
+
+	// measure object photometry
+	status = pmSourcePhotometry (&fitMag, &apMag, model, source->pixels, source->mask);
+
+	// subtract object, leave local sky
+	// XXX EAM : need to add 'sky' boolean option here 
+	model->params->data.F32[0] = 0;
+	pmSourceSubModel (source->pixels, source->mask, model, false);
+	model->params->data.F32[0] = sky;
+
+	if (status == FALSE) continue;
+
+	// XXX EAM : fix this to match sextractor types?
+	// type = pmSourceSextractType (source);
+
+	dmag = dPAR[1] / PAR[1];
+	psLineInit (line);
+	psLineAdd (line, "%5.2f", 0.0); // should be type
+	psLineAdd (line, "%11.3f", PAR[2]);
+	psLineAdd (line, "%11.3f", PAR[3]);
+	psLineAdd (line, "%9.4f", fitMag);
+	psLineAdd (line, "%9.4f", dmag);
+	psLineAdd (line, "%13.4f", PAR[0]);
+	psLineAdd (line, "%9.2f", 0.0); // should be FWHMx
+	psLineAdd (line, "%9.2f", 0.0); // should be FWHMy
+	psLineAdd (line, "%6.1f", 0.0); // should be Theta
+	psLineAdd (line, "%9.4f", 32.0); // should be MAG_ISO
+	psLineAdd (line, "%9.4f", apMag);
+	psLineAdd (line, "%4d\n", 0); // should be flags
+	fwrite (line->line, 1, line->Nline, f);
+    }
+    fclose (f);
+    return true;
+}
+
+// elixir-style pseudo FITS table (header + ascii list)
+bool pmSourcesWriteCMP (psImageData *imdata, char *filename, psArray *sources) {
+
+    int i, type, status;
+    psModel *model;
+    psF32 *PAR, *dPAR;
+    float sky, dmag, apMag, fitMag, lsky, Theta;
+
+    // create file, write-out header
+    unlink (filename);
+    psFits *fits = psFitsAlloc (filename);
+    psFitsWriteHeader (imdata->header, fits);
+    psFree (fits);
+
+    // re-open, add data to end of file
+    FILE *f = fopen (filename, "a+");
+    if (f == NULL) {
+	psLogMsg ("WriteSourceOBJ", 3, "can't open output file for output %s\n", filename);
+	return false;
+    }
+    fseek (f, 0, SEEK_END);
+
+    psLine *line = psLineAlloc (104);  // 104 is dophot-defined line length
+
+    // write sources with models first 
+    for (i = 0; i < sources->n; i++) {
+	psSource *source = (psSource *) sources->data[i];
+
+	// use the correct model
+	switch (source->type) {
+
+	    // use PSF model with stars
+	  case PS_SOURCE_PSFSTAR:
+	  case PS_SOURCE_SATSTAR:
+	  case PS_SOURCE_GOODSTAR:
+	    model = source->modelPSF;
+	    break;
+
+	    // use FLT model with galaxies
+	  case PS_SOURCE_GALAXY:
+	    model = source->modelFLT;
+	    break;
+	    
+	    // skip defects and poorly fitted stars or galaxies
+	  case PS_SOURCE_DEFECT:
+	  case PS_SOURCE_SATURATED:
+	  case PS_SOURCE_FAINTSTAR:
+	  case PS_SOURCE_POOR_FIT_PSF:
+	  case PS_SOURCE_POOR_FIT_GAL:
+	  case PS_SOURCE_FAIL_FIT_GAL:
+	  default:
+	    model = NULL;
+	}
+	if (model == NULL) continue;
+
+	PAR = model->params->data.F32;
+	dPAR = model->dparams->data.F32;
+
+	// save local sky for later
+	sky = model->params->data.F32[0];
+
+	// XXX EAM : get Theta from PAR[4,5,6]
+	Theta = 0.0;
+
+	// replace source flux
+	// XXX EAM : need to add 'sky' boolean option here 
+	model->params->data.F32[0] = 0;
+	pmSourceAddModel (source->pixels, source->mask, model, false);
+	model->params->data.F32[0] = sky;
+
+	// measure object photometry
+	status = pmSourcePhotometry (&fitMag, &apMag, model, source->pixels, source->mask);
+
+	// subtract object, leave local sky
+	// XXX EAM : need to add 'sky' boolean option here 
+	model->params->data.F32[0] = 0;
+	pmSourceSubModel (source->pixels, source->mask, model, false);
+	model->params->data.F32[0] = sky;
+
+	if (status == FALSE) continue;
+
+	dmag = dPAR[1] / PAR[1];
+	type = pmSourceDophotType (source);
+	lsky = (PAR[0] < 1.0) ? 0.0 : log10(PAR[0]);
+
+	psLineInit (line);
+	psLineAdd (line, "%6.1f ", PAR[2]);
+	psLineAdd (line, "%6.1f ", PAR[3]);
+	psLineAdd (line, "%6.3f ", fitMag + 25.0);
+	psLineAdd (line, "%03d ", (int)(1000*dmag));
+	psLineAdd (line, "%2d ",   type);
+	psLineAdd (line, "%3.1f ", lsky);
+	psLineAdd (line, "%6.3f ", 32.0);    // should be 'Mgal
+	psLineAdd (line, "%6.3f ", apMag + 25.0);
+	psLineAdd (line, "%6.2f ", PAR[4]);  // should be 'FHWM x'
+	psLineAdd (line, "%6.2f ", PAR[5]);  // should be 'FHWM y'
+	psLineAdd (line, "%5.1f\n", Theta);   // should be theta
+	fwrite (line->line, 1, line->Nline, f);
+    }
+    fclose (f);
+    return true;
+}
+
+// elixir-style FITS table output (header + table in 1st extension)
+bool pmSourcesWriteCMF (psImageData *imdata, char *filename, psArray *sources) {
+
+    // write imdata->header to file
+
+    // write FITS table to file (use autocode tools)
+
+    return true;
+}
+
+/***** Text Output Methods *****/
+
+// write the peaks to an output file
+bool psPeaksWriteText (psArray *sources, char *filename) {
+
+    int i;
+    FILE *f;
+    psPeak *peak;
+	
+    f = fopen (filename, "w");
+    if (f == NULL) {
+	psLogMsg ("psPeaksWriteText", 3, "can't open output file for peaks%s\n", filename);
+	return false;
+    }
+
+    for (i = 0; i < sources->n; i++) {
+	peak = (psPeak *) sources->data[i];
+	if (peak == NULL) continue;
+	fprintf (f, "%5d %5d  %7.1f\n", 
+		 peak->x, peak->y, peak->counts); 
+    }
+    fclose (f);
+    return true;
+}
+
+// write the moments to an output file 
+bool psMomentsWriteText (psArray *sources, char *filename) {
+
+    int i;
+    FILE *f;
+    psSource *source;
+	
+    f = fopen (filename, "w");
+    if (f == NULL) {
+	psLogMsg ("psMomentsWriteText", 3, "can't open output file for moments%s\n", filename);
+	return false;
+    }
+
+    for (i = 0; i < sources->n; i++) {
+	source = (psSource *) sources->data[i];
+	if (source->moments == NULL) continue;
+	fprintf (f, "%5d %5d  %7.1f  %7.1f %7.1f  %6.3f %6.3f  %8.1f %7.1f %7.1f %7.1f  %4d %2d\n", 
+		 source->peak->x, source->peak->y, source->peak->counts, 
+		 source->moments->x, source->moments->y, 
+		 source->moments->Sx, source->moments->Sy, 
+		 source->moments->Sum, source->moments->Peak, 
+		 source->moments->Sky, source->moments->SN, 
+		 source->moments->nPixels, source->type); 
+    }
+    fclose (f);
+    return true;
+}
+
+// dump the sources to an output file
+bool psModelWritePSFs (psArray *sources, char *filename) {
+
+    double dP, flux;
+    int i, j;
+    FILE *f;
+    psVector *params;
+    psVector *dparams;
+    psModel  *model;
+
+    f = fopen (filename, "w");
+    if (f == NULL) {
+	psLogMsg ("psModelWritePSFs", 3, "can't open output file for moments%s\n", filename);
+	return false;
+    }
+
+    // write sources with models first 
+    for (i = 0; i < sources->n; i++) {
+	psSource *source = (psSource *) sources->data[i];
+	model = (psModel  *) source->modelPSF;
+	if (model == NULL) continue;
+
+	// valid source types for this function
+	if (source->type == PS_SOURCE_SATSTAR) goto valid;
+	if (source->type == PS_SOURCE_PSFSTAR) goto valid;
+	if (source->type == PS_SOURCE_GOODSTAR) goto valid;
+	continue;
+
+    valid:
+	params = model->params;
+	dparams = model->dparams;
+	
+	dP = 0;
+	dP += PS_SQR(dparams[0].data.F32[2]);
+	dP += PS_SQR(dparams[0].data.F32[3]);
+	dP = sqrt (dP);
+	
+	psModelFlux modelFluxFunc = psModelFlux_GetFunction (model->type);
+	flux = modelFluxFunc (params);
+
+	fprintf (f, "%7.1f %7.1f  %5.1f %7.3f  %7.4f %7.4f  ", 
+		 params[0].data.F32[2], 
+		 params[0].data.F32[3], 
+		 params[0].data.F32[0], 
+		 -2.5*log10(flux), 
+		 (dparams[0].data.F32[1]/params[0].data.F32[1]),
+		 dP);
+	for (j = 0; j < model->params->n - 4; j++) {
+	    fprintf (f, "%9.6f ", params[0].data.F32[j+4]);
+	}
+	fprintf (f, " : ");
+	for (j = 0; j < model->params->n - 4; j++) {
+	    fprintf (f, "%9.6f ", dparams[0].data.F32[j+4]);
+	}
+	fprintf (f, ": %2d %7.3f %7.3f %7.2f %4d %2d\n", 
+		 source[0].type, 
+		 log10(model[0].chisq/model[0].nDOF), 
+		 source[0].moments->SN, 
+		 model[0].radius,
+		 model[0].nDOF, 
+		 model[0].nIter);
+    }
+    fclose (f);
+    return true;
+}
+
+// dump the sources to an output file
+bool psModelWriteFLTs (psArray *sources, char *filename) {
+
+    double dP;
+    int i, j;
+    FILE *f;
+    psVector *params;
+    psVector *dparams;
+    psModel  *model;
+
+    f = fopen (filename, "w");
+    if (f == NULL) {
+	psLogMsg ("psModelWriteFLTs", 3, "can't open output file for moments%s\n", filename);
+	return false;
+    }
+
+    // write sources with models first 
+    for (i = 0; i < sources->n; i++) {
+	psSource *source = (psSource *) sources->data[i];
+	model = (psModel  *) source->modelFLT;
+	if (model == NULL) continue;
+	if (source->type == PS_SOURCE_GALAXY) goto valid;
+	continue;
+
+    valid:
+	params = model->params;
+	dparams = model->dparams;
+
+	// XXX these are hardwired for SGAUSS : this should be pushed into the
+	// model functions as an abstract function 
+	dP = 0;
+	dP += PS_SQR(dparams[0].data.F32[4] / params[0].data.F32[4]);
+	dP += PS_SQR(dparams[0].data.F32[5] / params[0].data.F32[5]);
+	dP += PS_SQR(dparams[0].data.F32[7] / params[0].data.F32[7]);
+	dP = sqrt (dP);
+
+	fprintf (f, "%7.1f %7.1f  %5.1f %7.3f  %7.4f %7.4f  ", 
+		 params[0].data.F32[2], 
+		 params[0].data.F32[3], 
+		 params[0].data.F32[0], 
+		 -2.5*log10(params[0].data.F32[1]), 
+		 (dparams[0].data.F32[1]/params[0].data.F32[1]),
+		 dP);
+
+	for (j = 4; j < model->params->n; j++) {
+	    fprintf (f, "%9.6f ", params[0].data.F32[j]);
+	}
+	fprintf (f, " : ");
+	for (j = 4; j < model->params->n; j++) {
+	    fprintf (f, "%9.6f ", dparams[0].data.F32[j]);
+	}
+	fprintf (f, ": %2d %7.3f %7.3f %7.2f %4d %2d\n", 
+		 source[0].type, log10(model[0].chisq/model[0].nDOF), 
+		 source[0].moments->SN, 
+		 model[0].radius, 
+		 model[0].nDOF, 
+		 model[0].nIter);
+    }
+    fclose (f);
+    return true;
+}
+
+// dump the sources to an output file
+bool psModelWriteNULLs (psArray *sources, char *filename) 
+{
+
+    int i;
+    FILE *f;
+
+    f = fopen (filename, "w");
+    if (f == NULL) {
+	psLogMsg ("DumpObjects", 3, "can't open output file for moments%s\n", filename);
+	return false;
+    }
+
+    psMoments *empty = pmMomentsAlloc ();
+
+    // write sources with models first 
+    for (i = 0; i < sources->n; i++) {
+	psSource *source = (psSource *) sources->data[i];
+
+	// skip these sources (in PSF or FLT)
+	if (source->type == PS_SOURCE_GALAXY) continue;
+	if (source->type == PS_SOURCE_PSFSTAR) continue;
+	if (source->type == PS_SOURCE_SATSTAR) continue;
+	if (source->type == PS_SOURCE_GOODSTAR) continue;
+
+	if (source->moments == NULL) {
+	  fprintf (f, "%5d %5d  %7.1f  %7.1f %7.1f  %6.3f %6.3f  %8.1f %7.1f %7.1f %7.1f  %4d %2d\n", 
+		   source->peak->x, source->peak->y, source->peak->counts, 
+		   empty->x, empty->y, 
+		   empty->Sx, empty->Sy, 
+		   empty->Sum, empty->Peak, 
+		   empty->Sky, empty->SN, 
+		   empty->nPixels, source->type); 
+	} else {
+	  fprintf (f, "%5d %5d  %7.1f  %7.1f %7.1f  %6.3f %6.3f  %8.1f %7.1f %7.1f %7.1f  %4d %2d\n", 
+		   source->peak->x, source->peak->y, source->peak->counts, 
+		   source->moments->x, source->moments->y, 
+		   source->moments->Sx, source->moments->Sy, 
+		   source->moments->Sum, source->moments->Peak, 
+		   source->moments->Sky, source->moments->SN, 
+		   source->moments->nPixels, source->type); 
+	}
+    }
+    fclose (f);
+    return true;
+}
+
+// translations between psphot object types and dophot object types
+int pmSourceDophotType (psSource *source) {
+
+    switch (source->type) {
+
+      case PS_SOURCE_DEFECT:
+      case PS_SOURCE_SATURATED:
+	return (8);
+
+      case PS_SOURCE_SATSTAR:
+	return (10);
+
+      case PS_SOURCE_PSFSTAR:
+      case PS_SOURCE_GOODSTAR:
+	return (1);
+
+      case PS_SOURCE_POOR_FIT_PSF:
+	return (7);
+
+      case PS_SOURCE_FAIL_FIT_PSF:
+      case PS_SOURCE_FAINTSTAR:
+	return (4);
+
+      case PS_SOURCE_GALAXY:
+      case PS_SOURCE_FAINT_GALAXY:
+      case PS_SOURCE_DROP_GALAXY:
+      case PS_SOURCE_FAIL_FIT_GAL:
+      case PS_SOURCE_POOR_FIT_GAL:
+	return (2);
+
+      case PS_SOURCE_OTHER:
+      default:
+	return (0);
+    }
+    return (0);
+}
+
+int psphotSaveImage (psMetadata *header, psImage *image, char *filename) {
+
+    unlink (filename);
+    psFits *fits = psFitsAlloc (filename);
+    psFitsWriteImage (fits, header, image, 0, NULL);
+    psFree (fits);
+    return (TRUE);
+}
Index: /tags/eam_tags/psphot_dev_01/psphot/src/psphotSetup.c
===================================================================
--- /tags/eam_tags/psphot_dev_01/psphot/src/psphotSetup.c	(revision 22334)
+++ /tags/eam_tags/psphot_dev_01/psphot/src/psphotSetup.c	(revision 22334)
@@ -0,0 +1,96 @@
+# include "psphot.h"
+
+// load the image
+// load or construct the mask
+//   apply X/Y border info
+//   apply the SATURATE mask?
+// load or construct the noise 
+
+psImageData *psphotSetup (psMetadata *config) 
+{
+    psMetadata *header = NULL;
+    psImage *image = NULL;
+    psImage *noise = NULL;
+    psImage *mask = NULL;
+    psRegion region = {0,0,0,0};	// a region representing the entire array
+    bool status;
+
+    psTimerStart ("psphot");
+
+    // setup header template, specify COMMENT and HISTORY as MULTI
+    // XXX EAM - shouldn't this be part of the basic FITS concept?
+    header = psMetadataAlloc ();
+    psMetadataAdd (header, PS_LIST_HEAD, "COMMENT", PS_META_MULTI, "folder for comment", NULL);
+    psMetadataAdd (header, PS_LIST_HEAD, "HISTORY", PS_META_MULTI, "folder for history", NULL);
+
+    // read header and image data from INPUT 
+    char *input = psMetadataLookupPtr (&status, config, "IMAGE");
+    psFits *file = psFitsAlloc (input);
+    header = psFitsReadHeader (header, file);
+    image = psFitsReadImage (NULL, file, region, 0);
+    psFree (file);
+    // psFree (input);
+    // XXX EAM - does the current API expect to get a copy or a view?
+
+    // grab these values from the approrpiate location (image header if necessary)
+    float RDNOISE  = pmConfigLookupF32 (&status, config, header, "RDNOISE");
+    float GAIN     = pmConfigLookupF32 (&status, config, header, "GAIN");
+    psMetadataAdd (config, PS_LIST_TAIL, "RDNOISE", PS_META_F32 | PS_META_REPLACE, "read noise value used", RDNOISE);
+    psMetadataAdd (config, PS_LIST_TAIL, "GAIN",    PS_META_F32 | PS_META_REPLACE, "gain value used", GAIN);
+
+    // load the noise image if it is supplied, otherwise build from input
+    char *noiseName = psMetadataLookupPtr (&status, config, "NOISE");
+    if (status == true) {
+	file = psFitsAlloc (noiseName);
+	noise = psFitsReadImage  (NULL, file, region, 0);
+	psFree (file);
+	// psFree (noiseName); XXX - see psFree (input)
+    } else {
+	// build the noise image from the input image, RDNOISE, and GAIN
+	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);
+	    }
+	}
+    }
+
+    // load the mask if specified, otherwise construct 
+    char *maskName = psMetadataLookupPtr (&status, config, "MASK");
+    if (status == true) {
+	// XXX EAM require / convert mask to psU8?
+	file = psFitsAlloc (maskName);
+	mask  = psFitsReadImage  (NULL, file, region, 0);
+	psFree (file);
+	// psFree (maskName); XXX - see psFree (input)
+    } else {
+	mask = psImageAlloc (image->numCols, image->numRows, PS_TYPE_U8);
+    }
+
+    // mask the excluded outer pixels
+    float XMIN  = psMetadataLookupF32 (&status, config, "XMIN");
+    float XMAX  = psMetadataLookupF32 (&status, config, "XMAX");
+    float YMIN  = psMetadataLookupF32 (&status, config, "YMIN");
+    float YMAX  = psMetadataLookupF32 (&status, config, "YMAX");
+    psRegion *keep = psRegionAlloc (XMIN, XMAX, YMIN, YMAX);
+    keep = psRegionForImage (keep, image, keep);
+    psImageKeepRegion (mask, keep, OR, PSPHOT_MASK_INVALID);
+    psFree (keep);
+
+    // mask the saturated pixels
+    // XXX EAM does the mask need to grow?
+    float SATURATION = psMetadataLookupF32 (&status, config, "SATURATION");
+    for (int i = 0; i < image->numRows; i++) {
+	for (int j = 0; j < image->numCols; j++) {
+	    if (image->data.F32[i][j] >= SATURATION) {
+		mask->data.U8[i][j] |= PSPHOT_MASK_SATURATED;
+	    }
+	}
+    }
+
+    psLogMsg ("psphot", 3, "load data: %f sec\n", psTimerMark ("psphot"));
+
+    // save the image data & return it
+    psImageData *imdata = psImageDataAlloc(image, noise, mask, header);
+    return (imdata);
+}
Index: /tags/eam_tags/psphot_dev_01/psphot/src/psphotSortBySN.c
===================================================================
--- /tags/eam_tags/psphot_dev_01/psphot/src/psphotSortBySN.c	(revision 22334)
+++ /tags/eam_tags/psphot_dev_01/psphot/src/psphotSortBySN.c	(revision 22334)
@@ -0,0 +1,19 @@
+# include "psphot.h"
+
+// sort by SN (descending)
+int psphotSortBySN (const void **a, const void **b)
+{
+    psSource *A = *(psSource **)a;
+    psSource *B = *(psSource **)b;
+
+    psF32 fA = (A->moments == NULL) ? 0 : A->moments->SN;
+    psF32 fB = (B->moments == NULL) ? 0 : B->moments->SN;
+    if (isnan (fA)) fA = 0;
+    if (isnan (fB)) fB = 0;
+
+    psF32 diff = fA - fB;
+    if (diff > FLT_EPSILON) return (-1);
+    if (diff < FLT_EPSILON) return (+1);
+    return (0);
+}
+
Index: /tags/eam_tags/psphot_dev_01/psphot/src/psphotSourceStats.c
===================================================================
--- /tags/eam_tags/psphot_dev_01/psphot/src/psphotSourceStats.c	(revision 22334)
+++ /tags/eam_tags/psphot_dev_01/psphot/src/psphotSourceStats.c	(revision 22334)
@@ -0,0 +1,54 @@
+# include "psphot.h"
+
+psArray *psphotSourceStats (psImageData *imdata, psMetadata *config, psArray *peaks) 
+{
+    bool     status  = false;
+    psArray *sources = NULL;
+
+    psTimerStart ("psphot");
+
+    // determine properties (sky, moments) of initial sources
+    float INNER    = psMetadataLookupF32 (&status, config, "SKY_INNER_RADIUS");
+    float OUTER    = psMetadataLookupF32 (&status, config, "SKY_OUTER_RADIUS");
+    float RADIUS   = psMetadataLookupF32 (&status, config, "PSF_MOMENTS_RADIUS");
+
+    sources = psArrayAlloc (peaks->n);
+    sources->n = 0;
+
+    for (int i = 0; i < peaks->n; i++) {
+
+	// create a new source, add peak
+	psSource *source = pmSourceAlloc();
+	source->peak = (psPeak *)psMemCopy(peaks->data[i]);
+
+	// allocate image, noise, mask arrays for each peak (square of radius OUTER)
+	// XXX EAM : this has 12 leaks (v.5)
+	pmSourceDefinePixels (source, imdata, source->peak->x, source->peak->y, OUTER);
+
+	// measure a local sky value
+	// XXX EAM : this should use ROBUST not SAMPLE median, but it is broken
+	status = pmSourceLocalSky_EAM (source, PS_STAT_SAMPLE_MEDIAN, INNER);
+	if (!status) {
+	  psFree (source);
+	  continue;
+	}
+
+	// measure basic source moments
+	status = pmSourceMoments (source, RADIUS);
+	if (!status) {
+	  psFree (source);
+	  continue;
+	}
+	
+	// add to the source array
+	psArrayAdd (sources, 100, source);
+	psFree (source);
+    }
+
+    psLogMsg ("psphot", 3, "moments: %f sec\n", psTimerMark ("psphot"));
+    psLogMsg ("psphot", 3, "%d sources\n", sources->n);
+
+    return (sources);
+}
+
+// XXX EAM : filter out bad peaks (eg, no valid pixels available for sky)
Index: /tags/eam_tags/psphot_dev_01/psphot/src/psphotSubtractPSF.c
===================================================================
--- /tags/eam_tags/psphot_dev_01/psphot/src/psphotSubtractPSF.c	(revision 22334)
+++ /tags/eam_tags/psphot_dev_01/psphot/src/psphotSubtractPSF.c	(revision 22334)
@@ -0,0 +1,36 @@
+# include "psphot.h"
+
+// subtract model fits from image (only PSFSTAR && BRIGHTSTAR)
+// need to control application of this with some thresholding
+
+bool psphotSubtractPSF (psSource *source)
+{ 
+  float sky;
+  psModel *model;
+  psImage *pixels;
+
+  // only subtract successful fits
+  if (source->type == PS_SOURCE_SATSTAR) goto valid;
+  if (source->type == PS_SOURCE_PSFSTAR) goto valid;
+  if (source->type == PS_SOURCE_GOODSTAR) goto valid;
+  return (false);
+
+valid:
+  model = source->modelPSF;
+  if (model == NULL) {
+    psLogMsg ("psphot.subtract_psf_source", 2, "missing model for %f, %f\n", source->moments->x, source->moments->y);
+    return (false); 
+  }	    
+
+  pixels = source->pixels;
+
+  // subtract object, leave local sky
+  sky = model->params->data.F32[0];
+  model->params->data.F32[0] = 0;
+  pmSourceSubModel (pixels, source->mask, model, false);
+  model->params->data.F32[0] = sky;
+
+  // XXX EAM : amplify the noise matrix a la dophot?
+
+  return (true);
+}
Index: /tags/eam_tags/psphot_dev_01/psphot/src/pspsf.c
===================================================================
--- /tags/eam_tags/psphot_dev_01/psphot/src/pspsf.c	(revision 22334)
+++ /tags/eam_tags/psphot_dev_01/psphot/src/pspsf.c	(revision 22334)
@@ -0,0 +1,367 @@
+# include "psphot.h"
+
+static void pmPSFFree (pmPSF *psf) {
+
+  if (psf == NULL) return;
+
+  psFree (psf->params);
+  return;
+}
+
+// a PSF always has 4 parameters fewer than the equivalent model
+pmPSF *pmPSFAlloc (psModelType type) {
+
+    int Nparams;
+
+    pmPSF *psf = (pmPSF *) psAlloc(sizeof(pmPSF));
+
+    psf->type = type;
+    psf->chisq = 0.0;
+
+    Nparams = psModelParameterCount (type);
+    if (!Nparams) {
+	psError(PS_ERR_UNKNOWN, true, "Undefined psModelType");
+	return(NULL);
+    }      
+
+    psf->params = psArrayAlloc (Nparams - 4);
+    for (int i = 0; i < psf->params->n; i++) {
+	// we need a way to set the allowed range of orders
+	psf->params->data[i] = Polynomial2DAlloc(1, 1, PS_POLYNOMIAL_ORD);
+    }
+
+    p_psMemSetDeallocator(psf, (psFreeFcn) pmPSFFree);
+    return(psf);
+}
+
+static void pmPSF_TestFree (pmPSF_Test *test) {
+
+  if (test == NULL) return;
+
+  psFree (test->psf);
+  psFree (test->sources);
+  psFree (test->modelFLT);
+  psFree (test->modelPSF);
+  psFree (test->metric);
+  psFree (test->fitMag);
+  psFree (test->mask);
+  return;
+}
+
+pmPSF_Test *pmPSF_TestAlloc (psArray *sources, char *modelName) {
+
+    pmPSF_Test *test = (pmPSF_Test *) psAlloc(sizeof(pmPSF_Test));
+
+    // XXX probably need to increment ref counter
+    test->modelType   = psModelSetType (modelName);
+    test->psf         = pmPSFAlloc (test->modelType);
+    test->sources     = psMemCopy(sources);
+    test->modelFLT    = psArrayAlloc (sources->n);
+    test->modelPSF    = psArrayAlloc (sources->n);
+    test->metric      = psVectorAlloc (sources->n, PS_TYPE_F64);
+    test->fitMag      = psVectorAlloc (sources->n, PS_TYPE_F64);
+    test->mask        = psVectorAlloc (sources->n, PS_TYPE_U8);
+    test->ApResid     = 0;
+    test->dApResid    = 0;
+    test->skyBias     = 0;
+
+    for (int i = 0; i < test->modelFLT->n; i++) {
+	test->mask->data.U8[i]  = 0;
+	test->modelFLT->data[i] = NULL;
+	test->modelPSF->data[i] = NULL;
+	test->metric->data.F64[i] = 0;
+	test->fitMag->data.F64[i] = 0;
+    }	
+
+    p_psMemSetDeallocator(test, (psFreeFcn) pmPSF_TestFree);
+    return (test);
+}
+
+// test->mask values indicate reason source was rejected:
+// 1: outlier in psf polynomial fit
+// 2: flt model failed to converge
+// 3: psf model failed to converge
+// 4: invalid source photometry
+pmPSF_Test *pmPSF_TestModel (psArray *sources, char *modelName, float RADIUS) 
+{
+    bool status;
+    float obsMag;
+    float fitMag;
+    float x;
+    float y;
+    int Nflt = 0;
+    int Npsf = 0;
+
+    pmPSF_Test *test = pmPSF_TestAlloc (sources, modelName);
+
+    // stage 1:  fit an independent model (freeModel) to all sources
+    psTimerStart ("fit");
+    for (int i = 0; i < test->sources->n; i++) {
+
+	psSource *source = test->sources->data[i];
+	psModel  *model  = pmSourceModelGuess (source, test->modelType); 
+	x = source->peak->x;
+	y = source->peak->y;
+
+	// set temporary object mask and fit object
+	// fit model as FLT, not PSF
+	psImageKeepCircle (source->mask, x, y, RADIUS, OR, 0x80);
+	status = pmSourceFitModel (source, model, false);
+	psImageKeepCircle (source->mask, x, y, RADIUS, AND, 0x7f);
+
+	// exclude the poor fits
+	if (!status) {
+	  test->mask->data.U8[i] = 2;
+	  psFree (model);
+	  continue;
+	}
+	test->modelFLT->data[i] = model;
+	Nflt ++;
+    }
+    psLogMsg ("psphot.psftest", 4, "fit flt:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
+    psTrace ("psphot.psftest", 3, "keeping %d of %d PSF candidates (FLT)\n", Nflt, sources->n);
+
+    // make this optional? 
+    // DumpModelFits (test->modelFLT, "modelsFLT.dat");
+
+    // stage 2: construct a psf (pmPSF) from this collection of model fits
+    pmPSFFromModels (test->psf, test->modelFLT, test->mask);
+
+    // stage 3: refit with fixed shape parameters
+    psTimerStart ("fit");
+    for (int i = 0; i < test->sources->n; i++) {
+	// masked for: bad model fit, outlier in parameters
+	if (test->mask->data.U8[i]) continue;
+
+	psSource *source = test->sources->data[i];
+	psModel  *modelFLT = test->modelFLT->data[i];
+
+	// set shape for this model based on PSF
+	psModel *modelPSF = psModelFromPSF (modelFLT, test->psf); 
+	x = source->peak->x;
+	y = source->peak->y;
+
+	psImageKeepCircle (source->mask, x, y, RADIUS, OR, 0x80);
+	status = pmSourceFitModel (source, modelPSF, true);
+
+	// skip poor fits
+	if (!status) {
+	    test->mask->data.U8[i] = 3;
+	    psFree (modelPSF);
+	    goto next_source;
+	}
+
+	// otherwise, save the resulting model
+	test->modelPSF->data[i] = modelPSF;
+
+	// XXX : use a different aperture radius from the fit radius?
+	// XXX : use a different estimator for the local sky?
+	// XXX : pass 'source' as input?
+	if (!pmSourcePhotometry (&fitMag, &obsMag, modelPSF, source->pixels, source->mask)) {
+	    test->mask->data.U8[i] = 4;
+	    goto next_source;
+	}	    
+
+	test->metric->data.F64[i] = obsMag - fitMag;
+	test->fitMag->data.F64[i] = fitMag;
+	Npsf ++;
+
+    next_source:
+	psImageKeepCircle (source->mask, x, y, RADIUS, AND, 0x7f);
+
+    }
+    psLogMsg ("psphot.psftest", 4, "fit psf:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
+    psTrace ("psphot.psftest", 3, "keeping %d of %d PSF candidates (PSF)\n", Npsf, sources->n);
+
+    // make this optional
+    // DumpModelFits (test->modelPSF, "modelsPSF.dat");
+
+    // XXX this function wants aperture radius from pmSourcePhotometry
+    pmPSFMetricModel (test, RADIUS);
+    psLogMsg ("psphot.pspsf", 3, "test model %s, ap-fit: %f +/- %f, sky bias: %f\n", 
+	      modelName, test->ApResid, test->dApResid, test->skyBias);
+
+    return (test);
+}
+
+// input: an array of psModels, pre-allocated psf
+// some of the array entries may be NULL, ignore them
+bool pmPSFFromModels (pmPSF *psf, psArray *models, psVector *mask) {
+
+    // construct the fit vectors from the collection of objects
+    // use the mask to ignore missing fits 
+    psVector *x  = psVectorAlloc (models->n, PS_TYPE_F64);
+    psVector *y  = psVectorAlloc (models->n, PS_TYPE_F64);
+    psVector *z  = psVectorAlloc (models->n, PS_TYPE_F64);
+    psVector *dz = psVectorAlloc (models->n, PS_TYPE_F64);
+
+    for (int i = 0; i < models->n; i++) {
+	psModel *model = models->data[i];
+	if (model == NULL) continue;
+
+	// XXX EAM : this is fragile: x and y must be stored consistently at 2,3
+	x->data.F64[i] = model->params->data.F32[2];
+	y->data.F64[i] = model->params->data.F32[3];
+    }
+
+    // we are doing a robust fit.  after each pass, we drop points which are 
+    // more deviant than three sigma. 
+    // mask is currently updated for each pass. this is inconsistent?
+
+    for (int i = 0; i < psf->params->n; i++) {
+	for (int j = 0; j < models->n; j++) {
+	    psModel *model = models->data[j];
+	    if (model == NULL) continue;
+	    z->data.F64[j] = model->params->data.F32[i + 4];
+	    dz->data.F64[j] = 1;
+	    // XXX EAM : need to use actual errors?
+	    // XXX EAM : this is fragile: psf(Nparams) = model(Nparams) - 4
+	}
+
+	psf->params->data[i] = RobustFit2D (psf->params->data[i], mask, x, y, z, dz);
+	// psTrace ("psphot.psftest", 3, "keeping %d of %d PSF candidates (PSF param %d)\n", Nkeep, mask->n, i);
+	// psPolynomial2DDump (psf->params->data[i]);
+    }
+    
+    psFree (x);
+    psFree (y);
+    psFree (z);
+    psFree (dz);
+    return (true);
+}
+
+psModel *psModelFromPSF (psModel *modelFLT, pmPSF *psf) {
+
+    // need to define the relationship between the modelFLT and modelPSF ?
+    
+    // find function used to set the model parameters
+    psModelFromPSFFunc modelFromPSFFunc = psModelFromPSFFunc_GetFunction (psf->type);
+
+    // do we need a different model for the PSF vs FLT version?
+    psModel *modelPSF = psModelAlloc (psf->type);
+
+    // set model parameters for this source based on PSF information
+    modelFromPSFFunc (modelPSF, modelFLT, psf);
+
+    return (modelPSF);
+}
+
+bool pmPSFMetricModel (pmPSF_Test *test, float RADIUS) {
+
+  float dBin;
+  int   nKeep, nSkip;
+
+  // the measured (aperture - fit) magnitudes (dA == test->metric)
+  //   depend on both the true ap-fit (dAo) and the bias in the sky measurement:
+  //     dA = dAo + dsky/flux
+  //   where flux is the flux of the star
+  // we fit this trend to find the infinite flux aperture correction (dAo),
+  //   the nominal sky bias (dsky), and the error on dAo
+  // the values of dA are contaminated by stars with close neighbors in the aperture
+  //   we use an outlier rejection to avoid this bias
+
+  // rflux = ten(0.4*fitMag);
+  psVector *rflux = psVectorAlloc (test->sources->n, PS_TYPE_F64);
+  for (int i = 0; i < test->sources->n; i++) {
+    if (test->mask->data.U8[i]) continue;
+    rflux->data.F64[i] = pow(10.0, 0.4*test->fitMag->data.F64[i]);
+  }
+
+  // find min and max of (1/flux):
+  psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
+  psVectorStats (stats, rflux, NULL, test->mask, 0xff);
+  
+  // build binned versions of rflux, metric
+  dBin = (stats->max - stats->min) / 10.0;
+  psVector *rfBin = psVectorCreate (stats->min, stats->max, dBin, PS_TYPE_F64);
+  psVector *daBin = psVectorAlloc (rfBin->n, PS_TYPE_F64);
+  psVector *maskB = psVectorAlloc (rfBin->n, PS_TYPE_U8);
+  psFree (stats);
+
+  psTrace ("psphot.metricmodel", 3, "rflux max: %g, min: %g, delta: %g\n", stats->max, stats->min, dBin);
+
+  // group data in daBin bins, measure lower 50% mean
+  for (int i = 0; i < daBin->n; i++) {
+
+    psVector *tmp = psVectorAlloc (test->sources->n, PS_TYPE_F64);
+    tmp->n = 0;
+
+    // accumulate data within bin range
+    for (int j = 0; j < test->sources->n; j++) {
+      // masked for: bad model fit, outlier in parameters
+      if (test->mask->data.U8[j]) continue;
+    
+      // skip points with extreme dA values
+      if (fabs(test->metric->data.F64[j]) > 0.5) continue;
+
+      // skip points outside of this bin
+      if (rflux->data.F64[j] < rfBin->data.F64[i] - 0.5*dBin) continue;
+      if (rflux->data.F64[j] > rfBin->data.F64[i] + 0.5*dBin) continue;
+
+      tmp->data.F64[tmp->n] = test->metric->data.F64[j];
+      tmp->n ++;
+    }
+
+    // is this a valid point?
+    maskB->data.U8[i] = 0;
+    if (tmp->n < 2) {
+      maskB->data.U8[i] = 1;
+      psFree (tmp);
+      continue;
+    } 
+
+    // dA values are contaminated with low outliers 
+    // measure statistics only on upper 50% of points
+    // this would be easier if we could sort in reverse:
+    //
+    // psVectorSort (tmp, tmp);
+    // tmp->n = 0.5*tmp->n;
+    // stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
+    // psVectorStats (stats, tmp, NULL, NULL, 0);
+    // psTrace ("psphot.metricmodel", 4, "rfBin %d (%g): %d pts, %g\n", i, rfBin->data.F64[i], tmp->n, stats->sampleMedian);
+
+    psVectorSort (tmp, tmp);
+    nKeep = 0.5*tmp->n;
+    nSkip = tmp->n - nKeep;
+
+    psVector *tmp2 = psVectorAlloc (nKeep, PS_TYPE_F64);
+    for (int j = 0; j < tmp2->n; j++) {
+      tmp2->data.F64[j] = tmp->data.F64[j + nSkip];
+    }
+
+    stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
+    psVectorStats (stats, tmp2, NULL, NULL, 0);
+    psTrace ("psphot.metricmodel", 4, "rfBin %d (%g): %d pts, %g\n", i, rfBin->data.F64[i], tmp->n, stats->sampleMedian);
+
+    daBin->data.F64[i] = stats->sampleMedian;
+
+    psFree (stats);
+    psFree (tmp);
+    psFree (tmp2);
+  }
+
+  // linear fit to rfBin, daBin
+  psPolynomial1D *poly = Polynomial1DAlloc (1, PS_POLYNOMIAL_ORD);
+  poly = VectorFitPolynomial1DOrd_EAM (poly, maskB, rfBin, daBin, NULL);
+
+  psVector *daBinFit = Polynomial1DEvalVector_EAM (poly, rfBin);
+  psVector *daResid  = (psVector *) psBinaryOp (NULL, (void *) daBin, "-", (void *) daBinFit);
+
+  stats = psStatsAlloc (PS_STAT_CLIPPED_STDEV);
+  stats = psVectorStats (stats, daResid, NULL, maskB, 1);
+
+  test->ApResid = poly->coeff[0];
+  test->dApResid = stats->clippedStdev;
+  test->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
+
+  psFree (rflux);
+  psFree (rfBin);
+  psFree (daBin);
+  psFree (maskB);
+  psFree (daBinFit);
+  psFree (daResid);
+  psFree (poly);
+  psFree (stats);
+
+  return true;
+}