# include "relastro.h"

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

// list of all images associated with a mosaic
static int   *Nmosaic_own_images; // number of images for this mosaic
static int   *Amosaic_own_images; // size of allocated array
static int   **mosaic_own_images; // array of arrays: mosaic -> images

// list of mosaic associated with each image  
static int    Nmosaic_for_images; // number of images (for internal checks)
static int    *mosaic_for_images; // array of: image -> mosaic

Mosaic *getmosaics (int *N) {
  *N = Nmosaic;
  return (mosaic);
}

// find mosaic frames (unique time periods & photcode name matches mosaic) 
void initMosaics (Image *image, int Nimage) {

  int i, j, found, NMOSAIC;
  unsigned int start, stop;

  Nmosaic = 0;
  NMOSAIC = 10;
  ALLOCATE (mosaic, Mosaic, NMOSAIC);

  ALLOCATE (Nmosaic_own_images, int, NMOSAIC);
  ALLOCATE (Amosaic_own_images, int, NMOSAIC);
  ALLOCATE (mosaic_own_images, int *, NMOSAIC);

  /* find the mosaic images (coords.ctype = DIS); generate list of unique mosaics */
  for (i = 0; i < Nimage; i++) {
    if (strcmp(&image[i].coords.ctype[4], "-DIS")) continue;

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

    /* a new mosaic, define ranges */
    mosaic[Nmosaic].start = start;
    mosaic[Nmosaic].stop  = stop;
    mosaic[Nmosaic].Mcal  = 0.0;
    mosaic[Nmosaic].dMcal = 0.0;
    mosaic[Nmosaic].Xm    = 0.0;
    mosaic[Nmosaic].code  = image[i].code;
    mosaic[Nmosaic].secz  = image[i].secz;
    mosaic[Nmosaic].coords = image[i].coords;

    // init the mosaic_own_images array data
    Nmosaic_own_images[Nmosaic] = 0;
    Amosaic_own_images[Nmosaic] = 10;
    ALLOCATE (mosaic_own_images[Nmosaic], int, Amosaic_own_images[Nmosaic]);

    Nmosaic ++;
    if (Nmosaic == NMOSAIC) {
      NMOSAIC += 10;
      REALLOCATE (mosaic, Mosaic, NMOSAIC);
      REALLOCATE (mosaic_own_images, int *, NMOSAIC);
      REALLOCATE (Nmosaic_own_images, int, NMOSAIC);
      REALLOCATE (Amosaic_own_images, int, NMOSAIC);
    }
  }

  // array to store image->mosaic index
  Nmosaic_for_images = Nimage;
  ALLOCATE (mosaic_for_images, int, Nmosaic_for_images);

  /* now assign the WRP images to these mosaics */
  for (i = 0; i < Nimage; i++) {
    mosaic_for_images[i] = -1; // default value for no mosaic found

    if (strcmp(&image[i].coords.ctype[4], "-WRP")) continue;

    /* set image time range */
    start = image[i].tzero - MAX(0.05*image[i].trate*image[i].NY, 1);
    stop  = image[i].tzero + MAX(1.05*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;
      break;
    }
    /* if no matching mosaic exists, skip this image */
    if (!found) continue;

    // mosaic corresponding to this image
    mosaic_for_images[i] = j;

    // add image to mosaic_own_image list 
    mosaic_own_images[j][Nmosaic_own_images[j]] = i;
    Nmosaic_own_images[j] ++;
    if (Nmosaic_own_images[j] == Amosaic_own_images[j]) {
      Amosaic_own_images[j] += 10;
      REALLOCATE (mosaic_own_images[j], int, Amosaic_own_images[j]);
    }
    assert (Nmosaic_own_images[j] < Amosaic_own_images[j]);
  }

  return;
}

// return StarData values for detections in the specified image, converting coordinates from the
// chip positions: X,Y -> L,M -> P,Q -> R,D
StarData *getMosaicRaw (Catalog *catalog, int Ncatalog, int mos, int *Nstars) {

  int i, j, im, Nraw, Nnew;
  StarData *raw, *new;

  Nraw = 0;
  ALLOCATE (raw, StarData, 1);

  // loop over the images owned by this mosaic
  for (i = 0; i < Nmosaic_own_images[mos]; i++) {

    im = mosaic_own_images[mos][i];
    
    // retrieve stars for this chip, applying chip and mosaic astrometry
    // this function does the reverse-lookup for the mosaic corresponding to this image
    new = getImageRaw (catalog, Ncatalog, im, &Nnew, MODE_MOSAIC);
    
    // merge new and raw
    REALLOCATE (raw, StarData, Nraw + Nnew);
    for (j = 0; j < Nnew; j++) {
      raw[Nraw+j] = new[j];
    }
    Nraw += Nnew;

    free (new);
  }

  *Nstars = Nraw;
  return (raw);
}

// return StarData values for averages positions in the specified image, converting coordinates from
// the sky positions: R,D -> P,Q -> L,M -> X,Y
StarData *getMosaicRef (Catalog *catalog, int Ncatalog, int mos, int *Nstars) {

  int i, j, im, Nref, Nnew;
  StarData *ref, *new;
  
  Nref = 0;
  ALLOCATE (ref, StarData, 1);

  for (i = 0; i < Nmosaic_own_images[mos]; i++) {

    im = mosaic_own_images[mos][i];
    
    // retrieve stars for this chip, applying chip and mosaic astrometry
    // this function does the reverse-lookup for the mosaic corresponding to this image
    new = getImageRef (catalog, Ncatalog, im, &Nnew, MODE_MOSAIC);
    
    // merge new and ref
    REALLOCATE (ref, StarData, Nref + Nnew);
    for (j = 0; j < Nnew; j++) {
      ref[Nref+j] = new[j];
    }
    Nref += Nnew;

    free (new);
  }

  *Nstars = Nref;
  return (ref);
}

Mosaic *getMosaicForImage (int im) {

  int mos;

  if (im < 0) abort();
  if (im >= Nmosaic_for_images) abort();

  // search for the mosaic that 
  mos = mosaic_for_images[im];
  if (mos < 0) return NULL;

  return &mosaic[mos];
}
