Changeset 40291 for trunk/Ohana/src/opihi/lib.data/starfuncs.c

Timestamp:

Dec 13, 2017, 10:53:48 AM (9 years ago)

Author:

eugene

Message:

merge EAM development branch changes for DR2 into trunk (add PS1_V6 dvo format; change Mcal to McalPSF, McalAPER; change opihi int vectors to 64bit)

Location:

trunk/Ohana

Files:

• . (modified) (1 prop)
• src/opihi/lib.data/starfuncs.c (modified) (1 diff)

trunk/Ohana/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;