Index: branches/eam_branches/relphot.20210521/src/MosaicOps.c
===================================================================
--- branches/eam_branches/relphot.20210521/src/MosaicOps.c	(revision 41602)
+++ branches/eam_branches/relphot.20210521/src/MosaicOps.c	(revision 41603)
@@ -50,104 +50,4 @@
 # undef COMPARE
 
-}
-
-/* find mosaic frames (unique time periods & photcode name matches mosaic) */
-void initMosaics_old (Image *image, off_t Nimage) {
-
-  off_t i, j, status, found, NMOSAIC, *MosaicN_IMAGE;
-  unsigned int start, stop;
-  char *pname;
-
-  if (!MOSAIC_ZEROPT) return;
-
-  Nmosaic = 0;
-  NMOSAIC = 10;
-  ALLOCATE (mosaic, Mosaic, NMOSAIC);
-
-  ALLOCATE (MosaicToImage, off_t *, NMOSAIC);
-  ALLOCATE (MosaicN_Image, off_t,   NMOSAIC);
-  ALLOCATE (MosaicN_IMAGE, off_t,   NMOSAIC);
-
-  ALLOCATE (ImageToMosaic, off_t, Nimage); // mosaic to which image belongs
-
-  /* a 'mosaic' in relphot is (unlike relastro) a virtual concept: there is no 
-   * entry in the image table that represents this mosaic.  Instead, it is an
-   * internal construct that defines a group of related images 
-   */
-
-  /* generate list of unique mosaics */
-  for (i = 0; i < Nimage; i++) {
-    ImageToMosaic[i] = -1;
-
-    /* select valid mosaic images by photcode */
-    pname = GetPhotcodeNamebyCode (image[i].photcode);
-    status = strncmp (pname, MOSAICNAME, strlen (MOSAICNAME));
-    if (status) continue;
-
-    /* set image time range */
-    start = image[i].tzero - MAX(0.01*image[i].trate*image[i].NY, 1);
-    stop  = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1);
-
-    /* find existing mosaic with this time range */
-    found = FALSE;
-    for (j = 0; !found && (j < Nmosaic); j++) { 
-      if (stop  < mosaic[j].start) continue;
-      if (start > mosaic[j].stop)  continue;
-      found = TRUE;
-
-      // add reference from image to mosaic
-      ImageToMosaic[i] = j;
-
-      /* add image to mosaic image list */
-      MosaicToImage[j][MosaicN_Image[j]] = i;
-      MosaicN_Image[j] ++;
-      if (MosaicN_Image[j] == MosaicN_IMAGE[j]) {
-	MosaicN_IMAGE[j] += 10;
-	REALLOCATE (MosaicToImage[j], off_t, MosaicN_IMAGE[j]);
-      }
-
-    }
-    if (found) continue;
-    
-    /* a new mosaic, define ranges -- preserve the original values incase this image is not used */
-    mosaic[Nmosaic].start     = start;
-    mosaic[Nmosaic].stop      = stop;
-    mosaic[Nmosaic].McalPSF   = 0.0; // note : at the end, mosaic.Mcal is added back to the input images
-    mosaic[Nmosaic].McalAPER  = 0.0; // note : mosaic stores only offsets relative to the original image values
-    mosaic[Nmosaic].dMcal     = 0.0; // note : at the end, mosaic.Mcal is added back to the input images
-    mosaic[Nmosaic].dMsys     = 0.0;
-    mosaic[Nmosaic].McalChiSq = 0.0;// NAN or 0.0?
-    mosaic[Nmosaic].flags     = image[i].flags;
-    mosaic[Nmosaic].secz      = image[i].secz;
-    mosaic[Nmosaic].photcode  = GetPhotcodeEquivCodebyCode (image[i].photcode);
-
-    // XXX do we need to do something about flag consistency across a mosaic?
-
-    /* add image to mosaic image list */
-    MosaicN_IMAGE[Nmosaic] = 10;
-    MosaicN_Image[Nmosaic] = 1;
-    ALLOCATE (MosaicToImage[Nmosaic], off_t, MosaicN_IMAGE[Nmosaic]);
-    MosaicToImage[Nmosaic][0] = i;
-
-    // add reference from image to mosaic
-    ImageToMosaic[i] = Nmosaic;
-    
-    Nmosaic ++;
-    if (Nmosaic == NMOSAIC) {
-      NMOSAIC += 10;
-      REALLOCATE (mosaic,  Mosaic,  NMOSAIC);
-      REALLOCATE (MosaicToImage,  off_t *, NMOSAIC);
-      REALLOCATE (MosaicN_Image, off_t,   NMOSAIC);
-      REALLOCATE (MosaicN_IMAGE, off_t,   NMOSAIC);
-    }
-  }
-
-  // free this or not?
-  free (MosaicN_IMAGE);
-
-  initMosaicGrid (image, Nimage);
-
-  fprintf (stderr, "matched %d images to %d mosaics\n", (int) Nimage, (int) Nmosaic);
-  return;
 }
 
@@ -318,5 +218,5 @@
   free (startTimesMosaic);
 
-  initMosaicGrid (subset, Nsubset);
+  initMosaicMcal (subset, Nsubset);
 
   fprintf (stderr, "matched %d images to %d mosaics, %d simple chips not matched to mosaics\n", (int) (Nsubset - Nsimple), (int) Nmosaic, (int) Nsimple);
@@ -460,5 +360,5 @@
 
   if (mergeMcal) {
-    initMosaicGrid (image, Nimage);
+    initMosaicMcal (image, Nimage);
   }
 
@@ -511,4 +411,7 @@
 }
 
+// this function sets mosaic->coords to the median of the individual chips.  This
+// coordinate frame is used by the parallel region analysis to assign exposures (mosaics)
+// to specific machines by a single center (rather than individual chips)
 void setMosaicCenters (Image *image, off_t Nimage) {
   OHANA_UNUSED_PARAM(Nimage);
@@ -599,95 +502,24 @@
     mosaic[i].coords.cdelt1 = 1.0 / 3600.0;
     mosaic[i].coords.cdelt2 = 1.0 / 3600.0;
-
-    mosaic[i].McalPSF   = 0.0; 
-    mosaic[i].McalAPER  = 0.0; 
-    mosaic[i].dMcal     = 0.0; 
-    mosaic[i].dMsys     = 0.0;
-    mosaic[i].McalChiSq = 0.0; 
   }
   return;
 }
 
-void initMosaicGrid (Image *image, off_t Nimage) {
+void initMosaicMcal (Image *image, off_t Nimage) {
   OHANA_UNUSED_PARAM(Nimage);
 
-  /* find max dR, dD range for all mosaics */
-  /* define mosaic.coords to cover dR, dD */
-  /* send results to initGridBins */
-
-  off_t i, j, m, NX, NY;
-  int dXmax, dYmax;
-  double dS, dX, dY;
-  double R, D, Rmin, Rmax, Dmin, Dmax;
   double McalPSF, McalAPER, dMcal, McalChiSq;
 
   fprintf (stderr, "*** moving Mcal from image.Mcal to mosaic.Mcal ***\n");
 
-  dXmax = dYmax = 0.0;
-  for (i = 0; i < Nmosaic; i++) {
-    Dmin = Rmin =  360.0;
-    Dmax = Rmax = -360.0;
-    dS = 0.0;
+  for (off_t i = 0; i < Nmosaic; i++) {
     McalPSF = McalAPER = dMcal = McalChiSq = 0;
-    for (j = 0; j < MosaicN_Image[i]; j++) {
-      m = MosaicToImage[i][j];
-      NX = image[m].NX;
-      NY = image[m].NY;
-      dS += hypot (image[m].coords.cdelt1*image[m].coords.pc1_1, image[m].coords.cdelt1*image[m].coords.pc2_1);
-
-      OhanaProjection proj = GetProjection (image[m].coords.ctype);
-      if ((proj == PROJ_WRP) && !image[m].coords.mosaic) {
-	XY_to_LM (&R, &D, 0.0, 0.0, &image[m].coords);
-      } else {
-	XY_to_RD (&R, &D, 0.0, 0.0, &image[m].coords);
-      }
-      Rmin = MIN (Rmin, R);
-      Rmax = MAX (Rmax, R);
-      Dmin = MIN (Dmin, D);
-      Dmax = MAX (Dmax, D);
-
-      if ((proj == PROJ_WRP) && !image[m].coords.mosaic) {
-	XY_to_LM (&R, &D, (double) NX, 0.0, &image[m].coords);
-      } else {
-	XY_to_RD (&R, &D, (double) NX, 0.0, &image[m].coords);
-      }
-      Rmin = MIN (Rmin, R);
-      Rmax = MAX (Rmax, R);
-      Dmin = MIN (Dmin, D);
-      Dmax = MAX (Dmax, D);
-
-      if ((proj == PROJ_WRP) && !image[m].coords.mosaic) {
-	XY_to_LM (&R, &D, (double) NX, (double) NY, &image[m].coords);
-      } else {
-	XY_to_RD (&R, &D, (double) NX, (double) NY, &image[m].coords);
-      }
-      Rmin = MIN (Rmin, R);
-      Rmax = MAX (Rmax, R);
-      Dmin = MIN (Dmin, D);
-      Dmax = MAX (Dmax, D);
-
-      if ((proj == PROJ_WRP) && !image[m].coords.mosaic) {
-	XY_to_LM (&R, &D, 0.0, (double) NY, &image[m].coords);
-      } else {
-	XY_to_RD (&R, &D, 0.0, (double) NY, &image[m].coords);
-      }
-      Rmin = MIN (Rmin, R);
-      Rmax = MAX (Rmax, R);
-      Dmin = MIN (Dmin, D);
-      Dmax = MAX (Dmax, D);
-
-      // XXX : this probably does not handle mosaics at RA = 0,360 well
-
-      /* we are using mosaic.Mcal, not image.Mcal. reset image.Mcal */
-
-      // XXX: how does this work with UBERCAL?  We want to keep the Mcal values supplied by ubercal, but 
-      // solve for a single offset for each exposure (Mosaic.Mcal).  
-      // we also want to keep the flat-field terms for each exposure (regardless of ubercal or not)
-      // if it helps, note that ubercal uses a single zp per exposure, so the mean of those values is the same as the value
+    for (off_t j = 0; j < MosaicN_Image[i]; j++) {
+      off_t m = MosaicToImage[i][j];
 
       if (!isfinite(image[m].McalPSF)) {
-	image[m].McalPSF = 0.0;
-	image[m].McalAPER = 0.0;
-	image[m].dMcal = 0.0;
+	image[m].McalPSF   = 0.0;
+	image[m].McalAPER  = 0.0;
+	image[m].dMcal     = 0.0;
 	image[m].McalChiSq = 0.0;
 	fprintf (stderr, "warning: resetting NAN value for Mcal %s\n", image[m].name);
@@ -703,15 +535,5 @@
       image[m].dMcal     = NAN;
       image[m].McalChiSq = NAN;
-
-    }
-    dS /= MosaicN_Image[i];
-
-    InitCoords (&mosaic[i].coords, "DEC--TAN");
-    mosaic[i].coords.crval1 = Rmin;
-    mosaic[i].coords.crval2 = Dmin;
-    mosaic[i].coords.cdelt1 = dS;
-    mosaic[i].coords.cdelt2 = dS;
-
-    RD_to_XY (&dX, &dY, Rmax, Dmax, &mosaic[i].coords);
+    }
 
     mosaic[i].McalPSF   = McalPSF / MosaicN_Image[i];
@@ -720,9 +542,4 @@
     mosaic[i].McalChiSq = McalChiSq / MosaicN_Image[i];
   }
-  if (!USE_GRID) return;
-
-  dXmax = MAX (dXmax, dX);
-  dYmax = MAX (dYmax, dY);
-  initGrid (dXmax, dYmax);
   return;
 }
@@ -1231,5 +1048,5 @@
       continue;
     }
-    float Mgrid = getMgrid (m, c);
+    float Mgrid = getMgridTiny (&catalog[c].measureT[m]);
     if (isnan(Mgrid)) {
       info->Ngrid ++;
@@ -1382,5 +1199,5 @@
     float Mcal     = getMcal  (m, c, MAG_CLASS_PSF);
     // float Mgrp     = getMgrp  (m, c, catalog[c].measureT[m].airmass, NULL);
-    float Mgrid    = getMgrid (m, c);
+    float Mgrid    = getMgridTiny (&catalog[c].measureT[m]);
     float MrelPSF  = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP);
     float Mflat    = getMflat (m, c, flatcorr, catalog);
@@ -1481,4 +1298,5 @@
 	   summary.Nskip);
   free (threadinfo);
+  // XXX SetMmosInfoFree (&summary);
 
   npass_output ++;
@@ -1811,5 +1629,5 @@
       Mgrp = getMgrp  (m, c, catalog[c].measureT[m].airmass, NULL);
       if (isnan(Mgrp)) continue;
-      Mgrid = getMgrid (m, c);
+      Mgrid = getMgridTiny (&catalog[c].measureT[m]);
       if (isnan(Mgrid)) continue;
       Mrel = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP);
