Index: trunk/psphot/src/pspsf.c
===================================================================
--- trunk/psphot/src/pspsf.c	(revision 4954)
+++ trunk/psphot/src/pspsf.c	(revision 4977)
@@ -1,39 +1,12 @@
 # include "psphot.h"
 
-static void pmPSFFree (pmPSF *psf) {
-
-  if (psf == NULL) return;
-
-  psFree (psf->params);
-  return;
-}
-
-// a PSF always has 4 parameters fewer than the equivalent model
-pmPSF *pmPSFAlloc (pmModelType type) {
-
-    int Nparams;
-
-    pmPSF *psf = (pmPSF *) psAlloc(sizeof(pmPSF));
-
-    psf->type = type;
-    psf->chisq = 0.0;
-
-    Nparams = pmModelParameterCount (type);
-    if (!Nparams) {
-	psError(PS_ERR_UNKNOWN, true, "Undefined pmModelType");
-	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);
-    }
-
-    psMemSetDeallocator(psf, (psFreeFunc) pmPSFFree);
-    return(psf);
-}
-
-static void pmPSF_TestFree (pmPSF_Test *test) {
+// ********  pmPSFtry functions  **************************************************
+// * pmPSFtry holds a single pmPSF model test, with the input sources, the freely
+// * fitted version of the model, the pmPSF fit to the fitted model parameters, 
+// * and the PSF fits to the source. It also includes the statistics from the 
+// * fits, both the individual sources, and the collection
+
+// free a pmPSFtry structure
+static void pmPSFtryFree (pmPSFtry *test) {
 
   if (test == NULL) return;
@@ -49,7 +22,8 @@
 }
 
-pmPSF_Test *pmPSF_TestAlloc (psArray *sources, char *modelName) {
-
-    pmPSF_Test *test = (pmPSF_Test *) psAlloc(sizeof(pmPSF_Test));
+// allocate a pmPSFtry based on the desired sources and the model (identified by name)
+pmPSFtry *pmPSFtryAlloc (psArray *sources, char *modelName) {
+
+    pmPSFtry *test = (pmPSFtry *) psAlloc(sizeof(pmPSFtry));
 
     // XXX probably need to increment ref counter
@@ -74,14 +48,23 @@
     }	
 
-    psMemSetDeallocator(test, (psFreeFunc) pmPSF_TestFree);
+    psMemSetDeallocator(test, (psFreeFunc) pmPSFtryFree);
     return (test);
 }
 
-// test->mask values indicate reason source was rejected:
+// build a pmPSFtry for the given model:
+// - fit each source with the free-floating model
+// - construct the pmPSF from the collection of models
+// - fit each source with the PSF-parameter models 
+// - measure the pmPSF quality metric (dApResid)
+
+// sources used in for pmPSFtry may be masked by the analysis
+// mask values indicate the reason the source was rejected:
 // 1: outlier in psf polynomial fit
 // 2: flt model failed to converge
 // 3: psf model failed to converge
 // 4: invalid source photometry
-pmPSF_Test *pmPSF_TestModel (psArray *sources, char *modelName, float RADIUS) 
+// XXX EAM : use an enum for these values?
+
+pmPSFtry *pmPSFtryModel (psArray *sources, char *modelName, float RADIUS) 
 {
     bool status;
@@ -93,12 +76,12 @@
     int Npsf = 0;
 
-    pmPSF_Test *test = pmPSF_TestAlloc (sources, modelName);
+    pmPSFtry *try = pmPSFtryAlloc (sources, modelName);
 
     // stage 1:  fit an independent model (freeModel) to all sources
     psTimerStart ("fit");
-    for (int i = 0; i < test->sources->n; i++) {
-
-	pmSource *source = test->sources->data[i];
-	pmModel  *model  = pmSourceModelGuess (source, test->modelType); 
+    for (int i = 0; i < try->sources->n; i++) {
+
+	pmSource *source = try->sources->data[i];
+	pmModel  *model  = pmSourceModelGuess (source, try->modelType); 
 	x = source->peak->x;
 	y = source->peak->y;
@@ -112,31 +95,31 @@
 	// exclude the poor fits
 	if (!status) {
-	  test->mask->data.U8[i] = 2;
+	  try->mask->data.U8[i] = 2;
 	  psFree (model);
 	  continue;
 	}
-	test->modelFLT->data[i] = model;
+	try->modelFLT->data[i] = model;
 	Nflt ++;
     }
-    psLogMsg ("psphot.psftest", 4, "fit flt:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
-    psTrace ("psphot.psftest", 3, "keeping %d of %d PSF candidates (FLT)\n", Nflt, sources->n);
+    psLogMsg ("psphot.psftry", 4, "fit flt:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
+    psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (FLT)\n", Nflt, sources->n);
 
     // make this optional? 
-    // DumpModelFits (test->modelFLT, "modelsFLT.dat");
+    // DumpModelFits (try->modelFLT, "modelsFLT.dat");
 
     // stage 2: construct a psf (pmPSF) from this collection of model fits
-    pmPSFFromModels (test->psf, test->modelFLT, test->mask);
+    pmPSFFromModels (try->psf, try->modelFLT, try->mask);
 
     // stage 3: refit with fixed shape parameters
     psTimerStart ("fit");
-    for (int i = 0; i < test->sources->n; i++) {
+    for (int i = 0; i < try->sources->n; i++) {
 	// masked for: bad model fit, outlier in parameters
-	if (test->mask->data.U8[i]) continue;
-
-	pmSource *source = test->sources->data[i];
-	pmModel  *modelFLT = test->modelFLT->data[i];
+	if (try->mask->data.U8[i]) continue;
+
+	pmSource *source = try->sources->data[i];
+	pmModel  *modelFLT = try->modelFLT->data[i];
 
 	// set shape for this model based on PSF
-	pmModel *modelPSF = pmModelFromPSF (modelFLT, test->psf); 
+	pmModel *modelPSF = pmModelFromPSF (modelFLT, try->psf); 
 	x = source->peak->x;
 	y = source->peak->y;
@@ -147,5 +130,5 @@
 	// skip poor fits
 	if (!status) {
-	    test->mask->data.U8[i] = 3;
+	    try->mask->data.U8[i] = 3;
 	    psFree (modelPSF);
 	    goto next_source;
@@ -153,5 +136,5 @@
 
 	// otherwise, save the resulting model
-	test->modelPSF->data[i] = modelPSF;
+	try->modelPSF->data[i] = modelPSF;
 
 	// XXX : use a different aperture radius from the fit radius?
@@ -159,10 +142,10 @@
 	// XXX : pass 'source' as input?
 	if (!pmSourcePhotometry (&fitMag, &obsMag, modelPSF, source->pixels, source->mask)) {
-	    test->mask->data.U8[i] = 4;
+	    try->mask->data.U8[i] = 4;
 	    goto next_source;
 	}	    
 
-	test->metric->data.F64[i] = obsMag - fitMag;
-	test->fitMag->data.F64[i] = fitMag;
+	try->metric->data.F64[i] = obsMag - fitMag;
+	try->fitMag->data.F64[i] = fitMag;
 	Npsf ++;
 
@@ -171,86 +154,25 @@
 
     }
-    psLogMsg ("psphot.psftest", 4, "fit psf:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
-    psTrace ("psphot.psftest", 3, "keeping %d of %d PSF candidates (PSF)\n", Npsf, sources->n);
+    psLogMsg ("psphot.psftry", 4, "fit psf:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
+    psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (PSF)\n", Npsf, sources->n);
 
     // make this optional
-    // DumpModelFits (test->modelPSF, "modelsPSF.dat");
-
-    // XXX this function wants aperture radius from pmSourcePhotometry
-    pmPSFMetricModel (test, RADIUS);
-    psLogMsg ("psphot.pspsf", 3, "test model %s, ap-fit: %f +/- %f, sky bias: %f\n", 
-	      modelName, test->ApResid, test->dApResid, test->skyBias);
-
-    return (test);
-}
-
-// input: an array of pmModels, pre-allocated psf
-// some of the array entries may be NULL, ignore them
-bool pmPSFFromModels (pmPSF *psf, psArray *models, psVector *mask) {
-
-    // construct the fit vectors from the collection of objects
-    // use the mask to ignore missing fits 
-    psVector *x  = psVectorAlloc (models->n, PS_TYPE_F64);
-    psVector *y  = psVectorAlloc (models->n, PS_TYPE_F64);
-    psVector *z  = psVectorAlloc (models->n, PS_TYPE_F64);
-    psVector *dz = psVectorAlloc (models->n, PS_TYPE_F64);
-
-    for (int i = 0; i < models->n; i++) {
-	pmModel *model = models->data[i];
-	if (model == NULL) continue;
-
-	// XXX EAM : this is fragile: x and y must be stored consistently at 2,3
-	x->data.F64[i] = model->params->data.F32[2];
-	y->data.F64[i] = model->params->data.F32[3];
-    }
-
-    // we are doing a robust fit.  after each pass, we drop points which are 
-    // more deviant than three sigma. 
-    // mask is currently updated for each pass. this is inconsistent?
-
-    for (int i = 0; i < psf->params->n; i++) {
-	for (int j = 0; j < models->n; j++) {
-	    pmModel *model = models->data[j];
-	    if (model == NULL) continue;
-	    z->data.F64[j] = model->params->data.F32[i + 4];
-	    dz->data.F64[j] = 1;
-	    // XXX EAM : need to use actual errors?
-	    // XXX EAM : this is fragile: psf(Nparams) = model(Nparams) - 4
-	}
-
-	psf->params->data[i] = RobustFit2D (psf->params->data[i], mask, x, y, z, dz);
-	// psTrace ("psphot.psftest", 3, "keeping %d of %d PSF candidates (PSF param %d)\n", Nkeep, mask->n, i);
-	// psPolynomial2DDump (psf->params->data[i]);
-    }
-    
-    psFree (x);
-    psFree (y);
-    psFree (z);
-    psFree (dz);
-    return (true);
-}
-
-pmModel *pmModelFromPSF (pmModel *modelFLT, pmPSF *psf) {
-
-    // need to define the relationship between the modelFLT and modelPSF ?
-    
-    // find function used to set the model parameters
-    pmModelFromPSFFunc modelFromPSFFunc = pmModelFromPSFFunc_GetFunction (psf->type);
-
-    // do we need a different model for the PSF vs FLT version?
-    pmModel *modelPSF = pmModelAlloc (psf->type);
-
-    // set model parameters for this source based on PSF information
-    modelFromPSFFunc (modelPSF, modelFLT, psf);
-
-    return (modelPSF);
-}
-
-bool pmPSFMetricModel (pmPSF_Test *test, float RADIUS) {
+    // DumpModelFits (try->modelPSF, "modelsPSF.dat");
+
+    // XXX this function wants aperture radius for pmSourcePhotometry
+    pmPSFtryMetric (try, RADIUS);
+    psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f, sky bias: %f\n", 
+	      modelName, try->ApResid, try->dApResid, try->skyBias);
+
+    return (try);
+}
+
+
+bool pmPSFtryMetric (pmPSFtry *try, float RADIUS) {
 
   float dBin;
   int   nKeep, nSkip;
 
-  // the measured (aperture - fit) magnitudes (dA == test->metric)
+  // the measured (aperture - fit) magnitudes (dA == try->metric)
   //   depend on both the true ap-fit (dAo) and the bias in the sky measurement:
   //     dA = dAo + dsky/flux
@@ -262,13 +184,13 @@
 
   // rflux = ten(0.4*fitMag);
-  psVector *rflux = psVectorAlloc (test->sources->n, PS_TYPE_F64);
-  for (int i = 0; i < test->sources->n; i++) {
-    if (test->mask->data.U8[i]) continue;
-    rflux->data.F64[i] = pow(10.0, 0.4*test->fitMag->data.F64[i]);
+  psVector *rflux = psVectorAlloc (try->sources->n, PS_TYPE_F64);
+  for (int i = 0; i < try->sources->n; i++) {
+    if (try->mask->data.U8[i]) continue;
+    rflux->data.F64[i] = pow(10.0, 0.4*try->fitMag->data.F64[i]);
   }
 
   // find min and max of (1/flux):
   psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
-  psVectorStats (stats, rflux, NULL, test->mask, 0xff);
+  psVectorStats (stats, rflux, NULL, try->mask, 0xff);
   
   // build binned versions of rflux, metric
@@ -284,14 +206,14 @@
   for (int i = 0; i < daBin->n; i++) {
 
-    psVector *tmp = psVectorAlloc (test->sources->n, PS_TYPE_F64);
+    psVector *tmp = psVectorAlloc (try->sources->n, PS_TYPE_F64);
     tmp->n = 0;
 
     // accumulate data within bin range
-    for (int j = 0; j < test->sources->n; j++) {
+    for (int j = 0; j < try->sources->n; j++) {
       // masked for: bad model fit, outlier in parameters
-      if (test->mask->data.U8[j]) continue;
+      if (try->mask->data.U8[j]) continue;
     
       // skip points with extreme dA values
-      if (fabs(test->metric->data.F64[j]) > 0.5) continue;
+      if (fabs(try->metric->data.F64[j]) > 0.5) continue;
 
       // skip points outside of this bin
@@ -299,5 +221,5 @@
       if (rflux->data.F64[j] > rfBin->data.F64[i] + 0.5*dBin) continue;
 
-      tmp->data.F64[tmp->n] = test->metric->data.F64[j];
+      tmp->data.F64[tmp->n] = try->metric->data.F64[j];
       tmp->n ++;
     }
@@ -342,8 +264,13 @@
 
   // linear fit to rfBin, daBin
-  psPolynomial1D *poly = Polynomial1DAlloc (1, PS_POLYNOMIAL_ORD);
+  psPolynomial1D *poly = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 1);
+
+  // XXX EAM : this is the intended API (cycle 7? cycle 8?) 
+  // poly = psVectorFitPolynomial1D (poly, maskB, 1, daBin, NULL, rfBin);
+
+  // XXX EAM : replace this when the above version is implemented
   poly = VectorFitPolynomial1DOrd_EAM (poly, maskB, rfBin, daBin, NULL);
 
-  psVector *daBinFit = Polynomial1DEvalVector_EAM (poly, rfBin);
+  psVector *daBinFit = psPolynomial1DEvalVector (poly, rfBin);
   psVector *daResid  = (psVector *) psBinaryOp (NULL, (void *) daBin, "-", (void *) daBinFit);
 
@@ -351,7 +278,7 @@
   stats = psVectorStats (stats, daResid, NULL, maskB, 1);
 
-  test->ApResid = poly->coeff[0];
-  test->dApResid = stats->clippedStdev;
-  test->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
+  try->ApResid = poly->coeff[0];
+  try->dApResid = stats->clippedStdev;
+  try->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
 
   psFree (rflux);