Changeset 41647

Timestamp:

Jun 7, 2021, 4:14:19 PM (5 years ago)

Author:

eugene

Message:

major work to revise the grid correction : do not use astrometry to define grid, just subdivide each photcode

Location:

trunk/Ohana/src/relphot

Files:

• . (modified) (1 prop)
• Makefile (modified) (4 diffs)
• include/relphot.h (modified) (28 diffs)
• src/BrightCatalog.c (modified) (4 diffs)
• src/ConfigInit.c (modified) (3 diffs)
• src/GridIO.c (copied) (copied from branches/eam_branches/relphot.20210521/src/GridIO.c )
• src/GridOps.c (modified) (1 diff)
• src/ImageOps.c (modified) (22 diffs)
• src/MagResidSave.c (modified) (5 diffs)
• src/MosaicOps.c (modified) (36 diffs)
• src/RepairWarpMeasures.c (deleted)
• src/ResetOps.c (copied) (copied from branches/eam_branches/relphot.20210521/src/ResetOps.c )
• src/Shutdown.c (modified) (1 diff)
• src/StackImageMaps.c (deleted)
• src/StarOps.c (modified) (19 diffs)
• src/TGroupOps.c (modified) (27 diffs)
• src/WarpImageMaps.c (deleted)
• src/ZeroPointModes.c (modified) (6 diffs)
• src/args.c (modified) (10 diffs)
• src/assign_images.c (modified) (2 diffs)
• src/bcatalog.c (modified) (4 diffs)
• src/extra.c (modified) (3 diffs)
• src/global_stats.c (modified) (4 diffs)
• src/initialize.c (modified) (2 diffs)
• src/launch_region_hosts.c (modified) (1 diff)
• src/load_catalogs.c (modified) (1 diff)
• src/load_images.c (modified) (2 diffs)
• src/plot_scatter.c (modified) (2 diffs)
• src/reload_catalogs.c (modified) (6 diffs)
• src/relphot.c (modified) (4 diffs)
• src/relphot_client.c (modified) (4 diffs)
• src/relphot_images.c (modified) (12 diffs)
• src/relphot_objects.c (modified) (2 diffs)
• src/relphot_parallel_images.c (modified) (8 diffs)
• src/relphot_parallel_regions.c (modified) (2 diffs)
• src/select_images.c (modified) (1 diff)
• src/setMrelCatalog.c (modified) (16 diffs)
• src/setMrelFinal.c (modified) (6 diffs)
• test (modified) (1 prop)
• test/relphot.grid.dvo (copied) (copied from branches/eam_branches/relphot.20210521/test/relphot.grid.dvo )
• test/relphot.samples.cmd (copied) (copied from branches/eam_branches/relphot.20210521/test/relphot.samples.cmd )
• test/sample.gdb (copied) (copied from branches/eam_branches/relphot.20210521/test/sample.gdb )

trunk/Ohana/src/relphot/Makefile

@@ -30 +30 @@
 $(SRC)/ConfigInit.$(ARCH).o      \
 $(SRC)/GridOps.$(ARCH).o         \
+$(SRC)/GridIO.$(ARCH).o  \
+$(SRC)/ResetOps.$(ARCH).o        \
 $(SRC)/ImageOps.$(ARCH).o        \
 $(SRC)/ImageSubset.$(ARCH).o     \
@@ -39 +41 @@
 $(SRC)/Shutdown.$(ARCH).o        \
 $(SRC)/StarOps.$(ARCH).o         \
-$(SRC)/WarpImageMaps.$(ARCH).o   \
-$(SRC)/StackImageMaps.$(ARCH).o  \
-$(SRC)/RepairWarpMeasures.$(ARCH).o      \
-$(SRC)/myIndex.$(ARCH).o         \
 $(SRC)/psps_ids.$(ARCH).o        \
 $(SRC)/args.$(ARCH).o            \
@@ -93 +91 @@
 $(SRC)/ConfigInit.$(ARCH).o      \
 $(SRC)/GridOps.$(ARCH).o         \
+$(SRC)/GridIO.$(ARCH).o  \
+$(SRC)/ResetOps.$(ARCH).o        \
 $(SRC)/ImageOps.$(ARCH).o        \
 $(SRC)/ImageSubset.$(ARCH).o     \
@@ -102 +102 @@
 $(SRC)/Shutdown.$(ARCH).o        \
 $(SRC)/StarOps.$(ARCH).o         \
-$(SRC)/WarpImageMaps.$(ARCH).o   \
-$(SRC)/StackImageMaps.$(ARCH).o  \
-$(SRC)/RepairWarpMeasures.$(ARCH).o      \
-$(SRC)/myIndex.$(ARCH).o         \
 $(SRC)/psps_ids.$(ARCH).o        \
 $(SRC)/bcatalog.$(ARCH).o        \

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

@@ -5 +5 @@
 # include <pthread.h>
 
-/* # define GRID_V1 */
-# define GRID_V2
-# define NO_IMAGE -100
-
 # define ID_SECF_STACK_PRIMARY 0x00004000
 
 # define NBOOTSTRAP 100
+
+// MEAS_BAD is used in ImageOps.c, MosaicOps.c, TGroupOps.c to skip specific measurements.
+// ID_MEAS_AREA is raised for measurements outside the user-specified region of the chips
+// ID_MEAS_NOCAL is raised for measurements outside the valid time range or for inactive photcodes
+// Note: these bits are only raised on the temporary MeasureTiny elements and are not saved
+// to the database files.
+# define MEAS_BAD (ID_MEAS_AREA | ID_MEAS_NOCAL)
 
 # ifndef MAX_INT
@@ -53 +56 @@
 
 typedef enum {
-    MODE_NONE   = 0,
-    MODE_LOAD   = 1,
-    MODE_UPDATE = 2,
+    MODE_NONE           = 0,
+    MODE_LOAD           = 1,
+    MODE_UPDATE         = 2,
     MODE_UPDATE_OBJECTS = 3,
-    MODE_SYNTH_PHOT = 4,
+    MODE_SYNTH_PHOT     = 4,
 } ModeType;
 
@@ -70 +73 @@
   PS1_w = 6,
 } PS1filters;
+
+typedef enum {
+  GRID_MEAN = 0,
+  GRID_STDEV = 1,
+  GRID_NPTS = 2,
+} GridOutputMode;
 
 typedef struct {
@@ -99 +108 @@
   char skipCal;               // if TRUE, this mosaic is incomplete and should not be calibrated
   char inTGroup;
-  Coords coords;
+  Coords coords;              // this is only used to set the mosaic center for assign_images used by region hosts
 } Mosaic; 
 
@@ -140 +149 @@
   int      nCode;
 } TGTimes;
+
+typedef struct {
+  unsigned short photcode;
+  float **Mgrid; // grid of average corrections
+  float **dMgrid; // grid of stdev of corrections
+    int **nMgrid; // grid of number of stars for corrections
+  int Nx; // bin = int(Xchip * (Nx / NxChip))
+  int Ny; // bin = int(Ychip * (Ny / NyChip))
+  float dX; // bin = int(Xchip * dX), dX = Nx / NxChip
+  float dY; // bin = int(Ychip * dY), dY = Ny / NyChip
+  // NxChip, NyChip = 4900,4900 for now
+} GridCorrectionType;
 
 typedef enum {
@@ -366 +387 @@
 int    SKY_DEPTH;  /** XXX EAM : depth of catalog tables, fix usage */
 char  *SYNTH_ZERO_POINTS;
+char  *GRID_MEANFILE;
 
 // globals for parallel region operations
@@ -389 +411 @@
 double IMAGE_OFFSET;
 double NIGHT_SCATTER;
-double NIGHT_OFFSET;
+double NIGHT_CHISQ;
+double MOSAIC_SCATTER;
+double MOSAIC_CHISQ;
 double STAR_SCATTER;
 double STAR_CHISQ;
@@ -413 +437 @@
 int    MOSAIC_ZEROPT;
 int    TGROUP_ZEROPT;
+int    GRID_ZEROPT;
+
+int    GRID_BIN_HSC;
+int    GRID_BIN_CFH;
+int    GRID_BIN_GPC1;
+int    GRID_BIN_GPC2;
 
 int    FREEZE_IMAGES;
@@ -424 +454 @@
 int    MANUAL_ITERATION;
 int    NLOOP;
-int    NGRID;
 int    RESET;
 int    RESET_ZEROPTS;
+int    RESET_FLATCORR;
 int    REPAIR_WARPS;
 int    PRESERVE_PS1;
+int    REQUIRE_PSFFIT;
 int    UPDATE;
 int    SAVE_IMAGE_UPDATES;
@@ -436 +467 @@
 char   MOSAICNAME[256];
 char   STATMODE[32];
-int    STAR_BAD;
-int    MEAS_BAD;
 int    STAR_TOOFEW;
 int    GRID_TOOFEW;
 int    IMAGE_TOOFEW;
 double IMAGE_GOOD_FRACTION;
-int    IMAGE_BAD;
 int    CALIBRATE_STACKS_AND_WARPS;
-int    USE_GRID;
+
 int    KEEP_UBERCAL;
 char  *OUTROOT;
@@ -457 +485 @@
 char  *PhotcodeList;
 
-int    RELPHOT_GRID_X;
-int    RELPHOT_GRID_Y;
-int    RELPHOT_GRID_BINNING;
-
 int      *photseclist;
 int      Nphotcodes;
@@ -490 +514 @@
 char *UserCatalog;
 
+enum {GRID_ZPT_MODE_NONE, GRID_ZPT_MODE_ALL, };
+int GRID_ZPT_MODE;
+
 enum {TGROUP_ZPT_MODE_NONE, TGROUP_ZPT_MODE_GOOD_NIGHT, TGROUP_ZPT_MODE_ALL, };
 int TGROUP_ZPT_MODE;
@@ -504 +531 @@
 
 RelphotStages STAGES;
-
-# ifdef GRID_V1
-int setGridMeasure (off_t meas, int cat, double X, double Y);
-# endif
-
-# ifdef GRID_V2
-int setGridMeasure (off_t meas, int cat, double X, double Y, int ccdnum);
-# endif
 
 /*** relphot prototypes ***/
@@ -523 +542 @@
 int           bcatalog            PROTO((Catalog *subcatalog, Catalog *catalog, int Ncat));
 void          clean_images        PROTO((void));
-void          clean_measures      PROTO((Catalog *catalog, int Ncatalog, int final, FlatCorrectionTable *flatcorr));
+void          clean_measures      PROTO((Catalog *catalog, int Ncatalog, int final));
 void          clean_mosaics       PROTO((void));
 void          clean_tgroups       PROTO((void));
@@ -549 +568 @@
 int SetSignals (void);
 
-void          freeGridBins        PROTO((int Ncatalog));
+void          freeGridBins        PROTO((void));
 void          freeImageBins       PROTO((int Ncatalog, int doImageList));
 void          freeMosaicBins      PROTO((int Ncatalog, int doMosaicList));
@@ -560 +579 @@
 
 float         getMcal             PROTO((off_t meas, int cat, dvoMagClassType class));
-float         getMflat            PROTO((off_t meas, int cat, FlatCorrectionTable *flatcorr, Catalog *catalog));
-float         getMgrid            PROTO((off_t meas, int cat));
+float         getMflat            PROTO((off_t meas, int cat, Catalog *catalog));
+float         getMgridTiny        PROTO((MeasureTiny *measure));
+float         getMgrid            PROTO((Measure *measure));
 float         getMmos             PROTO((off_t meas, int cat));
 float         getMgrp             PROTO((off_t meas, int cat, float airmass, float *dZpt));
@@ -570 +590 @@
 Image        *getimages           PROTO((off_t *N, off_t **LineNumber));
 ImageSubset  *getimages_subset    PROTO((off_t *N));
-void          global_stats        PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr, int nloop));
-void          initGrid            PROTO((int dX, int dY));
-void          initGridBins        PROTO((Catalog *catalog, int Ncatalog));
+void          global_stats        PROTO((Catalog *catalog, int Ncatalog, int nloop));
+void          initGrid            PROTO((void));
+void          setMflatFromGrid    PROTO((Catalog *catalog));
+void          initGridBins        PROTO((void));
+GridCorrectionType *getGridCorrByCode PROTO((int code));
+GridCorrectionType *getGridCorrNext   PROTO((int *Nlast));
+int           GridCorrectionSave  PROTO((void));
+int           GridCorrectionSaveFile  PROTO((char *filename, GridOutputMode mode));
+void          GridCorrectionLoad  PROTO((char *filename));
 void          initImageBins       PROTO((Catalog *catalog, int Ncatalog, int doImageList));
 void          initImagesSubset    PROTO((ImageSubset *input, off_t *line_number, off_t N));
 void          initImages          PROTO((Image *input, off_t *LineNumber, off_t N));
 void          initMosaicBins      PROTO((Catalog *catalog, int Ncatalog, int doMosaicList));
-void          initMosaicGrid      PROTO((Image *image, off_t Nimage));
+void          initMosaicMcal      PROTO((Image *image, off_t Nimage));
 void          initMosaics         PROTO((Image *subset, off_t Nsubset, Image *image, char *inSubset, off_t Nimage));
 void          initTGroups         PROTO((Image *subset, off_t Nsubset));
@@ -621 +647 @@
 void          plot_chisq          PROTO((Catalog *catalog, int Ncatalog));
 void          plot_defaults       PROTO((Graphdata *graphdata));
-void          plot_grid           PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr));
 void          plot_images         PROTO((void));
 void          plot_list           PROTO((Graphdata *graphdata, double *xlist, double *ylist, int N, char *label, char *format, ...) OHANA_FORMAT(printf, 6, 7) );
 void          plot_mosaic_fields  PROTO((Catalog *catalog));
 void          plot_mosaics        PROTO((void));
-void          plot_scatter        PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr));
+void          plot_scatter        PROTO((Catalog *catalog, int Ncatalog));
 void          plot_star_coords    PROTO((Catalog *catalog, int Ncatalog));
 void          plot_stars          PROTO((Catalog *catalog, int Ncatalog));
@@ -635 +660 @@
 int           open_graph          PROTO((int N));
 
-void          reload_catalogs     PROTO((SkyList *skylist, FlatCorrectionTable *flatcorr, int hostID, char *hostpath));
+void          reload_catalogs     PROTO((SkyList *skylist, int hostID, char *hostpath));
 int           reload_catalogs_parallel PROTO((SkyList *sky));
 int           reload_images       PROTO((FITS_DB *db));
 int           setExclusions       PROTO((Catalog *catalog, int Ncatalog, int verbose));
-void          setMgrid            PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr));
+void          setMgrid            PROTO((Catalog *catalog, int Ncatalog));
+void          resetMgrid          PROTO((void));
+
 void          setMcalFromMosaics  PROTO((void));
 void          setMcalFromTGroups  PROTO((void));
-int           setMcalOutput       PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr));
-
-void          setMcal             PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr));
-int           setMmos             PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr));
-int           setMgrp             PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr));
-
-int           setMrel             PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr));
-void          setMrelFinal        PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr, int simpleAverage));
-int           setMrelOutput       PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr));
-
-void          setMcalTest         PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr));
+int           setMcalOutput       PROTO((Catalog *catalog, int Ncatalog));
+
+void          setMcal             PROTO((Catalog *catalog));
+int           setMmos             PROTO((Catalog *catalog));
+int           setMgrp             PROTO((Catalog *catalog));
+
+int           setMrel             PROTO((Catalog *catalog, int Ncatalog));
+void          setMrelFinal        PROTO((Catalog *catalog, int simpleAverage));
+int           setMrelOutput       PROTO((Catalog *catalog, int Ncatalog));
+
+void          setMcalTest         PROTO((Catalog *catalog));
 
 int           setMave             PROTO((Catalog *catalog, int Ncatalog));
 void          set_ZP              PROTO((double ZERO));
 int           setrefcode          PROTO((Image *image, off_t Nimage)); 
-void          skip_measurements   PROTO((Catalog *catalog, int pass, FlatCorrectionTable *flatcorr));
+void          skip_measurements   PROTO((Catalog *catalog, int pass));
 void          sortA               PROTO((double *X, int N));
 void          sortB               PROTO((double *X, double *Y, int N));
@@ -663 +690 @@
 void          sortD               PROTO((double *X, double *Y, double *Z, int N));
 StatType      statsImageM         PROTO((Catalog *catalog));
-StatType      statsImageN         PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr));
+StatType      statsImageN         PROTO((Catalog *catalog));
 StatType      statsImageX         PROTO((Catalog *catalog));
 StatType      statsImagedM        PROTO((Catalog *catalog));
 StatType      statsMosaicM        PROTO((Catalog *catalog));
-StatType      statsMosaicN        PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr));
+StatType      statsMosaicN        PROTO((Catalog *catalog));
 StatType      statsMosaicX        PROTO((Catalog *catalog));
 StatType      statsMosaicdM       PROTO((Catalog *catalog));
@@ -674 +701 @@
 StatType      statsTGroupdM       PROTO((Catalog *catalog));
 
-StatType      statsStarN          PROTO((Catalog *catalog, int Ncatalog, int Nsec, int seccode, FlatCorrectionTable *flatcorr));
+StatType      statsStarN          PROTO((Catalog *catalog, int Ncatalog, int Nsec, int seccode));
 StatType      statsStarS          PROTO((Catalog *catalog, int Ncatalog, int Nsec));
 StatType      statsStarX          PROTO((Catalog *catalog, int Ncatalog, int Nsec));
@@ -720 +747 @@
 
 int print_measure_set_alt (Average *average, SecFilt *secfilt, Measure *measure);
-int setMrel_catalog_alt (Catalog *catalog, int Nc, int isSetMrelFinal, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt);
-int setMrelAverageExposure (Catalog *catalog, int cat, off_t ave, int Nsecfilt, int isSetMrelFinal, FlatCorrectionTable *flatcorr, SetMrelInfo *results);
-int setMrelAverageStack (Catalog *catalog, int cat, off_t ave, int Nsecfilt, FlatCorrectionTable *flatcorr);
-int setMrelAverageForcedWarp (Catalog *catalog, int cat, off_t ave, int Nsecfilt, FlatCorrectionTable *flatcorr, SetMrelInfo *results);
+
+// in setMrelCatalog.c:
+int setMrelCatalog (Catalog *catalog, int Nc, int isSetMrelFinal, SetMrelInfo *results, int Nsecfilt);
+int setMrelAverageExposure (Catalog *catalog, int cat, off_t ave, int Nsecfilt, int isSetMrelFinal, SetMrelInfo *results);
+int setMrelAverageStack (Catalog *catalog, int cat, off_t ave, int Nsecfilt);
+int setMrelAverageForcedWarp (Catalog *catalog, int cat, off_t ave, int Nsecfilt, SetMrelInfo *results);
 
 int setGlobalObjStats (Average *average, Measure *measure);
@@ -735 +764 @@
 StatDataSet *StatDataSetAlloc (int Nsecfilt, int Nmax);
 
-int init_synthetic_mags ();
+int init_synthetic_mags (void);
 int add_synthetic_mags (AverageTiny *average, SecFilt *secfilt, MeasureTiny *measure, off_t *Nmeasure, off_t *Nm);
 
@@ -796 +825 @@
 
 int SynthZeroPointsLoad (char *filename);
-SynthZeroPoints *SynthZeroPointsGet ();
+SynthZeroPoints *SynthZeroPointsGet (void);
 
 int relphot_synthphot (SkyList *sky, int hostID, char *hostpath);
@@ -815 +844 @@
 
 void put_astrom_table (AstromOffsetTable *myTable);
-void free_astrom_table ();
-
-int RepairWarpMeasures (Catalog *catalog);
-
-int FindWarpGroups (void);
-int GetWarpSeq (Image *image, int obstime, unsigned short photcode, double Rave, double Dave, float X, float Y);
-void FreeWarpGroups (void);
-
-myIndexType *myIndexInit ();
-int myIndexFree (myIndexType *myIndex);
-int myIndexUpdateLimits (myIndexType *myIndex, int value);
-int myIndexSetRange (myIndexType *myIndex);
-int myIndexSetEntry (myIndexType *myIndex, int value, int entry);
-int myIndexGetEntry (myIndexType *myIndex, int value);
+void free_astrom_table (void);
 
 uint64_t CreatePSPSObjectID(double ra, double dec);
@@ -834 +850 @@
 uint64_t CreatePSPSDetectionID(double tobs, int ccdid, int detID);
 
-int MakeStackIndex (void);
-void FreeStackGroups (void);
-int GetStackSeq (Image *image, double Rstk, double Dstk, unsigned short photcode, float X, float Y);
 void sort_by_ra (double *R, double *D, int *I, int *S, int N);
-
-void RepairWarpMeasuresOpenLogfile ();
-void RepairWarpMeasuresCloseLogfile ();
 
 void dump_tgroups (Catalog *catalog, int Npass);
@@ -849 +859 @@
 int GetZptIteration (void);
 
-void SetZeroPointModes (void);
+void SetZeroPointModes (Catalog *catalog, int Ncatalog);
 int UseStandardOLS (ZptFitModeType mode);
 int GetActivePhotcodeIndex (int photcode);
@@ -857 +867 @@
 
 int MagResidSave(char *filename, Catalog *catalog);
+
+void ResetAverageAndMeasure (Catalog *catalog);
+void ResetAverageObjects (Catalog *catalog);
+void ResetImages (Image *subset, off_t Nsubset);
+void ResetMeasureZeroPoints (MeasureTiny *measure, off_t Nmeasure, off_t Ncat);
+void ResetAverageActivePhotcodes (SecFilt *secfilt);
+
+

trunk/Ohana/src/relphot/src/BrightCatalog.c

@@ -12 +12 @@
   myAssert (!strcmp(type, #TYPE), "wrong column type");
 
-// XXX double check the header extname values for each table
-// XXX make sure we free things as we can
-// XXX make sure we close files as we can
-// XXX handle free and close on error return as well
 BrightCatalog *BrightCatalogLoad(char *filename) {
 
@@ -622 +618 @@
 int BrightCatalogSplitFree (CatalogSplitter *catalogs) {
 
-  // XXX don't free the catalogs : 
-  // for (i = 0; i < NCATALOG; i++) {
-  //   free (catalogs->catalog[i].averageT);
-  //   free (catalogs->catalog[i].measureT);
-  //   free (catalogs->catalog[i].secfilt);
-  // }
-  // free (catalogs->catalog);
-
+  // don't free the catalogs : free elsewhere
   free (catalogs->catIDs);
   free (catalogs->NAVERAGE);
@@ -652 +641 @@
     catIDmax = MAX(catIDmax, bcatalog->average[i].catID);
   }
-  // XXX validate the measure ID range here
     
   int maxIDold = catalogs->maxID;
@@ -746 +734 @@
 
   // assign the measures to the corresponding catalog
-  // XXX what about averef and related links?  Do I need them?
   for (i = 0; i < bcatalog->Nmeasure; i++) {
     int ID = bcatalog->measure[i].catID;

trunk/Ohana/src/relphot/src/ConfigInit.c

@@ -26 +26 @@
   GetConfig (config, "SIGMA_LIM",              "%lf", 0, &SIGMA_LIM);
 
+  DefConfig ("STAR_SCATTER",                    "%lf",  0.1, STAR_SCATTER);
+  DefConfig ("STAR_CHISQ",                      "%lf", 10.0, STAR_CHISQ);
 
-  GetConfig (config, "STAR_SCATTER",           "%lf", 0, &STAR_SCATTER);
+  DefConfig ("NIGHT_SCATTER",                   "%lf",  0.1, NIGHT_SCATTER);
+  DefConfig ("NIGHT_CHISQ",                     "%lf", 10.0, NIGHT_CHISQ);
+
+  DefConfig ("MOSAIC_SCATTER",                  "%lf",  0.1, MOSAIC_SCATTER);
+  DefConfig ("MOSAIC_CHISQ",                    "%lf", 10.0, MOSAIC_CHISQ);
 
   GetConfig (config, "IMAGE_SCATTER",          "%lf", 0, &IMAGE_SCATTER);
   GetConfig (config, "IMAGE_OFFSET",           "%lf", 0, &IMAGE_OFFSET);
 
-  GetConfig (config, "STAR_CHISQ",             "%lf", 0, &STAR_CHISQ);
   GetConfig (config, "GRID_TOOFEW",            "%d",  0, &GRID_TOOFEW);
   GetConfig (config, "STAR_TOOFEW",            "%d",  0, &STAR_TOOFEW);
@@ -54 +59 @@
 
   DefConfig ("RELPHOT_IMFIT_SYS_SIGMA_LIM", "%lf", 0.01, IMFIT_SYS_SIGMA_LIM);
-  DefConfig ("NIGHT_SCATTER",               "%lf", 0.05, NIGHT_SCATTER);
-  DefConfig ("NIGHT_OFFSET",                "%lf", 0.50, NIGHT_OFFSET);
   DefConfig ("SKY_DEPTH",                   "%d",     2, SKY_DEPTH);
 
@@ -63 +66 @@
 
   GetConfig (config, "ZERO_PT",                "%lf", 0, &ZERO_POINT);
-
-  GetConfig (config, "RELPHOT_GRID_X",         "%d",  0, &RELPHOT_GRID_X);
-  GetConfig (config, "RELPHOT_GRID_Y",         "%d",  0, &RELPHOT_GRID_Y);
-  GetConfig (config, "RELPHOT_GRID_BINNING",   "%d",  0, &RELPHOT_GRID_BINNING);
   GetConfig (config, "CAMERA_CONFIG",          "%s",  0, CameraConfig);
   GetConfig (config, "MOSAICNAME",             "%s",  0, MOSAICNAME);

trunk/Ohana/src/relphot/src/GridOps.c

@@ -1 +1 @@
 # include "relphot.h"
 
-enum {
-    GRID_FITTED,
-    GRID_FROZEN,
-    GRID_REFERENCE,
-};
-
-static int     Ngrid;  // number of grid elements (gridX * gridY)
-static float   *gridM; // magnitude offset for this grid cell 
-static float   *gridS; // stdev of the magnitude offset for this grid cell 
-static int     *gridN; // number of stars used to measure the magnitude offset for this grid cell 
-static int     *gridV; // data mode for this cell: fitted, frozen, reference
-static int      gridX; // number of grid elements in X direction
-static int      gridY; // number of grid elements in Y direction
-
-static int    **bin;   // link from catalog, measure to grid element
-static int    **Xmeas; // grid x-coordinate for a measurement
-static int    **Ymeas; // grid y-coordinate for a measurement
-
-static int    **clist; // link from measurement on a cell to catalog containing measurement
-static off_t  **mlist; // link from measurement on a cell to measurement in a catalog
-static off_t   *Nlist; // number of measurements for each grid cell
-static off_t   *NLIST; // allocated number of measurements for each grid cell
-
-static struct {
-  int Nchip;
-  int Mx, My;  /* mosaic size in chips */
-  int Nx, Ny;  /* chip size in pixels */
-  int *Fx, *Fy;  /* chip flip */
-  int *Ox, *Oy;  /* chip offset */
-  int *valid;
-  char **ccdname;
-} camera;    
-
-static int *ccdnum;
-static char *config;
-
-void initGrid (int dX, int dY) {
-  OHANA_UNUSED_PARAM(dX);
-  OHANA_UNUSED_PARAM(dY);
-
-  int i, N, ccdnum_max, refX, refY, refBin;
-  char *p, field[64], line[256];
-  int *Fx, *Fy;  /* chip flip */
-  int *Ox, *Oy;  /* chip offset */
-  char **ccdname;
-
-  /* load camera config file */
-  config = LoadConfigFile (CameraConfig);
-  if (config == (char *) NULL) {
-    fprintf (stderr, "ERROR: can't find camera config file %s\n", CameraConfig);
-    exit (1);
-  }
-
-  /* load basic mosaic parameters */
-  ScanConfig (config, "NCCD", "%d", 1, &camera.Nchip);
-  ScanConfig (config, "MOSAIC_X", "%d", 1, &camera.Mx);
-  ScanConfig (config, "MOSAIC_Y", "%d", 1, &camera.My);
-  ScanConfig (config, "NAXIS1", "%d", 1, &camera.Nx);
-  ScanConfig (config, "NAXIS2", "%d", 1, &camera.Ny);
-
-  ScanConfig (config, "REFCELL.X", "%d", 1, &refX);
-  ScanConfig (config, "REFCELL.Y", "%d", 1, &refY);
-
-  // temporary storage
-  ALLOCATE (Fx,      int,    camera.Nchip);
-  ALLOCATE (Fy,      int,    camera.Nchip);
-  ALLOCATE (Ox,      int,    camera.Nchip);
-  ALLOCATE (Oy,      int,    camera.Nchip);
-  ALLOCATE (ccdname, char *, camera.Nchip);
-  ALLOCATE (ccdnum,  int,    camera.Nchip);
-
-  /* load per-chip parameters */
-  for (i = 0; i < camera.Nchip; i++) {
-    ALLOCATE (ccdname[i], char, 256);
-    sprintf (field, "CCD.%d", i);
-    ScanConfig (config, field, "%s", 1, line);
-    // XXX get error status!
-
-    sscanf (line, "%s %d %d %d %d", ccdname[i], &Ox[i], &Oy[i], &Fx[i], &Fy[i]);
+/* 
+   We define a 'grid correction', essentially the flat-field correction, as a correction per 
+   photcode subdivided into an NxN array.  The dimensions of the chips corresponding to a
+   photcode will need to be added to the photcode table.  This means a schema update,
+   which I detest.  For now (2021.05.16), I will hard-wire the GPC1 / GPC2 chip size and
+   worry about HSC & Megacam in the future.
+
+   We have a collection of photcodes, with photcodeID limited by design to 64k (unsigned short).
+   
+   Thus we can generate an array of pointers to the grid correction structures and access them 
+   by photcode.
+
+
+ */
+
+static GridCorrectionType **GridCorr = NULL;
+static int                 NGridCorr = 0;
+
+# if (0)
+/* PS1 / PS2 values */
+# define NX_CHIP 4900
+# define NY_CHIP 4900
+# define NX_BIN 16
+# define NY_BIN 16
+# endif
+
+/* test values */
+# define NX_CHIP_DEFAULT 1000
+# define NY_CHIP_DEFAULT 1000
+# define NX_BIN_DEFAULT 2
+# define NY_BIN_DEFAULT 2
+
+void initGridBins (void) {
+
+  // allocate the full possible range of GridCorrectionType pointers, Nphotcode
+  // loop over all images to find actual existing photcodes
+  // generate initial grid values for each existing photcode
+
+  if (!GRID_ZEROPT) return; // skip if we are ignoring the grid correction
+
+  if (GridCorr) return;
+
+  PhotCodeData *photcodes = GetPhotcodeTable();
+  if (!photcodes) return;
+
+  // we have photcodes->Ncodes actually loaded.  loop over them and allocate
+  // array only as large as the max photcodes->code[i].code
+
+  int maxCode = 0;
+  for (int i = 0; i < photcodes->Ncode; i++) {
+    maxCode = MAX(maxCode, photcodes->code[i].code);
+    myAssert (maxCode < 0x10000, "oops");
+  }
+
+  NGridCorr = maxCode + 1;
+  ALLOCATE (GridCorr, GridCorrectionType *, NGridCorr);
+  for (int i = 0; i < NGridCorr; i++) {
+    GridCorr[i] = NULL;
+  }
+
+  fprintf (stderr, "Generating grid corrections for %d photcodes\n", NGridCorr);
+  return;
+}
+
+void freeGridBins() {
+
+  if (!GridCorr) return;
+
+  for (int code = 0; code < NGridCorr; code++) {
+    if (!GridCorr[code]) continue;
+
+    for (int ix = 0; ix < GridCorr[code]->Nx; ix++) {
+      FREE (GridCorr[code]-> Mgrid[ix]);
+      if (GridCorr[code]->dMgrid) FREE (GridCorr[code]->dMgrid[ix]);
+      if (GridCorr[code]->nMgrid) FREE (GridCorr[code]->nMgrid[ix]);
+    }      
+
+    FREE (GridCorr[code]-> Mgrid);
+    FREE (GridCorr[code]->dMgrid);
+    FREE (GridCorr[code]->nMgrid);
+
+    FREE(GridCorr[code]);
+  }
+
+  FREE (GridCorr);
+  GridCorr = NULL;
+
+  return;
+}
+
+GridCorrectionType *getGridCorrNext (int *Nlast) {
+
+  if (GridCorr == NULL) return NULL;
+
+  if (*Nlast >= NGridCorr) return NULL;
+
+  if (*Nlast < 0) *Nlast = -1;
+
+  GridCorrectionType *result = NULL;
+  for (int i = *Nlast + 1; i < NGridCorr; i++) {
+    if (GridCorr[i] == NULL) continue;
+    result = GridCorr[i];
+    *Nlast = i;
+    break;
+  }
+  return result;
+}
+
+GridCorrectionType *getGridCorrByCode (int code) {
+
+  if (GridCorr == NULL) return NULL;
+
+  if (code >= NGridCorr) return NULL;
+  if (code <          0) return NULL;
+
+  ALLOCATE(GridCorr[code], GridCorrectionType, 1)
+  GridCorrectionType *result = GridCorr[code];
+  return result;
+}
+
+/* for GPC1, we have 60 chips, 5 filters, 16x16 grid cells = 2MB of memory for this stuff
+   if we go to 64x64 grid cells (~75 pixels), then it is still only 29MB */
+void initGrid (void) {
+
+  if (!GRID_ZEROPT) return;
+
+  off_t Nimages = 0;
+  Image *images = getimages (&Nimages, NULL);
+
+  int NGridReal = 0;
+
+  for (int i = 0; i < Nimages; i++) {
+    int code = images[i].photcode;
+    myAssert (code >= 0, "oops");
+    myAssert (code < NGridCorr, "oops");
     
-    p = ccdname[i];
-    while (!isdigit(*p) && *p) p++;
-    if (*p == 0) continue;
-    ccdnum[i] = atoi (p);
-  }
-
-  /* we now have the parameters loaded into a minimal length list; reshuffle to the a list of length 0 - MAX(ccdnum) */
-  ccdnum_max = 0;
-  for (i = 0; i < camera.Nchip; i++) {
-    ccdnum_max = MAX(ccdnum_max, ccdnum[i]);
-  }
-  ccdnum_max ++;
-
-  if (ccdnum_max < camera.Nchip) {
-    fprintf (stderr, "problem with camera config: duplicate ccd IDs\n");
-    exit (1);
-  }
-  
-  if (ccdnum_max > 0x1000) {
-    fprintf (stderr, "problem with camera config: absurd max ccd ID number %d\n", ccdnum_max);
-    exit (1);
-  }
-
-  ALLOCATE (camera.Fx,      int,    ccdnum_max);
-  ALLOCATE (camera.Fy,      int,    ccdnum_max);
-  ALLOCATE (camera.Ox,      int,    ccdnum_max);
-  ALLOCATE (camera.Oy,      int,    ccdnum_max);
-  ALLOCATE (camera.valid,   int,    ccdnum_max);
-  ALLOCATE (camera.ccdname, char *, ccdnum_max);
-
-  for (i = 0; i < ccdnum_max; i++) {
-    camera.valid[i] = FALSE;
-  }
-
-  for (i = 0; i < camera.Nchip; i++) {
-    N = ccdnum[i];
-    camera.Fx[N] = Fx[i];
-    camera.Fy[N] = Fy[i];
-    camera.Ox[N] = Ox[i];
-    camera.Oy[N] = Oy[i];
-    camera.ccdname[N] = ccdname[i];
-    camera.valid[N] = TRUE;
-  }
-
-  /* define mosaic 2d correction grid: 
-   * GRID_X is the number of grid pixels per chip in the X direction */
-  gridX = RELPHOT_GRID_X * camera.Mx;
-  gridY = RELPHOT_GRID_Y * camera.My;
-  Ngrid = gridX * gridY;
-
-  ALLOCATE (gridM, float, Ngrid);
-  ALLOCATE (gridS, float, Ngrid);
-  ALLOCATE (gridN, int,   Ngrid);
-  ALLOCATE (gridV, int,   Ngrid);
-
-  // the grid bins may have one of three possible states: fitted, frozen, reference
-  // set the initial values to indicate that the bins are frozen 
-  for (i = 0; i < Ngrid; i++) {
-    gridM[i] = 0.0;
-    gridS[i] = 0.0;
-    gridN[i] = 0;
-    gridV[i] = GRID_FROZEN;
-  }
-
-  // refBin is the index of the grid cell which is kept at 0.0 (all others are relative to this)
-  refBin = refX + refY*gridX;
-  gridV[refBin] = GRID_REFERENCE;
-}
-
-void initGridBins (Catalog *catalog, int Ncatalog) {
-
-  int i, j;
-
-  if (!USE_GRID) return;
-
-  /* define cat,meas -> grid pointers */
-  ALLOCATE (bin, int *, Ncatalog);
-  ALLOCATE (Xmeas, int *, Ncatalog);
-  ALLOCATE (Ymeas, int *, Ncatalog);
-  for (i = 0; i < Ncatalog; i++) {
-    ALLOCATE (bin[i],   int, MAX (catalog[i].Nmeasure, 1));
-    ALLOCATE (Xmeas[i], int, MAX (catalog[i].Nmeasure, 1));
-    ALLOCATE (Ymeas[i], int, MAX (catalog[i].Nmeasure, 1));
-    for (j = 0; j < catalog[i].Nmeasure; j++) bin[i][j] = -1;
-  }
-
-  /* define grid -> cat,meas pointers */
-  ALLOCATE (Nlist, off_t,   Ngrid);
-  ALLOCATE (NLIST, off_t,   Ngrid);
-  ALLOCATE (clist, int *,   Ngrid);
-  ALLOCATE (mlist, off_t *, Ngrid);
-
-  for (i = 0; i < Ngrid; i++) {
-    Nlist[i] = 0;
-    NLIST[i] = 100;
-    ALLOCATE (clist[i], int,   NLIST[i]);
-    ALLOCATE (mlist[i], off_t, NLIST[i]);
-  }
-}
-
-void freeGridBins (int Ncatalog) {
-
-  int i;
-
-  if (!USE_GRID) return;
-
-  /* define cat,meas -> grid pointers */
-  for (i = 0; i < Ncatalog; i++) {
-    free (bin[i]);
-    free (Xmeas[i]);
-    free (Ymeas[i]);
-  }
-  free (bin);
-  free (Xmeas);
-  free (Ymeas);
-
-  /* define grid -> cat,meas pointers */
-  for (i = 0; i < Ngrid; i++) {
-    free (clist[i]);
-    free (mlist[i]);
-  }
-  free (Nlist);
-  free (NLIST);
-  free (clist);
-  free (mlist);
-}
-
-int setGridMeasure (off_t meas, int cat, double X, double Y, int ccdnum) {
-
-  int ix, iy, Cx, Cy, i;
-  double x, y;
-
-  /* X, Y are chip coords on chip ccdnum */
-
-  /* normalize X & Y */
-  x = X;
-  if (camera.Fx[ccdnum]) x = camera.Nx - X;
-  y = Y;
-  if (camera.Fy[ccdnum]) y = camera.Ny - Y;
-
-  /* grid coords on the chip */
-  Cx = MIN (MAX ((x / camera.Nx) * RELPHOT_GRID_X, 0), RELPHOT_GRID_X - 1);
-  Cy = MIN (MAX ((y / camera.Ny) * RELPHOT_GRID_Y, 0), RELPHOT_GRID_Y - 1);
-  
-  /* coordinates in the grid */
-  ix = Cx + camera.Ox[ccdnum]*RELPHOT_GRID_X;
-  iy = Cy + camera.Oy[ccdnum]*RELPHOT_GRID_Y;
-
-  i = ix + iy*gridX;
-
-  bin[cat][meas] = i;
-  Xmeas[cat][meas] = x + camera.Ox[ccdnum]*camera.Nx;
-  Ymeas[cat][meas] = y + camera.Oy[ccdnum]*camera.Ny;
-  clist[i][Nlist[i]] = cat;
-  mlist[i][Nlist[i]] = meas;
-
-  // for the moment, add up the total grid count
-  gridN[i]++;
-
-  Nlist[i] ++;
-  if (Nlist[i] == NLIST[i]) {
-    NLIST[i] += 100;
-    REALLOCATE (clist[i], int,   NLIST[i]);
-    REALLOCATE (mlist[i], off_t, NLIST[i]);
-  }     
-  return (TRUE);
-
-  fprintf (stderr, "error: star out of grid\n");
-  exit (1);
-}
-
-int showGridCount() {
-
-    int ix, iy, i;
-
-    for (iy = 0; iy < gridY; iy++) {
-        for (ix = 0; ix < gridX; ix++) {
-            i = ix + iy*gridX;
-            fprintf (stderr, "%3d  ", gridN[i]);
-        }
-        fprintf (stderr, "\n");
+    // valid photcodes values (code) are in range 1 <= code < 0x10000
+    // photcode == 0 are e.g., PHU (mosaic) images, and should be ignored here
+    if (!code) continue;
+
+    int NX_CHIP = NX_CHIP_DEFAULT;
+    int NY_CHIP = NY_CHIP_DEFAULT;
+    int NX_BIN  = NX_BIN_DEFAULT;
+    int NY_BIN  = NY_BIN_DEFAULT;
+    if (isGPC1chip(code)) { NX_CHIP = 4900; NY_CHIP = 4900; NX_BIN = NY_BIN = GRID_BIN_GPC1; }
+    if (isGPC2chip(code)) { NX_CHIP = 4900; NY_CHIP = 4900; NX_BIN = NY_BIN = GRID_BIN_GPC2; }
+    if (isHSCchip(code))  { NX_CHIP = 2100; NY_CHIP = 4200; NX_BIN = NY_BIN = GRID_BIN_HSC;  }
+    if (isCFHchip(code))  { NX_CHIP = 2100; NY_CHIP = 4200; NX_BIN = NY_BIN = GRID_BIN_CFH;  }
+
+    if (GridCorr[code]) continue; // already created this one
+
+    NGridReal ++;
+
+    ALLOCATE(GridCorr[code], GridCorrectionType, 1); 
+    GridCorr[code]->photcode = code;
+    GridCorr[code]->Nx = NX_BIN;
+    GridCorr[code]->Ny = NY_BIN;
+    GridCorr[code]->dX = NX_BIN / (float) NX_CHIP;
+    GridCorr[code]->dY = NY_BIN / (float) NY_CHIP;
+
+    // we are normally accessing this array randomly, so there is no advantage to
+    // doing Mgrid[y][x] vs Mgrid[x][y]
+    ALLOCATE (GridCorr[code]-> Mgrid, float *, NX_BIN);
+    ALLOCATE (GridCorr[code]->dMgrid, float *, NX_BIN);
+    ALLOCATE (GridCorr[code]->nMgrid,   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);
+    }      
+  }
+  resetMgrid(); // start with values of 0
+  fprintf (stderr, "Init grid corrections for %d photcodes\n", NGridReal);
+}
+
+// reset the values in the arrays to 0
+void resetMgrid () {
+
+  if (!GRID_ZEROPT) return;
+  if (!GridCorr) return;
+
+  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++) {
+        GridCorr[code]-> Mgrid[ix][iy] = 0.0;
+        GridCorr[code]->dMgrid[ix][iy] = 0.0;
+        GridCorr[code]->nMgrid[ix][iy] =   0;
+      }
+    }      
+  }
+}
+
+void setMgrid (Catalog *catalog, int Ncatalog) {
+
+  // check if we are actually doing this step
+  if (!GRID_ZEROPT) return;
+  if (GRID_ZPT_MODE == GRID_ZPT_MODE_NONE) return;
+
+  resetMgrid(); // start with values of 0
+
+  // loop over all measurements, accumulate Sum (in Mgrid), Sum2 (in dMgrid), and Npts
+
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+
+  for (int nc = 0; nc < Ncatalog; nc++) {
+    for (int na = 0; na < catalog[nc].Naverage; na++) {
+
+      int nm = catalog[nc].averageT[na].measureOffset;
+      for (int k = 0; k < catalog[nc].averageT[na].Nmeasure; k++, nm++) {
+
+        // skip measurements marked by AREA or TIME
+        if (catalog[nc].measureT[nm].dbFlags & MEAS_BAD) continue;
+
+        float Mcal = getMcal  (nm, nc, MAG_CLASS_PSF);
+        if (isnan(Mcal)) continue;
+
+        float Mgrp = getMgrp  (nm, nc, catalog[nc].measureT[nm].airmass, NULL);
+        if (isnan(Mgrp)) continue;
+
+        float Mmos  = getMmos  (nm, nc);
+        if (isnan(Mmos)) continue;
+
+        float Mflat = getMflat (nm, nc, catalog);
+        if (isnan(Mflat)) continue;
+
+        // Mrel* is the average magnitude for this star.  For PS1 stacks, we have too much
+        // PSF variability.  We need to calibrate the PSF magnitudes separately from the
+        // Aperture-like magnitues.  (We have an option to use the kron magnitudes or the
+        // other apertures here).  I basically need to do this analysis separately for each
+        // magnitude type
+    
+        float MrelPSF = getMrel  (catalog, nm, nc, MAG_CLASS_PSF, MAG_SRC_CHP);
+        if (isnan(MrelPSF)) continue;
+      
+        float MsysPSF = PhotSysTiny (&catalog[nc].measureT[nm], &catalog[nc].averageT[na], &catalog[nc].secfilt[na*Nsecfilt], MAG_CLASS_PSF);
+        if (isnan(MsysPSF)) continue;
+
+        float Moff =  Mcal + Mgrp + Mmos + Mflat;
+
+        // Msys = Mrel + Moff + Mgrid
+        // thus Mgrid = Msys - Mrel - Moff
+
+        int code = catalog[nc].measureT[nm].photcode;
+        if (code <= 0) continue;
+        if (code >= NGridCorr) continue; // does not match one of our image, skip
+        
+        GridCorrectionType *grid = GridCorr[code];
+        if (!grid) continue; // does not match one of our images, skip
+
+        // edge effects could cause some positions to be slightly out of range
+        // probably should trap extreme outliers
+        int ix = MIN(MAX(0, (int)(catalog[nc].measureT[nm].Xccd * grid->dX)), grid->Nx - 1);
+        int iy = MIN(MAX(0, (int)(catalog[nc].measureT[nm].Yccd * grid->dY)), grid->Ny - 1);
+
+        float dM = MsysPSF - MrelPSF - Moff;
+        
+        // XXX by the time we get here, we should have already mostly fixed up the zero
+        // points.  reject measurements which are way off
+        if (fabs(dM) > 0.5) continue;
+
+        grid-> Mgrid[ix][iy] += dM;
+        grid->dMgrid[ix][iy] += dM*dM;
+        grid->nMgrid[ix][iy] ++;
+      }
     }
-    return (TRUE);
-}
-
-float getMgrid (off_t meas, int cat) {
-
-  int i;
-  float value;
-
-  if (!USE_GRID) return (0);
-  i = bin[cat][meas];
-  if (i == -1) return (NAN);
-
-  // during the grid annealing process, we skip over grid cells until they have enough
-  // valid stars to be fitted
-  if (gridV[i] == GRID_FROZEN) {
-    return (NAN);
-  }
-
-  value = gridM[i];
-  return (value);
-}
-
-/* direct (non-iterative) solution for Mgrid values for all grid bins */
-void setMgridDirect (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) {
-
-  int **gotstar, **gridmeas;
-  int i, j, k, Ngood, Nbad, Nmos, Ncal, Nrel, Nsys, Ngrp;
-  double **A, **B, *Mjx, *Wjx;
-  float Msys, Mcal, Mmos, Mgrp, Merr, Wsys;
-  double Mj, Wj;
-  
-  if (!USE_GRID) return;
-
-  ALLOCATE (A, double *, Ngrid);
-  ALLOCATE (B, double *, Ngrid);
-  for (i = 0; i < Ngrid; i++) {
-    ALLOCATE (A[i], double, Ngrid);
-    ALLOCATE (B[i], double, 1);
-    memset (A[i], 0, Ngrid*sizeof(double));
-    memset (B[i], 0, sizeof(double));
-  }
-
-  Ngood = Nbad = Ncal = Nmos = Nrel = Nsys = Ngrp = 0;
-
-  ALLOCATE (gotstar, int *, Ncatalog);
-  for (i = 0; i < Ncatalog; i++) {
-    ALLOCATE (gotstar[i], int, catalog[i].Naverage);
-  }  
-
-  // set up gridmeas table : grid index for each measurement
-  ALLOCATE (gridmeas, int *, Ncatalog);
-  for (i = 0; i < Ncatalog; i++) {
-    ALLOCATE (gridmeas[i], int, catalog[i].Nmeasure);
-    for (j = 0; j < catalog[i].Nmeasure; j++) {
-      gridmeas[i][j] = -1;
-    }
-  }  
-  for (i = 0; i < Ngrid; i++) {
-    for (j = 0; j < Nlist[i]; j++) {
-      int m, c;
-      m = mlist[i][j];
-      c = clist[i][j];
-      gridmeas[c][m] = i;
-    }
-  }
-
-  // as we loop over the grid cells, we need to accumulate the values of Wjx for each star
-  ALLOCATE (Wjx, double, Ngrid);
-  ALLOCATE (Mjx, double, Ngrid);
-
-  // accumulate the elements of the matrix equation.  We have an equation of the form: Ax = B
-  // where x is the vector of grid cell values G_x (x = 0 - Ngrid), A is an Ngrid x Ngrid matrix,
-  // and B is an Ngrid vector.  For a cell A(x,y), we need the following elements from each
-  // star which touches grid cell (x):
-  // 
-  // Mj  : sum over all measurements of Msys / dMsys^2 
-  // Wj  : sum over all measurements of 1.0 / dMsys^2 
-  // Mjx : sum over all measurements which touch cell (x) of Msys / dMsys^2
-  // Wjx : sum over all measurements which touch cell (x) of 1.0 / dMsys^2
-
-  // Mjx and Wjx can be calculated by summing over all measurements which touch the cell 
-  // Mj requires looping over stars which touch (x)
-
-  // this is tricky because we need to know both the measurements which touch a cell
-  // and the stars for which any measurement touches a cell.  we need to accumulate
-  // sums for each star which touches as cell on both bases.
-
-  int Nsecfilt = GetPhotcodeNsecfilt ();
-  int thisCode = photcodes[0][0].code;
-  int Nsec = GetPhotcodeNsec(thisCode);
-
-  for (i = 0; i < Ngrid; i++) {
-    
-    for (j = 0; j < Ncatalog; j++) {
-      memset (gotstar[j], 0, catalog[j].Naverage*sizeof(int));
-    }  
-
-    // we are looping over the stars, but doing so by looping over the set of measurements:
-    // every star which touches this grid cell has a measurement in Nlist[i]
-    for (j = 0; j < Nlist[i]; j++) {
-      
-      int mx, c, n, m0, Npts;
-
-      mx = mlist[i][j];
-      c  = clist[i][j];
-      n  = catalog[c].measureT[mx].averef;
-      
-      // if we have already visited this star, skip the stuff below
-      if (gotstar[c][n]) continue;
-      gotstar[c][n] = TRUE;
-
-      // skip stars marked as BAD
-      if (catalog[c].secfilt[n*Nsecfilt+Nsec].flags & STAR_BAD) {
-        Nrel ++;
-        continue;
-      }
-
-      m0 = catalog[c].average[n].measureOffset;
-
-      // we accumuate an entry for each cell
-      memset (Wjx, 0, Ngrid*sizeof(double));
-      memset (Mjx, 0, Ngrid*sizeof(double));
-      Npts = Mj = Wj = 0.0;
-
-      // if we have not yet visited this star, accumulate the Mj, Wj entries
-      for (k = 0; k < catalog[c].average[n].Nmeasure; k++) {
-
-        int m, Ng;
-
-        m = m0 + k;
-
-        // skip measurements marked as BAD
-        if (catalog[c].measureT[m].dbFlags & MEAS_BAD) {
-          Nbad ++;
+  }
+
+  // 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;
         }
 
-        // skip images marked as BAD
-        Mcal = getMcal  (m, c, MAG_CLASS_PSF);
-        if (isnan(Mcal)) {
-          Ncal ++;
-          continue;
-        }
-
-        // skip mosaics marked as BAD
-        Mmos = getMmos  (m, c);
-        if (isnan(Mmos)) {
-          Nmos ++;
-          continue;
-        }
-
-        // skip mosaics marked as BAD
-        Mgrp = getMgrp  (m, c, catalog[c].measureT[m].airmass, NULL);
-        if (isnan(Mgrp)) {
-          Ngrp ++;
-          continue;
-        }
-
-        // select the color- and airmass-corrected observed magnitude for this star
-        // XXX need to be able to turn off the color-correction until initial average mags are found
-        Msys = PhotCatTiny (&catalog[c].measureT[m], MAG_CLASS_PSF);
-        if (isnan(Msys)) {
-          Nsys++;
-          continue;
-        }
-
-        // mag-error for this measurement
-        Merr =  MAX (catalog[c].measureT[m].dM, MIN_ERROR);
-
-        // disable Wsys for now
-        Wsys = TRUE ? 1.0 : 1.0 / SQ(Merr);
-
-        Ng = gridmeas[c][m];
-        if (Ng == -1) continue;  // skip measurements which do not touch any cell
-
-        Mj += Msys * Wsys;  // we are only including measurements touching this cell
-        Wj += Wsys;  // we are only including measurements touching this cell
-        Npts ++;
-        Ngood ++;
-
-        Mjx[Ng] += Msys * Wsys;  // we are only including measurements touching cell (x)
-        Wjx[Ng] += Wsys;  // we are only including measurements touching cell (x)
-      }
-
-      // some stars will not have any valid measurements, skip these
-      if (Npts == 0) continue;
-
-      B[i][0] += Mj*Wjx[i]/Wj - Mjx[i];
-      A[i][i] -= Wjx[i];
-      for (k = 0; k < Ngrid; k++) {
-        A[i][k] += Wjx[i]*Wjx[k]/Wj;
-        // fprintf (stderr, "%3.0f ", Wjx[k]);
+        float Mgrid  = GridCorr[code]-> Mgrid[ix][iy] / GridCorr[code]->nMgrid[ix][iy]; // average Mgrid
+        float Mgrid2 = GridCorr[code]-> dMgrid[ix][iy] / GridCorr[code]->nMgrid[ix][iy]; // average Mgrid^2
+
+        float r = GridCorr[code]->nMgrid[ix][iy] / (float) (GridCorr[code]->nMgrid[ix][iy] - 1.0); // pop -> sample stdev
+
+        GridCorr[code]-> Mgrid[ix][iy] = Mgrid;
+        GridCorr[code]->dMgrid[ix][iy] = sqrt(r*(Mgrid2 - Mgrid*Mgrid)); // sample stdev
+        // fprintf (stderr, "grid code %d, %d x %d : %f +/- %f : %d\n", code, ix, iy, GridCorr[code]-> Mgrid[ix][iy], GridCorr[code]->dMgrid[ix][iy], GridCorr[code]->nMgrid[ix][iy]);
       }
     }
   }
-
-  if (1) {
-
-    FILE *f;
-    Header theader;
-    Matrix matrix;
-
-    /* we are writing to this file */
-    f = fopen ("matrix.fits", "w");
-    if (f == (FILE *) NULL) { 
-      fprintf (stderr, "cannot open matrix.fits for output\n");
-      return;
-    }
-
-    /* save grid mag values */
-    gfits_init_header (&theader);
-    theader.Naxes = 2;
-    theader.Naxis[0] = Ngrid;
-    theader.Naxis[1] = Ngrid;
-    theader.bitpix   = -32;
-    gfits_create_Theader (&theader, "IMAGE");
-    gfits_modify (&theader, "EXTNAME", "%s", 1, "MATRIX");
-    gfits_create_matrix  (&theader, &matrix);
-    for (i = 0; i < Ngrid; i++) {
-      for (j = 0; j < Ngrid; j++) {
-        gfits_set_matrix_value (&matrix, i, j, (double) A[i][j]);
-      }
-    }
-    gfits_fwrite_header (f, &theader);
-    gfits_fwrite_matrix (f, &matrix);
-    gfits_free_matrix (&matrix);
-
-    gfits_modify (&theader, "EXTNAME", "%s", 1, "TRPOSE");
-    gfits_create_matrix  (&theader, &matrix);
-    for (i = 0; i < Ngrid; i++) {
-      for (j = 0; j < Ngrid; j++) {
-        gfits_set_matrix_value (&matrix, i, j, (double) A[j][i]);
-      }
-    }
-    gfits_fwrite_header (f, &theader);
-    gfits_fwrite_matrix (f, &matrix);
-    gfits_free_matrix (&matrix);
-    fclose (f);
-  }
-
-  dgaussjordan (A, B, Ngrid, 1);
-
-  fprintf (stderr, "grid cells fitted (Ngood: %d, Nbad: %d, Nmos: %d, Ncal: %d, Nrel: %d, Nsys: %d)\n", Ngood, Nbad, Nmos, Ncal, Nrel, Nsys);
-
-  for (i = 0; i < Ngrid; i++) {
-    gridM[i] = B[i][0];
-    gridS[i] = sqrt(A[i][i]);
-    gridN[i] = Ngood;
-  }
-
-  free (Wjx);
-  free (Mjx);
-
-  for (i = 0; i < Ngrid; i++) {
-    free (A[i]);
-  }
-  free (A);
-  free (B);
-
-  for (i = 0; i < Ncatalog; i++) {
-    free (gotstar[i]);
-    free (gridmeas[i]);
-  }  
-  free (gotstar);
-  free (gridmeas);
-}
-
-/* determine Mgrid values for all grid bins */
-void setMgrid (Catalog *catalog, FlatCorrectionTable *flatcorr) {
-
-  int i, j, m, c, n, N, Nmax, Nbad, Nmos, Ngrp, Ncal, Nrel, Nsys, Nfit;
-  double *list, *dlist;
-  float Msys, Mrel, Mcal, Mmos, Mgrp;
-
-  StatType stats;
-  liststats_setmode (&stats, "INNER_WTMEAN");
+  return;
+}
+
+float getMgrid (Measure *measure) {
+
+  if (!GRID_ZEROPT) return 0.0;
+  if (GRID_ZPT_MODE == GRID_ZPT_MODE_NONE) return 0.0;
+
+  int code = measure->photcode;
+  if (code <= 0) return 0.0;
+  if (code >= NGridCorr) return 0.0; // does not match one of our image, skip
+        
+  GridCorrectionType *grid = GridCorr[code];
+  if (!grid) return 0.0; // does not match one of our images, skip
   
-  if (!USE_GRID) return;
-
-  int Nsecfilt = GetPhotcodeNsecfilt ();
-
-  Nmax = Nlist[0];
-  for (i = 0; i < Ngrid; i++) {
-    Nmax = MAX (Nmax, Nlist[i]);
-  }
-  ALLOCATE (list, double, Nmax);
-  ALLOCATE (dlist, double, Nmax);
-
-  Nbad = Ncal = Nmos = Ngrp = Nrel = Nsys = Nfit = 0;
-
-  for (i = 0; i < Ngrid; i++) {
-
-    N = 0;
-    for (j = 0; j < Nlist[i]; j++) {
-      
-      m = mlist[i][j];
-      c = clist[i][j];
-      
-      if (catalog[c].measureT[m].dbFlags & MEAS_BAD) {
-        Nbad ++;
-        continue;
-      }
-      Mcal = getMcal  (m, c, MAG_CLASS_PSF);
-      if (isnan(Mcal)) {
-        Ncal ++;
-        continue;
-      }
-      Mmos = getMmos  (m, c);
-      if (isnan(Mmos)) {
-        Nmos ++;
-        continue;
-      }
-      Mgrp = getMgrp  (m, c, catalog[c].measureT[m].airmass, NULL);
-      if (isnan(Mgrp)) {
-        Ngrp ++;
-        continue;
-      }
-      Mrel  = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP);
-      if (isnan(Mrel)) {
-        Nrel ++;
-        continue;
-      }
-      
-      n = catalog[c].measureT[m].averef;
-      Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF);
-      if (isnan(Msys)) {
-        Nsys++;
-        continue;
-      }
-      list[N] = Msys - Mrel - Mcal - Mmos - Mgrp;
-      dlist[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
-      N++;
-    }
-
-    // the reference Cell is forced to have a value of 0.0, and is never changed to GRID_FITTED
-    if (gridV[i] == GRID_REFERENCE) {
-      gridM[i] = 0.0;
-      gridS[i] = 0.0;
-      gridN[i] = N;
-      continue;
-    }
-
-    // until we have enough valid measurements on this grid cell, skip it
-    if (N < GRID_TOOFEW) {
-      gridV[i] = GRID_FROZEN;
-      continue;
-    }
-
-    liststats (list, dlist, NULL, N, &stats);
-    gridM[i] = stats.mean;
-    gridS[i] = stats.sigma;
-    gridN[i] = N;
-    gridV[i] = GRID_FITTED;
-    Nfit++;
-  }
-
-  fprintf (stderr, "%d of %d grid cells fitted (+ reference cell) (Nbad: %d, Nmos: %d, Ncal: %d, Nrel: %d, Nsys: %d)\n", Nfit, Ngrid, Nbad, Nmos, Ncal, Nrel, Nsys);
-
-  free (list);
-  free (dlist);
-}
-
-void plot_grid (Catalog *catalog, FlatCorrectionTable *flatcorr) {
-
-  int i, j, m, c, n, N, Narea;
-  float Msys, Mrel, Mcal, Mmos, Mgrp;
-  double *xlist, *Mlist, *dlist, *ylist;
-  Graphdata graphdata;
-
-  if (!USE_GRID) return;
-
-  int Nsecfilt = GetPhotcodeNsecfilt ();
-
-  N = 0;
-  for (i = 0; i < Ngrid; i++) 
-    N += Nlist[i];
-
-  ALLOCATE (xlist, double, N);
-  ALLOCATE (ylist, double, N);
-  ALLOCATE (Mlist, double, N);
-  ALLOCATE (dlist, double, N);
-
-  Narea = 0;
-  N = 0;
-  for (i = 0; i < Ngrid; i++) {
-    for (j = 0; j < Nlist[i]; j++) {
-      
-      m = mlist[i][j];
-      c = clist[i][j];
-      
-      if (catalog[c].measureT[m].dbFlags & MEAS_BAD) {
-        Narea ++;
-        continue;
-      }
-      Mcal  = getMcal  (m, c, MAG_CLASS_PSF);
-      if (isnan(Mcal)) continue;
-      Mmos  = getMmos  (m, c);
-      if (isnan(Mmos)) continue;
-      Mgrp  = getMgrp  (m, c, catalog[c].measureT[m].airmass, NULL);
-      if (isnan(Mgrp)) continue;
-      Mrel  = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP);
-      if (isnan(Mrel)) continue;
-
-      n = catalog[c].measureT[m].averef;
-      Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF);
-
-      xlist[N] = Xmeas[c][m];
-      ylist[N] = Ymeas[c][m];
-      Mlist[N] = Msys - Mrel - Mcal - Mmos - Mgrp;
-      dlist[N] = Msys - Mrel - Mcal - Mmos - Mgrp - gridM[i];
-      N++;
-    }
-  }
-
-  fprintf (stderr, "skipped %d meas for area\n", Narea);
-
-  plot_defaults (&graphdata);
-  graphdata.ymin = PlotdMmin;
-  graphdata.ymax = PlotdMmax;
-  plot_list (&graphdata, xlist, Mlist, N, "X vs dM raw", "%s.XdM.png", OUTROOT);
-  plot_list (&graphdata, xlist, dlist, N, "X vs dM corrected", "%s.XdMf.png", OUTROOT);
-  plot_list (&graphdata, ylist, dlist, N, "Y vs dM corrected", "%s.YdMf.png", OUTROOT);
-
-  plot_defaults (&graphdata);
-  plot_list (&graphdata, xlist, ylist, N, "X vs Y", "%s.XY.png", OUTROOT);
-
-  free (ylist);
-  free (xlist);
-  free (Mlist);
-  free (dlist);
-
-}
-
-void dump_grid () { 
-
-  off_t i, Nimage;
-  int j, Nbytes, Nformat;
-  FILE *f;
-  Header header, theader;
-  Matrix matrix;
-  Mosaic *refmosaic;
-  char *filename;
-  char formatline[32], key[32], value[64];
-    
-  Nbytes = strlen (OUTROOT) + 6;
-  ALLOCATE (filename, char, Nbytes);
-  snprintf (filename, Nbytes, "%s.fits", OUTROOT);
-
-  /* select reference mosaic image */
-  // off_t *imlist = SelectRefMosaic (&refmosaic, &Nimage); return value ignored
-  SelectRefMosaic (&refmosaic, &Nimage);
-
-  /* we are writing to this file */
-  f = fopen (filename, "w");
-  if (f == (FILE *) NULL) { 
-    fprintf (stderr, "cannot open %s for output\n", filename);
-    free (filename);
-    return;
-  }
-
-  /* create empty phu */
-  gfits_init_header (&header);
-  header.extend = TRUE;
-  gfits_create_header (&header);
-  gfits_create_matrix (&header, &matrix);
-  gfits_modify (&header, "NEXTEND", OFF_T_FMT, 1,  Nimage + 3);
-  gfits_modify (&header, "FILTER", "%s", 1, photcodes[0][0].name);  // XXXX note that this expects a single photcode, enforced in initialize.d
-  gfits_modify_alt (&header, "COMMENT", "%S", 1, "Mosaic Photometry Grid Analysis");
-
-  // we need to add lines to the PHU to identify the camera and format; these are used by the ipp config system
-  // Note that config must have been loaded (and not freed) above.
-  ScanConfig (config, "NFORMAT", "%d", 1, &Nformat);
-  for (i = 1; i <= Nformat; i++) {
-      ScanConfig (config, "FORMAT", "%s", i, formatline);
-      sscanf (formatline, "%s %s", key, value);
-      gfits_modify (&header, key, "%s", 1, value);
-  }
-
-  gfits_fwrite_header (f, &header);
-  gfits_fwrite_matrix (f, &matrix);
-  gfits_free_matrix (&matrix);
-
-  /* save grid mag values */
-  gfits_init_header (&theader);
-  theader.Naxes = 2;
-  theader.Naxis[0] = gridX;
-  theader.Naxis[1] = gridY;
-  theader.bitpix   = -32;
-  gfits_create_Theader (&theader, "IMAGE");
-  gfits_modify (&theader, "FILTER", "%s", 1, photcodes[0][0].name);
-  gfits_modify (&theader, "EXTNAME", "%s", 1, "MAG_OFFSET");
-  gfits_create_matrix  (&theader, &matrix);
-  for (i = 0; i < gridX; i++) {
-    for (j = 0; j < gridY; j++) {
-      gfits_set_matrix_value (&matrix, i, j, (double) gridM[i + j*gridX]);
-    }
-  }
-  write_coords (&theader, &refmosaic[0].coords);
-  gfits_fwrite_header (f, &theader);
-  gfits_fwrite_matrix (f, &matrix);
-  gfits_free_matrix (&matrix);
-
-  /* save grid Nmeas values */
-  gfits_modify (&theader, "EXTNAME", "%s", 1, "NMEAS");
-  gfits_modify (&theader, "FILTER", "%s", 1, photcodes[0][0].name);
-  gfits_create_matrix  (&theader, &matrix);
-  for (i = 0; i < gridX; i++) {
-    for (j = 0; j < gridY; j++) {
-      gfits_set_matrix_value (&matrix, i, j, (double) gridN[i + j*gridX]);
-    }
-  }
-  write_coords (&theader, &refmosaic[0].coords);
-  gfits_fwrite_header (f, &theader);
-  gfits_fwrite_matrix (f, &matrix);
-  gfits_free_matrix (&matrix);
-
-  /* save grid sigma values */
-  gfits_modify (&theader, "EXTNAME", "%s", 1, "SIGMA");
-  gfits_modify (&theader, "FILTER", "%s", 1, photcodes[0][0].name);
-  gfits_create_matrix  (&theader, &matrix);
-  for (i = 0; i < gridX; i++) {
-    for (j = 0; j < gridY; j++) {
-      gfits_set_matrix_value (&matrix, i, j, (double) gridS[i + j*gridX]);
-    }
-  }
-  write_coords (&theader, &refmosaic[0].coords);
-  gfits_fwrite_header (f, &theader);
-  gfits_fwrite_matrix (f, &matrix);
-  gfits_free_matrix (&matrix);
-
-  /* calculate value for each CCD pixel, write out CCD images */
-  /* grid pixels are tied to detector pixels, but are flipped to match focal plane */
-  for (i = 0; i < camera.Nchip; i++) {
-    int N, ix, iy, x, y, X, Y, bin;
-
-    N = ccdnum[i];
-
-    gfits_modify (&theader, "EXTNAME", "%s", 1, camera.ccdname[N]);
-    gfits_modify (&theader, "FILTER", "%s", 1, photcodes[0][0].name);
-    gfits_modify (&theader, "NX", "%d", 1, camera.Nx);
-    gfits_modify (&theader, "NY", "%d", 1, camera.Ny);
-      
-    theader.Naxis[0] = RELPHOT_GRID_X;
-    theader.Naxis[1] = RELPHOT_GRID_Y;
-    gfits_modify (&theader, "NAXIS1", "%d", 1, RELPHOT_GRID_X);
-    gfits_modify (&theader, "NAXIS2", "%d", 1, RELPHOT_GRID_Y);
-    gfits_create_matrix  (&theader, &matrix);
-
-    for (Y = 0; Y < RELPHOT_GRID_Y; Y++) {
-      for (X = 0; X < RELPHOT_GRID_X; X++) {
-              
-        /* normalize X & Y */
-        x = X;
-        if (camera.Fx[N]) x = RELPHOT_GRID_X - X - 1;
-        y = Y;
-        if (camera.Fy[N]) y = RELPHOT_GRID_Y - Y - 1;
-              
-        /* coordinates in the grid */
-        ix = x + camera.Ox[N]*RELPHOT_GRID_X;
-        iy = y + camera.Oy[N]*RELPHOT_GRID_Y;
-              
-        bin = ix + iy*gridX;
-        gfits_set_matrix_value (&matrix, X, Y, (double) gridM[bin]);
-      }
-    }
-    gfits_fwrite_header (f, &theader);
-    gfits_fwrite_matrix (f, &matrix);
-    gfits_free_matrix (&matrix);
-  }
-
-  free (filename);
-}
-
-void InterpolateGrid (float *buffer, int Nx, int Ny, Coords *ccd, Coords *gcoords) {
-
-  int i, j;
-  double x, y, r, d, X, Y, dx, dy;
-  double V00, V01, V10, V11;
-  double wV00, wV01, wV10, wV11;
-  double dV00, dV01, dV10, dV11;
-  double v1, v2, value;
-  int ix, iy, N;
-
-  for (i = 0; i < Nx; i++) {
-    for (j = 0; j < Ny; j++) {
-      x = i * RELPHOT_GRID_BINNING / 2;
-      y = j * RELPHOT_GRID_BINNING / 2;
-      XY_to_RD (&r, &d, x, y, ccd);
-      RD_to_XY (&X, &Y, r, d, gcoords);
-
-      X = X / RELPHOT_GRID_BINNING;
-      Y = Y / RELPHOT_GRID_BINNING;
-
-      ix = (int) X;
-      dx = X - ix;
-      iy = (int) Y;
-      dy = Y - iy;
-
-      if (ix < 0) continue;
-      if (iy < 0) continue;
-      if (ix >= gridX) continue;
-      if (iy >= gridY) continue;
-
-      N = ix + iy*gridX;
-      V00 = gridM[N];
-      V10 = gridM[N + 1];
-      V01 = gridM[N + gridX];
-      V11 = gridM[N + gridX + 1];
-
-      dV00 = gridS[N];
-      dV10 = gridS[N + 1];
-      dV01 = gridS[N + gridX];
-      dV11 = gridS[N + gridX + 1];
-
-      wV00 = (dV00 == 0) ? 0.0 : 1 / SQ(dV00);
-      wV01 = (dV01 == 0) ? 0.0 : 1 / SQ(dV01);
-      wV10 = (dV10 == 0) ? 0.0 : 1 / SQ(dV10);
-      wV11 = (dV11 == 0) ? 0.0 : 1 / SQ(dV11);
-
-      v1 = wV00*V00*(1 + dx*dy - dx - dy) +
-        wV10*V10*(dx - dx*dy) +
-        wV01*V01*(dy - dx*dy) +
-        wV11*V11*(dx*dy);
-
-      v2 = wV00*(1 + dx*dy - dx - dy) +
-        wV10*(dx - dx*dy) +
-        wV01*(dy - dx*dy) +
-        wV11*(dx*dy);
-
-      value = v1 / v2;
-      buffer[j*Nx + i] = value;
-    }
-  }
-}
+  // edge effects could cause some positions to be slightly out of range
+  // probably should trap extreme outliers
+  int ix = MIN(MAX(0, (int)(measure->Xccd * grid->dX)), grid->Nx - 1);
+  int iy = MIN(MAX(0, (int)(measure->Yccd * grid->dY)), grid->Ny - 1);
+
+  float Mgrid = grid-> Mgrid[ix][iy];
+  return Mgrid;
+}
+
+float getMgridTiny (MeasureTiny *measure) {
+
+  if (!GRID_ZEROPT) return 0.0;
+  if (GRID_ZPT_MODE == GRID_ZPT_MODE_NONE) return 0.0;
+
+  int code = measure->photcode;
+  if (code <= 0) return 0.0;
+  if (code >= NGridCorr) return 0.0; // does not match one of our image, skip
+        
+  GridCorrectionType *grid = GridCorr[code];
+  if (!grid) return 0.0; // does not match one of our images, skip
+  
+  // edge effects could cause some positions to be slightly out of range
+  // probably should trap extreme outliers
+  int ix = MIN(MAX(0, (int)(measure->Xccd * grid->dX)), grid->Nx - 1);
+  int iy = MIN(MAX(0, (int)(measure->Yccd * grid->dY)), grid->Ny - 1);
+
+  float Mgrid = grid-> Mgrid[ix][iy];
+  return Mgrid;
+}
+
+// for historical reasons, Mflat and Mgrid have opposite signs
+void setMflatFromGrid (Catalog *catalog) {
+  if (!GRID_ZEROPT) return;
+  for (off_t j = 0; j < catalog->Nmeasure; j++) {
+    float Mgrid = getMgrid (&catalog->measure[j]);
+    catalog->measure[j].Mflat += Mgrid;
+  }
+}
+

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

@@ -319 +319 @@
 }
 
-int findCCD (off_t idx, off_t meas, int cat, MeasureTiny *measure) {
-
-  int ccdnum;
-  double X, Y;
-  char *pname, *filter, *p, base[256];
-
-  /* identify the ccd on the basis of the photcode name */
-  pname = GetPhotcodeNamebyCode (image[idx].photcode);
-
-  // skip measurements which do not match one of the requested photcodes (
-  // (do we not already exclude in bcatalog -- maybe needed for reload_objects?
-  int Ns = GetActivePhotcodeIndex (measure[0].photcode);
-  if (Ns < 0) return FALSE;
-
-  filter = photcodes[Ns][0].name;
-  sprintf (base, "%s.%s.", MOSAICNAME, filter);
-  if (strncmp (pname, base, strlen (base))) return (FALSE);
-  p = pname + strlen(base);
-
-  /* p is pointing at the DETECTOR ID part of the photcode: CAMERA.FILTER.DETECTOR 
-     for now, we will let this be of the form SSSSnn where SSS is an arbitrary string
-     and nn is an integer -> ccdnum */
-      
-  /// XXX we have an inconsistency here wrt dvo.layout and dvo.photcode
-
-  while (!isdigit(*p) && *p) p++;
-  if (*p == 0) return (FALSE);
-  ccdnum = atoi (p);
-
-  /* ccdnum is an integer, but not necessarily a sequence number.  when the camera layout is
-     constructed, there will be null values for undefined ccdnums */
-
-  // old code to add this measurement to the grid cell for this chip 
-  // ave = measureT[0].averef;
-  // ra  = catalog[cat].averageT[ave].R - measureT[0].dR / 3600.0;
-  // dec = catalog[cat].averageT[ave].D - measureT[0].dD / 3600.0;
-  // RD_to_XY (&X, &Y, ra, dec, &image[i].coords);
-
-  // XXX we can now use these values (but need to be careful about old formats)
-  X = measure[0].Xccd;
-  Y = measure[0].Yccd;
-  setGridMeasure (meas, cat, X, Y, ccdnum);
-
-  return (TRUE);
-}
-
 void matchImage (Catalog *catalog, off_t meas, int cat, int doImageList) {
 
   off_t idx, ID;
-  int status;
   MeasureTiny *measure;
   
@@ -381 +334 @@
   catalog[cat].measureT[meas].myDet = TRUE;
 
-  if (USE_GRID) {
-    status = findCCD (idx, meas, cat, measure);
-    if (!status) {
-      if (VERBOSE2) fprintf (stderr, "failed to determine CCD for "OFF_T_FMT", %d\n",  meas, cat);
-      return;
-    }
-  }
-
   // index for (catalog, measure) -> image
   MeasureToImage[cat][meas] = idx;
@@ -440 +385 @@
   if (i == -1) return NAN;
 
-  if (image[i].flags & IMAGE_BAD) return NAN;  
+  // allow context dependence?: in relphot_images, we want to exclude detections from poor images
+  // but in the final pass, we want to allow even bad detections 
+  if (image[i].flags & (ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW)) return NAN;  
 
   switch (class) {
@@ -464 +411 @@
   if (i == -1) return (1000);
 
-  if (image[i].flags & IMAGE_BAD) return (1000);  
   distance = image[i].ubercalDist; // was dummy3 in structure
   return (distance);
@@ -552 +498 @@
 }
 
-float getMflat (off_t meas, int cat, FlatCorrectionTable *flatcorr, Catalog *catalog) {
-  OHANA_UNUSED_PARAM(flatcorr);
-
-  off_t i = MeasureToImage[cat][meas];
-  if (i == -1) return (NAN);
-
-  float offset = catalog[cat].measureT[meas].Mflat;
-  if (!isfinite(offset)) {
-    offset = 0.0;
-  }
-
-  // to do this, I need to pass in the catalog and flatcorr pointers
-  // int flat_id = image[i].photom_map_id;
-  // if (flat_id) {
-  //   offset = FlatCorrectionOffset (flatcorr, flat_id, catalog[cat].measureT[meas].Xccd, catalog[cat].measureT[meas].Yccd);
-  // }
-
-  return (offset);
-}
-
 Coords *getCoords (off_t meas, int cat) {
 
@@ -582 +508 @@
 
 /* determine Mcal values for all images */
-void setMcal (Catalog *catalog, FlatCorrectionTable *flatcorr) {
+void setMcal (Catalog *catalog) {
 
   off_t i, j, m, c, n;
@@ -690 +616 @@
           continue;
       }
-      float Mgrid = getMgrid (m, c);
+      float Mgrid = getMgridTiny (&catalog[c].measureT[m]);
       if (isnan(Mgrid)) {
           Ngrid++;
@@ -711 +637 @@
       // apply that offset as well here for this image (in other words, each detection is
       // being compared to the model, excluding the zero point, Mcal.  The model includes
-      // the flat-correction.  NOTE the sign of Mflat (Image.Mcal = Measure.Mcal - Mflat)
-
-      float Mflat = getMflat (m, c, flatcorr, catalog);
+      // the flat-correction.  NOTE the sign of Mflat (Image.Mcal = Measure.Mcal + Mflat)
+      // this was inconsistent w.r.t. PhotRel pre-r41606
+
+      float Mflat = getMflat (m, c, catalog);
 
       n = catalog[c].measureT[m].averef;
@@ -722 +649 @@
       }
 
-      float Moff = Mmos + Mgrp + Mgrid - Mflat;
+      float Moff = Mmos + Mgrp + Mgrid + Mflat;
 
       PhotCode *code = GetPhotcodebyCode (catalog[c].measureT[m].photcode);
@@ -763 +690 @@
     if (mark) {
       image[i].flags |= ID_IMAGE_PHOTOM_FEW;
+      image[i].McalPSF    = 0.0;
+      image[i].McalAPER   = 0.0;
+      image[i].McalChiSq  = NAN;
       image[i].dMcal      = NAN;
-      image[i].McalChiSq  = NAN;
+      image[i].dMagSys    = NAN;
       image[i].nFitPhotom = 0;
       Nfew ++;
@@ -835 +765 @@
 
 /* determine McalTEST values for all images -- this is not used to set measure.Mcal, but only to test */
-void setMcalTest (Catalog *catalog, FlatCorrectionTable *flatcorr) {
+void setMcalTest (Catalog *catalog) {
 
   off_t i, j, m, c, n;
@@ -889 +819 @@
           continue;
       }
-      float Mgrid = getMgrid (m, c);
+      float Mgrid = getMgridTiny (&catalog[c].measureT[m]);
       if (isnan(Mgrid)) {
           Ngrid++;
@@ -905 +835 @@
       // apply that offset as well here for this image (in other words, each detection is
       // being compared to the model, excluding the zero point, Mcal.  The model includes
-      // the flat-correction.  NOTE the sign of Mflat (Image.Mcal = Measure.Mcal - Mflat)
-
-      float Mflat = getMflat (m, c, flatcorr, catalog);
+      // the flat-correction.  NOTE the sign of Mflat (Image.Mcal = Measure.Mcal + Mflat)
+
+      float Mflat = getMflat (m, c, catalog);
 
       // get the PSF magnitude for thie measurement, with airmass slope applied
@@ -917 +847 @@
       }
 
-      float Moff = Mmos + Mgrp + Mgrid - Mflat;
+      float Moff = Mmos + Mgrp + Mgrid + Mflat;
 
       psfStars.alldata-> yVector[Nref] = MsysPSF - MrelPSF - Moff;
@@ -968 +898 @@
   // measure stats for Mcal and dMcal
   for (i = N = 0; i < Nimage; i++) {
-    if (image[i].flags & IMAGE_BAD) continue;
+    if (image[i].flags & (ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW)) continue;
 
     if (FREEZE_IMAGES && isMosaicChip(image[i].photcode)) continue;
@@ -996 +926 @@
 
     mark = FALSE;
-    image[i].flags &= ~ID_IMAGE_PHOTOM_POOR;
-    mark = (image[i].dMcal > MaxScatter) || (fabs(image[i].McalPSF - MedOffset) > MaxOffset);
+
+    if (image[i].flags & ID_IMAGE_PHOTOM_FEW) {
+      mark = TRUE;
+    }
+    if (image[i].dMcal > MaxScatter) {
+      mark = TRUE;
+    }
+    if (fabs(image[i].McalPSF - MedOffset) > MaxOffset) {
+      mark = TRUE;
+    }
+
     if (mark) { 
       Nmark ++;
@@ -1146 +1085 @@
 }
 
-StatType statsImageN (Catalog *catalog, FlatCorrectionTable *flatcorr) {
+StatType statsImageN (Catalog *catalog) {
 
   off_t i, j, m, c, n, N;
@@ -1163 +1102 @@
   n = 0;
   for (i = 0; i < Nimage; i++) {
-    if (image[i].flags & IMAGE_BAD)  continue;
+    if (image[i].flags & (ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW)) continue;
 
     if (FREEZE_IMAGES && isMosaicChip(image[i].photcode)) continue;
@@ -1179 +1118 @@
       Mgrp  = getMgrp  (m, c, catalog[c].measureT[m].airmass, NULL);
       if (isnan(Mgrp)) continue;
-      Mgrid = getMgrid (m, c);
+      Mgrid = getMgridTiny (&catalog[c].measureT[m]);
       if (isnan(Mgrid)) continue;
       N++;
@@ -1212 +1151 @@
   for (i = 0; i < Nimage; i++) {
 
-    if (image[i].flags & IMAGE_BAD)  continue;
+    if (image[i].flags & (ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW)) continue;
 
     if (FREEZE_IMAGES && isMosaicChip(image[i].photcode)) continue;
@@ -1245 +1184 @@
   for (i = 0; i < Nimage; i++) {
 
-    if (image[i].flags & IMAGE_BAD)  continue;
+    if (image[i].flags & (ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW)) continue;
 
     if (FREEZE_IMAGES && isMosaicChip(image[i].photcode)) continue;
@@ -1278 +1217 @@
   for (i = 0; i < Nimage; i++) {
 
-    if (image[i].flags & IMAGE_BAD)  continue;
+    if (image[i].flags & (ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW)) continue;
 
     if (FREEZE_IMAGES && isMosaicChip(image[i].photcode)) continue;

trunk/Ohana/src/relphot/src/MagResidSave.c

@@ -192 +192 @@
   CHECK_STATUS (!status, "ERROR: problem closing meanmags file %s\n", filename);
 
+  FitDataSetFree (&psfStars);
+
   return TRUE;
 }
@@ -213 +215 @@
     off_t c = ImageToCatalog[j];
       
-    // XXX drop this?
     if (catalog[c].measureT[m].dbFlags & MEAS_BAD) continue;
 
@@ -225 +226 @@
     if (isnan(Mgrp)) continue;
 
-    float Mgrid = getMgrid (m, c);
+    float Mgrid = getMgridTiny (&catalog[c].measureT[m]);
     if (isnan(Mgrid)) continue;
 
@@ -231 +232 @@
     // apply that offset as well here for this image (in other words, each detection is
     // being compared to the model, excluding the zero point, Mcal.  The model includes
-    // the flat-correction.  NOTE the sign of Mflat (Image.Mcal = Measure.Mcal - Mflat)
-
-    // float Mflat = getMflat (m, c, flatcorr, catalog);
+    // the flat-correction.  NOTE the sign of Mflat (Image.Mcal = Measure.Mcal + Mflat)
+
+    float Mflat = getMflat (m, c, catalog);
 
     // MrelPSF is the average magnitude for this star.  
@@ -246 +247 @@
     if (isnan(MsysPSF)) continue;
 
-    // XXX skip Mflat for now
-    // float Moff = Mcal + Mmos + Mgrp + Mgrid - Mflat;
-
-    float Moff = Mcal + Mmos + Mgrp + Mgrid;
+    float Moff = Mcal + Mmos + Mgrp + Mgrid + Mflat;
 
     // Mrel is the true average apparent magnitude of this star

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

@@ -1 +1 @@
 # include "relphot.h"
 off_t findMosaic (unsigned int *startTimes, off_t Nmosaic, unsigned int start);
-int save_test_mosaic_measures (FILE *fout, int Nmos, Catalog *catalog, FlatCorrectionTable *flatcorr);
+int save_test_mosaic_measures (FILE *fout, int Nmos, Catalog *catalog);
 
 // see discussion in ImagesOps.c re: IDX_T
@@ -50 +50 @@
 # undef COMPARE
 
-}
-
-/* find mosaic frames (unique time periods & photcode name matches mosaic) */
-void initMosaics_old (Image *image, off_t Nimage) {
-
-  off_t i, j, status, found, NMOSAIC, *MosaicN_IMAGE;
-  unsigned int start, stop;
-  char *pname;
-
-  if (!MOSAIC_ZEROPT) return;
-
-  Nmosaic = 0;
-  NMOSAIC = 10;
-  ALLOCATE (mosaic, Mosaic, NMOSAIC);
-
-  ALLOCATE (MosaicToImage, off_t *, NMOSAIC);
-  ALLOCATE (MosaicN_Image, off_t,   NMOSAIC);
-  ALLOCATE (MosaicN_IMAGE, off_t,   NMOSAIC);
-
-  ALLOCATE (ImageToMosaic, off_t, Nimage); // mosaic to which image belongs
-
-  /* a 'mosaic' in relphot is (unlike relastro) a virtual concept: there is no 
-   * entry in the image table that represents this mosaic.  Instead, it is an
-   * internal construct that defines a group of related images 
-   */
-
-  /* generate list of unique mosaics */
-  for (i = 0; i < Nimage; i++) {
-    ImageToMosaic[i] = -1;
-
-    /* select valid mosaic images by photcode */
-    pname = GetPhotcodeNamebyCode (image[i].photcode);
-    status = strncmp (pname, MOSAICNAME, strlen (MOSAICNAME));
-    if (status) continue;
-
-    /* set image time range */
-    start = image[i].tzero - MAX(0.01*image[i].trate*image[i].NY, 1);
-    stop  = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1);
-
-    /* find existing mosaic with this time range */
-    found = FALSE;
-    for (j = 0; !found && (j < Nmosaic); j++) { 
-      if (stop  < mosaic[j].start) continue;
-      if (start > mosaic[j].stop)  continue;
-      found = TRUE;
-
-      // add reference from image to mosaic
-      ImageToMosaic[i] = j;
-
-      /* add image to mosaic image list */
-      MosaicToImage[j][MosaicN_Image[j]] = i;
-      MosaicN_Image[j] ++;
-      if (MosaicN_Image[j] == MosaicN_IMAGE[j]) {
-        MosaicN_IMAGE[j] += 10;
-        REALLOCATE (MosaicToImage[j], off_t, MosaicN_IMAGE[j]);
-      }
-
-    }
-    if (found) continue;
-    
-    /* a new mosaic, define ranges -- preserve the original values incase this image is not used */
-    mosaic[Nmosaic].start     = start;
-    mosaic[Nmosaic].stop      = stop;
-    mosaic[Nmosaic].McalPSF   = 0.0; // note : at the end, mosaic.Mcal is added back to the input images
-    mosaic[Nmosaic].McalAPER  = 0.0; // note : mosaic stores only offsets relative to the original image values
-    mosaic[Nmosaic].dMcal     = 0.0; // note : at the end, mosaic.Mcal is added back to the input images
-    mosaic[Nmosaic].dMsys     = 0.0;
-    mosaic[Nmosaic].McalChiSq = 0.0;// NAN or 0.0?
-    mosaic[Nmosaic].flags     = image[i].flags;
-    mosaic[Nmosaic].secz      = image[i].secz;
-    mosaic[Nmosaic].photcode  = GetPhotcodeEquivCodebyCode (image[i].photcode);
-
-    // XXX do we need to do something about flag consistency across a mosaic?
-
-    /* add image to mosaic image list */
-    MosaicN_IMAGE[Nmosaic] = 10;
-    MosaicN_Image[Nmosaic] = 1;
-    ALLOCATE (MosaicToImage[Nmosaic], off_t, MosaicN_IMAGE[Nmosaic]);
-    MosaicToImage[Nmosaic][0] = i;
-
-    // add reference from image to mosaic
-    ImageToMosaic[i] = Nmosaic;
-    
-    Nmosaic ++;
-    if (Nmosaic == NMOSAIC) {
-      NMOSAIC += 10;
-      REALLOCATE (mosaic,  Mosaic,  NMOSAIC);
-      REALLOCATE (MosaicToImage,  off_t *, NMOSAIC);
-      REALLOCATE (MosaicN_Image, off_t,   NMOSAIC);
-      REALLOCATE (MosaicN_IMAGE, off_t,   NMOSAIC);
-    }
-  }
-
-  // free this or not?
-  free (MosaicN_IMAGE);
-
-  initMosaicGrid (image, Nimage);
-
-  fprintf (stderr, "matched %d images to %d mosaics\n", (int) Nimage, (int) Nmosaic);
-  return;
 }
 
@@ -318 +218 @@
   free (startTimesMosaic);
 
-  initMosaicGrid (subset, Nsubset);
+  initMosaicMcal (subset, Nsubset);
 
   fprintf (stderr, "matched %d images to %d mosaics, %d simple chips not matched to mosaics\n", (int) (Nsubset - Nsimple), (int) Nmosaic, (int) Nsimple);
@@ -460 +360 @@
 
   if (mergeMcal) {
-    initMosaicGrid (image, Nimage);
+    initMosaicMcal (image, Nimage);
   }
 
@@ -511 +411 @@
 }
 
+// this function sets mosaic->coords to the median of the individual chips.  This
+// coordinate frame is used by the parallel region analysis to assign exposures (mosaics)
+// to specific machines by a single center (rather than individual chips)
 void setMosaicCenters (Image *image, off_t Nimage) {
   OHANA_UNUSED_PARAM(Nimage);
@@ -599 +502 @@
     mosaic[i].coords.cdelt1 = 1.0 / 3600.0;
     mosaic[i].coords.cdelt2 = 1.0 / 3600.0;
-
-    mosaic[i].McalPSF   = 0.0; 
-    mosaic[i].McalAPER  = 0.0; 
-    mosaic[i].dMcal     = 0.0; 
-    mosaic[i].dMsys     = 0.0;
-    mosaic[i].McalChiSq = 0.0; 
   }
   return;
 }
 
-void initMosaicGrid (Image *image, off_t Nimage) {
+void initMosaicMcal (Image *image, off_t Nimage) {
   OHANA_UNUSED_PARAM(Nimage);
 
-  /* find max dR, dD range for all mosaics */
-  /* define mosaic.coords to cover dR, dD */
-  /* send results to initGridBins */
-
-  off_t i, j, m, NX, NY;
-  int dXmax, dYmax;
-  double dS, dX, dY;
-  double R, D, Rmin, Rmax, Dmin, Dmax;
   double McalPSF, McalAPER, dMcal, McalChiSq;
 
   fprintf (stderr, "*** moving Mcal from image.Mcal to mosaic.Mcal ***\n");
 
-  dXmax = dYmax = 0.0;
-  for (i = 0; i < Nmosaic; i++) {
-    Dmin = Rmin =  360.0;
-    Dmax = Rmax = -360.0;
-    dS = 0.0;
+  for (off_t i = 0; i < Nmosaic; i++) {
     McalPSF = McalAPER = dMcal = McalChiSq = 0;
-    for (j = 0; j < MosaicN_Image[i]; j++) {
-      m = MosaicToImage[i][j];
-      NX = image[m].NX;
-      NY = image[m].NY;
-      dS += hypot (image[m].coords.cdelt1*image[m].coords.pc1_1, image[m].coords.cdelt1*image[m].coords.pc2_1);
-
-      OhanaProjection proj = GetProjection (image[m].coords.ctype);
-      if ((proj == PROJ_WRP) && !image[m].coords.mosaic) {
-        XY_to_LM (&R, &D, 0.0, 0.0, &image[m].coords);
-      } else {
-        XY_to_RD (&R, &D, 0.0, 0.0, &image[m].coords);
-      }
-      Rmin = MIN (Rmin, R);
-      Rmax = MAX (Rmax, R);
-      Dmin = MIN (Dmin, D);
-      Dmax = MAX (Dmax, D);
-
-      if ((proj == PROJ_WRP) && !image[m].coords.mosaic) {
-        XY_to_LM (&R, &D, (double) NX, 0.0, &image[m].coords);
-      } else {
-        XY_to_RD (&R, &D, (double) NX, 0.0, &image[m].coords);
-      }
-      Rmin = MIN (Rmin, R);
-      Rmax = MAX (Rmax, R);
-      Dmin = MIN (Dmin, D);
-      Dmax = MAX (Dmax, D);
-
-      if ((proj == PROJ_WRP) && !image[m].coords.mosaic) {
-        XY_to_LM (&R, &D, (double) NX, (double) NY, &image[m].coords);
-      } else {
-        XY_to_RD (&R, &D, (double) NX, (double) NY, &image[m].coords);
-      }
-      Rmin = MIN (Rmin, R);
-      Rmax = MAX (Rmax, R);
-      Dmin = MIN (Dmin, D);
-      Dmax = MAX (Dmax, D);
-
-      if ((proj == PROJ_WRP) && !image[m].coords.mosaic) {
-        XY_to_LM (&R, &D, 0.0, (double) NY, &image[m].coords);
-      } else {
-        XY_to_RD (&R, &D, 0.0, (double) NY, &image[m].coords);
-      }
-      Rmin = MIN (Rmin, R);
-      Rmax = MAX (Rmax, R);
-      Dmin = MIN (Dmin, D);
-      Dmax = MAX (Dmax, D);
-
-      // XXX : this probably does not handle mosaics at RA = 0,360 well
-
-      /* we are using mosaic.Mcal, not image.Mcal. reset image.Mcal */
-
-      // XXX: how does this work with UBERCAL?  We want to keep the Mcal values supplied by ubercal, but 
-      // solve for a single offset for each exposure (Mosaic.Mcal).  
-      // we also want to keep the flat-field terms for each exposure (regardless of ubercal or not)
-      // if it helps, note that ubercal uses a single zp per exposure, so the mean of those values is the same as the value
+    for (off_t j = 0; j < MosaicN_Image[i]; j++) {
+      off_t m = MosaicToImage[i][j];
 
       if (!isfinite(image[m].McalPSF)) {
-        image[m].McalPSF = 0.0;
-        image[m].McalAPER = 0.0;
-        image[m].dMcal = 0.0;
+        image[m].McalPSF   = 0.0;
+        image[m].McalAPER  = 0.0;
+        image[m].dMcal     = 0.0;
         image[m].McalChiSq = 0.0;
         fprintf (stderr, "warning: resetting NAN value for Mcal %s\n", image[m].name);
@@ -703 +535 @@
       image[m].dMcal     = NAN;
       image[m].McalChiSq = NAN;
-
-    }
-    dS /= MosaicN_Image[i];
-
-    InitCoords (&mosaic[i].coords, "DEC--TAN");
-    mosaic[i].coords.crval1 = Rmin;
-    mosaic[i].coords.crval2 = Dmin;
-    mosaic[i].coords.cdelt1 = dS;
-    mosaic[i].coords.cdelt2 = dS;
-
-    RD_to_XY (&dX, &dY, Rmax, Dmax, &mosaic[i].coords);
+    }
 
     mosaic[i].McalPSF   = McalPSF / MosaicN_Image[i];
@@ -720 +542 @@
     mosaic[i].McalChiSq = McalChiSq / MosaicN_Image[i];
   }
-  if (!USE_GRID) return;
-
-  dXmax = MAX (dXmax, dX);
-  dYmax = MAX (dYmax, dY);
-  initGrid (dXmax, dYmax);
   return;
 }
@@ -963 +780 @@
   if (i == -1) return (0.0);
 
-  if (mosaic[i].flags & IMAGE_BAD) return (NAN);  
+  // if the mosaic cannot be calibrated (too few measurements, skip it)
+  // XXX: if the mosaic is bad (few) but the chips are OK, I can used them.
+  // XXX: if (mosaic[i].flags & ID_IMAGE_PHOTOM_FEW) return (NAN);  
   value = mosaic[i].McalPSF;
   return (value);
@@ -1001 +820 @@
   Catalog *catalog;
   Image *image;
-  FlatCorrectionTable *flatcorr;
   SetMmosInfo info;
 } ThreadInfo;
 
-int setMmos_mosaic (Mosaic *mosaic, off_t Nmos, Image *image, Catalog *catalog, SetMmosInfo *info, FlatCorrectionTable *flatcorr);
+int setMmos_mosaic (Mosaic *mosaic, off_t Nmos, Image *image, Catalog *catalog, SetMmosInfo *info);
 void *setMmos_worker (void *data);
-int setMmos_threaded (Catalog *catalog, FlatCorrectionTable *flatcorr);
+int setMmos_threaded (Catalog *catalog);
 
 void SetMmosInfoInit (SetMmosInfo *info, off_t Nmax, int allocLists) {
@@ -1073 +891 @@
 }
 
-int setMmos (Catalog *catalog, FlatCorrectionTable *flatcorr) {
+int setMmos (Catalog *catalog) {
 
   off_t i, N, Nmax;
@@ -1086 +904 @@
   // so do not run setMmos in threaded mode if PLOTSTUFF is set
   if (NTHREADS && !PLOTSTUFF) {
-    int status = setMmos_threaded (catalog, flatcorr);
+    int status = setMmos_threaded (catalog);
     return status;
   }
@@ -1106 +924 @@
 
   for (i = 0; i < Nmosaic; i++) {
-    setMmos_mosaic (&mosaic[i], i, image, catalog, &info, flatcorr);
+    setMmos_mosaic (&mosaic[i], i, image, catalog, &info);
   }
   SetMmosInfoFree (&info);
@@ -1123 +941 @@
   
 // 'mosaic' is a pointer to the current mosaic of interest (Nmos)
-int setMmos_mosaic (Mosaic *myMosaic, off_t Nmos, Image *image, Catalog *catalog, SetMmosInfo *info, FlatCorrectionTable *flatcorr) {
+int setMmos_mosaic (Mosaic *myMosaic, off_t Nmos, Image *image, Catalog *catalog, SetMmosInfo *info) {
 
   off_t j, NimageReal;
@@ -1204 +1022 @@
     snprintf (filename, 64, "test.%05d.%02d.dat", (int) Nmos, npass_output);
     fout = fopen (filename, "w");
-    save_test_mosaic_measures (fout, Nmos, catalog, flatcorr);
+    save_test_mosaic_measures (fout, Nmos, catalog);
     fclose (fout);
   }
@@ -1231 +1049 @@
       continue;
     }
-    float Mgrid = getMgrid (m, c);
+    float Mgrid = getMgridTiny (&catalog[c].measureT[m]);
     if (isnan(Mgrid)) {
       info->Ngrid ++;
@@ -1253 +1071 @@
     // apply that offset as well here for this image (in other words, each detection is
     // being compared to the model, excluding the zero point, Mcal.  The model includes
-    // the flat-correction.  NOTE the sign of Mflat (Image.Mcal = Measure.Mcal - Mflat)
-
-    float Mflat = getMflat (m, c, flatcorr, catalog);
+    // the flat-correction.  NOTE the sign of Mflat (Image.Mcal = Measure.Mcal + Mflat).
+    // this was inconsistent pre r41606
+
+    float Mflat = getMflat (m, c, catalog);
 
     off_t n = catalog[c].measureT[m].averef;
@@ -1278 +1097 @@
     assert (Nbright >= 0);
 
-    float Moff =  Mcal + Mgrp + Mgrid - Mflat;
+    float Moff =  Mcal + Mgrp + Mgrid + Mflat;
 
     psfStars->alldata-> yVector[N] = MsysPSF - MrelPSF - Moff;
@@ -1370 +1189 @@
 }
   
-int save_test_mosaic_measures (FILE *fout, int Nmos, Catalog *catalog, FlatCorrectionTable *flatcorr) {
+int save_test_mosaic_measures (FILE *fout, int Nmos, Catalog *catalog) {
 
   int Nsecfilt = GetPhotcodeNsecfilt ();
@@ -1382 +1201 @@
     float Mcal     = getMcal  (m, c, MAG_CLASS_PSF);
     // float Mgrp     = getMgrp  (m, c, catalog[c].measureT[m].airmass, NULL);
-    float Mgrid    = getMgrid (m, c);
+    float Mgrid    = getMgridTiny (&catalog[c].measureT[m]);
     float MrelPSF  = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP);
-    float Mflat    = getMflat (m, c, flatcorr, catalog);
+    float Mflat    = getMflat (m, c, catalog);
 
     off_t n = catalog[c].measureT[m].averef;
     float MsysPSF = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF);
     
-    float delta = MsysPSF - MrelPSF - Mcal - Mgrid + Mflat;
+    float delta = MsysPSF - MrelPSF - Mcal - Mgrid - Mflat;
 
     int isBad = (catalog[c].measureT[m].dbFlags & MEAS_BAD);
@@ -1398 +1217 @@
 }
 
-int setMmos_threaded (Catalog *catalog, FlatCorrectionTable *flatcorr) {
+int setMmos_threaded (Catalog *catalog) {
 
   int i;
@@ -1435 +1254 @@
     threadinfo[i].catalog  =  catalog;
     threadinfo[i].image    =    image;
-    threadinfo[i].flatcorr = flatcorr;
 
     // we do NOT allocate the arrays here, we only supply basic info (Nmax, Nloop) used in
@@ -1481 +1299 @@
            summary.Nskip);
   free (threadinfo);
+  // XXX SetMmosInfoFree (&summary);
 
   npass_output ++;
@@ -1503 +1322 @@
 
     Catalog *catalog = threadinfo->catalog;
-    FlatCorrectionTable *flatcorr = threadinfo->flatcorr;
     Image *image = threadinfo->image;
 
-    setMmos_mosaic (&mosaic[i], i, image, catalog, &results, flatcorr);
+    setMmos_mosaic (&mosaic[i], i, image, catalog, &results);
     SetMmosInfoAccum (&threadinfo->info, &results);
   }
@@ -1512 +1330 @@
   SetMmosInfoFree (&results);
   return NULL;
-}
-
-// When we rationalize the images/mosaics, we are driving the negative cloud images back
-// to 0.0.  At the same time, we make a guess to the effective impact on all other images,
-// driven by the coupling of common stars. 
-// XXX this function is no longer used because we force significantly negative clouds to 0.0
-int rationalize_mosaics (Catalog *catalog, int Ncatalog) {
-
-  double *imageOffset, **starOffset;
-  int **starNcount, *seclist;
-  int **Slist, *NSlist, *NSLIST;
-  int i, j, k, m, nMos, Ns, found;
-
-  off_t Nimage;
-  Image *image;
-
-  // set a test value for now
-  float CLOUD_TOLERANCE = 0.025;
-
-  if (!MOSAIC_ZEROPT) return (FALSE);
-  if (FREEZE_MOSAICS) return (FALSE);
-
-  image = getimages (&Nimage, NULL);
-
-  ALLOCATE (imageOffset, double, Nmosaic);
-  ALLOCATE ( Slist, int *, Nmosaic); // array of calibrated star indexes on this mosaic
-  ALLOCATE (NSlist, int,   Nmosaic); // number of stars on mosaic
-  ALLOCATE (NSLIST, int,   Nmosaic); // number of Slist entries allocated
-  memset (Slist, 0, Nmosaic*sizeof(int *));
-
-  // find the images / mosaics with negative clouds and save their offset
-  for (i = 0; i < Nmosaic; i++) {
-  
-    NSlist[i] =   0;
-    NSLIST[i] = 100;
-    ALLOCATE (Slist[i], int, NSLIST[i]);
-
-    imageOffset[i] = 0.0;
-
-    if (VERBOSE2 && (fabs(mosaic[i].McalPSF) < CLOUD_TOLERANCE)) {
-      fprintf (stderr, "cloud-free: %s : %f\n", image[MosaicToImage[i][0]].name, mosaic[i].McalPSF);
-    } 
-    if (VERBOSE2 && (mosaic[i].McalPSF < -CLOUD_TOLERANCE)) {
-      imageOffset[i] = -mosaic[i].McalPSF; 
-      // NOTE the negative sign: down below, we are going to add in the negative of McalPSF
-      // to this image, and the propagated mean values for other images
-      fprintf (stderr, "anti-clouds: %s : %f\n", image[MosaicToImage[i][0]].name, mosaic[i].McalPSF);
-    } 
-    if (VERBOSE2 && (mosaic[i].McalPSF > CLOUD_TOLERANCE)) {
-      fprintf (stderr, "cloudy    : %s : %f\n", image[MosaicToImage[i][0]].name, mosaic[i].McalPSF);
-    } 
-  }
-
-  int Nsecfilt = GetPhotcodeNsecfilt ();
-  ALLOCATE (seclist, int, Nphotcodes);
-  for (Ns = 0; Ns < Nphotcodes; Ns ++) { 
-    int thisCode = photcodes[Ns][0].code;
-    seclist[Ns] = GetPhotcodeNsec(thisCode);
-  }
-
-  // allocate an array for star offsets
-  int Nstars = 0;
-  int NSTARS = 1000;
-  ALLOCATE (starOffset, double *, NSTARS);
-  ALLOCATE (starNcount, int *,    NSTARS);
-  memset (starOffset, 0, NSTARS*sizeof(double *));
-  memset (starNcount, 0, NSTARS*sizeof(int *));
-
-  // find the mean offset for each star
-  for (i = 0; i < Ncatalog; i++) {
-    for (j = 0; j < catalog[i].Naverage; j++) {
-      ALLOCATE (starOffset[Nstars], double, Nphotcodes);
-      ALLOCATE (starNcount[Nstars], int,    Nphotcodes);
-      memset (starOffset[Nstars], 0, Nphotcodes*sizeof(double));
-      memset (starNcount[Nstars], 0, Nphotcodes*sizeof(int));
-      
-      m = catalog[i].averageT[j].measureOffset;
-
-      // determine the mosaic for each measurement
-      for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) {
-
-        // skip unused measurements
-        if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue;
-
-        // skip unused measurements
-
-        // XXX replace with:
-        // Ns = GetActivePhotcodeIndex (photcode);
-        // if (Ns < 0) continue;
-
-        int Nsec;
-        int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode);
-        found = FALSE;
-        for (Ns = 0; !found && (Ns < Nphotcodes); Ns ++) { 
-          if (ecode == photcodes[Ns][0].code) {
-            found = TRUE;
-            break;
-          }
-        }
-        if (!found) continue;
-        Nsec = seclist[Ns];
-
-        // bad stars for this secfilt
-        if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;  
-
-        // skip REF measurements (not tied to an image)
-        if (getImageEntry (m, i) < 0) continue;
-
-        // find the source of this measurement (skip unassigned measurements)
-        nMos = MeasureToMosaic[i][m];
-        if (nMos == -1) continue;
-
-        // XXX how can this not be true?  
-        if (mosaic[nMos].photcode != ecode) {
-          fprintf (stderr, "*");
-        }
-
-        assert (Ns >= 0);
-        assert (Ns < Nphotcodes);
-
-        // accumulate the offsets from the negative cloud images (others have 0.0 value)
-        starOffset[Nstars][Ns] += imageOffset[nMos];
-        starNcount[Nstars][Ns] ++;
-        
-        // record the mosaic->star reference
-        Slist[nMos][NSlist[nMos]] = Nstars;
-        NSlist[nMos] ++;
-        if (NSlist[nMos] == NSLIST[nMos]) {
-          NSLIST[nMos] += 100;
-          REALLOCATE (Slist[nMos], int, NSLIST[nMos]);
-        }         
-      }
-      Nstars ++;
-      if (Nstars == NSTARS) {
-        NSTARS += 1000;
-        REALLOCATE (starOffset, double *, NSTARS);
-        REALLOCATE (starNcount, int *,    NSTARS);
-        memset (&starOffset[NSTARS-1000], 0, 1000*sizeof(double *));
-        memset (&starNcount[NSTARS-1000], 0, 1000*sizeof(int *));
-      }
-    }
-  }
-
-  // find the mean offset of the images without negative clouds
-  for (i = 0; i < Nmosaic; i++) {
-
-    found = FALSE;
-    for (Ns = 0; !found && (Ns < Nphotcodes); Ns ++) { 
-      if (mosaic[i].photcode == photcodes[Ns][0].code) {
-        found = TRUE;
-        break;
-      }
-    }
-    if (!found) {
-      fprintf (stderr, "invalid photcode for mosaic?\n");
-      abort();
-    }
-
-    // a negative cloud image (cloud: McalPSF > 0; anti-clouds: McalPSF < 0; imageOffset = -McalPSF)
-    if (imageOffset[i] > 0.0) continue;
-
-    // we need to actually have cross-references to count
-    if (NSlist[i] < 2) continue;
-
-    float dM = 0.0;
-    for (j = 0; j < NSlist[i]; j++) {
-      k = Slist[i][j];
-      if (starNcount[k][Ns] > 1) {
-        dM += (starOffset[k][Ns] / starNcount[k][Ns]);
-      }
-    }
-    imageOffset[i] = dM / NSlist[i];
-    // fprintf (stderr, "adjust image: %s : (%f %d) : %f\n", image[MosaicToImage[i][0]].name, dM, NSlist[i], imageOffset[i]);
-  }
-
-  // for (i = 0; i < Nmosaic; i++) {
-  //   fprintf (stderr, "correction: %s : %f\n", image[MosaicToImage[i][0]].name, imageOffset[i]);
-  // } 
-
-  // apply all offset values to the mosaics
-  // find the images / mosaics with negative clouds and save their offset
-  for (i = 0; i < Nmosaic; i++) {
-    mosaic[i].McalPSF += imageOffset[i];
-  }
-
-  for (i = 0; i < Nstars; i++) {
-    free (starOffset[i]);
-    free (starNcount[i]);
-  }
-  free (starOffset);
-  free (starNcount);
-
-  free (seclist);
-  free (NSlist);
-  free (NSLIST);
-  for (i = 0; i < Nmosaic; i++){
-    free (Slist[i]);
-  }
-  free (Slist);
-  free (imageOffset);
-  return (TRUE);
 }
 
@@ -1731 +1348 @@
   n = 0;
   for (i = 0; i < Nmosaic; i++) {
-    if (mosaic[i].flags & IMAGE_BAD) continue;
+    if (mosaic[i].flags & (ID_IMAGE_MOSAIC_POOR | ID_IMAGE_PHOTOM_FEW)) continue;
     if (mosaic[i].skipCal) continue;
     if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue;
@@ -1764 +1381 @@
   n = 0;
   for (i = 0; i < Nmosaic; i++) {
-    if (mosaic[i].flags & IMAGE_BAD) continue;
+    if (mosaic[i].flags & (ID_IMAGE_MOSAIC_POOR | ID_IMAGE_PHOTOM_FEW)) continue;
     if (mosaic[i].skipCal) continue;
     if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue;
@@ -1781 +1398 @@
 }
 
-StatType statsMosaicN (Catalog *catalog, FlatCorrectionTable *flatcorr) {
+StatType statsMosaicN (Catalog *catalog) {
 
   off_t i, j, m, c, n, N;
@@ -1797 +1414 @@
   n = 0;
   for (i = 0; i < Nmosaic; i++) {
-    if (mosaic[i].flags & IMAGE_BAD)  continue;
+    if (mosaic[i].flags & (ID_IMAGE_MOSAIC_POOR | ID_IMAGE_PHOTOM_FEW))  continue;
     if (mosaic[i].skipCal) continue;
     if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue;
@@ -1811 +1428 @@
       Mgrp = getMgrp  (m, c, catalog[c].measureT[m].airmass, NULL);
       if (isnan(Mgrp)) continue;
-      Mgrid = getMgrid (m, c);
+      Mgrid = getMgridTiny (&catalog[c].measureT[m]);
       if (isnan(Mgrid)) continue;
       Mrel = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP);
@@ -1847 +1464 @@
   n = 0;
   for (i = 0; i < Nmosaic; i++) {
-    if (mosaic[i].flags & IMAGE_BAD) continue;
+    if (mosaic[i].flags & (ID_IMAGE_MOSAIC_POOR | ID_IMAGE_PHOTOM_FEW)) continue;
     if (mosaic[i].skipCal) continue;
     if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue;
@@ -1866 +1483 @@
 /* mark mosaic if: abs(Mcal - <Mcal>) too large, dMcal too large */
 
-static float MinMaxChiSq = NAN;
-static float MinMaxScatter = NAN;
+static float MinChiSqLim = NAN;
+static float MaxChiSqLim = NAN;
+static float MinScatterLim = NAN;
+static float MaxScatterLim = NAN;
 
 void clean_mosaics () {
@@ -1877 +1496 @@
 
   if (VERBOSE) fprintf (stderr, "marking poor mosaics\n");
+
+  if (isnan (MaxChiSqLim))   MaxChiSqLim   = MOSAIC_CHISQ;
+  if (isnan (MaxScatterLim)) MaxScatterLim = MOSAIC_SCATTER;
 
   ALLOCATE (mlist, double, Nmosaic);
@@ -1899 +1521 @@
   liststats (mlist, NULL, NULL, N, &stats);
   float ChiSqUpper90 = stats.Upper90; 
-  if (isnan (MinMaxChiSq)) MinMaxChiSq = 2.0*stats.median;
-  float MaxChiSq = MAX (MinMaxChiSq, ChiSqUpper90);
+  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, N, &stats);
   float ScatterUpper90 = stats.Upper90; 
-  if (isnan (MinMaxScatter)) MinMaxScatter = 2.0*stats.median;
-  float MaxScatter = MAX (MinMaxScatter, ScatterUpper90);
-
-  fprintf (stderr, "MOSAIC: Max ChiSq: %f, Max Scatter: %f | ChiSquare Upper 90: %f, Scatter Upper 90: %f\n", MaxChiSq, MaxScatter, ChiSqUpper90, ScatterUpper90);
+  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, "MOSAIC: ChiSqLimit: %f, ScatterLimit: %f | ChiSquare Upper 90: %f, Scatter Upper 90: %f\n", ChiSqLimit, ScatterLimit, ChiSqUpper90, ScatterUpper90);
   
   int Ntotal = 0, Npoor = 0, Nmark = 0, Nscatter = 0, Nfew = 0, Nchisq = 0;
@@ -1930 +1552 @@
       Nfew ++;
     }
-    if (mosaic[i].stdev > MaxScatter) {
+    if (mosaic[i].stdev > ScatterLimit) {
       mark = TRUE;
       Nscatter ++;
     }
-    if (mosaic[i].McalChiSq > MaxChiSq) {
+    if (mosaic[i].McalChiSq > ChiSqLimit) {
       mark = TRUE;
       Nchisq ++;
@@ -1979 +1601 @@
       c = MosaicToCatalog[i][j];
       
-      if (catalog[c].measureT[m].dbFlags & (ID_MEAS_AREA | ID_MEAS_NOCAL)) continue;
+      if (catalog[c].measureT[m].dbFlags & MEAS_BAD) continue;
 
       // ave = catalog[c].measureT[m].averef;

trunk/Ohana/src/relphot/src/Shutdown.c

@@ -10 +10 @@
 int Shutdown (char *format, ...) {  
   va_list argp;
-  char *formatplus;
+  char formatplus[1024];
   
-  ALLOCATE (formatplus, char, strlen(format));
-  strcpy (formatplus, format);
-  strcat (formatplus, "\n");
+  // we cannot allocate this with the ohana allocation tools
+  // because it could deadlock
+  snprintf (formatplus, 1024, "%s\n", format);
 
   va_start (argp, format);
   vfprintf (stderr, formatplus, argp);
-  free (formatplus);
   va_end (argp);
 

trunk/Ohana/src/relphot/src/StarOps.c

@@ -1 +1 @@
 # include "relphot.h"
+
+# define STAR_BAD (ID_SECF_STAR_FEW | ID_SECF_STAR_POOR)
 
 enum {THREAD_RUN, THREAD_DONE};
@@ -8 +10 @@
   Catalog *catalog;
   int Ncatalog;
-  FlatCorrectionTable *flatcorr;
   SetMrelInfo summary;
 } ThreadInfo;
@@ -14 +15 @@
 // ** internal functions:
 void *setMrel_worker (void *data);
-int setMrel_threaded (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr);
-int setMrel_catalog (Catalog *catalog, int Nc, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt);
+int setMrel_threaded (Catalog *catalog, int Ncatalog);
 int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure);
 
@@ -292 +292 @@
 
 // setMrel and setMrelOutput are extremely similar, but have slightly different implications:
-// * setMrel uses the internal Tiny structures only
-// * setMrelOutput skips stars for which there are too few good measurements 
-// * setMrelOutput is meant to be called repeatedly, relaxing the criteria for 'good' on each pass
+// * both call setMrelCatalog, but setMrel set isSetMrelFinal to FALSE
+// * setMrel only uses the internal Tiny structures
+// * setMrel skips stars for which there are too few good measurements
+// -> (setMrelOutput assigns an average mag for all objects, even those with only 1 measurement)
+// * setMrelOutput sets psfQF, psfQFperf, extNsigma, rank info, etc
+// * setMrelOutput sets average Map
 // * setMrelOutput updates 2MASS average flags
 // * setMrelOutput updates average EXT flags (PS1 and 2MASS)
-
-int setMrel (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) {
+// * setMrelOutput sets average stack & warp photometry (if requested), setMrel always skips those
+// * setMrelOutput is NOT threaded it reads/writes catalogs one at a time
+
+int setMrel (Catalog *catalog, int Ncatalog) {
 
   int i;
 
   if (NTHREADS) {
-    int status = setMrel_threaded (catalog, Ncatalog, flatcorr);
+    int status = setMrel_threaded (catalog, Ncatalog);
     return status;
   }
@@ -314 +319 @@
 
   for (i = 0; i < Ncatalog; i++) {
-    setMrel_catalog_alt (catalog, i, FALSE, flatcorr, &results, Nsecfilt);
+    setMrelCatalog (catalog, i, FALSE, &results, Nsecfilt);
     SetMrelInfoAccum (&summary, &results);
   }
@@ -325 +330 @@
 }
 
-// setMrelOutput is NOT threaded because we read/write catalogs one at a time (for now)
-int setMrelOutput (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) {
+int setMrelOutput (Catalog *catalog, int Ncatalog) {
 
   int i;
@@ -337 +341 @@
 
   for (i = 0; i < Ncatalog; i++) {
-    setMrel_catalog_alt (catalog, i, TRUE, flatcorr, &results, Nsecfilt);
+    setMrelCatalog (catalog, i, TRUE, &results, Nsecfilt);
     SetMrelInfoAccum (&summary, &results);
   }
@@ -347 +351 @@
 }
 
-int setMrel_threaded (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) {
+int setMrel_threaded (Catalog *catalog, int Ncatalog) {
 
   int i;
@@ -374 +378 @@
     threadinfo[i].catalog  =  catalog;
     threadinfo[i].Ncatalog = Ncatalog;
-    threadinfo[i].flatcorr = flatcorr;
     SetMrelInfoInit (&threadinfo[i].summary, FALSE);
     pthread_create (&threads[i], NULL, setMrel_worker, &threadinfo[i]);
@@ -435 +438 @@
 
     Catalog *catalog = threadinfo->catalog;
-    FlatCorrectionTable *flatcorr = threadinfo->flatcorr;
-
-    setMrel_catalog_alt (catalog, i, FALSE, flatcorr, &results, Nsecfilt);
+
+    setMrelCatalog (catalog, i, FALSE, &results, Nsecfilt);
     SetMrelInfoAccum (&threadinfo->summary, &results);
   }
@@ -444 +446 @@
   return NULL;
 }
-
-# define SKIP_THIS_MEAS(REASON) {                               \
-    catalog[Nc].measureT[m].dbFlags |= ID_MEAS_SKIP_PHOTOM;     \
-    if (catalog[Nc].measure) {                                  \
-      catalog[Nc].measure [m].dbFlags |= ID_MEAS_SKIP_PHOTOM;   \
-    }                                                           \
-    results->REASON ++;                                         \
-    continue; }
-
 
 int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure) {
@@ -473 +466 @@
 }
 
-/* set measure.Mcal for all measures except ID_MEAS_NOCAL and ID_IMAGE_PHOTOM_NOCAL */
-int setMcalOutput (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) {
+/* set measure.Mcal for all measures except ID_MEAS_NOCAL */
+int setMcalOutput (Catalog *catalog, int Ncatalog) {
 
   int i;
   off_t j, k, m;
 
-  MEAS_BAD = ID_MEAS_NOCAL;
-  IMAGE_BAD = ID_IMAGE_PHOTOM_NOCAL;
-
   for (i = 0; i < Ncatalog; i++) {
     for (j = 0; j < catalog[i].Naverage; j++) {
@@ -487 +477 @@
       m = catalog[i].averageT[j].measureOffset;
       for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) {
-        if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue;
+        if (catalog[i].measureT[m].dbFlags & ID_MEAS_NOCAL) continue;
         float McalPSF  = getMcal  (m, i, MAG_CLASS_PSF);
         float McalAPER = getMcal  (m, i, MAG_CLASS_APER);
         if (isnan(McalPSF)) continue;
-        float Mmos  = getMmos  (m, i);
-        if (isnan(Mmos)) continue;
-        float Mgrp  = getMgrp  (m, i, catalog[i].measureT[m].airmass, NULL);
-        if (isnan(Mgrp)) continue;
-        float Mgrid = getMgrid (m, i);
-        if (isnan(Mgrid)) continue;
+
+        // Mmos, Mgrp have already been transfered to Mcal and should be zero here
+        // float Mmos  = getMmos  (m, i);
+        // if (isnan(Mmos)) continue;
+        // float Mgrp  = getMgrp  (m, i, catalog[i].measureT[m].airmass, NULL);
+        // if (isnan(Mgrp)) continue;
+
+        // Mgrid is applied to Mflat and should be ignored here
+        // float Mgrid = getMgridTiny (&catalog[i].measureT[m]);
+        // if (isnan(Mgrid)) continue;
 
         // Note that this operation is setting measure->McalAPER to image->McalAPER only
@@ -509 +503 @@
 
         // set the output calibration
-        // float Moff = Mmos + Mgrp + Mgrid;
-        float Moff = Mmos + Mgrid; // Mgrp has already been applied to the image.Mcal
-        catalog[i].measure[m].McalPSF  = McalPSF + Moff;
-        catalog[i].measure[m].McalAPER = useStackAper ? McalAPER + Moff : McalPSF + Moff;
+        catalog[i].measure[m].McalPSF  = McalPSF;
+        catalog[i].measure[m].McalAPER = useStackAper ? McalAPER : McalPSF;
 
         if (catalog[i].measureT[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) {
@@ -584 +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;
@@ -597 +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++) {
-        if ( catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD ) continue;
+    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;
         if (isnan(Mchisq)) continue;
         xlist[Ntot] = Mchisq;
         slist[Ntot] = catalog[i].secfilt[Nsecfilt*j+Nsec].dMpsfChp;
-        dlist[Ntot] = 1;
         Ntot ++;
       }
     }
-  
-    liststats (xlist, dlist, NULL, Ntot, &stats);
-    MaxChisq = MAX (STAR_CHISQ, 2*stats.median);
-
-    liststats (slist, dlist, NULL, Ntot, &stats);
-    MaxScatter = MAX (STAR_SCATTER, 2*stats.median);
-    fprintf (stderr, "Max Scatter: %f, Max Chisq: %f\n", MaxScatter, MaxChisq);
+
+    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;
+    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_OBJ_POOR;
+          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_OBJ_POOR;
+          catalog[i].secfilt[Nsecfilt*j+Nsec].flags &= ~ID_SECF_STAR_POOR;
         }
         Nave ++;
       }
     }
-    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);
-}
-
-// clean_measures examines the stats for a single star.  It measures the INNER 50% mean
-// and sigma, it then re-measures the mean and sigma using all stars with NSIGMA_CLIP (3)
-// sigma of the INNER 50% mean.  it then flags any measurements which are more than
-// NSIGMA_REJECT (5) sigma of the mean
-
-// this function only operations on the PSF magnitudes
-
-// XXX EAM clean_measures is no longer used -- we rely on IRLS to down-weight outliers
-
-# define NSIGMA_CLIP 3.0
-# define NSIGMA_REJECT 5.0
-void clean_measures (Catalog *catalog, int Ncatalog, int final, FlatCorrectionTable *flatcorr) {
-
-  off_t j, k, m, Nmax, Ndel, Nave;
-  int i, N, image_bad, TOOFEW;
-  off_t *ilist;
-  double *tlist, *list, *dlist;
-  float Msys, Mcal, Mmos, Mgrp, Mgrid;
-  int Ncal, Nmos, Ngrp, Ngrid, Nfew;
-
-  int Nsecfilt = GetPhotcodeNsecfilt ();
-
-  if (VERBOSE) fprintf (stderr, "marking poor measures\n");
-  /* Nmeasure is now different, need to reallocate */
-  Nmax = 0;
-  for (i = 0; i < Ncatalog; i++) {
-    for (j = 0; j < catalog[i].Naverage; j++) {
-      Nmax = MAX (Nmax, catalog[i].averageT[j].Nmeasure);
-    }
-  }
-  ALLOCATE (list, double, MAX (1, Nmax));
-  ALLOCATE (dlist, double, MAX (1, Nmax));
-  ALLOCATE (ilist, off_t, Nmax);
-  ALLOCATE (tlist, double, Nmax);
-  
-  /* it makes no sense to mark 3-sigma outliers with <5 measurements */
-  TOOFEW = MAX (5, STAR_TOOFEW);
-
-  // stats structures for inner and full stats
-  StatType instats, stats;
-  liststats_setmode (&instats, "INNER_MEAN");
-  liststats_setmode (&stats, "MEAN");
-
-  Ndel = Nave = 0;
-  Ncal = Nmos = Ngrp = Ngrid = Nfew = 0;
-  for (i = 0; i < Ncatalog; i++) {
-    for (j = 0; j < catalog[i].Naverage; j++) {
-
-      int Ns;
-      for (Ns = 0; Ns < Nphotcodes; Ns++) {
-        
-        int thisCode = photcodes[Ns][0].code;
-        int Nsec = GetPhotcodeNsec(thisCode);
-        
-        /* on final processing, skip stars already measured */
-        if (final && catalog[i].found_t[Nsecfilt*j + Nsec]) continue;  
-
-        /* skip bad stars to prevent them from becoming good (on inner sample) */
-        if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;  
-
-        /* accumulate list of valid measurements */
-        m = catalog[i].averageT[j].measureOffset;
-        N = 0;
-        for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) {
-          // skip measurements that do not match the current photcode
-          int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode);
-          if (ecode != thisCode) { continue; }
-
-          // SKIP gpc1 stack data
-          if (isGPC1stack(catalog[i].measureT[m].photcode)) continue;
-          
-          // SKIP gpc1 forced-warp data
-          if (isGPC1warp(catalog[i].measureT[m].photcode)) continue;
-          
-          // NOTE: we do not skip MEAS_BAD because this measurement is just an internal assessment of the outliers
-          Mcal  = getMcal  (m, i, MAG_CLASS_PSF);
-          if (isnan(Mcal)) { Ncal ++; continue; }
-          Mmos  = getMmos  (m, i);
-          if (isnan(Mmos)) { Nmos ++; continue; }
-          Mgrp  = getMgrp  (m, i, catalog[i].measureT[m].airmass, NULL);
-          if (isnan(Mgrp)) { Ngrp ++; continue; }
-          Mgrid = getMgrid (m, i);
-          if (isnan(Mgrid)) { Ngrid ++; continue; }
-
-          Msys = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt], MAG_CLASS_PSF);
-          list[N] = Msys - Mcal - Mmos - Mgrp - Mgrid;
-          dlist[N] = MAX (catalog[i].measureT[m].dM, MIN_ERROR);
-          N++;
-        }
-        if (N <= TOOFEW) { 
-          Nfew ++; 
-          continue; 
-        }
-
-        /* 3-sigma clip based on stats of inner 50% */
-
-        // calculated mean of inner 50%
-        liststats (list, dlist, NULL, N, &instats);
-        instats.sigma = MAX (MIN_ERROR, instats.sigma); /* if measurements agree too well, sigma -> 0.0 */
-
-        // ignore entries > 3sigma from inner mean
-        for (k = m = 0; k < N; k++) {
-          if (fabs (list[k] - instats.median) < NSIGMA_CLIP*instats.sigma) {
-            list[m] = list[k];
-            m++;
-          }
-        }
-        // recalculate the mean & sigma of the accepted measurements
-        liststats (list, dlist, NULL, m, &stats);
-        stats.sigma = MAX (MIN_ERROR, stats.sigma);
-
-        /* apply to list of all relevant measurements, including IMAGE_POOR & IMAGE_FEW */
-        image_bad = IMAGE_BAD;
-        IMAGE_BAD = ID_IMAGE_PHOTOM_NOCAL;
-        m = catalog[i].averageT[j].measureOffset;
-        N = 0;
-        for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) {
-          // skip measurements that do not match the current photcode
-          int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode);
-          if (ecode != thisCode) { continue; }
-
-          // SKIP gpc1 stack data
-          if (isGPC1stack(catalog[i].measureT[m].photcode)) continue;
-          
-          // SKIP gpc1 forced-warp data
-          if (isGPC1warp(catalog[i].measureT[m].photcode)) continue;
-          
-          // NOTE: we do not skip MEAS_BAD because this measurement is just an internal assessment of the outliers
-          Mcal  = getMcal  (m, i, MAG_CLASS_PSF);
-          if (isnan(Mcal)) continue;
-          Mmos  = getMmos  (m, i);
-          if (isnan(Mmos)) continue;
-          Mgrid = getMgrid (m, i);
-          if (isnan(Mgrid)) continue;
-
-          Msys = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt], MAG_CLASS_PSF);
-          list[N] = Msys - Mcal - Mmos - Mgrid;
-          dlist[N] = MAX (catalog[i].measureT[m].dM, MIN_ERROR);
-          ilist[N] = m;
-          N++;
-          Nave ++;
-        }
-        if (N < TOOFEW) continue;
-
-        /* mark bad measures (> 5 sigma deviant) */
-        for (k = 0; k < N; k++) {
-          // treat the scatter of the star as a systematic term.  this is a bit of an
-          // over-estimage (a perfect Gauss distribution with perfect errors would have
-          // mySigma = sqrt(2) too large)
-          float mySigma = hypot (stats.sigma, dlist[k]);
-          if (fabs (list[k] - stats.median) > NSIGMA_REJECT*mySigma) {
-            catalog[i].measureT[ilist[k]].dbFlags |= ID_MEAS_POOR_PHOTOM;
-            if (final) {
-              // for the final pass, we have a duplicate set of values in measure and measureT
-              catalog[i].measure[ilist[k]].dbFlags |= ID_MEAS_POOR_PHOTOM;
-            }
-            Ndel ++;
-          } else {
-            catalog[i].measureT[ilist[k]].dbFlags &= ~ID_MEAS_POOR_PHOTOM;
-            if (final) {
-              // for the final pass, we have a duplicate set of values in measure and measureT
-              catalog[i].measure[ilist[k]].dbFlags &= ~ID_MEAS_POOR_PHOTOM;
-            }
-          }
-        }
-        IMAGE_BAD = image_bad;
-      }
-    }
-  }
-  if (VERBOSE) fprintf (stderr, OFF_T_FMT" measures marked poor, "OFF_T_FMT" total\n", Ndel, Nave);
-
-  free (list);
-  free (dlist);
-  free (ilist);
-  free (tlist);
-}
-
-StatType statsStarN (Catalog *catalog, int Ncatalog, int Nsec, int seccode, FlatCorrectionTable *flatcorr) {
+}
+
+StatType statsStarN (Catalog *catalog, int Ncatalog, int Nsec, int seccode) {
 
   off_t j, k, m, Ntot;
@@ -869 +692 @@
         Mmos = getMmos  (m, i);
         if (isnan(Mmos)) { continue; }
-        Mgrid = getMgrid (m, i);
+        Mgrid = getMgridTiny (&catalog[i].measureT[m]);
         if (isnan(Mgrid)) { continue;}
         N++;
@@ -1105 +928 @@
 
       float MrelPSF  = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); // average magnitude
-      // float Mgrid    = getMgrid (m, c);                    // camera offset (deprecated?)
-      // float Mflat    = getMflat (m, c, flatcorr, catalog); // flat-field correction
+      // float Mgrid    = getMgridTiny (&catalog[c].measureT[m]);                    // camera offset (deprecated?)
+      // float Mflat    = getMflat (m, c, catalog); // flat-field correction
 
       off_t n = catalog[c].measureT[m].averef;
@@ -1121 +944 @@
   fclose (fout);
 }
+
+// maybe make this a macro
+float getMflat (off_t meas, int cat, Catalog *catalog) {
+  float offset = catalog[cat].measureT[meas].Mflat;
+  if (!isfinite(offset)) {
+    offset = 0.0;
+  }
+  return (offset);
+}

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

@@ -2 +2 @@
 
 void assignMosaicsToTGroups (void);
-int save_test_night_measures (FILE *fout, TGroup *myTGroup, Catalog *catalog, FlatCorrectionTable *flatcorr);
+int save_test_night_measures (FILE *fout, TGroup *myTGroup, Catalog *catalog);
 
 // tgroupTimes carries the times of the tgroups (initially, just the photometric nights)
@@ -647 +647 @@
   Catalog *catalog;
   Image *image;
-  FlatCorrectionTable *flatcorr;
   SetMgrpInfo info;
 } ThreadInfo;
 
-int   setMgrp_tgroup (TGroup *tgroup, off_t Nmos, Image *image, Catalog *catalog, SetMgrpInfo *info, FlatCorrectionTable *flatcorr);
+int   setMgrp_tgroup (TGroup *tgroup, off_t Nmos, Image *image, Catalog *catalog, SetMgrpInfo *info);
 void *setMgrp_worker (void *data);
-int   setMgrp_threaded (Catalog *catalog, FlatCorrectionTable *flatcorr);
+int   setMgrp_threaded (Catalog *catalog);
 
 void SetMgrpInfoInit (SetMgrpInfo *info, off_t Nmax, int allocLists) {
@@ -727 +726 @@
 }
 
-int setMgrp (Catalog *catalog, FlatCorrectionTable *flatcorr) {
+int setMgrp (Catalog *catalog) {
 
   off_t N;
@@ -738 +737 @@
   // so do not run setMgrp in threaded mode if PLOTSTUFF is set
   if (NTHREADS && !PLOTSTUFF) {
-    int status = setMgrp_threaded (catalog, flatcorr);
+    int status = setMgrp_threaded (catalog);
     return status;
   }
@@ -758 +757 @@
     TGroup *tgroup = tgroupTimes[i][0].byCode;
     for (off_t j = 0; j < tgroupTimes[i][0].nCode; j++) {
-      setMgrp_tgroup (&tgroup[j], i, image, catalog, &info, flatcorr);
+      setMgrp_tgroup (&tgroup[j], i, image, catalog, &info);
       // fprintf (stderr, "TGROUP time %.0f photcode %d Mcal: %f, dK: %f, dMcal: %f, sigma: %f, chisq: %f, %d of %d, limiting negative clouds to %f\n",
       // ohana_sec_to_mjd(tgroupTimes[i][0].start), tgroup[j].photcode, tgroup[j].McalPSF, tgroup[j].dKlam, tgroup[j].dMcal, tgroup[j].stdev, tgroup[j].McalChiSq, tgroup[j].nFitPhotom, tgroup[j].nValPhotom, CLOUD_TOLERANCE);
@@ -781 +780 @@
 
 // 'tgroup' is a pointer to the current tgroup of interest (entry Ngrp)
-int setMgrp_tgroup (TGroup *myTGroup, off_t Ngrp, Image *image, Catalog *catalog, SetMgrpInfo *info, FlatCorrectionTable *flatcorr) {
+int setMgrp_tgroup (TGroup *myTGroup, off_t Ngrp, Image *image, Catalog *catalog, SetMgrpInfo *info) {
 
   off_t j;
@@ -826 +825 @@
     fout = fopen (filename, "w");
     fprintf (fout, "# tgroup %s (%.0f)\n", ohana_sec_to_date(tgroup->start), mjdStart);
-    save_test_night_measures (fout, myTGroup, catalog, flatcorr);
+    save_test_night_measures (fout, myTGroup, catalog);
     fclose (fout);
   }
@@ -853 +852 @@
       continue;
     }
-    float Mgrid = getMgrid (m, c);
+    float Mgrid = getMgridTiny (&catalog[c].measureT[m]);
     if (isnan(Mgrid)) {
       info->Ngrid ++;
@@ -875 +874 @@
     // apply that offset as well here for this image (in other words, each detection is
     // being compared to the model, excluding the zero point, Mcal.  The model includes
-    // the flat-correction.  NOTE the sign of Mflat (Image.Mcal = Measure.Mcal - Mflat)
-
-    float Mflat = getMflat (m, c, flatcorr, catalog);
+    // the flat-correction.  NOTE the sign of Mflat (Image.Mcal = Measure.Mcal + Mflat)
+    // this was inconsistent with PhotRel pre-r41606
+
+    float Mflat = getMflat (m, c, catalog);
 
     off_t n = catalog[c].measureT[m].averef;
@@ -900 +900 @@
     assert (Nbright >= 0);
 
-    float Moff =  Mcal + Mmos + Mgrid - Mflat;
+    float Moff =  Mcal + Mmos + Mgrid + Mflat;
 
     psfStars->alldata-> yVector[N] = MsysPSF - MrelPSF - Moff;
@@ -1000 +1000 @@
 }
   
-int save_test_night_measures (FILE *fout, TGroup *myTGroup, Catalog *catalog, FlatCorrectionTable *flatcorr) {
+int save_test_night_measures (FILE *fout, TGroup *myTGroup, Catalog *catalog) {
 
   int Nsecfilt = GetPhotcodeNsecfilt ();
@@ -1013 +1013 @@
     float Mmos     = getMmos  (m, c);                    // mosaic zero point
 
-    float Mgrid    = getMgrid (m, c);                    // camera offset (deprecated?)
+    float Mgrid    = getMgridTiny (&catalog[c].measureT[m]);                    // camera offset (deprecated?)
     float MrelPSF  = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); // average magnitude
-    float Mflat    = getMflat (m, c, flatcorr, catalog); // flat-field correction
+    float Mflat    = getMflat (m, c, catalog); // flat-field correction
 
     off_t n = catalog[c].measureT[m].averef;
@@ -1023 +1023 @@
 
     // for historical reasons, Mflat is defined with the wrong sign
-    float delta = MsysPSF - MrelPSF - Mcal - Mmos - Mgrid + Mflat - myTGroup->McalPSF;
+    float delta = MsysPSF - MrelPSF - Mcal - Mmos - Mgrid - Mflat - myTGroup->McalPSF;
 
     int isBad = (catalog[c].measureT[m].dbFlags & MEAS_BAD);
@@ -1037 +1037 @@
 }
 
-int setMgrp_threaded (Catalog *catalog, FlatCorrectionTable *flatcorr) {
+int setMgrp_threaded (Catalog *catalog) {
 
   int i;
@@ -1075 +1075 @@
     threadinfo[i].catalog  =  catalog;
     threadinfo[i].image    =    image;
-    threadinfo[i].flatcorr = flatcorr;
 
     // we do NOT allocate the arrays here, we only supply basic info used in the threads
@@ -1124 +1123 @@
   free (threadinfo);
 
+  // XXX not allocated SetMgrpInfoFree (&summary);
+
   npass_output ++;
 
@@ -1145 +1146 @@
 
     Catalog *catalog = threadinfo->catalog;
-    FlatCorrectionTable *flatcorr = threadinfo->flatcorr;
     Image *image = threadinfo->image;
 
     TGroup *tgroup = tgroupTimes[i][0].byCode;
     for (int j = 0; j < tgroupTimes[i][0].nCode; j++) {
-      setMgrp_tgroup (&tgroup[j], i, image, catalog, &results, flatcorr);
+      setMgrp_tgroup (&tgroup[j], i, image, catalog, &results);
       SetMgrpInfoAccum (&threadinfo->info, &results);
     }
@@ -1172 +1172 @@
   ALLOCATE (dlist, double, NtgroupTimes*Nphotcodes);
 
-  int NIGHT_BAD = ID_IMAGE_NIGHT_POOR | ID_IMAGE_PHOTOM_FEW;
-
   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].flags & NIGHT_BAD) continue;
+      if (tgroup[j].flags & (ID_IMAGE_NIGHT_POOR | ID_IMAGE_PHOTOM_FEW)) continue;
       list[n] = tgroup[j].McalPSF;
       dlist[n] = 1;
@@ -1209 +1207 @@
     TGroup *tgroup = tgroupTimes[i][0].byCode;
     for (int j = 0; j < tgroupTimes[i][0].nCode; j++) {
-      if (tgroup[j].flags & IMAGE_BAD) continue;
+      if (tgroup[j].flags & (ID_IMAGE_NIGHT_POOR | ID_IMAGE_PHOTOM_FEW)) continue;
       list[n] = tgroup[j].dMcal;
       dlist[n] = 1;
@@ -1240 +1238 @@
     TGroup *tgroup = tgroupTimes[i][0].byCode;
     for (int j = 0; j < tgroupTimes[i][0].nCode; j++) {
-      if (tgroup[j].flags & IMAGE_BAD) continue;
+      if (tgroup[j].flags & (ID_IMAGE_NIGHT_POOR | ID_IMAGE_PHOTOM_FEW)) continue;
       list[n] = tgroup[j].McalChiSq;
       dlist[n] = 1;
@@ -1267 +1265 @@
  */
 
-static float MinMaxChiSq = NAN;
-static float MinMaxScatter = NAN;
+static float MinChiSqLim = NAN;
+static float MaxChiSqLim = NAN;
+static float MinScatterLim = NAN;
+static float MaxScatterLim = NAN;
 
 void clean_tgroups () {
@@ -1277 +1277 @@
 
   if (VERBOSE) fprintf (stderr, "marking poor tgroups\n");
+
+  if (isnan (MaxChiSqLim))   MaxChiSqLim   = NIGHT_CHISQ;
+  if (isnan (MaxScatterLim)) MaxScatterLim = NIGHT_SCATTER;
 
   ALLOCATE (mlist, double, NtgroupTimes*Nphotcodes);
@@ -1300 +1303 @@
   liststats (mlist, NULL, NULL, N, &stats);
   float ChiSqUpper90 = stats.Upper90; 
-  if (isnan (MinMaxChiSq)) MinMaxChiSq = 2.0*stats.median;
-  float MaxChiSq = MAX (MinMaxChiSq, ChiSqUpper90);
+  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, N, &stats);
   float ScatterUpper90 = stats.Upper90; 
-  if (isnan (MinMaxScatter)) MinMaxScatter = 2.0*stats.median;
-  float MaxScatter = MAX (MinMaxScatter, ScatterUpper90);
-
-  fprintf (stderr, "TGROUPS: Max ChiSq: %f, Max Scatter: %f | ChiSquare Upper 90: %f, Scatter Upper 90: %f\n", MaxChiSq, MaxScatter, ChiSqUpper90, ScatterUpper90);
+  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, "TGROUPS: ChiSqLimit: %f, ScatterLimit: %f | ChiSquare Upper 90: %f, Scatter Upper 90: %f\n", ChiSqLimit, ScatterLimit, ChiSqUpper90, ScatterUpper90);
   
   int Ntotal = 0, Npoor = 0, Nmark = 0, Nscatter = 0, Nchisq = 0, NfewNights = 0, NfewExp = 0;
@@ -1337 +1340 @@
       }
       // scatter too large
-      if (tgroup[j].stdev > MaxScatter) {
+      if (tgroup[j].stdev > ScatterLimit) {
         mark = TRUE;
         Nscatter ++;
       }
       // chisq too large
-      if (tgroup[j].McalChiSq > MaxChiSq) {
+      if (tgroup[j].McalChiSq > ChiSqLimit) {
         mark = TRUE;
         Nchisq ++;
@@ -1394 +1397 @@
         c = tgroup[i].catalog[j];
       
-        if (catalog[c].measureT[m].dbFlags & (ID_MEAS_AREA | ID_MEAS_NOCAL)) continue;
+        if (catalog[c].measureT[m].dbFlags & MEAS_BAD) continue;
 
         // ave = catalog[c].measureT[m].averef;
@@ -1501 +1504 @@
 
         float MrelPSF  = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP); // average magnitude
-        // float Mgrid    = getMgrid (m, c);                    // camera offset (deprecated?)
-        // float Mflat    = getMflat (m, c, flatcorr, catalog); // flat-field correction
+        // float Mgrid    = getMgridTiny (&catalog[c].measureT[m]); // camera offset (deprecated?)
+        float Mflat    = getMflat (m, c, catalog); // flat-field correction
 
         off_t n = catalog[c].measureT[m].averef;
@@ -1510 +1513 @@
 
         // for historical reasons, Mflat is defined with the wrong sign
-        float delta = MsysPSF - MrelPSF - Mcal - Mmos - Mgrp;
+        float delta = MsysPSF - MrelPSF - Mcal - Mmos - Mgrp - Mflat;
 
         fprintf (fout, "%f %f : %f %f : %f %f %f  : %f %f\n", catalog[c].averageT[n].R, catalog[c].averageT[n].D, MsysPSF, MrelPSF, Mcal, Mmos, Mgrp, catalog[c].measureT[m].airmass, delta);

trunk/Ohana/src/relphot/src/ZeroPointModes.c

@@ -22 +22 @@
 }
 
-void SetZeroPointModes (void) {
+void SetZeroPointModes (Catalog *catalog, int Ncatalog) {
 
   if (TGROUP_ZEROPT && MOSAIC_ZEROPT) {
@@ -49 +49 @@
       clean_mosaics(); // do this on every pass 
     }
-    if ((CurrentLoop >= 18) && (CurrentLoop <= 999)) {
+    if ((CurrentLoop >= 18) && (CurrentLoop < 21)) {
       // after iterating a few times on the TGroup & Mosaic zero points:
       // * identify the good / bad Mosaics
@@ -59 +59 @@
       clean_mosaics(); // do this on every pass 
       clean_images(); // do this on every pass 
+      clean_stars(catalog, Ncatalog); // do this on every pass 
+    }
+    if ((CurrentLoop >= 21) && (CurrentLoop <= 999)) {
+      // after iterating a few times on the TGroup & Mosaic zero points:
+      // * identify the good / bad Mosaics
+      // * fit the chips from the bad Mosaics
+      TGROUP_ZPT_MODE = TGROUP_ZPT_MODE_GOOD_NIGHT; // stop fitting the bad nights
+      MOSAIC_ZPT_MODE = MOSAIC_ZPT_MODE_BAD_NIGHT_GOOD_MOSAIC; // only fit good mosaics on bad nights
+      IMAGE_ZPT_MODE  = IMAGE_ZPT_MODE_BAD_NIGHT_BAD_MOSAIC; // only fit bad mosaics on bad nights
+      GRID_ZPT_MODE  = GRID_ZPT_MODE_ALL; // only fit bad mosaics on bad nights
+      clean_tgroups(); // do this on every pass 
+      clean_mosaics(); // do this on every pass 
+      clean_images(); // do this on every pass 
+      clean_stars(catalog, Ncatalog); // do this on every pass 
     }
   }
@@ -68 +82 @@
       MOSAIC_ZPT_MODE = MOSAIC_ZPT_MODE_NONE;
     }
-    if ((CurrentLoop > 4) && (CurrentLoop <= 999)) {
+    if ((CurrentLoop > 4) && (CurrentLoop <= 8)) {
       // after iterating a few times on the TGroup zero points:
       // * identify the photometric nights
@@ -77 +91 @@
       clean_tgroups(); // do this on every pass or just sometimes?
     }
+    if ((CurrentLoop > 8) && (CurrentLoop <= 999)) {
+      // after iterating a few times on the TGroup zero points:
+      // * identify the photometric nights
+      // * fit the images from the non-photometric nights
+      TGROUP_ZPT_MODE = TGROUP_ZPT_MODE_GOOD_NIGHT; // stop fitting the bad nights
+      MOSAIC_ZPT_MODE = MOSAIC_ZPT_MODE_NONE;
+      IMAGE_ZPT_MODE  = IMAGE_ZPT_MODE_BAD_NIGHT;
+      GRID_ZPT_MODE  = GRID_ZPT_MODE_ALL; // only fit bad mosaics on bad nights
+      clean_tgroups(); // do this on every pass or just sometimes?
+    }
   }
 
   if (!TGROUP_ZEROPT && MOSAIC_ZEROPT) {
-    if (CurrentLoop <= 4) {
+    if (CurrentLoop < 6) {
       TGROUP_ZPT_MODE = TGROUP_ZPT_MODE_NONE;
       MOSAIC_ZPT_MODE = MOSAIC_ZPT_MODE_ALL;
       IMAGE_ZPT_MODE  = IMAGE_ZPT_MODE_NONE;
     }
-    if ((CurrentLoop > 4) && (CurrentLoop <= 999)) {
+    if ((CurrentLoop >= 6) && (CurrentLoop < 12)) {
       // after iterating a few times on the Mosaic zero points:
       // * identify the good / bad mosaics
@@ -94 +118 @@
       clean_mosaics(); // do this on every pass or just sometimes?
     }
+    if ((CurrentLoop >= 12) && (CurrentLoop < 18)) {
+      // after iterating a few times on the Mosaic zero points:
+      // * identify the good / bad mosaics
+      // * fit the images from the bad mosaics
+      TGROUP_ZPT_MODE = TGROUP_ZPT_MODE_NONE; // stop fitting the bad nights
+      MOSAIC_ZPT_MODE = MOSAIC_ZPT_MODE_GOOD_MOSAIC;
+      IMAGE_ZPT_MODE  = IMAGE_ZPT_MODE_BAD_MOSAIC;
+      clean_mosaics(); // do this on every pass or just sometimes?
+      clean_images(); // do this on every pass 
+      clean_stars(catalog, Ncatalog); // do this on every pass 
+    }
+    if ((CurrentLoop >= 18) && (CurrentLoop < 999)) {
+      // after iterating a few times on the Mosaic zero points:
+      // * identify the good / bad mosaics
+      // * fit the images from the bad mosaics
+      TGROUP_ZPT_MODE = TGROUP_ZPT_MODE_NONE; // stop fitting the bad nights
+      MOSAIC_ZPT_MODE = MOSAIC_ZPT_MODE_GOOD_MOSAIC;
+      IMAGE_ZPT_MODE  = IMAGE_ZPT_MODE_BAD_MOSAIC;
+      GRID_ZPT_MODE  = GRID_ZPT_MODE_ALL; // only fit bad mosaics on bad nights
+      clean_mosaics(); // do this on every pass or just sometimes?
+      clean_images(); // do this on every pass 
+      clean_stars(catalog, Ncatalog); // do this on every pass 
+    }
   }
 

trunk/Ohana/src/relphot/src/args.c

@@ -198 +198 @@
   }
 
-  NGRID = 8;
-  if ((N = get_argument (argc, argv, "-ngrid"))) {
-    remove_argument (N, &argc, argv);
-    NGRID = atof (argv[N]);
-    remove_argument (N, &argc, argv);
-  }
-
   TEST_IMAGE1 = -1;
   if ((N = get_argument (argc, argv, "-test-image1"))) {
@@ -227 +220 @@
     RESET = TRUE;
   }
-
   RESET_ZEROPTS = FALSE;
   if ((N = get_argument (argc, argv, "-reset-zpts"))) {
@@ -233 +225 @@
     RESET_ZEROPTS = TRUE;
   }
+  RESET_FLATCORR = FALSE;
+  if ((N = get_argument (argc, argv, "-reset-flat"))) {
+    remove_argument (N, &argc, argv);
+    RESET_FLATCORR = TRUE;
+  }
 
   REPAIR_WARPS = FALSE;
@@ -244 +241 @@
     remove_argument (N, &argc, argv);
     PRESERVE_PS1 = TRUE;
+  }
+  REQUIRE_PSFFIT = FALSE;
+  if ((N = get_argument (argc, argv, "-require-psffit"))) {
+    remove_argument (N, &argc, argv);
+    REQUIRE_PSFFIT = TRUE;
   }
 
@@ -387 +389 @@
   }
 
-  // USE_GRID is not valid for all cases, probably should be its own mode...
-  USE_GRID = FALSE;
+  // GRID_ZEROPT is not valid for all cases, probably should be its own mode...
+  GRID_ZEROPT = FALSE;
+  GRID_MEANFILE = NULL; // this is set by GridCorrectionSave()
+  GRID_ZPT_MODE = GRID_ZPT_MODE_NONE; // start with grid off
   if ((N = get_argument (argc, argv, "-grid"))) {
     remove_argument (N, &argc, argv);
-    USE_GRID = TRUE;
+    GRID_ZEROPT = TRUE;
+  }
+  if ((N = get_argument (argc, argv, "-grid-meanfile"))) {
+    GRID_ZEROPT = TRUE;
+    remove_argument (N, &argc, argv);
+    GRID_MEANFILE = strcreate (argv[N]); // used by relphot -averages and relphot -apply-offsets
+    remove_argument (N, &argc, argv);
+  }
+
+  GRID_BIN_GPC1 = 16;
+  if ((N = get_argument (argc, argv, "-grid-bin-gpc1"))) {
+    remove_argument (N, &argc, argv);
+    GRID_BIN_GPC1 = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  GRID_BIN_GPC2 = 16;
+  if ((N = get_argument (argc, argv, "-grid-bin-gpc2"))) {
+    remove_argument (N, &argc, argv);
+    GRID_BIN_GPC2 = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  GRID_BIN_HSC = 8;
+  if ((N = get_argument (argc, argv, "-grid-bin-hsc"))) {
+    remove_argument (N, &argc, argv);
+    GRID_BIN_HSC = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  GRID_BIN_CFH = 8;
+  if ((N = get_argument (argc, argv, "-grid-bin-cfh"))) {
+    remove_argument (N, &argc, argv);
+    GRID_BIN_CFH = atof (argv[N]);
+    remove_argument (N, &argc, argv);
   }
 
@@ -555 +590 @@
     case SYNTH_PHOT:
     case UPDATE_AVERAGES:
+      // note: initialize sets PhotcodeList to the full set of average photcodes
+      // if the mode is UPDATE_AVERAGES
       if (argc != 1) relphot_usage();
       break;
@@ -608 +645 @@
   FREE (REGION_FILE);
   FREE (SYNTH_ZERO_POINTS);
+  FREE (GRID_MEANFILE);
   FREE (IMAGE_TABLE);
   FREE (CATDIR);
@@ -763 +801 @@
     RESET = TRUE;
   }
-
   RESET_ZEROPTS = FALSE;
   if ((N = get_argument (argc, argv, "-reset-zpts"))) {
     remove_argument (N, &argc, argv);
     RESET_ZEROPTS = TRUE;
+  }
+  RESET_FLATCORR = FALSE;
+  if ((N = get_argument (argc, argv, "-reset-flat"))) {
+    remove_argument (N, &argc, argv);
+    RESET_FLATCORR = TRUE;
   }
 
@@ -864 +906 @@
   }
 
+  // load the GridCorrection file
+  GRID_ZEROPT = FALSE;
+  GRID_MEANFILE = NULL;
+  if ((N = get_argument (argc, argv, "-grid"))) {
+    GRID_ZEROPT = TRUE;
+    remove_argument (N, &argc, argv);
+    GRID_MEANFILE = strcreate (argv[N]);
+    remove_argument (N, &argc, argv);
+    GRID_ZPT_MODE  = GRID_ZPT_MODE_ALL; // since we are applying the grid, need to active this
+  }
+
   ImagSelect = FALSE;
   if ((N = get_argument (argc, argv, "-instmag"))) {
@@ -922 +975 @@
 
   FREE (SYNTH_ZERO_POINTS);
+  FREE (GRID_MEANFILE);
   FREE (BOUNDARY_TREE);
   FREE (UPDATE_CATFORMAT);

trunk/Ohana/src/relphot/src/assign_images.c

@@ -36 +36 @@
   if (MOSAIC_ZEROPT) {
     makeMosaics (image, Nimage, FALSE);
-
-    // center coords and Mcal, dMcal, Mchisq for the mosaics
+    MARKTIME("set mosaic coordinates and Mcal values: %f sec\n", dtime);
+
+    // center coords and zero Mcal, dMcal, Mchisq for the mosaics
     setMosaicCenters (image, Nimage);
     MARKTIME("set mosaic coordinates and Mcal values: %f sec\n", dtime);
@@ -124 +125 @@
       Mosaic *mosaic = getMosaicForImage (j); 
       if (mosaic) {
-        Rc = mosaic->coords.crval1;
-        Dc = mosaic->coords.crval2;
-        // NOTE : have defined mosaic Rc,Dc to choose the side of 0,360 on which most of the
-        // chips are located.  but, for host assignment, we rationalize to 0.0 - 360.0
-        Rc = ohana_normalize_angle_to_midpoint (Rc, 180.0);
+        Rc = mosaic->coords.crval1;
+        Dc = mosaic->coords.crval2;
+        // NOTE : have defined mosaic Rc,Dc to choose the side of 0,360 on which most of the
+        // chips are located.  but, for host assignment, we rationalize to 0.0 - 360.0
+        Rc = ohana_normalize_angle_to_midpoint (Rc, 180.0);
       }
     }

trunk/Ohana/src/relphot/src/bcatalog.c

@@ -35 +35 @@
   ALLOCATE (Nvalid, int, Nsecfilt);
 
-  // flags used by the photometry analysis (excluding UBERCAL)
-  unsigned int PHOTOM_FLAGS = 
-    ID_MEAS_NOCAL          | // detection ignored for this analysis (photcode, time range)
-    ID_MEAS_POOR_PHOTOM    | // detection is photometry outlier 
-    ID_MEAS_SKIP_PHOTOM    | // detection was ignored for photometry measurement 
-    ID_MEAS_AREA           | // detetion was outside acceptable area of device
-    ID_MEAS_SYNTH_MAG      | // magnitude is synthetic
-    ID_MEAS_STACK_PRIMARY  | // this stack measurement is in the primary skycell
-    ID_MEAS_STACK_PHOT_SRC | // this measurement supplied the stack photometry
-    ID_MEAS_PHOTOM_PSF     | // this measurement is used for the mean psf mag
-    ID_MEAS_PHOTOM_APER    | // this measurement is used for the mean ap mag
-    ID_MEAS_PHOTOM_KRON    ; // this measurement is used for the mean kron mag
-    
   /* exclude stars not in range or with too few measurements */
   for (i = 0; i < catalog[0].Naverage; i++) {
@@ -61 +48 @@
     }
 
-    // reset the calculated average magnitudes (does not affect ubercal-tied measurements or images)
-    if (RESET) {
-      int Ns;
-
-      for (Ns = 0; Ns < Nphotcodes; Ns++) {
-
-        int thisCode = photcodes[Ns][0].code;
-        int Nsec = GetPhotcodeNsec(thisCode);
-
-        dvo_secfilt_init (&subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec], SECFILT_RESET_ALL);
-      }
-    }
+    // reset the calculated average magnitudes for active photcode only
+    ResetAverageActivePhotcodes (&subcatalog[0].secfilt[Nsecfilt*Naverage]);
 
     Nm = 0;
@@ -111 +88 @@
       if (isnan(catalog[0].measure[offset].Mkron)) { Nnan ++; continue; }
 
-      // require 0x01 in photFlags (fitted with a PSF)
-      if ((catalog[0].measure[offset].photFlags & 0x01) == 0) { Nbad ++; continue; }
+      // require 0x01 in photFlags (fitted with a PSF) -- add to the test data
+      if (REQUIRE_PSFFIT && (catalog[0].measure[offset].photFlags & 0x01) == 0) { Nbad ++; continue; }
 
       // very loose cut on PSF - Kron 
@@ -156 +133 @@
       subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_SKIP_PHOTOM;
       subcatalog[0].measureT[Nmeasure].averef = Naverage;
-      if (RESET) { 
-        // only reset Mcal for measures with a matching image
-        // do not reset Mcal for ubercal images unless explicitly requested
-        if (subcatalog[0].measureT[Nmeasure].dbFlags & ID_MEAS_PHOTOM_UBERCAL) {
-          if (!KEEP_UBERCAL) {
-            subcatalog[0].measureT[Nmeasure].McalPSF  = 0.0;
-            subcatalog[0].measureT[Nmeasure].McalAPER = 0.0;
-            subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_PHOTOM_UBERCAL;
-          } 
-        } else {
-          if (getImageEntry (Nmeasure, Ncat) >= 0) {
-            subcatalog[0].measureT[Nmeasure].McalPSF  = 0.0;
-            subcatalog[0].measureT[Nmeasure].McalAPER = 0.0;
-          }
-        }
-        subcatalog[0].measureT[Nmeasure].dbFlags &= ~PHOTOM_FLAGS;
-      }
+      ResetMeasureZeroPoints (&subcatalog[0].measureT[Nmeasure], Nmeasure, Ncat);
+
       Nmeasure ++;
       Nm ++;

trunk/Ohana/src/relphot/src/extra.c

@@ -1 +1 @@
 # include "relphot.h"
+# define DISABLE_CORNER_GPC1 0
 
 int isMosaicChip (int photcode) {
@@ -14 +15 @@
 int whichGPC1filter (int photcode) {
 
-# if (1)
+# if (DISABLE_CORNER_GPC1)
   // disable the corner chips:
   if (((photcode >  10000) && (photcode <  10077)) || (photcode == 4100)) return PS1_g; // g-band 
@@ -37 +38 @@
 int isGPC1chip (int photcode) {
 
-# if (1)
+# if (DISABLE_CORNER_GPC1)
   if (((photcode >  10000) && (photcode <  10077)) || (photcode == 4100)) return TRUE; // g-band
   if (((photcode >  10100) && (photcode <  10177)) || (photcode == 4200)) return TRUE; // r-band

trunk/Ohana/src/relphot/src/global_stats.c

@@ -1 +1 @@
 # include "relphot.h"
 
-void global_stats (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr, int nloop) {
+void global_stats (Catalog *catalog, int Ncatalog, int nloop) {
 
   StatType stN, stX, stS, imN, imX, imM, imD, msM, msX, msN, msD, tgM, tgX, tgD;
@@ -33 +33 @@
     int seccode = photcodes[Ns][0].code;
 
-    stN = statsStarN (catalog, Ncatalog, Nsec, seccode, flatcorr);
+    stN = statsStarN (catalog, Ncatalog, Nsec, seccode);
     stX = statsStarX (catalog, Ncatalog, Nsec);
     stS = statsStarS (catalog, Ncatalog, Nsec);
@@ -43 +43 @@
   }
   
-  imN = statsImageN (catalog, flatcorr);
+  imN = statsImageN (catalog);
   imX = statsImageX (catalog);
   imM = statsImageM (catalog);
   imD = statsImagedM (catalog);
   
-  msN = statsMosaicN (catalog, flatcorr);
+  msN = statsMosaicN (catalog);
   msM = statsMosaicM (catalog);
   msD = statsMosaicdM (catalog);
@@ -58 +58 @@
   
   fprintf (stderr, "meas / image:   %7.0f  %7.0f  %7.0f  %7.0f  %7.0f  %6d\n",   imN.median, imN.mean, imN.sigma, imN.min, imN.max, imN.Nmeas);
-  fprintf (stderr, "Mcal image:     %7.4f  %7.4f  %7.4f  %7.4f  %7.4f  %6d\n",   imM.median, imM.mean, imM.sigma, imM.min, imM.max, imM.Nmeas);
+  fprintf (stderr, "Mcal  image:    %7.4f  %7.4f  %7.4f  %7.4f  %7.4f  %6d\n",   imM.median, imM.mean, imM.sigma, imM.min, imM.max, imM.Nmeas);
   fprintf (stderr, "dMcal image:    %7.4f  %7.4f  %7.4f  %7.4f  %7.4f  %6d\n",   imD.median, imD.mean, imD.sigma, imD.min, imD.max, imD.Nmeas);
   fprintf (stderr, "chisq image:    %7.1f  %7.1f  %7.1f  %7.1f  %7.1f  %6d\n",   imX.median, imX.mean, imX.sigma, imX.min, imX.max, imX.Nmeas);
 
   fprintf (stderr, "meas / mosaic:  %7.0f  %7.0f  %7.0f  %7.0f  %7.0f  %6d\n",   msN.median, msN.mean, msN.sigma, msN.min, msN.max, msN.Nmeas);
-  fprintf (stderr, "Mcal mosaic:    %7.4f  %7.4f  %7.4f  %7.4f  %7.4f  %6d\n",   msM.median, msM.mean, msM.sigma, msM.min, msM.max, msM.Nmeas);
+  fprintf (stderr, "Mcal  mosaic:   %7.4f  %7.4f  %7.4f  %7.4f  %7.4f  %6d\n",   msM.median, msM.mean, msM.sigma, msM.min, msM.max, msM.Nmeas);
   fprintf (stderr, "dMcal mosaic:   %7.4f  %7.4f  %7.4f  %7.4f  %7.4f  %6d\n",   msD.median, msD.mean, msD.sigma, msD.min, msD.max, msD.Nmeas);
   fprintf (stderr, "chisq mosaic:   %7.1f  %7.1f  %7.1f  %7.1f  %7.1f  %6d\n",   msX.median, msX.mean, msX.sigma, msX.min, msX.max, msX.Nmeas);

trunk/Ohana/src/relphot/src/initialize.c

@@ -27 +27 @@
   }    
     
-  if (USE_GRID && (Nphotcodes > 1)) {
-    fprintf (stderr, "grid correction analysis currently can only operate on a single photcode\n");
-    exit (1);
-  }
-
-  IMAGE_BAD = ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP;
-  STAR_BAD  = ID_OBJ_POOR | ID_OBJ_FEW;
-  MEAS_BAD  = ID_MEAS_NOCAL | ID_MEAS_POOR_PHOTOM | ID_MEAS_SKIP_PHOTOM | ID_MEAS_AREA;
-
   if (SHOW_PARAMS) {
     int Ns;
@@ -67 +58 @@
 
     fprintf (stderr, "VERBOSE: %d, PLOTSTUFF: %d\n", VERBOSE, PLOTSTUFF);
-    fprintf (stderr, "GRID_X: %d, GRID_Y: %d, BINNING: %d == Nmx: %d, Nmy: %d\n", 
-             RELPHOT_GRID_X, 
-             RELPHOT_GRID_Y, 
-             RELPHOT_GRID_BINNING, 
-             (int)(RELPHOT_GRID_X/RELPHOT_GRID_BINNING), (int)(RELPHOT_GRID_Y/RELPHOT_GRID_BINNING));
 
     fprintf (stderr, "STAR_SCATTER           %lf\n", STAR_SCATTER);
+    fprintf (stderr, "STAR_CHISQ             %lf\n", STAR_CHISQ);
+
+    fprintf (stderr, "MOSAIC_SCATTER         %lf\n", MOSAIC_SCATTER);
+    fprintf (stderr, "MOSAIC_CHISQ           %lf\n", MOSAIC_CHISQ);
+
+    fprintf (stderr, "NIGHT_SCATTER          %lf\n", NIGHT_SCATTER);
+    fprintf (stderr, "NIGHT_CHISQ            %lf\n", NIGHT_CHISQ);
+
     fprintf (stderr, "IMAGE_SCATTER          %lf\n", IMAGE_SCATTER);
     fprintf (stderr, "IMAGE_OFFSET           %lf\n", IMAGE_OFFSET);

trunk/Ohana/src/relphot/src/launch_region_hosts.c

@@ -62 +62 @@
     if (RESET)                      strextend (&command, "-reset");
     if (RESET_ZEROPTS)              strextend (&command, "-reset-zpts");
+    if (RESET_FLATCORR)             strextend (&command, "-reset-flat");
     if (!KEEP_UBERCAL)              strextend (&command, "-reset-ubercal");
     if (DophotSelect)               strextend (&command, "-dophot %d", DophotValue);

trunk/Ohana/src/relphot/src/load_catalogs.c

@@ -165 +165 @@
     if (RESET)               { strextend (&command, "-reset"); }
     if (RESET_ZEROPTS)       { strextend (&command, "-reset-zpts"); }
+    if (RESET_FLATCORR)      { strextend (&command, "-reset-flat"); }
     if (!KEEP_UBERCAL)       { strextend (&command, "-reset-ubercal"); }
     if (DophotSelect)        { strextend (&command, "-dophot %d", DophotValue); }

trunk/Ohana/src/relphot/src/load_images.c

@@ -57 +57 @@
   // reset image values as needed.  always allow 'few' images to succeed, if possible (new
   // images / detections may have been added
-  for (off_t i = 0; i < Nsubset; i++) {
-    // reset these two regardless (we will re-determine)
-    subset[i].flags &= ~(ID_IMAGE_PHOTOM_NOCAL | ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_SKIP | ID_IMAGE_PHOTOM_FEW);
-    if (RESET) {
-      if (RESET_ZEROPTS) {
-        if (!KEEP_UBERCAL || !(subset[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) {
-          subset[i].McalPSF  = 0.0;
-          subset[i].McalAPER = 0.0;
-          subset[i].dMcal    = NAN;
-          subset[i].flags   &= ~ID_IMAGE_PHOTOM_UBERCAL;
-        }
-      }
-      subset[i].flags &= ~(ID_IMAGE_MOSAIC_POOR | ID_IMAGE_NIGHT_POOR);
-      subset[i].flags &= ~(ID_IMAGE_IMAGE_PHOTCAL | ID_IMAGE_MOSAIC_PHOTCAL | ID_IMAGE_TGROUP_PHOTCAL);
-      subset[i].ubercalDist = 1000;
-    }
-  }
+  ResetImages (subset, Nsubset);
 
   // XXX consider allowing relphot to skip the AstroMap load (do we need precise astrometry here?)
@@ -115 +99 @@
   return TRUE;
 }
+
+// This function re-creates the image subset array using the array of LineNumbers saved
+// initially by load_images.  This function is only called by relphot_images.c after the
+// loops have been completed (if any).  I am not certain this function is doing anything
+// at all : the values in the image subset are copied back to the master table in
+// save_images_updates() so the values copied from the master to the subset should already
+// be there.
 
 int reload_images (FITS_DB *db) {

trunk/Ohana/src/relphot/src/plot_scatter.c

@@ -1 +1 @@
 # include "relphot.h"
    
-void plot_scatter (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) {
+void plot_scatter (Catalog *catalog, int Ncatalog) {
 
   off_t i, j, k, m, N, Ntot;
@@ -28 +28 @@
     N = 0;
     for (i = 0; i < Ncatalog; i++) {
-        for (j = 0; j < catalog[i].Naverage; j++) {
+      for (j = 0; j < catalog[i].Naverage; j++) {
 
-            /* calculate the average value for a single star */
-            if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;  
-            m = catalog[i].averageT[j].measureOffset;
+        /* calculate the average value for a single star */
+        // if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;  
+        m = catalog[i].averageT[j].measureOffset;
 
-            for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) {
-                // skip measurements that do not match the current photcode
-                int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode);
-                if (ecode != thisCode) { continue; }
+        for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) {
+          // skip measurements that do not match the current photcode
+          int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode);
+          if (ecode != thisCode) { continue; }
 
-                if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue;
-                Mcal = getMcal  (m, i, MAG_CLASS_PSF);
-                if (isnan(Mcal)) continue;
-                Mmos = getMmos  (m, i);
-                if (isnan(Mmos)) continue;
-                Mgrid = getMgrid (m, i);
-                if (isnan(Mgrid)) continue;
+          if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue;
+          Mcal = getMcal  (m, i, MAG_CLASS_PSF);
+          if (isnan(Mcal)) continue;
+          Mmos = getMmos  (m, i);
+          if (isnan(Mmos)) continue;
+          Mgrid = getMgridTiny (&catalog[i].measureT[m]);
+          if (isnan(Mgrid)) continue;
 
-                Mrel = catalog[i].secfilt[Nsecfilt*j+Nsec].MpsfChp;
-                if (isnan(Mrel)) continue;
+          Mrel = catalog[i].secfilt[Nsecfilt*j+Nsec].MpsfChp;
+          if (isnan(Mrel)) continue;
 
-                xlist[N] = Mrel;
-                ylist[N] = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt], MAG_CLASS_PSF) - Mcal - Mmos - Mgrid - Mrel;
-                ilist[N] = PhotInstTiny (&catalog[i].measureT[m], MAG_CLASS_PSF);
-                N++;
-            }
+          xlist[N] = Mrel;
+          ylist[N] = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt], MAG_CLASS_PSF) - Mcal - Mmos - Mgrid - Mrel;
+          ilist[N] = PhotInstTiny (&catalog[i].measureT[m], MAG_CLASS_PSF);
+          N++;
         }
+      }
     }
 

trunk/Ohana/src/relphot/src/reload_catalogs.c

@@ -9 +9 @@
     gettimeofday (&start, (void *) NULL);
 
-void reload_catalogs (SkyList *skylist, FlatCorrectionTable *flatcorr, int hostID, char *hostpath) {
+void reload_catalogs (SkyList *skylist, int hostID, char *hostpath) {
 
   int i;
@@ -17 +17 @@
 
   struct timeval start, stop;
-  double dtime = 0.0;
-  double time1 = 0.0;
-  double time2 = 0.0;
-  double time3a = 0.0;
-  double time3b = 0.0;
-  double time4 = 0.0;
-  double time5 = 0.0;
-  double time6 = 0.0;
-  double time7 = 0.0;
+  double dtime = 0.0, time1 = 0.0, time2 = 0.0, time3 = 0.0;
+  double time4 = 0.0, time5 = 0.0, time6 = 0.0, time7 = 0.0, time8 = 0.0;
 
   // XXX need to decide how to determine PARALLEL mode...
@@ -81 +74 @@
 
     populate_tiny_values(&catalog, DVO_TV_MEASURE | DVO_TV_AVERAGE);
-    TIMESTAMP(time3a);
+    TIMESTAMP(time3);
 
     // XXX need to worry about the image subset data
     initImageBins  (&catalog, 1, FALSE);
-    initMosaicBins (&catalog, 1, FALSE);
-    initTGroupBins (&catalog, 1);
-
-    initGridBins   (&catalog, 1); 
-    TIMESTAMP(time3b);
+    TIMESTAMP(time4);
 
     findImages (&catalog, 1, FALSE);
-    findMosaics (&catalog, 1, FALSE);
-    findTGroups (&catalog, 1);
-
-    TIMESTAMP(time4);
+    TIMESTAMP(time5);
 
     initMrel (&catalog, 1);
-    setMrelFinal (&catalog, flatcorr, FALSE);
-    TIMESTAMP(time5);
+
+    setMrelFinal (&catalog, FALSE);
+    TIMESTAMP(time6);
 
     // modify the output format as desired (ignore current format on disk)
@@ -119 +106 @@
     free_tiny_values(&catalog);
     dvo_catalog_free (&catalog);
-    TIMESTAMP(time6);
-
-    // XXX freeTgroupBins here: 
+    TIMESTAMP(time7);
 
     freeImageBins (1, FALSE);
-    freeMosaicBins (1, FALSE);
-    freeTGroupBins (1);
-    freeGridBins (1);
-    TIMESTAMP(time7);
+    TIMESTAMP(time8);
 
     if (hostID) {
@@ -134 +116 @@
   }
 
-  fprintf (stderr, "time step 1  %10.3f sec : find catalog\n", time1);
-  fprintf (stderr, "time step 2  %10.3f sec : load catalog\n", time2);
-  fprintf (stderr, "time step 3  %10.3f sec : make tiny values\n", time3a);
-  fprintf (stderr, "time step 4  %10.3f sec : init imbins\n",  time3b);
-  fprintf (stderr, "time step 5  %10.3f sec : find images\n",  time4);
-  fprintf (stderr, "time step 6  %10.3f sec : set Mrel\n",     time5);
-  fprintf (stderr, "time step 7  %10.3f sec : save catalog\n", time6);
-  fprintf (stderr, "time step 8  %10.3f sec : free catalog\n", time7);
+  fprintf (stderr, "time step 1  %10.3f sec : find catalog\n",     time1);
+  fprintf (stderr, "time step 2  %10.3f sec : load catalog\n",     time2);
+  fprintf (stderr, "time step 3  %10.3f sec : make tiny values\n", time3);
+  fprintf (stderr, "time step 4  %10.3f sec : init imbins\n",      time4);
+  fprintf (stderr, "time step 5  %10.3f sec : find images\n",      time5);
+  fprintf (stderr, "time step 6  %10.3f sec : set Mrel\n",         time6);
+  fprintf (stderr, "time step 7  %10.3f sec : save catalog\n",     time7);
+  fprintf (stderr, "time step 8  %10.3f sec : free catalog\n",     time8);
 }
 
@@ -231 +213 @@
     char *command = NULL;
     strextend (&command, "relphot_client %s -update-catalogs %s", PhotcodeList, imageFile);
-    strextend (&command, "-hostID %d -D CATDIR %s -hostdir %s", group->hosts[i][0].hostID, CATDIR, group->hosts[i][0].pathname);
+    strextend (&command, "-hostID %d", group->hosts[i][0].hostID);
+    strextend (&command, "-D CATDIR %s", CATDIR);
+    strextend (&command, "-hostdir %s", group->hosts[i][0].pathname);
     strextend (&command, "-region %f %f %f %f", UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax);
     strextend (&command, "-statmode %s", STATMODE);
     strextend (&command, "-D CAMERA %s", CAMERA);
-    strextend (&command, "-D STAR_TOOFEW %d -minerror %f", STAR_TOOFEW, MIN_ERROR);
-
-    // options & configs which affect relphot_client -update-catalogs
-    // VERBOSE, VERBOSE2
-    // RESET
-    // RESET_ZEROPTS
-    // TimeSelect
-    // AreaSelect
-    // STATMODE
-    // STAR_TOOFEW
-    // MIN_ERROR
-    
+    strextend (&command, "-D STAR_TOOFEW %d", STAR_TOOFEW);
+    strextend (&command, "-minerror %f", MIN_ERROR);
+
     if (VERBOSE)           { strextend (&command, "-v"); }
     if (VERBOSE2)          { strextend (&command, "-vv"); }
     if (RESET)             { strextend (&command, "-reset"); }
     if (RESET_ZEROPTS)     { strextend (&command, "-reset-zpts"); }
+    if (RESET_FLATCORR)    { strextend (&command, "-reset-flat"); }
+    if (!KEEP_UBERCAL)     { strextend (&command, "-reset-ubercal"); }
     if (PRESERVE_PS1)      { strextend (&command, "-preserve-ps1"); }
+    if (IS_DIFF_DB)        { strextend (&command, "-is-diff-db"); }
     if (UPDATE)            { strextend (&command, "-update"); }
-    if (IS_DIFF_DB)        { strextend (&command, "-is-diff-db"); }
-    if (!KEEP_UBERCAL)     { strextend (&command, "-reset-ubercal"); }
     if (UPDATE_CATFORMAT)  { strextend (&command, "-update-catformat %s", UPDATE_CATFORMAT); }
     if (BOUNDARY_TREE)     { strextend (&command, "-boundary-tree %s", BOUNDARY_TREE); }
     if (SYNTH_ZERO_POINTS) { strextend (&command, "-synthphot-zpts %s", SYNTH_ZERO_POINTS); }
+    if (GRID_ZEROPT)       { strextend (&command, "-grid %s", GRID_MEANFILE); } 
     if (USE_BASIC_CHECK)   {  strextend (&command, "-basic-image-search"); }
+    // if (USE_ALL_IMAGES)    { strextend (&command, "-use-all-images"); }
     if (USE_MCAL_PSF_FOR_STACK_APER) { strextend (&command, "-use-mcal-psf-for-stack-aper"); }
 

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

@@ -34 +34 @@
   switch (mode) {
     case UPDATE_IMAGES:
+      // calculate zero points for images (may group by exposure [mosaic], night [tgroups]; may measure flat-correction)
       // IF CALLED WITH NLOOP == 0, DOES NOT LOAD bcatalog, just loads image table and generates ImageSubset.dat
+      // If called with -update, calls reload_catalogs() just like apply_offsets
       relphot_images (skylist);
       relphot_free (sky, skylist);
@@ -52 +54 @@
 
     case PARALLEL_REGIONS:
-      // run image updates in parallel across multiple remote machines
+      // calculate zero points for images (may group my exposure [mosaic], night [tgroups]; may measure flat-correction)
+      // equivalent to relphot_images, but run in parallel across multiple remote machines
       relphot_parallel_regions (sky);
       relphot_free (sky, skylist);
@@ -69 +72 @@
         exit (2);
       }
-      reload_catalogs (skylist, NULL, 0, NULL);
+      // we do not need to load grid corrections here.  the value of
+      // GRID_MEANFILE, if specified, is passed to the clients
+      reload_catalogs (skylist, 0, NULL);
       relphot_free (sky, skylist);
       exit (0);
@@ -84 +89 @@
 
     I would like to merge the functionality of relphot_images (Nloop == 0), relphot_objects, and apply_offsets
+
+    relphot_objects vs reload_objects (apply offsets):
+
+    relphot_objects 
+    
     
 ***/

trunk/Ohana/src/relphot/src/relphot_client.c

@@ -80 +80 @@
       initImagesSubset (image, NULL, Nimage);
 
-      // load the flat-field correction table from CATDIR
-      FlatCorrectionTable *flatcorr = NULL;
-      reload_catalogs (skylist, flatcorr, HOST_ID, HOSTDIR);
+      // load grid corrections here
+      GridCorrectionLoad (GRID_MEANFILE);
+
+      reload_catalogs (skylist, HOST_ID, HOSTDIR);
       freeImages ((char *)image);
       free (image);
+      freeGridBins ();
       client_logger_message ("updated catalogs\n");
       relphot_client_free (sky, skylist);
@@ -90 +92 @@
     }
 
+/* deprecated (absorbed into MODE_UPDATE)
     case MODE_UPDATE_OBJECTS: {
       // take the current set of detections and set the mean magnitudes
@@ -97 +100 @@
       break;
     }
+*/
 
     case MODE_SYNTH_PHOT:
@@ -104 +108 @@
       break;
 
+    case MODE_UPDATE_OBJECTS:
     default:
       fprintf (stderr, "impossible!");

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

@@ -38 +38 @@
   MARKTIME("-- load image data: %f sec\n", dtime);
 
-  /* load regions and images based on specified sky patch and/or catalog */
+  // load regions and images based on specified sky patch and/or catalog 
+  // NOTE: load_images transfers zero points from images to mosaics (initMosaicMcal)
+  // and to tgroups (initTGroupMcal), resetting the image zero points to 0.0.
+
+  // (if both tgroups and mosaics are used, the mosaics will hold the value and tgroups
+  // will get values of 0.0.  This is not clearly the right way to go, but may not matter
+  // if we reset everything).
+
   load_images (&db, skylist, &UserPatch, FALSE, USE_ALL_IMAGES);
   MARKTIME("-- load images: %f sec\n", dtime);
@@ -45 +52 @@
   if (!UPDATE) dvo_image_unlock (&db); 
 
-  // load the flat correction table (if defined)
-  // we no longer need to do this: the flatcorrection is now saved in measure.Mflat
-# if (0)
-  char flatcorrfile[256];
-  sprintf (flatcorrfile, "%s/flatcorr.fits", CATDIR);
-  FlatCorrectionTable *flatcorr = FlatCorrectionLoad (flatcorrfile, VERBOSE);
-# endif
-  FlatCorrectionTable *flatcorr = NULL;
+  initGridBins (); // allocates the empty array of corrections (elements are built below)
 
   if (CALIBRATE_STACKS_AND_WARPS || (NLOOP > 0)) {
+    initGrid (); // allocate grid correction entries for existing photcodes
+
     /* load catalog data from region files (hostID is 0 since we are not a client */
     catalog = load_catalogs (skylist, &Ncatalog, 0, NULL, NULL);
@@ -68 +70 @@
     initTGroupBins (catalog, Ncatalog);
 
-    initGridBins   (catalog, Ncatalog);
     initMrel (catalog, Ncatalog);
 
@@ -74 +75 @@
     MARKTIME("-- set up image indexes: %f sec\n", dtime);
 
-    findMosaics (catalog, Ncatalog, TRUE);  /* also sets Grid values */
+    findMosaics (catalog, Ncatalog, TRUE);
     MARKTIME("-- set up mosaic indexes: %f sec\n", dtime);
 
@@ -82 +83 @@
     SAVEPLOT = FALSE;
 
+    // mark ID_MEAS_AREA (user-selected chip region) and ID_MEAS_NOCAL (time range & active photcode)
     setExclusions (catalog, Ncatalog, TRUE);
 
-    global_stats (catalog, Ncatalog, flatcorr, 0);
+    global_stats (catalog, Ncatalog, 0);
 
     if (PLOTSTUFF) {
@@ -91 +93 @@
     }
 
-    // if we are measuring the flat-field correction grid, we need to perform a number of iterations first:
-    if (USE_GRID) {
-      int star_toofew;
-
-      // until we finish the grid analysis, do not reject stars out-of-hand based on ID_OBJ_FEW
-      // XXX this is kind of poor: need to have a better distinctions about STAR_BAD in setMrel vs getMrel
-      star_toofew = STAR_TOOFEW;
-      STAR_TOOFEW = 0;
-      for (i = 0; i < NGRID; i++) {
-        STAR_BAD = ID_OBJ_POOR;
-        setMrel  (catalog, Ncatalog, flatcorr);
-        STAR_BAD  = ID_OBJ_POOR | ID_OBJ_FEW;
-        setMgrid (catalog, flatcorr);
-      }
-      STAR_BAD  = ID_OBJ_POOR | ID_OBJ_FEW;
-      STAR_TOOFEW = star_toofew;
-    }
-
     /* determine fit values */
     for (i = 0; i < NLOOP; i++) {
       SetZptIteration (i);
 
-      setMrel  (catalog, Ncatalog, flatcorr); // threaded
-
-      setMcal  (catalog, flatcorr);
-      setMmos  (catalog, flatcorr);
-      setMgrp  (catalog, flatcorr);
-
-      setMgrid (catalog, flatcorr);
+      setMrel  (catalog, Ncatalog); // threaded (calls setMrelCatalog)
+
+      setMcal  (catalog);
+      setMmos  (catalog);
+      setMgrp  (catalog);
+
+      setMgrid (catalog, Ncatalog);
       MARKTIME("-- set Mrel, Mcal, Mmos, Mgrid : %f sec\n", dtime);
     
       if (PLOTSTUFF) {
-        plot_scatter (catalog, Ncatalog, flatcorr); 
-        plot_grid (catalog, flatcorr); 
+        plot_scatter (catalog, Ncatalog); 
         plot_mosaics ();
         plot_images ();
@@ -131 +114 @@
       }
 
-      global_stats (catalog, Ncatalog, flatcorr, i);
-      SetZeroPointModes ();
+      global_stats (catalog, Ncatalog, i);
+      SetZeroPointModes (catalog, Ncatalog);
 
       MARKTIME("-- finished loop %d: %f sec\n", i, dtime);
@@ -139 +122 @@
 
     if (PLOTSTUFF) {
-      plot_scatter (catalog, Ncatalog, flatcorr); 
-      plot_grid (catalog, flatcorr); 
+      plot_scatter (catalog, Ncatalog); 
       plot_mosaics ();
       plot_images ();
@@ -147 +129 @@
     }
   
-    if (USE_GRID) dump_grid ();
+    // if (GRID_ZEROPT) dump_grid ();
 
     MARKTIME("-- finalize Mcal values: %f sec\n", dtime);
@@ -169 +151 @@
     freeMosaicBins (Ncatalog, TRUE);
     freeTGroupBins (Ncatalog);
-    freeGridBins (Ncatalog);
+
+    GridCorrectionSave ();
 
     // end of if (NLOOP > 0) block : this loop determines the offsets per chip 
@@ -179 +162 @@
   }
 
+  // once the zero points have been calculated, they are reassigned to the
+  // individual images.  At this point, the mosaic and tgroup values should not 
+  // be used.
+  freeMosaics ();
+  freeTGroups ();
+  MOSAIC_ZEROPT = FALSE;
+  TGROUP_ZEROPT = FALSE;
+
+  // only change the real database files if -update is requested
+  if (!UPDATE) {
+    dvo_image_unlock (&db); 
+    freeImages (db.ftable.buffer);
+    gfits_db_free (&db);
+    freeGridBins ();
+    return TRUE;
+  }
+  
   reload_images (&db);
 
-  // only change the real database files if -update is requested
-  if (!UPDATE) exit (0);
-  
   /* Load catalog data from region files, update Mrel include all data.  In a parallel
      context, this function writes the image parameters as a subset table for the remote
      clients */
   if (!CALIBRATE_STACKS_AND_WARPS) {
-    reload_catalogs (skylist, flatcorr, 0, NULL);
+    reload_catalogs (skylist, 0, NULL); // calls setMrelFinal which setMrelOutput which calls setMrelCatalog
     MARKTIME("-- updated all catalogs: %f sec\n", dtime);
   }
 
   if (UPDATE_CATFORMAT) {
-    // ensure the db format is updated
+    // ensure the Image db format is updated
     db.format = dvo_catalog_catformat (UPDATE_CATFORMAT);
     gfits_modify (&db.header, "FORMAT", "%s", 1, UPDATE_CATFORMAT);
@@ -204 +201 @@
 
   freeImages (db.ftable.buffer);
-  freeMosaics ();
+  freeGridBins ();
+
   gfits_db_free (&db);
 

trunk/Ohana/src/relphot/src/relphot_objects.c

@@ -6 +6 @@
 // has been correctly propagated to the measurement, the average will respect that value.
 
-int relphot_objects_parallel (SkyList *sky);
-
 int relphot_objects (SkyList *skylist, int hostID, char *hostpath) {
 
-  off_t i, j, k;
-  int Nsecfilt;
-  struct stat filestat;
-
-  Catalog catalog;
+  FITS_DB db;
 
   INITTIME;
 
-  // XXX need to decide how to determine PARALLEL mode...
-  if (PARALLEL && !hostID) {
-    relphot_objects_parallel (skylist);
-    return TRUE;
-  }
-
-  FITS_DB db;
-      
   set_db (&db);
   gfits_db_init (&db);
@@ -41 +27 @@
   MARKTIME("loaded images: %f sec\n", dtime);
 
-  // load the ZP corrections here
-  if (SYNTH_ZERO_POINTS) SynthZeroPointsLoad (SYNTH_ZERO_POINTS);
+  // load grid corrections here (specified by -grid-meanfile)
+  GridCorrectionLoad (GRID_MEANFILE);
 
-  if (REPAIR_WARPS) { 
-    FindWarpGroups ();
-    MARKTIME("setup warp groups: %f sec\n", dtime);
-    MakeStackIndex ();
-    MARKTIME("setup stack index: %f sec\n", dtime);
-    RepairWarpMeasuresOpenLogfile ();
-  }
-
-  // load data from each region file, only use bright stars
-  for (i = 0; i < skylist[0].Nregions; i++) {
-
-    // does this host ID match the desired location for the table?
-    if (!HostTableTestHost(skylist[0].regions[i], hostID)) continue;
-
-    // set up the basic catalog info
-    char hostfile[1024];
-    snprintf (hostfile, 1024, "%s/%s.cpt", hostpath, skylist[0].regions[i]->name);
-
-    dvo_catalog_init (&catalog, TRUE);
-    catalog.filename    = hostID ? hostfile : skylist[0].filename[i];
-
-    // only update existing db tables
-    int status = stat (catalog.filename, &filestat);
-    if ((status == -1) && (errno == ENOENT)) {
-      if (VERBOSE) fprintf (stderr, "no file %s, skipping\n", catalog.filename);
-      continue;
-    }
-
-    catalog.catflags    = DVO_LOAD_AVERAGE | DVO_LOAD_MEASURE | DVO_LOAD_MISSING | DVO_LOAD_SECFILT;
-    // if we are going to update the format, we should update all tables
-    if (UPDATE_CATFORMAT) {
-      catalog.catflags |= DVO_LOAD_LENSING | DVO_LOAD_LENSOBJ | DVO_LOAD_STARPAR | DVO_LOAD_GALPHOT;
-    }
-    catalog.Nsecfilt    = GetPhotcodeNsecfilt ();
-
-    if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE, "w")) {
-      fprintf (stderr, "ERROR: failure reading catalog %s, skipping\n", catalog.filename);
-      continue;
-    }
-    if (!catalog.Naverage_disk) {
-      if (VERBOSE) fprintf (stderr, "no data in %s, skipping\n", catalog.filename);
-      dvo_catalog_unlock (&catalog);
-      dvo_catalog_free (&catalog);
-      continue;
-    }
-
-    // reset the calculated average magnitudes (does not affect ubercal-tied measurements or images)
-    if (RESET) {
-      Nsecfilt = catalog.Nsecfilt;
-
-      DVOAverageFlags photomBits = 
-        ID_OBJ_EXT             | // extended in our data (eg, PS)
-        ID_OBJ_EXT_ALT         | // extended in external data (eg, 2MASS)
-        ID_OBJ_GOOD            | // good-quality measurement in our data (eg,PS)
-        ID_OBJ_GOOD_ALT        | // good-quality measurement in  external data (eg, 2MASS)
-        ID_OBJ_GOOD_STACK      | // good-quality object in the stack (> 1 good stack)
-        ID_OBJ_BEST_STACK      | // the primary stack measurement are the best measurements
-        ID_OBJ_SUSPECT_STACK   | // suspect object in the stack (> 1 good or suspect stack, < 2 good)
-        ID_OBJ_BAD_STACK;        // good-quality object in the stack (> 1 good stack)
-      
-      for (j = 0; j < catalog.Naverage; j++) {
-        catalog.average[j].flags    &= ~photomBits; // reset all except astrometry bits
-        catalog.average[j].psfQF     = NAN;     // reset (will be re-calculated here)
-        catalog.average[j].psfQFperf = NAN;     // reset (will be re-calculated here)
-        catalog.average[j].stargal   = NAN;     // reset (will be re-calculated here)
-        catalog.average[j].photFlagsUpper = 0;  // reset (will be re-calculated here)
-        catalog.average[j].photFlagsLower = 0;  // reset (will be re-calculated here)
-        catalog.average[j].NwarpOK        = 0;  // reset (will be re-calculated here)
-        for (k = 0; k < Nsecfilt; k++) {
-          dvo_secfilt_init (&catalog.secfilt[j*Nsecfilt + k], SECFILT_RESET_ALL);
-        }
-      }
-    }
-
-    populate_tiny_values(&catalog, DVO_TV_MEASURE | DVO_TV_AVERAGE);
-
-    // XXX need to worry about the image subset data
-    initImageBins  (&catalog, 1, FALSE);
-    initMosaicBins (&catalog, 1, FALSE);
-    initGridBins   (&catalog, 1); 
-
-    findImages (&catalog, 1, FALSE);
-    findMosaics (&catalog, 1, FALSE);
-
-    // update the detection coordinates using the new image parameters
-    if (REPAIR_WARPS) RepairWarpMeasures (&catalog);
-
-    initMrel (&catalog, 1);
-    setMrelFinal (&catalog, NULL, TRUE);
-
-    // uncomment for extra verbosity
-    // MARKTIME("setMrelFinal for "OFF_T_FMT" average "OFF_T_FMT" measure : %f sec\n", catalog.Naverage, catalog.Nmeasure, dtime);
-    // gettimeofday (&startTimer, NULL); // reset timer
-
-    // XXX if we want to have options for setting warp, chip, stack independently, we need to init only the desired ones
-    // NOTE flatcorr == NULL, but it should have been applied already by setphot
-
-    if (!UPDATE) {
-      dvo_catalog_unlock (&catalog);
-      free_tiny_values(&catalog);
-      dvo_catalog_free (&catalog);
-      freeImageBins (1, FALSE);
-      freeMosaicBins (1, FALSE);
-      freeGridBins (1);
-      continue;
-    }
-    
-    if (VERBOSE) fprintf (stderr, "saving catalog %s\n", catalog.filename);
-    
-    char history[128];
-    struct timeval now;
-    gettimeofday (&now, (void *) NULL);
-    char *moddate = ohana_sec_to_date (now.tv_sec);
-    gfits_modify (&catalog.header, "RELPHOT", "%s", 1, moddate);      
-    if (REPAIR_WARPS) {
-      snprintf (history, 128, "repair warp measure.imageID: %s", moddate);
-      gfits_modify_alt (&catalog.header, "HISTORY", "%S", 0, history); // adds a new entry
-    }
-    free (moddate);
-
-    // we can optionally convert output format here
-    // but it would be better to define a dvo crawler program to do this
-    // catalog.catformat = DVO_FORMAT_PS1_V1;
-
-    // modify the output format as desired (ignore current format on disk)
-    if (UPDATE_CATFORMAT) {
-      catalog.catformat = dvo_catalog_catformat (UPDATE_CATFORMAT);
-    }
-    
-    SetProtect (TRUE);
-    if (!dvo_catalog_save (&catalog, VERBOSE)) { fprintf (stderr, "ERROR: failed to save %s\n", catalog.filename); exit (1); }
-    if (!dvo_catalog_unlock (&catalog)) { fprintf (stderr, "ERROR: failed to unlock %s\n", catalog.filename); exit (1); }
-    SetProtect (FALSE);
-
-    free_tiny_values(&catalog);
-    dvo_catalog_free (&catalog);
-
-    freeImageBins (1, FALSE);
-    freeMosaicBins (1, FALSE);
-    freeGridBins (1);
-  }
-
-  if (REPAIR_WARPS) {
-    FreeWarpGroups ();
-    FreeStackGroups ();
-    RepairWarpMeasuresCloseLogfile ();
-  }
+  reload_catalogs (skylist, hostID, hostpath);
 
   freeImages(db.ftable.buffer);
+  freeGridBins ();
+
   gfits_db_free (&db);
 
   return (TRUE);
 }
-
-// CATDIR is supplied globally
-# define DEBUG 1
-int relphot_objects_parallel (SkyList *sky) {
-
-  // launch the setphot_client jobs to the parallel hosts
-
-  // load the list of hosts
-  HostTable *table = HostTableLoad (CATDIR, sky->hosts);
-  if (!table) {
-    fprintf (stderr, "ERROR: failure reading Host Table %s for database %s\n", sky->hosts, CATDIR);
-    exit (1);
-  }    
-
-  if (BOUNDARY_TREE) {
-    char *tmppath = abspath(BOUNDARY_TREE, DVO_MAX_PATH);
-    free (BOUNDARY_TREE);
-    BOUNDARY_TREE = tmppath;
-  }
-
-  int i;
-  for (i = 0; i < table->Nhosts; i++) {
-
-    // ensure that the paths are absolute path names
-    char *tmppath = abspath (table->hosts[i].pathname, DVO_MAX_PATH);
-    free (table->hosts[i].pathname);
-    table->hosts[i].pathname = tmppath;
-
-    // options / arguments that can affect relphot_client -update-objects:
-    // VERBOSE, VERBOSE2
-    // KEEP_UBERCAL
-    // RESET (-reset)
-    // TimeSelect -time
-    // DophotSelect
-    // (note that psfQF is applied rigidly at 0.85, as is the galaxy test)
-    // MAG_LIM
-    // SIGMA_LIM
-    // ImagSelect, ImagMin, ImagMax
-    // MaxDensityUse, MaxDensityValue
-
-    char *command = NULL;
-    strextend (&command, "relphot_client -update-objects");
-    strextend (&command, "-hostID %d", table->hosts[i].hostID);
-    strextend (&command, "-D CATDIR %s", CATDIR);
-    strextend (&command, "-hostdir %s", table->hosts[i].pathname);
-    strextend (&command, "-region %f %f %f %f", UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax);
-    strextend (&command, "-statmode %s", STATMODE);
-    strextend (&command, "-D CAMERA %s", CAMERA);
-    strextend (&command, "-D STAR_TOOFEW %d", STAR_TOOFEW);
-    strextend (&command, "-minerror %f", MIN_ERROR);
-
-    if (VERBOSE)           { strextend (&command, "-v"); }
-    if (VERBOSE2)          { strextend (&command, "-vv"); }
-    if (RESET)             { strextend (&command, "-reset"); }
-    if (RESET_ZEROPTS)     { strextend (&command, "-reset-zpts"); }
-    if (REPAIR_WARPS)      { strextend (&command, "-repair-warps"); }
-    if (PRESERVE_PS1)      { strextend (&command, "-preserve-ps1"); }
-    if (IS_DIFF_DB)        { strextend (&command, "-is-diff-db"); }
-    if (UPDATE)            { strextend (&command, "-update"); }
-    if (!KEEP_UBERCAL)     { strextend (&command, "-reset-ubercal"); }
-    if (UPDATE_CATFORMAT)  { strextend (&command, "-update-catformat %s", UPDATE_CATFORMAT); }
-    if (BOUNDARY_TREE)     { strextend (&command, "-boundary-tree %s", BOUNDARY_TREE); }
-    if (SYNTH_ZERO_POINTS) { strextend (&command, "-synthphot-zpts %s", SYNTH_ZERO_POINTS); }
-    if (USE_BASIC_CHECK)   { strextend (&command, "-basic-image-search"); }
-    if (USE_ALL_IMAGES)    { strextend (&command, "-use-all-images"); }
-    if (USE_MCAL_PSF_FOR_STACK_APER) { strextend (&command, "-use-mcal-psf-for-stack-aper"); }
-
-    if (!(STAGES & STAGE_CHIP))  { strextend (&command, "-skip-chip"); }
-    if (!(STAGES & STAGE_WARP))  { strextend (&command, "-skip-warp"); }
-    if (!(STAGES & STAGE_STACK)) { strextend (&command, "-skip-stack"); }
-
-    // if (SET_MREL_VERSION != 1) { strextend (&command, "-set-mrel-version %d", SET_MREL_VERSION); } 
-
-    fprintf (stderr, "command: %s\n", command);
-
-    if (PARALLEL_MANUAL) {
-      free (command);
-      continue;
-    }
-
-    if (PARALLEL_SERIAL) {
-      int status = system (command);
-      if (status) {
-        fprintf (stderr, "ERROR running relphot_client\n");
-        exit (2);
-      }
-    } else {
-      // launch the job on the remote machine (no handshake)
-      int errorInfo = 0;
-      int pid = rconnect ("ssh", table->hosts[i].hostname, command, table->hosts[i].stdio, &errorInfo, FALSE);
-      if (!pid) {
-        if (DEBUG) fprintf (stderr, "failure to start %s (error %d)\n", table->hosts[i].hostname, errorInfo);
-        exit (1);
-      }
-      table->hosts[i].pid = pid; // save for future reference
-    }
-    free (command);
-  }
-
-  if (PARALLEL_MANUAL) {
-    fprintf (stderr, "run the relphot_client commands above.  when these are done, hit return\n");
-    getchar();
-  }
-  if (!PARALLEL_MANUAL && !PARALLEL_SERIAL) {
-    HostTableWaitJobsGetIO (table, __FILE__, __LINE__, VERBOSE);
-  }
-
-  FreeHostTable (table);
-  return TRUE;
-}      

trunk/Ohana/src/relphot/src/relphot_parallel_images.c

@@ -8 +8 @@
  * update the unowned detections for owned objects from neighbor regions
  * update the unowned objects for owned detections 
+
+ * TGROUP and GRID corrections are not calculated.
  */
 
@@ -39 +41 @@
 
   initImages (image, NULL, Nimage);
-
-  // load the flat correction table (if defined)
-  // we no longer need to do this: the flatcorrection is now saved in measure.Mflat
-# if (0)
-  char flatcorrfile[256];
-  sprintf (flatcorrfile, "%s/flatcorr.fits", CATDIR);
-  FlatCorrectionTable *flatcorr = FlatCorrectionLoad (flatcorrfile, VERBOSE);
-# endif
-  FlatCorrectionTable *flatcorr = NULL;
 
   // UserPatch.Rmin,Rmax may have range from a few degrees < 0.0 to few degrees > 360.0.
@@ -88 +81 @@
   setExclusions (catalog, Ncatalog, TRUE);
 
-  global_stats (catalog, Ncatalog, flatcorr, 0);
+  global_stats (catalog, Ncatalog, 0);
 
   if (PLOTSTUFF) {
@@ -101 +94 @@
 
     // set the mean stellar mags given the measurements and the image calibrations
-    setMrel  (catalog, Ncatalog, flatcorr);
+    setMrel  (catalog, Ncatalog);
 
     // share mean mags for objects at the boundary (number of unowned meas > 0)
@@ -112 +105 @@
 
     // set the image (Mcal) and mosaic (Mmos) zero point offsets given the mean mags and measurements
-    setMcal  (catalog, flatcorr);
-    setMmos  (catalog, flatcorr);
-    // setMgrp  (catalog, FALSE, flatcorr); XXX think this through: this may not make sense for tgroups
+    setMcal  (catalog);
+    setMmos  (catalog);
     MARKTIME("-- set Mrel, Mcal, Mmos, Mgrid : %f sec\n", dtime);
     
@@ -124 +116 @@
     client_logger_message ("slurped image data : loop %d \n", i);
 
-    global_stats (catalog, Ncatalog, flatcorr, i);
+    global_stats (catalog, Ncatalog, i);
 
-    SetZeroPointModes ();
+    SetZeroPointModes (catalog, Ncatalog);
     MARKTIME("-- finished loop %d: %f sec\n", i, dtime);
   }
@@ -147 +139 @@
 
   setMcalFromMosaics (); // copy per-mosaic calibrations to the images
-  // XXX keep? setMcalFromTGroups (); // copy per-tgroup calibrations to the images
 
   share_image_mags (regionHosts, -1);
@@ -161 +152 @@
   freeImageBins(Ncatalog, TRUE);
   freeMosaicBins (Ncatalog, TRUE);
-  freeGridBins (Ncatalog);
   freeImages((char *)image);
   free (image);

trunk/Ohana/src/relphot/src/relphot_parallel_regions.c

@@ -14 +14 @@
     exit (2);
   }
-
-  // load the flat correction table (if defined)
-# if (0)
-  char flatcorrfile[256];
-  sprintf (flatcorrfile, "%s/flatcorr.fits", CATDIR);
-  FlatCorrectionTable *flatcorr = FlatCorrectionLoad (flatcorrfile, VERBOSE);
-# endif
-  FlatCorrectionTable *flatcorr = NULL;
 
   // register database handle with shutdown procedure
@@ -81 +73 @@
 
   /* update catalogs (in parallel) */
-  reload_catalogs (skylist, flatcorr, 0, NULL);
+  reload_catalogs (skylist, 0, NULL);
 
   // save the changes to the image parameters

trunk/Ohana/src/relphot/src/select_images.c

@@ -103 +103 @@
     if (!(i % 300000)) fprintf (stderr, ".");
 
-    // XXX I am just going for force this for the moment:
-    // goto found_it;
-
-    /* exclude images by photcode (we are not doing this, but I'm not sure why) */
+    // only include active photcodes in the analysis
     if (timage[i].photcode) {
       int Ns = GetActivePhotcodeIndex (timage[i].photcode);

trunk/Ohana/src/relphot/src/setMrelCatalog.c

@@ -60 +60 @@
 int useOLS = TRUE; // chosen based on iteration below
 
-int setMrel_catalog_alt (Catalog *catalog, int Nc, int isSetMrelFinal, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt) {
+int setMrelCatalog (Catalog *catalog, int Nc, int isSetMrelFinal, SetMrelInfo *results, int Nsecfilt) {
 
   off_t j;
@@ -75 +75 @@
 
     if (STAGES & STAGE_CHIP) {
-      setMrelAverageExposure (&catalog[Nc], Nc, j, Nsecfilt, isSetMrelFinal, flatcorr, results);
+      setMrelAverageExposure (&catalog[Nc], Nc, j, Nsecfilt, isSetMrelFinal, results);
     }
 
     // only apply Stack operation on setMrelFinal in first pass 
     if (isSetMrelFinal && (STAGES & STAGE_STACK) && !IS_DIFF_DB) {
-      setMrelAverageStack (&catalog[Nc], Nc, j, Nsecfilt, flatcorr);
+      setMrelAverageStack (&catalog[Nc], Nc, j, Nsecfilt);
     }
 
     // only measure force-warp mean values if issetMrelFinal (make it optional?)
     if (isSetMrelFinal && (STAGES & STAGE_WARP)) {
-      setMrelAverageForcedWarp (&catalog[Nc], Nc, j, Nsecfilt, flatcorr, results);
+      setMrelAverageForcedWarp (&catalog[Nc], Nc, j, Nsecfilt, results);
     }
   }
@@ -113 +113 @@
 
 // set mean of chip measurements (selected by photcode range for now):
-int setMrelAverageExposure (Catalog *catalog, int cat, off_t ave, int Nsecfilt, int isSetMrelFinal, FlatCorrectionTable *flatcorr, SetMrelInfo *results) {
+int setMrelAverageExposure (Catalog *catalog, int cat, off_t ave, int Nsecfilt, int isSetMrelFinal, SetMrelInfo *results) {
 
   int Nsec;
@@ -225 +225 @@
     }
 
+    /* some things to note here:
+
+       Mcal, Mmos, Mgrid are only relevant (non-zero) during the relphot_images and
+       relphot_parallel_images analysis steps. For the final output steps (setMrelFinal),
+       these are turned off
+
+       During the image analysis, Mflat is the prior stored value while Mgrid is the new
+       measurement.  In the final output step, Mgrid is transferred to Mflat and zero
+
+       On the final calculation of Mrel, the value of Mgrid has been applied to Mflat, but
+       it has not been removed from the GridOps.c structures.  In order to avoid
+       double-counting, we need to skip Mgrid on the final calculation.
+    */
+
     // ** Choose the calibration (depends on the mode : do I have an image reference or not?) 
     float Mcal = 0.0, Mmos = 0.0, Mgrid = 0.0, Mflat = 0.0, Mgrp = 0.0;
@@ -233 +247 @@
       // data for which the associated image has not been loaded (probably because of
       // overlaps).  Msys + measure.Mcal is our best guess of the true magnitude
-      Mmos = Mgrid = 0;
       Mcal = measureT[k].McalPSF; // check that this is zero for loaded REF value
     } else {
-      // getMcal returns image[].Mcal; note the flat-field correction is stored in measure.Mflat
+      // getMcal returns image[].Mcal
       Mcal  = getMcal (meas, cat, MAG_CLASS_PSF);
       if (isnan(Mcal))  SKIP_THIS_MEAS(Ncal);
-      Mmos  = getMmos  (meas, cat);
-      if (isnan(Mmos))  SKIP_THIS_MEAS(Nmos);
-      Mgrp  = getMgrp  (meas, cat, measureT[k].airmass, &dMgrp);
-      if (isnan(Mgrp))  SKIP_THIS_MEAS(Ngrp);
-      Mgrid = getMgrid (meas, cat);
-      if (isnan(Mgrid)) SKIP_THIS_MEAS(Ngrid);
+
+      // the flat-field correction is stored in measure.Mflat
       Mflat = isnan (measureT[k].Mflat) ? 0.0 : measureT[k].Mflat;
+
+      // see note above re: final output vs image analysis
+      if (!isSetMrelFinal) {
+        Mmos  = getMmos  (meas, cat);
+        if (isnan(Mmos))  SKIP_THIS_MEAS(Nmos);
+        Mgrp  = getMgrp  (meas, cat, measureT[k].airmass, &dMgrp);
+        if (isnan(Mgrp))  SKIP_THIS_MEAS(Ngrp);
+        Mgrid = getMgridTiny (&measureT[k]); 
+        if (isnan(Mgrid)) SKIP_THIS_MEAS(Ngrid);
+      }
     }
 
@@ -269 +288 @@
     float dMap = NAN;
 
-    // XXX I am not convinced the sign on Mflat here is correct
-    float Moff =  Mflat + Mcal + Mmos + Mgrp + Mgrid;
+    // This definition is consistent with PhotRel: Mrel = Msys - Mcal - Mflat
+    float Moff =  Mcal + Mmos + Mgrp + Mgrid + Mflat;
 
     if (isSetMrelFinal) {
@@ -349 +368 @@
   for (Nsec = 0; Nsec < Nsecfilt; Nsec++) {
 
-    // if we detected this object in PS1, or do not request -preserve-ps1, keep the mean photometry values
+    // -preserve-ps1 means keep an existing average PS1 value
+    // if we did not detect this object in PS1, or do not request -preserve-ps1, keep the mean photometry values
     if (!PRESERVE_PS1 || !(secfilt[Nsec].flags & ID_SECF_HAS_PS1)) {
       dvo_secfilt_init (&secfilt[Nsec], SECFILT_RESET_CHIP); // this does not reset astrometry or STACK bits
@@ -364 +384 @@
     }
 
+    // if -preserve-ps1 is set and this object has PS1 data, skip the rest of the steps:
+    // NOTE: This test only applies if this bit was set by an earlier analysis
     if (PRESERVE_PS1 && (secfilt[Nsec].flags & ID_SECF_HAS_PS1)) continue;
-    // if -preserve-ps1 is set and this object has PS1 data, skip the rest of the steps:
 
     // XXX hardwired grizy = (01234) JHK = (567) w = (8)
@@ -416 +437 @@
     int NrankingPSF = useOLS ? magStatsByRanking (&results->psfData[Nsec], psfstats) : magStatsByRankingIRLS (&results->psfData[Nsec], psfstats);
     if (NrankingPSF < Nminmeas) { 
-      secfilt[Nsec].flags |= ID_OBJ_FEW;
+      secfilt[Nsec].flags |= ID_SECF_STAR_FEW;
     } else {
       secfilt[Nsec].MpsfChp  = psfstats->mean;
       secfilt[Nsec].dMpsfChp = psfstats->error;
       secfilt[Nsec].Mchisq   = (psfstats->Nmeas > 1) ? psfstats->chisq : NAN;
+      secfilt[Nsec].flags &= ~ID_SECF_STAR_FEW;
     }
 
@@ -552 +574 @@
 // measure->McalAPER to image->McalPSF for chips and warps, but not stacks
 
-int setMrelAverageStack (Catalog *catalog, int cat, off_t ave, int Nsecfilt, FlatCorrectionTable *flatcorr) {
+int setMrelAverageStack (Catalog *catalog, int cat, off_t ave, int Nsecfilt) {
 
   // we are guaranteed to have average, measure, secfilt
@@ -562 +584 @@
   SecFilt *secfilt = &catalog[0].secfilt[ave*Nsecfilt];
 
-  float McalPSF = 0, McalAPER = 0, Mmos = 0, Mgrid = 0;
+  float McalPSF = 0, McalAPER = 0;
 
   // set the primary projection cell and skycell for this coordinate
@@ -604 +626 @@
 
     // reset all stack-related values for this secfilt:
-    secfilt[Nsec].stackBestOff = -1;
-    secfilt[Nsec].stackPrmryOff = -1;
-    secfilt[Nsec].FpsfStk   = NAN;
-    secfilt[Nsec].dFpsfStk  = NAN;
-    secfilt[Nsec].FkronStk  = NAN;
-    secfilt[Nsec].dFkronStk = NAN;
-    secfilt[Nsec].FapStk    = NAN;
-    secfilt[Nsec].dFapStk   = NAN;
-    secfilt[Nsec].MpsfStk   = NAN;
-    secfilt[Nsec].MkronStk  = NAN;
-    secfilt[Nsec].MapStk    = NAN;
-    secfilt[Nsec].Nstack    = 0;
-    secfilt[Nsec].NstackDet = 0;
-    secfilt[Nsec].flags    &= ~ID_SECF_STACK_FLAGS; // reset the stack flags
+    dvo_secfilt_init (&secfilt[Nsec], SECFILT_RESET_STACK);
 
     for (off_t k = 0; k < Nmeasure; k++) {
@@ -638 +647 @@
       haveStack = TRUE;
       haveStackObject = TRUE;
-
-      // if (measure[k].dbFlags & MEAS_BAD) SKIP_THIS_MEAS_STACK(Nbad); 
 
       int isPrimary = FALSE;
@@ -726 +733 @@
       McalPSF  = measure[meas].McalPSF; // check that this is zero for loaded REF value
       McalAPER = McalPSF; // check that this is zero for loaded REF value
-      Mmos     = 0.0;
-      Mgrid    = 0.0;
     } else {
       McalPSF   = getMcal  (measSeq, cat, MAG_CLASS_PSF);
       McalAPER  = USE_MCAL_PSF_FOR_STACK_APER ? getMcal  (measSeq, cat, MAG_CLASS_PSF) : getMcal  (measSeq, cat, MAG_CLASS_KRON);
-      Mmos      = getMmos  (measSeq, cat);
-      Mgrid     = getMgrid (measSeq, cat);
-      // XXX can Mmos and Mgrid exist for stacks?
     }
     
@@ -769 +771 @@
     // get the zero point for the selected image
     // Use a different zero point for the PSF-like and APERTURE-like magnitudes
-    float zpPSF  = PhotZeroPoint (&measure[meas], &average[0], &secfilt[0]) - (McalPSF  + Mmos + Mgrid);
-    float zpAPER = PhotZeroPoint (&measure[meas], &average[0], &secfilt[0]) - (McalAPER + Mmos + Mgrid);
+    float zpPSF  = PhotZeroPoint (&measure[meas], &average[0], &secfilt[0]) - McalPSF;
+    float zpAPER = PhotZeroPoint (&measure[meas], &average[0], &secfilt[0]) - McalAPER;
 
     // zpFactor to go from instrumental flux to Janskies
@@ -884 +886 @@
 // chips and warps, but not stacks
 
-int setMrelAverageForcedWarp (Catalog *catalog, int cat, off_t ave, int Nsecfilt, FlatCorrectionTable *flatcorr, SetMrelInfo *results) {
-  OHANA_UNUSED_PARAM(flatcorr);
+int setMrelAverageForcedWarp (Catalog *catalog, int cat, off_t ave, int Nsecfilt, SetMrelInfo *results) {
 
   // we are guaranteed to have average, measure, secfilt

trunk/Ohana/src/relphot/src/setMrelFinal.c

@@ -5 +5 @@
 // output dbFlags values
 
-void setMrelFinal (Catalog *catalog, FlatCorrectionTable *flatcorr, int simpleAverage) {
+void setMrelFinal (Catalog *catalog, int simpleAverage) {
 
   off_t i;
@@ -11 +11 @@
   int Nsecfilt = GetPhotcodeNsecfilt ();
 
-  /*** RESET photometry flags and Mcal values ***/
-  if (RESET) {
-    // flags used by the photometry analysis (excluding UBERCAL)
-    unsigned int PHOTOM_FLAGS = 
-      ID_MEAS_NOCAL |          // detection ignored for this analysis (photcode, time range)
-      ID_MEAS_POOR_PHOTOM |    // detection is photometry outlier 
-      ID_MEAS_SKIP_PHOTOM |    // detection was ignored for photometry measurement 
-      ID_MEAS_AREA |           // detetion was outside acceptable area of device
-      ID_MEAS_SYNTH_MAG |      // magnitude is synthetic
-      ID_MEAS_STACK_PRIMARY |  // this stack measurement is in the primary skycell
-      ID_MEAS_STACK_PHOT_SRC;  // this measurement supplied the stack photometry
-
-    // ID_MEAS_PHOTOM_UBERCAL -- externally-supplied zero point from ubercal analysis
-    // this is set by 'setphot', do not reset here
-
-    // flags used by the photometry analysis (excluding UBERCAL)
-    // unsigned int secfiltFlags = 
-    //   ID_PHOTOM_PASS_0 | // average measured at pass 0
-    //   ID_PHOTOM_PASS_1 | // average measured at pass 1
-    //   ID_PHOTOM_PASS_2 | // average measured at pass 2
-    //   ID_PHOTOM_PASS_3 | // average measured at pass 3
-    //   ID_PHOTOM_PASS_4 | // average measured at pass 3
-    //   ID_SECF_USE_SYNTH | // average measured at pass 3
-    //   ID_SECF_USE_UBERCAL | // average measured at pass 3
-    //   ID_SECF_OBJ_EXT; // average measured at pass 3
-    
-    // XXX I should really deprecate the concept of applying the average 
-    // calculation to a limited set of photcodes.  
-    // for now, just do all photcodes here (
-    // as it stands, only stacks are limited by photcode; mean exp and forced warp 
-    // are applied to all Nsecfilt
-
-    for (i = 0; i < catalog[0].Naverage; i++) {
-      int Ns;
-      for (Ns = 0; Ns < Nsecfilt; Ns++) {
-        
-        off_t N = Nsecfilt*i+Ns;
-        dvo_secfilt_init (&catalog[0].secfilt[N], SECFILT_RESET_ALL);
-
-        off_t m = catalog[0].average[i].measureOffset;
-        off_t j;
-        for (j = 0; j < catalog[0].average[i].Nmeasure; j++, m++) {
-        
-          /* select measurements by time */
-          if (TimeSelect) {
-            if (catalog[0].measure[m].t < TSTART) continue;
-            if (catalog[0].measure[m].t > TSTOP) continue;
-          }
-        
-          // only reset Mcal for measures with a matching image
-          // do not reset Mcal for ubercal images unless explicitly requested
-
-          if (catalog[0].measure[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) {
-            if (!KEEP_UBERCAL) {
-              catalog[0].measure[m].McalPSF  = 0.0;
-              catalog[0].measure[m].McalAPER = 0.0;
-              catalog[0].measure[m].dbFlags &= ~ID_MEAS_PHOTOM_UBERCAL;
-            } 
-          } else {
-            if (RESET_ZEROPTS && (getImageEntry (m, 0) >= 0)) {
-              catalog[0].measure[m].McalPSF  = 0.0;
-              catalog[0].measure[m].McalAPER = 0.0;
-            }
-          }
-          catalog[0].measure[m].dbFlags &= ~PHOTOM_FLAGS;
-        }
-      }
-    }
-  }
-
-  // this sets flags in the measureT element, not the measure element
+  ResetAverageAndMeasure (catalog);
+  
+  // This sets ID_MEAS flags in the measureT element, not the measure element
+  // These bits are thus ephemeral and not saved.  
   setExclusions (catalog, 1, VERBOSE);  /* mark by area */
 
@@ -92 +25 @@
   }
 
-  for (i = 0; i < catalog[0].Naverage; i++) {
-    catalog[0].average[i].psfQF     = NAN;      // force recalculation below
-    catalog[0].average[i].psfQFperf = NAN;      // force recalculation below
-    catalog[0].average[i].stargal   = NAN;      // force recalculation below
-    catalog[0].average[i].photFlagsUpper = 0;   // reset (will be re-calculated here)
-    catalog[0].average[i].photFlagsLower = 0;   // reset (will be re-calculated here)
-    catalog[0].average[i].NwarpOK        = 0;   // reset (will be re-calculated here)
-  }
-
-  // XXX make this optional? (do not clean for -averages?)
-  // XXX : for the moment, disable clean_measures
-  if (!simpleAverage && FALSE) clean_measures (catalog, 1, TRUE, flatcorr);    /* mark outliers ID_MEAS_POOR_PHOTOM */
-
   ALLOCATE (catalog[0].measureRank, char, catalog[0].Nmeasure);
   setMeasureRank (catalog);
-  setMrelOutput (catalog, 1, flatcorr); // sets the values secfilt.MpsfChp = <measure.M + measure.Mflat - image.Mcal>
-  setMcalOutput (catalog, 1, flatcorr); // sets measure.Mcal = image.Mcal
 
-  /* clear ID_OBJ_POOR, ID_OBJ_FEW values before writing ??? */
-  /* ID_MEAS_NOCAL is an internal bit, so it should be cleared */
-  // XXX eventually: move these internal flags to an internal-only catalog->FOO element (like found)
-  for (i = 0; i < catalog[0].Naverage; i++) {
-    catalog[0].average[i].flags &= ~ID_OBJ_FEW;
-    catalog[0].average[i].flags &= ~ID_OBJ_POOR;
-    off_t j;
-    for (j = 0; j < Nsecfilt; j++) {
-        catalog[0].secfilt[i*Nsecfilt+j].flags &= ~ID_OBJ_FEW; 
-        catalog[0].secfilt[i*Nsecfilt+j].flags &= ~ID_OBJ_POOR; 
-    }
-    off_t m = catalog[0].average[i].measureOffset;
-    for (j = 0; j < catalog[0].average[i].Nmeasure; j++, m++) {
-      catalog[0].measure[m].dbFlags &= ~ID_MEAS_NOCAL;
-    }
-  }
+  setMflatFromGrid (catalog); // Mgrid is used to set Mflat; Mgrid is the ignored in setMrelOutput / setMrelCatalog
+  setMrelOutput (catalog, 1); // sets the values secfilt.MpsfChp = <measure.M - image.Mcal - measure.Mflat>
+  setMcalOutput (catalog, 1); // sets measure.Mcal = image.Mcal
 }
-
-/* ID_MEAS_SKIP marks measurements which were not used to calculate Mrel */
-// For each star & photcode, we are selecting the measurements to use in the average
-// photometry value.  This function is called several times, lowering the bar to
-// acceptance on each pass. The rejections at the different passes are:
-
-// 0 : only use measurements thought to be GOOD (photflags not POOR or BAD)
-// 1 : accept measurements thought to be POOR based on photflags
-// 2 : accept the measurements marked as outliers; accept images marked as outliers
-// 3 : accept measurements thought to be BAD based on photflags
-// 4 : accept the measurements outside of the instrumental magnitude limits (eg, SAT)
-
-// XXX if I'm reading this correctly, passes 0, 1, 2 are identical & accept all
-// non-outliers with Imag range
 
 // int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure) {
@@ -161 +52 @@
 // }
 
-int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure);
-
-void skip_measurements (Catalog *catalog, int pass, FlatCorrectionTable *flatcorr) {
-
-  off_t i, k, m;
-  off_t Ntot, Ntry, Nkeep, Nskip;
-  float mag;
-
-  int Nsecfilt = GetPhotcodeNsecfilt ();
-
-  Ntot = Ntry = Nskip = Nkeep = 0;
-
-  /* allow measures from images marked POOR and FEW */
-  if (pass >= 3) IMAGE_BAD = ID_IMAGE_PHOTOM_NOCAL;
-  
-  /* allow measures marked as outliers (POOR) and off image region (AREA) */
-  if (pass >= 3) {
-    MEAS_BAD = ID_MEAS_NOCAL | ID_MEAS_SKIP_PHOTOM;
-  } else {
-    MEAS_BAD  = ID_MEAS_NOCAL | ID_MEAS_POOR_PHOTOM | ID_MEAS_SKIP_PHOTOM | ID_MEAS_AREA;
-  }
-
-  /* mark measures which should be ignored on this pass */
-  for (i = 0; i < catalog[0].Naverage; i++) {
-    Ntot += catalog[0].average[i].Nmeasure;
-
-    if (FALSE && (catalog[0].average[i].objID == 0x46a4) && (catalog[0].average[i].catID == 0xf40e)) {
-      fprintf (stderr, "test obj\n");
-      print_measure_set (&catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt], catalog[0].measure);
-    }
-
-    // mark measurements for each secfilt separately
-    int Ns;
-    for (Ns = 0; Ns < Nphotcodes; Ns++) {
-      
-      int thisCode = photcodes[Ns][0].code;
-      int Nsec = GetPhotcodeNsec(thisCode);
-      
-      /* star/photcodes already calibrated */
-      if (catalog[0].found_t[Nsecfilt*i+Nsec]) continue;  
-
-      m = catalog[0].average[i].measureOffset;
-      for (k = 0; k < catalog[0].average[i].Nmeasure; k++, m++) {
-        Ntry++;
-
-        // skip measurements not related to this photcode
-        PhotCode *code = GetPhotcodebyCode (catalog[0].measure[m].photcode);
-        if (!code) continue;
-        if (code->equiv != thisCode) continue;
-
-        // skip measurements by time range (mark as skipped for this secfilt)
-        // XXX note that this is a bit dangerous : some objects may never get calibrated 
-        if (TimeSelect) {
-          if (catalog[0].measure[m].t < TSTART) goto skip;
-          if (catalog[0].measure[m].t > TSTOP) goto skip;
-        }
-        
-        // we now have a measurement relavant to this photcode (and time range)
-
-        // clear SKIP for all measures at first 
-        catalog[0].measureT[m].dbFlags &= ~ID_MEAS_SKIP_PHOTOM;
-        catalog[0].measure [m].dbFlags &= ~ID_MEAS_SKIP_PHOTOM;
-
-        // skip measurements from BAD images and mosaics (not REF mags)
-        // do NOT skip measurements without a matching image (REF mags)
-        off_t Nim = getImageEntry (m, 0);
-        if (Nim > -1) {
-          if (isnan(getMcal (m, 0, MAG_CLASS_PSF))) goto skip;
-          if (isnan(getMmos (m, 0))) goto skip;
-        }
-
-        // PASS 4 : skip measurements by inst mag limit (not REF mags)
-        if ((pass < 4) && ImagSelect) {
-          if (Nim > -1) {
-            mag = PhotInst (&catalog[0].measure[m], MAG_CLASS_PSF);
-            if (mag < ImagMin) goto skip;
-            if (mag > ImagMax) goto skip;
-          }
-        }
-
-        // PASS 3 : accept bad measurements (eg, SAT, CR), internal outliers
-        if (pass < 3) {
-          if (catalog[0].measure[m].photFlags & code->photomBadMask) goto skip;
-          if ((catalog[0].measure[m].photcode > 10000) && (catalog[0].measure[m].photcode < 10500)) {
-            if (catalog[0].measure[m].psfQF < 0.85) goto skip;
-          }
-        }
-        
-        // PASS 2 : (no additional cut)
-
-        // PASS 1 : accept poor measurements as well (eg, POOR FIT, etc)
-        if ((pass < 1) && (catalog[0].measure[m].photFlags & code->photomPoorMask)) goto skip;
-
-        // PASS 0 : skip poor, outlier, bad measurements as well (eg, POOR FIT, etc)
-
-        Nkeep ++;
-        
-        continue;
-
-      skip:
-        catalog[0].measure [m].dbFlags |= ID_MEAS_SKIP_PHOTOM;
-        catalog[0].measureT[m].dbFlags |= ID_MEAS_SKIP_PHOTOM;
-        Nskip ++;
-      }
-    }
-  }
-  if (VERBOSE) fprintf (stderr, "pass %d, Ntot: "OFF_T_FMT", Ntry: "OFF_T_FMT", Nskip: "OFF_T_FMT", Nkeep: "OFF_T_FMT"\n",
-                        pass, Ntot, Ntry, Nskip, Nkeep);
-}
-
+// This function is only called for the final output step.  By this point, we have
+// propagated the mosaic and tgroup flags to each image.
 void setMeasureRank (Catalog *catalog) {
 
   int i;
 
-  // make these global
-  int IMAGE_BAD  = ID_IMAGE_PHOTOM_NOCAL;
-  int IMAGE_POOR = ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP;
-  int MEAS_BAD   = ID_MEAS_NOCAL        | ID_MEAS_SKIP_PHOTOM;
-  int MEAS_POOR  = ID_MEAS_POOR_PHOTOM  | ID_MEAS_AREA;
-
   Measure     *measure     = catalog[0].measure;
-  MeasureTiny *measureT    = catalog[0].measureT;
   char        *measureRank = catalog[0].measureRank;
 
@@ -289 +65 @@
     measureRank[i] = 11; // start at a low rank
 
-    // clear SKIP for all measures at first 
-    measureT[i].dbFlags &= ~ID_MEAS_SKIP_PHOTOM;
-    measure [i].dbFlags &= ~ID_MEAS_SKIP_PHOTOM;
-
-    // skip measurements without a valid photcode
+    // measurements without a valid photcode have lowest rank (should not be used anyway)
     PhotCode *code = GetPhotcodebyCode (measure[i].photcode);
     if (!code) continue;
@@ -299 +71 @@
     // measurements outside time range have poor rank
     if (TimeSelect) {
-      if (measure[i].t < TSTART)                     { measureRank[i] = 10; continue; }
-      if (measure[i].t > TSTOP)                      { measureRank[i] = 10; continue; }
+      if (measure[i].t < TSTART)                                        { measureRank[i] = 10; continue; }
+      if (measure[i].t > TSTOP)                                         { measureRank[i] = 10; continue; }
     }
     
-    // measurements from BAD images and mosaics (not REF mags)
-    off_t Nim = getImageEntry (i, 0);
+    // Nim < 0 for REF mags, imageFlags have bits for IMAGE, MOSAIC, NIGHT
+    off_t Nim      = getImageEntry (i, 0);
     int imageFlags = getImageFlags (i, 0);
-    int mosaicFlags = getMosaicFlags (i, 0);
 
     if (Nim > -1) {
-      if (imageFlags  & IMAGE_BAD)                   { measureRank[i] = 9; continue; }
-
-      if (mosaicFlags & IMAGE_BAD)                   { measureRank[i] = 9; continue; }
-    
       // measurements ranked by inst mag limit (not REF mags)
       if (ImagSelect) {
         float mag = PhotInst (&measure[i], MAG_CLASS_PSF);
-        if (mag < ImagMin)                           { measureRank[i] = 8; continue; }
-        if (mag > ImagMax)                           { measureRank[i] = 8; continue; }
+        if (mag < ImagMin)                                              { measureRank[i] = 9; continue; }
+        if (mag > ImagMax)                                              { measureRank[i] = 9; continue; }
       }
     }
 
-    // RANK 1 : BAD photFlags
-    if (measure[i].dbFlags & MEAS_BAD)               { measureRank[i] = 7; continue; }
-
-    // RANK 3 : BAD photFlags (eg, SAT, CR), internal outliers
-    if (measure[i].photFlags & code->photomBadMask)  { measureRank[i] = 6; continue; }
-
+    // RANK 8 : Poor image
     if (Nim > -1) {
-      if (imageFlags  & IMAGE_POOR)                  { measureRank[i] = 5; continue; }
-      if (mosaicFlags & IMAGE_POOR)                  { measureRank[i] = 5; continue; }
+      if (imageFlags & ID_IMAGE_PHOTOM_POOR)                            { measureRank[i] = 8; continue; }
     }
 
-    // RANK 2 : psfQF value
-    if (!isfinite(measure[i].psfQF) || measure[i].psfQF < 0.85) { measureRank[i] = 4; continue; }
+    // RANK 7 : BAD photFlags (eg, SAT, CR), internal outliers
+    if (measure[i].photFlags & code->photomBadMask)                     { measureRank[i] = 7; continue; }
+
+    // RANK 6 : bad psfQF value
+    if (!isfinite(measure[i].psfQF) || measure[i].psfQF < 0.85)         { measureRank[i] = 6; continue; }
         
-    // RANK 1 : POOR photFlags
-    if (measure[i].dbFlags & MEAS_POOR)              { measureRank[i] = 3; continue; }
+    // RANK 5 : not in valid chip region 
+    if (measure[i].dbFlags & ID_MEAS_AREA)                              { measureRank[i] = 5; continue; }
 
-    // RANK 1 : POOR photFlags
-    if (measure[i].photFlags & code->photomPoorMask) { measureRank[i] = 2; continue; }
+    // RANK 4 : POOR photFlags
+    if (measure[i].photFlags & code->photomPoorMask)                    { measureRank[i] = 4; continue; }
     
-    // RANK 1 : psfQFperfect value
-    if (!isfinite(measure[i].psfQFperf) || measure[i].psfQFperf < 0.85) { measureRank[i] = 1; continue; }
+    // RANK 3 : bad psfQFperfect value
+    if (!isfinite(measure[i].psfQFperf) || measure[i].psfQFperf < 0.85) { measureRank[i] = 3; continue; }
         
+    // RANK 2 : Poor mosaic
+    if (Nim > -1) {
+      if (imageFlags & ID_IMAGE_MOSAIC_POOR)                            { measureRank[i] = 2; continue; }
+    }
+
+    // RANK 1 : Poor night
+    if (Nim > -1) {
+      if (imageFlags & ID_IMAGE_NIGHT_POOR)                             { measureRank[i] = 1; continue; }
+    }
     // RANK 0 : perfect measurement:
     measureRank[i] = 0;