Index: trunk/Ohana/src/opihi/dvo/photometry.c
===================================================================
--- trunk/Ohana/src/opihi/dvo/photometry.c	(revision 4420)
+++ trunk/Ohana/src/opihi/dvo/photometry.c	(revision 4585)
@@ -3,25 +3,28 @@
 /* match code to measure  */
 # define TESTCODE(C,M) \
-  switch (C[0].type) { \
-  case PHOT_DEP: \
-  case PHOT_REF: \
-    if (C[0].code != M.source) continue; \
-    break; \
-  case PHOT_PRI: \
-  case PHOT_SEC: \
-    if (C[0].code != GetPhotcodeEquivCodebyCode (M.source)) continue; \
-    break; \
-  default: \
-    break; \
-} 
+  if ((C != NULL) && ApplySelections[SelectionParam]) { \
+    switch (C[0].type) { \
+    case PHOT_DEP: \
+    case PHOT_REF: \
+      if (C[0].code != M.source) continue; \
+      break; \
+    case PHOT_PRI: \
+    case PHOT_SEC: \
+      if (C[0].code != GetPhotcodeEquivCodebyCode (M.source)) continue; \
+      break; \
+    default: \
+      break; \
+  } }
 
 /* exclusions based on measure.params  */
 # define TESTMEASURE(M) \
-  if (TimeSelect && (M.t < tzero)) continue; \
-  if (TimeSelect && (M.t > tend)) continue; \
-  if (ErrSelect  && (M.dM > ErrValue)) continue; \
-  if (TypeSelect && (TypeValue != GetMeasureTypeCode (&M))) continue; \
-  if (iMagSelect && (PhotInst (&M) < iMagMin)) continue; \
-  if (FlagSelect && (M.flags != FlagValue)) continue;
+  if (ApplySelections[SelectionParam]) { \
+    if (TimeSelect && (M.t < tzero)) continue; \
+    if (TimeSelect && (M.t > tend)) continue; \
+    if (ErrSelect  && (M.dM > ErrValue)) continue; \
+    if (TypeSelect && (TypeValue != GetMeasureTypeCode (&M))) continue; \
+    if (iMagSelect && (PhotInst (&M) < iMagMin)) continue; \
+    if (FlagSelect && (M.flags != FlagValue)) continue; \
+  }
 
 # define SETMAG(MOUT,MEAS,MODE) \
@@ -34,6 +37,8 @@
   if (MODE == MAG_AVE)  MOUT = PhotRel  (&MEAS, average, secfilt); \
   if (MODE == MAG_REF)  MOUT = PhotCal  (&MEAS, average, secfilt, measure, GetPhotcodeEquivbyCode (MEAS.source)); \
-  if (MagSelect && (MOUT > MagMax)) continue; \
-  if (MagSelect && (MOUT < MagMin)) continue;
+  if (ApplySelections[SelectionParam]) { \
+    if (MagSelect && (MOUT > MagMax)) continue; \
+    if (MagSelect && (MOUT < MagMin)) continue; \
+  }
 
 /* selection criteria */
@@ -51,4 +56,8 @@
 static double TypefracValue;
 
+/* apply selections or not */
+static int ApplySelections[4];
+static int SelectionParam;
+
 /* applied to Average quantities */
 static int PhotcodeSelect;
@@ -69,10 +78,9 @@
 static int TimeFormat;
 
-/* db image table */
-static Image *image = NULL;
-static int *subset = NULL;
-static int Nimage = 0;
-static int Nsubset = 0;
-static Coords mosaic;
+int GetTimeSelection (time_t *tz, time_t *te) {
+  *tz = tzero;
+  *te = tend;
+  return (TimeSelect);
+}
 
 int GetPhotcodeInfo (char *string, PhotCode **Code, int *Mode) {
@@ -156,4 +164,13 @@
 }
  
+int SetSelectionParam (int param) {
+  SelectionParam = param;
+  return (TRUE);
+}
+
+int GetSelectionParam () {
+  return (SelectionParam);
+}
+
 int GetMagMode (char *string) {
 
@@ -196,33 +213,10 @@
     fprintf (stderr, "value may be one of the following:\n");
     fprintf (stderr, " ra dR dec dD mag dmag Mrel Mcal photcode time fwhm dophot xccd yccd xmosaic ymosaic flags\n");
+    fprintf (stderr, "value may also be a valid photcode\n");
+    fprintf (stderr, "photcodes or 'mag' may have optional magnitude type: mag,[Minst, Mcat, Msys, Mrel, Mcal]\n");
   }
   return (param);
 }
   
-int TestPhotSelections (PhotCode **code) {
-
-  int NeedPhotcode;
-
-  /* if i've supplied a photcode (code != NULL), i'm not allowed to restrict it */
-  if ((code[0] != NULL) && PhotcodeSelect) return (FALSE);
-
-  /* if I have an average or ensemble restriction, I must have a photcode */
-  if (ChiSelect || NphotSelect || NcodeSelect || FWHMSelect) {
-    if (code[0] != NULL) return (TRUE);
-    if (PhotcodeSelect) {
-      code[0] = PhotcodeValue;
-      return (TRUE);
-    }
-    return (FALSE);
-  }
-
-  /* if I have a photcode restriction, apply it to the code */
-  if (PhotcodeSelect) {
-    code[0] = PhotcodeValue;
-    return (TRUE);
-  }
-  return (TRUE);
-}
-
 int GetAverageParam (char *parname) {
 
@@ -233,4 +227,5 @@
   if (!strcasecmp (parname, "dec"))   param = AVE_DEC;
   if (!strcasecmp (parname, "dmag"))  param = AVE_dMAG;
+  if (!strcasecmp (parname, "mag"))   param = AVE_MAG;
   if (!strcasecmp (parname, "Nmeas")) param = AVE_NMEAS;
   if (!strcasecmp (parname, "Nmiss")) param = AVE_NMISS;
@@ -251,4 +246,45 @@
 }
 
+/* I've set some selections - if these require a photcode, check if I set one */
+int TestPhotSelections (PhotCode **code, int param) {
+
+  int NeedPhotcode;
+
+  /* if i've supplied a photcode (code != NULL), i'm not allowed to restrict it */
+  if (code[0] != NULL) {
+    if (PhotcodeSelect) {
+      fprintf (stderr, "photcode selection rules violated: cannot restrict photcode with a photcode\n");
+      return (FALSE);
+    } else {
+      return (TRUE);
+    }
+  }
+
+  /* if I have an average or ensemble restriction, I must have a photcode */
+  NeedPhotcode = FALSE;
+  NeedPhotcode |= ChiSelect;
+  NeedPhotcode |= NphotSelect;
+  NeedPhotcode |= NcodeSelect;
+  NeedPhotcode |= ErrSelect;
+  NeedPhotcode |= TypeSelect;
+  NeedPhotcode |= TypefracSelect;
+  
+  /* for measure tests, supply MEAS_ZERO */
+  NeedPhotcode |= (param == AVE_Xm);
+  NeedPhotcode |= (param == AVE_MAG);
+  NeedPhotcode |= (param == AVE_dMAG);
+  NeedPhotcode |= (param == AVE_TYPE);
+  NeedPhotcode |= (param == AVE_NPHOT);
+  NeedPhotcode |= (param == AVE_NCODE);
+
+  if (NeedPhotcode || PhotcodeSelect) {
+    if (!PhotcodeSelect) return (FALSE);
+    code[0] = PhotcodeValue;
+    if (code[0][0].type == PHOT_DEP) return (FALSE);
+    if (code[0][0].type == PHOT_REF) return (FALSE);
+  }
+  return (TRUE);
+}
+
 void GetAverageParamHelp () {
   fprintf (stderr, "value may be one of the following:\n");
@@ -259,44 +295,4 @@
 }
 
-/* load images based on parameters and region, etc */
-int SetImageSelection (int mode, Graphdata *graphsky, int RegionSelect) {
-
-  image = NULL;
-  subset = NULL;
-  
-  switch (mode) {
-    case MEAS_XCCD:
-    case MEAS_YCCD:
-      break;
-    case MEAS_XMOSAIC: 
-    case MEAS_YMOSAIC:
-      /* mosaic defines a frame with 0,0 at the mosaic center, and 1 arcsec / pixel */
-      mosaic.crpix1 = mosaic.crpix2 = 0.0;
-      mosaic.cdelt1 = mosaic.cdelt2 = 1.0 / 3600;
-      mosaic.pc1_1  = mosaic.pc2_2  = 1.0;
-      mosaic.pc1_2  = mosaic.pc2_1  = 0.0;
-      mosaic.Npolyterms = 0;
-      strcpy (mosaic.ctype, "RA---SIN");
-      break;
-    default:
-      return (TRUE);
-  }
-
-  if ((image = LoadImages (&Nimage)) == NULL) return (FALSE);
-  BuildChipMatch (image, Nimage);
-  image_subset (image, Nimage, &subset, &Nsubset, graphsky, RegionSelect, tzero, (double) tend - tzero, TimeSelect);
-  sort_image_subset (image, subset, Nsubset);
-  return (TRUE);
-}
-
-/* free loaded images */
-void FreeImageSelection () {
-  if (image != NULL) free (image);
-  if (subset != NULL) free (subset);
-  image = NULL;
-  subset = NULL;
-  return;
-}
-
 /* (re)load photcodes from photcode table */
 int InitPhotcodes () {
@@ -318,7 +314,7 @@
 /* remove standard photometry filtering options, set selections */
 /* not all functions respect all selections... */
-int SetPhotSelections (int *argc, char **argv) {
-
-  int N;
+int SetPhotSelections (int *argc, char **argv, int Nparams) {
+
+  int i, N;
   double trange;
 
@@ -367,10 +363,19 @@
   }
 
-  /* select on value of flag (MEASURE ONLY) */
-  ChiSelect = FALSE;
-  if (N = get_argument (*argc, argv, "-chisq")) {
-    ChiSelect = TRUE;
-    remove_argument (N, argc, argv);
-    ChiLimit = atof (argv[N]);
+  /* select on value of Chisq (AVERAGE ONLY) */
+  for (i = 0; i < 4; i++) ApplySelections[i] = TRUE;
+  if (N = get_argument (*argc, argv, "-apply")) {
+    remove_argument (N, argc, argv);
+    if (strlen(argv[N]) != Nparams) {
+      fprintf (stderr, "-apply selection must define all parameter choices\n");
+      return (FALSE);
+    }
+    for (i = 0; i < Nparams; i++) {
+      if (toupper(argv[N][i]) == 'Y') {
+	ApplySelections[i] = TRUE;
+      } else {
+	ApplySelections[i] = FALSE;
+      }
+    }
     remove_argument (N, argc, argv);
   }
@@ -384,5 +389,5 @@
     if (PhotcodeValue == NULL) {
       fprintf (stderr, "photcode not found in photcode table\n");
-      PhotcodeSelect = FALSE;
+      return (FALSE);;
     }
     remove_argument (N, argc, argv);
@@ -422,4 +427,13 @@
   }
 
+  /* select on value of Chisq (AVERAGE ONLY) */
+  ChiSelect = FALSE;
+  if (N = get_argument (*argc, argv, "-chisq")) {
+    ChiSelect = TRUE;
+    remove_argument (N, argc, argv);
+    ChiLimit = atof (argv[N]);
+    remove_argument (N, argc, argv);
+  }
+
   /* select on measurement type: 1,2,3 (AVERAGE ONLY) */
   TypeSelect = FALSE;
@@ -508,104 +522,15 @@
   ALLOCATE (list, double, NLIST);
 
-  /* chisq, fwhm, Nphot, Ncode */
-  /* if (!TestAverage (code, average, secfilt, measure)) return (list); */
+  /* check selections based on averages & ensembles: chisq, Nphot, etc */
+  if (!TestAverage (code, average, secfilt, measure)) return (list); 
 
   /* look for measures */
   for (i = 0; i < average[0].Nm; i++) {
-    if (code != NULL) {
-      TESTCODE (code, measure[i]);
-    }
-    /* exclusions based on measure.params  */
-    TESTMEASURE (measure[i]);
-
-    /* need to have the mag value to test MagSelect */ 
-    SETMAG (M, measure[i], mode);
+    TESTCODE (code, measure[i]);  /* skip measurements not matching photcode */
+    TESTMEASURE (measure[i]);     /* exclusions based on measure.params  */
+    SETMAG (M, measure[i], mode); /* set appropriate magnitude (also does MagSelect) */ 
 
     /* assign value */
-    value = 0;
-    switch (param) {
-      case MEAS_MAG: /* magnitudes are already determined above */
-	value = M;
-	break;
-      case MEAS_RA: /* OK */
-	value = average[0].R - measure[i].dR / 360000.0;
-	break;
-      case MEAS_DEC: /* OK */
-	value = average[0].D - measure[i].dD / 360000.0;
-	break;
-      case MEAS_dMAG: /* OK */
-	value = 0.001*measure[i].dM;
-	break;
-      case MEAS_AIRMASS: /* OK */
-	value = 0.001*measure[i].airmass;
-	break;
-      case MEAS_EXPTIME: /* OK */
-	value = pow (10.0, measure[i].dt * 0.0004);
-	break;
-      case MEAS_PHOTCODE: /* OK */
-	value = measure[i].source;
-	break;
-      case MEAS_TIME: /* OK */
-	value = TimeValue (measure[i].t, TimeReference, TimeFormat);
-	break;
-      case MEAS_dR: /* OK */
-	value = 0.01*measure[i].dR;
-	break;
-      case MEAS_dD: /* OK */
-	value = 0.01*measure[i].dD;
-	break;
-      case MEAS_FWHM: /* OK */
-	value = 0.01*measure[i].FWx;
-	break;
-      case MEAS_DOPHOT: /* OK */
-	value = measure[i].dophot;
-	break;
-      case MEAS_FLAGS: /* ? */
-	value = measure[i].flags;
-	break;
-      case MEAS_XCCD: /* OK */
-	ra  = average[0].R - measure[i].dR / 360000.0;
-	dec = average[0].D - measure[i].dD / 360000.0;
-	m = match_image_subset (image, subset, Nsubset, measure[i].t, measure[i].source);
-	if (m == -1) break;
-	if (!FindMosaicForImage (image, Nimage, m)) break;
-	RD_to_XY (&x, &y, ra, dec, &image[m].coords);
-	value = x;
-	break;
-      case MEAS_YCCD: /* OK */
-	ra  = average[0].R - measure[i].dR / 360000.0;
-	dec = average[0].D - measure[i].dD / 360000.0;
-	m = match_image_subset (image, subset, Nsubset, measure[i].t, measure[i].source);
-	if (m == -1) break;
-	if (!FindMosaicForImage (image, Nimage, m)) break;
-	RD_to_XY (&x, &y, ra, dec, &image[m].coords);
-	value = y;
-	break;
-      case MEAS_XMOSAIC: /* OK */
-	ra  = average[0].R - measure[i].dR / 360000.0;
-	dec = average[0].D - measure[i].dD / 360000.0;
-	m = match_image_subset (image, subset, Nsubset, measure[i].t, measure[i].source);
-	if (m == -1) break;
-	/* this should use the mosaic assoicated with the image */
-	/* the solution here depends on image.coords having crref1,2 set to match the boresite center */
-	/* XXX set the mosaic.coords to have a unity WRP transformation in which x,y = ra,dec */
-	mosaic.crval1 = image[m].coords.crval1;
-	mosaic.crval2 = image[m].coords.crval2;
-	RD_to_XY (&x, &y, ra, dec, &mosaic);
-	value = x;
-	break;
-      case MEAS_YMOSAIC: /* OK */
-	ra  = average[0].R - measure[i].dR / 360000.0;
-	dec = average[0].D - measure[i].dD / 360000.0;
-	m = match_image_subset (image, subset, Nsubset, measure[i].t, measure[i].source);
-	if (m == -1) break;
-	/* same comments as above */
-	mosaic.crval1 = image[m].coords.crval1;
-	mosaic.crval2 = image[m].coords.crval2;
-	RD_to_XY (&x, &y, ra, dec, &mosaic);
-	value = y;
-	break;
-    }
-    list[Nlist] = value;
+    list[Nlist] = GetMeasure (param, &average[0], &measure[i], M);
     Nlist ++;
   }
@@ -622,19 +547,8 @@
   value = NO_MAG;
 
-  /* need valid code if param is: */
-  switch (param) {
-    case AVE_MAG:
-    case AVE_dMAG:
-    case AVE_Xm:
-    case AVE_TYPE:
-    case AVE_NCODE:
-    case AVE_NPHOT: 
-      if (code == NULL) return (NO_MAG);
-      break;
-    default:
-      break;
-  }
-
-  /* chisq, fwhm, Nphot, Ncode */
+  /* this function requires code set for certain value of param.  
+     use TestPhotSelectionsAverage to validate code/param choices */
+
+  /* filter by average quantities (eg, chisq, Nphot, etc) */
   if (!TestAverage (code, average, secfilt, measure)) return (NO_MAG);
 
@@ -728,5 +642,5 @@
 }
 
-/* determine the representative dophot type for this photcode */
+/* determine the representative dophot type for this photcode (must be PRI/SEC) */
 int DetermineTypeCode (Average *average, Measure *measure, int code) {
 
@@ -785,25 +699,13 @@
 }
 
+/* test if this average object meets the specified selection criteria.
+   for photcode-dependent quantities, only test for PRI/SEC photcodes */
 int TestAverage (PhotCode *code, Average *average, SecFilt *secfilt, Measure *measure) {
 
-  int i, Nm, Type;
+  int i, Nm, Type, Select;
   double fwhm, typefrac, dM, Xm;
 
-  /* exclusions based on average.params  */
-  if (ChiSelect) {
-    /* some average.params need a valid photcode for reference */
-    if (code == NULL) return (FALSE);
-    Xm = PhotXm (code, average, secfilt);
-    if (Xm == -1) return (FALSE);
-    if (Xm > ChiLimit) return (FALSE);
-  }
-  
-  /* for ErrSelect, check average errors */
-  if (ErrSelect) {
-    if (code == NULL) return (FALSE);  /* need a valid code */
-    dM = iPhotdM (code, average, secfilt);
-    if (dM > ErrValue) return (NO_MAG);
-  }
-  
+  if (!ApplySelections[SelectionParam]) return (TRUE);
+
   /* pass objects with more than FWHMfrac points with FWHM above / below FWHMvalue */ 
   if (FWHMSelect) {
@@ -827,4 +729,28 @@
   }
 
+  /* all selections below require a valid photcode */
+  Select = ChiSelect || ErrSelect || NcodeSelect || NphotSelect || TypeSelect || TypefracSelect;
+  if (!Select) return (TRUE);
+
+  /* must have a valid code of some kind */
+  if (code == NULL) return (FALSE);
+
+  /* only PRI/SEC photcodes apply the filter */
+  if (code[0].type == PHOT_DEP) return (TRUE);
+  if (code[0].type == PHOT_REF) return (TRUE);
+
+  /* exclusions based on average.params  */
+  if (ChiSelect) {
+    Xm = PhotXm (code, average, secfilt);
+    if (Xm == -1) return (FALSE);
+    if (Xm > ChiLimit) return (FALSE);
+  }
+  
+  /* for ErrSelect, check average errors */
+  if (ErrSelect) {
+    dM = iPhotdM (code, average, secfilt);
+    if (dM > ErrValue) return (NO_MAG);
+  }
+  
   /* for NcodeSelect, count Nmeas for appropriate photcode */
   if (NcodeSelect) {
@@ -937,10 +863,11 @@
 }
 
+/* extract delta-mag pairs applying specified selections */
 double *ExtractDMag (PhotCode **code, int *mode, Average *average, SecFilt *secfilt, Measure *measure, int *nlist) {
 
-  int i, j, A1, A2, N1, N2, Nlist, NLIST;
-  double *M1, *M2, *list;
-
-  /* check for special case of measure-measure */
+  int i, j, A1, A2, N1, N2, Np, Nlist, NLIST;
+  double *M1, *M2, *list, NoMag;
+
+  /* check for special case of measure-measure - this is needed to drop self-matches */
   A1 = ((mode[0] == MAG_AVE) || (mode[0] == MAG_REF));
   A2 = ((mode[1] == MAG_AVE) || (mode[1] == MAG_REF));
@@ -955,4 +882,5 @@
   ALLOCATE (list, double, NLIST);
   M1 = M2 = NULL;
+  NoMag = NO_MAG * 0.001;
 
   /* one of the two is an average, must do independently */
@@ -960,29 +888,34 @@
   if (N1 == 0) goto skip;
   
+  Np = GetSelectionParam ();
+  SetSelectionParam (Np + 1);
   M2 = ExtractMagnitudes (code[1], mode[1], average, secfilt, measure, &N2);
   if (N2 == 0) goto skip;
 
-  /* average mags can return NO_MAG : skip them? */
+  /* magnitudes may be NO_MAG : set delta to NO_MAG */
   for (i = 0; i < N1; i++) {
     for (j = 0; j < N2; j++) {
-      list[Nlist] = M1[i] - M2[j];
+      if ((M1[i] == NoMag) || (M2[j] == NoMag)) {
+	list[Nlist] = NoMag;
+      } else {
+	list[Nlist] = M1[i] - M2[j];
+      }
       Nlist ++;
     }
   }
-  *nlist = Nlist;
-  return (list);
 
  skip: 
   if (M1 != NULL) free (M1);
   if (M2 != NULL) free (M2);
-  M1 = M2 = NULL;
+  *nlist = Nlist;
   return (list);
 }
   
-/* extract a list of delta-measure-mags from the specified average entry based on the pre-set selections */
+/* extract a list of delta-measure-mags from the specified average entry based on the 
+   pre-set selections - does not return self-matched measurements */
 double *ExtractMeasuresDMag (PhotCode **code, int *mode, Average *average, SecFilt *secfilt, Measure *measure, int *nlist) {
 
-  int i, j, Nlist, NLIST;
-  double *list, M1, M2, value;
+  int i, j, Np0, Np1, Nlist, NLIST;
+  double *list, M1, M2, value, NoMag;
   
   *nlist = 0; 
@@ -990,22 +923,21 @@
   NLIST = MAX (1, average[0].Nm*average[0].Nm);
   ALLOCATE (list, double, NLIST);
-
-  /* must have two code values */
+  NoMag = NO_MAG * 0.001;
+
+  /* must have two code values - drop this test? this is programming case, not a user case */
   if (code == NULL) return (list);
   if (code[0] == NULL) return (list);
   if (code[1] == NULL) return (list);
 
-  /* exclude based on chisq for both codes  */
-  if (ChiSelect) {
-    value = PhotXm (code[0], average, secfilt);
-    if (value == -1) return (list);
-    if (value > ChiLimit) return (list);
-    value = PhotXm (code[1], average, secfilt);
-    if (value == -1) return (list);
-    if (value > ChiLimit) return (list);
-  }
+  /* exclude based on average parameters for both codes  */
+  if (!TestAverage (code[0], average, secfilt, measure)) return (list);
+  if (!TestAverage (code[1], average, secfilt, measure)) return (list);
+
+  Np0 = GetSelectionParam ();
+  Np1 = Np0 + 1;
 
   /* loop twice over all measures */
   for (i = 0; i < average[0].Nm; i++) {
+    SetSelectionParam (Np0);
     TESTCODE (code[0], measure[i]);
     TESTMEASURE (measure[i]);
@@ -1013,8 +945,13 @@
     for (j = 0; j < average[0].Nm; j++) {
       if (i == j) continue;
+      SetSelectionParam (Np1);
       TESTCODE (code[1], measure[j]);
       TESTMEASURE (measure[j]);
       SETMAG(M2, measure[j], mode[1]);
-      list[Nlist] = M1 - M2;
+      if ((M1 == NoMag) || (M2 == NoMag)) {
+	list[Nlist] = NoMag;
+      } else {
+	list[Nlist] = M1 - M2;
+      }
       Nlist ++;
     }
@@ -1024,4 +961,5 @@
 }
 
+/* extract a measurement list matching the number of dmag entries */
 double *ExtractByDMag (PhotCode **code, int *mode, Average *average, SecFilt *secfilt, Measure *measure, int *nlist, int param) {
 
@@ -1037,4 +975,5 @@
   }
 
+  /* one of the two entries results in a single element. extract the other */
   if (A1) {
     list = ExtractMeasures (code[1], mode[1], average, secfilt, measure, &N1, param);
@@ -1049,9 +988,8 @@
 double *ExtractMeasuresByDMag (PhotCode **code, int *mode, int use_first, Average *average, SecFilt *secfilt, Measure *measure, int *nlist, int param) {
 
-  int i, j, k, m, n, Nlist, NLIST;
+  int i, j, k, n, Nlist, NLIST;
   time_t t;
   int np, found;
   double *list, M1, M2, value;
-  double ra, dec, x, y;
   
   *nlist = 0; 
@@ -1082,85 +1020,5 @@
 
       /* assign value */
-      value = 0;
-      switch (param) {
-	case MEAS_MAG: /* magnitudes are already determined above */
-	  value = use_first ? M1 : M2;
-	  break;
-	case MEAS_RA: /* OK */
-	  value = average[0].R - measure[n].dR / 360000.0;
-	  break;
-	case MEAS_DEC: /* OK */
-	  value = average[0].D - measure[n].dD / 360000.0;
-	  break;
-	case MEAS_dMAG: /* OK */
-	  value = 0.001*measure[n].dM;
-	  break;
-	case MEAS_AIRMASS: /* OK */
-	  value = 0.001*measure[n].airmass;
-	  break;
-	case MEAS_EXPTIME: /* OK */
-	  value = pow (10.0, measure[n].dt * 0.0004);
-	  break;
-	case MEAS_PHOTCODE: /* OK */
-	  value = measure[n].source;
-	  break;
-	case MEAS_TIME: /* OK */
-	  value = TimeValue (measure[n].t, TimeReference, TimeFormat);
-	  break;
-	case MEAS_dR: /* OK */
-	  value = 0.01*measure[n].dR;
-	  break;
-	case MEAS_dD: /* OK */
-	  value = 0.01*measure[n].dD;
-	  break;
-	case MEAS_FWHM: /* OK */
-	  value = 0.01*measure[n].FWx;
-	  break;
-	case MEAS_DOPHOT: /* OK */
-	  value = measure[n].dophot;
-	  break;
-	case MEAS_FLAGS: /* ? */
-	  value = measure[n].flags;
-	  break;
-	case MEAS_XCCD: /* OK */
-	  ra  = average[0].R - measure[n].dR / 360000.0;
-	  dec = average[0].D - measure[n].dD / 360000.0;
-	  m = match_image_subset (image, subset, Nsubset, measure[n].t, measure[n].source);
-	  if (m == -1) break;
-	  if (!FindMosaicForImage (image, Nimage, m)) break;
-	  RD_to_XY (&x, &y, ra, dec, &image[m].coords);
-	  value = x;
-	  break;
-	case MEAS_YCCD: /* OK */
-	  ra  = average[0].R - measure[n].dR / 360000.0;
-	  dec = average[0].D - measure[n].dD / 360000.0;
-	  m = match_image_subset (image, subset, Nsubset, measure[n].t, measure[n].source);
-	  if (m == -1) break;
-	  if (!FindMosaicForImage (image, Nimage, m)) break;
-	  RD_to_XY (&x, &y, ra, dec, &image[m].coords);
-	  value = y;
-	  break;
-	case MEAS_XMOSAIC: /* OK */
-	  ra  = average[0].R - measure[n].dR / 360000.0;
-	  dec = average[0].D - measure[n].dD / 360000.0;
-	  m = match_image_subset (image, subset, Nsubset, measure[n].t, measure[n].source);
-	  if (m == -1) break;
-	  mosaic.crval1 = image[m].coords.crval1;
-	  mosaic.crval2 = image[m].coords.crval2;
-	  RD_to_XY (&x, &y, ra, dec, &mosaic);
-	  value = x;
-	  break;
-	case MEAS_YMOSAIC: /* OK */
-	  ra  = average[0].R - measure[n].dR / 360000.0;
-	  dec = average[0].D - measure[n].dD / 360000.0;
-	  m = match_image_subset (image, subset, Nsubset, measure[n].t, measure[n].source);
-	  if (m == -1) break;
-	  mosaic.crval1 = image[m].coords.crval1;
-	  mosaic.crval2 = image[m].coords.crval2;
-	  RD_to_XY (&x, &y, ra, dec, &mosaic);
-	  value = y;
-	  break;
-      }
-      list[Nlist] = value;
+      list[Nlist] = GetMeasure (param, &average[0], &measure[n], (use_first ? M1 : M2));
       Nlist ++;
     }
@@ -1170,2 +1028,91 @@
 }
 
+double GetMeasure (int param, Average *average, Measure *measure, double mag) {
+
+  int m;
+  double ra, dec, x, y;
+  double value;
+  Image *image;
+  Coords *mosaic;
+
+  value = 0;
+  switch (param) {
+    case MEAS_MAG: /* magnitudes are already determined above */
+      value = mag;
+      break;
+    case MEAS_RA: /* OK */
+      value = average[0].R - measure[0].dR / 360000.0;
+      break;
+    case MEAS_DEC: /* OK */
+      value = average[0].D - measure[0].dD / 360000.0;
+      break;
+    case MEAS_dMAG: /* OK */
+      value = 0.001*measure[0].dM;
+      break;
+    case MEAS_AIRMASS: /* OK */
+      value = 0.001*measure[0].airmass;
+      break;
+    case MEAS_EXPTIME: /* OK */
+      value = pow (10.0, measure[0].dt * 0.0004);
+      break;
+    case MEAS_PHOTCODE: /* OK */
+      value = measure[0].source;
+      break;
+    case MEAS_TIME: /* OK */
+      value = TimeValue (measure[0].t, TimeReference, TimeFormat);
+      break;
+    case MEAS_dR: /* OK */
+      value = 0.01*measure[0].dR;
+      break;
+    case MEAS_dD: /* OK */
+      value = 0.01*measure[0].dD;
+      break;
+    case MEAS_FWHM: /* OK */
+      value = 0.01*measure[0].FWx;
+      break;
+    case MEAS_DOPHOT: /* OK */
+      value = measure[0].dophot;
+      break;
+    case MEAS_FLAGS: /* ? */
+      value = measure[0].flags;
+      break;
+    case MEAS_XCCD: /* OK */
+      ra  = average[0].R - measure[0].dR / 360000.0;
+      dec = average[0].D - measure[0].dD / 360000.0;
+      image = MatchImage (measure[0].t, measure[0].source);
+      if (image == NULL) break;
+      RD_to_XY (&x, &y, ra, dec, &image[0].coords);
+      value = x;
+      break;
+    case MEAS_YCCD: /* OK */
+      ra  = average[0].R - measure[0].dR / 360000.0;
+      dec = average[0].D - measure[0].dD / 360000.0;
+      image = MatchImage (measure[0].t, measure[0].source);
+      if (image == NULL) break;
+      RD_to_XY (&x, &y, ra, dec, &image[0].coords);
+      value = y;
+      break;
+    case MEAS_XMOSAIC: /* OK */
+      ra  = average[0].R - measure[0].dR / 360000.0;
+      dec = average[0].D - measure[0].dD / 360000.0;
+      mosaic = MatchMosaic (measure[0].t, measure[0].source);
+      if (mosaic == NULL) break;
+      RD_to_XY (&x, &y, ra, dec, mosaic);
+      value = x;
+      break;
+    case MEAS_YMOSAIC: /* OK */
+      ra  = average[0].R - measure[0].dR / 360000.0;
+      dec = average[0].D - measure[0].dD / 360000.0;
+      mosaic = MatchMosaic (measure[0].t, measure[0].source);
+      if (mosaic == NULL) break;
+      RD_to_XY (&x, &y, ra, dec, mosaic);
+      value = y;
+      break;
+  }
+  return (value);
+}
+
+/** the mosaic entries do not use the registered mosaic found 
+    by MatchImage (via FindMosaicForImage).  Rather, they use
+    a coordinate frame saved by SetImageSelection 
+**/