Changeset 41662

Timestamp:

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

Author:

eugene

Message:

add code to rationalize the zero points on each pass

Location:

trunk/Ohana/src/relphot

Files:

• include/relphot.h (modified) (1 diff)
• src/ImageOps.c (modified) (1 diff)
• src/MosaicOps.c (modified) (1 diff)
• src/TGroupOps.c (modified) (1 diff)
• src/relphot_images.c (modified) (1 diff)

trunk/Ohana/src/relphot/include/relphot.h

@@ -874 +874 @@
 void ResetAverageActivePhotcodes (SecFilt *secfilt);
 
-
+void   rationalize_zeropoints (int Niter);
+double get_median_zpt_images (short photcode);
+void   set_median_zpt_images (short photcode, double zpt);
+double get_median_zpt_mosaics (short photcode);
+void   set_median_zpt_mosaics (short photcode, double zpt);
+double get_median_zpt_tgroups (short photcode);
+void   set_median_zpt_tgroups (short photcode, double zpt);

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) {

trunk/Ohana/src/relphot/src/MosaicOps.c

@@ -1574 +1574 @@
 }
 
+double get_median_zpt_mosaics (short photcode) {
+
+  double *mlist;
+
+  if (!MOSAIC_ZEROPT) return NAN;
+
+  ALLOCATE (mlist, double, Nmosaic);
+
+  int N = 0;
+  for (int i = 0; i < Nmosaic; i++) {
+    if (mosaic[i].photcode != photcode) continue;
+    if (!(mosaic[i].flags & ID_IMAGE_MOSAIC_PHOTCAL)) continue;
+    mlist[N] = mosaic[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_mosaics (short photcode, double zpt) {
+
+  if (!MOSAIC_ZEROPT) return;
+  if (!isfinite(zpt)) return; // do not break the zero points
+
+  for (int i = 0; i < Nmosaic; i++) {
+    if (mosaic[i].photcode != photcode) continue;
+    mosaic[i].McalPSF -= zpt;
+    mosaic[i].McalAPER -= zpt;
+  }
+
+  return; 
+}
+
 void plot_mosaic_fields (Catalog *catalog) {
 

trunk/Ohana/src/relphot/src/TGroupOps.c

@@ -1365 +1365 @@
 }
 
+double get_median_zpt_tgroups (short photcode) {
+
+  double *mlist;
+
+  if (!TGROUP_ZEROPT) return NAN;
+
+  ALLOCATE (mlist, double, NtgroupTimes);
+
+  int N = 0;
+  for (int i = 0; i < NtgroupTimes; i++) {
+    TGroup *tgroup = tgroupTimes[i][0].byCode;
+    for (int j = 0; j < tgroupTimes[i][0].nCode; j++) {
+      if (tgroup[j].photcode != photcode) continue;
+      if (!(tgroup[j].flags & ID_IMAGE_TGROUP_PHOTCAL)) continue;
+      mlist[N] = tgroup[j].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_tgroups (short photcode, double zpt) {
+
+  if (!TGROUP_ZEROPT) return;
+  if (!isfinite(zpt)) return; // do not break the zero points
+
+  for (int i = 0; i < NtgroupTimes; i++) {
+    TGroup *tgroup = tgroupTimes[i][0].byCode;
+    for (int j = 0; j < tgroupTimes[i][0].nCode; j++) {
+      if (tgroup[j].photcode != photcode) continue;
+      tgroup[j].McalPSF -= zpt;
+      tgroup[j].McalAPER -= zpt;
+    }
+  }
+  return; 
+}
+
 void plot_tgroup_fields (Catalog *catalog) {
 

trunk/Ohana/src/relphot/src/relphot_images.c

@@ -102 +102 @@
       setMmos  (catalog);
       setMgrp  (catalog);
+
+      rationalize_zeropoints (i);
 
       setMgrid (catalog, Ncatalog);