Changeset 41662 for trunk/Ohana/src/relphot/src/ImageOps.c

Timestamp:

Jun 14, 2021, 11:02:20 AM (5 years ago)

Author:

eugene

Message:

add code to rationalize the zero points on each pass

File:

• trunk/Ohana/src/relphot/src/ImageOps.c (modified) (1 diff)

trunk/Ohana/src/relphot/src/ImageOps.c

@@ -951 +951 @@
 }
 
+// find the median zero point and subtract it (for each active average photcode)
+void rationalize_zeropoints (int Niter) {
+
+  if (VERBOSE) fprintf (stderr, "rationalize zero points\n");
+
+  // if we are calculating zero points for tgroups,
+  // rationalize based on the tgroups
+  if (TGROUP_ZEROPT) {
+    for (int ic = 0; ic < Nphotcodes; ic++) {
+      double zpt = get_median_zpt_tgroups (photcodes[ic][0].code);
+      fprintf (stderr, "rationalize zero points by tgroup for %s (%d) : %f\n", photcodes[ic][0].name, photcodes[ic][0].code, zpt);
+      if (isnan(zpt)) continue;
+      set_median_zpt_tgroups (photcodes[ic][0].code, zpt);
+      set_median_zpt_mosaics (photcodes[ic][0].code, zpt);
+      set_median_zpt_images  (photcodes[ic][0].code, zpt);
+    }
+    return;
+  }
+
+  // if we are calculating zero points for mosaics, but not tgroups,
+  // rationalize based on the mosaics
+  if (MOSAIC_ZEROPT) {
+    for (int ic = 0; ic < Nphotcodes; ic++) {
+      double zpt = get_median_zpt_mosaics (photcodes[ic][0].code);
+      fprintf (stderr, "rationalize zero points by mosaic for %s (%d) : %f\n", photcodes[ic][0].name, photcodes[ic][0].code, zpt);
+      if (isnan(zpt)) continue;
+      set_median_zpt_tgroups (photcodes[ic][0].code, zpt);
+      set_median_zpt_mosaics (photcodes[ic][0].code, zpt);
+      set_median_zpt_images  (photcodes[ic][0].code, zpt);
+    }
+    return;
+  }    
+
+  // otherwise, rationalize based on the images
+  for (int ic = 0; ic < Nphotcodes; ic++) {
+    double zpt = get_median_zpt_images (photcodes[ic][0].code);
+    fprintf (stderr, "rationalize zero points by image for %s (%d) : %f\n", photcodes[ic][0].name, photcodes[ic][0].code, zpt);
+    if (isnan(zpt)) continue;
+    set_median_zpt_tgroups (photcodes[ic][0].code, zpt);
+    set_median_zpt_mosaics (photcodes[ic][0].code, zpt);
+    set_median_zpt_images  (photcodes[ic][0].code, zpt);
+  }
+}
+
+double get_median_zpt_images (short photcode) {
+
+  double *mlist;
+
+  ALLOCATE (mlist, double, Nimage);
+
+  int N = 0;
+  for (int i = 0; i < Nimage; i++) {
+    int mycode = GetPhotcodeEquivCodebyCode (image[i].photcode);
+    if (mycode != photcode) continue;
+    if (!(image[i].flags & ID_IMAGE_IMAGE_PHOTCAL)) continue;
+    mlist[N] = image[i].McalPSF;
+    N++;
+  }
+
+  if (N == 0) {
+    free (mlist);
+    return NAN;
+  }
+
+  StatType stats;
+  liststats_setmode (&stats, "MEAN");
+
+  liststats (mlist, NULL, NULL, N, &stats);
+  free (mlist);
+
+  return stats.median; 
+}
+
+void set_median_zpt_images (short photcode, double zpt) {
+
+  if (!isfinite(zpt)) return; // do not break the zero points
+
+  for (int i = 0; i < Nimage; i++) {
+    int mycode = GetPhotcodeEquivCodebyCode (image[i].photcode);
+    if (mycode != photcode) continue;
+    image[i].McalPSF -= zpt;
+    image[i].McalAPER -= zpt;
+  }
+
+  return; 
+}
+
 static int setMcal_init_done = FALSE;
 void plot_setMcal (double *list, int Npts) {