Changeset 41637 for branches/eam_branches/relphot.20210521/src/StarOps.c

Timestamp:

Jun 4, 2021, 10:46:36 AM (5 years ago)

Author:

eugene

Message:

add user-control over night, mosaic, image chisq & scatter cuts

File:

• branches/eam_branches/relphot.20210521/src/StarOps.c (modified) (4 diffs)

branches/eam_branches/relphot.20210521/src/StarOps.c

@@ -576 +576 @@
 }
 
+static float MinChiSqLim = NAN;
+static float MaxChiSqLim = NAN;
+static float MinScatterLim = NAN;
+static float MaxScatterLim = NAN;
+
 void clean_stars (Catalog *catalog, int Ncatalog) {
 
-  int i, j, Ndel, Nave, Ntot, mark, Ns, Nscat, Nchi, Nnan;
-  float dM;
-  double MaxScatter, MaxChisq;
-  double *xlist, *slist, *dlist;
+  int Ndel, Nave, Ntot, Nscat, Nchi, Nnan;
+  double *xlist, *slist;
+
+  if (isnan (MaxChiSqLim))   MaxChiSqLim   = STAR_CHISQ;
+  if (isnan (MaxScatterLim)) MaxScatterLim = STAR_SCATTER;
 
   StatType stats;
@@ -589 +595 @@
 
   /* find Mchisq median -> ChiSq lim must be > median */
-  for (i = Ntot = 0; i < Ncatalog; i++) {
+  for (int i = Ntot = 0; i < Ncatalog; i++) {
     Ntot += catalog[i].Naverage; 
   }
   ALLOCATE (xlist, double, Ntot);
   ALLOCATE (slist, double, Ntot);
-  ALLOCATE (dlist, double, Ntot);
 
   int Nsecfilt = GetPhotcodeNsecfilt ();
 
   // eliminate bad stars using the stats for a single secfilt at a time
-  for (Ns = 0; Ns < Nphotcodes; Ns ++) {
+  for (int Ns = 0; Ns < Nphotcodes; Ns ++) {
     
     int thisCode = photcodes[Ns][0].code;
     int Nsec = GetPhotcodeNsec(thisCode);
 
-    for (i = Ntot = 0; i < Ncatalog; i++) {
-      for (j = 0; j < catalog[i].Naverage; j++) {
+    for (int i = Ntot = 0; i < Ncatalog; i++) {
+      for (int j = 0; j < catalog[i].Naverage; j++) {
         if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;
         float Mchisq = catalog[i].secfilt[Nsecfilt*j+Nsec].Mchisq;
@@ -611 +616 @@
         xlist[Ntot] = Mchisq;
         slist[Ntot] = catalog[i].secfilt[Nsecfilt*j+Nsec].dMpsfChp;
-        dlist[Ntot] = 1;
         Ntot ++;
       }
     }
 
-    // XXX the limits for MaxChisq and MaxScatter should be
-    // user-defined.  This uses 
-    liststats (xlist, dlist, NULL, Ntot, &stats);
-    float ChiSqUpper90 = stats.Upper90;
-    MaxChisq = MAX (STAR_CHISQ, ChiSqUpper90);
-
-    liststats (slist, dlist, NULL, Ntot, &stats);
+    liststats (xlist, NULL, NULL, Ntot, &stats);
+    float ChiSqUpper90 = stats.Upper90; 
+    if (isnan (MinChiSqLim)) MinChiSqLim = 2.0*stats.median;             // chi-square cut cannot fall below this value (even if this is > MaxChiSqLim)
+    float ChiSqLimit = MAX(MinChiSqLim, MIN(MaxChiSqLim, ChiSqUpper90)); // chi-square cut should be between MinChiSqLim and MaxChiSqLim
+
+    liststats (slist, NULL, NULL, Ntot, &stats);
     float ScatterUpper90 = stats.Upper90;
-    MaxScatter = MAX (STAR_SCATTER, ScatterUpper90);
-    fprintf (stderr, "Max Scatter: %f, Max Chisq: %f\n", MaxScatter, MaxChisq);
+    if (isnan (MinScatterLim)) MinScatterLim = 2.0*stats.median;         // scatter cut cannto fall below this value (even if this is > MaxScatterLim)
+    float ScatterLimit = MAX(MinScatterLim, MIN(MaxScatterLim, ScatterUpper90));
+
+    fprintf (stderr, "STARS: ChiSqLimit: %f, ScatterLimit: %f | ChiSquare Upper 90: %f, Scatter Upper 90: %f\n", ChiSqLimit, ScatterLimit, ChiSqUpper90, ScatterUpper90);
 
     Ndel = Nave = Nscat = Nnan = Nchi = 0;
-    for (i = 0; i < Ncatalog; i++) {
-      for (j = 0; j < catalog[i].Naverage; j++) {
-        dM = catalog[i].secfilt[Nsecfilt*j+Nsec].dMpsfChp;
+    for (int i = 0; i < Ncatalog; i++) {
+      for (int j = 0; j < catalog[i].Naverage; j++) {
+        float dM = catalog[i].secfilt[Nsecfilt*j+Nsec].dMpsfChp;
         float Mchisq = catalog[i].secfilt[Nsecfilt*j+Nsec].Mchisq;
-        mark = (dM > MaxScatter) || (isnan(Mchisq)) || (Mchisq > MaxChisq);
+        int mark = (dM > ScatterLimit) || (isnan(Mchisq)) || (Mchisq > ChiSqLimit);
         if (mark) {
           catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_STAR_POOR;
           Ndel ++;
-          if (dM > MaxScatter)   { Nscat ++; }
-          if (isnan(Mchisq))     { Nnan ++; }
-          if (Mchisq > MaxChisq) { Nchi ++; }
+          if (dM > ScatterLimit)   { Nscat ++; }
+          if (isnan(Mchisq))       { Nnan ++; }
+          if (Mchisq > ChiSqLimit) { Nchi ++; }
         } else {
           catalog[i].secfilt[Nsecfilt*j+Nsec].flags &= ~ID_SECF_STAR_POOR;
@@ -645 +650 @@
       }
     }
-    fprintf (stderr, "%d stars marked variable (%d scat, %d nan, %d chi), %d total\n", Ndel, Nscat, Nnan, Nchi, Nave);
+    fprintf (stderr, "%d of %d stars marked variable (%d scat, %d nan, %d chi)\n", Ndel, Nave, Nscat, Nnan, Nchi);
   }
   free (xlist);
   free (slist);
-  free (dlist);
 }