Changeset 31261

Timestamp:

Apr 10, 2011, 4:42:40 PM (15 years ago)

Author:

eugene

Message:

modifications to relphot to support multi-filter analysis (not yet working)

Location:

branches/eam_branches/ipp-20110404/Ohana/src/relphot

Files:

• include/relphot.h (modified) (1 diff)
• src/GridOps.c (modified) (10 diffs)
• src/ImageOps.c (modified) (5 diffs)
• src/MosaicOps.c (modified) (5 diffs)
• src/StarOps.c (modified) (19 diffs)
• src/bcatalog.c (modified) (7 diffs)
• src/initialize.c (modified) (2 diffs)
• src/plot_scatter.c (modified) (1 diff)
• src/select_images.c (modified) (1 diff)
• src/setExclusions.c (modified) (2 diffs)
• src/setMrelFinal.c (modified) (3 diffs)

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

@@ -85 +85 @@
 int    RELPHOT_GRID_BINNING;
 
-PhotCode      *photcode;
-int            PhotNsec;
-int            PhotSec;
+int      *photseclist;
+int      Nphotcodes;
+PhotCode **photcodes;
+// int            PhotSec;
+// int            PhotNsec;
 
 PhotCode      *refPhotcode;

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

@@ -410 +410 @@
         // 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]);
         if (isnan(Msys)) {
@@ -528 +527 @@
   if (!USE_GRID) return;
 
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+
   Nmax = Nlist[0];
   for (i = 0; i < Ngrid; i++) {
@@ -566 +567 @@
       
       n = catalog[c].measure[m].averef;
-      Msys = PhotSys (&catalog[c].measure[m], &catalog[c].average[n], &catalog[c].secfilt[n*PhotNsec]);
+      Msys = PhotSys (&catalog[c].measure[m], &catalog[c].average[n], &catalog[c].secfilt[n*Nsecfilt]);
       if (isnan(Msys)) {
         Nsys++;
@@ -612 +613 @@
 
   if (!USE_GRID) return;
+
+  int Nsecfilt = GetPhotcodeNsecfilt ();
 
   N = 0;
@@ -642 +645 @@
 
       n = catalog[c].measure[m].averef;
-      Msys = PhotSys (&catalog[c].measure[m], &catalog[c].average[n], &catalog[c].secfilt[n*PhotNsec]);
+      Msys = PhotSys (&catalog[c].measure[m], &catalog[c].average[n], &catalog[c].secfilt[n*Nsecfilt]);
 
       xlist[N] = Xmeas[c][m];
@@ -704 +707 @@
   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 +730 @@
   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 +745 @@
   /* 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 +759 @@
   /* 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 +779 @@
 
     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);

branches/eam_branches/ipp-20110404/Ohana/src/relphot/src/ImageOps.c

@@ -128 +128 @@
 }
 
-/* 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;
+  
+ 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;
+      found = FALSE;
+      for (Ns = 0; !found && (Ns < Nphotcodes); Ns++) {
+        if (ecode == photcodes[Ns][0].code) found = TRUE;
+      }
+      if (!found) continue; 
       matchImage (catalog, j, i);
     }
@@ -145 +149 @@
 int findCCD (off_t idx, off_t meas, int cat, Measure *measure) {
 
-  int ccdnum;
+  int ccdnum, found, Ns, ecode;
   double X, Y;
   char *pname, *filter, *p, base[256];
@@ -151 +155 @@
   /* 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);
@@ -307 +320 @@
 
   if (FREEZE_IMAGES) return;
+
+  int Nsecfilt = GetPhotcodeNsecfilt ();
 
   if (PoorImages) {
@@ -360 +375 @@
       
       n = catalog[c].measure[m].averef;
-      Msys = PhotSys (&catalog[c].measure[m], &catalog[c].average[n], &catalog[c].secfilt[n*PhotNsec]);
+      Msys = PhotSys (&catalog[c].measure[m], &catalog[c].average[n], &catalog[c].secfilt[n*Nsecfilt]);
       if (isnan(Msys)) {
         Nsys++;

branches/eam_branches/ipp-20110404/Ohana/src/relphot/src/MosaicOps.c

@@ -262 +262 @@
 int findMosaics (Catalog *catalog, int Ncatalog) {
   
-  int i, ecode;
+  int i, ecode, found, Ns;
   off_t j;
 
@@ -274 +274 @@
       }
       ecode = GetPhotcodeEquivCodebyCode (catalog[i].measure[j].photcode);
-      if (photcode[0].code != ecode) continue;
+      found = FALSE;
+      for (Ns = 0; !found && (Ns < Nphotcodes); Ns++) {
+        if (ecode == photcodes[Ns][0].code) found = TRUE;
+      }
+      if (!found) continue;
       matchMosaics (catalog, j, i);
     }
@@ -387 +391 @@
 
   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;
@@ -397 +401 @@
 
   image = getimages (&N);
+
+  Nsecfilt = GetPhotcodeNsecfilt ();
 
   if (PoorImages) {
@@ -448 +454 @@
       
       n = catalog[c].measure[m].averef;
-      Msys = PhotSys (&catalog[c].measure[m], &catalog[c].average[n], &catalog[c].secfilt[n*PhotNsec]);
+      Msys = PhotSys (&catalog[c].measure[m], &catalog[c].average[n], &catalog[c].secfilt[n*Nsecfilt]);
       if (isnan(Msys)) {
         Nsys++;

branches/eam_branches/ipp-20110404/Ohana/src/relphot/src/StarOps.c

@@ -29 +29 @@
 float getMrel (Catalog *catalog, off_t meas, int cat) {
 
+  int Nsec, Nsecfilt, code;
   int ave;
   float value;
 
   ave = catalog[cat].measure[meas].averef;
+  code = catalog[cat].measure[meas].photcode;
+
+  // XXX DEP : this should use a flag associated with secfilt
   if (catalog[cat].average[ave].flags & STAR_BAD) return (NAN);  
  
-  value = catalog[cat].secfilt[PhotNsec*ave+PhotSec].M;
+  Nsec = GetPhotcodeNsec(code);
+  Nsecfilt = GetPhotcodeNsecfilt ();
+
+  value = catalog[cat].secfilt[Nsecfilt*ave+Nsec].M;
   return (value);
 }
@@ -45 +52 @@
   float Msys, Mcal, Mmos, Mgrid;
   StatType stats;
-
+  int Nsec, Nsecfilt, ecode;
+
+  Nsecfilt = GetPhotcodeNsecfilt ();
   Nfew = Nsys = Nbad = Ncal = Nmos = Ngrid = 0;
 
@@ -51 +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;
+        Nsec = GetPhotcodeNsec(thisCode);
+
+        /* calculate the average value for a single star */
+        // XXX this flag should be set by secfilt, not average
+        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++) {
+          // skip measurements that do not match the current photcode
+          ecode = GetPhotcodeEquivCodebyCode (catalog[i].measure[m].photcode);
+          if (ecode != thisCode) { continue; }
+
+          if (catalog[i].measure[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 = PhotSys (&catalog[i].measure[m], &catalog[i].average[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].measure[m].dM, MIN_ERROR);
+
+          // tie down reference photometry if the -refcode (code) option is selected
+          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;
+            }
           }
-        }
-
-        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;
+        // XXX DEP this makes no sense: need a separate flag for each secfilt
+        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[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;
+      }
     }
   }
@@ -129 +151 @@
   double *list, *dlist;
   StatType stats;
+  int Nsec, Nsecfilt, ecode;
+
+  Nsecfilt = GetPhotcodeNsecfilt ();
 
   /* Nmeasure is now different, need to reallocate */
@@ -142 +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].average[j].measureOffset;
+        for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
+          // skip measurements that do not match the current photcode
+          ecode = GetPhotcodeEquivCodebyCode (catalog[i].measure[m].photcode);
+          if (ecode != thisCode) { continue; }
+
+          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*Nsecfilt]);
+          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[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;
+      }
     }
   }
@@ -209 +244 @@
   ALLOCATE (dlist, double, MAX (1, Nmax));
 
-  Nsecfilt = catalog[0].Nsecfilt;
+  Nsecfilt = GetPhotcodeNsecfilt ();
 
 # define PSFQUALSTATS 1
@@ -227 +262 @@
           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;
 
@@ -346 +381 @@
 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;
@@ -361 +396 @@
   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].average[j].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].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);
+  }
   free (xlist);
   free (slist);
@@ -416 +462 @@
   int Ncal, Nmos, Ngrid, Nfew;
 
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+
   if (VERBOSE) fprintf (stderr, "marking poor measures\n");
   /* Nmeasure is now different, need to reallocate */
@@ -437 +485 @@
     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;
+
+        /* skip bad stars to prevent them from becoming good (on inner sample) */
+        if (catalog[i].average[j].flags & STAR_BAD) 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++) {
+          // skip measurements that do not match the current photcode
+          int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measure[m].photcode);
+          if (ecode != thisCode) { continue; }
+
+          /* 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*Nsecfilt]);
+          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++) {
+          // skip measurements that do not match the current photcode
+          int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measure[m].photcode);
+          if (ecode != thisCode) { continue; }
+
+          /* 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*Nsecfilt]);
+          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;
+      }
     }
   }
@@ -567 +629 @@
 }
 
+// XXXX do this for all secfilt at once or one at a time?
 StatType statsStarX (Catalog *catalog, int Ncatalog) {
 
@@ -574 +637 @@
   StatType stats;
 
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+
   Ntot = 0;
   for (i = 0; i < Ncatalog; i++) {
@@ -586 +651 @@
     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;  
-
-      Xm = catalog[i].secfilt[PhotNsec*j+PhotSec].Xm;
-      if (Xm == NAN_S_SHORT) continue;
-      list[n] = pow (10.0, 0.01*Xm);
-      dlist[n] = 1;
-      n++;
+      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 */
+        if (catalog[i].average[j].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);
+        dlist[n] = 1;
+        n++;
+      }
     }
   }
@@ -611 +683 @@
   StatType stats;
 
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+
   Ntot = 0;
   for (i = 0; i < Ncatalog; i++) {
@@ -623 +697 @@
     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;  
-
-      dM = catalog[i].secfilt[PhotNsec*j+PhotSec].dM;
-      list[n] = dM;
-      dlist[n] = 1;
-      n++;
+      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 */
+        if (catalog[i].average[j].flags & STAR_BAD) continue;  
+
+        dM = catalog[i].secfilt[Nsecfilt*j+Nsec].dM;
+        list[n] = dM;
+        dlist[n] = 1;
+        n++;
+      }
     }
   }
@@ -651 +732 @@
   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].average[j].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);
@@ -678 +767 @@
   Graphdata graphdata;
 
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+
   Ntotal = 0;
   for (i = 0; i < Ncatalog; i++) Ntotal += catalog[i].Naverage;
@@ -684 +775 @@
   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].average[j].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);

branches/eam_branches/ipp-20110404/Ohana/src/relphot/src/bcatalog.c

@@ -6 +6 @@
   
   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] */
@@ -21 +19 @@
   NMEASURE = 1000;
   ALLOCATE (subcatalog[0].average, Average, NAVERAGE);
-  ALLOCATE (subcatalog[0].secfilt, SecFilt, NAVERAGE*PhotNsec);
+  ALLOCATE (subcatalog[0].secfilt, SecFilt, NAVERAGE*Nsecfilt);
   ALLOCATE (subcatalog[0].measure, Measure, NMEASURE);
   Nmeasure = Naverage = 0;
@@ -34 +32 @@
     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];
+    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].average[Naverage].flags &= ~ID_STAR_FEW;
+        subcatalog[0].average[Naverage].flags &= ~ID_STAR_POOR;
+      }
     }
 
@@ -54 +59 @@
       /* 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 */
@@ -129 +141 @@
       NAVERAGE += 50;
       REALLOCATE (subcatalog[0].average, Average, NAVERAGE);
-      REALLOCATE (subcatalog[0].secfilt, SecFilt, NAVERAGE*PhotNsec);
+      REALLOCATE (subcatalog[0].secfilt, SecFilt, NAVERAGE*Nsecfilt);
     }
   }
   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].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
@@ -161 +172 @@
   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);
@@ -205 +218 @@
   // allocate the output data 
   ALLOCATE (tmpcatalog.average, Average, Nmax);
-  ALLOCATE (tmpcatalog.secfilt, SecFilt, Nmax * PhotNsec);
+  ALLOCATE (tmpcatalog.secfilt, SecFilt, Nmax * Nsecfilt);
   ALLOCATE (tmpcatalog.measure, Measure, NMEASURE);
 

branches/eam_branches/ipp-20110404/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);

branches/eam_branches/ipp-20110404/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].average[j].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].measure[m].photcode);
+                if (ecode != thisCode) { continue; }
+
+                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;
+
+                Mrel = catalog[i].secfilt[Nsecfilt*j+Nsec].M;
+                xlist[N] = Mrel;
+                ylist[N] = PhotSys  (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]) - Mcal - Mmos - Mgrid - Mrel;
+                ilist[N] = PhotInst (&catalog[i].measure[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);

branches/eam_branches/ipp-20110404/Ohana/src/relphot/src/select_images.c

@@ -119 +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 */

branches/eam_branches/ipp-20110404/Ohana/src/relphot/src/setExclusions.c

@@ -4 +4 @@
 
   off_t i, j, k, m, Narea, Nnocal, Ngood;
-  int ecode;
+  int ecode, found, Ns;
   Coords *coords;
   double r, d, x, y;
@@ -16 +16 @@
         /* select measurements by photcode */
         ecode = GetPhotcodeEquivCodebyCode (catalog[i].measure[m].photcode);
-        if (ecode != photcode[0].code) goto mark_nocal;
+        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 */

branches/eam_branches/ipp-20110404/Ohana/src/relphot/src/setMrelFinal.c

@@ -9 +9 @@
   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;
       }
     }
@@ -68 +77 @@
 
   off_t i, k, m;
-  int ecode;
+  int ecode, found, Ns;
   off_t Ntot, Ntry, Nkeep, Nskip;
   float mag;
@@ -95 +104 @@
       /* 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 */