Index: /branches/eam_branches/ipp-20120405/Ohana/src/opihi/cmd.astro/Makefile
===================================================================
--- /branches/eam_branches/ipp-20120405/Ohana/src/opihi/cmd.astro/Makefile	(revision 33738)
+++ /branches/eam_branches/ipp-20120405/Ohana/src/opihi/cmd.astro/Makefile	(revision 33739)
@@ -46,4 +46,5 @@
 $(SRC)/ringflux.$(ARCH).o	   \
 $(SRC)/petrosian.$(ARCH).o	   \
+$(SRC)/kronflux.$(ARCH).o	   \
 $(SRC)/multifit.$(ARCH).o	   \
 $(SRC)/objload.$(ARCH).o	   \
Index: /branches/eam_branches/ipp-20120405/Ohana/src/opihi/cmd.astro/init.c
===================================================================
--- /branches/eam_branches/ipp-20120405/Ohana/src/opihi/cmd.astro/init.c	(revision 33738)
+++ /branches/eam_branches/ipp-20120405/Ohana/src/opihi/cmd.astro/init.c	(revision 33739)
@@ -32,4 +32,5 @@
 int ringflux                PROTO((int, char **));
 int petrosian               PROTO((int, char **));
+int kronflux                PROTO((int, char **));
 int mkgauss                 PROTO((int, char **));
 int mksersic                PROTO((int, char **));
@@ -88,4 +89,5 @@
   {1, "ringflux",    ringflux,     "mean flux in a ring"},
   {1, "petrosian",   petrosian,    "petrosian parameters given radial bins"},
+  {1, "kronflux",    kronflux,     "measure kronflux stats"},
   {1, "multifit",    multifit,     "fit multi-order spectrum"},
   {1, "objload",     objload,      "plot obj data on Ximage "},
Index: /branches/eam_branches/ipp-20120405/Ohana/src/opihi/cmd.astro/kronflux.c
===================================================================
--- /branches/eam_branches/ipp-20120405/Ohana/src/opihi/cmd.astro/kronflux.c	(revision 33739)
+++ /branches/eam_branches/ipp-20120405/Ohana/src/opihi/cmd.astro/kronflux.c	(revision 33739)
@@ -0,0 +1,173 @@
+# include "astro.h"
+int kronflux_stats (Matrix *matrix, int X, int Y, int radius, int maxradius, int border, float alpha);
+
+int VERBOSE = TRUE;
+
+// manual calculation of the Kron parameters  
+int kronflux (int argc, char **argv) {
+
+  int x, y, N, radius, border, maxradius;
+  float alpha;
+  Buffer *buf;
+
+  radius = 5;
+  if ((N = get_argument (argc, argv, "-radius"))) {
+    remove_argument (N, &argc, argv);
+    radius  = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  maxradius = 100;
+  if ((N = get_argument (argc, argv, "-maxradius"))) {
+    remove_argument (N, &argc, argv);
+    maxradius  = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  border = 5;
+  if ((N = get_argument (argc, argv, "-border"))) {
+    remove_argument (N, &argc, argv);
+    border  = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  
+  alpha = 2.5;
+  if ((N = get_argument (argc, argv, "-alpha"))) {
+    remove_argument (N, &argc, argv);
+    alpha  = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  
+  VERBOSE = TRUE;
+  if ((N = get_argument (argc, argv, "-q"))) {
+    remove_argument (N, &argc, argv);
+    VERBOSE = FALSE;
+  }
+  
+  if (argc != 4) {
+    gprint (GP_ERR, "USAGE: kronflux (buffer) x y [-radius N] [-border N] [-alpha F] [-q]\n");
+    gprint (GP_ERR, " -radius is the radius of the initial window (default 5)\n");
+    gprint (GP_ERR, " -border is the size of the sky annulus (default 5)\n");
+    gprint (GP_ERR, " -alpha  is the scale factor for kron radius vs 1st radial moment (default 2.5)\n");
+    return (FALSE);
+  }
+  if ((buf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
+  x = atof (argv[2]);
+  y = atof (argv[3]);
+
+  kronflux_stats (&buf[0].matrix, x, y, radius, maxradius, border, alpha);
+
+  return TRUE;
+}
+
+// I'd like to be able to skip the sky measurement
+// if Nborder > 0, calculate local sky as the median of the pixels in the region between radius and border
+int kronflux_stats (Matrix *matrix, int X, int Y, int radius, int maxradius, int border, float alpha) {
+
+  int i, j;
+
+  double gain;
+  char *string = get_variable ("GAIN");
+  if (string == (char *) NULL) {
+    if (VERBOSE) gprint (GP_ERR, "assuming a value of 1.0\n");
+    gain = 1.0;
+  } else {
+    gain = atof (string);
+  }
+
+  float *data = (float *) matrix->buffer;
+  int Nx = matrix->Naxis[0];
+  int Ny = matrix->Naxis[1];
+  
+  double sky = 0.0;
+  double dsky2 = 0.0;
+  if (border > 0) {
+    int radius2 = radius + border;
+
+    // allocate enough space for all pixels in the 2*radius2 + 1 square
+    int NRING = SQ(2*radius2 + 1);
+    double *ring = NULL;
+    ALLOCATE (ring, double, NRING);
+    bzero (ring, sizeof(double)*NRING);
+
+    int Nring = 0;  
+    for (j = -radius2; j <= radius2; j++) {
+      int iy = j + Y;
+      if (iy < 0) continue;
+      if (iy >= Ny) continue;
+      for (i = -radius2; i <= radius2; i++) {
+	if (hypot(i, j) <= radius) continue;
+	int ix = i + X;
+	if (ix < 0) continue;
+	if (ix >= Nx) continue;
+	ring[Nring] = data[ix + iy*Nx];
+	Nring ++;
+      }
+    }
+    dsort (ring, Nring);
+
+    int Npts = 0;
+    for (i = 0.25*Nring; i < 0.75*Nring; i++, Npts += 1.0) {
+      sky += ring[i];
+      dsky2 += ring[i]*ring[i];
+    }
+    sky = sky / Npts;
+    dsky2 = dsky2 / Npts - sky*sky;
+    free (ring);
+  }
+
+  // find the 1st radial moment
+  float flux1 = 0.0;
+  float Srf = 0.0;
+  for (j = -radius; j <= radius; j++) {
+    int iy = j + Y;
+    if (iy < 0) continue;
+    if (iy >= Ny) continue;
+    for (i = -radius; i <= radius; i++) {
+      float r = hypot(i, j);
+      if (r > radius) continue;
+      int ix = i + X;
+      if (ix < 0) continue;
+      if (ix >= Nx) continue;
+      float value = data[ix + iy*Nx] - sky;
+	
+      flux1 += value;
+      Srf += value * r;
+    }
+  }
+
+  float r1 = Srf / flux1;
+  float rK = r1 * alpha;
+
+  rK = MIN(MAX(rK, 5.0), maxradius);
+
+  // measure the flux within A*r1
+  float fluxK = 0.0;
+  for (j = -rK; j <= rK; j++) {
+    int iy = j + Y;
+    if (iy < 0) continue;
+    if (iy >= Ny) continue;
+    for (i = -rK; i <= rK; i++) {
+      float r = hypot(i, j);
+      if (r > rK) continue;
+      int ix = i + X;
+      if (ix < 0) continue;
+      if (ix >= Nx) continue;
+      float value = data[ix + iy*Nx] - sky;
+	
+      fluxK += value;
+    }
+  }
+
+  float magK = -2.5*log10(fluxK);
+  set_variable ("KronFlux", fluxK);
+  set_variable ("KronRadius", rK);
+  set_variable ("KronMag", magK);
+  set_variable ("KronF1", flux1);
+  set_variable ("KronR1", r1);
+  
+  if (VERBOSE) gprint (GP_LOG, "%f %f %f %f %f\n", magK, fluxK, rK, flux1, r1);
+
+  return TRUE;
+}
+