Index: trunk/psModules/src/objects/pmPSFtry.c
===================================================================
--- trunk/psModules/src/objects/pmPSFtry.c	(revision 13034)
+++ trunk/psModules/src/objects/pmPSFtry.c	(revision 13064)
@@ -5,6 +5,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.36 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-26 01:20:29 $
+ *  @version $Revision: 1.37 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-04-27 22:14:08 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -352,4 +352,9 @@
         x->data.F64[i] = source->modelEXT->params->data.F32[PM_PAR_XPOS];
         y->data.F64[i] = source->modelEXT->params->data.F32[PM_PAR_YPOS];
+
+	// weight by the error on the source flux
+	if (dz) {
+	    dz->data.F64[i] = source->modelEXT->dparams->data.F32[PM_PAR_I0];
+	}
     }
 
@@ -360,14 +365,50 @@
     psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
 
-    // skip the unfitted parameters (X, Y, Io, Sky)
+    // The shape parameters (SXX, SXY, SYY) are strongly coupled.  We have to handle them very
+    // carefully.  First, we convert them to the Ellipse Polarization terms (E0, E1, E2) for
+    // each source and fit this set of parameters.  These values are less tightly coupled, but
+    // are still inter-related.  The fitted values do a good job of constraining the major axis
+    // and the position angle, but the minor axis is weakly measured.  When we apply the PSF
+    // model to construct a source model, we convert the fitted values of E0,E1,E2 to the shape
+    // parameters, with the constraint that the minor axis must be greater than a minimum
+    // threshold.
+
+    // convert the measured source shape paramters to polarization terms
+    psVector *e0   = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
+    psVector *e1   = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
+    psVector *e2   = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
+    for (int i = 0; i < psfTry->sources->n; i++) {
+	pmSource *source = psfTry->sources->data[i];
+	if (source->modelEXT == NULL) continue;
+
+	psEllipsePol pol = pmPSF_ModelToFit (source->modelEXT->params->data.F32);
+
+	e0->data.F64[i] = pol.e0;
+	e1->data.F64[i] = pol.e1;
+	e2->data.F64[i] = pol.e2;
+    }
+
+    // we run 'clipIter' cycles clipping in each of x and y, with only one iteration each.
+    // This way, the parameters masked by one of the fits will be applied to the others
+    for (int i = 0; i < stats->clipIter; i++) {
+        psVectorClipFitPolynomial2D (psf->params_NEW->data[PM_PAR_E0], stats, psfTry->mask, 0xff, e0, dz, x, y);
+        psVectorClipFitPolynomial2D (psf->params_NEW->data[PM_PAR_E1], stats, psfTry->mask, 0xff, e1, dz, x, y);
+        psVectorClipFitPolynomial2D (psf->params_NEW->data[PM_PAR_E2], stats, psfTry->mask, 0xff, e2, dz, x, y);
+    }
+
+    // skip the unfitted parameters (X, Y, Io, Sky) and the shape parameters (SXX, SYY, SXY)
     for (int i = 0; i < psf->params_NEW->n; i++) {
-        if (i == PM_PAR_SKY)
-            continue;
-        if (i == PM_PAR_I0)
-            continue;
-        if (i == PM_PAR_XPOS)
-            continue;
-        if (i == PM_PAR_YPOS)
-            continue;
+	switch (i) {
+	  case PM_PAR_SKY:
+	  case PM_PAR_I0:
+	  case PM_PAR_XPOS:
+	  case PM_PAR_YPOS:            
+	  case PM_PAR_SXX:            
+	  case PM_PAR_SYY:            
+	  case PM_PAR_SXY:            
+	    continue;
+	  default:
+	    break;
+	}
 
         // select the per-object fitted data for this parameter
@@ -376,19 +417,10 @@
             if (source->modelEXT == NULL)
                 continue;
-
             z->data.F64[j] = source->modelEXT->params->data.F32[i];
-            if (applyWeight) {
-                dz->data.F64[j] = source->modelEXT->dparams->data.F32[i];
-            }
-
-            // for SXY, we actually fit SXY / (SXX^-2  + SYY^-2)
-            if (i == PM_PAR_SXY) {
-                z->data.F64[j] = pmPSF_SXYfromModel (source->modelEXT->params->data.F32);
-            }
         }
 
         // fit the collection of measured parameters to the PSF 2D model
         // the mask is carried from previous steps and updated with this operation
-        // the weight is either the parameter error or NULL, depending on 'applyWeights'
+        // the weight is either the flux error or NULL, depending on 'applyWeights'
         if (!psVectorClipFitPolynomial2D(psf->params_NEW->data[i], stats, psfTry->mask, 0xff, z, NULL, x, y)) {
             psError(PS_ERR_UNKNOWN, false, "failed to build psf model for parameter %d", i);