Changeset 41387 for branches/eam_branches/relphot.20200519/src/TGroupOps.c

Timestamp:

Jul 7, 2020, 9:30:05 AM (6 years ago)

Author:

eugene

Message:

updates to get tgroups working

File:

• branches/eam_branches/relphot.20200519/src/TGroupOps.c (modified) (22 diffs)

branches/eam_branches/relphot.20200519/src/TGroupOps.c

@@ -36 +36 @@
 // load the photometric nights from a file : these will be used to define the tgroups 
 // XXX call this in args.c where TGROUP_ZEROPT is set
+// XXX cannot call this in args.c : the Nphotcodes have not yet been loaded.
 void loadTGroups (char *filename) {
 
@@ -73 +74 @@
       tgroup[Nt].McalPSF   = 0.0; // note : at the end, tgroup.Mcal is added back to the input images
       tgroup[Nt].McalAPER  = 0.0; // note : tgroup stores only offsets relative to the original image values
+      tgroup[Nt].dKlam     = 0.0; 
       tgroup[Nt].dMcal     = 0.0; // note : at the end, tgroup.Mcal is added back to the input images
       tgroup[Nt].dMsys     = 0.0;
@@ -111 +113 @@
 
     // ignore non-GPC1 images
-    if (!isGPC1chip(subset[i].photcode)) {
-      Nsimple ++;
+    // if (!isGPC1chip(subset[i].photcode)) {
+    //   Nsimple ++;
+    //   continue;
+    // }
+
+    // ignore mosaic images (photcode == 0)
+    if (!subset[i].photcode) {
       continue;
     }
@@ -138 +145 @@
 
   fprintf (stderr, "matched %d images to %d tgroups, %d simple chips not matched to tgroups\n", (int) (Nsubset - Nsimple), (int) Ntgroup, (int) Nsimple);
+  return;
+}
+
+void freeTGroups (void) {
+
+  if (!TGROUP_ZEROPT) return;
+
+  free (ImageToTGroup);
+  free (tgroupTimes);
+
+  for (int i = 0; i < Ntgroup; i++) {
+    free (tgroup[i].image);
+    free (tgroup[i].measure);
+    free (tgroup[i].catalog);
+  }
+  free (tgroup);
   return;
 }
@@ -221 +244 @@
   // if start is not in this range, return -1
   if (start < tgroupTimes[Nlo]) return (-1);
-  if (start > tgroupTimes[Nhi]) return (-1);
+  if (start > tgroupTimes[Nhi] + 86399) return (-1); // 1-day range for tgroup for now
 
   while (Nhi - Nlo > 4) {
@@ -259 +282 @@
     for (j = 0; j < tgroup[i].Nimage; j++) {
       im = tgroup[i].image[j];
-      image[im].McalPSF    += tgroup[i].McalPSF;
-      image[im].McalAPER   += tgroup[i].McalAPER;
+      double Mgrp = (TGROUP_FIT_AIRMASS) ? tgroup[i].McalPSF + tgroup[i].dKlam*(image[im].secz - 1.0) : tgroup[i].McalPSF;
+      image[im].McalPSF    += Mgrp;
+      image[im].McalAPER   += Mgrp;
       image[im].dMcal       = tgroup[i].dMcal;
       image[im].McalChiSq   = tgroup[i].McalChiSq;
@@ -267 +291 @@
       image[im].flags      |= (tgroup[i].flags & ID_IMAGE_PHOTOM_FEW);
       image[im].flags      |= (tgroup[i].flags & ID_IMAGE_PHOTOM_POOR);
+
+      fprintf (stderr, "TG to IMAGE: %f -> %f (%d)\n", tgroup[i].McalPSF, image[im].McalPSF, tgroup[i].nFitPhotom);
     }
     tgroup[i].McalPSF  = 0.0;
     tgroup[i].McalAPER = 0.0;
+    tgroup[i].dKlam    = 0.0;
   }      
 }
@@ -287 +314 @@
 }
 
-void freeTGroupBins (int Ncatalog, int doTGroupList) {
+void freeTGroupBins (int Ncatalog) {
 
   off_t i;
@@ -329 +356 @@
     }
   }
-  // fprintf (stderr, "matched %d detections to tgroups\n", Nmatch);
+  fprintf (stderr, "Matched %d detections to tgroups\n", Nmatch);
+  for (i = 0; i < Ntgroup; i++) {
+    fprintf (stderr, "  TGROUP time %d, Nmeasure: %d, Nimage: %d\n", tgroup[i].start, (int) tgroup[i].Nmeasure, (int) tgroup[i].Nimage);
+  }
   return (TRUE);
 }
@@ -380 +410 @@
 }
 
-float getMgrp (off_t meas, int cat) {
+float getMgrp (off_t meas, int cat, float airmass) {
 
   off_t i;
@@ -392 +422 @@
 
   // if (tgroup[i].flags & IMAGE_BAD) return (NAN);  
-  value = tgroup[i].McalPSF;
+  value = tgroup[i].McalPSF + tgroup[i].dKlam*(airmass - 1.0);
   return (value);
 }
@@ -419 +449 @@
   int PoorImages;
   double *psfMagList;
+  double *psfAirList;
   double *psfErrList;
   double *brightMagList;
@@ -456 +487 @@
   if (allocLists) {
     ALLOCATE (info->psfMagList, double, Nmax);
+    ALLOCATE (info->psfAirList, double, Nmax);
     ALLOCATE (info->psfErrList, double, Nmax);
     ALLOCATE (info->kronMagList, double, Nmax);
@@ -466 +498 @@
 void SetMgrpInfoFree (SetMgrpInfo *info) {
   free (info->psfMagList   );
+  free (info->psfAirList   );
   free (info->psfErrList   );
   free (info->kronMagList  );
@@ -522 +555 @@
   image = getimages (&N, NULL);
 
-  fprintf (stderr, "limiting negative clouds to %f\n", CLOUD_TOLERANCE);
-
   if (PoorImages) {
     // XXX use bad stars and measurements for PoorImages? or not?
@@ -543 +574 @@
   for (i = 0; i < Ntgroup; i++) {
     setMgrp_tgroup (&tgroup[i], i, image, catalog, &info, flatcorr);
+    fprintf (stderr, "TGROUP Mcal: %f, dK: %f, limiting negative clouds to %f\n", tgroup[i].McalPSF, tgroup[i].dKlam, CLOUD_TOLERANCE);
   }
   SetMgrpInfoFree (&info);
@@ -575 +607 @@
 
   double *psfMagList    = info->psfMagList;
+  double *psfAirList    = info->psfAirList;
   double *psfErrList    = info->psfErrList;
   double *kronMagList   = info->kronMagList;
@@ -593 +626 @@
 
   int testImage = FALSE;
-  // testImage |= (abs(myTGroup[0].start - 1324104046) < 10);
+  testImage |= (abs(myTGroup[0].start - 1245283200) < 10);
   // testImage |= (abs(myTGroup[0].start - 1324103823) < 10);
 
@@ -633 +666 @@
       int isBad = (catalog[c].measureT[m].dbFlags & MEAS_BAD);
 
-      fprintf (fout, "%f %f : %f %f %f %f %f  : %f %d\n", catalog[c].averageT[n].R, catalog[c].averageT[n].D, MsysPSF, MrelPSF, Mcal, Mgrid, Mflat, delta, isBad);
+      fprintf (fout, "%f %f : %f %f %f %f %f  : %f %f %d\n", catalog[c].averageT[n].R, catalog[c].averageT[n].D, MsysPSF, MrelPSF, Mcal, Mgrid, Mflat, catalog[c].measureT[m].airmass, delta, isBad);
     }
 
@@ -697 +730 @@
 
     psfMagList[N]  = MsysPSF - MrelPSF - Mcal - Mgrid + Mflat;
+    psfAirList[N]  = (catalog[c].measureT[m].airmass - 1.0);
     psfErrList[N]  = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
     kronMagList[N] = MsysKron - MrelKron - Mcal - Mgrid + Mflat;
@@ -725 +759 @@
   }
 
-  liststats (psfMagList, psfErrList, NULL, N, &stats);
+  double dk = 0.0;
+  if (TGROUP_FIT_AIRMASS) {
+    liststats_fit1d (psfMagList, psfErrList, psfAirList, N, &stats, &dk);
+  } else {
+    liststats (psfMagList, psfErrList, NULL, N, &stats);
+  }
   if (VERBOSE2 && info->PoorImages) fprintf (stderr, "Mgrp: %f %f %d %d\n", stats.mean, stats.sigma, stats.Nmeas, N);
 
@@ -733 +772 @@
   myTGroup[0].McalChiSq  = stats.chisq;
   myTGroup[0].nFitPhotom = N;
+  myTGroup[0].dKlam      = dk;
 
   // drop this for tgroups?