Index: /trunk/Ohana/src/relastro/src/FitChip.c
===================================================================
--- /trunk/Ohana/src/relastro/src/FitChip.c	(revision 17208)
+++ /trunk/Ohana/src/relastro/src/FitChip.c	(revision 17209)
@@ -1,9 +1,21 @@
 # include "relastro.h"
-# define SCATTER_MAX_ERROR 0.05
-// XXX make this a user parameters
+
+// XXX make these user parameters
+# define FIT_CHIP_MAX_ERROR 0.05
+# define FIT_CHIP_NITER     3
+# define FIT_CHIP_NSIGMA    3.0
+
+// XXX we should test if the fit is sufficiently constrained across the chip, or if the
+// new positions deviate too much from the old positions
+
+// XXX add visualization tools: per chip residual plots
+
+// XXX save the fit[0].Npts value in the image table?
+
+// XXX save measurements of the fit quality (scatter, chisq) in the image table
 
 void FitChip (StarData *raw, StarData *ref, int Nmatch, Coords *coords) {
 
-  int i, Nscatter, Niter;
+  int i, Nscatter, Niter, skip;
   CoordFit *fit;
   double dL, dM, dR, dRmax, *values;
@@ -11,9 +23,10 @@
   ALLOCATE (values, double, Nmatch);
 
-  for (Niter = 0; Niter < 3; Niter ++) {
+  for (Niter = 0; Niter < FIT_CHIP_NITER; Niter ++) {
+
     // measure the scatter distribution (use only the bright end detections)
     for (i = Nscatter = 0; i < Nmatch; i++) {
       if (raw[i].mask) continue;
-      if (raw[i].dMag > SCATTER_MAX_ERROR) continue;
+      if (raw[i].dMag > FIT_CHIP_MAX_ERROR) continue;
 
       dL = raw[i].L - ref[i].L;
@@ -27,6 +40,7 @@
     // for a 2D Gaussian, 40% of the points are within 1 sigma; dRmax is ~ 3 sigma
     dsort (values, Nscatter);
-    dRmax = 3.0*values[(int)(0.40*Nscatter)];
+    dRmax = FIT_CHIP_NSIGMA*values[(int)(0.40*Nscatter)];
 
+    // fit the requested order polynomial
     if (CHIPORDER > 0) {
       coords[0].Npolyterms = CHIPORDER;
@@ -34,30 +48,44 @@
     fit = fit_init (coords[0].Npolyterms);
 
+    // generate the fit matches
     for (i = 0; i < Nmatch; i++) {
       if (raw[i].mask) continue;
 
-      // require radius of XXX arcsec
+      // only keep objects within dRmax
       dL = raw[i].L - ref[i].L;
       dM = raw[i].M - ref[i].M;
       dR = hypot (dL, dM);
-
-      // XXX ??? why are we not applying the mask?
-      // the detection is not included in the fit
-      if (dR > dRmax) {
-	continue;
-	raw[i].mask = TRUE;
-      }
+      if (dR > dRmax) continue;
     
       fit_add (fit, raw[i].X, raw[i].Y, ref[i].L, ref[i].M, raw[i].dPos);
     }
-    fprintf (stderr, "scatter limit: %f based on %d detections; using %d of %d for fit\n",
-	     dRmax, Nscatter, fit[0].Npts, Nmatch);
 
-    // XXX this limit should depend on the order of the fit (see psastro)
-    if (fit[0].Npts < 25) {
+    // check if the fit has enough data points for the polynomial order
+    skip = FALSE;
+    switch (coords[0].Npolyterms) {
+      case 0:
+      case 1:
+	skip = (fit[0].Npts < 8);
+	break;
+      case 2:
+	skip = (fit[0].Npts < 11);
+	break;
+      case 3:
+	skip = (fit[0].Npts < 15);
+	break;
+      default:
+	fprintf (stderr, "invalid chip order %d\n", coords[0].Npolyterms);
+	abort ();
+    }
+    if (skip) {
+      fprintf (stderr, "insufficient measurements (%d) for requested order (%d)\n", fit[0].Npts, coords[0].Npolyterms);
       fit_free (fit);
       free (values);
       return;
     }
+
+    fprintf (stderr, "scatter limit: %f based on %d detections; using %d of %d for fit\n", dRmax, Nscatter, fit[0].Npts, Nmatch);
+
+    // measure the fit, update the coords & object coordinates
     fit_eval (fit);
     fit_apply_coords (fit, coords);
@@ -82,4 +110,21 @@
    P,Q -> L,M (polynomial transformation : DIS)
    L,M -> X,Y (polynomial transformation : WRP)
+
+Some details about this function:
+
+- the initial value of the clipping radius is set based on the distribution of the dR
+  values for the bright sources.
+
+- NITER clipping passes are performed
+
+- the per-star astrometric errors are included in the fit.  The photcode table controls
+  whether only the reported positional errors are used, or if the magnitudes errors are
+  scaled to determine the error, and if there is a systematic floor for all sources.
+
+- input detections are pre-filtered on the basis of the photFlags, etc, in bcatalog
+
+- outlying detections are clipped in the iterative passes, but the clipped detections are
+  not recorded
+
 */
 
@@ -96,16 +141,2 @@
   fclose (f);
 */
-
-/*** XXX this function can be improved in 4 important ways:
-
-1) the initial value of the clipping radius could be set based on the distribution of the
-   value (eg, inner 50%) (better to use dL, dM rather than dP, dQ)
-
-2) multiple clipping passes could be performed
-
-3) the per-star astrometric errors could be included in the fit (non-systematic terms at
-   least, or scaled values of the mag error, or simply the mag error itself)
-
-4) we could ignore input detections on the basis of the photFlags
-
-*/