Changeset 4585

Timestamp:

Jul 20, 2005, 4:27:25 AM (21 years ago)

Author:

eugene

Message:

substantial cleanup of photometry extraction functions

Location:

trunk/Ohana/src/opihi

Files:

• dvo/ImageOps.c (modified) (1 diff)
• dvo/ImageSelection.c (added)
• dvo/Makefile (modified) (1 diff)
• dvo/avextract.c (modified) (4 diffs)
• dvo/calextract.c (modified) (5 diffs)
• dvo/calmextract.c (modified) (5 diffs)
• dvo/ccd.c (modified) (7 diffs)
• dvo/cmd.c (modified) (7 diffs)
• dvo/ddmags.c (modified) (6 diffs)
• dvo/dmagaves.c (modified) (7 diffs)
• dvo/dmagmeas.c (modified) (5 diffs)
• dvo/dmags.c (modified) (7 diffs)
• dvo/elixir.c (modified) (1 diff)
• dvo/fitcolors.c (modified) (5 diffs)
• dvo/imdata.c (modified) (8 diffs)
• dvo/mextract.c (modified) (4 diffs)
• dvo/photometry.c (modified) (28 diffs)
• dvo/region_list.c (modified) (3 diffs)
• include/dvo1.h (modified) (1 diff)

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

@@ -2 +2 @@
 
 void image_subset (Image *image, int Nimage, int **Subset, int *Nsubset,
-              Graphdata *graph, int RegionSelect, 
-              unsigned long int tzero, double trange, int TimeSelect) 
+                   Graphdata *graph, int RegionSelect, 
+                   unsigned long int tzero, double trange, int TimeSelect) 
 {
 

trunk/Ohana/src/opihi/dvo/Makefile

@@ -28 +28 @@
 $(SDIR)/init.$(ARCH).o                  \
 $(SDIR)/ImageOps.$(ARCH).o              \
+$(SDIR)/ImageSelection.$(ARCH).o        \
 $(SDIR)/LoadImages.$(ARCH).o            \
 $(SDIR)/aregion.$(ARCH).o               \

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

@@ -3 +3 @@
 int avextract (int argc, char **argv) {
   
+  int i, j, m, N, NPTS, param;
+  int Ngraph, Nsec, Nregions, mode;
+  char filename[256], catdir[256], *RegionName, *RegionList, *p;
+
   Catalog catalog;
-  Graphdata graphsky;
   RegionFile *regions;
   PhotCode *code;
   Vector *vec;
 
-  char filename[256], catdir[256], *p;
-  int i, j, m, N, NPTS, param;
-  int Ngraph, Nsec, Nregions, mode;
-
+  /* defaults */
   regions = NULL;
   catalog.average = NULL; 
   catalog.secfilt = NULL;
   catalog.measure = NULL;
+  RegionName = NULL;
+  RegionList = NULL;
+  code = NULL;
+  mode = MAG_AVE;
 
   /* load photcode information */
@@ -22 +26 @@
   Nsec = GetPhotcodeNsecfilt ();
 
-  /* load data about plot windows */
-  Ngraph = 0;
-  if (!GetGraphData (&graphsky, NULL, &Ngraph)) goto escape;
- 
   /* find CATDIR in config system */
-  VarConfig ("CATDIR", "%s", catdir);
-
-  /* grab value of photcode, if available */
-  code = NULL;
-  mode = MAG_AVE;
-  if (N = get_argument (argc, argv, "-photcode")) {
-    remove_argument (N, &argc, argv);
-    code = GetPhotcodebyName (argv[N]);
-    if (code == NULL) {
-      fprintf (stderr, "photcode not found in photcode table\n");
-      goto usage;
-    }
-    remove_argument (N, &argc, argv);
-  }
+  if (VarConfig ("CATDIR", "%s", catdir) == NULL) goto escape;
 
   /* interpret command-line options */
-  if (!SetPhotSelections (&argc, argv)) goto usage;
-  if (argc != 3) { goto usage; }
+  SetSelectionParam (0);
+  if (!SetRegionSelection (&argc, argv, &RegionName, &RegionList)) goto escape;
+  if (!SetPhotSelections (&argc, argv, 1)) goto usage;
 
-  /* identify selection */
-  param = GetAverageParam (argv[2]);
+  /* interpret required command-line arguments: mextract (value) */
+  if (argc != 2) { goto usage; }
+  param = GetAverageParam (argv[1]);
   if (param == AVE_ZERO) {
-    if (!GetPhotcodeInfo (argv[2], &code, &mode)) {
+    if (!GetPhotcodeInfo (argv[1], &code, &mode)) {
       GetAverageParamHelp ();
       goto escape;
     }
     param = AVE_MAG;
-    for (p = argv[2]; *p != 0; p++) { 
+    for (p = argv[1]; *p != 0; p++) { 
       if (*p == '.') *p = ':';
     }
   } 
-  if ((vec = SelectVector (argv[2], ANYVECTOR, TRUE)) == NULL) goto escape;
+  if (!TestPhotSelections (&code, param)) goto escape;
 
-  /* Need valid code for some params: */
-  switch (param) {
-    case AVE_MAG:
-    case AVE_dMAG:
-    case AVE_Xm:
-    case AVE_TYPE:
-    case AVE_NCODE:
-    case AVE_NPHOT: 
-      if (code == NULL) goto escape;
-      if (code[0].type == PHOT_DEP) goto escape;
-      if (code[0].type == PHOT_REF) goto escape;
-      break;
-    default:
-      break;
-  }
-
-  /* determine region-file names */
-  if (!strcmp (argv[1], "all")) {
-    float Radius;
-    Radius = MAX (fabs(graphsky.xmax), fabs(graphsky.ymax));
-    regions = find_regions (graphsky.coords.crval1, graphsky.coords.crval2, Radius, &Nregions);
-    if (!SetImageSelection (param, &graphsky, TRUE)) goto escape;
-    /* this is not really used by avextract -- no parameters require images */
-  } else {
-    Nregions = 1;
-    ALLOCATE (regions, RegionFile, 1);
-    strcpy (regions[0].name, argv[1]);
-    if (!SetImageSelection (param, NULL, FALSE)) goto escape;
-  }
+  /* load region corresponding to selection above */
+  if ((regions = SelectRegions (RegionName, RegionList, &Nregions)) == NULL) goto escape;
 
   /* create storage vector */
   N = 0;
-  NPTS = 1000;
-  REALLOCATE (vec[0].elements, float, NPTS);
+  NPTS = 1;
+  if ((vec = SelectVector (argv[1], ANYVECTOR, TRUE)) == NULL) goto escape;
 
   for (i = 0; i < Nregions; i++) {
@@ -115 +77 @@
       vec[0].elements[N] = ExtractAverages (code, mode, &catalog.average[j], &catalog.secfilt[j*Nsec], &catalog.measure[m], param);
       N++;
-      if (N == NPTS - 1) {
-        NPTS += 2000;
-        REALLOCATE (vec[0].elements, float, NPTS);
-      }
+      CHECK_REALLOCATE (vec[0].elements, float, NPTS, N, 2000);
     }
     if (catalog.average != NULL) free (catalog.average);
@@ -131 +90 @@
 
   if (regions != NULL) free (regions);
-  FreeImageSelection ();
   return (TRUE);
 
 usage:
-  fprintf (stderr, "USAGE: avextract (from) (value) [options]\n");
-  fprintf (stderr, "  from: cpt name or 'all'\n");
+  fprintf (stderr, "USAGE: avextract (value) [options]\n");
   fprintf (stderr, "  value: average.parameter or photcode\n");
-  return (FALSE);
 
 escape:
-  fprintf (stderr, "error in avextract\n");
   if (regions != NULL) free (regions);
-  FreeImageSelection ();
   if (catalog.average != NULL) free (catalog.average);
   if (catalog.secfilt != NULL) free (catalog.secfilt);
   if (catalog.measure != NULL) free (catalog.measure);
+  if (RegionName != NULL) free (RegionName);
+  if (RegionList != NULL) free (RegionList);
   return (FALSE);
 }

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

@@ -5 +5 @@
 int calextract (int argc, char **argv) {
   
+  int i, j, m, N, Nr, Nregions, mode[2];
+  int Nsec, NSTAR;
+  char filename[256], catdir[256], *RegionName, *RegionList;
+  double Radius, M1, M2, dM2, color;
+
   Catalog catalog;
-  Graphdata graphsky;
   RegionFile *regions;
   PhotCode *code[2];
   Vector **vec;
-
-  int i, j, m, N, Nr, Nregions, mode[2];
-  int RegionList, Ngraph, Nsec, NSTAR;
-  char filename[256], catdir[256], RegionListFile[256];
-  double Radius, M1, M2, dM2, color;
 
   /* these need to be freed in the end */
@@ -21 +20 @@
   catalog.secfilt = NULL;
   catalog.measure = NULL;
+  RegionName = NULL;
+  RegionList = NULL;
 
   /* load photcode information */
@@ -26 +27 @@
   Nsec = GetPhotcodeNsecfilt ();
 
-  /* load data about plot windows */
-  Ngraph = 0;
-  if (!GetGraphData (&graphsky, NULL, &Ngraph)) goto escape;
-
   /* find CATDIR in config system */
-  VarConfig ("CATDIR", "%s", catdir);
-
-  /* specify catalog files */
-  RegionList = FALSE;
-  if (N = get_argument (argc, argv, "-list")) {
-    RegionList = TRUE;
-    remove_argument (N, &argc, argv);
-    strcpy (RegionListFile, argv[N]);
-    remove_argument (N, &argc, argv);
-  }
+  if (VarConfig ("CATDIR", "%s", catdir) == NULL) goto escape;
 
   /* command line arguments */
-  if (!SetPhotSelections (&argc, argv)) goto usage;
+  SetSelectionParam (0);
+  if (!SetRegionSelection (&argc, argv, &RegionName, &RegionList)) goto escape;
+  if (!SetPhotSelections (&argc, argv, 2)) goto usage;
+
   if (argc != 4) goto usage;
   if (strcmp (argv[2], "-")) goto usage;
@@ -66 +57 @@
   if ((vec[Nd2] = SelectVector ("cal:dm2",      ANYVECTOR, TRUE)) == NULL) goto escape;
 
-  /* select region files */
-  if (RegionList) {
-    regions = region_list (RegionListFile, &Nregions);
-  } else {
-    Radius = MAX (fabs(graphsky.xmax), fabs(graphsky.ymax));
-    regions = find_regions (graphsky.coords.crval1, graphsky.coords.crval2, Radius, &Nregions);
-  }
+  /* load region corresponding to selection above */
+  if ((regions = SelectRegions (RegionName, RegionList, &Nregions)) == NULL) goto escape;
 
   for (Nr = 0; Nr < Nregions; Nr++) {
@@ -150 +136 @@
  escape:
   
+  if (RegionName != NULL) free (RegionName);
+  if (RegionList != NULL) free (RegionList);
   if (regions != NULL) free (regions);  
   for (i = 0; i < NVEC; i++) {

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

@@ -6 +6 @@
 int calmextract (int argc, char **argv) {
   
+  int i, j, k, m, N, N1, Nr, mode[2];
+  int NSTAR, Nstar, Nsec, Nregions;
+  char filename[256], catdir[256], *RegionName, *RegionList;
+  double *M1, M2, dM2, color, Radius;
+
   Catalog catalog;
-  Graphdata graphsky;
   RegionFile *regions;
   PhotCode *code[2];
   Vector **vec;
-
-  int i, j, k, m, N, N1, Nr, mode[2];
-  int Ngraph, NSTAR, Nstar, Nsec, Nregions, RegionList;
-  char filename[256], catdir[256], RegionListFile[256];
-  double *M1, M2, dM2, color, Radius;
 
   /* these need to be freed in the end */
@@ -22 +21 @@
   catalog.secfilt = NULL;
   catalog.measure = NULL;
+  RegionName = NULL;
+  RegionList = NULL;
 
   /* load photcode information */
-  if (!InitPhotcodes ()) return (FALSE);
+  if (!InitPhotcodes ()) goto escape;
   Nsec = GetPhotcodeNsecfilt ();
 
-  /* load data about plot windows */
-  Ngraph = 0;
-  if (!GetGraphData (&graphsky, NULL, &Ngraph)) goto escape;
-
   /* find CATDIR in config system */
-  VarConfig ("CATDIR", "%s", catdir);
-
-  /* specify catalog files */
-  RegionList = FALSE;
-  if (N = get_argument (argc, argv, "-list")) {
-    RegionList = TRUE;
-    remove_argument (N, &argc, argv);
-    strcpy (RegionListFile, argv[N]);
-    remove_argument (N, &argc, argv);
-  }
+  if (VarConfig ("CATDIR", "%s", catdir) == NULL) goto escape;
 
   /* command line arguments */
-  if (!SetPhotSelections (&argc, argv)) goto usage;
+  if (!SetRegionSelection (&argc, argv, &RegionName, &RegionList)) goto escape;
+  if (!SetPhotSelections (&argc, argv, 2)) goto usage;
+
+  /* interpret required command-line arguments: calmextract F1 - F2 */
   if (argc != 4) goto usage;
   if (strcmp (argv[2], "-")) goto usage;
   if (!GetPhotcodeInfo (argv[1], &code[0], &mode[0])) goto usage;
   if (!GetPhotcodeInfo (argv[3], &code[1], &mode[1])) goto usage;
-  /* code.type must be PHOT_REF */
-
-  /* test PhotSelections: is photcode specified if needed? */
-  if (!TestPhotSelections (&code)) {
-    fprintf (stderr, "photcode selection rules violated\n");
-    fprintf (stderr, "average and ensemble restrictions require a corresponding photcode\n");
-    goto escape;
-  }
-
+  if (!TestPhotSelections (&code[0], MEAS_ZERO)) goto escape;
 
   /* returned vectors are dmag, mag, color, time, airmass, ra, dec, x, y, exptime */
@@ -86 +69 @@
   }
 
-  /* select region files */
-  if (RegionList) {
-    regions = region_list (RegionListFile, &Nregions);
-    if (!SetImageSelection (MEAS_XMOSAIC, NULL, FALSE)) goto escape;
-  } else {
-    Radius = MAX (fabs(graphsky.xmax), fabs(graphsky.ymax));
-    regions = find_regions (graphsky.coords.crval1, graphsky.coords.crval2, Radius, &Nregions);
-    if (!SetImageSelection (MEAS_XMOSAIC, &graphsky, TRUE)) goto escape;
-  }
+  /* load region corresponding to selection above */
+  if ((regions = SelectRegions (RegionName, RegionList, &Nregions)) == NULL) goto escape;
+  if (!SetImageSelection (MEAS_XMOSAIC)) goto escape;
 
   for (Nr = 0; Nr < Nregions; Nr++) {
@@ -177 +154 @@
     vec[i][0].Nelements = N;
   }
+  FreeImageSelection ();
   return (TRUE);
   
@@ -185 +163 @@
  escape:
   
+  FreeImageSelection ();
   if (regions != NULL) free (regions);  
   for (i = 0; i < NVEC; i++) {

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

@@ -3 +3 @@
 int ccd (int argc, char **argv) {
   
-  Catalog catalog;
-  Graphdata graphsky;
-  RegionFile *regions;
-  PhotCode *code[4];
-  Vector *xvec, *yvec;
-
-  char filename[256], catdir[256], *CPTfile;
+  char filename[256], catdir[256], *RegionName, *RegionList;
   double Radius;
   double *M1, *M2;
   int i, m, k, Npts, NPTS, N;
   int N1, N2, i1, i2, mode[4];
-  int Ngraph, Nsec, Nregions, UseAverages;
+  int Ngraph, Nsec, Nregions;
+
+  Catalog catalog;
+  RegionFile *regions;
+  PhotCode *code[4];
+  Vector *xvec, *yvec;
 
   /* defaults */
@@ -21 +20 @@
   catalog.secfilt = NULL;
   catalog.measure = NULL;
+  RegionName = NULL;
+  RegionList = NULL;
 
   /* load photcode information */
@@ -26 +27 @@
   Nsec = GetPhotcodeNsecfilt ();
 
-  /* load data about plot windows */
-  Ngraph = 0;
-  if (!GetGraphData (&graphsky, NULL, &Ngraph)) return (FALSE);
- 
   /* find CATDIR in config system */
-  VarConfig ("CATDIR", "%s", catdir);
-
-  /* check for CPT selection */
-  CPTfile = NULL;
-  if (N = get_argument (argc, argv, "-cpt")) {
-    remove_argument (N, &argc, argv);
-    CPTfile = strcreate (argv[N]);
-    remove_argument (N, &argc, argv);
-  }    
+  if (VarConfig ("CATDIR", "%s", catdir) == NULL) goto escape;
 
   /* interpret command-line options */
-  if (!SetPhotSelections (&argc, argv)) goto usage;
+  if (!SetRegionSelection (&argc, argv, &RegionName, &RegionList)) goto escape;
+  if (!SetPhotSelections (&argc, argv, 4)) goto usage;
+
+  /* interpret command-line options */
   if (argc != 8) goto usage;
   if (strcmp (argv[2], "-")) goto usage;
   if (strcmp (argv[4], ":")) goto usage;
   if (strcmp (argv[6], "-")) goto usage;
-
-  /* interpret command-line options */
   if (!GetPhotcodeInfo (argv[1], &code[0], &mode[0])) return (FALSE);
   if (!GetPhotcodeInfo (argv[3], &code[1], &mode[1])) return (FALSE);
   if (!GetPhotcodeInfo (argv[5], &code[2], &mode[2])) return (FALSE);
   if (!GetPhotcodeInfo (argv[7], &code[3], &mode[3])) return (FALSE);
-  if ((mode[2] == MAG_AVE) || (mode[2] == MAG_REF)) UseAverages = TRUE;
+  if (!TestPhotSelections (&code[0], MEAS_ZERO)) goto escape;
 
-  /* determine region-file names */
-  if (CPTfile == NULL) {
-    float Radius;
-    Radius = MAX (fabs(graphsky.xmax), fabs(graphsky.ymax));
-    regions = find_regions (graphsky.coords.crval1, graphsky.coords.crval2, Radius, &Nregions);
-  } else {
-    Nregions = 1;
-    ALLOCATE (regions, RegionFile, 1);
-    strcpy (regions[0].name, CPTfile);
-  }
+  /* load region corresponding to selection above */
+  if ((regions = SelectRegions (RegionName, RegionList, &Nregions)) == NULL) goto escape;
 
   /* init vectors to save data */
@@ -94 +76 @@
       m = catalog.average[i].offset;
 
+      SetSelectionParam (0);
       M1 = ExtractDMag (&code[0], &mode[0], &catalog.average[i], &catalog.secfilt[i*Nsec], &catalog.measure[m], &N1);
       if (N1 == 0) goto skip;
 
+      SetSelectionParam (2);
       M2 = ExtractDMag (&code[2], &mode[2], &catalog.average[i], &catalog.secfilt[i*Nsec], &catalog.measure[m], &N2);
       if (N2 == 0) goto skip;
@@ -105 +89 @@
           yvec[0].elements[Npts] = M2[i2];
           Npts++;
-          if (Npts == NPTS) {
+          if (Npts >= NPTS) {
             NPTS += 2000;
             REALLOCATE (xvec[0].elements, float, NPTS);
@@ -128 +112 @@
 
 usage:
-  fprintf (stderr, "USAGE: ccd F - F : measure.param\n");
+  fprintf (stderr, "USAGE: ccd F - F : F - F\n");
   return (FALSE);
 
@@ -136 +120 @@
   if (catalog.secfilt != NULL) free (catalog.secfilt);
   if (catalog.measure != NULL) free (catalog.measure);
+  if (RegionName != NULL) free (RegionName);
+  if (RegionList != NULL) free (RegionList);
   return (FALSE);
 }

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

@@ -3 +3 @@
 int cmd (int argc, char **argv) { /* really need to think about upper limits & how to represent them */
   
+  char filename[256], catdir[256], *RegionName, *RegionList;
+  double Radius, *M1, *M3;
+  int i, j, m, i1, i3, N1, N3, N;
+  int Npts, NPTS, mode[3];
+  int Ngraph, Nsec, Nregions;
+
   Catalog catalog;
   Graphdata graphsky;
@@ -9 +15 @@
   Vector *xvec, *yvec;
 
-  char filename[256], catdir[256], *CPTfile;
-  double Radius, *M1, *M3;
-  int i, j, m, i1, i3, N1, N3, N;
-  int Npts, NPTS, mode[3];
-  int Ngraph, Nsec, Nregions;
-
   /* defaults */
   regions = NULL;
@@ -20 +20 @@
   catalog.secfilt = NULL;
   catalog.measure = NULL;
+  RegionName = NULL;
+  RegionList = NULL;
 
   /* load photcode information */
@@ -25 +27 @@
   Nsec = GetPhotcodeNsecfilt ();
 
-  /* load data about plot windows */
-  Ngraph = 0;
-  if (!GetGraphData (&graphsky, NULL, &Ngraph)) return (FALSE);
- 
   /* find CATDIR in config system */
-  VarConfig ("CATDIR", "%s", catdir);
-
-  /* check for CPT selection */
-  CPTfile = NULL;
-  if (N = get_argument (argc, argv, "-cpt")) {
-    remove_argument (N, &argc, argv);
-    CPTfile = strcreate (argv[N]);
-    remove_argument (N, &argc, argv);
-  }    
+  if (VarConfig ("CATDIR", "%s", catdir) == NULL) goto escape;
 
   /* interpret command-line options */
-  if (!SetPhotSelections (&argc, argv)) goto usage;
+  if (!SetRegionSelection (&argc, argv, &RegionName, &RegionList)) goto escape;
+  if (!SetPhotSelections (&argc, argv, 3)) goto usage;
+
+  /* interpret command-line options */
   if (argc != 6) { goto usage; }
   if (strcmp (argv[2], "-")) goto usage;
   if (strcmp (argv[4], ":")) goto usage;
-
-  /* interpret command-line options */
   if (!GetPhotcodeInfo (argv[1], &code[0], &mode[0])) return (FALSE);
   if (!GetPhotcodeInfo (argv[3], &code[1], &mode[1])) return (FALSE);
   if (!GetPhotcodeInfo (argv[5], &code[2], &mode[2])) return (FALSE);
+  if (!TestPhotSelections (&code[0], MEAS_ZERO)) goto escape;
 
-  /* determine region-file names */
-  if (CPTfile == NULL) {
-    float Radius;
-    Radius = MAX (fabs(graphsky.xmax), fabs(graphsky.ymax));
-    regions = find_regions (graphsky.coords.crval1, graphsky.coords.crval2, Radius, &Nregions);
-  } else {
-    Nregions = 1;
-    ALLOCATE (regions, RegionFile, 1);
-    strcpy (regions[0].name, CPTfile);
-  }
+  /* load region corresponding to selection above */
+  if ((regions = SelectRegions (RegionName, RegionList, &Nregions)) == NULL) goto escape;
 
   /* init vectors to save data */
@@ -90 +74 @@
       m = catalog.average[i].offset;
 
+      SetSelectionParam (0);
       M1 = ExtractDMag (code, mode, &catalog.average[i], &catalog.secfilt[i*Nsec], &catalog.measure[m], &N1);
       if (N1 == 0) goto skip;
 
+      SetSelectionParam (2);
       M3 = ExtractMagnitudes (code[2], mode[2], &catalog.average[i], &catalog.secfilt[i*Nsec], &catalog.measure[m], &N3);
       if (N3 == 0) goto skip;
@@ -101 +87 @@
           yvec[0].elements[Npts] = M3[i3];
           Npts++;
-          if (Npts == NPTS) {
+          if (Npts >= NPTS) {
             NPTS += 2000;
             REALLOCATE (xvec[0].elements, float, NPTS);
@@ -132 +118 @@
   if (catalog.secfilt != NULL) free (catalog.secfilt);
   if (catalog.measure != NULL) free (catalog.measure);
+  if (RegionName != NULL) free (RegionName);
+  if (RegionList != NULL) free (RegionList);
   return (FALSE);
 }

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

@@ -3 +3 @@
 int ddmags (int argc, char **argv) {
   
-  Catalog catalog;
-  Graphdata graphsky;
-  RegionFile *regions;
-  PhotCode *code[4];
-  Vector *xvec, *yvec;
-
-  char filename[256], catdir[256];
+  char filename[256], catdir[256], *RegionName, *RegionList;
   double Radius;
   double *M1, *M2;
   int i, m, k, Npts, NPTS;
   int N1, N2, i1, i2, mode[4];
-  int Ngraph, Nsec, Nregions, UseAverages;
+  int Ngraph, Nsec, Nregions;
+
+  Catalog catalog;
+  RegionFile *regions;
+  PhotCode *code[4];
+  Vector *xvec, *yvec;
 
   /* defaults */
@@ -21 +20 @@
   catalog.secfilt = NULL;
   catalog.measure = NULL;
-
-  /* find CATDIR in config system */
-  VarConfig ("CATDIR", "%s", catdir);
+  RegionName = NULL;
+  RegionList = NULL;
 
   /* load photcode information */
@@ -29 +27 @@
   Nsec = GetPhotcodeNsecfilt ();
 
-  /* load data about plot windows */
-  Ngraph = 0;
-  if (!GetGraphData (&graphsky, NULL, &Ngraph)) return (FALSE);
- 
-  /* find catalog files which overlap this region */
-  Radius = MAX (fabs(graphsky.xmax), fabs(graphsky.ymax));
-  regions = find_regions (graphsky.coords.crval1, graphsky.coords.crval2, Radius, &Nregions);
+  /* find CATDIR in config system */
+  if (VarConfig ("CATDIR", "%s", catdir) == NULL) goto escape;
 
   /* interpret command-line options */
-  if (!SetPhotSelections (&argc, argv)) goto usage;
+  if (!SetRegionSelection (&argc, argv, &RegionName, &RegionList)) goto escape;
+  if (!SetPhotSelections (&argc, argv, 4)) goto usage;
+
+  /* interpret command-line options */
   if (argc != 8) goto usage;
   if (strcmp (argv[2], "-")) goto usage;
   if (strcmp (argv[4], ":")) goto usage;
   if (strcmp (argv[6], "-")) goto usage;
-
-  /* interpret command-line options */
   if (!GetPhotcodeInfo (argv[1], &code[0], &mode[0])) return (FALSE);
   if (!GetPhotcodeInfo (argv[3], &code[1], &mode[1])) return (FALSE);
   if (!GetPhotcodeInfo (argv[5], &code[2], &mode[2])) return (FALSE);
   if (!GetPhotcodeInfo (argv[7], &code[3], &mode[3])) return (FALSE);
-  if ((mode[2] == MAG_AVE) || (mode[2] == MAG_REF)) UseAverages = TRUE;
+  if (!TestPhotSelections (&code[0], MEAS_ZERO)) goto escape;
+
+  /* load region corresponding to selection above */
+  if ((regions = SelectRegions (RegionName, RegionList, &Nregions)) == NULL) goto escape;
 
   /* init vectors to save data */
@@ -79 +76 @@
       m = catalog.average[i].offset;
 
+      SetSelectionParam (0);
       M1 = ExtractDMag (&code[0], &mode[0], &catalog.average[i], &catalog.secfilt[i*Nsec], &catalog.measure[m], &N1);
       if (N1 == 0) goto skip;
 
+      SetSelectionParam (2);
       M2 = ExtractDMag (&code[2], &mode[2], &catalog.average[i], &catalog.secfilt[i*Nsec], &catalog.measure[m], &N2);
       if (N2 == 0) goto skip;
@@ -90 +89 @@
           yvec[0].elements[Npts] = M2[i2];
           Npts++;
-          if (Npts == NPTS) {
+          if (Npts >= NPTS) {
             NPTS += 2000;
             REALLOCATE (xvec[0].elements, float, NPTS);
@@ -120 +119 @@
   if (catalog.secfilt != NULL) free (catalog.secfilt);
   if (catalog.measure != NULL) free (catalog.measure);
+  if (RegionName != NULL) free (RegionName);
+  if (RegionList != NULL) free (RegionList);
   return (FALSE);
 }

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

@@ -3 +3 @@
 int dmagaves (int argc, char **argv) {
   
-  Catalog catalog;
-  Graphdata graphsky;
-  RegionFile *regions;
-  PhotCode *code[3];
-  Vector *xvec, *yvec;
-
-  char filename[256], catdir[256], *CPTfile;
+  char filename[256], catdir[256], *RegionName, *RegionList;
   double Radius, *M1, M2;
   int i, j, k, m, N1, N;
   int Npts, NPTS, param, mode[3];
   int Ngraph, Nsec, Nregions;
+
+  Catalog catalog;
+  RegionFile *regions;
+  PhotCode *code[3];
+  Vector *xvec, *yvec;
 
   /* defaults */
@@ -20 +19 @@
   catalog.secfilt = NULL;
   catalog.measure = NULL;
+  RegionName = NULL;
+  RegionList = NULL;
 
   /* load photcode information */
@@ -25 +26 @@
   Nsec = GetPhotcodeNsecfilt ();
 
-  /* load data about plot windows */
-  Ngraph = 0;
-  if (!GetGraphData (&graphsky, NULL, &Ngraph)) return (FALSE);
- 
   /* find CATDIR in config system */
-  VarConfig ("CATDIR", "%s", catdir);
-
-  /* check for CPT selection */
-  CPTfile = NULL;
-  if (N = get_argument (argc, argv, "-cpt")) {
-    remove_argument (N, &argc, argv);
-    CPTfile = strcreate (argv[N]);
-    remove_argument (N, &argc, argv);
-  }    
+  if (VarConfig ("CATDIR", "%s", catdir) == NULL) goto escape;
 
   /* interpret command-line options */
-  if (!SetPhotSelections (&argc, argv)) goto usage;
+  if (!SetRegionSelection (&argc, argv, &RegionName, &RegionList)) goto escape;
+  if (!SetPhotSelections (&argc, argv, 3)) goto usage;
+
+  /* interpret command-line options: dmagaves F1 - F2 : (value) */
   if (argc != 6) { goto usage; }
   if (strcmp (argv[2], "-")) goto usage;
   if (strcmp (argv[4], ":")) goto usage;
+  if (!GetPhotcodeInfo (argv[1], &code[0], &mode[0])) goto usage;
+  if (!GetPhotcodeInfo (argv[3], &code[1], &mode[1])) goto usage;
+  if ((param = GetAverageParam (argv[5])) == AVE_ZERO) goto usage;
+  if (!TestPhotSelections (&code[2], param)) goto escape;
 
-  /* interpret command-line options */
-  if (!GetPhotcodeInfo (argv[1], &code[0], &mode[0])) return (FALSE);
-  if (!GetPhotcodeInfo (argv[3], &code[1], &mode[1])) return (FALSE);
-  if ((param = GetAverageParam (argv[5])) == AVE_ZERO) {
-    GetAverageParamHelp ();
-    return (FALSE);
-  }
-  code[2] = GetPhotcodeEquivbyCode (code[0][0].code);
-  mode[2] = mode[0];
-
-  /* test PhotSelections: is photcode specified if needed? */
-  for (i = 0; i < 3; i++) {
-    if (!TestPhotSelections (&code[i])) {
-      fprintf (stderr, "photcode selection rules violated\n");
-      fprintf (stderr, "average and ensemble restrictions require a corresponding photcode\n");
-      goto escape;
-    }
-  }
-
-  /* determine region-file names */
-  if (CPTfile == NULL) {
-    float Radius;
-    Radius = MAX (fabs(graphsky.xmax), fabs(graphsky.ymax));
-    regions = find_regions (graphsky.coords.crval1, graphsky.coords.crval2, Radius, &Nregions);
-  } else {
-    Nregions = 1;
-    ALLOCATE (regions, RegionFile, 1);
-    strcpy (regions[0].name, CPTfile);
-  }
+  /* load region corresponding to selection above */
+  if ((regions = SelectRegions (RegionName, RegionList, &Nregions)) == NULL) goto escape;
 
   /* init vectors to save data */
@@ -103 +72 @@
       m = catalog.average[i].offset;
 
+      SetSelectionParam (0);
       M1 = ExtractDMag (code, mode, &catalog.average[i], &catalog.secfilt[i*Nsec], &catalog.measure[m], &N1);
       if (N1 == 0) goto skip;
 
+      SetSelectionParam (2);
       M2 = ExtractAverages (code[2], mode[2], &catalog.average[i], &catalog.secfilt[i*Nsec], &catalog.measure[m], param);
 
@@ -112 +83 @@
         yvec[0].elements[Npts] = M2;
         Npts++;
-        if (Npts == NPTS) {
+        if (Npts >= NPTS) {
           NPTS += 2000;
           REALLOCATE (xvec[0].elements, float, NPTS);
@@ -128 +99 @@
     catalog.measure = (Measure *) NULL;
   }
-  FreeImageSelection ();
   if (regions != NULL) free (regions);
   xvec[0].Nelements = yvec[0].Nelements = Npts;
@@ -138 +108 @@
 
 escape:
-  FreeImageSelection ();
   if (regions != NULL) free (regions);
   if (catalog.average != NULL) free (catalog.average);
   if (catalog.secfilt != NULL) free (catalog.secfilt);
   if (catalog.measure != NULL) free (catalog.measure);
+  if (RegionName != NULL) free (RegionName);
+  if (RegionList != NULL) free (RegionList);
   return (FALSE);
 }

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

@@ -1 +1 @@
 # include "dvo1.h"
-double *ExtractByDMag (PhotCode **code, int *mode, Average *average, SecFilt *secfilt, Measure *measure, int *nlist, int param);
 
 int dmagmeas (int argc, char **argv) {
   
+  char filename[256], catdir[256], *RegionName, *RegionList;
+  double Radius, *M1, *M3;
+  int i, j, k, m, i1, i3, N1, N3, N;
+  int Npts, NPTS, param, mode[3];
+  int Ngraph, Nsec, Nregions;
+
   Catalog catalog;
-  Graphdata graphsky;
   RegionFile *regions;
   PhotCode *code[3];
   Vector *xvec, *yvec;
-
-  char filename[256], catdir[256], *CPTfile;
-  double Radius, *M1, *M3;
-  int i, j, k, m, N1, N3, N;
-  int Npts, NPTS, param, mode[3];
-  int Ngraph, Nsec, Nregions;
 
   /* defaults */
@@ -21 +19 @@
   catalog.secfilt = NULL;
   catalog.measure = NULL;
+  RegionName = NULL;
+  RegionList = NULL;
 
   /* load photcode information */
@@ -26 +26 @@
   Nsec = GetPhotcodeNsecfilt ();
 
-  /* load data about plot windows */
-  Ngraph = 0;
-  if (!GetGraphData (&graphsky, NULL, &Ngraph)) return (FALSE);
- 
   /* find CATDIR in config system */
-  VarConfig ("CATDIR", "%s", catdir);
-
-  /* check for CPT selection */
-  CPTfile = NULL;
-  if (N = get_argument (argc, argv, "-cpt")) {
-    remove_argument (N, &argc, argv);
-    CPTfile = strcreate (argv[N]);
-    remove_argument (N, &argc, argv);
-  }    
+  if (VarConfig ("CATDIR", "%s", catdir) == NULL) goto escape;
 
   /* interpret command-line options */
-  if (!SetPhotSelections (&argc, argv)) goto usage;
+  if (!SetRegionSelection (&argc, argv, &RegionName, &RegionList)) goto escape;
+  if (!SetPhotSelections (&argc, argv, 3)) goto usage;
+
+  /* interpret command-line arguments: dmagmeas F1 - F2 : (value) */
   if (argc != 6) { goto usage; }
   if (strcmp (argv[2], "-")) goto usage;
   if (strcmp (argv[4], ":")) goto usage;
-
-  /* interpret command-line options */
   if (!GetPhotcodeInfo (argv[1], &code[0], &mode[0])) goto usage;
   if (!GetPhotcodeInfo (argv[3], &code[1], &mode[1])) goto usage;
   if ((param = GetMeasureParam (argv[5])) == MEAS_ZERO) goto usage;
+  if (!TestPhotSelections (&code[2], MEAS_ZERO)) goto escape;
 
-  /* test PhotSelections: is photcode specified if needed? */
-  for (i = 0; i < 2; i++) {
-    if (!TestPhotSelections (&code[i])) {
-      fprintf (stderr, "photcode selection rules violated\n");
-      fprintf (stderr, "average and ensemble restrictions require a corresponding photcode\n");
-      goto escape;
-    }
-  }
-
-  /* determine region-file names */
-  if (CPTfile == NULL) {
-    float Radius;
-    Radius = MAX (fabs(graphsky.xmax), fabs(graphsky.ymax));
-    regions = find_regions (graphsky.coords.crval1, graphsky.coords.crval2, Radius, &Nregions);
-  } else {
-    Nregions = 1;
-    ALLOCATE (regions, RegionFile, 1);
-    strcpy (regions[0].name, CPTfile);
-  }
-  if (!SetImageSelection (param, &graphsky, TRUE)) return (FALSE);
+  /* load region corresponding to selection above */
+  if ((regions = SelectRegions (RegionName, RegionList, &Nregions)) == NULL) goto escape;
+  if (!SetImageSelection (param)) goto escape;
 
   /* init vectors to save data */
@@ -101 +74 @@
       m = catalog.average[i].offset;
 
+      SetSelectionParam (0);
       M1 = ExtractDMag (code, mode, &catalog.average[i], &catalog.secfilt[i*Nsec], &catalog.measure[m], &N1);
       if (N1 == 0) goto skip;
 
-      M3 = ExtractByDMag (code, mode, &catalog.average[i], &catalog.secfilt[i*Nsec], &catalog.measure[m], &N3, param);
-      if (N3 != N1) { fprintf (stderr, "mismatch: programming error\n"); goto escape; }
+      SetSelectionParam (2);
+      M3 = ExtractMeasures (code[2], mode[2], &catalog.average[i], &catalog.secfilt[i*Nsec], &catalog.measure[m], &N3, param);
+      if (N3 == 0) goto skip;
 
-      for (k = 0; k < N1; k++) {
-        xvec[0].elements[Npts] = M1[k];
-        yvec[0].elements[Npts] = M3[k];
-        Npts++;
-        if (Npts == NPTS) {
-          NPTS += 2000;
-          REALLOCATE (xvec[0].elements, float, NPTS);
-          REALLOCATE (yvec[0].elements, float, NPTS);
+      for (i1 = 0; i1 < N1; i1++) {
+        for (i3 = 0; i3 < N3; i3++) {
+          xvec[0].elements[Npts] = M1[i1];
+          yvec[0].elements[Npts] = M3[i3];
+          Npts++;
+          if (Npts >= NPTS) {
+            NPTS += 2000;
+            REALLOCATE (xvec[0].elements, float, NPTS);
+            REALLOCATE (yvec[0].elements, float, NPTS);
+          }
         }
       }
@@ -143 +120 @@
   if (catalog.secfilt != NULL) free (catalog.secfilt);
   if (catalog.measure != NULL) free (catalog.measure);
+  if (RegionName != NULL) free (RegionName);
+  if (RegionList != NULL) free (RegionList);
   return (FALSE);
 }

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

@@ -3 +3 @@
 int dmags (int argc, char **argv) {
   
-  Catalog catalog;
-  Graphdata graphsky;
-  RegionFile *regions;
-  PhotCode *code[3];
-  Vector *xvec, *yvec;
-
-  char filename[256], catdir[256], *CPTfile;
+  char filename[256], catdir[256], *RegionName, *RegionList;
   double Radius, *M1, *M3;
   int i, j, m, i1, i3, N1, N3, N;
   int Npts, NPTS, mode[3];
-  int Ngraph, Nsec, Nregions;
+  int Nsec, Nregions;
+
+  PhotCode *code[3];
+  Catalog catalog;
+  RegionFile *regions;
+  Vector *xvec, *yvec;
 
   /* defaults */
@@ -20 +19 @@
   catalog.secfilt = NULL;
   catalog.measure = NULL;
+  RegionName = NULL;
+  RegionList = NULL;
 
   /* load photcode information */
@@ -25 +26 @@
   Nsec = GetPhotcodeNsecfilt ();
 
-  /* load data about plot windows */
-  Ngraph = 0;
-  if (!GetGraphData (&graphsky, NULL, &Ngraph)) return (FALSE);
- 
   /* find CATDIR in config system */
-  VarConfig ("CATDIR", "%s", catdir);
-
-  /* check for CPT selection */
-  CPTfile = NULL;
-  if (N = get_argument (argc, argv, "-cpt")) {
-    remove_argument (N, &argc, argv);
-    CPTfile = strcreate (argv[N]);
-    remove_argument (N, &argc, argv);
-  }    
+  if (VarConfig ("CATDIR", "%s", catdir) == NULL) goto escape;
 
   /* interpret command-line options */
-  if (!SetPhotSelections (&argc, argv)) goto usage;
+  if (!SetRegionSelection (&argc, argv, &RegionName, &RegionList)) goto escape;
+  if (!SetPhotSelections (&argc, argv, 3)) goto usage;
+
+  /* interpret required command-line arguments: dmags F1 - F2 : F3 */
   if (argc != 6) { goto usage; }
   if (strcmp (argv[2], "-")) goto usage;
   if (strcmp (argv[4], ":")) goto usage;
-
-  /* interpret command-line options */
   if (!GetPhotcodeInfo (argv[1], &code[0], &mode[0])) goto usage;
   if (!GetPhotcodeInfo (argv[3], &code[1], &mode[1])) goto usage;
   if (!GetPhotcodeInfo (argv[5], &code[2], &mode[2])) goto usage;
+  if (!TestPhotSelections (&code[0], MEAS_ZERO)) goto escape;
 
-  /* test PhotSelections: is photcode specified if needed? */
-  for (i = 0; i < 3; i++) {
-    if (!TestPhotSelections (&code[i])) {
-      fprintf (stderr, "photcode selection rules violated\n");
-      fprintf (stderr, "average and ensemble restrictions require a corresponding photcode\n");
-      goto escape;
-    }
-  }
-
-  /* determine region-file names */
-  if (CPTfile == NULL) {
-    float Radius;
-    Radius = MAX (fabs(graphsky.xmax), fabs(graphsky.ymax));
-    regions = find_regions (graphsky.coords.crval1, graphsky.coords.crval2, Radius, &Nregions);
-  } else {
-    Nregions = 1;
-    ALLOCATE (regions, RegionFile, 1);
-    strcpy (regions[0].name, CPTfile);
-  }
+  /* load region corresponding to selection above */
+  if ((regions = SelectRegions (RegionName, RegionList, &Nregions)) == NULL) goto escape;
 
   /* init vectors to save data */
@@ -99 +73 @@
       m = catalog.average[i].offset;
 
+      SetSelectionParam (0);
       M1 = ExtractDMag (code, mode, &catalog.average[i], &catalog.secfilt[i*Nsec], &catalog.measure[m], &N1);
       if (N1 == 0) goto skip;
 
+      SetSelectionParam (2);
       M3 = ExtractMagnitudes (code[2], mode[2], &catalog.average[i], &catalog.secfilt[i*Nsec], &catalog.measure[m], &N3);
       if (N3 == 0) goto skip;
@@ -110 +86 @@
           yvec[0].elements[Npts] = M3[i3];
           Npts++;
-          if (Npts == NPTS) {
+          if (Npts >= NPTS) {
             NPTS += 2000;
             REALLOCATE (xvec[0].elements, float, NPTS);
@@ -134 +110 @@
 usage:
   fprintf (stderr, "USAGE: dmags F - F : F\n");
-  fprintf (stderr, "  F : any photcodes with matched qualifiers:\n");
+  fprintf (stderr, "    F : any photcodes with matched qualifiers:\n");
   fprintf (stderr, "    pri: F:inst, F:cat, F:sys, F:rel, F:cal, F:ave, F:ref\n");
   fprintf (stderr, "    sec: F:inst, F:cat, F:sys, F:rel, F:cal, F:ave, F:ref\n");
   fprintf (stderr, "    dep: F:inst, F:cat, F:sys, F:rel, F:cal\n");
   fprintf (stderr, "    ref: F:cat\n");
-  return (FALSE);
 
 escape:
@@ -146 +121 @@
   if (catalog.secfilt != NULL) free (catalog.secfilt);
   if (catalog.measure != NULL) free (catalog.measure);
+  if (RegionName != NULL) free (RegionName);
+  if (RegionList != NULL) free (RegionList);
   return (FALSE);
 }

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

@@ -129 +129 @@
 
   mode = S_IRUSR | S_IWUSR | S_IRGRP | S_IWGRP | S_IROTH | S_IWOTH;
-  chmod (f, mode);
+  chmod (fifo, mode);
   fclearlockfile (fifo, f, LCK_XCLD, &state);
 

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

@@ -13 +13 @@
   int i, k, m, NP1, NP2, NP, Np, Npts, NPTS;
   int N1, N2, i1, i2, mode[4];
-  int Nsec, Ngraph, Nregions, Ncatalog, status;
-  char catdir[256], filename[256];
+  int Nsec, Nregions, Ncatalog, status;
+  char catdir[256], filename[256], *RegionName, *RegionList;
   char *cmd, *outcmd, *camera;
   double Radius;
   double *M1, *M2;
   float *out;
+
   Catalog *catalog;
-  Graphdata graphsky;
   RegionFile *regions;
   PhotCode **codelist, *tcode, *code[4];
@@ -32 +32 @@
   xvec = yvec = NULL;
 
-  /* find CATDIR in config system */
-  VarConfig ("CATDIR", "%s", catdir);
-
   /* load photcode information */
   if (!InitPhotcodes ()) goto escape;
   Nsec = GetPhotcodeNsecfilt ();
 
-  /*** determine appropriate catalog files (find_regions) ***/
-  /* load data about plot windows */
-  Ngraph = 0;
-  if (!GetGraphData (&graphsky, NULL, &Ngraph)) return (FALSE);
-  /* find catalog files which overlap this region */
-  Radius = MAX (fabs(graphsky.xmax), fabs(graphsky.ymax));
-  regions = find_regions (graphsky.coords.crval1, graphsky.coords.crval2, Radius, &Nregions);
+  /* find CATDIR in config system */
+  if (VarConfig ("CATDIR", "%s", catdir) == NULL) goto escape;
 
   /* interpret command-line options */
-  /*** determine any restrictions (time, ?) ***/
-  if (!SetPhotSelections (&argc, argv)) goto usage;
+  if (!SetRegionSelection (&argc, argv, &RegionName, &RegionList)) goto escape;
+  if (!SetPhotSelections (&argc, argv, 1)) goto usage;
+
+  /* interpret command-line options */
   if (argc != 4) goto usage;
-
-  /* vectors to save data */
-  Npts = 0;
-  NPTS = 1;
-  if ((xvec = SelectVector ("tmp_x", ANYVECTOR, TRUE)) == NULL) goto escape;
-  if ((yvec = SelectVector ("tmp_y", ANYVECTOR, TRUE)) == NULL) goto escape;
 
   /* determine relevant photcodes, colors */
@@ -108 +96 @@
     out[i] = -1;
   }
+
+  /* load region corresponding to selection above */
+  if ((regions = SelectRegions (RegionName, RegionList, &Nregions)) == NULL) goto escape;
 
   /* loop over regions, extract data for each region */
@@ -131 +122 @@
   fprintf (stderr, "using %d regions\n", Nregions);
 
+  /* vectors to save data */
+  Npts = 0;
+  NPTS = 1;
+  if ((xvec = SelectVector ("tmp_x", ANYVECTOR, TRUE)) == NULL) goto escape;
+  if ((yvec = SelectVector ("tmp_y", ANYVECTOR, TRUE)) == NULL) goto escape;
+
   /*** generate the color-color vectors for the pairs ***/
   /* loop over chip photcode pairs */
@@ -159 +156 @@
               xvec[0].elements[Npts] = M2[i2];
               Npts++;
-              if (Npts == NPTS) {
+              if (Npts >= NPTS) {
                 NPTS += 2000;
                 REALLOCATE (xvec[0].elements, float, NPTS);

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

@@ -135 +135 @@
         vec[0].elements[N] = catalog.average[n].R - catalog.measure[i].dR / 360000.0;
         N++;
-        if (N == NPTS - 1) {
-          NPTS += 1000;
-          REALLOCATE (vec[0].elements, float, NPTS);
-        }
+        CHECK_REALLOCATE (vec[0].elements, float, NPTS, N, 1000);
       }
       break;
@@ -147 +144 @@
         vec[0].elements[N] = catalog.average[n].D - catalog.measure[i].dD / 360000.0;
         N++;
-        if (N == NPTS - 1) {
-          NPTS += 1000;
-          REALLOCATE (vec[0].elements, float, NPTS);
-        }
+        CHECK_REALLOCATE (vec[0].elements, float, NPTS, N, 1000);
       }
       break;
@@ -158 +152 @@
         vec[0].elements[N] = catalog.measure[i].M / 1000.0;
         N++;
-        if (N == NPTS - 1) {
-          NPTS += 1000;
-          REALLOCATE (vec[0].elements, float, NPTS);
-        }
+        CHECK_REALLOCATE (vec[0].elements, float, NPTS, N, 1000);
       }
       break;
@@ -169 +160 @@
         vec[0].elements[N] = catalog.measure[i].dM / 1000.0;
         N++;
-        if (N == NPTS - 1) {
-          NPTS += 1000;
-          REALLOCATE (vec[0].elements, float, NPTS);
-        }
+        CHECK_REALLOCATE (vec[0].elements, float, NPTS, N, 1000);
       }
       break;
@@ -180 +168 @@
         vec[0].elements[N] = catalog.measure[i].Mcal / 1000.0;
         N++;
-        if (N == NPTS - 1) {
-          NPTS += 1000;
-          REALLOCATE (vec[0].elements, float, NPTS);
-        }
+        CHECK_REALLOCATE (vec[0].elements, float, NPTS, N, 1000);
       }
       break;
@@ -192 +177 @@
         vec[0].elements[N] = catalog.average[n].M / 1000.0;
         N++;
-        if (N == NPTS - 1) {
-          NPTS += 1000;
-          REALLOCATE (vec[0].elements, float, NPTS);
-        }
       }
       break;
@@ -203 +184 @@
         vec[0].elements[N] = catalog.measure[i].source;
         N++;
-        if (N == NPTS - 1) {
-          NPTS += 1000;
-          REALLOCATE (vec[0].elements, float, NPTS);
-        }
+        CHECK_REALLOCATE (vec[0].elements, float, NPTS, N, 1000);
       }
       break;
@@ -214 +192 @@
         vec[0].elements[N] = TimeValue (catalog.measure[i].t, TimeReference, TimeFormat);
         N++;
-        if (N == NPTS - 1) {
-          NPTS += 1000;
-          REALLOCATE (vec[0].elements, float, NPTS);
-        }
+        CHECK_REALLOCATE (vec[0].elements, float, NPTS, N, 1000);
       }
       break;

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

@@ -4 +4 @@
   
   int i, j, k, m, N, N1, NPTS;
-  int param, mode, Nregions, Ngraph, Nsec;
-  char filename[128], catdir[256], *p;
-  double *M1;
+  int param, mode, Nregions, Nsec;
+  char filename[128], catdir[256], *RegionName, *RegionList, *p;
+  double Radius, *M1;
 
   PhotCode *code;
   Catalog catalog;
   RegionFile *regions;
-  Graphdata graphsky;
   Vector *vec;
 
+  /* defaults */
   regions = NULL;
   catalog.average = (Average *) NULL; 
   catalog.secfilt = (SecFilt *) NULL;
   catalog.measure = (Measure *) NULL;
+  RegionName = NULL;
+  RegionList = NULL;
+  code = NULL;
+  mode = MAG_REL;
 
   /* load photcode information */
@@ -23 +27 @@
   Nsec = GetPhotcodeNsecfilt ();
 
-  /* load data about plot windows */
-  Ngraph = 0;
-  if (!GetGraphData (&graphsky, NULL, &Ngraph)) goto escape;
- 
   /* find CATDIR in config system */
-  VarConfig ("CATDIR", "%s", catdir);
-
-  /* specify catalog files
-  RegionList = FALSE;
-  if (N = get_argument (argc, argv, "-list")) {
-    RegionList = TRUE;
-    remove_argument (N, &argc, argv);
-    strcpy (RegionListFile, argv[N]);
-    remove_argument (N, &argc, argv);
-  } */
+  if (VarConfig ("CATDIR", "%s", catdir) == NULL) goto escape;
 
   /* interpret command-line options */
-  if (!SetPhotSelections (&argc, argv)) goto usage;
-  if (argc != 3) goto usage;
+  SetSelectionParam (0);
+  if (!SetRegionSelection (&argc, argv, &RegionName, &RegionList)) goto escape;
+  if (!SetPhotSelections (&argc, argv, 1)) goto usage;
 
-  /* identify selection */
-  code = NULL;
-  mode = MAG_REL;
-  param = GetMeasureParam (argv[2]);
+  /* interpret required command-line arguments: mextract (value) */
+  if (argc != 2) goto usage;
+  param = GetMeasureParam (argv[1]);
   if (param == MEAS_ZERO) {
-    if (!GetPhotcodeInfo (argv[2], &code, &mode)) {
+    if (!GetPhotcodeInfo (argv[1], &code, &mode)) {
       GetMeasureParam ("help");
-      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");
       goto escape;
     }
     param = MEAS_MAG;
-    for (p = argv[2]; *p != 0; p++) { 
+    for (p = argv[1]; *p != 0; p++) { 
       if (*p == '.') *p = ':';
     }
-  } 
-  /* test PhotSelections: is photcode specified if needed? */
-  if (!TestPhotSelections (&code)) {
-    fprintf (stderr, "photcode selection rules violated\n");
-    fprintf (stderr, "average and ensemble restrictions require a corresponding photcode\n");
-    goto escape;
   }
+  if (!TestPhotSelections (&code, MEAS_ZERO)) goto escape;
 
-  if ((vec = SelectVector (argv[2], ANYVECTOR, TRUE)) == NULL) goto escape;
-
-  /* determine region-file names */
-  if (!strcmp (argv[1], "all")) {
-    float Radius;
-    Radius = MAX (fabs(graphsky.xmax), fabs(graphsky.ymax));
-    regions = find_regions (graphsky.coords.crval1, graphsky.coords.crval2, Radius, &Nregions);
-    /* I don't have a tool to select images by the overlapping region */
-    if (!SetImageSelection (param, &graphsky, TRUE)) goto escape;
-  } else {
-    Nregions = 1;
-    ALLOCATE (regions, RegionFile, 1);
-    strcpy (regions[0].name, argv[1]);
-    if (!SetImageSelection (param, NULL, FALSE)) goto escape;
-  }
+  /* load region corresponding to selection above */
+  if ((regions = SelectRegions (RegionName, RegionList, &Nregions)) == NULL) goto escape;
+  if (!SetImageSelection (param)) goto escape;
 
   /* create storage vector */
-  NPTS = 1000;
-  REALLOCATE (vec[0].elements, float, NPTS);
-  vec[0].Nelements = N = 0;
+  N = 0;
+  NPTS = 1;
+  if ((vec = SelectVector (argv[1], ANYVECTOR, TRUE)) == NULL) goto escape;
 
   for (i = 0; i < Nregions; i++) {
@@ -110 +82 @@
         vec[0].elements[N] = M1[k];
         N++;
-        if (N == NPTS - 1) {
-          NPTS += 2000;
-          REALLOCATE (vec[0].elements, float, NPTS);
-        }
+        CHECK_REALLOCATE (vec[0].elements, float, NPTS, N, 2000);
       }
       if (M1 != NULL) free (M1);
@@ -132 +101 @@
 
 usage:
-  fprintf (stderr, "USAGE: mextract (from) (value) [options]\n");
-  fprintf (stderr, "  from: cpt name or 'all'\n");
+  fprintf (stderr, "USAGE: mextract (value) [options]\n");
   fprintf (stderr, "  value: measure.parameter or photcode\n");
   return (FALSE);
 
 escape:
+  FreeImageSelection ();
   if (regions != NULL) free (regions);
-  FreeImageSelection ();
   if (catalog.average != NULL) free (catalog.average);
   if (catalog.secfilt != NULL) free (catalog.secfilt);
   if (catalog.measure != NULL) free (catalog.measure);
+  if (RegionName != NULL) free (RegionName);
+  if (RegionList != NULL) free (RegionList);
   return (FALSE);
 }

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

@@ -3 +3 @@
 /* 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 +37 @@
   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 +56 @@
 static double TypefracValue;
 
+/* apply selections or not */
+static int ApplySelections[4];
+static int SelectionParam;
+
 /* applied to Average quantities */
 static int PhotcodeSelect;
@@ -69 +78 @@
 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 +164 @@
 }
  
+int SetSelectionParam (int param) {
+  SelectionParam = param;
+  return (TRUE);
+}
+
+int GetSelectionParam () {
+  return (SelectionParam);
+}
+
 int GetMagMode (char *string) {
 
@@ -196 +213 @@
     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 +227 @@
   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 +246 @@
 }
 
+/* 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 +295 @@
 }
 
-/* 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 +314 @@
 /* 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 +363 @@
   }
 
-  /* 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 +389 @@
     if (PhotcodeValue == NULL) {
       fprintf (stderr, "photcode not found in photcode table\n");
-      PhotcodeSelect = FALSE;
+      return (FALSE);;
     }
     remove_argument (N, argc, argv);
@@ -422 +427 @@
   }
 
+  /* 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 +522 @@
   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 +547 @@
   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 +642 @@
 }
 
-/* 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 +699 @@
 }
 
+/* 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 +729 @@
   }
 
+  /* 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 +863 @@
 }
 
+/* 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 +882 @@
   ALLOCATE (list, double, NLIST);
   M1 = M2 = NULL;
+  NoMag = NO_MAG * 0.001;
 
   /* one of the two is an average, must do independently */
@@ -960 +888 @@
   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 +923 @@
   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 +945 @@
     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 +961 @@
 }
 
+/* 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 +975 @@
   }
 
+  /* 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 +988 @@
 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 +1020 @@
 
       /* 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 +1028 @@
 }
 
+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 
+**/

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

@@ -1 +1 @@
 # include "dvo1.h"
+
+static int RegionSelect;
+
+int GetRegionSelection () {
+  return (RegionSelect);
+}
+
+int SetRegionSelection (int *argc, char **argv, char **RegionName, char **RegionList) {
+  
+  int N;
+
+  RegionSelect = FALSE;
+
+  /* check for Region selection */
+  *RegionName = NULL;
+  if (N = get_argument (*argc, argv, "-cpt")) {
+    remove_argument (N, argc, argv);
+    *RegionName = strcreate (argv[N]);
+    remove_argument (N, argc, argv);
+    return (TRUE);
+  }    
+  /* check for Region list */
+  *RegionList = NULL;
+  if (N = get_argument (*argc, argv, "-cptlist")) {
+    remove_argument (N, argc, argv);
+    *RegionList = strcreate (argv[N]);
+    remove_argument (N, argc, argv);
+    return (TRUE);
+  } 
+  if ((*RegionName == NULL) && (*RegionList == NULL)) {
+    RegionSelect = TRUE;
+    return (TRUE);
+  }
+
+  fprintf (stderr, "-cpt and -cptlist are incompatible\n");
+  free (*RegionName);
+  free (*RegionList);
+  *RegionName = NULL;
+  *RegionList = NULL;
+  return (FALSE);
+}
+
+RegionFile *SelectRegions (char *RegionName, char *RegionList, int *nregions) {
+
+  double Radius;
+  int Ngraph;
+  Graphdata graphsky;
+  RegionFile *regions;
+
+  /* determine region-file names */
+  if (RegionName != NULL) {
+    *nregions = 1;
+    ALLOCATE (regions, RegionFile, 1);
+    strcpy (regions[0].name, RegionName);
+    free (RegionName);
+    RegionName = NULL;
+    return (regions);
+  } 
+  if (RegionList != NULL) {
+    regions = region_list (RegionList, nregions);
+    free (RegionList);
+    RegionList = NULL;
+    return (regions);
+  }
+
+  Ngraph = 0;
+  if (!GetGraphData (&graphsky, NULL, &Ngraph)) {
+    fprintf (stderr, "region display not available\n");
+    return (NULL);
+  }
+
+  Radius = MAX (fabs(graphsky.xmax), fabs(graphsky.ymax));
+  regions = find_regions (graphsky.coords.crval1, graphsky.coords.crval2, Radius, nregions);
+  return (regions);
+}
 
 /* returns a list of region files names from file */
@@ -12 +87 @@
     fprintf (stderr, "ERROR: can't find region list file %s\n", filename);
     *Nregions = 0;
-    return ((RegionFile *) NULL);
+    return (NULL);
   }
   
@@ -28 +103 @@
 
   *Nregions = nregion;
-
   return (regions);
-
 }
 

trunk/Ohana/src/opihi/include/dvo1.h

@@ -40 +40 @@
 
 /*** photometry functions ***/
-int           SetPhotSelections     PROTO((int *argc, char **argv));
-int           SetImageSelection     PROTO((int mode, Graphdata *graphsky, int state));
-void          FreeImageSelection    PROTO(());
+double       *ExtractByDMag         PROTO((PhotCode **code, int *mode, Average *average, SecFilt *secfilt, Measure *measure, int *nlist, int param));
+double       *ExtractDMag           PROTO((PhotCode **code, int *mode, Average *average, SecFilt *secfilt, Measure *measure, int *nlist));
+double       *ExtractMagnitudes     PROTO((PhotCode *code, int mode, Average *average, SecFilt *secfilt, Measure *measure, int *n));
+double       *ExtractMeasures       PROTO((PhotCode *code, int mode, Average *average, SecFilt *secfilt, Measure *measure, int *nlist, int param));
+double       *ExtractMeasuresByDMag PROTO((PhotCode **code, int *mode, int use_first, Average *average, SecFilt *secfilt, Measure *measure, int *nlist, int param));
+double       *ExtractMeasuresDMag   PROTO((PhotCode **code, int *mode, Average *average, SecFilt *secfilt, Measure *measure, int *nlist));
+Image        *MatchImage            PROTO((unsigned int time, short int source));
+Coords       *MatchMosaic           PROTO((unsigned int time, short int source));
+int           DetermineTypeCode     PROTO((Average *average, Measure *measure, int code));
+double        DetermineTypefrac     PROTO((Average *average, Measure *measure, PhotCode *code));
+double        ExtractAverages       PROTO((PhotCode *code, int mode, Average *average, SecFilt *secfilt, Measure *measure, int param));
+void          FreeImageSelection    PROTO(());
+int           GetAverageParam       PROTO((char *parname));
+void          GetAverageParamHelp   PROTO(());
+int           GetMagMode            PROTO((char *string));
+double        GetMeasure            PROTO((int param, Average *average, Measure *measure, double mag));
+int           GetMeasureParam       PROTO((char *parname));
+int           GetMeasureTypeCode    PROTO((Measure *measure));
+int           GetPhotcodeInfo       PROTO((char *string, PhotCode **Code, int *Mode));
+int           GetSelectionParam     PROTO(());
+int           InitPhotcodes         PROTO(());
+int           Quality               PROTO((Measure *measure, int IsDophot));
+int           SelectMags            PROTO((int Nphot, int Tphot, int Ns, Average *average, Measure *measure, SecFilt *secfilt, int UL));
+int           SetImageSelection     PROTO((int mode));
+int           SetPhotSelections     PROTO((int *argc, char **argv, int Nparams));
+int           SetSelectionParam     PROTO((int param));
+int           TestAverage           PROTO((PhotCode *code, Average *average, SecFilt *secfilt, Measure *measure));
+int           TestPhotSelections    PROTO((PhotCode **code, int param));
+void          print_value           PROTO((FILE *f, double value, short int ival));
 
-int           InitPhotcodes         PROTO(());
-int           GetPhotcodeInfo       PROTO((char *string, PhotCode **Code, int *Mode));
-
-int           GetAverageParam       PROTO((char *parname));
-int           GetMeasureParam       PROTO((char *parname));
-int           GetMagMode            PROTO((char *string));
-int           SelectMags            PROTO((int Nphot, int Tphot, int Ns, Average *average, Measure *measure, SecFilt *secfilt, int UL));
-void          print_value           PROTO((FILE *f, double value, short int ival));
-int           GetMeasureTypeCode    PROTO((Measure *measure));
-int           DetermineTypeCode     PROTO((Average *average, Measure *measure, int Nphot));
-double        ExtractAverages       PROTO((PhotCode *code, int mode, Average *ave, SecFilt *sec, Measure *meas, int param));
-double       *ExtractMagnitudes     PROTO((PhotCode *code, int mode, Average *ave, SecFilt *sec, Measure *meas, int *n));
-double       *ExtractMeasures       PROTO((PhotCode *code, int mode, Average *ave, SecFilt *sec, Measure *meas, int *nlist, int param));
-double       *ExtractDMag           PROTO((PhotCode **code, int *mode, Average *average, SecFilt *secfilt, Measure *measure, int *nlist));
-double       *ExtractMeasuresDMag   PROTO((PhotCode **code, int *mode, Average *average, SecFilt *secfilt, Measure *measure, int *nlist));
-double       *ExtractMeasuresByDMag PROTO((PhotCode **code, int *mode, int use_first, Average *average, SecFilt *secfilt, Measure *measure, int *nlist, int param));
-double        DetermineTypefrac     PROTO((Average *average, Measure *measure, PhotCode *code));
-void          GetAverageParamHelp   PROTO(());
-
+RegionFile   *SelectRegions         PROTO((char *RegionName, char *RegionList, int *nregions));
 
 # endif