Changeset 40202 for branches/eam_branches/ohana.20170822

Timestamp:

Nov 4, 2017, 1:12:06 PM (9 years ago)

Author:

eugene

Message:

add option to measure star stats in a rectangular box; fix long-standing bug with limits for all negative values

Location:

branches/eam_branches/ohana.20170822/src/opihi

Files:

• cmd.astro/star.c (modified) (4 diffs)
• include/data.h (modified) (1 diff)
• lib.data/graphtools.c (modified) (2 diffs)
• lib.data/starfuncs.c (modified) (1 diff)

branches/eam_branches/ohana.20170822/src/opihi/cmd.astro/star.c

@@ -3 +3 @@
 int star (int argc, char **argv) {
 
-  int x, y, N, dx, Nborder;
+  int x, y, N, Nborder;
   double max;
   Buffer *buf;
@@ -33 +33 @@
   }
   
+  int dx = 11;
+  int dy = 11;
+  int BOX = FALSE;
+  if ((N = get_argument (argc, argv, "-box"))) {
+    remove_argument (N, &argc, argv);
+    dx  = atoi(argv[N]);
+    remove_argument (N, &argc, argv);
+    dy  = atoi(argv[N]);
+    remove_argument (N, &argc, argv);
+    BOX = TRUE;
+  }
+
   if ((argc != 4) && (argc != 5)) {
-    gprint (GP_ERR, "USAGE: star (buffer) x y [dx] [-border N] [-sat cnts]\n");
+    gprint (GP_ERR, "USAGE: star (buffer) x y [dx] [-border N] [-sat cnts] [-box dx dy]\n");
     gprint (GP_ERR, " dx is the aperture diameter, but is adjusted up to the next odd number\n");
     return (FALSE);
@@ -40 +52 @@
   if ((buf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
 
-  dx = 11;
   x = atof (argv[2]);
   y = atof (argv[3]);
@@ -47 +58 @@
   }
 
-  get_aperture_stats (&buf[0].matrix, x, y, dx, Nborder, max, VERBOSE);
+  if (BOX) {
+    get_box_stats (&buf[0].matrix, x, y, dx, dy, Nborder, max, VERBOSE);
+  } else {
+    get_aperture_stats (&buf[0].matrix, x, y, dx, Nborder, max, VERBOSE);
+  }
   
   return (TRUE);

branches/eam_branches/ohana.20170822/src/opihi/include/data.h

@@ -162 +162 @@
 /* starfuncs.c */
 double get_aperture_stats (Matrix *matrix, int X, int Y, int Npix, int Nborder, double max, int VERBOSE);
+double get_box_stats (Matrix *matrix, int X, int Y, int dX, int dY, int Nborder, double max, int VERBOSE);
+
 int set_rough_radii (double Ra, double Ri, double Ro);
 int get_rough_star (float *data, int Nx, int Ny, int x, int y, opihi_flt *xc, opihi_flt *yc, opihi_flt *sx, opihi_flt *sy, opihi_flt *sxy, opihi_flt *zs, opihi_flt *zp, opihi_flt *sk);

branches/eam_branches/ohana.20170822/src/opihi/lib.data/graphtools.c

@@ -10 +10 @@
   if (xvec != NULL) {
     if (xvec->type == OPIHI_FLT) {
-      maxX = DBL_MIN;
+      maxX = -DBL_MAX;
       minX = DBL_MAX;
       for (i = 0; i < xvec[0].Nelements; i++) {
@@ -33 +33 @@
   if (yvec != NULL) {
     if (yvec->type == OPIHI_FLT) {
-      maxY = DBL_MIN;
+      maxY = -DBL_MAX;
       minY = DBL_MAX;
       for (i = 0; i < yvec[0].Nelements; i++) {

branches/eam_branches/ohana.20170822/src/opihi/lib.data/starfuncs.c

@@ -101 +101 @@
 }
 
+double get_box_stats (Matrix *matrix, int X, int Y, int dX, int dY, int Nborder, double max, int VERBOSE) {
+
+  double *ring;
+  double x, y, x2, y2, xy, I, sky, FWHMx, FWHMy, value, mag, Sxy;
+  int i, j, n, Nring, Nmax;
+  double Npts, gain, dsky2, dmag, peak, offset;
+  char *string;
+  
+  string = get_variable ("GAIN");
+  if (string == (char *) NULL) {
+    gprint (GP_ERR, "assuming a value of 1.0\n");
+    gain = 1.0;
+  } else {
+    gain = atof (string);
+  }
+  Nborder = MAX (1, Nborder);
+  Nborder = MIN (1000, Nborder);
+  
+  int dX2 = (int)(0.5*dX);
+  int dY2 = (int)(0.5*dY);
+  dX = 2 * dX2 + 1;
+  dY = 2 * dY2 + 1;
+
+  Nring = 2*Nborder*(dX + 2*Nborder) + 2*Nborder*(dY + 2*Nborder);
+  ALLOCATE (ring, double, Nring);
+  bzero (ring, sizeof(double)*Nring);
+
+  // get the pixels in the border regions:
+  // XXX gfits_get_matrix_value returns 0 for out-of-bounds pixels, but should return NAN
+  // and they should be skipped
+  n = 0;  
+  for (j = 0; j < Nborder; j++) {
+    for (i = X - dX2 - Nborder; i < X + dX2 + Nborder + 1; i++) {
+      value = gfits_get_matrix_value (matrix, i, (int)(Y - dY2 - j));
+      if (isfinite(value)) { ring[n] = value; n++; }
+      value = gfits_get_matrix_value (matrix, i, (int)(Y + dY2 + j));
+      if (isfinite(value)) { ring[n] = value; n++; }
+    }
+    for (i = Y - dY2; i < Y + dY2 + 1; i++) {
+      value = gfits_get_matrix_value (matrix, (int)(X - dX2 - j), i);
+      if (isfinite(value)) { ring[n] = value; n++; }
+      value = gfits_get_matrix_value (matrix, (int)(X + dX2 + j), i);
+      if (isfinite(value)) { ring[n] = value; n++; }
+    }
+  }
+  Nring = n;
+  dsort (ring, Nring);
+  for (Npts = sky = dsky2 = 0, 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);
+
+  float dx, dy;
+
+  peak = 0;
+  Npts = Nmax = 0;
+  x = y = x2 = y2 = xy = I = 0;
+  for (i = X - dX2; i < X + dX2 + 1; i++) {
+    for (j = Y - dY2; j < Y + dY2 + 1; j++) {
+      value = gfits_get_matrix_value (matrix, i, j);
+      if (!isfinite(value)) continue;
+      offset = value - sky;
+      dx = i - X;
+      dy = j - Y;
+      x  += dx*offset;
+      y  += dy*offset;
+      x2 += dx*dx*offset;
+      y2 += dy*dy*offset;
+      xy += dx*dy*offset;
+      I  += offset;
+      Npts ++;
+      if (value > max) {
+        Nmax ++;
+      }
+      if (value > peak) peak = value;
+    }
+  }
+
+  x = x / I;
+  y = y / I;
+  FWHMx = 2.355*sqrt (fabs(x2 / I - x*x));
+  FWHMy = 2.355*sqrt (fabs(y2 / I - y*y));
+  Sxy   = xy / I - x*y;
+  mag = -2.5*log10(I);
+
+  // flux_error = sqrt( I + Npts*dsky2 )
+  // dmag = 1.086 * flux_error / flux
+  dmag = 1.086 * sqrt (fabs(I + Npts*dsky2)) / (gain * I);
+  x = x + X;
+  y = y + Y;
+  
+  set_variable ("Xg", x);
+  set_variable ("Yg", y);
+  set_variable ("SXg", FWHMx);
+  set_variable ("SYg", FWHMy);
+  set_variable ("SXYg", Sxy);
+  set_variable ("Sg", sky);
+  set_variable ("dSg", sqrt (fabs (dsky2)));
+  set_variable ("Zg", mag);
+  set_variable ("dZg", dmag);
+  set_variable ("Zcg", I);
+  set_variable ("Zpk", peak);
+  set_int_variable ("Nsat", Nmax);
+  set_int_variable ("Npts", Npts);
+  
+  if (VERBOSE) gprint (GP_LOG, "%f %f %f %f %f %f %f %f\n", x, y, FWHMx, FWHMy, sky, I, mag, dmag);
+
+  return (mag);
+
+}
+
 static double Raper  =  5;
 static double Rinner = 10;