Changeset 31450

Timestamp:

May 5, 2011, 10:43:45 AM (15 years ago)

Author:

eugene

Message:

various improvements to relphot, related support

Location:

trunk/Ohana/src

Files:

• dvomerge/src/dvo_image_merge_dbs.c (modified) (1 diff)
• dvomerge/src/dvoverify.c (modified) (4 diffs)
• libdvo/include/dvo.h (modified) (4 diffs)
• libdvo/src/dvo_catalog_split.c (modified) (1 diff)
• libdvo/src/dvo_photcode_ops.c (modified) (1 diff)
• opihi/cmd.data/Makefile (modified) (1 diff)
• opihi/cmd.data/gridify.c (modified) (3 diffs)
• opihi/cmd.data/init.c (modified) (2 diffs)
• opihi/cmd.data/vshift.c (copied) (copied from branches/eam_branches/ipp-20110404/Ohana/src/opihi/cmd.data/vshift.c )
• opihi/dvo/Makefile (modified) (1 diff)
• opihi/dvo/avextract.c (modified) (3 diffs)
• opihi/dvo/init.c (modified) (2 diffs)
• opihi/dvo/objectcoverage.c (copied) (copied from branches/eam_branches/ipp-20110404/Ohana/src/opihi/dvo/objectcoverage.c )
• photdbc/src/make_subcatalog.c (modified) (2 diffs)
• relastro/src/UpdateObjects.c (modified) (1 diff)
• relphot/doc/allsky.txt (copied) (copied from branches/eam_branches/ipp-20110404/Ohana/src/relphot/doc/allsky.txt )
• relphot/include/relphot.h (modified) (4 diffs)
• relphot/src/GridOps.c (modified) (18 diffs)
• relphot/src/ImageOps.c (modified) (16 diffs)
• relphot/src/MosaicOps.c (modified) (19 diffs)
• relphot/src/StarOps.c (modified) (31 diffs)
• relphot/src/args.c (modified) (1 diff)
• relphot/src/bcatalog.c (modified) (10 diffs)
• relphot/src/global_stats.c (modified) (1 diff)
• relphot/src/initialize.c (modified) (3 diffs)
• relphot/src/load_images.c (modified) (1 diff)
• relphot/src/plot_scatter.c (modified) (1 diff)
• relphot/src/reload_catalogs.c (modified) (4 diffs)
• relphot/src/relphot.c (modified) (10 diffs)
• relphot/src/relphot_objects.c (modified) (3 diffs)
• relphot/src/select_images.c (modified) (7 diffs)
• relphot/src/setExclusions.c (modified) (3 diffs)
• relphot/src/setMrelFinal.c (modified) (6 diffs)

trunk/Ohana/src/dvomerge/src/dvo_image_merge_dbs.c

@@ -202 +202 @@
     if (newID == 0) {
       fprintf (stderr, "cannot find image ID "OFF_T_FMT"\n",  oldID);
-      exit (2);
+      // exit (2);
     }
     catalog[0].measure[i].imageID = newID;

trunk/Ohana/src/dvomerge/src/dvoverify.c

@@ -148 +148 @@
 int VerifyTableFile (char *filename) {
 
-  int status;
+  int status, Next;
   off_t Nbytes;
   Header header;
@@ -205 +205 @@
   // scan all extentions
   Nbytes = 0;
+  Next = -1;
   if (DEBUG) fprintf (stderr, "sizes: ("OFF_T_FMT" vs "OFF_T_FMT")\n", Nbytes, fileStats.st_size);
   while (Nbytes < fileStats.st_size) {
@@ -210 +211 @@
     // Check on the PHU
     if (!gfits_fread_header (file, &header)) {
-      fprintf (stderr, "unable to read PHU header for %s\n", filename);
+      if (Next == -1) {
+        fprintf (stderr, "unable to read PHU header for %s\n", filename);
+      } else {
+        fprintf (stderr, "unable to read header for %s, extension %d (or file has excess bytes)\n", filename, Next);
+      }
       fclose(file);
       return (FALSE);
@@ -243 +248 @@
       }
     }
+    Next ++;
   }
   if (DEBUG) fprintf (stderr, "file is good: %s\n", filename);

trunk/Ohana/src/libdvo/include/dvo.h

@@ -20 +20 @@
   DVO_FORMAT_PS1_DEV_2,
   DVO_FORMAT_PS1_DEV_3,
+  DVO_FORMAT_PS1_REF,
   DVO_FORMAT_PS1_V1,
   DVO_FORMAT_PS1_V2,
-  DVO_FORMAT_PS1_REF
 } DVOTableFormat;
 
@@ -210 +210 @@
 } PhotCodeData;
 
+// a reduced-subset structure for relphot
+typedef struct {
+  double         R;
+  double         D;
+  unsigned short Nmeasure;
+  int            measureOffset;
+  uint32_t       flags;
+} AverageTiny;
+
+// a reduced-subset structure for relphot
+typedef struct {
+  float          dR;
+  float          dD;
+  float          M;
+  float          Mcal;
+  float          dM;
+  float          airmass;
+  float          Xccd;
+  float          Yccd;
+  float          dt;
+  int            t;
+  unsigned int   averef;
+  unsigned int   imageID;
+  unsigned int   dbFlags;
+  unsigned short photcode;
+} MeasureTiny;
+
 /* a catalog contains this data */
 typedef struct Catalog {
@@ -226 +253 @@
   off_t Naves_disk, Nmeas_disk, Nmiss_disk, Nsecf_disk; /* current number of each component on disk */
   off_t Naves_off,  Nmeas_off,  Nmiss_off,  Nsecf_off;  /* index of first loaded data value */
+
+  // note that we use these for the full-sky relphot analysis
+  AverageTiny *averageT;
+  MeasureTiny *measureT; 
 
   /* the Nsecf_* values above are number of table rows (eg, Naverage*Nsecfilt) */
@@ -330 +361 @@
 int PhotColor (Average *average, SecFilt *secfilt, Measure *measure, int c1, int c2, double *color);
 
+float PhotInstTiny (MeasureTiny *measure);
+float PhotCatTiny (MeasureTiny *measure);
+float PhotAperTiny (MeasureTiny *measure);
+float PhotSysTiny (MeasureTiny *measure, AverageTiny *average, SecFilt *secfilt);
+float PhotRelTiny (MeasureTiny *measure, AverageTiny *average, SecFilt *secfilt);
+float PhotCalTiny (MeasureTiny *thisone, AverageTiny *average, SecFilt *secfilt, MeasureTiny *measure, PhotCode *code);
+float PhotAveTiny (PhotCode *code, AverageTiny *average, SecFilt *secfilt);
+float PhotRefTiny (PhotCode *code, AverageTiny *average, SecFilt *secfilt, MeasureTiny *measure);
+float PhotXmTiny (PhotCode *code, AverageTiny *average, SecFilt *secfilt);
+float PhotdMTiny (PhotCode *code, AverageTiny *average, SecFilt *secfilt);
+
+float PhotColorForCodeTiny (AverageTiny *average, SecFilt *secfilt, MeasureTiny *measure, PhotCode *code);
+int PhotColorTiny (AverageTiny *average, SecFilt *secfilt, MeasureTiny *measure, int c1, int c2, double *color);
+
+
 PhotCodeData *GetPhotcodeTable (void);
 void SetPhotcodeTable (PhotCodeData *);

trunk/Ohana/src/libdvo/src/dvo_catalog_split.c

@@ -448 +448 @@
     catalog[0].measure = FtableToMeasure (&ftable, &Nmeasure, &catalog[0].catformat);
     if (Nmeasure != Nrows) {
-      fprintf (stderr, "Warning: mismatch between Nmeasure in PHU and Table headers ("OFF_T_FMT" vs "OFF_T_FMT")\n",  Nmeasure,  Nrows);
+      // XXX this condition denotes the eof has been reached; not an error or a warning
+      // fprintf (stderr, "Warning: mismatch between Nmeasure in PHU and Table headers ("OFF_T_FMT" vs "OFF_T_FMT")\n",  Nmeasure,  Nrows);
     }
     gfits_free_header (&header);

trunk/Ohana/src/libdvo/src/dvo_photcode_ops.c

@@ -536 +536 @@
   return (TRUE);
 }
+
+/******** alternate photometry conversion functions using MeasureTiny and AverageTiny *********/
+float PhotInstTiny (MeasureTiny *measure) {
+
+  int Np;
+  float M;
+
+  Np = photcodes[0].hashcode[measure[0].photcode];
+  if (Np == -1) return (NAN);
+
+  if (photcodes[0].code[Np].type == PHOT_REF) {
+    M = measure[0].M;
+    return (M);
+  }
+
+  M = measure[0].M - measure[0].dt - ZERO_POINT;
+          
+  return (M);
+
+}
+
+float PhotCatTiny (MeasureTiny *measure) {
+
+  int Np;
+  float Mcat;
+  PhotCode *code;
+
+  Np = photcodes[0].hashcode[measure[0].photcode];
+  if (Np == -1) return (NAN);
+
+  if (photcodes[0].code[Np].type == PHOT_REF) {
+    Mcat = measure[0].M;
+    return (Mcat);
+  }
+  code = &photcodes[0].code[Np];
+  Mcat = measure[0].M - ZERO_POINT + code[0].K*(measure[0].airmass - 1.000) + SCALE*code[0].C;
+  
+  return (Mcat);
+}
+
+# if (0)
+float PhotAperTiny (MeasureTiny *measure) {
+
+  int Np;
+  float Mcat;
+  PhotCode *code;
+
+  Np = photcodes[0].hashcode[measure[0].photcode];
+  if (Np == -1) return (NAN);
+
+  if (photcodes[0].code[Np].type == PHOT_REF) {
+    Mcat = measure[0].Map;
+    return (Mcat);
+  }
+  code = &photcodes[0].code[Np];
+  Mcat = measure[0].Map - ZERO_POINT + code[0].K*(measure[0].airmass - 1.000) + SCALE*code[0].C;
+  
+  return (Mcat);
+}
+# endif
+
+float PhotSysTiny (MeasureTiny *measure, AverageTiny *average, SecFilt *secfilt) {
+
+  int i, Np;
+  float Mcat, Mcol, Msys, mc, Mc;
+  PhotCode *code;
+
+  Np = photcodes[0].hashcode[measure[0].photcode];
+  if (Np == -1) return (NAN);
+
+  if (photcodes[0].code[Np].type == PHOT_REF) {
+    Msys = measure[0].M;
+    return (Msys);
+  }
+  code = &photcodes[0].code[Np];
+  Mcat = measure[0].M - ZERO_POINT + code[0].K*(measure[0].airmass - 1.000) + SCALE*code[0].C;
+
+  /* for DEP, color must be made of PRI/SEC */
+  mc = PhotColorForCodeTiny (average, secfilt, NULL, code);
+  if (isnan(mc)) return (Mcat);
+  mc = mc - SCALE*code[0].dX;
+
+  Mc = mc;
+  Mcol = 0;
+  for (i = 0; i < code[0].Nc; i++) {
+    Mcol += code[0].X[i]*Mc;
+    Mc *= mc;
+  }
+  Msys = Mcat + Mcol;
+  return (Msys);
+}
+
+float PhotRelTiny (MeasureTiny *measure, AverageTiny *average, SecFilt *secfilt) {
+
+  int i, Np;
+  float Mcat, Mcol, Mrel, mc, Mc;
+  PhotCode *code;
+
+  Np = photcodes[0].hashcode[measure[0].photcode];
+  if (Np == -1) return (NAN);
+
+  if (photcodes[0].code[Np].type == PHOT_REF) {
+    Mcat = measure[0].M;
+    return (Mcat);
+  }
+  code = &photcodes[0].code[Np];
+  Mrel = measure[0].M - ZERO_POINT + code[0].K*(measure[0].airmass - 1.000) + SCALE*code[0].C - measure[0].Mcal;
+
+  /* for DEP, color must be made of PRI/SEC */
+  mc = PhotColorForCodeTiny (average, secfilt, NULL, code);
+  if (isnan(mc)) return (Mrel);
+  mc = mc - SCALE*code[0].dX;
+
+  Mc = mc;
+  Mcol = 0;
+  for (i = 0; i < code[0].Nc; i++) {
+    Mcol += code[0].X[i]*Mc;
+    Mc *= mc;    /* the 0.001 is needed for higher order terms to keep the units mag = mag^n */
+  }
+  Mrel += Mcol;
+  return (Mrel);
+}
+
+/* return calibrated magnitude from measure for given photcode */
+float PhotCalTiny (MeasureTiny *thisone, AverageTiny *average, SecFilt *secfilt, MeasureTiny *measure, PhotCode *code) {
+
+  int i, Np; 
+  float Mcal, Mrel, Mcol, mc, Mc;
+
+  if (code == NULL) return NAN;
+
+  /* code must be the matching PRI/SEC code for this measurement or an equivalent ALT */
+  Np = photcodes[0].hashcode[thisone[0].photcode];
+  if (Np == -1) return (NAN);
+
+  if (photcodes[0].code[Np].type == PHOT_REF) {
+    Mrel = thisone[0].M;
+    return (Mrel);
+  }
+  if (code[0].code != photcodes[0].code[Np].equiv) return (NAN);
+
+  Mcal = PhotRelTiny (thisone, average, secfilt) + SCALE*code[0].C;
+
+  mc = PhotColorForCodeTiny (average, secfilt, measure, code);
+  if (isnan(mc)) return (Mcal);
+  mc = mc - SCALE*code[0].dX;
+
+  Mc = mc;
+  Mcol = 0;
+  for (i = 0; i < code[0].Nc; i++) {
+    Mcol += code[0].X[i]*Mc;
+    Mc *= mc;
+  }
+  Mcal += Mcol;
+  return (Mcal);
+}
+
+/* color term may not use DEP magnitude */
+float PhotColorForCodeTiny (AverageTiny *average, SecFilt *secfilt, MeasureTiny *measure, PhotCode *code) {
+
+  int i, Ns1, Ns2, Ns;
+  float m1, m2, mc;
+  PhotCode *color;
+
+  m1 = m2 = NAN;
+
+  if (measure == NULL) {
+    Ns1 = photcodes[0].hashNsec[code[0].c1];
+    Ns2 = photcodes[0].hashNsec[code[0].c2];
+  
+    m1 = (Ns1 == -1) ? NAN : secfilt[Ns1].M;
+    m2 = (Ns2 == -1) ? NAN : secfilt[Ns2].M;
+    mc = (isnan(m1) || isnan(m2)) ? NAN : (m1 - m2);
+    return (mc);
+  }
+
+  /* find magnitude matching first color term */
+  color = GetPhotcodebyCode (code[0].c1);
+  if (color == NULL) return (NAN);
+  if (color[0].type == PHOT_REF) {
+    for (i = 0; (i < average[0].Nmeasure) && (isnan(m1)); i++) {
+      if (measure[i].photcode == color[0].code) {
+        m1 = measure[i].M;
+      }
+    }
+  } else {
+    Ns = photcodes[0].hashNsec[color[0].code];
+    m1 = (Ns == -1) ? NAN : secfilt[Ns].M;
+  }     
+
+  /* find magnitude matching second color term */
+  color = GetPhotcodebyCode (code[0].c2);
+  if (color == NULL) return (NAN);
+  if (color[0].type == PHOT_REF) {
+    for (i = 0; (i < average[0].Nmeasure) && (isnan(m2)); i++) {
+      if (measure[i].photcode == color[0].code) {
+        m2 = measure[i].M;
+      }
+    }
+  } else {
+    Ns = photcodes[0].hashNsec[color[0].code];
+    m2 = (Ns == -1) ? NAN : secfilt[Ns].M;
+  }     
+  mc = (isnan(m1) || isnan(m2)) ? NAN : (m1 - m2);
+  return (mc);
+}
+
+/* return calibrated magnitude from average/secfilt for given photcode */
+float PhotRefTiny (PhotCode *code, AverageTiny *average, SecFilt *secfilt, MeasureTiny *measure) {
+
+  int i, Ns;
+  float Mave, Mref, Mcol, mc;
+  double Mc;
+
+  if (code == NULL) return NAN;
+
+  Ns = photcodes[0].hashNsec[code[0].code];
+  Mave = (Ns == -1) ? NAN : secfilt[Ns].M;
+  Mref = Mave + SCALE*code[0].C;
+
+  mc = PhotColorForCodeTiny (average, secfilt, measure, code);
+  if (isnan(mc)) return (Mref);
+  mc = mc - SCALE*code[0].dX;
+
+  Mc = mc;
+  Mcol = 0;
+  for (i = 0; i < code[0].Nc; i++) {
+    Mcol += code[0].X[i]*Mc;
+    Mc *= mc;    /* the 0.001 is needed for higher order terms to keep the units mag = mag^n */
+  }
+  Mref += Mcol;
+  return (Mref);
+}
+
+/***/
+float PhotAveTiny (PhotCode *code, AverageTiny *average, SecFilt *secfilt) {
+
+  int Ns;
+  float Mave;
+
+  if (code == NULL) return NAN;
+
+  Ns = photcodes[0].hashNsec[code[0].code];
+  Mave = (Ns == -1) ? NAN : secfilt[Ns].M;
+  return (Mave);
+}
+
+float PhotdMTiny (PhotCode *code, AverageTiny *average, SecFilt *secfilt) {
+
+  int Ns;
+  float dM;
+
+  if (code == NULL) return NAN;
+
+  Ns = photcodes[0].hashNsec[code[0].code];
+  dM  = (Ns == -1) ? NAN : secfilt[Ns].dM;
+  return (dM);
+}
+
+// XXX return NAN or NAN_S_SHORT? (secfilt->Xm is short)
+float PhotXmTiny (PhotCode *code, AverageTiny *average, SecFilt *secfilt) {
+
+  int Ns;
+  short Mi;
+  float Xm;
+
+  if (code == NULL) return NAN;
+
+  Ns = photcodes[0].hashNsec[code[0].code];
+  Mi = (Ns == -1) ? NAN : secfilt[Ns].Xm;
+  Xm = (isnan(Mi)) ? -1.0 : pow (10.0, 0.01*Mi);
+  return (Xm);
+}
+
+/* given a photcode pair c1 & c2, return the color of this star (NaN if not found) */
+int PhotColorTiny (AverageTiny *average, SecFilt *secfilt, MeasureTiny *measure, int c1, int c2, double *color) {
+
+  int i, Ns;
+  double M1, M2, dM;
+  PhotCode *code;
+  
+  code = GetPhotcodebyCode (c1);
+  if (code == NULL) return (FALSE);
+  if (code[0].type == PHOT_REF) {
+    for (i = 0; i < average[0].Nmeasure; i++) {
+      if (measure[i].photcode == c1) {
+        M1 = measure[i].M;
+        goto filter1;
+      }
+    }   
+    return (FALSE);
+  } else {
+    Ns = photcodes[0].hashNsec[code[0].code];
+    M1 = (Ns == -1) ? NAN : secfilt[Ns].M;
+  }     
+
+filter1:
+  code = GetPhotcodebyCode (c2);
+  if (code == NULL) return (FALSE);
+  if (code[0].type == PHOT_REF) {
+    for (i = 0; i < average[0].Nmeasure; i++) {
+      if (measure[i].photcode == c2) {
+        M2 = measure[i].M;
+        goto filter2;
+      }
+    }   
+    return (FALSE);
+  } else {
+    Ns = photcodes[0].hashNsec[code[0].code];
+    M2 = (Ns == -1) ? NAN : secfilt[Ns].M;
+  }     
+  
+filter2:
+
+  dM = M1 - M2;
+  *color = dM;
+  
+  return (TRUE);
+}
+/***********************************************/
 
 // Create a map between the secfilt values from one photcode table to another

trunk/Ohana/src/opihi/cmd.data/Makefile

@@ -141 +141 @@
 $(SRC)/vpop.$(ARCH).o              \
 $(SRC)/vroll.$(ARCH).o             \
+$(SRC)/vshift.$(ARCH).o            \
 $(SRC)/vsmooth.$(ARCH).o           \
 $(SRC)/vstats.$(ARCH).o            \

trunk/Ohana/src/opihi/cmd.data/gridify.c

@@ -4 +4 @@
 
   int i, Nx, Ny, Xb, Yb, Normalize, N;
-  float Xmin, Xmax, dX, Ymin, Ymax, dY;
+  float Xmin, Xmax, dX, Ymin, Ymax, dY, initValue;
   float *buf, *val;
   int *Nval;
@@ -15 +15 @@
     remove_argument (N, &argc, argv);
     Normalize = FALSE;
+  }
+
+  initValue = 0.0;
+  if ((N = get_argument (argc, argv, "-init-value"))) {
+    remove_argument (N, &argc, argv);
+    initValue = atof(argv[N]);
+    remove_argument (N, &argc, argv);
   }
 
@@ -69 +76 @@
   buf = (float *) bf[0].matrix.buffer;
   for (i = 0; i < Nx*Ny; i++) {
+    buf[i] = initValue;
     if (Normalize) {
       if (Nval[i] == 0) {
-        buf[i] = 0;
         continue;
       }

trunk/Ohana/src/opihi/cmd.data/init.c

@@ -130 +130 @@
 int vstats           PROTO((int, char **));
 int vroll            PROTO((int, char **));
+int vshift           PROTO((int, char **));
 int vpop             PROTO((int, char **));
 int vsmooth          PROTO((int, char **));
@@ -274 +275 @@
   {1, "vmaxwell",     vmaxwell,         "fit a Maxwellian to a vector"},
   {1, "vpop",         vpop,             "remove first element of a vector"},
-  {1, "vroll",        vroll,            "roll vector elements"},
+  {1, "vroll",        vroll,            "roll vector elements by 1 entry"},
+  {1, "vshift",       vshift,           "shift vector elements by arbitrary amount"},
   {1, "vsmooth",      vsmooth,          "Gaussian smooth of a vector"},
   {1, "vstats",       vstats,           "statistics on a vector"},

trunk/Ohana/src/opihi/dvo/Makefile

@@ -80 +80 @@
 $(SRC)/mextract.$(ARCH).o               \
 $(SRC)/mmextract.$(ARCH).o              \
+$(SRC)/objectcoverage.$(ARCH).o         \
 $(SRC)/photcodes.$(ARCH).o              \
 $(SRC)/pmeasure.$(ARCH).o               \

trunk/Ohana/src/opihi/dvo/avextract.c

@@ -5 +5 @@
   off_t i, j, n, m;
   int N, Npts, NPTS, last, next, state, Nfields, Nreturn, Ncstack, Nstack;
-  int Nsecfilt, mode, VERBOSE;
+  int Nsecfilt, mode, VERBOSE, needMeasures;
   char **cstack, name[1024];
   void *Signal;
@@ -101 +101 @@
   }
 
+  // check the requested fields : are all average/secfilt entries, or do we need measures?
+  needMeasures = FALSE;
+  for (i = 0; !needMeasures && (i < Nfields); i++) {
+    if (fields[i].magMode == MAG_NONE) continue;
+    if (fields[i].photcode == NULL) continue; // assert this?
+    if (fields[i].photcode[0].type == PHOT_REF) needMeasures = TRUE;
+    if (fields[i].photcode[0].type == PHOT_DEP) needMeasures = TRUE;
+  }
+
   // grab data from all selected sky regions
   Signal = signal (SIGINT, handle_interrupt);
@@ -107 +116 @@
     /* lock, load, unlock catalog */
     catalog.filename = skylist[0].filename[i];
-    catalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_SECF;
+    catalog.catflags = LOAD_AVES | LOAD_SECF;
+    if (needMeasures) {
+      catalog.catflags |= LOAD_MEAS;
+    }
     catalog.Nsecfilt = 0;
 

trunk/Ohana/src/opihi/dvo/init.c

@@ -40 +40 @@
 int lightcurve      PROTO((int, char **));
 int mextract        PROTO((int, char **));
-int mmextract        PROTO((int, char **));
+int mmextract       PROTO((int, char **));
+int objectcoverage  PROTO((int, char **));
 int pcat            PROTO((int, char **));
 int photcodes       PROTO((int, char **));
@@ -92 +93 @@
   {1, "mextract",    mextract,     "extract measure data values"},
   {1, "mmextract",   mmextract,    "extract joined measurements"},
+  {1, "objectcoverage", objectcoverage, "plot catalog boundaries"},
   {1, "pcat",        skycat,       "plot catalog boundaries"},
   {1, "photcodes",   photcodes,    "list photometry codes"},

trunk/Ohana/src/photdbc/src/make_subcatalog.c

@@ -71 +71 @@
     }
 
+    // if the input catalog is an old type, generate the catID entries:
+    if (catalog[0].catformat < DVO_FORMAT_PS1_V1) {
+      subcatalog[0].average[Naverage].catID = catalog[0].catID;
+    }
+
     minMag   = 32;
     minSigma = 32;
@@ -127 +132 @@
       subcatalog[0].measure[Nmeasure]        = catalog[0].measure[offset];
       subcatalog[0].measure[Nmeasure].averef = Naverage;
+
+      // if the input catalog is an old type, generate the catID entries:
+      if (catalog[0].catformat < DVO_FORMAT_PS1_V1) {
+        subcatalog[0].measure[Nmeasure].catID = catalog[0].catID;
+      }
+
       Nmeasure ++;
       Nm ++;

trunk/Ohana/src/relastro/src/UpdateObjects.c

@@ -189 +189 @@
       coords.crval1 = R[0];
       coords.crval2 = D[0];
-      Tmean /= (float) N;
+
+      if (FIT_TARGET == TARGET_HIGH_SPEED) {
+          Tmean = 0.5*(Tmax - Tmin);
+      } else {
+          Tmean /= (float) N;
+      }
       
       XVERB = FALSE && (catalog[i].measure[m].dM < 0.01) && (N == 6) && (mode == FIT_PM_ONLY);

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

@@ -11 +11 @@
   unsigned int start;
   unsigned int stop;
+  short photcode;
   float Mcal;
   float dMcal;
@@ -85 +86 @@
 int    RELPHOT_GRID_BINNING;
 
-PhotCode      *photcode;
-int            PhotNsec;
-int            PhotSec;
+int      *photseclist;
+int      Nphotcodes;
+PhotCode **photcodes;
+// int            PhotSec;
+// int            PhotNsec;
 
 PhotCode      *refPhotcode;
+
+int MaxDensityUse;
+double MaxDensityValue;
 
 int AreaSelect;
@@ -213 +219 @@
 StatType      statsMosaicX        PROTO((Catalog *catalog));
 StatType      statsMosaicdM       PROTO((Catalog *catalog));
-StatType      statsStarN          PROTO((Catalog *catalog, int Ncatalog));
-StatType      statsStarS          PROTO((Catalog *catalog, int Ncatalog));
-StatType      statsStarX          PROTO((Catalog *catalog, int Ncatalog));
+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));
 void          wcatalog            PROTO((Catalog *catalog));
 void          wimages             PROTO((void));
@@ -223 +229 @@
 void relphot_usage (void);
 void relphot_help (int argc, char **argv);
+
+off_t getImageByID (off_t ID);
+
+int rationalize_mosaics ();
+int LimitDensityCatalog (Catalog *subcatalog, Catalog *catalog);
+
+int populate_tiny_values (Catalog *catalog);
+int free_tiny_values (Catalog *catalog);

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

@@ -348 +348 @@
   // 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++) {
     
@@ -362 +366 @@
       mx = mlist[i][j];
       c  = clist[i][j];
-      n  = catalog[c].measure[mx].averef;
+      n  = catalog[c].measureT[mx].averef;
       
       // if we have already visited this star, skip the stuff below
@@ -369 +373 @@
 
       // skip stars marked as BAD
-      if (catalog[c].average[n].flags & STAR_BAD) {
+      if (catalog[c].secfilt[n*Nsecfilt+Nsec].flags & STAR_BAD) {
         Nrel ++;
         continue;
@@ -389 +393 @@
 
         // skip measurements marked as BAD
-        if (catalog[c].measure[m].dbFlags & MEAS_BAD) {
+        if (catalog[c].measureT[m].dbFlags & MEAS_BAD) {
           Nbad ++;
           continue;
@@ -410 +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 = PhotSys (&catalog[c].measure[m], &catalog[c].average[n], &catalog[c].secfilt[n*PhotNsec]);
-        Msys = PhotCat (&catalog[c].measure[m]);
+        Msys = PhotCatTiny (&catalog[c].measureT[m]);
         if (isnan(Msys)) {
           Nsys++;
@@ -418 +421 @@
 
         // mag-error for this measurement
-        Merr =  MAX (catalog[c].measure[m].dM, MIN_ERROR);
+        Merr =  MAX (catalog[c].measureT[m].dM, MIN_ERROR);
 
         // Wsys = 1.0 / SQ(Merr);
@@ -528 +531 @@
   if (!USE_GRID) return;
 
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+
   Nmax = Nlist[0];
   for (i = 0; i < Ngrid; i++) {
@@ -545 +550 @@
       c = clist[i][j];
       
-      if (catalog[c].measure[m].dbFlags & MEAS_BAD) {
+      if (catalog[c].measureT[m].dbFlags & MEAS_BAD) {
         Nbad ++;
         continue;
@@ -565 +570 @@
       }
       
-      n = catalog[c].measure[m].averef;
-      Msys = PhotSys (&catalog[c].measure[m], &catalog[c].average[n], &catalog[c].secfilt[n*PhotNsec]);
+      n = catalog[c].measureT[m].averef;
+      Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt]);
       if (isnan(Msys)) {
         Nsys++;
@@ -572 +577 @@
       }
       list[N] = Msys - Mrel - Mcal - Mmos;
-      dlist[N] = MAX (catalog[c].measure[m].dM, MIN_ERROR);
+      dlist[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
       N++;
     }
@@ -613 +618 @@
   if (!USE_GRID) return;
 
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+
   N = 0;
   for (i = 0; i < Ngrid; i++) 
@@ -630 +637 @@
       c = clist[i][j];
       
-      if (catalog[c].measure[m].dbFlags & MEAS_BAD) {
+      if (catalog[c].measureT[m].dbFlags & MEAS_BAD) {
         Narea ++;
         continue;
@@ -641 +648 @@
       if (isnan(Mrel)) continue;
 
-      n = catalog[c].measure[m].averef;
-      Msys = PhotSys (&catalog[c].measure[m], &catalog[c].average[n], &catalog[c].secfilt[n*PhotNsec]);
+      n = catalog[c].measureT[m].averef;
+      Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt]);
 
       xlist[N] = Xmeas[c][m];
@@ -704 +711 @@
   gfits_create_matrix (&header, &matrix);
   gfits_modify (&header, "NEXTEND", OFF_T_FMT, 1,  Nimage + 3);
-  gfits_modify (&header, "FILTER", "%s", 1, photcode[0].name);
+  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");
 
@@ -727 +734 @@
   theader.bitpix   = -32;
   gfits_create_Theader (&theader, "IMAGE");
-  gfits_modify (&theader, "FILTER", "%s", 1, photcode[0].name);
+  gfits_modify (&theader, "FILTER", "%s", 1, photcodes[0][0].name);
   gfits_modify (&theader, "EXTNAME", "%s", 1, "MAG_OFFSET");
   gfits_create_matrix  (&theader, &matrix);
@@ -742 +749 @@
   /* save grid Nmeas values */
   gfits_modify (&theader, "EXTNAME", "%s", 1, "NMEAS");
-  gfits_modify (&theader, "FILTER", "%s", 1, photcode[0].name);
+  gfits_modify (&theader, "FILTER", "%s", 1, photcodes[0][0].name);
   gfits_create_matrix  (&theader, &matrix);
   for (i = 0; i < gridX; i++) {
@@ -756 +763 @@
   /* save grid sigma values */
   gfits_modify (&theader, "EXTNAME", "%s", 1, "SIGMA");
-  gfits_modify (&theader, "FILTER", "%s", 1, photcode[0].name);
+  gfits_modify (&theader, "FILTER", "%s", 1, photcodes[0][0].name);
   gfits_create_matrix  (&theader, &matrix);
   for (i = 0; i < gridX; i++) {
@@ -776 +783 @@
 
     gfits_modify (&theader, "EXTNAME", "%s", 1, camera.ccdname[N]);
-    gfits_modify (&theader, "FILTER", "%s", 1, photcode[0].name);
+    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);

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

@@ -1 +1 @@
 # include "relphot.h"
-
-# define USE_IMAGE_ID 1
 
 static off_t       **bin;     // link from catalog,measure to image          
@@ -12 +10 @@
 static off_t        Nimage;   // number of available images
 
-// if we search by image ID, we sort (imageIDs, imageIdx) by imageIDs to get a sorted
-// index
-
-# if USE_IMAGE_ID
+// to search by image ID, we sort (imageIDs, imageIdx) by imageIDs to get a sorted index
 static off_t        *imageIDs; // list of all image IDs
 static off_t        *imageIdx; // list of index for image IDs 
-# else
-static unsigned int *start;
-static unsigned int *stop;
-# endif
 
 Image *getimages (off_t *N) {
@@ -40 +31 @@
   Nimage = N;
 
-# if USE_IMAGE_ID
   ALLOCATE (imageIDs, off_t, Nimage);
   ALLOCATE (imageIdx, off_t, Nimage);
@@ -49 +39 @@
   }
   llsortpair (imageIDs, imageIdx, Nimage);
-# else
-  ALLOCATE (start,   unsigned, Nimage);
-  ALLOCATE (stop,    unsigned, Nimage);
-
-  for (i = 0; i < Nimage; i++) {
-    start[i] = image[i].tzero - MAX(0.01*image[i].trate*image[i].NY, 1);
-    stop[i]  = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1);
-  }
-# endif
 }
 
@@ -65 +46 @@
   // use bisection to find the specified image ID
 
-# if USE_IMAGE_ID
   off_t Nlo, Nhi, N;
 
@@ -82 +62 @@
       return (imageIdx[N]);
   }
-# endif
 
   return (-1);
@@ -128 +107 @@
 }
 
-/* select all image equivalent to the current photcode */
+/* select all image equivalent to the active photcode set */
 void findImages (Catalog *catalog, int Ncatalog) {
 
   off_t j;
-  int i, ecode;
-
-  for (i = 0; i < Ncatalog; i++) {
+  int i, ecode, Ns, found;
+  
+  int Nmatch = 0;
+ for (i = 0; i < Ncatalog; i++) {
     for (j = 0; j < catalog[i].Nmeasure; j++) {
-      ecode = GetPhotcodeEquivCodebyCode (catalog[i].measure[j].photcode);
-      if (photcode[0].code != ecode) continue;
+      ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[j].photcode);
+      found = FALSE;
+      for (Ns = 0; !found && (Ns < Nphotcodes); Ns++) {
+        if (ecode == photcodes[Ns][0].code) found = TRUE;
+      }
+      if (!found) continue; 
       matchImage (catalog, j, i);
-    }
-  }
-}
-
-int findCCD (off_t idx, off_t meas, int cat, Measure *measure) {
-
-  int ccdnum;
+      Nmatch ++;
+    }
+  }
+ // fprintf (stderr, "matched %d detections to images\n", Nmatch);
+}
+
+int findCCD (off_t idx, off_t meas, int cat, MeasureTiny *measure) {
+
+  int ccdnum, found, Ns, ecode;
   double X, Y;
   char *pname, *filter, *p, base[256];
@@ -151 +137 @@
   /* identify the ccd on the basis of the photcode name */
   pname = GetPhotcodeNamebyCode (image[idx].photcode);
-  filter = photcode[0].name;
+
+  // XXX this seems quite terrible...
+  ecode = GetPhotcodeEquivCodebyCode (measure[0].photcode);
+  found = FALSE;
+  for (Ns = 0; !found && (Ns < Nphotcodes); Ns++) {
+    if (ecode == photcodes[Ns][0].code) found = TRUE;
+  }
+  if (!found) return (FALSE); 
+
+  filter = photcodes[Ns][0].name;
   sprintf (base, "%s.%s.", MOSAICNAME, filter);
   if (strncmp (pname, base, strlen (base))) return (FALSE);
@@ -170 +165 @@
 
   // old code to add this measurement to the grid cell for this chip 
-  // ave = measure[0].averef;
-  // ra  = catalog[cat].average[ave].R - measure[0].dR / 3600.0;
-  // dec = catalog[cat].average[ave].D - measure[0].dD / 3600.0;
+  // 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);
 
@@ -183 +178 @@
 }
 
-# if USE_IMAGE_ID
 void matchImage (Catalog *catalog, off_t meas, int cat) {
 
   off_t idx, ID;
   int status;
-  Measure *measure;
+  MeasureTiny *measure;
   
-  measure = &catalog[cat].measure[meas];
+  measure = &catalog[cat].measureT[meas];
 
   ID = measure[0].imageID;
@@ -226 +220 @@
 }
 
-# else
-// this is the time-based match
-void matchImage (Catalog *catalog, off_t meas, int cat) {
-
-  off_t i;
-  int status;
-  Measure *measure;
-  
-  measure = &catalog[cat].measure[meas];
-  for (i = 0; i < Nimage; i++) {
-    if (image[0].photcode == -1) continue;
-    if (measure[0].photcode != image[i].photcode) continue;
-    if (measure[0].t < start[i]) continue;
-    if (measure[0].t > stop[i]) continue;
-    
-    if (USE_GRID) {
-      status = findCCD(i, meas, cat);
-      if (!status) continue;
-    }
-
-    // index for (catalog, measure) -> image
-    bin[cat][meas] = i;
-
-    // index for image, Nentry -> catalog
-    clist[i][Nlist[i]] = cat;
-
-    // index for image, Nentry -> measure
-    mlist[i][Nlist[i]] = meas;
-    Nlist[i] ++;
-
-    if (Nlist[i] == NLIST[i]) {
-      NLIST[i] += 100;
-      REALLOCATE (clist[i], off_t, NLIST[i]);
-      REALLOCATE (mlist[i], off_t, NLIST[i]);
-    }   
-    return;
-  }
-  /*  fprintf (stderr, "can't find source image for this measurement: %d (%d)\n", measure[0].t, measure[0].photcode); */
-}
-# endif
-
 off_t getImageEntry (off_t meas, int cat) {
 
@@ -307 +260 @@
 
   if (FREEZE_IMAGES) return;
+
+  int Nsecfilt = GetPhotcodeNsecfilt ();
 
   if (PoorImages) {
@@ -339 +294 @@
       c = clist[i][j];
       
-      if (catalog[c].measure[m].dbFlags & MEAS_BAD) {
+      if (catalog[c].measureT[m].dbFlags & MEAS_BAD) {
           Nbad++;
           continue;
@@ -359 +314 @@
       }
       
-      n = catalog[c].measure[m].averef;
-      Msys = PhotSys (&catalog[c].measure[m], &catalog[c].average[n], &catalog[c].secfilt[n*PhotNsec]);
+      n = catalog[c].measureT[m].averef;
+      Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt]);
       if (isnan(Msys)) {
         Nsys++;
@@ -366 +321 @@
       }
       list[N] = Msys - Mrel - Mmos - Mgrid;
-      dlist[N] = MAX (catalog[c].measure[m].dM, MIN_ERROR);
-      if (catalog[c].measure[m].dM < IMFIT_SYS_SIGMA_LIM) {
+      dlist[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
+      if (catalog[c].measureT[m].dM < IMFIT_SYS_SIGMA_LIM) {
         McalBright += list[N];
         McalBright2 += SQ(list[N]);
@@ -632 +587 @@
   return (stats);
 }
-

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

@@ -9 +9 @@
 static off_t   **imlist; /* mosaic -> image[] */
 static off_t   **bin;    /* catalog, measure -> mosaic */
+
+// list of mosaics associated with an image
+static off_t   *mosimage;
 
 // list of mosaic associated with each image  
@@ -33 +36 @@
   ALLOCATE (NIMLIST, off_t,   NMOSAIC);
 
+  ALLOCATE (mosimage, 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
@@ -40 +45 @@
   /* generate list of unique mosaics */
   for (i = 0; i < Nimage; i++) {
+    mosimage[i] = -1;
 
     /* select valid mosaic images by photcode */
@@ -56 +62 @@
       if (start > mosaic[j].stop)  continue;
       found = TRUE;
+
+      // add reference from image to mosaic
+      mosimage[i] = j;
 
       /* add image to mosaic image list */
@@ -77 +86 @@
     mosaic[Nmosaic].flags  = image[i].flags;
     mosaic[Nmosaic].secz  = image[i].secz;
+    mosaic[Nmosaic].photcode = GetPhotcodeEquivCodebyCode (image[i].photcode);
 
     /* add image to mosaic image list */
@@ -84 +94 @@
     imlist[Nmosaic][0] = i;
 
+    // add reference from image to mosaic
+    mosimage[i] = Nmosaic;
+    
     Nmosaic ++;
     if (Nmosaic == NMOSAIC) {
@@ -250 +263 @@
 int findMosaics (Catalog *catalog, int Ncatalog) {
   
-  int i, ecode;
+  int i, ecode, found, Ns;
   off_t j;
 
   if (!MOSAIC_ZEROPT) return (FALSE);
 
+  int Nmatch = 0;
   for (i = 0; i < Ncatalog; i++) {
     for (j = 0; j < catalog[i].Nmeasure; j++) {
       if (TimeSelect) {
-        if (catalog[i].measure[j].t < TSTART) continue;
-        if (catalog[i].measure[j].t > TSTOP) continue;
-      }
-      ecode = GetPhotcodeEquivCodebyCode (catalog[i].measure[j].photcode);
-      if (photcode[0].code != ecode) continue;
+        if (catalog[i].measureT[j].t < TSTART) continue;
+        if (catalog[i].measureT[j].t > TSTOP) continue;
+      }
+      ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[j].photcode);
+      found = FALSE;
+      for (Ns = 0; !found && (Ns < Nphotcodes); Ns++) {
+        if (ecode == photcodes[Ns][0].code) found = TRUE;
+      }
+      if (!found) continue;
       matchMosaics (catalog, j, i);
-    }
-  }
+      Nmatch ++;
+    }
+  }
+  // fprintf (stderr, "matched %d detections to mosaics\n", Nmatch);
   return (TRUE);
 }
@@ -271 +291 @@
 void matchMosaics (Catalog *catalog, off_t meas, int cat) {
 
-  int i;
-
-  for (i = 0; i < Nmosaic; i++) {
-    if (catalog[cat].measure[meas].t < mosaic[i].start) continue;
-    if (catalog[cat].measure[meas].t > mosaic[i].stop) continue;
+  off_t idx, ID, mosID;
+  MeasureTiny *measure;
+
+  measure = &catalog[cat].measureT[meas];
+
+  ID = measure[0].imageID;
+  idx = getImageByID (ID);
+  if (idx == -1) {
+    if (VERBOSE2) fprintf (stderr, "missed measurement "OFF_T_FMT", %d\n", meas, cat);
+    return;
+  }
+
+  mosID = mosimage[idx];
+  if (mosID < 0) {
+    Image *image = getimage(idx);
+    fprintf (stderr, "unmatched image %s\n", image[0].name);
+    return;
+  }
+
+  // test to check we got the right match:
+  {
+    Image *image = getimage(idx);
+    unsigned int imageStart = image[0].tzero - MAX(0.01*image[0].trate*image[0].NY, 1);
+    if (imageStart != mosaic[mosID].start) {
+      fprintf (stderr, "error in image to mosaic match\n");
+      abort();
+    }
+  }
+
+  bin[cat][meas] = mosID;
+
+  clist[mosID][Nlist[mosID]] = cat;
+  mlist[mosID][Nlist[mosID]] = meas;
+  Nlist[mosID] ++;
     
-# ifdef GRID_V1
-    if (USE_GRID) {
-      ave = catalog[cat].measure[meas].averef;
-      ra  = catalog[cat].average[ave].R - catalog[cat].measure[meas].dR / 3600.0;
-      dec = catalog[cat].average[ave].D - catalog[cat].measure[meas].dD / 3600.0;
-
-      /* X,Y always positive-definite in range 0,0 - dX, dY */
-      RD_to_XY (&X, &Y, ra, dec, &mosaic[i].coords);
-      setGridMeasure (meas, cat, X, Y);
-    }
-# endif
-
-    bin[cat][meas] = i;
-
-    clist[i][Nlist[i]] = cat;
-    mlist[i][Nlist[i]] = meas;
-    Nlist[i] ++;
-    
-    if (Nlist[i] == NLIST[i]) {
-      NLIST[i] += 100;
-      REALLOCATE (clist[i], int,   NLIST[i]);
-      REALLOCATE (mlist[i], off_t, NLIST[i]);
-    }   
-    return;
-  }
-  fprintf (stderr, "missed measurement\n");
+  if (Nlist[mosID] == NLIST[mosID]) {
+    NLIST[mosID] += 100;
+    REALLOCATE (clist[mosID], int,   NLIST[mosID]);
+    REALLOCATE (mlist[mosID], off_t, NLIST[mosID]);
+  }     
   return;
 }
@@ -324 +352 @@
 
   off_t i, j, m, c, N, Nmax;
-  int n, mark, bad, Nfew, Nbad, Ncal, Nrel, Ngrid, Nsys;
+  int n, mark, bad, Nfew, Nbad, Ncal, Nrel, Ngrid, Nsys, Nsecfilt;
   float Msys, Mrel, Mcal, Mgrid;
   double *list, *dlist, *Mlist, *dMlist;
@@ -334 +362 @@
 
   image = getimages (&N);
+
+  Nsecfilt = GetPhotcodeNsecfilt ();
 
   if (PoorImages) {
@@ -364 +394 @@
       c = clist[i][j];
       
-      if (catalog[c].measure[m].dbFlags & MEAS_BAD) {
+      if (catalog[c].measureT[m].dbFlags & MEAS_BAD) {
           Nbad ++;
           continue;
@@ -384 +414 @@
       }
       
-      n = catalog[c].measure[m].averef;
-      Msys = PhotSys (&catalog[c].measure[m], &catalog[c].average[n], &catalog[c].secfilt[n*PhotNsec]);
+      n = catalog[c].measureT[m].averef;
+      Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt]);
       if (isnan(Msys)) {
         Nsys++;
@@ -391 +421 @@
       }
       list[N]  = Msys - Mrel - Mcal - Mgrid;
-      dlist[N] = MAX (catalog[c].measure[m].dM, MIN_ERROR);
+      dlist[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
       Mlist[N] = Msys;
       dMlist[N] = list[N];
@@ -433 +463 @@
 }
   
+// 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. 
+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);
+
+  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].Mcal) < CLOUD_TOLERANCE)) {
+      fprintf (stderr, "cloud-free: %s : %f\n", image[imlist[i][0]].name, mosaic[i].Mcal);
+    } 
+    if (VERBOSE2 && (mosaic[i].Mcal < -CLOUD_TOLERANCE)) {
+      imageOffset[i] = -mosaic[i].Mcal; 
+      // NOTE the negative sign: down below, we are going to add in the negative of Mcal
+      // to this image, and the propagated mean values for other images
+      fprintf (stderr, "anti-clouds: %s : %f\n", image[imlist[i][0]].name, mosaic[i].Mcal);
+    } 
+    if (VERBOSE2 && (mosaic[i].Mcal > CLOUD_TOLERANCE)) {
+      fprintf (stderr, "cloudy    : %s : %f\n", image[imlist[i][0]].name, mosaic[i].Mcal);
+    } 
+  }
+
+  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
+        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 = bin[i][m];
+        if (nMos == -1) continue;
+
+        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: Mcal > 0; anti-clouds: Mcal < 0; imageOffset = -Mcal)
+    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++) {
+      Nstars = Slist[i][j];
+      if (starNcount[Nstars][Ns] > 1) {
+        dM += (starOffset[Nstars][Ns] / starNcount[Nstars][Ns]);
+      }
+    }
+    imageOffset[i] = dM / NSlist[i];
+    // fprintf (stderr, "adjust image: %s : (%f %d) : %f\n", image[imlist[i][0]].name, dM, NSlist[i], imageOffset[i]);
+  }
+
+  // for (i = 0; i < Nmosaic; i++) {
+  //   fprintf (stderr, "correction: %s : %f\n", image[imlist[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].Mcal += imageOffset[i];
+  }
+
+  for (i = 0; i < Nstars; i++) {
+    free (starOffset[i]);
+    free (starNcount[i]);
+  }
+  free (starOffset);
+  free (starNcount);
+
+  for (i = 0; i < Nmosaic; i++){
+    free (Slist[i]);
+  }
+  free (NSlist);
+  free (NSLIST);
+  free (imageOffset);
+  free (seclist);
+  return (TRUE);
+}
+
 StatType statsMosaicM (Catalog *catalog) {
 
@@ -524 +748 @@
     n++;
   }
-  fprintf (stderr, "Nmosaic: "OFF_T_FMT", n: "OFF_T_FMT"\n",  Nmosaic,  n);
+  // fprintf (stderr, "Nmosaic: "OFF_T_FMT", n: "OFF_T_FMT"\n",  Nmosaic,  n);
 
   liststats (list, dlist, n, &stats);
@@ -563 +787 @@
 void clean_mosaics () {
 
-  off_t i, N, mark, Nmark;
+  off_t i, N, mark, Nmark, Nscatter, Noffset;
   double *mlist, *slist, *dlist;
   double MaxOffset, MedOffset, MaxScatter;
@@ -592 +816 @@
   fprintf (stderr, "Mrel: %f, dMrel: %f, Max Scatter: %f, Max Offset: %f\n", MedOffset, stats.median, MaxScatter, MaxOffset);
   
-  Nmark = 0;
+  Nmark = Nscatter = Noffset = 0;
   for (i = 0; i < Nmosaic; i++) {
     mark = FALSE;
-    mark = (mosaic[i].dMcal > MaxScatter) || (fabs(mosaic[i].Mcal - MedOffset) > MaxOffset);
+    if (mosaic[i].dMcal > MaxScatter) {
+      mark = TRUE;
+      Nscatter ++;
+    }
+    if (fabs(mosaic[i].Mcal - MedOffset) > MaxOffset) {
+      mark = TRUE;
+      Noffset ++;
+    }
     if (mark) { 
       Nmark ++;
@@ -604 +835 @@
   }
 
-  fprintf (stderr, OFF_T_FMT" mosaics marked poor\n",  Nmark);
+  fprintf (stderr, OFF_T_FMT" mosaics marked poor ("OFF_T_FMT" scatter, "OFF_T_FMT" offset)\n",  Nmark, Nscatter, Noffset);
   initstats (STATMODE);
   free (mlist);
@@ -640 +871 @@
       c = clist[i][j];
       
-      if (catalog[c].measure[m].dbFlags & (ID_MEAS_AREA | ID_MEAS_NOCAL)) continue;
-
-      ave = catalog[c].measure[m].averef;
-      xlist[N] = catalog[c].average[ave].R - catalog[c].measure[m].dR / 3600.0;
-      ylist[N] = catalog[c].average[ave].D - catalog[c].measure[m].dD / 3600.0;
+      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;
       N++;
     }

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

@@ -5 +5 @@
 static double *dlist;
 
+// 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.  This array carries the impact of those offsets
+// on each star
+static double *Moffset;
+
 void initMrel (Catalog *catalog, int Ncatalog) {
 
@@ -12 +18 @@
   for (i = 0; i < Ncatalog; i++) {
     for (j = 0; j < catalog[i].Naverage; j++) {
-      Nmax = MAX (Nmax, catalog[i].average[j].Nmeasure);
-    }
-  }
-
-  ALLOCATE (list, double, MAX (1, Nmax));
-  ALLOCATE (dlist, double, MAX (1, Nmax));
+      Nmax = MAX (Nmax, catalog[i].averageT[j].Nmeasure);
+    }
+  }
+
+  ALLOCATE (list,    double, MAX (1, Nmax));
+  ALLOCATE (dlist,   double, MAX (1, Nmax));
+  ALLOCATE (Moffset, double, MAX (1, Nmax));
 }  
 
 float getMrel (Catalog *catalog, off_t meas, int cat) {
 
+  int Nsec, Nsecfilt, photcode;
   int ave;
   float value;
 
-  ave = catalog[cat].measure[meas].averef;
-  if (catalog[cat].average[ave].flags & STAR_BAD) return (NAN);  
+  ave = catalog[cat].measureT[meas].averef;
+  photcode = catalog[cat].measureT[meas].photcode;
+
+  int ecode = GetPhotcodeEquivCodebyCode (photcode);
+  Nsec = GetPhotcodeNsec(ecode);
+  Nsecfilt = GetPhotcodeNsecfilt ();
+
+  // is this star OK?
+  if (catalog[cat].secfilt[Nsecfilt*ave+Nsec].flags & STAR_BAD) return (NAN);  
  
-  value = catalog[cat].secfilt[PhotNsec*ave+PhotSec].M;
+  value = catalog[cat].secfilt[Nsecfilt*ave+Nsec].M;
   return (value);
 }
@@ -39 +54 @@
   StatType stats;
 
+  int Nsecfilt = GetPhotcodeNsecfilt ();
   Nfew = Nsys = Nbad = Ncal = Nmos = Ngrid = 0;
 
@@ -44 +60 @@
     for (j = 0; j < catalog[i].Naverage; j++) {
 
-      /* calculate the average value for a single star */
-      if (catalog[i].average[j].flags & STAR_BAD) continue;  
-      m = catalog[i].average[j].measureOffset;
-
-      N = 0;
-      for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
-        if (catalog[i].measure[m].dbFlags & MEAS_BAD) {
-          Nbad ++;
-          continue;
-        }
-        // XXX allow REF stars (no Image Entry) to be included in the calculation this
-        // should be optionally set, and should allow for REF stars to be downweighted by
-        // more than their reported errors.  how such info is carried is unclear...
-        if (getImageEntry (m, i) < 0) {
-          Mcal = Mmos = Mgrid = 0;
-        } else {
-          Mcal  = getMcal  (m, i);
-          if (isnan(Mcal)) { 
-            Ncal ++;
+      int Ns;
+      for (Ns = 0; Ns < Nphotcodes; Ns++) {
+
+        int thisCode = photcodes[Ns][0].code;
+        int Nsec = GetPhotcodeNsec(thisCode);
+
+        /* calculate the average value for a single star */
+
+        // skip bad stars
+        if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;
+        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; }
+
+          if (catalog[i].measureT[m].dbFlags & MEAS_BAD) {
+            Nbad ++;
             continue;
           }
-          Mmos  = getMmos  (m, i);
-          if (isnan(Mmos)) {
-            Nmos ++;
+          // XXX allow REF stars (no Image Entry) to be included in the calculation this
+          // should be optionally set, and should allow for REF stars to be downweighted by
+          // more than their reported errors.  how such info is carried is unclear...
+          if (getImageEntry (m, i) < 0) {
+            Mcal = Mmos = Mgrid = 0;
+          } else {
+            Mcal  = getMcal  (m, i);
+            if (isnan(Mcal)) { 
+              Ncal ++;
+              continue;
+            }
+            Mmos  = getMmos  (m, i);
+            if (isnan(Mmos)) {
+              Nmos ++;
+              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]);
+          if (isnan(Msys)) {
+            Nsys++;
             continue;
           }
-          Mgrid = getMgrid (m, i);
-          if (isnan(Mgrid)) {
-            Ngrid++;
-            continue;
+          list[N] = Msys - Mcal - Mmos - Mgrid;
+          dlist[N] = MAX (catalog[i].measureT[m].dM, MIN_ERROR);
+
+          // tie down reference photometry if the -refcode (code) option is selected
+          if (refPhotcode) {
+            if (GetPhotcodeEquivCodebyCode(catalog[i].measureT[m].photcode) == refPhotcode[0].equiv) {
+              // increase the weight by a factor of 100:
+              dlist[N] = 0.01*catalog[i].measureT[m].dM;
+            }
           }
-        }
-
-        Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
-        if (isnan(Msys)) {
-          Nsys++;
-          continue;
-        }
-        list[N] = Msys - Mcal - Mmos - Mgrid;
-        dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
-
-        if (refPhotcode) {
-          if (GetPhotcodeEquivCodebyCode(catalog[i].measure[m].photcode) == refPhotcode[0].equiv) {
-            // increase the weight by a factor of 100:
-            dlist[N] = 0.01*catalog[i].measure[m].dM;
-          }
-        }
-        N++;
-      }
-
-      // when performing the grid analysis, STAR_TOOFEW will be set to 1;
-
-      if (N <= STAR_TOOFEW) { /* too few measurements */
-        catalog[i].average[j].flags |= ID_STAR_FEW;
-        Nfew ++;
-      } else {
-        catalog[i].average[j].flags &= ~ID_STAR_FEW;
-      } 
-
-      liststats (list, dlist, N, &stats);
-
-      catalog[i].secfilt[PhotNsec*j+PhotSec].M  = stats.mean;
-      catalog[i].secfilt[PhotNsec*j+PhotSec].dM = stats.sigma;
-      catalog[i].secfilt[PhotNsec*j+PhotSec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NAN_S_SHORT;
+          N++;
+        }
+
+        // when performing the grid analysis, STAR_TOOFEW will be set to 1;
+        if (N <= STAR_TOOFEW) { /* too few measurements */
+          catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_STAR_FEW;
+          Nfew ++;
+        } else {
+          catalog[i].secfilt[Nsecfilt*j+Nsec].flags &= ~ID_STAR_FEW;
+        }       
+
+        liststats (list, dlist, N, &stats);
+
+        catalog[i].secfilt[Nsecfilt*j+Nsec].M  = stats.mean;
+        catalog[i].secfilt[Nsecfilt*j+Nsec].dM = stats.sigma;
+        catalog[i].secfilt[Nsecfilt*j+Nsec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NAN_S_SHORT;
+      }
     }
   }
@@ -122 +151 @@
   double *list, *dlist;
   StatType stats;
+  int Nsec, Nsecfilt, ecode;
+
+  Nsecfilt = GetPhotcodeNsecfilt ();
 
   /* Nmeasure is now different, need to reallocate */
@@ -127 +159 @@
   for (i = 0; i < Ncatalog; i++) {
     for (j = 0; j < catalog[i].Naverage; j++) {
-      Nmax = MAX (Nmax, catalog[i].average[j].Nmeasure);
+      Nmax = MAX (Nmax, catalog[i].averageT[j].Nmeasure);
     }
   }
@@ -135 +167 @@
   for (i = 0; i < Ncatalog; i++) {
     for (j = 0; j < catalog[i].Naverage; j++) {
-
       /* skip stars already calibrated */
       if (catalog[i].found[j]) continue;  
 
-      N = 0;
-      m = catalog[i].average[j].measureOffset;
-      for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
-        if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue;
-        // XXX allow REF stars (no Image Entry) to be included in the calculation this
-        // should be optionally set, and should allow for REF stars to be downweighted by
-        // more than their reported errors.  how such info is carried is unclear...
-        if (getImageEntry (m, i) < 0) {
-          Mcal = Mmos = Mgrid = 0;
-        } else {
-          Mcal  = getMcal  (m, i);
-          if (isnan(Mcal)) continue;
-          Mmos  = getMmos  (m, i);
-          if (isnan(Mmos)) continue;
-          Mgrid = getMgrid (m, i);
-          if (isnan(Mgrid)) continue;
-        }
-
-        Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
-        list[N] = Msys - Mcal - Mmos - Mgrid;
-        dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
-        N++;
-      }
-      if (N < 1) continue;
-
-      liststats (list, dlist, N, &stats);
-      if (mark) catalog[i].found[j] = TRUE;
-
-      /* use sigma or error in dM for output? */
-      catalog[i].secfilt[PhotNsec*j+PhotSec].M  = stats.mean;
-      catalog[i].secfilt[PhotNsec*j+PhotSec].dM = MAX (stats.error, stats.sigma);
-      catalog[i].secfilt[PhotNsec*j+PhotSec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NAN_S_SHORT;
+      int Ns;
+      for (Ns = 0; Ns < Nphotcodes; Ns++) {
+        int thisCode = photcodes[Ns][0].code;
+        Nsec = GetPhotcodeNsec(thisCode);
+
+        N = 0;
+        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
+          ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode);
+          if (ecode != thisCode) { continue; }
+
+          if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue;
+
+          // XXX allow REF stars (no Image Entry) to be included in the calculation this
+          // should be optionally set, and should allow for REF stars to be downweighted by
+          // more than their reported errors.  how such info is carried is unclear...
+          if (getImageEntry (m, i) < 0) {
+            Mcal = Mmos = Mgrid = 0;
+          } else {
+            Mcal  = getMcal  (m, i);
+            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]);
+          list[N] = Msys - Mcal - Mmos - Mgrid;
+          dlist[N] = MAX (catalog[i].measureT[m].dM, MIN_ERROR);
+          N++;
+        }
+        if (N < 1) continue;
+
+        liststats (list, dlist, N, &stats);
+        if (mark) catalog[i].found[j] = TRUE;
+
+        /* use sigma or error in dM for output? */
+        catalog[i].secfilt[Nsecfilt*j+Nsec].M  = stats.mean;
+        catalog[i].secfilt[Nsecfilt*j+Nsec].dM = MAX (stats.error, stats.sigma);
+        catalog[i].secfilt[Nsecfilt*j+Nsec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NAN_S_SHORT;
+      }
     }
   }
@@ -179 +221 @@
 }
 
-// for each average object, set the average mags based on existing equiv photometry
+// For each average object, set the average mags based on existing equiv photometry.
+// NOTE: this function operates on the real Measure & Average structures, not the
+// MeasureTiny & AverageTiny structures
 int setMave (Catalog *catalog, int Ncatalog) {
 
@@ -202 +246 @@
   ALLOCATE (dlist, double, MAX (1, Nmax));
 
-  Nsecfilt = catalog[0].Nsecfilt;
+  Nsecfilt = GetPhotcodeNsecfilt ();
 
 # define PSFQUALSTATS 1
@@ -220 +264 @@
           if (GetPhotcodeEquivCodebyCode (catalog[i].measure[m].photcode) != Nc) continue;
 
-          Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
+          Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
           if (isnan(Msys)) continue;
 
@@ -321 +365 @@
     for (j = 0; j < catalog[i].Naverage; j++) {
 
-      m = catalog[i].average[j].measureOffset;
-      for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
-        if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue;
+      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;
         Mcal  = getMcal  (m, i);
         if (isnan(Mcal)) continue;
@@ -330 +374 @@
         Mgrid = getMgrid (m, i);
         if (isnan(Mgrid)) continue;
+
+        // set the output calibration
         catalog[i].measure[m].Mcal = Mcal + Mmos + Mgrid;
       }
@@ -339 +385 @@
 void clean_stars (Catalog *catalog, int Ncatalog) {
 
-  int i, j, Ndel, Nave, Ntot, mark;
+  int i, j, Ndel, Nave, Ntot, mark, Ns;
   float dM, Xm;
   double Chisq, MaxScatter, MaxChisq;
@@ -354 +400 @@
   ALLOCATE (slist, double, Ntot);
   ALLOCATE (dlist, double, Ntot);
-  for (i = Ntot = 0; i < Ncatalog; i++) {
-    for (j = 0; j < catalog[i].Naverage; j++) {
-      if (catalog[i].average[j].flags & STAR_BAD) continue;
-      Xm = catalog[i].secfilt[PhotNsec*j+PhotSec].Xm;
-      if (Xm == -1) continue;
-      Chisq = pow (10.0, 0.01*Xm);
-      xlist[Ntot] = Chisq;
-      slist[Ntot] = catalog[i].secfilt[PhotNsec*j+PhotSec].dM;
-      dlist[Ntot] = 1;
-      Ntot ++;
-    }
-  }
+
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+
+  // eliminate bad stars using the stats for a single secfilt at a time
+  // XXX DEP replace average.flags with secfilt flags
+  for (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;
+        Xm = catalog[i].secfilt[Nsecfilt*j+Nsec].Xm;
+        if (Xm == -1) continue;
+        Chisq = pow (10.0, 0.01*Xm);
+        xlist[Ntot] = Chisq;
+        slist[Ntot] = catalog[i].secfilt[Nsecfilt*j+Nsec].dM;
+        dlist[Ntot] = 1;
+        Ntot ++;
+      }
+    }
   
-  initstats ("MEAN");
-  liststats (xlist, dlist, Ntot, &stats);
-  MaxChisq = MAX (STAR_CHISQ, 2*stats.median);
-  liststats (slist, dlist, Ntot, &stats);
-  MaxScatter = MAX (STAR_SCATTER, 2*stats.median);
-  fprintf (stderr, "Max Scatter: %f, Max Chisq: %f\n", MaxScatter, MaxChisq);
-
-  Ndel = Nave = 0;
-  for (i = 0; i < Ncatalog; i++) {
-    for (j = 0; j < catalog[i].Naverage; j++) {
-      dM = catalog[i].secfilt[PhotNsec*j+PhotSec].dM;
-      Xm = catalog[i].secfilt[PhotNsec*j+PhotSec].Xm;
-      Chisq = pow (10.0, 0.01*Xm);
-      mark = (dM > MaxScatter) || (Xm == NAN_S_SHORT) || (Chisq > MaxChisq);
-      if (mark) {
-        catalog[i].average[j].flags |= ID_STAR_POOR;
-        Ndel ++;
-      } else {
-        catalog[i].average[j].flags &= ~ID_STAR_POOR;
-      }
-      Nave ++;
-    }
-  }
-  fprintf (stderr, "%d stars marked variable, %d total\n", Ndel, Nave);
-  initstats (STATMODE);
+    initstats ("MEAN");
+    liststats (xlist, dlist, Ntot, &stats);
+    MaxChisq = MAX (STAR_CHISQ, 2*stats.median);
+    liststats (slist, dlist, Ntot, &stats);
+    MaxScatter = MAX (STAR_SCATTER, 2*stats.median);
+    fprintf (stderr, "Max Scatter: %f, Max Chisq: %f\n", MaxScatter, MaxChisq);
+
+    Ndel = Nave = 0;
+    for (i = 0; i < Ncatalog; i++) {
+      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);
+        if (mark) {
+          catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_STAR_POOR;
+          Ndel ++;
+        } else {
+          catalog[i].secfilt[Nsecfilt*j+Nsec].flags &= ~ID_STAR_POOR;
+        }
+        Nave ++;
+      }
+    }
+    fprintf (stderr, "%d stars marked variable, %d total\n", Ndel, Nave);
+    initstats (STATMODE);
+  }
   free (xlist);
   free (slist);
@@ -409 +466 @@
   int Ncal, Nmos, Ngrid, Nfew;
 
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+
   if (VERBOSE) fprintf (stderr, "marking poor measures\n");
   /* Nmeasure is now different, need to reallocate */
@@ -414 +473 @@
   for (i = 0; i < Ncatalog; i++) {
     for (j = 0; j < catalog[i].Naverage; j++) {
-      Nmax = MAX (Nmax, catalog[i].average[j].Nmeasure);
+      Nmax = MAX (Nmax, catalog[i].averageT[j].Nmeasure);
     }
   }
@@ -430 +489 @@
     for (j = 0; j < catalog[i].Naverage; j++) {
 
-      /* skip bad stars to prevent them from becoming good (on inner sample) */
-      if (catalog[i].average[j].flags & STAR_BAD) continue;  
-
-      /* on final processing, skip stars already measured */
-      if (final && catalog[i].found[j]) continue;  
-
-      /* accumulate list of valid measurements */
-      m = catalog[i].average[j].measureOffset;
-      N = 0;
-      for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
-        /* if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue; */
-        Mcal  = getMcal  (m, i);
-        if (isnan(Mcal)) { Ncal ++; continue; }
-        Mmos  = getMmos  (m, i);
-        if (isnan(Mmos)) { Nmos ++; continue; }
-        Mgrid = getMgrid (m, i);
-        if (isnan(Mgrid)) { Ngrid ++; continue; }
-
-        Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
-        list[N] = Msys - Mcal - Mmos - Mgrid;
-        dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
-        N++;
-      }
-      if (N <= TOOFEW) { Nfew ++; continue; }
-
-      /* 3-sigma clip based on stats of inner 50% */
-
-      // calculated mean of inner 50%
-      initstats ("INNER_MEAN");
-      liststats (list, dlist, N, &stats);
-      stats.sigma = MAX (MIN_ERROR, stats.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] - stats.median) < NSIGMA_CLIP*stats.sigma) {
-          list[m] = list[k];
-          m++;
-        }
-      }
-      // recalculate the mean & sigma of the accepted measurements
-      initstats ("MEAN");
-      liststats (list, dlist, 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].average[j].measureOffset;
-      N = 0;
-      for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
-        /* if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue; */
-        Mcal  = getMcal  (m, i);
-        if (isnan(Mcal)) continue;
-        Mmos  = getMmos  (m, i);
-        if (isnan(Mmos)) continue;
-        Mgrid = getMgrid (m, i);
-        if (isnan(Mgrid)) continue;
-
-        Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
-        list[N] = Msys - Mcal - Mmos - Mgrid;
-        dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
-        ilist[N] = m;
-        N++;
-        Nave ++;
-      }
-      if (N < TOOFEW) continue;
-
-      /* mark bad measures (> 3 sigma deviant) */
-      for (k = 0; k < N; k++) {
-        if (fabs (list[k] - stats.median) > NSIGMA_REJECT*stats.sigma) {
-          catalog[i].measure[ilist[k]].dbFlags |= ID_MEAS_POOR_PHOTOM;
-          Ndel ++;
-        }
-      }
-      IMAGE_BAD = image_bad;
+      int Ns;
+      for (Ns = 0; Ns < Nphotcodes; Ns++) {
+        
+        /* on final processing, skip stars already measured */
+        if (final && catalog[i].found[j]) continue;  
+
+        int thisCode = photcodes[Ns][0].code;
+        int Nsec = GetPhotcodeNsec(thisCode);
+        
+        /* 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; }
+
+          /* if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue; */
+          Mcal  = getMcal  (m, i);
+          if (isnan(Mcal)) { Ncal ++; continue; }
+          Mmos  = getMmos  (m, i);
+          if (isnan(Mmos)) { Nmos ++; 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]);
+          list[N] = Msys - Mcal - Mmos - 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%
+        initstats ("INNER_MEAN");
+        liststats (list, dlist, N, &stats);
+        stats.sigma = MAX (MIN_ERROR, stats.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] - stats.median) < NSIGMA_CLIP*stats.sigma) {
+            list[m] = list[k];
+            m++;
+          }
+        }
+        // recalculate the mean & sigma of the accepted measurements
+        initstats ("MEAN");
+        liststats (list, dlist, 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; }
+
+          /* if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue; */
+          Mcal  = getMcal  (m, i);
+          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]);
+          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 (> 3 sigma deviant) */
+        for (k = 0; k < N; k++) {
+          if (fabs (list[k] - stats.median) > NSIGMA_REJECT*stats.sigma) {
+            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 ++;
+          }
+        }
+        IMAGE_BAD = image_bad;
+      }
     }
   }
@@ -513 +591 @@
 }
 
-StatType statsStarN (Catalog *catalog, int Ncatalog) {
+StatType statsStarN (Catalog *catalog, int Ncatalog, int Nsec, int seccode) {
 
   off_t j, k, m, Ntot;
@@ -520 +598 @@
   float Mcal, Mmos, Mgrid;
   StatType stats;
+  // int N1, N2, N3, N4, N0;
+  // N1 = N2 = N3 = N4 = N0 = 0;
+
+  int Nsecfilt = GetPhotcodeNsecfilt ();
 
   Ntot = 0;
@@ -534 +616 @@
 
       /* calculate the average value for a single star */
-      if (catalog[i].average[j].flags & STAR_BAD) continue;  
-      m = catalog[i].average[j].measureOffset;
+      if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) { continue;  }
+      m = catalog[i].averageT[j].measureOffset;
 
       N = 0;
-      for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
+      for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) {
+        int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode);
+        if (ecode != seccode) { continue;}
         Mcal = getMcal  (m, i);
-        if (isnan(Mcal)) continue;
+        if (isnan(Mcal)) { continue;}
         Mmos = getMmos  (m, i);
-        if (isnan(Mmos)) continue;
+        if (isnan(Mmos)) { continue; }
         Mgrid = getMgrid (m, i);
-        if (isnan(Mgrid)) continue;
+        if (isnan(Mgrid)) { continue;}
         N++;
       }
@@ -554 +638 @@
   }
 
+  // fprintf (stderr, "N1: %d, N2: %d, N3: %d, N4: %d, N0: %d\n", N1, N2, N3, N4, N0);
   liststats (list, dlist, n, &stats);
   free (list);
@@ -560 +645 @@
 }
 
-StatType statsStarX (Catalog *catalog, int Ncatalog) {
+// stats for a single secfilt at a time
+StatType statsStarX (Catalog *catalog, int Ncatalog, int Nsec) {
 
   off_t j, Ntot;
@@ -567 +653 @@
   StatType stats;
 
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+
   Ntot = 0;
   for (i = 0; i < Ncatalog; i++) {
@@ -580 +668 @@
 
       /* calculate the average value for a single star */
-      if (catalog[i].average[j].flags & STAR_BAD) continue;  
-
-      Xm = catalog[i].secfilt[PhotNsec*j+PhotSec].Xm;
+      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);
@@ -596 +684 @@
 }
 
-StatType statsStarS (Catalog *catalog, int Ncatalog) {
+StatType statsStarS (Catalog *catalog, int Ncatalog, int Nsec) {
 
   int i, n;
@@ -604 +692 @@
   StatType stats;
 
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+
   Ntot = 0;
   for (i = 0; i < Ncatalog; i++) {
@@ -617 +707 @@
 
       /* calculate the average value for a single star */
-      if (catalog[i].average[j].flags & STAR_BAD) continue;  
-
-      dM = catalog[i].secfilt[PhotNsec*j+PhotSec].dM;
+      if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;  
+
+      dM = catalog[i].secfilt[Nsecfilt*j+Nsec].dM;
       list[n] = dM;
       dlist[n] = 1;
@@ -644 +734 @@
   ALLOCATE (Mlist, double, NBIN);
 
-  for (i = 0; i < NBIN; i++) xlist[i] = 0.00025*i;
-  bzero (Mlist, NBIN*sizeof(double));
-  for (i = 0; i < Ncatalog; i++) {
-    for (j = 0; j < catalog[i].Naverage; j++) {
-      if (catalog[i].average[j].flags & STAR_BAD) continue;  
-      dMrel = catalog[i].secfilt[PhotNsec*j+PhotSec].dM;
-      bin = dMrel / 0.00025;
-      bin = MAX (0, MIN (NBIN-1, bin));
-      Mlist[bin] += 1.0;
-    }
-  }
-
-  plot_defaults (&graphdata);
-  graphdata.style = 1;
-  plot_list (&graphdata, xlist, Mlist, NBIN, "dMrel hist", "%s.dMhist.png", OUTROOT);
-
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+
+  int Ns;
+  for (Ns = 0; Ns < Nphotcodes; Ns++) {
+
+    int thisCode = photcodes[Ns][0].code;
+    int Nsec = GetPhotcodeNsec(thisCode);
+
+    for (i = 0; i < NBIN; i++) xlist[i] = 0.00025*i;
+    bzero (Mlist, NBIN*sizeof(double));
+    for (i = 0; i < Ncatalog; i++) {
+      for (j = 0; j < catalog[i].Naverage; j++) {
+        if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;  
+        dMrel = catalog[i].secfilt[Nsecfilt*j+Nsec].dM;
+        bin = dMrel / 0.00025;
+        bin = MAX (0, MIN (NBIN-1, bin));
+        Mlist[bin] += 1.0;
+      }
+    }
+
+    plot_defaults (&graphdata);
+    graphdata.style = 1;
+    plot_list (&graphdata, xlist, Mlist, NBIN, "dMrel hist", "%s.dMhist.png", OUTROOT);
+  }
   free (xlist);
   free (Mlist);
@@ -671 +769 @@
   Graphdata graphdata;
 
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+
   Ntotal = 0;
   for (i = 0; i < Ncatalog; i++) Ntotal += catalog[i].Naverage;
@@ -677 +777 @@
   ALLOCATE (ylist, double, Ntotal);
 
-  N = 0;
-  for (i = 0; i < Ncatalog; i++) {
-    for (j = 0; j < catalog[i].Naverage; j++) {
-      if (catalog[i].average[j].flags & STAR_BAD) continue;
-      xlist[N] = catalog[i].secfilt[PhotNsec*j+PhotSec].M;
-      value    = catalog[i].secfilt[PhotNsec*j+PhotSec].Xm;
-      if (value == NAN_S_SHORT) continue;
-      ylist[N] = 0.01*value;
-      N++;
-    }
-  }
-
-  plot_defaults (&graphdata);
-  graphdata.ymin = -3.0;
-  plot_list (&graphdata, xlist, ylist, N, "chisq", "%s.chisq.png", OUTROOT);
+  int Ns;
+  for (Ns = 0; Ns < Nphotcodes; Ns++) {
+
+    int thisCode = photcodes[Ns][0].code;
+    int Nsec = GetPhotcodeNsec(thisCode);
+
+    N = 0;
+    for (i = 0; i < Ncatalog; i++) {
+      for (j = 0; j < catalog[i].Naverage; j++) {
+        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;
+        N++;
+      }
+    }
+
+    plot_defaults (&graphdata);
+    graphdata.ymin = -3.0;
+    plot_list (&graphdata, xlist, ylist, N, "chisq", "%s.chisq.png", OUTROOT);
+  }
+
   free (xlist);
   free (ylist);
@@ -716 +824 @@
   for (i = 0; i < Ncatalog; i++) {
     for (j = 0; j < catalog[i].Naverage; j++) {
-      xlist[N] = catalog[i].average[j].R;
-      ylist[N] = catalog[i].average[j].D;
+      xlist[N] = catalog[i].averageT[j].R;
+      ylist[N] = catalog[i].averageT[j].D;
       N++;
     }

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

@@ -128 +128 @@
   }
 
+  MaxDensityUse = FALSE;
+  if ((N = get_argument (argc, argv, "-max-density"))) {
+    remove_argument (N, &argc, argv);
+    MaxDensityValue = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+    MaxDensityUse = TRUE;
+  }
+
   SHOW_PARAMS = FALSE;
   if ((N = get_argument (argc, argv, "-params"))) {

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

@@ -1 +1 @@
 # include "relphot.h"
+
+extern double drand48();
+
+int CopyAverageTiny (AverageTiny *averageT, Average *average) {
+
+  averageT[0].R     = average[0].R;
+  averageT[0].D     = average[0].D;
+  averageT[0].flags = average[0].flags;
+  averageT[0].Nmeasure      = average[0].Nmeasure;
+  averageT[0].measureOffset = average[0].measureOffset;
+
+  // make Nmeasure & measureOffset optional?
+
+  return (TRUE);
+}
+
+int CopyMeasureTiny (MeasureTiny *measureT, Measure *measure) {
+
+  measureT[0].dR       = measure[0].dR;
+  measureT[0].dD       = measure[0].dD;
+  measureT[0].M        = measure[0].M;
+  measureT[0].Mcal     = measure[0].Mcal;
+  measureT[0].dM       = measure[0].dM;
+  measureT[0].airmass  = measure[0].airmass;
+  measureT[0].Xccd     = measure[0].Xccd;
+  measureT[0].Yccd     = measure[0].Yccd;
+  measureT[0].t        = measure[0].t;
+  measureT[0].dt       = measure[0].dt;
+  measureT[0].averef   = measure[0].averef;
+  measureT[0].imageID  = measure[0].imageID;
+  measureT[0].dbFlags  = measure[0].dbFlags;
+  measureT[0].photcode = measure[0].photcode;
+
+  return (TRUE);
+}
 
 int bcatalog (Catalog *subcatalog, Catalog *catalog) {
   
   off_t i, j, offset;
-  int ecode;
+  int ecode, found, Ns;
   off_t NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm;
   float mag;
   int Ncode, Ntime, Ndophot, Nmag, Nsigma, Nimag, Nfew, Ngalaxy, Npsfqf;
 
-  // XXX PhotNsec as a global is a bad idea; either get it from catalog
-  // or get it from:
-  // Nsecfilt = GetPhotcodeNsecfilt ();
-  // assert (catalog[0].Nsecfilt == Nsecfilt);
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+  assert (Nsecfilt == catalog[0].Nsecfilt);
 
   /* we are moving only the subset of measurements from catalog[0] to subcatalog[0] */
@@ -18 +51 @@
   NAVERAGE = 50;
   NMEASURE = 1000;
-  ALLOCATE (subcatalog[0].average, Average, NAVERAGE);
-  ALLOCATE (subcatalog[0].secfilt, SecFilt, NAVERAGE*PhotNsec);
-  ALLOCATE (subcatalog[0].measure, Measure, NMEASURE);
+  ALLOCATE (subcatalog[0].measureT, MeasureTiny, NMEASURE);
+  ALLOCATE (subcatalog[0].averageT, AverageTiny, NAVERAGE);
+  ALLOCATE (subcatalog[0].secfilt, SecFilt, NAVERAGE*Nsecfilt);
   Nmeasure = Naverage = 0;
 
@@ -30 +63 @@
 
     /* start with all stars good */
-    subcatalog[0].average[Naverage] = catalog[0].average[i];
-    subcatalog[0].average[Naverage].measureOffset = Nmeasure;
-    for (j = 0; j < PhotNsec; j++) {
-      subcatalog[0].secfilt[PhotNsec*Naverage+j] = catalog[0].secfilt[PhotNsec*i+j];
+    CopyAverageTiny (&subcatalog[0].averageT[Naverage], &catalog[0].average[i]);
+    subcatalog[0].averageT[Naverage].measureOffset = Nmeasure;
+
+    for (j = 0; j < Nsecfilt; j++) {
+      subcatalog[0].secfilt[Nsecfilt*Naverage+j] = catalog[0].secfilt[Nsecfilt*i+j];
     }
 
     if (RESET) {
-      subcatalog[0].secfilt[PhotNsec*Naverage+PhotSec].M  = NAN;
-      subcatalog[0].secfilt[PhotNsec*Naverage+PhotSec].dM = NAN;
-      subcatalog[0].average[Naverage].flags &= ~ID_STAR_FEW;
-      subcatalog[0].average[Naverage].flags &= ~ID_STAR_POOR;
+      int Ns;
+      for (Ns = 0; Ns < Nphotcodes; Ns++) {
+
+        int thisCode = photcodes[Ns][0].code;
+        int Nsec = GetPhotcodeNsec(thisCode);
+
+        subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].M  = NAN;
+        subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].dM = NAN;
+        subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].flags &= ~ID_STAR_FEW;
+        subcatalog[0].secfilt[Nsecfilt*Naverage+Nsec].flags &= ~ID_STAR_POOR;
+      }
     }
 
@@ -52 +93 @@
       /* select measurements by photcode */
       ecode = GetPhotcodeEquivCodebyCode (catalog[0].measure[offset].photcode);
-      if (ecode != photcode[0].code) { Ncode ++; continue; }
+      found = FALSE;
+      for (Ns = 0; !found && (Ns < Nphotcodes); Ns++) {
+        if (ecode == photcodes[Ns][0].code) found = TRUE;
+      }
+      if (!found) {
+        Ncode ++; 
+        continue; 
+      }
 
       /* select measurements by time */
@@ -61 +109 @@
 
       /* select measurements by quality */
-      // XXX ignore this criterion for REF measurements?
-      // XXX chnage this to select by bitflags
       if (DophotSelect && ((catalog[0].measure[offset].photFlags >> 16) != DophotValue)) { Ndophot ++; continue; }
 
@@ -70 +116 @@
       // check for galaxies
       if (!isnan(catalog[0].measure[offset].Map)) {
-          if (catalog[0].measure[offset].M - catalog[0].measure[offset].Map > 0.15) {
-              nEXT ++;
-          } else {
-              nPSF ++;
-          }
+        if (catalog[0].measure[offset].M - catalog[0].measure[offset].Map > 0.15) {
+          nEXT ++;
+        } else {
+          nPSF ++;
+        }
       }
 
@@ -91 +137 @@
       }
 
-      subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_SKIP_PHOTOM;
-      subcatalog[0].measure[Nmeasure]        = catalog[0].measure[offset];
-      subcatalog[0].measure[Nmeasure].averef = Naverage;
+      CopyMeasureTiny (&subcatalog[0].measureT[Nmeasure], &catalog[0].measure[offset]);
+      subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_SKIP_PHOTOM;
+      subcatalog[0].measureT[Nmeasure].averef = Naverage;
       if (RESET) { 
-        subcatalog[0].measure[Nmeasure].Mcal = 0;
-        subcatalog[0].measure[Nmeasure].dbFlags &= 0xff00;
-        subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_POOR_PHOTOM;
-        subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_AREA;
-        subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_NOCAL;
+        subcatalog[0].measureT[Nmeasure].Mcal = 0;
+        subcatalog[0].measureT[Nmeasure].dbFlags &= 0xff00;
+        subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_POOR_PHOTOM;
+        subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_AREA;
+        subcatalog[0].measureT[Nmeasure].dbFlags &= ~ID_MEAS_NOCAL;
       }
       Nmeasure ++;
@@ -105 +151 @@
       if (Nmeasure == NMEASURE) {
         NMEASURE += 1000;
-        REALLOCATE (subcatalog[0].measure, Measure, NMEASURE);
+        REALLOCATE (subcatalog[0].measureT, MeasureTiny, NMEASURE);
       }
     }
@@ -122 +168 @@
       continue; 
     }
-    subcatalog[0].average[Naverage].Nmeasure = Nm;
+    subcatalog[0].averageT[Naverage].Nmeasure = Nm;
     Naverage ++;
     if (Naverage == NAVERAGE) {
       NAVERAGE += 50;
-      REALLOCATE (subcatalog[0].average, Average, NAVERAGE);
-      REALLOCATE (subcatalog[0].secfilt, SecFilt, NAVERAGE*PhotNsec);
-    }
-  }
-  REALLOCATE (subcatalog[0].average, Average, MAX (Naverage, 1));
-  REALLOCATE (subcatalog[0].measure, Measure, MAX (Nmeasure, 1));
-  REALLOCATE (subcatalog[0].secfilt, SecFilt, PhotNsec*MAX (Naverage, 1));
+      REALLOCATE (subcatalog[0].averageT, AverageTiny, NAVERAGE);
+      REALLOCATE (subcatalog[0].secfilt,  SecFilt,     NAVERAGE*Nsecfilt);
+    }
+  }
+  REALLOCATE (subcatalog[0].averageT, AverageTiny, MAX (Naverage, 1));
+  REALLOCATE (subcatalog[0].measureT, MeasureTiny, MAX (Nmeasure, 1));
+  REALLOCATE (subcatalog[0].secfilt, SecFilt, Nsecfilt*MAX (Naverage, 1));
   subcatalog[0].Naverage = Naverage;
   subcatalog[0].Nmeasure = Nmeasure;
   subcatalog[0].Nsecfilt = catalog[0].Nsecfilt;
   subcatalog[0].Nsecf_mem = Naverage * catalog[0].Nsecfilt;
-  assert (PhotNsec == catalog[0].Nsecfilt);
+
+  // limit the total number of stars in the catalog
+  if (MaxDensityUse) {
+    LimitDensityCatalog (subcatalog, catalog);
+  }
 
   if (VERBOSE) {
@@ -147 +197 @@
   return (TRUE);
 }
+
+int LimitDensityCatalog (Catalog *subcatalog, Catalog *catalog) {
+
+  Catalog tmpcatalog;
+
+  double Rmin, Rmax, Dmin, Dmax;
+
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+
+  gfits_scan (&catalog[0].header, "RA0",  "%lf", 1, &Rmin);
+  gfits_scan (&catalog[0].header, "DEC0", "%lf", 1, &Dmin);
+  gfits_scan (&catalog[0].header, "RA1",  "%lf", 1, &Rmax);
+  gfits_scan (&catalog[0].header, "DEC1", "%lf", 1, &Dmax);
+
+  if (VERBOSE2) fprintf (stderr, "extracting from catalog covering region %f,%f to %f,%f\n", Rmin, Dmin, Rmax, Dmax);
+
+  float AREA = fabs(Dmax - Dmin) * fabs(Rmax - Rmin) * cos (0.5*RAD_DEG*(Dmax + Dmin));
+  assert (AREA > 0);
+
+  off_t Nmax = MaxDensityValue * AREA;
+  if (subcatalog[0].Naverage <= Nmax) {
+    if (VERBOSE) {
+      fprintf (stderr, "subcatalog has less than the max density\n");
+    }
+    return (TRUE);
+  }
+
+  off_t Naverage = subcatalog[0].Naverage;
+
+  // select a random subset of Nmax stars from subcatalog using Fisher-Yates
+
+  // we are going to select Nmax entries by generating a random-sorted index list
+  off_t *index, tmp, i, j, ave;
+  ALLOCATE (index, off_t, Naverage);
+  for (i = 0; i < Naverage; i++) {
+    index[i] = i;
+  }
+  for (i = 0; i < Naverage; i++) {
+    j = (Naverage - i) * drand48() + i; // a number between i and Naverage
+    tmp = index[j];
+    index[j] = index[i];
+    index[i] = tmp;
+  }
+
+  // count the number of measurements this selection will yield
+  off_t NMEASURE = 0;
+  for (i = 0; i < Nmax; i++) {
+    ave = index[i];
+    NMEASURE += subcatalog[0].averageT[ave].Nmeasure;
+  }
+
+  // allocate the output data 
+  ALLOCATE (tmpcatalog.averageT, AverageTiny, Nmax);
+  ALLOCATE (tmpcatalog.measureT, MeasureTiny, NMEASURE);
+  ALLOCATE (tmpcatalog.secfilt,  SecFilt, Nmax * Nsecfilt);
+
+  off_t Nmeasure = 0;
+
+  // copy the Nmax selected entries from subcatalog to tmpcatalog (adjusting links)
+  for (i = 0; i < Nmax; i++) {
+    ave = index[i];
+    tmpcatalog.averageT[i] = subcatalog[0].averageT[ave];
+    tmpcatalog.averageT[i].measureOffset = Nmeasure;
+    for (j = 0; j < tmpcatalog.averageT[i].Nmeasure; j++) {
+      off_t offset = subcatalog[0].averageT[ave].measureOffset + j;
+      tmpcatalog.measureT[Nmeasure] = subcatalog[0].measureT[offset];
+      tmpcatalog.measureT[Nmeasure].averef = i;
+      Nmeasure ++;
+    }
+  }
+
+  if (VERBOSE) {
+    fprintf (stderr, "limited to "OFF_T_FMT" of "OFF_T_FMT" stars ("OFF_T_FMT" of "OFF_T_FMT" measures) for catalog %s\n", 
+             Nmax, subcatalog[0].Naverage, Nmeasure, subcatalog[0].Nmeasure,  catalog[0].filename);
+  }
+
+  free (subcatalog[0].averageT);
+  free (subcatalog[0].measureT);
+  free (subcatalog[0].secfilt);
+
+  subcatalog[0].averageT = tmpcatalog.averageT;
+  subcatalog[0].measureT = tmpcatalog.measureT;
+  subcatalog[0].secfilt = tmpcatalog.secfilt;
+  subcatalog[0].Naverage = Nmax;
+  subcatalog[0].Nmeasure = Nmeasure;
+  subcatalog[0].Nsecfilt = catalog[0].Nsecfilt;
+  subcatalog[0].Nsecf_mem = Naverage * catalog[0].Nsecfilt;
+
+  return (TRUE);
+}
+
+// for the cases where we are not using a subset of the data, we still need to have a copy of these fields
+int populate_tiny_values (Catalog *catalog) {
+
+  off_t i;
+
+  ALLOCATE (catalog[0].measureT, MeasureTiny, catalog[0].Nmeasure);
+  ALLOCATE (catalog[0].averageT, AverageTiny, catalog[0].Naverage);
+
+  for (i = 0; i < catalog[0].Naverage; i++) {
+    CopyAverageTiny (&catalog[0].averageT[i], &catalog[0].average[i]);
+  }
+
+  for (i = 0; i < catalog[0].Nmeasure; i++) {
+    CopyMeasureTiny (&catalog[0].measureT[i], &catalog[0].measure[i]);
+  }
+
+  return (TRUE);
+}
+
+int free_tiny_values (Catalog *catalog) {
+
+  free (catalog[0].averageT);
+  free (catalog[0].measureT);
+  return (TRUE);
+}
+

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

@@ -7 +7 @@
   initstats ("MEAN");
 
-  stN = statsStarN (catalog, Ncatalog);
-  stX = statsStarX (catalog, Ncatalog);
-  stS = statsStarS (catalog, Ncatalog);
+  fprintf (stderr, "\n");
+  fprintf (stderr, "STATS            median     mean    sigma      min      max   Nmeas\n");
 
+  int Ns;
+  for (Ns = 0; Ns < Nphotcodes; Ns++) {
+
+    int thisCode = photcodes[Ns][0].code;
+    int Nsec = GetPhotcodeNsec(thisCode);
+    int seccode = photcodes[Ns][0].code;
+
+    stN = statsStarN (catalog, Ncatalog, Nsec, seccode);
+    stX = statsStarX (catalog, Ncatalog, Nsec);
+    stS = statsStarS (catalog, Ncatalog, Nsec);
+  
+    fprintf (stderr, "   --- stats for %s ---\n", photcodes[Ns][0].name);
+    fprintf (stderr, "meas / star:    %7.0f  %7.1f  %7.1f  %7.0f  %7.0f  %6d\n",   stN.median, stN.mean, stN.sigma, stN.min, stN.max, stN.Nmeas);
+    fprintf (stderr, "dMrel star:     %7.4f  %7.4f  %7.4f  %7.4f  %7.4f  %6d\n", stS.median, stS.mean, stS.sigma, stS.min, stS.max, stS.Nmeas);
+    fprintf (stderr, "chisq star:     %7.1f  %7.1f  %7.1f  %7.1f  %7.1f  %6d\n",   stX.median, stX.mean, stX.sigma, stX.min, stX.max, stX.Nmeas);
+  }
+  
   imN = statsImageN (catalog);
   imX = statsImageX (catalog);
   imM = statsImageM (catalog);
   imD = statsImagedM (catalog);
-
+  
   msN = statsMosaicN (catalog);
   msM = statsMosaicM (catalog);
   msD = statsMosaicdM (catalog);
   msX = statsMosaicX (catalog);
-
-  fprintf (stderr, "STATS            median     mean    sigma      min      max   Nmeas\n");
+  
   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, "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, "meas / star:    %7.0f  %7.1f  %7.1f  %7.0f  %7.0f  %6d\n",   stN.median, stN.mean, stN.sigma, stN.min, stN.max, stN.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, "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);
-  fprintf (stderr, "chisq star:     %7.1f  %7.1f  %7.1f  %7.1f  %7.1f  %6d\n",   stX.median, stX.mean, stX.sigma, stX.min, stX.max, stX.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, "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, "dMrel star:     %7.4f  %7.4f  %7.4f  %7.4f  %7.4f  %6d\n\n", stS.median, stS.mean, stS.sigma, stS.min, stS.max, stS.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);
+  
 
   initstats (STATMODE);

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

@@ -13 +13 @@
     N = UserPatchSelect ? 1 : 2;
 
+# if (0)
+    // XXX DEP
     if ((photcode = GetPhotcodebyName (argv[N])) == NULL) {
       fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", argv[N]);
@@ -21 +23 @@
       exit (1);
     }
+    // PhotSec is used to select the single average photcode being processed
     PhotSec = GetPhotcodeNsec (photcode[0].code);
+# endif
+
+    Nphotcodes = 0;
+    photcodes = NULL;
+    int NPHOTCODES = 10;
+    ALLOCATE (photcodes, PhotCode *, NPHOTCODES);
+
+    /* parse the comma-separated list of photcodesKeep */
+    char *myList = strcreate(argv[N]);
+    char *list = myList;
+    char *codename = NULL;
+    char *ptr = NULL;
+    while ((codename = strtok_r (list, ",", &ptr)) != NULL) {
+      list = NULL; // pass NULL on successive strtok_r calls
+      fprintf (stderr, "PHOTCODE LIST: %s\n", myList);
+      fprintf (stderr, "codename: %s\n", codename);
+      if ((photcodes[Nphotcodes] = GetPhotcodebyName (codename)) == NULL) {
+        fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", codename);
+        exit (1);
+      }
+      if (photcodes[Nphotcodes][0].type != PHOT_SEC) {
+          fprintf (stderr, "photcode %s is not an filter type (SEC)\n", codename);
+          exit (1);
+      }
+      Nphotcodes ++;
+      CHECK_REALLOCATE (photcodes, PhotCode *, NPHOTCODES, Nphotcodes, 10);
+    }
   }
-
-  PhotNsec = GetPhotcodeNsecfilt ();
+  // XXX DEP PhotNsec = GetPhotcodeNsecfilt ();
+  if (USE_GRID && (Nphotcodes > 1)) {
+    fprintf (stderr, "grid correction analysis currently can only operate on a single photcode\n");
+    exit (1);
+  }
 
   initstats (STATMODE);
@@ -55 +88 @@
     exit (0);
   }
+
+  // init the random seed
+  long A, B;
+  A = time(NULL);
+  for (B = 0; A == time(NULL); B++);
+  srand48(B);
 }
-

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

@@ -44 +44 @@
   // select the images which overlap the selected sky regions
   subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset);
-  MARKTIME("selected images: %f sec\n", dtime);
+  MARKTIME("selected %d overlapping images: %f sec\n", (int) Nsubset, dtime);
 
   // generate db->vtable from db->ftable based on the selection

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

@@ -18 +18 @@
   ALLOCATE (ilist, double, Ntot);
 
-  N = 0;
-  for (i = 0; i < Ncatalog; i++) {
-    for (j = 0; j < catalog[i].Naverage; j++) {
+  int Nsecfilt = GetPhotcodeNsecfilt ();
 
-      /* calculate the average value for a single star */
-      if (catalog[i].average[j].flags & STAR_BAD) continue;  
-      m = catalog[i].average[j].measureOffset;
+  int Ns;
+  for (Ns = 0; Ns < Nphotcodes; Ns++) {
 
-      for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
-        if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue;
-        Mcal = getMcal  (m, i);
-        if (isnan(Mcal)) continue;
-        Mmos = getMmos  (m, i);
-        if (isnan(Mmos)) continue;
-        Mgrid = getMgrid (m, i);
-        if (isnan(Mgrid)) continue;
+    int thisCode = photcodes[Ns][0].code;
+    int Nsec = GetPhotcodeNsec(thisCode);
 
-        Mrel = catalog[i].secfilt[PhotNsec*j+PhotSec].M;
-        xlist[N] = Mrel;
-        ylist[N] = PhotSys  (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]) - Mcal - Mmos - Mgrid - Mrel;
-        ilist[N] = PhotInst (&catalog[i].measure[m]);
-        N++;
-      }
+    N = 0;
+    for (i = 0; i < Ncatalog; i++) {
+        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].average[j].measureOffset;
+
+            for (k = 0; k < catalog[i].average[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);
+                if (isnan(Mcal)) continue;
+                Mmos = getMmos  (m, i);
+                if (isnan(Mmos)) continue;
+                Mgrid = getMgrid (m, i);
+                if (isnan(Mgrid)) continue;
+
+                Mrel = catalog[i].secfilt[Nsecfilt*j+Nsec].M;
+                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]);
+                N++;
+            }
+        }
     }
+
+    plot_defaults (&graphdata);
+    graphdata.xmin = PlotMmin;
+    graphdata.xmax = PlotMmax;
+    graphdata.ymin = PlotdMmin;
+    graphdata.ymax = PlotdMmax;
+    plot_list (&graphdata, xlist, ylist, N, "mag vs dmag", "%s.Mag.png", OUTROOT);
+
+    plot_defaults (&graphdata);
+    graphdata.ymin = PlotdMmin;
+    graphdata.ymax = PlotdMmax;
+    plot_list (&graphdata, ilist, ylist, N, "imag vs dmag", "%s.iMag.png", OUTROOT);
   }
-
-  plot_defaults (&graphdata);
-  graphdata.xmin = PlotMmin;
-  graphdata.xmax = PlotMmax;
-  graphdata.ymin = PlotdMmin;
-  graphdata.ymax = PlotdMmax;
-  plot_list (&graphdata, xlist, ylist, N, "mag vs dmag", "%s.Mag.png", OUTROOT);
-
-  plot_defaults (&graphdata);
-  graphdata.ymin = PlotdMmin;
-  graphdata.ymax = PlotdMmax;
-  plot_list (&graphdata, ilist, ylist, N, "imag vs dmag", "%s.iMag.png", OUTROOT);
   free (xlist);
   free (ylist);

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

@@ -1 +1 @@
 # include "relphot.h"
+
+# define TIMESTAMP(TIME) \
+    gettimeofday (&stop, (void *) NULL);        \
+    dtime = DTIME (stop, start);                \
+    TIME += dtime;                              \
+    gettimeofday (&start, (void *) NULL);
 
 void reload_catalogs (SkyList *skylist) {
@@ -8 +14 @@
   Catalog catalog;
 
+  struct timeval start, stop;
+  double dtime = 0.0;
+  double time1 = 0.0;
+  double time2 = 0.0;
+  double time3 = 0.0;
+  double time4 = 0.0;
+  double time5 = 0.0;
+  double time6 = 0.0;
+  double time7 = 0.0;
+
   if (VERBOSE) fprintf (stderr, "re-loading catalog data\n");
 
   /* load data from each region file */
   for (i = 0; i < skylist[0].Nregions; i++) {
+    gettimeofday (&start, (void *) NULL);
     catalog.filename = skylist[0].filename[i];
 
@@ -20 +37 @@
       continue;
     }
+    TIMESTAMP(time1);
 
     catalog.catformat = dvo_catalog_catformat (CATFORMAT);    // set the default catformat from config data
@@ -36 +54 @@
         continue;
     }
+    TIMESTAMP(time2);
+
+    populate_tiny_values(&catalog);
 
     initImageBins  (&catalog, 1);
     initMosaicBins (&catalog, 1);
-    initGridBins   (&catalog, 1);
+    initGridBins   (&catalog, 1); 
+    TIMESTAMP(time3);
 
-    findImages (&catalog, 1);
-    findMosaics (&catalog, 1);
+    findImages (&catalog, 1);  // FX
+    findMosaics (&catalog, 1); // 
+    TIMESTAMP(time4);
 
     setMrelFinal (&catalog);
+    TIMESTAMP(time5);
+
     dvo_catalog_save (&catalog, VERBOSE); 
     dvo_catalog_unlock (&catalog);
+
+    free_tiny_values(&catalog);
     dvo_catalog_free (&catalog);
+    TIMESTAMP(time6);
 
     freeImageBins (1);
     freeMosaicBins (1);
     freeGridBins (1);
+    TIMESTAMP(time7);
   }
+
+  fprintf (stderr, "time1 %f : find catalog\n", time1);
+  fprintf (stderr, "time2 %f : load catalog\n", time2);
+  fprintf (stderr, "time3 %f : init imbins\n",  time3);
+  fprintf (stderr, "time4 %f : find images\n",  time4);
+  fprintf (stderr, "time5 %f : set Mrel\n",     time5);
+  fprintf (stderr, "time6 %f : save catalog\n", time6);
+  fprintf (stderr, "time7 %f : free catalog\n", time7);
 }

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

@@ -47 +47 @@
     if (!dvo_image_load (&db, VERBOSE, FALSE)) Shutdown ("can't read image catalog %s", db.filename);
   }
-  MARKTIME("load image data: %f sec\n", dtime);
+  MARKTIME("-- load image data: %f sec\n", dtime);
 
   /* load regions and images based on specified sky patch */
@@ -54 +54 @@
   // XXX this is fairly lame: argv[1] is photcode if UserPatchSelect is true
   skylist = load_images (&db, argv[1], &UserPatch, UserPatchSelect);
-  MARKTIME("load images: %f sec\n", dtime);
+  MARKTIME("-- load images: %f sec\n", dtime);
 
   /* unlock, if we can (else, unlocked below) */
@@ -61 +61 @@
   /* load catalog data from region files */
   catalog = load_catalogs (skylist, &Ncatalog);
-  MARKTIME("load catalog data: %f sec\n", dtime);
+  MARKTIME("-- load catalog data: %f sec\n", dtime);
   
   /* add in a loop over the catalogs calling dvo_catalog_chipcoords */
@@ -67 +67 @@
   /* match measurements with images, mosaics */
   initImageBins  (catalog, Ncatalog);
-  MARKTIME("make image bins: %f sec\n", dtime);
+  MARKTIME("-- make image bins: %f sec\n", dtime);
 
   initMosaicBins (catalog, Ncatalog);
@@ -74 +74 @@
 
   findImages (catalog, Ncatalog);
-  MARKTIME("set up image indexes: %f sec\n", dtime);
+  MARKTIME("-- set up image indexes: %f sec\n", dtime);
 
   findMosaics (catalog, Ncatalog);  /* also sets Grid values */
-  MARKTIME("set up mosaic indexes: %f sec\n", dtime);
+  MARKTIME("-- set up mosaic indexes: %f sec\n", dtime);
 
   SAVEPLOT = FALSE;
 
   setExclusions (catalog, Ncatalog);
+
+  global_stats (catalog, Ncatalog);
 
   if (PLOTSTUFF) {
@@ -91 +93 @@
   if (USE_GRID) {
       int star_toofew;
-
-# if (USE_DIRECT)
-      // until we finish the grid analysis, do not reject stars out-of-hand based on ID_STAR_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;
-      STAR_BAD  = ID_STAR_POOR;
-
-      showGridCount ();
-      setMgridDirect (catalog, Ncatalog);
-
-      STAR_BAD  = ID_STAR_POOR | ID_STAR_FEW;
-      STAR_TOOFEW = star_toofew;
-
-      dump_grid ();
-      exit (0);
-
-# else
 
       // until we finish the grid analysis, do not reject stars out-of-hand based on ID_STAR_FEW
@@ -122 +106 @@
       STAR_BAD  = ID_STAR_POOR | ID_STAR_FEW;
       STAR_TOOFEW = star_toofew;
-# endif
   }
 
@@ -140 +123 @@
       plot_chisq (catalog, Ncatalog);
     }
-    if ((i == 1) || (i == 5) || (i ==  9) || (i == 13)) clean_measures (catalog, Ncatalog, FALSE); 
-    if ((i == 2) || (i == 6) || (i == 10) || (i == 14)) clean_stars (catalog, Ncatalog);
-    if ((i == 4) || (i == 8) || (i == 12) || (i == 16)) clean_mosaics ();
-    if ((i == 4) || (i == 8) || (i == 12) || (i == 16)) clean_images ();
+    if (i % 6 == 1) rationalize_mosaics (catalog, Ncatalog);
+    // if (i % 6 == 1) rationalize_images ();
+    if (i % 6 == 2) clean_measures (catalog, Ncatalog, FALSE); 
+    if (i % 6 == 3) clean_stars (catalog, Ncatalog);
+    if (i % 6 == 5) clean_mosaics ();
+    if (i % 6 == 5) clean_images ();
+
+    // if ((i == 1) || (i == 5) || (i ==  9) || (i == 13)) clean_measures (catalog, Ncatalog, FALSE); 
+    // if ((i == 2) || (i == 6) || (i == 10) || (i == 14)) clean_stars (catalog, Ncatalog);
+    // if ((i == 4) || (i == 8) || (i == 12) || (i == 16)) clean_mosaics ();
+    // if ((i == 4) || (i == 8) || (i == 12) || (i == 16)) clean_images ();
     global_stats (catalog, Ncatalog);
+    MARKTIME("-- finished loop %d: %f sec\n", i, dtime);
   }
 
-  SAVEPLOT = TRUE;
-  plot_scatter (catalog, Ncatalog); 
-  plot_grid (catalog); 
-  plot_mosaics ();
-  plot_images ();
-  plot_stars (catalog, Ncatalog);
-  plot_chisq (catalog, Ncatalog);
-
+  if (PLOTSTUFF) {
+    plot_scatter (catalog, Ncatalog); 
+    plot_grid (catalog); 
+    plot_mosaics ();
+    plot_images ();
+    plot_stars (catalog, Ncatalog);
+    plot_chisq (catalog, Ncatalog);
+  }
+  
   if (USE_GRID) dump_grid ();
   if (!UPDATE) exit (0);
@@ -161 +153 @@
   setMcal  (catalog, TRUE);
   setMmos  (catalog, TRUE);
+  MARKTIME("-- finalize Mcal values: %f sec\n", dtime);
 
   /* at this point, we have correct cal coeffs in the image/mosaic structures */
-  for (i = 0; i < Ncatalog; i++) dvo_catalog_free (&catalog[i]);
+  for (i = 0; i < Ncatalog; i++) {
+    free_tiny_values (&catalog[i]);
+    dvo_catalog_free (&catalog[i]);
+  }
   freeImageBins (Ncatalog);
   freeMosaicBins (Ncatalog);
@@ -170 +166 @@
   /* load catalog data from region files, update Mrel include all data */
   reload_catalogs (skylist);
+  MARKTIME("-- updated all catalogs: %f sec\n", dtime);
+
   setMcalFinal ();
-
   reload_images (&db);
   

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

@@ -3 +3 @@
 int relphot_objects () {
 
-  off_t i, j, k, m;
+  off_t i, j, k;
   int Nsecfilt;
 
@@ -38 +38 @@
     }
 
-    // XXX consider what gets reset (only PHOTOM flags)
+    // reset 
     if (RESET) {
       Nsecfilt = catalog.Nsecfilt;
@@ -51 +51 @@
           catalog.secfilt[j*Nsecfilt + k].Ncode = 0;
           catalog.secfilt[j*Nsecfilt + k].Nused = 0;
-        }
-        m = catalog.average[j].measureOffset;
-        for (k = 0; k < catalog.average[j].Nmeasure; k++) {
-          catalog.measure[m+k].dbFlags = 0;
-          catalog.measure[m+k].Mcal = 0;
+          // XXX reset the photometry flags for secfilt entries?
         }
       }

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

@@ -64 +64 @@
   DmaxSkyRegion = -90.0;
 
+  // FILE *ftest = fopen ("relphot.dump.dat", "w");
+
   /* compare with each region file */
   for (i = 0; i < skylist[0].Nregions; i++) { 
@@ -117 +119 @@
     /* exclude images by photcode */
     ecode = GetPhotcodeEquivCodebyCode (timage[i].photcode);
-    if (ecode != photcode[0].code) continue;
+    found = FALSE;
+    int Ns;
+    for (Ns = 0; !found && (Ns < Nphotcodes); Ns++) {
+        if (ecode == photcodes[Ns][0].code) found = TRUE;
+    }
+    if (!found) continue;
 
     /* exclude images by time */
@@ -130 +137 @@
       continue;
     }
+
+    // XXX temporary test : record center coords for each accepted and rejected chip/mosaic
+    // double RAo, DECo;
 
     /* define image corners - note the DIS images (mosaic phu) are special */
@@ -138 +148 @@
       Xi[3] = -0.5*timage[i].NX; Yi[3] = +0.5*timage[i].NY;
       Xi[4] = -0.5*timage[i].NX; Yi[4] = -0.5*timage[i].NY;
+      // XY_to_RD(&RAo, &DECo, 0.0, 0.0, &timage[i].coords);
     } else {
       Xi[0] = 0;            Yi[0] = 0;
@@ -144 +155 @@
       Xi[3] = 0;            Yi[3] = timage[i].NY;
       Xi[4] = 0;            Yi[4] = 0;
+      // XY_to_RD(&RAo, &DECo, 0.5*timage[i].NX, 0.5*timage[i].NY, &timage[i].coords);
     }
     found = FALSE;
 
     /* transform corners to ra,dec -- costs ~3sec for 3M images (pikake) */
-    double RminImage = 360.0;
-    double RmaxImage =   0.0;
+    double RminImage = RmidSkyRegion + 180.0;
+    double RmaxImage = RmidSkyRegion - 180.0;
     double DminImage = +90.0;
     double DmaxImage = -90.0;
@@ -218 +230 @@
 
   found_it:
+    // XXX We claim this is a good image: write to a test file
+    // fprintf (ftest, "%s : %lf %lf  : %f %f %d %x\n", timage[i].name, RAo, DECo, timage[i].Mcal, timage[i].dMcal, timage[i].nstar, timage[i].flags);
+
     image[nimage] = timage[i]; 
     /* always allow 'few' images to succeed, if possible */
@@ -239 +254 @@
   }
   MARKTIME("finish image selection: %f sec\n", dtime);
+
+  // fclose (ftest);
 
   if (VERBOSE) fprintf (stderr, "found "OFF_T_FMT" images\n", nimage);

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

@@ -1 +1 @@
 # include "relphot.h"
+
+// this function sets the NOCAL and AREA dbFlags bits for the MeasureTiny elements these
+// are used elsewhere (StarOps.c, ImageOps.c, MosaicOps.c, GridOps.c, etc) to skip bad
+// measurements.  The only exception is 'setMave' which is called by 'relphot_objects',
+// and uses the bits read from disk as the test
 
 int setExclusions (Catalog *catalog, int Ncatalog) {
 
   off_t i, j, k, m, Narea, Nnocal, Ngood;
-  int ecode;
+  int ecode, found, Ns;
   Coords *coords;
   double r, d, x, y;
@@ -11 +16 @@
   for (i = 0; i < Ncatalog; i++) {
     for (j = 0; j < catalog[i].Naverage; j++) {
-      m = catalog[i].average[j].measureOffset;
-      for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
+      m = catalog[i].averageT[j].measureOffset;
+      for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) {
 
         /* select measurements by photcode */
-        ecode = GetPhotcodeEquivCodebyCode (catalog[i].measure[m].photcode);
-        if (ecode != photcode[0].code) goto mark_nocal;
+        ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode);
+        found = FALSE;
+        for (Ns = 0; !found && (Ns < Nphotcodes); Ns++) {
+            if (ecode == photcodes[Ns][0].code) found = TRUE;
+        }
+        if (!found) goto mark_nocal;
         
         /* select measurements by time */
         if (TimeSelect) {
-          if (catalog[i].measure[m].t < TSTART) goto mark_nocal;
-          if (catalog[i].measure[m].t > TSTOP) goto mark_nocal;
+          if (catalog[i].measureT[m].t < TSTART) goto mark_nocal;
+          if (catalog[i].measureT[m].t > TSTOP) goto mark_nocal;
         }
 
         /* select measurements by mag limit */
         if (AreaSelect) {
-          r = catalog[i].average[j].R + catalog[i].measure[m].dR / 3600.0;
-          d = catalog[i].average[j].D + catalog[i].measure[m].dD / 3600.0;
+          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;
           if ((coords = getCoords (m, i)) == NULL) goto markbad;
           RD_to_XY (&x, &y, r, d, coords);
@@ -39 +48 @@
 
       markbad:
-        catalog[i].measure[m].dbFlags |= ID_MEAS_AREA;
+        catalog[i].measureT[m].dbFlags |= ID_MEAS_AREA;
         Narea ++;
         continue;
         
       mark_nocal:
-        catalog[i].measure[m].dbFlags |= ID_MEAS_NOCAL;
+        catalog[i].measureT[m].dbFlags |= ID_MEAS_NOCAL;
         Nnocal ++;
         continue;

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

@@ -1 +1 @@
 # include "relphot.h"
+
+// we've just reloaded the data from disk; we now need to apply the Image/Mosaic/Grid
+// calibrations determined by the rest of the program.  We also need to set the final
+// output dbFlags values
 
 void setMrelFinal (Catalog *catalog) {
@@ -9 +13 @@
   if (RESET) {
 
-    for (i = 0; i < catalog[0].Naverage; i++) {
-      catalog[0].secfilt[PhotNsec*i+PhotSec].M  = NAN;
-      catalog[0].secfilt[PhotNsec*i+PhotSec].dM = NAN;
-      catalog[0].secfilt[PhotNsec*i+PhotSec].Xm = NAN_S_SHORT;
+    int Nsecfilt = GetPhotcodeNsecfilt ();
 
-      m = catalog[0].average[i].measureOffset;
-      for (j = 0; j < catalog[0].average[i].Nmeasure; j++, m++) {
+    int Ns;
+    for (Ns = 0; Ns < Nphotcodes; Ns++) {
         
-        /* select measurements by photcode */
-        ecode = GetPhotcodeEquivCodebyCode (catalog[0].measure[m].photcode);
-        if (ecode != photcode[0].code) continue;
+      int thisCode = photcodes[Ns][0].code;
+      int Nsec = GetPhotcodeNsec(thisCode);
+
+      for (i = 0; i < catalog[0].Naverage; i++) {
+        catalog[0].secfilt[Nsecfilt*i+Nsec].M  = NAN;
+        catalog[0].secfilt[Nsecfilt*i+Nsec].dM = NAN;
+        catalog[0].secfilt[Nsecfilt*i+Nsec].Xm = NAN_S_SHORT;
+
+        m = catalog[0].average[i].measureOffset;
+        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;
+          // skip measurements that do not match the current photcode
+          ecode = GetPhotcodeEquivCodebyCode (catalog[0].measure[m].photcode);
+          if (ecode != thisCode) { continue; }
+
+          /* select measurements by time */
+          if (TimeSelect) {
+            if (catalog[0].measure[m].t < TSTART) continue;
+            if (catalog[0].measure[m].t > TSTOP) continue;
+          }
+        
+          catalog[0].measure[m].Mcal = 0;
+          catalog[0].measure[m].dbFlags &= 0xff00;
+          catalog[0].measure[m].dbFlags &= ~ID_MEAS_POOR_PHOTOM;
+          catalog[0].measure[m].dbFlags &= ~ID_MEAS_SKIP_PHOTOM;
+          catalog[0].measure[m].dbFlags &= ~ID_MEAS_AREA;
+          catalog[0].measure[m].dbFlags &= ~ID_MEAS_NOCAL;
         }
-        
-        catalog[0].measure[m].Mcal = 0;
-        catalog[0].measure[m].dbFlags &= 0xff00;
-        catalog[0].measure[m].dbFlags &= ~ID_MEAS_POOR_PHOTOM;
-        catalog[0].measure[m].dbFlags &= ~ID_MEAS_SKIP_PHOTOM;
-        catalog[0].measure[m].dbFlags &= ~ID_MEAS_AREA;
-        catalog[0].measure[m].dbFlags &= ~ID_MEAS_NOCAL;
       }
     }
   }
 
+  // this sets flags in the measureT element, not the measure element
   setExclusions (catalog, 1);  /* mark by area */
 
@@ -53 +67 @@
   setMcalOutput (catalog, 1);
 
-  /* clear ID_STAR_POOR, ID_STAR_FEW, ID_MEAS_NOCAL values before writing ??? */
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+
+  /* clear ID_STAR_POOR, ID_STAR_FEW values before writing ??? */
+  /* ID_MEAS_NOCAL is an internal bit, so it should be cleared */
   for (i = 0; i < catalog[0].Naverage; i++) {
     catalog[0].average[i].flags &= ~ID_STAR_FEW;
     catalog[0].average[i].flags &= ~ID_STAR_POOR;
+    for (j = 0; j < Nsecfilt; j++) {
+        catalog[0].secfilt[i*Nsecfilt+j].flags &= ~ID_STAR_FEW; 
+        catalog[0].secfilt[i*Nsecfilt+j].flags &= ~ID_STAR_POOR; 
+    }
     m = catalog[0].average[i].measureOffset;
     for (j = 0; j < catalog[0].average[i].Nmeasure; j++, m++) {
@@ -68 +89 @@
 
   off_t i, k, m;
-  int ecode;
+  int ecode, found, Ns;
   off_t Ntot, Ntry, Nkeep, Nskip;
   float mag;
@@ -90 +111 @@
 
       /* clear SKIP for all measures at first */
-      catalog[0].measure[m].dbFlags &= ~ID_MEAS_SKIP_PHOTOM;
+      catalog[0].measureT[m].dbFlags &= ~ID_MEAS_SKIP_PHOTOM;
 
       /** never use these measurements (wrong photcode, bad time range) */
       /* skipped via NOCAL, don't mark as skipped */
       ecode = GetPhotcodeEquivCodebyCode (catalog[0].measure[m].photcode);
-      if (ecode != photcode[0].code) continue;
+      found = FALSE;
+      for (Ns = 0; !found && (Ns < Nphotcodes); Ns++) {
+        if (ecode == photcodes[Ns][0].code) found = TRUE;
+      }
+      if (!found) continue;
 
       /* skip measurements by time range */
@@ -134 +159 @@
 
     skip:
-      catalog[0].measure[m].dbFlags |= ID_MEAS_SKIP_PHOTOM;
+      catalog[0].measure [m].dbFlags |= ID_MEAS_SKIP_PHOTOM;
+      catalog[0].measureT[m].dbFlags |= ID_MEAS_SKIP_PHOTOM;
       Nskip ++;
     }