Index: /trunk/psphot/Makefile
===================================================================
--- /trunk/psphot/Makefile	(revision 4114)
+++ /trunk/psphot/Makefile	(revision 4114)
@@ -0,0 +1,80 @@
+default: psphot
+help:
+	@echo "USAGE: make psphot"
+
+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)/pspsf.$(ARCH).o \
+$(SRC)/psphot-utils.$(ARCH).o \
+$(SRC)/psPolynomials.$(ARCH).o \
+$(SRC)/psUtils.$(ARCH).o \
+$(SRC)/setup.$(ARCH).o \
+$(SRC)/image_stats.$(ARCH).o \
+$(SRC)/find_peaks.$(ARCH).o \
+$(SRC)/source_moments.$(ARCH).o \
+$(SRC)/choose_psf_model.$(ARCH).o \
+$(SRC)/apply_psf_model.$(ARCH).o \
+$(SRC)/test_psf_scatter.$(ARCH).o \
+$(SRC)/mark_psf_sources.$(ARCH).o \
+$(SRC)/subtract_psf_sources.$(ARCH).o \
+$(SRC)/fit_galaxies.$(ARCH).o \
+$(SRC)/subtract_galaxies.$(ARCH).o
+
+psphot: $(BIN)/psphot.$(ARCH)
+
+$(BIN)/psphot.$(ARCH) : $(PSPHOT)
+
+$(PSPHOT): $(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: /trunk/psphot/doc/psphot.txt
===================================================================
--- /trunk/psphot/doc/psphot.txt	(revision 4114)
+++ /trunk/psphot/doc/psphot.txt	(revision 4114)
@@ -0,0 +1,144 @@
+
+2005.04.12 : Notes on psphot
+
+- 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: /trunk/psphot/src/apply_psf_model.c
===================================================================
--- /trunk/psphot/src/apply_psf_model.c	(revision 4114)
+++ /trunk/psphot/src/apply_psf_model.c	(revision 4114)
@@ -0,0 +1,65 @@
+# include "psphot.h"
+
+// fit psf model to all objects not already fitted (modelPSF == NULL)
+// run this function to a specific flux limit?
+
+bool apply_psf_model (psImage *image, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) 
+{ 
+    bool  status;
+    int   Nfit = 0;
+    int   Niter = 0;
+
+    psTimerStart ("psphot");
+
+    // we may set this differently here from the value used to mark likely saturated stars
+    float SATURATE    = psMetadataLookupF32 (&status, config, "SATURATE");
+    float FIT_NSIGMA  = psMetadataLookupF32 (&status, config, "FIT_NSIGMA");
+    float FIT_PADDING = psMetadataLookupF32 (&status, config, "FIT_PADDING");
+
+    float XBORDER  = psMetadataLookupF32 (&status, config, "XBORDER");
+    float YBORDER  = psMetadataLookupF32 (&status, config, "YBORDER");
+    psRegion *keep = psRegionAlloc (XBORDER, -XBORDER, YBORDER, -YBORDER);
+    keep           = psRegionForImage (keep, image, keep);
+
+    // set the object surfact-brightness limit for fitted pixels
+    float FLUX_LIMIT  = 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];
+	if (source->modelPSF != NULL) continue;
+
+	// skip non-astronomical objects (very likely defects)
+	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  *model  = pmSourceModelGuess (source, psf->type); 
+
+	// set PSF parameters for this model
+	model  = psModelFromPSF (model, psf);
+
+	// set the fit radius based on the object flux limit and the model
+	// FLUX_LIMIT should be set based on local sky model (not global median)
+	model->radius = modelRadius (model->params, FLUX_LIMIT) + FIT_PADDING;
+	if (isnan(model->radius)) {
+	  fprintf (stderr, "error in radius\n");
+	}
+	    
+	// set the valid/invalid pixel mask
+	pmSourceSetPixelCircle (source, image, model->radius);
+	pmSourceMaskRegion (source, keep);
+	pmSourceMaskSaturated (source, SATURATE);
+
+	// fit as PSF, not FLT (drop poor fits)
+	if (!pmSourceFitModel (source, model, true)) continue;
+	source->modelPSF = model;
+	Niter += model[0].nIter;
+	Nfit ++;
+    }
+    psLogMsg ("psphot", 3, "fit PSF models: %f sec for %d objects (%d total iterations)\n", psTimerMark ("psphot"), Nfit, Niter);
+    return (true);
+}
Index: /trunk/psphot/src/choose_psf_model.c
===================================================================
--- /trunk/psphot/src/choose_psf_model.c	(revision 4114)
+++ /trunk/psphot/src/choose_psf_model.c	(revision 4114)
@@ -0,0 +1,73 @@
+# include "psphot.h"
+
+// test PSF models and select best option
+
+pmPSF *choose_psf_model (psImage *image, psMetadata *config, psArray *sources) 
+{ 
+    bool        status;
+    char       *modelName;
+    char        key[64];
+    pmPSF_Test *test  = NULL;
+    psArray    *stars = NULL;
+
+    // select the possible PSF stars (pointers to original sources)
+    stars = psArrayAlloc (sources->n);
+    stars->n = 0;
+
+    for (int i = 0; i < sources->n; i++) {
+	psSource *source = sources->data[i];
+	if (source->type != PS_SOURCE_PSFSTAR) continue;
+	psArrayAdd (stars, 200, source);
+    }
+
+    // define model fit pixels
+    // XXX use a different radius - PSF_FIT_RADIUS?
+    float RADIUS   = psMetadataLookupF32 (&status, config, "FIT_RADIUS");
+
+    float XBORDER  = psMetadataLookupF32 (&status, config, "XBORDER");
+    float YBORDER  = psMetadataLookupF32 (&status, config, "YBORDER");
+    psRegion *keep = psRegionAlloc (XBORDER, -XBORDER, YBORDER, -YBORDER);
+    keep           = psRegionForImage (keep, image, keep);
+
+    for (int i = 0; i < stars->n; i++) {
+	psSource *source = stars->data[i];
+	pmSourceSetPixelCircle (source, image, RADIUS);
+	pmSourceMaskRegion (source, keep);
+	// XXX EAM - this is silly: reconsider how we use the masks
+	// XXX - we are needlessly re-defining the pixels here
+    }
+
+    // the user specifies 1 or more PSF model options
+    // run test on each model option
+    int      Ntest = psMetadataLookupS32 (&status, config, "PSF_MODEL_N");
+    psArray *tests = psArrayAlloc (Ntest);
+    for (int i = 0; i < tests->n; i++) { 
+	sprintf (key, "PSF_MODEL_%d", i);
+	modelName = psMetadataLookupPtr (&status, config, key);
+	tests->data[i] = pmPSF_TestModel (stars, modelName);
+    }
+
+    // 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->metricStats->clippedStdev;
+    for (int i = 1; i < Ntest; i++) {
+	test = tests->data[i];
+	int M = test->metricStats->clippedStdev;
+	if (M < bestM) {
+	    bestM = M;
+	    bestN = i;
+	}
+    }
+    test = tests->data[bestN];
+    modelName = psModelGetType (test->modelType);
+    psLogMsg ("psphot.pspsf", 3, "selected psf model %s, metric: %f +/- %f\n", modelName, test->metricStats->clippedMean, test->metricStats->clippedStdev);
+
+    // set source->model based on best psf model fit
+    for (int i = 0; i < test->sources->n; i++) {
+	psSource *source = test->sources->data[i];
+	source->modelPSF = test->modelPSF->data[i];
+    }
+    return (test->psf);
+}
Index: /trunk/psphot/src/find_peaks.c
===================================================================
--- /trunk/psphot/src/find_peaks.c	(revision 4114)
+++ /trunk/psphot/src/find_peaks.c	(revision 4114)
@@ -0,0 +1,46 @@
+# include "psphot.h"
+
+psArray *find_peaks (psImage *image, psMetadata *config, psStats *sky) 
+{
+    bool  status = false;
+    float NSIGMA;
+    float SIGMA;
+    float threshold;
+
+    // smooth the image 
+    // should we also subtract a super-binned image? (as an option)
+
+    psTimerStart ("psphot");
+
+    SIGMA  = psMetadataLookupF32 (&status, config, "SMOOTH_SIGMA");
+    NSIGMA = psMetadataLookupF32 (&status, config, "SMOOTH_NSIGMA");
+
+    psImage *smooth = psImageCopy (NULL, image, PS_TYPE_F32);
+    psImageSmooth (smooth, SIGMA, NSIGMA);
+
+    // XXX missing pslib functions:
+    // kernel = psKernelParts ();
+    // smooth = psImageConvolve (smooth, image, kernel, PS_PARTS);
+
+    psLogMsg ("psphot", 4, "smooth: %f sec\n", psTimerMark ("psphot"));
+
+
+    // find the peaks in the smoothed image 
+
+    psTimerStart ("psphot");
+
+    NSIGMA    = psMetadataLookupF32 (&status, config, "PEAK_NSIGMA");
+    threshold = NSIGMA*sky->sampleStdev + sky->sampleMean;
+    psLogMsg ("psphot", 3, "threshold: %f DN\n", threshold);
+
+    psArray *peaks = pmFindImagePeaks (smooth, threshold);
+    if (peaks == NULL) psAbort ("find peaks", "no peaks found");
+    // XXX Need to gracefully handle no source detections
+
+    // make this optional
+    DumpPeaks (peaks, "peaks.dat");
+
+    psLogMsg ("psphot", 3, "%d peaks: %f sec\n", peaks->n, psTimerMark ("psphot"));
+
+    return (peaks);
+}
Index: /trunk/psphot/src/fit_galaxies.c
===================================================================
--- /trunk/psphot/src/fit_galaxies.c	(revision 4114)
+++ /trunk/psphot/src/fit_galaxies.c	(revision 4114)
@@ -0,0 +1,40 @@
+# include "psphot.h"
+
+// fit selected galaxy model (GAUSS) to all bright objects of type GALAXY
+
+bool fit_galaxies (psImage *image, psMetadata *config, psArray *sources) 
+{ 
+    bool  status;
+    int   Nfit = 0;
+    int   Niter = 0;
+
+    float XBORDER  = psMetadataLookupF32 (&status, config, "XBORDER");
+    float YBORDER  = psMetadataLookupF32 (&status, config, "YBORDER");
+    psRegion *keep = psRegionAlloc (XBORDER, -XBORDER, YBORDER, -YBORDER);
+    keep           = psRegionForImage (keep, image, keep);
+
+    float RADIUS   = psMetadataLookupF32 (&status, config, "FIT_RADIUS");
+    float snFaint  = psMetadataLookupF32 (&status, config, "FAINT_SN_LIM");
+    psModelType modelType = psModelSetType ("PS_MODEL_PGAUSS");
+
+    psTimerStart ("psphot");
+    for (int i = 0; i < sources->n; i++) {
+	psSource *source = sources->data[i];
+	if (source->type != PS_SOURCE_GALAXY) continue;
+	if (source->moments->SN < snFaint) continue;
+
+	// need a better model guess and a better radius choice
+	pmSourceSetPixelCircle (source, image, RADIUS);
+	pmSourceMaskRegion (source, keep);
+	psModel  *model  = pmSourceModelGuess (source, modelType); 
+
+	// fit as FLT, not PSF (skip poor fits)
+	if (!pmSourceFitModel (source, model, false)) continue;
+	source->modelFLT = model;
+
+	Niter += model[0].nIter;
+	Nfit++;
+    }
+    psLogMsg ("psphot", 3, "fit galaxies: %f sec for %d galaxies (%d total iterations)\n", psTimerMark ("psphot"), Nfit, Niter);
+    return (true);
+}
Index: /trunk/psphot/src/image_stats.c
===================================================================
--- /trunk/psphot/src/image_stats.c	(revision 4114)
+++ /trunk/psphot/src/image_stats.c	(revision 4114)
@@ -0,0 +1,52 @@
+# include "psphot.h"
+
+psStats *image_stats (psImage *image, psMetadata *config) 
+{ 
+    psStats *stats = NULL;
+    psStats *sky   = NULL;
+
+    // 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;
+	psRandom *rnd    = psRandomAlloc (PS_RANDOM_TAUS, 0);
+	int    Nsubset   = psMetadataLookupF32 (&status, config, "NSUBSET");
+	double fSubset   = (double) Nsubset / (image->numRows*image->numCols);
+	psVector *subset = psVectorAlloc (Nsubset, PS_TYPE_F32);
+
+	subset->n = 0;
+	for (int i = 0; i < image->numRows; i++) {
+	    for (int j = 0; j < image->numCols; j++) {
+		double frnd = psRandomUniform (rnd);
+		if ((fSubset < 1.0) && (fSubset < frnd)) continue;
+		subset->data.F32[subset->n] = image->data.F32[j][i];
+		subset->n ++;
+		if (subset->n == Nsubset) goto got_subset;
+	    }
+	}
+
+    got_subset:
+	// 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?
+    }
+
+    // we store these values in mean,stdev 
+    bool status = false;
+    float NOISE = psMetadataLookupF32 (&status, config, "RDNOISE");
+    float GAIN  = psMetadataLookupF32 (&status, config, "GAIN");
+
+    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);
+
+    return (sky);
+}
Index: /trunk/psphot/src/mark_psf_sources.c
===================================================================
--- /trunk/psphot/src/mark_psf_sources.c	(revision 4114)
+++ /trunk/psphot/src/mark_psf_sources.c	(revision 4114)
@@ -0,0 +1,57 @@
+# 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)
+
+// any objects which is consistent with the PSF should have 
+// abs(dparams[5]) < N * dsxLine(mag) & abs(dparams[6]) < N * dsyLine(mag)
+// this should include a minimum dsx buffer for the brighter objects
+// saturated stars should fall outside (but are already IDed)
+// galaxies should be larger, cosmic rays smaller, but need to test?
+
+// need to demote PSFSTARs which don't meet criterion
+
+bool mark_psf_sources (psArray *sources, psMetadata *config)
+{ 
+    bool  status      = false;
+    float shapeNsigma = psMetadataLookupF32 (&status, config, "PSF_SHAPE_NSIGMA");
+    float snFaint     = psMetadataLookupF32 (&status, config, "FAINT_SN_LIM");
+    float dSX, dSY, SX, SY, SN;
+    float nSx, nSy;
+
+    for (int i = 0; i < sources->n; i++) {
+	psSource *source = sources->data[i];
+	if (source->type == PS_SOURCE_PSFSTAR) continue;
+	if (source->type == PS_SOURCE_SATSTAR) continue;
+	if (source->modelPSF == NULL) continue;
+
+	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];
+
+	nSx = dSX * SX * SN;
+	nSy = dSY * SY * SN;
+
+	// assign PS_SOURCE_BRIGHTSTAR to bright objects within PSF region of dparams[]
+	if ((fabs(nSx) < shapeNsigma) && (fabs(nSy) < shapeNsigma)) {
+	    if (SN > snFaint) {
+		source->type = PS_SOURCE_BRIGHTSTAR;
+	    } else {
+		source->type = PS_SOURCE_FAINTSTAR;
+	    }
+	}
+	    
+	if ((nSx >= shapeNsigma) || (nSy >= shapeNsigma)) {
+	    source->type = PS_SOURCE_GALAXY;
+	}
+	if ((nSx <= -shapeNsigma) || (nSy <= -shapeNsigma)) {
+	    source->type = PS_SOURCE_DEFECT;
+	}
+	    
+    }
+    return (true);
+}	
Index: /trunk/psphot/src/pixregion.c
===================================================================
--- /trunk/psphot/src/pixregion.c	(revision 4114)
+++ /trunk/psphot/src/pixregion.c	(revision 4114)
@@ -0,0 +1,14 @@
+# include "psphot.h"
+
+psRegion *define_valid_pixels (psImage *image, psMetadata *config)
+{
+
+    bool  status;
+    float XBORDER  = psMetadataLookupF32 (&status, config, "XBORDER");
+    float YBORDER  = psMetadataLookupF32 (&status, config, "YBORDER");
+
+    psRegion *pixregion = psRegionAlloc (XBORDER, -XBORDER, YBORDER, -YBORDER);
+    pixregion = psRegionForImage (pixregion, image, pixregion);
+
+    return (pixregion);
+}
Index: /trunk/psphot/src/psPolynomials.c
===================================================================
--- /trunk/psphot/src/psPolynomials.c	(revision 4114)
+++ /trunk/psphot/src/psPolynomials.c	(revision 4114)
@@ -0,0 +1,473 @@
+# include "psphot.h"
+
+// 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);
+}
+
+// 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));
+    // p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial1DFree);
+    // XXX EAM : me, being lazy
+
+    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;
+    }
+    return(newPoly);
+}
+
+// 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));
+    // p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial2DFree);
+    // XXX EAM : me, being lazy
+
+    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;
+        }
+    }
+    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 : 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 1d fitting
+psPolynomial1D* VectorFitPolynomial1DOrd_EAM(psPolynomial1D* myPoly,
+        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++) {
+	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(A);
+    psFree(B);
+    psFree(xSums);
+
+    psTrace(".psLib.dataManip.VectorFitPolynomial1DOrd", 4,
+            "---- VectorFitPolynomial1DOrd() begin ----\n");
+    return (myPoly);
+}
+
+// XXX EAM : test version of 2d fitting
+psPolynomial2D* VectorFitPolynomial2DOrd_EAM(psPolynomial2D* myPoly,
+        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++) {
+	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]);
+	}
+    }
+}    
+
+// 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]);
+    }
+}    
+
+psPolynomial2D* RobustFit2D(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, 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 : 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: /trunk/psphot/src/psUtils.c
===================================================================
--- /trunk/psphot/src/psUtils.c	(revision 4114)
+++ /trunk/psphot/src/psUtils.c	(revision 4114)
@@ -0,0 +1,262 @@
+# include "psphot.h"
+
+static psHash *timers = NULL;
+
+// set actual region based on image parameters
+// XXX EAM : this needs to be changes to use psRegion rather than psRegion*
+psRegion *psRegionForImage (psRegion *out, psImage *image, psRegion *in) {
+    
+    if (out == NULL) {
+	out = psRegionAlloc(in->x0, in->x1, in->y0, in->y1);
+    } else {
+	*out = *in;
+    }
+    if (out->x1 <= 0) {
+	out->x1 = image->numCols + out->x1;
+    }
+    if (out->y1 <= 0) {
+	out->y1 = image->numRows + out->y1;
+    }
+    return (out);
+}
+
+
+// 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);
+
+    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);
+    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 = get_argument (*argc, argv, "-v"))) {
+	remove_argument (N, argc, argv);
+	level = 3;
+    }
+    if ((N = get_argument (*argc, argv, "-vv"))) {
+	remove_argument (N, argc, argv);
+	level = 4;
+    }
+    if ((N = get_argument (*argc, argv, "-vvv"))) {
+	remove_argument (N, argc, argv);
+	level = 5;
+    }
+
+    if ((N = get_argument (*argc, argv, "-logfmt"))) {
+	if (*argc < N + 2) {
+	    psAbort ("psLogArguments", "USAGE: -logfmt (format)");
+	}
+	remove_argument (N, argc, argv);
+	psLogSetFormat (argv[N]); // XXX EAM : this function should return an error if the log format is invalid
+	remove_argument (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 = get_argument (*argc, argv, "-trace"))) {
+	if (*argc < N + 3) {
+	    psAbort ("psTraceArguments", "USAGE: -trace (facility) (level)");
+	}
+	remove_argument (N, argc, argv);
+	psTraceSetLevel (argv[N], atoi(argv[N+1]));
+	remove_argument (N, argc, argv);
+	remove_argument (N, argc, argv);
+    }
+    if ((N = get_argument (*argc, argv, "-trace-levels"))) {
+	psTracePrintLevels ();
+	exit (2);
+    }
+    return (TRUE);
+}
+
+// find the location of the specified argument
+int get_argument (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 remove_argument (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;
+	    return (value);
+	}	
+    }
+
+    //  ... or a value (F32?)
+    if (item->type == PS_META_F32) {
+	value = item->data.F32;
+	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);
+}
+    
Index: /trunk/psphot/src/psphot-utils.c
===================================================================
--- /trunk/psphot/src/psphot-utils.c	(revision 4114)
+++ /trunk/psphot/src/psphot-utils.c	(revision 4114)
@@ -0,0 +1,181 @@
+# include "psphot.h"
+
+// dump the peaks to an output file
+bool DumpPeaks (psArray *sources, char *filename) 
+{
+
+    int i;
+    FILE *f;
+    psPeak *peak;
+	
+    f = fopen (filename, "w");
+    if (f == NULL) {
+	psLogMsg ("DumpPeaks", 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;
+}
+
+// dump the moments to an output file 
+bool DumpMoments (psArray *sources, char *filename) 
+{
+
+    int i;
+    FILE *f;
+    psSource *source;
+	
+    f = fopen (filename, "w");
+    if (f == NULL) {
+	psLogMsg ("DumpMoments", 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 DumpModelPSF (psArray *sources, char *filename) 
+{
+
+    int i, j;
+    FILE *f;
+    psVector *params;
+    psVector *dparams;
+    psModel  *model;
+
+    f = fopen (filename, "w");
+    if (f == NULL) {
+	psLogMsg ("DumpObjects", 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;
+	params = model->params;
+	dparams = model->dparams;
+	fprintf (f, "%7.1f %7.1f  %5.1f %7.1f  %7.4f  ", 
+		 params[0].data.F32[2], params[0].data.F32[3], 
+		 params[0].data.F32[0], params[0].data.F32[1], (dparams[0].data.F32[1]/params[0].data.F32[1]));
+	for (j = 0; j < model->params->n - 4; j++) {
+	    fprintf (f, "%9.3g ", params[0].data.F32[j+4]);
+	}
+	for (j = 0; j < model->params->n - 4; j++) {
+	    fprintf (f, "%9.3g ", dparams[0].data.F32[j+4]);
+	}
+	fprintf (f, ": %2d %7.4f %7.4f %7.2f (%d %d)\n", source[0].type, log10(model[0].chisq), source[0].moments->SN, model[0].radius, model[0].nDOF, model[0].nIter);
+    }
+    fclose (f);
+    return true;
+}
+
+// dump the models to an output file
+bool DumpModelFits (psArray *models, char *filename) 
+{
+
+    int i, j;
+    FILE *f;
+    psVector *params;
+    psVector *dparams;
+    psModel  *model;
+
+    f = fopen (filename, "w");
+    if (f == NULL) {
+	psLogMsg ("DumpObjects", 3, "can't open output file for moments%s\n", filename);
+	return false;
+    }
+
+    // this loop is valid for GAUSS, PGAUSS:
+    for (i = 0; i < models->n; i++) {
+	model = (psModel *) models->data[i];
+	if (model == NULL) continue;
+	params = model->params;
+	dparams = model->dparams;
+	fprintf (f, "%7.1f %7.1f  %5.1f %7.1f  ", 
+		 params[0].data.F32[2], params[0].data.F32[3], 
+		 params[0].data.F32[0], params[0].data.F32[1]);
+	for (j = 0; j < model->params->n - 4; j++) {
+	    fprintf (f, "%9.3g ", params[0].data.F32[j+4]);
+	}
+	for (j = 0; j < model->params->n - 4; j++) {
+	    fprintf (f, "%9.3g ", dparams[0].data.F32[j+4]);
+	}
+	fprintf (f, ":  %7.1f (%d %d)\n", model[0].chisq, model[0].nDOF, model[0].nIter);
+    }
+    fclose (f);
+    return true;
+}
+
+void pmSourceMaskSaturated (psSource *source, float saturate) {
+
+    psImage *image = source->pixels;
+    psImage *mask  = source->mask;
+
+    for (int iy = 0; iy < image->numRows; iy++) {
+	for (int ix = 0; ix < image->numCols; ix++) {
+	    if (image->data.F32[iy][ix] < saturate) continue;
+	    mask->data.U8[iy][ix] |= 0x02;
+	}
+    }
+}
+
+void pmSourceMaskRegion (psSource *source, psRegion *region) {
+
+    psImage *image = source->pixels;
+    psImage *mask  = source->mask;
+
+    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:
+	    mask->data.U8[iy][ix] |= 0x02;
+	}
+    }
+}
+
+// create a subset of sources for the PS_SOURCE_PSFSTAR entries
+// XXX deprecated
+psArray *SelectPSFStars (psArray *sources) {
+
+    psArray  *stars;
+    psSource *source;
+    int i;
+
+    stars = psArrayAlloc (sources->n);
+    stars->n = 0;
+
+    for (i = 0; i < sources->n; i++) {
+	source = (psSource *) sources->data[i];
+	if (source->type != PS_SOURCE_PSFSTAR) continue;
+	psArrayAdd (stars, 200, source);
+    }
+
+    return (stars);
+}
+
Index: /trunk/psphot/src/psphot.c
===================================================================
--- /trunk/psphot/src/psphot.c	(revision 4114)
+++ /trunk/psphot/src/psphot.c	(revision 4114)
@@ -0,0 +1,64 @@
+# include "psphot.h"
+
+int main (int argc, char **argv) {
+
+    psMetadata *config  = NULL;		// user-provided configuration information
+    psMetadata *header  = NULL;		// input image header
+    psImage    *image   = NULL;		// input image data
+    psArray    *sources = NULL;
+    psArray    *peaks   = NULL;		// a list of pmPeaks 
+    pmPSF      *psf     = NULL;
+    psStats    *sky     = NULL;
+
+    psLogArguments (&argc, argv);
+    psTraceArguments (&argc, argv);
+    if (argc != 4) usage ();
+
+    // load input data (image and config)
+    setup (&image, &header, &config, argv);
+
+    // measure image stats for initial guess 
+    sky = image_stats (image, config);
+
+    // find the peaks in the image
+    peaks = find_peaks (image, config, sky);
+
+    sources = source_moments (image, config, peaks);
+
+    // use stellar objects SN > PSF_SN_LIM
+    psf = choose_psf_model (image, config, sources);
+
+    // test PSF on all except SATURATE and DEFECT (artifacts)
+    apply_psf_model (image, config, sources, psf, sky);
+
+    // XXX this test compares dSx,dSy distribution with model.
+    // XXX it has some problems
+    // test_psf_scatter (sources);
+
+    // identify PSF-like objects (< PSF_SHAPE_NSIGMA from dSx,dSy model)
+    // SN <= FAINT_SN_LIM : FAINTSTAR
+    // SN >  FAINT_SN_LIM : BRIGHTSTAR
+    mark_psf_sources (sources, config);
+
+    // subtract BRIGHTSTAR, PSFSTAR, SATSTAR
+    subtract_psf_sources (sources);
+
+    // fit_galaxies (image, config, sources);
+
+    // subtract_galaxies (sources, config);
+ 
+    // write out subtracted image
+    unlink (argv[2]);
+    psFits *fits = psFitsAlloc (argv[2]);
+    psFitsWriteImage (fits, NULL, image, 0, NULL);
+    psFree (fits);
+
+    DumpModelPSF (sources, "sources.dat");
+    exit (0);
+}
+
+int usage () {
+
+    fprintf (stderr, "USAGE: psphot (image.fits) (output.fits) (config)\n");
+    exit (2);
+}
Index: /trunk/psphot/src/psphot.h
===================================================================
--- /trunk/psphot/src/psphot.h	(revision 4114)
+++ /trunk/psphot/src/psphot.h	(revision 4114)
@@ -0,0 +1,75 @@
+# include <pslib.h>
+# include <pmObjects.h>
+# include <string.h>
+# include <strings.h>
+# include <unistd.h>
+# include <stdlib.h>
+
+// 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   *metric;
+  psStats    *metricStats;
+} pmPSF_Test;
+
+bool psTimerStart (char *name);
+psF64 psTimerMark (char *name);
+
+void psImageSmooth (psImage *image, float sigma, float Nsigma);
+psRegion *psRegionForImage (psRegion *out, psImage *image, psRegion *in);
+
+bool DumpPeaks (psArray *sources, char *filename);
+bool DumpMoments (psArray *sources, char *filename);
+bool DumpModelFits (psArray *models, char *filename);
+bool DumpModelPSF (psArray *sources, char *filename);
+
+int usage ();
+int get_argument (int argc, char **argv, char *arg);
+int remove_argument (int N, int *argc, char **argv);
+
+psArray *pmPeaksSubset(psArray *peaks, psF32 maxValue, const psRegion *valid);
+psArray *SelectPSFStars (psArray *sources);
+
+psF32 pmConfigLookupF32 (bool *status, psMetadata *config, psMetadata *header, char *name);
+
+psF32 Polynomial2DEval(const psPolynomial2D* myPoly, psF32 x, psF32 y);
+psImage *psBuildSums2D(psImage* sums,psF64 x,psF64 y,psS32 nXterm, psS32 nYterm);
+psPolynomial1D* VectorFitPolynomial1DOrd_EAM(psPolynomial1D* myPoly, const psVector* x, const psVector* y, const psVector* yErr);
+psPolynomial2D* VectorFitPolynomial2DOrd_EAM(psPolynomial2D* myPoly, 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(psPolynomial2D* poly, 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);
+
+psS32 psLogArguments (int *argc, char **argv);
+psS32 psTraceArguments (int *argc, char **argv);
+psVector *p_psGetRowVectorFromImage(psImage *image, psU32 row);
+psVector *psGetRowVectorFromImage(psImage *image, psU32 row);
+
+pmPSF *pmPSFAlloc (psModelType type);
+pmPSF_Test *pmPSF_TestAlloc (psArray *stars, char *modelName);
+pmPSF_Test *pmPSF_TestModel (psArray *stars, char *modelName);
+
+bool pmPSFFromModels (pmPSF *psf, psArray *models);
+psModel *psModelFromPSF (psModel *model, pmPSF *psf);
+void pmSourceMaskRegion (psSource *source, psRegion *region);
+void pmSourceMaskSaturated (psSource *source, float saturate);
+
+// top-level psphot functions
+bool setup (psImage **image, psMetadata **header, psMetadata **config, char **argv);
+psStats *image_stats (psImage *image, psMetadata *config);
+psArray *find_peaks (psImage *image, psMetadata *config, psStats *sky);
+psArray *source_moments (psImage *image, psMetadata *config, psArray *allpeaks);
+pmPSF *choose_psf_model (psImage *image, psMetadata *config, psArray *sources);
+bool apply_psf_model (psImage *image, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
+bool test_psf_scatter (psArray *sources);
+bool mark_psf_sources (psArray *sources, psMetadata *config);
+bool subtract_psf_sources (psArray *sources);
+bool fit_galaxies (psImage *image, psMetadata *config, psArray *sources);
+bool subtract_galaxies (psArray *sources, psMetadata *config);
Index: /trunk/psphot/src/pspsf.c
===================================================================
--- /trunk/psphot/src/pspsf.c	(revision 4114)
+++ /trunk/psphot/src/pspsf.c	(revision 4114)
@@ -0,0 +1,208 @@
+# include "psphot.h"
+
+// 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);
+    }
+    return(psf);
+}
+
+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     = sources;
+    test->modelFLT    = psArrayAlloc (sources->n);
+    test->modelPSF    = psArrayAlloc (sources->n);
+    test->metricStats = NULL;
+    test->metric      = psVectorAlloc (sources->n, PS_TYPE_F64);
+
+    for (int i = 0; i < test->modelFLT->n; i++) {
+	test->modelFLT->data[i] = NULL;
+	test->modelPSF->data[i] = NULL;
+    }	
+    return (test);
+}
+
+pmPSF_Test *pmPSF_TestModel (psArray *sources, char *modelName) {
+  
+    int Nflt = 0;
+    int Npsf = 0;
+    psVector *metricMask = psVectorAlloc (sources->n, PS_TYPE_U8);
+    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); 
+
+	// fit model as FLT, not PSF (skip poor fits)
+	if (!pmSourceFitModel (source, model, false)) continue;
+	test->modelFLT->data[i] = model;
+	Nflt ++;
+    }
+    psLogMsg ("psphot.psftest", 4, "fit flt:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
+    DumpModelFits (test->modelFLT, "modelsFLT.dat");
+
+    // stage 2: construct a psf (pmPSF) from this collection of model fits
+    pmPSFFromModels (test->psf, test->modelFLT);
+
+    // stage 3: refit with fixed shape parameters
+    psTimerStart ("fit");
+    for (int i = 0; i < test->sources->n; i++) {
+	psSource *source = test->sources->data[i];
+	psModel  *modelFLT = test->modelFLT->data[i];
+	if (modelFLT == NULL) continue;
+	psModel  *modelPSF = psModelFromPSF (modelFLT, test->psf); // set shape for this model
+
+	// fit model as PSF, not FLT (skip poor fits)
+	if (!pmSourceFitModel (source, modelPSF, true)) continue;
+	test->modelPSF->data[i] = modelPSF;
+	Npsf ++;
+    }
+    psLogMsg ("psphot.psftest", 4, "fit psf:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
+    DumpModelFits (test->modelPSF, "modelsPSF.dat");
+
+    // stage 4: measure PSF metric (Ap-Fit)
+    // should I be saving the aperture photometry for all sources?
+    // save the metric for all stars
+    for (int i = 0; i < test->sources->n; i++) {
+
+	float obsSum = 0;
+	float fitSum = 0;
+      
+	metricMask->data.U8[i] = 1;
+	test->metric->data.F64[i] = 0;
+
+	psModel *model = test->modelPSF->data[i];
+	if (model == NULL) continue;
+
+        psImage *image = ((psSource *)test->sources->data[i])->pixels;
+        psImage *mask  = ((psSource *)test->sources->data[i])->mask;
+
+	float sky = model->params->data.F32[0];
+
+	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;
+	    }
+	}
+	// this metric is Ap-Fit
+	// fprintf (stderr, "%d %f %f  %f\n", i, obsSum, fitSum, test->metric->data.F64[i]);
+	if ((fitSum > 0) && (obsSum > 0)) {
+	    metricMask->data.U8[i] = 0;
+	    test->metric->data.F64[i] = -2.5*log10(obsSum/fitSum);
+	}	
+    }
+
+    // XXX use robust stats, not clipped stats
+    psStats *stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
+    psVectorStats (stats, test->metric, NULL, metricMask, 1);
+    test->metricStats = stats;
+    psLogMsg ("psphot.pspsf", 3, "test model %s, metric: %f +/- %f\n", modelName, stats->clippedMean, stats->clippedStdev);
+    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) {
+
+    int n;
+
+    // construct the x and y vectors from the collection of objects
+    // count the number of valid model measurements
+    psVector *x = psVectorAlloc (models->n, PS_TYPE_F64);
+    psVector *y = psVectorAlloc (models->n, PS_TYPE_F64);
+    n = 0;
+    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[n] = model->params->data.F32[2];
+	y->data.F64[n] = model->params->data.F32[3];
+	n++;
+    }
+    x->n = y->n = n;
+
+    // we are doing a robust fit.  after each pass, we drop points which are 
+    // more deviant than three sigma.
+    psVector *z = psVectorAlloc (x->n, PS_TYPE_F64);
+    psVector *dz = psVectorAlloc (x->n, PS_TYPE_F64);
+    for (int i = 0; i < psf->params->n; i++) {
+	n = 0;
+	for (int j = 0; j < models->n; j++) {
+	    psModel *model = models->data[j];
+	    if (model == NULL) continue;
+	    z->data.F64[n] = model->params->data.F32[i + 4];
+	    dz->data.F64[n] = 1;
+	    n++;
+	    // XXX EAM : need to use actual errors?
+	    // XXX EAM : this is fragile: psf(Nparams) = model(Nparams) - 4
+	}
+	if (n != x->n) psAbort ("pmPSFFromModels", "programming error: mismatched number of NULL models\n");
+	z->n = n;
+	dz->n = n;
+
+	psf->params->data[i] = RobustFit2D (psf->params->data[i], x, y, z, dz);
+	// psPolynomial2DDump (psf->params->data[i]);
+    }
+    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);
+}
+
+// this was an attempt to measure ap-fit using image stats.  it was giving
+// strangely wrong answers.  
+# if (0)
+	psImage  *pixels = psImageCopy (NULL, source->pixels, source->pixels->type.type);
+	pixels->col0 = source->pixels->col0;
+	pixels->row0 = source->pixels->row0;
+	// psImageCopy does not retain the parent offset
+	// this copy is relative to the parent array.
+
+	pmSourceSubModel (pixels, source->mask, model, false);
+
+	// use a robust statistic or not?
+	psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN);
+	psImageStats (stats, pixels, NULL, 0);
+	// psImageStats (stats, pixels, source->mask, 0);
+	// XXX this should be the SUM, but we don't have such a statistic
+	// XXX given the mask, we need stat functions to return Npixels used
+# endif
Index: /trunk/psphot/src/setup.c
===================================================================
--- /trunk/psphot/src/setup.c	(revision 4114)
+++ /trunk/psphot/src/setup.c	(revision 4114)
@@ -0,0 +1,32 @@
+# include "psphot.h"
+
+bool setup (psImage **image, psMetadata **header, psMetadata **config, char **argv) 
+{
+    int       Nfail;
+    psFits   *file = NULL;		// FITS pointer to input image
+    psRegion  region = {0,0,0,0};	// a region representing the entire array
+
+    psTimerStart ("psphot");
+
+    *config = psMetadataParseConfig (NULL, &Nfail, argv[3], FALSE);
+
+    file   = psFitsAlloc (argv[1]);
+
+    *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);
+    *header = psFitsReadHeader (*header, file);
+    *image  = psFitsReadImage  (*image,  file, region, 0);
+    psFree (file);
+
+    // grab these values from the approrpiate location (image header if necessary)
+    bool  status   = false;
+    float RDNOISE  = pmConfigLookupF32 (&status, *config, *header, "RDNOISE");
+    float GAIN     = pmConfigLookupF32 (&status, *config, *header, "GAIN");
+    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);
+
+    psLogMsg ("psphot", 3, "load data: %f sec\n", psTimerMark ("psphot"));
+
+    return (true);
+}
Index: /trunk/psphot/src/source_moments.c
===================================================================
--- /trunk/psphot/src/source_moments.c	(revision 4114)
+++ /trunk/psphot/src/source_moments.c	(revision 4114)
@@ -0,0 +1,53 @@
+# include "psphot.h"
+
+psArray *source_moments (psImage *image, psMetadata *config, psArray *allpeaks) 
+{
+    bool     status  = false;
+    psArray *sources = NULL;
+    psArray *peaks   = NULL;
+
+    psTimerStart ("psphot");
+
+    // reduce the initial peaks to a subset of the better peaks
+    // is this the appropriate place for this culling?
+
+    float XBORDER  = psMetadataLookupF32 (&status, config, "XBORDER");
+    float YBORDER  = psMetadataLookupF32 (&status, config, "YBORDER");
+    float SATURATE = psMetadataLookupF32 (&status, config, "SATURATE");
+
+    psRegion *keep = psRegionAlloc (XBORDER, -XBORDER, YBORDER, -YBORDER);
+    keep           = psRegionForImage (keep, image, keep);
+    peaks          = pmPeaksSubset (allpeaks, 100000, keep);
+
+    // I need to exclude the saturated objects here, but include them again later...?
+    peaks          = pmPeaksSubset (allpeaks, 100000, keep);
+
+    // determine properties (sky, moments) of initial sources
+    float INNER    = psMetadataLookupF32 (&status, config, "INNER_RADIUS");
+    float OUTER    = psMetadataLookupF32 (&status, config, "OUTER_RADIUS");
+    float RADIUS   = psMetadataLookupF32 (&status, config, "PSF_MOMENTS_RADIUS");
+
+    sources = psArrayAlloc (peaks->n);
+    sources->n = 0;
+
+    // this should use ROBUST not SAMPLE median, but it is broken
+    for (int i = 0; i < peaks->n; i++) {
+	// provide the user arguments as metadata?
+	psSource *source  = pmSourceLocalSky (image, (psPeak *)peaks->data[i], PS_STAT_SAMPLE_MEDIAN, INNER, OUTER);
+	if (source == NULL) continue;
+	pmSourceMaskRegion (source, keep);
+	pmSourceMoments (source, RADIUS);
+	psArrayAdd (sources, 100, source);
+    }
+
+    // group into STAR, COSMIC, GALAXY, SATURATED
+    pmSourceRoughClass (sources, config);
+
+    // make this optional
+    DumpMoments (sources, "moments.dat");
+
+    psLogMsg ("psphot", 3, "moments: %f sec\n", psTimerMark ("psphot"));
+    psLogMsg ("psphot", 3, "%d sources\n", sources->n);
+
+    return (sources);
+}
Index: /trunk/psphot/src/subtract_galaxies.c
===================================================================
--- /trunk/psphot/src/subtract_galaxies.c	(revision 4114)
+++ /trunk/psphot/src/subtract_galaxies.c	(revision 4114)
@@ -0,0 +1,25 @@
+# include "psphot.h"
+
+// fit selected galaxy model (GAUSS) to all bright objects of type GALAXY
+
+bool subtract_galaxies (psArray *sources, psMetadata *config) 
+{ 
+    bool  status;
+
+    for (int i = 0; i < sources->n; i++) {
+	psSource *source = sources->data[i];
+	if (source->type != PS_SOURCE_GALAXY) continue;
+
+	psModel *model = source->modelFLT;
+	if (model == NULL) continue;
+
+	psImage  *pixels = source->pixels;
+
+	// subtract object, leave local sky
+	psF64 sky = model->params->data.F32[0];
+	model->params->data.F32[0] = 0;
+	pmSourceSubModel (pixels, source->mask, model, false);
+	model->params->data.F32[0] = sky;
+    }
+    return (true);
+}
Index: /trunk/psphot/src/subtract_psf_sources.c
===================================================================
--- /trunk/psphot/src/subtract_psf_sources.c	(revision 4114)
+++ /trunk/psphot/src/subtract_psf_sources.c	(revision 4114)
@@ -0,0 +1,39 @@
+# include "psphot.h"
+
+// subtract model fits from image (only PSFSTAR && BRIGHTSTAR)
+// need to control application of this with some thresholding
+
+bool subtract_psf_sources (psArray *sources)
+{ 
+    for (int i = 0; i < sources->n; i++) {
+	float sky;
+	psSource *source = sources->data[i];
+
+	if (source->type == PS_SOURCE_SATSTAR) {
+	  if (source->modelPSF == NULL) {
+	    fprintf (stderr, "missing model for %f, %f\n", source->moments->x, source->moments->y);
+	    continue;
+	  }
+	  fprintf (stderr, "subtracting model for %f, %f\n", source->moments->x, source->moments->y);
+	}	    
+
+
+	bool valid = false;
+	valid |= (source->type == PS_SOURCE_SATSTAR);
+	valid |= (source->type == PS_SOURCE_PSFSTAR);
+	valid |= (source->type == PS_SOURCE_BRIGHTSTAR);
+	if (!valid) continue;
+
+	psModel *model = source->modelPSF;
+	if (model == NULL) continue;
+
+	psImage  *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;
+    }
+    return(true);
+}
Index: /trunk/psphot/src/test_psf_scatter.c
===================================================================
--- /trunk/psphot/src/test_psf_scatter.c	(revision 4114)
+++ /trunk/psphot/src/test_psf_scatter.c	(revision 4114)
@@ -0,0 +1,81 @@
+# include "psphot.h"
+
+// measure the distribution of dparams[4,5] as a function of magnitude
+// use this as a test comparison with the theoretical distribution
+
+bool test_psf_scatter (psArray *sources)
+{
+    // analysis of PSF shape distribution
+    // select the PSF sources
+    int Npt = 0;
+    psVector *deltaSX = psVectorAlloc (100, PS_TYPE_F64);
+    psVector *deltaSY = psVectorAlloc (100, PS_TYPE_F64);
+    psVector *logSN   = psVectorAlloc (100, PS_TYPE_F64);
+    for (int i = 0; i < sources->n; i++) {
+	psSource *source = sources->data[i];
+	if (source->type != PS_SOURCE_PSFSTAR) continue;
+	if (source->modelPSF == NULL) continue;
+	deltaSX->data.F64[Npt] = source->modelPSF->dparams->data.F32[4];
+	deltaSY->data.F64[Npt] = source->modelPSF->dparams->data.F32[5];
+	logSN  ->data.F64[Npt] = log10(source->modelPSF->params->data.F32[1]/source->modelPSF->dparams->data.F32[1]);
+	Npt ++;
+	if (Npt == deltaSX->nalloc) {
+	    deltaSX = psVectorRealloc (deltaSX, deltaSX->nalloc + 100);
+	    deltaSY = psVectorRealloc (deltaSY, deltaSY->nalloc + 100);
+	    logSN   = psVectorRealloc (logSN,   logSN  ->nalloc + 100);
+	}		
+    }
+    if (Npt < 45) psAbort ("PSF analysis", "too few PSF stars for scatter analysis");
+    logSN->n = deltaSX->n = deltaSY->n = Npt;
+
+    // use a min number of bins or a min number of stars per bin?
+    // determine dynamic range of magnitudes
+    psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
+    psVectorStats (stats, logSN, NULL, NULL, 0);
+    float snMin = stats->min;
+    float snMax = stats->max;
+    float snDel = 25*(snMax - snMin) / logSN->n;
+    int   nBins = 1 + (int)((snMax - snMin) / snDel);
+
+    // generate a log-Sx, mag line and fit for appropriate mag bins
+    stats->options = PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV;
+    psVector *snBin  = psVectorAlloc (nBins, PS_TYPE_F64);
+    psVector *dsxBin = psVectorAlloc (nBins, PS_TYPE_F64);
+    psVector *dsyBin = psVectorAlloc (nBins, PS_TYPE_F64);
+    psVector *tsx    = psVectorAlloc (logSN->n, PS_TYPE_F64);
+    psVector *tsy    = psVectorAlloc (logSN->n, PS_TYPE_F64);
+    int bin = 0;
+    for (float sn = snMin; sn <= snMax - snDel; sn += snDel) {
+	tsx->n = tsy->n = 0;
+	for (int i = 0; i < logSN->n; i++) {
+	    if (logSN->data.F64[i] < sn) continue;
+	    if (logSN->data.F64[i] > sn + snDel) continue;
+	    tsx->data.F64[tsx->n] = deltaSX->data.F64[i];
+	    tsy->data.F64[tsy->n] = deltaSY->data.F64[i];
+	    tsx->n++;
+	    tsy->n++;
+	}
+	if (tsx->n < 5) continue;
+	snBin ->data.F64[bin] = sn + 0.5*snDel;
+	psVectorStats (stats, tsx, NULL, NULL, 0);
+	dsxBin->data.F64[bin] = log10(stats->clippedStdev);
+	psVectorStats (stats, tsy, NULL, NULL, 0);
+	dsyBin->data.F64[bin] = log10(stats->clippedStdev);
+	fprintf (stderr, "bin: %d, sn: %f, ndsx: %d, ndsy: %d, dsx: %f, dsy: %f\n",
+		 bin, sn + 0.5*snDel, tsx->n, tsy->n, dsxBin->data.F64[bin], dsyBin->data.F64[bin]);
+	bin ++;
+    }
+    dsxBin->n = bin;
+    dsyBin->n = bin;
+    snBin ->n = bin;
+
+    psPolynomial1D *dsxLine = psPolynomial1DAlloc (2, PS_POLYNOMIAL_ORD);
+    psVectorFitPolynomial1D (dsxLine, snBin, dsxBin, NULL);
+    psPolynomial1D *dsyLine = psPolynomial1DAlloc (2, PS_POLYNOMIAL_ORD);
+    psVectorFitPolynomial1D (dsyLine, snBin, dsyBin, NULL);
+    // these should have a slope of 1.0
+    psPolynomial1DDump (dsxLine);
+    psPolynomial1DDump (dsyLine);
+
+    return (true);
+}