# include "relphot.h"
void plot_setMcal (double *list, int Npts, StatType *stats, float clouds);
off_t findMosaic (unsigned int *startTimes, off_t Nmosaic, unsigned int start);

// see discussion in ImagesOps.c re: IDX_T

// array of mosaic definition structures
static off_t  Nmosaic;
static Mosaic *mosaic;

// relationships between the mosaics and their associated images
static off_t   *MosaicN_Image; // number of images associated with the given mosaic
static off_t  **MosaicToImage; // list of imagesa associated with the given mosaic

// mosaic index for given image : ImageToMosaic[ImageIndex] = MosaicIndex (ImageIndex : 0 < Nimage)
static off_t   *ImageToMosaic;

// elsewhere, we have loaded a set of catalogs with measures (catalog[cat].measure[meas])
// each mosaic has N_onMosaic[MosaicIndex] measurements
static off_t    *N_onMosaic;   // actual number of measurements on mosaic	 
static off_t    *N_ONMOSAIC;   // allocated number of measurements on mosaic   

// relationships between the measure,catalog set and the mosaics:
static off_t   **MeasureToMosaic; // Mosaic index from measure,catalog  : MeasureToMosaic[cat][meas] = MosaicIndex 
static off_t   **MosaicToCatalog; // catalog for given measure on mosaic : MosaicCatalog[MosaicIndex][i] = cat (i : 0 < NonMosaic[MosaicIndex])
static off_t   **MosaicToMeasure; // measure for given measure on mosaic : MosaicMeasure[MosaicIndex][i] = cat (i : 0 < NonMosaic[MosaicIndex])

// MeasureToMosaic was 'bin'
// MosaicToCatalog was 'clist'
// MosaicToMeasure was 'mlist'

// N_onMosaic was 'Nlist'
// N_ONMOSAIC was 'NLIST'

// ImageToMosaic was 'mosimage'

// MosaicN_image was 'Nimlist'
// MosaicToImage was 'imlist'

void sort_times (unsigned int *T, int N) {

# define SWAPFUNC(A,B){ unsigned int tmp; \
  tmp = T[A]; T[A] = T[B]; T[B] = tmp; \
}
# define COMPARE(A,B)(T[A] < T[B])

  OHANA_SORT (N, COMPARE, SWAPFUNC);

# undef SWAPFUNC
# 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].Mcal   = 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].Xm = 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;
}

/* find mosaic frames (unique time periods) (NOTE : we do NOT require matching photcodes...)
   this function will also identify the images NOT in the subset which belong to a selected mosaic
 */
void initMosaics (Image *subset, off_t Nsubset, Image *image, char *inSubset, off_t Nimage) {

  off_t i, j, found, NMOSAIC, *MosaicN_IMAGE;
  unsigned int start, stop, *startTimes, *startTimesMosaic;

  if (!MOSAIC_ZEROPT) return;

  INITTIME;

  /* 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 a list of all subset image start times
  ALLOCATE (startTimes, unsigned int, Nsubset);
  int Nmoschip = 0;
  for (i = 0; i < Nsubset; i++) {
    if (!isGPC1chip(subset[i].photcode)) continue;
    startTimes[Nmoschip] = subset[i].tzero;
    Nmoschip ++;
  }
  sort_times (startTimes, Nmoschip);
  MARKTIME("create array of all image obstimes: %f sec\n", dtime);
  
  Nmosaic = 0;
  NMOSAIC = 1000;
  ALLOCATE (startTimesMosaic, unsigned int, NMOSAIC);
  startTimesMosaic[0] = startTimes[0];

  // generate a list of the unique start times (these define the mosaics)
  for (i = 0; i < Nmoschip; i++) {
    myAssert (startTimes[i] >= startTimesMosaic[Nmosaic], "times out of order?");
    if (startTimes[i] == startTimesMosaic[Nmosaic]) continue;
    Nmosaic ++;
    if (Nmosaic >= NMOSAIC) {
      NMOSAIC += 1000;
      REALLOCATE (startTimesMosaic, unsigned int, NMOSAIC);
    }
    startTimesMosaic[Nmosaic] = startTimes[i];
  }
  Nmosaic ++;
  MARKTIME("create subset array of mosaic obstimes: %f sec\n", dtime);

  // now I have a list of uniq start times, and they are in order
  // create the mosaic arrays for these times
  ALLOCATE (mosaic, Mosaic, Nmosaic);

  ALLOCATE (MosaicToImage, off_t *, Nmosaic);
  ALLOCATE (MosaicN_Image, off_t,   Nmosaic);
  ALLOCATE (MosaicN_IMAGE, off_t,   Nmosaic);

  for (i = 0; i < Nmosaic; i++) {
    /* a new mosaic, define ranges */
    mosaic[i].start    = startTimesMosaic[i];
    mosaic[i].stop     = 0;
    mosaic[i].Mcal     = 0.0;
    mosaic[i].dMcal    = 0.0;
    mosaic[i].dMsys    = 0.0;
    mosaic[i].Xm       = 0.0;
    mosaic[i].flags    = 0;
    mosaic[i].secz     = NAN;
    mosaic[i].photcode = 0;
    mosaic[i].skipCal  = FALSE;
    
    memset (&mosaic[i].coords, 0, sizeof(Coords));

    MosaicN_IMAGE[i] = 10;
    MosaicN_Image[i] = 0;
    ALLOCATE (MosaicToImage[i], off_t, MosaicN_IMAGE[i]);
    MosaicToImage[i][0] = -1;
  }

  int Nskip, Nmark;
  Nskip = Nmark = 0;
  // find any mosaics (startTimesMosaic) which match unselected images 
  for (i = 0; i < Nimage; i++) {
    if (inSubset[i]) {
      Nskip ++;
      continue;
    }
    
    if (!isGPC1chip(image[i].photcode)) continue;

    /* set image time range */
    start = image[i].tzero;
    stop  = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1);

    /* find a matching mosaic */
    j = findMosaic(startTimesMosaic, Nmosaic, start);
    if (j != -1) {
      // mark this mosaic as bad
      mosaic[j].skipCal = TRUE;
      Nmark ++;
    }
  }
  fprintf (stderr, "%d total images, %d overlap skyregion & match selection criteria, %d do not overlap, but match overlapping mosaics\n", (int) Nimage, (int) Nskip, (int) Nmark);
  MARKTIME("find unselected images matching selected mosaics: %f sec\n", dtime);

  ALLOCATE (ImageToMosaic, off_t, Nsubset); // mosaic to which image belongs

  // assign each image to a mosaic
  int Nsimple = 0;
  for (i = 0; i < Nsubset; i++) {
    ImageToMosaic[i] = -1;

    if (!isGPC1chip(subset[i].photcode)) {
      Nsimple ++;
      continue;
    }

    start = subset[i].tzero;
    stop  = subset[i].tzero + MAX(1.01*subset[i].trate*subset[i].NY, 1);

    j = findMosaic(startTimesMosaic, Nmosaic, start);
    if (j == -1) {
      fprintf (stderr, "programming error? all subset images should belong to a mosaic\n");
      abort();
    }

    // add reference from image to mosaic
    ImageToMosaic[i] = j;

    // have we already found this mosaic?
    found = (MosaicN_Image[j] > 0);

    /* 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 */
    if (mosaic[j].start != start) { 
      fprintf (stderr, "error?\n");
      abort();
    }
    mosaic[j].stop     = stop;
    mosaic[j].Mcal     = 0.0;
    mosaic[j].dMcal    = 0.0;
    mosaic[j].Xm       = 0.0;
    mosaic[j].dMsys    = subset[i].flags;
    mosaic[j].flags    = subset[i].flags;
    mosaic[j].secz     = subset[i].secz;
    mosaic[j].photcode = GetPhotcodeEquivCodebyCode (subset[i].photcode);
  }

  // free this or not?
  free (MosaicN_IMAGE);
  free (startTimes);
  free (startTimesMosaic);

  initMosaicGrid (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);
  return;
}

void freeMosaics () {
  

}

/* find mosaic frames (unique time periods) (NOTE : require gpc1 chips, which is pretty limiting)
   if mergeMcal is TRUE, <image.Mcal> values will be saved on Mosaic.Mcal
 */
void makeMosaics (Image *image, off_t Nimage, int mergeMcal) {

  off_t i, j, found, NMOSAIC, *MosaicN_IMAGE;
  unsigned int start, stop, *startTimes, *startTimesMosaic;

  if (!MOSAIC_ZEROPT) return;

  INITTIME;

  /* 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 a list of all image start times
  ALLOCATE (startTimes, unsigned int, Nimage);
  int Nmoschip = 0;
  for (i = 0; i < Nimage; i++) {
    if (!isGPC1chip(image[i].photcode)) continue;
    startTimes[Nmoschip] = image[i].tzero;
    Nmoschip ++;
  }
  sort_times (startTimes, Nmoschip);
  MARKTIME("create array of all image obstimes: %f sec\n", dtime);
  
  Nmosaic = 0;
  NMOSAIC = 1000;
  ALLOCATE (startTimesMosaic, unsigned int, NMOSAIC);
  startTimesMosaic[0] = startTimes[0];

  // generate a list of the unique start times (these define the mosaics)
  for (i = 0; i < Nmoschip; i++) {
    myAssert (startTimes[i] >= startTimesMosaic[Nmosaic], "times out of order?");
    if (startTimes[i] == startTimesMosaic[Nmosaic]) continue;
    Nmosaic ++;
    if (Nmosaic >= NMOSAIC) {
      NMOSAIC += 1000;
      REALLOCATE (startTimesMosaic, unsigned int, NMOSAIC);
    }
    startTimesMosaic[Nmosaic] = startTimes[i];
  }
  Nmosaic ++;
  MARKTIME("create array of mosaic obstimes: %f sec\n", dtime);

  // now I have a list of uniq start times, and they are in order
  // create the mosaic arrays for these times
  ALLOCATE (mosaic, Mosaic, Nmosaic);

  ALLOCATE (MosaicToImage, off_t *, Nmosaic);
  ALLOCATE (MosaicN_Image, off_t,   Nmosaic);
  ALLOCATE (MosaicN_IMAGE, off_t,   Nmosaic);

  // init the mosaic array values
  for (i = 0; i < Nmosaic; i++) {
    mosaic[i].start    = startTimesMosaic[i];
    mosaic[i].stop     = 0;
    mosaic[i].Mcal     = 0.0;
    mosaic[i].dMcal    = 0.0;
    mosaic[i].dMsys    = 0.0;
    mosaic[i].Xm       = 0.0;
    mosaic[i].flags    = 0;
    mosaic[i].secz     = NAN;
    mosaic[i].photcode = 0;
    mosaic[i].skipCal  = FALSE;
    
    memset (&mosaic[i].coords, 0, sizeof(Coords));

    MosaicN_IMAGE[i] = 10;
    MosaicN_Image[i] = 0;
    ALLOCATE (MosaicToImage[i], off_t, MosaicN_IMAGE[i]);
    MosaicToImage[i][0] = -1;
  }

  ALLOCATE (ImageToMosaic, off_t, Nimage); // mosaic to which image belongs

  // assign each image to a mosaic
  int Nsimple = 0;
  for (i = 0; i < Nimage; i++) {
    ImageToMosaic[i] = -1;

    if (!isGPC1chip(image[i].photcode)) {
      Nsimple ++;
      continue;
    }

    start = image[i].tzero;
    stop  = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1);

    j = findMosaic(startTimesMosaic, Nmosaic, start);
    if (j == -1) {
      fprintf (stderr, "programming error? all image images should belong to a mosaic\n");
      abort();
    }

    // add reference from image to mosaic
    ImageToMosaic[i] = j;

    // have we already found this mosaic?
    found = (MosaicN_Image[j] > 0);

    /* 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 */
    if (mosaic[j].start != start) { 
      fprintf (stderr, "error?\n");
      abort();
    }
    mosaic[j].stop     = stop;
    mosaic[j].Mcal     = 0.0;
    mosaic[j].dMcal    = 0.0;
    mosaic[j].Xm       = 0.0;
    mosaic[j].dMsys    = image[i].flags;
    mosaic[j].flags    = image[i].flags;
    mosaic[j].secz     = image[i].secz;
    mosaic[j].photcode = GetPhotcodeEquivCodebyCode (image[i].photcode);
  }
  MARKTIME("assign images to mosaic: %f sec\n", dtime);

  // free this or not?
  free (MosaicN_IMAGE);
  free (startTimes);
  free (startTimesMosaic);

  if (mergeMcal) {
    initMosaicGrid (image, Nimage);
  }

  fprintf (stderr, "matched %d images to %d mosaics\n", (int) Nimage, (int) Nmosaic);
  return;
}

Mosaic *getMosaicForImage (off_t im) {

  if (im < 0) return NULL;
  if (!ImageToMosaic) return NULL;

  off_t m = ImageToMosaic[im];
  if (m < 0) return NULL;
  if (m >= Nmosaic) return NULL;

  return (&mosaic[m]);
}

// use bisection to find the overlapping mosaic (returns exact match)
// startTimes is a sorted, unique list of times
// start might not be in the list
off_t findMosaic (unsigned int *startTimes, off_t Nmosaic, unsigned int start) {

  off_t Nlo, Nhi, N;

  // find the last mosaic before start
  Nlo = 0; // first valid startTimes value
  Nhi = Nmosaic - 1; // last valid startTimes value

  // if start is not in this range, return -1
  if (start < startTimes[Nlo]) return (-1);
  if (start > startTimes[Nhi]) return (-1);

  while (Nhi - Nlo > 4) {
    N = 0.5*(Nlo + Nhi);
    if (startTimes[N] < start) {
      Nlo = MAX(N, 0);
    } else {
      Nhi = MIN(N, Nmosaic - 1);
    }
  }
  // we now have : startTimes[Nlo] < start <= starTimes[Nhi]

  // find a matched mosaic starting from Nlo, or return -1
  for (N = Nlo; N <= Nhi; N++) { 
    if (startTimes[N] == start) return N;
  }
  return (-1);
}

void setMosaicCenters (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, NC, Nc;
  double R, D, Rmid, Dmid;
  double *Rc, *Dc;

  NC = 100;
  ALLOCATE (Rc, double, NC);
  ALLOCATE (Dc, double, NC);

  for (i = 0; i < Nmosaic; i++) {
    Nc = 0;
    Rmid = Dmid = NAN;
    for (j = 0; j < MosaicN_Image[i]; j++) {
      m = MosaicToImage[i][j];

      NX = image[m].NX;
      NY = image[m].NY;
      if (!XY_to_RD (&R, &D, 0.5*NX, 0.5*NY, &image[m].coords)) continue;
      R = ohana_normalize_angle_to_midpoint (R, 180.0);

      // Exclude images with crazy astrometry
      // XXX NOTE : this is gpc1-specific
      { 
	double dP1 = hypot(image[m].coords.pc1_1, image[m].coords.pc1_2);
	double dP2 = hypot(image[m].coords.pc2_1, image[m].coords.pc2_2);
	if (fabs(dP1 - 1.0) > 0.02) continue;
	if (fabs(dP2 - 1.0) > 0.02) continue;

	double X00, Y00, X10, Y10, X01, Y01;
	XY_to_LM (&X00, &Y00, 0.0, 0.0, &image[m].coords);
	XY_to_LM (&X10, &Y10, image[m].NX, 0.0, &image[m].coords);
	XY_to_LM (&X01, &Y01, 0.0, image[m].NY, &image[m].coords);
	double dS0 = hypot ((X00 - X10), (Y00 - Y10));
	double dS1 = hypot ((X00 - X01), (Y00 - Y01));
	if (dS0 > 6000) continue;
	if (dS1 > 6500) continue;
      }	

      Rc[Nc] = R;
      Dc[Nc] = D;
      Nc ++;
      if (Nc >= NC) {
	NC += 100;
	REALLOCATE (Rc, double, NC);
	REALLOCATE (Dc, double, NC);
      }
    }

    if (Nc > 0) {
      dsort (Rc, Nc);
      if (Rc[Nc-1] - Rc[0] > 180.0) {
	// in our list, Rc is in the range 0.0 to 360.0.  
	// any mosaic which is close to the 0.0, 360.0 boundary may have some on
	// one side or the other.  count how many have values more than Rc[0] + 180.
	// if more than half are at the large end, re-normalize to that range
	int Nbig = 0;
	for (j = 1; j < Nc; j++) {
	  if (Rc[j] - Rc[0] > 180.0) Nbig ++;
	}
	if (Nbig  > 0.5*Nc) {
	  for (j = 0; j < Nc; j++) {
	    Rc[j] = ohana_normalize_angle_to_midpoint (Rc[j], 360.0);
	  }
	  dsort (Rc, Nc);
	} else if (Nbig > 0) {
	  for (j = 0; j < Nc; j++) {
	    Rc[j] = ohana_normalize_angle_to_midpoint (Rc[j], 0.0);
	  }
	  dsort (Rc, Nc);
	}
      }
      dsort (Dc, Nc);

      Rmid = Rc[(int)(0.5*Nc)];
      Dmid = Dc[(int)(0.5*Nc)];
    }

    InitCoords (&mosaic[i].coords, "DEC--TAN");
    mosaic[i].coords.crval1 = Rmid;
    mosaic[i].coords.crval2 = Dmid;
    mosaic[i].coords.cdelt1 = 1.0 / 3600.0;
    mosaic[i].coords.cdelt2 = 1.0 / 3600.0;

    mosaic[i].Mcal   = 0.0;
    mosaic[i].dMcal  = 0.0;
    mosaic[i].Xm     = 0.0;
  }
  return;
}

void initMosaicGrid (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 Mcal, dMcal, Xm;

  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;
    Mcal = dMcal = Xm = 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

      Mcal  += image[m].Mcal;
      dMcal += image[m].dMcal;
      Xm    += image[m].Xm;

      image[m].Mcal   = 0.0;
      image[m].dMcal  = NAN;
      image[m].Xm     = NAN_S_SHORT;
    }
    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].Mcal   = Mcal / MosaicN_Image[i];
    mosaic[i].dMcal  = dMcal / MosaicN_Image[i];
    mosaic[i].Xm     = Xm / MosaicN_Image[i];
  }
  if (!USE_GRID) return;

  dXmax = MAX (dXmax, dX);
  dYmax = MAX (dYmax, dY);
  initGrid (dXmax, dYmax);
  return;
}

void setMcalFinal () {

  off_t i, j, im, Nimage;
  Image *image;

  if (!MOSAIC_ZEROPT) return;

  image = getimages (&Nimage, NULL);

  fprintf (stderr, "*** return Mcal from mosaic.Mcal to image.Mcal ***\n");

  // copy the mosaic results to the images.  set the mosaic Mcal to 0.0 since we have moved its
  // impact to the images
  for (i = 0; i < Nmosaic; i++) {
    for (j = 0; j < MosaicN_Image[i]; j++) {
      im = MosaicToImage[i][j];
      image[im].Mcal += mosaic[i].Mcal;
      image[im].dMcal = mosaic[i].dMcal;
      image[im].Xm = mosaic[i].Xm;
      image[im].ubercalDist = mosaic[i].ubercalDist;
      image[im].dMagSys = mosaic[i].dMsys;
      image[im].nFitPhotom = mosaic[i].nFitPhotom;
      image[im].flags |= (mosaic[i].flags & ID_IMAGE_PHOTOM_FEW);
      image[im].flags |= (mosaic[i].flags & ID_IMAGE_PHOTOM_POOR);
    }
    mosaic[i].Mcal = 0.0;
  }      
}

void initMosaicBins (Catalog *catalog, int Ncatalog, int doMosaicList) {

  off_t i, j;

  /* measure -> mosaic */
  if (!MOSAIC_ZEROPT) return;

  ALLOCATE (MeasureToMosaic, off_t *, Ncatalog);
  for (i = 0; i < Ncatalog; i++) {
    ALLOCATE (MeasureToMosaic[i], off_t, MAX (catalog[i].Nmeasure, 1));
    for (j = 0; j < catalog[i].Nmeasure; j++) MeasureToMosaic[i][j] = -1;
  }

  if (doMosaicList) {
    /* mosaic -> measure */
    ALLOCATE (N_onMosaic, off_t,   Nmosaic);
    ALLOCATE (N_ONMOSAIC, off_t,   Nmosaic);
    ALLOCATE (MosaicToCatalog, off_t *, Nmosaic);
    ALLOCATE (MosaicToMeasure, off_t *, Nmosaic);

    for (i = 0; i < Nmosaic; i++) {
      N_onMosaic[i] = 0;
      N_ONMOSAIC[i] = 100;
      ALLOCATE (MosaicToCatalog[i], off_t, N_ONMOSAIC[i]);
      ALLOCATE (MosaicToMeasure[i], off_t, N_ONMOSAIC[i]);
    }
  }
}

void freeMosaicBins (int Ncatalog, int doMosaicList) {

  off_t i;

  /* measure -> mosaic */
  if (!MOSAIC_ZEROPT) return;

  for (i = 0; i < Ncatalog; i++) {
    free (MeasureToMosaic[i]);
  }
  free (MeasureToMosaic);

  if (doMosaicList) {
    /* mosaic -> measure */
    for (i = 0; i < Nmosaic; i++) {
      free (MosaicToCatalog[i]);
      free (MosaicToMeasure[i]);
    }
    free (N_onMosaic);
    free (N_ONMOSAIC);
    free (MosaicToCatalog);
    free (MosaicToMeasure);
  }
}

int findMosaics (Catalog *catalog, int Ncatalog, int doMosaicList) {
  
  int i, ecode, found, Ns;
  off_t j;

  if (!MOSAIC_ZEROPT) return (FALSE);
  // if we are calibrating by mosaic, redefine myDet == on one of my mosaics

  int Nmatch = 0;
  for (i = 0; i < Ncatalog; i++) {
    for (j = 0; j < catalog[i].Nmeasure; j++) {
      catalog[i].measureT[j].myDet = FALSE; // a detetion is not mine until proven otherwise

      if (TimeSelect) {
	if (catalog[i].measureT[j].t < TSTART) continue;
	if (catalog[i].measureT[j].t > TSTOP) continue;
      }
      ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[j].photcode);
      found = FALSE;
      for (Ns = 0; !found && (Ns < Nphotcodes); Ns++) {
	if (ecode == photcodes[Ns][0].code) found = TRUE;
      }
      if (!found) continue;
      matchMosaics (catalog, j, i, doMosaicList);
      Nmatch ++;
    }
  }
  // fprintf (stderr, "matched %d detections to mosaics\n", Nmatch);
  return (TRUE);
}

void matchMosaics (Catalog *catalog, off_t meas, int cat, int doMosaicList) {

  off_t idx, ID, mosID;
  MeasureTiny *measure;

  measure = &catalog[cat].measureT[meas];

  ID = measure[0].imageID;
  idx = getImageByID (ID);
  if (idx == -1) {
    if (VERBOSE2) fprintf (stderr, "missed measurement "OFF_T_FMT", %d\n", meas, cat);
    return;
  }

  mosID = ImageToMosaic[idx];
  if (mosID < 0) {
    // Image *image = getimage(idx);
    // fprintf (stderr, "unmatched image %s\n", image[0].name);
    // skip measurements from simple chips (not matched to a mosaic by definition)
    return;
  }

  // test to check we got the right match:
  {
    Image *image = getimage(idx);
    // XXX we are now matching with just tzero.  be careful for cameras with drift
    // unsigned int imageStart = image[0].tzero - MAX(0.01*image[0].trate*image[0].NY, 1);
    unsigned int imageStart = image[0].tzero;
    if (imageStart != mosaic[mosID].start) {
      fprintf (stderr, "error in image to mosaic match\n");
      abort();
    }
  }

  // this measurement is on one of my mosaics, mark it as mine.
  catalog[cat].measureT[meas].myDet = TRUE;
  MeasureToMosaic[cat][meas] = mosID;

  if (doMosaicList) {
    MosaicToCatalog[mosID][N_onMosaic[mosID]] = cat;
    MosaicToMeasure[mosID][N_onMosaic[mosID]] = meas;
    N_onMosaic[mosID] ++;
    
    if (N_onMosaic[mosID] == N_ONMOSAIC[mosID]) {
      N_ONMOSAIC[mosID] += 100;
      REALLOCATE (MosaicToCatalog[mosID], off_t, N_ONMOSAIC[mosID]);
      REALLOCATE (MosaicToMeasure[mosID], off_t, N_ONMOSAIC[mosID]);
    }	
  }
  return;
}

float getMmos (off_t meas, int cat) {

  off_t i;
  float value;

  if (!MOSAIC_ZEROPT) return (0);

  // unassigned measurements belong to simple chips
  i = MeasureToMosaic[cat][meas];
  if (i == -1) return (0.0);

  if (mosaic[i].flags & IMAGE_BAD) return (NAN);  
  value = mosaic[i].Mcal;
  return (value);
}

int getMosaicFlags (off_t meas, int cat) {

  if (!MOSAIC_ZEROPT) return (0);

  // unassigned measurements belong to simple chips
  int i = MeasureToMosaic[cat][meas];
  if (i == -1) return (0);

  return (mosaic[i].flags);
}

typedef struct {
  int Nfew;
  int Nbad;
  int Ncal;
  int Ngrid;
  int Nrel;
  int Nsys;
  int Nskip;
  off_t Nmax;
  int PoorImages;
  double *list;
  double *dlist;
  double *Mlist;
  double *dMlist;
} SetMmosInfo;

enum {THREAD_RUN, THREAD_DONE};

typedef struct {
  int entry;
  int state;
  Catalog *catalog;
  Image *image;
  FlatCorrectionTable *flatcorr;
  SetMmosInfo info;
} ThreadInfo;

int setMmos_mosaic (Mosaic *mosaic, off_t Nmos, Image *image, Catalog *catalog, SetMmosInfo *info, FlatCorrectionTable *flatcorr);
void *setMmos_worker (void *data);
int setMmos_threaded (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr);

void SetMmosInfoInit (SetMmosInfo *info, off_t Nmax, int allocLists, int PoorImages) {
  info->Nfew = 0;
  info->Nbad = 0;
  info->Ncal = 0;
  info->Nrel = 0;
  info->Ngrid = 0;
  info->Nskip = 0;
  info->Nsys = 0;

  info->Nmax = Nmax;
  info->PoorImages = PoorImages;

  if (allocLists) {
    ALLOCATE (info->list, double, Nmax);
    ALLOCATE (info->dlist, double, Nmax);
    ALLOCATE (info->Mlist, double, Nmax);
    ALLOCATE (info->dMlist, double, Nmax);
  }
}

void SetMmosInfoFree (SetMmosInfo *info) {
  free (info->list);
  free (info->dlist);
  free (info->Mlist);
  free (info->dMlist);
}

void SetMmosInfoAccum (SetMmosInfo *summary, SetMmosInfo *results) {
  summary->Nfew  += results->Nfew ;
  summary->Nsys  += results->Nsys ;
  summary->Nbad  += results->Nbad ;
  summary->Ncal  += results->Ncal ;
  summary->Nrel  += results->Nrel ;
  summary->Ngrid += results->Ngrid;
  summary->Nskip += results->Nskip;
}

static int npass_output = 0;

// mutex to lock setMmos_worker operations 
static pthread_mutex_t setMmos_mutex = PTHREAD_MUTEX_INITIALIZER;
static int nextMosaic = 0;

// we have an array of mosaics (mosaic, Nmosaic).  we need to hand out mosaics one at a time to
// the worker threads as they need
off_t getNextMosaicForThread () {

  pthread_mutex_lock (&setMmos_mutex);
  if (nextMosaic >= Nmosaic) {
    pthread_mutex_unlock (&setMmos_mutex);
    return (-1);
  }
  int thisMosaic = nextMosaic;
  nextMosaic ++;

  pthread_mutex_unlock (&setMmos_mutex);
  return (thisMosaic);
}

int setMmos (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) {

  off_t i, N, Nmax;
  Image *image;

  if (!MOSAIC_ZEROPT) return (FALSE);
  if (FREEZE_MOSAICS) return (FALSE);

  // plots cannot be done in a threaded context (the plot commands collide)
  // so do not run setMmos in threaded mode if PLOTSTUFF is set
  if (NTHREADS && !PLOTSTUFF) {
    int status = setMmos_threaded (catalog, PoorImages, flatcorr);
    return status;
  }

  image = getimages (&N, NULL);

  fprintf (stderr, "limiting negative clouds to %f\n", CLOUD_TOLERANCE);

  if (PoorImages) {
    // XXX use bad stars and measurements for PoorImages? or not?
    // IMAGE_BAD = STAR_BAD = MEAS_BAD = 0;
    IMAGE_BAD = 0;
  }

  Nmax = 0;
  for (i = 0; i < Nmosaic; i++) {
    Nmax = MAX (Nmax, N_onMosaic[i]);
  }

  SetMmosInfo info;
  SetMmosInfoInit (&info, Nmax, TRUE, PoorImages);

  // int savePlotDelay = PLOTDELAY;
  // if (PLOTSTUFF) PLOTDELAY = 0.0;

  for (i = 0; i < Nmosaic; i++) {
    setMmos_mosaic (&mosaic[i], i, image, catalog, &info, flatcorr);
  }
  SetMmosInfoFree (&info);

  if (PLOTSTUFF) {
    fprintf (stdout, "press return\n"); 
    if (fscanf (stdin, "%*c") != 1) fprintf (stderr, "\n");
  }

  npass_output ++;

  fprintf (stderr, "%d mosaics marked having too few measurements (Nbad: %d, Ncal: %d, Ngrid: %d, Nrel: %d, Nsys: %d)\n", info.Nfew, info.Nbad, info.Ncal, info.Ngrid, info.Nrel, info.Nsys);

  if (PoorImages) {
    IMAGE_BAD = ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP;
    STAR_BAD  = ID_OBJ_POOR | ID_OBJ_FEW;
    MEAS_BAD  = ID_MEAS_NOCAL | ID_MEAS_POOR_PHOTOM | ID_MEAS_SKIP_PHOTOM | ID_MEAS_AREA;
  }
  return (TRUE);
}
  
// 'mosaic' is a pointer to the current mosaic of interest (Nmos)
int setMmos_mosaic (Mosaic *myMosaic, off_t Nmos, Image *image, Catalog *catalog, SetMmosInfo *info, FlatCorrectionTable *flatcorr) {

  off_t j, NimageReal;

  StatType stats;
  // liststats_setmode (&stats, STATMODE);
  liststats_setmode (&stats, "INNER_WTMEAN");

  double *list   = info->list;
  double *dlist  = info->dlist;
  double *Mlist  = info->Mlist;
  double *dMlist = info->dMlist;

  assert (Nmos >= 0);
  assert (Nmos < Nmosaic);

  // Image *imageReal = getimages (&NimageReal, NULL); returned pointer is not used
  getimages (&NimageReal, NULL);

  /* on PoorImages run, skip good images */
  if (info->PoorImages) {
    int bad = myMosaic[0].flags & (ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_SKIP);
    if (!bad) return TRUE;
  }      

  // UBERCAL image: if this is an ubercal image, set minUbercalDist to 0:
  // we optionally do not recalibrate images with UBERCAL zero points 
  if (myMosaic[0].flags & ID_IMAGE_PHOTOM_UBERCAL) {
    myMosaic[0].ubercalDist = 0;
    // propagate ubercalDist to the images
    for (j = 0; j < MosaicN_Image[Nmos]; j++) {
      off_t im = MosaicToImage[Nmos][j];
      assert (im < NimageReal);
      assert (im >= 0);
      image[im].ubercalDist = myMosaic[0].ubercalDist;
      // fprintf (stderr, "%d %d %d\n", (int) i, (int) im, image[im].ubercalDist);
    }
    if (KEEP_UBERCAL) return TRUE;
  }

  int Nsecfilt = GetPhotcodeNsecfilt ();

  int minUbercalDist = 1000;

  // do not modify the calibration for mosaics with partial images loaded (skipCal TRUE)
  if (myMosaic[0].skipCal) {
    info->Nskip ++;
    return TRUE;
  }

  int testImage = FALSE;
  // testImage |= (abs(myMosaic[0].start - 1324104046) < 10);
  // testImage |= (abs(myMosaic[0].start - 1324103823) < 10);
  // testImage |= (abs(myMosaic[0].start - 1323003245) < 10);
  // testImage |= (abs(myMosaic[0].start - 1323003069) < 10);
  // testImage |= (abs(myMosaic[0].start - 1323003125) < 10);
  // testImage |= (abs(myMosaic[0].start - 1323003300) < 10);
  // testImage |= (abs(myMosaic[0].start - 1323003365) < 10);
  // testImage |= (abs(myMosaic[0].start - 1323003191) < 10);
  // testImage |= (abs(myMosaic[0].start - 1323003014) < 10);
  // testImage |= (abs(myMosaic[0].start - 1323003484) < 10);
  // testImage |= (abs(myMosaic[0].start - 1323003419) < 10);
  // testImage |= (abs(myMosaic[0].start - 1323002949) < 10);

  FILE *fout = NULL;
  if (testImage) {
    char filename[64];
    snprintf (filename, 64, "test.%05d.%02d.dat", (int) Nmos, npass_output);
    fout = fopen (filename, "w");
  }

  // number of stars to measure the bright-end scatter
  int Nbright = 0;

  int N = 0;
  for (j = 0; j < N_onMosaic[Nmos]; j++) {
    float Msys, Mrel, Mcal, Mgrid, Mflat;
      
    off_t m = MosaicToMeasure[Nmos][j];
    off_t c = MosaicToCatalog[Nmos][j];
      
    if (fout) {
      Mcal  = getMcal  (m, c, flatcorr, catalog);
      Mgrid = getMgrid (m, c);
      Mrel  = getMrel  (catalog, m, c);
      Mflat = getMflat (m, c, flatcorr, catalog);

      off_t n = catalog[c].measureT[m].averef;
      Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF);

      float delta = Msys - Mrel - Mcal - Mgrid + Mflat;

      int isBad = (catalog[c].measureT[m].dbFlags & MEAS_BAD);

      fprintf (fout, "%f %f : %f %f %f %f %f  : %f %d\n", catalog[c].averageT[n].R, catalog[c].averageT[n].D, Msys, Mrel, Mcal, Mgrid, Mflat, delta, isBad);
    }

    if (catalog[c].measureT[m].dbFlags & MEAS_BAD) {
      info->Nbad ++;
      continue;
    }
    Mcal  = getMcal  (m, c, flatcorr, catalog);
    if (isnan(Mcal)) {
      info->Ncal++;
      continue;
    }
    Mgrid = getMgrid (m, c);
    if (isnan(Mgrid)) {
      info->Ngrid ++;
      continue;
    }
    Mrel  = getMrel  (catalog, m, c);
    if (isnan(Mrel)) {
      info->Nrel ++;
      continue;
    }
      
    // image.Mcal is not supposed to include the flat-field correction, so we need to
    // apply that offset as well here for this image (in other words, each detection is
    // being compared to the model, excluding the zero point, Mcal.  The model includes
    // the flat-correction.  NOTE the sign of Mflat (Image.Mcal = Measure.Mcal - Mflat)

    Mflat = getMflat (m, c, flatcorr, catalog);

    off_t n = catalog[c].measureT[m].averef;
    Msys = PhotSysTiny (&catalog[c].measureT[m], &catalog[c].averageT[n], &catalog[c].secfilt[n*Nsecfilt], MAG_CLASS_PSF);
    if (isnan(Msys)) {
      info->Nsys++;
      continue;
    }

    PhotCode *code = GetPhotcodebyCode (catalog[c].measureT[m].photcode);
    if (!code) goto skip;
    if (code->equiv < 1) goto skip;
    int Nsec = GetPhotcodeNsec (code->equiv);
    if (Nsec == -1) goto skip;
    minUbercalDist = MIN (catalog[c].secfilt[n*Nsecfilt + Nsec].ubercalDist, minUbercalDist);

  skip:
    assert (N < info->Nmax);
    assert (N >= 0);
    assert (Nbright < info->Nmax);
    assert (Nbright >= 0);

    list[N]  = Msys - Mrel - Mcal - Mgrid + Mflat;
    dlist[N] = MAX (catalog[c].measureT[m].dM, MIN_ERROR);
    if (catalog[c].measureT[m].dM < IMFIT_SYS_SIGMA_LIM) {
      Mlist[Nbright] = list[N];
      dMlist[Nbright] = dlist[N];
      Nbright ++;
    }
    N++;
  }
  /* N_onMosaic[Nmos] is all measurements, N is good measurements */

  if (fout) {
    fclose (fout);
  }

  /* too few good measurements or too many bad measurements (skip in PoorImages run) */
  if (!info->PoorImages) {
    int mark = (N < IMAGE_TOOFEW) || (N < IMAGE_GOOD_FRACTION*N_onMosaic[Nmos]);
    if (mark) {
      if (VERBOSE2) { fprintf (stderr, "marked mosaic %s ("OFF_T_FMT"), (%d < %d) || (%d < %f*"OFF_T_FMT")\n", image[MosaicToImage[Nmos][0]].name,  Nmos,  N, IMAGE_TOOFEW,  N, IMAGE_GOOD_FRACTION,  N_onMosaic[Nmos]); }
      myMosaic[0].flags |= ID_IMAGE_PHOTOM_FEW;
      info->Nfew ++;
      if (testImage) {
	fprintf (stderr, "NOTE: *** marked test image poor : %d %d %d***\n", (int) N, (int) IMAGE_TOOFEW, (int) (IMAGE_GOOD_FRACTION*N_onMosaic[Nmos]));
      }
    } else {
      myMosaic[0].flags &= ~ID_IMAGE_PHOTOM_FEW;
    }
  }

  liststats (list, dlist, NULL, N, &stats);
  if (VERBOSE2 && info->PoorImages) fprintf (stderr, "Mmos: %f %f %d %d\n", stats.mean, stats.sigma, stats.Nmeas, N);

  myMosaic[0].Mcal  = stats.mean;
  myMosaic[0].dMcal = stats.error;
  myMosaic[0].nFitPhotom = N;
  myMosaic[0].Xm    = 100.0*log10(stats.chisq);

  if (testImage) {
    fprintf (stderr, "test image %d (%d) %f %f %d ... ", (int) Nmos, myMosaic[0].start, stats.mean, stats.error, myMosaic[0].nFitPhotom);
  }

  if (PLOTSTUFF) {
    fprintf (stderr, "Mmos: %6.3f %6.3f +/- %6.3f %5d %5d | %s\n", stats.mean, stats.median, stats.sigma, stats.Nmeas, N, image[MosaicToImage[Nmos][0]].name);
    plot_setMcal (list, N, &stats, CLOUD_TOLERANCE);
  }

  // bright end scatter
  liststats (Mlist, dMlist, NULL, Nbright, &stats);
  myMosaic[0].dMsys = stats.sigma;

  if (myMosaic[0].Mcal < -CLOUD_TOLERANCE) {
    myMosaic[0].Mcal = 0.0;
  }

  if (testImage) {
    fprintf (stderr, "%f %f  :  %f\n", myMosaic[0].Mcal, myMosaic[0].dMsys, pow (10.0, 0.01*myMosaic[0].Xm));
  }

  // minUbercalDist calculated here is the min value for any star owned by this image
  // since this particular image is tied to that star, bump its distance by 1
  myMosaic[0].ubercalDist = minUbercalDist + 1;

  // propagate ubercalDist to the images
  for (j = 0; j < MosaicN_Image[Nmos]; j++) {
    off_t im = MosaicToImage[Nmos][j];
    image[im].ubercalDist = myMosaic[0].ubercalDist;
    // fprintf (stderr, "%d %d %d\n", (int) i, (int) im, image[im].ubercalDist);
  }
  
  return TRUE;
}
  
int setMmos_threaded (Catalog *catalog, int PoorImages, FlatCorrectionTable *flatcorr) {

  int i;
  off_t N;

  Image *image = getimages (&N, NULL);

  fprintf (stderr, "limiting negative clouds to %f\n", CLOUD_TOLERANCE);

  if (PoorImages) {
    // XXX use bad stars and measurements for PoorImages? or not?
    // IMAGE_BAD = STAR_BAD = MEAS_BAD = 0;
    IMAGE_BAD = 0;
  }

  off_t Nmax = 0;
  for (i = 0; i < Nmosaic; i++) {
    Nmax = MAX (Nmax, N_onMosaic[i]);
  }

  SetMmosInfo summary;
  SetMmosInfoInit (&summary, Nmax, FALSE, PoorImages);

  pthread_attr_t attr;
  pthread_attr_init (&attr);
  pthread_attr_setdetachstate (&attr, PTHREAD_CREATE_DETACHED);
  
  pthread_t *threads;
  ALLOCATE (threads, pthread_t, NTHREADS);

  ThreadInfo *threadinfo;
  ALLOCATE (threadinfo, ThreadInfo, NTHREADS);

  // each time this function is called, we cycle through the available mosaics
  // make sure we start at 0
  nextMosaic = 0;;

  // launch N worker threads
  for (i = 0; i < NTHREADS; i++) {
    threadinfo[i].entry = i;
    threadinfo[i].state = THREAD_RUN;
    threadinfo[i].catalog  =  catalog;
    threadinfo[i].image    =    image;
    threadinfo[i].flatcorr = flatcorr;
    SetMmosInfoInit (&threadinfo[i].info, Nmax, FALSE, PoorImages);
    pthread_create (&threads[i], NULL, setMmos_worker, &threadinfo[i]);
  }
  pthread_attr_destroy (&attr);

  // wait until all threads have finished
  while (1) {
    int allDone = TRUE;
    for (i = 0; i < NTHREADS; i++) {
      if (threadinfo[i].state == THREAD_RUN) allDone = FALSE;
    }
    if (allDone) {
      break;
    }
    usleep (500000);
  }

  // all threads are done, free the threads array and grab the info
  free (threads);
  
  // report stats & summary from the threads
  for (i = 0; i < NTHREADS; i++) {
    fprintf (stderr, "setMmos thread %d : %d mosaics marked having too few measurements (Nbad: %d, Ncal: %d, Ngrid: %d, Nrel: %d, Nsys: %d), %d partials skipped\n", 
	     i, 
	     threadinfo[i].info.Nfew, 
	     threadinfo[i].info.Nbad, 
	     threadinfo[i].info.Ncal, 
	     threadinfo[i].info.Ngrid, 
	     threadinfo[i].info.Nrel, 
	     threadinfo[i].info.Nsys,
      	     threadinfo[i].info.Nskip);
    SetMmosInfoAccum (&summary, &threadinfo[i].info);
  }
  fprintf (stderr, "total : %d mosaics marked having too few measurements (Nbad: %d, Ncal: %d, Ngrid: %d, Nrel: %d, Nsys: %d), %d partials skipped\n", 
	   summary.Nfew, 
	   summary.Nbad, 
	   summary.Ncal, 
	   summary.Ngrid, 
	   summary.Nrel, 
	   summary.Nsys, 
	   summary.Nskip);
  free (threadinfo);

  npass_output ++;

  if (PoorImages) {
    IMAGE_BAD = ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP;
    STAR_BAD  = ID_OBJ_POOR | ID_OBJ_FEW;
    MEAS_BAD  = ID_MEAS_NOCAL | ID_MEAS_POOR_PHOTOM | ID_MEAS_SKIP_PHOTOM | ID_MEAS_AREA;
  }
  return TRUE;
}

void *setMmos_worker (void *data) {

  ThreadInfo *threadinfo = data;

  SetMmosInfo results;
  SetMmosInfoInit (&results, threadinfo->info.Nmax, TRUE, threadinfo->info.PoorImages); // allocate list, dlist arrays here

  while (1) {

    off_t i = getNextMosaicForThread();
    if (i == -1) {
      threadinfo->state = THREAD_DONE;
      return NULL;
    }

    Catalog *catalog = threadinfo->catalog;
    FlatCorrectionTable *flatcorr = threadinfo->flatcorr;
    Image *image = threadinfo->image;

    setMmos_mosaic (&mosaic[i], i, image, catalog, &results, flatcorr);
    SetMmosInfoAccum (&threadinfo->info, &results);
  }

  SetMmosInfoFree (&results);
  return NULL;
}

// When we rationalize the images/mosaics, we are driving the negative cloud images back
// to 0.0.  At the same time, we make a guess to the effective impact on all other images,
// driven by the coupling of common stars. 
// XXX this function is no longer used because we force significantly negative clouds to 0.0
int rationalize_mosaics (Catalog *catalog, int Ncatalog) {

  double *imageOffset, **starOffset;
  int **starNcount, *seclist;
  int **Slist, *NSlist, *NSLIST;
  int i, j, k, m, nMos, Ns, found;

  off_t Nimage;
  Image *image;

  // set a test value for now
  float CLOUD_TOLERANCE = 0.025;

  if (!MOSAIC_ZEROPT) return (FALSE);
  if (FREEZE_MOSAICS) return (FALSE);

  image = getimages (&Nimage, NULL);

  ALLOCATE (imageOffset, double, Nmosaic);
  ALLOCATE ( Slist, int *, Nmosaic); // array of calibrated star indexes on this mosaic
  ALLOCATE (NSlist, int,   Nmosaic); // number of stars on mosaic
  ALLOCATE (NSLIST, int,   Nmosaic); // number of Slist entries allocated
  memset (Slist, 0, Nmosaic*sizeof(int *));

  // find the images / mosaics with negative clouds and save their offset
  for (i = 0; i < Nmosaic; i++) {
  
    NSlist[i] =   0;
    NSLIST[i] = 100;
    ALLOCATE (Slist[i], int, NSLIST[i]);

    imageOffset[i] = 0.0;

    if (VERBOSE2 && (fabs(mosaic[i].Mcal) < CLOUD_TOLERANCE)) {
      fprintf (stderr, "cloud-free: %s : %f\n", image[MosaicToImage[i][0]].name, mosaic[i].Mcal);
    } 
    if (VERBOSE2 && (mosaic[i].Mcal < -CLOUD_TOLERANCE)) {
      imageOffset[i] = -mosaic[i].Mcal; 
      // NOTE the negative sign: down below, we are going to add in the negative of Mcal
      // to this image, and the propagated mean values for other images
      fprintf (stderr, "anti-clouds: %s : %f\n", image[MosaicToImage[i][0]].name, mosaic[i].Mcal);
    } 
    if (VERBOSE2 && (mosaic[i].Mcal > CLOUD_TOLERANCE)) {
      fprintf (stderr, "cloudy    : %s : %f\n", image[MosaicToImage[i][0]].name, mosaic[i].Mcal);
    } 
  }

  int Nsecfilt = GetPhotcodeNsecfilt ();
  ALLOCATE (seclist, int, Nphotcodes);
  for (Ns = 0; Ns < Nphotcodes; Ns ++) { 
    int thisCode = photcodes[Ns][0].code;
    seclist[Ns] = GetPhotcodeNsec(thisCode);
  }

  // allocate an array for star offsets
  int Nstars = 0;
  int NSTARS = 1000;
  ALLOCATE (starOffset, double *, NSTARS);
  ALLOCATE (starNcount, int *,    NSTARS);
  memset (starOffset, 0, NSTARS*sizeof(double *));
  memset (starNcount, 0, NSTARS*sizeof(int *));

  // find the mean offset for each star
  for (i = 0; i < Ncatalog; i++) {
    for (j = 0; j < catalog[i].Naverage; j++) {
      ALLOCATE (starOffset[Nstars], double, Nphotcodes);
      ALLOCATE (starNcount[Nstars], int,    Nphotcodes);
      memset (starOffset[Nstars], 0, Nphotcodes*sizeof(double));
      memset (starNcount[Nstars], 0, Nphotcodes*sizeof(int));
      
      m = catalog[i].averageT[j].measureOffset;

      // determine the mosaic for each measurement
      for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) {

	// skip unused measurements
	if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue;

	// skip unused measurements
	int Nsec;
	int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode);
	found = FALSE;
	for (Ns = 0; !found && (Ns < Nphotcodes); Ns ++) { 
	  if (ecode == photcodes[Ns][0].code) {
	    found = TRUE;
	    break;
	  }
	}
	if (!found) continue;
	Nsec = seclist[Ns];

	// bad stars for this secfilt
	if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;  

	// skip REF measurements (not tied to an image)
	if (getImageEntry (m, i) < 0) continue;

	// find the source of this measurement (skip unassigned measurements)
	nMos = MeasureToMosaic[i][m];
	if (nMos == -1) continue;

	if (mosaic[nMos].photcode != ecode) {
	  fprintf (stderr, "*");
	}

	assert (Ns >= 0);
	assert (Ns < Nphotcodes);

	// accumulate the offsets from the negative cloud images (others have 0.0 value)
	starOffset[Nstars][Ns] += imageOffset[nMos];
	starNcount[Nstars][Ns] ++;
	
	// record the mosaic->star reference
	Slist[nMos][NSlist[nMos]] = Nstars;
	NSlist[nMos] ++;
	if (NSlist[nMos] == NSLIST[nMos]) {
	  NSLIST[nMos] += 100;
	  REALLOCATE (Slist[nMos], int, NSLIST[nMos]);
	}	  
      }
      Nstars ++;
      if (Nstars == NSTARS) {
	NSTARS += 1000;
	REALLOCATE (starOffset, double *, NSTARS);
	REALLOCATE (starNcount, int *,    NSTARS);
	memset (&starOffset[NSTARS-1000], 0, 1000*sizeof(double *));
	memset (&starNcount[NSTARS-1000], 0, 1000*sizeof(int *));
      }
    }
  }

  // find the mean offset of the images without negative clouds
  for (i = 0; i < Nmosaic; i++) {

    found = FALSE;
    for (Ns = 0; !found && (Ns < Nphotcodes); Ns ++) { 
      if (mosaic[i].photcode == photcodes[Ns][0].code) {
	found = TRUE;
	break;
      }
    }
    if (!found) {
      fprintf (stderr, "invalid photcode for mosaic?\n");
      abort();
    }

    // a negative cloud image (cloud: Mcal > 0; anti-clouds: Mcal < 0; imageOffset = -Mcal)
    if (imageOffset[i] > 0.0) continue;

    // we need to actually have cross-references to count
    if (NSlist[i] < 2) continue;

    float dM = 0.0;
    for (j = 0; j < NSlist[i]; j++) {
      k = Slist[i][j];
      if (starNcount[k][Ns] > 1) {
	dM += (starOffset[k][Ns] / starNcount[k][Ns]);
      }
    }
    imageOffset[i] = dM / NSlist[i];
    // fprintf (stderr, "adjust image: %s : (%f %d) : %f\n", image[MosaicToImage[i][0]].name, dM, NSlist[i], imageOffset[i]);
  }

  // for (i = 0; i < Nmosaic; i++) {
  //   fprintf (stderr, "correction: %s : %f\n", image[MosaicToImage[i][0]].name, imageOffset[i]);
  // } 

  // apply all offset values to the mosaics
  // find the images / mosaics with negative clouds and save their offset
  for (i = 0; i < Nmosaic; i++) {
    mosaic[i].Mcal += imageOffset[i];
  }

  for (i = 0; i < Nstars; i++) {
    free (starOffset[i]);
    free (starNcount[i]);
  }
  free (starOffset);
  free (starNcount);

  free (seclist);
  free (NSlist);
  free (NSLIST);
  for (i = 0; i < Nmosaic; i++){
    free (Slist[i]);
  }
  free (Slist);
  free (imageOffset);
  return (TRUE);
}

StatType statsMosaicM (Catalog *catalog) {
  OHANA_UNUSED_PARAM(catalog);

  off_t i;
  int n;
  double *list, *dlist;
  StatType stats;

  bzero (&stats, sizeof (StatType));
  if (!MOSAIC_ZEROPT) return (stats);

  ALLOCATE (list, double, Nmosaic);
  ALLOCATE (dlist, double, Nmosaic);

  n = 0;
  for (i = 0; i < Nmosaic; i++) {
    if (mosaic[i].flags & IMAGE_BAD) continue;
    if (mosaic[i].skipCal) continue;
    if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue;

    list[n] = mosaic[i].Mcal;
    dlist[n] = 1;
    n++;
  }

  liststats_setmode (&stats, "MEAN");

  liststats (list, dlist, NULL, n, &stats);
  free (list);
  free (dlist);
  return (stats);
}

StatType statsMosaicdM (Catalog *catalog) {
  OHANA_UNUSED_PARAM(catalog);

  off_t i, n;
  double *list, *dlist;
  StatType stats;

  bzero (&stats, sizeof (StatType));
  if (!MOSAIC_ZEROPT) return (stats);
  if (FREEZE_MOSAICS) return (stats); 

  ALLOCATE (list, double, Nmosaic);
  ALLOCATE (dlist, double, Nmosaic);

  n = 0;
  for (i = 0; i < Nmosaic; i++) {
    if (mosaic[i].flags & IMAGE_BAD) continue;
    if (mosaic[i].skipCal) continue;
    if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue;

    list[n] = mosaic[i].dMcal;
    dlist[n] = 1;
    n++;
  }

  liststats_setmode (&stats, "MEAN");

  liststats (list, dlist, NULL, n, &stats);
  free (list);
  free (dlist);
  return (stats);
}

StatType statsMosaicN (Catalog *catalog, FlatCorrectionTable *flatcorr) {

  off_t i, j, m, c, n, N;
  double *list, *dlist;
  float Mcal, Mgrid, Mrel;
  StatType stats;

  bzero (&stats, sizeof (StatType));
  if (!MOSAIC_ZEROPT) return (stats);
  if (FREEZE_MOSAICS) return (stats); 

  ALLOCATE (list, double, Nmosaic);
  ALLOCATE (dlist, double, Nmosaic);

  n = 0;
  for (i = 0; i < Nmosaic; i++) {
    if (mosaic[i].flags & IMAGE_BAD)  continue;
    if (mosaic[i].skipCal) continue;
    if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue;

    N = 0;
    for (j = 0; j < N_onMosaic[i]; j++) {

      m = MosaicToMeasure[i][j];
      c = MosaicToCatalog[i][j];

      Mcal = getMcal  (m, c, flatcorr, catalog);
      if (isnan(Mcal)) continue;
      Mgrid = getMgrid (m, c);
      if (isnan(Mgrid)) continue;
      Mrel = getMrel  (catalog, m, c);
      if (isnan(Mrel)) continue;
      N++;
    }
    list[n] = N;
    dlist[n] = 1;
    n++;
  }
  // fprintf (stderr, "Nmosaic: "OFF_T_FMT", n: "OFF_T_FMT"\n",  Nmosaic,  n);

  liststats_setmode (&stats, "MEAN");

  liststats (list, dlist, NULL, n, &stats);
  free (list);
  free (dlist);
  return (stats);
}

StatType statsMosaicX (Catalog *catalog) {
  OHANA_UNUSED_PARAM(catalog);

  off_t i, n;
  double *list, *dlist;
  StatType stats;

  bzero (&stats, sizeof (StatType));
  if (!MOSAIC_ZEROPT) return (stats);
  if (FREEZE_MOSAICS) return (stats); 

  ALLOCATE (list, double, Nmosaic);
  ALLOCATE (dlist, double, Nmosaic);

  n = 0;
  for (i = 0; i < Nmosaic; i++) {
    if (mosaic[i].flags & IMAGE_BAD) continue;
    if (mosaic[i].skipCal) continue;
    if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue;

    list[n]  = pow (10.0, 0.01*mosaic[i].Xm);
    dlist[n] = 1;
    n++;
  }

  liststats_setmode (&stats, "MEAN");

  liststats (list, dlist, NULL, n, &stats);
  free (list);
  free (dlist);
  return (stats);
}

/* mark mosaic if: abs(Mcal - <Mcal>) too large, dMcal too large */
void clean_mosaics () {

  off_t i, N, mark, Nmark, Nscatter, Noffset;
  double *mlist, *slist, *dlist;
  double MaxOffset, MedOffset, MaxScatter;

  if (!MOSAIC_ZEROPT) return;
  if (FREEZE_MOSAICS) return;

  if (VERBOSE) fprintf (stderr, "marking poor mosaics\n");

  ALLOCATE (mlist, double, Nmosaic);
  ALLOCATE (slist, double, Nmosaic);
  ALLOCATE (dlist, double, Nmosaic);

  for (i = N = 0; i < Nmosaic; i++) {
    if (mosaic[i].flags & IMAGE_BAD) continue;
    if (mosaic[i].skipCal) continue;
    mlist[N] = mosaic[i].Mcal;
    slist[N] = mosaic[i].dMcal;
    dlist[N] = 1;
    N++;
  }

  StatType stats;
  liststats_setmode (&stats, "MEAN");

  liststats (mlist, dlist, NULL, N, &stats);
  MaxOffset = MAX (IMAGE_OFFSET, 2*stats.sigma);
  MedOffset = stats.median;

  liststats (slist, dlist, NULL, N, &stats);
  MaxScatter = MAX (IMAGE_SCATTER, 2*stats.median);
  fprintf (stderr, "Mrel: %f, dMrel: %f, Max Scatter: %f, Max Offset: %f\n", MedOffset, stats.median, MaxScatter, MaxOffset);
  
  Nmark = Nscatter = Noffset = 0;
  for (i = 0; i < Nmosaic; i++) {
    // if we are keeping ubercal sacrosanct, then we should not be allowed to break them...
    if (KEEP_UBERCAL && (mosaic[i].flags & ID_IMAGE_PHOTOM_UBERCAL)) continue;

    if (mosaic[i].flags & (ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP)) continue;
    if (mosaic[i].skipCal) continue;

    mark = FALSE;
    if (mosaic[i].dMcal > MaxScatter) {
      mark = TRUE;
      Nscatter ++;
    }
    if (fabs(mosaic[i].Mcal - MedOffset) > MaxOffset) {
      mark = TRUE;
      Noffset ++;
    }
    if (mark) { 
      Nmark ++;
      mosaic[i].flags |= ID_IMAGE_PHOTOM_POOR;
    } else {
      mosaic[i].flags &= ~ID_IMAGE_PHOTOM_POOR;
    }
  }

  fprintf (stderr, OFF_T_FMT" mosaics marked poor ("OFF_T_FMT" scatter, "OFF_T_FMT" offset)\n",  Nmark, Nscatter, Noffset);

  free (mlist);
  free (slist);
  free (dlist);
}

void plot_mosaic_fields (Catalog *catalog) {

  off_t i, j, m, c, N, Nimage;
  double *xlist, *ylist;
  char string[64];
  Graphdata graphdata;

  if (!MOSAIC_ZEROPT) return;

  // Image *image = getimages (&Nimage, NULL); returned value ignored
  getimages (&Nimage, NULL);

  N = 0;
  for (i = 0; i < Nmosaic; i++) {
    N = MAX (N, N_onMosaic[i]);
  }

  ALLOCATE (xlist, double, N);
  ALLOCATE (ylist, double, N);

  for (i = 0; i < Nmosaic; i++) {
    N = 0;
    for (j = 0; j < N_onMosaic[i]; j++) {
      
      m = MosaicToMeasure[i][j];
      c = MosaicToCatalog[i][j];
      
      if (catalog[c].measureT[m].dbFlags & (ID_MEAS_AREA | ID_MEAS_NOCAL)) continue;

      // ave = catalog[c].measureT[m].averef;
      xlist[N] = catalog[c].measureT[m].R;
      ylist[N] = catalog[c].measureT[m].D;
      N++;
    }
  
    sprintf (string, "Mosaic "OFF_T_FMT,  i);
    plot_defaults (&graphdata);
    plot_list (&graphdata, xlist, ylist, N, string, NULL);
  }

  free (ylist);
  free (xlist);
}

void plot_mosaics () {

  off_t i, bin;
  double *xlist, *Mlist, *dlist;
  Graphdata graphdata;

  if (!MOSAIC_ZEROPT) return;
  if (FREEZE_MOSAICS) return;

  ALLOCATE (xlist, double, Nmosaic);
  ALLOCATE (dlist, double, Nmosaic);
  ALLOCATE (Mlist, double, Nmosaic);

  for (i = 0; i < Nmosaic; i++) {
    Mlist[i] = mosaic[i].Mcal;
    dlist[i] = mosaic[i].dMcal;
    xlist[i] = mosaic[i].secz;
  }

  plot_defaults (&graphdata);
  graphdata.xmin = 0.95;
  graphdata.xmax = 2.50;
  graphdata.ymin = PlotdMmin;
  graphdata.ymax = PlotdMmax;
  plot_list (&graphdata, xlist, Mlist, Nmosaic, "airmass vs Mcal", "%s.airmass.png", OUTROOT);
  plot_defaults (&graphdata);
  graphdata.size = 1.5;
  graphdata.ptype = 7;
  plot_list (&graphdata, Mlist, dlist, Nmosaic, "Mcal vs dMcal", "%s.MdM.png", OUTROOT);

# define NBIN 200
  REALLOCATE (xlist, double, NBIN);
  REALLOCATE (Mlist, double, NBIN);

  /**** dMcal histgram ****/
  for (i = 0; i < NBIN; i++) xlist[i] = 0.00005*i;
  bzero (Mlist, NBIN*sizeof(double));
  for (i = 0; i < Nmosaic; i++) {
    bin = mosaic[i].dMcal / 0.00005;
    bin = MAX (0, MIN (NBIN - 1, bin));
    Mlist[bin] += 1.0;
  }
  plot_defaults (&graphdata);
  graphdata.style = 1;
  plot_list (&graphdata, xlist, Mlist, NBIN, "dMcal hist", "%s.dMcalhist.png", OUTROOT);

  free (dlist);
  free (xlist);
  free (Mlist);
}

off_t *SelectRefMosaic (Mosaic **refmosaic, off_t *Nimage) {

  off_t i, Imax, Nmax;

  Imax = 0;
  Nmax = MosaicN_Image[0];
  for (i = 0; i < Nmosaic; i++) {
    if (MosaicN_Image[i] > Nmax) {
      Imax = i;
      Nmax = MosaicN_Image[i];
    }
  }

  *refmosaic = &mosaic[Imax];
  *Nimage = Nmax;
  return (MosaicToImage[Imax]);
}
