Index: /branches/eam_branch_20080223/Ohana/src/relastro/src/FitChip.c
===================================================================
--- /branches/eam_branch_20080223/Ohana/src/relastro/src/FitChip.c	(revision 16791)
+++ /branches/eam_branch_20080223/Ohana/src/relastro/src/FitChip.c	(revision 16792)
@@ -1,22 +1,67 @@
 # include "relastro.h"
+# define SCATTER_MAX_ERROR 0.05
 
 void FitChip (StarData *raw, StarData *ref, int Nmatch, Coords *coords) {
 
-  int i;
+  int i, Nscatter, Niter;
   CoordFit *fit;
+  double dL, dM, dR, dRmax, *values;
 
-  fit = fit_init (coords[0].Npolyterms);
-  for (i = 0; i < Nmatch; i++) {
-    if (raw[i].mask) continue;
-    fit_add (fit, raw[i].X, raw[i].Y, ref[i].L, ref[i].M, 1.0);
+  ALLOCATE (values, double, Nmatch);
+
+  for (Niter = 0; Niter < 3; 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;
+
+      dL = raw[i].L - ref[i].L;
+      dM = raw[i].M - ref[i].M;
+      dR = hypot (dL, dM);
+    
+      values[Nscatter] = dR;
+      Nscatter++;
+    }
+
+    // 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)];
+
+    fit = fit_init (coords[0].Npolyterms);
+
+    for (i = 0; i < Nmatch; i++) {
+      if (raw[i].mask) continue;
+
+      // require radius of XXX arcsec
+      dL = raw[i].L - ref[i].L;
+      dM = raw[i].M - ref[i].M;
+      dR = hypot (dL, dM);
+
+      if (dR > dRmax) {
+	continue;
+	raw[i].mask = TRUE;
+      }
+    
+      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);
+
+    if (fit[0].Npts < 25) {
+      fit_free (fit);
+      free (values);
+      return;
+    }
+    fit_eval (fit);
+    fit_apply_coords (fit, coords);
+    fit_free (fit);
+
+    for (i = 0; i < Nmatch; i++) {
+      XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, coords);
+    }
   }
-  fit_eval (fit);
-  fit_apply_coords (fit, coords);
-  fit_free (fit);
 
-  // apply new coords to raw (X,Y -> L,M)
-  for (i = 0; i < Nmatch; i++) {
-    XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, coords);
-  }
+  free (values);
+  return;
 }
 
@@ -32,3 +77,28 @@
 */
 
-/* XXX I'm not using the errors at all : this could at least be done with the dMag values */
+/* example using fit_apply() :
+
+  f = fopen ("test3.dat", "w");
+
+  // apply new coords to raw (X,Y -> L,M)
+  for (i = 0; i < Nmatch; i++) {
+    fprintf (f, "%f %f  %f %f  ", raw[i].X, raw[i].Y, raw[i].L, raw[i].M);
+    fit_apply (newfit, &L1, &M1, raw[i].X - coords[0].crpix1, raw[i].Y - coords[0].crpix2);
+    fprintf (f, "%f %f\n", L1, M1);
+  }
+  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
+
+*/