Changeset 42818

Timestamp:

May 8, 2025, 4:07:43 PM (15 months ago)

Author:

eugene

Message:

merge from trunk: relphot grid improvements

Location:

branches/eam_branches/ipp-20230313/Ohana

Files:

• . (modified) (1 prop)
• src/addstar/src/grefstars.c (modified) (1 diff)
• src/relphot (modified) (1 prop)
• src/relphot/Makefile (modified) (1 diff)
• src/relphot/include/relphot.h (modified) (3 diffs)
• src/relphot/src/DumpAllMags.c (copied) (copied from trunk/Ohana/src/relphot/src/DumpAllMags.c )
• src/relphot/src/GridIO.c (modified) (1 diff)
• src/relphot/src/GridOps.c (modified) (8 diffs)
• src/relphot/src/MagResidSave.c (modified) (1 diff)
• src/relphot/src/relphot_images.c (modified) (1 diff)
• src/relphot/src/setMrelCatalog.c (modified) (1 diff)

branches/eam_branches/ipp-20230313/Ohana/src/addstar/src/grefstars.c

@@ -41 +41 @@
     catalog->average[N].R = ohana_normalize_angle (catalog->average[N].R);
 
+    catalog->measure[N].McalPSF  = 0.0;
+    catalog->measure[N].McalAPER = 0.0;
+    catalog->measure[N].dMcal    = 0.0;
+       
     catalog->measure[N].R = catalog->average[N].R;
     catalog->measure[N].D = catalog->average[N].D;

branches/eam_branches/ipp-20230313/Ohana/src/relphot/Makefile

@@ -82 +82 @@
 $(SRC)/setExclusions.$(ARCH).o   \
 $(SRC)/setMrelFinal.$(ARCH).o    \
-$(SRC)/setMrelCatalog.$(ARCH).o          \
-$(SRC)/BoundaryTreeOps.$(ARCH).o         \
+$(SRC)/setMrelCatalog.$(ARCH).o   \
+$(SRC)/BoundaryTreeOps.$(ARCH).o  \
+$(SRC)/DumpAllMags.$(ARCH).o      \
 $(SRC)/write_coords.$(ARCH).o
 

branches/eam_branches/ipp-20230313/Ohana/src/relphot/include/relphot.h

@@ -160 +160 @@
   float dY; // bin = int(Ychip * dY), dY = Ny / NyChip
   // NxChip, NyChip = 4900,4900 for now
+
+    int **nAlloc; // allocated vector length
+  double ***dMval; // values used to calculate corrections
 } GridCorrectionType;
 
@@ -615 +618 @@
 void          initGridBins        PROTO((void));
 GridCorrectionType *getGridCorrByCode PROTO((int code));
+GridCorrectionType *newGridCorrByCode PROTO((int code));
 GridCorrectionType *getGridCorrNext   PROTO((int *Nlast));
 int           GridCorrectionSave  PROTO((void));
@@ -902 +906 @@
 double get_median_zpt_tgroups (short photcode);
 void   set_median_zpt_tgroups (short photcode, double zpt);
+
+int DumpAllMags(char *filename, Catalog *catalog, int Ncatalog);

branches/eam_branches/ipp-20230313/Ohana/src/relphot/src/GridIO.c

@@ -151 +151 @@
 
     if (!gfits_scan (&header, "PHOTCODE", "%d", 1, &photcode))    myAbort ("failed to get PHOTCODE");
-    GridCorrectionType *GridCorr = getGridCorrByCode (photcode);
+    GridCorrectionType *GridCorr = newGridCorrByCode (photcode);
     GridCorr->photcode = photcode;
 

branches/eam_branches/ipp-20230313/Ohana/src/relphot/src/GridOps.c

@@ -76 +76 @@
       if (GridCorr[code]->dMgrid) FREE (GridCorr[code]->dMgrid[ix]);
       if (GridCorr[code]->nMgrid) FREE (GridCorr[code]->nMgrid[ix]);
+      if (GridCorr[code]->nAlloc) FREE (GridCorr[code]->nAlloc[ix]);
+
+      // dMval is a 2D array of vectors
+      for (int iy = 0; iy < GridCorr[code]->Ny; iy++) {
+        FREE (GridCorr[code]->dMval[ix][iy]);
+      }
+      FREE (GridCorr[code]->dMval[ix]);
     }      
 
@@ -81 +88 @@
     FREE (GridCorr[code]->dMgrid);
     FREE (GridCorr[code]->nMgrid);
+    FREE (GridCorr[code]->nAlloc);
+    FREE (GridCorr[code]-> dMval);
 
     FREE(GridCorr[code]);
@@ -109 +118 @@
 }
 
-GridCorrectionType *getGridCorrByCode (int code) {
+// generate a new GridCorr structure in the array
+GridCorrectionType *newGridCorrByCode (int code) {
 
   if (GridCorr == NULL) return NULL;
@@ -116 +126 @@
   if (code <          0) return NULL;
 
+  myAssert (GridCorr[code] == NULL, "oops, should not already be allocated!");
   ALLOCATE(GridCorr[code], GridCorrectionType, 1)
+  GridCorrectionType *result = GridCorr[code];
+  return result;
+}
+
+// return the GridCode structure corresponding to 'code'
+GridCorrectionType *getGridCorrByCode (int code) {
+
+  if (GridCorr == NULL) return NULL;
+
+  if (code >= NGridCorr) return NULL;
+  if (code <          0) return NULL;
+
   GridCorrectionType *result = GridCorr[code];
   return result;
@@ -163 +186 @@
     // we are normally accessing this array randomly, so there is no advantage to
     // doing Mgrid[y][x] vs Mgrid[x][y]
+    ALLOCATE (GridCorr[code]-> dMval, double **, NX_BIN);
     ALLOCATE (GridCorr[code]-> Mgrid, float *, NX_BIN);
     ALLOCATE (GridCorr[code]->dMgrid, float *, NX_BIN);
     ALLOCATE (GridCorr[code]->nMgrid,   int *, NX_BIN);
+    ALLOCATE (GridCorr[code]->nAlloc,   int *, NX_BIN);
     for (int ix = 0; ix < NX_BIN; ix++) {
       ALLOCATE (GridCorr[code]-> Mgrid[ix], float, NY_BIN);
       ALLOCATE (GridCorr[code]->dMgrid[ix], float, NY_BIN);
       ALLOCATE (GridCorr[code]->nMgrid[ix],   int, NY_BIN);
+      ALLOCATE (GridCorr[code]->nAlloc[ix],   int, NY_BIN);
+
+      // when we first allocate this array, set the elements to 0 (NULL)
+      // these arrays will be allocated or reallocated in resetMgrid
+      ALLOCATE_ZERO (GridCorr[code]-> dMval[ix], double *, NX_BIN);
     }      
   }
@@ -190 +220 @@
         GridCorr[code]->dMgrid[ix][iy] = 0.0;
         GridCorr[code]->nMgrid[ix][iy] =   0;
+        GridCorr[code]->nAlloc[ix][iy] = 100;
+
+        // allocate or reset vector length to default
+        if (GridCorr[code]->dMval[ix][iy]) {
+          REALLOCATE (GridCorr[code]->dMval[ix][iy], double, GridCorr[code]->nAlloc[ix][iy]);
+        } else {
+          ALLOCATE (GridCorr[code]->dMval[ix][iy], double, GridCorr[code]->nAlloc[ix][iy]);
+        }
       }
     }      
@@ -263 +301 @@
         if (fabs(dM) > 0.5) continue;
 
-        grid-> Mgrid[ix][iy] += dM;
-        grid->dMgrid[ix][iy] += dM*dM;
+        int ival = grid->nMgrid[ix][iy];
+        grid-> dMval[ix][iy][ival] = dM;
         grid->nMgrid[ix][iy] ++;
+        CHECK_REALLOCATE (grid-> dMval[ix][iy], double, grid->nAlloc[ix][iy], grid->nMgrid[ix][iy], 100);
+
+        /* XXX old code to calculate simple means:
+           grid-> Mgrid[ix][iy] += dM;
+           grid->dMgrid[ix][iy] += dM*dM;
+           grid->nMgrid[ix][iy] ++;
+        */
       }
     }
   }
+
+  StatType stats;
+  liststats_setmode (&stats, "MEDIAN");
+
+  // now calculate Mgrid, dMgrid, nMgrid from Sum, Sum, Npt
+  for (int code = 0; code < NGridCorr; code++) {
+    if (!GridCorr[code]) continue;
+
+    for (int ix = 0; ix < GridCorr[code]->Nx; ix++) {
+      for (int iy = 0; iy < GridCorr[code]->Ny; iy++) {
+
+        // cells without sufficient coverage stay at 0.0
+        if (GridCorr[code]->nMgrid[ix][iy] < 5) {
+          GridCorr[code]-> Mgrid[ix][iy] = 0.0;
+          GridCorr[code]->dMgrid[ix][iy] = 0.0;
+          continue;
+        }
+
+        liststats_init (&stats);
+
+        liststats (GridCorr[code]->dMval[ix][iy], NULL, NULL, GridCorr[code]->nMgrid[ix][iy], &stats);
+        double altSigma = (stats.Upper80 - stats.Lower20) / 1.6;  // 20% to 80% encompasses 60% of the values, corresponds to the range (-0.85 sigma : +0.85 sigma)
+
+        GridCorr[code]-> Mgrid[ix][iy] = stats.median;
+        GridCorr[code]->dMgrid[ix][iy] = altSigma; // sample stdev
+      }
+    }
+  }
+  return;
+}
+
+# if (0)
 
   // now calculate Mgrid, dMgrid, nMgrid from Sum, Sum, Npt
@@ -300 +377 @@
     // float GridAve = GridSum / GridCnt;
     // fprintf (stderr, "grid average: %f\n", GridAve);
-  }
-  return;
-}
+# endif
 
 float getMgrid (Measure *measure) {

branches/eam_branches/ipp-20230313/Ohana/src/relphot/src/MagResidSave.c

@@ -54 +54 @@
   gfits_create_table_header (&theader, "BINTABLE", "MAG_PSF_RESID");
 
-  gfits_define_bintable_column (&theader, "D", "RA",       "ra",                         "degree", 1.0, 0.0);
-  gfits_define_bintable_column (&theader, "D", "DEC",      "dec",                        "degree", 1.0, 0.0);
-  gfits_define_bintable_column (&theader, "E", "MAG_RESID", "magnitude residual", "magnitudes", 1.0, 0.0);
-  gfits_define_bintable_column (&theader, "E", "MAG_ERROR", "magnitude error", "magnitudes", 1.0, 0.0);
-  gfits_define_bintable_column (&theader, "E", "MAG_CHISQ", "magnitude resid chisq", "unitless", 1.0, 0.0);
-  gfits_define_bintable_column (&theader, "E", "MAG_MEDIAN", "median magnitude residual", "magnitudes", 1.0, 0.0);
+  gfits_define_bintable_column (&theader, "D", "RA",         "ra",                           "degree",     1.0, 0.0);
+  gfits_define_bintable_column (&theader, "D", "DEC",        "dec",                          "degree",     1.0, 0.0);
+  gfits_define_bintable_column (&theader, "E", "MAG_RESID",  "magnitude residual",           "magnitudes", 1.0, 0.0);
+  gfits_define_bintable_column (&theader, "E", "MAG_ERROR",  "magnitude error",              "magnitudes", 1.0, 0.0);
+  gfits_define_bintable_column (&theader, "E", "MAG_CHISQ",  "magnitude resid chisq",        "unitless",   1.0, 0.0);
+  gfits_define_bintable_column (&theader, "E", "MAG_MEDIAN", "median magnitude residual",    "magnitudes", 1.0, 0.0);
   gfits_define_bintable_column (&theader, "E", "MAG_ALTSIG", "68 percentile range as sigma", "magnitudes", 1.0, 0.0);
-  gfits_define_bintable_column (&theader, "E", "AIRMASS",   "", "unitless", 1.0, 0.0);
-  gfits_define_bintable_column (&theader, "D", "MJD",       "", "unitless", 1.0, 0.0);
-  gfits_define_bintable_column (&theader, "E", "EXPTIME",   "", "unitless", 1.0, 0.0);
-  gfits_define_bintable_column (&theader, "J", "PHOTCODE",  "", "unitless", 1.0, 0.0);
-  gfits_define_bintable_column (&theader, "J", "FLAGS",     "image flags", "unitless",   1.0, 0.0);
-  gfits_define_bintable_column (&theader, "J", "NFIT",      "number of stars fitted", NULL, 1.0, FT_BZERO_INT32);
-  gfits_define_bintable_column (&theader, "J", "NMEAS",     "number of stars on image", NULL, 1.0, FT_BZERO_INT32);
+  gfits_define_bintable_column (&theader, "E", "AIRMASS",    "",                             "unitless",   1.0, 0.0);
+  gfits_define_bintable_column (&theader, "D", "MJD",        "",                             "unitless",   1.0, 0.0);
+  gfits_define_bintable_column (&theader, "E", "EXPTIME",    "",                             "unitless",   1.0, 0.0);
+  gfits_define_bintable_column (&theader, "J", "PHOTCODE",   "",                             "unitless",   1.0, 0.0);
+  gfits_define_bintable_column (&theader, "J", "FLAGS",      "image flags",                  "unitless",   1.0, 0.0);
+  gfits_define_bintable_column (&theader, "J", "NFIT",       "number of stars fitted",       NULL,         1.0, 0.0);
+  gfits_define_bintable_column (&theader, "J", "NMEAS",      "number of stars on image",     NULL,         1.0, 0.0);
 
   // generate the output array that carries the data

branches/eam_branches/ipp-20230313/Ohana/src/relphot/src/relphot_images.c

@@ -147 +147 @@
     MagResidSave ("mag.resid.fits", catalog);
 
+    // DumpAllMags ("all.mags.fits", catalog, Ncatalog);
+
     save_images_updates (&db);
 

branches/eam_branches/ipp-20230313/Ohana/src/relphot/src/setMrelCatalog.c

@@ -412 +412 @@
       secfilt[Nsec].flags |= ID_SECF_HAS_CFH;
     }
+//  secfilt[Nsec].flags |= ID_SECF_AST_GAIA ; // for building reference catalog only, STAR_FEW need to comment out
+//  secfilt[Nsec].flags |= ID_SECF_PHO_SYNTH; // for building reference catalog only, STAR_FEW need to comment out
+//  secfilt[Nsec].flags |= ID_SECF_PHO_ATLAS; // for building reference catalog only, STAR_FEW need to comment out
+//  secfilt[Nsec].flags |= ID_SECF_PHO_PV3  ; // for building reference catalog only, STAR_FEW need to comment out
+//  secfilt[Nsec].flags |= ID_SECF_PHO_SOUTH; // for building reference catalog only, STAR_FEW need to comment out
 
     // if -preserve-ps1 is set and this object has PS1 data, skip the rest of the steps: