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);
+}
+
Index: /trunk/Ohana/src/opihi/mana/Makefile
===================================================================
--- /trunk/Ohana/src/opihi/mana/Makefile	(revision 2826)
+++ /trunk/Ohana/src/opihi/mana/Makefile	(revision 2827)
@@ -31,4 +31,5 @@
 
 manacmds = \
+$(SDIR)/rawstars.$(ARCH).o \
 $(SDIR)/findpeaks.$(ARCH).o 
 
Index: /trunk/Ohana/src/opihi/mana/init.c
===================================================================
--- /trunk/Ohana/src/opihi/mana/init.c	(revision 2826)
+++ /trunk/Ohana/src/opihi/mana/init.c	(revision 2827)
@@ -2,7 +2,9 @@
 
 int findpeaks	    PROTO((int, char **));
+int rawstars	    PROTO((int, char **));
 
 static Command cmds[] = {  
   {"findpeaks",   findpeaks,    "find image peaks"},
+  {"rawstars",    rawstars,     "find raw star stats"},
 }; 
 
Index: /trunk/Ohana/src/opihi/mana/rawstars.c
===================================================================
--- /trunk/Ohana/src/opihi/mana/rawstars.c	(revision 2827)
+++ /trunk/Ohana/src/opihi/mana/rawstars.c	(revision 2827)
@@ -0,0 +1,68 @@
+# include "mana.h"
+
+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);
+
+int rawstars (int argc, char **argv) {
+  
+  int i, x, y, Nx, Ny, Np;
+  float *v;
+  Vector *xp, *yp;
+  Vector *xc, *yc, *sx, *sy, *zs, *zp;
+  Buffer *buff;
+
+  if (argc < 4) goto usage;
+
+  if ((buff = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
+  if ((xp = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((yp = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if (xp[0].Nelements != yp[0].Nelements) {
+    fprintf (stderr, "vectors are not the same length\n");
+    return (FALSE);
+  }
+
+  /* output vectors */
+  if ((xc = SelectVector ("xc", ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((yc = SelectVector ("yc", ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((sx = SelectVector ("sx", ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((sy = SelectVector ("sy", ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((zs = SelectVector ("zs", ANYVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((zp = SelectVector ("zp", ANYVECTOR, TRUE)) == NULL) return (FALSE);
+
+  Nx = buff[0].matrix.Naxis[0];
+  Ny = buff[0].matrix.Naxis[1];
+  Np = xp[0].Nelements;
+
+  REALLOCATE (xc[0].elements, float, Np);
+  REALLOCATE (yc[0].elements, float, Np);
+  REALLOCATE (sx[0].elements, float, Np);
+  REALLOCATE (sy[0].elements, float, Np);
+  REALLOCATE (zs[0].elements, float, Np);
+  REALLOCATE (zp[0].elements, float, Np);
+  xc[0].Nelements = yc[0].Nelements = sx[0].Nelements = Np;
+  sy[0].Nelements = zs[0].Nelements = zp[0].Nelements = Np;
+
+  v = (float *) buff[0].matrix.buffer;
+  for (i = 0; i < Np; i++) {
+    x = xp[0].elements[i];
+    y = yp[0].elements[i];
+    if (x < 0) continue;
+    if (x >= Nx) continue;
+    if (y < 0) continue;
+    if (y >= Ny) continue;
+
+    get_rough_star (v, Nx, Ny, x, y, 
+		    &xc[0].elements[i], 
+		    &yc[0].elements[i], 
+		    &sx[0].elements[i], 
+		    &sy[0].elements[i], 
+		    &zs[0].elements[i], 
+		    &zp[0].elements[i]);
+  }
+
+  return (TRUE);
+
+ usage:
+  fprintf (stderr, "rawstars (buffer) (xp) (yp)\n");
+  return (FALSE);
+}
+