Index: trunk/Ohana/src/opihi/lib.data/starfuncs.c
===================================================================
--- trunk/Ohana/src/opihi/lib.data/starfuncs.c	(revision 2827)
+++ trunk/Ohana/src/opihi/lib.data/starfuncs.c	(revision 2843)
@@ -92,25 +92,40 @@
 }
 
+static double Raper  =  5;
+static double Rinner = 10;
+static double Router = 15;
+static double *sky = NULL;
+
+int set_rough_radii (double Ra, double Ri, double Ro) {
+
+  Raper = Ra;
+  Rinner = Ri;
+  Router = Ro;
+  if (sky == NULL) {
+    ALLOCATE (sky, double, SQ(2*Router + 1));
+  } else {
+    REALLOCATE (sky, double, SQ(2*Router + 1));
+  }
+  return (TRUE);
+}
+
 /* 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;
+		    float *xc, float *yc, 
+		    float *sx, float *sy, float *sxy,
+		    float *zs, float *zp, float *sk) {
+
+  double 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;
+  double Sx, Sy, Sx2, Sy2, Sxy, Sum;
   
   /* define circular boundaries */
-  Raper  = 15;  Ra2 = SQ(Raper);
-  Rinner = 25;  Ri2 = SQ(Rinner);
-  Router = 35;  Ro2 = SQ(Router);
+  Ra2 = SQ(Raper);
+  Ri2 = SQ(Rinner);
+  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);
@@ -119,4 +134,6 @@
   Ye = MIN (y + Router + 1, Ny);
 
+/* this sample uses a circular aperture */
+# if (0)
   Nsky = 0;  
   for (j = Ys; j < Ye; j++) {
@@ -135,5 +152,54 @@
   }
   fsky = fsky / Npts;
-  free (sky);
+# else
+
+/* this sample uses a square outer annulus, without loop if tests */
+  Nsky = 0;  
+  Xs = MAX (x - Router, 0);
+  Xe = MIN (x - Rinner + 1, Nx);
+  Ys = MAX (y - Rinner, 0);
+  Ye = MIN (y + Rinner + 1, Ny);
+  for (j = Ys; j < Ye; j++) {
+    off = j*Nx;
+    for (i = Xs; i < Xe; i++) { 
+      sky[Nsky] = data[i+off];
+      Nsky ++;
+    }
+  }
+  Xs = MAX (x + Rinner, 0);
+  Xe = MIN (x + Router + 1, Nx);
+  for (j = Ys; j < Ye; j++) {
+    off = j*Nx;
+    for (i = Xs; i < Xe; i++) { 
+      sky[Nsky] = data[i+off];
+      Nsky ++;
+    }
+  }
+  Xs = MAX (x - Rinner, 0);
+  Xe = MIN (x - Rinner + 1, Nx);
+  Ys = MAX (y - Router, 0);
+  Ye = MIN (y - Rinner + 1, Ny);
+  for (j = Ys; j < Ye; j++) {
+    off = j*Nx;
+    for (i = Xs; i < Xe; i++) { 
+      sky[Nsky] = data[i+off];
+      Nsky ++;
+    }
+  }
+  Ys = MAX (y + Rinner, 0);
+  Ye = MIN (y + Router + 1, Ny);
+  for (j = Ys; j < Ye; j++) {
+    off = j*Nx;
+    for (i = Xs; i < Xe; i++) { 
+      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;
+# endif
 
   /* boundaries for the star region */
@@ -145,5 +211,5 @@
   /** note that this will fail on negative flux objects */
   peak = Npts = 0;
-  Sx = Sy = Sx2 = Sy2 = Sum = 0;
+  Sx = Sy = Sx2 = Sy2 = Sxy = Sum = 0;
   for (j = Ys; j < Ye; j++) {
     off = j*Nx;
@@ -158,4 +224,5 @@
       Sx2 += Xc*Xc*value;
       Sy2 += Yc*Yc*value;
+      Sxy += Xc*Yc*value;
       Sum += value;
       Npts ++;
@@ -168,8 +235,10 @@
   *sx = sqrt (fabs (Sx2 / Sum - SQ(*xc)));
   *sy = sqrt (fabs (Sy2 / Sum - SQ(*yc)));
+  *sxy = Sxy / Sum;
   *xc += x;
   *yc += y;
   *zs = Sum;
   *zp = peak;
+  *sk = fsky;
   /* note sigma is rough: round-off errors can introduce errors */
   /* using values relative to x,y should minimize this effect */