Index: trunk/Ohana/src/opihi/lib.data/starfuncs.c
===================================================================
--- trunk/Ohana/src/opihi/lib.data/starfuncs.c	(revision 2826)
+++ trunk/Ohana/src/opihi/lib.data/starfuncs.c	(revision 2827)
@@ -92,2 +92,88 @@
 }
 
+/* use a circular aperture */
+int get_rough_star (float *data, int Nx, int Ny, int x, int y,
+		       float *xc, float *yc, 
+		       float *sx, float *sy, 
+		       float *zs, float *zp) {
+
+  double *sky;
+  double Raper, Rinner, Router, Ra2, Ri2, Ro2, rad2;
+  int i, j, Npts, Nsky;
+  int Xs, Xe, Ys, Ye, off, Xc, Yc;
+  double peak, fsky, value;
+  double Sx, Sy, Sx2, Sy2, Sum;
+  
+  /* define circular boundaries */
+  Raper  = 15;  Ra2 = SQ(Raper);
+  Rinner = 25;  Ri2 = SQ(Rinner);
+  Router = 35;  Ro2 = SQ(Router);
+
+  /* measure the sky level */
+  ALLOCATE (sky, double, SQ(2*Router + 1));
+
+  /* boundaries for the outer sky region */
+  Xs = MAX (x - Router, 0);
+  Xe = MIN (x + Router + 1, Nx);
+  Ys = MAX (y - Router, 0);
+  Ye = MIN (y + Router + 1, Ny);
+
+  Nsky = 0;  
+  for (j = Ys; j < Ye; j++) {
+    off = j*Nx;
+    for (i = Xs; i < Xe; i++) { 
+      rad2 = SQ(i - x) + SQ(j - y);
+      if (rad2 > Ro2) continue;
+      if (rad2 < Ri2) continue;
+      sky[Nsky] = data[i+off];
+      Nsky ++;
+    }
+  }
+  sort (sky, Nsky);
+  for (Npts = fsky = 0, i = 0.25*Nsky; i < 0.75*Nsky; i++, Npts += 1.0) {
+    fsky += sky[i];
+  }
+  fsky = fsky / Npts;
+  free (sky);
+
+  /* boundaries for the star region */
+  Xs = MAX (x - Raper, 0);
+  Xe = MIN (x + Raper + 1, Nx);
+  Ys = MAX (y - Raper, 0);
+  Ye = MIN (y + Raper + 1, Ny);
+
+  /** note that this will fail on negative flux objects */
+  peak = Npts = 0;
+  Sx = Sy = Sx2 = Sy2 = Sum = 0;
+  for (j = Ys; j < Ye; j++) {
+    off = j*Nx;
+    Yc = j - y;
+    for (i = Xs; i < Xe; i++) {
+      Xc = i - x;
+      rad2 = SQ(Xc) + SQ(Yc);
+      if (rad2 > Ro2) continue;
+      value = data[i+off] - fsky;
+      Sx  += Xc*value;
+      Sy  += Yc*value;
+      Sx2 += Xc*Xc*value;
+      Sy2 += Yc*Yc*value;
+      Sum += value;
+      Npts ++;
+      if (value > peak) peak = value;
+    }
+  }
+
+  *xc = Sx / Sum;
+  *yc = Sy / Sum;
+  *sx = sqrt (fabs (Sx2 / Sum - SQ(*xc)));
+  *sy = sqrt (fabs (Sy2 / Sum - SQ(*yc)));
+  *xc += x;
+  *yc += y;
+  *zs = Sum;
+  *zp = peak;
+  /* note sigma is rough: round-off errors can introduce errors */
+  /* using values relative to x,y should minimize this effect */
+
+  return (Npts);
+}
+