Index: /trunk/psphot/src/Makefile.am
===================================================================
--- /trunk/psphot/src/Makefile.am	(revision 6850)
+++ /trunk/psphot/src/Makefile.am	(revision 6851)
@@ -41,4 +41,5 @@
 	psphotModelTest.c	\
 	psphotFitSet.c		\
+	psphotWeightBias.c	\
 	psphotCleanup.c	   	
 
Index: /trunk/psphot/src/psphot.c
===================================================================
--- /trunk/psphot/src/psphot.c	(revision 6850)
+++ /trunk/psphot/src/psphot.c	(revision 6851)
@@ -3,6 +3,4 @@
 // XXX need a better structure for handling optional sequences
 int main (int argc, char **argv) {
-
-    // headtest (argc, argv);
 
     psTimerStart ("complete");
@@ -27,19 +25,2 @@
     exit (0);
 }
-
-/** I/O test code
-
-    psFits *fits = psFitsOpen (argv[1], "r");
-    psMetadata *header = psFitsReadHeader (NULL, fits);
-    psFitsClose (fits);
-
-    psMetadata *new = psMetadataCopy (NULL, header);
-    psMetadataConfigWrite (new, "test.cnf");
-
-    fits = psFitsOpen ("test.fits", "w");
-    psFitsWriteHeaderNotImage (fits, new);
-    psFitsClose (fits);
-
-    exit (0);
-
-**/
Index: /trunk/psphot/src/psphot.h
===================================================================
--- /trunk/psphot/src/psphot.h	(revision 6850)
+++ /trunk/psphot/src/psphot.h	(revision 6851)
@@ -54,4 +54,5 @@
 bool 		psphotDumpMoments (psMetadata *config, psArray *sources);
 psMetadata     *psphotDefineHeader (psMetadata *config);
+bool            psphotWeightBias (pmReadout *readout, psArray *sources, psMetadata *config, pmPSF *psf);
 
 // PSF / DBL / EXT evaluation functions
@@ -63,14 +64,14 @@
 bool 		psphotInitLimitsPSF (psMetadata *config);
 bool		psphotInitLimitsEXT (psMetadata *config);
-bool		psphotFitBlend (pmReadout *readout, pmSource *source);
-bool		psphotFitBlob (pmReadout *readout, pmSource *source, psArray *sources);
-bool		psphotFitPSF (pmReadout *readout, pmSource *source);
+bool		psphotFitBlend (pmReadout *readout, pmSource *source, pmPSF *psf);
+bool		psphotFitBlob (pmReadout *readout, pmSource *source, psArray *sources, pmPSF *psf);
+bool		psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf);
 pmModel	       *psphotFitEXT (pmReadout *readout, pmSource *source);
 psArray	       *psphotFitDBL (pmReadout *readout, pmSource *source);
 
+// functions to support simultaneous multi-source fitting
+bool 		psphotFitSet (pmSource *oneSrc, pmModel *oneModel, char *fitset, pmSourceFitMode mode);
+
 // XXX these can probably be dropped:
-
-// functions to support simultaneous multi-source fitting
-bool 		psphotFitSet (pmSource *oneSrc, pmModel *oneModel, char *fitset, bool PSF);
 
 # if (0)
Index: /trunk/psphot/src/psphotArguments.c
===================================================================
--- /trunk/psphot/src/psphotArguments.c	(revision 6850)
+++ /trunk/psphot/src/psphotArguments.c	(revision 6851)
@@ -10,4 +10,5 @@
 
     int N;
+    bool status;
 
     if (*argc == 1) usage ();
@@ -118,57 +119,7 @@
     }
 
-    // we load all input files onto a psArray, to be parsed later
-    psArray *input = psArrayAlloc (16);
-    input->n = 0;
-
-    // load the list of filenames the supplied file (may be a glob: "file*.fits")
-    if ((N = psArgumentGet (*argc, argv, "-file"))) {
-	glob_t globList;
-	psArgumentRemove (N, argc, argv);
-	globList.gl_offs = 0;
-	glob (argv[N], 0, NULL, &globList);
-	for (int i = 0; i < globList.gl_pathc; i++) {
-	    char *filename = psStringCopy (globList.gl_pathv[i]);
-	    psArrayAdd (input, 16, filename);
-	    psFree (filename);
-	}
-	psArgumentRemove (N, argc, argv);
-    }
-
-    // load the list from the supplied text file
-    if ((N = psArgumentGet (*argc, argv, "-list"))) {
-	int nItems;
-	char line[1024]; // XXX limits the list lines to 1024 chars
-	char word[1024];
-	char *filename;
-
-	psArgumentRemove (N, argc, argv);
-	FILE *f = fopen (argv[N], "r");
-	if (f == NULL) {
-	    psAbort ("psphot", "unable to open specified list file");
-	}
-	while (fgets (line, 1024, f) != NULL) {
-	    nItems = sscanf (line, "%s", word);
-	    switch (nItems) {
-	      case 0:
-		break;
-	      case 1:
-		filename = psStringCopy (word);
-		psArrayAdd (input, 16, filename);
-		psFree (filename);
-		break;
-	      default:
-		// rigid format, no comments allowed?
-		psAbort ("psphot", "error parsing input list file");
-		break;
-	    }
-	}
-	psArgumentRemove (N, argc, argv);
-    }		
-    if (input->n < 1) usage ();
-
-    // input list gets places as an array on the config->arguements list
-    psMetadataAddPtr (config->arguments, PS_LIST_TAIL, "INPUT",  PS_DATA_ARRAY, "", input);
-    psFree (input);
+    // the input file is a required argument; if not found, we will exit
+    status = pmConfigFileSetsMD (config->arguments, argc, argv, "INPUT", "-file", "-list");
+    if (!status) { usage ();}
 
     if (*argc != 2) usage ();
Index: /trunk/psphot/src/psphotBlendFit.c
===================================================================
--- /trunk/psphot/src/psphotBlendFit.c	(revision 6850)
+++ /trunk/psphot/src/psphotBlendFit.c	(revision 6851)
@@ -54,6 +54,6 @@
 	
 	// try fitting PSFs, then try extended sources
-	if (psphotFitBlend (readout, source)) continue;
-	if (psphotFitBlob (readout, source, sources)) continue;
+	if (psphotFitBlend (readout, source, psf)) continue;
+	if (psphotFitBlob (readout, source, sources, psf)) continue;
 
 	psTrace ("psphot.blend", 5, "failed fits\n");
Index: /trunk/psphot/src/psphotChoosePSF.c
===================================================================
--- /trunk/psphot/src/psphotChoosePSF.c	(revision 6850)
+++ /trunk/psphot/src/psphotChoosePSF.c	(revision 6851)
@@ -27,4 +27,9 @@
     int NSTARS = psMetadataLookupS32 (&status, config, "PSF_MAX_NSTARS");
     if (!status) NSTARS = PS_MIN (sources->n, 200);
+
+    // use poissonian errors or local-sky errors
+    bool POISSON_ERRORS = psMetadataLookupBool (&status, config, "POISSON_ERRORS");
+    if (!status) POISSON_ERRORS = true;
+    pmSourceFitModelInit (15, 0.1, POISSON_ERRORS);
 
     stars = psArrayAlloc (sources->n);
@@ -63,5 +68,5 @@
 	psMetadataItem *item = psListGetAndIncrement (iter);
 	modelName = item->data.V;
-	models->data[i] = pmPSFtryModel (stars, modelName, RADIUS);
+	models->data[i] = pmPSFtryModel (stars, modelName, RADIUS, POISSON_ERRORS);
     }
     psFree (iter);
Index: /trunk/psphot/src/psphotEvalPSF.c
===================================================================
--- /trunk/psphot/src/psphotEvalPSF.c	(revision 6850)
+++ /trunk/psphot/src/psphotEvalPSF.c	(revision 6851)
@@ -109,5 +109,5 @@
     dSX = model->dparams->data.F32[4];
     dSY = model->dparams->data.F32[5];
-    Chi = model->chisq / model->nDOF;
+    Chi = model->chisqNorm / model->nDOF;
 
     // swing of sigma_x,y in sigmas
Index: /trunk/psphot/src/psphotFitSet.c
===================================================================
--- /trunk/psphot/src/psphotFitSet.c	(revision 6850)
+++ /trunk/psphot/src/psphotFitSet.c	(revision 6851)
@@ -1,5 +1,5 @@
 # include "psphot.h"
 
-bool psphotFitSet (pmSource *source, pmModel *oneModel, char *fitset, bool PSF) {
+bool psphotFitSet (pmSource *source, pmModel *oneModel, char *fitset, pmSourceFitMode mode) {
 
     double x, y, Io;
@@ -24,5 +24,5 @@
     }
 
-    pmSourceFitSet (source, modelSet, PSF);
+    pmSourceFitSet (source, modelSet, mode);
 
     // write out positive object
Index: /trunk/psphot/src/psphotImageMedian.c
===================================================================
--- /trunk/psphot/src/psphotImageMedian.c	(revision 6850)
+++ /trunk/psphot/src/psphotImageMedian.c	(revision 6851)
@@ -1,11 +1,3 @@
 # include "psphot.h"
-double psImageClippedStats (psImage *image, psImage *mask, psU8 maskValue, double fmin, double fmax);
-void fsort (float *value, int N);
-
-// random number seed to select a fraction of the image pixels
-static psRandom *rnd;
-
-// use no more than MAX_SAMPLE_PIXELS pixels for each median box
-static int MAX_SAMPLE_PIXELS;
 
 // generate the median in NxN boxes, clipping heavily
@@ -15,19 +7,24 @@
     bool status;
     psRegion region;
-    psImage *model = NULL;
+    int MAX_SAMPLE_PIXELS;
 
     psTimerStart ("psphot");
 
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, "PSPHOT");
+
+    MAX_SAMPLE_PIXELS = psMetadataLookupF32 (&status, recipe, "IMSTATS_NPIX");
+    if (!status) MAX_SAMPLE_PIXELS = 1000;
+    psImageClippedStatsInit(MAX_SAMPLE_PIXELS);
+
+    // subtract this amount extra from the sky
+    float SKY_BIAS = psMetadataLookupF32 (&status, recipe, "SKY_BIAS");
+    if (!status) SKY_BIAS = 0;
+
     // find the currently selected readout
-    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, "PSPHOT");
-    pmFPAfile  *input   = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT");
-    pmReadout  *readout = pmFPAviewThisReadout (view, input->fpa);
+    pmReadout  *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT");
 
     psImage *image = readout->image;
     psImage *mask  = readout->mask;
-
-    rnd = psRandomAlloc (PS_RANDOM_TAUS, 0);
-    MAX_SAMPLE_PIXELS = psMetadataLookupF32 (&status, recipe, "IMSTATS_NPIX");
-    if (!status) MAX_SAMPLE_PIXELS = 1000;
 
     // dimensions of input & output image
@@ -49,6 +46,14 @@
     int ny = (Ny % DY) ? (int)(Ny / DY) + 1 : Ny / DY;
 
-    // select model pixels, from output background model file, or create
-    model = pmFPAfileReadoutImage (config->files, view, "PSPHOT.BACKMDL", nx, ny, PS_TYPE_F32);
+    // select model pixels (from output background model file, or create internal file)
+    pmReadout *model = pmFPAfileThisReadout (config->files, view, "PSPHOT.BACKMDL");
+    if (model == NULL) {
+	// select model pixels, from output background model file, or create
+	model = pmFPAfileCreateInternal (config->files, "PSPHOT.BACKMDL", nx, ny, PS_TYPE_F32);
+    } else {
+	// replace the supplied image data with an image of the desired size
+	psImageRecycle (model->image, nx, ny, PS_TYPE_F32);
+    }
+    psF32 **modelData = model->image->data.F32;
 
     // measure clipped median for subimages 
@@ -63,5 +68,6 @@
 	    psImage *submask = psImageSubset (mask, region);
 
-	    model->data.F32[iy][ix] = psImageClippedStats (subset, submask, 0xff, 0.25, 0.50);
+	    // XXX the value of the upper and lower cuts probably should be studied...
+	    modelData[iy][ix] = psImageClippedStats (subset, submask, 0xff, 0.25, 0.75) + SKY_BIAS;
 
 	    psFree (subset);
@@ -76,151 +82,36 @@
     psLogMsg ("psphot", 3, "build median image: %f sec\n", psTimerMark ("psphot"));
 
-    // linear interpolation to full-scale
+    // XXX temporarily until we get the pmFPAfile output finished
+    psphotSaveImage (NULL, model->image, "model.fits");
 
     // select background pixels, from output background file, or create
-    psImage *background = pmFPAfileReadoutImage (config->files, view, "PSPHOT.BACKGND", Nx, Ny, PS_TYPE_F32);
-
-    // XXX this code skips the initial pixels
-    for (int Iy = 0; Iy < ny-1; Iy ++) {
-	for (int Ix = 0; Ix < nx-1; Ix ++) {
-
-	    float V00 = model->data.F32[Iy+0][Ix+0];
-	    float V01 = model->data.F32[Iy+0][Ix+1];
-	    float V10 = model->data.F32[Iy+1][Ix+0];
-	    float V11 = model->data.F32[Iy+1][Ix+1];
-
-	    // a single binned pixel quad
-	    // (Xs,Ys) : (Xe,Ye) : binned pixel centers in unbinned coords
-	    // corresponding to (Ix,Iy), (Ix+1,Iy+1)
-	    int Xs = (Ix + 1)*DX - dx;
-	    int Ys = (Iy + 1)*DY - dy;
-	    int Xe = Xs + DX;
-	    int Ye = Ys + DY;
-
-	    for (int iy = Ys; (iy < Ye) && (iy < Ny); iy++) {
-		float Vxs = (V10 - V00)*(iy - Ys) / DY + V00;
-		float Vxe = (V11 - V01)*(iy - Ys) / DY + V01;
-		float dV = (Vxe - Vxs) / DX;
-		float V  = Vxs;
-		for (int ix = Xs; (ix < Xe) && (ix < Nx); ix++) {
-		    background->data.F32[iy][ix] = V;
-		    V += dV;
-		}
-	    }
-	}
+    pmReadout *background = pmFPAfileThisReadout (config->files, view, "PSPHOT.BACKGND");
+    if (background == NULL) {
+	background = pmFPAfileCreateInternal (config->files, "PSPHOT.BACKGND", Nx, Ny, PS_TYPE_F32);
+    } else {
+	// replace the supplied image data with an image of the desired size
+	// XXX should not have to do this: it should be allocated correctly
+	// psImageRecycle (background->image, Nx, Ny, PS_TYPE_F32);
     }
 
-    // side pixels
-    int Xs = DX - dx;
-    int Xe = nx*DX - dx;
-    for (int Iy = 0; Iy < ny - 1; Iy++) {
+    psF32 **backData = background->image->data.F32;
 
-	int Ys = (Iy + 1)*DY - dy;
-	int Ye = Ys + DY;
-
-	// leading edge
-	float V0 = model->data.F32[Iy+0][0];
-	float V1 = model->data.F32[Iy+1][0];
-	float dV = (V1 - V0) / DY;
-	float V = V0;
-	for (int iy = Ys; (iy < Ye) && (iy < Ny); iy++) {
-	    for (int ix = 0; ix < Xs; ix++) {
-		background->data.F32[iy][ix] = V;
-	    }
-	    V += dV;
-	}
-
-	// trailing edge
-	V0 = model->data.F32[Iy+0][nx-1];
-	V1 = model->data.F32[Iy+1][nx-1];
-	dV = (V1 - V0) / DY;
-	V = V0;
-	for (int iy = Ys; (iy < Ye) && (iy < Ny); iy++) {
-	    for (int ix = Xe; ix < Nx; ix++) {
-		background->data.F32[iy][ix] = V;
-	    }
-	    V += dV;
-	}
-    }
-
-    // top and bottom pixels
-    int Ys = DY - dy;
-    int Ye = ny*DY - dy;
-    for (int Ix = 0; Ix < nx - 1; Ix++) {
-
-	int Xs = (Ix + 1)*DX - dx;
-	int Xe = Xs + DX;
-
-	// top edge
-	float V0 = model->data.F32[0][Ix+0];
-	float V1 = model->data.F32[0][Ix+1];
-	float dV = (V1 - V0) / DX;
-	for (int iy = 0; iy < Ys; iy++) {
-	    float V = V0;
-	    for (int ix = Xs; (ix < Xe) && (ix < Nx); ix++) {
-		background->data.F32[iy][ix] = V;
-		V += dV;
-	    }
-	}
-
-	// bottom edge
-	V0 = model->data.F32[ny-1][Ix+0];
-	V1 = model->data.F32[ny-1][Ix+1];
-	dV = (V1 - V0) / DX;
-	for (int iy = Ye; iy < Ny; iy++) {
-	    float V = V0;
-	    for (int ix = Xs; (ix < Xe) && (ix < Nx); ix++) {
-		background->data.F32[iy][ix] = V;
-		V += dV;
-	    }
-	}
-    }
-
-    // the four corners
-    {
-	float V;
-	// 0,0
-	V = model->data.F32[0][0];
-	for (int iy = 0; iy < DY - dy; iy++) {
-	    for (int ix = 0; ix < DX - dx; ix++) {
-		background->data.F32[iy][ix] = V;
-	    }
-	}
-	// Nx,0
-	V = model->data.F32[0][nx-1];
-	for (int iy = 0; iy < DY - dy; iy++) {
-	    for (int ix = nx*DX - dx; ix < Nx; ix++) {
-		background->data.F32[iy][ix] = V;
-	    }
-	}
-	// 0,Ny
-	V = model->data.F32[ny-1][0];
-	for (int iy = ny*DY - dy; iy < Ny; iy++) {
-	    for (int ix = 0; ix < DX - dx; ix++) {
-		background->data.F32[iy][ix] = V;
-	    }
-	}
-	// Nx,Ny
-	V = model->data.F32[nx-1][ny-1];
-	for (int iy = ny*DY - dy; iy < Ny; iy++) {
-	    for (int ix = nx*DX - dx; ix < Nx; ix++) {
-		background->data.F32[iy][ix] = V;
-	    }
-	}
-    }	
+    // linear interpolation to full-scale
+    psImageUnbin (background->image, model->image, DX, DY, dx, dy);
     psLogMsg ("psphot", 3, "build resampled image: %f sec\n", psTimerMark ("psphot"));
 
-    psphotSaveImage (NULL, background, "back.fits");
+    // XXX temporarily until we get the pmFPAfile output finished...
+    psphotSaveImage (NULL, background->image, "back.fits");
 
     // back-sub image pixels, from output background file (don't create if not requested)
-    psImage *backSub = pmFPAfileReadoutImage (config->files, view, "PSPHOT.BACKSUB", 0, 0, PS_TYPE_F32);
+    pmReadout *backSub = pmFPAfileThisReadout (config->files, view, "PSPHOT.BACKSUB");
 
-    // subtract the background model
+    // subtract the background model (save in backSub, if requested)
     for (int j = 0; j < image->numRows; j++) {
 	for (int i = 0; i < image->numCols; i++) {
 	    if (!mask->data.U8[j][i]) {
-		image->data.F32[j][i] -= background->data.F32[j][i];
+		image->data.F32[j][i] -= backData[j][i];
 		if (backSub) {
-		    backSub->data.F32[j][i] = image->data.F32[j][i];
+		    backSub->image->data.F32[j][i] = image->data.F32[j][i];
 		}
 	    }
@@ -229,96 +120,7 @@
 
     psLogMsg ("psphot", 3, "subtracted background model: %f sec\n", psTimerMark ("psphot"));
-    psFree (rnd);
+    psImageClippedStatsCleanup();
 
-    // it is safe to free all of these: either they are saved on config->files or are local
-    psFree (model);
-    psFree (background);
-    psFree (backSub);
+    // the pmReadout selected in this function are all view on entries in config->files
     return true;
 }
-
-double psImageClippedStats (psImage *image, psImage *mask, psU8 maskValue, double fmin, double fmax) {
-
-    double value;
-    int nx = image->numCols;
-    int ny = image->numRows;
-    
-    if (nx*ny <= 0) return 0.0;
-
-    int Nsubset = PS_MIN (MAX_SAMPLE_PIXELS, nx*ny);
-    int Npixels = nx*ny;
-
-    psVector *values = psVectorAlloc (Nsubset, PS_TYPE_F32);
-
-    float min = values->data.F32[0];
-    float max = values->data.F32[0];
-
-    int n = 0;
-    for (int i = 0; i < Nsubset; i++) {
-	double frnd = psRandomUniform (rnd);
-	int pixel = Npixels * frnd;
-	int ix = pixel % nx;
-	int iy = pixel / nx;
-
-	if (mask->data.U8[iy][ix] & maskValue) continue;
-
-	value = image->data.F32[iy][ix];
-	min = PS_MIN (value, min);
-	max = PS_MIN (value, max);
-	values->data.F32[n] = value;
-	n++;
-    }
-    values->n = n;
-
-    int imin = fmin * n;
-    int imax = fmax * n;
-    int npts = imax - imin + 1;
-
-    fsort (values->data.F32, n);
-
-    value = 0;
-    for (int i = imin; (i <= imax) && (i < n); i++) {
-	value += values->data.F32[i];
-    }
-    value = value / npts;
-    
-    // XXX correct for selection bias??
-
-    psFree (values);
-    return value;
-}
-
-void fsort (float *value, int N) {
-
-    int l,j,ir,i;
-    float temp;
-  
-    if (N < 2) return;
-    l = N >> 1;
-    ir = N - 1;
-    for (;;) {
-	if (l > 0) {
-	    temp = value[--l];
-	}
-	else {
-	    temp = value[ir];
-	    value[ir] = value[0];
-	    if (--ir == 0) {
-		value[0] = temp;
-		return;
-	    }
-	}
-	i = l;
-	j = (l << 1) + 1;
-	while (j <= ir) {
-	    if (j < ir && value[j] < value[j+1]) ++j;
-	    if (temp < value[j]) {
-		value[i]=value[j];
-		j += (i=j) + 1;
-	    }
-	    else j = ir + 1;
-	}
-	value[i] = temp;
-    }
-}
-
Index: /trunk/psphot/src/psphotModelTest.c
===================================================================
--- /trunk/psphot/src/psphotModelTest.c	(revision 6850)
+++ /trunk/psphot/src/psphotModelTest.c	(revision 6851)
@@ -10,6 +10,12 @@
     char name[64];
     pmPSF *psf = NULL;
+    pmSourceFitMode fitMode;
 
     psMetadataItem *item  = NULL;
+
+    // use poissonian errors or local-sky errors
+    bool POISSON_ERRORS = psMetadataLookupBool (&status, recipe, "POISSON_ERRORS");
+    if (!status) POISSON_ERRORS = true;
+    pmSourceFitModelInit (15, 0.1, POISSON_ERRORS);
 
     // run model fitting tests on a single source
@@ -17,12 +23,13 @@
 
     // what fitting mode to use?
-    char *psfModeWord = psMetadataLookupStr (&status, recipe, "TEST_FIT_MODE");
+    char *fitModeWord = psMetadataLookupStr (&status, recipe, "TEST_FIT_MODE");
     if (!status) {
-	psfModeWord = DEFAULT_MODE;
+	fitModeWord = DEFAULT_MODE;
     }
-    bool psfMode = !strcasecmp (psfModeWord, "PSF");
+    fitMode = PM_SOURCE_FIT_EXT;
+    if (!strcasecmp (fitModeWord, "PSF")) fitMode = PM_SOURCE_FIT_PSF;
 
-    // in psfMode, psf sets the model type
-    if (psfMode) {
+    // in fitMode, psf sets the model type
+    if (fitMode == PM_SOURCE_FIT_PSF) {
 	char *psfFile = psMetadataLookupStr (&status, recipe, "PSF_INPUT_FILE");
 	if (!status) psAbort ("psphotModelTest", "PSF_INPUT_FILE not supplied");
@@ -128,5 +135,5 @@
     fprintf (stderr, "peak: %f @ (%f, %f)\n", source->moments->Sum*area, (double)source->peak->x, (double)source->peak->y);
 
-    if (psfMode) {
+    if (fitMode == PM_SOURCE_FIT_PSF) {
 	pmModel *modelPSF = pmModelFromPSF (model, psf);
 	psFree (model);
@@ -155,9 +162,9 @@
     char *fitset = psMetadataLookupStr (&status, recipe, "TEST_FIT_SET");
     if (status) {
-	status = psphotFitSet (source, model, fitset, psfMode);
+	status = psphotFitSet (source, model, fitset, fitMode);
 	exit (0);
     }
 
-    status = pmSourceFitModel (source, model, psfMode);
+    status = pmSourceFitModel (source, model, fitMode);
 
     // measure the source mags
Index: /trunk/psphot/src/psphotParseCamera.c
===================================================================
--- /trunk/psphot/src/psphotParseCamera.c	(revision 6850)
+++ /trunk/psphot/src/psphotParseCamera.c	(revision 6851)
@@ -4,86 +4,10 @@
 
     bool status;
-    psFits *fits = NULL;
-    psMetadata *phu = NULL;
-    psMetadata *format = NULL;
 
-    // psphot is supplied with a list of input images (may be only one image)
-    psArray *infiles = psMetadataLookupPtr(&status, config->arguments, "INPUT");
-    if (!status) psAbort (__func__, "missing INPUT entry");
-    if (infiles->n < 1) psAbort (__func__, "empty input list");
-
-    // if no camera has been specified, use the first image as a template for the rest.
-    if (config->camera == NULL) {
-	fits = psFitsOpen (infiles->data[0], "r");
-	phu = psFitsReadHeader (NULL, fits);
-	format = pmConfigCameraFormatFromHeader (config, phu);
-	psFitsClose (fits);
+    status = false;
+    pmFPAfile *input = pmFPAfileFromArgs (&status, config, "PSPHOT.INPUT", "INPUT");
+    if (!status) {
+	psAbort (__func__, "missing INPUT entry");
     }
-    // There's no point in continuing if we can't recognize what we've got
-    if (! config->camera) {
-	psErrorStackPrint(stderr, "Can't find camera configuration!\n");
-	exit(EXIT_FAILURE);
-    }
-
-    // files
-    config->files = psMetadataAlloc ();
-
-    // build the template fpa, set up the basic view
-    // XXX : 216 leaks in pmFPAConstruct / psFree (input)
-    pmFPA *input = pmFPAConstruct (config->camera);
-
-    // assign the I/O files (potentially) needed by psphot
-    // the output file is just a view to the file on config->files 
-    pmFPAfile *file = pmFPAfileDefine (config->files, format, input, "PSPHOT.INPUT");
-    if (!file) {
-	exit(EXIT_FAILURE);
-    }
-
-    for (int i = 0; i < infiles->n; i++) {
-	if (phu == NULL) {
-	    fits = psFitsOpen (infiles->data[i], "r");
-	    phu = psFitsReadHeader (NULL, fits);
-	    pmConfigValidateCameraFormat (format, phu);
-	    psFitsClose (fits);
-	}
-
-	// set the view to the corresponding entry for this phu
-	pmFPAview *view = pmFPAAddSource (input, phu, format);
-
-	// XXX is this the correct psMD to save the filename?
-	char *name = pmFPAfileNameFromRule (file->filextra, file, view);
-	psMetadataAddStr (file->names, PS_LIST_TAIL, name, 0, "", infiles->data[i]);
-
-	psFree (view);
-	psFree (name);
-	psFree (phu);
-	phu = NULL;
-    }
-
-    pmFPAfileDefine (config->files, format, input, "PSPHOT.OUTPUT");
-    pmFPAfileDefine (config->files, format, input, "PSPHOT.RESID");
-
-    // pmFPAfileDefine (config->files, format, input, "PSPHOT.PSF_INPUT");
-    // pmFPAfileDefine (config->files, format, input, "PSPHOT.PSF_OUTPUT");
-
-    // build the template fpa, set up the basic view
-    // supply the backgnd with a different camera?
-    // allow alternate format?
-    // pmFPAfileConstruct (config->files, format, config->camera, "PSPHOT.BACKSUB");
-    // pmFPAfileConstruct (config->files, format, config->camera, "PSPHOT.BACKGND");
-    // pmFPAfileConstruct (config->files, format, config->camera, "PSPHOT.BACKMDL");
-    // pmFPAfileConstruct (config->files, format, config->camera, "PSPHOT.PSF_SAMPLE");
-
-    // psphot is supplied with the output name
-    // this needs to be better: supply to all output (WRITE) files?
-    char *output = psMetadataLookupPtr(&status, config->arguments, "OUTPUT");
-
-    file = psMetadataLookupPtr (&status, config->files, "PSPHOT.OUTPUT");
-    if (!status) psAbort (__func__, "missing OUTPUT entry");
-    psMetadataAddStr (file->names, PS_LIST_TAIL, "OUTPUT", 0, "", output);
-
-    file = psMetadataLookupPtr (&status, config->files, "PSPHOT.RESID");
-    if (!status) psAbort (__func__, "missing OUTPUT entry");
-    psMetadataAddStr (file->names, PS_LIST_TAIL, "OUTPUT", 0, "", output);
 
     // recipe override values (command-line options):
@@ -96,4 +20,27 @@
     }
     psFree (iter);
+
+    // these calls bind the I/O handle to the specified fpa
+    pmFPAfileDefine (config->files, config->camera, input->fpa, "PSPHOT.OUTPUT");
+    pmFPAfileDefine (config->files, config->camera, input->fpa, "PSPHOT.RESID");
+
+    // pmFPAfileDefine (config->files, format, input, "PSPHOT.PSF_INPUT");
+    // pmFPAfileDefine (config->files, format, input, "PSPHOT.PSF_OUTPUT");
+
+    // these calls construct a new fpa for the I/O handle 
+    int DX = psMetadataLookupS32 (&status, recipe, "BACKGROUND_XBIN");
+    if (!status) {DX = 64;}
+    int DY = psMetadataLookupS32 (&status, recipe, "BACKGROUND_YBIN");
+    if (!status) {DY = 64;}
+    pmFPAfileFromFPA (config, input->fpa, DX, DY, "PSPHOT.BACKMDL");
+    pmFPAfileFromFPA (config, input->fpa,  1,  1, "PSPHOT.BACKGND");
+    pmFPAfileFromFPA (config, input->fpa,  1,  1, "PSPHOT.BACKSUB");
+
+    // pmFPAfileConstruct (config->files, format, config->camera, "PSPHOT.PSF_SAMPLE");
+
+    // psphot is supplied with the output name
+    // this needs to be better: supply to all output (WRITE) files?
+    char *output = psMetadataLookupPtr(&status, config->arguments, "OUTPUT");
+    pmFPAfileAddFileNames (config->files, "OUTPUT", output, PM_FPA_MODE_WRITE);
 
     // set default recipe values here
@@ -129,7 +76,4 @@
     # endif
 
-    psFree (input);
-    psFree (format);
-
     psTrace(__func__, 1, "Done with psphotParseCamera...\n");
     return true;
Index: /trunk/psphot/src/psphotReadout.c
===================================================================
--- /trunk/psphot/src/psphotReadout.c	(revision 6850)
+++ /trunk/psphot/src/psphotReadout.c	(revision 6851)
@@ -9,11 +9,14 @@
     bool         status;
 
+    // select the current recipe
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, "PSPHOT");
+
     // find the currently selected readout
-    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, "PSPHOT");
-    pmFPAfile  *input   = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT");
-    pmReadout  *readout = pmFPAviewThisReadout (view, input->fpa);
+    pmReadout  *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT");
 
     // XXX does this need to invoke I/O?
+    // XXX move this input the psphotImageLoop level?
     pmReadoutSetWeights (readout);
+    psphotSaveImage (NULL, readout->weight, "weight.fits");
 
     // I have a valid mask, now mask in the analysis region of interest
@@ -49,4 +52,6 @@
     psphotBlendFit (readout, sources, recipe, psf);
 
+    // psphotWeightBias (readout, sources, recipe, psf);
+
     // replace fitted sources
     psphotReplaceUnfit (sources);
@@ -75,4 +80,8 @@
     status = psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.PSF",     PS_DATA_UNKNOWN,  "psphot psf", psf);
 
+    // remove internal pmFPAfiles, if created
+    pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL");
+    pmFPAfileDropInternal (config->files, "PSPHOT.BACKGND");
+
     // free up the local copies of the data
     psFree (psf);
Index: /trunk/psphot/src/psphotSkyReplace.c
===================================================================
--- /trunk/psphot/src/psphotSkyReplace.c	(revision 6850)
+++ /trunk/psphot/src/psphotSkyReplace.c	(revision 6851)
@@ -5,31 +5,25 @@
 bool psphotSkyReplace (pmConfig *config, pmFPAview *view) {
 
-    bool status;
+    // find the currently selected readout
+    pmReadout *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT");
+    if (readout == NULL) psAbort ("psphot", "input not defined");
 
-    // find the currently selected readout
-    //  psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, "PSPHOT");
-    pmFPAfile *input   = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT");
-    pmReadout *readout = pmFPAviewThisReadout (view, input->fpa);
+    // select background pixels, from output background file, or create
+    pmReadout *background = pmFPAfileThisReadout (config->files, view, "PSPHOT.BACKGND");
+    if (background == NULL) psAbort ("psphot", "background not defined");
 
     // select the corresponding images
-    psImage *image = readout->image;
-    psImage *mask  = readout->mask;
-
-    // select background pixels, from output background file, or create
-    psImage *background = pmFPAfileReadoutImage (config->files, view, "PSPHOT.BACKGND", 0, 0, PS_TYPE_F32);
-    if (background == NULL) {
-	return false;
-    }
+    psF32 **image = readout->image->data.F32;
+    psU8  **mask  = readout->mask->data.U8;
+    psF32 **back  = background->image->data.F32;
 
     // replace the background model
-    for (int j = 0; j < image->numRows; j++) {
-	for (int i = 0; i < image->numCols; i++) {
-	    if (!mask->data.U8[j][i]) {
-		image->data.F32[j][i] += background->data.F32[j][i];
+    for (int j = 0; j < readout->image->numRows; j++) {
+	for (int i = 0; i < readout->image->numCols; i++) {
+	    if (!mask[j][i]) {
+		image[j][i] += back[j][i];
 	    }
 	}
     }
-    
-    psFree (background);
     return true;
 }
Index: /trunk/psphot/src/psphotSourceFits.c
===================================================================
--- /trunk/psphot/src/psphotSourceFits.c	(revision 6850)
+++ /trunk/psphot/src/psphotSourceFits.c	(revision 6851)
@@ -3,7 +3,8 @@
 // given a source with an existing modelPSF, attempt a full PSF fit, subtract if successful
 
-bool psphotFitPSF (pmReadout *readout, pmSource *source) { 
+bool psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf) { 
 
     float x, y;
+    double chiTrend;
 
     // save the PSF model from the Ensemble fit
@@ -18,6 +19,11 @@
     // fit PSF model (set/unset the pixel mask)
     psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "OR", PM_SOURCE_MASK_MARKED);
-    pmSourceFitModel (source, PSF, true);
+    pmSourceFitModel (source, PSF, PM_SOURCE_FIT_PSF);
     psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "AND", ~PM_SOURCE_MASK_MARKED);
+    
+    // correct model chisq for flux trend
+    chiTrend = psPolynomial1DEval (psf->ChiTrend, PSF->params->data.F32[1]);
+    PSF->chisqNorm = PSF->chisq / chiTrend;
+    // PSF->chisqNorm = PSF->chisq;
 
     // does the PSF model succeed?
@@ -60,8 +66,10 @@
 }
 
-bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *sources) { 
+bool psphotFitBlob (pmReadout *readout, pmSource *source, psArray *sources, pmPSF *psf) { 
 
     bool okEXT, okDBL;
     float chiEXT, chiDBL;
+    double chiTrend;
+    pmModel *ONE = NULL;
 
     // skip the source if we don't think it is extended
@@ -86,6 +94,19 @@
     okDBL  = psphotEvalDBL (tmpSrc, DBL->data[0]);
     okDBL &= psphotEvalDBL (tmpSrc, DBL->data[1]);
-    pmModel *ONE = DBL->data[0];
-    chiDBL = ONE->chisq / ONE->nDOF;
+
+    // correct first model chisqs for flux trend
+    ONE = DBL->data[0];
+    chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]);
+    ONE->chisqNorm = ONE->chisq / chiTrend;
+    // ONE->chisqNorm = ONE->chisq;
+
+    // save chisq for double-star/galaxy comparison
+    chiDBL = ONE->chisqNorm;
+
+    // correct second model chisqs for flux trend
+    ONE = DBL->data[1];
+    chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]);
+    ONE->chisqNorm = ONE->chisq / chiTrend;
+    // ONE->chisqNorm = ONE->chisq;
 
     psFree (tmpSrc); // XXX should I keep / save the flags set in the eval functions?
@@ -93,5 +114,5 @@
     if (okEXT && okDBL) {
 	psTrace ("psphot.blend", 5, "blob chisq: %f vs %f\n", chiEXT, chiDBL);
-	// XXX EAM : a bogus bias: need to examine this better
+	// XXX EAM : a bogus bias: need to examine this better 
 	if (3*chiEXT > chiDBL) goto keepDBL;
 	goto keepEXT;
@@ -188,7 +209,7 @@
     y = source->moments->y;
 
-    // fit EXT (not PSF) model (set/unset the pixel mask)
+    // fit PSF model (set/unset the pixel mask)
     psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "OR", PM_SOURCE_MASK_MARKED);
-    pmSourceFitSet (source, modelSet, true);
+    pmSourceFitSet (source, modelSet, PM_SOURCE_FIT_PSF);
     psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "AND", ~PM_SOURCE_MASK_MARKED);
 
@@ -210,5 +231,5 @@
     // fit EXT (not PSF) model (set/unset the pixel mask)
     psImageKeepCircle (source->mask, x, y, EXT->radiusTMP, "OR", PM_SOURCE_MASK_MARKED);
-    pmSourceFitModel (source, EXT, false);
+    pmSourceFitModel (source, EXT, PM_SOURCE_FIT_EXT);
     psImageKeepCircle (source->mask, x, y, EXT->radiusTMP, "AND", ~PM_SOURCE_MASK_MARKED);
 
@@ -216,5 +237,5 @@
 }
 
-bool psphotFitBlend (pmReadout *readout, pmSource *source) { 
+bool psphotFitBlend (pmReadout *readout, pmSource *source, pmPSF *psf) { 
 
     float x, y, dR;
@@ -222,5 +243,5 @@
     // if this source is not a possible blend, just fit as PSF
     if ((source->blends == NULL) || (source->mode & PM_SOURCE_MODE_SATSTAR)) {
-	bool status = psphotFitPSF (readout, source);
+	bool status = psphotFitPSF (readout, source, psf);
 	return status;
     }
@@ -273,6 +294,11 @@
     // fit PSF model (set/unset the pixel mask)
     psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "OR", PM_SOURCE_MASK_MARKED);
-    pmSourceFitSet (source, modelSet, true);
+    pmSourceFitSet (source, modelSet, PM_SOURCE_FIT_PSF);
     psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "AND", ~PM_SOURCE_MASK_MARKED);
+
+    // correct model chisq for flux trend
+    double chiTrend = psPolynomial1DEval (psf->ChiTrend, PSF->params->data.F32[1]);
+    PSF->chisqNorm = PSF->chisq / chiTrend;
+    // PSF->chisqNorm = PSF->chisq;
 
     // evaluate the blend objects, subtract if good, free otherwise
@@ -280,4 +306,9 @@
 	pmSource *src = sourceSet->data[i];
 	pmModel *model = modelSet->data[i];
+
+	// correct model chisq for flux trend
+	chiTrend = psPolynomial1DEval (psf->ChiTrend, model->params->data.F32[1]);
+	model->chisqNorm = model->chisq / chiTrend;
+	// model->chisqNorm = model->chisq;
 
 	// if we skip this one, free the corresponding blend entry model
Index: /trunk/psphot/src/psphotWeightBias.c
===================================================================
--- /trunk/psphot/src/psphotWeightBias.c	(revision 6851)
+++ /trunk/psphot/src/psphotWeightBias.c	(revision 6851)
@@ -0,0 +1,78 @@
+# include "psphot.h"
+
+// select objects fitted with PSF model
+// re-fit all of them with the non-poisson errors
+// only allow the normalization to vary 
+// XXX eventually, we should be able to do this with linear fitting...
+bool psphotWeightBias (pmReadout *readout, psArray *sources, psMetadata *config, pmPSF *psf) { 
+
+    int Nfit = 0;
+    bool status;
+    psF32 *PARp, *PARc;
+    psF32 x, y;
+
+    psTimerStart ("psphot");
+
+    // source analysis is done in S/N order (brightest first)
+    sources = psArraySort (sources, psphotSortBySN);
+    
+    // set fitting method to use non-poisson errors (local sky error)
+    pmSourceFitModelInit (15, 0.1, false);
+    
+    // option to limit analysis to a specific region
+    char *region = psMetadataLookupStr (&status, config, "ANALYSIS_REGION");
+    psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
+    if (psRegionIsNaN (AnalysisRegion)) psAbort ("psphot", "analysis region mis-defined");
+
+    FILE *f = fopen ("bias.dat", "w");
+    if (f == NULL) psAbort ("psphot", "can't open output file");
+
+    for (int i = 0; i < sources->n; i++) {
+	pmSource *source = sources->data[i];
+
+	// skip lower-quality objects
+	if (source->type != PM_SOURCE_TYPE_STAR) continue; 
+	if (source->mode &  PM_SOURCE_MODE_POOR) continue;
+	if (source->mode &  PM_SOURCE_MODE_FAIL) continue;
+	if (source->mode &  PM_SOURCE_MODE_PAIR) continue;
+	if (source->mode &  PM_SOURCE_MODE_BLEND) continue;
+	if (source->mode &  PM_SOURCE_MODE_SATSTAR) continue;
+
+	// if model is NULL, we don't have a starting guess
+	if (source->modelPSF == NULL) continue;
+
+	if (source->moments->x < AnalysisRegion.x0) continue;
+	if (source->moments->y < AnalysisRegion.y0) continue;
+	if (source->moments->x > AnalysisRegion.x1) continue;
+	if (source->moments->y > AnalysisRegion.y1) continue;
+
+	// replace object in image
+	pmModelAdd (source->pixels, source->mask, source->modelPSF, false, false);
+
+	// make a temporary model (we don't keep the result of this analysis)
+	pmModel *PSF = pmModelCopy (source->modelPSF);
+
+	// extend source radius as needed
+	psphotCheckRadiusPSF (readout, source, PSF);
+
+	x = PSF->params->data.F32[2];
+	y = PSF->params->data.F32[3];
+
+	// fit PSF model (set/unset the pixel mask)
+	psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "OR", PM_SOURCE_MASK_MARKED);
+	pmSourceFitModel (source, PSF, PM_SOURCE_FIT_NORM);
+	psImageKeepCircle (source->mask, x, y, PSF->radiusTMP, "AND", ~PM_SOURCE_MASK_MARKED);
+    
+	// re-subtract PSF for object, leave local sky
+	pmModelSub (source->pixels, source->mask, source->modelPSF, false, false);
+
+	PARp = source->modelPSF->params->data.F32;
+	PARc = PSF->params->data.F32;
+	fprintf (f, "%7.1f %7.1f %9.2f %9.2f %10.3f\n", PARp[2], PARp[3], PARp[1], PARc[1], source->moments->dSky);
+	Nfit ++;
+    }
+    fclose (f);
+    psLogMsg ("psphot", 3, "measure PSF weighting bias for %d objects: %f sec\n", Nfit, psTimerMark ("psphot"));
+
+    return (true);
+}