Changeset 37037

Timestamp:

Jul 17, 2014, 10:12:36 AM (12 years ago)

Author:

eugene

Message:

major re-work of mean photometry calculation to handle stack, forced warp data

Location:

trunk/Ohana/src/relphot

Files:

• Makefile (modified) (2 diffs)
• doc/mosaic.txt (modified) (1 diff)
• doc/pv2.txt (copied) (copied from branches/eam_branches/ipp-20140610/Ohana/src/relphot/doc/pv2.txt )
• include/relphot.h (modified) (6 diffs)
• src/BrightCatalog.c (modified) (13 diffs)
• src/GridOps.c (modified) (3 diffs)
• src/ImageMagIO.c (modified) (2 diffs)
• src/ImageOps.c (modified) (17 diffs)
• src/ImageSubset.c (modified) (2 diffs)
• src/MeanMagIO.c (modified) (6 diffs)
• src/MosaicOps.c (modified) (39 diffs)
• src/Shutdown.c (modified) (1 diff)
• src/StarOps.c (modified) (14 diffs)
• src/args.c (modified) (1 diff)
• src/assign_images.c (modified) (1 diff)
• src/bcatalog.c (modified) (10 diffs)
• src/extra.c (copied) (copied from branches/eam_branches/ipp-20140610/Ohana/src/relphot/src/extra.c )
• src/initialize.c (modified) (4 diffs)
• src/launch_region_hosts.c (modified) (1 diff)
• src/load_catalogs.c (modified) (5 diffs)
• src/plot_scatter.c (modified) (1 diff)
• src/relphot_client.c (modified) (1 diff)
• src/relphot_images.c (modified) (1 diff)
• src/relphot_parallel_images.c (modified) (2 diffs)
• src/relphot_parallel_regions.c (modified) (3 diffs)
• src/setExclusions.c (modified) (1 diff)
• src/setMrelCatalog.c (modified) (19 diffs)
• src/setMrelFinal.c (modified) (3 diffs)
• src/share_mean_mags.c (modified) (1 diff)

trunk/Ohana/src/relphot/Makefile

@@ -34 +34 @@
 $(SRC)/args.$(ARCH).o            \
 $(SRC)/help.$(ARCH).o            \
+$(SRC)/extra.$(ARCH).o           \
 $(SRC)/bcatalog.$(ARCH).o        \
 $(SRC)/global_stats.$(ARCH).o    \
@@ -82 +83 @@
 $(SRC)/args.$(ARCH).o            \
 $(SRC)/help.$(ARCH).o            \
+$(SRC)/extra.$(ARCH).o           \
 $(SRC)/synthetic_mags.$(ARCH).o  \
 $(SRC)/plotstuff.$(ARCH).o       \
 $(SRC)/liststats.$(ARCH).o       \
 $(SRC)/initialize.$(ARCH).o      \
+$(SRC)/syncfile.$(ARCH).o        \
 $(SRC)/load_catalogs.$(ARCH).o   \
 $(SRC)/reload_catalogs.$(ARCH).o \

trunk/Ohana/src/relphot/doc/mosaic.txt

@@ +1 @@
+
+
+2014.07.09
+
+* trying to make sense of mosaic.Mcal vs image.Mcal in the context of UBERCAL:
+
+  * relphot_images:
+    * load_images (convert raw FITS table to Image structure, select subset matching selection)
+      * initMosaics (associate mosaics to images (gpc1 photcodes only))
+        * initMosaicGrid (define spatial range of mosaics [only used by grid analysis]; set Mcal to <image.Mcal> and image.Mcal to 0.0)
+    --- calculate image or mosaic Mcal values
+    * setMcalFinal (set image.Mcal = mosaicMcal)
+
+  * relphot_parallel_regions:
+    * assign_images
+      * makeMosaics (equivalent to initMosaics, but works on full image table, not subset)
+      * setMosaicCenters (set Mcal to <image.Mcal> and image.Mcal to 0.0) <--- this is wrong!
+      (sends images ONLY, not mosaics, to remote machines)
+      (slurps back new image values, applies to db)
+      
+  * relphot_parallel_images:
+    * makeMosaics (associate mosaics to images, gpc1 photcodes only)
+      * initMosaicGrid (needed in this one)
+---
 
 For 'mosaic' zero points, I define the mosaics by grouping the images

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

@@ -4 +4 @@
 # include <signal.h>
 # include <pthread.h>
-
-# define MARKTIME(MSG,...) {                    \
-    gettimeofday (&stopTimer, (void *) NULL);   \
-    float dtime = DTIME (stopTimer, startTimer);        \
-    fprintf (stderr, MSG, __VA_ARGS__); }
-
-# define INITTIME \
-  struct timeval startTimer, stopTimer; \
-  gettimeofday (&startTimer, (void *) NULL);
 
 /* # define GRID_V1 */
@@ -119 +110 @@
   float M;
   float dM;
-  float Xm;
+  float Mchisq;
   int Nsec;
   unsigned int objID;
@@ -130 +121 @@
   float dMcal;
   float dMagSys;
-  float Xm;
+  short Xm;
   int nFitPhotom;
   int flags;
@@ -317 +308 @@
 int           findMosaics         PROTO((Catalog *catalog, int Ncatalog, int doMosaicList));
 
-void makeMosaics (Image *image, off_t Nimage);
+void makeMosaics (Image *image, off_t Nimage, int mergeMcal);
 Mosaic *getMosaicForImage (off_t im);
 void setMosaicCenters (Image *image, off_t Nimage);
@@ -360 +351 @@
 void          liststats_setmode   PROTO((StatType *stats, char *strmode));
 int           liststats           PROTO((double *value, double *dvalue, double *wvalue, int N, StatType *stats));
-Catalog      *load_catalogs       PROTO((SkyList *skylist, int *Ncatalog, int hostID, char *hostpath));
-Catalog      *load_catalogs_parallel PROTO((SkyList *sky, int *Ncatalog));
+Catalog      *load_catalogs       PROTO((SkyList *skylist, int *Ncatalog, int hostID, char *hostpath, char *syncfile));
+Catalog      *load_catalogs_parallel PROTO((SkyList *sky, int *Ncatalog, char *syncfile));
 
 SkyList      *load_images         PROTO((FITS_DB *db, char *regionName, SkyRegion *region));
@@ -509 +500 @@
 int relphot_parallel_images ();
 int relphot_parallel_regions ();
+
+// in extra.c
+int isGPC1chip  (int photcode);
+int isGPC1stack (int photcode);
+int isGPC1warp  (int photcode);
+
+

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

@@ -63 +63 @@
  
     // need to create and assign to flat-field correction
-    GET_COLUMN(dR,        "RA_OFF",     float);
-    GET_COLUMN(dD,        "DEC_OFF",    float);
+    GET_COLUMN(R,         "RA",         double);
+    GET_COLUMN(D,         "DEC",        double);
     GET_COLUMN(M,         "MAG_SYS",    float);
     GET_COLUMN(Mcal,      "MAG_CAL",    float);
@@ -85 +85 @@
     ALLOCATE (measure, MeasureTiny, Nrow);
     for (i = 0; i < Nrow; i++) {
-      measure[i].dR        = dR[i];
-      measure[i].dD        = dD[i];
+      measure[i].R         = R[i];
+      measure[i].D         = D[i];
       measure[i].M         = M[i];
       measure[i].Mcal      = Mcal[i];
@@ -105 +105 @@
     fprintf (stderr, "loaded data for %lld measure\n", (long long) Nrow);
 
-    free (dR      );
-    free (dD      );
+    free (R       );
+    free (D       );
     free (M       );
     free (Mcal    );
@@ -136 +136 @@
  
     // need to create and assign to flat-field correction
-    GET_COLUMN(R,             "RA",        double);
-    GET_COLUMN(D,             "DEC",       double);
+    GET_COLUMN(R,             "RA",          double);
+    GET_COLUMN(D,             "DEC",         double);
     GET_COLUMN(Nmeasure,      "NMEAS",       int);
     GET_COLUMN(measureOffset, "MEASURE_OFF", int);
@@ -182 +182 @@
  
     // need to create and assign to flat-field correction
-    GET_COLUMN(M,     "MAG",      float);
-    GET_COLUMN(dM,    "MAG_ERR",  float);
-    GET_COLUMN(Xm,    "MAG_CHI",  float);
-    GET_COLUMN(flags, "FLAGS",    int);
-    GET_COLUMN(Ncode, "NCODE",    short);
-    GET_COLUMN(Nused, "NUSED",    short);
-    GET_COLUMN(M_20,  "MAG_20",   short);
-    GET_COLUMN(M_80,  "MAG_80",   short);
+    GET_COLUMN(M,      "MAG",      float);
+    GET_COLUMN(dM,     "MAG_ERR",  float);
+    GET_COLUMN(Mchisq, "MAG_CHI",  float);
+    GET_COLUMN(flags,  "FLAGS",    int);
+    GET_COLUMN(Ncode,  "NCODE",    short);
+    GET_COLUMN(Nused,  "NUSED",    short);
+    GET_COLUMN(Mmin,   "MAG_MIN",  float);
+    GET_COLUMN(Mmax,   "MAG_MAX",  float);
     gfits_free_header (&theader);
     gfits_free_table  (&ftable);
@@ -196 +196 @@
     ALLOCATE (secfilt, SecFilt, Nrow);
     for (i = 0; i < Nrow; i++) {
-      secfilt[i].M     = M[i];         
-      secfilt[i].dM    = dM[i];
-      secfilt[i].Xm    = Xm[i];
-      secfilt[i].flags = flags[i];
-      secfilt[i].Ncode = Ncode[i];
-      secfilt[i].Nused = Nused[i];
-      secfilt[i].M_20  = M_20[i];
-      secfilt[i].M_80  = M_80[i];
+      secfilt[i].M      = M[i];         
+      secfilt[i].dM     = dM[i];
+      secfilt[i].Mchisq = Mchisq[i];
+      secfilt[i].flags  = flags[i];
+      secfilt[i].Ncode  = Ncode[i];
+      secfilt[i].Nused  = Nused[i];
+      secfilt[i].Mmin   = Mmin[i];
+      secfilt[i].Mmax   = Mmax[i];
     }
     fprintf (stderr, "loaded data for %lld secfilt\n", (long long) Nrow);
 
-    free (M    );
-    free (dM   );
-    free (Xm   );
-    free (flags);
-    free (Ncode);
-    free (Nused);
-    free (M_20 );
-    free (M_80 );
+    free (M     );
+    free (dM    );
+    free (Mchisq);
+    free (flags );
+    free (Ncode );
+    free (Nused );
+    free (Mmin  );
+    free (Mmax  );
     catalog->secfilt = secfilt;
     // assert Nsecfilt * Naverage = Nrow?
@@ -265 +265 @@
     gfits_create_table_header (&theader, "BINTABLE", "MEASURE_TINY");
 
-    gfits_define_bintable_column (&theader, "E", "RA_OFF",   "ra offset",                  "arcsec", 1.0, 0.0);
-    gfits_define_bintable_column (&theader, "E", "DEC_OFF",  "dec offset",                 "arcsec", 1.0, 0.0);
+    gfits_define_bintable_column (&theader, "D", "RA",       "ra",                         "degree", 1.0, 0.0);
+    gfits_define_bintable_column (&theader, "D", "DEC",      "dec",                        "degree", 1.0, 0.0);
     gfits_define_bintable_column (&theader, "E", "MAG_SYS",  "magnitude (sys)",             NULL,    1.0, 0.0);
     gfits_define_bintable_column (&theader, "E", "MAG_CAL",  "magnitude (cal)",             NULL,    1.0, 0.0);
@@ -274 +274 @@
     gfits_define_bintable_column (&theader, "E", "Y_CCD",    "ccd y coord",                "pix",    1.0, 0.0);
     gfits_define_bintable_column (&theader, "E", "EXPTIME",  "-2.5 * log (exposure time)", "sec",    1.0, 0.0);
-    gfits_define_bintable_column (&theader, "J", "TIME",     "time of exp",                "sec",    1.0, 1.0*0x8000);
-    gfits_define_bintable_column (&theader, "J", "AVE_REF",  "pointer to average table",    NULL,    1.0, 1.0*0x8000);
-    gfits_define_bintable_column (&theader, "J", "IMAGE_ID", "image",                       NULL,    1.0, 1.0*0x8000);
-    gfits_define_bintable_column (&theader, "J", "DB_FLAGS", "flags",                       NULL,    1.0, 1.0*0x8000);
-    gfits_define_bintable_column (&theader, "J", "PHOT_FLAGS", "photflags",                 NULL,    1.0, 1.0*0x8000);
-    gfits_define_bintable_column (&theader, "J", "CAT_ID",   "catalog",                     NULL,    1.0, 1.0*0x8000);
-    gfits_define_bintable_column (&theader, "I", "PHOTCODE", "photcode",                    NULL,    1.0, 1.0*0x80);
+    gfits_define_bintable_column (&theader, "J", "TIME",     "time of exp",                "sec",    1.0, FT_BZERO_INT32);
+    gfits_define_bintable_column (&theader, "J", "AVE_REF",  "pointer to average table",    NULL,    1.0, FT_BZERO_INT32);
+    gfits_define_bintable_column (&theader, "J", "IMAGE_ID", "image",                       NULL,    1.0, FT_BZERO_INT32);
+    gfits_define_bintable_column (&theader, "J", "DB_FLAGS", "flags",                       NULL,    1.0, FT_BZERO_INT32);
+    gfits_define_bintable_column (&theader, "J", "PHOT_FLAGS", "photflags",                 NULL,    1.0, FT_BZERO_INT32);
+    gfits_define_bintable_column (&theader, "J", "CAT_ID",   "catalog",                     NULL,    1.0, FT_BZERO_INT32);
+    gfits_define_bintable_column (&theader, "I", "PHOTCODE", "photcode",                    NULL,    1.0, FT_BZERO_INT16);
 
     // generate the output array that carries the data
@@ -286 +286 @@
 
     // create intermediate storage arrays
-    float *dR        ; ALLOCATE (dR       ,  float, catalog->Nmeasure);
-    float *dD        ; ALLOCATE (dD       ,  float, catalog->Nmeasure);
-    float *M         ; ALLOCATE (M        ,  float, catalog->Nmeasure);
-    float *Mcal      ; ALLOCATE (Mcal     ,  float, catalog->Nmeasure);
-    float *dM        ; ALLOCATE (dM       ,  float, catalog->Nmeasure);
-    float *airmass   ; ALLOCATE (airmass  ,  float, catalog->Nmeasure);
-    float *Xccd      ; ALLOCATE (Xccd     ,  float, catalog->Nmeasure);
-    float *Yccd      ; ALLOCATE (Yccd     ,  float, catalog->Nmeasure);
-    float *dt        ; ALLOCATE (dt       ,  float, catalog->Nmeasure);
-    int   *t         ; ALLOCATE (t        ,  int  , catalog->Nmeasure);
-    int   *averef    ; ALLOCATE (averef   ,  int  , catalog->Nmeasure);
-    int   *imageID   ; ALLOCATE (imageID  ,  int  , catalog->Nmeasure);
-    int   *dbFlags   ; ALLOCATE (dbFlags  ,  int  , catalog->Nmeasure);
-    int   *photFlags ; ALLOCATE (photFlags,  int  , catalog->Nmeasure);
-    int   *catID     ; ALLOCATE (catID    ,  int  , catalog->Nmeasure);
-    short *photcode  ; ALLOCATE (photcode ,  short, catalog->Nmeasure);
+    double *R         ; ALLOCATE (R        ,  double, catalog->Nmeasure);
+    double *D         ; ALLOCATE (D        ,  double, catalog->Nmeasure);
+    float  *M         ; ALLOCATE (M        ,  float,  catalog->Nmeasure);
+    float  *Mcal      ; ALLOCATE (Mcal     ,  float,  catalog->Nmeasure);
+    float  *dM        ; ALLOCATE (dM       ,  float,  catalog->Nmeasure);
+    float  *airmass   ; ALLOCATE (airmass  ,  float,  catalog->Nmeasure);
+    float  *Xccd      ; ALLOCATE (Xccd     ,  float,  catalog->Nmeasure);
+    float  *Yccd      ; ALLOCATE (Yccd     ,  float,  catalog->Nmeasure);
+    float  *dt        ; ALLOCATE (dt       ,  float,  catalog->Nmeasure);
+    int    *t         ; ALLOCATE (t        ,  int  ,  catalog->Nmeasure);
+    int    *averef    ; ALLOCATE (averef   ,  int  ,  catalog->Nmeasure);
+    int    *imageID   ; ALLOCATE (imageID  ,  int  ,  catalog->Nmeasure);
+    int    *dbFlags   ; ALLOCATE (dbFlags  ,  int  ,  catalog->Nmeasure);
+    int    *photFlags ; ALLOCATE (photFlags,  int  ,  catalog->Nmeasure);
+    int    *catID     ; ALLOCATE (catID    ,  int  ,  catalog->Nmeasure);
+    short  *photcode  ; ALLOCATE (photcode ,  short,  catalog->Nmeasure);
 
     // assign the storage arrays
     MeasureTiny *measure = catalog->measure;
     for (i = 0; i < catalog->Nmeasure; i++) {
-      dR[i]       = measure[i].dR       ;
-      dD[i]       = measure[i].dD       ;
+      R[i]        = measure[i].R       ;
+      D[i]        = measure[i].D       ;
       M[i]        = measure[i].M        ;
       Mcal[i]     = measure[i].Mcal     ;
@@ -325 +325 @@
 
     // add the columns to the output array
-    gfits_set_bintable_column (&theader, &ftable, "RA_OFF",     dR,        catalog->Nmeasure);
-    gfits_set_bintable_column (&theader, &ftable, "DEC_OFF",    dD,        catalog->Nmeasure);
+    gfits_set_bintable_column (&theader, &ftable, "RA",         R,         catalog->Nmeasure);
+    gfits_set_bintable_column (&theader, &ftable, "DEC",        D,         catalog->Nmeasure);
     gfits_set_bintable_column (&theader, &ftable, "MAG_SYS",    M,         catalog->Nmeasure);
     gfits_set_bintable_column (&theader, &ftable, "MAG_CAL",    Mcal,      catalog->Nmeasure);
@@ -342 +342 @@
     gfits_set_bintable_column (&theader, &ftable, "PHOTCODE",   photcode,  catalog->Nmeasure);
 
-    free (dR      );
-    free (dD      );
+    free (R       );
+    free (D       );
     free (M       );
     free (Mcal    );
@@ -428 +428 @@
     gfits_create_table_header (&theader, "BINTABLE", "SECFILT");
 
-    gfits_define_bintable_column (&theader, "E", "MAG",      "ra offset",                "arcsec", 1.0, 0.0);
-    gfits_define_bintable_column (&theader, "E", "MAG_ERR",  "dec offset",               "arcsec", 1.0, 0.0);
-    gfits_define_bintable_column (&theader, "E", "MAG_CHI",  "magnitude (sys)",           NULL,    1.0, 0.0);
-    gfits_define_bintable_column (&theader, "J", "FLAGS",    "magnitude (cal)",           NULL,    1.0, 0.0);
-    gfits_define_bintable_column (&theader, "I", "NCODE",    "magnitude (err)",           NULL,    1.0, 0.0);
-    gfits_define_bintable_column (&theader, "I", "NUSED",    "airmass",                   NULL,    1.0, 0.0);
-    gfits_define_bintable_column (&theader, "I", "MAG_20",   "ccd x coord",              "pix",    1.0, 0.0);
-    gfits_define_bintable_column (&theader, "I", "MAG_80",   "ccd y coord",              "pix",    1.0, 0.0);
+    gfits_define_bintable_column (&theader, "E", "MAG",      "",              "arcsec", 1.0, 0.0);
+    gfits_define_bintable_column (&theader, "E", "MAG_ERR",  "",              "arcsec", 1.0, 0.0);
+    gfits_define_bintable_column (&theader, "E", "MAG_CHI",  "",              NULL,    1.0, 0.0);
+    gfits_define_bintable_column (&theader, "J", "FLAGS",    "",              NULL,    1.0, 0.0);
+    gfits_define_bintable_column (&theader, "I", "NCODE",    "",              NULL,    1.0, 0.0);
+    gfits_define_bintable_column (&theader, "I", "NUSED",    "",              NULL,    1.0, 0.0);
+    gfits_define_bintable_column (&theader, "E", "MAG_MIN",  "min valid mag", "mag",   1.0, 0.0);
+    gfits_define_bintable_column (&theader, "E", "MAG_MAX",  "max valid mag", "mag",   1.0, 0.0);
 
     // generate the output array that carries the data
@@ -447 +447 @@
     float *M        ; ALLOCATE (M      ,  float, Nsec);
     float *dM       ; ALLOCATE (dM     ,  float, Nsec);
-    float *Xm       ; ALLOCATE (Xm     ,  float, Nsec);
+    float *Mchisq   ; ALLOCATE (Mchisq ,  float, Nsec);
     int   *flags    ; ALLOCATE (flags  ,  int,   Nsec);
     short *Ncode    ; ALLOCATE (Ncode  ,  short, Nsec);
     short *Nused    ; ALLOCATE (Nused  ,  short, Nsec);
-    short *M_20     ; ALLOCATE (M_20   ,  short, Nsec);
-    short *M_80     ; ALLOCATE (M_80   ,  short, Nsec);
+    float *Mmin     ; ALLOCATE (Mmin   ,  float, Nsec);
+    float *Mmax     ; ALLOCATE (Mmax   ,  float, Nsec);
 
     // assign the storage arrays
     SecFilt *secfilt = catalog->secfilt;
     for (i = 0; i < Nsec; i++) {
-      M     [i]       = secfilt[i]. M      ;
-      dM    [i]       = secfilt[i]. dM     ;
-      Xm    [i]       = secfilt[i]. Xm     ;
-      flags [i]       = secfilt[i]. flags  ;
-      Ncode [i]       = secfilt[i]. Ncode  ;
-      Nused [i]       = secfilt[i]. Nused  ;
-      M_20  [i]       = secfilt[i]. M_20   ;
-      M_80  [i]       = secfilt[i]. M_80   ;
+      M     [i]       = secfilt[i].M      ;
+      dM    [i]       = secfilt[i].dM     ;
+      Mchisq[i]       = secfilt[i].Mchisq ;
+      flags [i]       = secfilt[i].flags  ;
+      Ncode [i]       = secfilt[i].Ncode  ;
+      Nused [i]       = secfilt[i].Nused  ;
+      Mmin  [i]       = secfilt[i].Mmin   ;
+      Mmax  [i]       = secfilt[i].Mmax   ;
     }
 
     // add the columns to the output array
-    gfits_set_bintable_column (&theader, &ftable, "MAG",      M    , Nsec);
-    gfits_set_bintable_column (&theader, &ftable, "MAG_ERR",  dM   , Nsec);
-    gfits_set_bintable_column (&theader, &ftable, "MAG_CHI",  Xm   , Nsec);
-    gfits_set_bintable_column (&theader, &ftable, "FLAGS",    flags, Nsec);
-    gfits_set_bintable_column (&theader, &ftable, "NCODE",    Ncode, Nsec);
-    gfits_set_bintable_column (&theader, &ftable, "NUSED",    Nused, Nsec);
-    gfits_set_bintable_column (&theader, &ftable, "MAG_20",   M_20 , Nsec);
-    gfits_set_bintable_column (&theader, &ftable, "MAG_80",   M_80 , Nsec);
+    gfits_set_bintable_column (&theader, &ftable, "MAG",      M     , Nsec);
+    gfits_set_bintable_column (&theader, &ftable, "MAG_ERR",  dM    , Nsec);
+    gfits_set_bintable_column (&theader, &ftable, "MAG_CHI",  Mchisq, Nsec);
+    gfits_set_bintable_column (&theader, &ftable, "FLAGS",    flags , Nsec);
+    gfits_set_bintable_column (&theader, &ftable, "NCODE",    Ncode , Nsec);
+    gfits_set_bintable_column (&theader, &ftable, "NUSED",    Nused , Nsec);
+    gfits_set_bintable_column (&theader, &ftable, "MAG_MIN",  Mmin  , Nsec);
+    gfits_set_bintable_column (&theader, &ftable, "MAG_MAX",  Mmax  , Nsec);
 
     free (M      );
     free (dM     );
-    free (Xm     );
+    free (Mchisq );
     free (flags  );
     free (Ncode  );
     free (Nused  );
-    free (M_20   );
-    free (M_80   );
+    free (Mmin   );
+    free (Mmax   );
 
     gfits_fwrite_Theader (f, &theader);

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

@@ -414 +414 @@
         // 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]);
+        Msys = PhotCatTiny (&catalog[c].measureT[m], MAG_CLASS_PSF);
         if (isnan(Msys)) {
           Nsys++;
@@ -573 +573 @@
       
       n = catalog[c].measureT[m].averef;
-      Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt]);
+      Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF);
       if (isnan(Msys)) {
         Nsys++;
@@ -651 +651 @@
 
       n = catalog[c].measureT[m].averef;
-      Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt]);
+      Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF);
 
       xlist[N] = Xmeas[c][m];

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

@@ -61 +61 @@
   GET_COLUMN (dMcal,             "MCAL_ERR",       float);
   GET_COLUMN (dMagSys,           "MCAL_SYSERR",    float);
-  GET_COLUMN (Xm,                "CHISQ",          float);
   GET_COLUMN (nFitPhotom,        "NFIT",             int);
   GET_COLUMN (flags,             "FLAGS",            int);
   GET_COLUMN (ubercalDist,       "UDIST",            int);
   GET_COLUMN (imageID,           "ID",               int);
+  GET_COLUMN (Xm,                "CHISQ",          short);
 
   // free the memory associated with the FITS files
@@ -121 +121 @@
   gfits_create_table_header (&theader, "BINTABLE", "IMAGE_MAGS");
 
-  gfits_define_bintable_column (&theader, "E", "MCAL",           "cal offset", "magnitudes", 1.0, 0.0);
-  gfits_define_bintable_column (&theader, "E", "MCAL_ERR",       "cal error", "magnitudes", 1.0, 0.0);
-  gfits_define_bintable_column (&theader, "E", "MCAL_SYSERR",    "systematic error", "magnitudes", 1.0, 0.0);
-  gfits_define_bintable_column (&theader, "E", "CHISQ",          "cal chisq", "unitless", 1.0, 0.0);
-  gfits_define_bintable_column (&theader, "J", "NFIT",           "number of fitted stars", "unitless", 1.0, 0.0);
-  gfits_define_bintable_column (&theader, "J", "FLAGS",          "analysis flags", "unitless", 1.0, 0.0);
-  gfits_define_bintable_column (&theader, "J", "UDIST",          "distance to ubercal images", "images", 1.0, 0.0);
-  gfits_define_bintable_column (&theader, "J", "ID",             "image ID", "unitless", 1.0, 0.0);
+  gfits_define_bintable_column (&theader, "E", "MCAL",           "cal offset",                 "magnitudes", 1.0, 0.0);
+  gfits_define_bintable_column (&theader, "E", "MCAL_ERR",       "cal error",                  "magnitudes", 1.0, 0.0);
+  gfits_define_bintable_column (&theader, "E", "MCAL_SYSERR",    "systematic error",           "magnitudes", 1.0, 0.0);
+  gfits_define_bintable_column (&theader, "J", "NFIT",           "number of fitted stars",     "unitless",   1.0, 0.0);
+  gfits_define_bintable_column (&theader, "J", "FLAGS",          "analysis flags",             "unitless",   1.0, 0.0);
+  gfits_define_bintable_column (&theader, "J", "UDIST",          "distance to ubercal images", "images",     1.0, 0.0);
+  gfits_define_bintable_column (&theader, "J", "ID",             "image ID",                   "unitless",   1.0, 0.0);
+  gfits_define_bintable_column (&theader, "I", "CHISQ",          "cal chisq",                  "unitless",   1.0, FT_BZERO_INT16);
 
   // generate the output array that carries the data

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

@@ -531 +531 @@
   liststats_setmode (&stats, STATMODE);
 
-  if (FREEZE_IMAGES) return;
+  // FREEZE_IMAGES only applies to mosaic data (eg, gpc1)
 
   fprintf (stderr, "limiting negative clouds to %f\n", CLOUD_TOLERANCE);
@@ -561 +561 @@
       if (!bad) continue;
     }      
+
+    if (FREEZE_IMAGES && isGPC1chip(image[i].photcode)) continue;
 
     // UBERCAL image: if this is an ubercal image, set minUbercalDist to 0:
@@ -608 +610 @@
 
       n = catalog[c].measureT[m].averef;
-      Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt]);
+      Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF);
       if (isnan(Msys)) {
         Nsys++;
@@ -695 +697 @@
   double MaxOffset, MaxScatter, MedOffset;
 
-  if (FREEZE_IMAGES) return;
+  // if (FREEZE_IMAGES) return;
 
   if (VERBOSE) fprintf (stderr, "marking poor images\n");
@@ -706 +708 @@
   for (i = N = 0; i < Nimage; i++) {
     if (image[i].flags & IMAGE_BAD) continue;
+
+    if (FREEZE_IMAGES && isGPC1chip(image[i].photcode)) continue;
+
     mlist[N] = image[i].Mcal;
     slist[N] = image[i].dMcal;
@@ -797 +802 @@
   Graphdata graphdata;
 
-  if (FREEZE_IMAGES) return;
+  // if (FREEZE_IMAGES) return;
 
   ALLOCATE (xlist, double, Nimage);
@@ -810 +815 @@
   float mindMcal    = +100.0;
   float maxdMcal    = -100.0;
-  for (i = 0; i < Nimage; i++) {
-    Mlist[i] = image[i].Mcal;
-    dlist[i] = image[i].dMcal;
-    xlist[i] = image[i].secz;
+
+  int Nplot = 0;
+
+  for (i = 0; i < Nimage; i++) {
+
+    if (FREEZE_IMAGES && isGPC1chip(image[i].photcode)) continue;
+
+    Mlist[Nplot] = image[i].Mcal;
+    dlist[Nplot] = image[i].dMcal;
+    xlist[Nplot] = image[i].secz;
     minAirmass = MIN (image[i].secz, minAirmass);
     maxAirmass = MAX (image[i].secz, maxAirmass);
@@ -820 +831 @@
     mindMcal = MIN (image[i].dMcal, mindMcal);
     maxdMcal = MAX (image[i].dMcal, maxdMcal);
+    
+    Nplot ++;
   }
 
@@ -852 +865 @@
   for (i = 0; i < NBIN; i++) xlist[i] = 0.00025*i;
   bzero (Mlist, NBIN*sizeof(double));
-  for (i = 0; i < Nimage; i++) {
+
+  for (i = 0; i < Nimage; i++) {
+    if (FREEZE_IMAGES && isGPC1chip(image[i].photcode)) continue;
+
     bin = image[i].dMcal / 0.00025;
     bin = MAX (0, MIN (NBIN - 1, bin));
@@ -877 +893 @@
   StatType stats;
   bzero (&stats, sizeof (StatType));
-  if (FREEZE_IMAGES) return (stats);
+
+  // we no longer blindly apply FREEZE_IMAGES to all images, only to mosaics
+  // if (FREEZE_IMAGES) return (stats);
 
   ALLOCATE (list, double, Nimage);
@@ -885 +903 @@
   for (i = 0; i < Nimage; i++) {
     if (image[i].flags & IMAGE_BAD)  continue;
+
+    if (FREEZE_IMAGES && isGPC1chip(image[i].photcode)) continue;
 
     N = 0;
@@ -920 +940 @@
 
   bzero (&stats, sizeof (StatType));
-  if (FREEZE_IMAGES) return (stats);
+  // if (FREEZE_IMAGES) return (stats);
 
   ALLOCATE (list, double, Nimage);
@@ -929 +949 @@
 
     if (image[i].flags & IMAGE_BAD)  continue;
+
+    if (FREEZE_IMAGES && isGPC1chip(image[i].photcode)) continue;
 
     list[n] = pow (10.0, 0.01*image[i].Xm);
@@ -950 +972 @@
 
   bzero (&stats, sizeof (StatType));
-  if (FREEZE_IMAGES) return (stats);
+  // if (FREEZE_IMAGES) return (stats);
 
   ALLOCATE (list, double, Nimage);
@@ -959 +981 @@
 
     if (image[i].flags & IMAGE_BAD)  continue;
+
+    if (FREEZE_IMAGES && isGPC1chip(image[i].photcode)) continue;
 
     list[n] = image[i].Mcal;
@@ -980 +1004 @@
 
   bzero (&stats, sizeof (StatType));
-  if (FREEZE_IMAGES) return (stats);
+  // if (FREEZE_IMAGES) return (stats);
 
   ALLOCATE (list, double, Nimage);
@@ -989 +1013 @@
 
     if (image[i].flags & IMAGE_BAD)  continue;
+
+    if (FREEZE_IMAGES && isGPC1chip(image[i].photcode)) continue;
 
     list[n] = image[i].dMcal;

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

@@ -126 +126 @@
   gfits_create_table_header (&theader, "BINTABLE", "IMAGE_SUBSET");
 
-  gfits_define_bintable_column (&theader, "E", "MCAL", "zero point offset", "magnitudes", 1.0, 0.0);
-  gfits_define_bintable_column (&theader, "E", "MCAL_ERR", "zero point error", "magnitudes", 1.0, 0.0);
-
-  gfits_define_bintable_column (&theader, "J", "IMAGE_ID", "image ID", NULL, 1.0, 1.0*0x8000);
-  gfits_define_bintable_column (&theader, "J", "PHOTOM_MAP", "map", NULL, 1.0, 1.0*0x8000);
-  gfits_define_bintable_column (&theader, "J", "FLAGS", "flags", NULL, 1.0, 1.0*0x8000);
-
-  gfits_define_bintable_column (&theader, "J", "TESS_ID", "ID", NULL, 1.0, 0.0);
-  gfits_define_bintable_column (&theader, "J", "PROJ_ID", "ID", NULL, 1.0, 0.0);
-  gfits_define_bintable_column (&theader, "J", "SKYCELL_ID", "ID", NULL, 1.0, 0.0);
-
-  gfits_define_bintable_column (&theader, "J", "TZERO", "exposure start", NULL, 1.0, 1.0*0x8000);
-  gfits_define_bintable_column (&theader, "I", "TRATE", "tti rate", NULL, 1.0, 0.0);
-
-  gfits_define_bintable_column (&theader, "I", "UBERCAL_DIST", "ubercal distance", NULL, 1.0, 1.0*0x80);
+  gfits_define_bintable_column (&theader, "E", "MCAL",         "zero point offset", "magnitudes", 1.0, 0.0);
+  gfits_define_bintable_column (&theader, "E", "MCAL_ERR",     "zero point error",  "magnitudes", 1.0, 0.0);
+  gfits_define_bintable_column (&theader, "J", "IMAGE_ID",     "image ID",          NULL,         1.0, FT_BZERO_INT32);
+  gfits_define_bintable_column (&theader, "J", "PHOTOM_MAP",   "map",               NULL,         1.0, FT_BZERO_INT32);
+  gfits_define_bintable_column (&theader, "J", "FLAGS",        "flags",             NULL,         1.0, FT_BZERO_INT32);
+  gfits_define_bintable_column (&theader, "J", "TESS_ID",      "ID",                NULL,         1.0, 0.0);
+  gfits_define_bintable_column (&theader, "J", "PROJ_ID",      "ID",                NULL,         1.0, 0.0);
+  gfits_define_bintable_column (&theader, "J", "SKYCELL_ID",   "ID",                NULL,         1.0, 0.0);
+  gfits_define_bintable_column (&theader, "J", "TZERO",        "exposure start",    NULL,         1.0, FT_BZERO_INT32);
+  gfits_define_bintable_column (&theader, "I", "TRATE",        "tti rate",          NULL,         1.0, FT_BZERO_INT16);
+  gfits_define_bintable_column (&theader, "I", "UBERCAL_DIST", "ubercal distance",  NULL,         1.0, 0.0);
 
   // generate the output array that carries the data
@@ -152 +148 @@
 
   // create intermediate storage arrays
-  ALLOCATE (Mcal,      float,        Nimage);
-  ALLOCATE (dMcal,     float,        Nimage);
+  ALLOCATE (Mcal,      float,          Nimage);
+  ALLOCATE (dMcal,     float,          Nimage);
   ALLOCATE (imageID,   unsigned int,   Nimage);
   ALLOCATE (map,       unsigned int,   Nimage);

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

@@ -60 +60 @@
   GET_COLUMN (M,         "MAG",          float);
   GET_COLUMN (dM,        "MAG_ERR",      float);
-  GET_COLUMN (Xm,        "MAG_CHISQ",    float);
+  GET_COLUMN (Mchisq,    "MAG_CHISQ",    float);
   GET_COLUMN (Nsec,      "NSEC",         int);
   GET_COLUMN (objID,     "OBJ_ID",       int);
@@ -73 +73 @@
   ALLOCATE (meanmags, MeanMag, Nrow);
   for (i = 0; i < Nrow; i++) {
-    meanmags[i].M              = M    [i];
-    meanmags[i].dM             = dM   [i];
-    meanmags[i].Xm             = Xm   [i];
-    meanmags[i].Nsec           = Nsec [i];
-    meanmags[i].objID          = objID[i];
-    meanmags[i].catID          = catID[i];
+    meanmags[i].M              = M      [i];
+    meanmags[i].dM             = dM     [i];
+    meanmags[i].Mchisq         = Mchisq [i];
+    meanmags[i].Nsec           = Nsec   [i];
+    meanmags[i].objID          = objID  [i];
+    meanmags[i].catID          = catID  [i];
   }
   fprintf (stderr, "loaded data for %lld objects (* filters)\n", (long long) Nrow);
 
-  free (M    );
-  free (dM   );
-  free (Xm   );
-  free (Nsec );
-  free (objID);
-  free (catID);
+  free (M     );
+  free (dM    );
+  free (Mchisq);
+  free (Nsec  );
+  free (objID );
+  free (catID );
 
   *nmeanmags = Nrow;
@@ -118 +118 @@
   gfits_define_bintable_column (&theader, "E", "MAG_ERR",   "mean magnitude error", "magnitudes", 1.0, 0.0);
   gfits_define_bintable_column (&theader, "E", "MAG_CHISQ", "mean magnitude chisq", "unitless", 1.0, 0.0);
-  gfits_define_bintable_column (&theader, "J", "NSEC",      "secfilt sequence", NULL, 1.0, 1.0*0x8000);
-  gfits_define_bintable_column (&theader, "J", "OBJ_ID",    "object ID",        NULL, 1.0, 1.0*0x8000);
-  gfits_define_bintable_column (&theader, "J", "CAT_ID",    "catalog ID",       NULL, 1.0, 1.0*0x8000);
+  gfits_define_bintable_column (&theader, "J", "NSEC",      "secfilt sequence", NULL, 1.0, FT_BZERO_INT32);
+  gfits_define_bintable_column (&theader, "J", "OBJ_ID",    "object ID",        NULL, 1.0, FT_BZERO_INT32);
+  gfits_define_bintable_column (&theader, "J", "CAT_ID",    "catalog ID",       NULL, 1.0, FT_BZERO_INT32);
 
   // generate the output array that carries the data
   gfits_create_table (&theader, &ftable);
 
-  float *M, *dM, *Xm;
+  float *M, *dM, *Mchisq;
   int *Nsec;
   unsigned int *objID, *catID;
@@ -132 +132 @@
   ALLOCATE (M,         float,          Nmeanmags);
   ALLOCATE (dM,        float,          Nmeanmags);
-  ALLOCATE (Xm,        float,          Nmeanmags);
+  ALLOCATE (Mchisq,    float,          Nmeanmags);
   ALLOCATE (Nsec,      int,            Nmeanmags);
   ALLOCATE (objID,     unsigned int,   Nmeanmags);
@@ -139 +139 @@
   // assign the storage arrays
   for (i = 0; i < Nmeanmags; i++) {
-    M    [i]   = meanmags[i].M    ;
-    dM   [i]   = meanmags[i].dM   ;
-    Xm   [i]   = meanmags[i].Xm   ;
-    Nsec [i]   = meanmags[i].Nsec ;
-    objID[i]   = meanmags[i].objID;
-    catID[i]   = meanmags[i].catID;
+    M     [i]   = meanmags[i].M     ;
+    dM    [i]   = meanmags[i].dM    ;
+    Mchisq[i]   = meanmags[i].Mchisq;
+    Nsec  [i]   = meanmags[i].Nsec  ;
+    objID [i]   = meanmags[i].objID ;
+    catID [i]   = meanmags[i].catID ;
   }
 
@@ -150 +150 @@
   gfits_set_bintable_column (&theader, &ftable, "MAG",       M,         Nmeanmags);
   gfits_set_bintable_column (&theader, &ftable, "MAG_ERR",   dM,        Nmeanmags);
-  gfits_set_bintable_column (&theader, &ftable, "MAG_CHISQ", Xm,        Nmeanmags);
+  gfits_set_bintable_column (&theader, &ftable, "MAG_CHISQ", Mchisq,    Nmeanmags);
   gfits_set_bintable_column (&theader, &ftable, "NSEC",      Nsec,      Nmeanmags);
   gfits_set_bintable_column (&theader, &ftable, "OBJ_ID",    objID,     Nmeanmags);
   gfits_set_bintable_column (&theader, &ftable, "CAT_ID",    catID,     Nmeanmags);
 
-  free (M    );
-  free (dM   );
-  free (Xm   );
-  free (Nsec );
-  free (objID);
-  free (catID);
+  free (M     );
+  free (dM    );
+  free (Mchisq);
+  free (Nsec  );
+  free (objID );
+  free (catID );
 
   FILE *f = fopen (filename, "w");

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

@@ -1 +1 @@
 # include "relphot.h"
 void plot_setMcal (double *list, int Npts, StatType *stats, float clouds);
+off_t findMosaic (unsigned int *startTimes, off_t Nmosaic, unsigned int start);
 
 // see discussion in ImagesOps.c re: IDX_T
@@ -36 +37 @@
 // MosaicN_image was 'Nimlist'
 // MosaicToImage was 'imlist'
+
+void sort_times (unsigned int *T, int N) {
+
+# define SWAPFUNC(A,B){ unsigned int tmp; \
+  tmp = T[A]; T[A] = T[B]; T[B] = tmp; \
+}
+# define COMPARE(A,B)(T[A] < T[B])
+
+  OHANA_SORT (N, COMPARE, SWAPFUNC);
+
+# undef SWAPFUNC
+# undef COMPARE
+
+}
 
 /* find mosaic frames (unique time periods & photcode name matches mosaic) */
@@ -96 +111 @@
     
     /* 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].Mcal  = 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].Xm    = 0.0;
+    mosaic[Nmosaic].start  = start;
+    mosaic[Nmosaic].stop   = stop;
+    mosaic[Nmosaic].Mcal   = 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].Xm = 0.0;
     mosaic[Nmosaic].flags  = image[i].flags;
-    mosaic[Nmosaic].secz  = image[i].secz;
+    mosaic[Nmosaic].secz   = image[i].secz;
     mosaic[Nmosaic].photcode = GetPhotcodeEquivCodebyCode (image[i].photcode);
 
@@ -136 +151 @@
 }
 
-off_t findMosaic (unsigned int *startTimes, off_t Nmosaic, unsigned int start);
-
-void sort_times (unsigned int *T, int N) {
-
-# define SWAPFUNC(A,B){ unsigned int tmp; \
-  tmp = T[A]; T[A] = T[B]; T[B] = tmp; \
-}
-# define COMPARE(A,B)(T[A] < T[B])
-
-  OHANA_SORT (N, COMPARE, SWAPFUNC);
-
-# undef SWAPFUNC
-# undef COMPARE
-
-}
-
 /* find mosaic frames (unique time periods) (NOTE : we do NOT require matching photcodes...)
    this function will also identify the images NOT in the subset which belong to a selected mosaic
@@ -157 +156 @@
 void initMosaics (Image *subset, off_t Nsubset, Image *image, char *inSubset, off_t Nimage) {
 
-  off_t i, j, status, found, NMOSAIC, *MosaicN_IMAGE;
+  off_t i, j, found, NMOSAIC, *MosaicN_IMAGE;
   unsigned int start, stop, *startTimes, *startTimesMosaic;
-  char *pname;
 
   if (!MOSAIC_ZEROPT) return;
@@ -172 +170 @@
   // generate a list of all subset image start times
   ALLOCATE (startTimes, unsigned int, Nsubset);
+  int Nmoschip = 0;
   for (i = 0; i < Nsubset; i++) {
-    startTimes[i] = subset[i].tzero;
-  }
-  sort_times (startTimes, Nsubset);
+    if (!isGPC1chip(subset[i].photcode)) continue;
+    startTimes[Nmoschip] = subset[i].tzero;
+    Nmoschip ++;
+  }
+  sort_times (startTimes, Nmoschip);
   MARKTIME("create array of all image obstimes: %f sec\n", dtime);
   
@@ -184 +185 @@
 
   // generate a list of the unique start times (these define the mosaics)
-  for (i = 0; i < Nsubset; i++) {
-    if (startTimes[i] < startTimesMosaic[Nmosaic]) {
-      fprintf (stderr, "error?\n");
-      abort();
-    }
+  for (i = 0; i < Nmoschip; i++) {
+    myAssert (startTimes[i] >= startTimesMosaic[Nmosaic], "times out of order?");
     if (startTimes[i] == startTimesMosaic[Nmosaic]) continue;
     Nmosaic ++;
@@ -210 +208 @@
   for (i = 0; i < Nmosaic; i++) {
     /* a new mosaic, define ranges */
-    mosaic[i].start = startTimesMosaic[i];
-    mosaic[i].stop  = 0;
-    mosaic[i].Mcal  = 0.0;
-    mosaic[i].dMcal = 0.0;
-    mosaic[i].dMsys = 0.0;
-    mosaic[i].Xm    = 0.0;
-    mosaic[i].flags = 0;
-    mosaic[i].secz  = NAN;
+    mosaic[i].start    = startTimesMosaic[i];
+    mosaic[i].stop     = 0;
+    mosaic[i].Mcal     = 0.0;
+    mosaic[i].dMcal    = 0.0;
+    mosaic[i].dMsys    = 0.0;
+    mosaic[i].Xm       = 0.0;
+    mosaic[i].flags    = 0;
+    mosaic[i].secz     = NAN;
     mosaic[i].photcode = 0;
-    mosaic[i].skipCal = FALSE;
+    mosaic[i].skipCal  = FALSE;
     
     memset (&mosaic[i].coords, 0, sizeof(Coords));
@@ -238 +236 @@
     }
     
-    /* select valid mosaic images by photcode */
-    pname = GetPhotcodeNamebyCode (image[i].photcode);
-    if (!pname) continue;
-    status = strncmp (pname, MOSAICNAME, strlen (MOSAICNAME));
-    if (status) continue;
+    if (!isGPC1chip(image[i].photcode)) continue;
 
     /* set image time range */
@@ -262 +256 @@
 
   // assign each image to a mosaic
+  int Nsimple = 0;
   for (i = 0; i < Nsubset; i++) {
     ImageToMosaic[i] = -1;
 
-    /* select valid mosaic images by photcode */
-    pname = GetPhotcodeNamebyCode (subset[i].photcode);
-    status = strncmp (pname, MOSAICNAME, strlen (MOSAICNAME));
-    if (status) continue;
+    if (!isGPC1chip(subset[i].photcode)) {
+      Nsimple ++;
+      continue;
+    }
 
     start = subset[i].tzero;
@@ -299 +294 @@
       abort();
     }
-    mosaic[j].stop  = stop;
-    mosaic[j].Mcal  = 0.0;
-    mosaic[j].dMcal = 0.0;
-    mosaic[j].Xm    = 0.0;
-    mosaic[j].dMsys = subset[i].flags;
-    mosaic[j].flags = subset[i].flags;
-    mosaic[j].secz  = subset[i].secz;
+    mosaic[j].stop     = stop;
+    mosaic[j].Mcal     = 0.0;
+    mosaic[j].dMcal    = 0.0;
+    mosaic[j].Xm       = 0.0;
+    mosaic[j].dMsys    = subset[i].flags;
+    mosaic[j].flags    = subset[i].flags;
+    mosaic[j].secz     = subset[i].secz;
     mosaic[j].photcode = GetPhotcodeEquivCodebyCode (subset[i].photcode);
   }
@@ -316 +311 @@
   initMosaicGrid (subset, Nsubset);
 
-  fprintf (stderr, "matched %d images to %d mosaics\n", (int) Nsubset, (int) Nmosaic);
+  fprintf (stderr, "matched %d images to %d mosaics, %d simple chips not matched to mosaics\n", (int) (Nsubset - Nsimple), (int) Nmosaic, (int) Nsimple);
   return;
 }
 
-/* find mosaic frames (unique time periods) (NOTE : we do NOT require matching photcodes but we 
-   match images by MOSAICNAME.  this last point is weak: it forces a single camera at a time.
-   we can extend the logic to multiple cameras if we make list of MOSAICNAMES (better to assign a camera ID)
+/* find mosaic frames (unique time periods) (NOTE : require gpc1 chips, which is pretty limiting)
+   if mergeMcal is TRUE, <image.Mcal> values will be saved on Mosaic.Mcal
  */
-void makeMosaics (Image *image, off_t Nimage) {
-
-  off_t i, j, status, found, NMOSAIC, *MosaicN_IMAGE;
+void makeMosaics (Image *image, off_t Nimage, int mergeMcal) {
+
+  off_t i, j, found, NMOSAIC, *MosaicN_IMAGE;
   unsigned int start, stop, *startTimes, *startTimesMosaic;
-  char *pname;
 
   if (!MOSAIC_ZEROPT) return;
@@ -341 +334 @@
   // generate a list of all image start times
   ALLOCATE (startTimes, unsigned int, Nimage);
+  int Nmoschip = 0;
   for (i = 0; i < Nimage; i++) {
-    startTimes[i] = image[i].tzero;
-  }
-  sort_times (startTimes, Nimage);
+    if (!isGPC1chip(image[i].photcode)) continue;
+    startTimes[Nmoschip] = image[i].tzero;
+    Nmoschip ++;
+  }
+  sort_times (startTimes, Nmoschip);
   MARKTIME("create array of all image obstimes: %f sec\n", dtime);
   
@@ -353 +349 @@
 
   // generate a list of the unique start times (these define the mosaics)
-  for (i = 0; i < Nimage; i++) {
-    if (startTimes[i] < startTimesMosaic[Nmosaic]) {
-      fprintf (stderr, "error?\n");
-      abort();
-    }
+  for (i = 0; i < Nmoschip; i++) {
+    myAssert (startTimes[i] >= startTimesMosaic[Nmosaic], "times out of order?");
     if (startTimes[i] == startTimesMosaic[Nmosaic]) continue;
     Nmosaic ++;
@@ -379 +372 @@
   // init the mosaic array values
   for (i = 0; i < Nmosaic; i++) {
-    mosaic[i].start = startTimesMosaic[i];
-    mosaic[i].stop  = 0;
-    mosaic[i].Mcal  = 0.0;
-    mosaic[i].dMcal = 0.0;
-    mosaic[i].dMsys = 0.0;
-    mosaic[i].Xm    = 0.0;
-    mosaic[i].flags = 0;
-    mosaic[i].secz  = NAN;
+    mosaic[i].start    = startTimesMosaic[i];
+    mosaic[i].stop     = 0;
+    mosaic[i].Mcal     = 0.0;
+    mosaic[i].dMcal    = 0.0;
+    mosaic[i].dMsys    = 0.0;
+    mosaic[i].Xm       = 0.0;
+    mosaic[i].flags    = 0;
+    mosaic[i].secz     = NAN;
     mosaic[i].photcode = 0;
-    mosaic[i].skipCal = FALSE;
+    mosaic[i].skipCal  = FALSE;
     
     memset (&mosaic[i].coords, 0, sizeof(Coords));
@@ -401 +394 @@
 
   // assign each image to a mosaic
+  int Nsimple = 0;
   for (i = 0; i < Nimage; i++) {
     ImageToMosaic[i] = -1;
 
-    /* select valid mosaic images by photcode */
-    pname = GetPhotcodeNamebyCode (image[i].photcode);
-    if (!pname) continue;
-
-    status = strncmp (pname, MOSAICNAME, strlen (MOSAICNAME));
-    if (status) continue;
+    if (!isGPC1chip(image[i].photcode)) {
+      Nsimple ++;
+      continue;
+    }
 
     start = image[i].tzero;
@@ -440 +432 @@
       abort();
     }
-    mosaic[j].stop  = stop;
-    mosaic[j].Mcal  = 0.0;
-    mosaic[j].dMcal = 0.0;
-    mosaic[j].Xm    = 0.0;
-    mosaic[j].dMsys = image[i].flags;
-    mosaic[j].flags = image[i].flags;
-    mosaic[j].secz  = image[i].secz;
+    mosaic[j].stop     = stop;
+    mosaic[j].Mcal     = 0.0;
+    mosaic[j].dMcal    = 0.0;
+    mosaic[j].Xm       = 0.0;
+    mosaic[j].dMsys    = image[i].flags;
+    mosaic[j].flags    = image[i].flags;
+    mosaic[j].secz     = image[i].secz;
     mosaic[j].photcode = GetPhotcodeEquivCodebyCode (image[i].photcode);
   }
@@ -456 +448 @@
   free (startTimesMosaic);
 
+  if (mergeMcal) {
+    initMosaicGrid (image, Nimage);
+  }
+
   fprintf (stderr, "matched %d images to %d mosaics\n", (int) Nimage, (int) Nmosaic);
-
   return;
 }
@@ -513 +508 @@
   off_t i, j, m, NX, NY, NC, Nc;
   double R, D, Rmid, Dmid;
-  double Mcal, dMcal, Xm;
   double *Rc, *Dc;
 
@@ -523 +517 @@
     Nc = 0;
     Rmid = Dmid = NAN;
-    Mcal = dMcal = Xm = 0;
     for (j = 0; j < MosaicN_Image[i]; j++) {
       m = MosaicToImage[i][j];
@@ -563 +556 @@
         REALLOCATE (Dc, double, NC);
       }
-
-      Mcal  += image[m].Mcal;
-      dMcal += image[m].dMcal;
-      Xm    += image[m].Xm;
-
-      // for ubercal images, we (elsewhere) keep Mcal frozen
-
-      /* we are using mosaic.Mcal, not image.Mcal. reset image.Mcal */
-      image[m].Mcal  = 0.0;
-      image[m].dMcal = NAN;
-      image[m].Xm    = NAN_S_SHORT;
     }
 
@@ -617 +599 @@
     mosaic[i].coords.pc2_1  = 0.0;
 
-    mosaic[i].Mcal  = Mcal / MosaicN_Image[i];
-    mosaic[i].dMcal = dMcal / MosaicN_Image[i];
-    mosaic[i].Xm    = Xm / MosaicN_Image[i];
+    mosaic[i].Mcal   = 0.0;
+    mosaic[i].dMcal  = 0.0;
+    mosaic[i].Xm     = 0.0;
   }
   return;
@@ -635 +617 @@
   double R, D, Rmin, Rmax, Dmin, Dmax;
   double Mcal, dMcal, Xm;
+
+  fprintf (stderr, "*** moving Mcal from image.Mcal to mosaic.Mcal ***\n");
 
   dXmax = dYmax = 0.0;
@@ -667 +651 @@
       Dmin = MIN (Dmin, D);
       Dmax = MAX (Dmax, D);
+
+      /* 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
+
       Mcal  += image[m].Mcal;
       dMcal += image[m].dMcal;
       Xm    += image[m].Xm;
 
-      // 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
-
-      /* we are using mosaic.Mcal, not image.Mcal. reset image.Mcal */
-      image[m].Mcal  = 0.0;
-      image[m].dMcal = NAN;
-      image[m].Xm    = NAN_S_SHORT;
+      image[m].Mcal   = 0.0;
+      image[m].dMcal  = NAN;
+      image[m].Xm     = NAN_S_SHORT;
     }
     dS /= MosaicN_Image[i];
@@ -697 +681 @@
     RD_to_XY (&dX, &dY, Rmax, Dmax, &mosaic[i].coords);
 
-    mosaic[i].Mcal  = Mcal / MosaicN_Image[i];
-    mosaic[i].dMcal = dMcal / MosaicN_Image[i];
-    mosaic[i].Xm    = Xm / MosaicN_Image[i];
+    mosaic[i].Mcal   = Mcal / MosaicN_Image[i];
+    mosaic[i].dMcal  = dMcal / MosaicN_Image[i];
+    mosaic[i].Xm     = Xm / MosaicN_Image[i];
   }
   if (!USE_GRID) return;
@@ -709 +693 @@
 }
 
-// XXX : what about mosaics with skipCal == TRUE?
 void setMcalFinal () {
 
@@ -719 +702 @@
   image = getimages (&Nimage, NULL);
 
-  // XXX I think this is OK in the ubercal context, but probably need to skip UBERCAL
-  // images? (no need to update them)
+  fprintf (stderr, "*** return Mcal from mosaic.Mcal to image.Mcal ***\n");
 
   // copy the mosaic results to the images.  set the mosaic Mcal to 0.0 since we have moved its
@@ -840 +822 @@
   mosID = ImageToMosaic[idx];
   if (mosID < 0) {
-    Image *image = getimage(idx);
-    fprintf (stderr, "unmatched image %s\n", image[0].name);
+    // Image *image = getimage(idx);
+    // fprintf (stderr, "unmatched image %s\n", image[0].name);
+    // skip measurements from simple chips (not matched to a mosaic by definition)
     return;
   }
@@ -882 +865 @@
   if (!MOSAIC_ZEROPT) return (0);
 
+  // unassigned measurements belong to simple chips
   i = MeasureToMosaic[cat][meas];
-  if (i == -1) return (NAN);
+  if (i == -1) return (0.0);
 
   if (mosaic[i].flags & IMAGE_BAD) return (NAN);  
@@ -1127 +1111 @@
 
       off_t n = catalog[c].measureT[m].averef;
-      Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt]);
+      Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF);
 
       float delta = Msys - Mrel - Mcal - Mgrid + Mflat;
@@ -1164 +1148 @@
 
     off_t n = catalog[c].measureT[m].averef;
-    Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt]);
+    Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF);
     if (isnan(Msys)) {
       info->Nsys++;
@@ -1239 +1223 @@
 
   if (testImage) {
-    fprintf (stderr, "%f %f  :  %d %f\n", myMosaic[0].Mcal, myMosaic[0].dMsys, myMosaic[0].Xm, pow(10.0, 0.01*myMosaic[0].Xm));
+    fprintf (stderr, "%f %f  :  %f\n", myMosaic[0].Mcal, myMosaic[0].dMsys, pow (10.0, 0.01*myMosaic[0].Xm));
   }
 
@@ -1592 +1576 @@
     if (mosaic[i].flags & IMAGE_BAD) continue;
     if (mosaic[i].skipCal) continue;
+    if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue;
+
     list[n] = mosaic[i].Mcal;
     dlist[n] = 1;
@@ -1622 +1608 @@
     if (mosaic[i].flags & IMAGE_BAD) continue;
     if (mosaic[i].skipCal) continue;
+    if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue;
+
     list[n] = mosaic[i].dMcal;
     dlist[n] = 1;
@@ -1653 +1641 @@
     if (mosaic[i].flags & IMAGE_BAD)  continue;
     if (mosaic[i].skipCal) continue;
+    if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue;
 
     N = 0;
@@ -1699 +1688 @@
     if (mosaic[i].flags & IMAGE_BAD) continue;
     if (mosaic[i].skipCal) continue;
-    list[n] = pow(10.0, 0.01*mosaic[i].Xm);
+    if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue;
+
+    list[n]  = pow (10.0, 0.01*mosaic[i].Xm);
     dlist[n] = 1;
     n++;
@@ -1782 +1773 @@
 void plot_mosaic_fields (Catalog *catalog) {
 
-  off_t i, j, m, c, N, ave, Nimage;
+  off_t i, j, m, c, N, Nimage;
   double *xlist, *ylist;
-  // double Xmin, Xmax, Ymin, Ymax;
   char string[64];
   Graphdata graphdata;
@@ -1803 +1793 @@
   for (i = 0; i < Nmosaic; i++) {
     N = 0;
-    // Xmin = Ymin = +360.0;
-    // Xmax = Ymax = -360.0;
     for (j = 0; j < N_onMosaic[i]; j++) {
       
@@ -1812 +1800 @@
       if (catalog[c].measureT[m].dbFlags & (ID_MEAS_AREA | ID_MEAS_NOCAL)) continue;
 
-      ave = catalog[c].measureT[m].averef;
-      xlist[N] = catalog[c].averageT[ave].R - catalog[c].measureT[m].dR / 3600.0;
-      ylist[N] = catalog[c].averageT[ave].D - catalog[c].measureT[m].dD / 3600.0;
+      // ave = catalog[c].measureT[m].averef;
+      xlist[N] = catalog[c].measureT[m].R;
+      ylist[N] = catalog[c].measureT[m].D;
       N++;
     }

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

@@ -23 +23 @@
   SetProtect (TRUE);
   if (db) gfits_db_close (db);
-  fprintf (stderr, "ERROR: addstar halted\n");
+  fprintf (stderr, "ERROR: relphot halted\n");
   exit (1);
 }

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

@@ -153 +153 @@
   int i;
 
+  // fprintf (stderr, "this version of setMrel is invalid for now\n");
+  // exit (1);
+        
   int Nsecfilt = GetPhotcodeNsecfilt ();
 
@@ -407 +410 @@
 
       /* star/photcodes already calibrated */
-      if ( isSetMrelFinal && catalog[Nc].found[Nsecfilt*j+Nsec]) continue;  
+      if ( isSetMrelFinal && catalog[Nc].found_t[Nsecfilt*j+Nsec]) continue;  
       
       // skip bad stars
@@ -467 +470 @@
 
         // skip some absurd values NAN, < 0.0, > 30.0
-        Msys = PhotSysTiny (&catalog[Nc].measureT[m], &catalog[Nc].averageT[j], &catalog[Nc].secfilt[j*Nsecfilt]);
+        Msys = PhotSysTiny (&catalog[Nc].measureT[m], &catalog[Nc].averageT[j], &catalog[Nc].secfilt[j*Nsecfilt], MAG_CLASS_PSF);
         if (isnan(Msys)) SKIP_THIS_MEAS(Nsys);
         if (Msys <  0.0) SKIP_THIS_MEAS(Nsys);
@@ -670 +673 @@
       liststats (list, dlist, wlist, N, &stats);
 
-      catalog[Nc].secfilt[Nsecfilt*j+Nsec].M  = stats.mean;
-      catalog[Nc].secfilt[Nsecfilt*j+Nsec].dM = stats.error;
-      catalog[Nc].secfilt[Nsecfilt*j+Nsec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq + 1e-4) : NAN_S_SHORT;
+      catalog[Nc].secfilt[Nsecfilt*j+Nsec].M      = stats.mean;
+      catalog[Nc].secfilt[Nsecfilt*j+Nsec].dM     = stats.error;
+      catalog[Nc].secfilt[Nsecfilt*j+Nsec].Mchisq = (stats.Nmeas > 1) ? stats.chisq : NAN;
 
       // when running -averages, we have no information about the images, so we cannot set this
@@ -680 +683 @@
 
       if (isSetMrelFinal) {
-        catalog[Nc].found[Nsecfilt*j+Nsec] = TRUE;
+        catalog[Nc].found_t[Nsecfilt*j+Nsec] = TRUE;
 
         catalog[Nc].secfilt[Nsecfilt*j+Nsec].Mstdev = 1000.0*stats.sigma; // Mstdev is in millimags (not enough space for more precision)
@@ -950 +953 @@
 
   int i, j, Ndel, Nave, Ntot, mark, Ns, Nscat, Nchi, Nnan;
-  float dM, Xm;
-  double Chisq, MaxScatter, MaxChisq;
+  float dM;
+  double MaxScatter, MaxChisq;
   double *xlist, *slist, *dlist;
 
@@ -978 +981 @@
       for (j = 0; j < catalog[i].Naverage; j++) {
         if ( catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD ) continue;
-        Xm = catalog[i].secfilt[Nsecfilt*j+Nsec].Xm;
-        if (Xm == -1) continue;
-        Chisq = pow (10.0, 0.01*Xm);
-        xlist[Ntot] = Chisq;
+        float Mchisq = catalog[i].secfilt[Nsecfilt*j+Nsec].Mchisq;
+        if (isnan(Mchisq)) continue;
+        xlist[Ntot] = Mchisq;
         slist[Ntot] = catalog[i].secfilt[Nsecfilt*j+Nsec].dM;
         dlist[Ntot] = 1;
@@ -999 +1001 @@
       for (j = 0; j < catalog[i].Naverage; j++) {
         dM = catalog[i].secfilt[Nsecfilt*j+Nsec].dM;
-        Xm = catalog[i].secfilt[Nsecfilt*j+Nsec].Xm;
-        Chisq = pow (10.0, 0.01*Xm);
-        mark = (dM > MaxScatter) || (Xm == NAN_S_SHORT) || (Chisq > MaxChisq);
+        float Mchisq = catalog[i].secfilt[Nsecfilt*j+Nsec].Mchisq;
+        mark = (dM > MaxScatter) || (isnan(Mchisq)) || (Mchisq > MaxChisq);
         if (mark) {
           catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_STAR_POOR;
           Ndel ++;
-          if (dM > MaxScatter) { Nscat ++; }
-          if (Xm == NAN_S_SHORT) { Nnan ++; }
-          if (Chisq > MaxChisq) { Nchi ++; }
+          if (dM > MaxScatter)   { Nscat ++; }
+          if (isnan(Mchisq))     { Nnan ++; }
+          if (Mchisq > MaxChisq) { Nchi ++; }
         } else {
           catalog[i].secfilt[Nsecfilt*j+Nsec].flags &= ~ID_STAR_POOR;
@@ -1072 +1073 @@
         
         /* on final processing, skip stars already measured */
-        if (final && catalog[i].found[Nsecfilt*j + Nsec]) continue;  
+        if (final && catalog[i].found_t[Nsecfilt*j + Nsec]) continue;  
 
         /* skip bad stars to prevent them from becoming good (on inner sample) */
@@ -1093 +1094 @@
           if (isnan(Mgrid)) { Ngrid ++; continue; }
 
-          Msys = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt]);
+          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);
@@ -1138 +1139 @@
           if (isnan(Mgrid)) continue;
 
-          Msys = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt]);
+          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);
@@ -1233 +1234 @@
 
   off_t j, Ntot;
-  int i, n, Xm;
+  int i, n;
   double *list, *dlist;
   StatType stats;
@@ -1254 +1255 @@
       if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;  
         
-      Xm = catalog[i].secfilt[Nsecfilt*j+Nsec].Xm;
-      if (Xm == NAN_S_SHORT) continue;
-      list[n] = pow (10.0, 0.01*Xm);
+      float Mchisq = catalog[i].secfilt[Nsecfilt*j+Nsec].Mchisq;
+      if (isnan(Mchisq)) continue;
+      list[n] = Mchisq;
       dlist[n] = 1;
       n++;
@@ -1377 +1378 @@
         if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;
         xlist[N] = catalog[i].secfilt[Nsecfilt*j+Nsec].M;
-        value    = catalog[i].secfilt[Nsecfilt*j+Nsec].Xm;
-        if (value == NAN_S_SHORT) continue;
-        ylist[N] = 0.01*value;
+        value    = catalog[i].secfilt[Nsecfilt*j+Nsec].Mchisq;
+        if (isnan(value)) continue;
+        ylist[N] = value;
         N++;
       }

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

@@ -227 +227 @@
   }
 
-  SHOW_PARAMS = FALSE;
+  SHOW_PARAMS = TRUE;
   if ((N = get_argument (argc, argv, "-params"))) {
     remove_argument (N, &argc, argv);
-    SHOW_PARAMS = TRUE;
+    SHOW_PARAMS = FALSE;
   }
 

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

@@ -25 +25 @@
 
   if (MOSAIC_ZEROPT) {
-    makeMosaics (image, Nimage);
-
-    // center coords and Mcal, dMcal, Xm for the mosaics
+    makeMosaics (image, Nimage, FALSE);
+
+    // center coords and Mcal, dMcal, Mchisq for the mosaics
     setMosaicCenters (image, Nimage);
     MARKTIME("set mosaic coordinates and Mcal values: %f sec\n", dtime);

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

@@ -23 +23 @@
   Nmeasure = Naverage = 0;
 
+  int NmMax = 0;
+  int NmMin = catalog[0].Nmeasure;
+
   Ncode = Ntime = Ndophot = Nmag = Nsigma = Nimag = Nfew = Npsfqf = Ngalaxy = Nnan = Nbad = Npoor = 0;
 
@@ -55 +58 @@
     if (RESET) {
       int Ns;
+
+# if (0)      
       DVOSecfiltFlags secfiltBits = 
         ID_SECF_STAR_FEW | 
@@ -66 +71 @@
         ID_SECF_USE_UBERCAL |
         ID_SECF_OBJ_EXT;
+# endif
+
       for (Ns = 0; Ns < Nphotcodes; Ns++) {
 
@@ -71 +78 @@
         int Nsec = GetPhotcodeNsec(thisCode);
 
-        subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].M      = NAN;
-        subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Map    = NAN;
-        subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].dM     = NAN;
-        subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Mstdev = NAN_S_SHORT;
-        subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Xm     = NAN;
-        subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].M_20   = NAN_S_SHORT;
-        subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].M_80   = NAN_S_SHORT;
-        subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Ncode  = 0;
-        subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Nused  = 0;
-        subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].flags &= ~secfiltBits;
-        subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].ubercalDist = 1000;
+        fprintf (stderr, "really use dvo_secfilt_init?");
+        dvo_secfilt_init (&subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec]);
+
+        // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].M      = NAN;
+        // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Map    = NAN;
+        // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].dM     = NAN;
+        // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Mstdev = NAN_S_SHORT;
+        // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Mchisq = NAN;
+        // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].M_20   = NAN_S_SHORT;
+        // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].M_80   = NAN_S_SHORT;
+        // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Ncode  = 0;
+        // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].Nused  = 0;
+        // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].flags &= ~secfiltBits;
+        // subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].ubercalDist = 1000;
       }
     }
@@ -136 +146 @@
 
       /* select measurements by mag limit */
-      mag = PhotCat (&catalog[0].measure[offset]);
+      mag = PhotCat (&catalog[0].measure[offset], MAG_CLASS_PSF);
       if (mag > MAG_LIM) { Nmag ++; continue; } 
 
@@ -144 +154 @@
       /* select measurements by mag limit */
       if (ImagSelect) {
-        mag = PhotInst (&catalog[0].measure[offset]);
+        mag = PhotInst (&catalog[0].measure[offset], MAG_CLASS_PSF);
         if (mag < ImagMin) { Nimag ++; continue; }
         if (mag > ImagMax) { Nimag ++; continue; }
@@ -180 +190 @@
       }
       myAssert (Nmeasure < NMEASURE, "realloc failure");
-    }
+    } // end of catalog.average.Nmeasure loop
 
     // skip object if it is likely to be a galaxy
@@ -205 +215 @@
       continue; 
     }
+    NmMin = MIN (Nm, NmMin);
+    NmMax = MAX (Nm, NmMax);
 
     // for w-band photometry (& other cases?) convert gri photometry
@@ -235 +247 @@
   subcatalog[0].Nsecf_mem = Naverage * catalog[0].Nsecfilt;
 
+  if (VERBOSE) {
+    fprintf (stderr, "using "OFF_T_FMT" stars ("OFF_T_FMT" measures) of "OFF_T_FMT" ("OFF_T_FMT" measures) for catalog %s, %d < Nm < %d\n", 
+             subcatalog[0].Naverage,  subcatalog[0].Nmeasure,  catalog[0].Naverage, catalog[0].Nmeasure, catalog[0].filename, NmMin, NmMax);
+    fprintf (stderr, "rejections: %d stars have too few measures:\n   %d code, %d time, %d dophot, %d mag, %d sigma, %d imag, %d psfqf, %d Nnan, %d galaxies, %d bad, %d poor\n", 
+             Nfew, Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Npsfqf, Nnan, Ngalaxy, Nbad, Npoor);
+  }
+
   // limit the total number of stars in the catalog
   if (MaxDensityUse) {
@@ -240 +259 @@
   }
 
-  if (VERBOSE) {
-    fprintf (stderr, "using "OFF_T_FMT" stars ("OFF_T_FMT" measures) of "OFF_T_FMT" for catalog %s\n", 
-             subcatalog[0].Naverage,  subcatalog[0].Nmeasure,  i, catalog[0].filename);
-    fprintf (stderr, "rejections: %d code, %d time, %d dophot, %d mag, %d sigma, %d imag, %d few, %d psfqf, %d Nnan, %d galaxies, %d bad, %d poor\n", 
-             Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Nfew, Npsfqf, Nnan, Ngalaxy, Nbad, Npoor);
-  }
   return (TRUE);
 }

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

@@ -36 +36 @@
 
   if (SHOW_PARAMS) {
-    fprintf (stderr, "current parameter settings:\n");
+    int Ns;
+    fprintf (stderr, "subset selection criteria:\n");
+    fprintf (stderr, "  photcodes ");
+    for (Ns = 0; Ns < Nphotcodes; Ns++) {
+      if (Ns == Nphotcodes - 1) {
+        fprintf (stderr, "%s\n", photcodes[Ns][0].name);
+      } else {
+        fprintf (stderr, "%s, ", photcodes[Ns][0].name);
+      }
+    }
     if (TimeSelect) {
       fprintf (stderr, "TimeSelect: TRUE (%s - %s)\n", ohana_sec_to_date (TSTART), ohana_sec_to_date (TSTOP));
@@ -42 +51 @@
       fprintf (stderr, "TimeSelect: FALSE\n");
     }
+    if (DophotSelect) {
+      fprintf (stderr, "DophotSelect: TRUE (%d)\n", DophotValue);
+    } else {
+      fprintf (stderr, "DophotSelect: FALSE\n");
+    }
+    fprintf (stderr, "PSF_QF limit: 0.85 (hardwired)\n");
+
+    // fprintf (stderr, "Photom Bad Mask: 0x%08x, Photom Poor Mask: 0x%08x\n");
+
+    fprintf (stderr, "MAG_LIM: %f, SIGMA_LIM: %f\n", MAG_LIM, SIGMA_LIM);
+    fprintf (stderr, "INST_MAG_MIN: %f, INST_MAG_MAX: %f\n", ImagMin, ImagMax);
+
+    fprintf (stderr, "STAR_TOOFEW: %d\n", STAR_TOOFEW);
+
     fprintf (stderr, "VERBOSE: %d, PLOTSTUFF: %d\n", VERBOSE, PLOTSTUFF);
     fprintf (stderr, "GRID_X: %d, GRID_Y: %d, BINNING: %d == Nmx: %d, Nmy: %d\n", 
@@ -49 +72 @@
              (int)(RELPHOT_GRID_X/RELPHOT_GRID_BINNING), (int)(RELPHOT_GRID_Y/RELPHOT_GRID_BINNING));
 
-    fprintf (stderr, "MAG_LIM                %lf\n", MAG_LIM);
     fprintf (stderr, "STAR_SCATTER           %lf\n", STAR_SCATTER);
     fprintf (stderr, "IMAGE_SCATTER          %lf\n", IMAGE_SCATTER);
@@ -56 +78 @@
     fprintf (stderr, "GSCFILE                %s\n",  GSCFILE);
     fprintf (stderr, "CATDIR                 %s\n",  CATDIR);
-    exit (0);
   }
 

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

@@ -55 +55 @@
 
     char command[1024];
-    snprintf (command, 1024, "relphot %s -parallel-images %s -region-hosts %s -region-hostID %d -D CATDIR %s -region %f %f %f %f -statmode %s -D CAMERA %s -D STAR_TOOFEW %d -minerror %f", 
-              PhotcodeList, filename, REGION_FILE, host->hostID, CATDIR, host->RminCat, host->RmaxCat, host->DminCat, host->DmaxCat, STATMODE, CAMERA, STAR_TOOFEW, MIN_ERROR);
+    snprintf (command, 1024, "relphot %s", PhotcodeList);
+    strextend (command, "-parallel-images %s", filename);
+    strextend (command, "-region-hosts %s", REGION_FILE);
+    strextend (command, "-region-hostID %d", host->hostID);
+    strextend (command, "-D CATDIR %s", CATDIR);
+    strextend (command, "-region %f %f %f %f", host->RminCat, host->RmaxCat, host->DminCat, host->DmaxCat);
+    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);
+    strextend (command, "-cloud-limit %f", CLOUD_TOLERANCE);
 
     if (VERBOSE)             strextend (command, "-v");

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

@@ -12 +12 @@
 // if this function is called with a specified hostID, then only the fraction of the
 // database hosted by that hostID is loaded
-Catalog *load_catalogs (SkyList *skylist, int *Ncatalog, int hostID, char *hostpath) {
+Catalog *load_catalogs (SkyList *skylist, int *Ncatalog, int hostID, char *hostpath, char *syncfile) {
 
   off_t i, Nmeas, Nstar, Nmeas_total, Nstar_total;
@@ -19 +19 @@
   // XXX need to decide how to determine PARALLEL mode...
   if (PARALLEL && !hostID) {
-    catalog = load_catalogs_parallel (skylist, Ncatalog);
+    catalog = load_catalogs_parallel (skylist, Ncatalog, syncfile);
     return catalog;
   }
@@ -92 +92 @@
   // in regular relphot, we shutdown here; in relphot_client, we generate and return an empty table (for consistency)
 
+  // if we are running with parallel_images but not a parallel database, we need to
+  // release the lock so the next image host can proceed
+  if (!hostID && syncfile) {
+    update_sync_file (syncfile, 1);
+  }
+
   // XXX consider only returning the populated catalogs
   *Ncatalog = skylist[0].Nregions;
@@ -105 +111 @@
 // CATDIR is supplied globally
 # define DEBUG 1
-Catalog *load_catalogs_parallel (SkyList *sky, int *Ncatalog) {
+Catalog *load_catalogs_parallel (SkyList *sky, int *Ncatalog, char *syncfile) {
 
   char uniquer[12];
@@ -206 +212 @@
   }
 
+  // update syncfile here (save lots of I/O time):
+
+  // at this point, the remote relastro_client jobs are done loading their data.  in a
+  // parallel_images mode, the next image host can be launched while this image host now
+  // reads that
+
+  // NOTE: if I let all hosts load blindly, I saturate the data clients with too many
+  // relphot_client requests.  I need to have the master mediate this.  the master
+  // will not launch the next remote job until this one says it is done
+  if (syncfile) {
+    update_sync_file (syncfile, 1);
+  }
+
   // each host generates a BrightCatalog structure, with the measure, average, etc value
   // loaded into a single set of arrays (of MeasureTiny, AverageTiny, Secfilt).  I need to

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

@@ -51 +51 @@
 
                 xlist[N] = Mrel;
-                ylist[N] = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt]) - Mcal - Mmos - Mgrid - Mrel;
-                ilist[N] = PhotInstTiny (&catalog[i].measureT[m]);
+                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/relphot_client.c

@@ -39 +39 @@
     case MODE_LOAD: {
       int Ncatalog;
-      Catalog *catalog = load_catalogs (skylist, &Ncatalog, HOST_ID, HOSTDIR);
+      Catalog *catalog = load_catalogs (skylist, &Ncatalog, HOST_ID, HOSTDIR, NULL);
       if (!catalog) {
           fprintf (stderr, "ERROR loading catalogs from %s\n", CATDIR);

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

@@ -52 +52 @@
   if (NLOOP > 0) {
     /* load catalog data from region files (hostID is 0 since we are not a client */
-    catalog = load_catalogs (skylist, &Ncatalog, 0, NULL);
+    catalog = load_catalogs (skylist, &Ncatalog, 0, NULL, NULL);
     MARKTIME("-- load catalog data: %f sec\n", dtime);
   

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

@@ -38 +38 @@
   // unlink (IMAGE_TABLE); // XXX a bit risky, add some protection?
 
-  // XXX need to deal with mosaic vs image...
-  makeMosaics (image, Nimage);
+  makeMosaics (image, Nimage, TRUE);
 
   initImages (image, NULL, Nimage);
@@ -56 +55 @@
 
   /* load catalog data from region files (hostID is 0 since we are not a client */
-  catalog = load_catalogs (skylist, &Ncatalog, 0, NULL);
+  char *syncfile = make_filename (CATDIR, regionHosts->hosts[myHost].hostname, REGION_HOST_ID, "loadcat.sync");
+  catalog = load_catalogs (skylist, &Ncatalog, 0, NULL, syncfile);
   MARKTIME("-- load catalog data: %f sec\n", dtime);
+  free (syncfile);
+
   client_logger_message ("loaded catalog data\n");
 
-  // NOTE: if I let all hosts load blindly, I saturate the data clients with too many
-  // relphot_client requests.  I need to have the master mediate this.  the master
-  // will not launch the next remote job until this one says it is done
-  char *syncfile = make_filename (CATDIR, regionHosts->hosts[myHost].hostname, REGION_HOST_ID, "loadcat.sync");
-  update_sync_file (syncfile, 1);
-  
   // generate tables go from catID,objID -> catSeq,objSeq
   indexCatalogs (catalog, Ncatalog);

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

@@ -1 +1 @@
 # include "relphot.h"
+int save_images_backup (FITS_DB *db);
 
 int relphot_parallel_regions () {
@@ -63 +64 @@
   SkyList *skylist = SkyListByBounds (sky, -1, regionHosts->Rmin, regionHosts->Rmax, regionHosts->Dmin, regionHosts->Dmax);
 
+  // I have to save a copy because dvo_image_save and _unlock swap and free the data
+  save_images_backup (&db);
+
   /* update catalogs (in parallel) */
   reload_catalogs (skylist, flatcorr, 0, NULL);
@@ -74 +78 @@
 }
 
+int save_images_backup (FITS_DB *db) {
+
+  FITS_DB dbX;
+  char filename[1024];
+
+  INITTIME;
+
+  gfits_db_init (&dbX);
+  dbX.lockstate = LCK_XCLD;
+  dbX.timeout   = 60.0;
+  dbX.mode      = db->mode;
+  dbX.format    = db->format;
+
+  snprintf (filename, 1024, "%s.bck", ImageCat);
+  if (!gfits_db_lock (&dbX, filename)) {
+    fprintf (stderr, "can't lock backup image image catalog\n");
+    return (FALSE);
+  }
+    
+  // copy Images.dat data (all or only the subset / vtable elements?)  I think I need to dump
+  // the entire Image table, but with the updates in place I think this says: re-work the
+  // ftable/vtable usage in this program to be more sensible...
+  gfits_copy_header (&db->header,  &dbX.header);
+  gfits_copy_matrix (&db->matrix,  &dbX.matrix);
+  gfits_copy_header (&db->theader, &dbX.theader);
+  gfits_copy_ftable (&db->ftable,  &dbX.ftable);
+
+  dbX.ftable.header = &dbX.theader;
+
+  dvo_image_save (&dbX, VERBOSE);
+  dvo_image_unlock (&dbX); 
+  MARKTIME("-- save Image.dat.bck: %f sec\n", dtime);
+
+  return TRUE;
+}

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

@@ -35 +35 @@
         /* select measurements by mag limit */
         if (AreaSelect) {
-          r = catalog[i].averageT[j].R + catalog[i].measureT[m].dR / 3600.0;
-          d = catalog[i].averageT[j].D + catalog[i].measureT[m].dD / 3600.0;
+          r = catalog[i].measureT[m].R;
+          d = catalog[i].measureT[m].D;
           if ((coords = getCoords (m, i)) == NULL) goto markbad;
           RD_to_XY (&x, &y, r, d, coords);

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

@@ -1 +1 @@
 # include "relphot.h"
+# define UBERCAL_WEIGHT 100.0
 
 # define SKIP_THIS_MEAS(REASON) {                               \
@@ -13 +14 @@
     continue; }
 
+static float MagToFlux (float Mag) {
+  float Flux = pow(10.0, -0.4*(Mag));
+  return (Flux);
+}
+
 int print_measure_set_alt (Average *average, SecFilt *secfilt, Measure *measure) {
 
@@ -29 +35 @@
   return (TRUE);
 }
+
+int setMrelAverageForcedWarp (off_t measureOffset, int cat, int pass, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, off_t *found);
 
 int setMrel_catalog_alt (Catalog *catalog, int Nc, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt) {
@@ -52 +60 @@
     Average *average = catalog[Nc].average ? &catalog[Nc].average[j] : NULL;
     Measure *measure = catalog[Nc].measure ? &catalog[Nc].measure[m] : NULL;
-    setMrelAverageExposure (m, Nc, pass, flatcorr, results, average, &catalog[Nc].averageT[j], &catalog[Nc].secfilt[j*Nsecfilt], measure, &catalog[Nc].measureT[m], &catalog[Nc].found[Nsecfilt*j]);
+    setMrelAverageExposure (m, Nc, pass, flatcorr, results, average, &catalog[Nc].averageT[j], &catalog[Nc].secfilt[j*Nsecfilt], measure, &catalog[Nc].measureT[m], &catalog[Nc].found_t[Nsecfilt*j]);
 
     // only apply Stack operation on setMrelFinal in first pass 
@@ -59 +67 @@
       setGlobalObjStats (&catalog[Nc].average[j], &catalog[Nc].measure[m]);
     }
+
+    // only measure force-warp mean values if issetMrelFinal (make it optional?)
+    if (isSetMrelFinal) {
+      setMrelAverageForcedWarp (m, Nc, pass, results, &catalog[Nc].average[j], &catalog[Nc].secfilt[j*Nsecfilt], &catalog[Nc].measure[m], &catalog[Nc].foundWarp_t[Nsecfilt*j]);
+      setGlobalObjStats (&catalog[Nc].average[j], &catalog[Nc].measure[m]);
+    }
   }
   if (primaryCell) free (primaryCell);
@@ -64 +78 @@
 }
 
+// set mean of chip measurements (selected by photcode range for now):
 int setMrelAverageExposure (off_t measureOffset, int cat, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, AverageTiny *averageT, SecFilt *secfilt, Measure *measure, MeasureTiny *measureT, off_t *found) {
 
@@ -144 +159 @@
 
       // are we a PS1 exposure photcode?
-      if ((measureT[k].photcode >= 10000) && (measureT[k].photcode <= 10500)) {
-        NexpPS1 ++;
-      }
+      if (isGPC1chip(measureT[k].photcode)) NexpPS1 ++;
+
+      // SKIP gpc1 stack data
+      if (isGPC1stack(measureT[k].photcode)) continue;
+
+      // SKIP gpc1 forced-warp data
+      if (isGPC1warp(measureT[k].photcode)) continue;
 
       if (measureT[k].dbFlags & MEAS_BAD) SKIP_THIS_MEAS(Nbad); 
@@ -166 +185 @@
 
       // skip some absurd values NAN, < 0.0, > 30.0
-      Msys = PhotSysTiny (&measureT[k], &averageT[0], &secfilt[0]);
+      Msys = PhotSysTiny (&measureT[k], &averageT[0], &secfilt[0], MAG_CLASS_PSF);
       if (isnan(Msys)) SKIP_THIS_MEAS(Nsys);
       if (Msys <  0.0) SKIP_THIS_MEAS(Nsys);
@@ -197 +216 @@
           }
           havePS1 = TRUE;
-        }
-
-        // gpc1 stack data
-        if ((measure[k].photcode >= 11000) && (measure[k].photcode <= 11400)) {
-          continue;
         }
 
@@ -255 +269 @@
 
         // Map (Maplist) and Mkron (Mkronlist,dkronlist) are used to calculate mean mags per filter
-        float Map = PhotAperCat (&measure[k]);
+        float Map = PhotCat (&measure[k], MAG_CLASS_APER);
         if (!isnan(Map)) {
             Maplist[Nap] = Map - Mcal - Mmos - Mgrid;
-            daplist[Nap] = dMpsf;
-            waplist[Nap] = isUbercal ? 10.0 : 1.0;
+            daplist[Nap] = dMpsf;// XXX check on this...
+            waplist[Nap] = isUbercal ? UBERCAL_WEIGHT : 1.0;
             Nap ++;
         }
 
-        float Mkron = PhotKronCat (&measure[k]);
+        float Mkron = PhotCat (&measure[k], MAG_CLASS_KRON);
         if (!isnan(Mkron)) {
           Mkronlist[Nkron] = Mkron - Mcal - Mmos - Mgrid;
           dkronlist[Nkron] = measure[k].dMkron;
-          wkronlist[Nkron] = isUbercal ? 10.0 : 1.0;
+          wkronlist[Nkron] = isUbercal ? UBERCAL_WEIGHT : 1.0;
           Nkron ++;
         }
@@ -280 +294 @@
       Mpsflist[Npsf] = Msys - Mcal - Mmos - Mgrid;
       dpsflist[Npsf] = dMpsf;
-      wpsflist[Npsf] = isUbercal ? 10.0 : 1.0;
+      wpsflist[Npsf] = isUbercal ? UBERCAL_WEIGHT : 1.0;
 
       // NOTE: 
@@ -293 +307 @@
       if (refPhotcode) {
         if (code->code == refPhotcode->code) {
-          wpsflist[Npsf] = 100.0;
+          wpsflist[Npsf] = UBERCAL_WEIGHT;
         }
       }
@@ -304 +318 @@
     if (isSetMrelFinal && (pass == 0)) {
       if ((thisCode < 6) || (thisCode == 9)) {
-        secfilt[Nsec].Ncode = NexpPS1;
+        secfilt[Nsec].Ncode = NexpPS1; 
       } else {
         secfilt[Nsec].Ncode = Ncode; // 2MASS data if it exists
@@ -334 +348 @@
     liststats (Mpsflist, dpsflist, wpsflist, Npsf, psfstats);
 
-    secfilt[Nsec].M  = psfstats->mean;
-    secfilt[Nsec].dM = psfstats->error;
-    secfilt[Nsec].Xm = (psfstats->Nmeas > 1) ? 100.0*log10(psfstats->chisq + 1e-4) : NAN_S_SHORT;
+    secfilt[Nsec].M      = psfstats->mean;
+    secfilt[Nsec].dM     = psfstats->error;
+    secfilt[Nsec].Mchisq = (psfstats->Nmeas > 1) ? psfstats->chisq : NAN;
 
     // when running -averages, we have no information about the images, so we cannot set this
@@ -346 +360 @@
       found[Nsec] = TRUE;
 
-      secfilt[Nsec].Mstdev = 1000.0*psfstats->sigma; // Mstdev is in millimags (not enough space for more precision)
+      secfilt[Nsec].Mstdev = psfstats->sigma; // Mstdev is in millimags (not enough space for more precision)
       // secfilt[Nsec].Ncode = Ncode;
       secfilt[Nsec].Nused = psfstats->Nmeas;
 
-      secfilt[Nsec].M_80 = 1000 * psfstats->Upper80;
-      secfilt[Nsec].M_20 = 1000 * psfstats->Lower20;
+      secfilt[Nsec].Mmax = psfstats->max;
+      secfilt[Nsec].Mmin = psfstats->min;
 
       // NOTE : use the modified weight for apmags as well as psf mags
       liststats (Maplist, daplist, waplist, Nap, apstats);
       secfilt[Nsec].Map  = Nap > 0 ? apstats->mean : NAN; 
+      secfilt[Nsec].dMap  = Nap > 0 ? apstats->error : NAN; 
+      secfilt[Nsec].sMap  = Nap > 0 ? apstats->sigma : NAN; 
+      secfilt[Nsec].NusedAp  = Nap;
 
       liststats (Mkronlist, dkronlist, wkronlist, Nkron, kronstats);
       secfilt[Nsec].Mkron  = Nkron > 0 ? kronstats->mean  : NAN; 
       secfilt[Nsec].dMkron = Nkron > 0 ? kronstats->error : NAN; 
+      secfilt[Nsec].sMkron = Nkron > 0 ? kronstats->sigma : NAN; 
+      secfilt[Nsec].NusedKron  = Nkron;
 
       // NOTE: for 2MASS measurements, Next should be 1, as should N
@@ -472 +491 @@
 
       // only examine gpc1 stack data
-      // XXX this is absurdly hardwired (along with several photcode tests)
-      if (measure[k].photcode < 11000) continue;
-      if (measure[k].photcode > 11400) continue;
+      if (!isGPC1stack(measure[k].photcode)) continue;
 
       haveStack = TRUE;
@@ -493 +510 @@
       if (isnan(Mgrid)) SKIP_THIS_MEAS_STACK(Ngrid);
 
-      Msys = PhotSys (&measure[k], &average[0], &secfilt[0]);
+      Msys = PhotSys (&measure[k], &average[0], &secfilt[0], MAG_CLASS_PSF);
       if (isnan(Msys)) SKIP_THIS_MEAS_STACK(Nsys);
 
@@ -581 +598 @@
 
     // need to put in AB mag factor to get to Janskies (or uJy?)
-    secfilt[Nsec].FluxPSF   = zpFactor * measure[k].FluxPSF;  
-    secfilt[Nsec].dFluxPSF  = zpFactor * measure[k].dFluxPSF; 
-    secfilt[Nsec].FluxKron  = zpFactor * measure[k].FluxKron; 
-    secfilt[Nsec].dFluxKron = zpFactor * measure[k].dFluxKron;
+    secfilt[Nsec].FpsfStk   = zpFactor * measure[k].FluxPSF;  
+    secfilt[Nsec].dFpsfStk  = zpFactor * measure[k].dFluxPSF; 
+    secfilt[Nsec].FkronStk  = zpFactor * measure[k].FluxKron; 
+    secfilt[Nsec].dFkronStk = zpFactor * measure[k].dFluxKron;
+    secfilt[Nsec].FapStk    = zpFactor * measure[k].FluxAp;
+    secfilt[Nsec].dFapStk   = zpFactor * measure[k].dFluxAp;
+
+    secfilt[Nsec].MpsfStk   = (measure[k].FluxPSF  > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FpsfStk) : NAN;
+    secfilt[Nsec].MkronStk  = (measure[k].FluxKron > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FkronStk) : NAN;
+    secfilt[Nsec].MapStk    = (measure[k].FluxAp   > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FapStk) : NAN;
 
     secfilt[Nsec].stackDetectID = ID;
@@ -607 +630 @@
 }
 
+# undef SKIP_THIS_MEAS
+# define SKIP_THIS_MEAS(REASON) {                               \
+    measure[k].dbFlags |= ID_MEAS_SKIP_PHOTOM;  \
+    results->REASON ++;                                         \
+    continue; }
+
+// set mean of forced-warp measurements (selected by photcode range for now):
+// somewhat simplified relative to chip-photometry:
+// * no grid, no mosaic, no 2MASS, no SYNTH, no Ubercal, no flatcorr
+int setMrelAverageForcedWarp (off_t measureOffset, int cat, int pass, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, off_t *found) {
+
+  off_t k;
+  float Fsys = 0, Mcal= 0;
+
+  // we are measuring means for 3 types of FLUXes: psf, ap, kron.  I am using the psf mag
+  // error for the ap mags, but krons have their own errors.  it is an open question if I
+  // should be doing weighted or unweighted fits (this is a user option)
+  double *Fpsflist  = results->Mpsflist;
+  double *dpsflist  = results->dpsflist;
+  double *wpsflist  = results->wpsflist;
+
+  double *Faplist   = results->Maplist;
+  double *daplist   = results->daplist;
+  double *waplist   = results->waplist;
+
+  double *Fkronlist = results->Mkronlist;
+  double *dkronlist = results->dkronlist;
+  double *wkronlist = results->wkronlist;
+
+  StatType *psfstats  = &results->psfstats;
+  StatType *apstats   = &results->apstats;
+  StatType *kronstats = &results->kronstats;
+
+  // option for a test print
+  if (FALSE && (average[0].objID == 0x7146) && (average[0].catID == 0x49d8)) {
+    fprintf (stderr, "test obj\n");
+    print_measure_set_alt (average, secfilt, measure);
+  }
+
+  int Ns;
+  for (Ns = 0; Ns < Nphotcodes; Ns++) {
+
+    int thisCode = photcodes[Ns][0].code;
+    int Nsec = GetPhotcodeNsec(thisCode);
+
+    /* calculate the average mag in this SEC photcode for a single star */
+
+    /* star/photcodes already calibrated */
+    if (found[Nsec]) continue;  
+      
+    off_t meas = measureOffset;
+
+    int Nap = 0;
+    int Npsf = 0;
+    int Nkron = 0;
+    for (k = 0; k < average[0].Nmeasure; k++, meas++) {
+
+      // only examine gpc1 forced-warp data
+      if (!isGPC1warp(measure[k].photcode)) continue;
+
+      // skip measurements that do not match the current photcode
+      PhotCode *code = GetPhotcodebyCode (measure[k].photcode);
+      if (!code) continue;
+      if (code->equiv != thisCode) { continue; }
+
+      if (measure[k].dbFlags & MEAS_BAD) SKIP_THIS_MEAS(Nbad); 
+
+      if (getImageEntry (meas, cat) < 0) {
+        // measurements without an image are either external reference photometry or
+        // 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
+        Mcal = measure[k].Mcal; // check that this is zero for loaded REF value
+      } else {
+        // use getMcal not getMcal_alt?
+        Mcal  = getMcal_alt (meas, cat, NULL, measure[k].Xccd, measure[k].Yccd);
+        // Mcal  = getMcal (meas, cat);
+        if (isnan(Mcal))  SKIP_THIS_MEAS(Ncal);
+      }
+      float Fcal = MagToFlux(-Mcal);
+
+      // in the calculations below, 
+      // ...list gives the error per measurement, wlist gives the weight
+      // we can modify the error and weight in a few ways:
+      // 1) MIN_ERROR guarantees a floor
+      // 2) photomErrSys is added in quadrature as a sytematic error, set per photcode
+
+      // skip some absurd values NAN
+      // NOTE : I am using PhotCat not PhotSys for now since GPC1 chip-to-chip color terms
+      // are small (and not measured)
+      float Fpsf = PhotFluxCat (&measure[k], MAG_CLASS_PSF);
+      if (isnan(Fsys)) SKIP_THIS_MEAS(Nsys);
+      // if (Msys <  0.0) SKIP_THIS_MEAS(Nsys);
+      // if (Msys > 30.0) SKIP_THIS_MEAS(Nsys);
+
+      float dFpsf = PhotFluxCatErr (&measure[k], MAG_CLASS_PSF);
+      dFpsf = MAX (dFpsf, MIN_ERROR*Fpsf); // MIN_ERROR is a fractional error
+      Fpsflist[Npsf] = Fpsf * Fcal;
+      dpsflist[Npsf] = dFpsf * Fcal;
+      wpsflist[Npsf] = 1.0;
+      Npsf ++;
+
+      float Fap = PhotFluxCat (&measure[k], MAG_CLASS_APER);
+      float dFap = PhotFluxCatErr (&measure[k], MAG_CLASS_APER);
+      if (!isnan(Fap)) {
+        Faplist[Nap] = Fap * Fcal;
+        daplist[Nap] = dFap * Fcal;
+        waplist[Nap] = 1.0; // drop weight lists?
+        Nap ++;
+      }
+
+      float Fkron = PhotFluxCat (&measure[k], MAG_CLASS_KRON);
+      float dFkron = PhotFluxCatErr (&measure[k], MAG_CLASS_KRON);
+      if (!isnan(Fkron)) {
+        Fkronlist[Nkron] = Fkron * Fcal;
+        dkronlist[Nkron] = dFkron * Fcal;
+        wkronlist[Nkron] = 1.0;
+        Nkron ++;
+      }
+    }
+    if (Npsf < 1) continue;
+
+    found[Nsec] = TRUE;
+
+    liststats (Fpsflist, dpsflist, wpsflist, Npsf, psfstats);
+
+    secfilt[Nsec].FpsfWrp  = psfstats->mean;
+    secfilt[Nsec].dFpsfWrp = psfstats->error;
+    secfilt[Nsec].sFpsfWrp = psfstats->sigma; // Mstdev is in millimags (not enough space for more precision)
+    secfilt[Nsec].NusedWrp = psfstats->Nmeas;
+    secfilt[Nsec].MpsfWrp  = isnan(secfilt[Nsec].FpsfWrp) ? NAN : 8.9 - 2.5*log10(secfilt[Nsec].FpsfWrp); // 8.9 since flux is in Jy
+
+    // NOTE : use the modified weight for apmags as well as psf mags
+    liststats (Faplist, daplist, waplist, Nap, apstats);
+    secfilt[Nsec].FapWrp     = Nap > 0 ? apstats->mean : NAN; 
+    secfilt[Nsec].dFapWrp    = Nap > 0 ? apstats->error : NAN; 
+    secfilt[Nsec].sFapWrp    = Nap > 0 ? apstats->sigma : NAN; 
+    secfilt[Nsec].NusedApWrp = Nap;
+    secfilt[Nsec].MapWrp  = isnan(secfilt[Nsec].FapWrp) ? NAN : 8.9 - 2.5*log10(secfilt[Nsec].FapWrp); // 8.9 since flux is in Jy
+
+    liststats (Fkronlist, dkronlist, wkronlist, Nkron, kronstats);
+    secfilt[Nsec].FkronWrp     = Nkron > 0 ? kronstats->mean  : NAN; 
+    secfilt[Nsec].dFkronWrp    = Nkron > 0 ? kronstats->error : NAN; 
+    secfilt[Nsec].sFkronWrp    = Nkron > 0 ? kronstats->sigma : NAN; 
+    secfilt[Nsec].NusedKronWrp = Nkron;
+    secfilt[Nsec].MkronWrp  = isnan(secfilt[Nsec].FkronWrp) ? NAN : 8.9 - 2.5*log10(secfilt[Nsec].FkronWrp); // 8.9 since flux is in Jy
+  }
+  return (TRUE);
+}
+
 int setGlobalObjStats (Average *average, Measure *measure) {
 

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

@@ -85 +85 @@
 
   /* set catalog[0].found[i] = FALSE */
-  ALLOCATE (catalog[0].found, off_t, MAX (1, Nsecfilt*catalog[0].Naverage));
+  ALLOCATE (catalog[0].found_t, off_t, MAX (1, Nsecfilt*catalog[0].Naverage));
+  ALLOCATE (catalog[0].foundWarp_t, off_t, MAX (1, Nsecfilt*catalog[0].Naverage));
   for (i = 0; i < Nsecfilt*catalog[0].Naverage; i++) {
-    catalog[0].found[i] = FALSE;
+    catalog[0].found_t[i] = FALSE;
+    catalog[0].foundWarp_t[i] = FALSE;
   }
 
@@ -194 +196 @@
       
       /* star/photcodes already calibrated */
-      if (catalog[0].found[Nsecfilt*i+Nsec]) continue;  
+      if (catalog[0].found_t[Nsecfilt*i+Nsec]) continue;  
 
       m = catalog[0].average[i].measureOffset;
@@ -229 +231 @@
         if ((pass < 4) && ImagSelect) {
           if (Nim > -1) {
-            mag = PhotInst (&catalog[0].measure[m]);
+            mag = PhotInst (&catalog[0].measure[m], MAG_CLASS_PSF);
             if (mag < ImagMin) goto skip;
             if (mag > ImagMax) goto skip;

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

@@ -123 +123 @@
   meanmags->M  = secfilt->M;
   meanmags->dM = secfilt->dM;
-  meanmags->Xm = secfilt->Xm;
+  meanmags->Mchisq = secfilt->Mchisq;
   meanmags->Nsec = Nsec; // key to secfilt entry