Index: trunk/psphot/src/psphotEnsemblePSF.c
===================================================================
--- trunk/psphot/src/psphotEnsemblePSF.c	(revision 5980)
+++ trunk/psphot/src/psphotEnsemblePSF.c	(revision 5993)
@@ -1,49 +1,3 @@
 # include "psphot.h"
-
-psPolynomial2D *psImageBicubeFit (psImage *image, int x, int y) {
-
-    int ix = x - image->col0;
-    int iy = y - image->row0;
-
-    psF32 *Fm = &image->data.F32[iy - 1][ix];
-    psF32 *Fo = &image->data.F32[iy + 0][ix];
-    psF32 *Fp = &image->data.F32[iy + 1][ix];
-
-    double Fxm = Fm[-1] + Fo[-1] + Fp[-1];
-    double Fxp = Fm[+1] + Fo[+1] + Fp[+1];
-    double Fym = Fm[-1] + Fm[+0] + Fm[+1];
-    double Fyp = Fp[-1] + Fp[+0] + Fp[+1];
-    double Foo = Fym + Fyp + Fo[-1] + Fo[+0] + Fo[+1];
-
-    psPolynomial2D *poly = psPolynomial2DAlloc (2, 2, PS_POLYNOMIAL_ORD);
-    poly->mask[2][2] = 1;
-    poly->mask[1][2] = 1;
-    poly->mask[2][1] = 1;
-
-    poly->coeff[0][0] = Foo*(5.0/9.0) - (Fxp + Fxm)/3.0 - (Fyp + Fym)/3.0 ;
-
-    poly->coeff[1][0] = (Fxp - Fxm)/6.0;
-    poly->coeff[0][1] = (Fyp - Fym)/6.0;
-    
-    poly->coeff[2][0] = (Fxp + Fxm)/2.0 - Foo/3.0;
-    poly->coeff[0][2] = (Fyp + Fym)/2.0 - Foo/3.0;
-    
-    poly->coeff[1][1] = (Fp[+1] + Fm[-1] - Fm[+1] - Fp[-1])/4.0;
-    
-    return (poly);
-}
-
-psPlane psImageBicubeMin (psPolynomial2D *poly) {
-
-    psPlane min;
-
-    min.xErr = min.yErr = 0;
-
-    double det = 4*poly->coeff[2][0]*poly->coeff[0][2] - PS_SQR(poly->coeff[1][1]);
-
-    min.x = (poly->coeff[1][1]*poly->coeff[0][1] - 2*poly->coeff[0][2]*poly->coeff[1][0]) / det;
-    min.y = (poly->coeff[1][1]*poly->coeff[1][0] - 2*poly->coeff[2][0]*poly->coeff[0][1]) / det;
-    return (min);
-}
 
 bool psphotEnsemblePSF (eamReadout *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) { 
@@ -107,10 +61,10 @@
 	}
 
-	// XXX EAM : add option to use FLT or PSF form
+	// XXX EAM : add option to use EXT or PSF form
 	// use the source moments, etc to guess basic model parameters
-	pmModel *modelFLT = pmSourceModelGuess (inSource, psf->type); 
+	pmModel *modelEXT = pmSourceModelGuess (inSource, psf->type); 
 	if (inSource->mode &  PM_SOURCE_SATSTAR) {
-	    modelFLT->params->data.F32[2] = inSource->moments->x;
-	    modelFLT->params->data.F32[3] = inSource->moments->y;
+	    modelEXT->params->data.F32[2] = inSource->moments->x;
+	    modelEXT->params->data.F32[3] = inSource->moments->y;
 	} else {
 	    // peak-up on peak (for non-sat objects)
@@ -123,10 +77,10 @@
 
 	    psTrace ("psphotEnsemblePSF", 5, "peak coord: %f %f -> %f %f\n", 
-		     modelFLT->params->data.F32[2], modelFLT->params->data.F32[3], min.x + ix, min.y + iy);
+		     modelEXT->params->data.F32[2], modelEXT->params->data.F32[3], min.x + ix, min.y + iy);
 	    
 	    // if min point is too deviant, keep the old value
 	    if ((fabs(min.x) < 1.5) && (fabs(min.y) < 1.5)) {
-		modelFLT->params->data.F32[2] = min.x + ix;
-		modelFLT->params->data.F32[3] = min.y + iy;
+		modelEXT->params->data.F32[2] = min.x + ix;
+		modelEXT->params->data.F32[3] = min.y + iy;
 	    }
 	    psFree (bicube);
@@ -134,6 +88,6 @@
 
 	// set PSF parameters for this model
-	pmModel *model = pmModelFromPSF (modelFLT, psf);
-	psFree (modelFLT);
+	pmModel *model = pmModelFromPSF (modelEXT, psf);
+	psFree (modelEXT);
 
 	// save the original coords
@@ -185,9 +139,9 @@
 
 	// diagonal element (auto-cross-product)
-	f = psphotCrossProduct (Mi, Mi);
+	f = pmSourceCrossProduct (Mi, Mi);
 	psSparseMatrixElement (sparse, i, i, f);
 
 	// find the image x model value
-	f = psphotCrossProduct (Fi, Mi);
+	f = pmSourceCrossProduct (Fi, Mi);
 	psSparseVectorElement (sparse, i, f);
 
@@ -203,5 +157,5 @@
 	    
 	    // got an overlap; calculate cross-product and add to output array
-	    f = psphotCrossProduct (Mi, Mj);
+	    f = pmSourceCrossProduct (Mi, Mj);
 	    psSparseMatrixElement (sparse, j, i, f); 
 	}
@@ -235,42 +189,2 @@
     return true;
 }
-
-
-float psphotCrossProduct (pmSource *Mi, pmSource *Mj) {
-
-    int Xs, Xe, Ys, Ye;
-    int xi, xj, yi, yj;
-    int xIs, xJs, yIs, yJs;
-    int xIe, yIe;
-    float flux;
-
-    psImage *Pi = Mi->pixels;
-    psImage *Pj = Mj->pixels;
-
-    psImage *Ti = Mi->mask;
-    psImage *Tj = Mj->mask;
-
-    Xs = PS_MAX (Pi->col0, Pj->col0);
-    Xe = PS_MIN (Pi->col0 + Pi->numCols, Pj->col0 + Pj->numCols);
-    
-    Ys = PS_MAX (Pi->row0, Pj->row0);
-    Ye = PS_MIN (Pi->row0 + Pi->numRows, Pj->row0 + Pj->numRows);
-    
-    xIs = Xs - Pi->col0;
-    xJs = Xs - Pj->col0;
-    yIs = Ys - Pi->row0;
-    yJs = Ys - Pj->row0;
-
-    xIe = Xe - Pi->col0;
-    yIe = Ye - Pi->row0;
-
-    flux = 0;
-    for (yi = yIs, yj = yJs; yi < yIe; yi++, yj++) {
-	for (xi = xIs, xj = xJs; xi < xIe; xi++, xj++) {
-	    if (Ti->data.U8[yi][xi]) continue;
-	    if (Tj->data.U8[yj][xj]) continue;
-	    flux += Pi->data.F32[yi][xi] * Pj->data.F32[yj][xj];
-	}
-    }
-    return (flux);
-}