Index: /trunk/psphot/doc/psfmodel.txt
===================================================================
--- /trunk/psphot/doc/psfmodel.txt	(revision 9772)
+++ /trunk/psphot/doc/psfmodel.txt	(revision 9772)
@@ -0,0 +1,22 @@
+
+2006.10.27
+
+  I have been working to fix the PSF modeling in psphot.  The PSF
+  model consists of a flux model for an individual object using an
+  analytical model with a number of parameters.  For a collection of
+  PSF objects, the variation of the parameters as a function of
+  position are then themselves fitted with a model.  The PSF model for
+  a single object consists of a radial profile with a functional form
+  f(z) where a given value of z defines an elliptical coutour of the
+  form z = \frac{x^2}{2\sigma_x^2} + \frac{y^2}{2\sigma_y^2} +
+  \sigma_{xy}xy.
+
+
+
+  The term \sigma_{xy} is difficult to model as a simple function of x
+  and y (eg, a low-order polynomial).  
+
+  A better model can be constructed for
+  \frac{\sigma_{xy}}{(\sigma_x^{-2} + \sigma_y^{-2})^2}, which varies
+  like a^2 \sin 2\theta 
+
Index: /trunk/psphot/src/polyfitTest.c
===================================================================
--- /trunk/psphot/src/polyfitTest.c	(revision 9772)
+++ /trunk/psphot/src/polyfitTest.c	(revision 9772)
@@ -0,0 +1,93 @@
+# include "psphot.h"
+
+int main (int argc, char **argv) {
+
+    float X, Y, SX, SY, SXY;
+
+    if (argc != 3) {
+	fprintf (stderr, "USAGE: polyfitTest (input) (output)\n");
+	exit (2);
+    }
+
+    FILE *f = fopen (argv[1], "r");
+    if (f == NULL) exit (1);
+
+    psVector *x = psVectorAlloc (100, PS_TYPE_F32);
+    psVector *y = psVectorAlloc (100, PS_TYPE_F32);
+    psVector *sx = psVectorAlloc (100, PS_TYPE_F32);
+    psVector *sy = psVectorAlloc (100, PS_TYPE_F32);
+    psVector *sxy = psVectorAlloc (100, PS_TYPE_F32);
+
+    int i = 0;
+    while (fscanf (f, "%f %f %f %f %f %*f %*d", &X, &Y, &SX, &SY, &SXY) != EOF) {
+	x->data.F32[i] = X;
+	y->data.F32[i] = Y;
+	sx->data.F32[i] = SX;
+	sy->data.F32[i] = SY;
+	sxy->data.F32[i] = SXY;
+
+	psVectorExtend (x, 100, 1);
+	psVectorExtend (y, 100, 1);
+	psVectorExtend (sx, 100, 1);
+	psVectorExtend (sy, 100, 1);
+	psVectorExtend (sxy, 100, 1);
+
+	i++;
+    }    
+    fclose (f);
+
+    fprintf (stderr, "loaded %ld pts\n", x->n);
+
+    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
+    psVector *mask = psVectorAlloc (x->n, PS_TYPE_MASK);
+    mask->n = x->n;
+    psVectorInit (mask, 0);
+
+    psPolynomial2D *polySX = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 2, 2);
+    polySX->mask[2][1] = 1;
+    polySX->mask[2][2] = 1;
+    polySX->mask[1][2] = 1;
+
+    psVectorClipFitPolynomial2D (polySX, stats, mask, 0xff, sx, NULL, x, y);
+    fprintf (stderr, "stats: %f +/- %f\n", stats->robustMedian, stats->robustStdev);
+    fprintf (stderr, "SX: %g %g %g\n", polySX->coeff[0][0], polySX->coeff[1][0], polySX->coeff[0][1]);
+    fprintf (stderr, "SX: %g %g %g\n", polySX->coeff[0][2], polySX->coeff[1][1], polySX->coeff[2][0]);
+
+    psPolynomial2D *polySY = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 2, 2);
+    polySY->mask[2][1] = 1;
+    polySY->mask[2][2] = 1;
+    polySY->mask[1][2] = 1;
+
+    psVectorClipFitPolynomial2D (polySY, stats, mask, 0xff, sy, NULL, x, y);
+    fprintf (stderr, "stats: %f +/- %f\n", stats->robustMedian, stats->robustStdev);
+    fprintf (stderr, "SY: %g %g %g\n", polySY->coeff[0][0], polySY->coeff[1][0], polySY->coeff[0][1]);
+    fprintf (stderr, "SY: %g %g %g\n", polySY->coeff[0][2], polySY->coeff[1][1], polySY->coeff[2][0]);
+
+    psPolynomial2D *polySXY = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 2, 2);
+    polySXY->mask[2][1] = 1;
+    polySXY->mask[2][2] = 1;
+    polySXY->mask[1][2] = 1;
+
+    psVectorClipFitPolynomial2D (polySXY, stats, mask, 0xff, sxy, NULL, x, y);
+    fprintf (stderr, "stats: %f +/- %f\n", stats->robustMedian, stats->robustStdev);
+    fprintf (stderr, "SXY: %g %g %g\n", polySXY->coeff[0][0], polySXY->coeff[1][0], polySXY->coeff[0][1]);
+    fprintf (stderr, "SXY: %g %g %g\n", polySXY->coeff[0][2], polySXY->coeff[1][1], polySXY->coeff[2][0]);
+
+    psVector *sxf = psPolynomial2DEvalVector (polySX, x, y);
+    psVector *syf = psPolynomial2DEvalVector (polySY, x, y);
+    psVector *sxyf = psPolynomial2DEvalVector (polySXY, x, y);
+
+    f = fopen (argv[2], "w");
+    for (i = 0; i < x->n; i++) {
+	fprintf (f, "%f %f  %f %f %f  %f %f %f\n",
+		 x->data.F32[i], y->data.F32[i], 
+		 sx->data.F32[i], sy->data.F32[i], sxy->data.F32[i], 
+		 sxf->data.F32[i], syf->data.F32[i], sxyf->data.F32[i]);
+    }
+    fclose (f);
+    exit (0);
+}
+
+    // 0,0 1,0 2,0
+    // 0,1 1,1 2,1
+    // 0,2 1,2 2,2