Changeset 39122

Timestamp:

Nov 12, 2015, 5:41:54 PM (11 years ago)

Author:

eugene

Message:

do a better job fitting for the major & minor axis error value; handle edge cases; interpolate to get the galaxy mags

File:

• trunk/Ohana/src/addstar/src/loadgalphot_readstars.c (modified) (4 diffs)

trunk/Ohana/src/addstar/src/loadgalphot_readstars.c

@@ -237 +237 @@
     // I need to find the minimum position (interpolated) in this 2D space
 
+    // I want to use this flag to set a bit in secfilt.  is that possible?  or maybe galphot.dummy becomes flags
     if (!FitChisqMinimum (fit, &stars[i].galphot, &GAL_CHISQ[Nsample*i], &GAL_FLUX[Nsample*i], &GAL_FLUX_ERR[Nsample*i], Nsample, FR_MAJOR_MIN[i], FR_MAJOR_MAX[i], FR_MAJOR_DEL[i], FR_MINOR_MIN[i], FR_MINOR_MAX[i], FR_MINOR_DEL[i])) {
-      Nbad ++;
+      stars[i].flag = TRUE;
     }
 
@@ -298 +299 @@
   // here are the steps:
   // find the minimum chisq grid point (set X,Y values)
-  // choose the nearest 9 points
-  // fit 2D parabola (6 free parameters)
-  // use this to outlier-clip remaining points
-  // fit 2D parabola to the rest
+  // rescale chisq values so minimum is == Npix
+  // fit 2D parabola
 
   memset (mask, 0, Npts*sizeof(char));
@@ -324 +323 @@
   }
   if (!Nvalid) {
-    fprintf (stderr, "no good points\n");
     return FALSE; // we cannot do anything if there is no finite chisq
   }
   
-  if (Nvalid < 6) {
-    galphot->majorAxis = Xpt[iMin];
-    galphot->minorAxis = Ypt[iMin];
-    galphot->majorAxisErr = MajorDel;
-    galphot->minorAxisErr = MinorDel;
-    fprintf (stderr, "too few good points: %d\n", Nvalid);
-    return FALSE;
-  }
-
-  // re-fit all except the most extreme set
-  if (!fit2d (fit, Xpt, Ypt, chisq, chisqFit, mask, Npts)) {
-    // raise a flag?
-    galphot->majorAxis = Xpt[iMin];
-    galphot->minorAxis = Ypt[iMin];
-    galphot->majorAxisErr = MajorDel;
-    galphot->minorAxisErr = MinorDel;
-    fprintf (stderr, "failed fit\n");
-    return FALSE;
-  }
+  for (i = 0; i < Npts; i++) {
+    chisq[i] *= (galphot->Npix / chisqMin);
+  }
+
+  if (Nvalid < 6) goto bad_fit;
+
+# if (0)
+  {
+    FILE *f = fopen ("galphot.dump.dat", "w");
+    for (i = 0; i < Npts; i++) {
+      fprintf (f, "%d %d %f %f %f %f\n", i, mask[i], Xpt[i], Ypt[i], chisq[i], flux[i]);
+    }
+    fclose (f);
+  }
+# endif
+
+  // fit with 3 iterations, 5 sigma clipping (set above after fit2d_init)
+  // on failure (for any reason), the stars[i].flag is set to TRUE (XXX bad choice)
+  if (!fit2d (fit, Xpt, Ypt, chisq, chisqFit, mask, Npts)) goto bad_fit;
 
   // get X,Y for the chisq min from the fit:
@@ -353 +351 @@
   float Ymin = (fit->c10*fit->c11 - 2.0*fit->c20*fit->c01) * Det;
 
+  // allow the fitting minimum to be just outside the grid, but no further
+  int inRange = TRUE;
+  inRange = inRange && (Xmin > MajorMin - MajorDel);
+  inRange = inRange && (Xmin < MajorMax + MajorDel);
+  inRange = inRange && (Ymin > MinorMin - MinorDel);
+  inRange = inRange && (Ymin < MinorMax + MinorDel);
+
+  if (!inRange) goto bad_fit;
+
   // chisqMin @ (Xmin,Ymin)
-  chisqMin = 
-    fit->c00 + 
-    fit->c10*Xmin + fit->c20*Xmin*Xmin +
-    fit->c01*Ymin + fit->c02*Ymin*Ymin +
-    fit->c11*Xmin*Ymin;
-
-  float chisqOff = chisqMin + 1.0;
-
-  float A, B, C, Xoff, Yoff;
-
-  // Xoff @ chisqMin + 1.0
+  float A, B, C, Q;
+
+  // Xoff @ chisqMin + 1.0, delta Y = 0.0
   A = fit->c20;
-  B = fit->c10 + fit->c11*Ymin;
-  C = fit->c00 + fit->c01*Ymin + fit->c02*Ymin*Ymin - chisqOff;
-  Xoff = (-B + sqrt(B*B - 4*A*C)) / (2.0 * A);
-
-  // Xoff @ chisqMin + 1.0
+  B = fit->c10 + 2*fit->c20*Xmin + fit->c11*Ymin;
+  C = -1;
+  Q = B*B - 4*A*C;
+  if (Q < 0.0) goto bad_err;
+  
+  float dXmin = (-B + sqrt(Q)) / (2.0 * A);
+
+  // Yoff @ chisqMin + 1.0, delta X = 0.0
   A = fit->c02;
-  B = fit->c01 + fit->c11*Xmin;
-  C = fit->c00 + fit->c10*Xmin + fit->c20*Xmin*Xmin - chisqOff;
-  Yoff = (-B + sqrt(B*B - 4*A*C)) / (2.0 * A);
+  B = fit->c01 + 2*fit->c02*Ymin + fit->c11*Xmin;
+  C = -1;
+  Q = B*B - 4*A*C;
+  if (Q < 0.0) goto bad_err;
+
+  float dYmin = (-B + sqrt(Q)) / (2.0 * A);
 
   galphot->majorAxis = Xmin;
   galphot->minorAxis = Ymin;
-  galphot->majorAxisErr = fabs(Xoff - Xmin);
-  galphot->minorAxisErr = fabs(Yoff - Ymin);
+  galphot->majorAxisErr = dXmin;
+  galphot->minorAxisErr = dYmin;
+
+  // convert (Xmin, Ymin) to (iXmin, iYmin)
+  float iXminF = (Xmin - MajorMin) / MajorDel;
+  float iYminF = (Ymin - MinorMin) / MinorDel;
+
+  // interpolate in X and in Y
+  int   iXmin = floor(iXminF);
+  iXmin = MAX(0,MIN(Nx - 2, iXmin)); // force range of iXmin to be 0,Nx-1
+  float fXmin = iXminF - iXmin;
+
+  int   iYmin = floor(iYminF);
+  iYmin = MAX(0,MIN(Ny - 2, iYmin)); // force range of iYmin to be 0,Ny-1
+  float fYmin = iYminF - iYmin;
+
+  float V00 = flux[iXmin+0 + (iYmin+0)*Nx];
+  float V01 = flux[iXmin+0 + (iYmin+1)*Nx];
+  float V10 = flux[iXmin+1 + (iYmin+0)*Nx];
+  float V11 = flux[iXmin+1 + (iYmin+1)*Nx];
+
+  float Vx0 = V10*fXmin + V00*(1.0 - fXmin);
+  float Vx1 = V11*fXmin + V01*(1.0 - fXmin);
+
+  float value = fYmin * Vx1 + (1.0 - fYmin) * Vx0;
+
+  float dV00 = fluxErr[iXmin+0 + (iYmin+0)*Nx];
+  float dV01 = fluxErr[iXmin+0 + (iYmin+1)*Nx];
+  float dV10 = fluxErr[iXmin+1 + (iYmin+0)*Nx];
+  float dV11 = fluxErr[iXmin+1 + (iYmin+1)*Nx];
+
+  float dVx0 = dV10*fXmin + dV00*(1.0 - fXmin);
+  float dVx1 = dV11*fXmin + dV01*(1.0 - fXmin);
+
+  float dValue = fYmin * dVx1 + (1.0 - fYmin) * dVx0;
 
   // use bilinear interpolation to choose Flux @ (Xmin,Ymin)?
+  // float magRaw    = -2.5*log10(flux[iMin]) + ZeroPt; // correct for exptime?
+  // float magErrRaw = sqrt(fluxErr[iMin]) / flux[iMin];
+
+  galphot->mag    = -2.5*log10(value) + ZeroPt; // correct for exptime?
+  galphot->magErr = sqrt(dValue) / value;
+  galphot->chisq  = chisqMin;
+  
+  return TRUE;
+
+bad_fit:
+  galphot->majorAxis = Xpt[iMin];
+  galphot->minorAxis = Ypt[iMin];
+  galphot->majorAxisErr = MajorDel;
+  galphot->minorAxisErr = MinorDel;
+
   galphot->mag    = -2.5*log10(flux[iMin]) + ZeroPt; // correct for exptime?
   galphot->magErr = sqrt(fluxErr[iMin]) / flux[iMin];
   galphot->chisq  = chisqMin;
-  
-  return TRUE;
+  return FALSE;
+
+bad_err:
+  galphot->majorAxisErr = MajorDel;
+  galphot->minorAxisErr = MinorDel;
+  return FALSE;
 }