Index: trunk/psModules/src/objects/pmPSF.c
===================================================================
--- trunk/psModules/src/objects/pmPSF.c	(revision 9730)
+++ trunk/psModules/src/objects/pmPSF.c	(revision 9770)
@@ -6,6 +6,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-10-24 22:55:05 $
+ *  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-10-28 20:23:51 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -69,14 +69,15 @@
     psFree (psf->ApTrend);
     psFree (psf->growth);
-    psFree (psf->params);
+    psFree (psf->params_NEW);
     return;
 }
 
-
-
 /*****************************************************************************
-pmPSFAlloc (type): allocate a pmPSF.
-    NOTE: a PSF always has 4 parameters fewer than the equivalent model.
-      This is because those 4 parameters are
+ pmPSFAlloc (type): allocate a pmPSF.
+ 
+ NOTE: PSF model parameters which are not modeled on an image are set to NULL in psf->params.
+ 
+ These are normally:
+ 
  X-center
  Y-center
@@ -119,17 +120,24 @@
         return(NULL);
     }
-
-    psf->params = psArrayAlloc(Nparams - 4);
-
-    // the order of the PSF parameter (X,Y) fits is determined by the
-    // psfTrendMask polynomial (user-specified as in the recipe). the
-    // masks of psfTrendMask are applied to each parameter.
-    // if psfTrendMask is NULL, these polynomials are not allocated.
-    // in this case, the user must set them by hand (as in pmPSFfromMD)
-    // XXX should we drop the hard-wired '4' above and use NULL to identify
-    // the parameters which are not fitted.  these could be selected by
-    // testing for the value of PM_PAR_XPOS, etc.
+    psf->params_NEW = psArrayAlloc(Nparams);
+
+    // the order of the PSF parameter (X,Y) fits is determined by the psfTrendMask polynomial
+    // (user-specified as in the recipe). the masks of psfTrendMask are applied to each
+    // parameter.  if psfTrendMask is NULL, these polynomials are not allocated.  in this case,
+    // the user must set them by hand (as in pmPSFfromMD).  the parameters which are not fitted
+    // are left as NULL.  these are selected by testing for them by the named values (
+    // PM_PAR_XPOS, etc)
+
     if (psfTrendMask) {
-        for (int i = 0; i < psf->params->n; i++) {
+        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;
+
             psPolynomial2D *param = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, psfTrendMask->nX, psfTrendMask->nY);
             for (int nx = 0; nx < param->nX + 1; nx++) {
@@ -138,5 +146,5 @@
                 }
             }
-            psf->params->data[i] = param;
+            psf->params_NEW->data[i] = param;
         }
     }
@@ -166,4 +174,5 @@
     psVector *dz = psVectorAlloc (models->n, PS_TYPE_F64);
 
+    // construct the x,y terms
     for (int i = 0; i < models->n; i++) {
         pmModel *model = models->data[i];
@@ -171,8 +180,7 @@
             continue;
 
-        // XXX EAM : this is fragile: x and y must be stored consistently at 2,3
-        // note that the data is provided in the F64 array
-        x->data.F64[i] = model->params->data.F32[2];
-        y->data.F64[i] = model->params->data.F32[3];
+        // use F64 for polynomial fitting
+        x->data.F64[i] = model->params->data.F32[PM_PAR_XPOS];
+        y->data.F64[i] = model->params->data.F32[PM_PAR_YPOS];
     }
 
@@ -183,23 +191,60 @@
     psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
 
-    for (int i = 0; i < psf->params->n; i++) {
+    // skip the unfitted parameters (X, Y, Io, Sky)
+    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;
+
+        // select the per-object fitted data for this parameter
         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
+            z->data.F64[j] = model->params->data.F32[i];
+            dz->data.F64[j] = 1; // use the model-fitted error? or S/N?
+
+            // for SXY, we actually fit SXY / (SXX^-2  + SYY^-2)
+            if (i == PM_PAR_SXY) {
+                z->data.F64[j] = pmPSF_SXYfromModel (model->params->data.F32);
+            }
         }
-
-        // XXX EAM : this is the expected API (cycle 7? cycle 8?)
-        psf->params->data[i] = psVectorClipFitPolynomial2D(psf->params->data[i], stats, mask, 0xff, z, dz, x, y);
-
-        // XXX EAM : drop this when above is implemented...
-        // psf->params->data[i] = RobustFit2D (psf->params->data[i], mask, x, y, z, dz);
-
+        psf->params_NEW->data[i] = psVectorClipFitPolynomial2D(psf->params_NEW->data[i], stats, mask, 0xff, z, dz, x, y);
         // psTrace ("psphot.psftest", 3, "keeping %d of %d PSF candidates (PSF param %d)\n", Nkeep, mask->n, i);
-        // psPolynomial2DDump (psf->params->data[i]);
+
+        // XXX Test output
+        psPolynomial2D *poly = psf->params_NEW->data[i];
+        fprintf (stderr, "stats: %f +/- %f\n", stats->robustMedian, stats->robustStdev);
+        fprintf (stderr, "PO: %g %g %g\n", poly->coeff[0][0], poly->coeff[1][0], poly->coeff[0][1]);
+        fprintf (stderr, "PO: %g %g %g\n", poly->coeff[0][2], poly->coeff[1][1], poly->coeff[2][0]);
+    }
+
+    // XXX test dump of star parameters vs position (compare with fitted values)
+    {
+        FILE *f = fopen ("params.dat", "w");
+
+        for (int j = 0; j < models->n; j++)
+        {
+            pmModel *model = models->data[j];
+            if (model == NULL)
+                continue;
+
+            pmModel *modelPSF = pmModelFromPSF (model, psf);
+
+            fprintf (f, "%f %f : ", model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
+
+            for (int i = 0; i < psf->params_NEW->n; i++) {
+                if (psf->params_NEW->data[i] == NULL)
+                    continue;
+                fprintf (f, "%f %f : ", model->params->data.F32[i], modelPSF->params->data.F32[i]);
+            }
+            fprintf (f, "%f %d\n", model->chisq, model->nIter);
+        }
+        fclose (f);
     }
 
@@ -211,6 +256,4 @@
     return (true);
 }
-
-
 
 /*****************************************************************************
@@ -250,4 +293,33 @@
     }
     return true;
+}
+
+// the PSF models the \sigma_{xy} variation of the elliptical contour as a function of position in the image with a
+// polynomial.  an individual object has a contour of the form (x^2/2sx^2) + (y^2/2sy^2) + sxy*x*y
+// these are the values of the model->params.  the psf->params term for sxy is actually fitted
+// to sxy/(sxx^-2 + syy^-2)^2
+
+// input: model->param, output: psf->param[PM_PAR_SXY]
+double pmPSF_SXYfromModel (psF32 *modelPar)
+{
+
+    double SXX = modelPar[PM_PAR_SXX];
+    double SYY = modelPar[PM_PAR_SYY];
+    double SXY = modelPar[PM_PAR_SXY];
+
+    double par = SXY / PS_SQR(1.0 / PS_SQR(SXX) + 1.0 / PS_SQR(SYY));
+    return (par);
+}
+
+// input: fitted psf->param, output: model->param[PM_PAR_SXY]
+double pmPSF_SXYtoModel (psF32 *fittedPar)
+{
+
+    double SXX = fittedPar[PM_PAR_SXX];
+    double SYY = fittedPar[PM_PAR_SYY];
+    double fit = fittedPar[PM_PAR_SXY];
+
+    double SXY = fit * PS_SQR(1.0 / PS_SQR(SXX) + 1.0 / PS_SQR(SYY));
+    return SXY;
 }