Changeset 30509

Timestamp:

Feb 7, 2011, 4:45:12 PM (15 years ago)

Author:

eugene

Message:

some improvements to the visualization; add user config values for IMFIT_CLIP_NITER, IMFIT_CLIP_NSIGMA, IMFIT_SYS_SIGMA_LIM; calcualte dXpixSys and dYpixSys for images and dRsys (in arcsec) for measures; better tracking of why images are masked; better matching of images to mosaics; do not reset magnitude calibrations in relastro; use the measurement position errors it the fits (use error model from photcode table)

Location:

branches/eam_branches/ipp-20101205/Ohana/src/relastro

Files:

• include/relastro.h (modified) (3 diffs)
• src/ConfigInit.c (modified) (1 diff)
• src/FitChip.c (modified) (6 diffs)
• src/GetAstromError.c (modified) (1 diff)
• src/ImageOps.c (modified) (8 diffs)
• src/UpdateChips.c (modified) (5 diffs)
• src/UpdateObjects.c (modified) (1 diff)
• src/args.c (modified) (1 diff)
• src/relastro.c (modified) (1 diff)
• src/relastroVisual.c (modified) (20 diffs)
• src/select_images.c (modified) (4 diffs)

branches/eam_branches/ipp-20101205/Ohana/src/relastro/include/relastro.h

@@ -94 +94 @@
 
 double SIGMA_LIM;
-int SRC_MEAS_TOOFEW; //catalog objects wich fewer detections then this are ignored
+int    SRC_MEAS_TOOFEW; //catalog objects wich fewer detections then this are ignored
 double MIN_ERROR;
+
+int    IMFIT_CLIP_NITER; // number of clipping iterations to perform in FitChip
+double IMFIT_CLIP_NSIGMA; // number of sigma to clip in FitChip
+double IMFIT_SYS_SIGMA_LIM; // max dMag for objects used to measure systematic scatter
 
 double RADIUS; // match radius for high-speed objects
@@ -325 +329 @@
 int UpdateObjectOffsets (SkyList *skylist);
 
-int relastroVisualPlotRawRef(StarData *raw, StarData *ref, double dRmax, int numObj);
-int relastroVisualPlotScatter(double values[], double thresh, int npts);
-int relastroVisualPlotOutliers(Catalog *catalog, int offset, int Nmeasure, 
-                               StatType statsR, StatType statsD, double thresh);
-
-
+int relastroVisualPlotChipFit(StarData *raw, StarData *ref, double dRmax, int numObj);
+// int relastroVisualPlotRawRef(StarData *raw, StarData *ref, double dRmax, int numObj);
+// int relastroVisualPlotScatter(double values[], double thresh, int npts);
+// int relastroVisualPlotOutliers(Catalog *catalog, int offset, int Nmeasure, StatType statsR, StatType statsD, double thresh);
+void relastroSetVisual(int state);
 
 int FixProblemImages (SkyList *skylist);
@@ -352 +355 @@
 int createStarMapPoints();
 int checkStarMap(int N);
+
+int GetScatterRawRef(float *dLsig, float *dMsig, float *dRsig, int *nKeep, StarData *raw, StarData *ref, int Nstars, float SigmaLimit);

branches/eam_branches/ipp-20101205/Ohana/src/relastro/src/ConfigInit.c

@@ -19 +19 @@
   if (VERBOSE) fprintf (stderr, "loaded config file: %s\n", file);
 
-  GetConfig (config, "RELASTRO_SIGMA_LIM",     "%lf", 0, &SIGMA_LIM);
+  GetConfig (config, "RELASTRO_SIGMA_LIM",         "%lf", 0, &SIGMA_LIM); // exclude measurements on this basis
   GetConfig (config, "RELASTRO_SRC_MEAS_TOOFEW",   "%d",  0, &SRC_MEAS_TOOFEW);
+
+  if (!ScanConfig (config, "RELASTRO_IMFIT_CLIP_NITER",    "%d",  0, &IMFIT_CLIP_NITER))    IMFIT_CLIP_NITER  = 3;
+  if (!ScanConfig (config, "RELASTRO_IMFIT_CLIP_NSIGMA",   "%lf", 0, &IMFIT_CLIP_NSIGMA))   IMFIT_CLIP_NSIGMA = 3.0;
+  if (!ScanConfig (config, "RELASTRO_IMFIT_SYS_SIGMA_LIM", "%lf", 0, &IMFIT_SYS_SIGMA_LIM)) IMFIT_SYS_SIGMA_LIM = 0.01;
 
   // XXX these are used in relphot to identify poor stars / images -- define

branches/eam_branches/ipp-20101205/Ohana/src/relastro/src/FitChip.c

@@ -2 +2 @@
 # include "relastroVisual.h"
 
-// XXX make these user parameters
-# define FIT_CHIP_NITER     3
-# define FIT_CHIP_NSIGMA    3.0
-
-// XXX save the fit[0].Npts value in the image table?
-
-// XXX save measurements of the fit quality (scatter, chisq) in the image table
-
 int FitChip (StarData *raw, StarData *ref, int Nmatch, Image *image) {
 
-  int i, Nscatter, Niter, skip;
-  CoordFit *fit;
-  double dL, dM, dR, dRmax, *values;
-
-  ALLOCATE (values, double, Nmatch);
-  for (Niter = 0; Niter < FIT_CHIP_NITER; Niter ++) {
+  int i, NstatFull, Nstat, Niter, skip;
+  float dLsig, dMsig, dRsig;
+  float dLsigFull, dMsigFull, dRsigFull;
+  float dL, dM, dR, dRmax;
+
+  CoordFit *fit = NULL;
+
+  for (Niter = 0; Niter < IMFIT_CLIP_NITER; Niter ++) {
+    
+    // measure the scatter for the unmarked (good) measurements with dM < limit
+    // SIGMA_LIM here is redundant with ImageOps.c:915
+    GetScatterRawRef(&dLsig, &dMsig, &dRsig, &Nstat, raw, ref, Nmatch, SIGMA_LIM);
+    dRmax = IMFIT_CLIP_NSIGMA*dRsig;
+
+    // (a) skip unfittable points
+    // (b) mark good and bad points for fit (in and outliers)
 
     // measure the scatter distribution (use only the bright end detections)
-    for (i = Nscatter = 0; i < Nmatch; i++) {
+    for (i = 0; i < Nmatch; i++) {
       if (raw[i].mask) continue;
-      if (isnan(raw[i].dMag) || raw[i].dMag > SIGMA_LIM) continue;
+      if (isnan(raw[i].dMag)) { 
+        raw[i].mask |= 0x0004; 
+        continue; 
+      }
+      if (raw[i].dMag > SIGMA_LIM) { 
+        raw[i].mask |= 0x0008; 
+        continue; 
+      } // is this redundant with ImageOps.c:915?
 
       dL = raw[i].L - ref[i].L;
       dM = raw[i].M - ref[i].M;
       dR = hypot (dL, dM);
-
-      values[Nscatter] = dR;
-      Nscatter++;
-    }
-
-    if (Nscatter > 5) {
-      // for a 2D Gaussian, 40% of the points are within 1 sigma; dRmax is ~ 3 sigma
-      // XXX this test is not sensible for Nscatter < XXX (5?)
-      dsort (values, Nscatter);
-      dRmax = FIT_CHIP_NSIGMA*values[(int)(0.40*Nscatter)];
-      relastroVisualPlotScatter(values, dRmax, Nscatter);
-      relastroVisualPlotRawRef(raw, ref, dRmax, Nmatch);
-    } else {
-      dRmax = 0;
-    }
+      if (dR > dRmax) {
+        raw[i].mask |= 0x0010;
+        continue;
+      }
+    }
+
+    // figures to assess the fitting process:
+    // x vs dx, x vs dy, y vs dx, y vs dy : 
+    // residual vector field 
+    // dx vs mag, dy vs mag
+    relastroVisualPlotChipFit(raw, ref, dRmax, Nmatch);
 
     // fit the requested order polynomial
@@ -47 +52 @@
       image[0].coords.Npolyterms = CHIPORDER;
     }
+    if (fit) fit_free (fit);
     fit = fit_init (image[0].coords.Npolyterms);
 
     // generate the fit matches
     for (i = 0; i < Nmatch; i++) {
-      if (raw[i].mask) {
-        continue;
-      }
-      if (isnan(raw[i].dMag) || raw[i].dMag > SIGMA_LIM) {
-        continue;
-      }
-
-      // only keep objects within dRmax
-      dL = raw[i].L - ref[i].L;
-      dM = raw[i].M - ref[i].M;
-      dR = hypot (dL, dM);
-
-      // fprintf (stderr, "fit %f %f -> %f %f : %f %f (%f vs %f) wt: %f\n", raw[i].X, raw[i].Y, ref[i].L, ref[i].M, raw[i].L, raw[i].M, dR, dRmax, raw[i].dPos);
-
-      if ((dRmax > 0.0) && (dR > dRmax)) continue;
-
+      if (raw[i].mask) continue;
       fit_add (fit, raw[i].X, raw[i].Y, ref[i].L, ref[i].M, raw[i].dPos);
     }
@@ -90 +81 @@
       if (VERBOSE2) fprintf (stderr, "insufficient measurements (%d) for requested order (%d)\n", fit[0].Npts, image[0].coords.Npolyterms);
       fit_free (fit);
-      free (values);
       image[0].flags |= ID_IMAGE_ASTROM_FEW;
       return FALSE;
     }
-
-    // fprintf (stderr, "scatter limit: %f based on %d detections; using %d of %d for fit\n", dRmax, Nscatter, fit[0].Npts, Nmatch);
 
     // measure the fit, update the coords & object coordinates
@@ -110 +98 @@
     }
 
-    fit_free (fit);
-
     for (i = 0; i < Nmatch; i++) {
       XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, &image[0].coords);
@@ -118 +104 @@
   }
 
-  free (values);
+  // count mask classes
+  int nMask1 = 0;
+  int nMask2 = 0;
+  int nMask3 = 0;
+  int nMask4 = 0;
+  int nMask5 = 0;
+  for (i = 0; i < Nmatch; i++) {
+    if (!raw[i].mask) continue;
+    if (raw[i].mask & 0x01) nMask1 ++;
+    if (raw[i].mask & 0x02) nMask2 ++;
+    if (raw[i].mask & 0x04) nMask3 ++;
+    if (raw[i].mask & 0x08) nMask4 ++;
+    if (raw[i].mask & 0x10) nMask5 ++;
+  }
+
+  GetScatterRawRef(&dLsigFull, &dMsigFull, &dRsigFull, &NstatFull, raw, ref, Nmatch, SIGMA_LIM);
+  GetScatterRawRef(&dLsig, &dMsig, &dRsig, &Nstat, raw, ref, Nmatch, IMFIT_SYS_SIGMA_LIM);
+  if (VERBOSE) fprintf (stderr, "fit sigma: %f (%f, %f) : full: %f (%f, %f), scatter limit: %f (%d full, %d bright, %d fit, %d all) (%d %d %d %d %d)\n", dRsig, dLsig, dMsig, dRsigFull, dLsigFull, dMsigFull, dRmax, NstatFull, Nstat, fit[0].Npts, Nmatch, nMask1, nMask2, nMask3, nMask4, nMask5);
+
+  image[0].dXpixSys = dLsig;
+  image[0].dYpixSys = dMsig;
+  image[0].nFitAstrom = fit[0].Npts;
+
+  if (fit) fit_free (fit);
+
   return TRUE;
+}
+
+// measure the scatter distribution (limit selected objects by dMag)
+int GetScatterRawRef(float *dLsig, float *dMsig, float *dRsig, int *nKeep, StarData *raw, StarData *ref, int Nstars, float SigmaLimit) {
+
+  int i, Ns;
+  float dL, dM, dLsum, dLsum2, dMsum, dMsum2, *dR;
+
+  ALLOCATE (dR, float, Nstars);
+
+  Ns = 0;
+  dLsum = dLsum2 = dMsum = dMsum2 = 0.0;
+  for (i = 0; i < Nstars; i++) {
+    if (raw[i].mask) continue;
+    if (isnan(raw[i].dMag)) continue;
+    if ((SigmaLimit > 0.0) && (raw[i].dMag > SigmaLimit)) continue;
+    
+    dL = raw[i].L - ref[i].L;
+    dM = raw[i].M - ref[i].M;
+
+    dLsum  += dL;
+    dLsum2 += dL*dL;
+    dMsum  += dM;
+    dMsum2 += dM*dM;
+
+    dR[Ns] = hypot (dL, dM);
+    Ns++;
+  }
+
+  *dLsig = sqrt((dLsum2 / Ns) - SQ(dLsum/Ns));
+  *dMsig = sqrt((dMsum2 / Ns) - SQ(dMsum/Ns));
+  *nKeep = Ns;
+
+  if (Ns < 5) {
+    *dRsig = NAN;
+    free  (dR);
+    return (FALSE);
+  }
+
+  // for a 2D Gaussian, 40% of the points are within R = 1 sigma
+  fsort (dR, Ns);
+  *dRsig = dR[(int)(0.40*Ns)];
+  
+  free  (dR);
+  return (TRUE);
 }
 
@@ -162 +217 @@
   fclose (f);
 */
+

branches/eam_branches/ipp-20101205/Ohana/src/relastro/src/GetAstromError.c

@@ -1 +1 @@
 # include "relastro.h"
-# define WEIGHTED_ERRORS 0
+# define WEIGHTED_ERRORS 1
 
 float GetAstromError (Measure *measure, int mode) {

branches/eam_branches/ipp-20101205/Ohana/src/relastro/src/ImageOps.c

@@ -336 +336 @@
   off_t i, m, c, n, nPos;
   double X, Y, L, M, P, Q, R, D, dR, dD;
-  double dPos, DPOS_MAX_ASEC;
+  double dPos, dPosSys, DPOS_MAX_ASEC;
 
   Mosaic *mosaic;
@@ -359 +359 @@
   }
 
-  // accumulate the rms position offsets.  if this value, or any specific entry, is too
-  // large, we will reset the image to the original coords at the end of the analysis
-
+  // these are used to accumulate the rms position offsets.  if this value, or any
+  // specific entry, is too large, we will reset the image to the original coords at the
+  // end of the analysis
   dPos = 0.0;
   nPos = 0;
+
+  // convert the image systematic error in pixels to a value in arcsec
+  { 
+    double dLsig, dMsig;
+    double Ro, Do, Rx, Dx, dP0, dP1;
+    Coords *coords;
+
+    // these values are in pixels, but we to convert to arcsec
+    dLsig = image[0].dXpixSys;
+    dMsig = image[0].dYpixSys;
+
+    if (moscoords == NULL) {
+      coords = imcoords;
+    } else {
+      coords = moscoords;
+    }
+    XY_to_LM (&Ro, &Do, 0.0, 0.0, coords);
+    XY_to_LM (&Rx, &Dx, dLsig, 0.0, coords);
+    dP0 = 3600.0 * hypot(Rx - Ro, Dx - Do); // convert to arcsec
+    XY_to_LM (&Rx, &Dx, 0.0, dLsig, coords);
+    dP1 = 3600.0 * hypot(Rx - Ro, Dx - Do); // convert to arcsec
+    dPosSys = 0.5 * (dP0 + dP1);
+  }      
 
   for (i = 0; i < Nlist[im]; i++) {
@@ -418 +441 @@
     catalog[c].measure[m].dR = dR;
     catalog[c].measure[m].dD = dD;
-
+    
     if (catalog[c].measure[m].dR > +180.0*3600.0) {
       // average on high end of boundary, move star up
@@ -429 +452 @@
       catalog[c].measure[m].dR = 3600.0*(catalog[c].average[n].R - R);
     }
+
+    // set the systematic error for this image:
+    catalog[c].measure[m].dRsys = ToShortPixels(dPosSys);
   }
 
@@ -544 +570 @@
 
     // an object with only one detection provides no information about the image calibration
-    //XXX this is already taken care of in bcatalog
-    raw[i].mask = FALSE;
-    int mask = FALSE;
+    // XXX this is already taken care of in bcatalog
+    raw[i].mask = 0x0000;
     if (catalog[c].average[n].Nmeasure <= SRC_MEAS_TOOFEW) {
-      mask = TRUE;
+      raw[i].mask |= 0x0001;
     }
     if (!finite(catalog[c].measure[m].dR) || !finite(catalog[c].measure[m].dD)) {
-      mask = TRUE;
-    }
-
-    raw[i].mask = mask;
-
+      raw[i].mask |= 0x0002;
+    }
 
     switch (mode) {
@@ -725 +747 @@
     }
 
-    //examine results
-    relastroVisualPlotOutliers(catalog, catalog[0].average[j].measureOffset, 
-                               catalog[0].average[j].Nmeasure, 
-                               statsR, statsD, Ns);
+    // examine results
+    // relastroVisualPlotOutliers(catalog, catalog[0].average[j].measureOffset, catalog[0].average[j].Nmeasure, statsR, statsD, Ns);
   }
   
@@ -852 +872 @@
     }  //done rejecting outliers
 
-    //examine results
-    relastroVisualPlotOutliers(catalog, catalog[0].average[j].measureOffset, 
-                               catalog[0].average[j].Nmeasure, 
-                               statsR, statsD, Ns);
+    // examine results
+    // relastroVisualPlotOutliers(catalog, catalog[0].average[j].measureOffset, catalog[0].average[j].Nmeasure, statsR, statsD, Ns);
     
   } //done looping over objects
@@ -917 +935 @@
 
   /* select measurements by measurement error */
-  if ((SIGMA_LIM > 0) && (measure[0].dM > SIGMA_LIM)) return FALSE;
+  if ((SIGMA_LIM > 0) && (measure[0].dM > SIGMA_LIM)) {
+    return FALSE;
+  }
   
   /* select measurements by mag limit */

branches/eam_branches/ipp-20101205/Ohana/src/relastro/src/UpdateChips.c

@@ -8 +8 @@
 
   /* we can measure new image parameters for each non-mosaic chip independently */
-  off_t i, Nimage, Nraw, Nref;
+  off_t i, Nimage, Nraw, Nref, nFitAstr;
   Image *image;
   StarData *raw, *ref;
   Coords *oldCoords;
-
+  float dXpixSys, dYpixSys;
   double *Ro, *Do;
   char *mode;
@@ -19 +19 @@
 
   image = getimages (&Nimage);
-  // XXX BuildChipMatch (image, Nimage);
 
   // save fit results for summary plot
@@ -54 +53 @@
     assert (Nraw == Nref);
 
+    // save these in case of failure
     saveCoords (&image[i].coords, i);
+    dXpixSys = image[i].dXpixSys;
+    dYpixSys = image[i].dYpixSys;
+    nFitAstr = image[i].nFitAstrom;
 
     // FitChip does iterative, clipped fitting
@@ -60 +63 @@
     if (!FitChip (raw, ref, Nraw, &image[i])) {
       if (VERBOSE) fprintf (stderr, "reject fit for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT"\n", image[i].name,  i,  Nraw);
+
+      // restore status quo ante
       oldCoords = getCoords (i);
       memcpy (&image[i].coords, oldCoords, sizeof(Coords));
+      image[i].dXpixSys = dXpixSys;
+      image[i].dYpixSys = dYpixSys; 
+      image[i].nFitAstrom = nFitAstr;
+
       saveCenter (image, &Ro[i], &Do[i], i);
       mode[i] = 1;
@@ -72 +81 @@
     if (!checkStarMap (i)) {
       if (VERBOSE) fprintf (stderr, "fit diverges too much for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT"\n", image[i].name,  i,  Nraw);
+      // restore status quo ante
       oldCoords = getCoords (i);
       memcpy (&image[i].coords, oldCoords, sizeof(Coords));
+      image[i].dXpixSys = dXpixSys;
+      image[i].dYpixSys = dYpixSys; 
+      image[i].nFitAstrom = nFitAstr;
+
       saveCenter (image, &Ro[i], &Do[i], i);
       mode[i] = 2;

branches/eam_branches/ipp-20101205/Ohana/src/relastro/src/UpdateObjects.c

@@ -136 +136 @@
         dX[N] = GetAstromError (&catalog[i].measure[m], ERROR_MODE_RA);
         dY[N] = GetAstromError (&catalog[i].measure[m], ERROR_MODE_DEC);
+
+        // add systematic error in quadrature, if desired
+        // only do this after the fit has converged (or you will never improve the poor images)
+        // if (INCLUDE_SYS_ERR) {
+        // float dRsys = FromShortPixels(catalog[i].measure[m].dRsys);
+        // dX[N] = hypot(dX[N], dRsys);
+        // dY[N] = hypot(dY[N], dRsys);
+        // }
+
         // dX[N] = 0.1;
         // dY[N] = 0.1;

branches/eam_branches/ipp-20101205/Ohana/src/relastro/src/args.c

@@ -144 +144 @@
     VERBOSE = VERBOSE2 = TRUE;
     remove_argument (N, &argc, argv);
+  }
+
+  if ((N = get_argument (argc, argv, "-visual"))) {
+    remove_argument (N, &argc, argv);
+    relastroSetVisual(TRUE);
   }
 

branches/eam_branches/ipp-20101205/Ohana/src/relastro/src/relastro.c

@@ -78 +78 @@
         MARKTIME("update chips: %f sec\n", dtime);
       }
+      // create summary plots of the process
+      // relastroVisualSummaryChips();
       break;
 

branches/eam_branches/ipp-20101205/Ohana/src/relastro/src/relastroVisual.c

@@ -11 +11 @@
 #define KAPAY 700
 
-static int kapa = -1;
+static int kapa1 = -1;
 static int kapa2 = -1;
 static int kapa3 = -1;
+static int kapa4 = -1;
 
 static int isVisual = FALSE;
-static int plotRawRef = FALSE;
-static int plotScatter = FALSE;
-static int plotResid = FALSE;
-static int plotVector = FALSE;
+static int plotRawRef = TRUE;
+// static int plotScatter = TRUE;
+// static int plotResid = TRUE;
+// static int plotVector = TRUE;
 static int plotOutliers = TRUE;
 
@@ -29 +30 @@
       fprintf(stderr, "Failure to open kapa.\n");
       isVisual = 0;
-      return 0;
+      return FALSE;
     }
-    //    KapaResize (*kapid, KAPAX, KAPAY);
-  }
-  return 1;
+  }
+  return TRUE;
 }
 
@@ -49 +49 @@
     isVisual = 0;
   }
-  return 1;
-}
+  return TRUE;
+}
+
+void relastroSetVisual(int state) {
+  isVisual = state;
+}
+
 
 /** Size graphdata to encompass all points */
@@ -70 +75 @@
     graphdata->xmax = xhi;
     graphdata->ymax = yhi;
-    return 1;
-}
-
-static int residPlot(float x[], float y[],
-                     float xVec[], float yVec[],
-                     int npts, int *kapaID) {
-    if (!isVisual || !plotResid) return TRUE;
+    return TRUE;
+}
+
+/** 4-panel plot of x vs dx, y vs dx, x vs dy, y vs dy **/
+int relastroVisualRawRef(int *kapaID, float *x, float *y, float *dx, float *dy, float *dPos, int npts) {
+
     Graphdata graphdata;
     KapaSection section;
@@ -83 +87 @@
 
     KapaInitGraph(&graphdata);
-    KapaClearPlots(*kapaID);
+    KapaClearSections(*kapaID);
     KapaSetFont(*kapaID, "helvetica", 14);
 
@@ -89 +93 @@
     section.x = 0.0; section.y = 0.0;
     section.dx = .45, section.dy = .45;
+    section.bg = KapaColorByName("white");
     graphdata.ptype = 7;
     graphdata.style = 2;
+    graphdata.etype |= 0x01;
 
     KapaSetSection(*kapaID, &section);
-    if(!scaleGraphdata(x, xVec, &graphdata, npts)) return 0;
+    if(!scaleGraphdata(x, dx, &graphdata, npts)) return 0;
     KapaSetLimits(*kapaID, &graphdata);
     KapaBox(*kapaID, &graphdata);
@@ -100 +106 @@
     KapaPrepPlot(*kapaID, npts, &graphdata);
     KapaPlotVector(*kapaID, npts, x, "x");
-    KapaPlotVector(*kapaID, npts, xVec, "y");
+    KapaPlotVector(*kapaID, npts, dx, "y");
+    KapaPlotVector(*kapaID, npts, dPos, "dym");
+    KapaPlotVector(*kapaID, npts, dPos, "dyp");
 
     section.x = .5; section.y = 0; section.name="1";
     KapaSetSection(*kapaID, &section);
-    if(!scaleGraphdata(x, yVec, &graphdata, npts)) return 0;
+    if(!scaleGraphdata(x, dy, &graphdata, npts)) return 0;
     KapaSetLimits(*kapaID, &graphdata);
     KapaBox(*kapaID, &graphdata);
@@ -111 +119 @@
     KapaPrepPlot(*kapaID, npts, &graphdata);
     KapaPlotVector(*kapaID, npts, x, "x");
-    KapaPlotVector(*kapaID, npts, yVec, "y");
+    KapaPlotVector(*kapaID, npts, dy, "y");
+    KapaPlotVector(*kapaID, npts, dPos, "dym");
+    KapaPlotVector(*kapaID, npts, dPos, "dyp");
 
     section.x = .0; section.y = .5; section.name="2";
     KapaSetSection(*kapaID, &section);
-    if(!scaleGraphdata(y, xVec, &graphdata, npts)) return 0;;
+    if(!scaleGraphdata(y, dx, &graphdata, npts)) return 0;;
     KapaSetLimits(*kapaID, &graphdata);
     KapaBox(*kapaID, &graphdata);
@@ -122 +132 @@
     KapaPrepPlot(*kapaID, npts, &graphdata);
     KapaPlotVector(*kapaID, npts, y, "x");
-    KapaPlotVector(*kapaID, npts, xVec, "y");
+    KapaPlotVector(*kapaID, npts, dx, "y");
+    KapaPlotVector(*kapaID, npts, dPos, "dym");
+    KapaPlotVector(*kapaID, npts, dPos, "dyp");
 
     section.x = .5; section.y = .5; section.name="3";
     KapaSetSection(*kapaID, &section);
-    if(!scaleGraphdata(y, yVec, &graphdata, npts)) return 0;
+    if(!scaleGraphdata(y, dy, &graphdata, npts)) return 0;
     KapaSetLimits(*kapaID, &graphdata);
     KapaBox(*kapaID, &graphdata);
@@ -133 +145 @@
     KapaPrepPlot(*kapaID, npts, &graphdata);
     KapaPlotVector(*kapaID, npts, y, "x");
-    KapaPlotVector(*kapaID, npts, yVec, "y");
-
-    return 1;
-}
-
-
-/**Plot a vector field*/
-static int plotVectorField(float x[], float y[],
-                           float xVec[], float yVec[],
-                           int npts, int *kapaID, double maxVecLength) {
-
-    if(!plotVector) return 1;
+    KapaPlotVector(*kapaID, npts, dy, "y");
+    KapaPlotVector(*kapaID, npts, dPos, "dym");
+    KapaPlotVector(*kapaID, npts, dPos, "dyp");
+
+    return TRUE;
+}
+
+
+/** Plot a vector field (circles at vector origin, scaled lines giving vector directions */
+int relastroVisualVectorField(int *kapaID, float *x, float *y, float *dx, float *dy, int npts, double maxVecLength) {
 
     Graphdata graphdata;
-    float singleX[2], singleY[2];
+    float *xVec, *yVec;
     float vecScaleFactor;
     float graphSize;
     int i;
     char plotTitle[50];
-    sprintf(plotTitle, "Maximum Vector Size = %5.1e", maxVecLength);
-
+
+    if (!npts) return FALSE;
     if (!initWindow(kapaID)) return 0;
 
+    snprintf(plotTitle, 50, "Maximum Vector Size = %5.1e", maxVecLength);
+
     KapaInitGraph(&graphdata);
-    KapaClearPlots(*kapaID);
-    if(!scaleGraphdata(x, y, &graphdata, npts)) return 0;
+    KapaClearSections(*kapaID);
+    KapaSetFont(*kapaID, "helvetica", 14);
+
+    if (!scaleGraphdata(x, y, &graphdata, npts)) return 0;
 
     graphSize = graphdata.xmax - graphdata.xmin;
@@ -164 +178 @@
         graphSize = graphdata.ymax - graphdata.ymin;
     }
-
     vecScaleFactor = graphSize * 0.02 / (float)maxVecLength;
-#ifdef TESTING
-    fprintf(stderr, "GraphSize: %e\n maxVecLength: %e\n vecScaleFactor: %e\n",
-            graphSize, maxVecLength, vecScaleFactor);
-#endif
-
+
+    // fprintf(stderr, "GraphSize: %e\n maxVecLength: %e\n vecScaleFactor: %e\n", graphSize, maxVecLength, vecScaleFactor);
 
     KapaSetFont (*kapaID, "helvetica", 14);
@@ -179 +189 @@
     graphdata.ptype = 7;
     graphdata.style = 2;
+    graphdata.color = KapaColorByName("black");
     KapaPrepPlot(*kapaID, npts, &graphdata);
     KapaPlotVector(*kapaID, npts, x, "x");
     KapaPlotVector(*kapaID, npts, y, "y");
 
-    //plot each vector individually
-    graphdata.ptype = 0;
-    graphdata.style = 0;
+    ALLOCATE (xVec, float, 2*npts);
+    ALLOCATE (yVec, float, 2*npts);
+    for (i = 0; i < npts; i++) {
+      xVec[2*i + 0] = x[i];
+      yVec[2*i + 0] = y[i];
+      xVec[2*i + 1] = x[i] + dx[i] * vecScaleFactor;
+      yVec[2*i + 1] = y[i] + dy[i] * vecScaleFactor;
+    }
+
+    graphdata.ptype = 100; // line segements by point pair
+    graphdata.style = 2;
     graphdata.color = KapaColorByName("blue");
-    for(i = 0; i < npts; i++) {
-        singleX[0] = x[i];
-        singleY[0] = y[i];
-        singleX[1] = x[i] + xVec[i] * vecScaleFactor;
-        singleY[1] = y[i] + yVec[i] * vecScaleFactor;
-        KapaPrepPlot(*kapaID, 2, &graphdata);
-        KapaPlotVector(*kapaID, 2, singleX, "x");
-        KapaPlotVector(*kapaID, 2, singleY, "y");
-    }
-    return 1;
-}
-
-int relastroVisualPlotScatter(double values[], double thresh, int npts) {
-    float *x, *data;
-    int i;
-    float xline[2], yline[2];
-    if (!isVisual || !plotScatter) return 1;
-    if (!initWindow(&kapa2)) return 0;
-
-    ALLOCATE(x, float, npts);
-    ALLOCATE(data, float, npts);
-
-    for(i = 0; i < npts; i++) {
-        x[i] = i;
-        data[i] = (float) values[i];
-    }
+    KapaPrepPlot  (*kapaID, 2*npts, &graphdata);
+    KapaPlotVector(*kapaID, 2*npts, xVec, "x");
+    KapaPlotVector(*kapaID, 2*npts, yVec, "y");
+
+    free (xVec);
+    free (yVec);
+    return TRUE;
+}
+
+int relastroVisualPlotScatter(int *kapaID, float *dXfit, float *dYfit, float *mag, int npts) {
 
     Graphdata graphdata;
     KapaSection section;
-    section.x = 0; section.y = 0; section.dx = 1; section.dy = 1;
-    section.name = "junk";
+
+    if (!initWindow(kapaID)) return 0;
 
     KapaInitGraph(&graphdata);
-    KapaClearPlots(kapa2);
-    KapaSetSection(kapa2, &section);
-
-    graphdata.ptype = 0;
-    graphdata.style = 0;
-    graphdata.xmin = 0;
-    graphdata.xmax = npts;
-    graphdata.ymin = 0;
-    graphdata.ymax = data[npts-1];
-
-    KapaSetFont(kapa2, "helvetica", 14);
-    KapaSetLimits(kapa2, &graphdata);
-    KapaBox(kapa2, &graphdata);
-    KapaSendLabel( kapa2, "Object", KAPA_LABEL_XM);
-    KapaSendLabel( kapa2, "Offset(pixels)", KAPA_LABEL_YM);
-    KapaSendLabel( kapa2, "Astrometric Offset with fit cutoff",
-                   KAPA_LABEL_XP);
-    KapaPrepPlot(kapa2, npts, &graphdata);
-    KapaPlotVector(kapa2, npts, x, "x");
-    KapaPlotVector(kapa2, npts, data, "y");
-
-    graphdata.color = KapaColorByName("red");
-    KapaPrepPlot(kapa2, 2, &graphdata);
-    yline[0] = (float) thresh;
-    yline[1] = (float) thresh;
-    xline[0] = graphdata.xmin;
-    xline[1] = graphdata.xmax;
-    KapaPlotVector(kapa2, 2, xline, "x");
-    KapaPlotVector(kapa2, 2, yline, "y");
-
-    askUser(&plotScatter);
-    return 1;
-}
-
-
+    KapaClearSections(*kapaID);
+    KapaSetFont(*kapaID, "helvetica", 14);
+
+    section.name = "a";
+    section.x = 0.0; section.y = 0.0; section.dx = 1.0; section.dy = 0.5;
+    section.bg = KapaColorByName("white");
+    KapaSetSection(*kapaID, &section);
+
+    graphdata.ptype = 2;
+    graphdata.style = 2;
+
+    if(!scaleGraphdata(mag, dXfit, &graphdata, npts)) return 0;
+    KapaSetLimits(*kapaID, &graphdata);
+    KapaBox(*kapaID, &graphdata);
+    KapaSendLabel(*kapaID, "mag", KAPA_LABEL_XM);
+    KapaSendLabel(*kapaID, "dXfit", KAPA_LABEL_YM);
+    KapaPrepPlot(*kapaID, npts, &graphdata);
+    KapaPlotVector(*kapaID, npts, mag, "x");
+    KapaPlotVector(*kapaID, npts, dXfit, "y");
+
+    section.name = "b";
+    section.x = 0.0; section.y = 0.5; section.dx = 1.0; section.dy = 0.5;
+    section.bg = KapaColorByName("white");
+    KapaSetSection(*kapaID, &section);
+
+    graphdata.ptype = 2;
+    graphdata.style = 2;
+
+    if(!scaleGraphdata(mag, dYfit, &graphdata, npts)) return 0;
+    KapaSetLimits(*kapaID, &graphdata);
+    KapaBox(*kapaID, &graphdata);
+    KapaSendLabel(*kapaID, "mag", KAPA_LABEL_XM);
+    KapaSendLabel(*kapaID, "dYfit", KAPA_LABEL_YM);
+    KapaPrepPlot(*kapaID, npts, &graphdata);
+    KapaPlotVector(*kapaID, npts, mag, "x");
+    KapaPlotVector(*kapaID, npts, dYfit, "y");
+
+    return TRUE;
+}
 
 /** plot raw vs ref (L, M). Only those whose distance is < drMax are used in fit*/
-int relastroVisualPlotRawRef(StarData *raw, StarData *ref, double dRmax, int numObj) {
-
-  if( !isVisual || !plotRawRef) return 1;
-  if( !initWindow(&kapa)) return 0;
+int relastroVisualPlotFittedStars(int *kapaID, float *rawX, float *rawY, float *refX, float *refY, int numNoFit, float *rawXfit, float *rawYfit, float *refXfit, float *refYfit, int numFit) {
+
+  Graphdata graphdata;
+
+  if (!initWindow(kapaID)) return 0;
+
+  KapaInitGraph(&graphdata);
+  KapaClearPlots(*kapaID);
+
+  graphdata.ptype = 7;
+  graphdata.style = 2;
+
+  if(!scaleGraphdata(rawXfit, rawYfit, &graphdata, numFit)) {
+      fprintf(stderr, "Not enough finite points for plotting");
+      return 0;
+  }
+
+  KapaSetFont(*kapaID, "helvetica", 14);
+  KapaSetLimits(*kapaID, &graphdata);
+  KapaBox(*kapaID, &graphdata);
+  KapaSendLabel( *kapaID, "X", KAPA_LABEL_XM);
+  KapaSendLabel( *kapaID, "Y", KAPA_LABEL_YM);
+  KapaSendLabel( *kapaID, "orange, red, green, blue: (raw, ref), (nofit, fit)",
+                 KAPA_LABEL_XP);
+
+  graphdata.color = KapaColorByName("orange");
+  graphdata.size = 1;
+  KapaPrepPlot(*kapaID, numNoFit, &graphdata);
+  KapaPlotVector(*kapaID, numNoFit, rawX, "x");
+  KapaPlotVector(*kapaID, numNoFit, rawY, "y");
+
+  graphdata.color = KapaColorByName("red");
+  graphdata.size = 2;
+  KapaPrepPlot(*kapaID, numNoFit, &graphdata);
+  KapaPlotVector(*kapaID, numNoFit, refX, "x");
+  KapaPlotVector(*kapaID, numNoFit, refY, "y");
+
+  graphdata.color = KapaColorByName("green");
+  graphdata.size = 1;
+  KapaPrepPlot(*kapaID, numFit,  &graphdata);
+  KapaPlotVector(*kapaID, numFit, rawXfit, "x");
+  KapaPlotVector(*kapaID, numFit, rawYfit, "y");
+
+  graphdata.color = KapaColorByName("blue");
+  graphdata.size = 2;
+  KapaPrepPlot(*kapaID, numFit, &graphdata);
+  KapaPlotVector(*kapaID, numFit, refXfit, "x");
+  KapaPlotVector(*kapaID, numFit, refYfit, "y");
+
+  return TRUE;
+}
+
+/** create various plots using the data in raw and ref **/
+int relastroVisualPlotChipFit(StarData *raw, StarData *ref, double dRmax, int numObj) {
 
   float *rawX, *rawY,  *refX, *refY;
   float *rawXfit, *rawYfit, *refXfit, *refYfit;
   float *magRaw, *magRef, *magRawfit, *magReffit;
-  float *xVec, *yVec;
+  float *dXfit, *dYfit, *dPos;
   int numFit = 0, numNoFit = 0;
   double dL, dM, dR;
+
+  if (!isVisual) return TRUE;
 
   ALLOCATE(rawX,      float, numObj);
@@ -273 +332 @@
   ALLOCATE(refX,      float, numObj);
   ALLOCATE(refY,      float, numObj);
+  ALLOCATE(magRaw,    float, numObj);
+  ALLOCATE(magRef,    float, numObj);
+
   ALLOCATE(rawXfit,   float, numObj);
   ALLOCATE(rawYfit,   float, numObj);
   ALLOCATE(refXfit,   float, numObj);
   ALLOCATE(refYfit,   float, numObj);
-  ALLOCATE(magRaw,    float, numObj);
-  ALLOCATE(magRef,    float, numObj);
   ALLOCATE(magRawfit, float, numObj);
   ALLOCATE(magReffit, float, numObj);
+  ALLOCATE(dXfit,     float, numObj);
+  ALLOCATE(dYfit,     float, numObj);
+  ALLOCATE(dPos,      float, numObj);
 
   int i;
   for(i = 0; i < numObj; i++) {
-    if  (raw[i].mask) continue;
+    if (raw[i].mask) continue; // XXX 
 
     dL = raw[i].L - ref[i].L;
     dM = raw[i].M - ref[i].M;
     dR = hypot (dL, dM);
+
+    // XXX change the selection to a mask-based thing
     if (dR > dRmax) {
+      // UNFITTED values
       rawX[numNoFit] = raw[i].X;
       rawY[numNoFit] = raw[i].Y;
@@ -298 +364 @@
       numNoFit++;
     } else {
-      rawXfit[numFit] = raw[i].X;
-      rawYfit[numFit] = raw[i].Y;
-      refXfit[numFit] = ref[i].X;
-      refYfit[numFit] = ref[i].Y;
-      magRaw[numFit] = raw[i].Mag;
-      magRef[numFit] = ref[i].Mag;
+      // FITTED values
+      rawXfit[numFit]   = raw[i].X;
+      rawYfit[numFit]   = raw[i].Y;
+      refXfit[numFit]   = ref[i].X;
+      refYfit[numFit]   = ref[i].Y;
+      magRawfit[numFit] = raw[i].Mag;
+      magReffit[numFit] = ref[i].Mag;
+      dPos[numFit]      = ref[i].dPos;
+      dXfit[numFit]     = raw[i].X - ref[i].X;
+      dYfit[numFit]     = raw[i].Y - ref[i].Y;
       numFit++;
     }
@@ -310 +380 @@
   if (numFit == 0) return 0;
 
-  Graphdata graphdata;
-
-  KapaInitGraph(&graphdata);
-  KapaClearPlots(kapa);
-
-  graphdata.ptype = 7;
-  graphdata.style = 2;
-
-  if(!scaleGraphdata(rawXfit, rawYfit, &graphdata, numFit)) {
-      fprintf(stderr, "Not enough finite points for plotting");
-      return 0;
-  }
-
-  KapaSetFont(kapa, "helvetica", 14);
-  KapaSetLimits(kapa, &graphdata);
-  KapaBox(kapa, &graphdata);
-  KapaSendLabel( kapa, "X", KAPA_LABEL_XM);
-  KapaSendLabel( kapa, "Y", KAPA_LABEL_YM);
-  KapaSendLabel( kapa, "orange, red, green, blue: (raw, ref), (nofit, fit)",
-                 KAPA_LABEL_XP);
-
-  graphdata.color = KapaColorByName("orange");
-  graphdata.size = 1;
-  KapaPrepPlot(kapa, numNoFit, &graphdata);
-  KapaPlotVector(kapa, numNoFit, rawX, "x");
-  KapaPlotVector(kapa, numNoFit, rawY, "y");
-
-  graphdata.color = KapaColorByName("red");
-  graphdata.size = 2;
-  KapaPrepPlot(kapa, numNoFit, &graphdata);
-  KapaPlotVector(kapa, numNoFit, refX, "x");
-  KapaPlotVector(kapa, numNoFit, refY, "y");
-
-  graphdata.color = KapaColorByName("green");
-  graphdata.size = 1;
-  KapaPrepPlot(kapa, numFit,  &graphdata);
-  KapaPlotVector(kapa, numFit, rawXfit, "x");
-  KapaPlotVector(kapa, numFit, rawYfit, "y");
-
-  graphdata.color = KapaColorByName("blue");
-  graphdata.size = 2;
-  KapaPrepPlot(kapa, numFit, &graphdata);
-  KapaPlotVector(kapa, numFit, refXfit, "x");
-  KapaPlotVector(kapa, numFit, refYfit, "y");
-
-  ALLOCATE(xVec, float, numFit);
-  ALLOCATE(yVec, float, numFit);
-
-  //plot the fitted objects as vectors
-  for(i = 0; i < numFit; i++) {
-      xVec[i] = rawXfit[i] - refXfit[i];
-      yVec[i] = rawYfit[i] - refYfit[i];
-  }
-
-  plotVectorField(rawXfit, rawYfit, xVec, yVec, numFit, &kapa3, dRmax);
-  if(!residPlot(rawXfit, rawYfit, xVec, yVec, numFit, &kapa2)) {
-      fprintf(stderr, "Unable to plot residuals");
-      return 0;
-  }
-
-  FREE(xVec);
-  FREE(yVec);
+  // 4-panel plot of x vs dx, y vs dx, x vs dy, y vs dy 
+  relastroVisualRawRef(&kapa1, rawXfit, rawYfit, dXfit, dYfit, dPos, numFit);
+
+  // vector line plot for the fit
+  relastroVisualVectorField(&kapa2, rawXfit, rawYfit, dXfit, dYfit, numFit, dRmax);
+
+  // dXfit, dYfit vs mag
+  relastroVisualPlotScatter(&kapa3, dXfit, dYfit, magRawfit, numFit);
+
+  // plot the positions of the fitted stars on the chip
+  relastroVisualPlotFittedStars(&kapa4, rawX, rawY, refX, refY, numNoFit, rawXfit, rawYfit, refXfit, refYfit, numFit);
 
   askUser(&plotRawRef);
 
+  FREE(dXfit);
+  FREE(dYfit);
+  FREE(dPos);
   FREE(rawX);
   FREE(rawY);
@@ -388 +410 @@
   FREE(magReffit);
 
-  return 1;
+  return TRUE;
 }
 
@@ -403 +425 @@
   KapaSection section;
   
-  if (!isVisual || !plotOutliers) return 1;
-  if (!initWindow(&kapa)) return 0;
+  if (!isVisual || !plotOutliers) return TRUE;
+  if (!initWindow(&kapa1)) return 0;
   
   // populate vectors
@@ -459 +481 @@
   section.x = 0; section.y = 0; section.dx = 1; section.dy = 1;
   section.name = "junk";
+  section.bg = KapaColorByName("white");
   
   KapaInitGraph(&graphdata);
-  KapaClearPlots(kapa);
-  KapaSetFont(kapa, "helvetica", 14);
+  KapaClearPlots(kapa1);
+  KapaSetFont(kapa1, "helvetica", 14);
  
   graphdata.ptype = 7;
@@ -472 +495 @@
   graphdata.ymax = ymax;
 
-  KapaSetSection(kapa, &section);
-  KapaSetLimits(kapa, &graphdata);
-  KapaBox(kapa, &graphdata);
-
-  KapaSendLabel( kapa, "RA (arcsec)", KAPA_LABEL_XM);
-  KapaSendLabel( kapa, "Dec (arcsec)", KAPA_LABEL_YM);
-  KapaSendLabel( kapa, "Points flagged as outliers (red)",
+  KapaSetSection(kapa1, &section);
+  KapaSetLimits(kapa1, &graphdata);
+  KapaBox(kapa1, &graphdata);
+
+  KapaSendLabel( kapa1, "RA (arcsec)", KAPA_LABEL_XM);
+  KapaSendLabel( kapa1, "Dec (arcsec)", KAPA_LABEL_YM);
+  KapaSendLabel( kapa1, "Points flagged as outliers (red)",
                  KAPA_LABEL_XP);
 
   graphdata.color = KapaColorByName("green");
-  KapaPrepPlot(kapa, Nin, &graphdata);
-  KapaPlotVector(kapa, Nin, Rin, "x");
-  KapaPlotVector(kapa, Nin, Din, "y");
+  KapaPrepPlot(kapa1, Nin, &graphdata);
+  KapaPlotVector(kapa1, Nin, Rin, "x");
+  KapaPlotVector(kapa1, Nin, Din, "y");
 
   graphdata.color = KapaColorByName("red");
-  KapaPrepPlot(kapa, Nout, &graphdata);
-  KapaPlotVector(kapa, Nout, Rout, "x");
-  KapaPlotVector(kapa, Nout, Dout, "y");
+  KapaPrepPlot(kapa1, Nout, &graphdata);
+  KapaPlotVector(kapa1, Nout, Rout, "x");
+  KapaPlotVector(kapa1, Nout, Dout, "y");
 
   graphdata.color = KapaColorByName("black");
   graphdata.ptype = 0;
   graphdata.style = 0;
-  KapaPrepPlot(kapa, 100, &graphdata);
-  KapaPlotVector(kapa, 100, xCirc, "x");
-  KapaPlotVector(kapa, 100, yCirc, "y");
-
-
-  
+  KapaPrepPlot(kapa1, 100, &graphdata);
+  KapaPlotVector(kapa1, 100, xCirc, "x");
+  KapaPlotVector(kapa1, 100, yCirc, "y");
+
   askUser(&plotOutliers);
 
@@ -510 +531 @@
 }
   
-  
+# if (0)
+int relastroVisualSummaryChips() {
+
+  // plot the dXsys, dYsys histograms
+
+  // plot x vs dx, y vs dy, etc for all mosaics
+
+  // plot a map of median star scatter
+
+}
+# endif
+

branches/eam_branches/ipp-20101205/Ohana/src/relastro/src/select_images.c

@@ -33 +33 @@
   struct timeval start, stop;
   
-  double RmaxSkyRegion, RminSkyRegion, DminSkyRegion, DmaxSkyRegion;
+  double RmaxSkyRegion, RminSkyRegion, RmidSkyRegion, DminSkyRegion, DmaxSkyRegion;
 
   double *RmaxSky;
@@ -102 +102 @@
     DmaxSkyRegion = MAX(DmaxSkyRegion, skylist[0].regions[i][0].Dmax);
   }
+  RmidSkyRegion = 0.5*(RminSkyRegion + RmaxSkyRegion);
   MARKTIME("create sky region coords: %f sec\n", dtime);
 
@@ -169 +170 @@
     double DminImage = +90.0;
     double DmaxImage = -90.0;
+    int leftside = FALSE;
     for (j = 0; j < 5; j++) {
       XY_to_RD (&Ri[j], &Di[j], Xi[j], Yi[j], &timage[i].coords);
+      Ri[j] = ohana_normalize_angle (Ri[j]);
+      while (Ri[j] < RminSkyRegion) { Ri[j] += 360.0; }
+      while (Ri[j] > RmaxSkyRegion) { Ri[j] -= 360.0; }
+      if (j == 0) {
+        leftside = (Ri[j] < RmidSkyRegion);
+      } else {
+        if (  leftside && (Ri[j] > RmidSkyRegion + 90)) { Ri[j] -= 360.0; }
+        if (! leftside && (Ri[j] < RmidSkyRegion - 90)) { Ri[j] += 360.0; }
+      }
+
       RminImage = MIN(RminImage, Ri[j]);
       RmaxImage = MAX(RmaxImage, Ri[j]);
@@ -232 +244 @@
       if (RESET) {
         // XXX do we need / want to do this in relastro?
-        assignMcal (&image[nimage], (double *) NULL, -1);
-        image[nimage].dMcal = NAN;
+        // assignMcal (&image[nimage], (double *) NULL, -1);
+        // image[nimage].Mcal = NAN;
+        // image[nimage].dMcal = NAN;
         image[nimage].flags &= ~badImage;
       }